CN102514734B - 基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法 - Google Patents

Abstract

本发明公开了一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，根据Halo轨道参数确定探测器布局及指向方案，由地面站测得探测器位置与速度信息，在地面计算机进行积分，得到n天内轨道信息，选取m个点的轨道信息上传至探测器的星载计算机；星载计算机根据探测器轨道信息，通过三次样条插值得到当前时刻轨道信息；根据星载计算机获得的探测器轨道信息，计算由惯性参考系到轨道坐标系的转换矩阵，最终结合姿态确定系统得到的绝对姿态角，得到探测器的相对姿态角；本发明的优点为：采用批量上传数据结合插值方法，为自主姿态确定提供参考基准，且计算量小；采用定期上注轨道信息的方式，数据传输的负荷相应小，减低数据传输系统的负担。

Description

基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法

技术领域

本发明涉及航天器设计领域，具体来说，是一种基于日地系统Halo轨道探测器的构型与姿态指向的姿态递推方法。

背景技术

平动点是太阳和地球的引力平衡点，运行于该点的探测器可以保持其位置而几乎不用消耗燃料。同时，日地系L1点位于日地之间具有观测太阳活动的天然优势，用来监测太阳风有约一个小时的预警时间；若将探测器直接放置在该点附近的Halo轨道，即可满足观测任务，还可防止通讯信号被太阳风湮没。

到目前为止，NASA和ESA已成功发射了6颗平动点任务探测器：InternationalSun-Earth Explorer 3(ISEE-3)：该计划由NASA和ESA合作实施，并于1978年8月12日经由Delta#144火箭发射到日地系L1点附近的Halo轨道，用来探测太阳风和宇宙射线，并于1982年6月结束该项任务并离开Halo轨道。随后，该探测器又进行了地球磁尾探测。之后，该探测器被重新命名为International CometaryExplorer(ICE)以进行Giacobini-Zinner和Halley彗星探测；ICE预计将于2014返回地球。Interplanetary Physics Laboratory(又被称为Wind)：于1994年11月1日由Delta-II发射至L1点附近用来研究太阳风对地磁的影响。Wind任务是迄今为止飞行轨迹最为复杂的任务，创造了多项纪录：到1997年共进行了38次月球重力辅助飞掠，首次完成月球Back-flip轨迹，首次到达大幅值顺行轨道等。SolarHeliospheric Observatory(SOHO)：该计划由ESA和NASA联合开发进行太阳观测的任务，于1995年12月2日由Atlas-II-AS发射至L1点的Halo轨道。SOHO曾因姿态控制系统故障，一度与地球失去联系。后经过营救，SOHO不仅完成预定任务，随后还进行了扩展任务。Advanced Composition Explorer(ACE)：该探测器于1997年8月25日由Delta-II-7920火箭发射升空，用来研究太阳风和太阳粒子，是第一个运行在真正意义的Lissajou轨道的探测器。由于跟踪系统故障，ACE的运行状态一度受到干扰，后成功排除。Microwave Anisotropy probe(MAP)：于2001年6月30日由Delta-II-7920火箭发射至L2点的Halo轨道，进行宇宙背景辐射研究。L2点远离太阳风及宇宙射线的干扰，是进行深空观测的理想场所，MAP是第一颗运行于L2点的探测器。Genesis：于2001年8月28日由Delta-7326火箭发射至L1点的Halo轨道，其主要任务是进行太阳风取样返回。Genesis于2001年11月16日达到Halo轨道，在运行轨道五圈后，于2004年4月离开Halo轨道返回地球。

可以说，和地球同步轨道一样，日地系平动点轨道是宝贵的空间资源，是人类共同的财富。如何更好地开发和利用这些资源，将是21世纪航天领域的重点课题。

平动点轨道与地球轨道有很大区别，为了计算得到平动点轨道，需要大量时间和跟踪数据。大多数平动点任务基于深空网络(DSN)的支持，或者世界空间网络(USN)等。先进技术采用Celestial Navigator(CelNav)，CelNav是一个星载Kalman滤波器，可以处理单通道前馈Doppler测量量和星载姿态敏感器数据。此外还有VeryLong Baseline Interferometry(VLBI)测量，称为Delta Differenced One-WayRange(DDOR)，它实际是一个从类星体附近到探测器的角度度量，目前已经应用到一些深空探测任务中，如适用于LPOs。由于Halo轨道目前不具有解析解，因而不能通过轨道方程设计递推轨道。

由于缺乏这些深空探测网络及先进星载设备的支持，目前国内要完成对于平动点轨道任务的自主导航尚存在很多亟需解决的技术难题，例如，如何解决姿态定向问题、如何解决远距离通讯的时间延迟及探测器能源供应等。对于轨道，可通过地面站指令进行修正调整，然而姿态变化的频率相对轨道要大得多，要实现高精度指向，必须及时调整姿态，因而不可能通过地面站测控实时注入姿态参考基准(即轨道坐标系信息)，必须要求探测器具有一定的自主定姿与控制能力。国内现有星载计算机无法独立根据地面注入的轨道测控信息对微分方程组进行积分，也无法通过解析的递推公式进行星上轨道递推，因而无法为姿态确定于控制系统提供姿态基准。

发明内容

为了解决上述问题，本发明提出了一种简单易行的姿态自主确定方法，通过地面站对测控数据处理、定期批量上注至探测器，再由星载计算机完成插值计算轨道参数与姿态基准，最后通过姿态敏感器数据与星上插值得到的姿态基准数据共同得到姿态信息，可大量减少星载计算机的计算负荷，数传的负荷也相应较小。

本发明一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，通过以下步骤来完成：

步骤1：定义Halo轨道坐标系、日心旋转坐标系、地心惯性坐标系、探测器轨道坐标系及探测器本体坐标系，并得到各个坐标系之间的转换矩阵；

设日心、地心、地日质心、探测器分别为S、E、M、E.S、S/C，且设R_k(γ)为绕k轴转角为γ的旋转矩阵，k＝x，y，z；

Halo轨道坐标系F中，原点取于地-日质心，x轴指向地球，z轴指向地-日旋转角速度方向；归一化单位：长度单位[L]＝L_S-E，L_S-E表示太阳与地球之间的距离；质量单位[M]＝m_S+m_E，m_S，m_E分别表示太阳、地球的质量；时间单位[T]＝(L_S-E ³/G(m_S+m_E))^1/2，表示地球绕太阳旋转周期，则地心在F中的坐标为(1-μ00)，日心在F中的坐标为(-μ00)，其中

日心旋转坐标系FR中，原点取在S；x_R指向E；z_R沿着地球绕太阳旋转角速度的方向，y_R轴满足笛卡尔右手法则；

地心惯性坐标系Fi中，原点取在E；z_i轴垂直于地球赤道平面指向北极；x_i轴指向春分点，y_i轴满足笛卡尔法则；

探测器轨道坐标系FO中，原点取在探测器的质心；x_o指向E；y_o垂直于太阳、探测器、地球所在平面，且与探测器速度方向成锐角，z_o轴满足笛卡尔右手法则；

探测器本体坐标系FB中，相对于Halo轨道坐标系，采用3-2-1的旋转顺序得到；

从Halo轨道坐标系F到日心旋转坐标系FR的转换关系如下式，(x y z)^T为探测器的位置矢量r在Halo轨道坐标系F中的坐标分量，(x_R y_R z_R)^T为探测器位置矢量r在FR中的坐标分量，且满足(x_R y_R z_R)^T＝(x+μ y z)^T；

$F^{\underset{\‾}{x+\mathrm{\μ}}}\→\mathrm{FR}$

设Halo轨道周期内轨道角变量α₁，α₂，探测器和太阳的连线与太阳和地球连线的夹角为α₁，探测器、地球、太阳所在平面与日心旋转坐标系FR的x_R与z_R张成的平面夹角为α₂，探测器与地球连线和地球与太阳连线的夹角为β₁，则由地心惯性坐标系Fi到探测器本体坐标系FB坐标转换关系以及从Halo轨道坐标系F到日心旋转坐标系FR的坐标转换关系如下：

其中，α_se为地球赤道平面与黄道面的夹角，且α_se＝23.43°；FR与Fi在t_orbit时刻的夹角α_G＝α_G0+ω_E-St_orbit，α_G0为FR与Fi的初始夹角，ω_E-S为地球绕太阳旋转的角速度；

θ_r，ψ_r分别为探测器的滚转、俯仰和偏航角，根据探测器位置矢量r在FR坐标系下的坐标分量，可得到：

${\mathrm{\α}}_1=\arccos \frac{x_R}{r_{S-S/C}}$

${\mathrm{\α}}_2=\arctan \frac{y_R}{z_R}$

${\mathrm{\β}}_1=\arccos \frac{1-x_R}{r_{S-S/C}}$

其中，r_S-S/C为日心到探测器的距离， $r_{S-S/C}=\sqrt{{(x+\mathrm{\μ})}^2+y^2+z^2};$

根据上述的坐标转换关系可知各坐标系之间的转换矩阵为：

Fi到FR的转换矩阵：R_Ri＝R_z(α_G+π)R_x(α_se)

FR到FO的转换矩阵：R_OR＝R_y(π/2-β₁)R_z(α₂-π)R_y(-π/2)

FO到FB的转换矩阵：

步骤2：根据Halo轨道参数确定探测器的布局及姿态指向；

对于探测器本体坐标系FB，其相对于探测器的轨道坐标系FO采用3-2-1的旋转顺序得到；+x_b，+y_b，z_b分别表示x_b，y_b，z_b轴的正方向，-x_b，-y_b，-z_b分别表示x_b，y_b，z_b轴的负方向，且+x_b，+y_b，+z_b面分别表示法线指向+x_b，+y_b，+z_b方向的面，-x_b，-y_b，-z_b面分别表示法线指向-x_b，-y_b，-z_b方向的面；探测器的数传天线固定安装在探测器本体的+x_b面，即数传天线的轴线垂直于探测器本体的+x_b面，且方向指向探测器本体的+x_b方向；探测器的太阳帆板固定安装于探测器本体的±y_b面，太阳帆板的轴线垂直于探测器本体的±y_b面，法线与-x_b轴夹角为θ，θ为太阳帆板安装角，θ＝(α+β)，其中，α和β分别为步骤1中所述Halo轨道周期内轨道角变量α₁β₁的平均值，

T为探测器的Halo轨道周期；所述探测器的+x_b轴保持对地定向，同时探测器太阳帆板的法线指向太阳；

步骤3：由卫星地面站测得探测器位置与速度信息，在卫星地面站计算机上进行积分得到n天内的探测器轨道信息；

由卫星地面站测得探测器t₀时刻的位置矢量

与速度矢量

在地心惯性坐标系Fi中的坐标分量(x_i y_i z_i)和分别用

表示；令时刻t＝t₀+Δt，Δt为时间增量，Δt∈[0，86400n]s，n为已设定的常数，表示卫星地面站遥测探测器位置与速度的时间周期，n的值影响探测器星上递推轨道的精度。通过在卫星地面站计算机进行积分得到：

$r_t^i=\underset{{r_{t_0}}^i}{\overset{{r_{t_0+\mathrm{\Δt}}}^i}{\∫}}{\stackrel{\·}{r}}^i\left(\begin{array}{ccc}{\stackrel{\·}{x}}_i& {\stackrel{\·}{y}}_i& {\stackrel{\·}{z}}_i\end{array}\right)\mathrm{dt},$

为探测器位置矢量相对间的一次导数；

${\stackrel{\·}{r}}_t^i=\underset{{{\stackrel{\·}{r}}_{t_0}}^i}{\overset{{{\stackrel{\·}{r}}_{t_0+\mathrm{\Δ}}}^i}{\∫}}{\stackrel{\·\·}{r}}^i\left(\begin{array}{ccc}{\stackrel{\·\·}{x}}_i& {\stackrel{\·\·}{y}}_i& {\stackrel{\·\·}{z}}_i\end{array}\right)\mathrm{dt},$

为探测器位置矢量对时间的二次导数；

其中，r_t ⁱ与

分别为在t∈[t₀，t₀+86400n]s时刻由卫星地面站计算机积分得到的探测器的位置与速度矢量在地心惯性坐标系Fi中的坐标分量，即探测器的轨道信息。

步骤4：由步骤3得到的n天内的探测器轨道信息中选择m个数据点的数据批量上传到星载计算机；

由步骤3得到的n天内的探测器轨道信息r_t ⁱ、中，等时间间隔δt选取m个点的轨道信息

p＝0，1，…，m-1，将其通过转换矩阵R_Ri，由Fi中的坐标分量转换为探测器的位置矢量

与速度矢量

在日心旋转坐标系中的坐标分量

分别为：

${r_{t_0+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{r_{t_0+\mathrm{p\δt}}}^i$

${{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^i$

将

和

p＝0，1，…，m-1，

数据打包，批量上传至探测器上的星载计算机。

步骤5：通过探测器上的星载计算机，通过三次养条插值，计算t时刻探测器的位置与速度在日心旋转坐标系中的坐标分量的递推值；

探测器上的星载计算机根据步骤4中得到的探测器的位置矢量与速度矢量

在日心旋转坐标系中的坐标分量

p＝0，1，…，m-1，

通过三次样条插值，得到时刻t的探测器的位置与速度矢量在日心旋转坐标系FR中的坐标分量的递推值

步骤6：通过探测器上的星载计算机获得探测器所在轨道当前时刻的

根据探测器上的星载计算机计算得到的当前时刻t的探测器轨道信息

t∈[t₀，t₀+86400n]s

由 ${\stackrel{\‾}{r}}_t^R={\left(\begin{array}{ccc}{\stackrel{\‾}{x}}_R& {\stackrel{\‾}{y}}_R& {\stackrel{\‾}{z}}_R\end{array}\right)}^T,$ ${\stackrel{\‾}{\stackrel{\·}{r}}}_t^R={\left(\begin{array}{ccc}{\stackrel{\‾}{\stackrel{\·}{x}}}_R& {\stackrel{\‾}{\stackrel{\·}{y}}}_R& {\stackrel{\‾}{\stackrel{\·}{z}}}_R\end{array}\right)}^T,$ 可得到：

${\stackrel{\‾}{\mathrm{\α}}}_2=\arctan \frac{{\stackrel{\‾}{y}}_R}{{\stackrel{\‾}{z}}_R}$

${\stackrel{\‾}{\mathrm{\β}}}_1=\arccos \frac{1-{\stackrel{\‾}{x}}_R}{{\stackrel{\‾}{r}}_{S-S/C}}$

其中，

为星载计算机计算得到的日心到探测器的距离，

为星载计算机插值得到的探测器位置

在FR坐标系中的坐标分量，

为星载计算机插值得到的探测器速度

在FR坐标系中的坐标分量，此处上角标T表示向量的转置；

步骤7：通过探测器上的姿态确定系统得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵R_Bi；

通过姿态确定系统得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵

其中，

θ_i，ψ_i为探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态角；

步骤8：通过探测器上的星载计算机得到探测器的相对姿态角

θ_r，ψ_r；

将步骤6中计算得到

带入步骤1中的FR到FO的转换矩阵R_OR中，可得到：

${\stackrel{\‾}{R}}_{\mathrm{OR}}=R_y(\mathrm{\π}/2-{\stackrel{\‾}{\mathrm{\β}}}_1)R_z({\stackrel{\‾}{\mathrm{\α}}}_2-\mathrm{\π})R_y(-\mathrm{\π}/2)=\left(\begin{array}{ccc}\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1\\ 0& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2& \sin {\stackrel{\‾}{\mathrm{\α}}}_2\\ -\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1\end{array}\right)$

通过步骤7得到的R_Bi，并根据步骤1得到的坐标转换关系R_Bi＝R_BOR_ORR_Ri，R_Oi＝R_ORR_Ri，通过公式R_BO＝R_Bi·R_Oi ^T，得到探测器轨道坐标系FO相对于探测器本体坐标系FB的姿态矩阵R_BO，由此通过解算姿态矩阵R_BO得到探测器的相对姿态角

θ_r，ψ_r，即FB相对于FO的姿态角为：

θ_r＝-arcsin(t₁₃ ^BO)

${\mathrm{\ψ}}_r=\arctan \left(\frac{{r_{12}}^{\mathrm{BO}}}{{r_{11}}^{\mathrm{BO}}}\right)$

其中，r_ij ^BO，(i＝1，2，3；j＝1，2，3)表示R_BO的第i行的第j个元素；

步骤9：判断时间t是否小于(t₀+86400n)s，如果是，重复步骤5～8，否则，重复步骤3～8，直至探测器寿命终结。

本发明的优点在于：

1、本发明采用批量上传数据结合插值方法，完成星上自主轨道参数确定，从而为自主姿态确定提供参考基准；且采用插值方法可以实现，且计算量小，；

2、本发明采用定期上注轨道信息的方式，数据传输的负荷相应较小，可以大大减低数据传输系统的负担。

附图说明

图1为本发明基于日地系统Halo轨道探测器构型与姿态指向的递推方法流程图；

图2为探测器构型空间示意图；

图3为探测器构型局部示意图；

图4为探测器在L1点的Halo轨道姿态指向；

图5为探测器在L2点Halo轨道上的姿态指向。

具体实施方式

下面将结合附图对本发明作进一步的详细说明。

本发明基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，如图1所示，通过以下步骤来完成：

步骤1：定义Halo轨道坐标系、日心旋转坐标系、地心惯性坐标系、探测器轨道坐标系及探测器本体坐标系，并得到各个坐标系之间的转换矩阵；

为方便说明，设日心、地心、地日质心、探测器分别为S、E、M、E.S、S/C，且设R_k(γ)为绕k轴转角为γ的旋转矩阵，k＝x，y，z。

Halo轨道坐标系F中，原点取于地-日质心，x轴指向地球，z轴指向地-日旋转角速度方向。归一化单位(用[AU]表示)，长度单位[L]＝L_S-E，L_S-E表示太阳与地球之间的距离；质量单位[M]＝m_S+m_E，m_S，m_E分别表示太阳、地球的质量；时间单位[T]＝(L_S-E ³/G(m_S+m_E))^1/2，表示地球绕太阳旋转周期，则地心在F中的坐标为(1-μ00)，日心在F中的坐标为(-μ00)，其中

日心旋转坐标系FR中，原点取在S；x_R指向E；z_R沿着地球绕太阳旋转角速度方向，y_R轴满足笛卡尔右手法则确定。

地心惯性坐标系Fi中，原点取在E；z_i轴垂直于地球赤道平面指向北极；x_i轴指向春分点，y_i轴满足笛卡尔法则确定。

探测器轨道坐标系FO中，原点取在探测器的质心；x_o指向E；y_o垂直于太阳、探测器、地球所在平面，且与探测器速度方向成锐角，z_o轴满足笛卡尔右手法则。

探测器本体坐标系FB中，相对于Halo轨道坐标系，采用3-2-1的旋转顺序定义。

从Halo轨道坐标系F到日心旋转坐标系FR的转换关系如下式，(x y z)^T为探测器的位置矢量r在Halo轨道坐标系F中的坐标分量，(x_R y_R z_R)^T为探测器位置矢量r在FR中的坐标分量，且它们满足(x_R y_R z_R)^T＝(x+μ y z)^T。

$F^{\underset{\‾}{x+\mathrm{\μ}}}\→\mathrm{FR}$

设Halo轨道周期内的轨道角变量α₁，α₂，β₁，探测器和太阳的连线与太阳和地球连线的夹角为α₁，探测器、地球、太阳所在平面与日心旋转坐标系FR的x_R与z_R张成的平面夹角为α₂，探测器与地球连线和地球与太阳连线的夹角为β₁，则由地心惯性坐标系Fi到探测器本体坐标系FB坐标转换关系以及从Halo轨道坐标系F到日心旋转坐标系FR的坐标转换关系如下：

其中，α_se为地球赤道平面与黄道面的夹角，且α_se＝23.43°；FR与Fi在t_orbit时刻的夹角α_G＝α_G0+ω_E-St_orbit，α_G0为FR与Fi的初始夹角，ω_E-S为地球绕太阳旋转的角速度；

θ_r，ψ_r分别为探测器的滚转、俯仰和偏航角(相对姿态角)，根据探测器位置矢量r在FR坐标系下的坐标分量，可得到：

${\mathrm{\α}}_1=\arccos \frac{x_R}{r_{S-S/C}}$

${\mathrm{\α}}_2=\arctan \frac{y_R}{z_R}$

${\mathrm{\β}}_1=\arccos \frac{1-x_R}{r_{S-S/C}}$

其中，r_S-S/C为日心到探测器的距离， $r_{S-S/C}=\sqrt{{(x+\mathrm{\μ})}^2+y^2+z^2}.$

根据上述的坐标转换关系可知各坐标系之间的转换矩阵为：

Fi到FR的转换矩阵：R_Ri＝R_z(α_G+π)R_x(α_se)

FR到FO的转换矩阵：R_OR＝R_y(π/2-β₁)R_z(α₂-π)R_y(-π/2)

FO到FB的转换矩阵：

步骤2：根据Halo轨道参数确定探测器的布局及姿态指向；

如图3所示，对于探测器本体坐标系FB，其相对于探测器的轨道坐标系FO采用3-2-1的旋转顺序得到。+x_b，+y_b，+z_b分别表示x_b，y_b，z_b轴的正方向，-x_b，-y_b，-z_b分别表示x_b，y_b，z_b，轴的负方向，且+x_b，+y_b，+z_b面分别表示法线指向+x_b，+y_b，+z_b方向的面，-x_b，-y_b，-z_b面分别表示法线指向-x_b，-y_b，-z_b方向的面。根据上述探测器本体坐标系FB与步骤1中定义的探测器轨道坐标系FO，本发明中将探测器的数传天线固定安装在探测器本体的+x_b面，即数传天线的轴线垂直于探测器本体的+x_b面，且方向指向探测器本体的+x_b方向。探测器的太阳帆板固定安装于探测器本体的±y_b面，太阳帆板的轴线垂直于探测器本体的±y_b面，法线与-x_b轴夹角为θ，θ为太阳帆板安装角。由于探测器在平动点Halo轨道上运行时，沿地球与太阳的连线方向的轨道浮动相对于垂直地球与太阳的连线方向平面内运动的尺度很小，因此通过上述探测器的构型设计，探测器的+x_b和+y_b不会受到太阳光照，探测器的+x_b和+y_b面可固定为探测器的散热面，由于散热面的固定，大大降低了探测器的热控设计难度。

在探测器运行期间，探测器在L1点与L2的Halo轨道上的探测器姿态指向方案相同，均为：探测器的+x_b轴保持对地定向，同时探测器太阳帆板的法线指向太阳。其中，探测器在L1点的Halo轨道姿态指向如图4所示，探测器在L2点Halo轨道上的姿态指向如图5所示。令θ＝(α+β)，其中，α和β分别为步骤1中Halo轨道周期内的轨道角变量α₁，β₁的平均值，

T为探测器的Halo轨道周期；而其中Halo轨道角变量α₁，β₁可由地面计算机根据徐明等的微分修正算法(徐明，徐世杰，地-月系平动点及Halo轨道的应用研究，宇航学报，第27卷第4期，2006年7月)得到的轨道位置信息计算得到。获得探测器所在的Halo轨道及相应的α和β，就可以确定太阳帆板2的安装角度θ。

步骤3：由卫星地面站测得探测器位置与速度信息，在卫星地面站计算机上进行积分得到n天内的探测器轨道信息；

由卫星地面站测得探测器t₀时刻的位置矢量

与速度矢量

在地心惯性坐标系Fi中的坐标分量(x_i y_i z_i)和

分别用

表示。令时刻t＝t₀+Δt，Δt为时间增量，Δt∈[0，86400n]s，n为已设定的常数，表示卫星地面站遥测探测器位置与速度的时间周期，n的值影响探测器星上递推轨道的精度。通过在卫星地面站计算机进行积分得到：

$r_t^i=\underset{{r_{t_0}}^i}{\overset{{r_{t_0+\mathrm{\Δt}}}^i}{\∫}}{\stackrel{\·}{r}}^i\left(\begin{array}{ccc}{\stackrel{\·}{x}}_i& {\stackrel{\·}{y}}_i& {\stackrel{\·}{z}}_i\end{array}\right)\mathrm{dt},$

为探测器位置矢量相对间的一次导数；

${\stackrel{\·}{r}}_t^i=\underset{{{\stackrel{\·}{r}}_{t_0}}^i}{\overset{{{\stackrel{\·}{r}}_{t_0+\mathrm{\Δ}}}^i}{\∫}}{\stackrel{\·\·}{r}}^i\left(\begin{array}{ccc}{\stackrel{\·\·}{x}}_i& {\stackrel{\·\·}{y}}_i& {\stackrel{\·\·}{z}}_i\end{array}\right)\mathrm{dt},$ 为探测器位置矢量对时间的二次导数；

其中，r_t ⁱ与

分别为在t∈[t₀，t₀+86400n]s时刻由卫星地面站计算机积分得到的探测器的位置与速度矢量在地心惯性坐标系Fi中的坐标分量，即探测器的轨道信息。

步骤4：由步骤3得到的n天内的探测器轨道信息中选择m个数据点的数据批量上传到星载计算机；

由步骤3得到的n天内的探测器轨道信息r_t ⁱ、

中，等时间间隔δt选取m个点的轨道信息

p＝0，1，…，m-1，

将其通过转换矩阵R_Ri，由Fi中的坐标分量转换为探测器的位置矢量

与速度矢量

在日心旋转坐标系中的坐标分量

${r_{t_0+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{r_{t_0+\mathrm{p\δt}}}^i$

${{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^i$

将

和

p＝0，1，…，m-1，

数据打包，批量上传至探测器上的星载计算机。其中，m的值同样会影响探测器星上递推轨道的精度，在保证插值精度的前提下，m应该尽量小，从而减小数据传输量，本发明中取m＝10，可满足工程需要。

步骤5：通过探测器上的星载计算机，通过三次养条插值，计算t时刻探测器的位置与速度在日心旋转坐标系中的坐标分量的递推值；

探测器上的星载计算机根据步骤4中得到的探测器的位置矢量

与速度矢量

在日心旋转坐标系中的坐标分量

以及通过三次样条插值，得到时刻t的探测器的位置与速度矢量在日心旋转坐标系FR中的坐标分量的递推值

本发明采用三次样条插值，是由于通过三次样条插值计算后的插值结果具有光滑性，在工程中应用较为合适；且由于步骤4中得到的

为函数值，为

的一阶导数值，由此可以满足三次样条插值的连续性条件和边界条件，保证三次样条插值函数唯一存在；且由于Halo轨道目前没有准确的解析形式，然而通过计算机仿真得到的插值误差表明，采用三次样条插值得到的递推轨道的精度满足工程需要的精度要求。

步骤6：通过探测器上的星载计算机获得探测器所在轨道当前时刻的

根据探测器上的星载计算机计算得到的当前时刻t的探测器轨道信息

t∈[t₀，t₀+86400n]s

由 ${\stackrel{\‾}{r}}_t^R={\left(\begin{array}{ccc}{\stackrel{\‾}{x}}_R& {\stackrel{\‾}{y}}_R& {\stackrel{\‾}{z}}_R\end{array}\right)}^T,$ ${\stackrel{\‾}{\stackrel{\·}{r}}}_t^R={\left(\begin{array}{ccc}{\stackrel{\‾}{\stackrel{\·}{x}}}_R& {\stackrel{\‾}{\stackrel{\·}{y}}}_R& {\stackrel{\‾}{\stackrel{\·}{z}}}_R\end{array}\right)}^T,$ 可得到：

${\stackrel{\‾}{\mathrm{\α}}}_2=\arctan \frac{{\stackrel{\‾}{y}}_R}{{\stackrel{\‾}{z}}_R}$

${\stackrel{\‾}{\mathrm{\β}}}_1=\arccos \frac{1-{\stackrel{\‾}{x}}_R}{{\stackrel{\‾}{r}}_{S-S/C}}$

其中，

为星载计算机得到的日心到探测器的距离，

为星载计算机插值得到的探测器位置

在FR坐标系中的坐标分量，

为星载计算机插值得到的探测器速度

在FR坐标系中的坐标分量，此处上角标T表示向量的转置。

步骤7：通过探测器上的姿态确定系统得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵R_Bi；

由于姿态确定系统由两个或两个以上恒星敏感器及绝对姿态确定算法组成，因此通过姿态确定系统可以得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵

其中

θ_i，ψ_i为探测器的绝对姿态角，即探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态角。

步骤8：通过探测器上的星载计算机得到探测器的相对姿态角

θ_r，ψ_r；

将步骤6中计算得到

带入步骤1中的FR到FO的转换矩阵R_OR中，可得到从FR到FO的坐标转换矩阵为：

${\stackrel{\‾}{R}}_{\mathrm{OR}}=R_y(\mathrm{\π}/2-{\stackrel{\‾}{\mathrm{\β}}}_1)R_z({\stackrel{\‾}{\mathrm{\α}}}_2-\mathrm{\π})R_y(-\mathrm{\π}/2)=\left(\begin{array}{ccc}\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1\\ 0& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2& \sin {\stackrel{\‾}{\mathrm{\α}}}_2\\ -\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1\end{array}\right)$

通过步骤7得到的R_Bi根据步骤1得到的坐标转换关系R_Bi＝R_BOR_ORR_Ri，R_Oi＝R_ORR_Ri，通过公式R_BO＝R_Bi·R_Oi ^T，得到探测器轨道坐标系FO相对于探测器本体坐标系FB的姿态矩阵R_BO，由此通过解算姿态矩阵R_BO得到探测器的相对姿态角

θ_r，ψ_r，即FB相对于FO的姿态角为：

θ_r＝-arcsin(r₁₃ ^BO)

${\mathrm{\ψ}}_r=\arctan \left(\frac{{r_{12}}^{\mathrm{BO}}}{{r_{11}}^{\mathrm{BO}}}\right)$

其中，r_ij ^BO，(i＝1，2，3；j＝1，2，3)表示R_BO的第i行的第j个元素。

步骤9：判断时间t是否小于(t₀+86400n)s，如果是，重复步骤5～8，否则，重复步骤3～8，直至探测器寿命终结。

通过上述方法，可完成星上自主轨道参数确定，从而为自主姿态确定提供参考基准，且计算量小；并且定期上注轨道信息的这种方式，数据传输的负荷相应较小，可以大大减低数据传输系统的负担。

上述步骤3中，探测器上的星载计算机通过三次样条插值所得到的探测器位置矢量的递推值与真实位置矢量的相对误差随上注间隔天数变化而变化。表1中给出了星上轨道确定误差，即星上插值所得轨道信息对姿态指向的影响，可以看出，姿态指向精度的相对误差与间隔天数n的选择有关。通过仿真分析，选择20天进行一次数据上注，误差较小且周期较大，且插值得到的轨道位置和速度与名义轨道位置和速度相比，相对误差为0.002％，对姿态指向精度的影响不大于0.234％。

表1星上插值所得轨道信息对姿态指向的影响

Claims (4)

1.一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，其特征在于：通过以下步骤来完成：

步骤1：定义Halo轨道坐标系、日心旋转坐标系、地心惯性坐标系、探测器轨道坐标系及探测器本体坐标系，并得到各个坐标系之间的转换矩阵；

设日心、地心、地日质心、探测器分别为S、E、E.S、S/C，且设R_k(γ)为绕k轴转角为γ的旋转矩阵，k＝x,y,z，具体为：

$R_x\left(\mathrm{\γ}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\γ}& \sin \mathrm{\γ}\\ 0& -\sin \mathrm{\γ}& \cos \mathrm{\γ}\end{array}\right],$ $R_y\left(\mathrm{\γ}\right)=\left[\begin{array}{ccc}\cos \mathrm{\γ}& 0& -\sin \mathrm{\γ}\\ 0& 1& 0\\ \sin \mathrm{\γ}& 0& \cos \mathrm{\γ}\end{array}\right],$ $R_z\left(\mathrm{\γ}\right)=\left[\begin{array}{ccc}\cos \mathrm{\γ}& \sin \mathrm{\γ}& 0\\ -\sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right];$

Halo轨道坐标系F中，原点取于地日质心，x轴指向地球，z轴指向地-日旋转角速度方向；归一化单位：长度单位[L]＝L_S-E，L_S-E表示太阳与地球之间的距离；质量单位[M]＝m_S+m_E，m_S,m_E分别表示太阳、地球的质量；时间单位[T]＝(L_S-E ³/G(m_S+m_E))^1/2，表示地球绕太阳旋转周期，则地心在F中的坐标为(1-μ 0 0)^T，日心在F中的坐标为(-μ 0 0)^T，其中

日心旋转坐标系FR中，原点取在S；x_R指向E；z_R沿着地球绕太阳旋转角速度方向，y_R轴满足笛卡尔右手法则；

地心惯性坐标系Fi中，原点取在E；z_i轴垂直于地球赤道平面指向北极；x_i轴指向春分点，y_i轴满足笛卡尔法则；

探测器轨道坐标系FO中，原点取在探测器的质心；x_o指向E；y_o垂直于太阳、探测器、地球所在平面，且与探测器速度方向成锐角，z_o轴满足笛卡尔右手法则；

探测器本体坐标系FB中，相对于Halo轨道坐标系，采用3-2-1的旋转顺序得到；

从Halo轨道坐标系F到日心旋转坐标系FR的转换关系如下式，(x y z)^T为探测器的位置矢量r在Halo轨道坐标系F中的坐标分量，(x_R y_R z_R)^T为探测器位置矢量r在FR中的坐标分量，且满足(x_R y_R z_R)^T＝(x+μ y z)^T；

$F\stackrel{x+\mathrm{\μ}}{\‾}\→\mathrm{FR}$

设Halo轨道周期内的轨道角变量α₁，α₂，探测器和太阳的连线与太阳和地球连线的夹角为α₁，探测器、地球、太阳所在平面与日心旋转坐标系FR的x_R与z_R张成的平面夹角为α₂，探测器与地球连线和地球与太阳连线的夹角为β₁，则由地心惯性坐标系Fi到探测器本体坐标系FB坐标转换关系以及从Halo轨道坐标系F到日心旋转坐标系FR的坐标转换关系如下：

其中，α_se为地球赤道面与黄道面的夹角，且α_se＝23.43°；FR与Fi在t_orbit时刻的夹角α_G＝α_G0+ω_E-St_orbit，α_G0为FR与Fi的初始夹角，ω_E-S为地球绕太阳旋转的角速度；分别为探测器的滚转、俯仰和偏航角，根据探测器位置矢量r在FR坐标系下的坐标分量，可得到：

${\mathrm{\α}}_1=\arccos \frac{x_R}{r_{S-S/C}}$

${\mathrm{\α}}_2=\arctan \frac{y_R}{z_R}$

${\mathrm{\β}}_1=\arccos \frac{1-x_R}{r_{S-S/C}}$

其中，r_S-S/C为日心到探测器的距离，

根据上述的坐标转换关系可知各坐标系之间的转换矩阵为：

Fi到FR的转换矩阵：R_Ri＝R_z(α_G+π)R_x(α_se)

FR到FO的转换矩阵：R_OR＝R_y(π/2-β₁)R_z(α₂-π)R_y(-π/2)

FO到FB的转换矩阵：

步骤2：根据Halo轨道参数确定探测器的布局及姿态指向；

对于探测器本体坐标系FB，其相对于探测器的轨道坐标系FO采用3-2-1的旋转顺序得到；+x_b,+y_b,+z_b分别表示x_b,y_b,z_b轴的正方向，-x_b,-y_b,-z_b分别表示x_b,y_b,z_b轴的负方向，且+x_b,+y_b,+z_b面分别表示法线指向+x_b,+y_b,+z_b方向的面，-x_b,-y_b,-z_b面分别表示法线指向-x_b,-y_b,-z_b方向的面；探测器的数传天线固定安装在探测器本体的+x_b面，即数传天线的轴线垂直于探测器本体的+x_b面，且方向指向探测器本体的+x_b方向；探测器的太阳帆板固定安装于探测器本体的±y_b面，太阳帆板的轴线垂直于探测器本体的±y_b面，法线与-x_b轴夹角为θ，θ为太阳帆板安装角，θ＝(α+β)，其中，α和β分别为步骤1中所述Halo轨道周期内轨道角变量α₁与β₁的平均值，

T为探测器的Halo轨道周期；所述探测器的+x_b轴保持对地定向，同时探测器太阳帆板的法线指向太阳；

步骤3：由卫星地面站测得探测器位置与速度信息，在卫星地面站计算机上进行积分得到n天内的探测器轨道信息；

由卫星地面站测得探测器t₀时刻的位置矢量

与速度矢量

在地心惯性坐标系Fi中的坐标分量(x_i y_i z_i)和

分别用

表示；令时刻t＝t₀+Δt，Δt为时间增量，Δt∈[0,86400n]s，n为已设定的常数，表示卫星地面站遥测探测器位置与速度的时间周期，n的值影响探测器星上递推轨道的精度；通过在卫星地面站计算机进行积分得到：

为探测器位置矢量相对间的一次导数；

为探测器位置矢量对时间的二次导数；

其中，r_t ⁱ与

分别为在t∈[t₀,t₀+86400n]s时刻由卫星地面站计算机积分得到的探测器的位置与速度矢量在地心惯性坐标系Fi中的坐标分量，即探测器的轨道信息；

步骤4：由步骤3得到的n天内的探测器轨道信息中选择m个数据点的数据批量上传到星载计算机；

由步骤3得到的n天内的探测器轨道信息r_t ⁱ、

中，等时间间隔δt选取m个点的轨道信息

p＝0,1,…,m-1，

将其通过转换矩阵R_Ri，由Fi中的坐标分量转换为探测器的位置矢量

与速度矢量

在日心旋转坐标系中的坐标分量分别为：

${r_{t_0+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{r_{t_0+\mathrm{p\δt}}}^i$

${{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^R=R_{\mathrm{Ri}}{{\stackrel{\·}{r}}_{t+\mathrm{p\δt}}}^i$

将

和

p＝0,1,…,m-1，

数据打包，批量上传至探测器上的星载计算机；

步骤5：通过探测器上的星载计算机，通过三次养条插值，计算t时刻探测器的位置与速度在日心旋转坐标系中的坐标分量的递推值；

探测器上的星载计算机根据步骤4中得到的探测器的位置矢量

与速度矢量

在日心旋转坐标系中的坐标分量

以及

通过三次样条插值，得到时刻t的探测器的位置与速度矢量在日心旋转坐标系FR中的坐标分量的递推值

步骤6：通过探测器上的星载计算机获得探测器所在轨道当前时刻的

根据探测器上的星载计算机计算得到的当前时刻t的探测器轨道信息

t∈[t₀,t₀+86400n]s

由 ${{\stackrel{\‾}{r}}_t}^R={\left(\begin{array}{ccc}{\stackrel{\‾}{x}}_R& {\stackrel{\‾}{y}}_R& {\stackrel{\‾}{z}}_R\end{array}\right)}^T,$ ${{\stackrel{\stackrel{\‾}{\·}}{r}}_t}^R={\left(\begin{array}{ccc}{\stackrel{\‾}{\stackrel{\·}{x}}}_R& {\stackrel{\‾}{\stackrel{\·}{y}}}_R& {\stackrel{\‾}{\stackrel{\·}{z}}}_R\end{array}\right)}^T,$ 可得到：

${\stackrel{\‾}{\mathrm{\α}}}_2=\arctan \frac{{\stackrel{\‾}{y}}_R}{{\stackrel{\‾}{z}}_R}$

${\stackrel{\‾}{\mathrm{\β}}}_1=\arccos \frac{1-{\stackrel{\‾}{x}}_R}{{\stackrel{\‾}{r}}_{S-S/C}}$

其中，

为星载计算机计算得到的日心到探测器的距离， ${\left(\begin{array}{ccc}{\stackrel{\‾}{x}}_R& {\stackrel{\‾}{y}}_R& {\stackrel{\‾}{z}}_R\end{array}\right)}^T$ 为星载计算机插值得到的探测器位置

在FR坐标系中的坐标分量， ${\left(\begin{array}{ccc}{\stackrel{\‾}{\stackrel{\·}{x}}}_R& {\stackrel{\‾}{\stackrel{\·}{y}}}_R& {\stackrel{\‾}{\stackrel{\·}{z}}}_R\end{array}\right)}^T$ 为星载计算机插值得到的探测器速度在FR坐标系中的坐标分量，此处上角标T表示向量的转置；

步骤7：通过探测器上的姿态确定系统得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵R_Bi；

通过姿态确定系统得到探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态矩阵其中，为探测器本体坐标系FB相对于地心惯性坐标系Fi的姿态角；

步骤8：通过探测器上的星载计算机得到探测器的相对姿态角

将步骤6中计算得到

带入步骤1中的FR到FO的转换矩阵R_OR中，可得到：

${\stackrel{\‾}{R}}_{\mathrm{OR}}=R_y(\mathrm{\π}/2-{\stackrel{\‾}{\mathrm{\β}}}_1)R_z({\stackrel{\‾}{\mathrm{\α}}}_2-\mathrm{\π})R_y(-\mathrm{\π}/2)=\left(\begin{array}{ccc}\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\sin {\stackrel{\‾}{\mathrm{\β}}}_1\\ 0& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2& \sin {\stackrel{\‾}{\mathrm{\α}}}_2\\ -\sin {\stackrel{\‾}{\mathrm{\β}}}_1& -\sin {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1& -\cos {\stackrel{\‾}{\mathrm{\α}}}_2\cos {\stackrel{\‾}{\mathrm{\β}}}_1\end{array}\right)$

通过步骤7得到的R_Bi，并根据步骤1得到的坐标转换关系R_Bi＝R_BOR_ORR_Ri，R_Oi＝R_ORR_Ri，通过公式R_BO＝R_Bi·R_Oi ^T，得到探测器轨道坐标系FO相对于探测器本体坐标系FB的姿态矩阵R_BO，由此通过解算姿态矩阵R_BO得到探测器的相对姿态角

即FB相对于FO的姿态角为：

θ_r＝-arcsin(r₁₃ ^BO)

${\mathrm{\ψ}}_r=\arctan \left(\frac{{r_{12}}^{\mathrm{BO}}}{{r_{11}}^{\mathrm{BO}}}\right)$

其中，r_ij ^BO,(i＝1,2,3;j＝1,2,3)表示R_BO的第i行的第j个元素；

步骤9：判断时间t是否小于(t₀+86400n)s，如果是，重复步骤5～8，否则，重复步骤3～8，直至探测器寿命终结。

2.如权利要求1所述一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，其特征在于：步骤2中，探测器的+x_b和+y_b面为探测器的散热面。

3.如权利要求1所述一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，其特征在于：步骤3中，取n=20。

4.如权利要求1所述一种基于日地系统Halo轨道探测器构型与姿态指向的姿态递推方法，其特征在于：步骤4中，取m=10。

Images