[发明专利]基于多信息融合的航天器组合导航方法

权利要求书

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)=vx;vy;vz-μexr3[1-J2(Rer)2(7.5z2r2-1.5)]+μs(xs-x|rs-r|-xs|rs|)+μm(xm-x|rm-r|-xm|rm|)-μeyr3[1-J2(Rer)2(7.5z2r2-1.5)]+μs(ys-y|rs-r|-ys|rs|)+μm(ym-y|rm-r|-ym|rm|)-μezr3[1-J2(Rer)2(7.5z2r2-4.5)]+μs(zs-z|rs-r|-zs|rs|)+μm+(zm-z|rm-r|-zm|rm|)---(2)]]>

其中，μ_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.8SNRT15T10+7T505T10-T50T50(T10-T50)---(3)]]>

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

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

[Δφ+ΔN]·λ-n^·rϵ=n^·r+cδtsat+Δrel+vk---(4)]]>

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

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

ΔN=round(n^·(rE+r~)/λ-Δφ)---(5)]]>

其中，为卫星动力学方程提供的卫星位置预报值；所述步骤(2)中的GDOP值计算方法：

GDOP=trace(C)---(6)]]>

其中，trace表示求矩阵的迹，C为位置误差协方差矩阵，由下式算得

C={(HTH)-1HT}λ12σφ12000λ22σφ22000λ32σφ32{(HTH)-1HT}T---(7)]]>

其中，H为量测方程的系数矩阵，λ_i(i＝1，2，3)为第i颗脉冲星的脉冲辐射波长，σ_φi为第i颗脉冲星的相位测量精度；

所述步骤(3)中，以星光角距和地心距为观测量的量测方程为

Z1=arccos(-r·sr)+v1---(8)]]>

Z2=x2+y2+z2+v2---(9)]]>

其中，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²，σw2=(dRMS)2exp(-(i+j)/Tauc)]]>时，式(11)中两维随机函数Z(x_i，y_i)的相关函数为

Rij=(dRMS)exp(-i+jTauc)---(12)]]>

由此得到的两维随机过程，其均方差为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=ukrk+vk---(15)]]>

其中，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分配到四个子滤波器和主滤波器，分配原则如下：

Pi0-1(k)=βiPg0-1(k)X^i0(k)=X^g0(k)i=1,2,3,4,mQi0-1(k)=βiQg0-1(k)---(17)]]>

其中，i＝1，2，3，4表示四个子滤波器，m表示主滤波器，信息分配因子β_i满足信息守恒原理：Σi=14,mβi=1,]]>此处取β_m＝0，β₁＝β₂＝β₃＝β₄＝1/4，各子滤波器分别独立进行时间更新，并利用其对应传感器的量测信息进行量测更新，得四个子滤波器的局部估计值和估计误差协方差阵P_i(k)i＝1，2，3，4；在主滤波器中，将子滤波器输出的局部信息按下式进行信息融合，得到全局状态估计信息为

X^g=PgΣi=14Pi-1X^iPg-1=P1-1+P2-1+P3-1+P4-1---(18)]]>

即按照上式进行信息融合，输出全局状态估计和其协方差矩阵P_g，将输出的全局状态估计作为卫星轨道动力学方程预测卫星状态的初值，并反馈到卫星的轨道动力学模型，作为卫星轨道动力学预报下一时刻卫星状态的初值；

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

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

各子滤波器量测更新之后，进行故障检测和隔离，其方法是首先设计四个子系统的故障检测函数D_i(k)i＝1，2，3，4，

di(k)=zi(k)-Hi(k)X^i(k/k-1)---(19)]]>

Si(k)=Hi(k)Pi(k/k-1)HiT(k)+Ri(k)---(20)]]>

Di(k)=diT(k)Si-1(k)di(k)---(21)]]>

然后进行故障判断：D(k)＞T_D时，有故障，D(k)＜T_D时，无故障；T_D为预先设定的门限，可由误警率确定，当误警率P_fa＝α时，由P_fa＝P[λ_k＞T_D|H₀]＝α解出门限T_D。