0 概述

图像配准是将待配准图像与参考图像进行匹配。在此过程中, 参考图像保持不变, 待配准图像通过几何信息转换为参考图像对应的坐标系统^[1], 其一般分为特征提取、特征模式匹配、转换模型估计、重采样4个步骤。多模态图像配准方法可以分为基于区域配准和基于特征配准2类。

基于区域的配准方法通常需要优化全局目标函数, 对于同模态的图像, 可以最小化2个图像灰度之间差值的平方和(Sum of Squared Differences, SSD)^[2], 对于不同模态的图像, 可以使用互信息(Mutual Information, MI)^[3]和交叉累积残差熵(Cross-Cumulative Residual Entropy, CCRE)^[4], 其代表算法有梯度下降法、牛顿法和Levenberg-Marquardt法。由于MI和CCRE在每次迭代中都需要计算联合概率函数和边缘概率函数, 因此随着图像大小或模板窗口大小的增加, 算法时间复杂度也将增加。

基于特征的方法首先从参考图像和待配准图像中提取特征, 然后通过某些相似性标准进行匹配, 最后达到配准的目的。常用的图像特征包括点特征^[5]、线特征^[6]和区域特征^[7]等。局部不变特征作为计算机视觉领域的研究热点, 已被广泛用于图像配准。其中, 尺度不变特征变换(Scale-Invariant Feature Transform, SIFT)^[8]被证明是有效的特征描述符, 被应用于遥感图像配准^[9]。然而, 与光学图像不同, 合成孔径雷达(Synthetic Aperture Radar, SAR)图像反映的是目标场景的微波散射特性, 其成像机制非常复杂, 且具有斑点噪声^[6]。由于SIFT算法被设计用于检测低噪可见光图像结构特征, 因此该算法在SAR图像上的表现并不令人满意。为此, 研究人员提出了一些基于SIFT的改进方法:文献[10]提出极坐标SIFT(PSIFT)描述符, 以此来描述极光图像的特征; 文献[11]也提出了基于PSIFT描述符的图像配准方法。然而, 基于SIFT的这些特征在多模态遥感图像配准上的效果并不尽如人意。

文献[12]基于“结构和形状特征在不同模态图像中都会保留”这一假设, 采用结构特征进行多模态图像配准。方向相位一致直方图(Histogram of Orientated Phase Congruency, HOPC)可用于描述结构特征, 并已应用于遥感图像配准。但是, 使用HOPC作为Harris角点的描述符不能详细地表示其邻域角点强度的变化关系。与相位相比, 角点^[13]特征的计算复杂度相对较低, 且对图像光照、对比度以及旋转更具鲁棒性。因此, 本文提出一种改进的方向Harris角点直方图(HIOHC)结构特征来进行SAR图像和可见光图像配准。从每一副图像中独立地提取HIOHC描述符, 在此基础上, 使用夹角余弦作为相似度量(HIOHCac)匹配不同的角点, 最后通过选取的匹配点进行图像配准。

1 HIOHC描述符

本文提出基于角点信息的特征描述符HIOHC来改进Harris角点检测器。

1.1 Harris响应函数的二阶泰勒展开

角点是指在局部小窗口中观察图像特征, 并且小窗口在任何方向上移动都将导致图像平均灰度变化显著的点。设I为图像灰度, E(u, v)为由窗口移动(u, v)所引起的灰度变化, 则E(u, v)可定义如下:

$ E(u, v) = \sum\limits_{x, y} w (x, y){[I(x + u, y + v) - I(x, y)]^2} $

(1)

其中, w(x, y)是权值函数。

传统角点检测器使用一阶泰勒展开得到E(u, v), 然后阈值化响应值获得角点, 其响应值是不稳定的, 不利于构造特征描述符。为了解决这个问题, 本文在角点响应函数中引入二阶泰勒展开, 因此, E(u, v)可重新构造如下:

$ \begin{array}{l} E(u, v) = \sum\limits_{x, y} w (x, y)\left[ {{I_x}u + {I_y}v + \frac{1}{2}{I_{xx}}{u^2} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{I_{xy}}uv + \frac{1}{2}{I_{y, y}}{v^2}] = \left[ {\begin{array}{*{20}{c}} {{u^2}}&u&{uv}&{{v^2}} \end{array}} \right].\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\sum\limits_{x, y} w (x, y)\left[ {\begin{array}{*{20}{l}} A&B&C&D&E\\ B&F&G&H&I\\ C&G&J&K&L\\ D&H&K&M&N\\ E&I&L&N&O \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u^2}}\\ u\\ {uv}\\ v\\ {{v^2}} \end{array}} \right] = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{l}} {{u^2}}&u&{uv}&v&{{v^2}} \end{array}} \right]{\mathit{\pmb{P}}}\left[ {\begin{array}{*{20}{c}} {{u^2}}\\ u\\ {uv}\\ v\\ {{v^2}} \end{array}} \right] = {{\mathit{\pmb{X}}}^{\rm{T}}}{\mathit{\pmb{PY}}} \end{array} $

(2)

其中, P是5×5的对称矩阵, ${I_x} = \frac{{\partial I(x, y)}}{{\partial x}}$, ${\mathit{I}_{\rm{y}}}{\rm{ = }}$$\frac{{\partial I(x, y)}}{{\partial y}}$, ${I_{xx}} = \frac{{{\partial ^2}I(x, y)}}{{\partial {x^2}}}$, ${I_{yy}} = \frac{{{\partial ^2}I(x, y)}}{{\partial {y^2}}}$, ${I_{xy}} = $$\frac{{{\partial ^2}I(x, y)}}{{\partial x\partial y}}$, $A = \frac{1}{4}I_{xx}^2$, $B = \frac{1}{2}{I_{xx}}{I_x}$, $C = \frac{1}{2}{I_{\mathit{xx}}}{I_{\mathit{xy}}}$, $D = \frac{1}{2}{I_{\mathit{xx}}}{I_\mathit{y}}$, $E = \frac{1}{4}{I_{\mathit{xx}}}{I_{\mathit{Yy}}}$, $F = I_x^2$, $G = {I_x}{I_{xy}}$, H=I_xI_y, $I = \frac{1}{2}{I_x}{I_{yy}}$, J=I_xy², K=I_yI_xy, $L = \frac{1}{2}{I_{yy}}{I_{xy}}$, M=I_y², $N = \frac{1}{2}{I_{yy}}{I_y}$, $O = \frac{1}{4}I_{\mathit{yy}}^2$, Y=[y₁  y₂  y₃  y₄  y₅], X=VY, V是对应于矩阵P相应特征值的单位正交特征向量。

假设α₁、α₂、α₃、α₄、α₅是二次型E(u, v)的系数, 它们是P的特征值且降序排列。这意味着α₁是最大的特征值, α₂其次, 依此类推。因此, E(u, v)也可以改写为如下形式:

$ E(u, v) = \sum\limits_{i = 1}^5 {{\alpha _i}} y_i^2 $

此时, 本文考虑以下3种情况:

1) 如果图像窗口是平坦的, 则窗口沿任意方向移动仅导致很小的灰度变化, 并且所有特征值将非常小。

2) 如果图像窗口位于图像边缘上, 沿边缘方向的移动将导致较小的变化, 但垂直于边缘方向的移位将导致较大的变化。因此, α₁和α₂相对较大, 而α₃~α₅会非常小。

3) 如果图像窗口位于角点, 则沿任意方向的移动都将导致非常大的变化。因此, α₁~α₄相对较大, α₅会很小。

基于以上分析, 本文定义一个新的响应函数如下:

$ R = \sum\limits_{i = 1}^4 {{\alpha _i}} - {\alpha _5} $

传统的Harris角点响应函数仅将泰勒公式扩展到一阶项, 矩阵P仅为2×2的对称矩阵, 使角点响应图像仅在给定阈值下保持稳定(如图 1(b)所示), 将会对构造角点区域的特征描述符产生不利影响:减小阈值将产生许多噪声点, 而增加阈值将消除许多真正的角点。相反, 当使用本文方法时, 角点响应图像变得稳定, 这利于构造特征描述符(如图 1(c)所示)并且角点检测对噪声更加鲁棒。

1.2 Harris响应函数的权值

传统Harris检测器响应函数的权重是一个指示函数, 定义如下:

$ w(x, y) = \left\{ {\begin{array}{*{20}{l}} {1, \;\;\;(x, y)在窗口内}\\ {0, \;\;\;其他} \end{array}} \right. $

窗口内所有像素的贡献相同, 这不利于精确定位角点。因此, 通过结合像素的空间和灰度信息来构造不对称权重函数将是更好的选择。通过增加边缘点的权重以增加其贡献, 新的权重函数可以定义如下:

$ w(x, y) = \left\{ {\begin{array}{*{20}{l}} {\sqrt {{I_x}*{I_x} + {I_y}*{I_y}} \otimes G, (x, y)在窗口内}\\ {0, 其他} \end{array}} \right. $

(3)

其中, I_x、I_y分别表示x、y方向的导数, G是高斯核(窗口大小为3×3, 标准差为1.5), *表示Hadamard积, ⊗表示卷积。图 2和图 3展示了在真实数据上不同权值函数下的响应值。从图 2和图 3可以看出, 改进的权重函数可以检测到更精确的角点。

1.3 Harris响应函数的方向

传统Harris检测器仅考虑角点的幅值, 并未考虑角点方向, 不能有效地描述图像模板窗口的结构信息。因此, 本文提出Harris检测器的方向表示用以构造鲁棒的特征描述符。对于大小为3×3的给定模板窗口, 若中间像素是角点, 则其有8个方向, 定义方向码0到7用于表示不同的方向, 如图 4所示。

本文计算8个方向的E(u, v), 然后选择2个最小值的均值作为角点的方向, 其计算公式如下:

$ \begin{array}{l} \theta = \frac{{\vartheta (\min \{ U - \min \{ U\} \} ) + \vartheta \{ \min \{ U\} \} }}{2}\\ U = \left\{ {{E_0}, {E_1}, \cdots , {E_7}} \right\} \end{array} $

其中, i和j表示方向码, $\vartheta $(·)表示方向码的角度, U表示具有不同方向的E(u, v)构成的集合, $\begin{array}{*{20}{c}} {{\rm{min}}}\\ \cdot \end{array}${·}表示当前集合U中最小的方向码。图 5展示了图像的角点方向, 角点范围为(0°, 180°]。

1.4 结构特征描述符

本文提出的HIOHC结构特征描述符使用角点的幅度和方向来描述局部图像特征。其中, 描述符的角点方向和SIFT方法的梯度方向一样, 可以刻画异源图像的结构信息。HIOHC描述符的主要计算步骤为:首先通过计算原始图像(如图 6(a)所示)中的输入图像块(如图 6(b)所示)的角点幅度和方向(如图 6(c)所示)来创建HIOHC描述符。通过高斯函数对方向进行滤波, 然后使用幅度进行加权, 再使用L2范数标准化单元方向直方图以形成最终的HIOHC描述符(如图 6(d)所示)。描述符的计算过程如图 7所示。其中, 大方块为m=2的块, 其包括4个单元格(n=4), 最小的方块是像素。

HIOHC的详细计算方法如下:

1) 在图像中选择具有特定大小的模板窗口, 并计算模板窗口内每个像素的角点幅度和方向。

2) 将模板窗口划分为多个重叠块, 每个重叠块由m×m单元组成, 每个单元包含n像素×n像素。这些概念之间的关系如图 7所示。

3) 用高斯函数对每个像素的方向进行滤波, 然后将它们累加到方向直方图中, 并使用三线插值法对直方图进行加权。整体直方图由不同的单元直方图组成, 再由L2范数标准化, 得到最终的全局描述符。

4) 收集窗口中所有块的HIOHC描述符以形成向量, 该向量可用于特征点的匹配。

SAR图像保留可见光图像中的大部分结构特征, 如点特征、线特征、面特征等。由点可以构成线, 线可以构成面, 且角点对光照、旋转、对比度和尺度等变化不敏感。HIOHC特征能够在SAR与可见光图像中提取出这些相应的特征, 因此它可以用于表述SAR图像和可见光图像的局部结构信息。

1.5 图像配准

本文通过上述步骤获取图像的点特征及其特征描述符进行图像配准, 具体步骤如下:

1) 使用有理多项式系数^[14]模型粗略地配准参考图像和待配准图像。

2) 使用HIOHC检测参考图像中的角点, 然后均匀选择一定数量的感兴趣点作为控制点。

3) 使用HIOHC描述符检测待配准图像中以控制点为中心固定大小的搜索块中所有点的描述符, 采用夹角余弦作为相似度量, 选择其与参考图像中对应点相似度最高的点作为匹配。

4) 删除不匹配的控制点, 并使用最小二乘法计算仿射变换模型矩阵参数。不匹配的控制点是指若在待配准图像的搜索区域内找到2个及以上的匹配点, 则认为该匹配是错误。

5) 根据仿射变换模型将待配准图像中的所有像素映射到参考图像。

HIOHC描述符配准流程如图 8所示。

2 实验结果与分析 2.1 数据集

本文采用了10组遥感图像对(不同视角、不同轨道下的多种传感器所获得10幅SAR图像和10幅可见光图像)进行实验, 实验图像覆盖中国东部不同的地形地貌, 包括城市、郊区、河流和山区。实验数据集的特征如表 1所示。

2.2 配准精度

为了验证基于HIOHC配准方法的有效性, 本文采用基于HOPCncc配准方法^[12]、基于SIFT配准方法^[8]、基于SURF配准方法^[15]和手动配准方法进行对比。通过实验分析, 本文方法参数设置为:块重叠率为0.75, m=3, n=3。在手动配准中, 在参考图像和待配准图像上均匀地选择10个~30个控制点。基于SIFT配准方法首先使用SIFT算法提取参考和待配准图像中的点, 然后采用随机抽样一致算法^[16]估计仿射变换模式以实现图像配准。

本文以RMSE评价各种方法的配准精度, 结果如图 9所示。从图 9可以看出, 本文方法性能较优, HOPCncc方法其次, 手动配准方法再次。而由于基于SIFT方法和基于SURF方法失效, 其RMSE值为无穷大, 因此在图中未显示。由图 10可以观察到, 可见光图像和SAR图像具有较大的辐射差异^[17], SIFT算法不能同时有效地提取可见光图像和SAR图像中的高度重复特征(如图 10(a)所示), SURF算法也无法有效地处理SAR图像中辐射差异及噪声(如图 10(b)所示)。经计算, 本文方法的平均标准差为0.258 5, 优于基于HOPCncc方法的0.457 6和手动方法的0.444 4。因此, 本文方法比其他方法稳定。图 11为采用本文方法的部分配准结果。从图 11可以看出, 本文方法的配准结果对所有测试图像匹配准确度均较高。

3 结束语

本文提出一种新的结构特征描述符HIOHC, 用于SAR图像和可见光图像的配准。通过将响应函数从一阶项扩展到二阶项来改进Harris角点检测器, 使其对SAR图像噪声具有鲁棒性。同时综合考虑相邻像素的空间和灰度信息改进模型中的权值函数, 使得检测器更易找到角点的位置。实验结果表明, 与传统的Harris角点检测器相比, 改进后检测器的配准精度较高, 验证了HIOHC描述符的有效性和准确性。下一步将提高实时应用的计算效率, 并使用多种遥感图像进行测试。