CN102901519A - 一种基于探元指向角光学推扫卫星在轨分步几何定标方法 - Google Patents

Abstract

一种基于探元指向角光学推扫卫星在轨分步几何定标方法，包括基于探元指向角的在轨几何定标模型以及分步定标解算方法；该方法采用探元指向角模型作为相机定标模型，将定标参数分为内外定标参数，采用分步定标方法先解算外定标参数，然后解算内定标参数，能够避免由于定标参数之间存在强相关性而无法获得可靠的定标结果。本发明提供了一种有效的光学推扫卫星在轨几何定标方法，能够有效地对光学推扫式卫星进行在轨几何定标、改善影像产品几何质量，实践表明该定标方法可行、有效，定标结果稳定、可靠、精度高。

Description

一种基于探元指向角光学推扫卫星在轨分步几何定标方法

技术领域

本发明属于遥感影像几何处理领域，涉及一种基于探元指向角光学推扫卫星在轨分步几何定标方法。

背景技术

在轨几何定标是充分发挥高分辨率光学遥感卫星几何性能的关键环节，对提高卫星影像产品几何质量具有重要意义。由于受卫星发射过程中应力的释放、成像环境的改变以及器件老化等各种因素的影响，卫星影像几何成像参数发生改变，实验室定标结果不再适用，必须通过在轨几何定标重新获取这些参数的精确值，否则，卫星影像产品的几何质量难以保证，无法满足后续处理和应用的要求，使其应用价值大大降低。

目前，针对近景相机、航空相机的几何定标已经开展了大量的研究与实践工作，对于光学线阵推扫式卫星的几何定标则研究得较少，尚未形成成熟的理论与方法。由于星载线阵相机高轨、窄视场角以及线阵推扫等特殊的成像特性，使得几何成像参数之间存在强相关性，若采用近景相机、航空相机的定标模型与方法，必然难以获得稳定、可靠、高精度的定标结果，因此必须针对星载线阵相机的特点构建合适的定标模型与方法。

发明内容

本发明所要解决的问题是，针对光学线阵推扫式卫星，提供一种有效的在轨几何定标方法。

本发明的技术方案为一种基于探元指向角光学推扫卫星在轨分步几何定标方法，包括以下步骤：

步骤1，利用控制点数据，在待定标影像上进行控制点量测，获取控制点量测信息。

步骤2，利用辅助数据以及实验室定标参数构建光学线阵推扫式卫星基于探元指向角的几何定标模型；

步骤3，利用分步定标方法解算内外定标参数，获取定标结果；实现方式为，首先基于最小二乘平差解算外定标参数，恢复相机坐标系在空间中的姿态；然后在此基础上，基于最小二乘平差解算内定标参数，确定CCD各探元在相机坐标系下的指向角。

而且，步骤2所述光学线阵推扫式卫星基于探元指向角的几何定标模型如下，

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

式中，(X_g,Y_g,Z_g)与(X_gps,Y_gps,Z_gps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；

分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(B_X,B_Y,B_Z)_body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψ_x(s),ψ_y(s))代表探元s在相机坐标系下的指向角，s代表探元号；

在以上几何定标模型中，待定标参数分为外定标参数X_E和内定标参数X_I，外定标参数pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；内定标参数X_I＝(ax₀,ax₁,ax₂,ax₃,ay₀,ay₁,ay₂,ay₃)，ax₀,ax₁,ax₂,ax₃,ay₀,ay₁,ay₂,ay₃为内定标探元指向角模型的系数。

而且，步骤3包括以下子步骤，

步骤3.1，设在待定标影像上量测了K个均匀分布的高精度地面控制点作为定向点，记为RP_i，i＝1,2,3,...,K；各定向点对应的物方点和像方点分别记为RPG_i和RPM_i，物方点RPG_i的WGS84地心直角坐标为(X_i,Y_i,Z_i)，像方点RPM_i的影像坐标为(s_i,l_i)；

步骤3.2，令式（1）中：

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

式（1）转化为式（2）：

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

上式中，矢量 $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ 为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；a₁,b₁,c₁ a₂,b₂,c₂ a₃,b₃,c₃分别代表相机安装矩阵的9个元素；F(X_E,X_I)、G(X_E,X_I)分别为沿轨指向角残差与垂轨指向角残差；

步骤3.3，对外定标参数XE和内定标参数XI赋初值

步骤3.4，将内定标参数XI的当前值视为真值，将外定标参数XE视为待求的未知参数，将内定标参数XI和外定标参数XE的当前值

代入式（2），对每个定向点RPi，对式（2）进行线性化处理，建立误差方程（3），

V_i＝A_iX-L_i，权为P_i        （3）

其中 $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入式（2）计算得到的误差向量；A_i是误差方程式的系数矩阵；X代表外定标参数改正数dX_E=(dpitch,droll,dyaw)，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程；

按式（4）计算法方程系数矩阵，

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

                  （4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

上式中，矩阵 $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ 矩阵

矩阵 $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算X，如式（5），

X＝(A^TPA)^-1(A^TPL)        （5）

利用式（6）更新外定标参数XE的当前值，然后返回执行步骤4迭代计算，迭代停止后进入步骤3.4；

$X_E=X_E^o+X---\left(6\right)$

步骤3.5，解算内定标参数，将步骤3.4所得外定标参数XE的当前值视为真值，而内定标参数XI则视为待求的未知参数，将内定标参数XI和外定标参数XE的当前值

代入公式（2），对每个定向点RPi，对式（2）进行线性化处理，建立误差方程（7），

V_i＝B_iY-L_i  权为P_i        （7）

其中 $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入公式（2）计算得到的误差向量；B_i是误差方程式的系数矩阵；Y代表内定标参数改正数dX_I，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程；

按式（8）计算法方程系数矩阵；

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

                   （8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

上式中， $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算Y，如式（9）；

Y＝(B^TPB)^-1(B^TPL)        （9）

利用式（10）更新内定标参数XI的当前值，然后返回执行步骤5迭代计算，迭代停止后进入步骤3.6；

$X_I=X_I^o+Y---\left(10\right)$

步骤3.6，根据步骤3.4和步骤3.5所得外定标参数的当前值X_E和内定标参数X_I的当前值，更新相机参数文件。

本发明的优点在于：1.基于探元指向角的定标模型确定的不再是各探元在相机视场中的位置，而是其在相机坐标系下的指向角，一方面可以简化内定标模型，避免了相机内部各种畸变误差之间由于强相关性造成的过度参数化问题；同时，由于外定标参数恢复相机坐标系在空间中的姿态（三轴指向），内定标参数确定各探元在相机坐标系下的指向角，外定标参数为内定标参数提供参考基准，内定标参数在外定标参数提供的参考基准下进行确定，内外定标参数共同确定各探元在空间中的绝对指向，这样能够避免内外定标参数之间的相关性问题，提高内外定标参数解算稳定性；2.采用分步定标方法，先解算外定标参数恢复相机坐标系在空间中的姿态，为内定标参数的确定提供参考基准，然后在此基础上解算内定标参数确定各探元在相机坐标系（参考基准）下的指向角，能够避免定标参数同时解算由于强相关性造成的方程病态问题，有利于获得稳定、可靠、高精度的定标结果。

附图说明

图1为本发明实施例的流程示意图；

图2为探元指向角示意图。

具体实施方式

以下结合附图和实施例详细说明本发明具体实施方式。参见图1，实施例的流程可以分为三个步骤，每个步骤实施的具体方法、公式以及流程如下：

1.利用控制点数据，在待定标影像上量测控制点信息。为了保证定标结果的解算精度，对于所量测的控制点在数量以及分布上提供了建议：在待定标影像中，所量测的控制点在沿轨方向上尽量分布于较短的区域内，垂轨方向则应均匀覆盖整个CCD范围。控制点数量方面，在合理成本的前提下尽可能多。

2.利用姿轨、成像时间等辅助数据以及实验室定标参数，构建光学线阵推扫式卫星基于探元指向角的几何定标模型如式（1）：

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

式中，(X_g,Y_g,Z_g)与(X_gps,Y_gps,Z_gps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标，其中，后者由卫星上搭载的GPS获取。

分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(B_X,B_Y,B_Z)_body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；其中，

由星敏、陀螺通过组合定姿得到，(B_X,B_Y,B_Z)_body则在卫星发射前由实验室检校，(t)表示当前参数是一个随时间变化的量；(ψ_x(s),ψ_y(s))代表探元s在相机坐标系下的指向角，s代表探元号，如图2所示：相机坐标系X₁Y₁Z₁的原点记为O₁，探元V_Image的指向角为(ψ_x,ψ_y)；在该定标模型中，待定标参数分为外定标参数X_E和内定标参数X_I，其中，外定标参数

（pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角）用于补偿相机安装角误差，恢复相机坐标系在空间中的指向，为内定标参数的解算确定参考基准；内定标参数X_I＝(ax₀,ax₁,ax₂,ax₃,ay₀,ay₁,ay₂,ay₃)则用于补偿由于相机内部各种畸变导致的像点误差，确定CCD各探元在相机坐标系（参考基准）下的指向角。内定标参数与外定标参数两者共同恢复CCD各探元在空间中的绝对指向。与常规相机检校模型不同之处在于，本发明提出的基于探元指向角的定标模型确定的不再是各探元在相机视场中的位置，而是其在相机坐标系下的指向角，这样可以避免定标参数之间由于存在强相关性而无法获得稳定、可靠的定标结果。

3.利用分步定标方法解算内外定标参数，分步解算内外定标参数，外定标参数为内定标参数提供参考基准，内定标参数在外定标参数提供的参考基准下确定。其原理为：首先基于最小二乘平差解算外定标参数，恢复相机坐标系在空间中的姿态；然后在此基础上，基于最小二乘平差解算内定标参数，确定CCD各探元在相机坐标系下的指向角。

实施例的解算公式及步骤如下：

步骤3.1，假设在待定标影像上量测了K个均匀分布的高精度地面控制点作为定向点，记为RP_i，这里，i＝1,2,3,...,K；其对应的物方点和像方点分别记为RPG_i和RPM_i，物方点RPG_i的WGS84地心直角坐标为(X_i,Y_i,Z_i)，像方点RPM_i的影像坐标为(s_i,l_i)；

步骤3.2，令式（1）中：

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

式（1）可转化为式（2）：

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

上式中，矢量 $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ 为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；a₁,b₁,c₁ a₂,b₂,c₂ a₃,b₃,c₃分别代表相机安装矩阵的9个元素；F(X_E,X_I)、G(X_E,X_I)分别为沿轨指向角残差与垂轨指向角残差；

步骤3.3，对外定标参数X_E、内定标参数X_I赋初值

这里均为实验室定标值。

步骤3.4，将内定标参数X_I的当前值视为“真值”，将外定标参数X_E视为待求的未知参数。将它们的当前值

代入式（2），对每个定向点RP_i，对式（2）进行线性化处理，建立误差方程（3）（下标i代表利用定向点RP_i建立的误差方程）：

Vi＝A_iX-L_i，权为P_i            （3）

其中 $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入式（2）计算得到的误差向量；A_i是误差方程式的系数矩阵；X代表外定标参数改正数dX_E，即(dpitch,droll,dyaw)，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程。

按式（4）计算法方程系数矩阵，

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

                     （4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

上式中，矩阵 $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ 矩阵

矩阵 $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算X，如式（5）；

X＝(A^TPA)^-1(A^TPL)        （5）

利用式（6）更新外定标参数X_E的当前值，然后返回执行步骤4迭代计算。当外定标参数改正数X中3个参数(dpitch,droll,dyaw)均小于阈值10^-12时，迭代停止。

$X_E=X_E^o+X---\left(6\right)$

步骤3.5，解算内定标参数，与步骤4解算内定标参数实现方式一致，此时将外定标参数的当前值视为“真值”，而内定标参数则视为待求的未知参数，将内外定标参数的当前值

代入公式（2），此时

为步骤3.4所得外定标参数X_E的当前值。内定标参数X_I和外定标参数X_E的对每个定向点RP_i，对式（2）进行线性化处理，建立误差方程（7）（下标i代表利用定向点RP_i建立的误差方程）：具体解算公式和过程如下：

V_i＝B_iY-L_i  权为P_i        （7）

其中 $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入公式（2）计算得到的误差向量；B_i是误差方程式的系数矩阵；Y代表内定标参数改正数dX_I，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程。

按式（8）计算法方程系数矩阵；

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

                   （8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

上式中， $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算Y，如式（9）；

Y＝(B^TPB)^-1(B^TPL)        （9）

利用式（10）更新内定标参数X_I的当前值，然后返回执行步骤5迭代计算。当内定标参数改正数Y中8个参数(dax₀,dax₁,dax₂,dax₃,day₀,day₁,day₂,day₃)均小于阈值10^-12时，迭代停止。

$X_I=X_I^o+Y---\left(10\right)$

步骤3.6，根据得到的定标结果，即内外定标参数当前值，更新相机参数文件。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种基于探元指向角光学推扫卫星在轨分步几何定标方法，包括以下步骤：

步骤1，利用控制点数据，在待定标影像上进行控制点量测，获取控制点量测信息。

步骤2，利用辅助数据以及实验室定标参数构建光学线阵推扫式卫星基于探元指向角的几何定标模型；

步骤3，利用分步定标方法解算内外定标参数，获取定标结果；实现方式为，首先基于最小二乘平差解算外定标参数，恢复相机坐标系在空间中的姿态；然后在此基础上，基于最小二乘平差解算内定标参数，确定CCD各探元在相机坐标系下的指向角。

2.如权利要求1所述基于探元指向角光学推扫星在轨分步几何定标方法，其特征在于：步骤2所述光学线阵推扫式卫星基于探元指向角的几何定标模型如下，

$\left(\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ \tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\\ 1\end{array}\right)={\mathrm{\λR}}_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw}){(R_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}})---\left(1\right)$

$\left\{\begin{array}{c}{\mathrm{\ψ}}_x\left(s\right)={\mathrm{ax}}_0+{\mathrm{ax}}_1\×s+{\mathrm{ax}}_2\×s^2+{\mathrm{ax}}_3\×s^3\\ {\mathrm{\ψ}}_y\left(s\right)={\mathrm{ay}}_0+{\mathrm{ay}}_1\×s+{\mathrm{ay}}_2\×s^2+{\mathrm{ay}}_3\×s^3\end{array}\right.$

式中，(X_g,Y_g,Z_g)与(X_gps,Y_gps,Z_gps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标；

分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵；(B_X,B_Y,B_Z)_boay代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标；(t)表示当前参数是一个随时间变化的量；(ψ_x(s),ψ_y(s))代表探元s在相机坐标系下的指向角，s代表探元号；

在以上几何定标模型中，待定标参数分为外定标参数X_E和内定标参数X_I，外定标参数

pitch、roll、yaw分别为俯仰、翻滚以及偏航方向夹角；内定标参数X_I＝(ax₀,ax₁,ax₂,ax₃,ay₀,ay₁,ay₂,ay₃)，ax₀,ax₁,ax₂,ax₃,ay₀,ay₁,ay₂,ay₃为内定标探元指向角模型的系数。

3.如权利要求2所述基于探元指向角光学推扫星在轨分步几何定标方法，其特征在于：步骤3包括以下子步骤，

步骤3.1，设在待定标影像上量测了K个均匀分布的高精度地面控制点作为定向点，记为RP_i，i＝1,2,3,...,K；各定向点对应的物方点和像方点分别记为RPG_i和RPM_i，物方点RPG_i的WGS84地心直角坐标为(X_i,Y_i,Z_i)，像方点RPM_i的影像坐标为(s_i,l_i)；

步骤3.2，令式（1）中：

$\left\{\begin{array}{c}\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]\left(\begin{array}{c}{R}_{J2000}^{\mathrm{body}}{R}_{\mathrm{wgs}}^{J2000}{\left[\begin{array}{c}{X}_g-X_{\mathrm{gps}}\left(t\right)\\ Y_g-Y_{\mathrm{gps}}\left(t\right)\\ Z_g-Z_{\mathrm{gps}}\left(t\right)\end{array}\right]}_{\mathrm{wgs}}-{\left[\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right]}_{\mathrm{body}}\end{array}\right)\\ R_{\mathrm{body}}^{\mathrm{cam}}(\mathrm{pitch},\mathrm{roll},\mathrm{yaw})=\left[\begin{array}{c}{a}_1,b_1,c_1\\ a_2,b_2,c_2\\ a_3,b_3,c_3\end{array}\right]\end{array}\right.$

式（1）转化为式（2）：

$\left\{\begin{array}{c}F(X_E,X_I)=\frac{a_1\mathrm{Ux}+b_1\mathrm{Uy}+c_1\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_x\left(s\right)\right)\\ G(X_E,X_I)=\frac{a_2\mathrm{Ux}+b_2\mathrm{Uy}+c_2\mathrm{Uz}}{a_3\mathrm{Ux}+b_3\mathrm{Uy}+c_3\mathrm{Uz}}-\tan \left({\mathrm{\ψ}}_y\left(s\right)\right)\end{array}\right.---\left(2\right)$

上式中，矢量 $\left[\begin{array}{c}{U}_x\\ U_y\\ U_z\end{array}\right]$ 为物方矢量U，代表从相机投影中心到物方点的矢量在本体坐标系下的坐标；a₁,b₁,c₁ a₂,b₂,c₂ a₃,b₃,c₃分别代表相机安装矩阵的9个元素；F(X_E,X_I)、G(X_E,X_I)分别为沿轨指向角残差与垂轨指向角残差；

步骤3.3，对外定标参数XE和内定标参数XI赋初值

步骤3.4，将内定标参数X_I的当前值视为真值，将外定标参数X_E视为待求的未知参数，将内定标参数X_I和外定标参数X_E的当前值

代入式（2），对每个定向点RP_i，对式（2）进行线性化处理，建立误差方程（3），

V_i＝A_iX-L_i，权为P_i        （3）

其中 $A_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_E}\\ \frac{\∂G_i}{\∂X_E}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂\mathrm{pitch}}\frac{\∂F_i}{\∂\mathrm{roll}}\frac{\∂F_i}{\∂\mathrm{yaw}}\\ \frac{\∂G_i}{\∂\mathrm{pitch}}\frac{\∂G_i}{\∂\mathrm{roll}}\frac{\∂G_i}{\∂\mathrm{yaw}}\end{array}\right]$ $X=d{X}_E=\left[\begin{array}{c}\mathrm{dpitch}\\ \mathrm{droll}\\ \mathrm{dyaw}\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入式（2）计算得到的误差向量；Ai是误差方程式的系数矩阵；X代表外定标参数改正数dX_E=(dpitch,droll,dyaw)，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程；

按式（4）计算法方程系数矩阵，

$A^T\mathrm{PA}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{A}_i$

                        （4）

$A^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{A}_i^T{P}_i{L}_i$

上式中，矩阵 $A=\left[\begin{array}{c}{A}_1\\ ...\\ A_i\\ ...\\ A_K\end{array}\right],$ 矩阵矩阵 $L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算X，如式（5），

X＝(A^TPA)^-1(A^TPL)        （5）

利用式（6）更新外定标参数X_E的当前值，然后返回执行步骤4迭代计算，迭代停止后进入步骤3.4；

$X_E=X_E^o+X---\left(6\right)$

步骤3.5，解算内定标参数，将步骤3.4所得外定标参数X_E的当前值视为真值，而内定标参数X_I则视为待求的未知参数，将内定标参数X_I和外定标参数X_E的当前值代入公式（2），对每个定向点RP_i，对式（2）进行线性化处理，建立误差方程（7），

V_i＝B_iY-L_i  权为P_i            （7）

其中 $B_i=\left[\begin{array}{c}\frac{\∂F_i}{\∂X_I}\\ \frac{\∂G_i}{\∂X_I}\end{array}\right]=\left[\begin{array}{c}\frac{\∂F_i}{\∂{\mathrm{ax}}_0}\frac{\∂F_i}{\∂{\mathrm{ax}}_1}\frac{\∂F_i}{\∂{\mathrm{ax}}_2}\frac{\∂F_i}{\∂{\mathrm{ax}}_3}\frac{\∂F_i}{\∂{\mathrm{ay}}_0}\frac{\∂F_i}{\∂{\mathrm{ay}}_1}\frac{\∂F_i}{\∂a{y}_2}\frac{\∂F_i}{\∂{\mathrm{ay}}_3}\\ \frac{\∂G_i}{\∂{\mathrm{ax}}_0}\frac{\∂G_i}{\∂{\mathrm{ax}}_1}\frac{\∂G_i}{\∂{\mathrm{ax}}_2}\frac{\∂G_i}{\∂{\mathrm{ax}}_3}\frac{\∂G_i}{\∂{\mathrm{ay}}_0}\frac{\∂G_i}{\∂{\mathrm{ay}}_1}\frac{\∂G_i}{\∂{\mathrm{ay}}_2}\frac{\∂G_i}{\∂{\mathrm{ay}}_3}\end{array}\right]$ $Y={\mathrm{dX}}_I=\left[\begin{array}{c}{\mathrm{dax}}_0\\ {\mathrm{dax}}_1\\ {\mathrm{dax}}_2\\ {\mathrm{dax}}_3\\ {\mathrm{day}}_0\\ {\mathrm{day}}_1\\ {\mathrm{day}}_2\\ {\mathrm{day}}_3\end{array}\right]$ $L_i=\left[\begin{array}{c}F(X_E^o,X_I^o)\\ G(X_E^o,X_I^o)\end{array}\right]$

式中，L_i是利用内外定标参数当前值

代入公式（2）计算得到的误差向量；Bi是误差方程式的系数矩阵；Y代表内定标参数改正数dX_I，d代表改正数符号；P_i是当前定向点RP_i的像点量测精度对应的权；F_i和G_i分别为沿轨指向角残差F(X_E,X_I)、垂轨指向角残差G(X_E,X_I)的函数模型，微分后得到相应误差方程；

按式（8）计算法方程系数矩阵；

$B^T\mathrm{PB}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{B}_i$

                   （8）

$B^T\mathrm{PL}=\underset{i=1}{\overset{K}{\∑}}{B}_i^T{P}_i{L}_i$

上式中， $B=\left[\begin{array}{c}{B}_1\\ ...\\ B_i\\ ...\\ B_K\end{array}\right],$

$L=\left[\begin{array}{c}{L}_1\\ ...\\ L_i\\ ...\\ L_K\end{array}\right]$

利用最小二乘平差计算Y，如式（9）；

Y＝(B^TPB)^-1(B^TPL)        （9）

利用式（10）更新内定标参数XI的当前值，然后返回执行步骤5迭代计算，迭代停止后进入步骤3.6；

$X_I=X_I^o+Y---\left(10\right)$

步骤3.6，根据步骤3.4和步骤3.5所得外定标参数的当前值XE和内定标参数XI的当前值，更新相机参数文件。

Images