随着传感器技术的不断发展，人们通过卫星观测地球的手段越来越多，不同传感器获取的遥感数据各自具有独特优势，充分融合和利用多源遥感影像数据实现优势互补是一个重要的应用趋势。然而，多源影像数据融合利用的前提是影像间高精度配准，这个过程一般是通过在影像间匹配一定数量的同名像点，然后根据一定的数学模型，通过平差求解最佳的模型变换参数来实现的。因此，在异源遥感影像之间匹配获取可靠的同名点至关重要，而不同谱段、不同载荷获取的多源遥感影像之间匹配难度往往较大。

影像匹配是在两张或多张影像上识别同名点的过程^[1]。常用的遥感影像匹配方法可以分为两类：一类是基于特征的方法，其中尺度不变特征变换^[2]（scale-invariant feature transform，SIFT）算法是计算机视觉领域经典的特征匹配算法，广泛用于各类遥感影像匹配问题的解决过程。但是由于异源影像间的辐射差异大，传统的基于像素梯度的SIFT描述子无法获得稳定的同名点。为了克服影像间辐射差异的影响，许多学者在SIFT算法基础上进行改进，用于多源遥感影像匹配^[3-5]。Dellinger等^[6]和Xiang等^[7]用比值梯度取代灰度差梯度来对光学和合成孔径雷达（synthetic aperture radar，SAR）图像进行匹配，一定程度上减弱了斑点噪声的影响。叶沅鑫等^[8]使用局部相位描述符，在不同谱段的光学影像之间取得了比SIFT和加速稳健特征（speeded up robust features，SURF）更好的匹配效果。但重复特征的获取和表达仍是特征匹配需解决的问题。另一类是基于区域的方法，区域匹配方法一般通过滑动模板窗口来比较窗口之间的相似性从而确定同名点，此方法不受特征提取的影响。其中相关系数是影像匹配中使用最为广泛的相似性测度，互信息（mutual information，MI）^[9]方法利用了变量间的统计依赖关系，对图像的辐射差异不敏感，因此在多种类型的多模态图像匹配应用中取得了相当不错的效果^[10-12]。近年来，基于相位一致性^[13]构建的区域特征描述子在多源遥感影像匹配中研究较多，如叶沅鑫等^[14]结合相位一致性和方向梯度直方图构建的相位一致性方向直方图（histogram of orientated phase congruency, HOPC）描述子和Fan等^[15]利用多尺度加权相位一致性构建的相位一致性结构描述子（phase congruency structural descriptor, PCSD）等。这些基于区域的匹配方法能够在特征重复率明显降低的多源遥感影像间获取更多的可靠同名点。

鉴于卫星遥感平台的特性，相机的成像模型参数和位置信息很容易获取，可以利用已有的成像参数和地理参考信息来确定影像和影像之间的尺度、方向和定位差异，这是遥感影像匹配的重要初始信息。本文针对异源遥感影像的上述特点，在经过初步纠正的影像上提取均匀分布的角点作为待匹配点，然后计算图像在不同方向上的相位一致性特征，并将这些特征组合成一种新的局部描述子，使其能够克服多源图像间非线性辐射差异带来的匹配困难，以在多源图像匹配中获取数量更多的同名点。

1 基于方向相位特征的匹配方法 1.1 特征点检测

由于异源遥感影像之间存在特征偏移和缺失的问题，对参考影像和目标影像同时进行特征点提取很难保证特征点的重复性，因此本文仅对参考影像进行特征点提取，然后在目标影像上通过滑动窗口来搜索同名点。同时为了保证特征点分布的均匀性，采用分块的方式来控制特征角点的数量和分布，即先把参考影像划分为一定数量大小一致的格网区域，然后分别在每个区域内进行角点检测，并按响应值的大小从大到小选取一系列点作为待匹配点。

本文选择Harris^[16]算法来提取特征角点作为待匹配点，Harris算法认为一个特征角点在任意方向上滑动窗口都应该有较大的灰度变化。其算法原理如下：

$ \begin{array}{l} {E_{x, y}} = \sum\limits_{u, v} {{w_{u, v}}{{[{I_{x + u, y + v}} - {I_{u, v}}]}^2} = } \\ \sum\limits_{u, v} {{w_{u, v}}{{[xX + yY + o({x^2}, {y^2})]}^2}} \end{array} $

(1)

式中，E表示角点强度；(x, y）表示图像坐标；u, v分别是x和y方向的偏移量；X和Y分别是图像在x和y方向的一阶导数；w是高斯卷积核，用于抑制图像噪声的影响；o(x², y²)是关于x²和y²的极小值量，其计算式为：

$ {\rm Det} = \frac{{\partial I}}{{\partial x}} \otimes w + \frac{{\partial I}}{{\partial y}} \otimes w $

(2)

$ {\rm Tr} = (\frac{{\partial I}}{{\partial x}} \otimes w)(\frac{{\partial I}}{{\partial y}} \otimes w) - ((\frac{{\partial I}}{{\partial x}}\frac{{\partial I}}{{\partial y}}) \otimes w) $

(3)

$ r = {\rm Det} - c{\rm T{r^2}} $

(4)

式中，Det和Tr分别代表Harris矩阵的行列式和迹的值；$ \frac{{\partial I}}{{\partial x}}$、$ \frac{{\partial I}}{{\partial y}}$分别表示原图像I在x方向和y方向的梯度；r即为求得的图像的Harris响应值，r的大小代表了角点的强度；c为常数0.04；其余变量同前。

1.2 方向相位一致性特征

在同源遥感影像处理中，一般是通过灰度梯度信息来进行特征描述，但是对于异源遥感影像，由于传感器成像性质的不同，影像间往往存在显著的非线性辐射差异，这使得影像上同一特征梯度的大小和角度都发生了改变，基于图像梯度的匹配方法鲁棒性大大降低。为了克服这一缺陷，本文引入不同方向的相位一致性结构特征来构建一种新的特征描述，充分考虑图像在每个方向上的特征分布，实现多源遥感影像的有效匹配。

Oppenheim等^[17]证实图像的相位信息更适合用于描述图像特征，相位一致性作为一个无量纲量，它的大小与图像亮度无关。传统的相位一致性特征只能反映图像边缘的强度信息，对不同方向的边缘强度不具有选择性，本文通过选择不同角度的Log-Gabor^[18]滤波器分别对图像进行滤波，得到图像在不同方向上相位一致性能量图。其计算式如下：

$ \begin{array}{l} {\rm P}{{\rm C}_o}(x, y) = \\ \ \ \ \ \ \ \ \frac{{\sum\limits_n {{W_o}(x, y)\left\lfloor {{A_{no}}(x, y)\Delta {\Phi _{no}}(x, y) - {T_o}} \right\rfloor } }}{{\sum\limits_n {{A_{no}}(x, y) + \varepsilon } }} \end{array} $

(5)

式中，PC_o(x, y)即图像上坐标为(x, y)的点在方向o上的相位一致性能量，它的值域是0到1；n表示尺度；W_o是一个权函数；A_no是用方向为o的Log-Gabor滤波器在尺度n上得到的振幅值；T_o是噪声估计值；$\left\lfloor \right\rfloor $表示取正运算，大于0的值取自身，小于0的值令其等于0；ε是为了避免除数接近0而设置的一个极小的常数; ΔΦ_no是一种相位差分函数，其表达式为：

$ \begin{array}{l} \Delta {{\mathit \Phi} _{no}}(x, y) = \cos ({\varphi _{no}}(x, y) - \overline {{\varphi _o}} (x, y)) - \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \left| {sin({\varphi _{no}}(x, y) - \overline {{\varphi _o}} (x, y))} \right| \end{array} $

(6)

式中，φ_no是方向o和尺度n上的相位值；$ \overline {{\varphi _o}} $是方向o上所有尺度上的加权相位值；其他参数同前。

图 1是根据式（5）~（6）得到的在0、π/4、π/2和3π/4 4个方向上的相位一致性特征图像。从图 1中可以看出，不同方向的边缘在不同方向的滤波器上具有不同的响应强度。

对于一副大小为M×N的输入图像，如果构建其在O个方向上的相位特征图像，以O取4为例，可以得到一组大小为M×N的不同方向上相位特征图像（PC₁，PC₂，PC₃，PC₄），这样对于输入图像上任意一点(x, y)，该点的方向相位特征描述可以用V_p(x, y)来表示：

$ \begin{array}{l} {V_p}(x, y) = \\ \ \left[ {{\rm P}{{\rm C}_1}(x, y), {\rm P}{{\rm C}_2}(x, y), {\rm P}{{\rm C}_3}(x, y), {\rm P}{{\rm C}_4}(x, y)} \right] \end{array} $

(7)

1.3 构建特征描述向量

Fan等^[19]的研究表明，在异源图像匹配中，选择更多的区域构建的特征会更加稳定。因此，为了充分考虑特征点邻域的结构特征，本文提出了一种密集的方向相位特征描述子，通过对不同方向的特征图像进行组合, 使其能够抵抗影像间复杂的辐射差异，其构建方法如图 2所示。

在§1.2中已经得到了O个不同方向的相位特征图像，在此基础上，以待描述点P₀为中心，当给定半径R时，以P₀为中心确定一个模板窗口（长宽均为2R+1像素），在此模板窗口内构建点P₀的特征描述向量。从P₀出发以间距k在模板窗口内选取一组采样点（图 2中左侧红色区域）P₀, P₁, P₂…P_s-1，假设一共包含s个采样点，以每一个采样点为中心，累加其t×t邻域内每个像素的方向相位特征描述向量V_P，得到s个累加向量，将这s个向量组合成一个长度为s×O的向量，作为点P₀的特征描述向量。图 2中演示的是当R、k、t分别取值为4、2、3像素时的情况，其中t的取值应该保证相邻的采样点的累加邻域有一定的重叠，本文实验中k和t取该值并保持不变。

1.4 影像匹配

对于参考影像上的每一个待匹配特征点，首先按§1.3方法构建其特征描述向量，由影像的初始定位信息确定待匹配点在目标影像上的坐标，并根据影像的定位精度确定一个搜索半径，搜索半径的设置应该确保同名点位于搜索范围内且不宜过大；然后对目标影像搜索范围内的每一个点按§1.3方法构建其特征向量。搜索窗口内每一个点和待匹配点的特征向量间的归一化互相关系数（normalized cross correlation, NCC），把相关系数最大的点视作同名点，并给定一个阈值T剔除可能的错误匹配点，阈值T的选择与图像间质量差异相关; 最后利用合适的变换模型（如仿射变换模型等）和最小二乘算法剔除错误匹配点。

图 3是本文方法与互信息和归一化互相关系数在搜索窗口内的相关曲面的对比。从图 3中可以看出，本文方法在非线性辐射差异条件下有更平滑的相关曲面和显著的响应峰值。

为了使定位精度达到子像素，对相关函数进行峰值拟合，此处使用泰勒级数展开式作为拟合函数。在X₀=(x₀, y₀)^T处采用泰勒级数二阶展开得（高阶项已舍去）：

$ \begin{array}{c} f(X) = f({X_0}) + \frac{{\partial {f_T}}}{{\partial X}}(X - {X_0}) + \\ \ \ \ \ \ \ \ \ \ \ \frac{1}{2}{(X - {X_0})^T}\frac{{{\partial ^2}f}}{{\partial X}}(X - {X_0}) \end{array} $

(8)

式中，X₀就是插值中心的坐标，即匹配得到的同名点在目标影像上的行列号坐标；f是相似性测度函数；极大值点对应的X处的导数应该为0，所以令$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over X} = X - {X_0} $，然后对$ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over X} }$求导并令导数等于0可得插值中心的偏移量：

$ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over X} = - \frac{{{\partial ^2}{f^{ - 1}}}}{{\partial {X^2}}}\frac{{\partial f}}{{\partial X}} $

(9)

2 实验结果与分析

为了验证本文算法的有效性，分别选用3组不同类型的异源遥感数据进行匹配测试，具体的实验数据如表 1所示。

实验所用的影像数据是利用影像的初始有理多项式系数（rational polynomial coefficient, RPC）模型从原始遥感影像纠正得到的正射影像上裁剪得到的，通过纠正已经初步消除了影像间的旋转差异和较大的几何畸变。便于进一步的实验，且所有影像数据都已重采样到同一分辨率下，其定位精度在5个像素以内，为了验证本文算法的可靠性，本文在半径为15（像素）的矩形区域内进行搜索。

因为不同类型的遥感影像差异巨大，为了确定最佳的实验参数（方向数量O和模板半径R），分别以不同的参数对以上3组数据进行测试，绘制正确率变化曲线见图 4。

从图 4中可以发现，一般情况下，O的取值越大, 正确匹配率越高；同时，特征向量长度的增加会带来更大的计算量。因此，在保证效率和正确率的前提下，O取6进行不同模板半径下匹配效果的比较，同时为了体现本文算法的优点，选取经典的区域相关算法NCC、MI以及同样采用相位一致性的方法HOPC作为对比方法，实验结果见图 5。

从图 5中可以看出，在实验一和实验二中，由于成像波段不同带来的辐射性质上的差异，依赖于图像灰度相关性的NCC方法的匹配正确率明显降低，本文方法和HOPC算法都能保持较高的正确率；在实验三中，由于光学和SAR影像成像机制不同，影像之间不仅存在显著的非线性辐射差异，且SAR影像受斑点噪声影响产生了严重的图像退化，相关系数几乎无法得到正确匹配点。MI方法虽然在模板窗口增大时能够匹配到较多的点，但是其计算量会显著增加，且其在模板窗口较小时，准确率比较低，而本文方法相比于HOPC算法能获取到更多的匹配点。

图 6为4种方法的匹配效率对比。图 6中可以看出，本文算法在计算效率上具有明显优势，在模板窗口显著增大的同时也能保持优良的匹配性能。

图 7~9是在参数设置为O=6、R=50像素时得到的一组匹配效果图，并通过一组检查点拟合的仿射变换模型来计算匹配点的均方根误差(root mean square error，RMSE)，以此来评估算法的精度，结果见表 2。像点均方误差的计算公式如下:

$ {\rm{RMSE}} = \frac{1}{{{C_{{\rm{match}}}}}}\sum\limits_{i = 1}^{{C_{{\rm{match}}}}} {{{\left\| {H({x_1}^i,{y_1}^i) - ({x_2}^i,{y_2}^i)} \right\|}_2}} $

(10)

式中，C_match是正确匹配点的数量；H表示两张影像的真实变换模型，这里用检查点拟合的仿射变换模型替代；($x_1^i,y_1^i$)、($x_2^i,y_2^i$)是一对匹配点分别在两张影像上的坐标; $ {\left\| \right\|_2} $表示求欧氏距离。

从表 2数据分析可得，在第一组和第二组实验中，所用影像数据都是辐射质量较好的全色和红外多光谱数据，这两种类型的数据之间虽然存在一定的非线性辐射差异，但图像结构相对比较清晰、信噪比高，所以最终得到的匹配点的数量和可靠性也较高。第三组实验用的是成像机制差异较大的光学和SAR影像，且影像存在一定的时相差异，由于SAR影像受相干斑噪声干扰严重，且SAR图像的纹理特征相对于光学图像可识别度大大降低，尤其是在农田、空地等弱纹理区域，因此这两种类型的影像成功匹配的点大多位于道路和建筑分布密集区域，匹配精度也相对较低。

3 结语

本文从异源遥感影像的匹配难点出发，结合相位一致性方法在特征检测上的优良特性，用不同方向的相位特征图像组合来构建局部的特征描述向量，并通过滑动窗口和归一化互相关系数来搜索同名点。实验表明，本文提出的匹配方法在很大程度上能够克服多源遥感影像之间显著的非线性辐射差异带来的匹配困难，在实验数据上取得了比较高的精度和可靠性，尤其是能够在光学和SAR图像之间获取更多的同名点。后续需要对算法的效率和抗噪性进行改进，并结合一些粗匹配策略来缩小搜索范围，提高匹配速度。