一种基于相位一致性相关的多源遥感影像配准方法

范登科，潘 励，叶沅鑫

(武汉大学遥感信息工程学院，武汉 430079)

1 概述

图像配准作为一项基础性的影像几何处理技术，在医学、计算机视觉等学科领域中发挥着重要作用。尽管它在遥感中的研究和应用起步较晚，但是由于遥感影像在采集环境、观测角、地形起伏等多方面因素作用下所体现出的复杂性远远大于其他图像，再加上其本身海量及多源异构的特性，对配准技术在该领域的研究和应用提出新的挑战。目前，国内外众多学者针对遥感影像匹配问题的不同方面，提出了很多具有可行性和研究价值的配准策略和方法，也因此产生了多种将配准方法进行分类的方式，包括基于配准影像数据源类型的分类、基于影像间几何畸变类型的分类等。最常见的是基于匹配时所使用影像信息类型的分类方式［1]，该方式将配准分为特征匹配和区域匹配两类。

特征匹配是当前配准技术研究的热点方向，它通过计算影像的梯度、矩等强度或对象的几何形态信息来描述局部范围内的显著性特征，然后选择适用于该特征匹配的相似性测度建立同名点的对应关系，最终完成对几何畸变的改正。常见的用于遥感影像配准的点特征有:Harris、Forsterner、Morovac，线形或边缘特征有:Canny、Sobel，面状或区域特征有Hu不变矩、形状上下文。SIFT［2]是近几年图像配准研究领域兴起的一种具有尺度和旋转不变性描述的特征算子，它以其多尺度、稳定的点特征表达特性赢得许多研究学者的青睐，并由此衍生出许多基于SIFT的配准方法，如PCA-SIFT［3]、CSIFT［4]、SR-SIFT［5]等。然而，由于大多数特征算子都是以图像梯度的强度和方向为基础的，且很少顾及区域性信息，造成多源异构遥感影像配准的同名点匹配率和纠正精度不能满足实际应用生产的需要，因此特征匹配方法具有应用局限性和问题针对性。

区域匹配作为一种早期发展成熟的、方便易用的图像配准方法，被广泛应用于各类商业化图像处理软件中。它通过统计运算的方式对特定区域内的图像信息进行描述，并借助相关关系模型度量彼此间的相似程度，相对于特征匹配需要稳定且充足的特征作为配准前提条件来说，区域匹配具有很强的通用性。目前研究和应用最多的区域匹配方法有相关系数和互信息，它们分别是基于影像原始灰度和信息熵进行计算和匹配的。Gao使用相关系数作为匹配策略［6]，借助影像金字塔和二次多项式实现对Landsat GeoCover数据子像素精度的自动配准;Peter［7]同样以相关系数为核心匹配手段，完成对HyMap、AIRSAR等24组同源和非同源影像对的几何校正;Sahil［8]通过迭代计算TerraSAR-X和Ikonos影像间的互信息，研究分析了城镇区域异构图像配准的精度。由于非同源异构影像间在原始灰度上通常是非线性变化的，且信息熵的计算量比较大，导致区域匹配方法的推广应用受到一定的制约。

本文综合特征匹配和区域匹配的优势，提出一种适用于非同源异构遥感影像的配准方法，该方法以不同滤波尺度下位置稳定不变的角点作为控制点，借助金字塔分层映射搜索策略，通过计算领域和搜索区域内的相位一致性，引入相关系数作为相似性测度实现同名点的匹配。并最终给出利用可见光、近红外及SAR影像对该方法进行配准实验的验证和分析。

2 配准方法

2.1 控制点提取

为了保证提取足够多且位置相对稳定的控制点进行匹配和模型参数估计，需要使用通用性强、可靠性高的角点特征提取算子在待配准的输入影像上执行计算和检测。当前，Harris［9]和 FAST［10]被证明具有较高的角点提取重复率，适用于各类光谱组成结构的影像。本文选择Harris算子在输入影像上提取特征控制点，Harris通过局部图像信息的自相关方式，在各个方向都计算了灰度变化，其角点提取的原理是:首先使用标准偏离差为σ的高斯滤波模板与原始灰度影像I进行卷积运算，得到平滑后的影像¯I，然后计算¯I的一阶梯度矩阵M，公式为

其中k为介于0.04～0.06的常数，如果2个特征值都比较大，说明在该点处沿任意方向运动会引起判别函数R的急剧变化，判断为具有强响应的角点，否则为弱边缘或平滑域点，在使用阈值T对该值进行局部非最大抑制后，可以获得稳定的特征点。

由于使用不同标准偏差σ的高斯卷积模板会使平滑处理后的影像¯I变模糊，在特定滤波器尺度下检测的角点可能会在其他尺度的滤波影像中发生位置偏移甚至消失。为了提取不受高斯滤波尺度影响的稳定角点，需要使用临近尺度的高斯模板重新提取角点进行检测和判断。假设角点(x0，y0)是使用尺度(标准偏差值)为σ0的高斯滤波模板所提取的Harris特征点，构造尺度为 σ1=σ0-δ×σ0和 σ2=σ0+δ×σ0的高斯滤波模板对原始影像I进行卷积运算，提取的角点分别为(x1，y1)和(x2，y2)，当提取的角点满足条件

并且

时，认为(x0，y0)是一个不受卷积尺度影响的稳定角点，选为控制点，否则被排除。公式(3)、(4)中D代表欧式距离，即以角点(x1，y1)和(x2，y2)相对于(x0，y0)的位置偏移作为特征点的判定依据。图1展示了分别使用原始Harris算子和多尺度高斯滤波后Harris提取的特征点分布情况，可以看出，在加入了公式(3)、(4)的判别后，位于道路边缘和开阔地中心的点受平滑尺度变化的影响，位置变化显著，被判断为不稳定的特征点排除，保留的点为稳定的不受平滑滤波干扰的特征点。

2.2 同名点匹配

图1 Harris算子提取特征点分布示意

由于Harris特征算子提取的控制点为梯度变化显著的角点，因此在其领域内影像的纹理信息比较丰富，灰度变化显著，相关系数在小范围内容易形成变化剧烈的极值，可以选择作为同名点的匹配测度。然而，非同源异构遥感影像间的光谱组成存在明显差异，在灰度强度上呈现非线性变化，导致相关系数不能稳定描述相似度高的同名区域。针对这一问题，本文引入相位一致性对影像进行转换，得到相关性程度高的特征图像。

相位一致性是一种使用傅里叶谐波分量描述信号局部强度的特征，它利用特定方向θ不同尺度s的Log Gabor滤波器，将信号分解为频率域下与θ和s相对应的傅里叶谐波分量，通过加权叠合不同尺度的分量，获得具有原始信号强弱特性描述的响应。Morrone证明了相位一致性在描述图像特征信息上与人类神经视觉具有相一致的敏感度［11]。Wong指出，相位一致性特征具有抵抗图像对比度和亮度变化的能力，并利用其构造的最小和最大距特征，分别实现了遥感影像特征点的提取和同名点的匹配［12]。应用于图像二维离散信号，相位一致性的计算公式为

式中，PC(x，y)为点(x，y)处的相位一致性，θ为滤波方向，在区间［0，π]之间以固定步长变化，n为滤波器尺度s的最大值，Wθ(x，y)为频率扩展的权重系数，As，θ(x，y)代表使用沿 θ方向尺度为 s的 Log Gabor滤波器进行分解后谐波分量的振幅值，」表示仅取正值的运算，负值全部赋为0，T为噪声阈值，ε是一个避免除零的常数，Δφs，θ(x，y)可由分解后谐波分量的相位角计算获得，公式为

图2 Landsat 5 TM影像及其相位一致性特征图像

从图2可以看出，边缘和角点处的相位一致性明显大于平坦开阔地区，通常在角点处取得极值。同时，相位一致性特征影像反映的并不是二值化的边缘，而代表的是不同尺度和方向频率域卷积运算下傅里叶谐波分量的累积值。因此，以特征点为中心的窗口发生位置改变时，不会引起相关系数的剧烈变化，且由于角点处相位一致性的显著特征，使得极值比较容易获得。经计算对应图像的相关系数，图2中(a)与(c)之间的相关系数为-0.045，而(b)与(d)达到0.605，可见采用相位一致性特征度量影像间的相关性明显优于原始灰度。

遥感影像在采集过程中一般表现出大视角、宽视场的特性，成像的图幅范围大，使全局匹配的效率降低，计算复杂度增大，需要使用特定的搜索策略优化同名点匹配过程。基于金字塔结构的同名区域预测与分层索引映射是目前最受欢迎的、应用于遥感影像匹配中的搜索策略，方法是以一定比例逐级降低图像空间分辨率，缩小图幅范围。首先，在输入影像的顶层金字塔图像上提取特征点，并搜索参考影像对应尺度的金字塔，找出相关系数大的窗口中心作为匹配点，使用简单的二维平面转换关系模拟初始的几何变形，根据最小二乘准则删除残差大的误匹配点对;而后，利用不同尺度金字塔图像的映射关系，在下一层更高分辨率尺度的金字塔图像上获取同名区域，执行与上述一致的精确匹配，直到达到底层即原始分辨率影像为止。两层金字塔结构的同名区域映射关系及几何变换模型估计如图3所示。

2.3 模型参数估计与重采样

利用匹配获取的同名点对，依据最小二乘原理估计几何畸变模型的参数，剔除残差大的同名点对以保证模型估计精度在可接受的范围内，最终使用特定的插值方法实现对输入影像的重采样处理，输出在地理位置上与参照影像配准的校正影像。在这一过程中，全局性几何畸变模型的选择显得十分重要。

图3 基于金字塔结构的遥感影像同名区域索引和搜索策略

由于传感器姿态和侧视角的差异，采集获得的多源异构遥感影像在全局几何畸变上表现出不同形式。中分辨率光学影像的成像方式为推扫式，其传感器所处轨道高度及线阵扫描方式决定了影像几乎不受侧视角和地形起伏的影响，在全局影像上仅存在平移、旋转和尺度畸变，因此可以使用仿射模型描述待配准光学影像间的几何变形关系。仿射变换模型公式为

式中，(x，y)和(X，Y)分别为对应同名点在输入影像和参考影像的图像坐标，λ代表尺度缩放因子，θ表示影像的相对旋转角，(c，r)是图像在二维平面的相对平移量。因此仅需要两对同名点即可以求解该模型参数。而同分辨率SAR影像的采集使用主动传感方式，虽然轨道高度避免了地形起伏的影响，但是其成像图幅范围较小，并且存在一定的侧视角，因此对于光学影像与SAR影像以及不同时相SAR影像间的几何畸变，通常采用投影变形模型来描述，一次项投影模型公式为

式中包括a～h总计8个模型参数，需要4对同名点进行求解，分母部分描述了成像平面的倾斜变形。在最小二乘原理基础上对模型参数进行估计后，使用参考影像中的控制点坐标反算输入影像控制点的模型坐标，并与真实坐标求差后计算均方根误差RMSE，RMSE值的大小代表了该同名点对偏离模型估计的程度，删除RMSE值最大的多个点对，重新估计模型参数，迭代执行这一过程直到同名点对集合中最大的RMSE值小于预先设定的阈值，以获得满足纠正精度要求的模型参数。最终，双线性插值方法被用于重采样过程输出配准结果。

2.4 总体配准策略

本文提出的基于相位一致性相关的多源遥感影像配准方法可概括为以下7个步骤:(1)根据影像大小，设定特定的分辨率比例，分别构建输入和参考影像的多级金字塔图像;(2)在输入影像的顶层金字塔上使用多尺度Harris提取控制点，并计算其领域内窗口的相位一致性;(3)计算参照影像顶层金字塔的相位一致性特征图像，使用相关系数作为匹配测度全局搜索最大相关的窗口，以其中心作为匹配同名点;(4)用简单二维平面变换模型估计畸变改正关系，预测下一级金字塔图像的同名区域;(5)在下级金字塔图像的子区域内迭代执行步骤(2)～(4)，直到完成原始分辨率影像上同名点的匹配;(6)根据待配准影像的类型确定全局几何变换模型，采用最小二乘原理迭代解算和修正模型精度;(7)重采样输出配准后的结果影像。图4展示了配准策略的完整流程。

图4 配准方法流程

3 实验与分析

3.1 实验背景和数据

实验以Windows XP系统为平台，在Visual C++2008集成开发环境中设计开发程序实现上述配准方法，程序运行的硬件环境为Inter Core2 CPU 2.66 GHz，2 GB内存。实验数据选用3组遥感影像，每一组中输入影像和参考影像分别采集自异构光谱范围或不同传感器。前2组数据均由光学遥感影像组成，采集的区域为湖北省大东湖和梁子湖水域，其中第1组数据来自于同一传感器，为了体现空间位置上的差异性，手工对输入影像旋转10°，第2组数据来自不同传感器，存在微小旋转及10～20像素的平移变形。第3组数据为光学影像与SAR影像之间的配准，为了正确计算相关系数，对SAR影像进行重采样，确保其与光学影像的分辨率一致。3组实验数据的详细信息如表1所示。

在进行影像配准运算处理过程中，需要对相关参数进行初始化，具体项目及其初始值为:金字塔各级比例1∶3，金字塔级数3，原始分辨率下搜索区域大小200像素，相关系数阈值0.75，相关系数计算窗口大小13像素，RMSE阈值1.5。

3.2 配准精度分析

使用以上3组数据对本文配准方法进行实验验证，并对该方法的适用性和精度进行评价。表2列出了经实验获得的配准评价结果。

由表2可以看出，在配准精度方面，第3组数据的效果最好，最大残差和总体均方根误差指标都控制在1个像素以内，然而同名点对的误匹配率却达到了95%以上，仅有5对正确匹配的同名点用于模型参数估计，且它们分布集中于图5(e)、(f)中左下角的岛屿海岸线上，占影像中主体的陆地范围内并没有正确匹配的同名点对，这是因为大陆中地物类型丰富而海岛内单一，SAR影像之于前者表现出更多的噪声，而相位一致性的计算结果对噪声十分敏感，造成陆地中相位一致性的相关程度并不是十分显著;反之，海岛在SAR影像中纹理结构简单，海岸线轮廓特征明显，有利于相位一致性稳定描述局部图像信息。

表1 实验数据信息概览

表2 配准方法的实验评价结果

图5 3组实验数据影像及同名点分布情况

前2组实验数据的配准结果在精度和误匹配率上达到一致水平，尽管第2组数据的总体均方根中误差略微偏高，但其匹配得到的同名点对数目是第1组的3倍，是因为第2组数据的影像对均由近红外波段构成，相对于第1组的可见光-近红外组合，相关程度更高。然而，在匹配获得大量同名点对的同时，发生误匹配的点对数目也相应提升，且从2组数据对比情况来看，这一比例水平并不受光谱组成结构的影响，因此要进一步降低误匹配率，仍需要对匹配方法和策略进行改进。另一方面，由于第2组影像采集自不同时间，地物变化也会造成同名点的误匹配。从图5(a)～(d)中可以看出，正确匹配的同名点均匀分布于整个图幅范围内，且在对比度高的局部区域内相对集中，符合相位一致性描述亮度和对比度差异的不变特征。

综合分析3组实验数据的配准结果，本文方法更适用于对非同源异构的两景光学影像进行配准，对于光学与SAR影像的配准有一定局限性和区域选择性，同时将相关系数作为匹配测度引入到该方法中，对于全局仿射和投影变形有一定的抵抗能力。

3.3 对比验证

为了进一步验证本文方法的配准效果，引入传统基于灰度的相关系数(NCC)、尺度不变特征(SIFT)、改进尺度不变特征(SR-SIFT)3种方法分别对3组实验数据进行配准实验，统计各项实验结果的误匹配率和均方根中误差(RMSE)，得到的对比情况如表3所示。

表3 3种常用方法实验结果对比(“—”表示配准失败)

实验过程中设定RMSE阈值为1.0像素，以保证几何校正模型参数估计的准确性和一致性。当3种方法达到相同精度水平时，结果显示实验2的误匹配率最低，原因是两景近红外影像间灰度相关性程度最高，且与特征表述相关的梯度强度和方向变化较为稳定，与此同时，NCC和SR-SIFT方法的误匹配率明显低于本文方法。相反，3种方法对实验1数据进行配准的误匹配率高于或接近本文方法，可见，对于光谱组成差异越显著的影像来说，本文方法的适用性更强，效果更好。3种方法不能完成对实验3数据配准的事实更进一步验证了这一结论。因此，与3种配准方法进行对比分析得出的结论是，基于相位一致性相关的配准方法可以应对传统方法所不能解决的异构遥感影像间难以配准的问题。

4 结语

本文提出一种基于相位一致性相关的多源遥感影像配准方法，该方法以多级图像金字塔结构为索引和搜索策略，首先使用多尺度高斯滤波判别和遴选稳定的Harris特征点，在控制点领域和相关系数极值搜索区域内，引入相位一致性特征进行相似性度量，在降低误匹配概率的同时获得更多精确匹配的同名点对。经实验分析和对比验证，该方法适用于中分辨率多源异构光学-光学、光学-SAR遥感影像间的配准，精度满足实际生产过程中影像几何校正处理的需要，且高于其他基于区域匹配的配准方法。下一步的研究工作将专注于提高运算效率，改进算法以降低误匹配率，并适用于不同分辨率、地形起伏等复杂条件下的多源遥感影像配准。同时希望本文方法能为图像配准领域的学者提供研究基础和借鉴。

［1]Zitová B，Flusser J.Image registration methods:A survey［J].Image and Vision Computing，2003，21(11):977-1000.

［2]Lowe D G.Distinctive image features from scale-invariant keypoints［J].International Journal of Computer Vision，2004，60(2):91-110.

［3]Ke Y，Sukthankar R.PCA-SIFT:A more distinctive representation for local image descriptors［C].Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2004:506-513.

［4]Abdel-Hakim A E，Farag A A.CSIFT:A SIFT descriptor with color invariant characteristics［C].Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition，2006:1978-1983.

［5]Yi Z，Zhiguo C，Yang X.Multi-spectral remote image registration based on SIFT［J].Electronics Letters，2008，44(2):107-108.

［6]Gao F，Masek J，Wolfe R E.Automated registration and orthorectification package for Landsat and Landsat-like data processing［J].Journal of Applied Remote Sensing，2009，3(1).

［7]Bunting P，Labrosse F，Lucas R.A multi-resolution area-based technique for automatic multi-modal image registration［J].Image and Vision Computing，2010，28(8):1203-1219.

［8]Suri S，Reinartz P.Mutual-information-based registration of Terra-SAR-X and Ikonos imagery in urban area［J].IEEE Transactions on Geoscience and Remote Sensing，2010，48(2):939-949.

［9]Harris C，Stephens M K.A combined corner and edge detector［C].Proc.Alvey Vision Conf.，1988:147-152.

［10]Rosten E，Drummond T.Machine learning for high-speed corner detection［C].Proceedings of the European Conference on Computer Vision，2006:430-443.

［11]Morrone M C，Owens R A.Feature detection from local energy［J].Pattern Recognition Letters，1987，6(5):303-313.

［12]Wong A，Clausi D A.AISIR:Automated inter-sensor/inter-band satellite image registration using robust complex wavelet feature representations［J].Pattern Recognition Letters，2010，31(10):1160-1167.