一种大倾角航空影像与正射影像配准的方法

本发明属于倾角影像与正射影像配准。特别涉及一种大倾角航空影像与正射影像配准的方法。

背景技术：

1、相对于正射影像，倾角影像可以从多个角度观察地物，更加真实的反映地物的实际情况，极大地弥补了基于正射影像应用的不足。将倾角影像和正射影像进行融合应用，可以对地物信息进行快速关联、整合与分析，对遥感解译技术发展有重要意义。图像配准作为遥感图像数据融合的第一步，为遥感图像地球观测应用奠定了基础。图像配准是指将由不同传感器、不同视角和不同时间获取的同一场景的两个或多个图像进行几何对齐和配准的过程，图像配准的精度对后续应用起着至关重要的作用，如何实现高效率高精度的倾角航空影像配准是本领域技术人员所关注的问题。

2、图像配准一般可以分为基于区域的方法和基于特征的方法。基于区域的方法是依赖地理参考技术对图像进行粗略配准，消除图像对之间明显的平移和旋转差异。但是这种方法并不适应于大倾角航空影像，因为大倾角航空影像存在非常严重的地物几何变形，很难获得其地理坐标信息。基于特征的方法是通过评估图像结构特征信息实现图像配准，这种方法更具可控性和稳健性。目前，常见的配准方法一般适用于倾斜角度较小的影像配准，而对大倾角航空影像配准并不适用，而针对大倾角航空影像配准方面的专利技术和论文比较匮乏。

3、文献1(lowe d g.distinctive image features from scale-invariantkeypoints[j].international journal ofcomputervision,2004,60:91-110.)为尺度不变特征转换(sift)，该方法可以很好的应对尺度、旋转情况下的影像匹配，但是对于倾斜影像的特征点匹配只能提取很少量的特征，配准误差大。文献2(zhao x,zhu q,xiaox,et al.automatic matching method for aviation oblique images based onhomography transformation[j].journal ofcomputer applications,2015,35(6):1720.)提出一种大倾角航空倾斜影像自动配准方法h-sift，采用影像粗略外方位元素计算两幅影像a和b之间的单应变换矩阵，利用单应变换矩阵得到原始影像a的校正影像a'，然后对校正影像a'和影像b进行sift匹配，再将校正影像a'上匹配出的特征点反算到原始影像a上。这种方法严重依赖于初始匹配点的分布和准确率，匹配结果不稳定。文献3(morel jm,yu g.asift:anew framework for fully affine invariant image comparison[j].siamjournal on imaging sciences,2009,2(2):438-469.)提出了一种仿射不变的图像特征匹配算法(asift)，解决了sift应对倾斜影像的匹配问题。asift具有完全的仿射不变性且在大倾角影像配准中取得了不错的效果。但这种方法在大倾角航空影像配准精度低。文献4(赵朝贺,杨化超,张磊,等.一种大角度亚像素的图像自动配准算法[j].测绘通报,2014(8):30-35.)利用asift的完全放射不变性，提出一种基于asift和harris集成互补不变特征的大角度影像亚像素自动配准算法。但是，大倾角航空影像存在严重的地物几何变形、分辨率差异、影像旋转、地物相互遮挡等问题。基于sift和asift的方法采用高斯核构建尺度空间且依赖图像灰度信息和梯度特征的构建特征描述算子，不能很好地应对航空影像地物信息复杂的问题，提取丰富的结构特征。另外，灰度信息和梯度特征对图像噪声较敏感，构建的描述算子不具鲁棒性，导致配准成功率低且误差较大。

技术实现思路

1、本发明旨在克服现有技术缺陷，目的在于提供一种配准误差小、抗噪性能好、成功率高和结果稳定的大倾角航空图像与正射影像配准的方法。

2、为了实现上述目的，本发明采用的技术方案是：

3、步骤1、用asift算法构建大倾角航空影像s的模拟仿射影像s'

4、步骤1.1、参数采样

5、根据asift算法，在一个半球形空间中，通过相机在经度方向角度φ和纬度方向角度θ的方向参数模拟出大倾角航空影像s的视角；对经度方向角度φ和纬度方向角度θ分别进行离散化。其中：

6、纬度方向角度θ采用几何级数离散化方式：

7、tn＝an                                             (1)式(1)中：

8、tn表示相机在纬度方向的第n个采样级数下的倾斜度；

9、n表示相机在纬度方向的采样级数，n＝1,...,n；

10、a表示相机在纬度方向的采样间隔，a＞1。

11、经度方向角度φ采用算术级离散化方式：

12、

13、式(2)中：

14、表示相机在纬度方向的倾斜度tn下，在经度方向的第k个采样级数下的旋转角度；

15、b/tn表示相机在纬度方向的倾斜度tn下，在经度方向的采样间隔，b为一个角度常数；

16、k表示相机在在纬度方向的倾斜度tn下，在经度方向的采样级数，k＝1,...,k，

17、

18、步骤1.2、构建大倾角航空影像s的模拟仿射影像s'

19、设相机在纬度方向上的倾斜度为tn，在经度方向的旋转角度为大倾角航空影像s的模拟仿射影像的构建过程是：先对大倾角航空影像s进行双线插值旋转操作，得到旋转图像；再用反锯齿滤波器与旋转图像进行卷积，得到滤波之后的旋转图像；然后对滤波之后的旋转图像进行倾斜操作，得到相机在纬度方向上的倾斜度为tn和在经度方向的旋转角度为时大倾角航空影像s的模拟仿射影像

20、

21、式(3)中：

22、表示相机在纬度方向上的倾斜度为tn在经度方向的旋转角度为时的大倾角航空影像s的模拟仿射影像；

23、tilt表示对经反锯齿滤波器滤波后的旋转图像进行倾斜操作；

24、gδ表示反锯齿滤波器，其中δ表示反锯齿滤波器的核标准差，δ＝0.8tn；

25、*表示卷积操作；

26、rot表示对大倾角航空影像s进行双线插值旋转操作；

27、s表示大倾角航空影像；

28、表示相机在纬度方向的倾斜度tn下，在经度方向的第k个采样级数下的旋转角度；

29、tn表示相机在纬度方向的第n个采样级数下的倾斜度。

30、以此类推，得到大倾角航空影像s的其余的模拟仿射影像，则大倾角航空影像s的模拟仿射影像

31、为叙述简便起见：以下将“大倾角航空影像s的模拟仿射影像s'”简称为“模拟仿射影像s'”；将“相机在纬度方向上的倾斜度为tn和在经度方向的旋转角度为时大倾角航空影像s的模拟仿射影像”简称为“子模拟仿射影像”。

32、步骤2、对模拟仿射影像s'和正射影像r进行配准

33、步骤2.1、非线性尺度空间构建

34、将子模拟仿射影像的非线性尺度空间离散化为1个基层空间和l个子级空间，离散化后子模拟仿射影像的非线性尺度空间所对应的尺度为：

35、

36、式(4)中：

37、i表示子模拟仿射影像的非线性尺度空间的第i层，i＝0,1,2,....,l；

38、σi表示子模拟仿射影像的第i层非线性尺度空间的尺度；

39、σ0表示子模拟仿射影像非线性尺度空间的初始尺度，σ0为一个常数；

40、l表示子模拟仿射影像的子级空间的个数。

41、设子模拟仿射影像的非线性尺度空间为：其中：表示子模拟仿射影像的基层空间，依次表示子模拟仿射影像的第1个子级空间，……，第l个子级空间。

42、子模拟仿射影像的基层空间构建过程如下：

43、用高斯滤波器对子模拟仿射影像进行平滑滤波，高斯滤波器的高斯核标准差与子模拟仿射影像非线性尺度空间的初始尺度σ0相同，得到子模拟仿射影像的基层空间

44、子模拟仿射影像的第j个子级空间构建过程如下：

45、先以采样率1/σ0对子模拟仿射影像的第j-1层非线性尺度空间进行下采样，得到采样空间再用各向异性扩散滤波对采样空间进行滤波，得到子模拟仿射影像的第j个子级空间

46、

47、式(5)中：

48、表示在时间步长条件下得到的各项异性扩散滤波空间；

49、j表示子模拟仿射影像的子级空间的第j个，j＝1,2,...,l；

50、↓表示下采样操作；

51、表示对子模拟仿射影像的第j-1层非线性尺度空间进行下采样得到的采样空间；

52、t表示时间度量；

53、div表示散度操作符；

54、d(x,y,t)表示扩散系数，扩散系数取决于梯度范数

55、表示梯度操作符；

56、表示子模拟仿射影像的第j个子级空间；

57、id表示单位矩阵；

58、σj表示子模拟仿射影像的第j层非线性尺度空间的尺度；

59、σj-1表示子模拟仿射影像的第j-1层非线性尺度空间的尺度；

60、l表示各向异性扩散的方向轴；

61、m表示各向异性扩散的方向轴个数；

62、al表示采样空间沿着各向异性扩散的方向轴l的扩散系数矩阵，al为扩散系数d的离散形式。

63、根据式(5)，得到子模拟仿射影像的l个子级空间则子模拟仿射影像的非线性尺度空间为：以此类推，得到模拟仿射影像s'中的其余子模拟仿射影像的非线性尺度空间。

64、参照得到模拟仿射影像s'中的其余子模拟仿射影像的非线性尺度空间的方法，得到正射影像r的非线性尺度空间为：r0,r1,...,rl。

65、步骤2.2、特征点检测

66、为叙述方便，将相位一致性用“pc”表示。

67、对子模拟仿射影像进行特征点检测，所述特征点检测过程如下：

68、先用pc提取子模拟仿射影像的第i层非线性尺度空间影像的结构特征，得到子模拟仿射影像的第i层非线性尺度空间的pc特征图；再用x方向和y方向的sobel算子模板对子模拟仿射影像的第i层非线性尺度空间的pc特征图进行卷积，得到子模拟仿射影像的第i层非线性尺度空间的pc特征图在x和y方向上的一阶导数：

69、

70、

71、式(6)～(7)中：

72、i表示子模拟仿射影像的非线性尺度空间的第i层，i＝0,1,2,...,l；

73、表示子模拟仿射影像的第i层非线性尺度空间的pc特征图在x方向上的一阶导数；

74、表示子模拟仿射影像的第i层非线性尺度空间的pc特征图在y方向上的一阶导数；

75、pci表示子模拟仿射影像的第i层非线性尺度空间的pc特征图；

76、*表示卷积操作；

77、γx表示x方向的sobel算子模板；

78、γy表示y方向的sobel算子模板。

79、然后用得到的子模拟仿射影像的第i层非线性尺度空间的pc特征图在x和y方向上的一阶导数，构建子模拟仿射影像的第i层非线性尺度空间的特征点度量矩阵mi：

80、

81、式(8)中：

82、i表示子模拟仿射影像的非线性尺度空间的第i层，i＝0,1,2,....,l；

83、mi表示子模拟仿射影像的第i层非线性尺度空间的特征点度量矩阵；

84、表示高斯核，高斯核的标准差为其中σi表示子模拟仿射影像的第i层尺度空间的尺度；

85、*表示卷积操作；

86、表示子模拟仿射影像的第i层非线性尺度空间的pc特征图在x方向上的一阶导数；

87、表示子模拟仿射影像的第i层非线性尺度空间的pc特征图在y方向上的一阶导数。

88、最后，对子模拟仿射影像的第i层非线性尺度空间的特征点度量矩阵mi进行特征分解，得到子模拟仿射影像的第i层非线性尺度空间的特征点度量矩阵mi的第一特征值λ1和第二特征值λ2。子模拟仿射影像的第i层非线性尺度空间特征点强度图为：ri＝min(λ1,λ2)，满足ri＞r的点则为子模拟仿射影像的第i层非线性尺度空间的特征点，r是一个阈值常数。

89、在得到的子模拟仿射影像的每层非线性尺度空间的特征点中，如果不同层的非线性尺度空间在相同的位置存在2个以上的特征点，则这些特征点只保留一个；如果不同层的非线性尺度空间在不同的位置存在特征点，则这些特征点全部保留，得到子模拟仿射影像的特征点集合

90、以此类推，对模拟仿射影像s′中其余的子模拟仿射影像进行特征点检测。

91、将所有子模拟仿射影像的特征点集合进行组合，得到模拟仿射影像s′的特征点集合

92、参照得到模拟仿射影像s'的特征点集合的方法，得到正射影像r的特征点集合keyr＝{(x1′,y1′),...,(x'm,y'm)}。

93、步骤2.3、特征描述算子构建

94、对模拟仿射影像s'中的任一特征点(xc,yc)，提取以所述特征点(xc,yc)为中心的圆形区域ω的幅值magω和方向angω：

95、

96、

97、式(9)～(10)中：

98、ω表示以所述特征点(xc,yc)为中心的圆形区域，圆形区域ω的半径为dr；

99、magω表示以所述特征点(xc,yc)为中心的圆形区域ω的幅值；

100、表示以所述特征点(xc,yc)为中心的圆形区域ω的pc特征图在x方向上的一阶导数；

101、表示以所述特征点(xc,yc)为中心的圆形区域ω的pc特征图在y方向上的一阶导数；

102、angω表示以所述特征点(xc,yc)为中心的圆形区域ω的方向；

103、ξ表示一个极小的正实数。

104、由式(9)～(10)得到所述特征点(xc,yc)所在圆形区域ω的幅值magω和方向angω，将方向angω范围设置为[0,360°)，再以每10°为间隔划分为36个区间。统计每一区间的幅值特征和方向特征，得到方向直方图分布；选择方向直方图的峰值方向作为所述特征点(xc,yc)的主方向。

105、用gloh-like对数极坐标构建所述特征点(xc,yc)的特征描述算子：将所述特征点(xc,yc)所在圆形区域ω的极径方向划分nr份，极角方向划分为nθ份，得到(nr-1)×nθ+1个对数极坐标子区域。统计每个对数极坐标子区域的nbin个方向幅度特征和角度特征，得到所述特征点(xc,yc)的特征描述算子。以此类推，得到模拟仿射影像s'上的其余的特征点对应的特征描述算子。

106、每个所述特征点对应的特征描述算子的维度为nbin×((nr-1)×nθ+1)。

107、再将模拟仿射影像s'上所有特征点的特征描述算子进行组合，得到模拟仿射影像s'的特征描述算子dess'。

108、然后参照所述得到模拟仿射影像s'的特征描述算子dess'的方法，得到正射影像r的特征描述算子desr。

109、步骤2.4、特征配准和误差剔除

110、先对模拟仿射影像s′的特征描述算子dess′和正射影像r的特征描述算子desr采用最近邻距离比率算法进行初始配准，用随机抽样一致性算法进行快速异常值剔除，得到模拟仿射影像s'和正射影像r的初始配准对应点；再用快速一致性采样方法剔除初始配准中的错误配准，当所述对应点的残差小于3个像素，则为模拟仿射影像s′和正射影像r正确配准的同名点。

111、步骤3、图像校正

112、设模拟仿射影像s′中正确配准的特征点(xp,yp)和正射影像r中正确配准的特征点(x′p,y′p)为一对同名点，将模拟仿射影像s′中正确配准的特征点(xp,yp)按照对应的倾斜度t和经度方向旋转角度找到与大倾角影像s中对应的特征点得到大倾角影像s和正射影像r的同名点为和(x'p,y'p)。以此类推，得到大倾角影像s和正射影像r之间的所有同名点。

113、然后通过所有同名点构建大倾角影像s和正射影像r之间的仿射变换参数ψ，将大倾角影像s和正射影像r按照仿射变换参数ψ进行校正，完成配准过程。

114、由于采用上述技术方案，本发明与现有技术相比具有如下积极效果：

115、本发明在步骤1采用asift算法模拟大倾角航空影像的仿射变换视角，能更好地适应大角度旋转问题。本发明在步骤2.1中采用非线性扩散滤波器构建非线性尺度空间，和传统高斯金字塔尺度空间相比，非线性尺度空间能在平滑图像噪声的同时保留图像边界纹理特征，更好地应对大倾角影像存在的地物几何变形、分辨率差异、地物相互遮挡的问题。

116、本发明在步骤2.2中用pc取代传统梯度信息，将pc与sobel算子结合构建特征点度量矩阵，应对航空影像地物信息复杂的问题，提取更加丰富且稳定的航空影像结构特征和点特征；在步骤2.3中，用pc和sobel算子计算特征点的幅值和方向，然后构建gloh-like特征描述算子。和传统sift-like特征描述算子相比，基于pc的幅值和方向比基于梯度的幅值和方向对图像噪声和明亮度变化具有不变性，能更好地表征图像结构特征和方向信息；另外，圆形结构gloh-like特征描述算子比方形结构sift-like特征描述算子更具鲁棒性和稳定性，使得配准误差小，配准成功率高，配准结果稳定。

117、因此，本发明具有配准误差小、抗噪性能好、成功率高、结果稳定的特点。