一种应用星间跟踪的航天器星座分散化自主导航方法

技术领域

本发明涉及一种应用星间跟踪的航天器星座分散化自主导航方法，它是一种实现多航天 器星座分散化自主导航的信息处理方法。该方法在基于星间测量进行自主导航定轨的各种轨 道构型的多航天器任务中得到有效应用，经改进后亦可在星座星-地联合定轨系统中获得应 用。属于航天器自主导航技术领域。

背景技术

在当今科学探测、技术应用、乃至军事斗争中，航天活动正发挥着日益重要甚至不可替 代的作用。在各类航天计划中，使用多航天器构成一个整体空间系统的任务模式以其分布协 同化构型、多样灵活的功能组合、高任务效率和低风险等特点，可在更高的技术水平上满足 日益复杂多样的任务要求，是航天技术发展的重要趋势之一。目前已取得广泛成功的人造卫 星星座即属于多航天器任务的典型应用。包括通信卫星、导航卫星以及部分对地观测卫星均 采用了卫星组网构成星座的方式。

星座自主运行是指卫星在不依赖地面设施的情况下，自主确定星座状态和维持星座构型, 在轨完成飞行任务所要求的功能或操作。与以地面测控为主的传统模式相比，自主运行能大 大降低星座运行和管理成本、减小系统风险，是一种必然的发展趋势。自主导航为星座构形 控制提供测量数据，是卫星星座实现自主运行和控制的前提和基础，尤其对于导航星座来讲， 实现星座的自主导航不仅能够实现战时星座的自主生存，还肩负着为星座系统提供高精度广 播星历，从而提高用户的定位精度以及整个导航系统性能的重任。

从上世纪70年代开始,美国、俄罗斯和欧空局先后研究了多种卫星自主导航方案。目 前星座自主导航主要有两种技术途径：

（1）依靠单星自主导航实现星座自主轨道确定。这种方法依靠星座中每个单星独立完 成定轨，主要手段包括通过卫星导航定位或采用天文导航技术。前者实际上仍依赖GPS这 样的人造系统，严格地说不是完全自主的导航方式。后者则实现了完全自主导航测量，但目 前精度还相对较低。

（2）基于星间测量的星座自主导航。从原理上来说，基于星间测量的自主导航把卫星 星座成员成对地当作若干基线很长的重力测量仪，则星座成员相对运动变化体现的引力场信 息与绝对位置关联。通过测量星座成员卫星相互之间的相对运动状态，包括相对距离、相对 距离变化率和视线方位，可用来改进卫星的预报星历，从而提高星座整网导航定轨精度。美 国从20世纪80年代就开始研究GPS星座的自主运行问题,1984年Markley提出可以通过 测量星间矢量在惯性空间的投影确定两颗卫星的轨道，Ananda等随后公布了关于GPS自主 导航可行性的研究成果,1985年初美国空军空间系统部委托IBM开展关于自主导航算法的 深入研究。从2000年起,具备自主导航功能的GPS Block IIR系列进入全面测试阶段，其自 主导航的基本思想就是利用星间伪距测量数据,对地面控制中心注入的轨道预报数据进行 改进。但是截至目前有关GPS星座自主导航试验的具体数据仍未见披露。

在理论研究和实验方面，Psiaki指出，由于非中心引力的存在，卫星间的相对运动通过 绝对引力场的变化与位置相关，因此上述方法可以适用于各种地球卫星以及其他行星星座的 定轨。刘林、Hill等人的工作进一步表明，仅依靠星间相对测距进行自主导航的可观测性随 着星座成员卫星所处引力场不对称度的增大而提高。多个天体共同作用或具有较强非对称性 的引力场有利于绝对导航状态估计；反之，在引力场结构接近对称的情形下，仅仅依靠相对 测距只能构成空间的相对位置约束，无法测量星座的整体旋转。因此在星间测距基础上，陈 培提出加入基于星载多接收机载波相位的星间测向信息以提高导航估计性能；陈金平等提出 基于星敏感器测量卫星相对方位，进而确定星座相对惯性参考坐标系方位的轨道；熊凯则引 入X射线脉冲星观测获得了理论上更精确的星间方位信息。Yim等则表明仅根据星间方位 测量即可实现中心引力场中的完全自主定轨。基于星间测量的自主导航技术正呈现多种方案 同时发展的趋势，可以预见将成为星座自主导航的重要甚至首选方式。

作为基于星间测量的星座自主导航的关键技术，导航算法的设计必须考虑如下要点。首 先在导航信息来源方面，星间测量是多个航天器协同和并行的过程；第二，在系统构型上， 星座中航天器数目往往较多，而且星间链路可用性和拓扑结构具有时变的特点；第三，为了 完成导航信息的融合与分发，要求各成员航天器协同工作。

然而，由于导航状态估计算法结构的限制，当前的星座自主导航方案大多采用整网集中 定轨或群组分片集中定轨的方式，在每个群组中指定中心航天器，负责获取和存储群组各成 员航天器的观测信息，并调用批处理算法或卡尔曼滤波算法同时确定出群组中所有卫星的轨 道参数或导航状态。算法的集中化势必导致导航计算量和计算流程均集中于中心航天器，同 时增加了系统运行风险，而且不利于星座链接构型变化时的算法结构调整，亦不利于解决不 同航天器节点测量信息的非同步采样问题。随着星座成员数目的增加和构型多样化，上述问 题还会更加突出。

研究者们逐步认识到，采用分散化的算法方案是应对上述困难的有效途径。已提出的方 案将各成员航天器观测任务分散化，以顺序级联的方式进行全局状态的观测更新。虽然将观 测更新过程分解到相应的成员航天器中进行，但仍然将星座整网或群组看作一个整体进行导 航状态估计。相较于集中式算法，此类方法很好地处理了分布式观测的问题，但由于没有实 现滤波算法的彻底分散化，群组中相关的每个航天器需要依次对全局状态进行更新，存在单 个航天器计算量大、星间通信量反而有所增加、以及系统容错性能不高等缺点。

综上所述，分散化协同运行已成为基于星间测量的星座自主导航系统重要的发展趋势， 将构成其自主运行的关键环节，然而目前从算法结构上尚无完整实用的分散化方法。本发明 就是专门针对这一难点问题，基于星间测量和信息共享技术，提出在星间观测、状态估计、 故障检测以及系统重构等方面均按照分散化原则设计的自主导航系统和方法，实现系统功能 和算法运行的高度分散化，旨在为各类基于星间测量的星座自主导航系统提供一种有效的技 术方案。

发明内容

1、目的：

本发明针对航天器星座自主运行的需要，目的是提供一种应用星间跟踪的航天器星座分 散化自主导航方法。该方法可以较好地解决现有系统方案在算法结构上的不足。

2、技术方案：

一种应用星间跟踪的航天器星座分散化自主导航方法，方法实施的载体为由多个航天器 按照一定构型组成的星座。星座中每个航天器配置有星载计算机、星间相对距离测量设备、 星间相对速度测量设备、星间相对方位观测设备以及星间无线通信设备，具备进行导航计算、 星间测量及星间通信功能。每个航天器在星间网络中所处的位置平等，在计算功能上亦等同。 按照航天器与子滤波器一一对应的原则，将星座自主导航问题分解为对各成员航天器系统状 态的估计问题，各子滤波器负责对应航天器的导航估计。在星座整体导航滤波算法中对应一 个子滤波器。参见图1，该方法采用递推计算方式实现，记k(k=1,2,3...)为计算步序号，t_k为 对应的特征时刻，以一个计算更新周期[t_k,t_k+1]为例，该方法具体步骤如下：

步骤1：各子滤波器初始化；

步骤2：各子滤波器进行本地状态采样；

步骤3：各子滤波器进行时间更新；

步骤4：各航天器之间建立星间通信链路并保持跟踪。对于建立链路且相互跟踪成功航天器，

进入步骤5。对于未成功建立任何链路的航天器，进入步骤11；

步骤5：已建立星间链路的航天器进行星间跟踪观测。对于成功进行星间观测的航天器，根 据可用的星间观测确定本地相关观测模型，进入步骤6。对于未成功进行星间观测 的航天器，执行步骤11；

步骤6：经星间链路共享各子滤波器的状态采样信息；

步骤7：各子滤波器进行本地相关量测采样；

步骤8：各子滤波器进行量测更新；

步骤9：各子滤波器进行性能监控，判断其运行是否正常。若判断结果为正常，则执行步骤 10。否则执行步骤11；

步骤10：各子滤波器将步骤8的量测更新结果作为本地导航估计输出，返回步骤1，开始执 行下一个计算周期；

步骤11：各子滤波器将步骤3的时间更新结果作为本地导航估计输出，返回步骤1，开始执 行下一个计算周期。

其中，步骤1中所述的各子滤波器初始化，其实现方法为：

各子滤波器初始化是指确定各子滤波器在当前计算时刻t_k的本地系统状态估计初值及相应的误差协方差矩阵初值

对于整体算法的起始时刻，即t₀时刻，各子滤波器系统状态估计初值包括相应本地 航天器在惯性参照坐标系中的位置矢量估计初值和速度矢量估计初值

${\hat{X}}_0^{+}={\left(\begin{array}{c}{\hat{r}}_X\\ {\hat{v}}_X\end{array}\right)}_0---\left(1\right)$

设t₀时刻本地航天器的系统状态真实值为X₀，则相应的误差协方差矩阵初值按照 下式计算：

$P_{\mathrm{XX},0}^{+}=E\left\{\right[{\hat{X}}_0^{+}-X_0\left]{[{\hat{X}}_0^{+}-X_0]}^T\right\}---\left(2\right)$

若缺乏系统状态真实值X₀的必要信息，亦可根据工程经验确定

对于t_k(k=1,2,…)时刻，和则等于上一步计算时刻的估计输出。

其中，步骤2中所述的各子滤波器进行本地状态采样，其实现方法为：

各子滤波器依据t_k时刻本地状态估计初值及相应的误差协方差矩阵初值并行 地使用下面的对称采样算法计算相应的本地状态采样

$\left(\begin{array}{c}{\hat{X}}_k^{\left(0\right)}={\hat{X}}_k^{+}\\ {\hat{X}}_k^{\left(i\right)}={\hat{X}}_k^{+}+\sqrt{n+\tau }{\left(\sqrt{P_{\mathrm{XX},k}^{+}}\right)}_i^T\\ {\hat{X}}_k^{(n+i)}={\hat{X}}_k^{+}-\sqrt{n+\tau }{\left(\sqrt{P_{\mathrm{XX},k}^{+}}\right)}_i^T\end{array}\right),(i=1,...,n)---\left(3\right)$

其中采样矢量共2n+1个，带括号的上标表示采样矢量序号；n为系统状态维数；τ为状态 采样系数；当系统状态误差满足高斯分布时，选取n+τ=3。

其中，步骤3中所述的各子滤波器进行时间更新，其实现方法为：

首先定义子滤波器状态动力学模型f_x(·)。本发明主要关注以太阳系行星、矮行星或大 卫星为中心天体的航天器星座系统，状态动力学模型在相应中心天体惯性系中建立。与公式 （1）对应，导航系统状态矢量X包含航天器在相应中心天体惯性系下的位置矢量r以及速 度矢量v，导航系统状态动力学模型f_x(·)为：

$\stackrel{·}{X}=\left(\begin{array}{c}\stackrel{·}{r}\\ \stackrel{·}{v}\end{array}\right)=f_X(X,t)+w=\left(\begin{array}{c}v\\ a_{\mathrm{cen}}+a_{\mathrm{ns}}+a_{\mathrm{bd}}+a_{\mathrm{srp}}\end{array}\right)+w---\left(4\right)$

其中航天器受到中心天体质点引力加速度a_cen、中心天体非球形摄动加速度a_ns、太阳系主要 天体质点引力加速度a_bg、太阳光压摄动加速度a_srp以及加速度模型误差w影响。根据航天器 轨道动力学理论可完成各引力项的计算，模型误差建模为零均值高斯白噪声。

根据子滤波器状态动力学模型f_x(·)，建立相应的离散化状态模型F_x(·)：

$X_{k+1}=F_X(X_k,t_k)={\int }_{t_k}^{t_{k+1}}{f}_X(X_k,t_k)\mathrm{dt}---\left(5\right)$

接下来各子滤波器使用离散化状态模型F_x(·)并行地进行本地状态采样的时间更 新，得到t_k+1时刻各自的本地状态一步预测及本地状态误差协方差矩阵一步预测计算公式为：

$\left(\begin{array}{c}{\hat{X}}_{k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}{F}_X({\hat{X}}_k^{\left(j\right)},t_k)\\ P_{\mathrm{XX},k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-})}^T+Q_k\end{array}\right)---\left(6\right)$

其中j=0，...,2n；Q_k为系统状态模型噪声对应的协方差阵；W^(j)为状态采样权值，计算公式 为：

$\left(\begin{array}{c}{W}^{\left(0\right)}=\frac{\tau }{n+\tau }\\ W^{\left(i\right)}=\frac{1}{2(n+\tau )}\\ W^{(n+i)}=W^{\left(i\right)}\end{array}\right),(i=1,...,n)---\left(7\right)$

其中，步骤4中所述的各航天器之间建立星间通信链路并保持跟踪，其实现方法为：

星间通信链路的建立和保持通过星载空间通信及其链路捕获、跟踪、瞄准（ATP）系统 来完成。

首先利用各航天器星载通信发射端机产生星间通信信号，通过天线向满足可视条件的其 它航天器发射。后者使用天线和接收端机对星间通信信号进行捕获和确认，然后返回信标到 发射端，从而完成初步的链路锁定，建立通信链路。接下来发射端航天器根据目标航天器的 估计方位，驱动天线ATP伺服机构完成粗跟踪指向，然后提取通信信号的测角信息，导入 信号发射方向微调反馈控制回路，保持通信链路稳定精确指向。

其中，步骤5中所述的已建立星间链路的航天器进行星间跟踪观测，其实现方法为：

首先定义航天器星间跟踪观测量包括航天器间的相对距离，相对速度和导航计算坐标系 （本发明为中心天体惯性系）中的相对方位。

参见图2所示，以航天器A对航天器B的测量为例，假设其在惯性参考坐标系（记为i 系）中的位置矢量分别为和速度矢量分别为和相对视线矢量（即相对位置矢 量）为相对速度矢量为相对距离为ρ_AB，相对速度为相对方位单位矢量为

采用伪距式载波相位进行星间相对距离测量，利用无线电信号在空间定速传播的特性， 测量其发射时刻与接收时刻的时间差来确定相对距离：

ρ_AB＝cΔt_AB    (8)

式中c为电磁波传播速度，即光速；Δt_AB是测量信号的传播时间，由测距设备测定。

利用多普勒频移可测定相对速度，测量关系为

${\stackrel{·}{\rho }}_{\mathrm{AB}}=c·\frac{{\Phi }^2\cos \theta \pm \sqrt{1-{\Phi }^2{\sin }^2\theta }}{{\Phi }^2{\cos }^2\theta +1}---\left(9\right)$

其中Φ为目标天体辐射观测频率与实际频率之比；θ为星间视线方向与相对速度方向的夹 角；

设参考恒星视线方向为该数据由星历数据库给出，且在惯性空间中视线位置变化 极小。以其作为方位参考进行航天器间方位观测，使用恒星敏感器和星间跟踪测量设备，可 测量航天器间的相对视线方向与在惯性空间中的相对角偏差接下来按照下面的公 式可精确获得航天器间的相对方位在相应中心天体惯性系中的单位矢量。

$n_{\mathrm{AB}}^i=n_S^i+{\mathrm{\Delta n}}_{S,\mathrm{AB}}^i---\left(10\right)$

综合式（8）~（10），以航天器A对航天器B的测量为例，星间测量值包括：

$Z_{\mathrm{AB}}=\left(\begin{array}{c}{\rho }_{\mathrm{AB}}\\ {\stackrel{·}{\rho }}_{\mathrm{AB}}\\ n_{\mathrm{AB}}^i\end{array}\right)---\left(11\right)$

每一对可进行星间测量的航天器间均可得到一组星间测量值。对于某一个航天器，所有 与之相关的星间测量值组成其本地相关量测矢量Z_r，k+1。

其中，步骤6中所述的经星间链路共享各子滤波器的状态采样信息，其实现方法为：

经由星间链路，在各个测量相关的航天器间共享相应各子滤波器在步骤2中产生的状态 采样信息。对于每个子滤波器，在将本地状态采样上传至星间链路的同时，获得来自所 有与其存在星间测量的子滤波器的外部状态采样信息

其中，步骤7中所述的各子滤波器进行本地相关量测采样，其实现方法为：

首先定义观测模型。对于某个航天器对应的子滤波器，定义本地相关观测模型h_r(·)包 括该航天器和所有与其存在星间测量链路的航天器间的相对距离观测模型、相对速度观测模 型和相对方位观测模型。

根据步骤5中的变量定义，并参见图2，以航天器A和航天器B为例，每一个星间观测 量都至少同时与两个航天器的状态相关，相对距离观测模型为：

${\tilde{\rho }}_{\mathrm{AB}}=\left|\right|r_A^i-r_B^i\left|\right|+{ϵ}_{\rho ,\mathrm{AB}}---\left(12\right)$

其中“～”标记表示相应变量的测量值（下同），ε_ρ，AB表示航天器A和航天器B的相对距离测 量误差，包括测量时延、钟差以及随机误差。

相对速度观测模型可表示为速度矢量差在位置矢量差方向上的投影：

${\widetilde{\stackrel{·}{\rho }}}_{\mathrm{AB}}=(v_A^i-v_B^i)·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{ϵ}_{\stackrel{·}{\rho },\mathrm{AB}}---\left(13\right)$

其中表示航天器A和航天器B的相对速度测量误差。

相对方位观测模型则为：

${\tilde{n}}_{\mathrm{AB}}^m=C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{ϵ}_{n,\mathrm{AB}}---\left(14\right)$

其中表示i系向星间测量坐标系（m系）的姿态转换矩阵，由星载姿态确定系统测定； ε_n，AB表示航天器A和航天器B的相对方位测量误差。

记为航天器A和航天器B的星间观测矢量，式（12）~式（14）构成了一组星间观 测模型：

${\tilde{Z}}_{\mathrm{AB}}=h_{\mathrm{AB}}(X_A,X_B)+{ϵ}_{\mathrm{AB}}\stackrel{\Delta }{=}\left(\begin{array}{c}\left|\right|r_A^i-r_B^i\left|\right|\\ (v_A^i-v_B^i)·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\\ C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\end{array}\right)+\left(\begin{array}{c}{ϵ}_{\rho ,\mathrm{AB}}\\ {ϵ}_{\stackrel{·}{\rho },\mathrm{AB}}\\ {ϵ}_{n,\mathrm{AB}}\end{array}\right)---\left(15\right)$

对于某个航天器，完整的本地相关量测模型h_r(·)包括该航天器和所有与之存在星间测 量链路的航天器间的星间观测模型。

接下来，各子滤波器并行地使用各自的h_r(·)计算相应的本地相关量测采样矢量

${\hat{Z}}_{r,k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}{h}_r({\hat{X}}_k^{\left(j\right)},{\hat{Y}}_{r,k}^{\left(j\right)}),(j=0,...,2n)---\left(16\right)$

其中，步骤8中所述的各子滤波器进行量测更新，其实现方法为：

各子滤波器首先并行地计算相应的本地状态量测协方差矩阵P_XZr，h+1和本地量测协方差 矩阵P_ZrZr，k+1：

$\left(\begin{array}{c}{P}_{\mathrm{XZr},k+1}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})}^T\\ P_{\mathrm{ZrZr},k+1}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-}){({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-})}^T\end{array}\right)---\left(17\right)$

进而计算相应的增益矩阵K_k+1：

$K_{k+1}=P_{\mathrm{XZr},k+1}{P}_{\mathrm{Zr},k+1}^{-1}---\left(18\right)$

然后并行地计算t_k+1时刻各子滤波器相应的本地状态估计和本地状态估计误差协方 差阵

${\hat{X}}_{k+1}^{+}={\hat{X}}_{k+1}^{-}+K_{k+1}(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})---\left(19\right)$

$P_{\mathrm{XX},k+1}^{+}=P_{\mathrm{XX},k+1}^{-}-K_{k+1}{P}_{\mathrm{Zr},k+1}{K}_{k+1}^T---\left(20\right)$

其中，步骤9中所述的各子滤波器进行性能监控，判断滤波器运行是否正常，其实现方法为：

针对成员航天器可能出现测量或计算失效而造成算法故障的情况，承袭各成员航天器单 独估计自身状态的独立估计方式，每个子滤波器独立检测自身故障。故障检测算法采用基于 新息的经验卡方分布分析法，方法步骤如下。

首先通过下面的表达式计算t_k+1时刻的新息ε_k+1：

${ϵ}_{k+1}=Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-}---\left(21\right)$

然后定义下面的等价统计函数：

$\gamma =\frac{{ϵ}_{k+1}^T{P}_{\mathrm{XZr},k+1}^{-}{ϵ}_{k+1}}{l}---\left(22\right)$

式中，l为量测量的维数，统计量γ是最小值为零的非负数。理论上，如果滤波器模型准确， 且滤波未出现发散现象，γ将是一个标准的卡方分布。本算法中，设定γ的一个上限阈值γ_max作为滤波器发散的判据，当γ≤γ_max，则认为滤波器运行较好，且γ越小滤波性能越好；当 γ>γ_max，则认为滤波器出现故障。阈值的取值需要通过监视运行中的系统来确定，工程上 可以通过仿真实验按经验和需求取定上限阀值γ_max。

3、优点及功效：本发明的特点和优势在于：（1）与集中式UKF相比，本发明的分散化 算法利用了不同航天器状态解耦的性质，将集中式最优估计算法分模块并行运行，在数学本 质上与集中式算法等价，因此不会影响导航估计精度；（2）通过分散化的计算机制，合理地 平衡各成员航天器的导航计算负担，提高整体计算效率；（3）观测信息由相应航天器单独处 理，且不同航天器局部状态相互解耦，显著减少星间通信量；（4）算法结构不因成员航天器 数目和星间链路拓扑关系变化而改变，能够灵活应对星座构型的变化，亦可避免因航天器单 点失效引起整体导航计算失误的情况；（5）由于不同航天器状态解耦而观测相关的特性，便 于检测成员航天器的导航系统故障。总的来说，本发明在不牺牲导航精度的前提下显著提高 了星座自主导航算法的效率、并行性、灵活性以及容错性，为提升星座自主智能运行水平构 建了基础。

附图说明

图1为本发明的导航算法流程图。

图2为星间测量几何模型定义图。

图3（a）为航天器A、B、C集中式算法与本发明算法的位置估计误差对比图；

图3（b）为航天器D、E、F集中式算法与本发明算法的位置估计误差对比图；

图3（c）为航天器A、B、C集中式算法与本发明算法的速度估计误差对比图；

图3（d）为航天器D、E、F集中式算法与本发明算法的速度估计误差对比图。

图4（a）为本发明在星座构型变化时各航天器位置估计误差图；

图4（b）为本发明在星座构型变化时各航天器速度估计误差图。

图2中符号说明如下：

参见图2所示，在惯性直角坐标系O_iX_iY_iZ_i中，O_A航天器A的质心位置；O_B航天器B 的质心位置；表示航天器A在惯性参考坐标系（记为i系）中的位置矢量；表示航天 器B在惯性参考坐标系中的位置矢量；表示航天器A在惯性参考坐标系中的速度矢量； 表示航天器B在惯性参考坐标系中的速度矢量；表示航天器B相对于航天器A在惯 性参考坐标系中的相对位置矢量；表示航天器B相对于航天器A在惯性参考坐标系中 的相对速度矢量。

具体实施方式

下面结合附图和设定的仿真场景对本发明给予进一步说明。

见图1，一种应用星间跟踪的航天器星座分散化自主导航方法，以一个计算更新周期 [t_k,t_k+1]为例，具体方法步骤如下：

步骤1：各子滤波器初始化；

步骤2：各子滤波器进行本地状态采样；

步骤3：各子滤波器进行时间更新；

步骤4：各航天器之间建立星间通信链路并保持跟踪。对于建立链路且相互跟踪成功航天器， 进入步骤5。对于未成功建立任何链路的航天器，进入步骤11；

步骤5：已建立星间链路的航天器进行星间跟踪观测。对于成功进行星间观测的航天器，根 据可用的星间观测确定本地相关观测模型，进入步骤6。对于未成功进行星间观测 的航天器，执行步骤11；

步骤6：经星间链路共享各子滤波器的状态采样信息；

步骤7：各子滤波器进行本地相关量测采样；

步骤8：各子滤波器进行量测更新；

步骤9：各子滤波器进行性能监控，判断滤波器运行是否正常。若判断结果为正常，则执行 步骤10。否则执行步骤11；

步骤10：各子滤波器将步骤8的量测更新结果作为本地导航估计输出，返回步骤1，开始执 行下一个计算周期；

步骤11：各子滤波器将步骤3的时间更新结果作为本地导航估计输出，返回步骤1，开始执 行下一个计算周期。

其中，步骤1中所述各子滤波器初始化，其实现方法为：

各子滤波器初始化是指确定各子滤波器在当前计算时刻tk的本地系统状态估计初值 及相应的误差协方差矩阵初值

对于整体算法的起始时刻，即t₀时刻，各子滤波器系统状态估计初值包括相应的本 地航天器在惯性参照坐标系中的位置矢量估计初值和速度矢量估计初值

${\hat{X}}_0^{+}={\left(\begin{array}{c}{\hat{r}}_X\\ {\hat{v}}_X\end{array}\right)}_0---\left(1\right)$

设t₀时刻本地航天器的系统状态真实值为X₀，则状态估计误差协方差矩阵初值按 照下式计算：

$P_{\mathrm{XX},0}^{+}=E\left\{\right[{\hat{X}}_0^{+}-X_0\left]{[{\hat{X}}_0^{+}-X_0]}^T\right\}---\left(2\right)$

若缺乏系统状态真实值X₀的必要信息，亦可根据工程经验确定

对于t_k(k=1,2,…)时刻，和则等于上一步计算时刻的估计输出。

以包含24颗卫星的GPS星座为例，24个子滤波器将分别进行初始化，得到各自的状态 矢量估计初值和状态估计误差协方差矩阵初值。

其中，步骤2中所述各子滤波器进行本地状态采样，其实现方法为：

各子滤波器依据t_k时刻本地状态估计初值及相应的误差协方差矩阵初值并行 地使用下面的对称采样算法计算相应的本地状态采样

$\left(\begin{array}{c}{\hat{X}}_k^{\left(0\right)}={\hat{X}}_k^{+}\\ {\hat{X}}_k^{\left(i\right)}={\hat{X}}_k^{+}+\sqrt{n+\tau }{\left(\sqrt{P_{\mathrm{XX},k}^{+}}\right)}_i^T\\ {\hat{X}}_k^{(n+i)}={\hat{X}}_k^{+}-\sqrt{n+\tau }{\left(\sqrt{P_{\mathrm{XX},k}^{+}}\right)}_i^T\end{array}\right),(i=1,...,n)---\left(3\right)$

其中采样矢量共2n+1个，带括号的上标表示采样矢量序号；n为系统状态维数；τ为状态 采样系数；当系统状态误差满足高斯分布时，选取n+τ=3。

同样以包含24颗卫星的GPS星座为例，24个子滤波器将分别进行本地状态采样，得到 各自的状态采样矢量。本发明中n=6，因此每个子滤波器将产生13个状态采样矢量。

其中，步骤3中所述各子滤波器进行时间更新，其实现方法为：

首先定义子滤波器状态动力学模型f_x(·)。本发明主要关注以太阳系行星、矮行星或大 卫星为中心天体的航天器星座系统，状态动力学模型在相应中心天体惯性系中建立。与公式 （1）对应，导航系统状态矢量X包含航天器在相应中心天体惯性系下的位置矢量r以及速 度矢量v，导航系统状态动力学模型f_x(·)为：

$\stackrel{·}{X}=\left(\begin{array}{c}\stackrel{·}{r}\\ \stackrel{·}{v}\end{array}\right)=f_X(X,t)+w=\left(\begin{array}{c}v\\ a_{\mathrm{cen}}+a_{\mathrm{ns}}+a_{\mathrm{bd}}+a_{\mathrm{srp}}\end{array}\right)+w---\left(4\right)$

其中航天器受到中心天体质点引力加速度a_cen、中心天体非球形摄动加速度a_ns、太阳系主要 天体质点引力加速度a_bg、太阳光压摄动加速度a_srp以及加速度模型误差w影响。根据航天器 轨道动力学理论可完成各引力项的计算，模型误差建模为零均值高斯白噪声。

根据子滤波器状态动力学模型f_x(·)，建立相应的离散化状态模型F_x(·)：

$X_{k+1}=F_X(X_k,t_k)={\int }_{t_k}^{t_{k+1}}{f}_X(X_k,t_k)\mathrm{dt}---\left(5\right)$

接下来各子滤波器使用离散化状态模型F_x(·)并行地对本地状态采样进行时间更 新，得到t_k+1时刻各自的本地状态一步预测及本地状态误差协方差矩阵一步预测计算公式为：

$\left(\begin{array}{c}{\hat{X}}_{k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}{F}_X({\hat{X}}_k^{\left(j\right)},t_k)\\ P_{\mathrm{XX},k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-})}^T+Q_k\end{array}\right)---\left(6\right)$

其中j=0，...,2n；Q_k为系统状态模型噪声对应的协方差阵；W^(j)为状态采样权值，计算公式 为：

$\left(\begin{array}{c}{W}^{\left(0\right)}=\frac{\tau }{n+\tau }\\ W^{\left(i\right)}=\frac{1}{2(n+\tau )}\\ W^{(n+i)}=W^{\left(i\right)}\end{array}\right),(i=1,...,n)---\left(7\right)$

同样以包含24颗卫星的GPS星座为例，24个子滤波器将分别建立本地状态动力学模型， 并独立进行的状态时间更新。注意到，由于各子滤波器状态维数相同，因此按照公式（7）， 各子滤波器状态采样权值相同。

其中，步骤4中所述各航天器之间建立星间通信链路并保持跟踪，其实现方法为：

星间通信链路的建立和保持通过星载空间通信及其链路捕获、跟踪、瞄准（ATP）系统 来完成。

首先利用各航天器星载通信发射端机产生星间通信信号，通过天线向满足可视条件的其 它航天器发射。后者使用天线和接收端机对星间通信信号进行捕获和确认，然后返回信标到 发射端，从而完成初步的链路锁定，建立通信链路。接下来发射端航天器根据目标航天器的 估计方位，驱动天线ATP伺服机构完成粗跟踪指向，然后提取通信信号的测角信息，导入 信号发射方向微调反馈控制回路，保持通信链路稳定精确指向。

同样以包含24颗卫星的GPS星座为例，设每个卫星与其轨道面相邻且相位相差为1的 所有卫星分别建立星间链路，则每个卫星参与建立星间链路4条，整个星座包含星间链路共 48条。

其中，步骤5中所述已建立星间链路的航天器进行星间跟踪观测，其实现方法为：

参见图2所示，以航天器A对航天器B的测量为例，假设其在惯性参考坐标系（记为i 系）中的位置矢量分别为和速度矢量分别为和相对视线矢量（即相对位置矢 量）为相对速度矢量为相对距离为ρ_AB，相对速度为相对方位单位矢量为

采用伪距式载波相位进行星间相对距离测量，利用无线电信号在空间定速传播的特性， 测量其发射时刻与接收时刻的时间差来确定相对距离：

ρ_AB＝cΔt_AB    (8)

式中c为电磁波传播速度，即光速；Δt_AB是测量信号的传播时间，由测距设备测定。

利用多普勒频移可测定相对速度，测量关系为

${\stackrel{·}{\rho }}_{\mathrm{AB}}=c·\frac{{\Phi }^2\cos \theta \pm \sqrt{1-{\Phi }^2{\sin }^2\theta }}{{\Phi }^2{\cos }^2\theta +1}---\left(9\right)$

其中Φ为目标天体辐射观测频率与实际频率之比；θ为星间视线方向与相对速度方向的夹 角；

设参考恒星视线方向为该数据由星历数据库给出，且在惯性空间中视线位置变化 极小。以其作为方位参考进行航天器间方位观测，使用恒星敏感器和星间跟踪测量设备，可 测量航天器间的相对视线方向与在惯性空间中的相对角偏差接下来按照下面的公 式可精确获得航天器间的相对方位在相应中心天体惯性系中的单位矢量。

$n_{\mathrm{AB}}^i=n_S^i+{\mathrm{\Delta n}}_{S,\mathrm{AB}}^i---\left(10\right)$

综合式（8）~（10），以航天器A对航天器B的测量为例，星间测量值包括：

$Z_{\mathrm{AB}}=\left(\begin{array}{c}{\rho }_{\mathrm{AB}}\\ {\stackrel{·}{\rho }}_{\mathrm{AB}}\\ n_{\mathrm{AB}}^i\end{array}\right)---\left(11\right)$

每一对可进行星间测量的航天器间均可得到一组星间测量值。对于某一个航天器，所有 与之相关的星间测量值组成其本地相关量测矢量Z_r，k+1。

同样以包含24颗卫星的GPS星座为例，设每个卫星与其轨道面相邻且相位相差为1的 所有卫星分别建立星间链路，则每个卫星与其相邻的4颗卫星参与建立星间链路4条，可获 得4组星间测量值，整个星座共96组星间测量值。

其中，步骤6中所述经星间链路共享各子滤波器的状态采样信息，其实现方法为：

经由星间链路，在各个测量相关的航天器间共享相应各子滤波器在步骤2中产生的状态 采样信息。对于每个子滤波器，在将本地状态采样上传至星间链路的同时，获得来自所 有与其存在星间测量的子滤波器的外部状态采样信息

同样以包含24颗卫星的GPS星座为例，设每个卫星与其轨道面相邻且相位相差为1的 所有卫星分别建立星间链路，则每个卫星与其相邻的4颗卫星通过星间链路共享状态采样信 息。

其中，步骤7中所述各子滤波器进行本地相关量测采样，其实现方法为：

首先定义观测模型。对于某个航天器对应的子滤波器，定义本地相关观测模型h_r(·)包 括该航天器和所有与其存在星间测量链路的航天器间的相对距离观测模型、相对速度观测模 型和相对方位观测模型。

根据步骤5中的变量定义，并参见图2，以航天器A和航天器B为例，每一个星间观测 量都至少同时与两个航天器的状态相关，相对距离观测模型为：

${\tilde{\rho }}_{\mathrm{AB}}=\left|\right|r_A^i-r_B^i\left|\right|+{ϵ}_{\rho ,\mathrm{AB}}---\left(12\right)$

其中“～”标记表示相应变量的测量值（下同），ε_ρ，AB表示航天器A和航天器B的相对距离测 量误差，包括测量时延、钟差以及随机误差。

相对速度观测模型可表示为速度矢量差在位置矢量差方向上的投影：

${\widetilde{\stackrel{·}{\rho }}}_{\mathrm{AB}}=(v_A^i-v_B^i)·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{ϵ}_{\stackrel{·}{\rho },\mathrm{AB}}---\left(13\right)$

其中表示航天器A和航天器B的相对速度测量误差。

相对方位观测模型则为：

${\tilde{n}}_{\mathrm{AB}}^m=C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{ϵ}_{n,\mathrm{AB}}---\left(14\right)$

其中表示i系向星间测量坐标系（m系）的姿态转换矩阵，由星载姿态确定系统测定； ε_n，AB表示航天器A和航天器B的相对方位测量误差。

记为航天器A和航天器B的星间观测矢量，式（12）~式（14）构成了一组星间观 测模型：

${\tilde{Z}}_{\mathrm{AB}}=h_{\mathrm{AB}}(X_A,X_B)+{ϵ}_{\mathrm{AB}}\stackrel{\Delta }{=}\left(\begin{array}{c}\left|\right|r_A^i-r_B^i\left|\right|\\ (v_A^i-v_B^i)·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\\ C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\end{array}\right)+\left(\begin{array}{c}{ϵ}_{\rho ,\mathrm{AB}}\\ {ϵ}_{\stackrel{·}{\rho },\mathrm{AB}}\\ {ϵ}_{n,\mathrm{AB}}\end{array}\right)---\left(15\right)$

对于某个航天器，完整的本地相关量测模型h_r(·)包括该航天器和所有与之存在星间测 量链路的航天器间的星间观测模型。

接下来，各子滤波器并行地使用各自的h_r(·)计算相应的本地相关量测采样矢量

${\hat{Z}}_{r,k+1}^{-}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}{h}_r({\hat{X}}_k^{\left(j\right)},{\hat{Y}}_{r,k}^{\left(j\right)}),(j=0,...,2n)---\left(16\right)$

同样以包含24颗卫星的GPS星座为例，24个子滤波器将分别建立本地相关量测模型， 并独立进行本地相关量测采样。其中每个卫星参与建立星间链路4条，则相应的本地相关量 测模型涉及本地卫星和与之建立星间链路的4个卫星，包含4组相对距离观测模型、相对速 度观测模型和相对方位观测模型。

其中，步骤8中所述各子滤波器进行量测更新，其实现方法为：

各子滤波器首先并行地计算相应的本地状态量测协方差矩阵P_XZr，k+1和本地量测协方差 矩阵P_ZrZr，k+1：

$\left(\begin{array}{c}{P}_{\mathrm{XZr},k+1}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})}^T\\ P_{\mathrm{ZrZr},k+1}=\underset{j=0}{\overset{2n}{\Sigma }}{W}^{\left(j\right)}({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-}){({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-})}^T\end{array}\right)---\left(17\right)$

进而计算相应的增益矩阵K_k+1：

$K_{k+1}=P_{\mathrm{XZr},k+1}{P}_{\mathrm{Zr},k+1}^{-1}---\left(18\right)$

然后并行地计算t_k+1时刻各子滤波器相应的本地状态估计和本地状态估计误差协方 差阵

${\hat{X}}_{k+1}^{+}={\hat{X}}_{k+1}^{-}+K_{k+1}(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})---\left(19\right)$

$P_{\mathrm{XX},k+1}^{+}=P_{\mathrm{XX},k+1}^{-}-K_{k+1}{P}_{\mathrm{Zr},k+1}{K}_{k+1}^T---\left(20\right)$

同样以包含24颗卫星的GPS星座为例，24个子滤波器将分别进行量测更新。

其中，步骤9中所述各子滤波器进行性能监控，判断滤波器运行是否正常，其实现方法为：

针对成员航天器可能出现测量或计算失效而造成算法故障的情况，承袭各成员航天器单 独估计自身状态的独立估计方式，每个子滤波器独立检测自身故障。故障检测算法采用基于 新息的经验卡方分布分析法，方法步骤如下。

首先通过下面的表达式计算t_k+1时刻的新息ε_k+1：

${ϵ}_{k+1}=Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-}---\left(21\right)$

然后定义下面的等价统计函数：

$\gamma =\frac{{ϵ}_{k+1}^T{P}_{\mathrm{XZr},k+1}^{-}{ϵ}_{k+1}}{l}---\left(22\right)$

式中，l为量测量的维数，统计量γ是最小值为零的非负数。理论上，如果滤波器模型准确， 且滤波未出现发散现象，γ将是一个标准的卡方分布。本算法中，设定γ的一个上限阈值γ_max作为滤波器发散的判据，当γ≤γ_max，则认为滤波器运行较好，且γ越小滤波性能越好；当 γ>γ_max，则认为滤波器出现故障。阈值的取值需要通过监视运行中的系统来确定，工程上 可以通过仿真实验按经验和需求取定上限阈值γ_max。

同样以包含24颗卫星的GPS星座为例，24个子滤波器将分别计算新息和确定上限阈值， 然后独立进行性能监控。

使用以上方法进行数值仿真验算，仿真初始条件以GPS星座为参考设定，选取6个处 于不同轨道平面的GPS卫星进行自主导航计算仿真，GPS卫星PRN编号分别为07、25、29、 01、05和15，分别运行于GPS星座A、B、C、D、E和F轨道面。星间相对距离测量精度 设定为1m（1σ），相对速度测量精度设定为0.01m/s（1σ），相对方位测量精度设定为0.01° （1σ）。仿真时空基准选取J2000地球中心赤道惯性系，起始时刻设为2012年1月1日0 时（UTC）。仿真计算在MATLAB/Simulink环境中进行，数值积分算法采用4阶Runge-Kutta 法，更新步长设为5秒。

仿真设定两种场景。场景一为正常仿真模式，涉及所有6个卫星，每个卫星与前后各相 邻两个轨道面的共4个卫星建立星间链路，星座共形成12条星间链路，且各星间链路一直 保持正常跟踪和测量；场景二为星座构型变化模式，其初始设定与场景一相同。在20000秒 时刻，F轨卫星失效，星座成员由6个卫星变为5个卫星，所有与F轨卫星相关的星间链路 亦失效。在40000秒时刻，F轨卫星恢复，再次形成完整的6个卫星的构型。

图3（a）、图3（b）、图3（c）及图3（d）为仿真场景一情形下，传统集中式算法与本 发明提出的分散式算法的导航误差仿真对比图。其中图3（a）为航天器A、B、C的位置估 计误差对比图；图3（b）为航天器D、E、F的位置估计误差对比图；图3（c）为航天器A、 B、C的速度估计误差对比图；图3（d）为航天器D、E、F的速度估计误差对比图。从这 四个图可以看出，各航天器导航估计均可平稳收敛，两者精度相当。这验证了本发明提出的 分散式算法与传统集中式算法在本质上相一致的性质。

图4（a）和图4（b）为仿真场景二下，本发明算法的导航误差仿真结果图。其中图4 （a）为本发明在星座构型变化时各航天器位置估计误差图；图4（b）为本发明在星座构型 变化时各航天器速度估计误差图。从两个图可以看出，本发明设计的算法能够动态调整测量 信息以适应星座构型变化，从而保持导航估计的稳定。