CN104361590B - 一种控制点自适应分布的高分辨率遥感影像配准方法 - Google Patents

Abstract

本发明公开一种控制点自适应分布的高分辨率遥感影像配准方法，将多尺度J‑image引入到图像配准中，首先提出了一种基于分块策略的自适应控制点提取方法，实现了多尺度J‑image影像中的控制点提取，从而克服了控制点仅对个别方向高频信息敏感的局限。同时，定义了自适应控制点提取策略以约束控制点的分布。进而采用归一化的互信息测度NMI对控制点进行多尺度进行匹配，能够有效平滑配准函数。最后采用Delaunay三角形局部变换实现几何校正。下文首先介绍了所提出算法的基本原理及关键步骤，进而对三组不同类型的遥感影像进行了实验和分析，并分别与基于小波变换及NSCT的不同配准方法进行了比较实验。

Description

一种控制点自适应分布的高分辨率遥感影像配准方法

技术领域

本发明涉及一种控制点自适应分布的高分辨率遥感影像配准方法属于，遥感影像配准技术领域。

背景技术

1图像配准的基本框架

随着高分辨率遥感影像的广泛应用，针对高分辨率遥感影像的图像配准研究成为遥感应用研究领域的一个重要课题。图像配准的基本框架主要包含以下几个方面的内容：

(1)特征空间，是指在参考影像与待配准影像中提取的图像特征，主要包括灰度特征、点特征、边缘特征以及不变矩结构特征等。

(2)搜索空间，是指图像配准在几何较正时变换的范围及采用的方式。其中，变换范围包括全局、局部以及位移场三种类型。变换方式分为线性变和非线性变换，线性变换又可分为刚体变换、投影变换以及仿射变换三种类型。

(3)相似性度量，用于对变换结果的优劣进行评价的准则。常见的相似性度量有互信息、归一化互信息、欧氏距离、联合熵、梯度互相关性等。通常认为特征提取与相似性度量直接决定了配准算法的鲁棒性。

(4)搜索策略，用于在搜索空间中寻找最优的配准参数，并以相似性度量作为评判的依据。常见的搜索策略包括蚁群算法、遗传算法、抛物线法、Powell法、牛顿法等。

2图像配准的分类

目前图像配准方法还没有统一的分类规则，通常情况下主要可分为基于灰度的图像配准以及基于特征的图像配准。

基于灰度的配准方法直接利用影像的灰度值建立参考影像与待配准影像的空间变换关系，其基本思想为：首先对待配准影像进行几何变换，而后建立参考影像与待配准影像之间的相似性度量，即在统计影像灰度信息的基础上定义一个目标函数。将该目标函数作为评价配准结果优劣的判别准则以及几何配准参数最优化的目标函数，使得影像配准问题转化为求该目标函数的极值问题。最后，采用一系列最优化方法确定最优配准参数。

尽管基于灰度的配准方法直观而且易于实现，但计算相似性度量时通常以模板所在位置为中心，因而形成的峰面较为平缓，难以获得精确的配准结果。另一方面，由于存在光照条件变化、传感器类型不同等因素，多时相高分辨率遥感影像配准中常存在较大的灰度、旋转等差异，基于灰度的配准方法对这些干扰因素的鲁棒性较差。而基于特征的配准方法能够有效克服这些问题，因此目前高分辨率遥感影像配准主要采用基于特征的配准方法。

基于特征的配准方法利用图像的特征作为配准的依据，这些特征应具有对灰度变换及空间变换不敏感的特点，从而使方法具有更好的可重复性与稳定性，常见的特征包括角点、边缘特征等。基于特征的配准方法基本步骤为：

(1)对图像进行预处理，以降低或消除参考影像与待配准影像之间的几何形变以及灰度偏差，减少对配准过程的影响。

(2)选择或定义合适的特征，还可以采用SIFT等特征提取算子实现特征自动提取。若采用角点特征(控制点)进行匹配，还应使控制点在影像中标识出明显的地形地物，并保证一定的数量及合理的分布。

(3)选择或定义合适的匹配策略，将满足条件的控制点对作为一对同名点，获得所需的同名点集。

(4)根据同名点集建立参考影像与待配准影像间的映射关系，并根据坐标及参考影像对待配准影像进行重采样，获得最终的配准结果。

在多种基于特征的配准方法中，基于点特征的配准方法应用及研究最为广泛。其基本思想为首先提取图像中的控制点，进而基于控制点的局部灰度、梯度等特征以及控制点之间的几何关系获得同名点对的集合，从而实现图像间的完全配准。其优点在于只对控制点进行匹配，从而与计算所有像素相比能够显著减少计算量；控制点在提取过程中能够有效减少灰度差异、几何形变、噪声以及其他干扰因素的影响，因而具有很好的鲁棒性；基于控制点的相似性度量对位置变化敏感，能有效提高配准精度；同名点集可以直接用于计算参考影像与待配准影像间的空间变换关系。基于这些优点，基于点特征的影像配准方法通常能够实现快速准确的高精度图像配准，同时具有广泛的适用性，是目前主流的图像配准方法，同时也代表了图像配准技术的发展方向。

3存在问题

针对几何配准误差对变化检测精度的影响，目前已有很多学者展开了深入的研究。如Dai等人的研究表明，影像间错误的位置匹配会导致“过变化检测”。而Verbyla等人则认为，由于影像纠正模型中不可避免的存在残差，因此几乎没有完全精确的影像匹配，这就导致了检测结果中的“伪变化”是几何配准误差在变化检测中的必然产物。在Townshend等人对多时相NDVI(Normalized Difference Vegetation Index)影像进行差值比较的变化检测中，至少40％是由于亚像素级别的配准误差造成的。而Dai等人的研究成果表明，在中、低分辨率遥感影像变化检测中，只有配准精度达到五分之一像素才能够获得90％的变化检测精度。而高分辨率遥感影像变化检测更容易受到配准误差的影响，因此无论采用哪种具体变化检测方法，精确的几何配准都是成功的前提和保证。

随着遥感影像空间分辨率的不断提高，遥感影像提供了更加丰富的特征信息用以表达地物。另一方面，与中低分辨率遥感影像相比，高分辨率遥感影像在不同时间尺度(时间尺度也叫做时间分辨率，是指在同一位置的连续采集图像之间的时间间隔)中地物特征随着季节变化所产生的状态差异更加突出。同时，由于传感器内部成像***以及地形起伏、云层遮挡、卫星轨道变化、飞行器姿态不稳定等造成的非线性畸变等干扰因素对配准造成的影响也更加明显。针对这些问题，基于点特征的配准方法能够提供一条有效的解决途径。通过基于点特征配准方法基本思想及实现流程可以看出，特征空间中的控制点合理选择及控制点间的精确匹配是提高配准精度的关键因素。

根据高分辨率遥感影像的特点，采用多尺度下的控制点提取与匹配能够获得更加可靠的同名点集，因而越来越多的受到重视。目前常见的多尺度分析工具主要包括小波变换(Wavelet Transform)、轮廓波变换(Contourlet Transform)[106]、非下采样轮廓波变换(NSCT，Non-Subsampled Contourlet Transform)[107]等，例如I.Zavorin等人提出了一种基于小波金字塔的高分辨率遥感影像配准算法，能够在多源传感器影像配准中取得良好的效果。Z.G.Chen等选择NSCT进行多尺度控制点提取及同名点匹配，显著提高了配准精度。尽管多尺度分析工具有利于更加准确的提取图像的特征，但基于小波变换的方法通常会忽略45°方向子带信息，且不具有平移不变性；而基于轮廓波变换以及NSCT的方法同样存在对方向选择敏感的问题，即选择不同方向的高频子带会直接影响图像的配准精度，降低了算法的稳定性与鲁棒性。另一方面，在选择地面控制点时，应考虑如下几个因素：应选择场景中具有显著定位标志的位置，例如道路的交叉口、建筑物边界等；应确保同名点对应的地物不随时间的变化而变化；应使地面控制点合理的分布于整个场景中，以免在控制点密集区陷入局部最优。尤其在高分辨率遥感影像内，受到大气折射、地形起伏，地球曲率等诸多因素的影响，场景中通常存在不规则的局部形变。为保证配准精度及减少局部形变对配准造成的影响，内部纹理特征丰富、结构复杂的对象应采用较多的同名点对进行标记；而对于纹理特征单一、结构简单的局部区域提取少量同名点对即可。现有多尺度配准算法约束控制点分布时，通常选择在某一特定尺寸窗口内提取一定数量的控制点，因此难以保证同名点集的合理分布。最后，目前基于点特征的影像配准中所采用的控制点匹配策略通常只考虑了局部区域的灰度统计特征，而忽略了原始影像中的纹理、尺度等特征，容易受到噪声等因素的干扰，降低同名点集的可靠性，增加了误匹配的可能。

发明内容

发明目的：针对现有技术中存在的问题，本发明提供了一种控制点自适应分布的高分辨率遥感影像配准方法，主要包括如下三个步骤：(1)多尺度控制点提取。(2)基于NMI的同名点匹配。(3)基于Delaunay三角形局部变换的几何校正。

技术方案：一种控制点自适应分布的高分辨率遥感影像配准方法，将多尺度J-image引入到图像配准中，首先提出了一种基于分块策略的自适应控制点提取方法，实现了多尺度J-image影像中的控制点提取，从而克服了控制点仅对个别方向高频信息敏感的局限。同时，定义了自适应控制点提取策略以约束控制点的分布。进而采用归一化的互信息测度NMI对控制点进行多尺度进行匹配，能够有效平滑配准函数。最后采用Delaunay三角形局部变换实现几何校正。下文首先介绍了所提出算法的基本原理及关键步骤，进而对三组不同类型的遥感影像进行了实验和分析，并分别与基于小波变换及NSCT的不同配准方法进行了比较实验。

JSEG与NMI相结合的配准方法

多尺度控制点提取

选择J-image影像作为多尺度分析工具进行控制点提取与匹配。

1计算J-image影像序列

利用某一特定尺寸的窗口计算原始影像中每个像素对应的局部同质性指标J-value并作为该像素的像素值，可以获得单一尺度的J-image影像。J-value的计算过程如下：

首先采用Deng等人提出的量化方法(Deng Y,Kenney C,Moore M S,et al.Peergroup filtering and perceptual color image quantization[C]//Proceedings ofthe 1999 IEEE International Symposium on Circuits and Systems,1999,4:21-24.)对原始影像进行量化从而获得量化影像。令量化影像中每一个像素的位置z(x,y)为像素z的像素值，z(x,y)∈Z。Z为以像素z为中心的特定尺寸窗口内所有像素组成的集合。

在量化影像中，定义N为以z为中心窗口中的像素总数，则均值m：

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差的和S_W可定义为：

定义S_T为窗口中所有像素的总体方差：

则J-value为：

J＝(S_T-S_W)/S_W (2.4)

利用不同尺寸窗口分别计算不通尺度下z的J-value，并作为z的像素值，从而获得多尺度J-image影像序列。根据以上J-value的定义，在某一尺度的J-image中，某一像素的对应的J-value值越大，则越可能处于对象的边缘；反之，则可能位于对象的中心。同时，J-value综合了原始影像所包含的光谱信息、纹理信息与尺度信息，且对方向信息不敏感，因此也可作为控制点的描述子用于同名点的匹配。

2控制点提取及分布约束

根据J-value的定义及特点，提出了一种基于分块策略的自适应控制点提取方法，其核心思想是在纹理特征丰富的区域提取更多的控制点，实现流程如下：

Step1：将各个尺度的J-image影像等分为不同尺寸的子图像，某尺度中子图像大小与计算该尺度J-image影像的特定窗口尺寸大小相同，但保留角点，为N×N像素的正方形窗口。

Step2：在原始影像中利用公式(2.4)计算整幅影像的J-value值，定义为T_a。T_a反映了原始影像的总体均质程度。

Step3：在某一尺度J-image中，分别将每个子图像中心像素的J-valve与T_a进行比较，若大于T_a，则认为该子图像内部纹理特征丰富，结构复杂，应提取较多的控制点，取子图像中J-value最大的前K个像素作为控制点；反之，则提取少量控制点即可，取子图像中J-value最大的前L个像素作为控制点，并且满足K和L的值可根据影像实际尺寸及影像纹理的复杂程度人工设定。

Step4：在所有尺度J-image影像中重复Step3的操作，实现控制点的多尺度提取及分布约束。

基于NMI的控制点匹配

在完成多尺度下的控制点提取后，需要进一步利用相似性度量对每个控制点逐一判别获得同名点集。控制点间的相似性度量通常采用以控制点为中心的局部区域的灰度统计特性作为相似性测度。因此，基于局部图像灰度特征的控制点相似性度量可以采用一些通用的基于灰度特征的图像配准方法。

互信息(MI，Mutual Information)是目前图像配准中广泛应用的一种可回溯性相似性测度，具有无需预分割、鲁棒性强以及精度高等特点。由于无需进行预处理，即使在图像数据有部分缺失的应用场合中，依然能够取得良好的配准结果，同时也可用于多源遥感影像的配准中。另一方面，互信息相似性测度对图像用重叠区域的应用场合过于敏感，而当重叠区域减少时，意味着更少的像素参与了互信息的统计，从而降低了配准的可靠性。为克服这一局限，Studholme C等提出的归一化的互信息测度(NMI，Normalized MutualInformation)能够有效平滑配准函数，减少图像配准中互信息对重叠区域的敏感程度，使目标函数更加精确地反映配准参数与互信息量之间的关系，从而提高配准精度。因此，本文选择NMI作为控制点相似性度量提取各尺度下的同名点对，获得全局同名点集。

1计算NMI相似性度量

在某一尺度J-image中，将参考影像和待配准影像中的控制点所在区域作为两个矢量C和D，计算控制点间的NMI。所采用的窗口尺寸与计算当前尺度J-image所采用的特定窗口尺寸相同。首先用熵来描述互信息，熵是指平均意义上对信源来表征的总体特征，可通过对可能时间上的信息量加权平均获得。定义熵H及互信息I(C,D)如公式(2.5)、(2.6)所示：

I(C,D)＝H(C)+H(D)-H(C,D) (2.6)

其中，P_i表示某个随机变量所有可能中第i种可能情况发生的概率。则H(C)、H(D)、H(C,D)分别为参考影像和待配准影像中同名点对应的熵以及他们的联合熵。根据公式(2.5)可知：

其中P_C(c)和P_D(d)分别为C和D完全独立时的概率分布，P_CD为C和D的联合概率分布。NMI定义如下：

在参考影像及待配准影像的同一尺度J-image中，以控制点为中心，以计算当前尺度J-value的特定窗口为局部区域生成子图像，计算子图像间的NMI从而度量控制点间的相似性。需要指出的是，尽管计算NMI的过程仅对控制点所在区域的局部灰度特征进行了统计，但由于J-image影像中每个像素的灰度值为原始影像中每个像素对应的J-value值，因此根据J-value的定义可知，所计算的NMI相似性度量综合反映了控制点所在局部区域的灰度、纹理以及光谱特征的相似性，从而更加可靠。

2基于NMI的最大双向匹配策略

为实现控制点间的准确匹配，提出了一种基于NMI的最大双向匹配策略，匹配过程如下：

Step1：在参考影像与待配准影像的同一尺度J-image中，NMI最大且大于阈值T_NMI的控制点对作为一对同名点，进而得到同名点集1。

Step2：与Step1计算方向相反，将待配准影像中某一控制点分别与参考影像中所有控制点进行比较，得到同名点集2.

Step3：比较同名点集1和同名点集2，将相同的同名点对作为当前尺度下获得的最终同名点集。

Step4：重复Step1至Step3，对所有尺度J-image进行判别，并将获得的同名点对全部映射到原始影像中的相同位置。

Step5：在原始影像中采用经典的RANSAC(Random Sample Consensus)算法进一步消除误匹配现象，获得最终的同名点集。

基于Delaunay三角形的图像配准

由于高分辨率遥感影像中地物种类众多，纹理信息丰富，通常存在大量的局部形变，因而单纯采用基于全局多项式的映射函数难以保证配准的精度，而局部三角形变换采用局部映射函数，能够有效解决这一问题。与普通三角网相比，Rognant等人提出的Delaunay三角网具有拓扑关系始终保持不变等优点，即无论从何处开始构造，所有控制点构造的Delaunay三角网格都是一致的。因而，本文采用Delaunay三角网构造几何映射函数，对待配准的影像进行重采样。映射函数为：

其中，(x,y)为参考影像中局部三角形的坐标，(x',y')为待配准影像中对应局部三角形的坐标，由于构建三角形的三个控制点坐标已知，因而可直接求出(2.11)式中待定系数。

本发明采用上述技术方案，具有以下有益效果：与基于小波变换及NSCT等传统多尺度分析工具的配准方法相比，本发明采用J-image序列作为多尺度分析工具，克服了控制点提取时采用不同方向高频子带对配准结果造成的影响。与此同时，所提出的基于分块策略的自适应控制点提取策略实现了控制点的准确提取及分布约束。而基于NMI的多尺度同名点匹配策略综合利用了原始影像中光谱、纹理以及尺度特征，能够有效消除误匹配现象。实验证明这样的控制点提取及匹配策略保证了获得的同名点对根据局部区域的纹理及结构特征不同而合理的分布。最后，通过引入Delaunay三角形局部变换实现了图像配准。实验采用三组不同类型的高分辨率遥感影像，并与两种基于传统多尺度分析工具的配准方法进行了比较。实验证明所提出算法具有更高的配准精度，且尤其适用于场景中存在大量图像畸变的应用场合，能够为后续的变化检测提供可靠的图像配准预处理。

附图说明

图1为本发明实施例的方法流程图；

图2为本发明实施例的9×9像素窗口；

图3为本发明实施例的18×18像素窗口；

图4为本发明实施例的参考影像#1

图5为本发明实施例的待配准影像#2；

图6为本发明实施例的方法提取的84对同名点；

图7为基于小波变换方法提取的43对同名点；

图8为基于NSCT方法提取的53对同名点；

图9为本发明实施例的方法配准结果；

图10为基于小波变换方法配准结果；

图11为基于NSCT方法配准结果；

图12为本发明实施例的参考影像#3；

图13为本发明实施例的待配准影像#4；

图14为本发明实施例的方法提取的61对同名点；

图15为基于小波变换方法提取的43对同名点；

图16为基于NSCT方法提取的53对同名点；

图17为本发明实施例的方法配准结果；

图18为基于小波变换方法配准结果；

图19为基于NSCT方法配准结果；

图20为本发明实施例的参考影像#5；

图21为本发明实施例的待配准影像#6；

图22为本发明实施例的方法提取的30对同名点；

图23为基于小波变换方法提取的23对同名点；

图24为基于NSCT方法提取的26对同名点；

图25为本发明实施例的方法配准结果；

图26为基于小波变换方法配准结果；

图27为基于NSCT方法配准结果。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，控制点自适应分布的高分辨率遥感影像配准方法，主要包括如下三个步骤：(1)多尺度控制点提取。(2)基于NMI的同名点匹配。(3)基于Delaunay三角形局部变换的几何校正。

将多尺度J-image引入到图像配准中，首先提出了一种基于分块策略的自适应控制点提取方法，实现了多尺度J-image影像中的控制点提取，从而克服了控制点仅对个别方向高频信息敏感的局限。同时，定义了自适应控制点提取策略以约束控制点的分布。进而采用归一化的互信息测度NMI对控制点进行多尺度进行匹配，能够有效平滑配准函数。最后采用Delaunay三角形局部变换实现几何校正。下文首先介绍了所提出算法的基本原理及关键步骤，进而对三组不同类型的遥感影像进行了实验和分析，并分别与基于小波变换及NSCT的不同配准方法进行了比较实验。

JSEG与NMI相结合的配准方法

多尺度控制点提取

选择J-image影像作为多尺度分析工具进行控制点提取与匹配。

1计算J-image影像序列

利用某一特定尺寸的窗口计算原始影像中每个像素对应的局部同质性指标J-value并作为该像素的像素值，可以获得单一尺度的J-image影像。J-value的计算过程如下：

首先采用Deng等人提出的量化方法(Deng Y,Kenney C,Moore M S,et al.Peergroup filtering and perceptual color image quantization[C]//Proceedings ofthe 1999 IEEE International Symposium on Circuits and Systems,1999,4:21-24.)对原始影像进行量化从而获得量化影像。令量化影像中每一个像素的位置z(x,y)为像素z的像素值，z(x,y)∈Z。Z为以像素z为中心的特定尺寸窗口内所有像素组成的集合。例如，采用图2所示的窗口作为第一尺度，以像素z为中心的特定窗口尺寸为9×9像素。则在第二尺度中，采用以z为中心18×18像素窗口，如图3所示。同时，为进一步减少计算量，图3中只有“+”代表的像素用于计算当前尺度下像素z的J-value，计算各尺度下的J-value依此类推。为尽量保证各个方向的一致性，窗口中的角点被去除。

在量化影像中，定义N为以z为中心窗口中的像素总数，则均值m：

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差的和S_W可定义为：

定义ST为窗口中所有像素的总体方差：

则J-value为：

J＝(S_T-S_W)/S_W (2.4)

利用不同尺寸窗口分别计算不通尺度下z的J-value，并作为z的像素值，从而获得多尺度J-image影像序列。根据以上J-value的定义，在某一尺度的J-image中，某一像素的对应的J-value值越大，则越可能处于对象的边缘；反之，则可能位于对象的中心。同时，J-value综合了原始影像所包含的光谱信息、纹理信息与尺度信息，且对方向信息不敏感，因此也可作为控制点的描述子用于同名点的匹配。

2控制点提取及分布约束

根据J-value的定义及特点，提出了一种基于分块策略的自适应控制点提取方法，其核心思想是在纹理特征丰富的区域提取更多的控制点，实现流程如下：

Step1：将各个尺度的J-image影像等分为不同尺寸的子图像，某尺度中子图像大小与计算该尺度J-image影像的特定窗口尺寸大小相同，但保留角点，为N×N像素的正方形窗口。

Step2：在原始影像中利用公式(2.4)计算整幅影像的J-value值，定义为T_a。T_a反映了原始影像的总体均质程度。

Step3：在某一尺度J-image中，分别将每个子图像中心像素的J-valve与T_a进行比较，若大于T_a，则认为该子图像内部纹理特征丰富，结构复杂，应提取较多的控制点，取子图像中J-value最大的前K个像素作为控制点；反之，则提取少量控制点即可，取子图像中J-value最大的前L个像素作为控制点，并且满足K和L的值可根据影像实际尺寸及影像纹理的复杂程度人工设定。

Step4：在所有尺度J-image影像中重复Step3的操作，实现控制点的多尺度提取及分布约束。

基于NMI的控制点匹配

在完成多尺度下的控制点提取后，需要进一步利用相似性度量对每个控制点逐一判别获得同名点集。控制点间的相似性度量通常采用以控制点为中心的局部区域的灰度统计特性作为相似性测度。因此，基于局部图像灰度特征的控制点相似性度量可以采用一些通用的基于灰度特征的图像配准方法。

互信息(MI，Mutual Information)是目前图像配准中广泛应用的一种可回溯性相似性测度，具有无需预分割、鲁棒性强以及精度高等特点。由于无需进行预处理，即使在图像数据有部分缺失的应用场合中，依然能够取得良好的配准结果，同时也可用于多源遥感影像的配准中。另一方面，互信息相似性测度对图像用重叠区域的应用场合过于敏感，而当重叠区域减少时，意味着更少的像素参与了互信息的统计，从而降低了配准的可靠性。为克服这一局限，Studholme C等提出的归一化的互信息测度(NMI，Normalized MutualInformation)能够有效平滑配准函数，减少图像配准中互信息对重叠区域的敏感程度，使目标函数更加精确地反映配准参数与互信息量之间的关系，从而提高配准精度。因此，本文选择NMI作为控制点相似性度量提取各尺度下的同名点对，获得全局同名点集。

1计算NMI相似性度量

在某一尺度J-image中，将参考影像和待配准影像中的控制点所在区域作为两个矢量C和D，计算控制点间的NMI。所采用的窗口尺寸与计算当前尺度J-image所采用的特定窗口尺寸相同。首先用熵来描述互信息，熵是指平均意义上对信源来表征的总体特征，可通过对可能时间上的信息量加权平均获得。定义熵H及互信息I(C,D)如公式(2.5)、(2.6)所示：

I(C,D)＝H(C)+H(D)-H(C,D) (2.6)

其中，P_i表示某个随机变量所有可能中第i种可能情况发生的概率。则H(C)、H(D)、H(C,D)分别为参考影像和待配准影像中同名点对应的熵以及他们的联合熵。根据公式(2.5)可知：

其中P_C(c)和P_D(d)分别为C和D完全独立时的概率分布，P_CD为C和D的联合概率分布。NMI定义如下：

在参考影像及待配准影像的同一尺度J-image中，以控制点为中心，以计算当前尺度J-value的特定窗口为局部区域生成子图像，计算子图像间的NMI从而度量控制点间的相似性。需要指出的是，尽管计算NMI的过程仅对控制点所在区域的局部灰度特征进行了统计，但由于J-image影像中每个像素的灰度值为原始影像中每个像素对应的J-value值，因此根据J-value的定义可知，所计算的NMI相似性度量综合反映了控制点所在局部区域的灰度、纹理以及光谱特征的相似性，从而更加可靠。

2基于NMI的最大双向匹配策略

为实现控制点间的准确匹配，提出了一种基于NMI的最大双向匹配策略，匹配过程如下：

Step1：在参考影像与待配准影像的同一尺度J-image中，NMI最大且大于阈值T_NMI的控制点对作为一对同名点，进而得到同名点集1。

Step2：与Step1计算方向相反，将待配准影像中某一控制点分别与参考影像中所有控制点进行比较，得到同名点集2.

Step3：比较同名点集1和同名点集2，将相同的同名点对作为当前尺度下获得的最终同名点集。

Step4：重复Step1至Step3，对所有尺度J-image进行判别，并将获得的同名点对全部映射到原始影像中的相同位置。

Step5：在原始影像中采用经典的RANSAC(Random Sample Consensus)算法进一步消除误匹配现象，获得最终的同名点集。

基于Delaunay三角形的图像配准

由于高分辨率遥感影像中地物种类众多，纹理信息丰富，通常存在大量的局部形变，因而单纯采用基于全局多项式的映射函数难以保证配准的精度，而局部三角形变换采用局部映射函数，能够有效解决这一问题。与普通三角网相比，Rognant等人提出的Delaunay三角网具有拓扑关系始终保持不变等优点，即无论从何处开始构造，所有控制点构造的Delaunay三角网格都是一致的。因而，本文采用Delaunay三角网构造几何映射函数，对待配准的影像进行重采样。映射函数为：

其中，(x,y)为参考影像中局部三角形的坐标，(x',y')为待配准影像中对应局部三角形的坐标，由于构建三角形的三个控制点坐标已知，因而可直接求出(2.11)式中待定系数。

实验结果与分析

为了验证所提算法的有效性，本文分别对三组来自不同传感器的高分辨率遥感影像数据集进行了实验。另外，实验结果将分别与采用基于小波变换的配准方法(文献Zavorin I,Le Moigne J.Use of multi-resolution wavelet feature pyramids forautomatic registration of multi-sensor imagery[J].Image Processing,IEEETransactions on,2005,14(6):770-782.)以及基于NSCT的配准方法(文献Zhigang C,Fuchang Y,Fu S.Registration technique for high-resolution remote sensingimages based on non-subsampled contourlet transform[J].Acta Optica Sinica,2009,29(10):2744-2750.)进行比较。其中，基于小波变换的配准方法采用全局多项式拟合获得配准影像，基于NSCT的配准方法则采用Delaunay三角网构建映射函数。所提出的方法实验参数设定均通过多组实验选择最优值，其他方法涉及参数均参考原文献设定。

所采用的评价指标为同名点对数量、配准后影像与参考影像的相对相关系数(Comparative Correlation)ρ_corr、均方根误差(Root-Mean-Square Error)R_MSE以及归一化均方根误差(Normalized-Root-Mean-Square Error)NR_MSE，定义如下：

相对相关系数，

均方根误差，

归一化均方根误差，

其中，N表示配准后影像的像素个数，x_i为配准后影像中像素z_i(i＝1,2...N)的像素值，为z_i在参考影像中对应像素的像素值。在各评价指标中，相对相关系数越大，说明两幅影像越接近，理想值为1。均方根误差及归一化均方根误差反映了样本的离散程度，参数值越小，两幅图像越接近，则配准精度越高。

数据集1实验结果与分析

选择数据一组影像#1、#2作为数据集1进行实验，如图4，图5所示。参考影像#1与待配准影像#2分别为2009年3月及2012年2月获取的航空遥感DOM(Digital Orthophoto Map)数据，所在地区为中国江苏省南京市河海大学江宁校区，空间分辨率为0.5m，图像大小为512×512像素。

如图4，图5所示，由于两时相影像均为由多幅不同视角的航空影像合成的DOM影像，各局部区域视角差异明显，从而造成了大量局部形变，这一现象尤其在结构复杂的人造目标及高层建筑中更为突出，因此需要更多的控制点进行校正。在实验中，计算J-value的窗口被设定为20×20像素，10×10像素和5×5像素，即计算三个尺度的J-image影像。设定参数K＝6，L＝3，T_NMI＝0.85。基于小波变换及NSCT的配准方法同样采用三级分解。三种方法提取的同名点集中，同名点对分布如图6～图8所示。

三种方法所获得的配准结果如图9～图11所示。

如图4～图11所示，根据目视分析可以得出如下结论：(1)本发明提出的方法提取的同名点对最多，且依据局部区域的结构与纹理特征不同合理分布于图像内部。由于本文提出了基于分块策略的自适应控制点提取方法，对于内部纹理特征丰富、结构复杂的人造建筑等区域，大量控制点标记了对象的细节特征，从而减少了视角差异等因素导致的局部形变所造成的影响；对于内部纹理均匀、结构单一的道路、湖面等区域，少量控制点即准确标记了对象所在的位置。(2)由于基于小波变换与NSCT的配准方法在控制点提取及同名点匹配过程中都分别采用了单一的阈值进行判别，而忽略了图像中局部区域纹理特征与结构的复杂程度，因而提取到的控制点分布较为分散，难以保证最优。(3)基于小波变换的配准方法由于采用了全局多项式拟合几何校正方法，因而配准结果中对大量局部性形变区域的校正效果较差。(4)基于NSCT的配准方法采用了基于三角形变换的几何校正方法，因而对图像中的局部形变较为敏感。但少量控制点难以充分标记发生局部形变的区域在不同时相影像中的形状及结构差异，从而影响最终的配准精度。另一方面，若单纯改变阈值来增加提取控制点与同名点对的数量，又会造成计算量的指数级增长以及误匹配现象的增加。

为进一步定量分析算法的性能，精度评价结果如表1所示。通过表1可以看出：(1)本发明提出的方法提取到了更多的同名点对，且配准精度明显高于其他两种算法，与目视分析结果一致。(2)本发明方法较其他两种方法相比，由于采用的J-image序列作为多尺度分析工具，不存在高频子带的方向选择，因而不会忽略某一方向上的高频信息，从而提取的控制点具有更好的准确性。(3)对于局部形变较大的高分辨率遥感影像配准，采用基于三角形变换的几何校正方法优于全局多项式拟合方法。

表2.1本发明方法与其他方法的比较

数据集2实验结果与分析

数据集2采用空间分辨率为5m的SPOT 5全色—多光谱融合影像，参考影像#3、待配准影像#4(如图12，图13所示)，尺寸均为512×512像素。融合波段包括全色波段以及红、绿、近红外波段。#3与#4获取时间分别为2004年6月及2008年7月，所在地区为中国上海。

与数据集1相比，数据集2两幅影像空间分辨更低，视角差异很小，且地物种类丰富，背景更加复杂。在实验中，计算J-value的窗口被设定为10×10像素，8×8像素和5×5像素，即计算三个尺度的J-image影像。设定参数K＝6，L＝3，T_NMI＝0.9。基于小波变换及NSCT的配准方法同样采用三级分解。三种方法提取的同名点集中，同名点对分布如图14～图16所示。

三种方法所获得的配准结果如图17～图19所示。

精度评价结果如表2所示：

表2本发明与其他方法的比较

通过图14～图19及表2可以看出，目视分析结果与精度评价结果与数据集1基本相同，进一步验证了本文算法有效性。另外，数据集2中三种算法的配准精度较数据集1均有所提高，主要原因是数据集2中视角差异更小，从而造成的局部形变区域更少。通过比较表1及表2可以看出，基于全局多项式拟合的校正方法仅在受到局部形变影响较小的图像配准中能够取得较好的效果。而采用三角形变换的校正方法除了对同名点对的准确性要求较高外，同名点对的数量及分布也会对配准精度产生重要影响。

数据集3实验结果与分析

在前两组实验中，数据集1为航空遥感影像，数据集2为卫星遥感影像。为进一步验证所提出算法对高分辨率SAR影像配准的效果，选择ENVISAT卫星拍摄的ASAR(AdvancedSynthetic Aperture Radar)参考影像#5、待配准影像#6作为数据集3进行实验，如图20，图21所示。数据集3所在地为中国武汉地区，空间分辨率为10m，采集时间分别为2008年2月和2009年8月，尺寸为512×512像素。

比较影像#5、#6可以发现，数据集3图像数据主要由陆地、水体、湖泊组成构成。尽管影像采集时间间隔仅三个月，但水域边界已发生了显著变化，场景中缺乏具有显著定位标志的对象，同时存在大量相干斑噪声，这些因素都给地理位置的对准增加了困难。在实验中，计算J-value的窗口被设定为10×10像素，8×8像素和5×5像素，即计算三个尺度的J-image影像。设定参数K＝6，L＝3，T_NMI＝0.85。基于小波变换及NSCT的配准方法同样采用三级分解。三种方法提取的同名点集中同名点对分布如图22～图24所示。

三种方法所获得的配准结果如图25～图27所示。

表3本发明与其他方法的比较

通过图22～图27及表3可以看出，在高分辨率SAR影像配准中，本发明提出方法提取的同名点对数量最多且分布更加均匀，同时各精度指标也优于其他两种算法。同时，通过比较三种算法在三组数据集实验中的配准精度可以看出，SAR影像中大量存在的相干斑噪声没有显著降低三种算法的配准精度。这也进一步证明了基于多尺度分析工具的遥感影像配准方法对噪声干扰具有良好的鲁棒性。

Claims (1)

1.一种控制点自适应分布的高分辨率遥感影像配准方法，其特征在于：主要包括多尺度控制点提取、基于NMI的控制点匹配和基于Delaunay三角形的图像配准三个步骤；

多尺度控制点提取

选择J-image影像作为多尺度分析工具进行控制点提取与匹配；

1计算J-image影像序列

利用某一特定尺寸的窗口计算原始影像中每个像素对应的局部同质性指标J-value并作为该像素的像素值，获得单一尺度的J-image影像，J-value的计算过程如下：

首先采用量化方法对原始影像进行量化从而获得量化影像；令量化影像中每一个像素的位置z(x,y)为像素z的像素值，z(x,y)∈Z，Z为以像素z为中心的特定尺寸窗口内所有像素组成的集合；

在量化影像中，定义N为以z为中心窗口中的像素总数，则均值m：

$m=\frac{1}{N}\underset{z\∈Z}{\Σ}z(x,y)---\left(2.1\right)$

定义m_p为窗口中属于相同灰度级p的所有像素均值，Z_p为窗口中属于灰度级p的所有像素的集合，P为量化影像中的灰度级总数，则窗口中属于同一灰度级像素的方差的和S_W定义为：

$S_W=\underset{p=1}{\overset{P}{\Σ}}\underset{z\∈Z_p}{\Σ}\left|\right|z-m_p|{|}^2---\left(2.2\right)$

定义S_T为窗口中所有像素的总体方差：

$S_T=\underset{z\∈Z}{\Σ}\left|\right|z-m|{|}^2---\left(2.3\right)$

则J-value为：

J＝(S_T-S_W)/S_W (2.4)

利用不同尺寸窗口分别计算不同尺度下z的J-value，并作为z的像素值，从而获得多尺度J-image影像序列；根据以上J-value的定义，在某一尺度的J-image中，某一像素的对应的J-value值越大，则越可能处于对象的边缘；反之，则可能位于对象的中心；同时，J-value综合了原始影像所包含的光谱信息、纹理信息与尺度信息，且对方向信息不敏感，因此也作为控制点的描述子用于同名点的匹配；

2控制点提取及分布约束

基于分块策略的自适应控制点提取方法，实现流程如下：

Step1-1：将各个尺度的J-image影像等分为不同尺寸的子图像，某尺度中子图像大小与计算该尺度J-image影像的特定窗口尺寸大小相同，但保留角点，为N×N像素的正方形窗口；

Step1-2：在原始影像中利用公式(2.4)计算整幅影像的J-value值，定义为T_a；T_a反映了原始影像的总体均质程度；

Step1-3：在某一尺度J-image中，分别将每个子图像中心像素的J-valve与T_a进行比较，若大于T_a，则认为该子图像内部纹理特征丰富，结构复杂，应提取较多的控制点，取子图像中J-value最大的前K个像素作为控制点；反之，则提取少量控制点即可，取子图像中J-value最大的前L个像素作为控制点，并且满足K和L的值根据影像实际尺寸及影像纹理的复杂程度设定；

Step1-4：在所有尺度J-image影像中重复Step1-3的操作，实现控制点的多尺度提取及分布约束；

基于NMI的控制点匹配

1计算NMI相似性度量

在某一尺度J-image中，将参考影像和待配准影像中的控制点所在区域作为两个矢量C和D，计算控制点间的NMI；所采用的窗口尺寸与计算当前尺度J-image所采用的特定窗口尺寸相同；首先用熵来描述互信息；定义熵H及互信息I(C,D)如公式(2.5)、(2.6)所示：

H＝-Σ_iP_ilog₂P_i (2.5)

I(C,D)＝H(C)+H(D)-H(C,D) (2.6)

其中，P_i表示某个随机变量所有可能中第i种可能情况发生的概率；则H(C)、H(D)、H(C,D)分别为参考影像和待配准影像中同名点对应的熵以及他们的联合熵；根据公式(2.5)可知：

H(C)＝-∑_cP_C(c)log₂P_C(c) (2.7)

H(D)＝-∑_dP_D(d)log₂P_D(d) (2.8)

H(C,D)＝-∑_c,dP_CD(c,d)log₂P_CD(c,d) (2.9)

其中P_C(c)和P_D(d)分别为C和D完全独立时的概率分布，P_CD为C和D的联合概率分布；NMI定义如下：

$NMI=\frac{H\left(C\right)+H\left(D\right)}{H(C,D)}---\left(2.10\right)$

在参考影像及待配准影像的同一尺度J-image中，以控制点为中心，以计算当前尺度J-value的特定窗口为局部区域生成子图像，计算子图像间的NMI从而度量控制点间的相似性；

2基于NMI的最大双向匹配策略

为实现控制点间的准确匹配，提出了一种基于NMI的最大双向匹配策略，匹配过程如下：

Step1：在参考影像与待配准影像的同一尺度J-image中，NMI最大且大于阈值T_NMI的控制点对作为一对同名点，进而得到同名点集1；

Step2：与Step1计算方向相反，将待配准影像中某一控制点分别与参考影像中所有控制点进行比较，得到同名点集2；

Step3：比较同名点集1和同名点集2，将相同的同名点对作为当前尺度下获得的最终同名点集；

Step4：重复Step1至Step3，对所有尺度J-image进行判别，并将获得的同名点对全部映射到原始影像中的相同位置；

Step5：在原始影像中采用经典的RANSAC算法进一步消除误匹配现象，获得最终的同名点集；

基于Delaunay三角形的图像配准

采用Delaunay三角网构造几何映射函数，对待配准的影像进行重采样；映射函数为：

$\left\{\begin{array}{c}{x}^{\′}=a_0+a_{1}x+a_{2}y\\ y^{\′}=b_0+b_{1}x+b_{2}y\end{array}\right.---\left(2.11\right)$

其中，a₀，a₁，a₂，b₀，b₁，b₂为参考影像与待配准影像间映射函数中的待定系数，(x,y)为参考影像中局部三角形的坐标，(x',y')为待配准影像中对应局部三角形的坐标，由于构建三角形的三个控制点坐标已知，因而直接求出(2.11)式中待定系数。