CN103344245B - 火星进入段imu和甚高频无线电组合导航的ud-skf方法 - Google Patents

Abstract

一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，步骤如下：一、建立火星进入段的动力学方程；二、建立火星进入段的量测方程；三、对上述动力学方程（5）和量测方程（6）进行离散化；四、采用UD-SKF滤波算法，并输出导航信息。本发明所采用的UD-SKF滤波算法与传统的EKF相比，不仅将量测方程的偏差信息融合到滤波过程中，极大地改善了滤波效果，提高了导航精度，而且利用UD分解，减少了滤波过程中出现较大误差甚至发散的情况，增强滤波过程的稳定性，从而能对探测器进行实时高效地估计导航状态。

Description

火星进入段IMU和甚高频无线电组合导航的UD-SKF方法

技术领域

本发明涉及一种火星进入段组合导航的自主导航方法，具体涉及以一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法,属于航天自主导航技术领域.

背景技术

2012年8月，美国的火星探测器“好奇号”成功登陆火星，着陆点距离预定目标点2385米。但是，未来的火星探测导航任务要求探测器要达到定点着陆的要求，也即是着陆精度在1000米之内。在着陆之前，火星探测器会经历火星进入段、下降段和着陆段（简称EDL），其中火星进入段的导航对探测器的最终着陆精度影响最为显著。这主要是由于火星进入段气动加热环境恶劣，且探测器包裹在整流罩内，导致很多导航敏感器无法正常工作，只能主要依靠惯性测量单元（简称IMU）提供导航信息；其次由于火星和地球之间的通信大延迟，无法利用深空网实时提供较为准确的导航信息。因此，如何规划火星进入段的自主导航方案和自主导航滤波算法成为火星探测器能否定点着陆的首要问题。

目前，大量学者提出了比较新颖的火星进入段自主导航方案，较著名的是Levesque等提出的基于星间甚高频无线电测量通信的自主导航方案（参见Levesque,JF.Advancednavigationandguidanceforhigh-precisionplanetarylandingonMars.Ph.D.Thesis,DepartmentofElectricalEngineering,UniversitedeSherbrooke,2006.）。此导航方案解决了导航信息不足导致的不可观性的问题，但是初始误差、大气模型的不确定性、探测器参数和气动特性的不确定性、信标的位置偏差和信标与探测器之间的通信偏差等导致此自主导航方案也存在导航精度不一定满足未来火星探测任务的要求。

实时高效的自主导航滤波算法也是解决这类问题的有效途径。目前火星进入段导航方式采用的是对IMU的输出进行积分，递推得到相应的探测器状态，但是此种方法收惯性单元漂移和偏差的影响较大，随着时间推移导航误差积累较大。由于火星进入段动力学具有很强的非线性，因此对探测器进行实时估计的算法必须是非线性滤波算法。工程上应用最多的是扩展卡尔曼滤波（简称EKF），此方法易于工程实现，但是由于线性化产生的截断误差会导致滤波精度和稳定性降低（参见Lightsey,EG,Mogensen,AE,Burkhart,PD,etal.,Real-timenavigationformarsmissionsusingthemarsnetwork.JournalofSpacecraftandRockets,2008.45(3):519-533.）。有学者提出了无迹卡尔曼滤波（简称UKF），并应用于火星进入段自主导航中，估计精度较EKF高，但是计算量较EKF大（参见Levesque,JF.Advancednavigationandguidanceforhigh-precisionplanetarylandingonMars.Ph.D.Thesis,DepartmentofElectricalEngineering,UniversitedeSherbrooke,2006.）。但是这两类方法无法有效地处理上述问题的不确定性和偏差，即使采用将这些不确定性或偏差进行状态扩维，然后进行估计，也能取得较好地效果。但是这种处理方式无疑增加了计算量和计算复杂度，不利于对探测器的状态进行实时估计。

本发明在提出的火星进入段IMU和甚高频无线电组合导航的UD-SKF滤波将这些偏差的方差加入状态方差矩阵的计算中，但是不估计这些偏差。此自主导航滤波算法成功地考虑了这些偏差的影响，但是又避免估计它们，有效地提高了导航精度，同时也大量的减少了计算量和计算复杂度，有利于火星进入段对探测器的进入状态进行实时高效地估计。

发明内容

本发明的目的是提供一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法。针对现有的非线性自主导航滤波算法无法消除这些偏差带来较大误差，甚至导致发散的现象，或者进行状态扩维后，导致计算量和计算复杂度急剧增加的不利情况，本发明所提出的基于UD-SKF滤波算法成功地考虑了这些偏差的影响，但是又避免估计它们，有效地提高了导航精度，同时也大量的减少了计算量和计算复杂度，有利于火星进入段对探测器的进入状态进行实时高效地估计。

本发明涉及一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，它包括以下四个步骤（具体流程图参见图1）：

步骤一、建立火星进入段的动力学方程

在火星进入段，将探测器看作是一个无动力的质点。假设进入段的大气相对火星是静止的，同时忽略火星的自转，建立火星进入段的三自由度动力学方程：

$\stackrel{\·}{r}=v\sin \mathrm{\γ}$

$\stackrel{\·}{\mathrm{\θ}}=\frac{v\cos \mathrm{\γ}\sin \mathrm{\ψ}}{r\cos \mathrm{\λ}}$

$\stackrel{\·}{\mathrm{\λ}}=\frac{v\cos \mathrm{\γ}\cos \mathrm{\ψ}}{r}$

$\stackrel{\·}{v}=-D-g\left(r\right)\sin \mathrm{\γ}---\left(1\right)$

$\stackrel{\·}{\mathrm{\γ}}=(\frac{v}{r}-\frac{g\left(r\right)}{v})\cos \mathrm{\γ}+\frac{L}{v}\cos \mathrm{\φ}$

$\stackrel{\·}{\mathrm{\ψ}}=\frac{v}{r}\sin \mathrm{\ψ}\cos \mathrm{\γ}\tan \mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}\sin \mathrm{\φ}$

式中，r,v,θ,λ,γ,ψ分别表示探测器距离火星中心的半径距离（高度）、探测器的速度、经度、纬度、飞行路径角和航向角，火星重力加速度g(r)＝μ/r²，μ是火星重力加速度常数，D和L分别为阻力加速度和升力加速度：

$D=\frac{C_{D}S}{m}\stackrel{\‾}{q}---\left(2\right)$

$L=\frac{C_{L}S}{m}\stackrel{\‾}{q}---\left(3\right)$

式中C_D和C_L分别为阻力系数和升力系数，S为探测器参考面积，m为探测器质量，为动压，ρ是火星大气密度，此处采用指数大气密度模型：

$\mathrm{\ρ}={\mathrm{\ρ}}_0\exp \left\{\frac{r_0-r}{h_s}\right\}---\left(4\right)$

式中ρ₀是标称密度，r₀为火星径向基准位置，h_s为大气定标高度。

取状态向量为X=[rθλvγψ]^T，考虑到系统噪声向量w(t)，则动力学程（1）可简写为

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(5\right)$

式中，f(X(t),t)为系统非线性连续状态转移函数，w(t)为零均值的高斯白噪声。

步骤二、建立火星进入段的量测方程

本发明以基于星间甚高频无线电测量通信的自主导航方案为基础，也即是将探测器与火星在轨信标或火星表面固定信标之间的双向测距信息作为测量信息和IMU提供的测量信息一起组合为新的测量信息，建立火星进入段量测方程：

Z=h(X(t),t)+b+v(t)(6)

式中，

$h(X\left(t\right),t)={\left[\begin{array}{cc}{a}^v& \tilde{R}\end{array}\right]}^T$

a^v=[-D-D·L/D·sinφD·L/D·cosφ]^T(7)

$\tilde{R}={\[{\tilde{R}}_1,{\tilde{R}}_2,\·\·\·,{\tilde{R}}_m\]}^T$

v(t)=[η_aξ_R]^T

式中，a^v为速度系下的加速度，表示双向测距信息，维数m由有效的火星在轨信标和火星表面固定信标数目决定，b为常值偏差向量，v(t)为零均值的高斯白噪声，η_a为IMU测量信息的测量噪声向量，ξ_R为双向测距的测量噪声向量。

其中，IMU的输出是由加速度计测量得到的本体系三个轴向的线加速度，测量模型为

$\tilde{a}=a+b_a+{\mathrm{\η}}_a---\left(8\right)$

式中，表示各个轴的加速度测量值，a表示真实的加速度，b_a加速度计常值漂移，测量噪声η_a为零均值的高斯白噪声。

其中，双向测距指的是探测器与火星在轨信标或火星表面固定信标之间通过双向无线电通信测量得到的二者之间的距离，无线电通信双向测距示意图参见图2，可表示为:

$\tilde{R}=R+b_R+{\mathrm{\ξ}}_R---\left(9\right)$

$R=\sqrt{{(r_l-r_k)}^T(r_l-r_k)}---\left(10\right)$

式中，r_l和r_k分别表示火星质心惯性系下火星探测器和火星在轨信标(k=o)或火星表面固定信标(k=b)，b_R为双向测距的常值偏差，测距噪声ξ_R为零均值的高斯白噪声。

步骤三、对上述动力学方程（5）和量测方程（6）进行离散化

X_k+1=F(X_k)+w_k（11）

式中，k=1,2,3,…，F(X_k)为f(X(t),t)离散后的非线性状态转移函数，为h(X(t),t)离散后的非线性量测函数，w_k和v_k互不相关，且其方差矩阵分别为Q_k和R_k。

然后将方程（11）和（12）进行线性化，也即将式（11）中的非线性离散函数F(X_k)围绕估计值展开成泰勒级数，并略去二阶以上项，得线性化之后的动力学方程

X_k+1=Φ_k+1/kX_k+U_k+w_k（13）

式中，

${\mathrm{\Φ}}_{k+1/k}=\frac{\∂F\left(X_k\right)}{\∂X_k}{|}_{x_k={\hat{x}}_k}---\left(14\right)$

$U_k=F\left({\hat{X}}_k\right)-\frac{\∂F\left(X_k\right)}{\∂X_k}{|}_{x_k={\hat{x}}_k}\·{\hat{X}}_k---\left(15\right)$

然后，再将式（12）中的非线性离散函数围绕估计值和展开成泰勒级数，并略去二阶以上项，得线性化之后的增量量测方程

Z_k=H_kX_k+Y_k+v_k（16）

式中，

通过上述过程，就得到了线性化后的动力学方程和量测方程，式中U_k和Y_k为非随机的外作用项。

步骤四、采用UD-SKF滤波算法，并输出导航信息。

考虑式（6）中的常值偏差向量b未能精确建模，故施密特-卡尔曼滤波算法（SchmidtKalmanfilter,简称SKF）考虑将其方差融入到算法中，也即是考虑了偏差与状态的互协方差，但是不估计这些偏差。同时，在计算增益矩阵时采用UD分解以确保滤波的数值稳定性。

其中步骤四所采用的UD-SKF滤波算法实现步骤为：

1.将常值偏差向量b扩维进入状态向量，则动力学方程（13）和量测方程（16）变为

$\left[\begin{array}{c}{X}_{k+1}\\ b_{k+1}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\Φ}}_{k+1|k}& 0\\ 0& I\end{array}\right]\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+\left[\begin{array}{c}{w}_k\\ 0\end{array}\right]---\left(19\right)$

$Z_k=[H_k,I]\begin{array}{}\end{array}\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+v_k---\left(20\right)$

式中，我们考虑的常值偏差向量满足条件：b_k=b_k-1。其方差矩阵B₀为

B₀=Cov{b₀}=Cov{b_k}(21)

和偏差与状态的互协方差矩阵C_k为

$C_k=E\left\{{\tilde{X}}_k{b}^T\right\}=E\left\{(X_k-{\hat{X}}_k)b^T\right\},---\left(22\right)$

并且初始值为C₀=0。式中，为卡尔曼滤波第k步的状态估计。

同时，相应与动力学方程（19）和量测方程（20）第k步的误差方差阵为

式中，为C_k转置矩阵，P_k为状态X_k的误差方差阵

$P_k=E\left\{{\tilde{X}}_k{\tilde{X}}_k^T\right\}=E\left\{(X_k-{\hat{X}}_k){(X_k-{\hat{X}}_k)}^T\right\}---\left(24\right)$

然后，初始化状态向量X₀及误差方差阵P₀。

2.由第k步的状态估计可得，第k+1步的状态一步预测为

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1|k}{\hat{X}}_k---\left(25\right)$

3.第k+1步的一步预测误差方差矩阵为

则可得到状态和偏差的一步预测误差方差矩阵P_k+1|k和C_k+1|k

$P_{k+1|k}={\mathrm{\Φ}}_{k+1|k}{P}_k{\mathrm{\Φ}}_{k+1|k}^T+Q_k---\left(27\right)$

C_k+1|k=Φ_k+1|kC_k(28)

4.第k+1步的滤波增益矩阵为

其中，由于不需要估计偏差，故强制令偏差项的增益矩阵为零。

则状态的增益矩阵K_k+1为

$K_{k+1}=[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]{\mathrm{\Ω}}_{k+1}^{-1}---\left(30\right)$

${\mathrm{\Ω}}_{k+1}=H_{k+1}{P}_{k+1|k}{H}_{k+1}^T+C_{k+1|k}^T{H}_{k+1}^T+H_{k+1}{C}_{k+1|k}+B_0+R_{k+1}---\left(31\right)$

由于式（31）中的Ω_k+1在滤波过程中无法保持非负定性，此时式（30）中的增益矩阵K_k+1在求逆运算中会产生较大的误差，甚至会发散。为了保持滤波的数值稳定性，本发明采用UD分解来计算增益矩阵K_k+1，会起到很好的效果。UD分解的详细步骤如下：

1）将对称正定矩阵Ω_k分解成一个对角阵D和一个单位上三角阵U的形式：

Ω_k=UDU^T(32)

2）将式（32）代入是式（30）中，得到

$K_k{\mathrm{UDU}}^T=P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}---\left(33\right)$

${\mathrm{UDU}}^T{K}_k^T={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T---\left(34\right)$

3）令则有

${\mathrm{UDU}}^{T}X={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T---\left(35\right)$

4）可以通过以下的迭代过程得到

${\mathrm{UX}}_1={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T$

DX₂=X₁(36)

U^TX=X₂

5）最后令然后即可得到K_k。

5.第k+1步的状态估计为

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-H_{k+1}{\hat{X}}_{k+1/k})---\left(37\right)$

6.第k+1步的估计误差方差矩阵为

则状态的估计误差方差矩阵P_k+1为

$P_{k+1}=P_{k+1|k}-K_{k+1}{\mathrm{\Ω}}_{k+1}{K}_{k+1}^T---\left(39\right)$

和偏差与状态的互协方差矩阵为

C_k+1=C_k+1|k-K_k+1(H_k+1C_k+1|k+B₀)(40)

通过以上6步循环进行即可得到探测器的实时状态估计值，包含探测器的高度、速度、经度、纬度、飞行路径角和航向角。本发明的UD‐SKF滤波算法的原理图参见图3。

本发明共由步骤一、步骤二、步骤三和步骤四四个步骤组成，通过建立IMU和信标无线电组合导航的量测方程，然后利用UD‐SKF滤波算法消除测量信息中的误差的影响，并增强滤波的稳定性，达到高效实时估计探测器导航状态的目的。

本发明的原理是：由于火星进入段气动加热环境恶劣，且探测器包裹在整流罩内，导致很多导航敏感器无法正常工作，只能主要依靠惯性测量单元（简称IMU）提供导航信息，而仅依靠IMU提供的导航信息无法满足完探测器任务精确导航的要求。由于在火星进入段甚高频无线电可以进行通信，因此利用探测器与预先投放的火星在轨信标或火星表面固定信标之间的甚高频无线电双向通信测出探测器与信标之间的距离作为新的测量信息加入量测方程。经可观性理论分析和仿真验证表明，与三个信标进行双向通信测距即可满足完全可观的，且状态估计更加精确。但是由于IMU提供的测量信息带有常值偏差和双向测距存在通信偏差，这些偏差将给导航状态估计带来较大的偏差，甚至可能导致发散。因此，采用UD‐SKF滤波算法将这些偏差的方差融合到对状态的估计中而不去估计这些偏差，同时使用UD分解将来避免增益矩阵可能产生的较多误差，增强滤波过程的稳定性，从而能对探测器进行实时高效地估计导航状态。

其中，在步骤一中所述的“建立火星进入段的动力学方程”，是指令

$f(X\left(t\right),t)=\left[\begin{array}{c}v\sin \mathrm{\γ}\\ \frac{v\cos \mathrm{\γ}\sin \mathrm{\ψ}}{r\cos \mathrm{\λ}}\\ \frac{v\cos \mathrm{\γ}\cos \mathrm{\ψ}}{r}\\ -D-g\left(r\right)\sin \mathrm{\γ}\\ (\frac{v}{r}-\frac{g\left(r\right)}{v})\cos \mathrm{\γ}+\frac{L}{v}\cos \mathrm{\φ}\\ \frac{v}{r}\sin \mathrm{\ψ}\cos \mathrm{\γ}\tan \mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}\sin \mathrm{\φ}\end{array}\right]---\left(41\right)$

得到式（5）。

其中，在步骤二中所述的“建立火星进入段的量测方程”，是指令式（7）中的第一式为

$h(X\left(t\right),t)={[-D,-D\·L/D\·\sin \mathrm{\φ},D\·L/D\·\cos \mathrm{\φ},{\tilde{R}}_1,{\tilde{R}}_2,\·\·\·,{\tilde{R}}_m]}^T---\left(42\right)$

得到式（6）。

其中，在步骤三中所述的“对动力学方程（5）和量测方程（6）进行离散化”，所采用的方法为泰勒级数展开法。泰勒级数在数学上是一个无穷可微的函数f(x)的幂级数展开式：

$f\left(x\right)=\underset{n=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{f^{\left(n\right)}\left(a\right)}{n!}{(x-a)}^n---\left(43\right)$

式中，n!表示n的阶乘,而f⁽ⁿ⁾(a)表示函数f(x)在点x=a处的n阶导数。实际应用中，泰勒级数需要截断，只取有限项，因此会产生相应的截断误差。

本发明的优点是：本发明所采用的UD‐SKF滤波算法与传统的EKF相比，增加了少量的计算量，不仅将量测方程的偏差信息融合到滤波过程中，极大地改善了滤波效果，提高了导航精度，而且利用UD分解，减少了滤波过程中出现较大误差甚至发散的情况，增强滤波过程的稳定性，从而能对探测器进行实时高效地估计导航状态。

附图说明

图1为本发明的流程图。

图2为本发明中的测量信息——无线电通信双向测距示意图。

图3为本发明中的滤波方法——UD‐SKF滤波算法的原理图。

图4为本发明的UD‐SKF滤波算法与EKF导航误差对比图。

图中符号说明如下：

图2中“r_orbiter”表示在轨信标与火星中心的距离向量，“r_lander”表示探测器与火星中心的距离向量，“r_beacon”表示火星表面固定信标与火星中心的距离向量。

图4中“EKFErrors”表示EKF滤波算法的误差，“UD‐SKFErrors”表本发明的UD‐SKF滤波算法的误差，“taxis”表示时间轴，“AltitudeError”表示探测器的高度误差，“VelocityError”表示探测器的速度误差，“FPAError”表示探测器的飞行路径角误差，“LongitudeError”表示探测器的经度误差，“LatitudeError”表示探测器的纬度误差，“AzimuthError”表示探测器的航向角误差。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

本发明涉及一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，其计算流程图如图1所示和UD-SKF算法方法流程图如图2所示，它包括以下四个步骤：

步骤一、建立火星进入段的动力学方程

在火星进入段，将探测器看作是一个无动力的质点。假设进入段的大气相对火星是静止的，同时忽略火星的自转，建立火星进入段的三自由度动力学方程：

$\stackrel{\·}{r}=v\sin \mathrm{\γ}$

$\stackrel{\·}{\mathrm{\θ}}=\frac{v\cos \mathrm{\γ}\sin \mathrm{\ψ}}{r\cos \mathrm{\λ}}$

$\stackrel{\·}{\mathrm{\λ}}=\frac{v\cos \mathrm{\γ}\cos \mathrm{\ψ}}{r}$

$\stackrel{\·}{v}=-D-g\left(r\right)\sin \mathrm{\γ}---\left(1\right)$

$\stackrel{\·}{\mathrm{\γ}}=(\frac{v}{r}-\frac{g\left(r\right)}{v})\cos \mathrm{\γ}+\frac{L}{v}\cos \mathrm{\φ}$

$\stackrel{\·}{\mathrm{\ψ}}=\frac{v}{r}\sin \mathrm{\ψ}\cos \mathrm{\γ}\tan \mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}\sin \mathrm{\φ}$

式中，r,v,θ,λ,γ,ψ分别表示探测器距离火星中心的半径距离（高度）、探测器的速度、经度、纬度、飞行路径角和航向角，火星重力加速度g(r)＝μ/r²，μ是火星重力加速度常数（μ＝42828.29×10⁹m³/s²），D和L分别为阻力加速度和升力加速度：

$D=\frac{C_{D}S}{m}\stackrel{\‾}{q}---\left(2\right)$

$L=\frac{C_{L}S}{m}\stackrel{\‾}{q}---\left(3\right)$

式中C_D和C_L分别为阻力系数和升力系数，S为探测器参考面积，m为探测器质量，为动压，ρ是火星大气密度，此处采用指数大气密度模型：

$\mathrm{\ρ}={\mathrm{\ρ}}_0\exp \left\{\frac{r_0-r}{h_s}\right\}---\left(4\right)$

式中ρ₀是标称密度，r₀为火星径向基准位置(r₀＝7.5km)，h_s为大气定标高度(h_s＝7.5km)。

取状态向量为X=[rθλvγψ]^T，考虑到系统噪声向量w(t)，则动力学方程（1）可简写为

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(5\right)$

式中，f(X(t),t)为系统非线性连续状态转移函数，w(t)为零均值的高斯白噪声。

其中，本发明的仿真验证中的探测器的参考值和估计值如下表所示：

其中，本发明的仿真验证中的弹道系数（B=C_DS/m）、升阻比（L/D=C_L/C_D）和标称大气密度参考值和估计值分别如下表所示：

步骤二、建立火星进入段的量测方程

本发明以基于星间甚高频无线电测量通信的自主导航方案为基础，也即是将探测器与火星在轨信标或火星表面固定信标之间的双向测距信息作为测量信息和IMU提供的测量信息一起组合为新的测量信息，建立火星进入段量测方程：

Z=h(X(t),t)+b+v(t)(6)

式中，

$h(X\left(t\right),t)={\left[\begin{array}{cc}{a}^v& \tilde{R}\end{array}\right]}^T$

a^v=[-D-D·L/D·sinφD·L/D·cosφ]^T(7)

$\tilde{R}={\[{\tilde{R}}_1,{\tilde{R}}_2,\·\·\·,{\tilde{R}}_m\]}^T$

v(t)=[η_aξ_R]^T

式中，a^v为速度系下的加速度，表示双向测距信息，维数m由有效的火星在轨信标和火星表面固定信标数目决定，b为常值偏差向量，v(t)为零均值的高斯白噪声，η_a为IMU测量信息的测量噪声向量，ξ_R为双向测距的测量噪声向量。

其中，IMU的输出是由加速度计测量得到的本体系三个轴向的线加速度，测量模型为

$\tilde{a}=a+b_a+{\mathrm{\η}}_a---\left(8\right)$

式中，表示各个轴的加速度测量值，a表示真实的加速度，b_a加速度计常值漂移，测量噪声η_a为零均值的高斯白噪声。

其中，双向测距指的是探测器与火星在轨信标或火星表面固定信标之间通过双向无线电通信测量得到的二者之间的距离，无线电通信双向测距示意图参见图2，可表示为:

$\tilde{R}=R+b_R+{\mathrm{\ξ}}_R---\left(9\right)$

$R=\sqrt{{(r_l-r_k)}^T(r_l-r_k)}---\left(10\right)$

式中，r_l和r_k分别表示火星质心惯性系下火星探测器和火星在轨信标(k=o)或火星表面固定信标(k=b)，b_R为双向测距的常值偏差，测距噪声ξ_R为零均值的高斯白噪声。

其中，本发明的仿真验证中选择一个火星在轨信标和两个火星表面固定信标，具体参数参见下表：

其中，本发明的仿真验证中对于测量中IMU的加速度常值偏差取为0.003m/s²和通信偏差针对在轨信标和固定信标分别取为24m和45m。

步骤三、对上述动力学方程（5）和量测方程（6）进行离散化

X_k+1=F(X_k)+w_k（11）

式中，k=1,2,3,…，F(X_k)为f(X(t),t)离散后的非线性状态转移函数，为h(X(t),t)离散后的非线性量测函数，w_k和v_k互不相关，且其方差矩阵分别为Q_k和R_k。

然后将方程（11）和（12）进行线性化，也即将式（11）中的非线性离散函数F(X_k)围绕估计值展开成泰勒级数，并略去二阶以上项，得线性化之后的动力学方程

X_k+1=Φ_k+1/kX_k+U_k+w_k（13）

式中，

${\mathrm{\Φ}}_{k+1/k}=\frac{\∂F\left(X_k\right)}{\∂X_k}{|}_{x_k={\hat{x}}_k}---\left(14\right)$

$U_k=F\left({\hat{X}}_k\right)-\frac{\∂F\left(X_k\right)}{\∂X_k}{|}_{x_k={\hat{x}}_k}\·{\hat{X}}_k---\left(15\right)$

然后，再将式（12）中的非线性离散函数围绕估计值和展开成泰勒级数，并略去二阶以上项，得线性化之后的增量量测方程

Z_k=H_kX_k+Y_k+v_k（16）

式中，

通过上述过程，就得到了线性化后的动力学方程和量测方程，式中U_k和Y_k为非随机的外作用项。

步骤四、采用UD-SKF滤波算法，并输出导航信息。

考虑式（6）中的常值偏差向量b未能精确建模，故施密特-卡尔曼滤波算法（SchmidtKalmanfilter,简称SKF）考虑将其方差融入到算法中，也即是考虑了偏差与状态的互协方差，但是不估计这些偏差。同时，在计算增益矩阵时采用UD分解以确保滤波的数值稳定性。

其中步骤四所采用的UD-SKF滤波算法实现步骤为：

7.将常值偏差向量b扩维进入状态向量，则动力学方程（13）和量测方程（16）变为

$\left[\begin{array}{c}{X}_{k+1}\\ b_{k+1}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\Φ}}_{k+1|k}& 0\\ 0& I\end{array}\right]=\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+\left[\begin{array}{c}{w}_k\\ 0\end{array}\right]---\left(19\right)$

$Z_k=[H_k,I]\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+v_k---\left(20\right)$

式中，我们考虑的常值偏差向量满足条件：b_k=b_k-1。其方差矩阵B₀为

B₀=Cov{b₀}=Cov{b_k}(21)

和偏差与状态的互协方差矩阵C_k为

$C_k=E\left\{{\tilde{X}}_k{b}^T\right\}=E\left\{(X_k-{\hat{X}}_k)b^T\right\},---\left(22\right)$

并且初始值为C₀=0。式中，为卡尔曼滤波第k步的状态估计。

同时，相应与动力学方程（19）和量测方程（20）第k步的误差方差阵为

式中，为C_k转置矩阵，P_k为状态X_k的误差方差阵

$P_k=E\left\{{\tilde{X}}_k{\tilde{X}}_k^T\right\}=E\left\{(X_k-{\hat{X}}_k){(X_k-{\hat{X}}_k)}^T\right\}---\left(24\right)$

然后，初始化状态向量X₀及误差方差阵P₀。

8.由第k步的状态估计可得，第k+1步的状态一步预测为

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1|k}{\hat{X}}_k---\left(25\right)$

9.第k+1步的一步预测误差方差矩阵为

则可得到状态和偏差的一步预测误差方差矩阵P_k+1|k和C_k+1|k

$P_{k+1|k}={\mathrm{\Φ}}_{k+1|k}{P}_k{\mathrm{\Φ}}_{k+1|k}^T+Q_k---\left(27\right)$

C_k+1|k=Φ_k+1|kC_k(28)

10.第k+1步的滤波增益矩阵为

其中，由于不需要估计偏差，故强制令偏差项的增益矩阵为零。

则状态的增益矩阵K_k+1为

$K_{k+1}=[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]{\mathrm{\Ω}}_{k+1}^{-1}---\left(30\right)$

${\mathrm{\Ω}}_{k+1}=H_{k+1}{P}_{k+1|k}{H}_{k+1}^T+C_{k+1|k}^T{H}_{k+1}^T+H_{k+1}{C}_{k+1|k}+B_0+R_{k+1}---\left(31\right)$

由于式（31）中的Ω_k+1在滤波过程中无法保持非负定性，此时式（30）中的增益矩阵K_k+1在求逆运算中会产生较大的误差，甚至会发散。为了保持滤波的数值稳定性，本发明采用UD分解来计算增益矩阵K_k+1，会起到很好的效果。UD分解的详细步骤如下：

1）将对称正定矩阵Ω_k分解成一个对角阵D和一个单位上三角阵U的形式：

Ω_k=UDU^T(32)

2）将式（32）代入是式（30）中，得到

$K_k{\mathrm{UDU}}^T=P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}---\left(33\right)$

${\mathrm{UDU}}^T{K}_k^T={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T---\left(34\right)$

3）令则有

${\mathrm{UDU}}^{T}X={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T---\left(35\right)$

4）可以通过以下的迭代过程得到

${\mathrm{UX}}_1={[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}]}^T$

DX₂=X₁(36)

U^TX=X₂

5）最后令然后即可得到K_k。

11.第k+1步的状态估计为

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-H_{k+1}{\hat{X}}_{k+1/k})---\left(37\right)$

12.第k+1步的估计误差方差矩阵为

则状态的估计误差方差矩阵P_k+1为

$P_{k+1}=P_{k+1|k}-K_{k+1}{\mathrm{\Ω}}_{k+1}{K}_{k+1}^T---\left(39\right)$

和偏差与状态的互协方差矩阵为

C_k+1=C_k+1|k-K_k+1(H_k+1C_k+1|k+B₀)(40)

通过以上6步循环进行即可得到探测器的实时状态估计值，包含探测器的高度、速度、经度、纬度、飞行路径角和航向角。本发明的UD‐SKF滤波算法的原理图参见图3。

其中，在步骤一中所述的“建立火星进入段的动力学方程”，是指令

$f(X\left(t\right),t)=\left[\begin{array}{c}v\sin \mathrm{\γ}\\ \frac{v\cos \mathrm{\γ}\sin \mathrm{\ψ}}{r\cos \mathrm{\λ}}\\ \frac{v\cos \mathrm{\γ}\cos \mathrm{\ψ}}{r}\\ -D-g\left(r\right)\sin \mathrm{\γ}\\ (\frac{v}{r}-\frac{g\left(r\right)}{v})\cos \mathrm{\γ}+\frac{L}{v}\cos \mathrm{\φ}\\ \frac{v}{r}\sin \mathrm{\ψ}\cos \mathrm{\γ}\tan \mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}\sin \mathrm{\φ}\end{array}\right]---\left(41\right)$

得到式（5）。

其中，在步骤二中所述的“建立火星进入段的量测方程”，是指令式（7）中的第一式为

$h(X\left(t\right),t)={[-D,-D\·L/D\·\sin \mathrm{\φ},D\·L/D\·\cos \mathrm{\φ},{\tilde{R}}_1,{\tilde{R}}_2,{\tilde{R}}_3]}^T---\left(42\right)$

得到式（6）。由于是具体实施中采用三个信标，故式（42）中只有三个信标的距离信息。

其中，在步骤三中所述的“对动力学方程（5）和量测方程（6）进行离散化”，所采用的方法为泰勒级数展开法。泰勒级数在数学上是一个无穷可微的函数f(x)的幂级数展开式：

$f\left(x\right)=\underset{n=0}{\overset{\∞}{\mathrm{\Σ}}}\frac{f^{\left(n\right)}\left(a\right)}{n!}{(x-a)}^n---\left(43\right)$

式中，n!表示n的阶乘,而f⁽ⁿ⁾(a)表示函数f(x)在点x=a处的n阶导数。实际应用中，泰勒级数需要截断，只取有限项，因此会产生相应的截断误差。

从图4的UD‐SKF滤波算法与EKF导航误差对比图可以看出，由EKF滤波算法所估计的状态仅有航向角是收敛的，其余的状态都是发散的，而由UD‐SKF滤波算法所估计的状态除了速度的收敛速度较慢之外，其余的状态都很快地收敛到参考值。从二者的对比来看，明显UD‐SKF滤波算法要远远好于EKF滤波算法，且最终的滤波效果也好于EKF滤波算法（见下表）。

Claims (4)

1.一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，其特征在于：它包括以下四个步骤：

步骤一、建立火星进入段的动力学方程

在火星进入段，将探测器看作是一个无动力的质点，假设进入段的大气相对火星是静止的，同时忽略火星的自转，建立火星进入段的三自由度动力学方程：

$\stackrel{\·}{r}=vsin\mathrm{\γ}$

$\stackrel{\·}{\mathrm{\θ}}=\frac{vcos\mathrm{\γ}sin\mathrm{\ψ}}{rcos\mathrm{\λ}}$

$\begin{array}{c}\stackrel{\·}{\mathrm{\λ}}=\frac{vcos\mathrm{\γ}cos\mathrm{\ψ}}{r}\\ \stackrel{\·}{v}=-D-g\left(r\right)\sin \mathrm{\γ}\end{array}---\left(1\right)$

$\stackrel{\·}{\mathrm{\γ}}=\left(\frac{v}{r}-\frac{g\left(r\right)}{v}\right)cos\mathrm{\γ}+\frac{L}{v}cos\mathrm{\φ}$

$\stackrel{\·}{\mathrm{\ψ}}=\frac{v}{r}sin\mathrm{\ψ}cos\mathrm{\γ}tan\mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}sin\mathrm{\φ}$

式中，r,v,θ,λ,γ,ψ分别表示探测器距离火星中心的半径距离即高度、探测器的速度、经度、纬度、飞行路径角和航向角，火星重力加速度g(r)＝μ/r²，μ是火星重力加速度常数，D和L分别为阻力加速度和升力加速度：

$D=\frac{C_{D}S}{m}\stackrel{\‾}{q}---\left(2\right)$

$L=\frac{C_{L}S}{m}\stackrel{\‾}{q}---\left(3\right)$

式中C_D和C_L分别为阻力系数和升力系数，S为探测器参考面积，m为探测器质量，为动压，ρ是火星大气密度，此处采用指数大气密度模型：

$\mathrm{\ρ}={\mathrm{\ρ}}_0\exp \left\{\frac{r_0-r}{h_s}\right\}---\left(4\right)$

式中ρ₀是标称密度，r₀为火星径向基准位置，h_s为大气定标高度；

取状态向量为X＝[rθλvγψ]^T，考虑到系统噪声向量w(t)，则动力学程(1)简写为

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(5\right)$

式中，f(X(t),t)为系统非线性连续状态转移函数，w(t)为零均值的高斯白噪声；

步骤二、建立火星进入段的量测方程

以基于星间甚高频无线电测量通信的自主导航方案为基础，也即是将探测器与火星在轨信标或火星表面固定信标之间的双向测距信息作为测量信息和IMU提供的测量信息一起组合为新的测量信息，建立火星进入段量测方程：

Z＝h(X(t),t)+b+v(t)(6)

式中，

$h(X\left(t\right),t)={\left[\begin{array}{cc}{a}^v& \tilde{R}\end{array}\right]}^T$

$a^v={\left[\begin{array}{ccc}-D& -D\·L/D\·\sin \mathrm{\φ}& D\·L/D\·\cos \mathrm{\φ}\end{array}\right]}^T---\left(7\right)$

$\tilde{R}={\left[\begin{array}{cccc}{\tilde{R}}_1,& {\tilde{R}}_2,& \mathrm{...},& {\tilde{R}}_m\end{array}\right]}^T$

$v\left(t\right)={\left[\begin{array}{cc}{\mathrm{\η}}_a& {\mathrm{\ξ}}_R\end{array}\right]}^T$

式中，a^v为速度系下的加速度，表示双向测距信息，维数m由有效的火星在轨信标和火星表面固定信标数目决定，b为常值偏差向量，v(t)为零均值的高斯白噪声，η_a为IMU测量信息的测量噪声向量，ξ_R为双向测距的测量噪声向量；

其中，IMU的输出是由加速度计测量得到的本体系三个轴向的线加速度，测量模型为

$\tilde{a}=a+b_a+{\mathrm{\η}}_a---\left(8\right)$

式中，表示各个轴的加速度测量值，a表示真实的加速度，b_a加速度计常值漂移，测量噪声η_a为零均值的高斯白噪声；

其中，双向测距指的是探测器与火星在轨信标或火星表面固定信标之间通过双向无线电通信测量得到的二者之间的距离，表示为:

$\tilde{R}=R+b_R+{\mathrm{\ξ}}_R---\left(9\right)$

$R=\sqrt{{(r_l-r_k)}^T(r_l-r_k)}---\left(10\right)$

式中，r_l和r_k分别表示火星质心惯性系下火星探测器和火星在轨信标或火星表面固定信标，b_R为双向测距的常值偏差，测距噪声ξ_R为零均值的高斯白噪声；

步骤三、对上述动力学方程(5)和量测方程(6)进行离散化

X_k+1＝F(X_k)+w_k(11)

式中，k＝1,2,3,…，F(X_k)为f(X(t),t)离散后的非线性状态转移函数，为h(X(t),t)离散后的非线性量测函数，w_k和v_k互不相关，且其方差矩阵分别为Q_k和R_k；

然后将方程(11)和(12)进行线性化，也即将式(11)中的非线性离散函数F(X_k)围绕估计值展开成泰勒级数，并略去二阶以上项，得线性化之后的动力学方程

X_k+1＝Φ_k+1/kX_k+U_k+w_k(13)

式中，

${\mathrm{\Φ}}_{k+1/k}=\frac{\∂F\left(X_k\right)}{\∂X_k}{|}_{X_k={\hat{X}}_k}---\left(14\right)$

$U_k=F\left({\hat{X}}_k\right)-\frac{\∂F\left(X_k\right)}{\∂X_k}|{}_{X_k={\hat{X}}_k}\·{\hat{X}}_k---\left(15\right)$

然后，再将式(12)中的非线性离散函数围绕估计值和展开成泰勒级数，并略去二阶以上项，得线性化之后的增量量测方程：

Z_k＝H_kX_k+Y_k+v_k(16)

式中，

通过上述过程，就得到了线性化后的动力学方程和量测方程，式中U_k和Y_k为非随机的外作用项；

步骤四、采用UD-SKF滤波算法，并输出导航信息

考虑式(6)中的常值偏差向量b未能精确建模，故施密特-卡尔曼滤波算法即SKF考虑将其方差融入到算法中，也即是考虑了偏差与状态的互协方差，但是不估计这些偏差，同时，在计算增益矩阵时采用UD分解以确保滤波的数值稳定性；其中步骤四所采用的UD-SKF滤波算法实现步骤为：

步骤4.1.将常值偏差向量b扩维进入状态向量，则动力学方程(13)和量测方程(16)变为

$\left[\begin{array}{c}{X}_{k+1}\\ b_{k+1}\end{array}\right]=\left[\begin{array}{cc}{\mathrm{\Φ}}_{k+1|k}& 0\\ 0& I\end{array}\right]\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+\left[\begin{array}{c}{w}_k\\ 0\end{array}\right]---\left(19\right)$

$Z_k=\[H_k,I\]\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+v_k---\left(20\right)$

式中，我们考虑的常值偏差向量满足条件：b_k＝b_k-1。其方差矩阵B₀为

B₀＝Cov{b₀}＝Cov{b_k}(21)

和偏差与状态的互协方差矩阵C_k为

$C_k=E\left\{{\tilde{X}}_k{b}^T\right\}=E\left\{\left(X_k-{\hat{X}}_k\right)b^T\right\},---\left(22\right)$

并且初始值为C₀＝0，式中，为卡尔曼滤波第k步的状态估计；

同时，相应与动力学方程(19)和量测方程(20)第k步的误差方差阵为

式中，为C_k转置矩阵，P_k为状态X_k的误差方差阵

$P_k=E\left\{{\tilde{X}}_k{\tilde{X}}_k^T\right\}=E\left\{\left(X_k-{\hat{X}}_k\right){\left(X_k-{\hat{X}}_k\right)}^T\right\}---\left(24\right)$

然后，初始化状态向量X₀及误差方差阵P₀；

步骤4.2.由第k步的状态估计可得，第k+1步的状态一步预测为

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1|k}{\hat{X}}_k---\left(25\right)$

步骤4.3.第k+1步的一步预测误差方差矩阵为

则得到状态和偏差的一步预测误差方差矩阵P_k+1|k和C_k+1|k

$P_{k+1|k}={\mathrm{\Φ}}_{k+1|k}{P}_k{\mathrm{\Φ}}_{k+1|k}^T+Q_k---\left(27\right)$

C_k+1|k＝Φ_k+1|kC_k(28)

步骤4.4.第k+1步的滤波增益矩阵为

其中，由于不需要估计偏差，故强制令偏差项的增益矩阵为零；

则状态的增益矩阵K_k+1为

$K_{k+1}=\[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}\]{\mathrm{\Ω}}_{k+1}^{-1}---\left(30\right)$

${\mathrm{\Ω}}_{k+1}=H_{k+1}{P}_{k+1|k}{H}_{k+1}^T+C_{k+1|k}^T{H}_{k+1}^T+H_{k+1}{C}_{k+1|k}+B_0+R_{k+1}---\left(31\right)$

由于式(31)中的Ω_k+1在滤波过程中无法保持非负定性，此时式(30)中的增益矩阵K_k+1在求逆运算中会产生较大的误差，甚至会发散；为了保持滤波的数值稳定性，采用UD分解来计算增益矩阵K_k+1，会起到很好的效果；UD分解的详细步骤如下：

1)将对称正定矩阵Ω_k分解成一个对角阵D和一个单位上三角阵U的形式：

Ω_k＝UDU^T(32)

2)将式(32)代入是式(30)中，得到

$K_k{\mathrm{UDU}}^T=P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}---\left(33\right)$

${\mathrm{UDU}}^T{K}_k^T={\[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}\]}^T---\left(34\right)$

3)令则有

${\mathrm{UDU}}^{T}X={\[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}\]}^T---\left(35\right)$

4)通过以下的迭代过程得到

${\mathrm{UX}}_1={\[P_{k+1|k}{H}_{k+1}^T+C_{k+1|k}\]}^T$

DX₂＝X₁(36)

U^TX＝X₂

5)最后令然后即得到K_k；

步骤4.5.第k+1步的状态估计为

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-H_{k+1}{\hat{X}}_{k+1/k})---\left(37\right)$

步骤4.6.第k+1步的估计误差方差矩阵为

则状态的估计误差方差矩阵P_k+1为

$P_{k+1}=P_{k+1|k}-K_{k+1}{\mathrm{\Ω}}_{k+1}{K}_{k+1}^T---\left(39\right)$

和偏差与状态的互协方差矩阵为

C_k+1＝C_k+|k-K_k+1(H_k+1C_k+1|k+B₀)(40)

通过步骤4.1至步骤4.6循环进行即得到探测器的实时状态估计值，包含探测器的高度、速度、经度、纬度、飞行路径角和航向角；

通过上述四个步骤建立IMU和信标无线电组合导航的量测方程，然后利用UD-SKF滤波算法消除测量信息中的误差的影响，并增强滤波的稳定性，达到高效实时估计探测器导航状态的目的。

2.根据权利要求1所述的一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，其特征在于：在步骤一中所述的“建立火星进入段的动力学方程”，是指令

$f(X\left(t\right),t)=\left[\begin{array}{c}vsin\mathrm{\γ}\\ \frac{vcos\mathrm{\γ}\sin \mathrm{\ψ}}{rcos\mathrm{\λ}}\\ \frac{vcos\mathrm{\γ}\cos \mathrm{\ψ}}{r}\\ -D-g\left(r\right)\sin \mathrm{\γ}\\ \left(\frac{v}{r}-\frac{g\left(r\right)}{v}\right)\cos \mathrm{\γ}+\frac{L}{v}cos\mathrm{\φ}\\ \frac{v}{r}sin\mathrm{\ψ}\cos \mathrm{\γ}\tan \mathrm{\λ}+\frac{L}{v\cos \mathrm{\γ}}sin\mathrm{\φ}\end{array}\right]---\left(41\right)$

得到式(5)。

3.根据权利要求1所述的一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，其特征在于：在步骤二中所述的“建立火星进入段的量测方程”，是指令式(7)中的第一式为

$h\left(X\right(t),t)={\left[\begin{array}{cccc}-D,-D\·L/D\·\sin \mathrm{\φ},D\·L/D\·\cos \mathrm{\φ},{\tilde{R}}_1,& {\tilde{R}}_2,& \mathrm{...},& {\tilde{R}}_m\end{array}\right]}^T---\left(42\right)$

得到式(6)。

4.根据权利要求1所述的一种火星进入段IMU和甚高频无线电组合导航的UD-SKF方法，其特征在于：在步骤三中所述的“对动力学方程(5)和量测方程(6)进行离散化”，所采用的方法为泰勒级数展开法；泰勒级数在数学上是一个无穷可微的函数f(x)的幂级数展开式：

$f\left(x\right)=\underset{n=0}{\overset{\mathrm{\∞}}{\Σ}}\frac{f^{\left(n\right)}\left(a\right)}{n!}{(x-a)}^n---\left(43\right)$

式中，n！表示n的阶乘,而f⁽ⁿ⁾(a)表示函数f(x)在点x＝a处的n阶导数，实际应用中，泰勒级数需要截断，只取有限项，因此会产生相应的截断误差。