[发明专利]一种稳健的多源卫星遥感影像配准方法

摘要

本发明公开了一种基于互信息和分块随机采样一致性的多源遥感影像配准方法，其步骤：A、各层金字塔影像生成和特征点提取。B、在最高层金字塔影像，利用全局配准法，计算主、副影像间的仿射变换系数，估计影像间旋转角度和分辨率差异系数。C、准确预测同名点初始位置，对匹配窗口影像进行几何粗纠正，利用归一化互信息寻找同名点。D、利用二次多项式和分块RANSAC算法剔除错误的匹配点。E、重复步骤C和D直到原始影像层，基于线性橡皮拉伸法实现影像精确配准。该方法极大减少了同名点量测的人工编辑量，提高多源卫星遥感影像的数据处理自动化程度，具有显著的经济和社会效益。

主权项

1.一种稳健的多源卫星遥感影像配准方法，其步骤是：A、多源卫星遥感影像匹配的数据预处理(1)：多源卫星遥感影像匹配的数据预处理包括：金字塔影像的生成、特征点的提取：步骤是：a、金字塔影像的生成：采用3×3像元平均法，将测区所有影像生成3级金字塔影像,首先对原始影像的每3×3个像元计算其平均灰度值，并赋给第1级金字塔影像的对应像元，生成第一级金字塔影像，直至生成第2级金字塔影像；b、特征点的提取：将影像划分为均匀的格网，格网点数依据影像大小给定，将影像划分为30×30以内的格网，采用摄影测量中常用的Forstner特征提取算子在每个格网内提取一个特征点,格网内的特征信息不明显，将格网中心点作为特征点；B、特征点初始点位的预测(2)：为了减少同名点的搜索时间，提高影像匹配的成功率，对特征点初始点位进行准确预测：步骤如下：a、最高层金字塔影像匹配的同名点预测：在最高层金字塔影像上，首先，通过全局匹配算法，利用化算法求解待配准影像间的几何变换参数,在最高层金字塔影像上，视差变化比较平缓，利用仿射变换来近似表达影像间的几何关系,具体如下：x′＝p₁+p₂x+p₃y    (1)y′＝p₄+p₅x+p₆y式中（x′,y′)、（x,y)分别表示同名点在主影像和副影像上的像点坐标，p₁,p₂,p₃,p₄,p₅,p₆为求解的仿射变换参数；全局匹配算法是一个迭代寻优的过程，目标方程F取值最大，所得的参数即为所求，即p^=argmaxpF(R,S;p)---(2)]]>式中，F为目标方程，R为主影像，S为副影像，p为待求解的仿射变换系数向量(p₁,p₂,p₃,p₄,p₅,p₆)^T,表示使目标方程F最大的参数p的寻优过程，表示使主、副影像归一化互信息最大的仿射变换系数向量；归一化互信息NMI的计算公式如下：NMI(A,B)=H(A)+H(B)H(A,B)---(3)]]>式中，A，B分别代表待匹配的影像块，H(A)和H(B)分别是图像块A和B的平均互信息量，H(A,B)是两个图像块的相关平均信息量;A和B的平均互信息和相关平均信息量如下计算：H(A)=Σa-PA(a)logPA(a)]]>H(B)=Σb-PB(b)logPB(b)---(4)]]>H(A,B)=Σa,b-PA,B(a,b)logPA,B(a,b)]]>式中，P_A(a)和P_B(b)为边缘概率密度函数，P_A,B(a,b)为联合概率密度函数,由下式得到：PA(a)=ΣaPA,B(a,b)]]>PB(b)=ΣbPA,B(a,b)---(5)]]>PA,B(a,b)=h(a,b)Σa,bh(a,b)]]>式中，h表示影像块A和B的联合直方图,是一个二维矩阵：这里，h(a,b)是联合直方图h的元素，表示图像块A中具有灰度值a，在图像块B中有灰度值b的点对个数，M和N分别表示图像块A、图像块B的灰度级数；其他层金字塔影像匹配的同名点预测：在其他金字塔影像层，利用上层匹配结果，依据邻近点视差连续的特性，对同名点初始位置进行预测；步骤如下：1）计算影像间的旋转角度和分辨率差异系数：利用最高层金字塔影像匹配中全局匹配算法所解算出的6个仿射变换系数，对影像间的旋转角度和分辨率差异系数进行近似估计；对于影像间的旋转角度，首先计算旋转角度的两个估值α₁,α₂：α1=arctan(a2/a1)α2=arctan(-b1/b2)---(7)]]>式中，arctan代表三角函数中的反正切函数运算,将α₁,α₂的平均值作为两影像旋转角度的估值同理，首先计算影像间的分辨率差异系数的两个估值λ₁、λ₂：λ1=a12+a22λ2=a32+a42---(8)]]>将λ₁、λ₂的平均值作为两幅影像分辨率差异系数的估值2）对于上层匹配成功的点，直接投影得到在该层的初始位置；对于匹配失败的点，利用邻近点视差预测点位，考虑影像间的旋转角度和分辨率差异带来的影响；设p₁(x₁,y₁)为主影像上层匹配失败点，p₂(x₂,y₂)为其邻近的匹配成功点，p′₂(x′₂,y′₂)为p₂在副影像上的同名点坐标；首先计算p₁和p₂两点在x,y方向上的位移dx=x2-x1dy=y2-y1---(9)]]>然后计算旋转α角后的位移大小dx′=dxcosα-dysinαdy′=dxsinα+dycosα---(10)]]>再计算补偿分辨率差异后的位移大小dx′′=λdx′dy′′=λdy′---(11)]]>p₁点的同名点在副影像上的初始位置(x′₁,y′₁)为：x1′=x2′+dx′′y1′=y2′+dy′′---(12)]]>C、带几何粗纠正的归一化互信息匹配(3)：减少影像旋转角度和几何变形的不利影响，利用归一化互信息测度搜索同名点之前，对匹配窗口影像进行几何粗纠正：步骤如下：a、匹配窗口影像的几何粗纠正：直接对待匹配特征点的搜索窗口影像进行几何粗纠正，这样遍历搜索窗口中的像素进行匹配时，无需每次进行重采样来获得匹配窗口影像；待匹配特征点在副影像上的初始位置A为原点，建立局部坐标系A-x′y′和A-xy，其中x轴、y轴分别与副影像像素坐标系对应坐标轴相互平行，x′轴、y′轴分别与旋转补偿后的搜索窗口对应坐标轴相互平行，两个坐标系间的夹角大小近似为影像间的旋转角度α，利用公式(13)计算搜索窗口中的像素在局部坐标系A-xy中的坐标：x=(x′cosα-y′sinα)×λy=(x′sinα+y′cosα)×λ---(13)]]>加上初始点位A对应的像素坐标，得到搜索窗口中该像素在副影像上的像素坐标,通过双线性内插法，计算出对应的灰度值,遍历处理搜索窗口中每个像素，最终生成经过旋转角度和分辨率差异补偿后的搜索窗口影像；b、归一化互信息匹配：按上步方法进行匹配窗口影像的几何粗纠正，利用公式(3)计算归一化互信息，匹配窗口大小为61×61像素,最后，选取搜索窗口中归一化互信息值最大的像素点作为该特征点的同名点；D、基于分块RANSAC剔除错误匹配点(4)：遥感卫星运行在轨道上，地形起伏引起的像点位移相对小，待匹配影像的像点坐标之间存在系统性，通常用二次多项式模型描述,提出的分块RANSAC剔除错误匹配点法，将整幅遥感影像分成3×3个分块或4×4个分块或5×5个分块，然后对每块中的匹配点分别用二次多项式模型描述其几何约束关系，并辅之于RANSAC算法来自动探测和剔除错误匹配点；二次多项式模型具体形式为：xs=a0+a1xm+a2ym+a3xmym+a4xm2+a5ym2ys=b0+b1xm+b2ym+b3xmym+b4xm2+b5ym2---(14)]]>式中，(x_m,y_m),(x_s,y_s)分别为同名像点在主影像和副影像上的像素坐标，RANSAC算法的求解过程如下：1)设样本数k为无穷大，样本计数器t为0；2)从当前分块中匹配特征点中随机抽取6对点建立方程组，解算二次多项式的12个系数；3)计算当前分块中其余特征点经多项式变换后与其匹配点间的距离,距离小于给定阈值，该对同名点正确匹配，最后统计错误匹配点对所占的比例4)计算样本数值k=log(1-η)log(1-(1-ϵ6)t),]]>η取0.99；5)样本计数器t加1；6)k＜t时，停止循环；7)迭代终止后，选取有最多正确匹配点对的集合重新计算二次多项式系数，并统计各对同名点的拟合误差；E、基于线性橡皮拉伸法进行影像精确配准(5)：首先将所有主、副影像上的同名点对组成三角网，对每个三角形区域，通过最小二乘算法计算仿射变换系数，每个三角形区域利用，利用对应的仿射变换系数进行坐标转换，通过双线性内插计算灰度值，最终完成影像间的精确配准。