本文方法的总体流程设计如下：（1）特征点提取。利用分块的Harris算子在输入影像上提取出分布均匀的特征点。（2）局部几何纠正。以特征点为中心取一定大小的区域影像，采用遥感影像自带的RPC文件对该区域进行局部几何纠正，以消除参考影像与输入影像间的旋转和尺度差异。（3）结构特征提取和快速匹配。利用CFOG算子进行结构特征提取，并利用傅里叶变换将其变换到频率空间，采用相位相关作为相似性测度，提高模板匹配的计算效率，实现同名点快速匹配。（4）粗差剔除和几何校正。利用RPC模型作为几何约束条件，采用最小二乘原理进行参数求解和误匹配剔除，并进行几何纠正，最终实现光学和SAR影像的自动配准，实验的总流程图如图 1所示。其中，RMSE（root mean square error）为均方根误差。

图  1  本文方法流程

Figure 1.  Flowchart of the Proposed Method in this Paper

将输入影像划分为N×N的互不重叠的格网区域，并采用具有较高稳定性和鲁棒性的Harris算子^[4]提取特征点。传统和分块Harris算子提取特征点结果如图 2所示。

图  2  两种方法的Harris提取结果

Figure 2.  Harris Extraction Results by Two Methods

对影像进行分块的目的是能够在输入影像上提取出分布均匀的特征点，因为Harris算子提取特征点需要设定特定的阈值，直接通过某一阈值在整幅影像上提取特征点会造成检测到的特征点分布不均匀，甚至会产生特征点聚簇现象，如图 2（a）所示。对影像分块后，在影像的每一小块中采用Harris算子计算每个像素的特征值并进行排序，选择前k个特征值较大的点作为特征点则可以避免提取的特征点分布不均匀的现象，如图 2（b）所示。

局部几何校正是基于有理函数模型（rational function model，RFM）对输入影像以提取的Harris特征点为中心的区域影像进行几何校正。有理函数模型是将像点坐标（r，c）表示为以其对应地面的点坐标的（P，L，H）为自变量的3次多项式比值函数^[18]：

式中，像方坐标（r，c）与物方坐标（P，L，H）均为正则化后的坐标值，其值在（-1，1）之间。式（1）为有理函数模型的正变换形式，由影像坐标到地面点坐标的转换可以得到反变换形式。

根据有理函数模型，以特征点为中心，在输入影像上选取一定大小的局部影像，引入DEM作为高程基准面，对局部影像进行几何校正，消除输入影像与参考影像之间局部影像的旋转和尺度差异，形成局部影像间的粗配准。局部几何校正前后的对比图如图 3所示，从对比图可以看出，几何校正后，影像间的旋转和尺度差异被消除了。之后根据特征点在几何校正后影像上的地理坐标预测其同名点在参考影像上所对应的位置，选取预测区域。

图  3  几何校正前后对比

Figure 3.  Comparison Before and After Geometric Correction

通过局部影像粗纠正，可以消除影像间整体的几何形变（如尺度和旋转差异），因此，光学和SAR影像配准的主要难点是影像间的非线性辐射差异。考虑到光学影像与SAR影像间具有较为相似的几何结构信息，这里采用CFOG算子^[13]提取影像间的几何结构特征，并在频率域空间建立一种快速的相似性测度进行模板匹配。CFOG是在梯度直方图（histogram of oriented gradient，HOG）的基础上发展而来的一种描述几何结构特征的描述子，它能够较好地抵抗影像间的非线性辐射差异。CFOG的构建过程见图 4，包括方向梯度通道和三维高斯卷积两个部分。

图  4  CFOG构建过程

Figure 4.  Processing Chain of CFOG

给定一幅影像，需要计算影像的m层方向梯度通道，称之为g_o，1≤o≤m，每层方向梯度通道g_o（x，y）是位于点（x，y）处方向o的梯度幅值，若值大于0，则保持不变；否则，其值等于0。可以定义每层的方向梯度通道为：

式中，I表示影像；o是所划分的方向；$⌊⌋$表示：当值为正时，等于其本身；否则，为0。实际操作不需要对每层的方向梯度通道g_o进行分别计算，可以通过使用水平方向和垂直方向上的梯度幅值$(g_x,g_y)$计算每层的方向梯度通道特征值：

式中，θ表示划分的梯度方向；绝对值的目的是将梯度方向限制在[0，π]，可以较好地处理多模态遥感影像间梯度反向的情况^[19]。

在形成方向梯度通道后，通过使用三维高斯卷积核进行卷积来实现卷积特征通道，如：

式中，$\mathrm{*}$为卷积；σ为高斯卷积核的标准差，该高斯卷积核并不是严格意义上的三维空间上的三维高斯函数，而是一个在X和Y方向上的二维高斯核和一个在梯度方向的核[1 2 1]^T（以下简称为Z方向）。通过在Z方向进行卷积平滑了Z方向上的梯度，降低了两幅影像间因局部几何形变和灰度强度变化造成方向畸变的影响；通过对$g_o^{\sigma }$进行归一化操作，形成最终的方向梯度通道特征。

CFOG是在HOG的基础上发展而来的，CFOG和HOG都可以描述影像的几何结构属性，但相比于HOG，CFOG具有以下两个优点：（1）HOG是通过对局部区域的梯度幅值和方向进行三维线性插值构建的直方图，其计算过程非常费时，而CFOG通过对方向梯度信息进行卷积运算来代替三维线性插值，有效地提高了计算效率；（2）HOG是在一个相对稀疏的采样格网上进行特征描述构建，而CFOG是一种逐像素的结构特征描述符，能更为精确地描述影像的形状和几何结构属性。

由于CFOG是逐像素的三维结构特征描述符，数据量过大，致使采用传统的相似性度量方法（如NCC）进行匹配非常耗时。根据空间域中两个模板间的相关或卷积等于两个模板的傅里叶变换在频率域中的乘积，本文采用傅里叶变换方法将特征模板从空间域变换到频率域，并以相位相关^[20]为相似性测度加速模板匹配。

假设f（x，y）和g（x，y）分别为参考影像和输入影像大小为M×M的特征模板，经过傅里叶变换分别表示为F（u，v）和G（u，v），因为已经进行过局部几何纠正，所以特征模板间只存在平移关系（x₀，y₀），即g（x，y）=f（x-x₀，y-y₀），根据傅里叶变换的平移特性，在频率域下，两者的关系为：

则它们的归一化互功率谱可以表示为：

式中，*表示复共轭。根据平移理论，归一化互功率谱的相位等于相位间的相位差，通过对归一化互功率谱进行傅里叶反变换后可以得到一个二维脉冲函数$\delta (x-x_{\mathrm{0}},y-y_{\mathrm{0}})$。此函数在偏移位置处有明显的尖锐峰值，其他位置的值接近于0，通过峰值的位置可以确定特征模板之间的偏移量x₀和y₀，实现特征模板的快速匹配，如当大小为N×N的模板窗口在大小为M×M的搜索窗口上进行逐像素滑动匹配时，传统相似性度量（如NCC）的运算量为$\zeta (M^{\mathrm{2}}{N}^{\mathrm{2}})$，进行傅里叶变换到频率域后，采用相位相关方法的运算量为$\zeta ({(M+N)}^{\mathrm{2}}\mathrm{l}\mathrm{g}(M+N))$，根据上面运算量的对比，采用相位相关方法的计算效率有显著提升。

这里使用相似性图对匹配模板在搜索区域的相似性系数进行可视化表达，并将本文方法与MI和NCC比较，如图 5所示，蓝色表示相似性较低，红色表示相似性较高。

图  5  NCC、MI和本文方法的相似性图对比

Figure 5.  Comparison of Similarity Map of NCC, MI and Our Proposed Method

从图 5（c）~图 5（e）的3种相似性图可以看出，NCC方法的相似性图有多个峰值，表示未能匹配到正确的同名点；而MI虽匹配到一个峰值较高的同名点，但峰值的位置具有较大的偏离，表示匹配到的同名点的误差较大；本文方法的相似性图只具有一个较高的峰值，且峰值的位置位于中心处，说明了本文方法的匹配性能优于MI和NCC的匹配性能。

由于遮挡或阴影等影响，通过相位相关法得到的匹配结果不可避免地存在着错误同名点对，因此，需要对其匹配结果进行筛选和优化。本文依据最小二乘原理，通过影像之间的几何约束关系补偿有理函数模型的系统误差，得到正确匹配，提高配准精度。基于影像之间的几何约束关系补偿有理函数模型的系统误差可描述为：

式中，（x，y）为同名点在输入影像上的图像坐标；$(a_{\mathrm{0}},a_{\mathrm{1}},a_{\mathrm{2}})$和$(b_{\mathrm{0}},b_{\mathrm{1}},b_{\mathrm{2}})$为补偿参数；（P，L，H）为对应同名点在参考影像上的物方坐标。依据最小二乘原理解算补偿参数后，计算同名点的残差和RMSE，剔除残差较大的同名点。重新解算补偿参数，并对上述过程进行迭代，直至RMSE小于设定的某一阈值（如1个像素）。基于有理函数模型和像方误差补偿参数对输入影像进行几何纠正，实现输入影像与参考影像的自动配准。

本文选择了3组国产高分辨率卫星遥感影像作为实验数据验证本文方法的配准性能，如图 6所示。

图  6  实验数据和匹配点分布图

Figure 6.  Experimental Data and Matching Points

图 6中，黄色圆点表示本文匹配的正确同名点。参考影像为做过正射纠正的高分三号SAR影像，输入影像分别使用高分二号多光谱影像和资源三号全色影像，实验数据的详细信息如表 1所示。

基于实验数据，采用人工刺点的方式在参考影像和配准后的影像上选取50对同名点作为检查点，通过计算它们的RMSE来表示配准精度。本文将分别从匹配正确率、配准精度和计算效率等方面对本文方法进行分析与评价，并将传统的模板匹配方法MI、NCC，以及同样基于结构特征的HOPC描述符，结合上述的局部纠正，与本文方法进行对比。实验所用PC配置为：Inter（R）Xeon（R）CPU E5-1603 v3 2.8 GHz，内存为8 ‍GB，采用C++和GDAL源码库编程实现。本文方法配准的参数设置如下：将输入影像划分为25×20的块区域，在每块区域中提取一个Harris特征点，共计500个特征点，为了保证不同方法性能对比的客观性，所有匹配方法的模板窗口均设置为121×121像素，搜索区域窗口的大小设置为200×200像素。

实验1和2位于地形平坦、地物信息较为丰富的区域；实验3位于四川省眉山市，地势西高东低，起伏较大，且分辨率相差1.25倍，配准难度较大。配准过程中采用不同方法匹配的正确匹配点数量、匹配正确率、RMSE和计算速率见表 2。

由表 2可以看出，对于3组实验，NCC方法的匹配正确率都处于最低，这是因为NCC方法对影像间的灰度只具有线性不变性，当影像间的辐射差异较大，尤其是面对非线性的辐射差异时，NCC方法通常不能取得满意的效果^[21]。MI方法在匹配点获得的数量方面虽具有可观性，但经过误差剔除后，正确匹配点的数量却较少，3组实验中的匹配正确率最高仅有44.56%，对于大区域影像的配准而言较低，难以满足大区域影像配准的要求。HOPC特征在正确匹配点对和正确匹配率上取得了良好的效果，3组实验的正确匹配点数量分别为441对、413对和395对，匹配正确率也分别达到了92.26%、88.45%和83.68%，但HOPC的匹配耗时较长，表明匹配速率较低。

相较于MI、NCC及HOPC，本文方法在正确匹配点数量和匹配正确率上表现更佳，正确匹配点分布如图 6所示。由图 6可知，实验1和实验2的正确匹配点均超过了400对，且匹配正确率也均超过了90%，究其原因是本文使用遥感影像的几何结构特征构建匹配模板，相比于灰度信息，影像的几何结构特征表现更加稳健和鲁棒。相比于实验1和实验2，实验3在正确匹配点数量上和匹配正确率都有所降低，原因是实验3的数据地形起伏较大，局部几何畸变显著，增加了匹配的难度。另外，基于有理函数模型的像方误差补偿参数也难以对地形起伏较大的区域进行精确拟合，致使剔除了部分正确的匹配点，但本文方法的匹配正确率依然远高于MI、NCC及HOPC。

从匹配耗时可以看出，MI极为耗时，这是因为MI需要计算模板间的联合概率直方图，计算量较大，NCC的计算效率虽快于MI，但匹配仍然耗时。HOPC方法需要提取影像间的结构特征，再计算结构特征间的相关系数，比NCC方法更为耗时。本文通过傅里叶变换技术将匹配模板变换到频率域，加快了模板匹配的速度，因此，耗时较少。3组实验的配准结果如图 7所示。

图  7  实验配准棋盘格图

Figure 7.  Registration Checkboard Results of Experiment

从配准精度来看，图 7（a）实验1数据是建筑物较为丰富的高分辨率影像，且时相差异较小，因此，配准精度达到了子像素的精度。图 7（b）实验2影像中包含了大量的水系，由于CFOG是根据影像的几何结构特征构建特征描述符进行匹配，应用于水系，CFOG难以在水系中提取出明显的几何结构特征，致使匹配点数量和匹配正确率有所下降。图 7（c）实验3数据虽受到地形起伏的影响，本文方法的配准精度也达到了1.052像素。从实验的局部配准结果可以看出，虽然光学影像与SAR影像间存在显著的噪声和非线性辐射差异，但配准后的输入影像与参考影像依然叠合得非常好。对于3组实验数据，本文方法的配准精度都好于基于MI、NCC和HOPC的方法。

针对SAR影像与光学影像间几何形变和非线性辐射差异所造成的配准困难问题，本文引入能够描述影像结构特征的描述符——CFOG，并根据遥感影像间的几何约束关系，结合相位相关技术实现了大区域SAR影像与光学影像的快速自动配准。为了验证本文方法的性能和表现，选择了3组具有时相差异和尺度差异的国产高分辨率卫星遥感影像进行配准实验，并将其与结合局部校正的MI、NCC以及HOPC等方法在多个方面进行对比与分析。结果表明，本文方法在3组实验中的正确匹配点数量和匹配正确率都远高于传统模板匹配方法，配准精度也表现最佳，且本文方法借助于傅里叶变换技术实现快速模板匹配，大幅度提高了计算效率，表明了本文方法可以有效抵抗光学影像与SAR影像间的几何形变和辐射差异，从而获得较高的匹配正确率，实现大区域光学影像与SAR影像间的配准。

由于本文方法是通过遥感影像的姿态参数信息RPC对影像进行局部校正，消除影像间的旋转和尺度差异，因此，本文的实验前提是可以利用卫星影像的姿态参数信息对影像进行局部校正。另外，本文方法是利用影像的几何结构特征构建三维特征描述符，当影像的几何结构特征不够丰富时，其匹配性能会有所下降，并且构建的三维特征描述符也存在数据量过大的问题。后续将对这些问题进行进一步研究，以提高光学影像与SAR影像间配准的精度与速度。