CN104299228B - 一种基于精确点位预测模型的遥感影像密集匹配方法 - Google Patents

Abstract

本发明涉及一种基于精确点位预测模型的遥感影像密集匹配方法，属于遥感测量技术领域。本发明首先提取影像特征点作为初始匹配同名点；提取待匹配影像的特征点P_i；利用仿射变换模型和初始匹配同名点对所提取的待匹配特征点P_i进行同名点预测，以及相关系数匹配，得到的P_i同名点，并将得到同名点作为下一层加密匹配的已知点，重复特征点提取、点位预测和相关系数匹配过程，进行逐层加密匹配，直至获取满足需要的匹配点数量。本发明所采用的逐层加密匹配以及精确点位预测机制能够有效解决火星表面影像密集匹配的问题，且随着加密匹配过程中同名点的逐渐增多，点位预测精度逐步提高，最终可达到1～3个像素的高精度点位预测结果。

Description

一种基于精确点位预测模型的遥感影像密集匹配方法

技术领域

本发明涉及一种基于精确点位预测模型的遥感影像密集匹配方法，属于遥感测量技术领域。

背景技术

火星地形测绘是开展火星科学研究的基础，火星探测任务实施以来，与地形测绘相关的火星探测任务主要有海盗号火星探测器、火星全球勘测者轨道探测器、火星快车探测器、火星侦察轨道器。相比其它火星探测器，火星快车携带的高分辨率立体测绘相机(High Resolution Stereo Camera，HRSC)采用三线阵原理，可构成同轨立体，在影像分辨率及全球覆盖性方面具有优势。

目前国内行星测绘的研究对象主要集中在月球。王任享对月球三线阵CCD影像EFP光束法平差进行了研究与实践。李春来等人研究了嫦娥一号三线阵CCD数据摄影测量处理及全月球数字地形图生成技术。周杨等人研究了大数据量月面形貌的三维可视化方法。赵双明等人研究了立体影像与激光高度计数据的不一致性。赵葆常等人介绍了嫦娥二号卫星CCD立体相机设计方法。崔平远等人在研究火星精确着陆制导问题的基础上，分析了火星好奇号探测器的着陆精度及着陆区地形。

发明内容

本发明的目的是提供一种基于精确点位预测模型的遥感影像密集匹配方法，以有效解决火星表面影像密集匹配的问题。

本发明为解决上述技术问题而提供一种基于精确点位预测模型的遥感影像密集匹配方法，该密集匹配方法包括以下步骤：

1)基于严密几何模型对HRSC影像进行近似核线重采样；

2)对重采样后的HRSC影像利用SURF算法提取影像特征点作为初始匹配同名点；

3)将得到的初始匹配同名点作为已知点进行第一次加密匹配，提取待匹配影像的特征点P_i；

4)利用仿射变换模型对所提取的待匹配特征点P_i进行同名点预测，并确定搜索窗口范围；

5)以预测点为中心，在预测点搜索窗口范围内计算各点与待匹配特征点P_i的相关系数，从而确定待匹配特征点P_i的同名点；

6)将得到同名点作为下一层加密匹配的已知点，重复步骤3)至步骤5)的特征点提取、点位预测和相关系数匹配过程，进行逐层加密匹配，直至获取满足需要的匹配点数量，或者影像中已经不能再提取出特征点。

所述步骤3)是采用用Shi-Tomasi算子进行特征点提取的。

所述步骤4)中的仿射变换模型为：

x₂＝a₀+a₁x₁+a₂y₁

y₂＝b₀+b₁x₁+b₂y₁

其中(x₁，y₁)为待匹配点坐标，(x₂，y₂)是同名点预测坐标，a₀、a₁、a₂、b₀、b₁和b₂为模型参数。

所述步骤4)中同名点的预测过程为：

a.对待匹配的特征点P_i以设定的搜索半径为步长从已知点中搜索与P_i邻近的n个点；

b.利用得到的n个点按照最小二乘原理解算仿射变换模型中模型参数；

c.根据解算后的模型参数计算仿射变换模型中的同名点预测坐标，从而确定同名点预测值。

所述步骤4)中搜索窗口的范围由邻近已知点与待匹配特征点P_i的距离以及核线、地形约束条件设定。

所述步骤a是采用KD树算法对邻近已知点进行搜索。

所述步骤1)中基于严密几何模型构建过程如下：

A.利用GDAL开源图像处理库将HRSC Level 2级PDS格式影像数据转换为TIFF格式；

B.利用SPICE库获取HRSC相机的焦距、像元大小、像元位置参数；

C.将HRSC原始PDS格式数据导入ISIS系统，导出扫描行时间文件；

D.利用SPICE库及卫星影像扫描行时文件获取每一扫描行的位置、姿态数据；

E.利用获取的相机几何参数、位置、姿态等信息构建HRSC影像严密几何模型。

所述步骤6)中当影像中已经不能在提取出特征点时，可在待匹配影像上按照一定的格网间距选择一个特征相对明显的像点作为待匹配点，并对获取的待匹配点进行点位预测、相关系数匹配以获取密集匹配点。

本发明的有益效果是：本发明根据火星表面影像匹配纹理稀疏和地形连续的特点提出了一种基于精确点位预测模型的密集匹配方法，该方法首先HRSC影像利用SURF算法提取影像特征点作为初始匹配同名点；将得到的初始匹配同名点作为已知点进行第一次加密匹配，提取待匹配影像的特征点P_i；利用仿射变换模型对所提取的待匹配特征点P_i进行同名点预测，以及相关系数匹配，并将得到同名点作为下一层加密匹配的已知点，重复特征点提取、点位预测和相关系数匹配过程，进行逐层加密匹配，直至获取满足需要的匹配点数量，或者影像中已经不能再提取出特征点。并以火星快车HRSC影像对本发明进行了验证，结果表明本发明所采用的逐层加密匹配以及精确点位预测机制能够有效解决火星表面影像密集匹配的问题，且随着加密匹配过程中同名点的逐渐增多，点位预测精度逐步提高，最终可达到1～3个像素的高精度点位预测结果，在此基础上利用相关系数沿核线匹配可快速、精确命中同名点。

附图说明

图1是HRSC相机CCD排列示意图；

图2是利用SPICE库获取HRSC外方位元素示意图；

图3是同名点位精确预测示意图；

图4是本发明基于精确点位预测模型的遥感影像密集匹配方法的流程图；

图5是基于KD数的邻近已知点搜索及点位预测效率变化示意图；

图6-a是H3304轨影像匹配示意图；

图6-b是H3304轨影像DEM结果示意图；

图6-c是H8433轨影像匹配示意图；

图6-d是H8433轨影像DEM结果示意图；

图7-a是H9465轨影像匹配示意图；

图7-b是H9465轨影像DEM结果示意图；

图8-a是H5273轨影像采用本发明匹配方法后提取DEM结果示意图；

图8-b是H5273轨影像欧空局提取DEM结果示意图；

图8-c是本发明DOM与DEM叠加效果示意图；

图9-a是DOM量测点示意图；

图9-b是DEM量测点示意图；

图9-c是本发明与欧空局的DOM结果对比图；

图9-d是本发明与欧空局的DEM结果对比图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的说明。

火星快车是欧洲的首个火星探测任务，搭载的HRSC相机影像分辨率达10m，且影像数据已经基本覆盖火星全球。HRSC相机CCD像素数为5184，焦距为175mm，影像辐射分辨率为12位，可同时获取9个波段线阵影像，即5个全色波段和4个多光谱波段，其全色影像分辨率高于多光谱影像，通过融合处理可以获取高分辨率多光谱影像，相机焦平面排列如图1所示。

卫星影像严密几何模型构建是后序摄影测量处理的基础，火星快车HRSC影像严密几何模型构建方法如下：

步骤1利用GDAL(开源图像处理库)将HRSC Level 2级PDS格式影像数据转换为TIFF格式；

步骤2利用SPICE库获取HRSC相机的焦距、像元大小、像元位置等参数；

步骤3将HRSC原始PDS格式数据导入ISIS系统，导出扫描行时间文件；

步骤4利用SPICE库及卫星影像扫描行时文件获取每一扫描行的位置、姿态数据；

步骤5利用获取的相机几何参数、位置、姿态等信息构建HRSC影像严密几何模型。

图2是利用SPICE库获取HRSC影像外方位元素的示意图，时间转换为J2000历书时，摄影测量物方坐标系转换至火星地固坐标系，HRSC影像严密几何模型公式如下：

其中(x，y)为像点坐标，f为相机焦距，[X Y Z]^T为火星地面点坐标，[X_S Y_S Z_S]^T为摄站位置，为J2000坐标系到火星地固坐标系旋转矩阵，为星敏感器测定的探测器姿态，为卫星本体至星敏感器的旋转矩阵，为相机与卫星本体间的安装矩阵。

火星表面影像纹理信息贫乏，影像特征点提取困难，这是火星表面影像匹配的主要难点。但是与对地观测卫星影像相比，火星表面影像匹配也有一些有利因素，比如火星上不存在建筑物遮挡、移动目标、河流、湖泊，而且火星表面地形比较连续，可以在匹配时利用这些有利条件。结合上述火星表面影像特点分析，本发明提出了一种基于精确点位预测模型的遥感影像密集匹配方法。

本实施例以立体影像I₁与I₂为例，其影像匹配过程如图4所示，具体步骤如下：

1.HRSC影像近似核线重采样，

基于严密几何模型对影像I₁与I₂进行近似核线重采样，针对HRSC影像设计了近似核线重采样方法，而HRSC核线重采样实际上也是几何纠正过程。

2.SURF匹配初始点

利用SURF算法提取影响特征点，计算特征描述向量并匹配特征点，设初始匹配同名点数为N₀。

3.提取加密匹配用的特征点

利用SURF匹配得出的N₀个同名点作为已知点进行第一次加密匹配，在影像I₁上利用Shi-Tomasi算子提取特征点。

常用的角点检测算法主要有Moravec算子、算子、Harris算子、SIFT算子、SURF算子等，Shi和Tomasi在Harris算子的基础上研究发现，若两个特征中较小的一个大于最小阈值，则会得到强角点，Shi-Tomasi算子大部分情况下可以得到优于Harris算子的效果，并且便于控制点位分布，因此本发明在逐层加密匹配时选用Shi-Tomasi算子提取特征点。

4.点位预测

对于每一个待匹配的角点，以搜索半径R(如50个像素)为步长从N₀个同名点中搜索与之邻近的n个点，将这n个点按照最小二乘原理解算仿射变换模型参数并对同名点进行预测，设预测点位为Q′_i，同时按照邻近已知点与Q′_i的距离以及核线、地形约束条件设定搜索窗口范围。

5.计算相关系数并确定同名点

以预测点Q′_i为中心，在该点搜索窗口范围内计算各点与Q′_i的相关系数ρ_i，一般当相关系数大于0.8时即认为是同名点，以此确定同名点。

相关系数测度是数字摄影测量中比较经典的测度函数，其数学公式描述如下：设影像匹配的目标窗口灰度矩阵为G(g_i,j)(i＝1,2,…,m；j＝1,2,…,n)，m与n分别是矩阵G的行列数，其构成的影像窗口称为匹配窗口(matching window)或者相关窗口，且一般取为奇数(如15×15)，与G相应的灰度函数为g(x,y)，(x,y)∈D，将G中元素排成一行构成一个N＝m×n的目标向量X＝(x₁,x₂,…,x_N)。而搜索区灰度矩阵为G′(g′_i,j)(i＝1,2,…,k；j＝1,2,…,l)，k与l是矩阵G′的行列数，其构成的影像窗口称为搜索窗口(search window)，一般也取为奇数，与G′相对应的灰度函数为g′(x′,y′)，(x′,y′)∈D′。G′中任意一个m行n列的子块可以表示为：

其中INT表示取整运算，将G′_r,c中的元素排成一行可同样构成一个N＝m×n的搜索向量，记为Y＝(y₁,y,…,y_N)。对于图像灰度数据，相关系数匹配测度公式如下：

6.逐层加密匹配

重复步骤3值步骤5的特征点提取、点位预测和相关系数匹配过程，进行逐层加密匹配，直至获取满足需要的匹配点数量，或者影像中已经不能再提取出特征点，此时，为获取更多同名点(如用于DEM生成)可转至步骤7。

7.匹配格网点

在影响I1上按照一定的格网间距(如)选择一个特征相对明显的像点作为待匹配点，对这些格网点进行点位预测、相关系统匹配以获取密集匹配点。

当经过逐层加密匹配得出一定数量的同名点后，待匹配点与邻近已知点一般距离较小(如50～100像素)，可使用仿射变换模型建立左右影像间的变换关系：

其中(x₁，y₁)为左影像I₁上待匹配点坐标，(x₂，y₂)是右影像I₂上同名点预测坐标，a₀、a₁、a₂、b₀、b₁和b₂为模型参数。公式中仿射变换模型参数共有6个未知数，而每个同名点可以列出两个方程，因此使用3个同名点即可解算，一般使用4个均匀分布的已知点，当搜索范围内有多个已知点时则采用最小二乘平差求解。解算出变换参数后，即可利用仿射变换公式由待匹配点坐标(x₁，y₁)预测出其同名点坐标(x₂，y₂)。实际上当待匹配点周围有足够已知点时，预测点位的精度可以达到1～3个像素，这就是本发明期望达到的“精确预测”效果。

在预测出点位之后，还需要根据邻近已知点与待匹配点的距离再结合核线、地形等约束条件设定合理的搜索窗口大小。因此，点位预测模型可以描述为：

①点位估计：

②核线、地形约束设定搜索窗口:k×l

其中k和l分别为搜索窗口宽度和搜索窗口高度。

通常情况下(平地、丘陵地形)当与邻近已知点的距离分别为50，100，200，400,…时，搜索窗口大于依次设置为5x3，9x3，15x3，21x3，当为山地地形时，搜索窗口扩大一倍，这些参数由多次试验确定。

本发明仅需要少量邻近已知点即可构建点位预测模型用于初次加密匹配，只是点数较少时点位预测精度稍差，此时可以将搜索窗口增大，仍然可以在后续逐层加密匹配时获得大量同名点。为进一步说明逐层加密匹配时对初始已知点的依赖程度，以H8433轨影像为例，在H8433轨影像四角及中心人工采集5个同名点，以这五个同名点作为初始已知点构建仿射变换预测模型，利用Shi-Tomasi算子在左影像上提取5000个特征点，右影像上共匹配出4685个特征点，匹配成功率为93.7％，可见，本发明逐层加密匹配对SURF初始匹配点的数量要求并不高。

邻近已知点的搜索效率是影响本发明匹配计算效率的一个关键因素，本实施例采用KD数(K-Dimension Tree)算法实现邻近已知点的快速搜索，下面通过试验分析基于KD树算法的邻近已知点搜索及点位预测效率。选取H5273轨影像进行试验，已知点数量从2069增大值144101，预测点数从1万逐渐增大至20万，点位预测效率如图5所示。从图5中可知，邻近已知点的搜索以及点位预测效率与已知点数量有关，图中数据走势一显示出来其增长管线基本上是线性的，即使已知点数量为14万时，20万个待匹配点的同名点预测时间约为30秒钟，点位预测效率是比较高的。

下面选取火星快车HRSC影像Level 2级数据进行试验，分析匹配算法点位预测精度、抗粗差能力以自动提取DEM的效果，试验算法程序设计平台为Visual Studio 2010+Qt，试验硬件环境为Windows 8操作系统，CPU为Intel Core i5-2450，主频2.50GHz，内存4GB。

点位预测精度分析

通过试验分析预测点位与实际点位的误差，其中实际点位是在预测点位的基础上由相关系数普洱陪得出。选取H5273轨影像进行两组试验，第一组已知点数量为9455，第二组已知点数量为44147，使用不同的搜索窗口分析匹配成功率以及点位预测误差的最大值(xmax、ymax)、最小值(xmin、ymin)、均值(xmean、ymean)和中误差(xstd、ystd)等指标，试验结果如表1所示。

表1

对表1中的试验结果进行分析可以得到：

1)点位预测精度X方向中误差为1.2～1.5个像素，Y方向中误差为0.5～0.7个像素，大部分点位预测误差在2个像素以内，仅少量点预测误差在2～4个像素，因此，可以证明本发明能够达到同名点精确预测的效果。

2)匹配成功率的分析，由于点位预测精度较高，增大搜索窗口基本不影响匹配成功率，因此本发明在逐层加密匹配后期可以限定一个较小的搜索窗口，以提升匹配计算效率。

3)核线几何约束分析，点位预测误差在核线影像Y方向较小，试验中搜索窗口沿核线影像Y方向从9×3变化至9×9，但是并未影响匹配成功率，精确点位预测模型结合核线约束共同形成高精度点位预测效果。

粗差对点位预测的影响

基于本发明中采用的逐层加密与邻近已知点精确预测的匹配策略，即使已知点数据中有少量粗差点，实际上对同名点预测精度的影响并不大，且在逐层加密匹配过程中如果粗差点导致预测精度降低，其在后续的相关系数匹配过程中是难以找出同名像点的，这等于在逐层加密匹配时舍弃了上一层粗差点。例如，在H8433轨影像上人工加入一个粗差点，粗差点的误差幅度约20个像素，已知点本身数量是9321，在此基础上加密匹配得出18766个同名点，试验结果表明粗差点未影响后序加密匹配，因此本发明的点位预测机制以及逐层加密的匹配策略具有潜在的抵抗粗差的能力。

火星表面影像匹配及DEM自动提取试验

选取4轨影像进行DEM自动提取试验，试验数据基本信息及结果见表2。H3304与H8433轨由DEM逐层加密匹配两次得出，其结果如图6-a至图6-d所示；而H9466轨DEM首先经过两次逐层加密匹配，然后利用格网点进一步加密匹配得出，结果如图7-a和图7-b所示。

表2

由图6-a至6-b逐层加密匹配直接构建出的DEM结果可知，H3304轨DEM接近于10×10网格内一个匹配点，H8433轨DEM接近7×7各网内一个匹配点，可见H8433轨影像匹配结果更为密集，主要是该轨影像信息更为丰富。虽然H3304匹配点数密集程度不如H8433轨影像，但是其地形相对简单，因此即使匹配点数较少，也不影响地形构建效果，如图6-b所示。从中可以看出当地形简单或者DEM精度要求不高时，利用逐层加密匹配得出的特征点能够直接构建DEM。

图7-a、图7-b、图8-a和图8-b是逐层加密匹配并结合格网点构建出的DEM，从中可以得出，对于纹理相对稀疏区域，在逐层加密匹配的基础上进一步利用格网点匹配可以得到更为密集的匹配点。图8-a和图8-b是采用本发明方式提取DEM结果与欧空局结果对比(H5273轨)，可以看出，本发明基于格网点加密匹配出的地形更为精细，图8-c为DOM与DEM叠加显示效果，可见地形与影像吻合程度较好，H5237轨影像位于好奇号着陆区盖尔陨石坑，图8-c中红色椭圆区域即为着陆区。

由于目前欧空局并未发布H8433与H9465轨测绘产品，因此精度分析仅利用H5273轨影像(与H3304轨相比地形相对复杂，影像范围较大)，基于H5273轨提取的DEM纠正生成DOM，并将本发明DOM与DEM结果与欧空局结果进行对比，数据处理时火星椭球参数采用正球体定义，半径为3396.0km，精度统计结果如图9-a和图9-b所示。

H5273轨生成DOM分辨率为25m，DEM分辨率为50m，分析图9-c和图9-d中的数据可知，经度方向误差最大值为269m，最小值为4m，中误差为95m，约34个像素；纬度方向误差最大值为244m，最小值为93m，中误差为38m，约为1.5个像素；高程方向误差最大值为124m，最小值为4m，中误差为52m，约为1个像素(相对于DEM格网间距)。本发明试验结果与欧空局结果在平面精度上有明显系统偏差，且在精度方向偏差较大，约为2-3个像素，而纬度方向系统性偏差约为1个像素，高程方向未见明显系统性偏差。而产生系统性偏差的原因主要包括以下两个方面。

1)本发明试验所使用的相机位置、姿态数据是粗定姿、定轨数据，欧空局的精确定姿、定轨数据难以获取，虽然可经过摄影测量处理提高几何定位精度，但是其改善程度是有限的。

2)本发明试验所用的HRSC相机的几何模型参数不够精确，考虑到火星快车已经在轨运行10年，本发明获取的参数与真实值可能存在少量误差。

3)HRSC发布的Level2既产品数据列宽为5176像元，而欧空局公布的传感器像元数为5184，两者存在8个像元的偏差，对于此问题本发明也进行了试验分析，但是由于误差是各种因素的综合影响，仅将相机模型直接减去或者加上8个像元时结果未见到明显改善。

因此，结合上述分析，如果去除系统误差的影响，本发明的结果与欧空局结果偏差存在1-2个像素，即本发明所提供的基于精确点位预测模型的遥感影像密集匹配方法能够有效解决火星表面影像密集匹配问题，且随着逐层加密匹配过程中同名点的逐渐增多，点位预测精度逐步提高，最终可达到1～3个像素的高精度点位预测结果，在此基础上进行相关系数匹配可快速、精确命中同名像点。

Claims (8)

1.一种基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，该密集匹配方法包括以下步骤：

1)基于严密几何模型对HRSC影像进行近似核线重采样；

2)对重采样后的HRSC影像利用SURF算法提取影像特征点作为初始匹配同名点；

3)将得到的初始匹配同名点作为已知点进行第一次加密匹配，提取待匹配影像的特征点P_i；

4)利用仿射变换模型对所提取的待匹配特征点P_i进行同名点预测，并确定搜索窗口范围；

5)以预测点为中心，在预测点搜索窗口范围内计算各点与待匹配特征点P_i的相关系数，从而确定待匹配特征点P_i的同名点，所述各点为搜索窗口范围内的HRSC影像上的像点；

6)将得到的同名点作为下一层加密匹配的已知点，重复步骤3)至步骤5)的特征点提取、点位预测和相关系数匹配过程，进行逐层加密匹配，直至获取满足需要的匹配点数量，或者影像中已经不能再提取出特征点。

2.根据权利要求1所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤3)是采用Shi-Tomasi算子进行特征点提取的。

3.根据权利要求2所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤4)中的仿射变换模型为：

x₁＝G₀+n₁x₁+n₂y₂

y₁＝b₀+b₁x₁+b₂y₂

其中(x₁，y₁)为待匹配点坐标，(x₂，y₂)是同名点预测坐标，a₀、a₁、a₂、b₀、b₁和b₂为模型参数。

4.根据权利要求3所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤4)中同名点的预测过程为：

a.对待匹配的特征点P_i以设定的搜索半径为步长从已知点中搜索与P_i邻近的n个点；

b.利用得到的n个点按照最小二乘原理解算仿射变换模型中模型参数；

c.根据解算后的模型参数计算仿射变换模型中的同名点预测坐标，从而确定同名点预测值。

5.根据权利要求4所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤4)中搜索窗口的范围由邻近已知点与待匹配特征点P_i的距离以及核线、地形约束条件设定。

6.根据权利要求4所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤a是采用KD树算法对邻近已知点进行搜索。

7.根据权利要求1所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤1)中基于严密几何模型构建过程如下：

A.利用GDAL开源图像处理库将HRSC Level 2级PDS格式影像数据转换为TIFF格式；

B.利用SPICE库获取HRSC相机的焦距、像元大小、像元位置参数；

C.将HRSC原始PDS格式数据导入ISIS系统，导出扫描行时间文件；

D.利用SPICE库及卫星影像扫描行时文件获取每一扫描行的位置、姿态数据；

E.利用获取的相机几何参数、位置、姿态信息构建HRSC影像严密几何模型。

8.根据权利要求1-7中任一项所述的基于精确点位预测模型的遥感影像密集匹配方法，其特征在于，所述步骤6)中当影像中已经不能再提取出特征点时，可在待匹配影像上按照一定的格网间距选择一个特征相对明显的像点作为待匹配点，并对获取的待匹配点进行点位预测、相关系数匹配以获取密集匹配点。