基于双特征的丘陵山区耕地低空遥感图像配准算法

宋飞， 杨扬， 杨昆， 张愫， 毕东升

(1. 云南师范大学信息学院， 昆明 650500； 2. 西部资源环境地理信息技术教育部工程研究中心， 昆明 650500；3. 云南师范大学模式识别与人工智能实验室， 昆明 650500)

丘陵山区是中国南方主要的地貌结构[1]，其地面崎岖不平、岗拗交错，又常年受多雾、多雨气候条件的影响，导致耕地数量减少，水土流失加重，以及土壤环境污染加剧等现象。严重威胁着丘陵山区的农业生产、粮食安全和可持续发展[2-3]。在此严峻背景下，如何及时提取耕地分布信息，对于改善丘陵山区居民生活水平至关重要。

航天、航空遥感技术是目前提取耕地分布信息的重要方法[4]。但现有大、中型的航天、航空遥感技术存在气象影响因子多、重访周期过长、应急不及时和成本高等缺点。因此，采用小型无人机进行耕地分布信息的遥感获取，具有机动灵活、高时效、高分辨率、能耗低、成本低和受天气状况影响小等优点[5]。在小型无人机获取遥感图像(低空遥感图像)的过程中，由于受飞行器姿态、高度、速度等因素影响，经常会导致获取的遥感图像相对于地面目标实际位置发生挤压、扭曲、伸展、偏移以及图像重叠度不规则等一系列问题。消除这些因素所造成的同一场景(如不同视角的同一区域)不在同一坐标系(相同尺度的空间坐标系)的情形，是顺利实现航拍数据分析的前提和保证。

针对以上问题，遥感图像配准方法可以依据一定的相似性度量，对在不同条件下获取的两幅或多幅遥感图像进行对准、叠加的操作，使得图像在同一坐标系下获得最佳的重合。遥感图像配准方法主要分为2类：基于区域的配准方法和基于特征(例如点、线、面)的配准方法。由于点特征具有光照、旋转不变，可靠性高，且易提取性、计算简单、运算速度快，因此，选择图像中点特征研究的最多，典型的点特征提取算法有Harris算子[6]、SIFT(Scale-Invariant Feature Transform)算子[7]、SURF(Speeded Up Robust Features)算子[8]等。其中，Bay等[8]提出的SURF算子，不仅具有尺度、光照不变性，而且相比于SIFT算子计算量小，处理速度快而受到了领域内学者的青睐。

在基于特征点的图像配准方法研究领域，Brook 和Bendor[9]通过在SURF算子中加入拓扑关系提高了遥感图像的配准精度；Zhao等[10]提出了基于SIFT和区域互信息优化的无人机遥感图像方案，提高了配准的精度；Lei等[11]通过在SIFT算法中加入随机抽样一致性和最小二乘法提高了无人机遥感图像自动配准的精度；乔川等[12]针对传统图像配准方法难以对海洋、沙漠和草原等特征不明显区域航空遥感图像进行配准的问题，提出了一种基于地理位置信息的航空遥感图像配准算法；Myronenko和Song[13]提出了可以同时处理刚性和非刚性配准问题的CPD(Coherent Point Drift)算法，可以在牺牲些许配准精度的前提下，通过采用快速高斯变换[14]和矩阵低秩逼近[15]技术减少计算量来提升算法的配准速度。Wang等[16]提出了一种基于非对称高斯模型的非刚性配准算法(AGMReg)，在再生核希尔伯特空间(RKHS)决了正则化理论下的优化问题。Yang等[17]提出了一种基于欧氏距离的全局特征描述子和局部距离特征描述子的混合特征描述子的非刚性配准(GLMDTPS) 算法，大幅提升了非刚性点阵配准在处理冗余点和噪音等方面的鲁棒性。Ma等[18]提出的PRGLS通过匹配形状上下文(Shape Context，SC)描述子来获取一个二进制对应矩阵，如果分配高斯混合模型(Gaussian Mixture Model，GMM)的成员概率接近1，则表示匹配，从而来改善CPD的配准性能。

本文在综述图像配准方法研究现状的基础上，针对丘陵山区耕地多视角航拍图像中的尺度变化、几何畸变和图像重叠等问题，开展了图像特征提取、匹配等相关算法和技术的研究。主要研究工作以及贡献：首先，针对耕地复杂的空间分布和几何形状，本文构建了能够稳健描述航拍图像几何特征的双特征描述子(Dual Feature Descriptor，DFD)，此描述子由基于欧氏距离的全局特征和基于和向量的局部特征构成。在图像配准过程中，通过使用不确定性退火方案合理调节2个特征的强度。其次，采用特征点集之间欧氏距离度量的匹配策略最为常用，但匹配精度不高。为了提高低空遥感图像的匹配精度，本文以GMM为核心，结合构建的双特征描述子，得到了能够同时通过2种特征进行对应关系评估的双特征有限混合模型。最后，为了维护特征点集在进行空间变换更新时的局部与全局结构稳定，构建了基于L2E[19-21]的双约束(局部、全局约束项)能量方程来估计变换函数。

1 方 法

首先介绍本文算法的3个贡献：构建双特征描述子、基于双特征有限混合模型、基于L2E的双约束能量方程。其次，详述具体的算法步骤。最后，进行算法分析。假定D维上的2个点集YM×D=[y1,y2,…,ym}T、XN×D=[x1,x2,…,xn]T(待配准图像特征点集Y、参考图像特征点集X)，其中ym和xn分别表示为Y的第m个特征点和X的第n个特征点。

1.1 构建双特征描述子

(1)

式中：γ为一个权重参数，用来调节μmk与Lm之

图1 心型特征点集中某中心点与其最近相邻点(5个点)构成的局部片段Fig.1 A central point and its nearest neighboring points (five points) construct a small fragment of heart-shaped feature points

(2)

式中：I为M×M对角矩阵；RY为一个M×M的矩阵；定义rij表示RY中任意一项，yi和yj表示Y中任意一点，且

(3)

DFDY=GY+εLY

(4)

1.2 基于双特征有限混合模型

Myronenko和Song[13]在CPD方法中提出2个特征点集对应关系评估可以视为概率密度估计问题。本节针对特征点集配准问题采用GMM进行建模，结合双特征描述子对混合模型的参数进行赋值，从而改进对应关系评估的效果。双特征有限混合模型的推导过程，具体如下：

1) 构建GMM。将整个参考图像特征点集X拟合成一个GMM，使其高斯模型(Gaussian Model，GM)的中心与待配准图像特征点集Y保持一致，故GMM的概率密度函数(Probability Density Function，PDF)为

(5)

式中：P(m)=1/M为混合系数，描述了第m个高斯模型在混合模型中的权重。

(6)

(7)

2) 构建双特征有限混合模型(DFMM)。若使用式(7)进行对应关系评估，只考虑点集X到点集Y之间的欧氏距离，可能导致配准过程的鲁棒性不好，例如：xa和xb的局部结构不同，但xa与xb到ym的欧氏距离一致，导致pma与pmb相等。

为了解决以上问题，将双特征描述子式(4)代替方程式(7)中的欧氏距离，重写式(7)的形式得到DFMM，并将其定义如下：

(8)

1.3 基于L2E的双约束能量方程

从全局(欧氏距离)的角度结合局部特征(和向量)对点集空间分布和结构进行分析，为点集匹配问题提供了一条新颖而具有潜力的解决途径。但该方法不可避免的会产生误匹配，即异常点。因此，在空间变换更新中，不得不对变换函数T(ym)采取鲁棒的估计。本文构建基于L2E[19-21]的双约束能量方程来估计变换函数T(ym)，其表达式为

E(T,σ2)=L2E(T,σ2)+(T)

(9)

1) 构建L2E估计子。假设有参数模型P(z|θ)，本文目标是对L2距离进行最小化，即

(10)

(11)

全局结构约束项，具体能量方程形式为

(12)

(13)

式中：β为一个常数，用来控制空间平滑性。

局部结构约束项，核心思想是通过判断某图像在空间变换前后的局部相似性,来约束其在进行空间变换时的局部形变，方程形式如下：

(14)

1.4 图像配准过程

针对以上研究，本文提出了基于双特征的丘陵山区耕地低空遥感图像配准的算法——DFMM算法，其流程图如图2所示。

图2 图像配准流程图Fig.2 Flowchart of image registration

1.4.1 基于SURF的特征点提取

鉴于丘陵山区耕地背景环境复杂、光照因素等影响，本文使用尺度不变特征SURF算法分别提取待配准遥感图像Is和参考图像Ir的特征点集YM×D=[y1,y2,…,ym]T、XN×D=[x1,x2,…,xn]T，D为特征点集的维度，M=N。具体步骤如下：

1) 对参考图像与待配准图像分别进行SURF特征点检测与定位和特征描述子计算。

2) SURF特征点匹配。本文采用马氏距离进行相似性度量，其定义如下：

(15)

1.4.2 特征点匹配

1) 对应关系评估

通过式(8)得到一个M×N矩阵P，便是待配准图像特征点集Y和参考图像特征点集X之间的一组模糊对应关系。然后由矩阵P得到对应的目标点阵Xc，其方程为

Xc=PX

(16)

2) 空间变换更新

在RKHS[24]中对变换函数T(ym)进行建模。首先，选择一个正定核G(y,ym)(即高斯核)来定义一个RKHSH。然后，采用GRBF，得到变换函数T(ym)最优解表达式为

(17)

式中：系数cm为一个D×1维的待求解向量。于是，在一个无限维H空间中最小化问题归结于求解一个有限的M系数向量cm。然后，采用能量方程式(9)评估变换函数T(ym)，具体形式为

(18)

(19)

式中：H=exp{d[(UC-PX)(UC-PX)T]/(2σ2)}为M×1的向量，d(·)为矩阵的对角线；11×D所有1的1×D行向量；⊙为Hadamard乘积；⊗为Kronecker卷积；A=[2KRx-2(RY)TRX+2K·(RY)T-2K2I]；B=[2(RY)TRY-2KRY-2K·(RY)T+2K2I]；R类似于二进制矩阵，只需将“1”的每个条目替换为高斯模糊的向量范数。

1.4.3 图像转换

本节采用“反推法(the backward approach)[25]”来建立TPS(Thin-Plate Spline)变换模型，以避免输出图像发生孔和(或)重叠的现象(由于离散化和四舍五入所引起的配准错误)。从待配准图像得来的转换图像是由目标像素坐标决定的(即与参考图像有着同样的坐标元)，并且估计的TPS变换模型参考图像到浮动图像是“反推法”的关键理念，可以将其定义为

It(x,y)=Is(T(x,y))

(20)

式中:It和Is分别为转换图像和待配准图像。It与参考图像Is具有相同的大小，且T(x,y)为从参考图像到待配准图像的TPS变换估计模型。T(x,y)的方程式为

(21)

(ω1,ω2,…,ωn,c1,cx,cy)T=

(22)

1.5 算法分析

1.5.1 SURF提取丘陵山区耕地低空遥感图像特征点的优势

小型无人机在航拍丘陵山区耕地时，由于受地面崎岖不平和航拍视角等影响，会加剧航拍图像间产生非刚性几何畸变，增加图像配准的复杂度。因此，有效地选择图像特征点集(优于随机选取)，会减小配准的计算复杂度。为了精确、鲁棒、快速地获得特征点集，本文开展了提取丘陵山区耕地低空遥感图像特征点集的研究(见图3)，通过实验比较，SUSAN算子提取的轮廓信息会增加匹配的计算复杂度;Harris算子不具有尺度不变性，而且在角点定位方面存在偏差，可能导致匹配不准确；SURF算子对平移、旋转、尺度变化和噪声等具有良好的不变性，对视觉变化、仿射变化也保持一定程度的稳定性，因而能够均匀提取到丘陵山区耕地低空遥感图像的特征点。

1.5.2 特征点匹配性能评估

1) 单一特征与双特征性能比较。

将只使用单一特征(欧氏距离)的CPD算法和使用双特征的DFMM算法的特征点匹配过程作对比(见图4)，在迭代开始时(如：Iter=1),CPD方法最初评估所得的对应点集整体位于质心区域，而使用DFMM算法最初评估，其分布更接近真正的目标点集。随着迭代次数增加(如：Iter=5,Iter=10)，CPD算法评估得到的点集依旧整团分布于区域质心周围，相比之下，DFMM算法迭代更新得到的源点集位置更优。

2)L2E的鲁棒性。

参数估计一般采用最大似然估计(Maximum Likelihood Estimation，MLE)来完成。然而，当采用的参考图像与待配准图像近似得不够好，或者提取的特征点集中存在大量冗余点时,MLE将会产生严重的偏差。L2E是一个用于最小化密度之间L2距离的鲁棒估计子，其与相比，更适用于分析包含冗余点的大规模数据集。将L2E和最大似然估计MLE进行对比，评估L2E鲁棒性。

图3 提取图像特征Fig.3 Image feature extraction

图4 单一特征与混合特征(双特征)性能比较Fig.4 Performance comparison of single feature and mixed feature (dual-feature)

图5 L2E与MLE随着冗余点数量的改变的鲁棒性对比Fig.5 Variation of L2E and MLE with number of redundant point and their robustness comparison

在4附近存在一个局部最小值，并且随着冗余点的增加极值会变得更深。与此相反，MLE的偏差随着冗余点的增加而变大，效果的稳定性变差。

1.5.3 参数设置

1) 退火参数。通过实验分析，式(4)在开始迭代时ε=exp(0)=1随着迭代次数的增加，ε逐渐减小，最后ε=exp(12.5)≈3×10-6。实验中，考虑到效率问题，将式(1)初始化γ参数设置为0.9。

2) 平滑性控制。在GRBF中，β用来控制空间平滑性。将源图像的特征点的空间坐标归一化为[-1.5，1.5]，因此β设置为2。

3) 协方差。在式(7)和式(18)中，σ2对非凸性问题处理上作用明显，退火方程σ2=ρσ2，其中ρ=0.9，σ2分别初始化为1和0.05。

4) 正则化参数。式(18)中有2个正则化参数λ和η，分别代表了全局和局部结构约束的权重大小。在本算法中，将λ和η分别设置为2。

5) 最近邻近点。在式(1)中，K的设定基于区别局部结构描述子所需的最少的点数量[17]。例如，区别角(包含2个相邻点)和十字(包含4个相邻点)，至少需要4个相邻点。通过实验分析，二维和三维配准中K的默认值设为5。

2 实 验

首先，使用丘陵山区耕地的小型无人机航拍多视角遥感图像(低空遥感图像)，验证本文提出的基于双特征的丘陵山区耕地低空遥感图像配准算法——DFMM算法；然后，验证本文算法在其他复杂地形的配准效果。

2.1 丘陵山区耕地低空遥感图像配准实验

2.1.1 研究区域及数据

坡耕地是丘陵山区的典型代表[26-28]，是分布在山坡上6°～25°之间的地貌类型，被开垦后形成的耕地。根据坡耕地的坡度(地表单元陡缓的程度)，本文实验数据集主要划分成2大类：(i)选取60张(30组)坡度为6°～15°(包括15°)，不同角度变化的坡耕地小型无人机多视角遥感图像；(ii)选取60张(30组)坡度为15°～25°，不同角度变化的坡耕地小型无人机多视角遥感图像。

本实验所有小型无人机多视角遥感图像来自720云平台，区域分别为云南、贵州、四川、广西、甘肃及湖南丘陵山区的典型坡耕地。每组图像的航拍视角变化为俯角变化(30°～90°)以及水平视角变化(-90°～90°)。其大小范围为640 mm×450 mm～1 100 mm×850 mm，实验数据集如表1所示。

表1 实验数据

2.1.2 评估标准

在评估检验时，首先，手动的在待配准图像和参考图像上确定了至少20个地面真值。为了减小评估过程中产生的误差，所有地面真值选择的都是合理分布于容易识别的区域。然后，使用均方根误差(Root of Mean Square Error，RMSE)、平均绝对误差(Mean Absolute Error， MAE)来量化配准精度。统计中的相关配方和定义如下：

(23)

(24)

平均绝对误差通过计算选定地标与实际位置对应点之间的平均绝对差异，来检验算法的效果。

2.1.3 实验结果及分析

为了证明本文算法的可行性及其性能，将DFMM与SIFT[7]、SURF[8]、CPD[13]、AGMReg[16]、GLMDTPS[17]及PRGLS[18]算法进行了比较。结果表明，本文算法在进行坡耕地的小型无人机多视角遥感图像(低空遥感图像)配准过程中，具有较好的鲁棒性。在图6、图7中分别列举出了5组坡度为6°～15°、15°～25°坡耕地的小型无人机多视角遥感图像(低空遥感图像)的匹配结果。图6和图7前2行表示两幅取自不同视角的无人机航拍图像(待配准图像和参考图像 )。每种配准方法用2行进行展示，首行表示特征点的配准结果，第2行中表示给定一个5×5的棋盘格，交替显示转换图像和参考图像。数据集(i)和(ii)的RMSE、MAE评估结果如表2所示。

结果分析：首先，分析以上5种算法的主要特征，SURF、SIFT特征点集的提取是基于强度信息，或者更确切地说，是尺度空间极值；CPD是采用欧氏距离，利用GMM建立对应关系；GLMDTPS主要利用基于欧氏距离的全局特征描述子和局部距离特征描述子的双特征描述子建立对应关系；AGMReg是以非对称高斯模型为核心，使用软分配技术来恢复对应关系；PRGLS基于全局特征和基于SC描述子来获取一个二进制对应矩阵，来完成图像相关匹配。DFMM主要以GMM为核心，结合2个单一特征差异描述子(基于欧氏距离的全局特征和基于和向量的局部特征)构建的双特征描述子(式(4)),得到了改进后的GMM，使其能够同时通过2种特征进行对应关系评估(式(8))。

图6 5组数据集(i)的图像配准结果示例Fig.6 Examples of image registration results for five sets of data sets (i)

图7 5组数据集(ii)的图像配准结果示例Fig.7 Examples of image registration results for five sets of data sets (ii)

误 差数 据 集SURFSIFTCPDAGMRegGLMDTPSPRGLSDFMMRMSE(i)7.283 712.528 75.436 55.112 34.919 13.716 31.321 1(ii)6.062 711.556 63.257 63.190 43.106 52.483 21.012 7MAE(i)4.234 4 7.288 93.110 23.044 03.058 12.038 01.074 9(ii)5.441 1 6.008 02.239 42.170 22.103 81.504 50.645 2

总结这5种算法中使用的主要特征包括：①强度信息差异；②全局几何结构差异；③局部几何结构差异；④全局几何约束；⑤局部几何约束。从表2中的实验结果容易看到CPD、AGMReg、及PRGLS配准结果优于SURF、SIFT，其原因在于SURF、SIFT只采用特征①提取对应的特征点完成图像配准，而CPD、AGMReg、GLMDTPS及PRGLS分别还采用了特征②④、②④、②③及②③④，能够在一定程度上提高图像配准的精度。本文DFMM算法配准结果优于其他配准结果，其主要原因是采用了特征①②③④⑤，不仅能够提取对图像尺度和旋转变化具有不变性、对光照变化和图像变形具有较强适应性的特征点，而且还能有效地利用了图像特征点的局部信息和全局信息完成图像之间的配准。

2.2 其他复杂地形低空遥感图像配准实验

复杂地形地质结构特殊，时常发生滑坡、山体崩塌等自然地质灾害,因此实时监测复杂地形非常必要。由于小型无人机具有机动灵活、高时效、高分辨率、能耗低、成本低、受天气状况影响小等优点，可以广泛地运用于复杂地形的监测。本节通过实验，验证了DFMM算法在不同复杂地形小型无人机航拍的多视角遥感图像(低空遥感图像)中的配准效果，所有数据来自720云平台。实验结果表明(表3)，DFMM算法具有较好的配准效果。在图8分别列举出了4组不同复杂地形的对应关系(图8(a))、匹配结果(图8(b))，图8(a)第1行表示不同视角的无人机航拍图像(待配准图像和参考图像)，第2行表示匹配的对应关系；图8(b)第1行表示特征点的配准结果，第2行表示给定一个5×5的棋盘格，交替显示转换图像和参考图像。

表3 使用RMSE和MAE进行复杂地形低空遥感图像实验的评估实验结果

图8 4组图像配准结果示例Fig.8 Examples of four sets of image registration results

3 结 论

采用小型无人机提取丘陵山区耕地分布信息，具有机动灵活、高时效、高分辨率、能耗低、成本低、受天气状况影响小等优点，但地面的起伏变化和航拍视角变化会加剧图像间产生非刚性几何畸变及重叠度不规则等问题。针对以上问题，本文提出了一种基于双特征的丘陵山区耕地低空遥感图像配准算法——DFMM算法。

实验验证表明，该算法可以有效用于解决丘陵山区耕地小型无人机低空遥感图像间存在的非刚性几何畸变及重叠度不规则等问题。DFMM算法不仅能够稳健描述航拍图像几何特征，而且配准结果的鲁棒性和精度有明显提高。这种基于双特征的低空遥感配准算法，也适用于部分复杂地形小型无人机航拍的多视角遥感图像配准(低空遥感图像)，将会进一步研究能否通用到所有复杂地形。