CN113160052B - 基于非均匀精度的近海养殖区域图像拼接方法 - Google Patents

Abstract

本发明公开了一种基于非均匀精度的近海养殖区域图像拼接方法，该方法，步骤如下包括：图像预处理，完成图形匀光操作与网箱区域掩膜获取；图像子区域划分，根据图形重叠区域的特点将图像划分为多个子区域；子区域优化拼接，使用基于特征点的DLT算法和基于光度对准的网格优化算法，结合网箱掩膜搜索缝合线，对子区域进行高精度拼接；全景优化拼接，利用位姿信息将多个子区域投影到全景平面，实现子区域间的拼接，得到全景图像。本发明针对近海养殖区域的特点，在子区域内部使用高精度拼接算法，子区域之间使用低精度拼接算法，实现了针对近海养殖区域的自然的拼接。

Description

基于非均匀精度的近海养殖区域图像拼接方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种基于非均匀精度的近海养殖区域图像拼接方法。

背景技术

在传统的近海养殖产业当中，对于养殖情况的观测通常由人工巡视来实现。但是，由于养殖设施分布分散，需要通过船只来到达各个养殖设施，费时费力。随着科学技术的快速发展，无人机已经被广泛的用于生产生活当中。利用无人机对养殖区域进行低空航拍，并将航拍图像拼接成全景图来辅助工作人员对养殖区域进行观测，可以帮助工作人员方便快速地了解养殖区域的水产品生长情况，极大地减轻了工作负担，提高了工作效率。

目前，常用的航拍图像拼接方法有两种，一种是基于特征点的拼接方法，另一种是基于位姿的方法。基于特征点的方法通过检测两张图像中的特征点从而实现对图像之间的单应矩阵的求取。单应矩阵可以很好的反映出两张图像之间的像素点的关系，通过单应矩阵可以很方便的将目标图像上的像素投影到参考图像中，实现两张图像的拼接。常用的基于特征点的方法主要有SIFT(Scale Invariant Feature Transform)、SURF(Speeded UpRobust Features)和ORB(Oriented Fast and Rotated Brief)。SIFT经过计算后对于每一个特征点会生成一个128维的特征向量，然后会根据特征向量之间的距离来进行特征点的匹配。这种检测方法检测出的特征点具有很高的精度，而且对于图像的光照、尺度和旋转具有很强的稳健性，因此可以获得很优秀的拼接效果。但是该方法在计算特征向量时需要进行高斯金字塔与拉普拉斯金字塔的计算，所以速度较慢。SURF是一种基于SIFT的改进方法，它使用盒式滤波器等方式提高了运算速度，而且特征向量的维度也有所减少，减少了匹配的计算量。ORB结合了FAST(Features from Accelerated Segment Test)角点检测算法与BRIEF(Binary Robust Independent Elementary Features)描述子来实现对特征点的检测，ORB相对于另外两种算法速度更快，但是误匹配率较高。这些基于特征点的方法在对重叠区域为固定场景时的拼接会有很好的效果，但是在重叠区域为水域时，往往难以实现有效的拼接。因为水面的波动会造成特征点的特征描述子的值发生改变，使得特征点难以实现匹配。

另一种方法是基于位姿的拼接方法，该方法利用相机的投影原理来获取单应矩阵实现图像的拼接，这种方法的拼接速度很快，计算简单，而且不受图像纹理、内容因素的影响，除此之外，还能够很好的记录拍摄位置等信息，但是无人机搭载的传感器精度往往不是很高，在拼接时有较大的误差。由于水域对拼接误差的容忍度较高，基于位姿的方法在对水域拼接时可以得到较为自然的拼接效果。但是，在对非水域，有着较复杂纹理的养殖设施拼接之时会有较为明显的错位。

由于近海养殖区域拼接场景中存在大量水域，且在水面上有很多呈块状离散分布的养殖设施。这种情况下，基于特征点的方式往往难以成功的实现拼接，而基于位姿的方法会在水面养殖设施的拼接时出现较大的错位。除此之外，传统的缝合线搜索方法主要基于图像的纹理，梯度等来计算拼接对齐的能量，从而保证缝合边界的自然性。若将其应用到对网箱区域的拼接时，缝合线往往会穿过网箱，鱼类的游动使得最终拼接结果偏离真实情况。因此，缝合线应该在保证缝合边界自然的同时尽量避开网箱区域。

发明内容

本发明的目的是为了解决现有拼接方法中在此水域图像场景的缺陷，提供一种基于非均匀精度的近海养殖区域图像拼接方法。该方法首先将近海养殖区域图像分为多个可以用特征点匹配的子区域。在子区域内使用特征点实现子区域拼接的粗对齐，然后使用基于光度对准的网格变换进行精对齐，实现子区域内部的高精度拼接。进一步结合深度学习u-net模型来检测出网箱区域，u-net(Ronneberger O,Fischer P,Brox T.U-net:Convolutional networks for biomedical image segmentation[C]//InternationalConference on Medical image computing and computer-assistedintervention.Springer,Cham,2015:234-241.)是在2015年提出的一种成熟的深度学习模型，分割精度高，效果好，此操作可以使得缝合线搜索时避开网箱区域，保护了网箱结构的完整性。由于不同子区域之间相交区域多为水，对位姿带来的错位敏感度较低。因此，在不同子区域拼接时，使用基于位姿的方法来进行低精度的拼接。

本发明的目的可以通过采取如下技术方案达到：

一种基于非均匀精度的近海养殖区域图像拼接方法，所述图像拼接方法包括如下步骤：

S1、利用无人机对养殖区域进行低空航拍得到航拍图像，对航拍图像执行预处理操作，调整航拍图像的亮度并获取网箱区域的掩膜；

S2、将经过预处理的航拍图像分成多个可以根据特征点进行拼接的子区域，所述子区域内部的图像重叠区域为浮于水面的养殖设施或者船只，所述子区域之间的重叠区域为水；

S3、对子区域进行高精度的内部拼接，利用特征点和基于光度对准的网格变换进行精准对齐，并结合步骤S1预处理中获取的网箱掩膜得到绕开网箱的缝合线，得到子区域的区域图，经过S3步骤的处理，获得的局部区域图像不仅精度高，而且很好的保证了网箱的完整性，降低了鱼类游动带来的负面效果；

S4、将所有子区域的区域图利用位姿获取的单应矩阵投影到全景平面实现全景的低精度拼接，获取全景图，拼接出的全景图不仅保证了对于养殖设施等水面目标的高精度拼接，而且在全局上也保证了自然的观感。

进一步地，所述步骤S1中对航拍图像执行预处理操作包括：对航拍图像进行匀光处理来调整亮度与通过u-net深度学习模型获取网箱掩膜。

进一步地，所述步骤S2过程如下：

S201、使用SIFT算法来检测每张航拍图像中的特征点，并创建一个集合，并第一张航拍图像放入其中，作为第一个子区域；

S202、从第二张航拍图像开始，依次对每一张未划分子区域的航拍图像执行步骤S203；

S203、针对当前未划分子区域的航拍图像，对每一张划分好子区域的航拍图像，执行步骤S204，如果当前未划分子区域的航拍图像对当前所有划分好子区域的航拍图像都无法完成区域匹配，则新建一个集合，并将当前未划分子区域的航拍图像放入新建的集合中，作为一个新的子区域，然后执行S202；

S204、利用位姿计算当前未划分子区域的航拍图像与当前划分好子区域的航拍图像之间的单应矩阵，进一步判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的相交情况，如果相交则执行S205，否则执行步骤S203，此操作可以大大降低特征点匹配的工作量，提高区域划分的速度；

S205、使用特征点匹配判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的匹配情况，判定所属的子区域，如果成功匹配则将当前航拍图像放入当前划分好子区域的航拍图像所属子区域，并执行S202，否则执行S203，其中对于特征点匹配步骤如下：

S20501、使用最佳优先搜索(Best-Bin-First,BBF)算法进行两张航拍图像之间的特征点粗匹配，得到匹配点对的数量n_f；

S20502、使用随机抽样一致性(Random Sample Consensus,RANSAC)算法进行特征点的精匹配，得到匹配点对的数量n,如果满足下式：

n>8+0.3*n_f 公式(1)

则认为其匹配成功。

进一步地，所述步骤S204中使用位姿判断两张航拍图像是否相交的过程如下：

首先建立辅助坐标系，辅助平面与辅助图像，以第一张航拍图像的拍摄位置为坐标原点，以正北方向为Y轴正方向，正西方向为X轴正方向，正地方向为Z轴正方向建立辅助坐标系，辅助平面即为大地平面，辅助图像为相机在辅助坐标系原点，相机坐标系三轴与辅助坐标系三轴重合时，大地平面在相机中所形成的图像；

若第一张航拍图像拍摄位置的坐标为

x₀、y₀、z₀为在辅助坐标系X轴、Y轴、Z轴上的值；旋转为

θ₀、ψ₀为第一张图像拍摄姿态绕辅助坐标系的X轴、Y轴、Z轴的旋转角度；第二张航拍图像拍摄位置的坐标为

三轴的旋转为

x₁、y₁、z₁为在辅助坐标系X轴、Y轴、Z轴上的值；

θ₁、ψ₁为第二张图像拍摄姿态绕辅助坐标系的X轴、Y轴、Z轴的旋转角度，第一张航拍图像拍摄位置与相对于参考坐标系原点的位移矩阵t为:

第一张航拍图像拍摄位姿与相对于辅助坐标系的旋转矩阵R为：

其中：

由相机的成像原理可知，对于空间中的一点，存在如下关系：

其中K为相机的内参矩阵，可由相机标定获得，t，R分别为公式(2)与公式(3)位移矩阵与旋转矩阵，

为空间中的一点，X′、Y′、Z′为在辅助坐标系X轴、Y轴、Z轴上的值，

为空间中该点在相机成像平面中的像素位置，u为像素的横坐标，v为像素的纵坐标；

对于第一张航拍图像到辅助图像，变换关系可由下式表示：

I为单位矩阵，K^-1为K的逆矩阵，d为参考坐标系原点距离参考平面的距离，n^T为参考平面的法向量，

为辅助图像中像素点的位置，u_f、v_f为辅助图像中像素点的横坐标、纵坐标，

为

投影到第一张航拍图像中的像素点的位置，u₀、v₀为第一张航拍图像中像素点的横坐标、纵坐标；其中，为了简化描述，可得下式：

其中H为辅助图像到第一张航拍图像的变换矩阵，而

H₁＝H^-1 公式(10)

为第一张航拍图像到辅助图像之间的变换矩阵,H^-1为H的逆；

同理可得到第二张航拍图像到辅助图像之间的变换矩阵H₂，接着获取两张航拍图像之间的单应矩阵H_t；然后，以其中一张航拍图像为参考图像，将另一张航拍图像投影到参考图像之上，并查看投影后的四个顶点是否在参考图像范围内，如果四个顶点中有任意一个顶点在参考图像范围内，则认为两张航拍图像相交。

进一步地，所述步骤S3过程如下：

S301、利用特征点匹配来获取每两张航拍图像之间的匹配关系；

S302、利用最大生成树计算出参考图像以及每张航拍图像的拼接顺序，进一步计算出每张航拍图像到参考图像的单应矩阵，并将航拍图像投影到参考图像平面上；

S303、利用基于光度对准的网格变换处理每一张已投影到参考图像平面的航拍图像；

S304、利用对齐质量与网箱掩膜来寻找最佳缝合线，并缝合子区域图像。

进一步地，所述步骤S303中循环遍历每张航拍图像，每一次循环都固定所有与当前航拍图像相交的航拍图像，并对当前航拍图像做基于光度对准的网格优化，优化过程如下：

S30301、网格化图像并对图像进行抽样，具体如下：

因为在航拍时拍摄的区域并不是完全的平面，而且基于特征点的配准也存在误差，因此会对对齐效果带来负面的影响。因此，需要对其进行网格优化。在优化时，将当前遍历到的航拍图像所在区域划分为多个网格，网格的顶点为vertices,并对当前航拍图像按距离进行像素点抽样，抽样后的点表示为sample_point,每一个抽样点都可以利用周围的4个网格顶点的双线性插值来表示，具体如下：

N_s为该抽样点周围的网格顶点集合，w_j为第j个顶点的权重，vertices_j为第j个顶点的位置；

将当前遍历到的航拍图像进行直线检测，对每条直线进行像素点抽样操作，抽样点为line_point,并对抽样点做双线性插值：

其中line_point为直线中的抽样点,N′_l为每个抽样点周围的4个网格顶点集合，w_k为第k个网格顶点的权重，vertices_k为第k个网格顶点的位置；

S30302、求取优化能量，获取优化后网格顶点的位置，对于当前优化的航拍图像与任意一张与其相交的航拍图像，当前优化的航拍图像为目标图像，相交的航拍图像为参考图像，能量定义如下：

第一项是光度对准约束项，其能量E_pm如下所示：

E_pm＝∑_i‖I_t(p_i+τ(p_i))-I_S(p_i)‖² 公式(13)

其中I_t为目标图像的灰度图，I_S为参考图像的灰度图，p_i为重叠区域图像中第i个采样点，τ(p_i)为第i个采样点的偏移量；

第二项是梯度约束项，其能量E_pg如下所示：

E_pg＝∑_i′‖G_t(p′_i′+τ(p′_i′))-G_S(′_i′)‖² 公式(14)

其中G_t为目标图像的灰度图的梯度图，G_S为参考图像的灰度图的梯度图，p′_i为重叠区域图像的第i′个采样点，τ(p′_i)为第i′个采样点的偏移量；

其中，在公式(13)与公式(14)中，p_i+τ(p_i)、p_i、p′_i′、p′_i′+τ(p′_i′)都通过公式(11)的插值方式表示；

第三项是三角形约束项，其能量E_s如公式(15)所示：

其中

为变换前的第i_tr个三角形的顶点，

为优化后的第i_tr个三角形的顶点，N_tr为三角形的数量，其中三角形由划分的网格顶点获得；

与

是通过第i_tr个变换前的三角形顶点根据公式(16)计算得出的变换参数，

第四项是直线约束项，其能量E_c为：

其中

第li条线段的起点，

为第li条线段的终点，

为第li条线段上第lj个抽样点，

均由公式(12)的方式表示，且u_l为当前点在线段上的相对比例位置，N_l为直线的数量，

为第li条线段抽样点的数量；

最终的能量E为：

E＝αE_pm+βE_pg+γE_s+ηE_c 公式(18)

其中α、β、γ、η为能量系数，对于多图联合优化，针对当前航拍图像，与所有的相交航拍图像求取公式(18)所述的能量项，并作网格优化；最终最小化能量E_fe，具体为：

其中np为与当前航拍图像相交的航拍图像数量，E_je为公式(18)表述的与第j_e张匹配图优化的能量项，最小化公式(19)即可得到优化后的网格顶点位置；

每一次优化后重新映射航拍图像并判断当前航拍图像的网格顶点的位置变化量，如果变化量小于事先指定的优化阈值，将该航拍图像移出优化列表，当优化列表为空时，停止优化。经过优化后，图像之间的对齐质量会得到进一步的提升，从而可以搜索到更为自然的缝合线。

进一步地，所述步骤S304过程如下：

利用对齐质量与网箱掩膜来寻找最佳缝合线，对网箱区域依照网箱掩膜施加能量值，使得缝合线绕开网箱，对于相邻的两张优化后的航拍图像，缝合线的能量E_seam可由下式表示：

其中E_seam为缝合线的总能量，p′_s为航拍图像中的像素；N_o为两张航拍图像重合区域相邻像素点的集合，p_s、q_s为N_o中的任意两个相邻像素点；N_t为两张航拍图像重合区域中的目标掩膜区域相邻像素点的集合,p_t、q_t为Nx中的任意一对相邻像素点，α′、β′为能量系数，能量函数各部分组成部分定义如下：

为非重合区域的能量，其中λ_m为一个惩罚值，

为像素的标签，确定像素是来自第一张航拍图像还是第二张航拍图像，这一项对重叠区域不做处理，但是对非重叠区域施加一个能量惩罚，

E_a为相邻的两张优化后的航拍图像对齐质量能量，

为p_s、q_s两个像素点的标签，即第一张航拍图像或第二张航拍图像；

为p_s点在

图中的像素值，

为p_s点在

图中的像素值,

为q_s点在

图中的像素值,

为q_s点在

图中的像素值，E_t为网箱区域根据网箱掩膜赋予的能量，mask_OE为对于网箱区域惩罚的能量值；结合这三项能量利用马尔科夫随机场获取最佳的缝合线。相对于传统的单纯基于图像特征缝合线搜索方法，在此基础上搜索到的缝合线，不仅可以在一定程度上保证了缝合的自然效果，而且，有效的降低了缝合线穿过网箱的概率，方便了对于网箱内部鱼类的观测。

找到最佳缝合线后，缝合航拍图像，得到子区域的航拍图像拼接完成图，也即该子区域的子区域图。

进一步地，所述步骤S4过程如下：

S401、计算全景平面的坐标系以及每个子区域的区域图到全景平面的单应矩阵；

S402、修正子区域的区域图的坐标范围；

S403、将子区域的区域图投影到全景平面获取全景图。

本发明相对于现有技术具有如下的优点及效果：

1)本发明可以针对重叠区域的特点选择拼接方法，在子区域内部拼接时，由于相交区域存在养殖设施等水面目标，因此采用基于特征点匹配结合网格优化来进行高精度的子区域内部拼接。子区域之间的相交区域由于绝大部分为水，特征点难以成功实现匹配，且对误差的容忍度较高，因此采用低精度的基于姿态的方式来进行拼接。最终在保持水面目标拼接精度高的同时实现了全局自然的拼接。

2)本发明在缝合线搜索时结合u-net检测的网箱掩膜与图像对齐误差来计算能量，在进行子区域内部拼接时，使得缝合线尽量绕开了网箱所在的区域，提升了养殖网箱区域在图像缝合时的完整度。

附图说明

图1是本发明实施例公开的基于非均匀精度的近海养殖区域图像拼接方法的工作流程图；

图2是本发明实施例中进行图像局部分类的流程图；

图3是本发明实施例中进行局部拼接的流程图；

图4是本发明实施例中进行全景拼接的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

本实例公开了一种基于非均匀精度的近海养殖区域图像拼接方法，主要结合位姿与特征点对图像进行先子区域内部后子区域之间的拼接，如图1所示，具体流程如下：

S1、利用无人机对养殖区域进行低空航拍得到航拍图像，对航拍图像执行预处理操作，调整航拍图像的亮度并获取网箱区域的掩膜；

对于图像的预处理，这里包括航拍图像匀光与网箱掩膜获取。由于相机在拍照时的曝光不一致，使得最终获取的航拍图像的亮度不一致。航拍图像的亮度不仅会影响最终拼接结果的观感，也会对基于光度对准的网格优化带来不利的影响。因此，首先要对航拍图像做匀光处理。在处理时，以第一张航拍图像为参考图像，其他的航拍图像为变换图像来进行处理，其中变换图像与参考图像的关系如下：

其中I²是变换图像，I¹是参考图像，I_i″j″是输出图像的第i″行j″列的像素，

是输出图像的第i″行j″列的像素，

σ_I1、σ_I2分别为对应参考图像与变换图像的均值与方差。

其中：

这里m′,n′为当前图像的行数与列数，I_{i″′j″′}为当前图像的第i″′行j″′列的像素。

亮度调整完毕后，使用u-net来对航拍图像进行分割，获取网箱区域的网箱掩膜。u-net(Ronneberger O,Fischer P,Brox T.U-net:Convolutional networks forbiomedical image segmentation[C]//International Conference on Medical imagecomputing and computer-assisted intervention.Springer,Cham,2015:234-241.)是一个2015年提出的图像分割模型，此模型分割精度高，所需数据量少。分割过程主要如下：1)首先搭建u-net深度学习模型；2)标注航拍图像中的目标，作为训练模型的数据集，这里的标注目标为网箱；3)训练模型；4)检测网箱掩膜，训练完毕后，模型即可通过输入的航拍图像预测出其中的网箱区域，针对网箱区域，再根据灰度阈值来分离出网箱的框架与水域，得到不同目标的区域的位置。

S2、将经过预处理的航拍图像分成多个可以根据特征点进行拼接的子区域，所述子区域内部的航拍图像重叠区域为浮于水面的养殖设施或者船只，所述子区域之间的重叠区域为水，其流程如图2所示，具体如下：

S201、使用SIFT算法来检测每张航拍图像中的特征点，创建一个集合，将第一张航拍图像放入其中，作为第一个子区域；

S202、从第二张航拍图像开始，依次对每一张未划分子区域的航拍图像执行步骤S203；

S203、针对当前未划分子区域的航拍图像，对每一张划分好子区域的航拍图像，执行步骤S204，如果当前未划分子区域的航拍图像对当前所有划分好子区域的航拍图像都无法完成区域匹配，则新建一个集合，并将当前未划分子区域的航拍图像放入新建的集合中，作为一个新的子区域，然后执行S202；

S204、利用位姿计算当前未划分子区域的航拍图像与当前划分好子区域的航拍图像之间的单应矩阵，进一步判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的相交情况，如果相交则执行S205，否则执行步骤S203，其中利用位姿判断是否相交的过程如下：

首先建立辅助坐标系，辅助平面与辅助图像，以第一张航拍图像的拍摄位置为坐标原点，以正北方向为Y轴正方向，正西方向为X轴正方向，正地方向为Z轴正方向建立辅助坐标系，辅助平面即为大地平面，辅助图像为相机在辅助坐标系原点，相机坐标系三轴与辅助坐标系三轴重合时，大地平面在相机中所形成的图像；

若第一张航拍图像拍摄位置的坐标为

x₀、y₀、z₀为在辅助坐标系X轴、Y轴、Z轴上的值；旋转为

θ₀、ψ₀为第一张图像拍摄姿态绕辅助坐标系的X轴、Y轴、Z轴的旋转角度；第二张航拍图像拍摄位置的坐标为

三轴的旋转为

x₁、y₁、z₁为在辅助坐标系X轴、Y轴、Z轴上的值；

θ₁、ψ₁为绕辅助坐标系的X轴、Y轴、Z轴的旋转角度，第一张航拍图像拍摄位置与相对于参考坐标系原点的位移矩阵t为:

第一张航拍图像拍摄位姿与相对于辅助坐标系的旋转矩阵R为：

其中：

由相机的成像原理可知，对于空间中的一点，存在如下关系：

其中K为相机的内参矩阵，可由相机标定获得，t，R分别为公式(d)与公式(e)位移矩阵与旋转矩阵，

为空间中的一点，X′、Y′、Z′为在辅助坐标系X轴、Y轴、Z轴上的值，

为空间中该点在相机成像平面中的像素位置，u为像素的横坐标，v为像素的纵坐标；

对于第一张航拍图像到辅助图像，变换关系可由下式表示：

I为单位矩阵，K^-1为K的逆矩阵，d为参考坐标系原点距离参考平面的距离，n^T为参考平面的法向量，

为辅助图像中像素点的位置，u_f、v_f为辅助图像中像素点的横坐标、纵坐标，

为

投影到第一张航拍图像中的像素点的位置，u₀、v₀为第一张航拍图像中像素点的横坐标、纵坐标；其中，为了简化描述，可得下式：

其中H为辅助图像到第一张航拍图像的变换矩阵，而

H₁＝H^-1 公式(l)

为第一张航拍图像到辅助图像之间的变换矩阵,H^-1为H的逆；

同理可得到第二张航拍图像到辅助图像之间的变换矩阵H₂，接着获取两张航拍图像之间的单应矩阵H_t；然后，以其中一张航拍图像为参考图像，将另一张航拍图像投影到参考图像之上，并查看投影后的四个顶点是否在参考图像范围内，如果四个顶点中有任意一个顶点在参考图像范围内，则认为两张航拍图像相交。

S205、使用特征点匹配判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的匹配情况，判定所属的子区域，如果成功匹配则将当前航拍图像放入当前划分好子区域的航拍图像所属子区域，并执行S202，否则执行S203，其中对于特征点匹配步骤如下：

S20501、使用最佳优先搜索(Best-Bin-First,BBF)算法进行两张航拍图像之间的特征点粗匹配，得到匹配点对的数量n_f；

S20502、使用随机抽样一致性(Random Sample Consensus,RANSAC)算法进行特征点的精匹配，得到匹配点对的数量n,如果满足下式：

n>8+0.3*n_f 公式(m)

则认为其匹配成功。

S3、对子区域进行高精度的内部拼接，利用特征点和基于光度对准的网格变换进行精准对齐，并结合步骤S1预处理中获取的网箱掩膜得到绕开网箱的缝合线，得到子区域的区域图，其流程如图3所示，具体如下：

S301、利用特征点匹配来获取每两张航拍图像之间的匹配关系；

S302、利用最大生成树计算出参考图像以及每张航拍图像的拼接顺序，进一步计算出每张航拍图像到参考图像的单应矩阵，并将航拍图像投影到参考图像的平面上；

S303、利用基于光度对准的网格变换依次处理每一张已投影到参考图像的航循环遍历每张航拍图像，每一次循环都固定所有与遍历的当前航拍图像相交的航拍图像，并对当前航拍图像做基于光度对准的网格优化，优化过程如下：

循环遍历每张航拍图像，每一次循环都固定所有与当前航拍图像相交的航拍图像，并对当前航拍图像做基于光度对准的网格优化，优化过程如下：

S30301、网格化图像并对图像进行抽样，具体如下：

将当前遍历到的航拍图像所在区域划分为多个网格，网格的顶点为vertices,并对当前航拍图像按距离进行像素点抽样，抽样后的点表示为sample_point,每一个抽样点都可以利用周围的4个网格顶点的双线性插值来表示，具体如下：

N_s为该抽样点周围的网格顶点集合，w_j为第j个顶点的权重，vertices_j为第j个顶点的位置；

将当前遍历到的航拍图像进行直线检测，对每条直线进行像素点抽样操作，抽样点为line_point,并对抽样点做双线性插值：

其中line_point为直线中的抽样点,N′_l为每个抽样点周围的4个网格顶点集合，w_k为第k个网格顶点的权重，vertices_k为第k个网格顶点的位置；

S30302、求取优化能量，获取优化后网格顶点的位置，对于当前优化的航拍图像与任意一张与其相交的航拍图像，当前优化的航拍图像为目标图像，相交的航拍图像为参考图像，能量定义如下：

第一项是光度对准约束项，其能量E_pm如下所示：

E_pm＝∑_i‖I_t(p_i+τ(p_i))-I_S(p_i)‖² 公式(p)

其中I_t为目标图像的灰度图，I_S为参考图像的灰度图，p_i为重叠区域图像中第i个采样点，τ(p_i)为第i个采样点的偏移量；

第二项是梯度约束项，其能量E_pg如下所示：

E_pg＝∑_i′‖G_t(p′_i′+τ(p′_i′))-G_S(p′_i′)‖² 公式(q)

其中G_t为目标图像的灰度图的梯度图，G_S为参考图像的灰度图的梯度图，p′_i′为重叠区域图像的第i′个采样点，τ(p′_i′)为第i′个采样点的偏移量；

其中，在公式(p)与公式(q)中，p_i+τ(p_i)、p_i、p′_i′、p′_i′+v(p′_i′)都通过公式(n)的插值方式表示；

第三项是三角形约束项，其能量E_s如公式(r)所示：

其中

为变换前的第i_tr个三角形的顶点，

为优化后的第i_tr个三角形的顶点，N_tr为三角形的数量，其中三角形由划分的网格顶点获得；

与

是通过第i_tr个变换前的三角形顶点根据公式(s)计算得出的变换参数，

第四项是直线约束项，其能量E_c为：

其中

第li条线段的起点，

为第li条线段的终点，

为第li条线段上第lj个抽样点，

均由公式(o)的方式表示，且u_l为当前点在线段上的相对比例位置，N_l为直线的数量，

为第li条线段抽样点的数量；

最终的能量E为：

E＝αE_pm+βE_pg+γE_s+ηE_c 公式(u)

其中α、β、γ、η为能量系数，对于多图联合优化，针对当前航拍图像，与所有的相交航拍图像求取公式(u)所述的能量项，并作网格优化；最终最小化能量E_fe，具体为：

其中np为与当前航拍图像相交的航拍图像数量，

为公式(u)表述的与第j_e张匹配图优化的能量项，最小化公式(v)即可得到优化后的网格顶点位置；

每一次优化后重新映射航拍图像并判断当前航拍图像的网格顶点的位置变化量，如果变化量小于事先指定的优化阈值，将该航拍图像移出优化列表，当优化列表为空时，停止优化。

S304、利用对齐质量与网箱掩膜来寻找最佳缝合线，并缝合子区域图像，具体如下所示：

利用对齐质量与网箱掩膜来寻找最佳缝合线，对网箱区域依照网箱掩膜施加能量值，使得缝合线绕开网箱，对于相邻的两张优化后的航拍图像，缝合线的能量E_seam可由下式表示：

其中E_seam为缝合线的总能量，p′_s为航拍图像中的像素；N_o为两张航拍图像重合区域相邻像素点的集合，p_s、q_s为N_o中的任意两个相邻像素点；N_t为两张航拍图像重合区域中的目标掩膜区域相邻像素点的集合,p_t、q_t为N_t中的任意一对相邻像素点，α′、β′为能量系数，能量函数各部分组成部分定义如下：

为非重合区域的能量，其中λ_m为一个惩罚值，

为像素的标签，确定像素是来自第一张航拍图像还是第二张航拍图像，这一项对重叠区域不做处理，但是对非重叠区域施加一个能量惩罚，

E_a为相邻的两张优化后的航拍图像对齐质量能量，

为p_s、q_s两个像素点的标签，即第一张航拍图像或第二张航拍图像；

为p_s点在

图中的像素值，

为p_s点在

图中的像素值,

为q_s点在

图中的像素值,

为q_s点在

图中的像素值，E_t为网箱区域根据网箱掩膜赋予的能量，mask_OE为对于网箱区域惩罚的能量值；结合这三项能量利用马尔科夫随机场获取最佳的缝合线；

找到最佳缝合线后，缝合航拍图像，得到子区域的航拍图像拼接完成图，也即该子区域的子区域图。

S4、将所有子区域的区域图利用位姿获取的单应矩阵投影到全景平面实现全景的低精度拼接，获取自然的全景图，其流程如图4所示，具体如下：

S401、计算全景平面的坐标系以及每个子区域的区域图到全景平面的单应矩阵，全景平面对应坐标系的原点的经度，纬度，相对拍摄平面的高度可以由下式获得：

其中lan为所有航拍图像位置的纬度集合，lon为所有航拍图像的经度集合，hei为所有航拍图像相对拍摄平面的高度集合；

然后基于该中心位置坐标为原点建立空间直角坐标系；其中正西方向为X轴正方向，正北为Y轴正方向，正下为Z轴正方向；坐标系建立完毕之后，将大地平面作为最终的投影平面；随后便可以根据(k)式来计算每个子区域的区域图到全景平面的单应矩阵；

S402、修正子区域的区域图的坐标范围，对子区域的区域图进行坐标修正；因为在进行子区域内部拼接时，子区域里的参考图像会进行横向与纵向上的位移Δx,Δy，以使得所有的航拍图像都能被记录到参考平面上；因此，子区域的区域图在横向与纵向上的坐标范围不再是[0:width,0:height]，修正后的局部区域在横向与纵向上的坐标范围为：[-Δx:width-Δx,-Δy:height-Δy]；其中width，height为每个子区域的区域图的高与宽，Δx与Δy为子区域参考图像在横向与纵向上的位移；

S403、将子区域的区域图投影到全景平面获取全景图，在获取每个子区域的区域图到全景平面的单应矩阵后，便可以获取最终的全景图的大小，以及每个子区域的区域图投影后的位移[offset_x,offset_y]；

其中位移的获取方式具体为，由步骤S402可知每个子区域的区域图的范围为[-Δx:width-Δx,-Δy:height-Δy]，则其顶点依次左上，左下，右上，右下为：

然后利用下式：

求出变换后的顶点，其中

为变换前图像的顶点,x_v、y_v为顶点的横坐标与纵坐标，

为变换后坐标的顶点，u_v、v_v为变换后顶点的纵坐标与横坐标，H_v为子区域的区域图到全景平面的单应矩阵。因为所有的图像都要在最终的全景图中，因此要对所有投影后的图像做位移操作。具体如下：

offset_x＝max(abs(min(0,min(vertice_na[:,0])))…) 公式(e1)

offset_y＝max(abs(min(0,min(vertice_na[:,1])))…) 公式(f1)

其中vertice_na[:,0]为第na个子区域的区域图投影过后的顶点横坐标，vertice_na[:,1]为第na个子区域的区域图投影过后的顶点纵坐标，max为求取最大值的操作，abs为求取绝对值的操作，min为求取最小值的操作；

对于全景图拼接，也就是将子区域的区域图投影到全景图中；投影后的坐标可以表示为：

其中，H_p为每个子区域的区域图到全局的单应矩阵，

为子区域的区域图的每个像素点的坐标，x_a、y_a分别为每个像素点的横坐标与纵坐标，

为投影到全景平面后该像素点的坐标，u_p，v_p为投影后的横坐标与纵坐标,

为位移；投影完毕后，便可以通过最佳缝合线缝合图像获取最终的全景图。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (8)

1.一种基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述图像拼接方法包括如下步骤：

S1、利用无人机对养殖区域进行低空航拍得到航拍图像，对航拍图像执行预处理操作，调整航拍图像的亮度并获取网箱区域的掩膜；

S2、将经过预处理的航拍图像分成多个可以根据特征点进行拼接的子区域，所述子区域内部的图像重叠区域为浮于水面的养殖设施或者船只，所述子区域之间的重叠区域为水；

S3、对子区域进行高精度的内部拼接，利用特征点和基于光度对准的网格变换进行精准对齐，并结合步骤S1预处理中获取的网箱掩膜得到绕开网箱的缝合线，得到子区域的区域图；

S4、将所有子区域的区域图利用位姿获取的单应矩阵投影到全景平面实现全景的低精度拼接，获取全景图。

2.根据权利要求1所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S1中对航拍图像执行预处理操作包括：对航拍图像进行匀光处理来调整亮度与通过u-net深度学习模型获取网箱掩膜。

3.根据权利要求1所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S2过程如下：

S201、使用SIFT算法来检测每张航拍图像中的特征点，并创建一个集合，并第一张航拍图像放入其中，作为第一个子区域；

S202、从第二张航拍图像开始，依次对每一张未划分子区域的航拍图像执行步骤S203；

S203、针对当前未划分子区域的航拍图像，对每一张划分好子区域的航拍图像，执行步骤S204，如果当前未划分子区域的航拍图像对当前所有划分好子区域的航拍图像都无法完成区域匹配，则新建一个集合，并将当前未划分子区域的航拍图像放入新建的集合中，作为一个新的子区域，然后执行S202；

S204、利用位姿计算当前未划分子区域的航拍图像与当前划分好子区域的航拍图像之间的单应矩阵，进一步判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的相交情况，如果相交则执行S205，否则执行步骤S203；

S205、使用特征点匹配判断当前未划分子区域的航拍图像与当前划分好子区域的航拍图像的匹配情况，判定所属的子区域，如果成功匹配则将当前航拍图像放入当前划分好子区域的航拍图像所属子区域，并执行S202，否则执行S203，其中对于特征点匹配步骤如下：

S20501、使用最佳优先搜索算法进行两张航拍图像之间的特征点粗匹配，得到匹配点对的数量n_f；

S20502、使用随机抽样一致性算法进行特征点的精匹配，得到匹配点对的数量n，如果满足下式：

n>8+0.3*n_f 公式(1)

则认为其匹配成功。

4.根据权利要求3所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S204中使用位姿判断两张航拍图像是否相交的过程如下：

首先建立辅助坐标系，辅助平面与辅助图像，以第一张航拍图像的拍摄位置为坐标原点，以正北方向为Y轴正方向，正西方向为X轴正方向，正地方向为Z轴正方向建立辅助坐标系，辅助平面即为大地平面，辅助图像为相机在辅助坐标系原点，相机坐标系三轴与辅助坐标系三轴重合时，大地平面在相机中所形成的图像；

若第一张航拍图像拍摄位置的坐标为

x₀、y₀、z₀为在辅助坐标系X轴、Y轴、Z轴上的值，旋转为

θ₀、ψ₀为第一张图像拍摄姿态绕辅助坐标系的X轴、Y轴、Z轴的旋转角度；第二张航拍图像拍摄位置的坐标为

三轴的旋转为

x₁、y₁、z₁为在辅助坐标系X轴、Y轴、Z轴上的值，

θ₁、ψ₁为绕辅助坐标系的X轴、Y轴、Z轴的旋转角度，第一张航拍图像拍摄位置与相对于参考坐标系原点的位移矩阵t为:

第一张航拍图像拍摄位姿与相对于辅助坐标系的旋转矩阵R为：

其中：

由相机的成像原理可知，对于空间中的一点，存在如下关系：

其中K为相机的内参矩阵，可由相机标定获得，t、R分别为公式(2)的位移矩阵与公式(3)的旋转矩阵，

为空间中的一点，X′、Y′、Z′为在辅助坐标系X轴、Y轴、Z轴上的值，

为空间中该点在相机成像平面中的像素位置，u为像素的横坐标，v为像素的纵坐标；

对于第一张航拍图像到辅助图像，变换关系可由下式表示：

I为单位矩阵，K^-1为K的逆矩阵，d为参考坐标系原点距离参考平面的距离，n^T为参考平面的法向量，

为辅助图像中像素点的位置，u_f、v_f为辅助图像中像素点的横坐标、纵坐标，

为

投影到第一张航拍图像中的像素点的位置，u₀、v₀为第一张航拍图像中像素点的横坐标、纵坐标；其中，为了简化描述，可得下式：

其中H为辅助图像到第一张航拍图像的变换矩阵，而

H₁＝H^-1 公式(10)

为第一张航拍图像到辅助图像之间的变换矩阵,H^-1为H的逆；

同理可得到第二张航拍图像到辅助图像之间的变换矩阵H₂，接着获取两张航拍图像之间的单应矩阵H_t；然后，以其中一张航拍图像为参考图像，将另一张航拍图像投影到参考图像之上，并查看投影后的四个顶点是否在参考图像范围内，如果四个顶点中有任意一个顶点在参考图像范围内，则认为两张航拍图像相交。

5.根据权利要求1所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S3过程如下：

S301、利用特征点匹配来获取每两张航拍图像之间的匹配关系；

S302、利用最大生成树计算出参考图像以及每张航拍图像的拼接顺序，进一步计算出每张航拍图像到参考图像的单应矩阵，并将航拍图像投影到参考图像平面上；

S303、利用基于光度对准的网格变换处理每一张已投影到参考图像平面的航拍图像；

S304、利用对齐质量与网箱掩膜来寻找最佳缝合线，并缝合子区域图像。

6.根据权利要求5所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S303中循环遍历每张航拍图像，每一次循环都固定所有与当前航拍图像相交的航拍图像，并对当前航拍图像做基于光度对准的网格优化，优化过程如下：

S30301、网格化图像并对图像进行抽样，具体如下：

将当前遍历到的航拍图像所在区域划分为多个网格，网格的顶点为vertices,并对当前航拍图像按距离进行像素点抽样，抽样后的点表示为sample_point,每一个抽样点都可以利用周围的4个网格顶点的双线性插值来表示，具体如下：

N_s为该抽样点周围的网格顶点集合，w_j为第j个顶点的权重，vertices_j为第j个顶点的位置；

将当前遍历到的航拍图像进行直线检测，对每条直线进行像素点抽样操作，抽样点为line_point,并对抽样点做双线性插值：

其中line_point为直线中的抽样点,N_l′为每个抽样点周围的4个网格顶点集合，w_k为第k个网格顶点的权重，vertices_k为第k个网格顶点的位置；

S30302、求取优化能量，获取优化后网格顶点的位置，对于当前优化的航拍图像与任意一张与其相交的航拍图像，当前优化的航拍图像为目标图像，相交的航拍图像为参考图像，能量定义如下：

第一项是光度对准约束项，其能量E_pm如下所示：

E_pm＝∑_i‖I_t(p_i+τ(p_i))-I_S(p_i)‖² 公式(13)

其中I_t为目标图像的灰度图，I_S为参考图像的灰度图，p_i为重叠区域图像中第i个采样点，τ(p_i)为第i个采样点的偏移量；

第二项是梯度约束项，其能量E_pg如下所示：

E_pg＝∑_i′‖G_t(p′_i′+τ(p′_i′))-G_S(p′_i′)‖² 公式(14)

其中G_t为目标图像的灰度图的梯度图，G_S为参考图像的灰度图的梯度图，p_i′为重叠区域图像的第i′个采样点，τ(p_i′)为第i′个采样点的偏移量；

其中，在公式(13)与公式(14)中，p_i+τ(p_i)、p_i、p′_i′、p′_i′+τ(p′_i′)都通过公式(11)的插值方式表示；

第三项是三角形约束项，其能量E_s如公式(15)所示：

其中

为变换前的第i_tr个三角形的顶点，

为优化后的第i_tr个三角形的顶点，N_tr为三角形的数量，其中三角形由划分的网格顶点获得；

与

是通过第i_tr个变换前的三角形顶点根据公式(16)计算得出的变换参数，

第四项是直线约束项，其能量E_c为：

其中

第li条线段的起点，

为第li条线段的终点，

为第li条线段上第lj个抽样点，

均由公式(12)的方式表示，且ul为当前点在线段上的相对比例位置，N_l为直线的数量，

为第li条线段抽样点的数量；

最终的能量E为：

E＝αE_pm+βE_pg+γE_s+ηE_c 公式(18)

其中α、β、γ、η为能量系数，对于多图联合优化，针对当前航拍图像，与所有的相交航拍图像求取公式(18)的能量项，并作网格优化；最终最小化能量E_fe，计算公式如下：

其中np为与当前航拍图像相交的航拍图像数量，

为公式(18)表述的与第j_e张匹配图优化的能量项，最小化公式(19)即可得到优化后的网格顶点位置；

每一次优化后重新映射航拍图像并判断当前航拍图像的网格顶点的位置变化量，如果变化量小于事先指定的优化阈值，将该航拍图像移出优化列表，当优化列表为空时，停止优化。

7.根据权利要求5所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S304过程如下：

利用对齐质量与网箱掩膜来寻找最佳缝合线，对网箱区域依照网箱掩膜施加能量值，使得缝合线绕开网箱，对于相邻的两张优化后的航拍图像，缝合线的能量E_seam可由下式表示：

其中E_seam为缝合线的总能量，p_s′为航拍图像中的像素；N_o为两张航拍图像重合区域相邻像素点的集合，p_a、q_s为N_o中的任意两个相邻像素点；N_t为两张航拍图像重合区域中的目标掩膜区域相邻像素点的集合,p_t、q_t为N_t中的任意一对相邻像素点，α′、β′为能量系数，能量函数各部分组成部分定义如下：

为非重合区域的能量，其中λ_m为一个惩罚值，

为像素的标签，确定像素是来自第一张航拍图像还是第二张航拍图像，这一项对重叠区域不做处理，但是对非重叠区域施加一个能量惩罚，

E_a为相邻的两张优化后的航拍图像对齐质量能量，

为p_s、q_s两个像素点的标签，即第一张航拍图像或第二张航拍图像；

为p_s点在

图中的像素值，

为p_s点在

图中的像素值,

为q_s点在

图中的像素值,

为q_s点在

图中的像素值，E_t为网箱区域根据网箱掩膜赋予的能量，mask_OE为对于网箱区域惩罚的能量值；结合这三项能量利用马尔科夫随机场获取最佳的缝合线；

找到最佳缝合线后，缝合航拍图像，得到子区域的航拍图像拼接完成图，也即该子区域的子区域图。

8.根据权利要求1所述的基于非均匀精度的近海养殖区域图像拼接方法，其特征在于，所述步骤S4过程如下：

S401、计算全景平面的坐标系以及每个子区域的区域图到全景平面的单应矩阵；

S402、修正子区域的区域图的坐标范围；

S403、将子区域的区域图投影到全景平面获取全景图。

Images