CN101251926A - 一种基于局部轮廓协方差矩阵的遥感图像配准方法 - Google Patents

Abstract

本发明提出了一种基于局部轮廓协方差矩阵的遥感图像配准方法，将角点特征、局部分割图像和局部边缘轮廓结合起来，共同作为局部特征来实现控制点提取，并将局部轮廓协方差矩阵引入配准处理，利用局部轮廓协方差矩阵来实现旋转缩放因子的自动提取，结合局部窗口自适应选取技术来减小旋转缩放因子对配准处理的影响，提高配准处理的鲁棒性，在无人工参与条件下实现大旋转缩放模式下遥感图像间的精确配准。

Description

一种基于局部轮廓协方差矩阵的遥感图像配准方法

技术领域

本发明属于图像处理领域，涉及一种图像配准方法，具体涉及一种基于局部轮廓协方差矩阵的遥感图像配准处理方法。

背景技术

图像配准是图像处理的一个基础问题，是对取自不同时间、不同传感器或不同视角的同一场景的两幅或多幅图像进行匹配、叠加的过程。准确的说，图像配准的目标就是找到输入图像之间的最佳映射关系。

目前，图像配准处理方法主要可以分为三类：

1、基于全局灰度统计信息的配准处理方法；

基于全局灰度统计信息的配准处理方法直接利用图像本身具有的灰度统计特性来度量图像间的相似性，进而实现图像间的配准处理。该处理方法最大优点在于原理简单、易于实现，对于同类图像且仿射误差较小的情况下能够获得很好的处理效果。然而对于异类遥感图像，由于图像灰度特性存在较大的差异，其处理效果并不理想。同时，由于处理过程中采用迭代搜索，当仿射误差较大时，其处理速度急剧下降，甚至无法收敛。

2、基于相位相关技术的配准处理方法；

基于相位相关技术的配准处理方法的理论基础是傅立叶相位相关技术，利用频域信息进行相关处理，搜索最佳匹配，进而在不需要寻找控制点的情况下实现图像自动配准。该算法的优点在于对噪声不敏感，由于可以采用现成的FFT算法，因此，该配准处理算法在处理速度上具有一定的优势。在实际使用中，该算法需要待配准图像间具有较大的重叠区域，且图像的灰度特征相似，这限制了该处理算法在遥感图像配准过程中的使用。

3、基于图像特征的配准处理方法；

基于图像特征的配准处理方法是图像配准方法中的另一大类，如图1所示，该方法从各图象通道中读入各图像后，首先提取图像中的特征信息，利用提取的特征来进行局部相似性匹配，并通过特征的匹配关系建立图像之间的映射关系，最后通过图像重采样处理来获取匹配的图像对。

相比较而言，由于前两类处理方法过多的依赖于图像的灰度特性，制约了它们在遥感图像配准处理中的应用，因此，基于特征的配准处理方法成为遥感图像配准处理的常用方法，其核心在于建立图像特征间的映射关系。

在特征匹配关系建立方面，传统的方法是通过***部相关系数或交互信息等统计参数来衡量图像中两点是否一致，但该方法受旋转因子和缩放因子的影响很大，随着缩放因子和旋转因子的增大，其相关性迅速降低。利用仿射不变量进行相似性测试，则依赖于提取完整封闭的边缘轮廓，而这对于遥感图像尤其是雷达图像来说通常很难实现。利用点特征之间的聚集性来实现相似性匹配则要求各输入图像所提取的点特征具有很高的一致性，对于遥感图像来说，这同样很难实现。因此，如何有效的建立输入图像间的映射关系，成为遥感图像配准处理的难点。

发明内容

本发明提出了一种基于局部轮廓协方差矩阵的遥感图像配准方法，将角点特征、局部分割图像和局部边缘轮廓结合起来，共同作为局部特征来实现控制点提取，并将局部轮廓协方差矩阵引入配准处理，利用局部轮廓协方差矩阵来实现旋转缩放因子的自动提取，结合局部窗口自适应选取技术来减小旋转缩放因子对配准处理的影响，提高配准处理的鲁棒性，在无人工参与条件下实现大旋转缩放模式下遥感图像间的精确配准。

首先提取输入图像中的角点特征，以角点特征为参考，提取角点附近的局部轮廓特征，并利用局部轮廓特征的协方差矩阵提取输入图像间的旋转缩放因子，减小旋转缩放因子对特征匹配的影响，提高同名点的提取精度；然后，利用检测到的同名点来进行仿射模型参数估计；最后，根据估算出来的变换模型参数进行图像重采样，获得配准后的遥感图像对。从整体上来看，该处理流程主要包含特征的选择与提取、旋转缩放因子的自动提取、最优窗口尺寸选取、角点特征匹配、仿射模型参数估计及图像重采样五个部分。

一种基于局部轮廓协方差矩阵的遥感图像配准方法，包括如下步骤：

步骤1：从各图像通道中同时读入各图像，采用Harris角点检测算法实现图像角点检测，通过统计移动搜索窗口内各点的角响应函数值来判断是否存在角点；

如果一点的角响应函数CRF值是局部邻域内的最大值，则将该点设置为备选角点；否则表明该点不是角点。在判断完图像中所有点后，利用门限处理来提取最终角点。

步骤2：采用类间方差自动阈值法对参考图像及待配准图像进行局部阈值分割；利用提取得最佳阈值对局部图像进行分割，获得局部二值图像，并依据提取的局部二值图像获得局部边缘轮廓。

步骤3：旋转缩放因子的自动提取；采用基于局部轮廓协方差矩阵的旋转缩放因子自动提取方法。利用阈值分割所获得的二值图像，分别计算参考图像与待配准图像的协方差矩阵；依据参考图像与待配准图像的局部结构分布即可提取图像间的旋转因子和缩放因子。

步骤4：通过对比分析参考图像和待配准图像中选取窗口尺寸的比值和利用输入图像提取的缩放因子来选取分析窗口尺寸，提取最优窗口。

步骤5：利用提取的旋转缩放因子和最优窗口尺寸对步骤1中提取的边缘轮廓进行重采样，减小旋转缩放因子对特征匹配的影响；并对重采样后的边缘特征进行Freeman编码，并通过链码相关技术统计两线段间的相似性；利用点集间的相对位置关系和点集间的相互支持特性剔除错误匹配点。

步骤6：当完成特征点匹配后，利用最小二乘法估算仿射模型中的参数；在确定变换参数后，依据仿射参数对输入图像对进行重采样处理，可获得相互匹配的图像对。

所述步骤3中，针对局部分割图像展开旋转缩放因子的自动提取。通过阈值分割将图像分割为亮、暗两类图像，然后分别统计两部分的协方差矩阵及旋转缩放因子，确定图像间的对应关系，并提取相应的旋转缩放因子。

所述步骤3中，确定图像间的对应关系时，对比分析参考图像中亮、暗图像对应待配准图像中的亮、暗图像及参考图像中亮、暗图像对应待配准图像中的暗、亮图像这两种对应关系，选取旋转角度差和缩放因子的比例，选取旋转角度差更小，两轴旋转缩放因子的比值更加接近的对应关系作为两图像间的对应关系。

所述步骤4中，***部轮廓协方差矩阵提取的缩放因子，然后计算旋转因子和缩放因子计算出对比分析参考图像和待配准图像中选取窗口尺寸的比值，并对比分析待配准图像中选取不同窗口尺寸条件下两者的差值，如果两者差值达到最小，则获取该待配准窗口的尺寸作为最优窗口的尺寸；否则修正待配准图像的窗口尺寸后重新计算局部轮廓协方差矩阵提取的缩放因子，从而自适应校正分析窗口的尺寸。

本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的优点在于：

(1)本发明能够实现在大旋转缩放模式下遥感图像间的高精确自动配准，增强了遥感图像配准处理的鲁棒性，提高了遥感数据的处理效率，为实现遥感数据的自动处理奠定基础。

(2)本发明将角点特征、局部分割图像和局部边缘轮廓结合起来，共同作为局部特征来实现同名点的自动提取；

(3)本发明将局部轮廓协方差矩阵引入配准处理，实现参考图像与待配准图像间旋转缩放因子的自动提取，同时结合局部窗口自适应选取技术来减小旋转缩放因子对配准处理的影响，提高配准处理的鲁棒性。

附图说明

图1是基于图像特征的配准处理方法的流程图；

图2是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的步骤流程图；

图3是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的角点特征提取步骤流程图；

图4是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的边缘特征提取步骤流程图；

图5是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的最优窗口尺寸提取步骤流程图；

图6是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的角点特征匹配处理步骤流程图；

图7a是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的直角坐标系示意图；

图7b是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的本征坐标系示意图；

图8a是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的不同角点处的旋转系数的变化曲线；

图8b是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的不同角点处的缩放因子的变化曲线；

图9是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的最优窗口尺寸提取的分析结果图；

图10a是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的局部校正前线段间的相似性Dkl的变化曲线图；

图10b是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的局部校正前线段间的相似性Dkl的变化曲线图；

图11是本发明一种基于局部轮廓协方差矩阵的遥感图像配准方法的配准处理结果图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明主要针对大旋转缩放模式下遥感图像特征间难以建立映射关系这一难点，提出了一种基于局部轮廓协方差矩阵的遥感图像配准处理方法，该处理方法能够实现遥感图像间旋转缩放因子的自动提取，有效地增加了遥感图像特征间的匹配精度，增强了遥感图像配准处理的鲁棒性，在无人工参与条件下实现大旋转缩放模式下遥感图像间的精确配准。

  本方法首先提取输入图像中的角点特征，以角点特征为参考，提取角点附近的局部轮廓特征，并利用局部轮廓特征的协方差矩阵提取输入图像间的旋转缩放因子，减小旋转缩放因子对特征匹配的影响，提高同名点的提取精度；然后，利用检测到的同名点来进行仿射模型参数估计；最后，根据估算出来的变换模型参数进行图像重采样，获得配准后的遥感图像对。从整体上来看，该处理流程主要包含特征的选择与提取、旋转缩放因子的自动提取、最优窗口尺寸选取、角点特征匹配、仿射模型参数估计及图像重采样六个部分。

本方法一种基于局部轮廓协方差矩阵的遥感图像配准方法，如图2所示，包括如下步骤：

步骤一：特征的选择与提取。

特征的选择与提取是遥感图像配准处理的前提，其直接影响到配准处理的效果，甚至导致无法实现匹配处理。原则上讲，完整的边缘轮廓是一种很好的图像特征，其不仅具有很好的定位能力，而且能够很容易的建立起映射关系，因此，对于光学图像来说利用边缘轮廓能够很好的实现图像配准处理。然而对于雷达图像来说，强烈的斑点噪声使得很难提取完整的边缘轮廓，此时，若直接采用边缘特征进行图像配准无法获得较好的处理结果。为了有效的解决这一问题，在本发明中采用将角点特征、局部分割图像与局部边缘轮廓特征相结合，充分利用角点特征定位特性好和局部轮廓特征易于建立映射关系的特征，同时考虑到角点特征通常位于灰度突变的地方，因此可以有效降低局部轮廓特征的提取难度，提高配准处理的鲁棒性。

角点特征的提取：在本文中采用Harris角点检测算法实现图像角点检测。该算法通过统计移动搜索窗口内图像灰度在各方向上的梯度来判断是否存在角点。如果移动搜索窗口内图像灰度在任意方向上的灰度变化都很小，则搜索区域内不包含角点特征，如果移动搜索窗口内图像灰度在任意方向上的均出现明显的灰度变化，则搜索区域内的特征是角点。为有效判断图像灰度的变化，定义角响应函数CRF，其包含两种模式分别为Harris模式和Nobel模式：

$\mathrm{CRF}=\left\{\begin{array}{cc}\det M-k\·\mathrm{trace}\left(M\right)& \left(\mathrm{Harris}\right)\\ \frac{\mathrm{trace}\left(M\right)}{\det M}& \left(\mathrm{Nobel}\right)\end{array}\right.---\left(1\right)$

其中：M表示梯度矩阵，其表达式为：

$M=G\⊗\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]---\left(2\right)$

I_x表示x方向上的梯度；I_y表示y方向上的梯度；G为高斯模板；

表示卷积运算；det表示矩阵的行列式；trace表示矩阵的迹；k为常系数，一般取0.04。

判断每一点的CRF值是否是局部邻域内的最大值(Nobel角响应函数对应着最小值)。若是，则将该点设置为备选角点；若不是，则表明该点不是角点。最后利用门限处理来提取最终角点。

图3给出了角点提取流程图。

具体操作时，分别使用模版 $\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]$ 和模版 $\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]$ 在图像上沿方位向和距离向移动，计算每个位置上的梯度值，得到两个方向上的梯度图像，并计算每个位置对应的梯度乘积，这样实质上得到了三幅图，其像素值分别对应了I_x、I_y和I_xI_y；然后利用高斯模版分别与这三幅图像进行卷积处理，接着，根据公式(1)分别计算各点的角响应函数(CRF)值，并判断每一点的CRF值是否为其局部邻域内(如5×5)的最大值，若不是表明该点不是角点，若是，则将该点设置为备选角点，最后利用门限处理提取一定数量的最优点作为最后的结果。

步骤二：局部轮廓提取：在完成角点特征提取的基础上，进一步提取各角点的局部轮廓特征。根据角点的定义及角点特征的提取方法可知，角点特征通常位于灰度突变的地方，因此，在角点的局部区域存在着明显的灰度变化，因此，利用简单的图像分割方法就能够将局部图像进行而致划分，进而提取局部轮廓。在本发明中采用的类间方差自动阈值法来进行阈值分割。

如图4所示，根据图像的灰度范围{0，1，…，l-1}，依次选择阈值t将其分为两类：C₀：{0，1，…，t}和C₁：{t+1，t+2，…，l-1}。首先，***部图像的直方图对图像直方图归一化，各灰度值的概率分布为：

p_i＝n_i/N，p_i≥0， $\underset{i=0}{\overset{l-1}{\mathrm{\Σ}}}{p}_i=1---\left(3\right)$

其中：n_i表示第i个灰度值的像素点个数，N表示总像素点个数，p_i表示该灰度值的像素点所占概率。

此时，可以得到两类的类间方差σ_B ²为：

${\mathrm{\σ}}_B^2={\mathrm{\ω}}_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2={\mathrm{\ω}}_0{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_0)}^2---\left(4\right)$

其中：ω₀和ω₁分别表示C₀类和C₁类的出现概率，其计算公式为：

$\left\{\begin{array}{c}{\mathrm{\ω}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{p}_i\\ {\mathrm{\ω}}_1=\underset{i=t+1}{\overset{l-1}{\mathrm{\Σ}}}{p}_i\end{array}\right.---\left(5\right)$

μ₀、μ₁和μ_T分别表示C₀类、C₁类及全部灰度阶的均值，其计算公式为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}\frac{{\mathrm{ip}}_i}{{\mathrm{\ω}}_0}\\ {\mathrm{\μ}}_1=\underset{i=t+1}{\overset{l-1}{\mathrm{\Σ}}}\frac{{\mathrm{ip}}_i}{{\mathrm{\ω}}_1}\\ {\mathrm{\μ}}_T=\underset{i=0}{\overset{l-1}{\mathrm{\Σ}}}{\mathrm{ip}}_i\end{array}\right.---\left(6\right)$

设定不同的分割阈值，根据式(4)、式(5)和式(6)分别计算在各种阈值分割情况下两类图像的类间方差，并以类间方差最大为准则，可以获得最佳阈值t^*为：

$t^{*}\text{=}\underset{0\≤t\≤l-1}{\mathrm{Max}}{\mathrm{\σ}}_{B---\left(7\right)}^2$

并利用该阈值完成对图像的二值分割，最后通过对比分析二值图像中各点与其周围邻域中各点的像素值变化来提取保局部邻域的边缘特征，当存在像素值变化时表明存在边缘特征，当不存在像素值变化时，表明不存在边缘特征，最后，通过分析角点与边缘的相对位置关系提取离角点特征最近的边缘特征作为局部轮廓特征。

步骤三：旋转缩放因子的自动提取

当输入图像间存在较大的旋转缩放误差时，两图像间的局部相似性随旋转缩放误差的增大急剧降低，这使得特征匹配变得更加困难。而利用输入图像间的局部协方差矩阵能够有效地提取输入图像间旋转缩放因子，进而减小旋转缩放因子对特征匹配的影响，提高配准处理的鲁棒性。

任何一幅图像其自身都存在着这样一个本征坐标系，即通过这一坐标系，该图像的协方差矩阵为对角矩阵，而协方差矩阵的特征值则反映了图像中的能量分布。以星形目标为例来说明本征坐标系，如图7a所示，在直角坐标系中存在一星形目标，此时，对于该星形目标的本征坐标系如图7b所示，在本征坐标系内图像关于轴线对称，其协方差矩阵为对角矩阵，特征值反映了图像能量的分布。通过对比输入图像角点处局部轮廓的本征坐标系的关系可以提取配准图像和待配准图像间的旋转因子，比较协方差矩阵的特征值可以得到两图像之间的缩放因子。下面详细推导旋转缩放因子的提取方法。

假设I_i和T_j分别为参考图像与待配准图像中的对应角点，则参考图像角点I_i的局部邻域I₀和待配准图像角点T_j的局部邻域T₀间存在仿射变换f，使得两图像间满足：

T₀＝f(I₀)＝AI₀(8)

其中：A为2×2的仿射变换矩阵。

在提取局部尺度旋转因子之前先依据局部灰度分布对局部图像采取阈值分割处理，将灰度值大于阈值的点设置为1，小于阈值的点设置为0，形成特征图像，记作I和T。考虑到虽然多模式遥感图像的灰度差异很大，但各图像的基本结构特征是一致的，因此，其特征图像之间仍然满足仿射变换关系。

采用基于局部轮廓协方差矩阵的旋转缩放因子自动提取方法。利用阈值分割所获得的二值图像，分别计算参考图像I与待配准图像T的协方差矩阵∑_I和∑_T，两协方差矩阵间满足：分别计算参考图像I与待配准图像T的协方差矩阵∑_I和∑_T，两协方差矩阵间满足：

∑_T＝A∑_IA′(9)

其中： ${\mathrm{\Σ}}_I=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\stackrel{\→}{r}}_{\mathrm{ij}}{\stackrel{\→}{r}}_{\mathrm{ij}}^{\′};$ ${\mathrm{\Σ}}_T=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\stackrel{\→}{r}}_{\mathrm{tj}}{\stackrel{\→}{r}}_{\mathrm{tj}}^{\′};$ N为不为O的像素点个数； ${\stackrel{\→}{r}}_j=\left[\begin{array}{c}x\\ y\end{array}\right]$ 为每个点的位置矢量，上标′表示转置。

由于斜方差矩阵∑_I和∑_T均为正规矩阵，存在正交矩阵U₁、U₂和对角矩阵Λ₁、Λ₂使得：

${\mathrm{\Σ}}_I=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\stackrel{\→}{r}}_{\mathrm{ij}}{\stackrel{\→}{r}}_{\mathrm{ij}}^{\′}=U_1{\mathrm{\Λ}}_1{U}_1^{\′}$ ${\mathrm{\Σ}}_T=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\stackrel{\→}{r}}_{\mathrm{tj}}{\stackrel{\→}{r}}_{\mathrm{tj}}^{\′}=U_2{\mathrm{\Λ}}_2{U}_2^{\′}---\left(10\right)$

其中： $U_1=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_1& \sin {\mathrm{\θ}}_1\\ -\sin {\mathrm{\θ}}_1& \cos {\mathrm{\θ}}_1\end{array}\right];$ $U_2=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_2& \sin {\mathrm{\θ}}_2\\ -\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\end{array}\right];$

${\mathrm{\Λ}}_1=\left[\begin{array}{cc}{\mathrm{\λ}}_{11}& 0\\ 0& {\mathrm{\λ}}_{12}\end{array}\right],$ λ₁₁＞λ₁₂； ${\mathrm{\Λ}}_2=\left[\begin{array}{cc}{\mathrm{\λ}}_{21}& 0\\ 0& {\mathrm{\λ}}_{22}\end{array}\right],$ λ₂₁＞λ₂₂。

将上式代入式(9)有：

U₂Λ₂U′₂＝AU₁Λ₁U′₁A′(11)

令 $B_2=U_2{\mathrm{\Λ}}_2^{0.5},$ $B_1=U_1{\mathrm{\Λ}}_1^{0.5}$ 有：

$A=B_2{B}_1^{-1}---\left(12\right)$

这样，依据参考图像与待配准图像的局部结构分布即可提取图像间的旋转因子

和缩放因子λ。

λ_i＝(λ_2i/λ_1i)(13)

在具体实现旋转缩放因子提取时，需要先确定特征图像间的对应关系，确定图像间的对应关系时，对比分析参考图像中亮、暗图像对应待配准图像中的亮、暗图像及参考图像中亮、暗图像对应待配准图像中的暗、亮图像，这两种情况下旋转角度差和缩放因子的比例，选取旋转角度差更小，两轴旋转缩放因子λ₁和λ₂的比值更加接近的对应关系作为两图像间的对应关系。

如图8a所示，为本实例中输入图像的不同角点处估算的旋转角度的提取结果，图8b给出了各点的缩放因子提取结果。从图可见，提取的旋转因子的变化范围在1.2度以内，缩放因子的变化范围在0.14以内，它们各自都具有很高的一致性，反应了输入图像间的旋转缩放关系。

步骤四：最优窗口尺寸选取

当两幅图像间存在较大的旋转和尺度误差时，常用的局部匹配准则如相关系数、交叉相关系数、交互信息等的测试性能急剧下降，造成测试性能下降的原因主要有两个方面：

1：由于旋转因子和缩放因子造成的测试性能下降；

2、由于缩放因子的影响，在相同的窗口大小下，图像内容的差异造成测试性能下降。因此，当两幅图像间存在较大的旋转和缩放误差时，需要根据图像的内容，提取两幅图像间的旋转因子和缩放因子，自适应的调整分析窗口的大小，消除旋转缩放误差的影响。

自适应调整分析窗口尺寸的基本原则如下：当所选取的窗口尺寸为最优窗口尺寸时，两窗口内的图像内容一致，此时两图像局部轮廓协方差矩阵特征值的比值反映了两图像之间的缩放因子，同时两图像中分析窗口的比值同样反映了两图像间的缩放因子，理论上将，当达到最优窗口尺寸时，两比值是一致的，因此，可以通过对比分析参考图像和待配准图像中选取窗口尺寸的比值p和输入图像对局部轮廓协方差矩阵提取的缩放因子λ来选取分析窗口尺寸。如图5所示，具体分析方法如下：

***部轮廓协方差矩阵提取的缩放因子λ，然后计算旋转因子和缩放因子计算出对比分析参考图像和待配准图像中选取窗口尺寸的比值p，并对比分析待配准图像中选取不同窗口尺寸条件下λ和p的差值，如果λ和p的差值达到最小，则获取该待配准窗口的尺寸作为最优窗口的尺寸；否则修正待配准图像的窗口尺寸后重新计算局部轮廓协方差矩阵提取的缩放因子λ，从而自适应校正分析窗口的尺寸。

如果参考图像中局部窗口内的局部图像与待配准图像中局部窗口内的局部图像纹理特征相同，此时，参考图像和待配准图像中选取窗口尺寸的比值也正好反映了两图像间的缩放因子，因此，根据参考图像和待配准图像中选取窗口尺寸的比值所计算出的缩放因子与根据式(13)所提取的缩放因子将会一致。

如果参考图像中局部窗口内的局部图像与待配准图像中局部窗口内的局部图像纹理特征不相同，根据参考图像和待配准图像中选取窗口尺寸的比值所获取的缩放因子和根据式(13)所提取的缩放因子都不能正确的提取局部图像间的缩放因子，此时，两结果之间将会存在一定的偏差。

因此，可以通过比较两种情况下所提取的缩放因子的差值来判断是否获得最优窗口尺寸，即当差值达到最小时达到最优窗口尺寸。图9给出了最优窗口尺寸提取的分析结果，从图中可见有选出的最优窗口在23～25之间波动，通过最终的配准处理可知两图像间的理论最优窗口为23，通过比较可见提取的最优窗口尺寸具有高度的一致性，很好的反应了输入图像间的旋转缩放关系。

步骤五：角点特征匹配

利用局部边缘轮廓特征来实现角点特征的匹配。首先利用提取的旋转缩放因子和最优窗口尺寸对步骤一中提取的边缘轮廓进行重采样，减小旋转缩放因子对特征匹配的影响。

然后，根据各点间的空间位置关系对重采样后的局部边缘特征进行Freeman编码，经过曲线编码后两条长度为N的线段间的相似性度量公式如下式所示：

$D_{\mathrm{kl}}=\frac{1}{N}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\cos \frac{\mathrm{\π}}{4}(a_{k+j}-b_{l+j})---\left(14\right)$

其中：a_k+j和b_l+j分别表示两条线段进行Freeman编码后所得到的码值；D_kl是一个0到1之间的值，两条线段越匹配，其值越趋近于1，两条直线差别越大，其值约趋近于0，两条线段完全匹配时，其值为1。

统计两线段的相似性，当两线段间的相似性D_kl满足如下条件时，认为这对角点为同名点。

$D_{\mathrm{kl}}=\underset{k\∈A}{\max }\left(D_{\mathrm{kl}}\right)$ $D_{\mathrm{kl}}=\underset{l\∈B}{\max }\left(D_{\mathrm{kl}}\right)$ D_kl＞D_th(15)

其中：A表示参考图像中提取的角点集；B表示待配准图像中提取的角点集；D_th表示曲线相似度的阈值。

至此，通过上述步骤能够实现遥感图像间角点特征的粗匹配，建立了点与点之间的映射关系，然而由于噪声及图像灰度间的差异，不可避免的存在错误匹配点对，为了提高配准处理的精度，有必要在粗匹配的基础上进一步剔除错误的匹配点对。

为此，可进一步利用点集间的相对位置关系和点集间的相互支持特性剔除错误匹配点。首先，将上面所述初始匹配点对按照相似性测试的大小进行排序。考虑到4对控制点能够很好反应控制点件的空间结构，并将所有控制点进行按照排列顺序组合成由4对控制点构成的控制点组，顺序检索所有的控制点组。假设它们之间满足仿射变换模型，当这4对控制点对在空间上不存在共线现象时，根据这4对控制点对估算的仿射矩阵与实际仿射矩阵间的误差很小，利用估算的仿射矩阵对待配准图像中的控制点进行仿射变换，并以4对控制点对间的平均位置误差为指标来衡量是否为正确控制点对，当控制点对间的平均位置误差小于设定的阈值时，可以认为找到正确控制点对，停止搜索。当控制点的平均位置误差大于设定的阈值或者4对控制点对在空间上存在共线现象时，则认为是误匹配点，继续搜索，直到初始提取的控制点全部分析完毕。

在确定4对正确控制点对的基础上，测试其余的控制点对是否为正确的控制点对。顺序将待测试的控制点对引入正确控制点对，根据正确控制点对估算出仿射矩阵，重新计算控制点对间的平均位置误差，当平均位置误差大于设定阈值时，认为新加入的控制点对为错误控制点对，将其从初始匹配点对中删除。否则，认为新引入的控制点对为正确控制点对，继续搜索其余控制点对，提取所有正确控制点对。图6给出角点特征匹配处理的流程图。

如图10a、10b所示，图10a为未提取旋转缩放因子及最优窗口尺寸下的相似性D_kl的变化曲线，图10b为提取旋转缩放因子及最优窗口尺寸后的相似性D_kl的变化曲线。从结果中可见，由于能够消除局部图像间旋转缩放因子的影响，大大提高了角点特征匹配的鲁棒性。

步骤六：模型参数估算及图像重采样

当完成特征点匹配后，可以利用最小二乘法估算仿射模型中的参数。具体估计方法如下：

假设输入图像间共提取了步骤四中的N对控制点，分别为{((x_r1，y_r1)，(x_w1，y_w1))，((x_r2，y_r2)，(x_w2，y_w2))，…，((x_rN，y_rN)，(x_wN，y_wN))}，考虑到图像间满足仿射关系，此时，任意控制点对间满足如下关系式：

$\left[\begin{array}{c}{x}_{\mathrm{wi}}\\ y_{\mathrm{wi}}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ri}}\\ y_{\mathrm{ri}}\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]=\left[\begin{array}{cccccc}1& 0& x_{\mathrm{ri}}& y_{\mathrm{ri}}& 0& 0\\ 0& 1& 0& 0& x_{\mathrm{ri}}& y_{\mathrm{ri}}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]---\left(16\right)$

将提取得所有控制点均代入方程，可以得到如下方程组：

$Y=\left[\begin{array}{c}{x}_{w1}\\ y_{\mathrm{wl}}\\ \·\\ \·\\ \·\\ x_{\mathrm{wN}}\\ y_{\mathrm{wN}}\end{array}\right]=\left[\begin{array}{cccccc}1& 0& x_{r1}& y_{r1}& 0& 0\\ 0& 1& 0& 0& x_{r1}& y_{r1}\\ & & & \·& & \\ & & & \·& & \\ & & & \·& & \\ 1& 0& x_{\mathrm{rN}}& y_{\mathrm{rN}}& 0& 0\\ 0& 1& 0& 0& x_{\mathrm{rN}}& y_{\mathrm{rN}}\end{array}\right]=\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]=H\·\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]---\left(17\right)$

将各对控制点的位置坐标代入方程组后，可以得到仿射参数的估算值如下式所示：

$\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]={(H^T\·H)}^{-1}\·H^T\·Y---\left(18\right)$

其中：上标T表示转置，上标-1表示求逆。

在确定变换参数后，可以依据仿射参数对输入图像对进行重采样处理，即可获得相互匹配的图像对，在本处理算法中采用双线性内插法实现图像的重采样。

首先根据参考图中各点的位置，计算出式(17)中的系数矩阵H及矩阵H^T·H，接着将矩阵H^T.H分解为非奇异对角矩阵D和非奇异下三角行矩阵Λ，

$D=\left[\begin{array}{cccccc}{d}_1& 0& 0& 0& 0& 0\\ 0& d_2& 0& 0& 0& 0\\ 0& 0& d_3& 0& 0& 0\\ 0& 0& 0& d_4& 0& 0\\ 0& 0& 0& 0& d_5& 0\\ 0& 0& 0& 0& 0& d_6\end{array}\right],$ $\mathrm{\Λ}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ {\mathrm{\λ}}_{21}& 1& 0& 0& 0& 0\\ {\mathrm{\λ}}_{31}& {\mathrm{\λ}}_{32}& 1& 0& 0& 0\\ {\mathrm{\λ}}_{41}& {\mathrm{\λ}}_{42}& {\mathrm{\λ}}_{43}& 1& 0& 0\\ {\mathrm{\λ}}_{51}& {\mathrm{\λ}}_{52}& {\mathrm{\λ}}_{53}& {\mathrm{\λ}}_{54}& 1& 0\\ {\mathrm{\λ}}_{61}& {\mathrm{\λ}}_{62}& {\mathrm{\λ}}_{63}& {\mathrm{\λ}}_{64}& {\mathrm{\λ}}_{65}& 1\end{array}\right]---\left(19\right)$

其中： $x_{\mathrm{ij}}=h_{\mathrm{ij}}-\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}{x}_{\mathrm{ik}}{\mathrm{\λ}}_{\mathrm{jk}},$ ${\mathrm{\λ}}_{\mathrm{ij}}=\frac{x_{\mathrm{ij}}}{d_j},$ $d_i=h_{\mathrm{ii}}-\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}{x}_{\mathrm{ik}}{\mathrm{\λ}}_{\mathrm{jk}}\·$

分别求取矩阵D和矩阵Λ的逆矩阵D^-1和Λ^-1。

$D^{-1}=\left[\begin{array}{cccccc}\frac{1}{d_1}& 0& 0& 0& 0& 0\\ 0& \frac{1}{d_2}& 0& 0& 0& 0\\ 0& 0& \frac{1}{d_3}& 0& 0& 0\\ 0& 0& 0& \frac{1}{d_4}& 0& 0\\ 0& 0& 0& 0& \frac{1}{d_5}& 0\\ 0& 0& 0& 0& 0& \frac{1}{d_6}\end{array}\right],$ $\mathrm{\Λ}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ {\mathrm{\ρ}}_{21}& 1& 0& 0& 0& 0\\ {\mathrm{\ρ}}_{31}& {\mathrm{\ρ}}_{32}& 1& 0& 0& 0\\ {\mathrm{\ρ}}_{41}& {\mathrm{\ρ}}_{42}& {\mathrm{\ρ}}_{43}& 1& 0& 0\\ {\mathrm{\ρ}}_{51}& {\mathrm{\ρ}}_{52}& {\mathrm{\ρ}}_{53}& {\mathrm{\ρ}}_{54}& 1& 0\\ {\mathrm{\ρ}}_{61}& {\mathrm{\ρ}}_{62}& {\mathrm{\ρ}}_{63}& {\mathrm{\ρ}}_{64}& {\mathrm{\ρ}}_{65}& 1\end{array}\right]---\left(20\right)$

其中： $\left\{\begin{array}{c}{\mathrm{\ρ}}_{\mathrm{ij}}=-\underset{k=j}{\overset{i-1}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{ik}}{\mathrm{\ρ}}_{\mathrm{kj}}\\ {\mathrm{\ρ}}_{\mathrm{ii}}=1\end{array}\right.(j=\mathrm{1,2},\·\·\·,n-1;i=j+1,j+2,\·\·\·,n)$

这样，利用逆矩阵D^-1和Λ^-1可以求取矩阵H^T·H的逆矩阵，结合矩阵H及待配准图像中控制点对的位置，可以计算出仿射变换参数。

表1给出了配准处理结果，图11给出了控制点位置的均方根误差的变化曲线，其中参考图像X坐标和Y坐标表示参考图像中控制点的位置，待配准图像X坐标和Y坐标表示待配准图像中控制点的位置，重采样后X坐标和Y坐标表示利用提取的房舍参数对待配准图像中的控制点进行变换处理后得到的位置坐标，位置的均方根误差表示重采样后控制点的位置与参考图像中控制点位置间偏离的距离，从表1和图11可以看出，在存在25度旋转误差、1.4倍的缩放因子条件下，10对同名点经过校正处理后位置的均方根最大为0.94个像素单元，能够实现配准误差优于1个像素单元。

表1配准处理结果

序号	参考图像X坐标	参考图像Y坐标	待配准图像X坐标	待配准图像Y坐标	重采样后X坐标	重采样后Y坐标	位置的均方根误差
								0	652	856	943	564	651.78	855.95	0.23
1	620	679	989	390	619.96	678.87	0.14
								2	498	655	889	317	498.54	654.95	0.54
3	188	520	665	64	188.78	520.30	0.84
								4	298	758	664	325	298.11	757.20	0.81
5	45	874	386	324	45.88	873.67	0.94
								6	316	782	670	355	316.22	781.85	0.27
7	104	549	576	54	103.94	548.82	0.19
								8	442	442	928	100	442.21	441.87	0.25
9	447	591	869	237	446.64	590.91	0.37

Claims (4)

1、一种基于局部轮廓协方差矩阵的遥感图像配准方法，其特征在于，包括如下步骤：

步骤一：从各图像通道中同时读入各图像，采用Harris角点检测算法实现图像角点检测，通过统计移动搜索窗口内各点的角响应函数值来判断是否存在角点；

如果一点的角响应函数CRF值是局部邻域内的最大值，则将该点设置为备选角点；否则表明该点不是角点；在判断完图像中所有点后，利用门限处理来提取最终角点；

步骤二：采用类间方差自动阈值法对参考图像及待配准图像进行局部阈值分割；利用提取得最佳阈值对局部图像进行分割，获得局部二值图像，并依据提取的局部二值图像获得局部边缘轮廓；

步骤三：旋转缩放因子的自动提取；采用基于局部轮廓协方差矩阵的旋转缩放因子自动提取方法；利用阈值分割所获得的二值图像，分别计算参考图像与待配准图像的协方差矩阵；依据参考图像与待配准图像的局部结构分布即可提取图像间的旋转因子和缩放因子；

步骤四：通过对比分析参考图像和待配准图像中选取窗口尺寸的比值和利用输入图像提取的缩放因子来选取分析窗口尺寸，提取最优窗口；

步骤五：利用提取的旋转缩放因子和最优窗口尺寸对步骤一中提取的边缘轮廓进行重采样，减小旋转缩放因子对特征匹配的影响；并对重采样后的边缘特征进行Freeman编码，并通过链码相关技术统计两线段间的相似性；利用点集间的相对位置关系和点集间的相互支持特性剔除错误匹配点；

步骤六：当完成特征点匹配后，利用最小二乘法估算仿射模型中的参数；在确定变换参数后，依据仿射参数对输入图像对进行重采样处理，可获得相互匹配的图像对。

2、根据权利要求1所述一种基于局部轮廓协方差矩阵的遥感图像配准方法，其特征在于：

所述步骤三中，针对局部分割图像展开旋转缩放因子的自动提取；通过阈值分割将图像分割为亮、暗两类图像，然后分别统计两部分的协方差矩阵及旋转缩放因子，确定图像间的对应关系，并提取相应的旋转缩放因子。

3、根据权利要求1所述一种基于局部轮廓协方差矩阵的遥感图像配准方法，其特征在于：

所述步骤三中，确定图像间的对应关系时，对比分析参考图像中亮、暗图像对应待配准图像中的亮、暗图像及参考图像中亮、暗图像对应待配准图像中的暗、亮图像这两种对应关系，选取旋转角度差和缩放因子的比例，选取旋转角度差更小，两轴旋转缩放因子的比值更加接近的对应关系作为两图像间的对应关系。

4、根据权利要求1所述一种基于局部轮廓协方差矩阵的遥感图像配准方法，其特征在于：

所述步骤四中，***部轮廓协方差矩阵提取的缩放因子，然后计算旋转因子和缩放因子计算出对比分析参考图像和待配准图像中选取窗口尺寸的比值，并对比分析待配准图像中选取不同窗口尺寸条件下两者的差值，如果两者差值达到最小，则获取该待配准窗口的尺寸作为最优窗口的尺寸；否则修正待配准图像的窗口尺寸后重新计算局部轮廓协方差矩阵提取的缩放因子，从而自适应校正分析窗口的尺寸。

Images