CN101178312A - 基于多信息融合的航天器组合导航方法 - Google Patents

Abstract

一种基于多信息融合的航天器组合导航方法，属于航天器自主导航方法。该导航方法包括以下步骤：建立基于轨道动力学的卫星运动状态方程，建立基于X射线探测器的卫星自主导航子***，建立基于红外地平仪和星敏感器的卫星自主导航子***，建立基于雷达高度计的卫星自主导航子***，建立基于紫外三轴敏感器的卫星自主导航子***，对上述四个子***的输出进行信息融合处理，输出导航信息，采用x²检验法检测发生故障的***并隔离故障子***。本组合导航方法，能实现深空探测器的高精度自主导航，可靠性高，容错能力强。

Description

基于多信息融合的航天器组合导航方法

技术领域

本发明属于航天器组合导航技术领域，适用于绕地球飞行航天器的高精度自主导航，也可以用于深空探测器，以及行星际飞行航天器的高精度自主导航。

背景技术

航天器自主导航***可以在少依赖甚至最终不依赖于地面***支持的情况下，在轨实时确定航天器的位置和速度，实现自主运行。对于卫星***来说，自主导航有利于降低卫星对地面的依赖程度，提高***生存能力，例如战时，当地面***遭到敌方的破坏和干扰时，仍能完成轨道的确定和保持，这对军事卫星来讲，具有非常重要的意义。另外，自主导航还可以有效地减轻地面测控站的负担，降低地面支持成本，从而降低整个航天计划的研制费用。

目前航天器的自主导航方法主要包括GPS，星间链路，磁力计，雷达和天文导航方法。其中GPS和星间链路虽然可进行实时导航且达到较高的精度，但是航天器必须依赖其他航天器提供量测信息，不能算是完全自主的导航方法。使用雷达高度计进行卫星的自主导航，受到地形高度变化的影响，但是在海洋上方时能获得较好的导航性能。

天文导航是利用天体敏感器测得的天体(如行星、恒星、X射线脉冲星等)信息进行航天器位置解算的一种完全自主的定位导航方法，可同时提供姿态、位置、速度以及时间信息；适用于低轨卫星、高轨卫星和深空探测器，因而备受青睐。目前卫星的天文导航中，根据被观测对象属性的不同，使用的敏感器主要有星敏感器、红外地平仪、紫外三轴敏感器以及X射线探测器等。其中X射线探测器的观测对象为X射线源的脉冲辐射，获取脉冲到达时间信息；使用X射线探测器能够为卫星提供精度较高的导航信息。对于多种天体测量敏感器的组合，由于其被观测对象属性的差异，因而其导航属性存在互补的可能。

随着电子技术、导航技术和控制理论的发展，卫星上拥有众多的导航传感器模块，具备了构成多传感器组合导航信息***的可能条件。随着航空航天事业的发展，特别是现代化战争对导航定位的精度和可靠性提出的更高要求，单一导航***已难以满足。把两种或多种导航***组合起来，应用最优估计理论，形成最优组合导航***，有利于充分运用各导航***的信息进行信息互补和信息合作，使组合后的导航***在精度和可靠性方面都有所提高。同时这也对卫星组合导航***发生故障时的自动检测和故障隔离能力，提出了更高的要求。

发明内容

本发明的目的在于：利用多个传感器进行组合，将各个传感器的测量数据进行信息融合处理，提供了一种采用信息融合手段为卫星提供高精度自主导航参数的方法，以及提高卫星在传感器发生故障时的自动检测能力和故障隔离能力，重构导航***的方法。

本发明所采用的技术方案是：基于多信息融合的卫星组合导航***，利用X射线探测器、星敏感器、红外地平仪、雷达高度计和紫外三轴敏感器五种传感器以及卫星的轨道动力学模型构成了四个卫星自主导航子***，对四个子***输出的卫星位置、速度参数进行数据融合处理，得到卫星的位置和速度参数最优估计值，根据所述通过信息融合得到的卫星位置和速度参数对卫星进行导航。对子***故障的情况下，设计故障检测与隔离方法，隔离故障子***，排除故障的影响。

其中四个卫星自主导航子***分别是基于X射线探测器的卫星自主导航子***，基于红外地平仪和星敏感器的卫星自主导航子***，基于雷达高度计的卫星自主导航子***，基于紫外三轴敏感器的卫星自主导航子***，分别称为子***1，子***2，子***3，子***4。

具体包括以下步骤：

(1)建立基于轨道动力学的卫星运动状态方程；

(2)建立基于X射线探测器的卫星自主导航子***；

(3)建立基于红外地平仪和星敏感器的卫星自主导航子***；

(4)建立基于雷达高度计的卫星自主导航子***；

(5)建立基于紫外三轴敏感器的卫星自主导航子***；

(6)对上述四个子***的输出进行信息融合处理，输出导航信息；

(7)采用χ²检验法检测发生故障的***并隔离故障子***。

所述的步骤(1)中，建立包含了J₂项地球形状摄动和日月摄动的卫星轨道动力学方程，在J2000.0地心惯性系下，卫星轨道动力学方程为

                    X＝f(X，t)+W(t)                                           (1)

式中，X＝[x，y，z，v_x，v_y，v_z]^T为状态变量，x，y，z，v_x，v_y，v_z分别为卫星在X，Y，Z三个方向的位置和速度值；f(X，t)的形式根据卫星所设摄动的不同，仅考虑J2项的地球形状摄动和日、月的三体引力情况下，

$f(X,t)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\μ}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\μ}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\μ}}_m+(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

其中，μ_e，μ_s，μ_m分别为地球、太阳、月亮的引力常数，R_e为地球半径，J₂为地球引力系数，r_s＝[x_s，y_s，z_s]^T和r_m＝[x_m，y_m，z_m]^T分别表示太阳和月亮在地心惯性系下的位置，x_s，y_s，z_s为太阳在地心惯性系下的坐标值，x_m，y_m，z_m为月亮在地心惯性系下的坐标值，x，y，z，v_x，v_y，v_z分别为卫星在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为脉冲星旋转周期，T₅₀和T₁₀分别为脉冲波峰强度50％和10％时的脉冲宽度；

所述步骤(2)中的卫星位置与X射线脉冲星自转相位之间关系的计算方法为

$[\mathrm{\Δ\φ}+\mathrm{\ΔN}]\·\mathrm{\λ}-\hat{n}\·r_E=\hat{n}\·r+\mathrm{c\δ}{t}_{\mathrm{sat}}+{\mathrm{\Δ}}_{\mathrm{rel}}+v_k---\left(4\right)$

其中，Δφ＝φ_SSB-φ_sat，ΔN为卫星与太阳系质心SSB之间的整周模糊度值，λ为脉冲星的脉冲辐射的波长，

为脉冲星的视线方向，r_E为地球相对与SSB的位置矢量，r为卫星相对于地心的位置矢量，c为光速，δt_sat为卫星上的钟差，v_k为脉冲相位测量噪声，Δ_rel为相对论效应改正项，主要包括Roemer延迟改正、Shapiro延迟改正、太阳系行星总延迟改正；

所述步骤(2)中的整周模糊度的计算方法：

$\mathrm{\ΔN}=\mathrm{round}(\hat{n}\·(r_E+\tilde{r})/\mathrm{\λ}-\mathrm{\Δ\φ})---\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}{\mathrm{\λ}}_1^2{\mathrm{\σ}}_{\mathrm{\φ}1}^2& 0& 0\\ 0& {\mathrm{\λ}}_2^2{\mathrm{\σ}}_{\mathrm{\φ}2}^2& 0\\ 0& 0& {\mathrm{\λ}}_3^2{\mathrm{\σ}}_{\mathrm{\φ}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)$

其中，Z₁为星光角距观测量，Z₂为地心距观测量，r为卫星本体系下的地心矢量，s为卫星本体系下星光方向的单位矢量，v₁为星光角距测量噪声，v₂为地心距测量噪声，(x，y，z)为在地心惯性系下表示的卫星位置坐标值；

所述步骤(4)中，雷达高度计垂直安装于卫星的正下面，测量卫星平台到地表的距离，且大地水准面采用克莱饶椭球模型；

所述步骤(4)中，卫星星下点的地形高度h(x，y)由随机软件生成的随机地形提供，随机地形生成方法如下：

地形在高度方向上分解为地形基准平面高度h₀和在此平面上叠加的地形起伏Z(x，y)，即(x，y)点的实际海拔高度h(x，y)表示成：

                h(x，y)＝h₀+Z(x，y)                                 (10)

人工地形由二维一阶离散自回归过程产生：

          Z(x_i，y_i)＝a₁Z(x_i-Δx，y_i)+a₂Z(x_i，y_i-Δy)+a₃Z(x_i-Δx，y_i-Δy)+W(x_i，y_i)    (11)

其中，Z(x_i，y_i)为(x_i，y_i)点的地形起伏高度，Δx，Δy分别是沿x，y方向的采样间隔，W(x_i，y_i)为零均值白噪声序列，且W(x_i，y_i)～N(0，σ_w ²)；

根据地形的各向同性和平稳特性，令系数a₁＝a₂＝a，a₃＝b，且当a＝exp(-1/T_auc)，b＝-a²， ${\mathrm{\σ}}_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ 时，式(11)中两维随机函数Z(x_i，y_i)的相关函数为

$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(12\right)$

由此得到的两维随机过程，其均方差为d_RMS，按指数衰减，抗相关系数为T_auc，模拟地形在x方向和y方向的边界地形由两个一维随机过程分别产生，在x方向的地形边界上，y_i＝y₀，则有

           Z(x_i，y₀)＝aZ(x_i-Δx，y₀)+W(x_i，y₀)                     (13)

在y方向的地形边界上，x_i＝x₀，则有

        Z(x₀，y_i)＝aZ(x₀，y_i-Δy)+W(x₀，y_i)                 (14)

通过设定h₀、d_RMS、T_auc三个参数，利用软件编程产生不同的随机地形，将产生的随机地形h(x，y)叠加到大地水准面上；

所述步骤(5)中，以地心方向和地心距观测量的离散形式量测方程为

$z=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(15\right)$

其中，u_k为地心方向的单位矢量，r_k为地心距，v_k为量测噪声；

所述步骤(2)(3)(4)(5)中，四个子***的子滤波器，采用扩展卡尔曼滤波方法或者Unscented卡尔曼滤波方法，分别独立进行时间更新和测量更新；

所述步骤(6)中，对上述步骤(2)(3)(4)(5)四个子***的输出进行信息融合处理，采用的方法如下：

第i(i＝1，2，3，4)个子***的量测方程为：

           Z_i(k)＝H_i(k)X(k)+V_i(k)，i＝1，2，3，4                (16)

初始时刻全局状态估计为

其协方差阵为P_g0，根据信息守恒原理，将这一信息通过信息分配因子β_i分配到四个子滤波器和主滤波器，分配原则如下：

$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\mathrm{\β}}_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)={\mathrm{\β}}_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}{\mathrm{\Σ}}}{\mathrm{\β}}_i=1,$ 此处取β_m＝0，β₁＝β₂＝β₃＝β₄＝1/4，各子滤波器分别独立进行时间更新，并利用其对应传感器的量测信息进行量测更新，得四个子滤波器的局部估计值

和估计误差协方差阵P_i(k)i＝1，2，3，4；在主滤波器中，将子滤波器输出的局部信息按下式进行信息融合，得到全局状态估计信息为

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{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)$

即按照上式进行信息融合，输出全局状态估计

和其协方差矩阵P_g，将输出的全局状态估计

作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模型，作为卫星轨道动力学预报下一时刻卫星状态的初值；

所述步骤(7)中，χ²检验法进行***故障检测与隔离的方法如下：

在每个子滤波器中均加入故障检测和隔离的χ²方法，构成容错联邦滤波器，利用四个子***的故障检测模块，分别对基于X射线探测器的卫星自主导航子***、基于红外地平仪和星敏感器的卫星自主导航子***、基于雷达高度计的卫星自主导航子***、基于紫外三轴敏感器的卫星自主导航子***故障与否进行检测；

各子滤波器量测更新之后，进行故障检测和隔离，其方法是首先设计四个子***的故障检测函数D_i(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)＞T_D时，有故障，D(k)＜T_D时，无故障；T_D为预先设定的门限，可由误警率确定，当误警率P_fa＝α时，由P_fa＝P[λ_k＞T_D|H₀]＝α解出门限T_D。

本发明与现有技术相比的优点在于：

(1)本发明将多种传感器的量测信息融合得到精度更高的导航信息，具有为卫星提供高精度自主导航信息的优点；同时本发明采用多传感器冗余配置，并采用***级的故障检测与隔离方法，提高了***的可靠性和容错能力。

(2)本发明中的基于X射线探测器的卫星自主导航子***，具有提供十米量级的自主定轨精度能力，且适用于深空探测器的高精度自主导航需求；基于雷达高度计的卫星自主导航子***，采用雷达高度计垂直安装的方式，使用合理的海平面模型和模拟地形，提高了卫星的导航精度。

(3)本发明中的每个自主导航子***，均可独立构成卫星自主导航***，因而可在不同的情况下选用不同的子***进行组合，构建新的组合导航***，也可以只使用其中的一种子***，作为卫星的自主导航***。适合在不同导航精度要求下，对导航***成本以及重量、体积的控制。

(4)本发明适用于月球探测器和火星探测器飞行各个阶段的高精度自主导航。在月球探测器和火星探测器飞行的各个阶段，星敏感器、紫外三轴敏感器、X射线探测器都能够为导航***提供足够的导航观测量，从而为航天器提供高精度导航信息。在绕月飞行和绕火星飞行阶段，雷达高度计测距信息提供航天器与月球或火星表面之间的距离，从而确定绕月或绕火飞行航天器的轨道高度，为航天器提供高精度导航信息，实现深空探测飞行器完全自主式的导航。

附图说明

图1为本发明的组合导航***容错联邦滤波器结构图；

图2为本发明中基于X射线探测器的卫星自主导航子***的脉冲星导航示意图。

具体实施方式

基于多信息融合的卫星自主导航***由四个子***和一个主滤波器构成，如图1所示，本发明的具体实施方法如下：

一、建立基于轨道动力学的卫星运动状态方程

(1)建立含日月摄动的卫星轨道动力学方程

在J2000.0地心惯性系下，卫星轨道动力学方程(***状态方程)为

                    X＝f(X，t)+W(t)                                             (1)

式中，X＝[x，y，z，v_x，v_y，v_z]^T为状态变量，x，y，z，v_x，v_y，v_z分别为卫星在X，Y，Z三个方向的位置和速度值。f(X，t)的形式根据卫星所设摄动的不同，在仅考虑J2项的地球形状摄动和日、月的三体引力情况下，

$f(X,t)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\μ}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\μ}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\μ}}_m+(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

其中，μ_e，μ_s，μ_m分别为地球、太阳、月亮的引力常数，R_e为地球半径，J₂为地球引力系数，r_s＝[x_s，y_s，z_s]^T和r_m＝[x_m，y_m，z_m]^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为脉冲星旋转周期，T₅₀和T₁₀分别为脉冲波峰强度50％和10％时的脉冲宽度。计算候选脉冲星的品质因数值，Q值越大，脉冲星的导航属性越好，越可以用作导航脉冲星。

(3)建立以脉冲星自转相位(TOA)为观测量的量测方程

图2为绕地卫星的脉冲星导航示意图。其中，Sat、E、SSB分别表示卫星、地球和太阳系质心，为某颗脉冲星的视线方向，T为一个脉冲周期；r_E和r_{SSB_sat}分别表示地球、卫星相对于SSB的位置矢量，r表示卫星相对于地心的位置矢量。观测历元t时卫星处测得的脉冲相位值为φ_sat，太阳系质心的相位值为φ_SSB。

定义在太阳系质心的脉冲星钟模型为

$\mathrm{\Φ}\left(t\right)=\mathrm{\φ}\left(t_0\right)+f(t-t_0)+\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\frac{f^{\left(m\right)}{(t-t_0)}^{m+1}}{(m+1)!}---\left(4\right)$

式中，φ(t₀)为参考历元t₀时的脉冲星相位，f是脉冲星在t₀时的自转频率，f^(m)为f的m阶导数(一般取m＝1，2，3)。由(4)式预报可以得到t时太阳系质心处的脉冲相位值Φ_SSB，取其小数部分φ_SSB(0≤φ_SSB＜1)。

卫星与SSB之间的整周模糊度为ΔN，脉冲星的脉冲辐射的波长为λ(λ＝c·T，c为光速)，则测量方程为

$[\mathrm{\Δ\φ}+\mathrm{\ΔN}]\·\mathrm{\λ}-\hat{n}\·r_E=\hat{n}\·r+\mathrm{c\δ}{t}_{\mathrm{sat}}+{\mathrm{\Δ}}_{\mathrm{rel}}+v_k---\left(5\right)$

其中，Δφ＝φ_SSB-φ_sat，v_k为脉冲相位测量噪声。Δ_rel为相对论效应改正项，主要包括Roemer延迟改正、Shapiro延迟改正、太阳系行星总延迟改正。地球矢量r_E由JPL星历表提供。使用三颗脉冲星的观测量，则量测方程为

$\left[\begin{array}{c}({\mathrm{\Δ\φ}}_1+{\mathrm{\ΔN}}_1)\·{\mathrm{\λ}}_1-{\hat{n}}_1\·r_E\\ ({\mathrm{\Δ\φ}}_2+{\mathrm{\ΔN}}_2)\·{\mathrm{\λ}}_2-{\hat{n}}_2\·r_E\\ ({\mathrm{\Δ\φ}}_3+\mathrm{\Δ}{N}_3)\·{\mathrm{\λ}}_3-{\hat{n}}_3\·r_E\end{array}\right]=H\·\left[\begin{array}{c}r\\ \mathrm{c\δ}{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{\ΔN}=\mathrm{round}(\hat{n}\·(r_E+\tilde{r})/\mathrm{\λ}-\mathrm{\Δ\φ})---\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}{\mathrm{\λ}}_1^2{\mathrm{\σ}}_{\mathrm{\φ}1}^2& 0& 0\\ 0& {\mathrm{\λ}}_2^2{\mathrm{\σ}}_{\mathrm{\φ}2}^2& 0\\ 0& 0& {\mathrm{\λ}}_3^2{\mathrm{\σ}}_{\mathrm{\φ}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{{\mathrm{\σ}}_x^2+{\mathrm{\σ}}_y^2+{\mathrm{\σ}}_z^2+{\mathrm{\σ}}_t^2}---\left(9\right)$

其中，t_race表示求矩阵的迹，σ_x ²，σ_y ²，σ_z ²，σ_t ²为矩阵C的对角线元素，分别表示x，y，z轴的距离以及时间变量t的测量精度。

(6)子滤波器的滤波方法

子滤波器采用非线性卡尔曼滤波方法，可以考虑扩展卡尔曼滤波方法(EKF)或者Unscented卡尔曼滤波分别独立进行时间更新和测量更新。

扩展卡尔曼滤波处理的过程如下：

对状态方程和量测方程进行离散化，得到状态方程的离散型线性干扰方程以及离散型量测方程为

              δX_k+1＝φ_k+1，kδX_k+W_k                       (10)

              Z_k+1＝H_k+1X_k+1+V_k+1                           (11)

式中φ_k+1，k为状态转移矩阵，Z_k+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),$ W_k与V_k+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}={\mathrm{\φ}}_{k+1/k}{P}_{k/k}{\mathrm{\φ}}_{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为状态转移矩阵，K_k+1为卡尔曼增益系数，H_k+1为量测方程系数矩阵，P_k+1，k为最优预测估值误差协方差阵，P_k+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个采样点为

${\mathrm{\χ}}_0=\hat{X}$

${\mathrm{\χ}}_i=\hat{X}+{\left(\sqrt{(n+\mathrm{\λ})P_{k/k}}\right)}_i,i=\mathrm{1,2},...,n---\left(14\right)$

${\mathrm{\χ}}_i=\hat{X}-{\left(\sqrt{(n+\mathrm{\λ})P_{k/k}}\right)}_i,i=n+1,...,2n$

每个采样点对应的权重为

$W_0^{\left(m\right)}=\mathrm{\λ}/(n+\mathrm{\λ})$

$W_0^{\left(c\right)}=\mathrm{\λ}/(n+\mathrm{\λ})+(1-{\mathrm{\α}}^2+\mathrm{\β})---\left(15\right)$

$W_i^{\left(m\right)}=W_i^{\left(c\right)}=\mathrm{\λ}/\left\{2(n+\mathrm{\λ})\right\}i=\mathrm{1,2},...,2n$

其中，λ＝α²(n+k)-n是一个标量参数，常数α确定采样点围绕的分布特征，通常被设置为一个小的正数(1≥α≥10^-4)，常数k是一个标量参数，通常被设置为2或3-n。W_i ^(m)为均值加权所用权值，W_i ^(c)为协方差加权所用权值。Unscented卡尔曼滤波的基本设施过程如下：

计算采样点

${\mathrm{\χ}}_k=\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{x}}_k+\sqrt{n+\mathrm{\λ}}\sqrt{P_{k/k}}& {\hat{X}}_k-\sqrt{n+\mathrm{\λ}}\sqrt{P_{k/k}}\end{array}\right]---\left(16\right)$

时间更新

$\left\{\begin{array}{c}{\mathrm{\χ}}_{k+1/k}^{*}=F({\mathrm{\χ}}_k,k),{\hat{X}}_{k+1/k}=\underset{0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{\mathrm{\χ}}_{k+1/k}^{*}\\ P_{k+1/k}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{i,k+1/k}^{*}-{\hat{X}}_{k+1/k}]{[{\mathrm{\χ}}_{i,k+1/k}^{*}-{\hat{X}}_{k+1/k}]}^T+Q_{k+1}\\ {\mathrm{\χ}}_{k+1/k}=\left[\begin{array}{ccc}{\widehat{\mathrm{\χ}}}_{k+1/k}& {\hat{X}}_{k+1/k}+\mathrm{\γ}\sqrt{P_{k+1/k}}& {\hat{X}}_{k+1/k}-\mathrm{\γ}\sqrt{P_{k+1/k}}\end{array}\right]\\ z_{k+1/k}=H\left({\mathrm{\χ}}_{k+1/k}\right),{\hat{Z}}_{k+1/k}=\underset{0}{\overset{2n}{\mathrm{\Σ}}}{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}{\mathrm{\Σ}}}{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}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}[{\mathrm{\χ}}_{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,$ W_k与V_k+1互不相关。

利用状态方程(1)和量测方程(6)组成子滤波器1，对方程(1)和(6)分别进行离散化，采用扩展卡尔曼滤波方法或者Unscented卡尔曼滤波方法进行信息融合处理，输出子滤波器得到的卫星位置、速度信息。

三、建立基于红外地平仪和星敏感器的卫星自主导航子***

(7)建立以星光角距和地心距为观测量的量测方程

星敏感器观测导航恒星，可以确定星光在卫星本体坐标系中的方向；利用红外地平仪测得地心矢量在卫星本体坐标系中的方向；这样得到星光矢量和地心矢量之间的夹角，即星光角距。红外地平仪观测地球红外图像，获得地球半张角ρ，地球半径R_e为已知，从而可解算出地心距 $r=\frac{R_e}{\sin \mathrm{\ρ}},$ 从而以星光角距和地心距为观测量的量测方程为

$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)$

其中，Z₁为星光角距观测量，Z₂为地心距观测量，r为卫星本体系下的地心矢量，s为卫星本体系下星光方向的单位矢量，v₁为星光角距测量噪声，v₂为地心距测量噪声，(x，y，z)为在地心惯性系下表示的卫星位置坐标值。

(8)量测方程(19)、(20)构成的量测方程组，利用状态方程(2)和量测方程组构建如图1所示的子滤波器2，重复上述步骤(6)的扩展卡尔曼滤波方法或者Unscented卡尔曼滤波处理方法进行数据处理，输出卫星的位置、速度估计值。

四、建立基于雷达高度计的卫星自主导航子***

(9)建立以星光角距和轨道高度为观测量的量测方程

雷达高度计垂直按照在卫星的正下方，测得卫星平台到地表的距离，即卫星的轨道高度量测信息。大地水准面采用克莱饶椭球模型，并将海平面的起伏作为修正项加以考虑。克莱饶水准面可以表示为下式

饶水准面可以表示为下式

其中，_e为地心纬度，P₂为二次勒得让多项式，J₂为地球带谐项系数，R为地球平均半径值，参数α_e为

${\mathrm{\α}}_e=\frac{\left[\frac{2}{3}{J}_2(\frac{{\mathrm{\α}}_e}{R}+\frac{m}{2})\right]}{[1-\frac{2}{3}m]m}\≈0.003449963---\left(22\right)$

则以轨道高度为观测量的量测方程为

$Z_2=r-R[1-\frac{{2\mathrm{\α}}_e}{3}\left(0.5\right){\left(\frac{z}{r}\right)}^2-\frac{1}{2}]-h(x,y)+v_2---\left(23\right)$

式中，Z₂为雷达高度计测得的卫星轨道高度，r为卫星到地心的距离，v₂为测量噪声，h(x，y)为星下点的地形高度，由随机地形软件生成。以星光角距为观测量的量测方程同方程(19)。

(10)随机地形生成方法

建立模拟地形，模拟实际地形的变化。地形在高度方向上分解为地形基准平面高度h₀和在此平面上叠加的地形起伏Z(x，y)，即(x，y)点的实际海拔高度h(x，y)可表示成：

                  h(x，y)＝h₀+Z(x，y)                                   (24)

人工地形由二维一阶离散自回归过程产生：

          Z(x_i，y_i)＝a₁Z(x_i-Δx，y_i)+a₂Z(x_i，y_i-Δy)+a₃Z(x_i-Δx，y_i-Δy)+W(x_i，y_i)    (25)

其中，Z(x_i，y_i)为(x_i，y_i)点的地形起伏高度，Δx，Δy分别是沿x，y方向的采样间隔，W(x_i，y_i)为零均值白噪声序列，W(x_i，y_i)～N(0，σ_w ²)。

根据地形的各向同性和平稳特性，可令系数a₁＝a₂＝a，a₃＝b，且当a＝exp(-1/T_auc)，b＝-a²， ${\mathrm{\σ}}_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ 时，式(25)中两维随机函数Z(x_i，y_i)的相关函数为

$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(26\right)$

由此得到的两维随机过程，其均方差为d_RMS，按指数衰减，抗相关系数为T_auc。模拟地形在x方向和y方向的边界地形由两个一维随机过程分别产生。若在x方向的地形边界上，y_i＝y₀，则有

              Z(x_i，y₀)＝aZ(x_i-Δx，y₀)+W(x_i，y₀)                         (27)

若在y方向的地形边界上，x_i＝x₀，则有

              Z(x₀，y_i)＝aZ(x₀，y_i-Δy)+W(x₀，y_i)                (28)

通过设定h₀、d_RMS、T_auc三个参数，利用软件编程可以产生不同的随机地形。将产生的随机地形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)$

其中，u_k为地心方向的单位矢量，r_k为地心距，v_k为量测噪声。

利用状态方程(2)和量测方程(29)构建如图1所示的子滤波器4，重复上述步骤(6)的扩展卡尔曼滤波方法或者Unscented卡尔曼滤波处理方法进行数据融合，输出卫星的位置、速度估计值。

六、主滤波器的信息融合、故障检测与隔离方法

(13)主滤波器的信息融合处理

第i(i＝1，2，3，4)个子***的量测方程为：

          Z_i(k)＝H_i(k)X(k)+V_i(k)，(i＝1，2，3，4)              (30)

初始时刻全局状态估计为

其协方差阵为P_g0。根据信息守恒原理，将这一信息通过信息分配因子β_i分配到四个子滤波器和主滤波器，分配原则如下：

$\left\{\begin{array}{c}{P}_{i0}^{-1}\left(k\right)={\mathrm{\β}}_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)={\mathrm{\β}}_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}{\mathrm{\Σ}}}{\mathrm{\β}}_i=1,$ 此处取β_m＝0，β₁＝β₂＝β₃＝β₄＝1/4。各子滤波器分别独立进行时间更新，并利用其对应传感器的量测信息进行量测更新，得四个子滤波器的局部估计值

和估计误差协方差阵P_i(k)(i＝1，2，3，4)，具体措施如方程(12)或(17)和(18)。在主滤波器中，将子滤波器输出的局部信息按下式进行信息融合，得到全局状态估计信息为

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{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)式进行信息融合，输出全局状态估计

和其协方差矩阵P_g。将输出的全局状态估计

作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模型，作为卫星轨道动力学预报下一时刻卫星状态的初值。

联邦滤波中，主滤波器和四个子滤波器联合工作，信息分配只在初始时刻进行，子滤波器独立工作，不接收主滤波器的反馈信息。

(14)组合导航的故障检测与隔离方法

采用***级的故障检测与隔离方法，排除发生故障的传感器及相应子***。在每个子滤波器中均加入故障检测和隔离的χ²方法，构成容错联邦滤波器，其结构如图1所示。图1中FDI 1、FDI 2、FDI 3、FDI 4分别为四个子***的故障检测模块，分别对子***1、2、3、4故障与否进行检测，并决定是否隔离，以便主滤波器对正常的子***进行重组，使部分导航传感器故障导致子***失效后，综合导航***仍能继续工作，具容错性。

各子滤波器分别独立进行时间更新和量测更新。各子滤波器量测更新之后，进行故障检测和隔离。其方法是首先设计四个子***的故障检测函数D_i(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)＞T_D时，有故障；D(k)＜T_D时，无故障。T_D为预先设定的门限，可由误警率确定。当误警率P_fa＝α时，P_fa＝P[λ_k＞T_D|H₀]＝α，解出门限T_D。若判断出子***无故障，则把其滤波结果送到主滤波器；若子***故障，则滤波结果不送到主滤波器，并对其隔离，直至该子***恢复正常。

在各子滤波器计算出各自的局部估计，并进行故障检测后，由主滤波器把输出正常的子滤波器的公共状态进行信息融合，得到全局估计信息。

(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，v_x，v_y，v_z的最优估计。

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

Claims (1)

1.一种基于多信息融合的航天器组合导航方法，其特征在于：具体包括一下步骤：

(1)建立基于轨道动力学的卫星运动状态方程；

(2)建立基于X射线探测器的卫星自主导航子***；

(3)建立基于红外地平仪和星敏感器的卫星自主导航子***；

(4)建立基于雷达高度计的卫星自主导航子***；

(5)建立基于紫外三轴敏感器的卫星自主导航子***；

(6)对上述步骤(2)(3)(4)(5)四个子***的输出进行信息融合处理，输出导航信息；

(7)采用χ²检验法检测发生故障的***并隔离故障子***；

所述步骤(1)中，建立包含J₂项地球形状摄动和日月摄动的卫星轨道动力学方程，在J2000.0地心惯性系下，卫星轨道动力学方程为

              X＝f(X，t)+W(t)                               (1)

式中，X＝[x，y，z，v_x，v_y，v_z]^T为状态变量，x，y，z，v_x，v_y，v_z分别为卫星在X，Y，Z三个方向的位置和速度值；f(X，t)的形式根据卫星所设摄动的不同，仅考虑J2项的地球形状摄动和日、月的三体引力情况下，

$f(X,t)=\left[\begin{array}{c}{v}_x;v_y;v_z\\ -\frac{{\mathrm{\μ}}_{e}x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{x_s-x}{|r_s-r|}-\frac{x_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{x_m-x}{|r_m-r|}-\frac{x_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+{\mathrm{\μ}}_s(\frac{y_s-y}{|r_s-r|}-\frac{y_s}{\left|r_s\right|})+{\mathrm{\μ}}_m(\frac{y_m-y}{|r_m-r|}-\frac{y_m}{\left|r_m\right|})\\ -\frac{{\mathrm{\μ}}_{e}z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+{\mathrm{\μ}}_s(\frac{z_s-z}{|r_s-r|}-\frac{z_s}{\left|r_s\right|})+{\mathrm{\μ}}_m+(\frac{z_m-z}{|r_m-r|}-\frac{z_m}{\left|r_m\right|})\end{array}\right]---\left(2\right)$

其中，μ_e，μ_s，μ_m分别为地球、太阳、月亮的引力常数，R_e为地球半径，J₂为地球引力系数，r_s＝[x_s，y_s，z_s]^T和r_m＝[x_m，y_m，z_m]^T分别表示太阳和月亮在地心惯性系下的位置，x_s，y_s，z_s为太阳在地心惯性系下的坐标值，x_m，y_m，z_m为月亮在地心惯性系下的坐标值，x，y，z，v_s，v_y，v_z分别为卫星在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}-T_{50}}{T_{50}(T_{10}-T_{50})}}---\left(3\right)$

式中，SNR为脉冲星的信噪比，T为脉冲星旋转周期，T₅₀和T₁₀分别为脉冲波峰强度50％和10％时的脉冲宽度；

所述步骤(2)中的卫星位置与X射线脉冲星自转相位之间关系的计算方法为

$[\mathrm{\Δ\φ}+\mathrm{\ΔN}]\·\mathrm{\λ}-\hat{n}\·r_{\mathrm{\ϵ}}=\hat{n}\·r+\mathrm{c\δ}{t}_{\mathrm{sat}}+{\mathrm{\Δ}}_{\mathrm{rel}}+v_k---\left(4\right)$

其中，Δφ＝φ_SSB-φ_sat，ΔN为卫星与太阳系质心SSB之间的整周模糊度值，λ为脉冲星的脉冲辐射的波长，为脉冲星的视线方向，r_s为地球相对与SSB的位置矢量，r为卫星相对于地心的位置矢量，c为光速，δt_sat为卫星上的钟差，v_k为脉冲相位测量噪声，Δ_rel为相对论效应改正项，主要包括Roemer延迟改正、Shapiro延迟改正、太阳系行星总延迟改正；

所述步骤(2)中的整周模糊度的计算方法：

$\mathrm{\ΔN}=\mathrm{round}(\hat{n}\·(r_E+\tilde{r})/\mathrm{\λ}-\mathrm{\Δ\φ})---\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}{\mathrm{\λ}}_1^2{\mathrm{\σ}}_{\mathrm{\φ}1}^2& 0& 0\\ 0& {\mathrm{\λ}}_2^2{\mathrm{\σ}}_{\mathrm{\φ}2}^2& 0\\ 0& 0& {\mathrm{\λ}}_3^2{\mathrm{\σ}}_{\mathrm{\φ}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)$

其中，Z₁为星光角距观测量，Z₂为地心距观测量，r为卫星本体系下的地心矢量，s为卫星本体系下星光方向的单位矢量，v₁为星光角距测量噪声，v₂为地心距测量噪声，(x，y，z)为在地心惯性系下表示的卫星位置坐标值；

所述步骤(4)中，雷达高度计垂直安装于卫星的正下面，测量卫星平台到地表的距离，且大地水准面采用克莱饶椭球模型；

所述步骤(4)中，卫星星下点的地形高度h(x，y)由随机软件生成的随机地形提供，随机地形生成方法如下：

地形在高度方向上可以分解为地形基准平面高度h₀和在此平面上叠加的地形起伏Z(x，y)，即(x，y)点的实际海拔高度h(x，y)可表示成：

                 h(x，y)＝h₀+Z(x，y)                                (10)

人工地形由二维一阶离散自回归过程产生：

    Z(x_i，y_i)＝a₁Z(x_i-Δx，y_i)+a₂Z(x_i，y_i-Δy)+a₃Z(x_i-Δx，y_i-Δy)+W(x_i，y_i)  (11)

其中，Z(x_i，y_i)为(x_i，y_i)点的地形起伏高度，Δx，Δy分别是沿x，y方向的采样间隔，W(x_i，y_i)为零均值白噪声序列，且W(x_i，y_i)～N(0，σ_w ²)；

根据地形的各向同性和平稳特性，令系数a₁＝a₂＝a，a₃＝b，且当a＝exp(-1/T_auc)，b＝-a²， ${\mathrm{\σ}}_w^2={\left(d_{\mathrm{RMS}}\right)}^2\exp (-(i+j)/T_{\mathrm{auc}})$ 时，式(11)中两维随机函数Z(x_i，y_i)的相关函数为

$R_{\mathrm{ij}}=\left(d_{\mathrm{RMS}}\right)\exp (-\frac{i+j}{T_{\mathrm{auc}}})---\left(12\right)$

由此得到的两维随机过程，其均方差为d_RMS，按指数衰减，抗相关系数为T_auc，模拟地形在x方向和y方向的边界地形由两个一维随机过程分别产生，在x方向的地形边界上，y_i＝y₀，则有

       Z(x_i，y₀)＝aZ(x_i-Δx，y₀)+W(x_i，y₀)                         (13)

在y方向的地形边界上，x_i＝x₀，则有

       Z(x₀，y_i)＝aZ(x₀，y_i-Δy)+W(x₀，y_i)                         (14)

通过设定h₀、d_RMS、T_auc三个参数，利用软件编程产生不同的随机地形，将产生的随机地形h(x，y)叠加到大地水准面上；

所述步骤(5)中，以地心方向和地心距观测量的离散形式量测方程为

$z=\left[\begin{array}{c}{u}_k\\ r_k\end{array}\right]+v_k---\left(15\right)$

其中，u_k为地心方向的单位矢量，r_k为地心距，v_k为量测噪声；

所述步骤(2)(3)(4)(5)中，四个子***的子滤波器，采用扩展卡尔曼滤波方法或者Unscented卡尔曼滤波方法，分别独立进行时间更新和测量更新；

所述步骤(6)中，对上述步骤(2)(3)(4)(5)四个子***的输出进行信息融合处理，采用的方法如下：

第i(i＝1，2，3，4)个子***的量测方程为：

           Z_i(k)＝H_i(k)X(k)+V_i(k)，i＝1，2，3，4               (16)

初始时刻全局状态估计为

其协方差阵为P_g0，根据信息守恒原理，将这一信息通过信息分配因子β_i分配到四个子滤波器和主滤波器，分配原则如下：

$\left\{\begin{array}{cc}{P}_{i0}^{-1}\left(k\right)={\mathrm{\β}}_i{P}_{g0}^{-1}\left(k\right)& \\ {\hat{X}}_{i0}\left(k\right)={\hat{X}}_{g0}\left(k\right)& i=\mathrm{1,2,3,4},m\\ Q_{i0}^{-1}\left(k\right)={\mathrm{\β}}_i{Q}_{g0}^{-1}\left(k\right)& \end{array}\right.---\left(17\right)$

其中，i＝1，2，3，4表示四个子滤波器，m表示主滤波器，信息分配因子β_i满足信息守恒原理： $\underset{i=1}{\overset{4,m}{\mathrm{\Σ}}}{\mathrm{\β}}_i=1,$ 此处取β_m＝0，β₁＝β₂＝β₃＝β₄＝1/4，各子滤波器分别独立进行时间更新，并利用其对应传感器的量测信息进行量测更新，得四个子滤波器的局部估计值

和估计误差协方差阵P_i(k)i＝1，2，3，4；在主滤波器中，将子滤波器输出的局部信息按下式进行信息融合，得到全局状态估计信息为

$\left\{\begin{array}{c}{\hat{X}}_g=P_g\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{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)$

即按照上式进行信息融合，输出全局状态估计

和其协方差矩阵P_g，将输出的全局状态估计

作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模型，作为卫星轨道动力学预报下一时刻卫星状态的初值；

所述步骤(7)中，χ²检验法进行***故障检测与隔离的方法如下：

在每个子滤波器中均加入故障检测和隔离的χ²方法，构成容错联邦滤波器，利用四个子***的故障检测模块，分别对基于X射线探测器的卫星自主导航子***、基于红外地平仪和星敏感器的卫星自主导航子***、基于雷达高度计的卫星自主导航子***、基于紫外三轴敏感器的卫星自主导航子***故障与否进行检测；

各子滤波器量测更新之后，进行故障检测和隔离，其方法是首先设计四个子***的故障检测函数D_i(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)＞T_D时，有故障，D(k)＜T_D时，无故障；T_D为预先设定的门限，可由误警率确定，当误警率P_fa＝α时，由P_fa＝P[λ_k＞T_D|H₀]＝α解出门限T_D。

Images