一种面向水域的无人机遥感多光谱图像拼接方法及系统

1.本发明涉及一种面向水域的无人机遥感多光谱图像拼接方法及系统，属于水域图像处理技术领域。


背景技术：

2.中国专利(公开号cn 110503604 a)公开了一种基于高精度pos的航空面阵影像实时正射拼接方法，包括飞行前，进行相机检校、相机与pos系统的检校；获取每张原始影像的高精度pos数据，利用相机检校结果去除原始影像的图像畸变；将地表平均高程面作为物方平面；计算每张影像的投影矩阵，再通过地表平均高程，计算物方平面到原始像方平面的唯一确定的正定矩阵；确定统一分辨率，得到纠正影像平面到原始像方平面的变换矩阵；更新总纠正影像的范围，重新计算总纠正影像的宽和高；获得当前纠正影像在总纠正影像中的位置，进行影像拼接，对于重叠区部分直接使用影像dom上已包含的部分；迭代执行直到得到最终的总纠正影像。本发明提高了正射纠正影像拼接的精度，满足影像实时拼接的需求。
3.上述方案，运用地理坐标信息进行图像拼接，但由于需要提供高精度的gps信息，否则遥感图像获取的地理信息存在较大的误差，导致拼接图像之间位置错位，而无人机的稳定性和抗风能力以及gps测量设备性能，使得每张图像对应的经纬度pos信息精度难以保证，进而导致拼接精度将受到严重影响。


技术实现要素：

4.针对现有技术的缺陷，本发明的目的一在于提供一种根据经纬度pos信息进行图像的几何校正预处理，得到校正图像，并对校正图像添加地理坐标投影信息；然后根据地理坐标投影信息，进行图像的粗配准，得到配准图像；再根据配准图像，构建计算机图形学知识模型，计算出重叠区域位置，并抠取重叠区域图像，得到待拼接图像；进而构建基于主成分图像的尺度不变特征变换sift模型，提取待拼接图像的特征点；并以空间距离作为特征点的相似性评估指标，筛选出符合要求的匹配点；最后构建基于拉普拉斯laplacian金字塔的多分辨率融合模型，对待拼接图形进行融合，实现无人机遥感多光谱图像的拼接，方案科学合理，切实可行的面向水域的无人机遥感多光谱图像拼接方法。
5.针对现有技术的缺陷，本发明的目的二在于提供一种通过处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法，能够有效提高图像拼接精度，方案简单，实用，便于推广使用的计算机设备。
6.为实现上述目的之一，本发明的第一种技术方案为：
7.一种面向水域的无人机遥感多光谱图像拼接方法，
8.包括以下步骤：
9.第一步，接收图像数据，并合成多波段图像；
10.同时获取多波段图像对应的经纬度pos信息，并剔除经纬高数据异常值；
11.第二步，根据第一步中对应的经纬度pos信息进行图像的几何校正预处理，得到校
正图像；
12.第三步，基于图像对应的经纬和偏航角信息，计算地理坐标投影信息，并对第二步中的校正图像添加地理坐标投影信息；
13.第四步，按照拼接策略选取无人机图像序列，根据第三步中的地理坐标投影信息，进行图像的粗配准，得到平移变换矩阵和两幅配准图像；
14.第五步，根据第四步中的配准图像，利用计算机图形学知识，计算出重叠区域位置，并抠取重叠区域图像，得到待拼接图像；
15.第六步，利用主成分分析法pca对第五步中的待拼接图像提取第一主成分图像，并对第一主成分图像构建尺度不变特征模型sift，以获取第一主成分图像的特征点；
16.第七步，判断第六步中的特征点数量，若数量少于n个，将第四步中的平移变换矩阵作为后期的单应性变化矩阵；
17.否则采用查询索引kd树方式，建立一对多假设匹配集，并以空间距离作为第六步中特征点的相似性评估指标，筛选出符合要求的匹配点；并使用距离直方图约束全局相似属性来剔除误匹配点对，计算单应性变换矩阵；
18.第八步，对第四步中的两幅配准图像构建基于拉普拉斯laplacian金字塔的多分辨率融合模型，并结合第七步中的单应性变换矩阵，对第四步中的两幅配准图像进行融合，实现无人机遥感多光谱图像的拼接；
19.之后返回第四步，直到完成所有序列图像的拼接后退出。
20.本发明经过不断探索以及试验，根据经纬度pos信息进行图像的几何校正预处理，得到校正图像，并对校正图像添加地理坐标投影信息；然后根据地理坐标投影信息，进行图像的粗配准，得到配准图像；再根据配准图像，构建计算机图形学知识模型，计算出重叠区域位置，并抠取重叠区域图像，得到待拼接图像；进而构建基于主成分图像的尺度不变特征变换sift模型，提取待拼接图像的特征点；并以空间距离作为特征点的相似性评估指标，筛选出符合要求的匹配点；最后构建基于拉普拉斯laplacian金字塔的多分辨率融合模型，对待拼接图形进行融合，实现无人机遥感多光谱图像的拼接。
21.进而，本发明特征点融合pos信息的方案，能够有效融合了特征点方法配准精准度高的特点，以及位姿拼接法配准速度快的特点，拼接方案可作为实时配准的基础，并降低了对pos信息的精度要求；同时特征点结合经纬度pos信息可以降低大规模图像拼接中出现的误差累积，从而使得拼接图配准精度明显提高，最终得到的拼接融合图像不仅含有地理坐标信息，而且具有良好视觉效果，满足一般的无人机水质遥感监测时的图像拼接要求，方案科学合理，切实可行。
22.进一步，本发明适用于渔业养殖水域环境，在特征匹配算法中，针对养殖水域图像检测的特征点幅值偏小，方向不明确，以及在图像中存在许多具有局部相似属性的特征点，提供了更好的剔除方法，使得获得的单应性变化矩阵更准确。
23.更进一步，本发明能够处理长条状，规划航线等图像序列的拼接，且最终的拼接图具有地理坐标。
24.优选的，n＝4。
25.作为优选技术措施：
26.所述第一步中，多波段图像的合成方法如下：
27.接收无人机遥感的多光谱图像，并将同时刻拍摄的同一位置的单波段图像进行合成，得到多波段图像；
28.并根据获取单波段图像对应的经纬度pos信息剔除起飞和返航时的经纬高数据异常值；
29.所述经纬度pos信息包括经度lon、纬度lat、航高h、航向角γ、俯仰角α、滚转角β、地面分辨率∈以及测量系统的地理坐标系和投影坐标系；
30.地面分辨率∈的计算公式如下：
31.∈＝pixelsize*h/f
ꢀꢀ
(1)。
32.作为优选技术措施：
33.所述第二步中，几何校正的方法包括以下内容：
34.步骤21，建立图像坐标系o-xy，相机坐标系s-x
syszs
，机体坐标系p-x
pypzp
，地理坐标系e-xeyeze；
35.步骤22，对步骤21中的图像坐标系选择以像主点为原点，沿飞行方向为t轴的正方向，垂直于飞行方向为x轴的正方向时，根据相机与机体的相对位置关系及无人机的飞行姿态，进行图像坐标系下像点g(x,y)与地理坐标系下对应物点g(xe,ye,ze)的变化关系转换；
36.变化关系转换的计算公式如下：
37.[xe,ye,ze]
et
＝λr
eprpsrso
[x,y]
ot
＝λr
eprps
[x,y,-f]
ot
ꢀꢀꢀ
(2)
[0038]
其中，f为焦距，用于表示图像坐标系与相机坐标系在z轴方向的平移量；
[0039]
步骤23，由于相机的镜头光心与无人机的质心重合，且相机坐标轴系经过平移后与机体坐标系完全重合，即r
ps
为单位矩阵i；
[0040]
根据步骤22中的平移量，构建相机成像模型，其计算公式如下：
[0041]
[xe,ye,ze]
et
＝λr
ep
[x,y,-f]
ot
ꢀꢀꢀ
(3)
[0042]
其中，r
ep
＝r(h)r(γ)r(β)r(α)，λ为比例系数，即λ＝h/f；
[0043]
在正直摄影条件下，遥感图像中像点g(x,y)与校正后的像点g
′
(x
′
,y
′
)的满足关系式如下：
[0044]
(x
′
,y
′
,-f)
ot
＝r
ep
(x,y,-f)
ot
ꢀꢀꢀ
(4)
[0045]
其中，r(α)、r(β)、r(γ)、r(h)分别为基于俯仰角的校正旋转矩阵、基于滚动角的校正旋转矩阵、基于偏航角角的校正旋转矩阵、基于航高的校正矩阵；
[0046]
步骤24，对步骤23中的r(α)、r(β)、r(γ)、r(h)分别建立俯仰角校正矩阵、滚动角校正矩阵、偏航角校正矩阵、高度校正矩阵，得到校正数学模型；
[0047]
俯仰角校正矩阵：
[0048][0049]
滚动角校正矩阵：
[0050][0051]
偏航角校正矩阵：
[0052][0053][0054]
高度校正矩阵：
[0055]
其中，θ表示任一像素点和焦点的连线与中心视轴之间的夹角；
[0056]
步骤25，根据步骤24中的校正数学模型，计算图像四角坐标校正后的新坐标，通过旧坐标与新坐标求出变换矩阵，并采用双线性插值重采样方法计算出新的校正图像，完成图像的几何校正，并更新图像的地理坐标信息。
[0057]
作为优选技术措施：
[0058]
所述第三步中，地理坐标投影信息的计算方法如下：
[0059]
利用仿射矩阵参数对遥感图像坐标与地理坐标进行转换，其包括6个参数，分别为xe,x
pixel
,r
γ
,ye,y
pixel
,r
γ
，描述的是图像行列号和地理坐标之间的关系，
[0060]
其中，xe、ye表示图像左上角像元的地理投影坐标，x
pixel
、y
pixel
分别表示图像像元在经度、维度方向的地面分辨率，r
γ
表示图像旋转角度的正弦值；
[0061]
将图像绕中心o点顺时针旋转γ度，然后根据无人机记录的经纬度经过坐标投影计算得到o点坐标(x
oe
,y
oe
)，则图像左上角像点g坐标的计算公式如下：
[0062][0063]
其中w和h分别表示图像尺寸的宽和高；
[0064]
地面分辨率x
pixel
＝-y
pixel
＝-∈，此时r
γ
＝0；
[0065]
图像坐标系下任一点g(row,col)的地理坐标g(xe′
,ye′
)的计算公式如下：
[0066][0067]
作为优选技术措施：
[0068]
所述第四步中，粗配准的过程如下：
[0069]
在拼接的前几轮使用帧到帧拼接策略，后续选择拼接图到拼接图的拼接策略以完成图像的拼接；
[0070]
获取同一航线上相邻的图像i1和图像i2，图像左上角的地理坐标分别为g1(x
e1
,y
e1
)和g2(x
e2
,y
e2
)，以图像i1为基准，求得图像i2相对于图像i1的偏移量；
[0071]
所述偏移量的计算公式如下：
[0072][0073]
根据偏移量，将几何校正后的图像i1和图像i2分别绕o1和o2旋转γ度至航线方向，图像i2中的任一点在图像i1中的图像坐标的计算公式如下：
[0074][0075]
其中，h
rigid
表示平移变化矩阵；
[0076]
配准图像的仿射矩阵参数的左上角像元地理坐标g(x
enew
,y
enew
)的计算公式如下：
[0077][0078]
配准图像的地面分辨率x
pixel
、y
pixel
和r
γ
均保持不变。
[0079]
作为优选技术措施：
[0080]
所述第五步中，抠取重叠区域图像包括以下内容：
[0081]
步骤51，将图像i1和i2分别绕o1和o2旋转γ度至航线方向；
[0082]
步骤52，步骤51中的旋转完成后，根据计算机图形学知识，计算出图像i1四角构成的四边形i
1ai1bi1ci1d
与图像i2四角构成的四边形i
2ai2bi2ci2d
的多边形重叠区域；
[0083]
步骤53，在步骤52中的多边形重叠区域向外添加一定偏置δ，最终获得多边形重叠区域abcd的坐标；
[0084]
步骤54，利用步骤53中的多边形重叠区域abcd的坐标制作掩膜图像，从图像i1和i2中获得只含重叠区域的图像，记为图像i1′
和i2′
。
[0085]
作为优选技术措施：
[0086]
所述第六步中，构建第一主成分图像的尺度不变特征变换sift模型的方法如下：
[0087]
步骤61，输入大小为w
×
h，波段数为c的多光谱图像矩阵i，将i重组为w
×
h行c列矩阵i
reshape
，对i
reshape
每一列进行归一化得到i
′
reshape
；
[0088]i′
reshape
的计算公式如下：
[0089][0090]
步骤62，根据步骤61中的i
′
reshape
，计算协方差矩阵cov；
[0091]
协方差矩阵cov的计算公式如下:
[0092][0093]
步骤63，计算步骤62中的协方差矩阵cov最大特征值对应的特征向量vec；
[0094]
步骤64，对重组图像的i
reshape
提取步骤63中的特征向量vec的主成分图像，并重组为w
×
h大小的图像i
pca
；
[0095]
图像i
pca
的计算公式如下：
[0096]ipca
＝i
reshape
·
vec
ꢀꢀꢀꢀ
(20)
[0097]
步骤65，将步骤64中的主成分图像作为尺度不变特征变换sift特征点提取的输入，然后进行构建dog尺度空间、检测dog尺度空间的极值点，删除不稳定的特征点、对特征点方向进行赋值和生成特征点描述子；
[0098]
所述构建dog尺度空间，包括以下内容：
[0099]
通过对源图像进行高斯模糊和降采样得到图像金字塔，降采样的计算公式如下：
[0100]
o为整数且o≥1
ꢀꢀ
(21)
[0101]
其次对图像金字塔的每层图像使用不同的σ参数进行高斯模糊，得到的多张模糊过图像构成了高斯金字塔，其计算公式如下：
[0102][0103][0104]
然后将相邻的两个高斯空间的图像相减得到dog图像，其计算公式如下：
[0105][0106]
式中，down表示降采样，其中g0＝i；i(x,y)表示源图像，l(x,y,σ)表示原图像卷积后的高斯尺度空间，代表卷积运算，σ表示高斯卷积核的尺度因子；g(x,y,σ)表示高斯核函数，(m,n)代表卷积核的大小；k表示相邻尺度空间的比例因子，取
[0107]
检测dog尺度空间的极值点，包括以下内容：
[0108]
将每个像素点与同一尺度领域的8个像素点和上下相邻的尺度的18个像素点进行比较，只有当该像素点的dog值全部大于或小于对比的26个像素点的dog值时，该像素点是dog尺度空间上一个极值点；
[0109]
删除不稳定的特征点，包括以下内容：
[0110]
首先，为了得到极值点的精确位置，需要对离散空间进行像素差值，然后通过类比拟合三维二次函数以得到极值点精确位置，最后剔除低对比度的极值点；
[0111]
同时，对边缘响应大的极值点进行舍弃；
[0112]
对特征点方向进行赋值，包括以下内容：
[0113]
计算每个特征点的模m(x,y)和方向θ(x,y)的计算公式如下：
[0114][0115][0116]
式中，l(x,y)是特征点所在的尺度空间值；以特征点为中心，计算以3σ为半径的邻域梯度幅值和方向，利用直方图统计领域像素点的梯度方向分布，将梯度方向从0～360
°
分成36等份，把直方图中主峰值的所在方向确定为该特征点的主方向，若存在一个等份方向的峰值大于主峰值80％时，则将该等份方向作为该特征点的辅方向；
[0117]
生成特征点描述子，包括以下内容：
[0118]
将该区域的图像坐标轴旋转到与特征点梯度方向重合；
[0119]
然后取以特征点为中心的16
×
16大小的邻域窗口，计算每一个像素的梯度，距离特征点越近则权值越大，子区域的像素梯度按照σ＝d/2来进行高斯加权，权重计算公式如下；
[0120][0121]
再将该区域细分为4
×
4个子区域，通过对梯度值进行加权的方法统计每个子区域中8个方向的梯度直方图，形成8位的向量描述，这将构成一个4
×4×
8＝128维的描述向量，最后将这些向量归一化后，即为尺度不变特征变换sift特征点的描述符。
[0122]
作为优选技术措施：
[0123]
所述第七步中，匹配点的筛选方法如下：
[0124]
采用查询索引kd树方式，通过空间距离选择与图像i1中每个特征点最接近的图像i2中的n个特征匹配点形成一对多的假设匹配集，并采取空间距离作为特征点相似性的评估指标，空间距离包括欧氏距离d
p
和像素坐标距离dd的加权和，其计算公式如下：
[0125][0126][0127]
其中，e1和e2分别为两个特征点的特征描述子向量；
[0128]
g1和g2为像素点位置坐标；
[0129]
计算每个特征点与假设匹配集点的空间距离ds，其计算公式如下：
[0130]ds
＝α
·dd
(1-α)
·dp
ꢀꢀꢀꢀ
(14)
[0131]
其中，α为匹配点的像素坐标距离项影响因子；
[0132]
匹配筛选策略采用最近邻空间距离与次近邻空间距离之比，其计算公式如下：
[0133]
r＝min_fst/min_scd
ꢀꢀꢀ
(15)
[0134]
其中min_fst为最近空间距离，min_scd次近空间距离；当r《t，t是预先确定的比值，则两个角点互为匹配点对，实现对所有的特征点进行粗匹配；
[0135]
对于局部相似属性的特征点，匹配点的筛选过程如下：
[0136]
计算匹配点对间的距离，将距离值的最大最小值均匀的分为10个区间，每个区间的频率为p＝{p1,
…
,p
10
}，则峰值区间的频率为max(p)，对应的区间为第i个，在区间[i-1,i 1]中的匹配点对是正确匹配点对，该匹配点对集为寻找的精确匹配点对；
[0137]
再根据随机抽样一致ransac算法消除错误特征点对，从而计算出单应性变化矩阵，用单应性变换矩阵乘以图像得到待融合图像。
[0138]
作为优选技术措施：
[0139]
所述第八步中，基于拉普拉斯laplacian金字塔的多分辨率融合模型的构建方法如下：
[0140]
每两幅图像的融合首先建立两幅图像i1和i2的高斯金字塔g1、g2，然后建立对应的4层laplace金字塔图像lap1、lap2；
[0141]
制作一个和图像i1同样大小的mask掩膜图像i
mask
，这个掩模图像代表融合的位置，然后求mask图像的高斯金字塔g
mask
，其代表每个像素点的融合权重；
[0142]
在每个尺度即分辨率下，根据当前尺度的g
mask
将两幅图像的laplace金字塔图像lap1，lap2进行相加，最终得到拼接的laplace金字塔图像lap
fused
；
[0143]
以lap
fused
的最低分辨率图作为起始图，重构得到最高分辨率的拼接结果。
[0144]
为实现上述目的之一，本发明的第二种技术方案为：
[0145]
一种计算机设备，包括：
[0146]
一个或多个处理器；
[0147]
存储装置，用于存储一个或多个程序；
[0148]
当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
[0149]
本发明通过处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法，能够有效提高图像拼接精度，方案简单，实用，便于推广使用。
[0150]
与现有技术相比，本发明具有以下有益效果：
[0151]
本发明经过不断探索以及试验，根据经纬度pos信息进行图像的几何校正预处理，得到校正图像，并对校正图像添加地理坐标投影信息；然后根据地理坐标投影信息，进行图像的粗配准，得到配准图像；再根据配准图像，构建计算机图形学知识模型，计算出重叠区域位置，并抠取重叠区域图像，得到待拼接图像；进而构建基于主成分图像的尺度不变特征变换sift模型，提取待拼接图像的特征点；并以空间距离作为特征点的相似性评估指标，筛选出符合要求的匹配点；最后构建基于拉普拉斯laplacian金字塔的多分辨率融合模型，对待拼接图形进行融合，实现无人机遥感多光谱图像的拼接。
[0152]
进而，本发明特征点融合pos信息的方案，能够有效融合了特征点方法配准精准度高的特点，以及位姿拼接法配准速度快的特点，拼接方案可作为实时配准的基础，并降低了对pos信息的精度要求；同时特征点结合经纬度pos信息可以降低大规模图像拼接中出现的误差累积，从而使得拼接图配准精度明显提高，最终得到的拼接融合图像不仅含有地理坐标信息，而且具有良好视觉效果，满足一般的无人机水质遥感监测时的图像拼接要求，方案科学合理，切实可行。
[0153]
进一步，本发明适用于渔业养殖水域环境，在特征匹配算法中，针对养殖水域图像检测的特征点幅值偏小，方向不明确，以及在图像中存在许多具有局部相似属性的特征点，
提供了更好的剔除方法，使得获得的单应性变化矩阵更准确。
[0154]
更进一步，本发明能够处理长条状，规划航线等图像序列的拼接，且最终的拼接图具有地理坐标。
附图说明
[0155]
图1是本发明实施例方法的流程图；
[0156]
图2是本发明各种坐标系关系示意图；
[0157]
图3是本发明无人机俯仰角变化时沿无人机纵轴剖面的成像几何示意图；
[0158]
图4是本发明无人机俯仰角变化时三维的成像几何示意图；
[0159]
图5是本发明无人机俯仰角和滚动角变化时像平面位置几何关系示意图；
[0160]
图6是本发明无人机偏航角变化时三维的成像几何示意图；
[0161]
图7是本发明图像仿射矩阵参数计算示意图；
[0162]
图8是本发明图像根据地理坐标粗配准计算示意图；
[0163]
图9是本发明dog尺度空间的构建过程示意图；
[0164]
图10是本发明旋转关键点邻域坐标轴示意图；
[0165]
图11是本发明sift特征描述子构造示意图
[0166]
图12是本发明图像特征点匹配连线图；
[0167]
图13是本发明图像帧到帧的拼接策略示意图；
[0168]
图14是本发明图像拼接图到拼接图的拼接策略示意图；
[0169]
图15是本发明建立高斯金字塔和拉普拉斯金字塔laplace示意图。
具体实施方式
[0170]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0171]
相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。
[0172]
本发明面向水域的无人机遥感多光谱图像拼接方法的第一种具体实施例：
[0173]
一种面向水域的无人机遥感多光谱图像拼接方法，
[0174]
包括以下步骤：
[0175]
第一步，接收图像数据，并合成多波段图像；
[0176]
同时获取多波段图像对应的经纬度pos信息，并剔除经纬高数据异常值；
[0177]
第二步，根据第一步中对应的经纬度pos信息进行图像的几何校正预处理，得到校正图像；
[0178]
第三步，基于图像对应的经纬和偏航角信息，计算地理坐标投影信息，并对第二步中的校正图像添加地理坐标投影信息；
[0179]
第四步，根据第三步中的地理坐标投影信息，进行图像的粗配准，得到配准图像；
[0180]
第五步，根据第四步中的配准图像，构建计算机图形学知识模型，计算出重叠区域位置，并抠取重叠区域图像，得到待拼接图像；
[0181]
第六步，构建基于主成分图像的尺度不变特征变换sift模型，提取第五步中待拼接图像的特征点；
[0182]
第七步，采用查询索引kd树方式，建立一对多假设匹配集，并以空间距离作为第六步中特征点的相似性评估指标，筛选出符合要求的匹配点；并使用距离直方图约束全局相似属性来剔除误匹配点对，计算单应性变换矩阵；
[0183]
第八步，根据第七步中的匹配点，按照拼接策略选取无人机图像序列，构建基于拉普拉斯laplacian金字塔的多分辨率融合模型，对待拼接图形进行融合，实现无人机遥感多光谱图像的拼接。
[0184]
如图1、图2所示，本发明面向水域的无人机遥感多光谱图像拼接方法的第二种具体实施例：
[0185]
一种面向水域的无人机遥感多光谱图像拼接方法，
[0186]
包括以下步骤：
[0187]
步骤(1).接收数据合成多波段图像，获取图像对应的pos信息，并剔除经纬高数据异常值；
[0188]
步骤(2).根据对应的pos信息进行图像的几何校正预处理；
[0189]
步骤(3).基于图像对应的经纬和偏航角信息，计算地理坐标投影并给图像添加地理坐标投影信息；
[0190]
步骤(4).根据相互配准的两幅图像的地理坐标信息，进行图像的粗配准；
[0191]
步骤(5).根据粗配准结果，利用计算机图形学知识计算出重叠区域位置，并抠取重叠区域图像；
[0192]
步骤(6).提取待拼接的两幅图像的基于第一主成分图像的sift特征点；
[0193]
步骤(7).采用kd树索引方式，建立一对多假设匹配集，并以空间距离作为特征点相似性的评估指标，筛选出最佳匹配点，进一步使用距离直方图约束全局相似属性来剔除误匹配点对，之后计算单应性变换矩阵；
[0194]
步骤(8).按照拼接策略选取无人机图像序列，每两幅图像之间采用4层laplacian金字塔的多分辨率融合算法。
[0195]
步骤(1)的具体方式为：
[0196]
接收无人机遥感的多光谱图像，并将同时刻拍摄的同一位置的单波段图像进行合成，根据获取图像对应的经纬度pos信息剔除起飞和返航时的经纬高数据异常值，所属经纬度pos信息包括经(lon)维(lay)高(h)、航向角(γ)、俯仰角(α)和滚转角(β)，以及测量系统的地理坐标系和投影坐标系。为了得到地面分辨率(∈)，还需获取相机的焦距(f)和传感器的像元尺寸(pixelsize)，其计算公式如下。为了方便后续的几何校正，还需知道图像的长(w)和宽(h)；
[0197]
∈＝pixelsize*h/f
ꢀꢀꢀꢀꢀ
(1)
[0198]
步骤(2)的具体方式为：
[0199]
在几何校正的过程中,建立图像坐标系o-xy，相机坐标系s-x
syszs
，机体坐标系p-x
pypzp
，地理坐标系e-xeyeze。当图像坐标系选择以像主点为原点，沿飞行方向为y轴的正方
向，垂直于飞行方向为x轴的正方向时，根据相机与机体的相对位置关系及无人机的飞行姿态，像点g(x,y)与对应物点g(xe,ye,ze)的变化关系可通过式(2)转换，而图像坐标系与相机坐标系存在z轴方向的平移，平移量为焦距f。
[0200]
[xe,ye,ze]
et
＝λr
eprpsrso
[x,y]
ot
＝λr
eprps
[x,y,-f]
ot
ꢀꢀꢀ
(2)
[0201]
假设相机的镜头光心与无人机的质心重合，且相机坐标轴系经过平移后与机体坐标系完全重合，即r
ps
为单位矩阵i。因此相机成像模型可简化为：
[0202]
[xe,ye,ze]
et
＝λr
ep
[x,y,-f]
ot
ꢀꢀꢀ
(3)
[0203]
其中r
ep
＝r(h)r(γ)r(β)r(α)，λ为比例系数，即λ＝h/f。
[0204]
因此，在正直摄影条件下，遥感图像中像点g(x,y)与校正后的像点g
′
(x
′
,y
′
)满足：
[0205]
(x
′
,y
′
,-f)
ot
＝r
ep
(x,y,-f)
ot
ꢀꢀ
(4)
[0206]
其中，r(α)、r(β)、r(γ)、r(h)分别为基于俯仰角的校正旋转矩阵、基于滚动角的校正旋转矩阵、基于偏航角角的校正旋转矩阵、基于航高的校正矩阵。在无人机其他飞行姿态参数不变的情况下，仅考虑某一姿态变化时的情形，可对r(α)、r(β)、r(γ)、r(h)分别建立如图3、图4、图5和图6所示的几何关系，然后计算得到如下的校正数学模型。
[0207]
俯仰角校正矩阵：
[0208][0209]
滚动角校正矩阵：
[0210][0211]
偏航角校正矩阵：
[0212]
[0213][0214]
高度校正矩阵：
[0215]
最后，计算图像四角坐标校正后的新坐标，通过旧坐标与新坐标可求出变换矩阵，然后采用双线性插值重采样方法计算出新的校正图像，完成几何校正。
[0216]
步骤(3)的具体方式为：
[0217]
仿射矩阵参数可以用来对遥感图像坐标与地理坐标进行转换，其包括6个参数(xe,x
pixel
,r
γ
,ye,y
pixel
,r
γ
)，描述的是图像行列号和地理坐标之间的关系，其中xe、ye表示图像左上角像元的地理投影坐标，x
pixel
、y
pixel
分别表示图像像元在经度、维度方向的地面分辨率，r
γ
表示图像旋转角度的正弦值。如图7所示，将图像绕o点顺时针旋转γ度，然后根据无人机记录的经纬度经过坐标投影计算可得到o点坐标(x
oe
,y
oe
)，则图像左上角像点g坐标可通过式(13)获得，而地面分辨率x
pixel
＝-y
pixel
＝-∈，此时r
γ
＝0。其中w和h分别表示图像尺寸的宽和高。
[0218][0219]
而图像坐标系下任一点g(row,col)的地理坐标g(xe′
,ye′
)可以通过式(14)计算得到。
[0220][0221]
步骤(4)的具体方式为：
[0222]
如图8所示，假设同一航线上相邻的图像i1和图像i2，图像左上角的地理坐标分别为g1(x
e1
,y
e1
)和g2(x
e2
,y
e2
)，以图像i1为基准，可求得图像i2相对于图像i1的偏移量，公式如下：
[0223][0224]
将几何校正后的图像i1和图像i2分别绕o1和o2旋转γ度至航线方向，其仍是标准的正直摄影，视轴相互平行，且相机坐标系中心到大地平面的距离相等，因此图像i1和图像i2的拼接可认为是刚体平移变化，图像i2中的任一点可通过式(16)计算其在图像i1中的图像坐标。
[0225][0226]
其配准图像的仿射矩阵参数的左上角像元地理坐标g(x
enew
,y
enew
)可根据式(17)求得，配准结果图像的地面分辨率x
pixel
、y
pixel
和r
γ
均保持不变。
[0227][0228]
步骤(5)的具体方式为：
[0229]
将图像i1和i2分别绕o1和o2旋转γ度至航线方向后，然后可根据计算机图形学知识，计算出图像i1四角构成的四边形i
1ai1bi1ci1d
与图像i2四角构成的四边形i
2ai2bi2ci2d
的多边形重叠区域，在重叠区域向外添加一定偏置δ，最终获得多边形重叠区域abcd的坐标，如图8所示。之后，利用多边形重叠区域abcd的坐标制作掩膜图像，就可以从图像i1和i2中获得只含重叠区域的图像，记为图像i1′
和i2′
。
[0230]
步骤(6)的具体方式为：
[0231]
1).提取多光谱图像第一主成分图像的具体步骤如下：
[0232]
a)输入大小为w
×
h，波段数为c的多光谱图像矩阵i，将i重组为w
×
h行c列矩阵i
reshape
，对i
reshape
每一列进行归一化得到i
′
reshape
：
[0233][0234]
b)求i
′
reshape
的协方差矩阵cov:
[0235][0236]
c)求协方差矩阵cov最大特征值对应的特征向量vec；
[0237]
d)对重组图像i
reshape
提取特征向量vec的主成分图像，并重组为w
×
h大小的图像i
pca
。
[0238]ipca
＝i
reshape
·
vec
ꢀꢀꢀꢀ
(20)
[0239]
2).将pca提取的第一主成分图像作为sift特征点提取的输入，然后进行构建dog尺度空间、检测dog尺度空间的极值点，删除不稳定的特征点、对特征点方向进行赋值和生成特征点描述子五个小步骤，具体实施如下：
[0240]
构建dog尺度空间：
[0241]
通过对源图像进行高斯模糊和降采样得到图像金字塔，降采样过程可以用式(21)表示；其次对图像金字塔的每层图像使用不同的σ参数进行高斯模糊，得到的多张模糊过图像构成了高斯金字塔，可用式(22)和式(23)表示；然后将相邻的两个高斯空间的图像相减就得到了dog图像，可用式(24)表示，整个构建过程可以用图9表示。
[0242]
o为整数且o≥1
ꢀꢀ
(21)
[0243][0244][0245][0246]
式中，down表示降采样，其中g0＝i；i(x,y)表示源图像，l(x,y,σ)表示原图像卷积后的高斯尺度空间，代表卷积运算，σ表示高斯卷积核的尺度因子；g(x,y,σ)表示高斯核函数，具体如式(24)所示，(m,n)代表卷积核的大小；k表示相邻尺度空间的比例因子，一般取
[0247]
检测dog尺度空间的极值点：
[0248]
将每个像素点与同一尺度领域的8个像素点和上下相邻的尺度的18个像素点进行比较，只有当该像素点的dog值全部大于或小于对比的26个像素点的dog值时，即认为该点是dog尺度空间上一个极值点。
[0249]
删除不稳定的特征点：
[0250]
首先，为了得到极值点的精确位置，需要对离散空间进行像素差值，然后通过类比拟合三维二次函数以得到极值点精确位置，最后剔除低对比度的极值点；另外，对边缘响应比较大导致不太稳定的极值点也进行舍弃。
[0251]
对特征点方向进行赋值：
[0252]
为了使特征点具有旋转不变性，需要确定特征点的最大梯度方向，因此在高斯图像中，定义每个特征点的模m(x,y)和方向θ(x,y)如下：
[0253][0254][0255]
式中，l(x,y)是特征点所在的尺度空间值。以特征点为中心，计算以3σ为半径的邻域梯度幅值和方向，利用直方图统计领域像素点的梯度方向分布，将梯度方向从0～360
°
分成36等份，把直方图中主峰值的所在方向确定为该特征点的主方向，若存在一个等份方向的峰值大于主峰值80％时，则将该等份方向作为该特征点的辅方向。
[0256]
生成特征点描述子：
[0257]
为了确保描述子具有旋转不变性，需要将该区域的图像坐标轴旋转到与特征点梯度方向重合，如图10所示，红色箭头表示特征点梯度方向。
[0258]
然后取以特征点为中心的16
×
16大小的邻域窗口，计算每一个像素的梯度，距离特征点越近则权值越大，根据lowe的研究，子区域的像素梯度按照σ＝d/2来进行高斯加权，权重计算公式如下。
[0259][0260]
再将该区域细分为4
×
4个子区域，通过对梯度值进行加权的方法统计每个子区域中8个方向的梯度直方图，形成8位的向量描述，这将构成一个4
×4×
8＝128维的描述向量，最后将这些向量归一化后，即为sift特征点的描述符，如图11所示。
[0261]
步骤(7)的具体方式为：
[0262]
采用kd树索引方式，通过欧式距离选择与图像i1中每个特征点最接近的图像i2中的n个特征匹配点形成一对多的假设匹配集，采取空间距离作为特征点相似性的评估指标，假设两个特征点的特征描述向量分别为e1和e2，位置坐标为g1和g2，那么欧氏距离d
p
和像素坐标距离dd的计算公式如下：
[0263][0264][0265]
为了限制匹配点不在全局匹配中搜索匹配点，在计算空间距离时引入了匹配点间
的像素坐标距离项，影响因子为α，然后计算每个特征点与假设匹配集点的空间距离ds，其计算公式如下：
[0266]ds
＝α
·dd
(1-α)
·dp
ꢀꢀꢀꢀ
(14)
[0267]
匹配筛选策略采用最近邻空间距离与次近邻空间距离之比，即：
[0268]
r＝min_fst/min_scd
ꢀꢀꢀ
(15)
[0269]
其中min_fst为最近空间距离，min_scd次近空间距离。当r《t(t是预先确定的比值)，则两个角点互为匹配点对，实现对所有的特征点进行粗匹配。
[0270]
上面的粗匹配仅约束了全局匹配，但是对于局部相似属性的特征点，可以假设在粗匹配集中有许多正确的匹配对，错误匹配是随机分布的，那么匹配点对的距离直方图中将会出现一个峰值，则正确匹配将分布在峰值的周围。具体步骤可表示为：计算匹配点对间的距离，将距离值的最大最小值均匀的分为10个区间，每个区间的频率为p＝{p1,
…
,p
10
}，则峰值区间的频率为max(p)，对应的区间为第i个，在区间[i-1,i 1]中的匹配点对认为是正确匹配点对，该匹配点对集为寻找的精确匹配点对。最后再根据ransac算法进一步消除错误特征点对，从而计算出单应性变化矩阵，用单应性变换矩阵乘以图像得到待融合图像。根据上述方法，生成的特征匹配点对的连线如图12所示。
[0271]
步骤(8)的具体方式为：
[0272]
在拼接的前4轮使用帧到帧拼接策略，后续选择拼接图到拼接图的拼接策略以完成图像的拼接，帧到帧的拼接策略如图13所示，拼接图到拼接图的拼接策略如图14所示。每两幅图像的融合首先建立两幅图像i1和i2的高斯金字塔g1、g2，然后建立对应的4层laplace金字塔图像lap1、lap2。制作一个和图像i1同样大小的mask掩膜图像i
mask
，这个掩模图像代表了融合的位置，然后求mask图像的高斯金字塔g
mask
，其代表每个像素点的融合权重。在每个尺度(分辨率)下，根据当前尺度的g
mask
将两幅图像的laplace金字塔图像lap1，lap2进行相加，最终得到拼接的laplace金字塔图像lap
fused
。以lap
fused
的最低分辨率图作为起始图，重构得到最高分辨率的拼接结果，整个构建构成如图15所示。
[0273]
应用本发明方法的一种装置实施例：
[0274]
一种计算机设备，其包括：
[0275]
一个或多个处理器；
[0276]
存储装置，用于存储一个或多个程序；
[0277]
当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
[0278]
应用本发明方法的一种计算机介质实施例：
[0279]
一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述的一种面向水域的无人机遥感多光谱图像拼接方法。
[0280]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0281]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程
图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0282]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0283]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0284]
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。