CN101226062B - 一种星上实时计算环月轨道的方法 - Google Patents

Abstract

一种星上实时计算环月轨道的方法，涉及航天器实时计算环月轨道技术领域，包括以下步骤：(1)以定轨结果作为初值，用数值方法计算出实际轨道；(2)选取一条二体轨道作为参考轨道；(3)计算实际轨道相对于参考轨道的运动轨迹；(4)根据公知的两条二体轨道之间的相对运动规律，计算出一条相对于参考轨道运动轨迹比实际轨道相对于参考轨道运动轨迹更接近的二体轨道；步骤(1)至步骤(4)均在地面完成，然后将计算结果注入到星上；(5)卫星采用二体轨道计算方法，根据步骤(4)中算得的二体轨道实时计算出卫星的轨道根数。本发明使得星上在计算量很小的情况下能够实施计算出卫星轨道，同时不需要地面和卫星保持实时联系。

Description

一种星上实时计算环月轨道的方法

技术领域

本发明涉及航天器实时计算环月轨道技术领域。

背景技术

目前星上实时轨道计算一般只在近地轨道卫星上使用，所使用的方法是平根数轨道计算方法。国内航天器上采用注入19根数进行计算，国外采用SGP4模型进行计算，这两种方法类似，都是采用平均轨道根数及其长期项系数、长周期项系数、短周期项系数计算卫星瞬时轨道根数。由于地球引力场的J2项摄动和大气阻力摄动比其他项大的多，在使用平根数计算时只需要考虑地球引力的J2项和大气阻力项，因而公式较为简单，计算量不大。而月球引力场的J2项摄动的量级与其他高阶项接近，如果采用平根数法进行计算，则需要考虑20阶月球引力和地球引力摄动，即使经过简化后计算量仍然非常大，不适合星上计算机进行计算。

在月球卫星环月期间，定姿需要卫星在月球坐标系内实时的轨道根数，对卫星轨道计算精度提出的要求是位置误差小于1km。如果采用传统的平根数法，要做到较高精度的轨道计算，由于月球引力场的复杂性，计算时采用的月球引力模型就需要取到20阶左右的较高阶次，另外在较高精度要求下还需要考虑地球引力摄动的影响，从而导致计算量过大，不适合在星上计算。而数值解的方法尽管精度较高，但也存在计算量较大的问题。

发明内容

本发明的技术解决问题是：克服现有技术的不足之处，提供一种地面运算与星上运算相结合的方式，使得环月卫星能够在星上实施计算出环月轨道。

本发明的方法的技术解决方案是：一种星上实时计算环月轨道的方法，包括以下步骤：

(1)以定轨结果作为初值，用RKF78数值方法计算出实际轨道；

(2)选取一条二体轨道作为参考轨道；

(3)计算实际轨道相对于参考轨道的运动轨迹，采用三角函数来拟合实际轨道相对于参考轨道的运动；

(4)根据公知的两条二体轨道之间的相对运动规律，计算出一条相对于参考轨道运动轨迹比实际轨道相对于参考轨道运动轨迹更接近的二体轨道；

步骤(1)至步骤(4)均在地面完成，将一天分成12个时间段，分别重复利用步骤(1)至步骤(4)计算相应的二体轨道，然后将计算结果注入到星上；

(5)卫星采用二体轨道计算方法，根据步骤(4)中算得的二体轨道实时计算出卫星的轨道根数。

本发明与现有技术相比的有益效果是：本发明通过地上运算与星上运算相结合的方式，特别是选取了一条二体轨道作为参考轨道，使得星上在计算量很小的情况下能够实施计算出卫星轨道，同时不需要地面和卫星保持实时联系，这种轨道计算方法的位置误差小于1km。

附图说明

图1为本发明星上计算轨道的水平面内误差情况；

图2为本发明星上计算轨道的轨道面内误差情况。

具体实施方式

(1)在地面通过高精度的数值运算计算出月球卫星在未来一段时间内的月心惯性系轨道根数，以一定的步长h保存(h≤60s)，轨道计算的初值由定轨所得。因为摄动原因得到的结果随时间而变化，以σ_t表示此根数，并认为此根数表示了卫星的实际位置。计算过程需要考虑的摄动包括月球引力非球形引力摄动，其阶次应取到30阶以上，同时还需考虑地球引力、太阳引力等对月球卫星的摄动。计算方法采用高精度的数值解法。

例如，取卫星定轨的结果为：

起始时间(北京时)t0＝2007年11月7日12:0:0.000

卫星轨道根数为：

半长轴a＝1934km；

偏心率e＝0.00237；

倾角i＝87.41度；

升交点赤经Ω＝265.51度；

近地点幅角ω＝281.585；

平近点角M＝202。

精确轨道计算采用RKF78数值方法，考虑月球引力场为70×70阶、考虑地球引力摄动、太阳引力摄动。

(2)选取一条二体轨道作为参考轨道，使得上述计算出的轨道相对于该二体轨道在小范围内运动。可以选择开始时刻对应的二体轨道作为参考轨道，其轨道根数记为σ₀。参考轨道选择为：

半长轴a(m)	偏心率e	倾角i(度)	升交点赤经Ω(度)	近地点幅角ω(度)	平近点角M(度)
						1934	0.00237	87.41	265.51	281.585	202

(3)计算实际轨道相对于参考轨道的运动轨迹。利用坐标变换可以得到σ_t相对于σ₀的运动轨迹：

(x_So，y_So，z_So)^T

＝L_t[R_z(u₀)R_x(i₀)R_z(Ω₀)R_z(-Ω_t)R_x(-i_t)R_z(-u_t)(r_t，0，0)^T-(r₀，0，0)^T]

其中：

$L_t=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& -1\\ -1& 0& 0\end{array}\right]$ R_x，R_y，R_z表示旋转矩阵

下标t表示(1)中计算出的数值解轨道σ_t，下标0表示参考轨道，r₀、r_t分别表示σ₀、σ_t的月心距。

(4)对得到的相对运动结果x_So，y_So做最小二乘拟合，对x方向和y方向的相对运动采用的拟合函数分别为：

x＝c₀+c₁M+c₂sinM+c₃cosM

y＝b₁sin M+b₂cos M

其中，在一个拟合时间段内(a_r确定的周期)M取值从0开始，以的步长增加。

解法方程A^TAc＝A^Tx得到拟合公式的系数c＝(c₀，c₁，c₂，c₃)^T，其中

$x={(x_{{1S}_0},x_{{2S}_0},\·\·\·,x_{{\mathrm{nS}}_0})}^T,$ $A=\left(\begin{array}{cccc}1& M_0& \sin 0& \cos 0\\ 1& M_1& \sin M_1& \cos M_1\\ 1& M_2& \sin M_2& \cos M_2\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& M_{n-1}& \sin M_{n-1}& \cos M_{n-1}\end{array}\right)$

解法方程A^TAb＝A^Ty得到拟合公式的系数b＝(b₁，b₂)^T，其中，

$y={(y_{{1S}_0},y_{{2S}_0},\·\·\·,y_{{\mathrm{nS}}_0})}^T,$ $A=\left(\begin{array}{cc}\sin 0& \cos 0\\ \sin M_1& \cos M_1\\ \sin M_2& \cos M_2\\ .& .\\ .& .\\ .& .\\ \sin M_{n-1}& \cos M_{n-1}\end{array}\right)$

上面A矩阵中的M_k＝kn₀h，(k＝1，2...n-1)，

由拟合系数计算相对轨道根数：

$D=\frac{c_1{n}_0}{a_0},$ $\mathrm{\Δ}{M}^{\′}=\frac{c_0}{a_0},$ u₀₀＝ω₀+M₀₀

其中M₀₀为参考轨道σ₀在一个拟合时间段初始时刻的平近点角M

$\mathrm{\Δ}{e}_x=\frac{(c_2\cos u_{00}+c_3\sin u_{00})}{2a_0}$

$\mathrm{\Δ}{e}_y=\frac{(c_2\sin u_{00}-c_3\cos u_{00})}{2a_0}$

$\mathrm{\Δ}{i}_x=\frac{b_2\cos u_{00}-b_1\sin u_{00}}{a_0}$

$\mathrm{\Δ}{i}_y=\frac{b_1\cos u_{00}+b_2\sin u_{00}}{a_0}$

由相对轨道根数求拟合轨道根数：

$a_f=a_0-\frac{2a_{0}D}{3\sqrt{G_m/{a_0}^3}}={\left[\frac{G_m}{{(D+n_0)}^2}\right]}^{1/3}$

$e_f=\sqrt{{(e_0\cos {\mathrm{\ω}}_0+\mathrm{\Δ}{e}_x)}^2+{(e_0\sin {\mathrm{\ω}}_0+\mathrm{\Δ}{e}_y)}^2}$

$i_f=a\cos (\mathrm{\Δ}{i}_y\sin i_0+\sqrt{1-\mathrm{\Δ}{i_x}^2-\mathrm{\Δ}{i_y}^2}\cos i_0)$

$\left\{\begin{array}{c}\sin {\mathrm{\Ω}}_f=\frac{1}{\sin i_f}(\mathrm{\Δ}{i}_x\cos {\mathrm{\Ω}}_0-\mathrm{\Δ}{i}_y\sin {\mathrm{\Ω}}_0\cos i_0+\sqrt{1-\mathrm{\Δ}{i_x}^2-\mathrm{\Δ}{i_y}^2}\sin {\mathrm{\Ω}}_0\sin i_0)\\ \cos {\mathrm{\Ω}}_f=-\frac{1}{\sin i_f}(\mathrm{\Δ}{i}_x\sin {\mathrm{\Ω}}_0+\mathrm{\Δ}{i}_y\cos {\mathrm{\Ω}}_0\cos i_0-\sqrt{1-\mathrm{\Δ}{i_x}^2-\mathrm{\Δ}{i_y}^2}\cos {\mathrm{\Ω}}_0\sin i_0)\end{array}\right.$

$\left\{\begin{array}{c}\sin {\mathrm{\ω}}_f=e_0\sin {\mathrm{\ω}}_0+\mathrm{\Δ}{e}_y\\ \cos {\mathrm{\ω}}_f=e_0\cos {\mathrm{\ω}}_0+\mathrm{\Δ}{e}_x\end{array}\right.$

M_f0＝ω₀+M₀₀+ΔM′-ω_f

每组轨道的使用时间段取为1或2个平均轨道周期，对12个拟合时间段分别使用以上(1)至(4)中所述的方法进行拟合，得到12组轨道根数，在加上轨道参数初始时刻对应的星时t₀、使用时间段对应的轨道半长轴a₀两个参数，最终形成上行注入的74个参数，分别用以下符号表示：

t<sub>0</sub>(s)	初始参数注入时刻
		a<sub>0</sub>(m)	使用时间段对应的轨道半长轴
a<sub>j0</sub>(m)，(j＝1...12)	注入的一组二体轨道根数之半长轴
		e<sub>j0</sub>，(j＝1...12)	注入的一组二体轨道根数之偏心率
i<sub>j0</sub>(rad)，(j＝1...12)	注入的一组二体轨道根数之倾角
		Ω<sub>j0</sub>(rad)，(j＝1...12)	注入的一组二体轨道根数之升交点赤径
ω<sub>j0</sub>(rad)，(j＝1...12)	注入的一组二体轨道根数之近地点幅角
		M<sub>j0</sub>(rad)，(j＝1...12)	注入的一组二体轨道根数之初始平近点角

具体计算结果为：

序号	半长轴a(m)	偏心率e	倾角i(度)	升交点赤经Ω(度)	近地点幅角ω(度)	平近点角M(度)
							1	1934.20921	0.00229469041	87.4186294	265.512224	280.67202	202.915221
2	1934.18368	0.00230057926	87.4332619	265.509514	284.119813	199.833363
							3	1934.15943	0.00229811322	87.4483868	265.507274	288.184877	196.137375
4	1934.13817	0.0022926994	87.4638563	265.50566	292.864857	191.829531
							5	1934.12125	0.00229118924	87.4795892	265.504796	298.066838	187.00259
6	1934.10953	0.00230004874	87.4955511	265.504779	303.595317	181.851867
							7	1934.10345	0.00232363461	87.5117213	265.505683	309.173532	176.653944
8	1934.10295	0.00236270604	87.5280748	265.507569	314.493695	171.716408
							9	1934.10767	0.0024144502	87.5445687	265.510485	319.279916	167.314922
10	1934.11699	0.00247347148	87.5611424	265.514464	323.335734	163.64566
							11	1934.1301	0.00253354174	87.5777315	265.519521	326.557106	160.812299
12	1934.14614	0.00258919365	87.5942805	265.525643	328.91743	158.840992

(5)卫星采用二体轨道计算方法，根据(4)中算得的二体轨道实时计算出卫星的轨道根数。计算公式为：

$n_0=\sqrt{\frac{G_m}{{a_0}^3}};$ dt＝t-t₀；

$j=\mathrm{int}\left(\frac{\mathrm{dt}}{2\mathrm{\π}}{n}_0\right)+1;$ 确定使用第j组轨道根数

a＝a_j0

${\mathrm{\ω}}_0=\sqrt{\frac{G_m}{a^3}}$

e＝e_j0

i＝i_j0

Ω＝Ω_j0

ω＝ω_j0

$M=M_{j0}+{\mathrm{\ω}}_0\·[\mathrm{dt}-\frac{2\mathrm{\π}(j-1)}{n_0}];$ M＝mod(M，π)；

$E=E-\frac{E-e\sin E-M}{1-e\cos E}$ (以E＝M为初值，迭代5次)

$f={\tan 2}^{-1}(\sqrt{1-e^2}\sin E,\cos E-e)$

u＝ω+f

$\mathrm{\ρ}={\sin }^{-1}\left(\frac{R_m(1+e\cos f)}{a(1-e^2)}\right)$

r＝a(1-e cos E)

惯性系至轨道坐标系转换矩阵记为

其中

a₁₁＝-sin(u)cos(Ω)-cos(u)cos(i)sin(Ω)

a₁₂＝-sin(u)sin(Ω)+cos(u)cos(i)cos(Ω)

a₁₃＝cos(u)sin(i)

a₂₁＝-sin(i)sin(Ω)

a₂₂＝sin(i)cos(Ω)

a₂₃＝-cos(i)

a₃₁＝-cos(u)cos(Ω)+sin(u)cos(i)sin(Ω)

a₃₂＝-cos(u)sin(Ω)-sin(u)cos(i)cos(Ω)

a₃₃＝-sin(u)sin(i)

A_OI用作卫星相对轨道坐标系进行姿态确定的输入。

如图1、2所示，为星上使用地面注入的二体轨道进行计算与地面高精度计算结果相比的水平面内误差情况和轨道面内误差情况。图中×轴表示卫星飞行方向，y轴表示轨道负法线方向，z轴表示月心方向，三个轴的单位都为km，从图中可见星上计算卫星×方向误差小于200m，y方向误差小于400m，z方向误差小于300m。

Claims (4)

1.一种星上实时计算环月轨道的方法，其特征在于包括以下步骤：

(1)以定轨结果作为初值，用数值方法计算出实际轨道；

(2)选取一条二体轨道作为参考轨道；

(3)计算实际轨道相对于参考轨道的运动轨迹；

(4)根据公知的两条二体轨道之间的相对运动规律，计算出一条相对于参考轨道运动轨迹比实际轨道相对于参考轨道运动轨迹更接近的二体轨道；

步骤(1)至步骤(4)均在地面完成，然后将计算结果注入到星上；

(5)卫星采用二体轨道计算方法，根据步骤(4)中算得的二体轨道实时计算出卫星的轨道根数。

2.根据权利要求1所述的星上实时计算环月轨道的方法，其特征在于：所述步骤(1)中的数值方法为RKF78数值方法。

3.根据权利要求1所述的星上实时计算环月轨道的方法，其特征在于：所述步骤(3)中计算实际轨道相对于参考轨道的运动轨迹，是采用三角函数来拟合实际轨道相对于参考轨道的运动轨迹。

4.根据权利要求1所述的星上实时计算环月轨道的方法，其特征在于：将一天分成12个时间段，分别重复利用步骤(1)至步骤(4)计算相应的二体轨道，再注入星上计算相应的轨道根数。

Images