CN102878995A - 一种静止轨道卫星自主导航方法 - Google Patents

Abstract

本发明涉及一种静止轨道卫星自主导航方法，属于卫星自主导航研究领域。基于星敏感器和地球敏感器的kalman滤波算法，实时得到卫星相对于定点的位置偏差值；对位置偏差值用最小二乘方法进行数据处理，获得一天内的轨道平面内的平均轨道根数，再用平均滤波方法获得一天内的轨道平面外的平均轨道根数；以平均轨道根数作为星上轨道解析外推算法的输入值，外推计算一天内卫星轨道位置，提供连续导航定位信息，实现卫星自主导航功能。本发明的方法已经在中星2A上成功应用，经过在轨标定后自主轨道确定精度优于10km，该方法可以推广应用于所有要求具备自主功能的地球静止轨道卫星。

Description

一种静止轨道卫星自主导航方法

技术领域

本发明涉及一种静止轨道卫星自主导航方法，属于卫星自主导航研究领域，可以应用于静止轨道卫星自主导航任务。

背景技术

传统卫星在寿命周期内完全依赖地面测控***的支持，与之相比，中星2A卫星是我国第一颗具有自主生存能力的地球静止轨道通信卫星，其所采用的自主导航方法为该技术首次工程应用，所用的方法属于天文导航方法，这一大类方法从算法角度已经较为成熟，并也产生了多种新型敏感器，如空间六分仪等，如中星2A上用到的基于地球敏感器和星敏感器的导航滤波算法是比较适合目前发展趋势的方法，首先是方法较为简单(已经获得我国专利保护)，其次是方法所用到的部件属于星上常用姿态确定部件，能够长期在轨不间断工作，即导航方案不对星上资源产生新的需求；在原有导航滤波方法的基础上，中星2A针对静止轨道卫星的运行要求，尤其是轨道控制要求，创造性将星上轨道解析外推算法和导航滤波后的处理结果有机结合起来，既使得卫星自主获取了实时位置，同时也获得了适合于轨道控制的轨道信息。

发明内容

本发明的目的是为了提出一种静止轨道卫星自主导航方法，该方法能够给出静止轨道卫星实时的轨道信息。

本发明的目的是通过以下技术方案实现的。

本发明的一种静止轨道卫星自主导航方法，该方法的步骤为：

1)基于星敏感器和地球敏感器的kalman滤波算法，实时得到卫星相对于定点的位置偏差值；

2)对步骤1)得到的位置偏差值用最小二乘方法进行数据处理，获得一天内的轨道平面内的平均轨道根数：即半长轴、偏心率、近地点幅角和真近点角；对步骤1)得到的位置偏差值用平均滤波方法获得一天内的轨道平面外的平均轨道根数：即轨道倾角和升交点赤经；

3)以步骤2)得到的平均轨道根数作为星上轨道解析外推算法的输入值，外推计算一天内卫星轨道位置，提供连续导航定位信息，实现卫星自主导航功能。

有益效果

本发明的方法已经在中星2A上成功应用，经过在轨标定后自主轨道确定精度优于10km，该方法可以推广应用于所有要求具备自主功能的地球静止轨道卫星。

附图说明

图1为卫星自主导航测量原理示意图；

图2为轨道平面参数确定示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

实施例

1)卫星自主导航所需的星上测量部件为星敏感器和地球敏感器，基于星敏感器和地球敏感器的kalman滤波算法

以此标称轨道定点位置建立Hill方程，即卫星相对标称轨道定点位置的相对运动方程，作为kalman滤波算法的状态方程，方程形式如下：

$\stackrel{\·\·}{\mathrm{\α}}-2{\mathrm{\ω}}_0\stackrel{\·}{\mathrm{\γ}}=\frac{a_x}{r}$

$\stackrel{\·\·}{\mathrm{\β}}+{\mathrm{\ω}}_0^2\mathrm{\β}=\frac{a_y}{r}$

$\stackrel{\·\·}{\mathrm{\γ}}+2{\mathrm{\ω}}_0\stackrel{\·}{\mathrm{\α}}-3{\mathrm{\ω}}_0^2\mathrm{\γ}=\frac{a_z}{r}$

其中：

α为经度偏差，单位为rad；

β为纬度偏差，单位为rad；

γ为径向偏差，单位为rad；

r为静止轨道卫星理论半径，单位为km；

ω₀为标称静止轨道角速度，单位为rad/s；

a_x、a_y、a_z是外界加速度在标称轨道定点位置为原点的坐标系的分量；

标称轨道定点位置定义为：不考虑卫星漂移且轨道倾角为零的情况下的理论定点位置，与地球自转完全同步；

卫星自主导航测量原理如图1所示，包括地球、卫星和标称轨道定点位置；其中，

表示星敏感器给出的第一个方向矢量，

表示星敏感器给出的第二个方向矢量，

代表地球敏感器给出的地心方向矢量。

方向矢量

在卫星本体系下的坐标位置由星敏的安装矩阵决定；根据星敏感器的测量原理，其测量输出为方向矢量

在惯性坐标系中的方位坐标，分别记为：

和

具体表达为 ${\left({\stackrel{\→}{S}}_1\right)}_I=\left[\begin{array}{ccc}{S}_{1\mathrm{ix}}& S_{1\mathrm{iy}}& S_{1\mathrm{iz}}\end{array}\right],$ ${\left({\stackrel{\→}{S}}_2\right)}_I=\left[\begin{array}{ccc}{S}_{2\mathrm{ix}}& S_{2\mathrm{iy}}& S_{2\mathrm{iz}}\end{array}\right].$

红外地球敏感器的测量输出为滚动角φ_h和俯仰角θ_h。根据测量得到的滚动角和俯仰角，经过两次旋转变换即可得到在本体坐标系下卫星轨道坐标系原点到地心的位置矢量，记为：

$\stackrel{\→}{E}=\left\{\begin{array}{c}{E}_{b1}=-\sin \left({\mathrm{\θ}}_h\right)\\ E_{b2}=\sin \left({\mathrm{\φ}}_h\right)\cos \left({\mathrm{\θ}}_h\right)\\ E_{b3}=\cos \left({\mathrm{\φ}}_h\right)\cos \left({\mathrm{\θ}}_h\right)\end{array}\right.$

从标称轨道定点位置的指向地心的单位矢量记为

在卫星本体坐标系得到两个夹角η₁、η₂：

$\cos {\mathrm{\η}}_1=\stackrel{\→}{E}\·\stackrel{\→}{S_1}$

$\cos {\mathrm{\η}}_2=\stackrel{\→}{E}\·\stackrel{\→}{S_2}$

在惯性坐标系下计算得到另外两个夹角η₁′、η₂′：

$\cos {\mathrm{\η}}_1^{\′}=\stackrel{\→}{Z}\·\stackrel{\→}{{\left(S_1\right)}_I}$

$\cos {\mathrm{\η}}_2^{\′}=\stackrel{\→}{Z}\·\stackrel{\→}{{\left(S_2\right)}_I}$

可得到测量方程g：

g＝(cosη₁-cosη′₁ cosη₂-cosη′₂)^T

为此构建一个自主导航的模型，其对象方程为：

$\left\{\begin{array}{c}\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{bU}\\ g=h\left(x\right)\end{array}\right.$

其中：

$x={\left[\begin{array}{cccccc}\mathrm{\α}& \mathrm{\β}& \mathrm{\γ}& \stackrel{\·}{\mathrm{\α}}& \stackrel{\·}{\mathrm{\β}}& \stackrel{\·}{\mathrm{\γ}}\end{array}\right]}^T$

$A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& {2\mathrm{\ω}}_0\\ 0& -{\mathrm{\ω}}_0^2& 0& 0& 0& 0\\ 0& 0& 3{\mathrm{\ω}}_0^2& -2{\mathrm{\ω}}_0& 0& 0\end{array}\right]$

b＝[0 0 0 1 1 1]^T

U为外部输入。

将测量方程在参考轨道处线性化：

$H\left(x\right)=\frac{\∂h\left(x\right)}{\∂x}=\left[\begin{array}{cccccc}{h}_{11}& h_{12}& 0& 0& 0& 0\\ h_{21}& h_{22}& 0& 0& 0& 0\end{array}\right]$

其中：

h₁₁＝-sin(α_T+α)cos(δ_T-β)S_1ix+cos(α_T+α)cos(δ_T-β)S_1iy

h₁₂＝cos(α_T+α)sin(δ_T-β)S_1ix+sin(α_T+α)sin(δ_T-β)S_1iy-cos(δ_T-β)S_1iz

h₂₁＝-sin(α_T+α)cos(δ_T-β)S_2ix+cos(α_T+α)cos(δ_T-β)S_2iy

h₂₂＝cos(α_T+α)sin(δ_T-β)S_2ix+sin(α_T+α)sin(δ_T-β)S_2iy-cos(δ_T-β)S_2iz

α_T为标称轨道定点位置的赤经，δ_T为标称轨道定点位置的赤纬

根据上述状态方程和测量方程，采用扩展Kalman滤波方法，可获得当前卫星轨道相对于标称轨道定点位置的偏移量，即：

$x={\left[\begin{array}{cccccc}\mathrm{\α}& \mathrm{\β}& \mathrm{\γ}& \stackrel{\·}{\mathrm{\α}}& \stackrel{\·}{\mathrm{\β}}& \stackrel{\·}{\mathrm{\γ}}\end{array}\right]}^T$

2)对步骤1)得到的位置偏差值用最小二乘方法进行数据处理，获得一天内的轨道平面内的平均轨道根数：即半长轴、偏心率、近地点幅角和真近点角；对步骤1)得到的位置偏差值用平均滤波方法获得一天内的轨道平面外的平均轨道根数：即轨道倾角和升交点赤经；

根据步骤1)中kalman滤波后获得的真实卫星相对于标称轨道定点位置的相对位置偏差值，为了获得便于轨道控制的控制变量以及提高导航精度，取滤波结果中的东西偏差α和南北偏差β，基于静止轨道卫星的受摄漂移原理，通过最小二乘法和平均滤波方法获得卫星的平均轨道六根数，包括：半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角和真近点角，其中平均的意义为：利用一天的α和β数据拟合获得。

根据α和β，可以转换成卫星的赤经和赤纬，即赤经为α_T+α：，赤纬为δ_T-β，可以获得当前k时刻的地心单位矢量

其表达式为：

$\stackrel{\→}{i_k}=\left[\begin{array}{ccc}\cos ({\mathrm{\α}}_{\mathrm{Tk}}-{\mathrm{\α}}_k)\cos ({\mathrm{\δ}}_{\mathrm{Tk}}-{\mathrm{\β}}_k)& \sin ({\mathrm{\α}}_{\mathrm{Tk}}-{\mathrm{\α}}_k)\cos ({\mathrm{\δ}}_{\mathrm{Tk}}-{\mathrm{\β}}_k)& \sin ({\mathrm{\δ}}_{\mathrm{Tk}}-{\mathrm{\β}}_k)\end{array}\right].$

由若干组这样的单位矢量即可确定卫星的轨道平面参数，即轨道倾角i和升交点赤经Ω，采用的方法为平均滤波方法。

在一个轨道周期T内，首先每隔设定时间dt，可获得一个地心单位矢量，分别记为：(k＝1...N)，其中N＝T/dt，如图2所示，则通过叉乘可得一个倾角矢量：

$\stackrel{\&RightArrow;}{{i_1}^{\&prime;}}=\stackrel{\&RightArrow;}{i_1}\&times;\stackrel{\&RightArrow;}{i_2}...\stackrel{\&RightArrow;}{{i_k}^{\&prime;}}=\stackrel{\&RightArrow;}{i_k}\&times;\stackrel{\&RightArrow;}{i_{k+1}}(k<N)$

最终形成一个总的倾角矢量：

$\stackrel{\&RightArrow;}{i}=\frac{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\stackrel{\&RightArrow;}{{i_k}^{\&prime;}}}{N}$

如此在一个轨道周期内通过一系列倾角矢量得到一个总的倾角矢量，根据倾角矢量

在惯性坐标系中的XOY平面内投影定义：[sin i cos(Ω-90°)，sin i sin(Ω-90°)]，可解算出一组轨道倾角i和升交点赤经Ω。由于该倾角矢量是由一个轨道周期内的数据平均后获得，所以可以认为求出的i和Ω为上一个轨道周期内轨道倾角和升交点赤经的平均值，即一天数据可获得一组

和

采用递推最小二乘法对平面内轨道参数进行估计，即卫星东西方向运行规律可简化为：

α＝a₀+a₁t_orb+b₁cos(ω_et_orb)+b₂sin(ω_et_orb)

其中，t_orb一个轨道周期的时间变量，ω_e为静止轨道角速度，待估参数为a₀、a₁、b₁、b₂。

通过递推最小二乘法可实时估计参数a₀、a₁、b₁、b₂，由此在一个轨道周期结束时，可解算出卫星轨道平面内的平均轨道根数：

半长轴可由a₁直接解出：

其中r_s为静止轨道标称半径，a₁对应的单位为rad/天。

根据数学模型2esin(ω_et_orb+θ)，其中θ＝atan2(b₁，b₂)为相位偏移量，由b₁、b₂可解算出偏心率

即：

$\stackrel{\&OverBar;}{e}=0.5\sqrt{b_1^2+b_2^2}$

并且，令 ${\stackrel{\&RightArrow;}{i}}_{\mathrm{\&Omega;}}=\left[\begin{array}{ccc}\cos \left(\stackrel{\&OverBar;}{Q}\right)& \sin \left(\stackrel{\&OverBar;}{Q}\right)& 0\end{array}\right]$ 与当前的卫星地心单位矢量

进行内积，再求反余弦获得卫星幅角u，另外，真近点角

则：

$\stackrel{\&OverBar;}{\mathrm{\&omega;}}=\mathrm{\&alpha;}-\stackrel{\&OverBar;}{\mathrm{\&Omega;}}-\stackrel{\&OverBar;}{f}$

3)以步骤2)得到的六个平均轨道根数作为星上轨道解析外推算法的输入值，外推计算一天内卫星轨道位置，提供连续导航定位信息，实现卫星自主导航功能。

由于地球同步轨道卫星偏心率e和轨道倾角i量级都比较小，故在进行解析外推计算时需采用适用于0≤e＜1和0≤i＜180°的第二类六个无奇点根数：

σ＝[a，h，k，ξ，η，λ]^T

其表达是为：a为半长轴，h＝sin i cosΩ，k＝sin i sinΩ，

$\mathrm{\&lambda;}=M+\widetilde{\mathrm{\&omega;}},$ 其中 $\widetilde{\mathrm{\&omega;}}=\mathrm{\&omega;}+\mathrm{\&Omega;.}$

将步骤2)得到的六个平均轨道根数代入上式，令f≈M，可获得第二类无奇点根数的平均值： ${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_0={\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{a}}_0,& {\stackrel{\&OverBar;}{h}}_0,& {\stackrel{\&OverBar;}{k}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_0\end{array}\right]}^T,$ 其中的

即为作为星上轨道解析外推算法的初始拟平均根数，且根据步骤2)的结果每天更新一次，进而获得自主导航期间的实时轨道根数。

采用拟平均根数法构造的分析解可以获得瞬时轨道根数σ(t)：

$\mathrm{\&sigma;}\left(t\right)=\stackrel{\&OverBar;}{\mathrm{\&sigma;}}+{\mathrm{\&sigma;}}_s^{\left(1\right)}\left(t\right)$

$\stackrel{\&OverBar;}{\mathrm{\&sigma;}}\left(t\right)={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_0+(\mathrm{\&delta;}{\stackrel{\&OverBar;}{n}}_0+{\mathrm{\&sigma;}}_1)\left({t-t}_0\right)$

$\mathrm{\&delta;}{\stackrel{\&OverBar;}{n}}_0={\left(\begin{array}{cccccc}0& 0& 0& 0& 0& {\stackrel{\&OverBar;}{n}}_0\end{array}\right)}^T$

${\stackrel{\&OverBar;}{n}}_0={{\stackrel{\&OverBar;}{a}}_0}^{-3/2}$

其中：

σ(t)为瞬时轨道根数

为拟平均轨道根数

σ₁为一阶长期项系数，包含长期项和长周期项

两部分

为一阶短周期项

为初始拟平均根数

t₀为初始拟平均根数所对应的轨道历元

为轨道平均角速度

地球非球形长期项系数σ₁：

a₁＝0

$h_1=-{\stackrel{\&OverBar;}{k}}_0{\mathrm{\&Omega;}}_1$

$k_1={\stackrel{\&OverBar;}{h}}_0{\mathrm{\&Omega;}}_1$

${\mathrm{\&xi;}}_1=-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0{\mathrm{\&omega;}}_1$

${\mathrm{\&eta;}}_1={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0{\mathrm{\&omega;}}_1$

λ₁＝M₁+ω₁+Ω₁

其中Ω₁，ω₁和M₁由下式表达：

$M_1=\frac{{3J}_2}{{2p}^2}n(1-\frac{3}{2}{\sin }^{2}i)\sqrt{1-e^2}$

${\mathrm{\&omega;}}_1=\frac{3J_2}{2p^2}n(2-\frac{5}{2}{\sin }^{2}i)$

${\mathrm{\&Omega;}}_1=-\frac{{3J}_2}{{2p}^2}n\cos i,$

上述各式中的a，e，i均为

且有：

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{a}={\stackrel{\&OverBar;}{a}}_0,\stackrel{\&OverBar;}{e}={({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0^2+{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0^2)}^{1/2},\\ {\sin }^2\stackrel{\&OverBar;}{i}=({\stackrel{\&OverBar;}{h}}_0^2+{\stackrel{\&OverBar;}{k}}_0^2),\cos i={[1-({\stackrel{\&OverBar;}{h}}_0^2+{\stackrel{\&OverBar;}{k}}_0^2)]}^{1/2}\\ \stackrel{\&OverBar;}{n}={\stackrel{\&OverBar;}{a}}_0^{-3/2},\stackrel{\&OverBar;}{p}=\stackrel{\&OverBar;}{a}(1-e^2)={\stackrel{\&OverBar;}{a}}_0(1-e_0^2)\end{array}\right.$

地球非球形短周期项

$a_s^{\left(1\right)}\left(t\right)=\frac{{3J}_2}{2a}\left\{(\frac{2}{3}-{\sin }^{2}i)\right[{\left(\frac{a}{r}\right)}^3-{(1-e^2)}^{-3/2}]+{\left(\frac{a}{r}\right)}^3[(h^2-k^2)\cos 2u+2\mathrm{hk}\sin 2u\left]\right\}$

${h_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{{4a}^2}\right)(h\cos 2u+k\sin 2u)$

${k_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{{4a}^2}\right)(h\sin 2u-k\cos 2u)$

${{\mathrm{\&lambda;}}_s}^{\left(1\right)}=\left(\frac{{21J}_2}{4a^2}\right)(\mathrm{\&xi;}\sin u-\mathrm{\&eta;}\cos u)$

${{\mathrm{\&xi;}}_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{4a^2}\right)[2(\cos u+\mathrm{\&xi;})+(\mathrm{\&xi;}\cos 2u+\mathrm{\&eta;}\sin 2u)]$

${{\mathrm{\&eta;}}_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{4a^2}\right)[2(\sin u+\mathrm{\&eta;})+(\mathrm{\&xi;}\sin 2u-\mathrm{\&eta;}\cos 2u)]$

上述短周期项各式右端出现的根数均为

而u和a/r由下列相应公式进行迭代计算获得：

$\left\{\begin{array}{c}\tilde{u}-\mathrm{\&lambda;}=\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u}\\ \frac{a}{r}={(1-\mathrm{\&xi;}\cos \tilde{u}-\mathrm{\&eta;}\sin \tilde{u})}^{-1}\\ \sin u=\left(\frac{a}{r}\right)[(\sin \tilde{u}-\mathrm{\&eta;})-\frac{\mathrm{\&xi;}}{1+\sqrt{1-e^2}}(\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u})]\\ \cos u=\left(\frac{a}{r}\right)[(\cos \tilde{u}-\mathrm{\&xi;})+\frac{\mathrm{\&eta;}}{1+\sqrt{{1-e}^2}}(\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u})]\end{array}\right.$

其中

u＝f+ω，E是偏近点角。

地球非球形长周期项系数：

$a_l^{\left(1\right)}\left(t\right)=\left(\frac{{6J}_{\mathrm{2,2}}}{a}\right)\frac{\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})}{(1-\mathrm{\&alpha;})}$

$h_l^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_{\mathrm{2,2}}}{{2a}^2}\right)\frac{1}{1-\mathrm{\&alpha;}}[k\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})-h\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})]$

$k_l^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_{\mathrm{2,2}}}{{2a}^2}\right)\frac{1}{1-\mathrm{\&alpha;}}[h\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})+k\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})]$

${\mathrm{\&Delta;\&xi;}}_l^{\left(1\right)}\left(t\right)=O\left(e\right)$

${\mathrm{\&Delta;\&eta;}}_l^{\left(1\right)}\left(t\right)=O\left(e\right)$

${\mathrm{\&lambda;}}_l^{\left(1\right)}\left(t\right)=\left(\frac{{9J}_{\mathrm{2,2}}}{a^2}\right)\frac{\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})}{(1-\mathrm{\&alpha;})}$

上述各式中的a，h，k等量均为

各式中的谐系数J_2，2以及

等量的意义如下：

其中，n_e为标称轨道角速度，C_2，2，S_2，2是非归一化的谐系数，d为J2000.0起算的儒略日，S_G是格林尼治恒星时，λ_2，2＝-14.545°。

日月引力摄动部分：

长期项系数σ₁：

a₁＝0

$h_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i)\cos i\cos \mathrm{\&Omega;}+({\mathrm{AA}}_i+{\mathrm{BB}}_i)\sin \mathrm{\&Omega;}]n$

$k_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i)\cos i\sin \mathrm{\&Omega;}-({\mathrm{AA}}_i+{\mathrm{BB}}_i)\cos \mathrm{\&Omega;}]n$

${\mathrm{\&xi;}}_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[5\mathrm{\&xi;}\left(\mathrm{AB}\right)-\mathrm{\&eta;}(1+B^2-{4A}^2)]n$

${\mathrm{\&eta;}}_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[5\mathrm{\&eta;}\left(\mathrm{AB}\right)+\mathrm{\&xi;}(1+B^2-{4A}^2)]n$

${\mathrm{\&lambda;}}_1=\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[\frac{4}{3}-2(A^2+B^2)+\frac{\sin i}{1+\cos i}({\mathrm{AA}}_i+{\mathrm{BB}}_i)]n$

各式右端出现的a，e，i和n等量均为轨道根数

和

日、月的有关量β定义如下：

$\mathrm{\&beta;}=\left(\frac{\mathrm{GS}}{\mathrm{GM}}\right)/r^{\&prime;3}=\left(\frac{S}{M}\right)/r^{\&prime;3}=m^{\&prime;}/r^{\&prime;3}$

(其中M为地球质量，S为日或月质量，r′为日或月至地心的距离)

日、月位置量

可采用t0时刻的值。式中出现的几个量A，B，A_i，B_i，A_Ω，，B_Ω定义如下，其中的i′，Ω′，ω′，u′分别为日或月在惯性坐标系下的轨道根数：

$A=\frac{1}{4}\left\{(1-\cos i)\right[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$+(1+\cos i)[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\sin i\sin i^{\&prime;}[\cos \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\cos (\mathrm{\&omega;}+u^{\&prime;})]\}$

$B=-\frac{1}{4}\left\{(1-\cos i)\right[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$+(1+\cos i)[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\sin i\sin i^{\&prime;}[\sin \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\sin (\mathrm{\&omega;}+u^{\&prime;})]\}$

$A_i=\frac{\&PartialD;A}{\&PartialD;i}=\frac{1}{4}\left\{\sin i\right[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$-\sin i[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\cos i\sin i^{\&prime;}[\cos \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\cos (\mathrm{\&omega;}+u^{\&prime;})]\}$

$B_i=\frac{\&PartialD;B}{\&PartialD;i}=-\frac{1}{4}\left\{\sin i\right[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$-\sin i[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\cos i\sin i^{\&prime;}[\sin \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\sin (\mathrm{\&omega;}+u^{\&prime;})]\}$

$({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i=(A\frac{\&PartialD;A}{\&PartialD;\mathrm{\&Omega;}}+B\frac{\&PartialD;B}{\&PartialD;\mathrm{\&Omega;}})/\sin i$

$=-\frac{1}{8}\left\{\sin i\right[{(1+{\cos i}^{\&prime;})}^2\sin 2(\mathrm{\&theta;}-u^{\&prime;})+{(1-\cos i^{\&prime;})}^2\sin 2(\mathrm{\&theta;}+u^{\&prime;})+2{\sin }^2{i}^{\&prime;}\sin 2\mathrm{\&theta;}]$

$-2\cos i\sin i^{\&prime;}[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&theta;}-{2u}^{\&prime;})-(1-\cos i^{\&prime;})\sin (\mathrm{\&theta;}+{2u}^{\&prime;})-2\cos i^{\&prime;}\sin \mathrm{\&theta;}]\}$

短周期项

$a_s^{\left(1\right)}\left(t\right)=a\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(A^2+B^2-\frac{2}{3})({\left(\frac{r}{a}\right)}^2-1)+{\left(\frac{r}{a}\right)}^2((A^2-B^2)\cos 2f+\left(2\mathrm{AB}\right)\sin 2f)]$

式中变量含义同上。

Claims (4)

1.一种静止轨道卫星自主导航方法，其特征在于该方法的步骤为：

1)基于星敏感器和地球敏感器的kalman滤波算法，实时得到卫星相对于定点的位置偏差值；

2)对步骤1)得到的位置偏差值用最小二乘方法进行数据处理，获得一天内的轨道平面内的平均轨道根数：即半长轴、偏心率、近地点幅角和真近点角；对步骤1)得到的位置偏差值用平均滤波方法获得一天内的轨道平面外的平均轨道根数：即轨道倾角和升交点赤经；

3)以步骤2)得到的平均轨道根数作为星上轨道解析外推算法的输入值，外推计算一天内卫星轨道位置，提供连续导航定位信息，实现卫星自主导航功能。

2.根据权利要求1所述的一种静止轨道卫星自主导航方法，其特征在于：步骤1)的具体步骤为：

以此标称轨道定点位置建立Hill方程，即卫星相对标称轨道定点位置的相对运动方程，作为kalman滤波算法的状态方程，方程形式如下：

$\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}-2{\mathrm{\&omega;}}_0\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{a_x}{r}$

$\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&beta;}}+{\mathrm{\&omega;}}_0^2\mathrm{\&beta;}=\frac{a_y}{r}$

$\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}+2{\mathrm{\&omega;}}_0\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}-3{\mathrm{\&omega;}}_0^2\mathrm{\&gamma;}=\frac{a_z}{r}$

其中：

α为经度偏差，单位为rad；

β为纬度偏差，单位为rad；

γ为径向偏差，单位为rad；

f为静止轨道卫星理论半径，单位为km；

ω₀为标称静止轨道角速度，单位为rad/s；

a_x、a_y、a_z是外界加速度在标称轨道定点位置为原点的坐标系的分量；

标称轨道定点位置定义为：不考虑卫星漂移且轨道倾角为零的情况下的理论定点位置，与地球自转完全同步；

卫星自主导航测量要用到地球、卫星和标称轨道定点位置；其中，表示星敏感器给出的第一个方向矢量，

表示星敏感器给出的第二个方向矢量，

代表地球敏感器给出的地心方向矢量；

方向矢量

在卫星本体系下的坐标位置由星敏的安装矩阵决定；根据星敏感器的测量原理，其测量输出为方向矢量

在惯性坐标系中的方位坐标，分别记为：和

具体表达为 ${\left({\stackrel{\&RightArrow;}{S}}_1\right)}_I=\left[\begin{array}{ccc}{S}_{1\mathrm{ix}}& S_{1\mathrm{iy}}& S_{1\mathrm{iz}}\end{array}\right],$ ${\left({\stackrel{\&RightArrow;}{S}}_2\right)}_I=\left[\begin{array}{ccc}{S}_{2\mathrm{ix}}& S_{2\mathrm{iy}}& S_{2\mathrm{iz}}\end{array}\right];$

红外地球敏感器的测量输出为滚动角φ_h和俯仰角θ_h；根据测量得到的滚动角和俯仰角，经过两次旋转变换即可得到在本体坐标系下卫星轨道坐标系原点到地心的位置矢量，记为：

$\stackrel{\&RightArrow;}{E}=\left\{\begin{array}{c}{E}_{b1}=-\sin \left({\mathrm{\&theta;}}_h\right)\\ E_{b2}=\sin \left({\mathrm{\&phi;}}_h\right)\cos \left({\mathrm{\&theta;}}_h\right)\\ E_{b3}=\cos \left({\mathrm{\&phi;}}_h\right)\cos \left({\mathrm{\&theta;}}_h\right)\end{array}\right.$

从标称轨道定点位置的指向地心的单位矢量记为

在卫星本体坐标系得到两个夹角η₁、η₂：

$\cos {\mathrm{\&eta;}}_1=\stackrel{\&RightArrow;}{E}\&CenterDot;\stackrel{\&RightArrow;}{S_1}$

$\cos {\mathrm{\&eta;}}_2=\stackrel{\&RightArrow;}{E}\&CenterDot;\stackrel{\&RightArrow;}{S_2}$

在惯性坐标系下计算得到另外两个夹角η₁′、η₂′：

$\cos {\mathrm{\&eta;}}_1^{\&prime;}=\stackrel{\&RightArrow;}{Z}\&CenterDot;\stackrel{\&RightArrow;}{{\left(S_1\right)}_I}$

$\cos {\mathrm{\&eta;}}_2^{\&prime;}=\stackrel{\&RightArrow;}{Z}\&CenterDot;\stackrel{\&RightArrow;}{{\left(S_2\right)}_I}$

可得到测量方程g：

g＝(cosη₁-cosη′₁ cosη₂-cosη′₂)^T

为此构建一个自主导航的模型，其对象方程为：

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{bU}\\ g=h\left(x\right)\end{array}\right.$

其中：

$x={\left[\begin{array}{cccccc}\mathrm{\&alpha;}& \mathrm{\&beta;}& \mathrm{\&gamma;}& \stackrel{\&CenterDot;}{\mathrm{\&alpha;}}& \stackrel{\&CenterDot;}{\mathrm{\&beta;}}& \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}\end{array}\right]}^T$

$A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& {2\mathrm{\&omega;}}_0\\ 0& -{\mathrm{\&omega;}}_0^2& 0& 0& 0& 0\\ 0& 0& 3{\mathrm{\&omega;}}_0^2& -2{\mathrm{\&omega;}}_0& 0& 0\end{array}\right]$

b＝[0 0 0 1 1 1]^T

U为外部输入；

将测量方程在参考轨道处线性化：

$H\left(x\right)=\frac{\&PartialD;h\left(x\right)}{\&PartialD;x}=\left[\begin{array}{cccccc}{h}_{11}& h_{12}& 0& 0& 0& 0\\ h_{21}& h_{22}& 0& 0& 0& 0\end{array}\right]$

其中：

h₁₁＝-sin(α_T+α)cos(δ_T-β)S_1ix+cos(α_T+α)cos(δ_T-β)S_1iy

h₁₂＝cos(α_T+α)sin(δ_T-β)S_1ix+sin(α_T+α)sin(δ_T-β)S_1iy-cos(δ_T-β)S_1iz

h₂₁＝-sin(α_T+α)cos(δ_T-β)S_2ix+cos(α_T+α)cos(δ_T-β)S_2iy

h₂₂＝cos(α_T+α)sin(δ_T-β)S_2ix+sin(α_T+α)sin(δ_T-β)S_2iy-cos(δ_T-β)S_2iz

α_T为标称轨道定点位置的赤经，δ_T为标称轨道定点位置的赤纬

根据上述状态方程和测量方程，采用扩展Kalman滤波方法，可获得当前卫星轨道相对于标称轨道定点位置的偏移量，即：

$x={\left[\begin{array}{cccccc}\mathrm{\&alpha;}& \mathrm{\&beta;}& \mathrm{\&gamma;}& \stackrel{\&CenterDot;}{\mathrm{\&alpha;}}& \stackrel{\&CenterDot;}{\mathrm{\&beta;}}& \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}\end{array}\right]}^T.$

3.根据权利要求1所述的一种静止轨道卫星自主导航方法，其特征在于：步骤2)的具体步骤为：对步骤1)得到的位置偏差值用最小二乘方法进行数据处理，获得一天内的轨道平面内的平均轨道根数：即半长轴、偏心率、近地点幅角和真近点角；对步骤1)得到的位置偏差值用平均滤波方法获得一天内的轨道平面外的平均轨道根数：即轨道倾角和升交点赤经；

根据步骤1)中kalman滤波后获得的真实卫星相对于标称轨道定点位置的相对位置偏差值，为了获得便于轨道控制的控制变量以及提高导航精度，取滤波结果中的东西偏差α和南北偏差β，基于静止轨道卫星的受摄漂移原理，通过最小二乘法和平均滤波方法获得卫星的平均轨道六根数，包括：半长轴、偏心率、轨道倾角、升交点赤经、近地点幅角和真近点角，其中平均的意义为：利用一天的α和β数据拟合获得；

根据α和β，可以转换成卫星的赤经和赤纬，即赤经为α_T+α：，赤纬为δ_T-β，可以获得当前k时刻的地心单位矢量

其表达式为：

$\stackrel{\&RightArrow;}{i_k}=\left[\begin{array}{ccc}\cos ({\mathrm{\&alpha;}}_{\mathrm{Tk}}-{\mathrm{\&alpha;}}_k)\cos ({\mathrm{\&delta;}}_{\mathrm{Tk}}-{\mathrm{\&beta;}}_k)& \sin ({\mathrm{\&alpha;}}_{\mathrm{Tk}}-{\mathrm{\&alpha;}}_k)\cos ({\mathrm{\&delta;}}_{\mathrm{Tk}}-{\mathrm{\&beta;}}_k)& \sin ({\mathrm{\&delta;}}_{\mathrm{Tk}}-{\mathrm{\&beta;}}_k)\end{array}\right];$

由若干组这样的单位矢量即可确定卫星的轨道平面参数，即轨道倾角i和升交点赤经Ω，采用的方法为平均滤波方法；

在一个轨道周期T内，首先每隔设定时间dt，可获得一个地心单位矢量，分别记为：

(k＝1...N)，其中N＝T/dt，则通过叉乘可得一个倾角矢量：

$\stackrel{\&RightArrow;}{{i_1}^{\&prime;}}=\stackrel{\&RightArrow;}{i_1}\&times;\stackrel{\&RightArrow;}{i_2}...\stackrel{\&RightArrow;}{{i_k}^{\&prime;}}=\stackrel{\&RightArrow;}{i_k}\&times;\stackrel{\&RightArrow;}{i_{k+1}}(k<N)$

最终形成一个总的倾角矢量：

$\stackrel{\&RightArrow;}{i}=\frac{\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}\stackrel{\&RightArrow;}{{i_k}^{\&prime;}}}{N}$

如此在一个轨道周期内通过一系列倾角矢量得到一个总的倾角矢量，根据倾角矢量

在惯性坐标系中的XOY平面内投影定义：[sin i cos(Ω-90°)，sin i sin(Ω-90°)]，可解算出一组轨道倾角i和升交点赤经Ω；由于该倾角矢量是由一个轨道周期内的数据平均后获得，所以可以认为求出的i和Ω为上一个轨道周期内轨道倾角和升交点赤经的平均值，即一天数据可获得一组

和

采用递推最小二乘法对平面内轨道参数进行估计，即卫星东西方向运行规律可简化为：

α＝a₀+a₁t_orb+b₁cos(ω_et_orb)+b₂sin(ω_et_orb)

其中，t_orb一个轨道周期的时间变量，ω_e为静止轨道角速度，待估参数为a₀、a₁、b₁、b₂；

通过递推最小二乘法可实时估计参数a₀、a₁、b₁、b₂，由此在一个轨道周期结束时，可解算出卫星轨道平面内的平均轨道根数：

半长轴

可由a₁直接解出：其中r_s为静止轨道标称半径，a₁对应的单位为rad/天；

根据数学模型2esin(ω_et_orb+θ)，其中θ＝atan2(b₁，b₂)为相位偏移量，由b₁、b₂可解算出偏心率

即：

$\stackrel{\&OverBar;}{e}=0.5\sqrt{b_1^2+b_2^2}$

并且，令 ${\stackrel{\&RightArrow;}{i}}_{\mathrm{\&Omega;}}=\left[\begin{array}{ccc}\cos \left(\stackrel{\&OverBar;}{Q}\right)& \sin \left(\stackrel{\&OverBar;}{Q}\right)& 0\end{array}\right]$ 与当前的卫星地心单位矢量

进行内积，再求反余弦获得卫星幅角u，另外，真近点角

则：

$\stackrel{\&OverBar;}{\mathrm{\&omega;}}=\mathrm{\&alpha;}-\stackrel{\&OverBar;}{\mathrm{\&Omega;}}-\stackrel{\&OverBar;}{f}.$

4.根据权利要求1所述的一种静止轨道卫星自主导航方法，其特征在于：步骤3)的具体步骤为：以步骤2)得到的六个平均轨道根数作为星上轨道解析外推算法的输入值，外推计算一天内卫星轨道位置，提供连续导航定位信息，实现卫星自主导航功能；

由于地球同步轨道卫星偏心率e和轨道倾角i量级都比较小，故在进行解析外推计算时需采用适用于0≤e＜1和0≤i＜180°的第二类六个无奇点根数：

σ＝[a，h，k，ξ，η，λ]^T

其表达是为：a为半长轴，h＝sin i cosΩ，k＝sin i sinΩ，

$\mathrm{\&lambda;}=M+\widetilde{\mathrm{\&omega;}},$ 其中 $\widetilde{\mathrm{\&omega;}}=\mathrm{\&omega;}+\mathrm{\&Omega;};$

将步骤2)得到的六个平均轨道根数代入上式，令f≈M，可获得第二类无奇点根数的平均值： ${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_0={\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{a}}_0,& {\stackrel{\&OverBar;}{h}}_0,& {\stackrel{\&OverBar;}{k}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0,& {\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_0\end{array}\right]}^T,$ 其中的

即为作为星上轨道解析外推算法的初始拟平均根数，且根据步骤2)的结果每天更新一次，进而获得自主导航期间的实时轨道根数；

采用拟平均根数法构造的分析解可以获得瞬时轨道根数σ(t)：

$\mathrm{\&sigma;}\left(t\right)=\stackrel{\&OverBar;}{\mathrm{\&sigma;}}+{\mathrm{\&sigma;}}_s^{\left(1\right)}\left(t\right)$

$\stackrel{\&OverBar;}{\mathrm{\&sigma;}}\left(t\right)={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_0+(\mathrm{\&delta;}{\stackrel{\&OverBar;}{n}}_0+{\mathrm{\&sigma;}}_1)\left({t-t}_0\right)$

$\mathrm{\&delta;}{\stackrel{\&OverBar;}{n}}_0={\left(\begin{array}{cccccc}0& 0& 0& 0& 0& {\stackrel{\&OverBar;}{n}}_0\end{array}\right)}^T$

${\stackrel{\&OverBar;}{n}}_0={{\stackrel{\&OverBar;}{a}}_0}^{-3/2}$

其中：

σ(t)为瞬时轨道根数

为拟平均轨道根数

σ₁为一阶长期项系数，包含长期项和长周期项两部分

为一阶短周期项

为初始拟平均根数

t₀为初始拟平均根数所对应的轨道历元

为轨道平均角速度

地球非球形长期项系数σ₁：

a₁＝0

$h_1=-{\stackrel{\&OverBar;}{k}}_0{\mathrm{\&Omega;}}_1$

$k_1={\stackrel{\&OverBar;}{h}}_0{\mathrm{\&Omega;}}_1$

${\mathrm{\&xi;}}_1=-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0{\mathrm{\&omega;}}_1$

${\mathrm{\&eta;}}_1={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0{\mathrm{\&omega;}}_1$

λ₁＝M₁+ω₁+Ω₁

其中Ω₁，ω₁和M₁由下式表达：

$M_1=\frac{{3J}_2}{{2p}^2}n(1-\frac{3}{2}{\sin }^{2}i)\sqrt{1-e^2}$

${\mathrm{\&omega;}}_1=\frac{{3J}_2}{2p^2}n(2-\frac{5}{2}{\sin }^{2}i)$

${\mathrm{\&Omega;}}_1=-\frac{{3J}_2}{{2p}^2}n\cos i,$

上述各式中的a，e，i均为

且有：

$\left\{\begin{array}{c}\stackrel{\&OverBar;}{a}={\stackrel{\&OverBar;}{a}}_0,\stackrel{\&OverBar;}{e}={({\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_0^2+{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_0^2)}^{1/2},\\ {\sin }^2\stackrel{\&OverBar;}{i}=({\stackrel{\&OverBar;}{h}}_0^2+{\stackrel{\&OverBar;}{k}}_0^2),\cos i={[1-({\stackrel{\&OverBar;}{h}}_0^2+{\stackrel{\&OverBar;}{k}}_0^2)]}^{1/2}\\ \stackrel{\&OverBar;}{n}={\stackrel{\&OverBar;}{a}}_0^{-3/2},\stackrel{\&OverBar;}{p}=\stackrel{\&OverBar;}{a}(1-e^2)={\stackrel{\&OverBar;}{a}}_0(1-e_0^2)\end{array}\right.$

地球非球形短周期项

$a_s^{\left(1\right)}\left(t\right)=\frac{{3J}_2}{2a}\left\{(\frac{2}{3}-{\sin }^{2}i)\right[{\left(\frac{a}{r}\right)}^3-{(1-e^2)}^{-3/2}]+{\left(\frac{a}{r}\right)}^3[(h^2-k^2)\cos 2u+2\mathrm{hk}\sin 2u\left]\right\}$

${h_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{{4a}^2}\right)(h\cos 2u+k\sin 2u)$

${k_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{{4a}^2}\right)(h\sin 2u-k\cos 2u)$

${{\mathrm{\&lambda;}}_s}^{\left(1\right)}=\left(\frac{{21J}_2}{4a^2}\right)(\mathrm{\&xi;}\sin u-\mathrm{\&eta;}\cos u)$

${{\mathrm{\&xi;}}_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{4a^2}\right)[2(\cos u+\mathrm{\&xi;})+(\mathrm{\&xi;}\cos 2u+\mathrm{\&eta;}\sin 2u)]$

${{\mathrm{\&eta;}}_s}^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_2}{4a^2}\right)[2(\sin u+\mathrm{\&eta;})+(\mathrm{\&xi;}\sin 2u-\mathrm{\&eta;}\cos 2u)]$

上述短周期项各式右端出现的根数均为

而u和a/r由下列相应公式进行迭代计算获得：

$\left\{\begin{array}{c}\tilde{u}-\mathrm{\&lambda;}=\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u}\\ \frac{a}{r}={(1-\mathrm{\&xi;}\cos \tilde{u}-\mathrm{\&eta;}\sin \tilde{u})}^{-1}\\ \sin u=\left(\frac{a}{r}\right)[(\sin \tilde{u}-\mathrm{\&eta;})-\frac{\mathrm{\&xi;}}{1+\sqrt{1-e^2}}(\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u})]\\ \cos u=\left(\frac{a}{r}\right)[(\cos \tilde{u}-\mathrm{\&xi;})+\frac{\mathrm{\&eta;}}{1+\sqrt{{1-e}^2}}(\mathrm{\&xi;}\sin \tilde{u}-\mathrm{\&eta;}\cos \tilde{u})]\end{array}\right.$

其中

u＝f+ω，E是偏近点角；

地球非球形长周期项系数：

$a_l^{\left(1\right)}\left(t\right)=\left(\frac{{6J}_{\mathrm{2,2}}}{a}\right)\frac{\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})}{(1-\mathrm{\&alpha;})}$

$h_l^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_{\mathrm{2,2}}}{{2a}^2}\right)\frac{1}{1-\mathrm{\&alpha;}}[k\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})-h\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})]$

$k_l^{\left(1\right)}\left(t\right)=\left(\frac{{3J}_{\mathrm{2,2}}}{{2a}^2}\right)\frac{1}{1-\mathrm{\&alpha;}}[h\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})+k\cos (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})]$

${\mathrm{\&Delta;\&xi;}}_l^{\left(1\right)}\left(t\right)=O\left(e\right)$

${\mathrm{\&Delta;\&eta;}}_l^{\left(1\right)}\left(t\right)=O\left(e\right)$

${\mathrm{\&lambda;}}_l^{\left(1\right)}\left(t\right)=\left(\frac{{9J}_{\mathrm{2,2}}}{a^2}\right)\frac{\sin (2\mathrm{\&lambda;}-2\stackrel{\&OverBar;}{S})}{(1-\mathrm{\&alpha;})}$

上述各式中的a，h，k等量均为

各式中的谐系数J_2，2以及

等量的意义如下：

其中，n_e为标称轨道角速度，C_2，2，S_2，2是非归一化的谐系数，d为J2000.0起算的儒略日，S_G是格林尼治恒星时，λ_2，2＝-14.545°；

日月引力摄动部分：

长期项系数σ₁：

a₁＝0

$h_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i)\cos i\cos \mathrm{\&Omega;}+({\mathrm{AA}}_i+{\mathrm{BB}}_i)\sin \mathrm{\&Omega;}]n$

$k_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i)\cos i\sin \mathrm{\&Omega;}-({\mathrm{AA}}_i+{\mathrm{BB}}_i)\cos \mathrm{\&Omega;}]n$

${\mathrm{\&xi;}}_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[5\mathrm{\&xi;}\left(\mathrm{AB}\right)-\mathrm{\&eta;}(1+B^2-{4A}^2)]n$

${\mathrm{\&eta;}}_1=-\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[5\mathrm{\&eta;}\left(\mathrm{AB}\right)+\mathrm{\&xi;}(1+B^2-{4A}^2)]n$

${\mathrm{\&lambda;}}_1=\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[\frac{4}{3}-2(A^2+B^2)+\frac{\sin i}{1+\cos i}({\mathrm{AA}}_i+{\mathrm{BB}}_i)]n$

各式右端出现的a，e，i和n等量均为轨道根数

和

日、月的有关量β定义如下：

$\mathrm{\&beta;}=\left(\frac{\mathrm{GS}}{\mathrm{GM}}\right)/r^{\&prime;3}=\left(\frac{S}{M}\right)/r^{\&prime;3}=m^{\&prime;}/r^{\&prime;3}$

(其中M为地球质量，S为日或月质量，r′为日或月至地心的距离)

日、月位置量

可采用t0时刻的值；式中出现的几个量A，B，A_i，B_i，A_Ω，，B_Ω定义如下，其中的i′，Ω′，ω′，u′分别为日或月在惯性坐标系下的轨道根数：

$A=\frac{1}{4}\left\{(1-\cos i)\right[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$+(1+\cos i)[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\sin i\sin i^{\&prime;}[\cos \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\cos (\mathrm{\&omega;}+u^{\&prime;})]\}$

$B=-\frac{1}{4}\left\{(1-\cos i)\right[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$+(1+\cos i)[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\sin i\sin i^{\&prime;}[\sin \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\sin (\mathrm{\&omega;}+u^{\&prime;})]\}$

$A_i=\frac{\&PartialD;A}{\&PartialD;i}=\frac{1}{4}\left\{\sin i\right[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$-\sin i[(1+{\cos i}^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\cos (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\cos i\sin i^{\&prime;}[\cos \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\cos (\mathrm{\&omega;}+u^{\&prime;})]\}$

$B_i=\frac{\&PartialD;B}{\&PartialD;i}=-\frac{1}{4}\left\{\sin i\right[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}+u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}-\mathrm{\&theta;}-u^{\&prime;})]$

$-\sin i[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}-u^{\&prime;})+(1-\cos i^{\&prime;})\sin (\mathrm{\&omega;}+\mathrm{\&theta;}+u^{\&prime;})]$

$+2\cos i\sin i^{\&prime;}[\sin \left({\mathrm{\&omega;}-u}^{\&prime;}\right)-\sin (\mathrm{\&omega;}+u^{\&prime;})]\}$

$({\mathrm{AA}}_{\mathrm{\&Omega;}}+{\mathrm{BB}}_{\mathrm{\&Omega;}})/\sin i=(A\frac{\&PartialD;A}{\&PartialD;\mathrm{\&Omega;}}+B\frac{\&PartialD;B}{\&PartialD;\mathrm{\&Omega;}})/\sin i$

$=-\frac{1}{8}\left\{\sin i\right[{(1+{\cos i}^{\&prime;})}^2\sin 2(\mathrm{\&theta;}-u^{\&prime;})+{(1-\cos i^{\&prime;})}^2\sin 2(\mathrm{\&theta;}+u^{\&prime;})+2{\sin }^2{i}^{\&prime;}\sin 2\mathrm{\&theta;}]$

$-2\cos i\sin i^{\&prime;}[(1+{\cos i}^{\&prime;})\sin (\mathrm{\&theta;}-{2u}^{\&prime;})-(1-\cos i^{\&prime;})\sin (\mathrm{\&theta;}+{2u}^{\&prime;})-2\cos i^{\&prime;}\sin \mathrm{\&theta;}]\}$

短周期项

$a_s^{\left(1\right)}\left(t\right)=a\left(\frac{3}{2}{\mathrm{\&beta;a}}^3\right)[(A^2+B^2-\frac{2}{3})({\left(\frac{r}{a}\right)}^2-1)+{\left(\frac{r}{a}\right)}^2((A^2-B^2)\cos 2f+\left(2\mathrm{AB}\right)\sin 2f)]$

式中变量含义同上。

Images