航天器自主导航系统可以在少依赖甚至最终不依赖于地面系统支持的情况下，在轨实时 确定航天器的位置和速度，实现自主运行。对于卫星系统来说，自主导航有利于降低卫星对 地面的依赖程度，提高系统生存能力，例如战时，当地面系统遭到敌方的破坏和干扰时，仍 能完成轨道的确定和保持，这对军事卫星来讲，具有非常重要的意义。另外，自主导航还可 以有效地减轻地面测控站的负担，降低地面支持成本，从而降低整个航天计划的研制费用。
目前航天器的自主导航方法主要包括GPS，星间链路，磁力计，雷达和天文导航方法。 其中GPS和星间链路虽然可进行实时导航且达到较高的精度，但是航天器必须依赖其他航天 器提供量测信息，不能算是完全自主的导航方法。使用雷达高度计进行卫星的自主导航，受 到地形高度变化的影响，但是在海洋上方时能获得较好的导航性能。
天文导航是利用天体敏感器测得的天体(如行星、恒星、X射线脉冲星等)信息进行航 天器位置解算的一种完全自主的定位导航方法，可同时提供姿态、位置、速度以及时间信息； 适用于低轨卫星、高轨卫星和深空探测器，因而备受青睐。目前卫星的天文导航中，根据被 观测对象属性的不同，使用的敏感器主要有星敏感器、红外地平仪、紫外三轴敏感器以及X 射线探测器等。其中X射线探测器的观测对象为X射线源的脉冲辐射，获取脉冲到达时间信 息；使用X射线探测器能够为卫星提供精度较高的导航信息。对于多种天体测量敏感器的组 合，由于其被观测对象属性的差异，因而其导航属性存在互补的可能。
随着电子技术、导航技术和控制理论的发展，卫星上拥有众多的导航传感器模块，具备 了构成多传感器组合导航信息系统的可能条件。随着航空航天事业的发展，特别是现代化战 争对导航定位的精度和可靠性提出的更高要求，单一导航系统已难以满足。把两种或多种导 航系统组合起来，应用最优估计理论，形成最优组合导航系统，有利于充分运用各导航系统 的信息进行信息互补和信息合作，使组合后的导航系统在精度和可靠性方面都有所提高。同 时这也对卫星组合导航系统发生故障时的自动检测和故障隔离能力，提出了更高的要求。

本发明的目的在于：利用多个传感器进行组合，将各个传感器的测量数据进行信息融合 处理，提供了一种采用信息融合手段为卫星提供高精度自主导航参数的方法，以及提高卫星 在传感器发生故障时的自动检测能力和故障隔离能力，重构导航系统的方法。
本发明所采用的技术方案是：基于多信息融合的卫星组合导航系统，利用X射线探测器、 星敏感器、红外地平仪、雷达高度计和紫外三轴敏感器五种传感器以及卫星的轨道动力学模 型构成了四个卫星自主导航子系统，对四个子系统输出的卫星位置、速度参数进行数据融合 处理，得到卫星的位置和速度参数最优估计值，根据所述通过信息融合得到的卫星位置和速 度参数对卫星进行导航。对子系统故障的情况下，设计故障检测与隔离方法，隔离故障子系 统，排除故障的影响。
其中四个卫星自主导航子系统分别是基于X射线探测器的卫星自主导航子系统，基于红 外地平仪和星敏感器的卫星自主导航子系统，基于雷达高度计的卫星自主导航子系统，基于 紫外三轴敏感器的卫星自主导航子系统，分别称为子系统1，子系统2，子系统3，子系统4。
具体包括以下步骤：
(1)建立基于轨道动力学的卫星运动状态方程；
(2)建立基于X射线探测器的卫星自主导航子系统；
(3)建立基于红外地平仪和星敏感器的卫星自主导航子系统；
(4)建立基于雷达高度计的卫星自主导航子系统；
(5)建立基于紫外三轴敏感器的卫星自主导航子系统；
(6)对上述四个子系统的输出进行信息融合处理，输出导航信息；
(7)采用χ2检验法检测发生故障的系统并隔离故障子系统。
所述的步骤(1)中，建立包含了J2项地球形状摄动和日月摄动的卫星轨道动力学方程，在 J2000.0地心惯性系下，卫星轨道动力学方程为
                    X＝f(X，t)+W(t)                                           (1)
式中，X＝[x，y，z，vx，vy，vz]T为状态变量，x，y，z，vx，vy，vz分别为卫星在X，Y，Z三个方向的位置和速度值； f(X，t)的形式根据卫星所设摄动的不同，仅考虑J2项的地球形状摄动和日、月的三体引力情 况下，
$f(X,t)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mu }_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mu }_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mu }_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mu }_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mu }_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mu }_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mu }_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mu }_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mu }_m+(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$
其中，μe，μs，μm分别为地球、太阳、月亮的引力常数，Re为地球半径，J2为地球引力系数， rs＝[xs，ys，zs]T和rm＝[xm，ym，zm]T分别表示太阳和月亮在地心惯性系下的位置，xs，ys，zs为太阳在地心 惯性系下的坐标值，xm，ym，zm为月亮在地心惯性系下的坐标值，x，y，z，vx，vy，vz分别为卫星在X，Y，Z 三个方向的位置和速度值，r为地心到卫星的距离，W(t)为模型误差向量，代表了地球非球形 摄动的高阶项、太阳光压摄动大气阻力摄动等摄动力的影响，假设为零均值高斯白噪声；
所述步骤(2)，包括X射线脉冲星的品质因素计算方法，卫星位置与X射线脉冲星自转相 位之间关系的计算方法，整周模糊度计算方法及GDOP值计算方法；
所述步骤(2)中的X射线脉冲星的品质因数的计算方法为
$Q=754.8\mathrm{SNR}\sqrt{T}\frac{1}{{5T}_{10}+{7T}_{50}}\sqrt{\frac{{5T}_{10}-7T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$
式中，SNR为脉冲星的信噪比，T为脉冲星旋转周期，T50和T10分别为脉冲波峰强度50％ 和10％时的脉冲宽度；
所述步骤(2)中的卫星位置与X射线脉冲星自转相位之间关系的计算方法为
$[\mathrm{\Delta \phi }+\mathrm{\Delta N}]·\lambda -\hat{n}·r_E=\hat{n}·r+\mathrm{c\delta }{t}_{\mathrm{sat}}+{\Delta }_{\mathrm{rel}}+v_k---\left(4\right)$
其中，Δφ＝φSSB-φsat，ΔN为卫星与太阳系质心SSB之间的整周模糊度值，λ为脉冲星的脉 冲辐射的波长，为脉冲星的视线方向，rE为地球相对与SSB的位置矢量，r为卫星相对于地 心的位置矢量，c为光速，δtsat为卫星上的钟差，vk为脉冲相位测量噪声，Δrel为相对论效应改 正项，主要包括Roemer延迟改正、Shapiro延迟改正、太阳系行星总延迟改正；
所述步骤(2)中的整周模糊度的计算方法：
$\mathrm{\Delta N}=\mathrm{round}(\hat{n}·(r_E+\tilde{r})/\lambda -\mathrm{\Delta \phi })---\left(5\right)$
其中，为卫星动力学方程提供的卫星位置预报值；
所述步骤(2)中的GDOP值计算方法：
$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left(C\right)}---\left(6\right)$
其中，trace表示求矩阵的迹，C为位置误差协方差矩阵，由下式算得
$C=\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}\left[\begin{array}{ccc}{\lambda }_1^2{\sigma }_{\phi 1}^2& 0& 0\\ 0& {\lambda }_2^2{\sigma }_{\phi 2}^2& 0\\ 0& 0& {\lambda }_3^2{\sigma }_{\phi 3}^2\end{array}\right]{\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}}^T---\left(7\right)$
其中，H为量测方程的系数矩阵，λi(i＝1，2，3)为第i颗脉冲星的脉冲辐射波长，σφi为第i颗 脉冲星的相位测量精度；
所述步骤(3)中，以星光角距和地心距为观测量的量测方程为
$Z_1=\arccos (-\frac{r·s}{r})+v_1---\left(8\right)$
$Z_2=\sqrt{x^2+y^2+z_2}+v_2---\left(9\right)$
其中，Z1为星光角距观测量，Z2为地心距观测量，r为卫星本体系下的地心矢量，s为卫 星本体系下星光方向的单位矢量，v1为星光角距测量噪声，v2为地心距测量噪声，(x，y，z)为在 地心惯性系下表示的卫星位置坐标值；
所述步骤(4)中，雷达高度计垂直安装于卫星的正下面，测量卫星平台到地表的距离，且 大地水准面采用克莱饶椭球模型；
所述步骤(4)中，卫星星下点的地形高度h(x，y)由随机软件生成的随机地形提供，随机地 形生成方法如下：
地形在高度方向上分解为地形基准平面高度h0和在此平面上叠加的地形起伏Z(x，y)，即 (x，y)点的实际海拔高度h(x，y)表示成：
                h(x，y)＝h0+Z(x，y)                                 (10)
人工地形由二维一阶离散自回归过程产生：
          Z(xi，yi)＝a1Z(xi-Δx，yi)+a2Z(xi，yi-Δy)+a3Z(xi-Δx，yi-Δy)+W(xi，yi)    (11)
其中，Z(xi，yi)为(xi，yi)点的地形起伏高度，Δx，Δy分别是沿x，y方向的采样间隔，W(xi，yi) 为零均值白噪声序列，且W(xi，yi)～N(0，σw 2)；
根据地形的各向同性和平稳特性，令系数a1＝a2＝a，a3＝b，且当a＝exp(-1/Tauc)，b＝-a2， ${\sigma }_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$时，式(11)中两维随机函数Z(xi，yi)的相关函数为
$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(12\right)$
由此得到的两维随机过程，其均方差为dRMS，按指数衰减，抗相关系数为Tauc，模拟地 形在x方向和y方向的边界地形由两个一维随机过程分别产生，在x方向的地形边界上， yi＝y0，则有
           Z(xi，y0)＝aZ(xi-Δx，y0)+W(xi，y0)                     (13)
在y方向的地形边界上，xi＝x0，则有
        Z(x0，yi)＝aZ(x0，yi-Δy)+W(x0，yi)                 (14)
通过设定h0、dRMS、Tauc三个参数，利用软件编程产生不同的随机地形，将产生的随机地 形h(x，y)叠加到大地水准面上；
所述步骤(5)中，以地心方向和地心距观测量的离散形式量测方程为
$z=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(15\right)$
其中，uk为地心方向的单位矢量，rk为地心距，vk为量测噪声；
所述步骤(2)(3)(4)(5)中，四个子系统的子滤波器，采用扩展卡尔曼滤波方法或者 Unscented卡尔曼滤波方法，分别独立进行时间更新和测量更新；
所述步骤(6)中，对上述步骤(2)(3)(4)(5)四个子系统的输出进行信息融合处理，采用的 方法如下：
第i(i＝1，2，3，4)个子系统的量测方程为：
           Zi(k)＝Hi(k)X(k)+Vi(k)，i＝1，2，3，4                (16)
初始时刻全局状态估计为其协方差阵为Pg0，根据信息守恒原理，将这一信息通过信 息分配因子βi分配到四个子滤波器和主滤波器，分配原则如下：
$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\beta }_i{P}_{g0}^{-1}\left(k\right)\\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)\\ Q_{i0}^{-1}\left(k\right)={\beta }_i{Q}_{g0}^{-1}\left(k\right)\end{array}\right.,i=\mathrm{1,2,3,4},m---\left(17\right)$
其中，i＝1，2，3，4表示四个子滤波器，m表示主滤波器，信息分配因子βi满足信息守恒原理： $\underset{i=1}{\overset{4,m}{\Sigma }}{\beta }_i=1,$此处取βm＝0，β1＝β2＝β3＝β4＝1/4，各子滤波器分别独立进行时间更新，并利用其对应传 感器的量测信息进行量测更新，得四个子滤波器的局部估计值和估计误差协方差阵 Pi(k)i＝1，2，3，4；在主滤波器中，将子滤波器输出的局部信息按下式进行信息融合，得到全局状 态估计信息为
$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\Sigma }}{P}_i^{-1}{\hat{X}}_i\\ P_g^{-1}=P_1^{-1}+P_2^{-1}+P_3^{-1}+P_4^{-1}\end{array}\right.---\left(18\right)$
即按照上式进行信息融合，输出全局状态估计和其协方差矩阵Pg，将输出的全局状态 估计作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模型，作 为卫星轨道动力学预报下一时刻卫星状态的初值；
所述步骤(7)中，χ2检验法进行系统故障检测与隔离的方法如下：
在每个子滤波器中均加入故障检测和隔离的χ2方法，构成容错联邦滤波器，利用四个子 系统的故障检测模块，分别对基于X射线探测器的卫星自主导航子系统、基于红外地平仪和 星敏感器的卫星自主导航子系统、基于雷达高度计的卫星自主导航子系统、基于紫外三轴敏 感器的卫星自主导航子系统故障与否进行检测；
各子滤波器量测更新之后，进行故障检测和隔离，其方法是首先设计四个子系统的故障 检测函数Di(k)i＝1，2，3，4，
$d_i\left(k\right)=z_i\left(k\right)-H_i\left(k\right){\hat{X}}_i(k/k-1)---\left(19\right)$
$S_i\left(k\right)=H_i\left(k\right)P_i(k/k-1)H_i^T\left(k\right)+R_i\left(k\right)---\left(20\right)$
$D_i\left(k\right)=d_i^T\left(k\right)S_i^{-1}\left(k\right)d_i\left(k\right)---\left(21\right)$
然后进行故障判断：D(k)＞TD时，有故障，D(k)＜TD时，无故障；TD为预先设定的门限， 可由误警率确定，当误警率Pfa＝α时，由Pfa＝P[λk＞TD|H0]＝α解出门限TD。
本发明与现有技术相比的优点在于：
(1)本发明将多种传感器的量测信息融合得到精度更高的导航信息，具有为卫星提供 高精度自主导航信息的优点；同时本发明采用多传感器冗余配置，并采用系统级的故障检测 与隔离方法，提高了系统的可靠性和容错能力。
(2)本发明中的基于X射线探测器的卫星自主导航子系统，具有提供十米量级的自主 定轨精度能力，且适用于深空探测器的高精度自主导航需求；基于雷达高度计的卫星自主导 航子系统，采用雷达高度计垂直安装的方式，使用合理的海平面模型和模拟地形，提高了卫 星的导航精度。
(3)本发明中的每个自主导航子系统，均可独立构成卫星自主导航系统，因而可在不 同的情况下选用不同的子系统进行组合，构建新的组合导航系统，也可以只使用其中的一种 子系统，作为卫星的自主导航系统。适合在不同导航精度要求下，对导航系统成本以及重量、 体积的控制。
(4)本发明适用于月球探测器和火星探测器飞行各个阶段的高精度自主导航。在月球 探测器和火星探测器飞行的各个阶段，星敏感器、紫外三轴敏感器、X射线探测器都能够为 导航系统提供足够的导航观测量，从而为航天器提供高精度导航信息。在绕月飞行和绕火星 飞行阶段，雷达高度计测距信息提供航天器与月球或火星表面之间的距离，从而确定绕月或 绕火飞行航天器的轨道高度，为航天器提供高精度导航信息，实现深空探测飞行器完全自主 式的导航。
附图说明
图1为本发明的组合导航系统容错联邦滤波器结构图；
图2为本发明中基于X射线探测器的卫星自主导航子系统的脉冲星导航示意图。

基于多信息融合的卫星自主导航系统由四个子系统和一个主滤波器构成，如图1所示， 本发明的具体实施方法如下：
一、建立基于轨道动力学的卫星运动状态方程
(1)建立含日月摄动的卫星轨道动力学方程
在J2000.0地心惯性系下，卫星轨道动力学方程(系统状态方程)为
                    X＝f(X，t)+W(t)                                             (1)
式中，X＝[x，y，z，vx，vy，vz]T为状态变量，x，y，z，vx，vy，vz分别为卫星在X，Y，Z三个方向的位置和速度值。 f(X，t)的形式根据卫星所设摄动的不同，在仅考虑J2项的地球形状摄动和日、月的三体引力 情况下，
$f(X,t)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mu }_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mu }_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mu }_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mu }_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mu }_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mu }_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mu }_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mu }_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mu }_m+(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$
其中，μe，μs，μm分别为地球、太阳、月亮的引力常数，Re为地球半径，J2为地球引力系数， rs＝[xs，ys，zs]T和rm＝[xm，ym，zm]T分别表示太阳和月亮在地心惯性系下的位置，W(t)为模型误差向量， 代表了地球非球形摄动的高阶项、太阳光压摄动大气阻力摄动等摄动力的影响，假设为零均 值高斯白噪声。
二、建立基于X射线探测器的卫星自主导航子系统
(2)导航用X射线脉冲星的选择方法
导航用X射线脉冲星的选取原则如下：毫秒脉冲星为主要候选脉冲星；用于导航的脉冲 星自转周期稳定，具有自转跃变的脉冲星应排除在外；并且考虑脉冲星辐射信号的强弱、脉 冲轮廓形状、脉冲星所在方向的背景辐射，并选择基本参数测定精度高的脉冲星；所有脉冲 星尽量在空间均匀分布。脉冲星品质因数计算公式为：
$Q=754.8\mathrm{SNR}\sqrt{T}\frac{1}{{5T}_{10}+{7T}_{50}}\sqrt{\frac{{5T}_{10}-7T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$
式中，SNR为脉冲星的信噪比，T为脉冲星旋转周期，T50和T10分别为脉冲波峰强度50％和 10％时的脉冲宽度。计算候选脉冲星的品质因数值，Q值越大，脉冲星的导航属性越好，越可 以用作导航脉冲星。
(3)建立以脉冲星自转相位(TOA)为观测量的量测方程
图2为绕地卫星的脉冲星导航示意图。其中，Sat、E、SSB分别表示卫星、地球和太阳 系质心，为某颗脉冲星的视线方向，T为一个脉冲周期；rE和rSSB_sat分别表示地球、卫星相 对于SSB的位置矢量，r表示卫星相对于地心的位置矢量。观测历元t时卫星处测得的脉冲相 位值为φsat，太阳系质心的相位值为φSSB。
定义在太阳系质心的脉冲星钟模型为
$\Phi \left(t\right)=\phi \left(t_0\right)+f(t-t_0)+\underset{m=1}{\overset{M}{\Sigma }}\frac{f^{\left(m\right)}{(t-t_0)}^{m+1}}{(m+1)!}---\left(4\right)$
式中，φ(t0)为参考历元t0时的脉冲星相位，f是脉冲星在t0时的自转频率，f(m)为f的m阶 导数(一般取m＝1，2，3)。由(4)式预报可以得到t时太阳系质心处的脉冲相位值ΦSSB，取其小 数部分φSSB(0≤φSSB＜1)。
卫星与SSB之间的整周模糊度为ΔN，脉冲星的脉冲辐射的波长为λ(λ＝c·T，c为光速)， 则测量方程为
$[\mathrm{\Delta \phi }+\mathrm{\Delta N}]·\lambda -\hat{n}·r_E=\hat{n}·r+\mathrm{c\delta }{t}_{\mathrm{sat}}+{\Delta }_{\mathrm{rel}}+v_k---\left(5\right)$
其中，Δφ＝φSSB-φsat，vk为脉冲相位测量噪声。Δrel为相对论效应改正项，主要包括Roemer 延迟改正、Shapiro延迟改正、太阳系行星总延迟改正。地球矢量rE由JPL星历表提供。使 用三颗脉冲星的观测量，则量测方程为
$\left[\begin{array}{c}({\mathrm{\Delta \phi }}_1+{\mathrm{\Delta N}}_1)·{\lambda }_1-{\hat{n}}_1·r_E\\ ({\mathrm{\Delta \phi }}_2+{\mathrm{\Delta N}}_2)·{\lambda }_2-{\hat{n}}_2·r_E\\ ({\mathrm{\Delta \phi }}_3+\Delta N_3)·{\lambda }_3-{\hat{n}}_3·r_E\end{array}\right]=H·\left[\begin{array}{c}r\\ \mathrm{c\delta }{t}_{\mathrm{sat}}\end{array}\right]+V_k---\left(6\right)$
其中$H=\left[\begin{array}{cc}{\hat{n}}_1& 1\\ {\hat{n}}_2& 1\\ {\hat{n}}_3& 1\end{array}\right].$
(4)计算整周模糊度值
使用卫星轨道预报器提供的卫星位置预报值整周模糊度为ΔN为
$\mathrm{\Delta N}=\mathrm{round}(\hat{n}·(r_E+\tilde{r})/\lambda -\mathrm{\Delta \phi })---\left(7\right)$
(5)计算GDOP值
当有多颗X射线脉冲星可见时，选择几何精度因子(GDOP)最小的三颗脉冲星，作为导航 解算脉冲星。对于同时观测三颗脉冲星的情况，其GDOP值计算方法如下：
位置误差协方差矩阵为
$\mathrm{cov}\left(P\right)=\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}\left[\begin{array}{ccc}{\lambda }_1^2{\sigma }_{\phi 1}^2& 0& 0\\ 0& {\lambda }_2^2{\sigma }_{\phi 2}^2& 0\\ 0& 0& {\lambda }_3^2{\sigma }_{\phi 3}^2\end{array}\right]{\left\{{\left(H^{T}H\right)}^{-1}{H}^T\right\}}^T---\left(8\right)$
其中σφi(i＝1，2，3)为第i颗脉冲星的相位测量精度。则GDOP值为
$\mathrm{GDOP}=\sqrt{\mathrm{trace}\left(C\right)}=\sqrt{{\sigma }_x^2+{\sigma }_y^2+{\sigma }_z^2+{\sigma }_t^2}---\left(9\right)$
其中，trace表示求矩阵的迹，σx 2，σy 2，σz 2，σt 2为矩阵C的对角线元素，分别表示x，y，z轴的距离 以及时间变量t的测量精度。
(6)子滤波器的滤波方法
子滤波器采用非线性卡尔曼滤波方法，可以考虑扩展卡尔曼滤波方法(EKF)或者 Unscented卡尔曼滤波分别独立进行时间更新和测量更新。
扩展卡尔曼滤波处理的过程如下：
对状态方程和量测方程进行离散化，得到状态方程的离散型线性干扰方程以及离散型量 测方程为
              δXk+1＝φk+1，kδXk+Wk                       (10)
              Zk+1＝Hk+1Xk+1+Vk+1                           (11)
式中φk+1，k为状态转移矩阵，Zk+1为在三颗脉冲星视线方向上卫星与太阳系质心之间的距离 观测量。状态模型噪声的协方差阵为$E\left[W_k{W}_k^T\right]=Q\left(k\right),$量测噪声的协方差阵为$E\left[V_{k+1}{V}_{k+1}^T\right]=R\left(k\right),$Wk 与Vk+1互不相关。
利用状态方程(10)和量测方程(11)组成如图1所示的子滤波器1，采用扩展卡尔曼 滤波方法进行信息融合处理，输出子滤波器得到的卫星位置、速度信息。扩展卡尔曼滤波基 本方程如下：
$\left\{\begin{array}{c}{\hat{X}}_{k+1,k}={\hat{X}}_{k,k}+f[{\hat{X}}_{k,k},t_k]·T\\ {\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-H_{k+1}{X}_{k+1/k})\\ P_{k+1/k}={\phi }_{k+1/k}{P}_{k/k}{\phi }_{k+1/k}^T+Q_k\\ K_{k+1}=P_{k+1/k}{H}_{k+1}^T[H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R_{k+1}{]}^{-1}\\ P_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T\end{array}\right.---\left(12\right)$
式中，f表示卫星轨道动力学方程，T为卡尔曼滤波周期，φk+1，k为状态转移矩阵，Kk+1为 卡尔曼增益系数，Hk+1为量测方程系数矩阵，Pk+1，k为最优预测估值误差协方差阵，Pk+1，k+1为最 优滤波值误差协方差阵，为一步预测值，为卡尔曼滤波值。
Unscented卡尔曼滤波处理的过程如下：
对于系统方程和量测方程如下的离散系统
$\left\{\begin{array}{c}{X}_{k+1}=F(X_k,k)+W_k\\ Z_{k+1}=H\left(X_k\right)+V_k\end{array}\right.---\left(13\right)$
设系统状态变量的维数为n×1维，那么2n+1个采样点为
${\chi }_0=\hat{X}$
${\chi }_i=\hat{X}+{\left(\sqrt{(n+\lambda )P_{k/k}}\right)}_i,i=\mathrm{1,2},...,n---\left(14\right)$
${\chi }_i=\hat{X}-{\left(\sqrt{(n+\lambda )P_{k/k}}\right)}_i,i=n+1,...,2n$
每个采样点对应的权重为
$W_0^{\left(m\right)}=\lambda /(n+\lambda )$
$W_0^{\left(c\right)}=\lambda /(n+\lambda )+(1-{\alpha }^2+\beta )---\left(15\right)$
$W_i^{\left(m\right)}=W_i^{\left(c\right)}=\lambda /\left\{2(n+\lambda )\right\}i=\mathrm{1,2},...,2n$
其中，λ＝α2(n+k)-n是一个标量参数，常数α确定采样点围绕的分布特征，通常被 设置为一个小的正数(1≥α≥10-4)，常数k是一个标量参数，通常被设置为2或3-n。Wi (m)为 均值加权所用权值，Wi (c)为协方差加权所用权值。Unscented卡尔曼滤波的基本设施过程如 下：
计算采样点
${\chi }_k=\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{x}}_k+\sqrt{n+\lambda }\sqrt{P_{k/k}}& {\hat{X}}_k-\sqrt{n+\lambda }\sqrt{P_{k/k}}\end{array}\right]---\left(16\right)$
时间更新
$\left\{\begin{array}{c}{\chi }_{k+1/k}^{*}=F({\chi }_k,k),{\hat{X}}_{k+1/k}=\underset{0}{\overset{2n}{\Sigma }}{W}_i^{\left(m\right)}{\chi }_{k+1/k}^{*}\\ P_{k+1/k}=\underset{i=0}{\overset{2n}{\Sigma }}{W}_i^{\left(c\right)}[{\chi }_{i,k+1/k}^{*}-{\hat{X}}_{k+1/k}]{[{\chi }_{i,k+1/k}^{*}-{\hat{X}}_{k+1/k}]}^T+Q_{k+1}\\ {\chi }_{k+1/k}=\left[\begin{array}{ccc}{\hat{\chi }}_{k+1/k}& {\hat{X}}_{k+1/k}+\gamma \sqrt{P_{k+1/k}}& {\hat{X}}_{k+1/k}-\gamma \sqrt{P_{k+1/k}}\end{array}\right]\\ z_{k+1/k}=H\left({\chi }_{k+1/k}\right),{\hat{Z}}_{k+1/k}=\underset{0}{\overset{2n}{\Sigma }}{W}_i^{\left(m\right)}{z}_{k+1/k}\end{array}\right.---\left(17\right)$
量测更新
$\left\{\begin{array}{c}{P}_{{\hat{Z}}_{k+1},{\hat{Z}}_{k+1}}=\underset{i=0}{\overset{2n}{\Sigma }}{W}_i^{\left(c\right)}[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]{[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]}^T+R_{k+1}\\ R_{{\hat{X}}_{k+1},{\hat{Z}}_{k+1}}=\underset{i=0}{\overset{2n}{\Sigma }}{W}_i^{\left(c\right)}[{\chi }_{i,k+1/k}-{\hat{X}}_{k+1/k}]{[z_{i,k+1/k}-{\hat{Z}}_{k+1/k}]}^T+R_{k+1}\\ K_{k+1}=P_{{\hat{X}}_{k+1},{\hat{Z}}_{k+1}}{P}_{{\hat{Z}}_{k+1},{\hat{Z}}_{k+1}}^{-1},{\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-{\hat{Z}}_{k+1/k})\\ P_{k+1/k+1}=P_{k+1/k}-K_{k+1}{P}_{{\hat{Z}}_{k+1},{\hat{Z}}_{k+1}}{K}_{k+1}^T\end{array}\right.---\left(18\right)$
其中，状态模型噪声的协方差阵为$E\left[W_k{W}_k^T\right]=Q_k,$量测噪声的协方差阵为$E\left[V_{k+1}{V}_{k+1}^T\right]=R_k,$Wk与Vk+1 互不相关。
利用状态方程(1)和量测方程(6)组成子滤波器1，对方程(1)和(6)分别进行离 散化，采用扩展卡尔曼滤波方法或者Unscented卡尔曼滤波方法进行信息融合处理，输出子 滤波器得到的卫星位置、速度信息。
三、建立基于红外地平仪和星敏感器的卫星自主导航子系统
(7)建立以星光角距和地心距为观测量的量测方程
星敏感器观测导航恒星，可以确定星光在卫星本体坐标系中的方向；利用红外地平仪测 得地心矢量在卫星本体坐标系中的方向；这样得到星光矢量和地心矢量之间的夹角，即星光 角距。红外地平仪观测地球红外图像，获得地球半张角ρ，地球半径Re为已知，从而可解算 出地心距$r=\frac{R_e}{\sin \rho },$从而以星光角距和地心距为观测量的量测方程为
$Z_1=\arccos (-\frac{r·s}{r})+v_1---\left(19\right)$
$Z_2=\sqrt{x^2+y^2+z^2}+v_2---\left(20\right)$
其中，Z1为星光角距观测量，Z2为地心距观测量，r为卫星本体系下的地心矢量，s为卫 星本体系下星光方向的单位矢量，v1为星光角距测量噪声，v2为地心距测量噪声，(x，y，z)为在 地心惯性系下表示的卫星位置坐标值。
(8)量测方程(19)、(20)构成的量测方程组，利用状态方程(2)和量测方程组构建如 图1所示的子滤波器2，重复上述步骤(6)的扩展卡尔曼滤波方法或者Unscented卡尔曼滤波 处理方法进行数据处理，输出卫星的位置、速度估计值。
四、建立基于雷达高度计的卫星自主导航子系统
(9)建立以星光角距和轨道高度为观测量的量测方程
雷达高度计垂直按照在卫星的正下方，测得卫星平台到地表的距离，即卫星的轨道高度 量测信息。大地水准面采用克莱饶椭球模型，并将海平面的起伏作为修正项加以考虑。克莱 饶水准面可以表示为下式
饶水准面可以表示为下式

其中，e为地心纬度，P2为二次勒得让多项式，J2为地球带谐项系数，R为地球平均半 径值，参数αe为
${\alpha }_e=\frac{\left[\frac{2}{3}{J}_2(\frac{{\alpha }_e}{R}+\frac{m}{2})\right]}{[1-\frac{2}{3}m]m}\approx 0.003449963---\left(22\right)$
则以轨道高度为观测量的量测方程为
$Z_2=r-R[1-\frac{{2\alpha }_e}{3}\left(0.5\right){\left(\frac{z}{r}\right)}^2-\frac{1}{2}]-h(x,y)+v_2---\left(23\right)$
式中，Z2为雷达高度计测得的卫星轨道高度，r为卫星到地心的距离，v2为测量噪声，h(x，y) 为星下点的地形高度，由随机地形软件生成。以星光角距为观测量的量测方程同方程(19)。
(10)随机地形生成方法
建立模拟地形，模拟实际地形的变化。地形在高度方向上分解为地形基准平面高度h0和 在此平面上叠加的地形起伏Z(x，y)，即(x，y)点的实际海拔高度h(x，y)可表示成：
                  h(x，y)＝h0+Z(x，y)                                   (24)
人工地形由二维一阶离散自回归过程产生：
          Z(xi，yi)＝a1Z(xi-Δx，yi)+a2Z(xi，yi-Δy)+a3Z(xi-Δx，yi-Δy)+W(xi，yi)    (25)
其中，Z(xi，yi)为(xi，yi)点的地形起伏高度，Δx，Δy分别是沿x，y方向的采样间隔，W(xi，yi) 为零均值白噪声序列，W(xi，yi)～N(0，σw 2)。
根据地形的各向同性和平稳特性，可令系数a1＝a2＝a，a3＝b，且当a＝exp(-1/Tauc)，b＝-a2， ${\sigma }_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$时，式(25)中两维随机函数Z(xi，yi)的相关函数为
$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(26\right)$
由此得到的两维随机过程，其均方差为dRMS，按指数衰减，抗相关系数为Tauc。模拟地 形在x方向和y方向的边界地形由两个一维随机过程分别产生。若在x方向的地形边界上， yi＝y0，则有
              Z(xi，y0)＝aZ(xi-Δx，y0)+W(xi，y0)                         (27)
若在y方向的地形边界上，xi＝x0，则有
              Z(x0，yi)＝aZ(x0，yi-Δy)+W(x0，yi)                (28)
通过设定h0、dRMS、Tauc三个参数，利用软件编程可以产生不同的随机地形。将产生的随 机地形h(x，y)叠加到大地水准面上。
(11)量测方程(19)、(23)构成的量测方程组，利用状态方程(2)和量测方程组构建如 图1所示的子滤波器3，重复上述步骤(6)的扩展卡尔曼滤波方法或者Unscented卡尔曼滤波 处理方法进行数据融合，输出卫星的位置、速度估计值。 五、建立基于紫外三轴敏感器的卫星自主导航子系统
(12)建立基于紫外三轴敏感器的卫星自主导航子系统
紫外三轴敏感器敏感地球紫外图像，通过图像处理提取星体系中的地心方向和地心距信 息；利用卫星的姿态矩阵转换得到卫星在地心惯性系中的位置矢量，从而得到地心方向和地 心距观测量，其离散形式的量测方程为
$z=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(29\right)$
其中，uk为地心方向的单位矢量，rk为地心距，vk为量测噪声。
利用状态方程(2)和量测方程(29)构建如图1所示的子滤波器4，重复上述步骤(6) 的扩展卡尔曼滤波方法或者Unscented卡尔曼滤波处理方法进行数据融合，输出卫星的位置、 速度估计值。
六、主滤波器的信息融合、故障检测与隔离方法
(13)主滤波器的信息融合处理
第i(i＝1，2，3，4)个子系统的量测方程为：
          Zi(k)＝Hi(k)X(k)+Vi(k)，(i＝1，2，3，4)              (30)
初始时刻全局状态估计为其协方差阵为Pg0。根据信息守恒原理，将这一信息通过 信息分配因子βi分配到四个子滤波器和主滤波器，分配原则如下：
$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\beta }_i{P}_{g0}^{-1}\left(k\right)\\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)\\ Q_{i0}^{-1}\left(k\right)={\beta }_i{Q}_{g0}^{-1}\left(k\right)\end{array}\right.(i=\mathrm{1,2,3,1},m)---\left(31\right)$
其中，i＝1，2，3，4表示四个子滤波器，m表示主滤波器，信息分配因子βi满足信息守恒原理： $\underset{i=1}{\overset{4,m}{\Sigma }}{\beta }_i=1,$此处取βm＝0，β1＝β2＝β3＝β4＝1/4。各子滤波器分别独立进行时间更新，并利用其对应传 感器的量测信息进行量测更新，得四个子滤波器的局部估计值和估计误差协方差阵 Pi(k)(i＝1，2，3，4)，具体措施如方程(12)或(17)和(18)。在主滤波器中，将子滤波器输出的 局部信息按下式进行信息融合，得到全局状态估计信息为
$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\Sigma }}{P}_i^{-1}{\hat{X}}_i\\ P_g^{-1}=P_1^{-1}+P_2^{-1}+P_3^{-1}+P_4^{-1}\end{array}\right.---\left(32\right)$
即按照(32)式进行信息融合，输出全局状态估计和其协方差矩阵Pg。将输出的全局 状态估计作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模 型，作为卫星轨道动力学预报下一时刻卫星状态的初值。
联邦滤波中，主滤波器和四个子滤波器联合工作，信息分配只在初始时刻进行，子滤波 器独立工作，不接收主滤波器的反馈信息。
(14)组合导航的故障检测与隔离方法
采用系统级的故障检测与隔离方法，排除发生故障的传感器及相应子系统。在每个子滤 波器中均加入故障检测和隔离的χ2方法，构成容错联邦滤波器，其结构如图1所示。图1中 FDI 1、FDI 2、FDI 3、FDI 4分别为四个子系统的故障检测模块，分别对子系统1、2、3、4 故障与否进行检测，并决定是否隔离，以便主滤波器对正常的子系统进行重组，使部分导航 传感器故障导致子系统失效后，综合导航系统仍能继续工作，具容错性。
各子滤波器分别独立进行时间更新和量测更新。各子滤波器量测更新之后，进行故障检 测和隔离。其方法是首先设计四个子系统的故障检测函数Di(k)(i＝1，2，3，4)：
$d_i\left(k\right)=z_i\left(k\right)-H_i\left(k\right){\hat{X}}_i(k/k-1)---\left(33\right)$
$S_i\left(k\right)=H_i\left(k\right)P_i(k/k-1)H_i^T\left(k\right)+R_i\left(k\right)---\left(34\right)$
$D_i\left(k\right)=d_i^T\left(k\right)S_i^{-1}\left(k\right)d_i\left(k\right)---\left(35\right)$
然后进行故障判断：D(k)＞TD时，有故障；D(k)＜TD时，无故障。TD为预先设定的门限， 可由误警率确定。当误警率Pfa＝α时，Pfa＝P[λk＞TD|H0]＝α，解出门限TD。若判断出子系统无故 障，则把其滤波结果送到主滤波器；若子系统故障，则滤波结果不送到主滤波器，并对其隔 离，直至该子系统恢复正常。
在各子滤波器计算出各自的局部估计，并进行故障检测后，由主滤波器把输出正常的子 滤波器的公共状态进行信息融合，得到全局估计信息。
(15)输出位置、速度信息。
按照上述步骤(1)-(15)进行计算机仿真，建立轨道动力学方程和量测方程，利用扩展卡 尔曼滤波方法或者Unscented卡尔曼滤波方法对四个子系统的观测量进行滤波处理，再采用 信息融合手段得到状态量的全局状态估计${\hat{X}}_g={\left[\begin{array}{cccccc}\hat{x}& \hat{y}& \hat{z}& {\hat{v}}_x& {\hat{v}}_y& {\hat{v}}_z\end{array}\right]}^T$和其协方差矩阵 $P_g={\left[\begin{array}{cccccc}{p}_x& p_y& p_z& p_{v_s}& p_{v_y}& p_{v_z}\end{array}\right]}^T,$其中 分别是在地心惯性系下对卫星在X、Y、Z 三个方向的位置、速度量x，y，z，vx，vy，vz的最优估计。
本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。