火星动力下降段多源信息组合导航的迭代SKF方法

摘要

一种火星动力下降段多源信息组合导航的迭代SKF方法，步骤如下：一、利用火星动力下降段的动力学方程;二、建立火星动力下降段的量测方程；三、离散化动力学方程和量测方程并线性化之后得动力学方程和量测方程；四、利用迭代SKF滤波算法输出导航信息。通过以上四个步骤，构造动力学方程和量测方程，然后利用迭代SKF滤波算法消除测量信息中误差的影响，确保滤波算法的稳定性，达到高效实时估计探测器导航状态的目的。它有效地修正了量测方程中偏差对滤波的影响，而且利用迭代方法，修正了由于泰勒级数带来的截断误差带来的滤波误差，提高了导航精度，增强了滤波过程的稳定性，从而能对探测器进行实时高效地估计导航状态。

主权项

1.一种火星动力下降段多源信息组合导航的迭代SKF方法，其特征在于：该方法步骤如下：步骤一、利用火星动力下降段的动力学方程火星动力下降段的情况比较复杂，建立精确的动力学模型非常困难，在考虑IMU输出的基础之上，利用其构造动力下降段的动力学方程：$\stackrel{·}{r}=v$$\stackrel{·}{v}=C_b^i(\tilde{a}-b_a)+C_g^{i}g-C_b^i{\eta }_a---\left(1\right)$$\stackrel{·}{\Omega }=K(\tilde{\omega }-b_{\omega })-K{\eta }_{\omega }$式中，r=[x y z]^T表示着陆点固连系下的位置向量，v=[v_x v_y v_z]^T表示着陆点固连系下的速度向量，Ω=[σ θ ψ]^T表示着陆点固连系下的姿态角，控制量σ为探测器的滚转角，θ为探测器的经度，ψ为探测器的航向角；表示火星中心惯性系到探测器本体坐标系的转换矩阵，表示火星地理坐标系到火星中心惯性系的转换矩阵；表示IMU中由加速度计输出的本体坐标系下的三个轴向的线加速度，b_a表示加速度计的常值漂移，η_a表示加速度计的噪声；表示IMU中由陀螺输出的本体坐标系下的三个轴向的瞬时旋转角速度，b_ω表示陀螺的常值漂移，η_ω表示陀螺的噪声；g表示地理坐标系的火星重力加速度；K表示姿态运动学矩阵$K=\frac{1}{\cos \theta }\left[\begin{array}{ccc}\cos \theta & \sin \theta \sin \sigma & \sin \theta \cos \sigma \\ 0& \cos \theta \cos \sigma & -\cos \theta \sin \sigma \\ 0& \sin \sigma & \cos \sigma \end{array}\right]---\left(2\right)$取状态向量为X=[r^T v^T Ω^T]^T，则动力下降段的动力学方程（1）简化为$\stackrel{·}{X}=f\left(X\right)+w---\left(3\right)$式中，f(X)为系统非线性连续状态转移函数，系统噪声$w={\left[\begin{array}{ccc}{0}_{3\times 1}& -C_b^i{\eta }_a^T& -{\mathrm{K\eta }}_{\omega }^T\end{array}\right]}^T$为零均值的高斯白噪声；步骤二、火星动力下降段的量测方程固定在本体坐标系上的微型传感器MCAV输出的是探测器的高度和三个轴的速率，则以其测量信息为基础建立火星动力下降段的量测方程：Z=h(X)+b+v_m  (4)式中，h(X)=[r_z v^b]^T，  (5)r_z为探测器的Z轴坐标，即探测器距离火星表面的高度，v^b为本体坐标系下的速度，其表达式为$v^b=C_i^{b}v---\left(6\right)$式中，表示探测器本体坐标系到火星中心惯性系的转换矩阵，b为常值偏差向量，测量噪声v_m为零均值的高斯白噪声，且与系统噪声w不相关；步骤三、离散化上述动力学方程（3）和量测方程（4），得X_k+1=F(X_k)+w_k  （7）式中，k=1,2,3,…，F(X_k)为f(X)离散后的非线性状态转移函数，为h(X)离散后的非线性量测函数，w_k和v_k互不相关，且其方差矩阵分别为Q_k和R_k；将式（7）中的非线性离散函数F(X_k)围绕估计值展开成泰勒级数，并略去二阶以上项，得线性化之后的动力学方程$X_{k+1}={\Phi }_{k+1|k}{X}_k+U_k+w_k---\left(9\right)$式中，${\Phi }_{k+1|k}=\frac{\partial F\left(X_k\right)}{{\partial X}_k}{|}_{X_k={\hat{X}}_k}---\left(10\right)$$U_k=F\left({\hat{X}}_k\right)-\frac{\partial F\left(X_k\right)}{{\partial F}_k}{|}_{X_k={\hat{X}}_k}·{\hat{X}}_k---\left(11\right)$然后，再将式（8）中的非线性离散函数围绕估计值和展开成泰勒级数，并略去二阶以上项，得线性化之后的量测方程Z_k=H_kX_k+Y_k+v_k  （12）式中，通过上述过程，就得到了线性化后的动力学方程和量测方程，式中U_k和Y_k为非随机的外作用项；步骤四、迭代SKF滤波算法及导航信息输出由于式（4）中的常值偏差b未能精确知道，故施密特-卡尔曼滤波算法即SKF在不估计这些偏差的基础之上，考虑将其方差融入到滤波算法中，也即是通过考虑偏差与状态的互协方差，增加滤波的精度；同时，由于面对的是非线性滤波，对量测方程进行泰勒级数展开时会产生截断误差，因此采用迭代的方法来减少对非线性量测方程进行线性化产生的误差对滤波过程的影响，从而达到减少发散现象，提高滤波精度，确保滤波的数值稳定性；所述的迭代SKF滤波算法实现步骤如下：1.扩维状态向量X_k，即将常值偏差向量b加入，则动力学方程（9）和量测方程（12）变为$\left[\begin{array}{c}{X}_{k+1}\\ b_{k+1}\end{array}\right]=\left[\begin{array}{cc}{\Phi }_{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(15\right)$$Z_k=[H_k,I]\left[\begin{array}{c}{X}_k\\ b_k\end{array}\right]+v_k---\left(16\right)$式中，常值偏差向量满足条件：b_k+1=b_k，且其方差矩阵B₀满足B₀=Cov{b₀}=Cov{b_k}  (17)和偏差与状态的互协方差矩阵C_k满足$C_k=E\left\{{\tilde{X}}_k{b}_k^T\right\}=E\left\{(X_k-{\hat{X}}_k)b_k^T\right\},---\left(18\right)$并且初始值为C₀=0；式中，为卡尔曼滤波第k步的状态估计；相应的与动力学方程（15）和量测方程（16）第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(20\right)$开始滤波计算时,需要初始化状态向量和误差方差阵，并设初始化状态向量为X₀及误差方差阵为P₀；2.时间更新过程由第k步的状态估计能得，第k+1步的状态一步预测为${\hat{X}}_{k+1|k}={\Phi }_{k+1|k}{\hat{X}}_k,---\left(21\right)$且第k+1步的一步预测误差方差矩阵为则得到状态和偏差的一步预测误差方差矩阵P_k+1|k和C_k+1|k$P_{k+1|k}={\Phi }_{k+1|k}{P}_k{\Phi }_{k+1|k}^T+Q_k---\left(23\right)$$C_{k+1|k}={\Phi }_{k+1|k}{C}_k---\left(24\right)$3.量测更新过程使用迭代方法进行量测更新：当i=1,2,3,…时，按如下步骤进行循环计算：1)计算第i步的滤波增益矩阵为其中，由于不需要估计偏差，故强制令偏差项的增益矩阵为零；则状态的增益矩阵为$K_{k+1}^i=[P_{k+1|k}{\left(H_{k+1}^i\right)}^T+C_{k+1|k}]{\left({\Omega }_{k+1}^i\right)}^{-1}---\left(26\right)$${\Omega }_{k+1}^i=H_{k+1}^i{P}_{k+1|k}{\left(H_{k+1}^i\right)}^T+C_{k+1|k}^T{\left(H_{k+1}^i\right)}^T+H_{k+1}^i{C}_{k+1|k}+B_0+R_{k+1}---\left(27\right)$2)计算第i步量测信息残差3)计算第i步的状态估计${\hat{X}}_{k+1}^i={\hat{X}}_{k+1|k}+K_{k+1}^i\{{\tilde{Z}}_{k+1}^i-H_{k+1}^i{\hat{X}}_{k+1|k}\}---\left(29\right)$式中，使用迭代方法进行量测更新过程中，当所得状态向量的估计值满足向量的2范数满足的条件—其中阈值设为ε_limit:${\left|\right|{\hat{X}}_{k+1}^i-{\hat{X}}_{k+1}^{i-1}\left|\right|}_2<{ϵ}_{\mathrm{limit}}---\left(31\right)$时即结束循环计算；4.利用量测更新中的状态估计值及相应的参数，计算第k+1步的估计误差方差矩阵为则状态的估计误差方差矩阵P_k+1为$P_{k+1}=P_{k+1|k}-K_{k+1}^i{\Omega }_{k+1}^i{\left(K_{k+1}^i\right)}^T---\left(33\right)$和偏差与状态的互协方差矩阵C_k+1为$C_{k+1}=C_{k+1|k}-K_{k+1}^i(H_{k+1}^i{C}_{k+1|k}+B_0)---\left(34\right)$通过以上4步循环进行即得到火星动力下降段探测器的实时状态估计值，包含探测器的位置向量、速度向量和姿态角；通过以上四个步骤，通过构造动力学方程和建立多源信息组合导航的量测方程，然后利用迭代SKF滤波算法消除测量信息中误差的影响，并确保滤波算法的稳定性，达到高效实时估计探测器导航状态的目的。