CN103823981B - 一种数字高程模型辅助的卫星影像区域网平差方法 - Google Patents

Abstract

本发明提供一种数字高程模型辅助的卫星影像区域网平差方法。其结合已有的DEM（数字高程模型）数据，在卫星影像以RFM模型（有理函数模型）进行区域网平差时提供高程控制，大大减少了的卫星影像在平差时对稠密分布的地面控制点的依赖。本方法在对弱交会的高分辨率下视（全色和多光谱）卫星影像进行平差时，与传统的非高程控制平差方法相比，其结果精度提升尤为明显。本发明仅使用少量的物方平面高程控制点就可以使一个较大的区域网中的卫星影像获得极高的接边精度和较为均匀的绝对精度。这种方法非常适用于难以大量获取高精度物方控制点的卫星影像用户。

Description

一种数字高程模型辅助的卫星影像区域网平差方法

技术领域

本发明属于测绘科学与技术领域，涉及一种数字高程模型辅助的卫星影像区域网平差方法，主要应用于高分辨率卫星遥感影像的自动空中三角测量、高精度正射影像生产等领域。

背景技术

区域网平差是摄影测量中影像定向参数全局优化的方法。目前最常用的方法是光束法区域网平差。例如：在使用框幅相机的摄影测量生产中，使用严格共线方程模型进行光束法区域网平差，在这一过程中，像点都会通过共线方程模拟摄影中的光束，在平差迭代过程中，当同名点的像点对应的光线以“最小二乘”的标准下达到最优交会时，光束法平差结束。这种流程也可用于卫星遥感影像的全局优化中，卫星影像的RFM模型（有理函数模型）可以使用“像方仿射纠正模型”进行优化，像点通过组成RFM模型的方程来模拟推扫式摄影中的光束，通过区域网平差获得最佳的“像方仿射纠正参数”，从而使连接点的每个光束达到最优交会。但是，卫星遥感影像的几何成像方式同框幅相机航空摄影测量相比，有两点主要差别，一是成像方式，二是影像间基高比较小。由于卫星遥感影像以推扫的方式成像，在进行下视影像的区域网平差时，相邻的航带影像的扫描线之间是近似共面的，如果某一航带的影像中没有控制点，即使使用连接点与其他航带进行连接，它仍然在扫描面上有三个较强的自由度（两个平移自由度和一个旋转自由度），从而导致无法解出稳定的像方仿射纠正参数。在进行非完全下视影像的平差时，虽然能解出稳定结果，但由于卫星影像极小的基高比，其高程误差累积很快。因此，在使用光束法区域网平差手段对卫星影像进行处理时，需要稠密分布的控制点（每景影像9-15个点）。然而，有时大批量的获取高精度的控制点的成本是高昂的，而大量使用精度不高的控制点反而会导致正射影像产品接边精度的下降。

发明内容

本发明的目的是提供一种数字高程模型辅助的卫星影像区域网平差方法。克服在缺乏控制点时，使用光束法区域网平差对卫星影像进行处理时，无法得到可靠结果或结果精度不均匀的问题。

为了解决上述技术问题，本发明提出一种结合DEM（数字高程模型）高程控制的卫星影像区域网平差方法。将连接点在DEM中内插得到的高程值作为控制数据，列出虚拟观测方程，纳入区域网平差中。这种方法可以在控制点稀疏分布（平均每10景影像1个物方控制点）的情况下，使整个区域网达到较为均匀的物方精度和较高的相对精度。因此，该方法在卫星遥感影像的自动空中三角测量，高精度正射影像生成等领域有良好的应用前景。

一种数字高程模型辅助的卫星影像区域网平差方法，包括以下步骤，

步骤1、数据准备；

准备好区域网中影像的连接点匹配结果、控制点物方坐标、控制点像方坐标、对应区域的DEM数据以及将DEM由其本身坐标系转为WGS-84下的经纬度坐标系的参数；

步骤2、获取虚拟高程观测值；

将同名连接点进行关联，结合DEM数据利用“多像迭代反投影方法”，得到连接点的物方平面坐标，并以物方平面坐标对应的DEM高程数据作为该连接点的虚拟高程观测值h；

步骤3、高程观测值定权；

高程观测值的权值的计算公式为：式中σ₀为单位权中误差；σ_h为步骤2中获取的高程观测值的中误差，σ_h取决于DEM本身的精度和区域网中卫星影像的直接定位精度；

步骤4、组建法方程并求解；

通过步骤2和步骤3，连接点都被转化为高程控制点，根据有理函数成像模型和像方仿射纠正模型列出光束法区域网平差的误差方程组并组建法方程组，利用最小二乘法求法方程组，得到成像模型纠正参数改正值的最小二乘解、每个误差方程的残差，然后更新所有模型纠正参数，计算出当前的残差的单位权中误差估计值；

步骤5、迭代进行步骤2至步骤4，直到平差结果的残差中误差的改变量小于给定的阈值，平差结束，得到每张影像的仿射纠正参数和所有连接点的物方坐标。

所述步骤2中的“多像迭代反投影方法”利用不同分辨率的影像上的同名像点具有不同的精度，以得到最优的DEM内插高程值，其具体步骤为：

步骤2.1、求解DEM的平均高程值，作为物方点高程的初始值；

步骤2.2、将物方点高程设为定值，将有理函数成像模型和像方仿射纠正模型针对连接点的物方平面坐标线性化，列出该物方点对应的所有像点的误差方程组，通过最小二乘法迭代求解出物方点当前的平面坐标；

步骤2.3、在步骤2.2求解的物方点平面位置处，内插出DEM的高程，比较这一高程与物方点高程的差绝对值，如果该值大于阈值，则将该高程内插值作为连接点高程，再进行步骤2.2，如果该值小于阈值，则该点迭代结束，内插出的高程即该点“多像迭代反投影”法的结果。

所述的高程观测值的中误差大小σ_h通过以下公式计算：

${\mathrm{\σ}}_h=\sqrt{{\mathrm{\σ}}_H^2+(G\·{\mathrm{\σ}}_P)}$

式中，σ_H为DEM本身的中误差；σ_P为卫星影像的平面“直接定位误差”；G为常数，其意义为卫星影像的平面直接定位误差方向上测区坡度正切值的标准差，基于地面坡度的方向无关性通过统计DEM的x方向的坡度正切值的平方平均数作为G的估值。

所述σ_P的估算方法的步骤为：

步骤3.1、利用步骤2中的“多像迭代反投影法”求出所有平高控制点的平面坐标解算值；

步骤3.2、对每个控制点，求出其平面坐标解算值与其平面坐标观测值的欧氏距离；

步骤3.3、求出这一欧氏距离的平方平均值，该值即卫星影像的平面“直接定位误差”估计值σ_P。

所述的统计DEM的x方向的坡度正切值的平方平均数的计算方法为：

首先：用下式统计出DEM每个格网点上的x方向坡度正切值：

$\tan \mathrm{Gra}{d}_{i,j}=\frac{H_{i+1,j}-H_{i,j}}{X_{i+1,j}-V_{i,j}}$

式中，Grad_i,j为DEM上第j行第i列的格网点在x方向上的坡度，X_i,j为该点x方向的坐标，H_i,j为该点的高程；

然后：用下式得到G的估值，即x方向上的坡度正切值的平方平均：

$G=\sqrt{\frac{1}{m\·n}\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}\tan \mathrm{Gra}{d}_{i,j}}$

式中，m和n分别为DEM的列数和行数。

与现有技术相比，本发明具有以下优点和有益效果：

1、本发明减少了推扫式卫星影像的区域网平差过程中需要的控制点数量，利用了DEM的辅助，大大提高了每个控制点的作用范围。

2、本发明中的“多像最小二乘迭代测图”方法考虑了影像分辨率对迭代结果的影响，得到的结果更为合理。

3、本发明中的“高程观测值定权”步骤，结合影像RFM模型参数误差，给“高程观测值”赋予了合理的权值，保证了在影像初始定向参数精度较低时，平差结果的可靠性。

附图说明

图1为本发明的流程图；

图2为本发明中多像最小二乘迭代测图方法的原理示意图；

图3为本发明中高程观测值定权原理示意图。

具体实施方式

下面结合附图，对本发明作进一步详细的说明。

本发明提供的技术方案是，一种数字高程模型辅助的卫星影像区域网平差方法，如图1所示，包括以下步骤，

步骤1、数据准备；

准备好区域网中影像的连接点匹配结果、控制点物方坐标、控制点像方坐标、对应区域的DEM数据以及将DEM由其本身坐标系转为WGS-84下的经纬度坐标系的参数；

步骤2、获取虚拟高程观测值；

将同名连接点进行关联，结合DEM数据利用“多像迭代反投影方法”，得到连接点的物方平面坐标，并以物方平面坐标对应的DEM高程数据作为该连接点的虚拟高程观测值h；

步骤3、高程观测值定权；

高程观测值的权值的计算公式为：式中σ₀为单位权中误差；σ_h为步骤2中获取的高程观测值的中误差，σ_h取决于DEM本身的精度和区域网中卫星影像的直接定位精度；

步骤4、组建法方程并求解；

通过步骤2和步骤3，连接点都被转化为高程控制点，根据有理函数成像模型和像方仿射纠正模型列出光束法区域网平差的误差方程组并组建法方程组，利用最小二乘法求法方程组，得到成像模型纠正参数改正值的最小二乘解、每个误差方程的残差，然后更新所有模型纠正参数，计算出当前的残差的单位权中误差估计值；

步骤5、迭代进行步骤2至步骤4，直到平差结果的残差中误差的改变量小于给定的阈值，平差结束，得到每张影像的仿射纠正参数和所有连接点的物方坐标。

所述步骤2中的“多像迭代反投影方法”利用不同分辨率的影像上的同名像点具有不同的精度，以得到最优的DEM内插高程值，其具体步骤为：

步骤2.1、求解DEM的平均高程值，作为物方点高程的初始值；

步骤2.2、将物方点高程设为定值，将有理函数成像模型和像方仿射纠正模型针对连接点的物方平面坐标线性化，列出该物方点对应的所有像点的误差方程组，通过最小二乘法迭代求解出物方点当前的平面坐标；

步骤2.3、在步骤2.2求解的物方点平面位置处，内插出DEM的高程，比较这一高程与物方点高程的差绝对值，如果该值大于阈值，则将该高程内插值作为连接点高程，再进行步骤2.2，如果该值小于阈值，则该点迭代结束，内插出的高程即该点“多像迭代反投影”法的结果。

所述的高程观测值的中误差大小σ_h通过以下公式计算：

${\mathrm{\σ}}_h=\sqrt{{\mathrm{\σ}}_H^2+(G\·{\mathrm{\σ}}_P)}$

式中，σ_H为DEM本身的中误差；σ_P为卫星影像的平面“直接定位误差”；G为常数，其意义为卫星影像的平面直接定位误差方向上测区坡度正切值的标准差，基于地面坡度的方向无关性通过统计DEM的x方向的坡度正切值的平方平均数作为G的估值。

所述σ_P的估算方法的步骤为：

步骤3.1、利用步骤2中的“多像迭代反投影法”求出所有平高控制点的平面坐标解算值；

步骤3.2、对每个控制点，求出其平面坐标解算值与其平面坐标观测值的欧氏距离；

步骤3.3、求出这一欧氏距离的平方平均值，该值即卫星影像的平面“直接定位误差”估计值σ_P。

所述的统计DEM的x方向的坡度正切值的平方平均数的计算方法为：

首先：用下式统计出DEM每个格网点上的x方向坡度正切值：

$\tan \mathrm{Gra}{d}_{i,j}=\frac{H_{i+1,j}-H_{i,j}}{X_{i+1,j}-V_{i,j}}$

式中，Grad_i,j为DEM上第j行第i列的格网点在x方向上的坡度，X_i,j为该点x方向的坐标，H_i,j为该点的高程；

然后：用下式得到G的估值，即x方向上的坡度正切值的平方平均：

$G=\sqrt{\frac{1}{m\·n}\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}\tan \mathrm{Gra}{d}_{i,j}}$

式中，m和n分别为DEM的列数和行数。

本发明采用的原理：

物方点(X_i,Y_i,Z_i)卫星影像j上的像点(x_i,j,y_i,j)之间的关系可以用RFM模型进行估算：

$\begin{array}{c}{x}_{i,j}^{\left(0\right)}={\mathrm{Fx}}_j(X_i,Y_i,Z_i)\\ y_{i,j}^{\left(0\right)}={\mathrm{Fy}}_i(X_i,Y_i,Z_i)\end{array}---\left(1\right)$

上式中，Fx和Fy是分母和分子均为三元三次多项式的分式，为使用原始RFM参数解算得到的像方坐标，由于RFM参数存在误差，其值与像点的实际位置(x_i,j,y_i,j)存在差异，以“像方仿射变换模型”为例，它们之间的关系可以表达为：

$\begin{array}{c}{x}_{i,j}={\mathrm{Px}}_j(x_{i,j}^{\left(0\right)},y_{i,j}^{\left(0\right)})=a_{0j}+a_{1j}{x}_{i,j}^{\left(0\right)}+a_{2j}{y}_{i,j}^{\left(0\right)}\\ y_{i,j}={\mathrm{Py}}_j(x_{i,j}^{\left(0\right)},y_{i,j}^{\left(0\right)})=b_{0j}+b_{1j}{x}_{i,j}^{\left(0\right)}+b_{2j}{y}_{i,j}^{\left(0\right)}\end{array}---\left(2\right)$

平差的目的则为解算得到每一张影像a₀,a₁,a₂,b₀,b₁,b₂6个参数，从而实现影像接边和定向的最优化。平差解算的流程见图1.

本方法涉及的“多像最小二乘迭代测图”步骤原理如图2所示。图中，S1和S2是同名点x1和x2拍摄瞬间的摄影中心，r1和r2是同名点x1和x2的量测先验中误差，L1和L2是同名点对应的光线。Le是同名光线的“等效光线”，对任意高程值，迭代得到的物方点均位于“等效光线”上。根据有理函数模型（RFM），像点对应的物方点解算的具体步骤如下：

①设起算高程值H0为物方点的高程（Z＝H₀），如果是第一次多像迭代测图，可以使用测区的平均高作为H0，如果不是第一次多像迭代测图，则使用上一次的高程观测值作为H0。

②根据式(2)，利用最小二乘法求解物方点的平面坐标(X,Y),线性化的方程组如下所示：

$\left\{\begin{array}{c}{x}_j=x_j^{\left(0\right)}+\frac{{\∂\mathrm{Px}}_j}{\∂X}\mathrm{dX}+\frac{{\∂\mathrm{Px}}_j}{\∂Y}\mathrm{dY}\\ y_j=y_j^{\left(0\right)}+\frac{{\∂\mathrm{Py}}_j}{\∂X}\mathrm{dX}+\frac{{\∂\mathrm{Py}}_j}{\∂Y}\mathrm{dY}\end{array}---\left(3\right)\right.$

式中j为同名物方点的像点编号和对应影像的编号,和为式(2)中的像方仿射变换函数和由于该方程是一个非线性方程，需要使用迭代解法。最终解出物方点如图中位置P0。

③在DEM中内插出P0点平面坐标对应的高程值H1，作为物方点的新高程值(Z＝H₁)。

④重复进行步骤②和③，直到第n次解算的物方坐标P_n和第n+1次结果P_n+1之间的距离小于给定的阈值。从而得到“多像迭代测图”的结果，并将解算得到的高程值Hn作为该点（编号设为i）的“高程虚拟观测值”。

本方法涉及的“高程观测值定权”步骤，其原理如图3所示。图中的光线L0是实际光线，Ln是平差之前，相片外方位元素含有误差情况下的光线，如果在“多像迭代测图中”，L0和Ln为“等效光线”。P0是实际光线与实际地面的交点，其坐标为实际物方坐标，Pn是误差光线与DEM表面的交点，其坐标为解算坐标。

为了简化误差模型，将光线近似看做平行直线，并且假设在L0和Ln之间，DEM坡度为常量，DEM高程与实际地面高程的误差为常量。因此，高程误差与平面误差之间的关系为：

$\begin{array}{c}\▿h={\mathrm{\δ}}_H+\▿P\·\tan \mathrm{grad}\\ {\▿P}_0=\▿P+\▿h\·\cot \mathrm{\θ}\end{array}---\left(4\right)$

式中为高程虚拟观测值的误差，δ_H为物方点邻域内DEM的高程误差，为物方点实际位置和通过DEM迭代计算得到的位置之间的平面误差，grad为物方点邻域中DEM在方向上的坡度。为地面为水平面时的物方平面误差，在每一步平差迭代中，区域网中的可以视作一个唯一常量。θ为光线中偏下的方向向量与方向向量的夹角。

对于正下视的卫星影像组成的立体像对来说，等效光线的θ角非常接近直角，可以近似认为cotθ≈0。因此得到高程误差的估算公式：

$\▿h={\mathrm{\δ}}_H+{\▿P}_0\·\mathrm{tam}\mathrm{grad}---\left(5\right)$

如果认为DEM不含有系统误差，仅含有随机误差，则式(5)中，δ_DEM和tangrad均可视作符合期望为正态分布的随机变量。因此根据误差传播定律，得到高程中误差的估计公式：

${\mathrm{\σ}}_h=\sqrt{{\mathrm{\σ}}_H^2+{\▿P}_0^2\·{\mathrm{\σ}}_{\tan \left(\mathrm{grad}\right)}^2}---\left(6\right)$

式(6)中，σ_H为高程观测值的中误差估计值；σ_DEM为DEM数据的高程值中误差，该值可以通过查阅DEM相关的资料和国家标准得到。

σ_tan(grad)为tangrad值的中误差，其值与地形起伏的程度有关，可以根据测区的地形灵活设置，如平原地区设为0.3-0.5，丘陵和山区0.5-1.5。

平面误差可以通过控制点进行估算。如估算测区中所有控制点的DEM迭代计算得到的位置和控制点坐标控制值的平面差值的中误差，将其作为。

因此，在每一步迭代的过程中，可以根据式(6)，估算出高程虚拟观测值的中误差，然后通过单位权中误差计算出权值。

Claims (5)

1.一种数字高程模型辅助的卫星影像区域网平差方法，其特征在于：包括以下步骤，

步骤1、数据准备；

准备好区域网中影像的连接点匹配结果、控制点物方坐标、控制点像方坐标、对应区域的DEM数据以及将DEM由其本身坐标系转为WGS-84下的经纬度坐标系的参数；

步骤2、获取虚拟高程观测值；

将同名连接点进行关联，结合DEM数据利用“多像迭代反投影方法”，得到连接点的物方平面坐标，并以物方平面坐标对应的DEM高程数据作为该连接点的虚拟高程观测值h；

步骤3、高程观测值定权；

高程观测值的权值的计算公式为：式中σ₀为单位权中误差；σ_h为步骤2中获取的高程观测值的中误差，σ_h取决于DEM本身的精度和区域网中卫星影像的直接定位精度；

步骤4、组建法方程并求解；

通过步骤2和步骤3，连接点都被转化为高程控制点，根据有理函数成像模型和像方仿射纠正模型列出光束法区域网平差的误差方程组并组建法方程组，利用最小二乘法求法方程组，得到成像模型纠正参数改正值的最小二乘解、每个误差方程的残差，然后更新所有模型纠正参数，计算出当前的残差的单位权中误差估计值；

步骤5、迭代进行步骤2至步骤4，直到平差结果的残差中误差的改变量小于给定的阈值，平差结束，得到每张影像的仿射纠正参数和所有连接点的物方坐标。

2.根据权利要求1所述的一种数字高程模型辅助的卫星影像区域网平差方法，其特征在于：所述步骤2中的“多像迭代反投影方法”利用不同分辨率的影像上的同名像点具有不同的精度，以得到最优的DEM内插高程值，其具体步骤为：

步骤2.1、求解DEM的平均高程值，作为物方点高程的初始值；

步骤2.2、将物方点高程设为定值，将有理函数成像模型和像方仿射纠正模型针对连接点的物方平面坐标线性化，列出该物方点对应的所有像点的误差方程组，通过最小二乘法迭代求解出物方点当前的平面坐标；

步骤2.3、在步骤2.2求解的物方点平面位置处，内插出DEM的高程，比较这一高程与物方点高程的差绝对值，如果该值大于阈值，则将该高程内插值作为连接点高程，再进行步骤2.2，如果该值小于阈值，则该点迭代结束，内插出的高程即该点“多像迭代反投影方法”的结果。

3.根据权利要求1或2所述的一种数字高程模型辅助的卫星影像区域网平差方法，其特征在于：所述的高程观测值的中误差大小σ_h通过以下公式计算：

${\mathrm{\σ}}_h=\sqrt{{\mathrm{\σ}}_H^2+(G\·{\mathrm{\σ}}_P)}$

式中，σ_H为DEM本身的中误差；σ_P为卫星影像的平面“直接定位误差”；G为常数，其意义为卫星影像的平面直接定位误差方向上测区坡度正切值的标准差，基于地面坡度的方向无关性通过统计DEM的x方向的坡度正切值的平方平均数作为G的估值。

4.根据权利要求3所述的一种数字高程模型辅助的卫星影像区域网平差方法，其特征在于：所述σ_P的估算方法的步骤为：

步骤3.1、利用步骤2中的“多像迭代反投影方法”求出所有平高控制点的平面坐标解算值；

步骤3.2、对每个控制点，求出其平面坐标解算值与其平面坐标观测值的欧氏距离；

步骤3.3、求出这一欧氏距离的平方平均值，该值即卫星影像的平面“直接定位误差”估计值σ_P。

5.根据权利要求3所述的一种数字高程模型辅助的卫星影像区域网平差方法，其特征在于：所述的统计DEM的x方向的坡度正切值的平方平均数的计算方法为：

首先：用下式统计出DEM每个格网点上的x方向坡度正切值：

${\mathrm{tanGrad}}_{i,j}=\frac{H_{i+1,j}-H_{i,j}}{X_{i+1,j}-X_{i,j}}$

式中，Grad_i,j为DEM上第j行第i列的格网点在x方向上的坡度，X_i,j为该点x方向的坐标，H_i,j为该点的高程；

然后：用下式得到G的估值，即x方向上的坡度正切值的平方平均：

$G=\sqrt{\frac{1}{m\·n}\underset{i=1}{\overset{m-1}{\Σ}}\underset{j=0}{\overset{n-1}{\Σ}}{\mathrm{tanGrad}}_{i,j}}$

式中，m和n分别为DEM的列数和行数。