异源遥感影像匹配的实质是识别不同传感器、不同时间拍摄的两幅或多幅影像同名点的过程^[1]。随着多源传感技术和摄影测量学的迅速发展，异源影像匹配逐渐成为研究热点，并被广泛应用于图像融合、变化检测、三维建模、目标识别、灾害评估、违章建筑检测和土地资源变化监测等多个领域。

专家学者对遥感影像匹配展开了广泛研究，提出了基于互信息相似性匹配、局部相似性匹配及形状属性相似性等^[2-4]多模态影像匹配方法，通过提高多模态影像的相似性来提升匹配效果，然而都不同程度上受图像轮廓或形状依赖约束。有学者相继提出合成孔径雷达尺度不变特征变换（synthetic aperture radar scale invariant feature transform，SAR-SIFT）^[5]、位置尺度定向不变特征变换（position scale orientation-SIFT，PSO-SIFT）^[6]、多模式匹配等方法^[7-9]，虽然这些方法对遥感影像具有较好的匹配效果，但仍然依赖影像梯度特征信息。梯度特征对光照差异、对比度差异和非线性辐射畸变较为敏感，会导致匹配困难。有学者将相位一致性模型应用于异源影像匹配，提出了相位一致性直方图^[10]、局部相位特征描述^[11]、辐射不变特征变换（radiation-variation insensitive feature transform，RIFT）^[12]、Log-Gabor滤波新型扩展相位一致性^[13]、迁移优化^[14]及数据驱动灰狼算法优化^[15]等方法，利用相位特征及优化算法在异源影像中表现出了良好的效果，但相位一致性直方图和局部相位特征描述等方法受地理信息的影响较大且计算成本较高，RIFT方法不支持影像尺度差异，新型扩展相位一致性算法又对影像几何结构信息较为依赖。因此，这些方法虽然一定程度上改善和优化了异源匹配问题，但光照差异、对比度差异和非线性辐射畸变等匹配难题依然有待研究。深度卷积特征模型和相似性监督等^[16]基于深度学习的方法也被用于异源影像匹配，对异源影像匹配的部分环节进行监督匹配，不断丰富匹配方法。但深度学习方法的计算复杂度和迁移性仍有待提高，这些方法多集中在医学领域多模态匹配研究中，针对复杂性更大的遥感领域研究相对较少。

综上所述，从经典特征方法到相位一致性模型再到深度学习匹配方法都取得了较快发展。随着技术需求的日益突破，对异源匹配的要求也在不断提高。但异源匹配依然面临一些挑战，主要存在以下两方面问题：（1）异源图像因拍摄模式的不同和成像机理的差异导致特征检测难度较大，间接增加了匹配难度；（2）异源影像显著的光照差异、对比度差异和非线性辐射差异对影像的梯度和方向更加敏感，降低了特征描述的稳健性，容易导致匹配失败。

因此，本文从异源特征检测和稳健特征描述的角度来改善异源影像匹配困难的问题，提出了一种顾及各向异性加权力矩与绝对相位方向梯度特征直方图（histogram of absolute phase consistency gradients，HAPCG）的匹配方法，从而弥补异源匹配特征检测和传统影像特征梯度对异源影像光照差异、对比度差异和非线性辐射畸变较为敏感的缺陷，实现异源影像的稳健匹配。

1 绝对相位方向梯度直方图匹配方法

本文提出的匹配方法流程如图 1所示，主要包括4个步骤：（1）通过各向异性滤波扩算计算非线性尺度影像空间，然后借助相位一致性模型中的Log-Gabor函数奇对称滤波器生成最大矩和最小矩，构建各向异性加权力矩方程，求解获取各向异性加权力矩图。（2）在各向异性加权力矩空间中通过Harris函数获取极值点。（3）将相位一致性特征值作为影像梯度特征，将Log-Gabor函数的奇对称滤波器计算的相位特征方向进行扩展，生成绝对相位一致性方向特征，并联合对数极坐标描述框架来计算特征向量，重复迭代计算得到最终的HAPCG描述符。（4）采用欧氏距离作为匹配测度，通过HAPCG描述符最近邻和次近邻距离之比得到匹配同名点，并利用快速样本共识算法完成误匹配剔除。HAPCG算法的代码下载链接为：https://skyearth.org/publication/project/HAPCG/。

图  1  本文方法流程图

Figure  1.  Flowchart of Our Proposed Method

下载: 全尺寸图片 幻灯片

1.1   各向异性加权力矩构建

异源影像（尤其红外影像/电子地图/深度图）之间存在的光照、对比度差异和非线性辐射畸变等问题，进一步提升了特征点的识别难度。相位一致性模型中的最大最小矩可用于提取影像的边缘和角点特征。因此，本文将异源影像在各向异性滤波中扩散，并通过构建加权力矩方程来生成各向异性加权力矩图并用于检测特征点。

各向异性滤波是一种非线性滤波，可以更好地保留图像边缘信息，方便影像结构特征提取，一定程度上增加特征点的丰富程度。因此，采用文献[17]优化的图像非线性扩散方法，建立各向异性扩散方程为：

式中，div表示散度操作符；t表示时间度量的尺度值；$\mathrm{1}/(\mathrm{1}+{\left|\nabla L_{\sigma }\right|}^{\mathrm{2}}/k^{\mathrm{2}})$表示扩散系数；k表示对比度因子；$L$为差分图像；I表示当前图像；${\mathit{A}}_{\mathrm{1}}(L^n)$和${\mathit{A}}_{\mathrm{2}}(L^n)$分别表示编码图像$L^n$行方向和列方向扩散系数矩阵；$L^{n+\mathrm{1}}$表示扩散后的结果；$\tau$表示扩散的时间步长。其中在扩散计算时，需要将尺度值转换为时间值，即$t_n=\mathrm{1}/\mathrm{2}{\sigma }_n^{\mathrm{2}}$，σ${}_n$表示尺度，层数n属于[0，1，2…S]，S表示最大图层数目。

完成各向异性影像扩散后，进行相位一致性计算得到相位图^[18]，即P图。给了更好描述边缘特征，为每个方向o计算一个独立的相位图映射P（θ_o），θ是方向o的角度。根据力矩分析算法将最小矩对应的轴称为主轴，主轴通常表示特征的方向信息。最大矩的轴垂直于主轴，反映特征的显著性。本文根据最大矩和最小矩构建了加权力矩方程，归一化后生成各向异性化加权力矩结果，计算如下：

式中，W表示最终的各向异性加权力矩结果；${\mathit{M}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$表示影像的相位一致性最大矩；${\mathit{M}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$表示影像的相位一致性最小矩；$\mathrm{\Delta }\Phi$表示影像权重系数（取值在-1~5之间）；A、B和C是相位力矩计算的中间量。

1.2   特征点提取

通过各向异性力矩图进行Harris特征点提取^[19-21]，进行重复性检测，将过滤后的结果保留为特征点。为了进一步对比基于各向异性加权力矩图进行特征点检测的效果，对一组深度图和光学影像分别进行了高斯影像尺度空间的特征提取和各向异性力矩图的特征提取，结果如图 2所示。由图 2可知，各向异性加权力矩图能够显著提升检测得到的特征点数量。

图  2  高斯影像尺度空间的特征点提取和各向异性加权力矩图特征点提取结果对比

Figure  2.  Comparison of Results of Feature Points Extraction in Gaussian Image Scale Space and Anisotropic Weighted Moment Map

下载: 全尺寸图片 幻灯片

1.3   绝对相位方向特征描述

如图 3所示，绝对相位方向特征描述主要包括以下3个步骤：（1）计算相位一致性模型，将相位一致性特征值转换为影像梯度特征；（2）对Log-Gabor函数奇对称滤波器计算的相位特征方向进行扩展，生成绝对相位一致性方向特征；（3）联合对数极坐标描述框架计算特征向量，重复迭代计算得到HAPCG描述符。

图  3  绝对相位方向特征描述子流程图

Figure  3.  Flowchart of Absolute Phase Orientation Feature Descriptor

下载: 全尺寸图片 幻灯片

相位一致性是一种基于频率域的特征检测方法，能够较好地在影像傅里叶谐波分量叠加最大相位处检测出影像角点和边缘等特征，且几乎不受信号特征幅值大小干扰，因此对影像光照和对比度具有很好的不变性。

1.3.1   相位梯度幅值

在异源影像中考虑到诸如SIFT等传统匹配方法的梯度幅值和梯度方向对光照差异显著、对比度显著等更加敏感，相位一致性模型又具有显著优势。因此，本文将相位一致性模型计算的相位特征幅值特征来替代影像梯度幅值特征。

文献[20]利用Log-Gabor滤波改进了相位一致性的计算模型，并将相位一致性扩展到二维空间，计算如下：

式中，C（x，y）表示相位一致性模型的相位特征结果，本文将其作为影像梯度幅值；w_o（x，y）为频率扩展的权重因子；${\mathit{A}}_{s\mathit{o}}(x,y)$为像点（x，y）在Log-Gabor滤波器尺度s和方向o上的振幅；$\xi$是一个避免分母为零的极小常数；$⌊\cdot ⌋$防止值为负数；T为噪声阈值；$\mathrm{\Delta }{\Phi }_{s\mathit{o}}(x,y)$为二维相位偏差函数。

1.3.2   绝对相位一致性特征梯度方向

相位一致性模型只能计算影像的特征值，仅采用相位特征幅值进行特征描述不能实现稳健匹配。文献[11]通过Log-Gabor函数的奇对称滤波器构建相位一致性的方向特征，虽然取得了较好的效果，但受地理坐标信息影响较大。

特征方向表示影像特征变化最剧烈的方向，对于描述子的构建非常重要。Log-Gabor函数的奇对称滤波器${\mathit{C}}_{s\mathit{o}}(x,y)$的卷积结果表示影像在某个方向的能量变化。因此，可以用Log-Gabor函数的奇对称滤波器${\mathit{C}}_{s\mathit{o}}(x,y)$构建相位一致性的方向特征。

由于Log-Gabor奇对称滤波器在计算时需要将多个方向的奇对称滤波器结果分别投影到水平方向（X）和垂直方向（Y）获取水平方向的能量A和垂直方向的能量C，计算如下：

式中，$\mathit{\Phi }$表示相位一致性方向特征；${\mathit{C}}_{s\mathit{o}}(\theta )$表示方向θ上的奇对称滤波器卷积结果；$\xi$是一个极小值，防止分母为零。

然而，由于计算结果是多个方向的卷积值，导致相位方向特征值存在反方向，难以正确描述影像特征信息。因此，本文对相位一致性模型进行扩展，通过求取绝对值的方式将卷积结果全部归算到正方向上，之后通过非负常数项优化绝对相位一致性方向特征，获得绝对相位一致性方向特征${\mathit{A}}_{\mathit{o}}$，计算如下：

式中，$\mathit{\Phi }[i][j]$表示$\mathit{\Phi }$中的元素；$\mathrm{\Delta }\omega$表示一个常数（文中设置为90）。

1.3.3   对数极坐标描述

本文以特征点为中心选取固定的圆形邻域，并计算其绝对相位一致性梯度特征和方向特征，生成绝对相位一致性方向梯度直方图。将此直方图以每10°为间隔，均匀划分为36等份，统计每一等份的绝对相位一致性梯度特征和方向特征，选择直方图的峰值方向作为特征点的主方向。

获取特征点主方向后，利用对数极坐标框架来构建描述符^[20]。为了更加详细描述特征点邻域内的特征信息，首先将圆形邻域进行等分，形成了41个子区域的对极数坐标网格，各个子区域面积近似一致，然后统计每个子区域的像素点的8方向梯度特征和方向特征，最后生成包含328维特征向量的HAPCG描述子。

1.4   匹配与粗差剔除

本文采用欧氏距离作为相似性测度进行初始匹配来获取异源遥感影像的初始匹配，利用快速样本共识（fast sample consensus，FSC）方法来剔除初始匹配中的错误匹配^[21]，并对获取的同名点进行精度评估。

2 实验结果与分析

为验证HAPCG匹配方法的效果，将其与当前4种比较先进的SIFT、PSO-SIFT^[7]、Log-Gabor直方图描述符（Log-Gabor histogram descriptor，LGHD）^[22]和RIFT^[10]方法进行了比较，采用同名点匹配数量（number of correct matches，NCM）和均方根误差（root mean square error，RMSE）两种指标来衡量。参考已有研究，将特征点提取阈值设置为0.4，影像尺度差异设置为1.6，邻域窗口设置为42像素，其他4种匹配方法的参数也相应调整到最优阶段。本文采用仿射变换作为约束条件，将匹配误差在3像素内的结果认定为正确匹配，其中匹配同名点数目不能少于4对^[10]。上述方法都在Matlab R2018a下实现，实验平台处理器是Intel（R）Core（TM）i7-9750H CPU 2.59 GHz，RAM是16 GB，Windows10 X64操作系统。

2.1   影像数据

本文将18组具有光照差异、对比度差异、尺度差异、旋转差异、位移差异以及综合差异（同时具有对比度、尺度、旋转等多种差异情况）的异源遥感影像作为数据源（见图 4）。其中，图 4（a）的第一、第三组影像是从美国国家航空航天局（National Aeronautics and Space Administration，NASA）数据网站中获取，尺寸分别为600×600像素和500×500像素，图 4（a）的第二组影像是从Google Earth中获取，尺寸为500×500像素；图 4（b）的第一组和第三组影像从Google Earth中获取，尺寸都为500×500像素，第二组影像由无人机航摄获取，尺寸为640×512像素；图 4（c）的3组影像都是由无人机航摄获取，尺寸分别是600×450像素、600×450像素和600×337像素；图 4（d）的3组影像是从Landsat和Sentinel系列卫星影像中获取，第一、第三组影像的尺寸都是500×500像素，第二组影像的尺寸是588×606像素；图 4（e）的第一、第二组影像是从Google Earth中获取，尺寸都是500×500像素，第三组影像是由无人机航摄获取，尺寸为505×329像素；图 4（f）的3组影像是从NASA数据网站中获取，第一组影像的尺寸是600×600像素，第二组影像的尺寸是500×500像素，第三组影像的尺寸是516×500像素。

图  4  异源遥感影像数据

Figure  4.  Heterologous Remote Sensing Image Data

下载: 全尺寸图片 幻灯片

2.2   实验结果

为了验证本文提出的HAPCG方法的鲁棒性和匹配精度，分别从方法本身改进、其他方法对比和定量检测3个方面进行验证。

首先，将HAPCG方法分别与直接利用Sobel算子计算的影像梯度匹配方法、采用未改进的相位一致性匹配方法进行对比。3种方法均在各向异性加权力矩图的基础上进行后续计算，结果如图 5所示。由图 5可知，在异源影像中采用3种方法构建的不同方向特征匹配结果差异显著，其中本文提出的绝对相位一致性方向特征匹配结果最佳。Sobel算子计算的影像方向特征虽然能较好地体现影像的结构特征，但由于影像存在明显的光照和对比度差异，难以获取到正确的相似性特征，导致匹配失败。图 5（b）和图 5（c）均采用将相位一致性模型的方向特征来代替影像梯度方向特征，虽然实现了匹配，但本文方法计算的相位方向特征更加丰富，相似性程度更高，因此能得到更加稳健的匹配结果，识别的正确同名匹配点数相比直接采用相位一致性方向特征识别的结果，其效果提升了约3倍。

图  5  不同方法的匹配结果对比

Figure  5.  Comparison of Matching Results Using Different Methods

下载: 全尺寸图片 幻灯片

为进一步验证本文提出的绝对相位方向特征匹配方法的效果，将其分别与SIFT、PSO-SIFT、LGHD、RIFT 4种匹配方法进行对比，选择6组异源遥感影像进行匹配对比实验，结果如图 6所示。由图 6可知，SIFT方法在夜间灯光影像和可见光影像以及大尺度差异的高光谱影像和可见光影像中能成功匹配，但在电子地图、红外影像与可见光影像匹配中效果较差（图 6（a））。而PSO-SIFT方法不能实现全部匹配，在红外与可见光影像中匹配效果最好（图 6（b）），但在大对比度影像中匹配效果不佳，甚至失败。LGHD方法虽然全部匹配成功，但在大对比度的导航地图和可见光影像中检测特征点难度较大，匹配效果不佳（图 6（c））。尽管SIFT/PSO-SIFT/LGHD 3种方法在部分影像中实现了成功匹配，但识别的匹配同名点数量较少，大部分都低于50对。RIFT方法在对比度差异、旋转差异和位移差异中具有较好的性能，但因其不支持尺度差异，因此在具有尺度差异的数据中表现不佳（图 6（d））。本文提出的HAPCG方法在18组匹配实验中均实现了成功匹配，且能识别到充足的匹配同名点，NCM数目都在100对以上，图 6（e）是6种情况中各自第一组像对的匹配结果。综上所述，HAPCG方法能在异源遥感影像中实现稳健匹配，与其他方法对比，能识别到的匹配同名点数量最多。同时，为进一步展示HAPCG方法的匹配效果，图 7展示了另外12组异源遥感影像的匹配结果。

图  6  5种方法匹配结果对比

Figure  6.  Comparison of Matching Results of Five Methods

下载: 全尺寸图片 幻灯片

图  7  12组影像的HAPCG方法的匹配结果

Figure  7.  Matching Results of HAPCG Method for 12 Groups of Images

下载: 全尺寸图片 幻灯片

匹配算法优劣也可以通过同名点的RMSE来进行定量检测，计算如下：

式中，N表示选取的同名点个数；（x_i，y_i）表示参考图像中与待匹配图像中的同名点坐标；（$x_i^{″}$，$y_i^{″}$）是第i个待匹配图像中同名点在经过匹配对应关系转换后的坐标。

表 1统计了5种方法在18组异源遥感影像匹配测试中的NCM和RMSE两个指标，定量地对比分析HAPCG方法的性能。由表 1可知，实验的18组数据在同名点对数量方面，HAPCG方法明显多于SIFT、PSO-SIFT和LGHD方法。SIFT方法效果最差，18组影像中只匹配成功了10组，且NCM大部分都低于50。PSO-SIFT方法优于SIFT，但在电子地图与光学影像匹配中效果较差；LGHD方法在18组影像中可以成功匹配14对异源影像，但在具有综合差异的异源影像中效果不佳，且NCM和RMSE均低于HAPCG方法。RIFT方法在对比度差异、旋转差异和位移差异中效果较好，也能适应较小尺度差异，但在尺度差异显著和综合差异匹配中性能受到较大影响，且在平均NCM及RMSE方面不及HAPCG方法。HAPCG的平均NCM约是SIFT方法的4倍，PSO-SIFT方法的3.5倍，LGHD方法的2倍，且HAPCG在18组影像中的平均RMSE是1.83像素，整体RMSE精度要优于其他4种方法。

表  1  5种方法在18组异源遥感影像中的定量匹配结果

Table  1.  Quantitative Matching Results of Five Methods in 18 Groups of Heterologous Remote Sensing Images

差异条件	影像	指标	SIFT	PSO-SIFT	LGHD	RIFT	HAPCG
光照差异	第一组	NCM	54	43	41		260
		RMSE/像素	2.458	2.044	1.965		1.726
	第二组	NCM	78	117	400	365	484
		RMSE/像素	2.306	1.924	1.997	2.046	1.931
	第三组	NCM	18	2	275	312	304
		RMSE/像素	1.360		2.070	2.047	1.893
对比度差异	第一组	NCM	3	2	44	594	465
		RMSE/像素			1.939	1.876	1.742
	第二组	NCM	6	39	7	62	163
		RMSE/像素		1.874		1.903	1.885
	第三组	NCM	5	31	31	332	408
		RMSE/像素		2.046	1.487	2.034	1.860
尺度差异	第一组	NCM	23	37	13	50	288
		RMSE/像素	1.867	1.906	1.318	1.761	1.755
	第二组	NCM	25	42	8		243
		RMSE/像素	2.025	2.080	2.028		1.804
	第三组	NCM	6	158	7		180
		RMSE/像素		1.878			1.928
旋转差异	第一组	NCM	5	203	295	213	424
		RMSE/像素		1.687	1.965	1.939	1.750
	第二组	NCM	3	234	593	237	302
		RMSE/像素		1.782	2.010	1.912	1.855
	第三组	NCM	52	198	192	670	686
		RMSE/像素	1.675	1.734	1.819	1.935	1.873
位移差异	第一组	NCM	3	15	25	460	379
		RMSE/像素		2.093	1.806	1.86	1.796
	第二组	NCM	80	179	192	608	810
		RMSE/像素	1.479	1.598	1.929	1.849	1.768
	第三组	NCM	228	125	95	207	540
		RMSE/像素	1.669	1.767	1.918	1.924	1.750
综合差异	第一组	NCM	90	33	56		148
		RMSE/像素	1.994	2.034	2.060		1.988
	第二组	NCM	5	1	9		118
		RMSE/像素					1.695
	第三组	NCM	75	36	40		102
		RMSE/像素	2.401	2.031	2.012		1.944

下载: 导出CSV 

| 显示表格

通过计算匹配时间成本来进一步检验HAPCG方法的性能，表 2统计了5种匹配方法在18组异源遥感影像中的时间消耗情况。由表 2可知，5种匹配方法的匹配时间成本具有显著差异。尽管SIFT方法的匹配时间最短，但其识别到的同名点大部分都少于50对，且18组影像中只匹配成功了10对，成功率仅有55.56%，因此SIFT方法的耗时不能真实反映异源遥感影像匹配的时间成本。PSO-SIFT方法和LGHD方法虽然能满足大部分像对的匹配任务，但是LGHD匹配时间成本较高，严重影响了匹配效率，PSO-SIFT的匹配时间虽然优于LGHD，但匹配时间成本仍高于RIFT和HAPCG方法，其匹配效率相对较低。RIFT方法虽然在获取较多同名点的同时，时间成本优于PSO-SIFT和LGHD方法，但其匹配适用场景受到较大约束，尤其在尺度差异较大和综合差异情况下表现明显。HAPCG方法能获取较佳的匹配成本，在保证识别出丰富的同名点对的同时，依然能实现较高效率的计算，匹配时间成本总体上优于另外4种方法。

表  2  5种方法在18组异源遥感影像中匹配时间的统计结果

Table  2.  Statistical Results of Matching Time of Five Methods in 18 Groups of Heterologous Remote Sensing Images

差异条件	影像	耗时/s
		SIFT	PSO-SIFT	LGHD	RIFT	HAPCG
光照差异	第一组	26.830	122.07	81.11	13.14	13.56
	第二组	4.080	10.84	36.07	11.78	9.30
	第三组	4.035		33.72	12.78	7.94
对比度差异	第一组	2.920		20.70	11.61	8.05
	第二组	7.340	24.87	43.55	14.04	10.92
	第三组	3.390	10.95	20.76	13.68	10.84
尺度差异	第一组	3.490	10.82	28.26	16.25	8.15
	第二组	5.100	12.30	36.16	15.14	7.70
	第三组	4.790	14.60	48.87	11.53	8.27
旋转差异	第一组	6.050	14.35	57.12	11.5	9.71
	第二组	11.470	24.08	147.37	12.77	13.67
	第三组	4.730	12.61	19.35	14.26	9.05
位移差异	第一组	2.560	5.33	19.95	13.09	6.99
	第二组	3.670	20.35	29.85	14.48	10.75
	第三组	4.860	11.60	17.10	10.41	6.30
综合差异	第一组	15.310	45.50	105.11	14.32	13.84
	第二组	5.480	10.27	53.11	13.4	9.48
	第三组	5.440	11.88	53.71	13.41	10.03

下载: 导出CSV 

| 显示表格

3 结语

为解决异源遥感影像因光照、对比度差异及非线性辐射畸变所导致的匹配困难问题，本文提出了一种顾及各向异性加权力矩和绝对相位方向特征描述的稳健匹配方法。通过18组具有光照差异、对比图差异、尺度差异、旋转差异、位移差异和综合差异的异源影像进行了实验验证，并与SIFT、PSO-SIFT、LGHD和RIFT方法进行了对比。实验结果表明，所提HAPCG方法具有较好的尺度不变性，在保证匹配精度的前提下能获取更多的正确同名点，其平均NCM和平均RMSE均优于其他4种方法。HAPCG方法在异源遥感影像中能实现稳健、鲁棒匹配，可为多源影像数据融合等任务提供稳定的数据支撑，尤其在红外影像、电子地图、多光谱影像等数据类型中效果更为突出。同时，本文方法在特征点提取阈值设置上可以更加灵活，取值在0~1之间。当需获取更多特征点时调小阈值即可；反之，则调大阈值，常规情况下采用0.4即可满足匹配需求。