高分辨率卫星遥感定位的组合岭估计光束法平差方法

技术领域

本方法涉及一种高分辨率卫星遥感定位的光束法平差法，尤其是涉及一种高分辨率卫星遥感定位的组合岭估计光束法平差方法。

背景技术

在空中三角测量中，光束法平差在理论上是最严密的方法，因为这种方法直接使用了像点坐标观测值，但是实验结果表明光束法平差并不是精度最高的，主要是由于还存在残余的系统误差，解决的办法主要是自检校法。李德仁等学者将光束法平差应用于高分辨率卫星影像定位中，但是由于法方程的病态性导致计算的结果不甚理想，因为未知外方位线角元素强相关，现有的解决方法主要有四种：增加虚拟观测方程、合并强相关项、线角分求法、岭估计。在处理高分辨率遥感影像定位法方程病态问题时岭估计特别是广义岭估计近来受到越来越多的重视。但是广义岭估计首先需要获取最小二乘解来解算岭参数，而当法方程系数阵严重病态时，可能无法获得最小二乘解，因为法方程系数阵无法求逆。本方法针对该问题结合文献的方法提出一种组合岭估计方法进行单线阵高分辨率影像的光束法平差，取得了良好的效果，可达0.5m的平面定位精度和0.7m的高程定位精度。

方法内容

本方法的目的就是为了克服上述现有技术存在的缺陷而提供的一种能提高定位精度的高分辨率卫星遥感定位的组合岭估计光束法平差方法

本方法的目的可以通过以下技术方案来实现：

一种高分辨率卫星遥感定位的组合岭估计光束法平差方法，其特征在于，包括以下步骤：

1)读取像对数据、控制数据、传感数据，并设定限差与迭代次数；

2)数据预处理；

3)通过L曲线法计算岭参数；

4)将岭参数通过组合岭估计计算残差

5)将残差作为广义岭估计的岭参数矩阵，应用广义岭估计计算平差

6)判断是否中的元素值是否超过限差或迭代次数未超过设定值，若为是，返回步骤3)，若为否，执行步骤7)；

7)打印计算结果并进行精度分析。

所述的步骤2)数据预处理包括空间后方交会、空间前方交会。

所述的步骤3)中的岭参数计算过程如下：

通过L曲线法计算岭参数K值，即选择L曲线上曲率最大的点相对应的K值作为所求的参数，根据岭估计的估计准则可得和都是岭参数K的函数，选择不同的K值并以为横坐标、为纵坐标画图，得到L曲线，则L曲线上的点的曲率q的计算公式为：

$q=\frac{{\hat{\rho }}^{\prime }{\hat{\eta }}^{\prime \prime }-{\hat{\rho }}^{\prime \prime }{\hat{\eta }}^{\prime }}{{[{\left({\hat{\rho }}^{\prime }\right)}^2+{\left({\hat{\eta }}^{\prime }\right)}^2]}^{3/2}}$

对上式求最大值，就可以得到最大曲率所对应的K值即为岭参数。

所述的4)中的残差计算过程如下：将步骤3)计算得到岭参数K代入可计算得到估计的残差

所述的5)中的平差计算过程如下：将步骤4中的残差作为广义岭估计的岭参数矩阵并将其带入可求得平差

与现有技术相比，本方法具有以下优点：

1、解决了由于法方程系数严重病态而无法求得最小二乘解，从而不能进行广义岭估计的问题；

2、提高了定位精度，用本方法进行单线阵高分辨率影像的光束法平差，取得了良好的效果，可达0.5m的平面定位精度和0.7m的高程定位精度。

附图说明

图1为本发明的流程图；

图2为前方交会误差分布图(放大20倍)；

图3为基于岭估计的光束法平差误差分布图(放大500倍)；

图4为基于组合岭估计的光束法平差误差分布图(放大600倍)。

具体实施方式

下面结合附图和具体实施例对本方法进行详细说明。

如图1所示，一种高分辨率卫星遥感定位的组合岭估计光束法平差方法，包括以下步骤：

1)读取像对数据、控制数据、传感数据，并设定限差与迭代次数；

2)对数据进行空间后方交会、空间前方交会处理；

3)通过L曲线法计算岭参数；

4)将岭参数通过组合岭估计计算残差

5)将残差作为广义岭估计的岭参数矩阵，应用广义岭估计计算平差

6)判断是否中的元素值是否超过限差或迭代次数未超过设定值，若为是，返回步骤3)，若为否，执行步骤7)；

7)打印计算结果并进行精度分析。

单线阵光束法平差模型

现有大部分商用高分辨率卫星影像(QuickBird、IKONOS、SPOT等)都是单线阵沿飞行方向推扫成像。与航片不同的是，卫星影像每一行的外方位元素都不一样，不过每一扫描行的影像与被摄物体之间仍然满足共线方程。

假设飞行方向为像坐标x轴方向，扫描行方向为像坐标y轴方向，若第i扫描行的外方位元素为X_si，Y_si，Z_si，ω_si，κ_si，(X，Y，Z)为地面点的物方空间坐标；(X_si，Y_si，Z_si)为第i行投影中心的物方空间坐标。其共线方程的形式：

$0=-f\frac{a_1(X-X_{\mathrm{si}})+b_1(Y-Y_{\mathrm{si}})+c_1(Z-Z_{\mathrm{si}})}{a_3(X-X_{\mathrm{si}})+b_3(Y-Y_{\mathrm{si}})+c_3(Z-Z_{\mathrm{si}})}---\left(1\right)$

$y_i=-f\frac{a_2(X-X_{\mathrm{si}})+b_2(Y-Y_{\mathrm{si}})+c_2(Z-Z_{\mathrm{si}})}{a_3(X-X_{\mathrm{si}})+b_3(Y-Y_{\mathrm{si}})+c_3(Z-Z_{\mathrm{si}})}$

设外方位角元素分别为ω_i，κ_i，由于卫星飞行相对稳定，可以认为外方位元素为时间t的函数(一般认为取2次足够保证精度)，故并可用下述数学模型进行描述：

$X_{\mathrm{si}}=X_{S0}+{\stackrel{·}{X}}_S·t+{\stackrel{··}{X}}_S·t^2+···$

$Y_{\mathrm{si}}=Y_{S0}+{\stackrel{·}{Y}}_S·t+{\stackrel{··}{Y}}_S·t^2+···$

$Z_{\mathrm{si}}=Z_{S0}+{\stackrel{·}{Z}}_S·t+{\stackrel{··}{Z}}_S·t^2+···$

${\omega }_i={\omega }_0+\stackrel{·}{\omega }·t+\stackrel{··}{\omega }·t^2+···$

${\kappa }_i={\kappa }_0+\stackrel{·}{\kappa }·t+\stackrel{··}{\kappa }·t^2+···$

式中(X_S0，Y_S0，Z_S0，ω₀，κ₀)为起始扫描行的外方位元素；为外方位元素的一阶变率，为外方位元素的二阶变率。

将(1)式按泰勒公式展开至一次项，得误差方程：

$-a_{11}d{X}_i-a_{12}d{Y}_i-a_{13}d{Z}_i-l_x$

$-a_{21}d{X}_i-a_{22}d{Y}_i-a_{23}d{Z}_i{l}_y$

令：

$L_i=\left(\begin{array}{c}{l}_{x_i}\\ l_{y_i}\end{array}\right)$$V_i=\left(\begin{array}{c}{V}_{x_i}\\ V_{y_i}\end{array}\right)$

则(2)式可表示为：

V_i＝A_iX_i-L_i              (3)

式中a₁₁，a₁₂，…，a₂₆的具体表达式为：

$\left(\begin{array}{c}{a}_{11}=(a_{1}f+a_{3}x)/\bar{Z}\\ a_{12}=(b_{1}f+b_{3}x)/\bar{Z}\\ a_{13}=(c_{1}f+c_{3}x)/\bar{Z}\\ a_{14}=-(f+x^2/f)·b_2\\ a_{14}=-(f+x^2/f)·\sin \kappa \\ a_{16}=0\end{array}\right)$$\left(\begin{array}{c}{a}_{21}=(a_2·f)/\bar{Z}\\ a_{22}=(b_2·f)/\bar{Z}\\ a_{23}=(c_2·f)/\bar{Z}\\ a_{24}=f{b}_1+x{b}_3\\ a_{24}=-f\cos \kappa \\ a_{26}=-x\end{array}\right)$

其中：$\bar{X}=a_1(X-X_{S0})+b_1(Y-Y_{S0})+c_1(Z-Z_{S0})$

$\bar{Y}=a_2(X-X_{S0})+b_2(Y-Y_{S0})+c_2(Z-Z_{S0})$

$\bar{Z}=a_3(X-X_{S0})+b_3(Y-Y_{S0})+c_3(Z-Z_{S0})$

X＝(A^TPA)^-1A^TPL

基于组合岭估计的光束法平差

由于线角元素高度相关，法方程系数A^TPA严重病态，故采用有偏估计，有偏估计中比较实用的是广义岭估计。

广义岭估计

用表示X的广义岭估计，则：

$\hat{X}={(A^T\mathrm{PA}+\mathrm{KI})}^{-1}{A}^T\mathrm{PL}---\left(4\right)$

解算广义岭估计值的关键是确定岭参数K，其基本原则是寻找使均方差达到最小的对角阵K。而岭估计中K是个正常数。

组合岭估计

由于广义岭估计求解K需要首先算出最小二乘解，但是当法方程严重病态时往往最小二乘解无解，本文提出应用L曲线法先计算出岭估计参数K值，即选择L曲线上曲率最大的点相对应的K值作为所求的参数。根据岭估计的估计准则知道和都是岭参数K的函数，选择不同的K值并以为横坐标、为纵坐标画图，得到L曲线。令

则L曲线上的点的曲率q的计算公式为：

$q=\frac{{\hat{\rho }}^{\prime }{\hat{\eta }}^{\prime \prime }-{\hat{\rho }}^{\prime \prime }{\hat{\eta }}^{\prime }}{{[{\left({\hat{\rho }}^{\prime }\right)}^2+{\left({\hat{\eta }}^{\prime }\right)}^2]}^{3/2}}$

对上式求最大值，就可以得到最大曲率所对应的K值即为所求值。

然后应用(4)式先估计的残差作为为广义岭估计的K阵。最后应用广义岭估计进行平差：

$\tilde{X}={(A^T\mathrm{PA}+\hat{x}I)}^{-1}{A}^T\mathrm{PL}$

实施例

本方法采用QuickBird卫星的某地区异轨立体像对，为Basic产品。影像分辨率为0.7m，标称地面定位精度为23m，覆盖范围大约为18km×18km。控制点坐标采用基于某市虚拟参考网络VRS(Virtual reference station)的GPSRTK((real time kinematic)测得，精度在5cm以内，共有69个测点。所有测点的高度变化在100m以内。从中选择任选23个已知点对立体定位结果进行检验。采用VC+编制程序进行光束法平差解算得到的结果如表1所示：

表1各种计算方法的误差比较

为了更清楚的展示误差的分布，利用matlab绘制了误差图，前方交会误差分布图如图2，基于岭估计的光束法平差误差分布图如图3，而基于组合岭估计的光束法平差误差分布图如图4，图2清楚的表明前方交会解算的结果具有很强的系统误差，而基于岭估计和组合岭估计的光束法平差计算出的结果则表现出杂乱无章的误差分布，亦即没有明显的系统误差，并且本文提出的基于组合岭估计方法能获取更高精度的加密点坐标。另外控制点的分布对加密点的坐标也有很大的影响，加密点周围控制点多的比周围控制点少的精度高。