CN111354040B - 一种基于Partial EIV模型的光学卫星影像区域网平差方法 - Google Patents

Abstract

本发明公开了一种基于Partial EIV模型的光学卫星影像区域网平差方法，将Partial EIV模型应用于光学卫星影像区域网平差，并推导基于该模型的影像几何定位方法。首先构建光学卫星影像区域网平差模型，根据Partial EIV模型构建平差目标函数，然后依据Lagrange条件极值原理，推导未知参数的混合总体最小二乘估计方法，最后利用资源三号光学卫星影像进行试验与分析。

Description

一种基于Partial EIV模型的光学卫星影像区域网平差方法

技术领域

本发明涉及遥感影像处理领域，具体为一种基于Partial EIV模型的光学卫星影像区域网平差方法。

背景技术

近年来，随着国产立体测绘卫星“天绘一号”[1-2]和“资源三号”[3]的成功发射，提高国产光学卫星影像的几何定位精度成为摄影测量工作者的研究热点。根据其成像模型的不同可分为基于严格几何成像模型[4-7]和基于有理函数模型(RFM)[8-11]。有理函数模型数学形式简单，使用方便，具有一般性和保密性，已广泛应用于光学卫星影像几何处理。文献[8]于2003年首次将附加像方补偿的RFM模型应用于光学卫星影像几何处理。文献[12]利用有理函数模型对SPOT5-HRS影像进行区域网平差实验，结果表明利用少量地面控制点可以较好的补偿RPC参数中的系统误差。文献[13]针对RPC参数受到系统误差的影响，采用基于像方改正的RPC定位方法，实现了系统误差改正与几何定位精度的提升。文献[14]研究了RPC系统误差改正方法及其参数重新生成技术，并利用QuickBird卫星影像进行了实验分析。文献[15]采用最小二乘配置理论，将像方的系统误差分为随机与非随机两部分，该方法具有更高的定位精度与稳定性。文献[16]使用RFM模型对覆盖全国的26 000多景资源三号立体影像进行无控区域网平差，平面和高程中误差均达到4m以内。

虽然以上方法在一些实际应用中都取得了较好的定位结果，但大都采用最小二乘方法进行未知参数估计，其前提是认为误差函数模型为Gauss-Markov(GM)模型，即假定误差的函数模型已知和非随机性，且观测方程中仅包含随机误差[17]。而该假设在光学卫星影像区域网平差中并不严格成立，因为观测方程用照泰勒公式展开成线性形式后，在系统误差未知参数对应的系数矩阵中，不但含有非随机量，而且还含有像点观测值等随机量。尽管在大多实际工程应用中采用最小二乘方法进行未知参数估计可以解决问题，但是其在理论上并不充分严密。

变量误差模型(Errors-in-variables，EIV)不但顾及了观测向量的随机误差，而且还考虑到系数矩阵也可能包含随机误差，采用总体最小二乘方法(Total leastsquares，TLS)，可以较好的估计未知参数[18-19]。但在实际应用中[20]，许多情况下系数矩阵只有部分含有随机误差，如线性回归模型和坐标转换问题，这类情况可以采用PartialEIV模型[21-23]进行未知参数估计，相比较于EIV模型，该模型更适合于解决实际应用问题。文献[21]推导了Partial EIV模型的混合总体最小二乘方法(Mixed total leastsquares，M-TLS)，但是该方法存在收敛速度较慢、迭代次数过多等问题。文献[24]在此基础上推导了适合数据量较大的M-TLS方法。文献[25]进一步利用Partial EIV模型解决了附有不等式约束的EIV模型问题。Partial EIV模型是EIV模型的更为通用的表达式，只将系数矩阵中随机元素部分提取出来，并作为未知参数进行估计，提高了整体计算效率。

发明内容

目前，Partial EIV模型在摄影测量方面的研究成果相当有限，本发明将PartialEIV模型应用于光学卫星影像区域网平差，并推导基于该模型的影像几何定位方法。首先构建光学卫星影像区域网平差模型，根据Partial EIV模型构建平差目标函数，然后依据Lagrange条件极值原理，推导未知参数的混合总体最小二乘估计方法，最后利用资源三号光学卫星影像进行试验与分析。

为实现上述目的，本发明提供如下技术方案：一种基于Partial EIV模型的光学卫星影像区域网平差方法，包括以下步骤：

步骤1，根据影像连接点和对应地面点的关系建立有理函数模型；

步骤2，在有理函数模型的基础上，构建区域网平差模型，解算整个测区中所有影像的仿射变换系数和所有加密点的地面坐标；

步骤3，基于P-EIV模型补偿平差方程中的观测误差和系数矩阵误差，利用Lagrange条件极值原理，推导基于混合总体最小二乘的光束法区域网平差计算方法。

与现有技术相比，本发明提供了一种基于Partial EIV模型的光学卫星影像区域网平差方法，具备以下有益效果：

与现有技术相比，本发明在光学卫星影像区域网平差中，顾及像点观测方程的系数矩阵中存在随机误差和非随机误差，提出了基于Partial EIV模型的混合总体最小二乘的光学卫星影像定位方法。与传统最小二乘方法相比，Partial EIV模型提高了定位精度。

附图说明

图1为本发明的流程框架图；

图2为法国地区地面控制点分布图；

图3为为本发明试验区域控制但布设方案；

图4为本发明法国地区方案1检查点误差分布(EIV模型)；

图5(a)-(c)为法国地区方案2-4检查点误差分布；

表1为法国地区ZY-3影像各方案平差精度。

具体实施方式

参照图1、2、3、5，本发明提供一种技术方案：一种基于Partial EIV模型的光学卫星影像区域网平差方法，包括以下步骤：

步骤1，根据影像连接点和对应地面点的关系建立有理函数模型；

步骤2，在有理函数模型的基础上，构建区域网平差模型，解算整个测区中所有影像的仿射变换系数和所有加密点的地面坐标；

步骤3，基于P-EIV模型补偿平差方程中的观测误差和系数矩阵误差，利用Lagrange条件极值原理，推导基于混合总体最小二乘的光束法区域网平差计算方法。

1.本实施例中，所述步骤1，根据影像连接点和对应地面点的关系建立有理函数模型，该有理函数模型RFM是将像方影像坐标(r,c)表示为以对应物方地面三维空间坐标(P,L,H)为自变量的多项式的比值，其本质是对严密模型的高度数学拟合，可以实现与其一样的几何定位精度；

(1)其形式如下所示：

式中，(r_n,c_n)(P_n,L_n,H_n)分别表示像点坐标和地面点坐标正则化后的坐标，其值在(-1～+1)之间；正则化可以提高RPC模型中各系数求解的稳定性，同时有助于减少计算过程中由于数据级数差别过大而引起的数据舍入误差；

(2)正则化是通过平移和缩放实现的，正则化公式如下所示：

式中，r₀,c₀,P₀,L₀,H₀为正则化的平移参数；r_s,c_s,P_s,L_s,H_s为正则化的缩放系数；

(3)利用附加像方补偿方法对有理函数模型的系统误差进行改正，由式(1)可得：

式中，分别表示光学卫星影像行、列方向系统误差改正数，(a₀a₁,a₂,b₀,b₁,b₂)表示系统误差的补偿系数。

本实施例中，基于RFM模型的区域网平差是利用少量地面控制点解算整个测区中所有影像的仿射变换系数和所有加密点的地面坐标，是摄影测量应用中重要的步骤；光学卫星影像的RPC参数存在较大的系统误差，通过影像之间连接点的约束关系对RFM模型的系统误差进行补偿，进而提高目标几何定位精度；

(4)根据式(3)进一步写成如下：

式中，

(5)对式(4)进行线性化后，

进一步得到基于RFM模型的区域网平差模型的误差方程的矩阵形式：

V＝At+Bx-L,P

式中，V＝[v_r v_c]^T表示像点行、列坐标残差向量； 分别表示相应的系数矩阵；x＝[ΔP ΔL ΔH]^T表示连接点对应地面点的空间大地坐标改正向量；t表示系统误差补偿模型参数改正向量；L表示表示观测值向量；P表示相应的权阵；

(6)注意到系数矩阵A由像点观测值(r,c)组成，因而矩阵A中包含观测误差，得到：

通常情况下，地面点坐标未知数x的个数要远大于系统误差补偿未知数t的个数，消去x后可得t的解为：

t＝[A^TPA-A^TPB(B^TPB)^-1(B^TPA)]^-1·[A^TPL-A^TPB(B^TPB)^-1(B^TPL)]

区域网平差是一个迭代逐步趋近的过程，当两次计算的结果小于给定的阈值时，迭代计算过程结束。

本实施例中，Partial EIV模型总体最小二乘方法是在传统最小二乘的基础上进一步发展得到的平差方法，是EIV模型的更为严密的估计方法，EIV模型同时顾及了观测向量与系数矩中包含的随机误差；

(7)其数学模型表示为：

L＝(A-E_A)X+e_L

(8)其随机模型为：

式中，vec(·)表示矩阵的拉直变换，表示单位权中误差；

(9)在基于RFM模型的区域网平差方法中，系数矩阵A的元素并非全都是随机的，部分是由0或者1等常数组成。针对EIV模型中系数矩阵由一些非随机元素与随机元素组成，对系数矩阵A进一步处理，将系数矩阵中的随机元素提取出来，进而EIV模型可以写成PartialEIV模型，其数学模型表示为：

(10)其随机模型为：

式中，表示n×m系数矩阵真值，h表示一个确定的(n×m)×1常数向量，由中非随机元素和零组成，B表示一个给定的(n×m)×s矩阵，s表示系数矩阵中随机元素的个数，L表示n×1的观测向量，e_L表示观测向量误差；a表示系数矩阵中随机元素组成的s×1列向量，表示其真值，e_a表示a的随机误差，Q_L表示观测量误差的协因数阵，Q_a表示系数矩阵误差的协因数阵，表示单位权方差；I_n表示n×n单位矩阵，表示Kronecker积；

(11)进一步构造平差目标函数：

(12)对式(11)求偏导，整理后得到参数估计为：

式中，

(13)改正数估计：

式中，

(14)单位权中误差估计：

迭代计算步骤如下：

1)给定向量A,Q_L,Q_a,L,h,B,a，利用最小二乘计算初值

2)根据式(12)计算和

3)重复步骤2)直至小于给定的阈值，此处设为1e-10，输出 和

但是针对待解未知数较多的情况下，采用上述方法求解，效率并不高。本发明采用一种高效的Partial EIV模型解法[22]。

(15)根据Lagrange条件极值定理，对式(9)进行重新构造目标函数如下：

式中，λ和u均为Lagrange因子；

(16)为使式(15)最小，根据条件极值定理可得到方程组：

(17)由式(16)得到：

(18)

(19)

(20)

(21)

式中，

(22)

式中，

迭代计算步骤如下：

1)给定向量A,Q_L,Q_a,L,h,B,a，利用最小二乘计算初值

2)根据式(20)和式(21)计算e_a和e_L；

3)计算和

4)再根据式(22)计算

5)重复步骤(2)-(4)直至小于给定的阈值，此处设为1e-10，输出

利用上述Partial EIV模型对RPC的系统误差的改正参数进行估计，并对像点观测量进行更新及改正之后，回代式(6)求解地面坐标。

具体实施例，为验证P-EIV模型应用于光学卫星影像区域网平差的有效性，采用资源三号卫星影像对该方法进行验证。每组试验主要分为三部分。

(1)无地面控制点直接定位，分别采用最小二乘(GM模型)和总体最小二乘(EIV模型)方法进行实验验证。地面控制点作为检查点进行精度评价。

(2)稀少地面控制点参与平差，分别采用最小二乘(GM模型)、总体最小二乘(EIV模型)和混合总体最小二乘(P-EIV模型)方法进行进行平差实验。没有参与平差的剩余地面控制点作为检查点。

(3)大量地面控制点参与平差，分别采用上述三种方法进行平差实验。

最后在进行定位精度评价之前，需要先将检查点的三维空间坐标转换到空间直角坐标系下，再进行比较与分析。

实验数据：试验数据采用法国Sainte-Maxime地区的资源三号卫星影像，该影像成像于2014年，影像范围内有12个地面控制点(如图2所示)。像方及物方坐标均由国际摄影测量与遥感协会ISPRS(International Society of Photogrammetry and Remote Sensing)提供。

实验方案：根据上述基本思路，设计了4种试验方案，如图3所示。对不同地面控制点布设情况下利用最小二乘(GM模型)、总体最小二乘(EIV模型)和混合总体最小二乘(P-EIV模型)方法进行区域网平差实验，并对实验结果进行了比较与分析。

结果与分析：利用上述ZY-3卫星影像对4种方案分别进行试验并得出相应的结果，试验中各方案检查点的最大误差及中误差如表1所示，试验中方案1利用GM模型和EIV模型直接定位后检查点的误差分布如图4所示。方案2-4中利用P-EIV模型区域网平差后检查点的误差分布如图5所示。

方案1为直接定位(前方交会)结果，可以看出直接定位结果均有明显系统误差；方案2-4为经区域网平差后的结果，均在不同程度上补偿了系统误差。从方案2的结果可以看出，在测区中央加入1个地面控制点后，法国地区ZY-3影像采用GM模型、EIV模型和P-EIV模型定位精度在XYZ方向分别提升为3.952、2.679和3.719m；3.457、2.169和3.189m；3.424、2.157和3.178m。从方案3的结果可以看出，在测区四角分别布设1个地面控制点时，比较P-EIV模型和GM模型的结果，法国地区ZY-3影像定位精度在XYZ方向分别提升了0.337、0.544和0.543m。从方案4的结果可以看出，在测区布设一个中央和四角地面控制点后，比较P-EIV模型和GM模型的结果，法国地区ZY-3影像定位精度在XYZ方向分别提升了0.430、0.412和0.572m。此外，比较方案2-4的结果可以看出，P-EIV模型与EIV模型的平差结果几乎完全一致，但是都比传统最小二乘方法几何定位精度高，整体提升约7.81％，并且P-EIV模型比EIV模型迭代次数较少，收敛速度较快，有利于提升整体解算效率。

以上所述，仅为发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (3)

1.一种基于Partial EIV模型的光学卫星影像区域网平差方法，其特征在于包括以下步骤：

步骤1，根据影像连接点和对应地面点的关系建立有理函数模型；

步骤2，在有理函数模型的基础上，构建区域网平差模型，解算整个测区中所有影像的仿射变换系数和所有加密点的地面坐标；

步骤3，基于PartialEIV模型补偿平差方程中的观测误差和系数矩阵误差，利用Lagrange条件极值原理，推导基于混合总体最小二乘的光束法区域网平差计算方法；

Partial EIV模型总体型同时顾及了观测向量与系数矩中包含的随机误差；包括以下公式步骤：

(1)其数学模型表示为：

L＝(A-E_A)X+e_L

(8)其随机模型为：

式中，vec(·)表示矩阵的拉直变换，

表示单位权中误差；

(9)在基于RFM模型的区域网平差方法中，系数矩阵A的元素并非全都是随机的，部分是由0或者1的常数组成。针对EIV模型中系数矩阵由一些非随机元素与随机元素组成，对系数矩阵A进一步处理，将系数矩阵中的随机元素提取出来，进而EIV模型可以写成Partial EIV模型，其数学模型表示为：

(10)其随机模型为：

式中，

表示n×m系数矩阵真值，h表示一个确定的(n×m)×1常数向量，由

中非随机元素和零组成，B表示一个给定的(n×m)×s矩阵，s表示系数矩阵

中随机元素的个数，L表示n×1的观测向量，e_L表示观测向量误差；a表示系数矩阵

中随机元素组成的s×1列向量，

表示其真值，e_a表示a的随机误差，Q_L表示观测量误差的协因数阵，Q_a表示系数矩阵误差的协因数阵，

表示单位权方差；I_n表示n×n单位矩阵，

表示Kronecker积；

(11)进一步构造平差目标函数：

(12)对式(11)求偏导，整理后得到参数估计为：

式中，

(13)改正数估计：

式中，

(14)单位权中误差估计：

(15)根据Lagrange条件极值定理，对式(9)进行重新构造目标函数如下：

式中，λ和u均为Lagrange因子；

(16)为使式(15)最小，根据条件极值定理可得到方程组：

(17)由式(16)得到：

(18)

(19)

(20)

(21)

式中，

(22) 

式中，

迭代计算步骤如下：

1)给定向量A,Q_L,Q_a,L,h,B,a，利用最小二乘计算初值

2)根据式(20)和式(21)计算e_a和e_L；

3)计算

和

4)再根据式(22)计算

5)重复步骤2)-4)直至

小于给定的阈值，此处设为1e-10，输出

利用上述Partial EIV模型对RPC的系统误差的改正参数进行估计，并对像点观测量进行更新及改正之后，回代式(6)求解地面坐标。

2.根据权利要求1所述的一种基于Partial EIV模型的光学卫星影像区域网平差方法，其特征在于：所述步骤1，根据影像连接点和对应地面点的关系建立有理函数模型中，该有理函数模型RFM是将像方影像坐标(r,c)表示为以对应物方地面三维空间坐标(P,L,H)为自变量的多项式的比值，实现与其一样的几何定位精度；包括以下公式步骤：

(2)其形式如下所示：

式中，(r_n,c_n)(P_n,L_n,H_n)分别表示像点坐标和地面点坐标正则化后的坐标，其值在(-1～+1)之间；正则化提高RPC模型中各系数求解的稳定性，同时有助于减少计算过程中由于数据级数差别过大而引起的数据舍入误差；

(3)正则化是通过平移和缩放实现的，正则化公式如下所示：

式中，r₀,c₀,P₀,L₀,H₀为正则化的平移参数；r_s,c_s,P_s,L_s,H_s为正则化的缩放系数；

(4)利用附加像方补偿方法对有理函数模型的系统误差进行改正，由式(1)可得：

式中，

分别表示光学卫星影像行、列方向系统误差改正数，(a₀,a₁,a₂,b₀,b₁,b₂)表示系统误差的补偿系数。

3.根据权利要求2所述的一种基于Partial EIV模型的光学卫星影像区域网平差方法，其特征在于：基于RFM模型的区域网平差是利用少量地面控制点解算整个测区中所有影像的仿射变换系数和所有加密点的地面坐标；光学卫星影像的RPC参数存在系统误差，通过影像之间连接点的约束关系对RFM模型的系统误差进行补偿，进而提高目标几何定位精度；包括以下公式步骤：

(5)根据式(3)进一步写成如下：

式中，

(6)对式(4)进行线性化后，

进一步得到基于RFM模型的区域网平差模型的误差方程的矩阵形式：

V＝At+Bx-L,P

式中，V＝[v_r v_c]^T表示像点行、列坐标残差向量；

分别表示相应的系数矩阵；x＝[ΔPΔLΔH]^T表示连接点对应地面点的空间大地坐标改正向量；t表示系统误差补偿模型参数改正向量；L表示观测值向量；P表示相应的权阵；

(7)注意到系数矩阵A由像点观测值(r,c)组成，因而矩阵A中包含观测误差，得到：

其中，地面点坐标未知数x的个数要远大于系统误差补偿未知数t的个数，消去x后可得t的解为：

t＝[A^TPA-A^TPB(B^TPB)^-1(B^TPA)]^-1·[A^TPL-A^TPB(B^TPB)^-1(B^TPL)]

区域网平差是一个迭代逐步趋近的过程，当两次计算的结果小于给定的阈值时，迭代计算过程结束。

Images