CN110826407B - 一种高分辨率卫星广义像对的立体匹配方法 - Google Patents

Abstract

该发明公开了一种高分辨率卫星广义像对的立体匹配方法，属于高分辨率卫星广义像对的立体匹配方法，利用高分辨率卫星广义像对进行立体匹配。本发明用图像处理算法进行特征提取和特征点匹配，摆脱了传统遥感软件需要大量手动修改校准的过程，即可以提高精度，也可以节约时间成本。同时本方法适用于广义立体像对，将得到的同名点代入有理函数前方交会模型，得出地面点的三维坐标。

Description

一种高分辨率卫星广义像对的立体匹配方法

技术领域

本发明属于高分辨率卫星广义像对的立体匹配方法，利用高分辨率卫星广义像对进行立体匹配。

背景技术

随着遥感和航天技术的高速发展，高分辨卫星影像成为获取地球表面信息的重要工具。随着卫星数量和分辨率的提高，利用高分辨卫星影像获取数字高程模型(DigitalElevation Model,DEM) 或者对地面目标进行三维重建成为可能。但是由于政治、成本等因素影响，利用严格立体像对去获取这些信息面临着巨大的困境。如果能利用具有重叠区域的影像，但又不考虑影像的来源，分辨率等因素，再利用遥感原理和图像处理算法构建广义上的遥感像对，就可以克服传统方法的弊端。

发明内容

本发明提出了一种高分辨率卫星广义像对的立体匹配方法。该方法使用有理函数模型进行影像定位，用图像处理算法进行特征提取和特征点匹配，将得到的同名点代入有理函数前方交会模型，得出地面点三维坐标，并对有理函数模型进行误差补偿，得出合适的三维坐标。

本发明技术方案为一种高分辨率卫星广义像对的立体匹配方法，该方法包括：

步骤1：选取合适的不同来源的遥感广义立体像对，并对其进行预处理；

选取的两幅遥感影像数据作为广义立体像对需要满足以下条件：(1)两幅影像有相同的覆盖区域；(2)两幅影像的获取时间基本一致，影像覆盖区域基本没有变化；(3)影像分辨率较好； (4)影像云量和积雪较少；然后对获取的影像进行线性拉伸或直方图均衡化；

步骤2：对获取的高分辨率广义立体像对的左右影像利用有理函数模型进行影像定位；

步骤3：利用SURF提取算法对左右影像进行特征点提取，然后对两幅影像的特征点进行匹配得到同名点对；

步骤4：利用左右影像同名点对，通过有理函数前方交会模型获取高程。

进一步的，所述步骤2的具体方法为：

由于本发明使用的高分辨卫星数据都是商业卫星，商家不提供成像时传感器参数，因此只能建立步骤1预处理后的广义立体像对之间的有理函数模型，为后续的定位，计算高程；

所述有理函数模型(rational function model,RFM)描述了像点坐标(r,c)和地面点三维坐标(X,Y,Z)一一对应关系，如式1：

其中，p₁、p₂、p₃、p₄为RFM的四个多项式函数；(r_n,c_n)和(X_n,Y_n,Z_n)分别表示像点坐标 (r,c)和地面点三维坐标(X,Y,Z)经平移和缩放后的归一化坐标，取值介于-1.0～+1.0之间；变换公式为：

其中(X₀,Y₀,Z₀,r₀,c₀)为标准化平移参数，(X_S,Y_S,Z_S,r_S,c_S)为标准化比例参数；RFM采用归一化坐标是为了减少计算过程中由于数据数量级差别过大引入的舍入误差；每种卫星影像中包含了RPC参数：行平移参数(Line Offset)，列平移参数(Sample Offset)，纬度平移参数(Latitude Offset)，经度平移参数(Longitude Offset)，高程平移参数(Height Offset)，行比例参数(Line Scale)，列比例参数(Sample Scale)，纬度比例参数(Latitude Scale)，经度比例参数(Longitude Scale)，高程比例参数(Height Scale)；

p₁、p₂、p₃、p₄多项式的形式为式(4)：

其中，a₀、a₁、a₂.......a₁₇、a₁₈、a₁₉为RPC参数中的多项式计算系数，(X,Y,Z)为地面点三维坐标，将像点坐标对代入有理函数模型，其中RPC参数由卫星影像提供，计算出这些像点的大地坐标。

进一步的，所述步骤3的具体方法为：

利用图像处理算法对立体像对的左右两幅影像进行特征点提取和匹配，使用快速鲁棒特征 SURF(Speeded-Up-Robust Features)算法提取到两幅图像的特征点，然后利用GMS算法精匹配不同高分辨卫星影像；

步骤3.1：将积分图像中的任意一点的值定义为从图像的左上角到这个点的所构成的矩形区域内所有的像素点的灰度值之和，将复杂的计算转化为几个矩形之间的四则运算，减少积分图像进行运算的时间；

步骤3.2：确定Box滤波器，并设定Box滤波器窗口的大小；

步骤3.3：确定极值点；

通过每幅图像与box滤波器进行多次卷积，生成的多幅图像构成金字塔；选取当前点的图像，以及上下相邻的图像，然后将周围8个点，上下各9个点与自身比较，如果上述26个点与自身的不等号恒成立，该点默认为特征点；

步骤3.4：对特征点进行数学量化，为其分配方向；

步骤3.4.1：对特征点所在金字塔那一层的图像进行Haar小波变换，得到一个圆心角为60°的扇形区域；小波周长取4s，s为特征点所在的尺度值；将得到的扇形区域绕着同一方向旋转，会得到6个不同的方向，对6个不同的方向做矢量加法；在得到的6个向量中，将模值最大的扇形区域里的矢量方向当作是特征点的方向；

步骤3.4.2：生成特征描述符；以特征点为中心，选取特征点所在金字塔那一层图像的20s×20s 矩形，再把这个矩形平均分成4×4＝16份，得到面积为5s×5s的小区域；然后对小区域进行Haar 小波变换，s为大于1的整数；

步骤3.5：匹配同名点；使用BF匹配进行初步匹配，然后，对粗匹配得到的像素点对使用GMS 算法进行精匹配。

进一步的，所述步骤4的具体方法为：利用左右相片的RFM模型，再根据同名点开始求解地面点的三维坐标，推导过程如下：

则有

f'_X(X)、f'_Y(Y)、f'_Z(Z)为对应函数的倒数，将

代入到公式(1)中，整理得

令

则公式(9)可以改写成：

将公式(10)按照泰勒定理进行线性化，变成只有一次项的结构：

其中：

表示公式(10)按公式(11)格式利用泰勒定理进行线性化之后的余项之和，ΔX 表示自变量X的增量，ΔY表示自变量Y的增量，ΔZ表示自变量Z的增量；

由公式(11)可以得出误差方程，如公式(12)所示：

其中v_r和v_c表示行和列的改正数，在计算时一般设为0，得到三个增量ΔX、ΔY、ΔZ，加到原来的值上面，用来重新计算，减小误差；

以上各式中的偏导数的表达式为：

根据左右相片的同名点坐标(r_l,c_l)和(r_r,c_r)，可以列出以下四个误差方程：

上式可以改写为：

V＝AΔ-l (23)

其中：

Δ的最小二乘解为:

Δ＝[ΔZ ΔY ΔX]^T＝(A^TA)^-1A^Tl (24)

地面点三维坐标(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)作为初始值采用最小二乘法对V＝AΔ-l进行迭代计算得到对应点的三维坐标；

(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)这样得到：将两幅影像中提供的RPC参数文件的给定参数平均值作为起始坐标 (X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)，即

本发明用图像处理算法进行特征提取和特征点匹配，摆脱了传统遥感软件需要大量手动修改校准的过程，即可以提高精度，也可以节约时间成本。同时本方法适用于广义立体像对，将得到的同名点代入有理函数前方交会模型，得出地面点的三维坐标。目前，直接购买处理好的立体像对价格是十分昂贵的，本发明在一定程度上可以节约资金。

附图说明

图1为高分一号左右影像图；

图2为高分二号左右影像图；

图3为Pleiades左右影像图；

图4为QuickBird和WorldView-2左右影像图；

图5为WorldView-3左右影像图；

图6为地面点三维坐标计算流程图；

图7为高分一号像对实验结果图；(a)为高分一号影像DEM，(b)为高分一号研究区域三维散点图；

图8为高分二号像对实验结果图；(a)为高分二号影像DEM，(b)为高分二号研究区域三维散点图；

图9为Pleiades像对实验结果图；(a)为Pleiades影像DEM，(b)为Pleiades研究区域三维散点图；

图10为WorldView-3立体像对实验结果图；(a)为WorldView-3影像DEM，(b)为WorldView-3 研究区域三维散点图；

图11为未纠正的RFM模型得到的高程误差图；

图12为纠正的RFM模型得到的高程误差图；

图13为高分一号像对三维渲染结果图；

图14为高分二号像对三维渲染结果图；

图15为Pleiades像对三维渲染结果图；

图16为WorldView-3像对三维渲染结果图；

图17为本发明总流程图；

具体实施方式

本发明使用了7种不同来源的影像：高分一号，高分二号，Pleiades 1A，Pleiades1B，QuickBird， WorldView-2，WorldView-3左，WorldView-3右，进行特征点的提取、匹配，然后利用有理函数前方交会模型计算高程，并绘制DEM。本发明实验流程如图17所示。

1)选择遥感图像

按照广义立体像对的要求选取遥感图像，图1-5为所选取的高分一号，高分二号，Pleiades 1A， Pleiades 1B，QuickBird，WorldView-2，WorldView-3左，WorldView-3右立体像对。

2)有理函数模型定位精度分析

本发明所使用影像数据RPC归一化参数，如表1和2：

表1前四种RPC归一化参数

表2后四种RPC归一化参数

根据遥感图像信息，可以得到有理函数模型参数，本文采用多种传感器卫星影像，将像点坐标代入有理函数模型，其中RPC参数卫星影像已经提供，然后计算出这些像点的大地坐标，并事先在Google Earth上选取对应的点作为检查点，计算均方误差和各向最大误差。表3-12为各个卫星影像RFM系数拟合误差。

表3高分1号左卫星影像RFM系数拟合误差

表4高分一号右卫星影像RFM系数拟合误差

表5高分二号A卫星影像RFM系数拟合误差

表6高分二号B卫星影像RFM系数拟合误差

表7 Pleiades 1A卫星影像RFM系数拟合误差

表8 Pleiades 1B卫星影像RFM系数拟合误差

表9 QuickBird卫星影像RFM系数拟合误差

表10 WorldView-2卫星影像RFM系数拟合误差

表11 WorldView-3左卫星影像RFM系数拟合误差

表12 WorldView-3右卫星影像RFM系数拟合误差

使用三阶多项式形式的有理函数模型去拟合传感器成像模型，能得到很好的精度。首先不管是各向残差还是方差的数量级基本上在小数点后三位以上。容易观察到，分母不相等且不为1 的时候，为有理函数模型最完整的形式，精度最高。其次，分母相等但是不为1时，解算模型时可能会产生很强的相关性，降低模型精度，精度要低于分母不相等这种情况。最后一种情况是分母都为1，这时有理函数模型退化成多项式模型，精度也是三种情形里最低的。

另外，从10张表中可以看出，从高分一号、高分二号、Pleiades、QuickBird、WorldView-2到WorldView-3影像的分辨率越来越高，相应的影像RFM系数拟合误差越来越小，但是其中高分一二号影像的拟合误差较大的原因还包括影像对应的地面区域地形起伏较大，后面几种影像数据拍摄的都是平坦的市区，这也证明了上面的分析结果是正确的。

3)特征提取和匹配

采用SURF特征提取算法提取广义立体像对中的特征点，完成两幅遥感影像的粗匹配。SURF 算法发现的任何一个关键点，都可以把它所在区域划分为4个小的区域，形成的16个区域的方向和Harr小波响应组成一个独特的64维特征向量。利用SURF算法对特征的描述进行粗匹配，可以得到粗略的像素点对。

然后，对粗匹配得到的这些像素点对使用GMS算法进行精匹配。选取最后匹配较好的同名点对，进行地面点三维坐标的获取。

4)地面点三维坐标获取

本发明利用左右相片的RFM模型，再根据同名点开始求解地面点的三维坐标，由于QuickBird 和WorldView-2提取的同名点过少，不参与后续计算。

计算步骤如下：

(1)计算地面坐标初始值(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)，并分别用左右影像的标准化参数将其转化为对应的标准化坐标

及

(2)分别利用

和

计算式中各个偏导数以及

和

并根据同名点坐标(r_l,c_l)，(r_r,c_r)组成误差方程；

(3)对误差方程并求解，得到坐标改正数ΔZ⁽ⁱ⁾,ΔY⁽ⁱ⁾,ΔX⁽ⁱ⁾(i＝0,1,...)。如果改正数超出容许范围，则修正当前的地面坐标值，令

X⁽ⁱ⁺¹⁾＝X⁽ⁱ⁾+ΔX⁽ⁱ⁾，Y⁽ⁱ⁺¹⁾＝Y⁽ⁱ⁾+ΔY⁽ⁱ⁾，Z⁽ⁱ⁺¹⁾＝Z⁽ⁱ⁾+ΔZ⁽ⁱ⁾，并计算X⁽ⁱ⁺¹⁾，Y⁽ⁱ⁺¹⁾，Z^{(i +1)}的左右片标准化坐标

和

转回第(2)步迭代进行。否则跳出循环，最终的X⁽ⁱ⁺¹⁾，Y⁽ⁱ⁺¹⁾，Z⁽ⁱ⁺¹⁾即为地面点空间坐标。

图6为地面点三维坐标计算流程。将高分一号像对、高分二号像对、Pleiades像对、QuickBird 和WorldView-2像对、WorldView-3立体像对求得特征点和同名点，再用上述方法求得对应的地面点高程，生成如下各图的三维散点图，利用这些信息可以求取DEM，各区域(图中矩形区域)的经纬度分别为：高分一号为103.59°～103.6°W，37.81°～37.82°N。高分二号为103.2°～103.21°W，31.425°～31.435°N。Pleiades为102.95°～102.96°W，31.5°～31.51°N，WorldView-3为43.17°～ 43.18°E，22.95°～22.97°S。图7-10为各广义立体像对实验结果。

5)RFM模型精度补偿由于没有实际的地面控制点，测量到的控制点计算的RPC参数精度不够高。物方补偿模型认为RFM模型定位的系统误差主要是由其物方坐标系与控制点坐标系不一致造成的。本节采用 Pleiades影像对对RFM模型的精度特性进行了验证，在Google Earth上查询对应点的高程用以精度验证。然后引进改正模型补偿RFM模型系统误差。

首先，利用RFM模型计算出各点的地面坐标，并统计出各点的定位误差。结果表明，RFM模型定位的最大误差分别为：高分一号像对10m，高分二号像对9m，Pleiades像对6.7m，WorldView-3立体像对7.2m。图11为未纠正的RFM模型得到的高程误差。

然后引进改正模型补偿RFM模型系统误差，物方补偿模型认为RFM模型定位的系统误差主要是由其物方坐标系与控制点坐标系不一致造成的，这种不一致可通过引入空间相似变换参数进行改正，认为两个坐标系之间存在空间的平移旋转和缩放，即:

其中，λ为尺度常量，M为旋转矩阵，由两坐标系间轴线偏差角(Φ,Ω,K)构成，(X₀,Y₀,Z₀)为平移参数。空间相似变换中共有7个独立参数，至少需要3个地面控制点才能解算。实际上，控制点坐标系和RFM模型的物方坐标系吻合是很好的，因此尺度参数一般为1，小范围情况下旋转矩阵M一般认为是单位矩阵，这时只需要一个平高控制点就可以确定参数。改正后结果如图12。

改正后高程最大误差为：高分一号像对5.2m，高分二号像对5.4m，Pleiades像对4.2m， WorldView-3立体像对5.2m。对每种像对随机取100个同名点计算出的高程，与Google Earth 上面的对应点对比，高程差在3米以下的点分别为83、78、76、85，可以认为本文求出的高程精度为83％、78％、76％、85％。图13-16为使用本文得到的DEM和原始图像进行实际地形渲染结果图，其中红色和蓝色矩形框为本文截取的实验影像区域。

Claims (2)

1.一种高分辨率卫星广义像对的立体匹配方法，该方法包括：

步骤1：选取合适的不同来源的遥感广义立体像对，并对其进行预处理；

选取的两幅遥感影像数据作为广义立体像对需要满足以下条件：(1)两幅影像有相同的覆盖区域；(2)两幅影像的获取时间基本一致，影像覆盖区域基本没有变化；(3)影像分辨率较好；(4)影像云量和积雪较少；然后对获取的影像进行线性拉伸或直方图均衡化；

步骤2：对获取的高分辨率广义立体像对的左右影像利用有理函数模型进行影像定位；

步骤3：利用SURF提取算法对左右影像进行特征点提取，然后对两幅影像的特征点进行匹配得到同名点对；

步骤4：利用左右影像同名点对，通过有理函数前方交会模型获取高程；

所述步骤2的具体方法为：

建立步骤1预处理后的广义立体像对之间的有理函数模型，为后续的定位，计算高程；

所述有理函数模型(rational function model,RFM)描述了像点坐标(r,c)和地面点三维坐标(X,Y,Z)一一对应关系，如式1：

其中，p₁、p₂、p₃、p₄为RFM的四个多项式函数；(r_n,c_n)和(X_n,Y_n,Z_n)分别表示像点坐标(r,c)和地面点三维坐标(X,Y,Z)经平移和缩放后的归一化坐标，取值介于-1.0～+1.0之间；变换公式为：

其中(X₀,Y₀,Z₀,r₀,c₀)为标准化平移参数，(X_S,Y_S,Z_S,r_S,c_S)为标准化比例参数；每种卫星影像中包含了RPC参数：行平移参数(Line Offset)，列平移参数(Sample Offset)，纬度平移参数(Latitude Offset)，经度平移参数(Longitude Offset)，高程平移参数(HeightOffset)，行比例参数(Line Scale)，列比例参数(Sample Scale)，纬度比例参数(LatitudeScale)，经度比例参数(Longitude Scale)，高程比例参数(Height Scale)；

p₁、p₂、p₃、p₄多项式的形式为式(4)：

其中，a₀、a₁、a₂.......a₁₇、a₁₈、a₁₉为RPC参数中的多项式计算系数，(X,Y,Z)为地面点三维坐标，将像点坐标对代入有理函数模型，其中RPC参数由卫星影像提供，计算出这些像点的大地坐标；

所述步骤4的具体方法为：利用左右相片的RFM模型，再根据同名点开始求解地面点的三维坐标，推导过程如下：

则有

f′_X(X)、f′_Y(Y)、f′_Z(Z)为对应函数的倒数，将

代入到公式(1)中，整理得

令

则公式(9)可以改写成：

将公式(10)按照泰勒定理进行线性化，变成只有一次项的结构：

其中：

表示公式(10)按公式(11)格式利用泰勒定理进行线性化之后的余项之和，ΔX表示自变量X的增量，ΔY表示自变量Y的增量，ΔZ表示自变量Z的增量；

由公式(11)可以得出误差方程，如公式(12)所示：

其中v_r和v_c表示行和列的改正数，在计算时一般设为0，得到三个增量ΔX、ΔY、ΔZ，加到原来的值上面，用来重新计算，减小误差；

以上各式中的偏导数的表达式为：

根据左右相片的同名点坐标(r_l,c_l)和(r_r,c_r)，可以列出以下四个误差方程：

上式可以改写为：

V＝AΔ-l (23)

其中：

Δ的最小二乘解为:

Δ＝[ΔZ ΔY ΔX]^T＝(A^TA)^-1A^Tl (24)

地面点三维坐标(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)作为初始值采用最小二乘法对V＝AΔ-l进行迭代计算得到对应点的三维坐标；

(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)这样得到：将两幅影像中提供的RPC参数文件的给定参数平均值作为起始坐标(X⁽⁰⁾,Y⁽⁰⁾,Z⁽⁰⁾)，即

2.如权利要求1所述的一种高分辨率卫星广义像对的立体匹配方法，其特征在于所述步骤3的具体方法为：

利用图像处理算法对立体像对的左右两幅影像进行特征点提取和匹配，使用快速鲁棒特征SURF(Speeded-Up-Robust Features)算法提取到两幅图像的特征点，然后利用GMS算法精匹配不同高分辨卫星影像；

步骤3.1：将积分图像中的任意一点的值定义为从图像的左上角到这个点的所构成的矩形区域内所有的像素点的灰度值之和，将复杂的计算转化为几个矩形之间的四则运算，减少积分图像进行运算的时间；

步骤3.2：确定Box滤波器，并设定Box滤波器窗口的大小；

步骤3.3：确定极值点；

通过每幅图像与box滤波器进行多次卷积，生成的多幅图像构成金字塔；选取当前点的图像，以及上下相邻的图像，然后将周围8个点，上下各9个点与自身比较，如果上述26个点与自身的不等号恒成立，该点默认为特征点；

步骤3.4：对特征点进行数学量化，为其分配方向；

步骤3.4.1：对特征点所在金字塔那一层的图像进行Haar小波变换，得到一个圆心角为60°的扇形区域；小波周长取4s，s为特征点所在的尺度值；将得到的扇形区域绕着同一方向旋转，会得到6个不同的方向，对6个不同的方向做矢量加法；在得到的6个向量中，将模值最大的扇形区域里的矢量方向当作是特征点的方向；

步骤3.4.2：生成特征描述符；以特征点为中心，选取特征点所在金字塔那一层图像的20s×20s矩形，再把这个矩形平均分成4×4＝16份，得到面积为5s×5s的小区域；然后对小区域进行Haar小波变换，s为大于1的整数；

步骤3.5：匹配同名点；使用BF匹配进行初步匹配，然后，对粗匹配得到的像素点对使用GMS算法进行精匹配。

Images