CN102693542A - 一种影像特征匹配方法 - Google Patents

Abstract

本发明涉及一种影像特征匹配方法，对获取的CCD影像进行预处理；利用SURF算子对预处理后的CCD影像提取特征点，利用准核线约束条件和最小欧氏距离条件进行影像匹配，获取同名点信息；根据获得的同名点信息建立CCD影像间的仿射变形关系；利用Harris角点提取算子提取参考影像上特征点，并利用仿射变换投影至搜索影像上得到待匹配点；在待匹配点周围一邻域内，利用积分影像计算特征点与邻域内各点之间的相关系数，取极值点作为其同名点；将两次匹配的结果综合作为最终同名点信息；本发明的方法可以利用在深空环境下获得的深空星体表面影像进行影像匹配而获得高精度的CCD影像同名点信息，实现特征匹配。

Description

一种影像特征匹配方法

技术领域

本发明属于摄影测量与遥感技术中的行星摄影测量技术领域，涉及一种适合于星球表面DEM生成的CCD影像特征匹配方法。

背景技术

影像匹配技术，是建立影像间的对应关系。基于特征的影像匹配方法是比较影像中提取的特征之间的相似性以实现同名点的提取的方法，影像特征主要包括点、线、面特征。这类方法通常包括三个步骤：（1）特征提取：利用兴趣算子提取影像上的特征；（2）候选特征的确定：通过对提取的特征的属性进行比较获得待配准的候选特征；（3）特征对应：对提取的特征进行对应获得匹配结果。该方法的优点受噪声、亮度、变化及各种几何变形影响小，如果再加上合适的约束条件，匹配速度和效率会有明显提高，而且在区域网平差、数字高程模型（DigitalElevation Model，DEM）的自动生成等后续处理中匹配的结果可以再次利用。

由于深空星体表面的特征与地球特征相比，没有大量的地物，因此对于获取的影像上的特征较少，在利用现有的常规摄影测量特征点提取算子提取影像上的特征点时，获得的点数较少，不利于后续月球表面DEM的生成，因此，本发明针对这一问题，研究了我国的探月卫星嫦娥一号获得的CCD影像进行影像匹配而获得较多特征点的方法。

发明内容

本发明的目的是提供一种影像匹配方法，以解决现有的影像匹配方法针对深空星体影像存在的提取特征点较少不利于星体的DEM生成问题。

为实现上述目的，本发明的影像特征匹配方法步骤如下：

（1）对获取的CCD影像进行预处理；

（2）利用SURF算子对预处理后的CCD影像提取特征点，利用准核线约束条件和最小欧氏距离条件进行影像匹配，获取同名点信息；

（3）根据获得的同名点信息建立CCD影像间的仿射变形关系；

（4）利用Harris角点提取算子提取参考影像上特征点，并利用仿射变换投影至搜索影像上得到待匹配点；

(5)在待匹配点周围一邻域内，利用积分影像计算特征点与邻域内各点之间的相关系数，取极值点作为其同名点；

(6)将步骤(2)和步骤(5)匹配的结果共同作为最终同名点信息。

进一步的，所述步骤(1)中使用Wallis滤波器对CCD影像进行预处理。

进一步的，所述步骤(2)利用SURF算子提取特征点是计算原始影像的积分影像，然后建立尺度空间，使用Hessian矩阵的行列式同时来确定每个尺度上特征点的位置和尺度，当Hessian矩阵的行列式局部取最大时，检测出的像点做为特征点。

进一步的，所述步骤(2)利用准核线约束条件和最小欧氏距离条件进行影像匹配的过程如下：

(21)根据参考影像构像模型，对于参考影像上的每个特征点建立投影光线；

(22)在投影光线上根据预先设定的月面点的近似高程Z₀以及高度误差ΔZ，结合搜索影像的构像模型在搜索影像上确定准核线；

(23)以准核线为中心，在其周围按照参与y视差扩展设置一邻域；

(24)计算邻域内部搜索影像上的特征点与参考影像上特征点之间的欧氏距离和正则化交叉相关系数NCC；

(25)对欧氏距离和NCC进行对比，如果两者均达到极值，则该点与参考影像上的特征点位一对同名像点。

进一步的，所述步骤(24)中相关系数NCC的计算中，设I₀为参考影像，p₀为I₀上的一个像点，W为点p₀处的相关窗口，窗口大小为M×N，NCC可以表示为：

$\mathrm{NCC}(p_0,Z)=\frac{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)\×(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}{\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)}\·\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}}$

其中： ${\stackrel{\‾}{I}}_0=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_0\left(p_0\right)$ 和 ${\stackrel{\‾}{I}}_1=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_1\left(p_1\left(Z\right)\right)$ 为窗口内的灰度均值，I₀(p₀)与I₁(p₁(Z))分别为参考影像点p₀处与搜索影像点p₁处的灰度值。

进一步的，所述步骤（2）完成后，在步骤（3）开始之前对同名点进行粗差剔除，并进行最小二乘优化。

进一步的，所述对同名点进行粗差剔除是应用RANSAC剔除粗差匹配点，过程如下：

（221）建立两幅影像上的同名点关系；

（222）任意选取3对同名点，建立两幅影像之间的仿射变换模型；

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{a}_0\\ b_0\end{array}\right]+\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]---\left(8\right)$

其中，(x，y)和(x′,y′)分布为同名点在两幅影像上坐标，(a_i,b_i)(i=0,1,2)为仿射变换参数。

（223）计算其余同名点对的匹配误差，根据事先设定的阈值，判断是否为内点，并统计内点个数；

（224）通过选取不同的同名点对，重复以上（222）—（223）步骤，直至所有点都经过一次计算；

（225）选取内点数量最多的内点集作为最佳匹配点集。

进一步的，所述步骤（4）中利用Harris角点提取算子提取参考影像上特征点的过程如下：

（41）计算图像在x和y方向的梯度g_x,g_y；

（42）计算自相关函数矩阵 $M=G\left(s\right)*\left[\begin{array}{cc}{g}_x^2& g_x{g}_y\\ g_x{g}_y& g_y^2\end{array}\right],$ 其中，G(s)为高斯滤波器；

（43）计算每一个像点的兴趣值R=det(M)-k(tr(M)²)，det(M)和tr(M)分别表示矩阵的行列式值和迹，k为常数；

（44）在局部范围内利用兴趣值进行局部非最大值抑制，寻找极值点作为特征点。

进一步的，所述步骤（5）中利用积分影像计算特征点与邻域内各点之间的相关系数是以待匹配点为中心，定义一窗口，采用积分图像计算窗口内部的各点的相关系数，

$\mathrm{\ρ}(r,c)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}\·g_{i+r,j+c}^{\′})-\frac{1}{m\×n}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}}{\sqrt{(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^2-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2)(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′2}-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′}\right)}^2)}}$

其中，g_i，j分别表示参考影像上行列分别为i，j的点的灰度值，g′_i，j分别表示搜索影像上行列分别为i，j的点的灰度值，m×n为在影像上定义的窗口大小。

本发明的影像特征匹配方法，可以保证在深空环境下，进行影像匹配而获得高精度的CCD影像同名点信息。本方法中首先利用SURF算子提取特征点进行初次匹配以获得待匹配CCD影像的同名点，然后对该结果进行多尺度约束而获得较可靠的同名点信息，这些点的匹配精度可以到达优于1个像素；但其获取的特征点分布不是十分均匀，不利于星体表面的DEM生成，因此本方法利用初次匹配提取的同名点建立的影像间的仿射变形关系，在此约束条件下利用Harris算子提取的特征点进行二次匹配而获得同名点信息；两次匹配的结果共同构成最终的匹配同名点信息，其分布也较为均匀，有利于DEM生成。

附图说明

图1本发明实施例CE-1线阵CCD影像特征匹配方法流程图；

图2是积分图像原理图；

图3是SURF尺度空间的创建方式示意图；

图4是残余y视差图；

图5是山脉地区匹配后同名点信息的效果图；

图6是月海地区匹配后同名点信息的效果图；

图7（a）是图6原始的参考影像，（b）为经过预处理后的影像；

图8（a）是图6原始的搜索影像，（b）为经过预处理后的影像；

图9是在图7上利用SURF提取特征点后经过初次匹配的结果；

图10是利用初次匹配的结果提取Harris角点进行匹配的结果。

具体实施方式

如图1所示，以嫦娥一号CCD影像特征匹配方法来进行详细说明。

1.嫦娥一号的载荷

嫦娥一号卫星上搭载光学传感器是一个大面阵CCD探测器，其获取的三线阵影像是分别读取面阵探测器上垂直于飞行方向上的第11行、第512行和第1013行，分别作为前视、正视和后视的CCD阵列。当卫星飞行时，这三个CCD阵列可以连续的扫描月面产生三条相互重叠的航带图像，即前视、后视和下视影像。嫦娥一号卫星上搭载的激光高度计的任务是激光测距，它有两个科学目标：(1)测量卫星到月球表面星下点的距离；(2)提供月球表面数字高程模型的径向高度数据。

2.影像预处理

为了在提取点特征和进行影像匹配时提高点特征的数量及匹配结果的可靠性、精度，因此首先对嫦娥一号CE-1获取的CCD影像进行预处理，影像预处理采用Wallis滤波。

Wallis滤波器可以在增强原始影像的反差的同时压制噪声，而且它可以大大增强影像中不同尺度的影像纹理模式。Wallis滤波的原理是将局部影像的灰度均值和方差映射到给定的灰度均值和方差，从而使反差大的影像区域反差减小，反差小的影像区域反差增大。Wallis滤波器在对影像处理过程中引入了平滑算子，所以能够抑制一部分噪声，提高影像的信噪比。

3.基于SURF算子的特征点提取

SURF（Speeded Up Robust Features）算子与摄影测量中常用的角点提取算子相比较，在各方面都有具有较好的稳定性，而且它们主要提取圆状特征点，这与CE-1线阵CCD月球影像的基本特征非常相似。

SURF使用Hessian矩阵来检测特征点，用箱式滤波器（Box Filters）来代替二阶高斯滤波，用积分图像（Integral Image）来加速卷积以提高计算速度，并减少了局部影像特征描述符的维数，更适用于实时处理。提取SURF特征是通过计算原始影像的积分影像，然后建立尺度空间并在每个尺度上检测特征点。

（1）、积分图像

如图2所示，积分图像定义为：积分图像上任一点(x，y)的值I_∑(x，y)可以表示成原始影像I上以影像左上角为原点，以(x，y)为右下角所组成的矩形区域内所有点的灰度值之和，即：

$I_{\mathrm{\Σ}}(x,y)=\underset{i=0}{\overset{i\≤x}{\mathrm{\Σ}}}\underset{j=0}{\overset{j\≤y}{\mathrm{\Σ}}}I(x,y)---\left(1\right)$

（2）、构建SURF尺度空间

为了在不同的尺度上寻找局部极值点，通常使用影像金字塔构建尺度空间，如图3所示。SURF在生成尺度空间金字塔时，是在保持原始图像不变的前提下，采用不同大小的滤波器得到的。由于减少了类似SIFT算子的亚采样过程，因此尺度空间的构建效率与影像的大小无关，这样可以提高计算速度，缩短计算时间。

（3）、Hessian矩阵

SURF算子的特征点检测是通过Hessian矩阵进行的，使用Hessian矩阵的行列式同时来确定特征点的位置和尺度，当Hessian矩阵的行列式局部取最大时，检测出的像点做为特征点，而其实际是一个小区域，即斑状结构。

影像上任一点(x，y)在尺度σ时的Hessian矩阵H(x，y，σ)定义为：

$H(x,y,\mathrm{\σ})=\left[\begin{array}{cc}{L}_{\mathrm{xx}}(x,y,\mathrm{\σ})& L_{\mathrm{xy}}(x,y,\mathrm{\σ})\\ L_{\mathrm{xx}}(x,y,\mathrm{\σ})& L_{\mathrm{yy}}(x,y,\mathrm{\σ})\end{array}\right]---\left(2\right)$

式中：L_xx(x，y，σ)、L_xy(x，y，σ)、L_yy(x，y，σ)分别是影像在点(x，y)处与高斯二阶导数

的卷积。其中g(x，y，σ)表示影像上任一点(x，y)尺度可变高斯函数，即： $g(x,y,\mathrm{\σ})=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)}$

求解高斯函数二阶导数时，SURF使用方框滤波器来近似高斯二阶导数，其原因在于使用积分图像可以迅速计算出均值滤波的图像卷积。将方框滤波器与图像卷积的结果表示为D_xx、D_xy、D_yy，那么Hessian矩阵的行列式可以表示为

det(H)=D_xxD_yy-(ωD_xy)²    （3）

其中ω是用来平衡高斯核的近似时产生的能量差异，可以利用一个常数代替，通常取ω=0.9。

（4）、局部最大值点定位

在建立尺度空间后，通过计算3×3×3局部区域内各点Hessian矩阵的行列式值，将比邻近的26个点值都大的点作为特征点。

（5）、特征点主方向的确定

检测出特征点后，在其周围选取一个以6σ为半径的圆形区域，在此区域内使用积分图像计算Harr小波在x，y轴方向上的响应，并赋以高斯权重系数，然后将圆形区域等分成6份，分别计算以圆心角为π/3的扇形窗口内x，y方向小波的响应矢量和，选择最大值所在的矢量方向作为该特征点的主方向。

（6）、特征点描述子特征向量生成

在确定特征点的主方向后，选取特征点周围的一个正方形区域，将该区域划分成4×4的子区域，在每一个子区域内，计算5×5内部点的相对于主方向的水平、垂直方向的Haar小波响应值，记做dx，dy，并赋以相应的权值，然后将每个子区域的响应、响应的绝对值相加构成四维矢量V=[∑dx，∑dy，∑|dx|,∑|dy|]。这样就形成了一个64维的向量，此向量就是描述该特征点的描述子特征向量。

SURF描述了特征点局部区域的梯度信息的空间分布情况，由于它将每个子块的梯度信息作为一个整体，因此其对噪声不敏感。

4.基于多尺度约束的影像匹配原理

在利用SURF算子提取特征点后，本文提出了结合准核线约束和最小欧氏距离进行匹配的方法。该方法的流程可以表示为：

（1）利用SURF算子提取CE-1线阵CCD影像上的特征点；

（2）根据参考影像构像模型，对于参考影像上的每个特征点建立投影光线；

CE-1传感器的成像模型是描述CE-1线阵CCD影像像点和对应月面点空间位置关系的数学表达式，是进行月球影像测图、月面点坐标确定及月面DEM生成的基础。常用的成像模型主要分为两类，一类是严格物理模型，一类是通用模型。本方法采用通用模型——有理函数模型描述成像状态。

有理函数模型（Rational Function Model，简称RFM）是将CCD影像上像点坐标(r，c)表示为以相应月面点空间坐标(X，Y，Z)为自变量的多项式的比值，即：

$\left\{\begin{array}{c}{r}_n=\frac{p_1(X_n,Y_n,Z_n)}{p_2(X_n,Y_n,Z_n)}\\ c_n=\frac{p_3(X_n,Y_n,Z_n)}{p_4(X_n,Y_n,Z_n)}\end{array}\right.---\left(4\right)$

其中，（r_n，c_n）和（X_n，Y_n，Z_n）分别表示像点坐标(r，c)和月面点坐标(X，Y，Z)经平移和缩放后的标准化坐标。RFM采用标准化坐标的目的是减少计算过程中由于数据数量级差别过大引入的舍入误差。p_i(X_n,Y_n,Z_n)（i＝1，2，3，4）为多项式，通常多项式最高次数不超过3次，多项式的系数称为有理函数系数。

由于在月球上无法类似进行对地观测任务时布设控制点，因此本文采用CE-1获得的激光测高数据作为辅助数据使用计算定向参数，即有理函数系数。但是激光测高数据本身的获取过程中不可避免的会含有一定的误差，而且将激光测高数据投影到像片上时，其对应影像上的像点坐标的获得也存在误差，因此如果直接将这些激光测高数据在计算定向参数时作为控制点使用时，必然导致CCD影像的构像模型中的定向参数的求解结果不可靠。因此，在定向参数的求解过程中，把激光测高数据考虑成辅助约束条件使用。

在线阵CCD影像的定向参数的计算过程中，采用选权迭代的方法来剔除粗差，该方法的基本思想是：计算从最小二乘平差方法开始，在每次平差以后，根据观测值的残差和其他有关参数，按所选择的权函数计算每个观测值在下一次迭代计算中的权值。如果权函数选择合适，并且粗差可定位，则含粗差的观测值的权将越来越小，直至趋近于零。当迭代中止时，相应的残差将直接指出含粗差的数值，而平差结果将不再受粗差的影响。

在选权迭代中，权函数的选择也是多种多样的，本文采用丹麦法来确定权函数：

式中f(v)称为权函数，它是余差v的函数，P₀为权因子，σ为单位权中误差。

利用选权迭代法进行运算的第一步仍是进行经典的最小二乘法，从而获得各定向参数的改正数v_i，然后按照丹麦法将每个观测值重新赋以新的权值。利用这些新的权值再进行最小二乘平差，如此反复迭代直至收敛为止，进而得到有理函数系数的值，得到待匹配影像的成像模型。

（3）在投影光线上根据预先设定的月面点的近似高程Z₀以及高度误差ΔZ，结合搜索影像的构像模型在搜索影像上确定准核线；

近似高程和高程误差的选择不仅决定在匹配过程中的计算时间和速度，而且匹配结果将影响到待生成的DEM的精度。在本文中，利用嫦娥一号获得的激光测高数据构成三角网后，按照内插的方法获取特征点的近似高程。进行影像匹配时，利用参考影像上待匹配点的像点坐标和已知的影像定向参数元素，确定像点对应的月面空间点落在激光测高数据三角网的哪一个三角形内，通过该三角形的3个顶点的高程来确定近似高程Z₀和高程误差ΔZ。

（4）以准核线为中心，在其周围按照残余y视差扩展设置一邻域；

在准核线的搜索过程中，由于投影光线的确定是建立在构像模型的基础上，而构像模型中影像的定向参数是根据后方交会得到的，其值并不是严格意义下的真值，因此，正确的匹配点有可能不在通过上述流程计算得到的准核线上。如图4所示，对于月面点P，在CCD影像的获取时，它位于两幅影像上的同名像点为p和p′，而利用计算得到的定向参数，利用构像模型计算出点P在获取的CCD影像上的实际位置为p₁和p′₁。如果定向参数不含有误差，那么这两对点p和p₁、p′和p′₁应该是对应重合的。但实际上，由于定向参数解算过程中的线性化等方法导致定向参数含有一定的误差，那么p和p₁、p′和p′₁存在一定的偏差。在每张像片上，定义像点在像片y轴方向上的差别为残余y视差，即dy和d′y。

根据有理函数模型可以估算出像点存在的残余y视差的大小，从而按照残余y视差将准核线沿垂直于准核线方向扩展到一组平行线束，这组线束的宽度就是的两倍的残余y视差，那么同名点的搜索就可以沿着这组平行线束进行。

（5）计算邻域内部搜索影像上的特征点与参考影像上特征点之间的欧氏距离和正则化交叉相关系数NCC（Normalized Cross－Correlation）；

设I₀为参考影像，p₀为I₀上的一个像点,W为点p₀处的相关窗口，窗口大小为M×N，NCC可以表示为：

$\mathrm{NCC}(p_0,Z)=\frac{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)\×(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}{\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)}\·\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}}---\left(6\right)$

其中： ${\stackrel{\‾}{I}}_0=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_0\left(p_0\right)$ 和 ${\stackrel{\‾}{I}}_1=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_1\left(p_1\left(Z\right)\right)$ 为窗口内的灰度均值。

式（6）中，I₀(p₀)与I₁(p₁(Z))分别为参考影像点p₀处与搜索影像点p₁处的灰度值。NCC表示的是匹配过程中，考察一个立体像对的同名点的相似程度的尺度。

对于分别位于立体像对P₁和P₂上的特征点A和特征点B，定义D_Ai和D_Bi为这两点的描述子向量的第i个分量，则点A和点B之间的距离可以表示为：

$D(A,B)=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(D_{\mathrm{Ai}}-D_{\mathrm{Bi}})}^2}---\left(7\right)$

（6）对欧氏距离和NCC进行对比，如果两者均达到极值，则认为该点与参考影像上的特征点为一对同名像点；进而对所有提取的特征点计算得到最终的匹配结果。

5.粗差点剔除

应用RANSAC剔除粗差匹配点的方法如下：

（1）利用前述方法建立两幅影像上的同名点关系；

（2）任意选取3对同名点，建立两幅影像之间的仿射变换模型；

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{a}_0\\ b_0\end{array}\right]+\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]---\left(8\right)$

其中，(x，y)和(x′,y′)分布为同名点在两幅影像上坐标，(a_i,b_i)(i=0,1,2)为仿射变换参数。

（3）计算其余同名点对的匹配误差，根据事先设定的阈值，判断是否为内点，并统计内点个数；

（4）通过选取不同的同名点对，重复以上（2）—（3）步骤，直至所有点都经过一次计算；

（5）选取内点数量最多的内点集作为最佳匹配点集，由此可以计算得到基于准核线约束和SURF算子进行匹配的误匹配率为

6.最小二乘优化

在利用前述方法完成了匹配后，即可得到线阵CCD影像上的同名像点，下面采用最小二乘匹配进一步优化匹配结果。

(1)、最小二乘匹配原理

设g₀(x，y)为参考影像I₀上以p₀为中心的一个影像窗口，定义为参考影像窗口，g_i(x，y)(i=1,2)为搜索影像上以p_i为中心的一个影像窗口，定义为搜索影像窗口，(x，y)为影像的扫描坐标系坐标。在最小二乘平差中，g₀(x，y)被认为是g_i(x，y)的一个观测。而它们之间的辐射变形可以在计算相关系数时进行处理，这里可以不考虑，因此，它们之间的误差方程式为：

$\begin{array}{c}{v}_i(x,y)=\frac{{\∂g}_i}{{\∂x}_i}({\mathrm{da}}_0+x_0\·{\mathrm{da}}_1+y_0\·{\mathrm{da}}_2)+\frac{{\∂g}_i}{{\∂y}_i}({\mathrm{db}}_0+x_0\·{\mathrm{db}}_1+y_0\·{\mathrm{db}}_2)\\ +g_i(T_{\mathrm{Gx}}(x_0,y_0),T_{\mathrm{Gy}}(x_0,y_0))-g_0(x,y)\end{array}---\left(9\right)$

其中，a_i,b_i(i=0,1,2)是参考窗口和搜索窗口之间存在的仿射变换系数，T_G代表了窗口间存在的几何变形。

(2)、准核线几何约束条件

将几何约束条件加入到已有的平差模型中，可以提高匹配的精度与可靠性。

准核线方程可以表示为：

y_i=m_ix_i+n_i    （10）

如果我们已经得到像点p_i的近似坐标

则对上式线性化可以得到几何观测方程为：

$v_i(x,y)={\mathrm{dy}}_i-m_i{\mathrm{dx}}_i+y_i^0-m_i{x}_i^0-n_i---\left(11\right)$

灰度观测方程（9）与几何约束条件方程(11)共同组成了一个平差系统，其最小二乘解的结果即为第一次匹配的同名点结果。

7.结合Harris点特征的影像匹配

利用前文所述方法尽管提取了较多的特征点进行了影像匹配，但对于后续的摄影测量应用（如DEM生成）而言其分布仍不够均匀，不能满足实际需求，经过利用多种算子提取影像的特征点试验研究发现，Harris算子提取的特征点较为均匀，有利于后续作业。因此在前文方法的基础上，进一步提出了结合Harris点特征的影像匹配方法。该方法的主要思想是，首先利用SURF算子结合准核线约束的方法得到影像同名点信息；其次，利用Harris点特征提取算子提取参考影像的点特征；第三，利用这些同名点信息建立原始影像间的仿射变换关系，将提取的特征点按照仿射变换关系投影到搜索影像上，得到待匹配点信息；最后，以待匹配点为中心的局部范围内利用相关系数测度进行影像匹配，并利用最小二乘进行高精度匹配，最终得到高精度的匹配结果。其计算步骤为：

（1）、利用Harris特征点提取算子提取影像特征点

Harris算子的其原理是如果图像上某一点向任一方向作微小偏移都会引起灰度的较大变化，这样的点就定义为角点。其计算步骤为：

（a）计算图像在x和y方向的梯度g_x,g_y；

（b）计算自相关函数矩阵

$M=G\left(s\right)*\left[\begin{array}{cc}{g}_x^2& g_x{g}_y\\ g_x{g}_y& g_y^2\end{array}\right]---\left(12\right)$

其中，G(s)为高斯滤波器。

（c）计算每一个像点的兴趣值

R=det(M)-k(tr(M)²)        （13）

det(M)和tr(M)分别表示矩阵的行列式值和迹，k通常为常数。

（d）在局部范围内利用兴趣值进行局部非最大值抑制，寻找极值点作为特征点。

(2)、建立待匹配影像间的仿射变换关系

利用基于SURF和准核线约束的匹配方法得到的同名点，按照仿射变换模型（同（8）式）列出误差方程式，用最小二乘法计算6个仿射变换参数，得出待匹配影像间的对应关系。

(3)、根据仿射变换参数及特征点在参考影像上的位置，计算其对应像点在搜索影像上的大概位置。

(4)、相关系数计算

以待匹配点为中心，定义一窗口，计算窗口内部的各点的相关系数。

$\mathrm{\ρ}(r,c)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}\·g_{i+r,j+c}^{\′})-\frac{1}{m\×n}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}}{\sqrt{(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^2-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2)(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′2}-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′}\right)}^2)}}---\left(14\right),$

其中，g_i，j分别表示参考影像上行列分别为i，j的点的灰度值，g′_i，j分别表示搜索影像上行列分别为i，j的点的灰度值，m×n为在影像上定义的窗口大小。

在计算相关系数的过程中需要算出图像窗口的灰度和及灰度的平方和，实际上计算量较大，会影响匹配的速度。由于在进行SURF算子提取特征点时，生成了积分影像，而积分影像上计算灰度和时较快，因此计算相关系数时采用了积分图像的关系进行了快速计算。

(5)、利用最小二乘匹配的思想进行优化得到同名点信息。

(6)、结合之前SURF和准核线约束匹配是提取的同名点信息，得到最终匹配结果。

8、下面结合具体实验及实验数据说明本发明的效果。

选择嫦娥一号获得的两个具有月球典型形貌特征的区域进行了影像匹配实验。图5和图6分别显示了匹配后同名点信息的效果图，其中图5为山脉地区，图6为月海地区。从图中可以看出，利用本方法进行影像匹配后的点位分布较均匀，能较好的覆盖重叠区域，有利于后续的DEM生成。

图7和图8分别为图6地区的原始影像和经过预处理（wallis滤波）后的影像，图7为参考影像（其中（a）是原始的参考影像，而（b）为经过预处理后的影像），图8为搜索影像，（其中（a）是原始的搜索影像，而（b）为经过预处理后的影像）。

图9为在图7基础上利用SURF提取特征点后经过初次匹配的结果，首先在参考影像（左像）上利用SURF算子提取特征点，然后再搜索影像（右像）上进行匹配得到同名点信息。

图10为利用初次匹配的结果提取Harris角点进行匹配的结果，首先在参考影像（左像）上利用Harris算子提取特征点，然后再搜索影像（右像）上进行匹配得到同名点信息。将图9和图10同名点数据结合即为图6所示最终结果。

Claims (9)

1.一种影像特征匹配方法，其特征在于，该方法的步骤如下：

（1）对获取的CCD影像进行预处理；

（2）利用SURF算子对预处理后的CCD影像提取特征点，利用准核线约束条件和最小欧氏距离条件进行影像匹配，获取同名点信息；

（3）根据获得的同名点信息建立CCD影像间的仿射变形关系；

（4）利用Harris角点提取算子提取参考影像上特征点，并利用仿射变换投影至搜索影像上得到待匹配点；

（5）在待匹配点周围一邻域内，利用积分影像计算特征点与邻域内各点之间的相关系数，取极值点作为其同名点；

（6）将步骤（2）和步骤（5）匹配的结果共同作为最终同名点信息。

2.根据权利要求1所述的方法，其特征在于：所述步骤（1）中使用Wallis滤波器对CCD影像进行预处理。

3.根据权利要求1所述的方法，其特征在于：所述步骤（2）利用SURF算子提取特征点是计算原始影像的积分影像，然后建立尺度空间，使用Hessian矩阵的行列式同时来确定每个尺度上特征点的位置和尺度，当Hessian矩阵的行列式局部取最大时，检测出的像点做为特征点。

4.根据权利要求2所述的方法，其特征在于：所述步骤（2）利用准核线约束条件和最小欧氏距离条件进行影像匹配的过程如下：

（21）根据参考影像构像模型，对于参考影像上的每个特征点建立投影光线；

（22）在投影光线上根据预先设定的月面点的近似高程Z₀以及高度误差ΔZ，结合搜索影像的构像模型在搜索影像上确定准核线；

（23）以准核线为中心，在其周围按照参与y视差扩展设置一邻域；

（24）计算邻域内部搜索影像上的特征点与参考影像上特征点之间的欧氏距离和正则化交叉相关系数NCC；

（25）对欧氏距离和NCC进行对比，如果两者均达到极值，则该点与参考影像上的特征点位一对同名像点。

5.根据权利要求4所述的方法，其特征在于，所述步骤（24）中相关系数NCC的计算中，设I₀为参考影像，p₀为I₀上的一个像点，W为点p₀处的相关窗口，窗口大小为M×N，NCC可以表示为：

$\mathrm{NCC}(p_0,Z)=\frac{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)\×(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}{\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_0\left(p_0\right)-{\stackrel{\‾}{I}}_0)}\·\sqrt{\underset{p\∈W}{\mathrm{\Σ}}(I_1\left(p_1\left(Z\right)\right)-{\stackrel{\‾}{I}}_1)}}$

其中： ${\stackrel{\‾}{I}}_0=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_0\left(p_0\right)$ 和 ${\stackrel{\‾}{I}}_1=\frac{1}{M\×N}\underset{p\∈W}{\mathrm{\Σ}}{I}_1\left(p_1\left(Z\right)\right)$ 为窗口内的灰度均值，I₀(p₀)与I₁(p₁(Z))分别为参考影像点p₀处与搜索影像点p₁处的灰度值。

6.根据权利要求1-5中任一项所述的方法，其特征在于，所述步骤(2)完成后，在步骤(3)开始之前对同名点进行粗差剔除，并进行最小二乘优化。

7.根据权利要求6所述的方法，其特征在于，所述对同名点进行粗差剔除是应用RANSAC剔除粗差匹配点，过程如下：

(221)建立两幅影像上的同名点关系；

(222)任意选取3对同名点，建立两幅影像之间的仿射变换模型；

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{a}_0\\ b_0\end{array}\right]+\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]---\left(8\right)$

其中，(x，y)和(x′，y′)分布为同名点在两幅影像上坐标，(a_i，b_i)(i＝0，1，2)为仿射变换参数。

(223)计算其余同名点对的匹配误差，根据事先设定的阈值，判断是否为内点，并统计内点个数；

(224)通过选取不同的同名点对，重复以上(222)-(223)步骤，直至所有点都经过一次计算；

(225)选取内点数量最多的内点集作为最佳匹配点集。

8.根据权利要求6所述的方法，其特征在于，所述步骤(4)中利用Harris角点提取算子提取参考影像上特征点的过程如下：

(41)计算图像在x和y方向的梯度g_x，g_y；

（42）计算自相关函数矩阵 $M=G\left(s\right)*\left[\begin{array}{cc}{g}_x^2& g_x{g}_y\\ g_x{g}_y& g_y^2\end{array}\right],$ 其中，G(s)为高斯滤波器；

（43）计算每一个像点的兴趣值R=det(M)-k(tr(M)²)，det(M)和tr(M)分别表示矩阵的行列式值和迹，k为常数；

（44）在局部范围内利用兴趣值进行局部非最大值抑制，寻找极值点作为特征点。

9.根据权利要求8所述的方法，其特征在于，所述步骤（5）中利用积分影像计算特征点与邻域内各点之间的相关系数是，以待匹配点为中心，定义一窗口，采用积分图像计算窗口内部的各点的相关系数，

$\mathrm{\ρ}(r,c)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}(g_{i,j}\·g_{i+r,j+c}^{\′})-\frac{1}{m\×n}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i+r,j+c}^{\′}}{\sqrt{(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^2-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}\right)}^2)(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′2}-\frac{1}{m\×n}{\left(\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_{i,j}^{\′}\right)}^2)}}$

其中，g_i，j分别表示参考影像上行列分别为i，j的点的灰度值，g′_i，j分别表示搜索影像上行列分别为i，j的点的灰度值，m×n为在影像上定义的窗口大小。

Images