CN102230969A - 一种卫星星座星间链路的长时间自主维持方法 - Google Patents

Abstract

一种卫星星座星间链路的长时间自主维持方法，用于保持两颗卫星(主星和副星)的相对位置。该方法在副星上累积基于GPS测量的主星和副星的原始轨道点位数据，在完成数据累积后，使用解析形式的低阶轨道摄动模型对主星和副星在相同时间序列下的采样数据进行轨道改进，得到改进后的轨道平根数，再分别预报两颗卫星的位置，差分后得到卫星间相对位置。该方法主要应用于星间距离保持在100km左右的中低轨道星座的卫星中，长时间预报可达三天时间，三天中两星的绝对位置误差均可保持在800米以内、星间相对位置预报误差可保持在50米以内，能够满足航天任务中小卫星编队飞行的队形保持和星间数传天线的指向控制等领域的长时间不连续依赖于GPS的需求。

Description

一种卫星星座星间链路的长时间自主维持方法

技术领域

本发明涉及一种卫星星座星间链路的长时间自主维持方法，属于卫星通信领域。

背景技术

航天技术的发展，使得人类对太空的探索逐步的深入和细化，航天交会对接、小卫星编队、星间链路的建立、空间站的维护等任务使得相对测量技术成为航天工作者关注的热点。

以GPS为代表的全球卫星导航定位系统，具有全球覆盖、全天时、全天候工作的特性，可实时的为用户提供高精度时间位置速度等信息的能力。但是基于测量设备的导航系统的数据稳定性受到系统稳定性的限制，如GPS接收机的正常工作依赖于GPS系统的稳定性，其抗干扰能力差，已无法满足航天任务对导航定位和相对位置数据的长时间连续稳定输出的要求。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种卫星星座星间链路的长时间自主维持方法，通过在星上使用解析法定轨的轨道改进和预报模型对星间相对位置进行长时间预测。

本发明的技术解决方案是：

一种卫星星座星间链路的长时间自主维持方法，所述卫星星座包括主星和副星，步骤如下：

(1)通过主星和副星的GPS接收机分别实时获得主星和副星的原始轨道点位数据，所述原始轨道点位数据包括时间、位置和速度信息；

(2)在主星和副星的原始轨道点位数据同时有效的条件下，对所述主星和副星的原始轨道点位数据进行等间隔采样，采样间隔在30s到90s之间；

(3)当步骤(2)中得到的主星和副星的采样数据积累的时间超过一个轨道周期时，以当前采样时刻为参考时间，同时对主星和副星累积的采样数据进行轨道改进，分别得到改进后的主星和副星的轨道平根数；

(4)根据步骤(3)中得到的改进后的主星和副星的轨道平根数分别对主星和副星进行轨道预报并通过转换分别得到主星和副星的绝对位置；

(5)将步骤(4)中得到的主星和副星的绝对位置进行差分，得到主星和副星之间的相对位置，所述相对位置的模即为主星和副星之间的距离。

所述步骤(3)中对主星和副星累积的采样数据进行轨道改进通过如下步骤进行：

(2.1)解析形式的低阶轨道摄动模型采用地球重力场模型的二阶带谐项和不大于四阶的田谐项；

(2.2)在步骤(2.1)的条件下，以当前采样时刻t₀为参考时间，计算t₀对应的轨道平根数

并将

作为初始轨道参数；

(2.3)采用解析形式的低阶轨道摄动模型对参考时间t₀的轨道平根数

进行轨道外推，即轨道预报，得到任一时刻的轨道平根数

(2.4)将步骤(2.3)中得到的所述任一时刻的轨道平根数

转换成轨道瞬时根数σ_j后，再转换成轨道系下的位置

和速度，即为预报轨道点位数据；

(2.5)计算步骤(2.4)中得到的位置对轨道瞬时根数σ_j的偏导数矩阵；

(2.6)将所述原始轨道点位和所述预报轨道点位之差作为右端项、将一个轨道周期内的所有采样时刻的偏导数矩阵作为系数矩阵，建立条件方程组，将条件方程组法化之后，采用基于最小二乘估计原理的高斯法求解线性观测方程组，得到轨道平根数的增量，所述线性观测方程组即为建立的条件方程组；

(2.7)根据步骤(2.6)中得到的轨道平根数的增量，修正上一次的轨道平根数，完成轨道改进的一次迭代过程；之后按照步骤(2.3)～(2.6)的过程继续迭代，直到原始轨道点位与预报轨道点位之间的残差收敛之后，得到轨道改进后的参考时间t₀的轨道平根数

所述主星和副星运行在高度为300～1200km的轨道上。

所述主星和副星之间的距离在0～150km的范围内。

本发明与现有技术相比的有益效果是：

(1)本发明方法选用基于解析法定轨作为星间相对位置长时间预报的核心，通过GPS接收机实时获取主星和副星的原始轨道点位数据，通过其轨道改进模型在轨实时得到改进后的轨道平根数，以此分别对主副星进行长时间的轨道预报，将两者的差分作为相对位置；该方法完全在星上实现，不需要地面的支持，因此能够实现卫星编队和星间链路的自主维持；

(2)在理论上使用轨道摄动理论分析解析法定轨技术的预报模型，对其预报的两颗卫星(主星和副星)的星间相对位置误差进行理论分析，有效提高了星间相对位置长期预报的精度；

(3)本发明方法不但适用于基于GPS定位的卫星，还可以扩展到基于GLONASS、Galileo、北斗定位的卫星。

附图说明

图1为本发明方法的流程框图。

具体实施方式

本发明提供了一种卫星星座星间链路的长时间自主维持方法，所述卫星星座包括主星和副星，如图1所示，步骤如下：

(1)通过主星和副星的GPS接收机分别实时获得主星和副星的原始轨道点位数据，所述原始轨道点位数据包括时间、位置和速度信息；

(2)在主星和副星的原始轨道点位数据同时有效的条件下，对所述主星和副星的原始轨道点位数据进行等间隔采样，采样间隔在30s到90s之间；主星或者副星的原始轨道点位数据有可能会有毛刺，对于主星或者副星其一或者全部都有毛刺的情况，视为无效。若某一个原始轨道点位数据无效，但是刚好这个点是待采样的点位，那么后移1s采样。

(3)当步骤(2)中得到的主星和副星的采样数据积累的时间超过一个轨道周期时，以当前采样时刻为参考时间，同时对主星和副星累积的采样数据进行轨道改进，分别得到改进后的主星和副星的轨道平根数；主星和副星组成卫星星座，主星和副星的相对位置是近似固定的，因此主星和副星的轨道周期一致。

(4)根据步骤(3)中得到的改进后的主星和副星的轨道平根数分别对主星和副星进行轨道预报并通过转换分别得到主星和副星的绝对位置和绝对速度；

(5)将步骤(4)中得到的主星和副星的绝对位置进行差分，得到主星和副星之间的相对位置，所述相对位置的模即为主星和副星之间的距离。

所述步骤(3)中对主星和副星累积的采样数据进行轨道改进通过如下步骤进行：

(2.1)解析形式的低阶轨道摄动模型采用地球重力场模型的二阶带谐项和不大于四阶的田谐项；

(2.2)在步骤(2.1)的条件下，以当前采样时刻t₀为参考时间，计算t₀对应的轨道平根数

并将

作为初始轨道参数；

(2.3)采用解析形式的低阶轨道摄动模型对参考时间t₀的轨道平根数

进行轨道外推，即轨道预报，得到任一时刻的轨道平根数

(2.4)将步骤(2.3)中得到的所述任一时刻的轨道平根数

转换成轨道瞬时根数σ_j后，再转换成轨道系下的位置

和速度，即为预报轨道点位数据；

(2.5)计算步骤(2.4)中得到的位置对轨道瞬时根数σ_j的偏导数矩阵；

(2.6)将所述原始轨道点位和所述预报轨道点位之差作为右端项、将一个轨道周期内的所有采样时刻的偏导数矩阵作为系数矩阵，建立条件方程组，将条件方程组法化之后，采用基于最小二乘估计原理的高斯法求解线性观测方程组，得到轨道平根数的增量，所述线性观测方程组即为建立的条件方程组；

(2.7)根据步骤(2.6)中得到的轨道平根数的增量，修正上一次的轨道平根数，完成轨道改进的一次迭代过程；之后按照步骤(2.3)～(2.6)的过程继续迭代，直到原始轨道点位与预报轨道点位之间的残差收敛之后，得到轨道改进后的参考时间t₀的轨道平根数

所述主星和副星运行在高度为300～1200km的轨道上。

所述主星和副星之间的距离在0～150km的范围内。

首先对相对位置长期预报误差的相关性的原理进行分析。

基于轨道摄动模型的在轨长期预报方法主要分为积分法与分析法两种。积分法主要采用RK4方法和高精度的轨道摄动模型，其优点为可选取各种高精度的模型，且不受模型阶数影响，因而对单颗卫星位置预报精度比较高；但是计算量大，维持复杂，对计算设备要求也高，不利于在轨计算。针对单颗卫星位置预报精度，分析法满足在轨计算的设备要求，但其计算精度较低；然而当两颗卫星距离相近时，由于误差的相关性，在星间相对位置的计算中有较高的预报精度。

下面用分析法来分析星间相对位置误差。分析法的原理是基于小参数幂级数解的方法，其摄动方程可表示为：

$\stackrel{\stackrel{\·\·}{\→}}{r}={\stackrel{\→}{F}}_0+{\stackrel{\→}{F}}_{\mathrm{\ϵ}}$

其中：

为卫星运动加速度，

为卫星所受到的中心引力，

为卫星所受到的摄动力。

通过常数变易法处理，上式可转化为相应的摄动方程：

$\frac{d\stackrel{\→}{\mathrm{\θ}}}{\mathrm{dt}}=f_s(\stackrel{\→}{\mathrm{\θ}},t,\mathrm{\ϵ})$

其中：

为卫星轨道根数向量，f_s为六维向量函数，ε为摄动力和中心引力之比。

利用小参数幂级数解的方法解上述摄动方程即可得到解析解：

$\stackrel{\→}{\mathrm{\θ}}\left(t\right)={\stackrel{\→}{\mathrm{\θ}}}_0\left(t\right)+\underset{i=1}{\overset{6}{\mathrm{\Σ}}}\mathrm{\Δ}{\stackrel{\→}{\mathrm{\θ}}}_i(t,{\mathrm{\ϵ}}_i)$

对于较长弧段，在分析法中引入了拟平均根数法，原理是分别将长期摄动项和长周期摄动项引入了参考轨道。分析法通常只考虑到地球引力场的低阶带谐项，对于单颗卫星一天轨道预报，其精度在2.0公里左右。

对于星间相对位置，利用分析法我们可将卫星相对位置简单表示为：

$\mathrm{\Δ}\stackrel{\→}{r}={\stackrel{\→}{r}}_m-{\stackrel{\→}{r}}_s=[{\stackrel{\→}{r}}_m\left({\stackrel{\→}{\mathrm{\θ}}}_0\right)-{\stackrel{\→}{r}}_s\left({\stackrel{\→}{\mathrm{\θ}}}_0\right)]+[{\stackrel{\→}{r}}_m\left(\mathrm{\ϵ}\right)-{\stackrel{\→}{r}}_s\left(\mathrm{\ϵ}\right)]$

其中：

为卫星空间位置，下标m和s表示主副星，

为星间相对位置矢量。通过理论分析可知，由于卫星高度1100km一致、卫星之间的相对位置100km始终保持在1°(100km/7470km＝0.767°)以内，根据小参数摄动量级公式可知，卫星轨道根数摄动量级约为2ε·t，主要由1°左右的重力场空间差异导致，一天弧长相对位置预报精度大致可达到0.1公里((3×π/180)×2.0＝0.1)，远小于单颗卫星预报精度。

本发明是基于解析法定轨的星间相对位置长时间预报方法，相对于两个小时以内的短期预报而言，长时间预报的长度为三天。本发明直接使用解析法定轨模型对单星进行绝对位置长期预报，在此基础上对两星的点位作位置差分得到相对位置结果。

本发明方法的技术要点是采用基于解析法定轨模型在星上进行星间相对位置的长期预报。该方法主要适应于星间距离在0～150km的两星之间的中低轨道卫星的长时间的相对位置保持，此时两星所受的大气阻力基本相同，长期预报中相对位置的误差主要由两星绝对距离下的地球引力场差异决定，因为两星相距在0～150km，对应地球引力场相位相差在1°左右，通过理论分析可知，由此引起的长期相对位置误差在100m/天的量级。由于直接使用解析法定轨模型对单星进行绝对位置长期预报的误差在2km/天，说明长期保持星间距离在0～150km的两星作长时间相对位置预报时，之间的位置误差存在高度的相关性。

解析法定轨用于相对位置长期预报的核心如下：

(一)基本定义：

轨道计算过程中采用坐标系为轨道坐标系：以地球质心为坐标原点，瞬时真赤道面为参考平面，X轴指向2000.0平春分点在真赤道上的投影。

采用轨道根数为第一类无奇点根数σ＝(a，i，Ω，e，ω，M)^T：

a，i，Ω，ξ＝ecosω，η＝-esinω，λ＝ω+M

其中a，i，Ω，e，ω，M为开普勒根数，分别表示轨道半长径、轨道倾角、升交点赤经、轨道偏心率、近地点辐角、平近点角。轨道改进结果为指定时间时刻的轨道根数和卫星角速度的一阶变化率。轨道预报的结果为指定时刻的地固坐标系下的卫星坐标。

(二)地固系与轨道系之间的坐标转换：

WGS84坐标系是建立在WGS-84参考椭球体上的地固坐标系：以地球质心为坐标原点，与地心和国际习用原点CIO连线正交的平面为参考平面，X轴指向参考平面与格林尼治子午面的交线方向。

(注：WGS-84参考椭球体a＝6378136.3m，f＝1/298.257223563)

$\stackrel{\→}{r}=R_z(-{\mathrm{\θ}}_G){\left(\mathrm{EP}\right)}^T{\stackrel{\→}{r}}^{\′}$

其中：

$R_z(-{\mathrm{\θ}}_G)=\left(\begin{array}{ccc}\cos {\mathrm{\θ}}_G& -\sin {\mathrm{\θ}}_G& 0\\ \mathrm{ins}{\mathrm{\θ}}_G& \cos {\mathrm{\θ}}_G& 0\\ 0& 0& 1\end{array}\right)$

θ_G为轨道坐标系中的格林尼治恒星时角，GMST为格林尼治平恒星时，μ为从

对应有时刻到轨道坐标系时间时刻的赤经岁差；

θ_G＝GMST+μ＝280°.460618375+360°.9856122882×d

d为自2000年1月1日12时(对应的儒略日为2451545.0)起算的儒略日；

(EP)＝R_y(-x_p)·R_x(-y_p)为极移矩阵；

x_p，y_p分别为极移量。

由于极移量一般不超过0″.5，它的量级小于2.4×10-6(约15米)，而且在实际测量过程中无法得到实时的极移量，我们将忽略其影响。

(三)轨道系中位置速度与轨道根数之间的相互转换，(为公知)：

利用

计算某时刻的轨道根数σ：

$a=\frac{1}{\frac{2}{\left|\stackrel{\→}{r}\right|}-\stackrel{\stackrel{\·}{\→}}{r}\·\stackrel{\stackrel{\·}{\→}}{r}}$

$\stackrel{\→}{r}\×\stackrel{\stackrel{\·}{\→}}{r}=\sqrt{p}\stackrel{\→}{W}=\sqrt{p}\left(\begin{array}{c}\sin i\sin \mathrm{\Ω}\\ -\sin i\cos \mathrm{\Ω}\\ \cos i\end{array}\right)\⇒p=a(1-e^2),i,\mathrm{\Ω}$

$\stackrel{\→}{\mathrm{\Ω}}=\left(\begin{array}{c}\cos \mathrm{\Ω}\\ \sin \mathrm{\Ω}\\ 0\end{array}\right)$ ${\stackrel{\→}{\mathrm{\Ω}}}^{\′}=\left(\begin{array}{c}-\cos i\sin \mathrm{\Ω}\\ \cos i\cos \mathrm{\Ω}\\ \sin i\end{array}\right)$

$\begin{array}{cc}\cos u=\stackrel{\→}{r}\·\stackrel{\→}{\mathrm{\Ω}}/r& \sin u=\stackrel{\→}{r}\·{\stackrel{\→}{\mathrm{\Ω}}}^{\′}/r\\ \mathrm{\η}-\sin u=\stackrel{\stackrel{\·}{\→}}{r}\·\stackrel{\→}{\mathrm{\Ω}}\·\sqrt{p}& \mathrm{\ξ}+\cos u=\stackrel{\stackrel{\·}{\→}}{r}\·{\stackrel{\→}{\mathrm{\Ω}}}^{\′}\·\sqrt{p}\end{array}\⇒\mathrm{\ξ},\mathrm{\η},u$

$\mathrm{\ψ}=u-2{\tan }^{-1}\left(\frac{\mathrm{\ξ}\sin u+\mathrm{\η}\cos u}{1+\sqrt{1-{\mathrm{\ξ}}^2-{\mathrm{\η}}^2}+\mathrm{\ξ}\cos u-\mathrm{\η}\sin u}\right)$

λ＝ψ-ξsinψ-ηcosψ

利用某时刻的瞬时根数σ计算某时刻的卫星坐标、速度

$\stackrel{\→}{r}=r(\cos u\stackrel{\→}{\mathrm{\Ω}}+\sin u{\stackrel{\→}{\mathrm{\Ω}}}^{\′})$

$\stackrel{\stackrel{\·}{\→}}{r}=\frac{1}{\sqrt{p}}[(\mathrm{\η}-\sin u)\stackrel{\→}{\mathrm{\Ω}}+(\mathrm{\ξ}+\cos u){\stackrel{\→}{\mathrm{\Ω}}}^{\′}]$

r＝a(1-ξcosψ+ηsinψ)

$\cos u=\frac{a}{r}[\cos \mathrm{\ψ}-\mathrm{\ξ}+\mathrm{\α\η}(\mathrm{\λ}-\mathrm{\ψ})]$

$\sin u=\frac{a}{r}[\sin \mathrm{\ψ}+\mathrm{\η}+\mathrm{\α\ξ}(\mathrm{\λ}-\mathrm{\ψ})]$

ψ用迭代法求解：λ＝ψ-(ξsinψ+ηcosψ)

p＝a(1-ξ²-η²)

$\mathrm{\α}=\frac{1}{1+\sqrt{1-{\mathrm{\ξ}}^2-{\mathrm{\η}}^2}}$

(四)轨道系中轨道瞬时根数与轨道平根数之间的相互转换：

利用t₀时刻的瞬时根数σ₀计算t₀时刻的平根数

${\stackrel{\→}{\mathrm{\σ}}}_0={\mathrm{\σ}}_0-\mathrm{\Δ}{\mathrm{\σ}}_s$

σ₀为t₀时刻的第一类无奇点形式的轨道瞬时根数且σ₀＝(a，i，Ω，ξ，η，λ)^T，为已知量。

Δσ_s为一阶短周期项且Δσ_s＝(Δa_s，Δi_s，ΔΩ_s，Δξ_s，Δη_s，Δλ_s)^T，其计算公式如下(忽略J₂e～10^-6以及小于该量级的项)：

$\mathrm{\Δ}{a}_s=+\frac{3J_2}{2a}{\sin }^{2}i\cos 2\mathrm{\λ}$

$\mathrm{\Δ}{i}_s=+\frac{3J_2}{4a^2}\cos i\sin i\cos 2\mathrm{\λ}-\frac{3J_{22}}{2a^2\mathrm{\α}}\sin i\cos \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{22})\right]$

$-\frac{15J_{41}}{16a^4\mathrm{\α}}(4-7{\sin }^{2}i)\·\cos i\·\sin (\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{41})$

$+\frac{45J_{42}}{16a^4\mathrm{\α}}(6-7{\sin }^{2}i)\·\sin i\·\cos \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{42})\right]$

$+\frac{315J_{43}}{8a^4\mathrm{\α}}{\sin }^{2}i\·\cos i\·\sin \left[3(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{43})\right]$

$-\frac{315J_{44}}{8a^4\mathrm{\α}}{\sin }^{3}i\·\cos \left[4(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{44})\right]$

$\mathrm{\Δ}{\mathrm{\Ω}}_s=+\frac{3J_2}{4a^2}\cos i\sin 2\mathrm{\λ}+\frac{3J_{22}}{2a^2\mathrm{\α}}\cos i\sin \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{22})\right]$

$-\frac{15J_{41}}{16a^4\mathrm{\α}}(4-29{\sin }^{2}i+28{\sin }^{4}i)\cos (\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{41})/\sin i$

$-\frac{45J_{42}}{8a^4\mathrm{\α}}(3-7{\sin }^{2}i)\cos i\·\sin \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{42})\right]$

$+\frac{105J_{43}}{8a^4\mathrm{\α}}(3-4{\sin }^{2}i)\sin i\·\cos \left[3(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{43})\right]$

$+\frac{315J_{44}}{8a^4\mathrm{\α}}{\sin }^{2}i\·\cos i\·\sin \left[4(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{44})\right]$

$\mathrm{\Δ}{\mathrm{\ξ}}_s=+\frac{3J_2}{8a^2}(4-5{\sin }^{2}i)\cos \mathrm{\λ}+\frac{{7J}_2}{8a^2}{\sin }^{2}i\cos 3\mathrm{\λ}$

$+\frac{3J_{31}}{8a^3\mathrm{\α}}(4-15{\sin }^{2}i)\cos i\cos (\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{31})$

$+\frac{15J_{32}}{8a^3\mathrm{\α}}(2-3{\sin }^{2}i)\sin i\sin \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{32})\right]$

$-\frac{15J_{33}}{4a^2\mathrm{\α}}{\sin }^{2}i\cos i\cos \left[3(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{33})\right]$

$\mathrm{\Δ}{\mathrm{\η}}_s=-\frac{3J_2}{8a^2}(4-7{\sin }^{2}i)\sin \mathrm{\λ}-\frac{7J_2}{8a^2}{\sin }^{2}i\sin 3\mathrm{\λ}$

$+\frac{3J_{31}}{8a^3\mathrm{\α}}(4-5{\sin }^{2}i)\sin (\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{31})$

$-\frac{15J_{32}}{4a^3\mathrm{\α}}\sin i\cos i\cos \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{32})\right]$

$-\frac{15J_{33}}{4a^3\mathrm{\α}}{\sin }^{2}i\sin \left[3(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{33})\right]$

$\mathrm{\Δ}{\mathrm{\λ}}_s=-\frac{3J_2}{8a^2}(2-5{\sin }^{2}i)\sin 2\mathrm{\λ}+\frac{3J_{22}}{2a^2\mathrm{\α}}(4{\sin }^{2}i-1)\sin \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{22})\right]$

$+\frac{15J_{41}}{16a^4\mathrm{\α}}(4-69{\sin }^{2}i+98{\sin }^{4}i)\·\mathrm{ctgi}\·\cos (\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{41})$

$+\frac{45J_{42}}{16a^4\mathrm{\α}}(6-50{\sin }^{2}i+49{\sin }^{4}i)\·\sin \left[2(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{42})\right]$

$-\frac{105J_{43}}{8a^4\mathrm{\α}}(3-14{\sin }^{2}i)\·\sin i\·\cos i\·\cos \left[3(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{43})\right]$

$-\frac{315J_{43}}{16a^4\mathrm{\α}}(2-7{\sin }^{2}i)\·\mathrm{si}{n}^{2}i\·\sin \left[4(\mathrm{\Ω}-{\mathrm{\θ}}_G-{\mathrm{\λ}}_{44})\right]$

其中：为地球自转角速度n_e与卫星平均运动角速度n之比。J₂₂，J₃₁，J₃₂，J₃₃，J₄₁，J₄₂，J₄₃，J₄₄为田谐系数，J₂为带谐系数且J₂＝1.082636×10^-3，

J₂₂＝-1.77116×10^-6，λ₂₂＝-0.260468

J₃₁＝-2.17121×10^-6，λ₃₁＝0.1113511

J₃₂＝-3.70416×10^-7，λ₃₂＝-0.385165

J₃₃＝-2.23802×10^-7，λ₃₃＝0.420156

J₄₁＝-6.70506×10^-7，λ₄₁＝-2.392299

J₄₂＝-1.68650×10^-7，λ₄₂＝0.5493703

J₄₃＝-6.24996×10^-8，λ₄₃＝-0.038085

J₄₄＝-7.85530×10^-9，λ₄₄＝0.464363

(五)轨道平根数的预报：

利用t₀时刻的平根数

计算t_i时刻的平根数

${\stackrel{\‾}{\mathrm{\σ}}}_0={\stackrel{\‾}{\mathrm{\σ}}}_i+(\mathrm{\δ}+{\mathrm{\σ}}_1+{\mathrm{\σ}}_2)T$

其中σ₁，σ₂分别为一阶、二阶长期项系数，当σ为λ时，δ＝1；当σ为其它根数时，δ＝0。将公式展开如下：

$T=\stackrel{\‾}{n}(t_i-t_0)$

$\stackrel{\‾}{i}={\stackrel{\‾}{i}}_0+(i_{20}+i_{21}+i_{22})T$

$\stackrel{\‾}{\mathrm{\Ω}}={\stackrel{\‾}{\mathrm{\Ω}}}_0+({\mathrm{\Ω}}_1+{\mathrm{\Ω}}_{20}+{\mathrm{\Ω}}_{21}+{\mathrm{\Ω}}_{22})T$

$\stackrel{\‾}{\mathrm{\ξ}}={\stackrel{\‾}{\mathrm{\ξ}}}_0\cos \left({\mathrm{\ω}}_{1}T\right)+{\stackrel{\‾}{\mathrm{\η}}}_0\sin \left({\mathrm{\ω}}_{1}T\right)+({\mathrm{\ξ}}_{20}+{\mathrm{\ξ}}_{21}+{\mathrm{\ξ}}_{22})T$

$\stackrel{\‾}{\mathrm{\η}}={\stackrel{\‾}{\mathrm{\η}}}_0\cos \left({\mathrm{\ω}}_{1}T\right)+{\stackrel{\‾}{\mathrm{\ξ}}}_0\sin \left({\mathrm{\ω}}_{1}T\right)+({\mathrm{\η}}_{20}+{\mathrm{\η}}_{21}+{\mathrm{\η}}_{22})T$

$\stackrel{\‾}{\mathrm{\λ}}={\stackrel{\‾}{\mathrm{\λ}}}_0+(1+{\mathrm{\λ}}_1+{\mathrm{\λ}}_{20}+{\mathrm{\λ}}_{21}+{\mathrm{\λ}}_{22})T+\frac{1}{2}{\stackrel{\·}{\stackrel{\‾}{n}}}_0{(t_i-t_0)}^2$

$\stackrel{\‾}{n}={\stackrel{\‾}{n}}_0+{\stackrel{\·}{\stackrel{\‾}{n}}}_0(t_i-t_0)$

因此，在计算当中，需要计算

的时候，均可以

用

相应代替。

其中σ₁为一阶长期项系数，σ₂₀，σ₂₁，σ₂₂为地球形状摄动和大气阻力产生的二阶长期项系数，且σ₂＝σ₂₀+σ₂₁+σ₂₂；

一阶长期项系数σ₁的计算公式

${\mathrm{\Ω}}_1=-\frac{3J_2}{2{\stackrel{\‾}{p}}^2}\cos {\stackrel{\‾}{i}}_0$

${\mathrm{\ω}}_1=+\frac{3J_2}{2{\stackrel{\‾}{p}}^2}(2-\frac{5}{2}{\sin }^2{\stackrel{\‾}{i}}_0)$

${\mathrm{\λ}}_1=M_1+{\mathrm{\ω}}_1=+\frac{3J_2}{2{\stackrel{\‾}{p}}^2}(1-\frac{3}{2}{\sin }^2{\stackrel{\‾}{i}}_0)\sqrt{1-e^2}+{\mathrm{\ω}}_1$

二阶长期项系数σ₂的计算公式(忽略

J₃e，J₄e，F_sune～10^-9)

1.项部分σ₂₀

$\left\{\begin{array}{c}{a}_{20}=0\\ i_{20}=0\\ {\mathrm{\Ω}}_{20}=-\frac{3J_2^2}{8{\stackrel{\‾}{a}}^4}\cos {\stackrel{\‾}{i}}_0(15-19{\sin }^2{\stackrel{\‾}{i}}_0)\\ {\mathrm{\ξ}}_{20}=0\\ {\mathrm{\η}}_{20}=0\\ {\mathrm{\λ}}_{20}=+\frac{3J_2^2}{32{\stackrel{\‾}{a}}^4}(216-526{\sin }^2{\stackrel{\‾}{i}}_0+341{\sin }^4{\stackrel{\‾}{i}}_0)\end{array}\right.$

2.J_n(n≥3)部分σ₂₁

J₃＝-2.540×10^-6，J₄＝-1.619×10^-6

$\left\{\begin{array}{c}{a}_{21}=0\\ i_{21}=0\\ {\mathrm{\Ω}}_{21}=+\frac{15J_4}{16{\stackrel{\‾}{a}}^4}\cos {\stackrel{\‾}{i}}_0(4-7\mathrm{si}{n}^2{\stackrel{\‾}{i}}_0)\\ {\mathrm{\ξ}}_{21}=-\frac{3J_3}{8{\stackrel{\‾}{a}}^3}\sin {\stackrel{\‾}{i}}_0(4-5\mathrm{si}{n}^2{\stackrel{\‾}{i}}_0)\\ {\mathrm{\η}}_{21}=0\\ {\mathrm{\λ}}_{21}=-\frac{15J_4}{32{\stackrel{\‾}{a}}^4}(16-62{\sin }^2{\stackrel{\‾}{i}}_0+49{\sin }^4{\stackrel{\‾}{i}}_0)\end{array}\right.$

3.大气阻力摄动σ₂₂

$\left\{\begin{array}{c}\frac{1}{{\stackrel{\‾}{a}}_0}{a}_{22}=-{\stackrel{\‾}{a}}_0\left(C_D\frac{s}{m}{\mathrm{\ρ}}_{\mathrm{\Π}}\right)I_0^{*}\left(v\right)\\ e_{22}=-{\stackrel{\‾}{a}}_0\left(C_D\frac{s}{m}{\mathrm{\ρ}}_{\mathrm{\Π}}\right)I_1^{*}\left(v\right)\\ i_{22}={\mathrm{\Ω}}_{22}={\mathrm{\ω}}_{22}=0\\ M_{22}=-\frac{3T}{4{\stackrel{\‾}{a}}_0}{a}_{22}\\ {\mathrm{\ξ}}_{22}=+\cos {\stackrel{\‾}{\mathrm{\ω}}}_{1/2}{e}_{22}\\ {\mathrm{\η}}_{22}=-\sin {\stackrel{\‾}{\mathrm{\ω}}}_{1/2}{e}_{22}\\ {\mathrm{\λ}}_{22}=M_{22}=-\frac{3T}{4{\stackrel{\‾}{a}}_0}{a}_{22}\end{array}\right.$

其中：C_D是阻力系数，

是卫星的面质比(s是对阻力而言的有效截面积)，ρ_∏是卫星轨道近地点处的大气密度，H是该处的标高。

$v=\frac{\mathrm{ae}}{H}$ $I_m^{*}\left(v\right)=\frac{1}{\mathrm{\π}}{\∫}_0^{\mathrm{\π}}\cos \mathrm{mEexp}[-v(1-\cos E)]\mathrm{dE}$ T＝n(t-t₀)

(六)卫星的位置、速度对轨道根数的偏导数矩阵：

计算卫星的位置、速度对t₀时刻的轨道瞬时根数的偏导数

微系数

计算时忽略其中含偏心率e的项

$\left\{\begin{array}{c}\frac{\∂\stackrel{\→}{r}}{\∂i_0}=\stackrel{\→}{\mathrm{\Ω}}\×\stackrel{\→}{r}+{\stackrel{\→}{\mathrm{\Δ}}}_i\\ \frac{\∂\stackrel{\→}{r}}{\∂{\mathrm{\Ω}}_0}=\stackrel{\→}{k}\×\stackrel{\→}{r}\\ \frac{\∂\stackrel{\→}{r}}{\∂{\mathrm{\ξ}}_0}=\cos \left[\mathrm{\ω}(t-t_0)\right]\·({\mathrm{\ξ}}_1\stackrel{\→}{r}+{\mathrm{\ξ}}_2\stackrel{\stackrel{\·}{\→}}{r})-\sin \left[\mathrm{\ω}(t-t_0)\right]\·({\mathrm{\η}}_1\stackrel{\→}{r}+{\mathrm{\η}}_2\stackrel{\stackrel{\·}{\→}}{r})\\ \frac{\∂\stackrel{\→}{r}}{\∂{\mathrm{\η}}_0}=\cos \left[\mathrm{\ω}(t-t_0)\right]\·({\mathrm{\η}}_1\stackrel{\→}{r}+{\mathrm{\η}}_2\stackrel{\stackrel{\·}{\→}}{r})+\sin \left[\mathrm{\ω}(t-t_0)\right]\·({\mathrm{\ξ}}_1\stackrel{\→}{r}+{\mathrm{\ξ}}_2\stackrel{\stackrel{\·}{\→}}{r})\\ \frac{\∂\stackrel{\→}{r}}{\∂{\mathrm{\λ}}_0}=\frac{\stackrel{\stackrel{\·}{\→}}{r}}{n}\\ \frac{\∂\stackrel{\→}{r}}{\∂n_0}=\frac{1}{n}[-\frac{2}{3}\stackrel{\→}{r}+(t-t_0)\stackrel{\stackrel{\·}{\→}}{r}]\\ \frac{\∂\stackrel{\→}{r}}{\∂{\stackrel{\·}{n}}_0}=\frac{1}{n}[-\frac{2}{2}\stackrel{\→}{r}+\frac{(t-t_0)}{2}\stackrel{\stackrel{\·}{\→}}{r}](t-t_0)\end{array}\right.$

其中：

${\stackrel{\→}{\mathrm{\Ω}}}_1={\left(\begin{array}{ccc}-\sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\end{array}\right)}^T$ $\stackrel{\→}{k}={\left(\begin{array}{ccc}0& 0& 1\end{array}\right)}^T$

$J_2^{*}=\frac{3J_2}{2a^2},$ J₂＝1.082636×10^-3

$\left\{\begin{array}{c}{\mathrm{\Ω}}_i=J_2^{*}T\sin i\\ {\mathrm{\ω}}_i=-5\cos i{\mathrm{\Ω}}_i\\ M_i=-3\cos i{\mathrm{\Ω}}_i\end{array}\right.$

${\stackrel{\→}{\mathrm{\Δ}}}_i=-{\mathrm{\Ω}}_i(4-5{\sin }^{2}i)(\stackrel{\→}{k}\×\stackrel{\→}{r})-{\mathrm{\ω}}_i\sin i({\stackrel{\→}{\mathrm{\Ω}}}_1\×\stackrel{\→}{r})+M_i\left(\frac{\stackrel{\stackrel{\·}{\→}}{r}}{n}\right)$

ξ₁＝-cosu  η₁＝sinu

u＝ω+f；ω为近地点幅角，f为真近点角，n＝a^-3/2，

注：计算有摄偏导数中的

时可以按下式作近似计算：

为

的对时间的导数。

(七)条件方程：

建立的条件方程为：

${\stackrel{\→}{r}}_i-{\stackrel{\→}{r}}_{\mathrm{ic}}=\mathrm{\Σ}\frac{\∂{\stackrel{\→}{r}}_{\mathrm{ic}}}{\∂{\mathrm{\σ}}_0}\mathrm{\Δ}{\mathrm{\σ}}_0+{\stackrel{\→}{v}}_i,i=\mathrm{1,2},......,N$

其中，

为预报轨道点位数据中的位置数据，

为原始轨道点位数据中的位置数据，

为条件方程的残差，

收敛判别条件：

|δ_k|＜ε₁或|δ_k-δ_k-1|＜ε₂

k表示迭代的次数，ε₁，ε₂是给定的收敛判别标准。

(八)原始采样数据累积

a.按照主副星数据同时有效的时标对原始数据进行等间隔采样，在采样后的处理中按照单颗用户星进行描述，主星与副星的采样数据在后续的处理中是相同的；

b.在采样过程中，从GPS接收机原始观测数据包中解析出当前时间t_i的GPS时间和对应的用户星在地固系(WGS84)下的位置

和速度

c.将当前时间t_i的位置和速度从地固系转换到轨道系【参见(二)】；

d.采样存储时间时间和对应的轨道系下的位置；

e.将当前时间的轨道系下的位置和速度转换成第一类无奇点形式的轨道瞬时根数σ_i【参见(三)】；

f.将当前时间的轨道系下第一类无奇点形式的轨道瞬时根数σ_i转换成轨道

平根数

【参见(四)】；

g.在轨道改进后得到的轨道平根数有效的条件下，将轨道系下的轨道平根数外推到当前时间的地固系下的位置速度，与输入的通道板的位置速度进行比较：当位置单点超差时，将输入量作为野值剔除；当连续剔除30个时间后，认为轨道平根数异常，在数据累积后重新进行轨道改进获取新的轨道平根数。

(九)轨道改进

a.解析形式的低阶轨道摄动模型采用地球重力场模型的二阶带谐项和不大于四阶的田谐项；

b.以当前时刻t₀为参考时间，以当前时间对应的轨道平根数

为初始轨道参数；

c.采用解析形式的低阶轨道摄动模型对参考时间t₀的轨道平根数

进行轨

道外推，即轨道预报，得到任一时刻的轨道平根数【参见(五)】；

d.将得到的所述任一时刻的轨道平根数

转换成轨道瞬时根数σ_j后，再转换成轨道系下的位置

和速度，即为预报轨道点位数据【参见(三)】；

e.计算得到的位置

对轨道瞬时根数σ_j的偏导数矩阵【参见(六)】；

f.将所述原始轨道点位和所述预报轨道点位之差作为右端项、将一个轨道周期内的所有采样时刻的偏导数矩阵作为系数矩阵，建立条件方程组，将条件方程组法化之后，采用基于最小二乘估计原理的高斯法求解线性观测方程组，得到轨道平根数的增量，所述线性观测方程组即为建立的条件方程组【参见(七)】；

g.根据得到的轨道平根数的增量，修正上一次的轨道平根数，完成轨道改进的一次迭代过程；之后按照步骤c～f的过程继续迭代，直到原始轨道点位与预报轨道点位之间的残差收敛之后，得到轨道改进后的参考时间t₀的轨道平根数

(十)轨道长时间预报和星间差分

a.采用解析形式的低阶轨道摄动模型对轨道改进后的参考t₀时间的轨道改进后的轨道平根数

进行轨道外推，得到之后任一预报时刻t_k的轨道平

根数

【参见(五)】；

b.将任一预报时刻t_k的轨道平根数

转换成瞬时根数σ_k后【参见(四)】，

再转换成轨道系下的位置和速度【参见(三)】；

c.将任一预报时刻t_k的位置和速度从轨道系转换到地固系【参见(二)】，完成轨道预报；

d.计算出任一预报时刻t_k下的主星和副星的地固系下的位置和速度，两者之间的差分作为星间相对位置，从而完成星间相对位置的长时间预报。

上述流程如图1所示。

在以XX-8为背景工程实践中，使用解析法定轨的星间相对位置预报算法运行在以VC33为DSP的单片机上，VC33的主频为33MHz，以1Hz的频率进行数据处理时占用其30％的机时，算法耗时最大的轨道改进过程的一次迭代需要50秒，一般3～5次迭代即可完成轨道改进，最大不超过10分钟；用于存储主星和副星的采样数据的缓存区在SRAM中最大占用12k×32bit。由此可知该算法在嵌入式平台中占用的运行时间和存储空间的资源较小，有利于推广到小卫星编队飞行、星间链路维持等型号任务中。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (4)

1.一种卫星星座星间链路的长时间自主维持方法，所述卫星星座包括主星和副星，其特征在于步骤如下：

(1)通过主星和副星的GPS接收机分别实时获得主星和副星的原始轨道点位数据，所述原始轨道点位数据包括时间、位置和速度信息；

(2)在主星和副星的原始轨道点位数据同时有效的条件下，对所述主星和副星的原始轨道点位数据进行等间隔采样，采样间隔在30s到90s之间；

(3)当步骤(2)中得到的主星和副星的采样数据积累的时间超过一个轨道周期时，以当前采样时刻为参考时间，同时对主星和副星累积的采样数据进行轨道改进，分别得到改进后的主星和副星的轨道平根数；

(4)根据步骤(3)中得到的改进后的主星和副星的轨道平根数分别对主星和副星进行轨道预报并通过转换分别得到主星和副星的绝对位置；

(5)将步骤(4)中得到的主星和副星的绝对位置进行差分，得到主星和副星之间的相对位置，所述相对位置的模即为主星和副星之间的距离。

2.根据权利要求1所述的一种卫星星座星间链路的长时间自主维持方法，其特征在于：所述步骤(3)中对主星和副星累积的采样数据进行轨道改进通过如下步骤进行：

(2.1)解析形式的低阶轨道摄动模型采用地球重力场模型的二阶带谐项和不大于四阶的田谐项；

(2.2)在步骤(2.1)的条件下，以当前采样时刻t₀为参考时间，计算t₀对应的轨道平根数并将

作为初始轨道参数；

(2.3)采用解析形式的低阶轨道摄动模型对参考时间t₀的轨道平根数

进行轨道外推，即轨道预报，得到任一时刻的轨道平根数

(2.4)将步骤(2.3)中得到的所述任一时刻的轨道平根数

转换成轨道瞬时根数σ_j后，再转换成轨道系下的位置

和速度，即为预报轨道点位数据；

(2.5)计算步骤(2.4)中得到的位置

对轨道瞬时根数σ_j的偏导数矩阵；

(2.6)将所述原始轨道点位和所述预报轨道点位之差作为右端项、将一个轨道周期内的所有采样时刻的偏导数矩阵作为系数矩阵，建立条件方程组，将条件方程组法化之后，采用基于最小二乘估计原理的高斯法求解线性观测方程组，得到轨道平根数的增量，所述线性观测方程组即为建立的条件方程组；

(2.7)根据步骤(2.6)中得到的轨道平根数的增量，修正上一次的轨道平根数，完成轨道改进的一次迭代过程；之后按照步骤(2.3)～(2.6)的过程继续迭代，直到原始轨道点位与预报轨道点位之间的残差收敛之后，得到轨道改进后的参考时间t₀的轨道平根数

3.根据权利要求1所述的一种卫星星座星间链路的长时间自主维持方法，其特征在于：所述主星和副星运行在高度为300～1200km的轨道上。

4.根据权利要求1所述的一种卫星星座星间链路维持方法，其特征在于：所述主星和副星之间的距离在0～150km的范围内。

Images