CN102855628B - 多源多时相高分辨率卫星遥感影像自动匹配方法 - Google Patents

Abstract

本发明公开了一种多源多时相高分辨率卫星遥感影像自动匹配方法，包括步骤：步骤一，根据参考影像对待匹配影像进行纠正和初步配准得到立体影像对；步骤二，从立体影像对中截取影像块，对各影像块进行如下操作：计算影像块的梯度从而得到其对应的梯度图像，并确定所得梯度图像的边缘区域；对影像块梯度图像中的边缘区域和非边缘区域加权，计算加权后影像块梯度图像的非线性亮度相关，得到匹配点；对所得匹配点进行粗差剔除，获取可靠的匹配点；步骤三，根据所得匹配点完成待匹配影像与参考影像的匹配。本方法采用非线性加权相关的方法，考虑了亮度的非线性变化，当影像亮度发生较大非线性变化时，也可以获得稳健的影像匹配结果。

Description

多源多时相高分辨率卫星遥感影像自动匹配方法

技术领域

本发明属于摄影测量与遥感领域，尤其涉及一种多源多时相高分辨率卫星遥感影像自动匹配方法。

背景技术

地面分辨率高于5米的卫星遥感影像为高分辨率卫星遥感影像，其具有覆盖范围广、成像周期短等优点，近年来在中小比例尺DEM(Digital Elevation Model，数字高程模型)、DOM（Digital Orthophoto Map，数字正射影像）制作、地形测绘、变化检测、地图更新等领域得到了广泛的应用。而快速、自动、稳定的多源多时相高分辨率卫星遥感影像匹配技术，是实现新影像纠正、对地定位及后继应用的前提。多源多时相影像匹配的关键技术在于能够克服影像间不同程度的几何和辐射差异,尤其是非线性差异，获取可靠的、高精度的匹配点。

在摄影测量领域，自最小二乘匹配算法提出以来，影像匹配都以像素亮度的线性相似度（CC，Correlation coefficient）作为匹配测度，更多的工作被集中在影像间几何关系的使用。在计算机视觉领域，大量的特征提取方法和描述算子被相继提出用以对影像进行匹配。这些特征提取方法和描述算子致力于发现立体像对在特定几何变形下的不变量，包括缩放不变、旋转不变、仿射不变等。如SIFT（Scale-invariant feature transform，尺度不变特征转换)算法主要考虑立体像对的缩放不变、旋转不变，并顾及仿射不变；SURF（Speeded-Up RobustFeatures，快速稳健特征）则重点考虑立体像对的缩放不变，顾及旋转不变。同摄影测量相似，在计算机视觉领域中的影像匹配通常也是认为影像亮度是线性变化的，或仅仅采用轻微的抑制手段，如抑制局部最大。

总之，无论在摄影测量领域还是计算机视觉领域，目前还未考虑“当影像亮度发生较大非线性变化时，如何获得稳健的影像匹配结果”这一问题。

发明内容

本发明的目的在于提供一种考虑了影像亮度非线性变化的、多源多时相高分辨率卫星遥感影像自动匹配方法。

为达到上述目的，本发明的技术方案为：

多源多时相高分辨率卫星遥感影像自动匹配方法，包括步骤：

步骤一，根据参考影像对待匹配影像进行纠正和初步配准得到立体影像对，所述的立体影像对包括参考影像和待匹配影像；

步骤二，从立体影像对中截取影像块，对各影像块进行如下操作：

计算影像块的梯度从而得到其对应的梯度图像，并确定所得梯度图像的边缘区域；

对影像块梯度图像中的边缘区域和非边缘区域加权，计算加权后影像块梯度图像的非线性亮度相关，得到匹配点；

对所得匹配点进行粗差剔除，获取可靠的匹配点；

步骤三，根据所得匹配点完成待匹配影像与参考影像的匹配。

步骤一中根据参考影像间的几何关系对待匹配影像进行纠正和初步配准。

步骤二中确定影像块梯度图像的边缘区域进一步包括子步骤：

2-1设定影像块梯度图像中边缘像素占总像素的比例K；

2-2统计影像块梯度图像中所有像素的灰度得到像素灰度直方图；

2-3根据预设K值和像素灰度直方图得到影像块梯度图像的边缘像素阈值g₀；

2-4逐一比较影像块梯度图像中各像素的灰度绝对值和边缘像素阈值g₀的大小，灰度绝对值大于边缘像素阈值g₀的像素则为边缘像素，所有边缘像素的集合即为边缘区域。K值一般为5～10%。

为了增强匹配的稳健性，在确定影像块梯度图像的边缘区域后，还有影像块梯度图像预处理步骤，具体为：将梯度图像中的非边缘区域滤除，将梯度图像中像素灰度绝对值小于某预设值δ的边缘像素滤除，对梯度图像中亮度高于某预设亮度值的边缘区域，采用“抑制局部最大”的方法进行抑制。作为优选，上述某预设值为0<δ≤5。

步骤二中对影像块梯度图像的边缘和非边缘区域加权，所述的边缘区域和非边缘区域的权值为：

将非边缘区域的权值衡定为1，边缘区域的权值为W，W根据关注边缘区域变化的程度来取值，其取值范围为10～100，在W的取值范围内，关注边缘区域变化程度越高W取值越大。

步骤二中计算加权后影像块梯度图像的非线性亮度相关，具体为：

对影像块梯度图像进行加权，得到非线性亮度相关ρ：ρ＝k₁ρ_x+k₂ρ_y+k₃ρ_-x+k₄ρ_-y，其中：

k₁、k₂、k₃、k₄为归一化常数，k₁、k₂的比值根据参考影像块梯度和的比值来确定，m、n分别为影像块梯度图像的长和宽，i、j为影像块梯度图像中像素点的二维坐标，dx_ij、dy_ij分别表示像素（i、j）在x方向、y方向的梯度，所述的x方向、y方向分别表示从左往右、从上往下方向，作为优选，k₁、k₂的值不小于0.4；k₃、k₄之和小于0.2；

ρ_x、ρ_y、ρ_-x、ρ_-y分别表示影像块梯度图像x方向、y方向、-x方向、-y方向梯度的相关系数，所述的x方向、y方向、-x方向、-y方向分别表示从左往右、从上往下、从右向左、从下往上方向；ρ_x、ρ_y、ρ_-x、ρ_-y根据边缘区域信息进行非线性加权获得：

${\mathrm{\&rho;}}_p=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{d}p\&CenterDot;w_{\mathrm{ij}}\&CenterDot;\stackrel{\&OverBar;}{d}{p}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}}}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}p\&CenterDot;w_{\mathrm{ij}})}^2}\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}{p}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}})}^2}}$

ρ_p表示影像块p方向梯度的相关系数，m、n分别为影像块梯度图像的长和宽，分别表示参考影像块、待匹配影像块p方向的重心化梯度，w_ij为参考影像块中像素（i，j）的权值；w_ij为待匹配影像中像素（i，j）的权值，p可取x、y、-x、-y。

步骤二中粗差剔除是根据是否已知待匹配影像的定向参数来自动选择粗差剔除几何模型。优选方案为：若不知待匹配影像的定向参数，则自动选择采用二次多项式进行粗差剔除；若已知待匹配影像的定向参数，则根据RPC（Rational Polynomial coefficients，有理多项式系数）参数进行粗差剔除。

传统的多源遥感影像匹配方法采用了线性相关的方法，受到不同时期、不同光照条件下亮度的非线性影响，从而导致匹配失败。本方法采用非线性加权相关的方法，考虑了亮度的非线性变化，以边缘的改变作为真实的变化，可以解决传统方法无法解决的问题。当影像亮度发生较大非线性变化时，采用本发明匹配方法也可以获得稳健的影像匹配结果。

附图说明

图1为本发明一种具体实施的流程图；

图2为第一组试验中试验区域的匹配结果，其中，图（a）为ADS40参考影像的整幅图像，图（b）为ADS40参考影像的试验区域，图（c）为待匹配CBERS-02B影像的整幅图像，图（d）为待匹配CBERS-02B影像的试验区域；

图3为第二组试验中试验区域的匹配结果，其中，图（a）为15幅胶片航空参考影像，图（b）为胶片航空参考影像的试验区域，图（c）为待匹配SPOT5HRS影像的整幅图像，图（d）为待匹配SPOT5HRS影像的试验区域。

具体实施方式

本发明提供了一种多源多时相高分辨率卫星遥感影像的自动匹配方法，首先，根据已知多源影像间的几何关系对待匹配影像进行纠正和初步配准；然后，从纠正和初步配准后的新影像中截取若干影像块，计算各影像块的梯度得到各影像块对应的梯度图像，并根据梯度图像确定各影像块的边缘；接着，对各影像块中的边缘和非边缘区域进行加权，计算加权后的非线性亮度相关；最后，完成粗差剔除，获取可靠的匹配点。流程图参见图1。

（1）待匹配影像的初步配准

若已知待匹配影像的定向参数，则利用待匹配影像的定位参数预测同名点，然后根据仿射变换对待匹配影像和参考影像进行采样，完成待匹配影像和参考影像的初步配准；若不知待匹配影像的定向参数，则在待匹配影像和参考影像上手工选择3对同名点，然后根据仿射变换对待匹配影像和参考影像进行采样，完成待匹配影像和参考影像的初步配准。

（2）计算梯度图像

从立体影像对中截取影像块，得到立体影像块（文中将立体影像块均简称为影像块），计算影像块中每一像素（i，j）不同方向的梯度dx_ij、dy_ij、-dx_ij、-dy_ij，dx_ij、dy_ij、-dx_ij、-dy_ij分别表示像素（i，j）从左往右、从上往下、从右向左、从下往上方向的梯度，这里的立体影像对是对待匹配影像进行纠正和初步配准后得到的立体影像对，其包括参考影像和待匹配影像；然后，对上述4个方向所对应的4幅影像块梯度图像向下重采样，以消除初步配准的几何误差，可取2×2或3×3向下采样。

（3）梯度图像的边缘区域

在梯度图像中，像素灰度绝对值较大的像素被认为是边缘像素，而边缘像素的集合即为边缘区域，本发明中以此思想为基础来确定预处理后的梯度图像中的边缘区域。

首先找到边缘像素阈值g₀，边缘像素阈值g₀通过统计像素灰度直方图得到，具体为：根据实际情况和经验，设定影像块梯度图像中边缘像素占总像素的比例K，一般地，边缘像素占整个图像像素的比例是较小的，因此，K取值为5～15%较合适。统计向下重采样后的影像块梯度图像中所有像素的灰度值得到像素灰度直方图，根据像素灰度直方图和预设K值得到该影像块梯度图像的边缘像素阈值g₀。

确定边缘像素阈值g₀后，逐一比较该影像块梯度图像中各像素的灰度绝对值和边缘像素阈值g₀的大小，灰度绝对值大于边缘像素阈值g₀的像素则为边缘像素，否则为非边缘像素。所有边缘像素的集合即为边缘区域。

（4）梯度图像预处理

为了增强匹配的稳健性，在确定影像块梯度图像的边缘区域后，对影像块梯度图像进行预处理，具体为：

将影像块梯度图像中像素灰度绝对值小于2的非边缘像素滤除。对于梯度图像中亮度高于某预设亮度值的高亮边缘区域，采用“抑制局部最大”的方法进行抑制，例如，当边缘区域中像素梯度值大于50时，重新设定该像素梯度值为50。

（5）对边缘和非边缘区域加权

将非边缘区域的权值衡定为1，对边缘区域给予较大的权值，边缘区域的权值W确定了图像亮度特征与边缘特征之间的关系，W在10～100间取值，W的取值和目的有关，如果想更加关注于边缘的变化，W可以取较大值；如果不想关注边缘的变化，则可以取较小值。

对边缘和非边缘像素的加权可以采用以下公式表示：

其中，（i、j）表示梯度图像中各像素点；g(i,j)为像素（i，j）的像素灰度绝对值；w_ij为像素（i，j）的权值。

（6）计算非线性亮度相关，得到匹配点

对影像块进行非线性亮度相关是通过加权相关完成，各影像块的非线性亮度相关ρ为：ρ＝k₁ρ_x+k₂ρ_y+k₃ρ_-x+k₄ρ_-y，

其中，ρ_x、ρ_y、ρ_-x、ρ_-y分别表示影像块梯度图像x方向、y方向、-x方向、-y方向梯度的相关系数，所述的x方向、y方向、-x方向、-y方向分别表示从左往右、从上往下、从右向左、从下往上方向；ρ_x、ρ_y、ρ_-x、ρ_-y根据边缘区域信息进行非线性加权获得，下面以ρ_x为例来说明ρ_x、ρ_y、ρ_-x、ρ_-y的计算方法：

${\mathrm{\&rho;}}_x=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{d}x\&CenterDot;w_{\mathrm{ij}}\&CenterDot;\stackrel{\&OverBar;}{d}{x}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}}}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}x\&CenterDot;w_{\mathrm{ij}})}^2}\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}{x}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}})}^2}}---\left(8\right)$

m、n分别为影像块梯度图像的长和宽；表示参考影像块的重心化梯度；为待匹配影像块的重心化梯度；w_ij为参考影像块中像素（i，j）的权值；w′_ij为待匹配影像中像素（i，j）的权值；

重心化梯度和的计算公式如下：

$\stackrel{\&OverBar;}{d}x=d{x}_{\mathrm{ij}}-\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{j=n}{\mathrm{\&Sigma;}}}{\mathrm{dx}}_{\mathrm{ij}}---\left(9\right)$

${\stackrel{\&OverBar;}{d}}^{\&prime;}x={{\mathrm{dx}}_{\mathrm{ij}}}^{\&prime;}-\frac{1}{\mathrm{mn}}\underset{i=1}{\overset{i=m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{j=n}{\mathrm{\&Sigma;}}}{{\mathrm{dx}}_{\mathrm{ij}}}^{\&prime;}---\left(10\right)$

k₁、k₂、k₃、k₄为归一化常数，k₁、k₂的比值根据参考影像块梯度和的比值来确定，m、n分别为参考影像块梯度图像的长和宽，i、j为参考影像块中像素点的二维坐标，dx_ij、dy_ij分别表示像素（i、j）在x方向、y方向的梯度，所述的x方向、y方向分别表示从左往右、从上往下方向，作为优选，k₁、k₂的值不小于0.4；k₃、k₄之和小于0.2，作为一种优选取值，可取k₃=k₄=0.08。

（7）粗差剔除

具体实施时，根据是否已知待匹配影像的定向参数来自动选择粗差剔除几何模型：

若不知待匹配影像的定向参数，则自动选择采用二次多项式进行粗差剔除：

x＝e₀+e₁·x₁+e₂·y₁+e₃x₁ ²+e₄y₁ ²+e₅x₁y₁   (11)

y＝f₀+f₁·x₁+f₂·y₁+f₃x₁ ²+f₄y₁ ²+f₅x₁y

其中，e₀,e₁,e₂,e₃,e₄,e₅,f₀,f₁,f₂,f₃,f₄,f₅为二次多项式系数，通过求解方差获得；x₁,y₁是参考影像上的像素坐标，x,y是待匹配影像上与x₁,y₁对应的匹配坐标；

若已知待匹配影像的定向参数，则根据RPC（Rational Polynomial coefficients，有理多项式系数）参数与仿射变换（Affine transformation）（以下简写为RPC+A）进行粗差剔除：

y＝e₀+e₁·sample+e₂·line

x＝f₀+f₁·sample+f₂·line    (12)

其中，e₀,e₁,e₂,f₀,f₁,f₂为仿射变换参数，通过求解方程计算得到；line,sample是由地理参考信息逐点计算出的待匹配影像的像素坐标，而匹配获取的待测影像扫描坐标x,y作为观测值。

下面将采用验证试验来说明本发明的有益效果。

为验证本发明在不同传感器光学遥感影像的匹配能力，选用两套控制点数据和对应区域的新影像进行试验。第一组试验中控制点选自航空ADS40数码条带式扫描影像，地面坐标已知，定位精度约0.5m；对应的待匹配遥感影像为中国高分辨率卫星CBERS02B影像，地面分辨率2.3m，见图2。第二组试验中控制点选自胶片框幅式航摄影像，地面坐标已知，定位精度约1m，对应的待匹配遥感影像为法国SPOT5HRS影像，地面分辨率5m，见图3。从图中可见两组试验中的试验区域的影像都具有很大的几何形变和亮度变化，且地物有一定的变化。

第一组试验中，首先从ADS40影像中均匀提取了46个控制点，并通过原始定位参数投影到CBERS影像上。由于控制点影像分辨率高，这里选择对控制点影像进行重采样。然后采用本发明方法对这46个控制点进行匹配。初步获得35个匹配点，经过RPC粗差剔除，最终获得31个正确匹配的控制点，部分匹配点见图2十字丝所示，定位精度见表1所示。

为了验证密集控制点和稀疏控制点两种情形下的匹配和定位效果，在第二组试验中，SPOT5影像的试验区密集选取了706个控制点。初步匹配获取了244对同名点，剔除粗差后最终保留了182对匹配点，部分匹配点见图3十字丝。匹配前后的定位精度见表1。表1反映的只是几何方程的像方内符合精度（内精度），为验证在有地理参考下的实际对地定位精度（外精度），分别采用手工精确量测的数个控制点作为检查点作验证。精度如表2所示。

表1第一组试验像空间匹配精度比较(*匹配前的几何误差)

表2对地目标定位精度比较

通过表1、2实验结果看出，经过剔粗差后，SPOT5像方平面中误差约0.6像素，CBERS像方平面中误差为1.5像素。该结果验证了在不同传感器卫星，不同控制点数量的情况下，

本发明方法在多源多时相影像匹配中的可靠性和较高的匹配精度。SPOT5和CBERS检查点残差分别为4.4m和3.8m，明显高于原始辅助数据的定位精度，可满足1：50000地形图的测图精度要求，而无需人工布设和量测新的控制点。此外，无论是采用大量控制点还是稀疏控制点，本发明都可适应不同控制点数量下粗差剔除和精确定位的需要。

Claims (6)

1.一种多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于，包括步骤：

步骤一，根据参考影像对待匹配影像进行纠正和初步配准得到立体影像对，所述的立体影像对包括参考影像和待匹配影像；

步骤二，从立体影像对中截取影像块，对各影像块进行如下操作：

计算影像块的梯度从而得到其对应的梯度图像，并确定所得梯度图像的边缘区域；

对影像块梯度图像中的边缘区域和非边缘区域加权，计算加权后影像块梯度图像的非线性亮度相关ρ＝k₁ρ_x+k₂ρ_y+k₃ρ_-x+k₄ρ_-y，得到匹配点，其中：

影像块梯度图像中非边缘区域的权值衡定为1，边缘区域的权值为W，W根据关注边缘区域变化的程度来取值，其取值范围为10～100，在W的取值范围内，关注边缘区域变化程度越高W取值越大；

计算加权后影像块梯度图像的非线性亮度相关具体为：

对影像块梯度图像进行加权，得到非线性亮度相关ρ：ρ＝k₁ρ_x+k₂ρ_y+k₃ρ_-x+k₄ρ_-y，k₁、k₂、k₃、k₄为归一化常数，k₁、k₂的比值根据参考影像块梯度和的比值来确定，m、n分别为影像块梯度图像的长和宽，i、j为影像块梯度图像中像素点的二维坐标，dx_ij、dy_ij分别表示像素(i、j)在x方向、y方向的梯度，所述的x方向、y方向分别表示从左往右、从上往下方向；k3、k4之和小于0.2；

ρ_x、ρ_y、ρ_-x、ρ_-y分别表示影像块梯度图像x方向、y方向、-x方向、-y方向梯度的相关系数，所述的x方向、y方向、-x方向、-y方向分别表示从左往右、从上往下、从右向左、从下往上方向；ρ_x、ρ_y、ρ_-x、ρ_-y根据边缘区域信息进行非线性加权获得：

${\mathrm{\&rho;}}_p=\frac{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{d}p\&CenterDot;w_{\mathrm{ij}}\&CenterDot;\stackrel{\&OverBar;}{d}{p}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}}}{\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}p\&CenterDot;w_{\mathrm{ij}})}^2}\sqrt{\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{(\stackrel{\&OverBar;}{d}{p}^{\&prime;}\&CenterDot;{w^{\&prime;}}_{\mathrm{ij}})}^2}}$

ρ_p表示影像块p方向梯度的相关系数，m、n分别为影像块梯度图像的长和宽，分别表示参考影像块、待匹配影像块p方向的重心化梯度，w_ij为参考影像块中像素(i，j)的权值，w′_ij为待匹配影像中像素(i，j)的权值，p可取x、y、-x、-y；

对所得匹配点进行粗差剔除，获取可靠的匹配点；

步骤三，根据所得匹配点完成待匹配影像与参考影像的匹配。

2.根据权利要求1所述的多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于：

所述的纠正和初步配准是根据参考影像间的几何关系进行。

3.根据权利要求1所述的多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于：

所述的确定影像块梯度图像的边缘区域进一步包括子步骤：

2-1设定影像块梯度图像中边缘像素占总像素的比例K；

2-2统计影像块梯度图像中所有像素的灰度得到像素灰度直方图；

2-3根据预设K值和像素灰度直方图得到影像块梯度图像的边缘像素阈值g₀；

2-4逐一比较影像块梯度图像中各像素的灰度绝对值和边缘像素阈值g₀的大小，灰度绝对值大于边缘像素阈值g₀的像素则为边缘像素，所有边缘像素的集合即为边缘区域。

4.根据权利要求3所述的多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于：

所述的影像块梯度图像中边缘像素占总像素的比例K值为5～10％。

5.根据权利要求1所述的多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于：

所述的对所得匹配点进行粗差剔除根据是否已知待匹配影像的定向参数来自动选择粗差剔除几何模型。

6.根据权利要求5所述的多源多时相高分辨率卫星遥感影像自动匹配方法，其特征在于：

所述的根据是否已知待匹配影像的定向参数来自动选择粗差剔除几何模型具体为：若不知待匹配影像的定向参数，则自动选择采用二次多项式进行粗差剔除；若已知待匹配影像的定向参数，则根据有理多项式参数进行粗差剔除。