CN103514362B - 基于模型误差补偿的两行根数生成方法 - Google Patents

Abstract

本发明公开了一种基于模型误差补偿的两行根数生成方法。该方法首先将观测到的位置和速度矢量由地心赤道惯性坐标系转化为真赤道平春分点坐标系，然后计算相应的瞬时轨道根数并将其作为两行根数迭代计算的初始估计值；利用初始估计值分别预报相应采样时刻的位置和速度向量，计算其与实际观测采样时刻的位置和速度向量的差值并判断是否满足终止条件；通过数值微分法计算简化普适摄动轨道预报函数在采样时刻关于两行根数的偏导数矩阵；通过求解补偿最小二乘问题得到两行根数的修正量；计算新的两行根数，并进行重复迭代直至满足精度要求。本发明不仅能够统一处理单点拟合和采样拟合的情况，而且拟合得到的两行根数用于轨道预报时平均误差更小。

Description

基于模型误差补偿的两行根数生成方法

技术领域

本发明属于空间目标轨道计算技术，特别是一种基于模型误差补偿的两行根数生成方法。

背景技术

两行轨道根数(TLE)是美国战略司令部发布空间目标编目数据采用的格式，是一种基于开普勒根数的平均轨道元，已被广泛用于描述空间目标（包括小卫星以及空间碎片等）的飞行轨道。对两行根数而言，其是用特定方法去掉了周期扰动项的平均轨道根数，它必须与特定的轨道预报模型一起使用来预报空间目标某一时刻的位置和速度，对近地目标而言该预报模型就是北美防空司令部(NORAD)提供的SGP4模型。作为目前应用最为广泛的一类解析模型，它考虑了地球非球形J2、J3、J4项摄动、大气摄动及日月引力摄动等摄动因素的影响，在轨道预报精度和计算效率之间作了较好的平衡。

目前，大部分空间目标的两行根数会在Goddard空间飞行中心定期公布更新，只要将TLE代入已经公开的SGP4模型,即可求出任意时刻空间目标的位置和速度，但是TLE生成算法至今仍未公开，这使得基于TLE的进一步应用受到较大限制。早在1967年德国的Walter（1.Walter H G.Conversion of osculating obbital elements into mean elements[J].The Astronomical Journal,1967,72(8):994-997）就在Brouwer和Kozai的轨道理论意义下研究了由瞬时轨道根数转换为相应的平均轨道根数问题，提出了一个迭代算法并证明了其收敛性。在SGP4模型出现后，作为一种特殊的且被广泛应用的平均轨道根数，两行根数的生成算法引起了广泛的研究。其中，Jochim等（2.Jochim E F,Gill E,Montenbruck O,et al.GPS based on-board and onground orbit operations for small satellites[J].Acta Astronautica,1996,39(9):917-922）以空间目标的星载GPS数据作为观测输人实时生成该目标的TLE，但该方法仅适用于合作目标。Lee（3.Byoung S L.NORAD TLEconversion from osculating orbital elements[J].Journal of Astronomy and SpaceSciences,2002,19(4):395-402）提出利用空间目标的瞬时密切轨道数据拟合出该目标的TLE，但是只研究了单点拟合的方法，且未考虑大气阻力因子的影响；Dwight（4.Dwight EA.Computing NORAD mean orbital elements from a state vector[D].Ohio:Air ForceInstitute of Technology,Wright-Patterson Air Force Base,1994）考虑了由位置和速度矢量生成两行根数的问题，提出了牛顿法和直接迭代两种方法，但是该方法并未解决在偏心率和轨道倾角很小时出现奇点导致无法收敛的问题；刘光明等（5.刘光明，文援兰，廖瑛.基于无奇异变换的双行轨道根数生成算法[J].***工程与电子技术，2011，33(5):1104-1107）提出了基于无奇异变换的空间目标TLE生成算法，针对TLE采样拟合过程中可能出现奇点的问题引入无奇异轨道根数代替开普勒根数并给出了目标状态矢量关于两行根数偏导数矩阵的近似解析表达式，但其拟合的目标函数仅考虑了位置向量。胡敏等（6.胡敏，曾国强.平均轨道根数与密切轨道根数的互换[J].飞行器测控学报,2012,31(2):77-81）基于TLE和SGP4模型从正反两方面总结了平均轨道根数与密切轨道根数的互换问题,仅用单点拟合法进行了计算。

由于SGP4模型本身是有误差的，这影响了两行根数拟合的精度，针对当前研究存在的问题，本发明在考虑SGP4模型误差的基础上提出了两行根数拟合的误差补偿优化模型，将单点拟合和采样拟合的情况统一处理，使得拟合得到的两行根数用于轨道预报时平均误差更小。

发明内容

发明目的：本发明的目的在于提供一种基于简化普适摄动轨道预报模型（SGP4）的空间目标两行根数生成方法。

技术方案：实现本发明目的的技术解决方案为：一种基于模型误差补偿的两行根数生成方法，步骤如下：

步骤1、将观测到的位置和速度矢量由地心赤道惯性（J2000）坐标系转化为真赤道平春分点（TEME）坐标系；

步骤2、利用采样初始时刻的位置和速度矢量计算相应的瞬时轨道根数并将其作为两行根数迭代计算的初始估计值；

步骤3、利用简化普适摄动轨道预报模型SGP4（取B^*=0）和两行根数的初始估计值分别预报相应采样时刻的位置和速度向量，计算其与实际观测采样时刻的位置和速度向量的差值并判断是否满足终止条件。若满足，算法终止；若不满足，转步骤4；

步骤4、通过数值微分法计算简化普适摄动轨道预报函数在相应采样时刻关于两行根数的偏导数矩阵；

步骤5、选取合适的平滑因子和正则矩阵，通过求解补偿最小二乘问题得到两行根数的修正量；

步骤6、计算新的两行根数，并进行重复迭代，直至满足精度要求。

有益效果：本发明与现有技术相比，其显著优点是：1）在考虑了简化普适摄动轨道预报模型的模型误差基础上，能够统一处理单点拟合和采样拟合的情况；2）拟合的目标函数同时兼顾了位置和速度向量，使得拟合得到的两行根数用于轨道预报时平均误差更小。

附图说明

图1是基于模型误差补偿的两行根数生成方法流程图。

具体实施方式

下面结合说明书附图对本发明进行进一步详述：

如图1所示，本发明涉及一种基于模型误差补偿的两行根数生成方法，具体步骤如下：

步骤1、将观测到的位置和速度矢量Y₀,Y₁,…,Y_n由地心赤道惯性（J2000）坐标系转化为真赤道平春分点（TEME）坐标系，即作如下变换：

${\left[\begin{array}{c}x\\ y\\ z\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·}{z}\end{array}\right]}_{\mathrm{TEME}}=R_z^{-1}\left(\mathrm{\Δ\φ}\cos \stackrel{\‾}{\mathrm{\ϵ}}\right)\·N\·P\·{\left[\begin{array}{c}x\\ y\\ z\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·}{z}\end{array}\right]}_{J2000},$

式中，表示t时刻观测到的位置和速度矢量且Y_i=Y(t_i),i=0,1,…,n，R_z为坐标旋转矩阵，N为章动矩阵，P为岁差矩阵，Δφ为黄金章动，为平黄赤交角。

步骤2、利用采样初始时刻t₀时的位置r=(x,y,z)^T和速度计算相应的瞬时轨道根数，将其作为两行根数迭代计算的初始估计值X₀=(e₀,i₀,Ω₀,ω₀,M₀,n₀)^T；具体为：

步骤2-1、计算中间变量如下：

W_x=h_x/||h||，W_y=h_y/||h||,W_x=h_x/||h||，

其中， $G{M}_{\⊕}=398600.4415k{m}^3/s^2$ 且

$h=\left(\begin{array}{c}{h}_x\\ h_y\\ h_z\end{array}\right)=\left(\begin{array}{c}y\stackrel{\·}{z}-z\stackrel{\·}{y}\\ z\stackrel{\·}{x}-x\stackrel{\·}{z}\\ x\stackrel{\·}{y}-y\stackrel{\·}{x}\end{array}\right),$ $\left|\right|h\left|\right|=\sqrt{h_x^2+h_y^2+h_z^2},$ $\left|\right|r\left|\right|=\sqrt{x^2+y^2+z^2},$ $\left|\right|\stackrel{\·}{r}\left|\right|=\sqrt{{\stackrel{\·}{x}}^2+{\stackrel{\·}{y}}^2+{\stackrel{\·}{z}}^2};$

步骤2-2、计算轨道倾角和升交点赤经如下：

$i_0=\arctan \left(\frac{\sqrt{{W_x}^2+{W_y}^2}}{W_z}\right),$ ${\mathrm{\Ω}}_0=\arctan (-\frac{W_x}{W_z});$

步骤2-3、计算平均角速度和偏心率如下：

$n_0=\sqrt{\frac{G{M}_{\⊕}}{a^3}},$ $e_0=\sqrt{1-\frac{p}{a}};$

步骤2-4、计算平近点角如下：

M₀=E₀-esin(E₀),

其中，

$E_0=\arctan \left(\frac{(x\stackrel{\·}{x}+y\stackrel{\·}{y}+z\stackrel{\·}{z})/\left(a^2{n}_0\right)}{1-\left|\right|r\left|\right|/a}\right);$

步骤2-5、计算近地点辐角如下：

${\mathrm{\ω}}_0=\arctan \left(\frac{z}{-x{W}_y+y{W}_x}\right)-\arctan \left(\frac{\sqrt{1-e_0^2}\sin \left(E_0\right)}{\cos \left(E_0\right)-e_0}\right).$

步骤3、利用简化普适摄动轨道预报模型SGP4（取B^*=0）和两行根数的当前估计值分别预报相应采样时刻的位置和速度向量，并计算其与实际观测采样时刻的位置和速度向量的差值；具体为：

步骤3-1、利用SGP4模型的预报函数F(X₀,t₀,t)计算t₀,t₁,…,t_n时刻的位置和速度向量即

$\tilde{y}\left(t\right)=F(X_0,t_0,t),$ t=t₀,t₁,…,t_n.

其中，X₀为t₀时刻两行根数的初始估计值。

步骤3-2、计算当前根数预报的位置和速度向量与实际观测的位置和速度向量在采样时刻t₀,t₁,…,t_n的差值f₀,f₁,…,f_n，即

$f\left(t\right)=Y\left(t\right)-\tilde{Y}\left(t\right),$ t=t₀,t₁,…,t_n.

步骤3-3、判断差值大小是否满足终止条件。若满足，算法终止；若不满足，转步骤4；

步骤4、通过数值微分法计算简化普适摄动轨道预报函数在采样时刻t₀,t₁,…,t_n关于两行根数的偏导数矩阵A_i=A(t_i),i=0,1,…,n，即

$A\left(t_i\right)=\frac{\∂F(X,t_0,t)}{\∂X}{|}_{t=t_i}={\left[\begin{array}{cccccc}\frac{\∂y_1}{\∂x_1}& \frac{\∂y_1}{\∂x_2}& \frac{\∂y_1}{\∂x_3}& \frac{\∂y_1}{\∂x_4}& \frac{\∂y_1}{\∂x_5}& \frac{\∂y_1}{\∂x_6}\\ \frac{\∂y_2}{\∂x_1}& \frac{\∂y_2}{\∂x_2}& \frac{\∂y_2}{\∂x_3}& \frac{\∂y_2}{\∂x_4}& \frac{\∂y_2}{\∂x_5}& \frac{\∂y_2}{\∂x_6}\\ \frac{\∂y_3}{\∂x_1}& \frac{\∂y_3}{\∂x_2}& \frac{\∂y_3}{\∂x_3}& \frac{\∂y_3}{\∂x_4}& \frac{\∂y_3}{\∂x_5}& \frac{\∂y_3}{\∂x_6}\\ \frac{\∂y_4}{\∂x_1}& \frac{\∂y_4}{\∂x_2}& \frac{\∂y_4}{\∂x_3}& \frac{\∂y_4}{\∂x_4}& \frac{\∂y_4}{\∂x_5}& \frac{\∂y_4}{\∂x_6}\\ \frac{\∂y_5}{\∂x_1}& \frac{\∂y_5}{\∂x_2}& \frac{\∂y_5}{\∂x_3}& \frac{\∂y_5}{\∂x_4}& \frac{\∂y_5}{\∂x_5}& \frac{\∂y_5}{\∂x_6}\\ \frac{\∂y_6}{\∂x_1}& \frac{\∂y_6}{\∂x_2}& \frac{\∂y_6}{\∂x_3}& \frac{\∂y_6}{\∂x_4}& \frac{\∂y_6}{\∂x_5}& \frac{\∂y_6}{\∂x_6}\end{array}\right]}_{t=t_i},i=\mathrm{0,1},\·\·\·,n.$

其中，

$F(X,t_0,t)=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\\ y_5\\ y_6\end{array}\right],$ $X=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\end{array}\right],$

$\frac{\∂y_i}{\∂x_j}=\frac{y_i(x_1,\·\·\·,x_j+\mathrm{\Δ}{x}_j/2,\·\·\·,x_6,t)-y_i(x_1,\·\·\·,x_j-\mathrm{\Δ}{x}_j/2,\·\·\·,x_6,t)}{\mathrm{\Δ}{x}_j},i,j=\mathrm{1,2},\·\·\·,6.$

Δx_j可根据精度要求选取，一般可取Δx_j=10^-3x_j。

步骤5、选取合适的平滑因子和正则矩阵，通过求解补偿最小二乘问题得到两行根数的修正量；具体为：

步骤5-1、选取平滑因子α和正则矩阵R=P^TP，并计算矩阵S=(I+αR)^-1，其中

α可根据经验尝试选取，取值范围一般为0<α<1，I为相应阶数的单位矩阵。

步骤5-2、计算两行根数的修正向量：

ΔX₀=(A^T(I-S)A)^-1A^T(I-S)f.

其中，

$f=\left[\begin{array}{c}{f}_0\\ f_1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_n\end{array}\right],$ $A=\left[\begin{array}{c}{A}_0\\ A_1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ A_n\end{array}\right],$

步骤6、计算新的两行根数，并进行重复迭代，直至满足精度要求。

$X_0=X_0^{*}+\mathrm{\&Delta;}{X}_0,$ 可令 $X_0^{*}=X_0$ 进行迭代求解。

由上可知，本发明提出了关于两行根数更一般的优化计算模型，既适用于单点拟合也适用于采样拟合，相比现有的方法适用范围更广。既考虑了SGP4的模型误差，也考虑了速度矢量的偏差作为优化目标，相比现有方法得到的两行根数，其用于轨道预报的平均精度更高。

Claims (1)

1.一种基于模型误差补偿的两行根数生成方法，其特征在于，用于空间目标轨道计算技术，包括如下步骤：

1)将观测到的位置和速度向量均由地心赤道惯性坐标系转化为真赤道平春分点坐标系；

2)将采样初始时刻的位置和速度向量转换为瞬时轨道根数，将其作为两行根数迭代计算的初始估计值；

3)利用简化普适摄动轨道预报模型和两行根数的初始估计值分别预报相应采样时刻的位置和速度向量，计算其与实际观测采样时刻的位置和速度向量的差值并判断是否满足终止条件；若满足，算法终止；若不满足，转步骤4)；

4)通过数值微分法计算简化普适摄动轨道预报函数在采样时刻关于两行根数的偏导数矩阵；

5)选取合适的平滑因子和正则矩阵，通过求解补偿最小二乘问题得到两行根数的修正量；

步骤5)包括以下步骤：

S5.1：选取平滑因子α和正则矩阵R＝P^TP，并计算矩阵S＝(I+αR)^-1，其中

α可根据经验尝试选取，取值范围为0＜α＜1，I为相应阶数的单位矩阵；

S5.2：计算两行根数的修正向量：

ΔX₀＝(A^T(I-S)A)^-1A^T(I-S)f

其中，

$f=\left[\begin{array}{c}{f}_0\\ f_1\\ .\\ .\\ .\\ f_n\end{array}\right],A=\left[\begin{array}{c}{A}_0\\ A_1\\ .\\ .\\ .\\ A_n\end{array}\right]$

其中，A_i为两行根数的偏导数矩阵，且i＝0,1，…，n，A为由A₁，…，A_n组成的矩阵，f_i是当前根数预报的位置和速度向量与实际观测的位置和速度向量在采样时刻t₀,t₁,…,t_n的差值，且i＝0,1，…，n，f为由f₁，…，f_n组成的矩阵；

6)计算新的两行根数，并进行重复迭代，直至满足精度要求；迭代采用如下公式：

令进行反复迭代求解

其中，X₀为两行根数迭代计算的初始估计值。