CN101949703B - 一种捷联惯性/卫星组合导航滤波方法 - Google Patents

Abstract

一种捷联惯性/卫星组合导航滤波方法。首先，建立包含模型误差的捷联惯性导航系统(SINS)/全球导航卫星系统(GNSS)组合导航滤波非线性数学模型；其次，对传统的二阶插值滤波进行改进，并将改进后的二阶插值滤波与预测滤波相结合，利用预测滤波在线估计出的模型误差修正改进二阶插值滤波器的状态，实现对真实误差状态量——SINS/GNSS组合导航的姿态误差、速度误差及位置误差的精确估计；最后利用估计出的误差状态量计算出更加准确的SINS/GNSS组合导航的位置、速度和姿态。

Description

一种捷联惯性/卫星组合导航滤波方法

技术领域

本发明涉及一种捷联惯性导航系统(Strapdown Inertial Navigation System，SINS)与全球导航卫星系统(Global Navigation Satellite System，GNSS)组合导航的滤波方法，可用于航空、航海、陆地等领域的SINS/GNSS组合导航系统。

背景技术

捷联惯性导航系统(SINS)的基本原理是根据牛顿提出的相对惯性空间的力学定律，利用陀螺仪、加速度计测量载体相对惯性空间的线运动和角运动参数，在给定的运动初始条件下，由计算机进行积分运算，连续、实时地提供位置、速度和姿态信息。SINS完全依靠自身的惯性敏感元件，不依赖任何外界信息测量导航参数，因此，具有隐蔽性好、不受气候条件限制、无信号丢失以及不受干扰等优点，是一种完全自主式、全天候的导航系统。但是，SINS也有其自身的不足。由于捷联解算中的积分原理，惯性器件的误差会导致导航误差随时间积累，因而纯惯性导航系统难以满足远程、长时间运动载体的高精度导航要求。全球导航卫星系统(GNSS)是在无线电技术的基础上，伴随航天技术发展而形成的一种天基无线电导航系统。它的优点是定位精度高，导航误差不随时间积累，可全天时、全天候工作。但是，GNSS难以直接提供姿态信息，并存在数据更新率低、易受电磁干扰等缺点。若将SINS与GNSS组合起来，能够实现二者的优势互补，显著提高导航系统的综合性能。目前，SINS/GNSS组合导航系统已广泛应用于航空、航海等领域，是一种较为理想的组合导航系统。

SINS/GNSS组合导航系统通常采用滤波技术对SINS和GNSS的数据进行融合，获得优于单一子系统的位置、速度和姿态精度。实现SINS/GNSS的高精度组合滤波面临的两大问题，一是SINS/GNSS组合系统为非线性系统，二是系统噪声和量测噪声为非高斯噪声。因此，滤波器的设计以及滤波方法的选择对SINS/GNSS组合导航系统的精度起着至关重要的作用。

在非线性滤波方面，传统非线性滤波方法扩展卡尔曼滤波(Extended KalmanFiltering，EKF)存在模型线性化截断误差问题，且假设噪声为高斯白噪声；迭代扩展卡尔曼滤波(Iterated Extended Kalman Filtering，IEKF)和Unscented卡尔曼滤波(Unscented Kalman Filtering，UKF)是EKF的改进方法，其中UKF无需对非线性模型进行线性化，IEKF通过多次迭代运算能够减小模型线性化截断误差从而提高了估计精度，但这两种方法同样假设噪声为高斯白噪声，且滤波算法耗时较长，难以在实时导航中应用。二阶插值非线性滤波(Second-order Divided Difference Filtering，DD2)采用多维Stirling插值方法来替代EKF非线性函数导数的计算，并获得优于EKF的估计精度，但仍无法克服插值近似引起的模型误差对滤波估计精度的不利影响。因此提高SINS/GNSS组合导航系统的滤波精度和实时性有着十分重要的意义。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种可以达到较高适用精度要求的基于改进二阶插值滤波与预测滤波相结合的捷联惯性/卫星组合导航滤波方法。

本发明的技术解决方案为：一种捷联惯性/卫星组合导航滤波方法，具体步骤如下：

(1)建立包含模型误差的捷联惯性导航系统/全球导航卫星系统组合导航滤波非线性数学模型；

(2)对传统二阶插值滤波进行改进，并基于步骤(1)建立的数学模型对状态进行估计，然后利用预测滤波估计出的模型误差修正改进二阶插值滤波的状态估计，得到修正后的姿态误差、位置误差和速度误差；

(3)利用以上步骤(2)得到的修正后的姿态误差、位置误差和速度误差对捷联惯性导航系统捷联解算出的姿态、位置和速度进行补偿，获得更加准确的导航信息，即补偿后的姿态、位置和速度信息；

(4)将以上步骤(3)中补偿后的位置、速度和姿态作为下一导航时刻的初始值，不断重复以上步骤(2)和步骤(3)，直至捷联惯性导航系统/全球导航卫星系统组合导航结束。

本发明与现有技术相比的优点在于：本发明克服了现有SINS/GNSS组合滤波估计方法的不足，对传统二阶插值滤波进行改进，利用精度优于

的平滑值

替换

计算获得更高精度的状态估计

然后将基于最小模型误差准则的预测滤波思想引入到改进二阶插值滤波估计中，建立包含模型误差的SINS/GNSS组合导航滤波的非线性数学模型，利用预测滤波在线估计出未知模型误差并以此修正改进二阶插值滤波的状态变量，提高状态估计的精度，进而提高SINS捷联解算的位置、速度和姿态精度。该方法克服了传统二阶插值滤波插值近似引起的模型误差对滤波估计精度的不利影响，同时可以在线实时进行，从而能够实现实时、准确的SINS/GNSS组合导航，为载体提供更高精度的位置、速度和姿态信息。

附图说明

图1为本发明采用改进二阶插值滤波与预测滤波的SINS/GNSS组合导航流程图；

图2为导航坐标系与载体坐标系之间的关系示意图，

θ、γ分别为SINS的航向角、俯仰角和横滚角。图中Ox_ny_nz_n为导航坐标系，取为东北天地理坐标系，Ox_by_bz_b为载体坐标系。其中，图2a表示从导航坐标系Ox_ny_nz_n绕z_n轴逆时针旋转

至Ox_n1y_n1z_n；图2b表示Ox_n1y_n1z_n绕x_n1轴逆时针旋转θ至Ox_n1y_n2z_n1，即Ox_n1y_bz_n1；图2c表示Ox_n1y_n2z_n1绕y_b(y_n2)轴逆时针旋转γ至Ox_n2y_n2z_n2，即Ox_by_bz_b；通过以上三次旋转，可实现导航坐标系Ox_ny_nz_n到载体坐标系Ox_by_bz_b的转换；

图3为导航坐标系Ox_ny_nz_n与计算导航坐标系Ox_n′y_n′z_n′之间的关系示意图。

具体实施方式

如图1所示，本发明的具体实施方法如下：

1、建立包含模型误差的SINS/GNSS组合导航滤波非线性数学模型

针对SINS/GNSS组合导航系统具有非线性的特点，建立SINS的非线性误差方程。在此基础上，建立包含模型误差的SINS/GNSS组合导航滤波非线性数学模型。

①SINS非线性误差方程

采用东北天地理坐标系作为导航坐标系。导航坐标系Ox_ny_nz_n和载体坐标系Ox_by_bz_b的定义如图2a所示。SINS的非线性误差方程如下：

a)姿态误差方程

SINS的姿态误差角指SINS数学平台坐标系与导航坐标系之间的失准角。在SINS的姿态计算过程中，四元数法因计算简单且精度高而被广泛采用。下面给出采用加性四元数误差来描述失准角的SINS姿态误差方程。

$\mathrm{\δ}\stackrel{\·}{Q}=\frac{1}{2}\mathrm{M\δQ}+\frac{1}{2}(N\left({\hat{Q}}_b^n\right)\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ib}}^b-Y\left({\hat{Q}}_b^n\right)\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n)---\left(1\right)$

式中为姿态误差；

和

分别为载体坐标系到导航坐标系的真实四元数和计算四元数；

为陀螺的理论输出，ω_x、ω_y、ω_z分别代表载体坐标系中x、y、z轴上陀螺的理论输出，

表示陀螺的测量误差；

为真实导航坐标系相对惯性坐标系的旋转角速度在导航坐标系的投影，ω_E、ω_N、ω_U分别为

在真实导航坐标系中东向、北向、天向轴上的投影，

为的计算误差；M是由ω_x、ω_y、ω_z、ω_E、ω_N、ω_U组成的矩阵，

和

是由

组成的矩阵，M，

和

的具体定义如下

$M\≡\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_x& -{\mathrm{\ω}}_y& -{\mathrm{\ω}}_z\\ {\mathrm{\ω}}_x& 0& {\mathrm{\ω}}_z& -{\mathrm{\ω}}_y\\ {\mathrm{\ω}}_y& -{\mathrm{\ω}}_z& 0& {\mathrm{\ω}}_x\\ {\mathrm{\ω}}_z& {\mathrm{\ω}}_y& -{\mathrm{\ω}}_x& 0\end{array}\right]-\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_E& -{\mathrm{\ω}}_N& -{\mathrm{\ω}}_U\\ {\mathrm{\ω}}_E& 0& -{\mathrm{\ω}}_U& {\mathrm{\ω}}_N\\ {\mathrm{\ω}}_N& {\mathrm{\ω}}_U& 0& -{\mathrm{\ω}}_E\\ {\mathrm{\ω}}_U& -{\mathrm{\ω}}_N& {\mathrm{\ω}}_E& 0\end{array}\right],$ $N\left({\hat{Q}}_b^n\right)\≡\left[\begin{array}{ccc}-{\hat{q}}_1& -{\hat{q}}_2& -{\hat{q}}_3\\ {\hat{q}}_0& -{\hat{q}}_3& {\hat{q}}_2\\ {\hat{q}}_3& {\hat{q}}_0& -{\hat{q}}_1\\ -{\hat{q}}_2& {\hat{q}}_1& {\hat{q}}_0\end{array}\right],$

$Y\left({\hat{Q}}_b^n\right)\≡\left[\begin{array}{ccc}-{\hat{q}}_1& -{\hat{q}}_2& -{\hat{q}}_3\\ {\hat{q}}_0& {\hat{q}}_3& -{\hat{q}}_2\\ -{\hat{q}}_3& {\hat{q}}_0& {\hat{q}}_1\\ {\hat{q}}_2& -{\hat{q}}_1& {\hat{q}}_0\end{array}\right].$

b)速度误差方程

$\mathrm{\δ}\stackrel{\·}{V}=-2\left[{\hat{C}}_b^n{\hat{f}}^b\right]Y^T\left({\hat{Q}}_b^n\right)\mathrm{\δQ}+2{\hat{C}}_b^n{\hat{f}}^b{\left({\hat{Q}}_b^n\right)}^T\mathrm{\δQ}-Y^T\left(\mathrm{\δQ}\right)N\left(\mathrm{\δQ}\right){\hat{f}}^b$ (2)

$-({2\mathrm{\δ\ω}}_{\mathrm{ie}}^n+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{en}}^n)V^n-({2\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\mathrm{\δV}$

式中Vⁿ＝[V_E V_N V_U]^T为SINS解算出的载体相对地球的速度在导航坐标系中的投影，V_E、V_n、V_U分别代表东向速度、北向速度、天向速度；δV＝[δV_E δV_n δV_U]^T为速度误差，δV_E、δV_N、δV_U分别代表东向速度误差、北向速度误差、天向速度误差；

为载体坐标系到导航坐标系的方向余弦矩阵，也称为姿态矩阵；

为

的计算值，

为

的计算误差；

为加速度计的测量值，f^b为载体真实加速度，为加速度计的测量误差，

分别载体坐标系x、y、z轴上加速度计的测量误差；为地球坐标系相对惯性坐标系的旋转角速度在导航坐标系的投影；

为

的计算误差；

为导航坐标系相对地球坐标系的旋转角速度在导航坐标系的投影；

为

的计算误差；

是δQ的非线性函数，具体表达式为

$Y^T\left(\mathrm{\δQ}\right)N\left(\mathrm{\δQ}\right){\hat{f}}^b=$

$\left[\begin{array}{ccc}{{\mathrm{\δq}}_0}^2+{\mathrm{\δ}{q}_1}^2-{\mathrm{\δ}{q}_2}^2-{{\mathrm{\δq}}_3}^2& 2({\mathrm{\δq}}_1\mathrm{\δ}{q}_2-{\mathrm{\δq}}_0\mathrm{\δ}{q}_3)& 2({\mathrm{\δq}}_1\mathrm{\δ}{q}_3+{\mathrm{\δq}}_0\mathrm{\δ}{q}_2)\\ 2({\mathrm{\δq}}_0\mathrm{\δ}{q}_3+{\mathrm{\δq}}_1\mathrm{\δ}{q}_2)& {{\mathrm{\δq}}_0}^2-{{\mathrm{\δq}}_1}^2+{{\mathrm{\δq}}_2}^2-{{\mathrm{\δq}}_3}^2& 2({\mathrm{\δq}}_2{\mathrm{\δq}}_3-{\mathrm{\δq}}_0{\mathrm{\δq}}_1)\\ 2({\mathrm{\δq}}_1{\mathrm{\δq}}_3-\mathrm{\δ}{q}_0\mathrm{\δ}{q}_2)& 2({\mathrm{\δq}}_0{\mathrm{\δq}}_1+{\mathrm{\δq}}_2\mathrm{\δ}{q}_3)& {{\mathrm{\δq}}_0}^2-{{\mathrm{\δq}}_1}^2-{{\mathrm{\δq}}_2}^2+{{\mathrm{\δq}}_3}^2\end{array}\right]{\hat{f}}^b.$

c)位置误差方程

位置误差方程的矩阵形式为

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{L}\\ \mathrm{\δ}\stackrel{\·}{\mathrm{\λ}}\\ \mathrm{\δ}\stackrel{\·}{h}\end{array}\right]=\left[\begin{array}{ccc}0& 0& -\stackrel{\·}{L}/(R_M+h)\\ \stackrel{\·}{\mathrm{\λ}}\tan L& 0& -\stackrel{\·}{\mathrm{\λ}}/(R_N+h)\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\δL}\\ \mathrm{\δ\λ}\\ \mathrm{\δh}\end{array}\right]+\left[\begin{array}{ccc}0& 1/(R_M+h)& 0\\ \sec L/(R_N+h)& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\mathrm{\δV}}_E\\ \mathrm{\δ}{V}_N\\ \mathrm{\δ}{V}_U\end{array}\right]---\left(3\right)$

式中δL、δλ、δh为位置误差，分别代表纬度误差、经度误差、高度误差；R_M和R_N分别为沿子午圈和卯酉圈的主曲率半径；L、λ、h分别为纬度、经度和高度，有

R_M＝R_e(1-2e+3e sin²L)，R_N＝R_e(1+e sin²L)，其中椭圆度e＝1/298.257，地球椭球长半轴R_e＝6378137m。

②SINS/GNSS组合导航滤波非线性数学模型

SINS/GNSS组合导航滤波的非线性数学模型包括状态方程和量测方程。

a)状态方程

综合以上SINS的姿态误差方程、速度误差方程和位置误差方程，并将SINS的陀螺漂移进行建模扩充为状态变量，考虑模型误差，可以得到如下系统状态方程：

$\stackrel{\·}{x}=f\left(x\right)+G^{d}d+G^{w}w---\left(4\right)$

式中x＝[δL δλ δh δV_E δV_N δV_U δq₀ δq₁ δq₂ δq₃ ε_x ε_y ε_z]^T为状态变量，ε_x、ε_y、ε_z分别为载体坐标系中x、y、z轴上陀螺常值漂移；d为模型误差；G^d为模型误差分布矩阵；G^w为系统噪声输入矩阵；w为系统噪声，Q为方差阵。

将式(4)表示成矩阵形式：

其中

$F_S={\left[\begin{array}{c}{0}_{3\×3}\\ 0_{3\×3}\\ \frac{1}{2}N\left({\hat{Q}}_b^n\right)\end{array}\right]}_{10\×3}$

系数矩阵F_N可写成如下形式：

其中

$F_{11}=\left[\begin{array}{ccc}0& 0& 0\\ \frac{V_E}{R_N+h}\sec L\tan L& 0& 0\\ 0& 0& 0\end{array}\right],$ $F_{12}=\left[\begin{array}{ccc}0& 1/(R_M+h)& 0\\ \sec L/(R_N+h)& 0& 0\\ 0& 0& 1\end{array}\right],$

$F_{21}=\left[\begin{array}{ccc}(\frac{V_E}{R_N+h}{\sec }^{2}L+2{\mathrm{\ω}}_{\mathrm{ie}}\cos L)V_N+2{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·V_U& 0& 0\\ -(\frac{V_E}{R_N+h}{\sec }^{2}L+2{\mathrm{\ω}}_{\mathrm{ie}}\cos L)V_E& 0& 0\\ -2{\mathrm{\ω}}_{\mathrm{ie}}\sin L\·V_E& 0& 0\end{array}\right],$

$F_{22}=\left[\begin{array}{ccc}{V}_N\tan L/(R_N+h)-V_U/(R_N+h)& (2{\mathrm{\ω}}_{\mathrm{ie}}+\stackrel{\·}{\mathrm{\λ}})\sin L& -({2\mathrm{\ω}}_{\mathrm{ie}}+\stackrel{\·}{\mathrm{\λ}})\cos L\\ -2({\mathrm{\ω}}_{\mathrm{ie}}+\stackrel{\·}{\mathrm{\λ}})\sin L& -V_U/(R_M+h)& -\stackrel{\·}{L}\\ 2({\mathrm{\ω}}_{\mathrm{ie}}+\stackrel{\·}{\mathrm{\λ}})\cos L& 2\stackrel{\·}{L}& 0\end{array}\right],$

$F_{23}=-2[{\hat{C}}_b^n+{\hat{f}}^b]\·Y^T\left(\hat{Q}\right)+2{\hat{C}}_b^n{\hat{f}}^b{\hat{Q}}^T,$

$F_{33}=\frac{1}{2}\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_x& -{\mathrm{\ω}}_y& -{\mathrm{\ω}}_z\\ {\mathrm{\ω}}_x& 0& {\mathrm{\ω}}_z& -{\mathrm{\ω}}_y\\ {\mathrm{\ω}}_y& -{\mathrm{\ω}}_z& 0& {\mathrm{\ω}}_x\\ {\mathrm{\ω}}_z& {\mathrm{\ω}}_y& -{\mathrm{\ω}}_x& 0\end{array}\right]-\frac{1}{2}\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_E& -{\mathrm{\ω}}_N& -{\mathrm{\ω}}_U\\ {\mathrm{\ω}}_E& 0& -{\mathrm{\ω}}_U& {\mathrm{\ω}}_N\\ {\mathrm{\ω}}_N& {\mathrm{\ω}}_U& 0& -{\mathrm{\ω}}_E\\ {\mathrm{\ω}}_U& -{\mathrm{\ω}}_N& {\mathrm{\ω}}_E& 0\end{array}\right].$

其中，ω_ie为地球自转角速度。

q(x，t)为非线性部分，其表达式为

$q(x,t)={\left[\begin{array}{c}{0}_{3\×1}\\ -Y^T\left(\mathrm{\δQ}\right)N\left(\mathrm{\δQ}\right){\hat{f}}^b\\ 0_{7\×1}\end{array}\right]}_{13\×1}$

将

简记为N′，即

$N^{\′}=-Y^T\left(\mathrm{\δQ}\right)N\left(\mathrm{\δQ}\right){\hat{f}}^b={\left[\begin{array}{c}{N}_1^{\′}\\ N_2^{\′}\\ N_3^{\′}\end{array}\right]}_{3\×1}$

式(4)中，G^d、G^w和d的具体表达式分别如下

$d={\left[\begin{array}{c}{d}_x\\ d_y\\ d_z\end{array}\right]}_{3\×1},$ $G^d={\left[\begin{array}{c}{0}_{3\×3}\\ C_b^n\\ 0_{7\×3}\end{array}\right]}_{13\×3},$ $G^w={\left[\begin{array}{c}{0}_{6\×3}\\ \frac{1}{2}N\left({\hat{Q}}_b^n\right)\\ 0_{3\×3}\end{array}\right]}_{13\×3}.$

其中，d_x、d_y、d_z代表模型误差分量，包含载体坐标系中x、y、z轴上的加速度计误差等未知模型误差量。

b)量测方程

取SINS捷联解算出的位置和速度与GNSS的速度和位置之差作为量测值，SINS/GNSS组合导航滤波量测方程为

y＝H(x)+v    (5)

式中y＝[δL′δλ′δh′δV′_E δV′_N δV′_U]^T，δL′，δλ′和δh′分别表示SINS与GNSS输出的纬度、经度和高度之差；δV′_E，δV′_N和δV′_U分别表示SINS与GNSS输出的东向、北向和天向速度之差；量测噪声

v_δL′、v_δλ′、v_δH′、

分别代表GNSS的纬度、经度、高度、东向速度、北向速度、天向速度的量测噪声；量测噪声方差阵R根据GNSS的位置、速度噪声水平选取；H(x)的具体表达式为H(x)＝[(R_M+h)δL(R_N+h)cosL·δλ δh δV_E δV_N δV_U]^T。

2、对传统二阶插值滤波进行改进，并基于1中建立的数学模型对状态进行估计

a)计算出t_k-1时刻协方差阵

和系统噪声方差阵Q_k-1均方根的各列向量

$s_{x,k-1}^i=\sqrt{3}{\left(\sqrt{P_{k-1}^x}\right)}_{i}i=\mathrm{1,2},...,13---\left(6\right)$

$s_{w,k-1}^i=\sqrt{3}{\left(\sqrt{Q_{k-1}}\right)}_{i}i=\mathrm{1,2,3}---\left(7\right)$

式中，下标k-1表示t_k-1时刻；

表示矩阵

的平方根矩阵中的第i列；

表示矩阵Q_k-1的平方根矩阵中的第i列。

b)计算系统状态一步预测

和状态一步预测协方差阵

${\hat{x}}_{k/k-1}=\frac{-13}{3}f({\hat{x}}_{k-1},w_{k-1})+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})]$

$+\frac{1}{6}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)]---\left(8\right)$

$P_{k/k-1}^x=\frac{1}{12}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}{[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})-f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})]}^2$

$+\frac{1}{12}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)-f{({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}{[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})-2f({\hat{x}}_{k-1},w_{k-1})]}^2---\left(9\right)$

$+\frac{1}{18}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)-2f({\hat{x}}_{k-1},w_{k-1})]}^2$

式中，下标k/k-1表示从t_k-1到t_k时刻的预测；

为t_k-1时刻的状态估计，w_k-1为系统噪声w在t_k-1时刻的值。

c)计算出t_k时刻状态协方差阵

和量测噪声方差阵R_k均方根的各列向量

$s_{x,k}^i=\sqrt{3}{\left(\sqrt{P_k^x}\right)}_{i}i=\mathrm{1,2},...,13---\left(10\right)$

$s_{v,k}^i=\sqrt{3}{\left(\sqrt{R_k}\right)}_{i}i=\mathrm{1,2},...,6---\left(11\right)$

式中，下标k表示t_k时刻；

表示矩阵

的平方根矩阵中的第i列；

表示矩阵R_k的平方根矩阵中的第i列。

d)量测修正更新，即计算量测一步预测

量测协方差阵

协方差阵滤波增益K_k、状态估计

和状态协方差阵

${\hat{y}}_{k/k-1}=\frac{16}{3}H({\hat{x}}_{k/k-1},v_k)+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)+H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]$

$+\frac{1}{6}\underset{i=1}{\overset{6}{\mathrm{\Σ}}}[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)+H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)]---\left(12\right)$

$P_k^y=\frac{1}{12}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)-H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]}^2$

$+\frac{1}{12}\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)-H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)+H({\hat{x}}_{k/k-1}-x_{x,k}^i,v_k)-2H({\hat{x}}_{k/k-1},v_k)]}^2---\left(13\right)$

$+\frac{1}{18}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)+H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)-2H({\hat{x}}_{k/k-1},v_k)]}^2$

$P_k^{\mathrm{xy}}=\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}{s}_{x,k}^i{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)-H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]}^T---\left(14\right)$

$K_k=P_k^{\mathrm{xy}}{\left[P_k^y\right]}^{-1}---\left(15\right)$

${\hat{x}}_k={\hat{x}}_{k/k-1}+K_k[y_k-{\hat{y}}_{k/k-1}]---\left(16\right)$

$P_k^x=P_{k/k-1}^x-K_k{P}_k^y{K}_k^T---\left(17\right)$

式中y_k和v_k分别为量测变量y和量测噪声v在t_k时刻的值。

e)状态修正更新。当t_k时刻的量测值y_k未到来时，

和

分别是状态x_k-1和x_k的最优估值。但是，当y_k到来时，

要优于

若以平滑值

取代式(8)中的

将使滤波估计的精度更高。因此，本发明为进一步提高状态估计

的精度，首先计算

然后用

替换

重新计算

和

${\hat{x}}_{k-1/k}={\hat{x}}_{k-1}+P_{k-1}^x{\left[f({\hat{x}}_k,w_k)\right]}^2{\left[P_{k/k-1}^x\right]}^{-1}[{\hat{x}}_k-{\hat{x}}_{k/k-1}]---\left(18\right)$

${\hat{x}}_{k/k-1}=\frac{-13}{3}f({\hat{x}}_{k-1/k},w_{k-1})+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\Σ}}}[f({\hat{x}}_{k-1/k}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1/k}-s_{x,k-1}^i,w_{k-1})]$ (19)

$+\frac{1}{6}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}[f({\hat{x}}_{k-1/k},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1/k},w_{k-1}-s_{w,k-1}^i)]$

${\hat{x}}_k={\hat{x}}_{k/k-1}+K_k[y_k-{\hat{y}}_{k/k-1}]---\left(20\right)$

综合以上a)至e)步骤，即为改进后的二阶插值滤波算法。

3、利用预测滤波估计模型误差d

a)由t_k-1时刻系统的状态估值计算系统的输出估值

令采样时间间隔为Δt(可根据SINS/GNSS组合导航的滤波周期来定)，根据量测方程(5)式，量测输出估值

的计算公式为

${\hat{y}}_{k-1}=H({\hat{x}}_{k-1},t_{k-1})---\left(21\right)$

b)估计出t_k-1时刻的模型误差

模型误差

的计算公式为

${\hat{d}}_{k-1}=-{\left\{{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)U\left({\hat{x}}_{k-1}\right)\right]}^T{R}^{-1}\right[\mathrm{\Λ}\left(\mathrm{\Δt}\right)U\left({\hat{x}}_{k-1}\right)]+W\}}^{-1}---\left(22\right)$

$\·{\left[\mathrm{\Λ}\left(\mathrm{\Δt}\right)U\left({\hat{x}}_{k-1}\right)\right]}^T{R}^{-1}[Z({\hat{x}}_{k-1},\mathrm{\Δt})-y_k+{\hat{y}}_{k-1}]$

式中参数

为6维列向量，其各分量为

$Z_i({\hat{x}}_{k-1},\mathrm{\Δt})=\underset{j=1}{\overset{p_i}{\mathrm{\Σ}}}\frac{{\mathrm{\Δt}}^j}{j!}{L}_f^j\left(H_i\right)i=1,...6---\left(23\right)$

式中H_i为H(x)的第i个分量；p_i为模型误差向量d的任何分量第一次出现在h_i的微分中的最低阶数；

的计算公式为

$L_f^0\left(H_i\right)=H_i$

$L_f^j\left(H_i\right)=\frac{\∂L_f^{j-1}\left(H_i\right)}{\∂\hat{x}}f({\hat{x}}_{k-1},t_{k-1})j\≥1---\left(24\right)$

求得

$Z_1=\mathrm{\Δt}(\left[\begin{array}{ccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x)$

$+\frac{{\mathrm{\Δt}}^2}{2}(\frac{1}{R_M+h}\·(\left[\begin{array}{ccccccccccccc}0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x+N_2^{\′}))$

$Z_2=\mathrm{\Δt}(\left[\begin{array}{ccccccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x)$

$+\frac{{\mathrm{\Δt}}^2}{2}(\frac{V_E\sec L\tan L}{R_N+h}(\left[\begin{array}{ccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x)$

$+\frac{\sec L}{R_N+h}(\left[\begin{array}{ccccccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x+N_1^{\′}))$

$Z_3=\mathrm{\Δt}(\left[\begin{array}{ccccccccccccc}0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x)$

$+\frac{{\mathrm{\Δt}}^2}{2}(\left[\begin{array}{ccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N$

$\·\left[\begin{array}{ccccccccccccc}0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right]F_N\·x+N_3^{\′}\left)\right)$

Z₄＝Δt([0 0 0 1 0 0 0 0 0 0 0 0 0]F_N·x+N′₁)

Z₅＝Δt([0 0 0 0 1 0 0 0 0 0 0 0 0]F_N·x+N′₂)

Z₆＝Δt([0 0 0 0 0 1 0 0 0 0 0 0 0]F_N·x+N′₃)

Λ(Δt)∈R^6×6为对角阵，其对角元素计算公式为

${\mathrm{\λ}}_{\mathrm{ii}}=\frac{{\mathrm{\Δt}}^{p_i}}{p_i!},i=1,...,6---\left(25\right)$

求得 $\mathrm{\Λ}\left(\mathrm{\Δt}\right)=\mathrm{diag}\left[\begin{array}{cccccc}\mathrm{\Δt}& \mathrm{\Δt}& \mathrm{\Δt}& \frac{{\left(\mathrm{\Δt}\right)}^2}{2}& \frac{{\left(\mathrm{\Δt}\right)}^2}{2}& \frac{{\left(\mathrm{\Δt}\right)}^2}{2}\end{array}\right]$

为灵敏度矩阵，各行元素可表示为

$U_i\left({\hat{x}}_{k-1}\right)=\left\{L_{g_1}\right[L_f^{p_i-1}\left(H_i\right)],...,L_{g_q}[L_f^{p_i-1}\left(H_i\right)\left]\right\}q=\mathrm{1,2,3};i=\mathrm{1,2}...,6---\left(26\right)$

式中，g_j为G^d的第j列，标量函数

关于向量场

的一阶李导数记为

其表达式为

$L_{g_j}\left[L_f^{p_i-1}\left(H_i\right)\right]=\frac{{\∂L}_f^{p_i-1}\left(H_i\right)}{\∂\hat{x}}{g}_j,j=\mathrm{1,2,3}---\left(27\right)$

求得

U₁₁＝C₁₁，U₁₂＝C₁₂，U₁₃＝C₁₃，U₂₁＝C₂₁，U₂₂＝C₂₂，U₂₃＝C₂₃，

U₃₁＝C₃₁，U₃₂＝C₃₂，U₃₃＝C₃₃，U₄₁＝C₂₁/R_M，U₄₂＝C₂₂/R_M，U₄₃＝C₂₃/R_M，

U₅₁＝C₁₁/R_N/cosL，U₅₂＝C₁₂/R_N/cosL，U₅₃＝C₁₃/R_N/cosL，

U₆₁＝C₃₁，U₆₂＝C₃₂，U₆₆＝C₃₃。

其中C_ij(i，j＝1，2，3)为

的第i行第j列元素，即

4、利用预测滤波估计出的模型误差

修正改进后的二阶插值滤波算法的状态估计，得到修正后的姿态误差、位置误差和速度误差

a)在改进二阶插值滤波

的计算式(18)中增加修正项

以修正模型误差引起的状态估计误差

${\hat{x}}_{k-1/k}={\hat{x}}_{k-1}+P_{k-1}^x{\left[f({\hat{x}}_k,w_k)\right]}^2{\left[P_{k/k-1}^x\right]}^{-1}[{\hat{x}}_k-{\hat{x}}_{k/k-1}]+G_{k-1}^d{\hat{d}}_{k-1}\mathrm{\Δt}---\left(28\right)$

式中，

为模型误差分布矩阵G^d在t_k-1时刻的值。

b)将步骤a)中计算出的

代入式(19)中计算

然后利用式(20)计算

最终计算得到的

即为最终修正后的状态变量，包括姿态误差(δq₀、δq₁、δq₂、δq₃)、位置误差(δL、δλ、δh)、速度误差(δV_E、δV_N、δV_U)等状态变量。

5、利用修正后的姿态误差、位置误差和速度误差修正SINS捷联解算出的姿态、位置和速度。将修正后的姿态、位置和速度作为下一导航时刻的初始值。

①姿态修正

利用修正后的姿态误差(δq₀、δq₁、δq₂、δq₃)计算出计算导航坐标系到真实导航坐标系间的转移矩阵然后对SINS计算出的姿态矩阵

(即上述的)进行校正，得到

最后根据

计算出更为准确的姿态——航向角

俯仰角θ和横滚角γ。

a)计算真实导航坐标系与计算导航坐标系间的转移矩阵

${\left[\begin{array}{ccc}{\mathrm{\φ}}_E& {\mathrm{\φ}}_N& {\mathrm{\φ}}_U\end{array}\right]}^T=-2{\left[Y\left({\hat{Q}}_b^n\right)\right]}^T\mathrm{\δQ}---\left(29\right)$

式中φ_E、φ_N和φ_U为导航坐标系东向、北向和天向轴上的数学平台失准角。

由φ_E、φ_N和φ_U可得

的表达式为：

$C_{n^{\′}}^n=\left[\begin{array}{ccc}1& -{\mathrm{\φ}}_U& {\mathrm{\φ}}_N\\ {\mathrm{\φ}}_U& 1& -{\mathrm{\φ}}_E\\ -{\mathrm{\φ}}_N& {\mathrm{\φ}}_E& 1\end{array}\right]---\left(30\right)$

b)计算姿态矩阵

利用

对SINS解算得到的

进行如下校正：

$C_b^n=C_{n^{\′}}^n{C}_b^{n^{\′}}---\left(31\right)$

式中

即为修正后的姿态矩阵。

c)计算航向角

俯仰角θ和横滚角γ航向角

俯仰角θ和横滚角γ的定义如图2a、图2b和图2c所示。

将(31)式计算出的

记为

$C_b^n=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]---\left(32\right)$

式中，T_ij(i＝1，2，3，j＝1，2，3)表示

的第i行第j列元素。

又因为

因此，由(32)式和(33)式，可以确定航向角

俯仰角θ和横滚角γ的主值，即

θ_主＝arcsin(T₃₂)

若航向角

俯仰角θ和横滚角γ的取值范围分别定义为[0，2π]、

[-π，+π]。那么，

θ和γ的真值可由下式确定：

θ＝θ_主    (35)

由(35)式确定的

θ和γ即为经过修正后的航向角、俯仰角和横滚角。

②速度修正

由SINS捷联解算出的东向速度V_E、北向速度V_N和天向速度V_U分别减去速度误差δV_E、δV_N和δV_U，得到修正后的速度：

$\left[\begin{array}{c}{V_E}^{\′}\\ {V_N}^{\′}\\ {V_U}^{\′}\end{array}\right]=\left[\begin{array}{c}{V}_E-\mathrm{\δ}{V}_E\\ V_N-\mathrm{\δ}{V}_N\\ V_U-\mathrm{\δ}{V}_U\end{array}\right]---\left(36\right)$

式中V′_E、V′_N和V′_U为修正后的东向速度、北向速度和天向速度。

③位置修正

由SINS捷联解算出的纬度L、经度λ和高度h分别减去位置误差δL、δλ和δh，得到修正后的位置：

$\left[\begin{array}{c}{L}^{\′}\\ {\mathrm{\λ}}^{\′}\\ h^{\′}\end{array}\right]=\left[\begin{array}{c}L-\mathrm{\δL}\\ \mathrm{\λ}-\mathrm{\δ\λ}\\ h-\mathrm{\δh}\end{array}\right]---\left(37\right)$

式中L′、λ′和h′为修正后的纬度、经度和高度。

本发明未详细阐述部分属于本领域公知技术。

Claims (4)

1.一种捷联惯性/卫星组合导航滤波方法，其特征在于实现步骤如下：

(1)建立包含模型误差的捷联惯性导航系统/全球导航卫星系统组合导航滤波非线性数学模型；

(2)对传统二阶插值滤波进行改进，并基于步骤(1)建立的数学模型对状态进行估计，然后利用预测滤波估计出的模型误差修正改进二阶插值滤波的状态估计，得到修正后的姿态误差、位置误差和速度误差；

(3)利用以上步骤(2)得到的修正后的姿态误差、位置误差和速度误差对捷联惯性导航系统捷联解算出的姿态、位置和速度进行补偿，获得更加准确的导航信息，即补偿后的姿态、位置和速度信息；

(4)将以上步骤(3)中补偿后的位置、速度和姿态作为下一导航时刻的初始值，不断重复以上步骤(2)和步骤(3)，直至捷联惯性导航系统/全球导航卫星系统组合导航结束；

所述的步骤(1)中，包含模型误差的捷联惯性导航系统/全球导航卫星系统组合导航滤波的非线性数学模型包括状态方程和量测方程，其中状态方程为：

$\stackrel{\·}{x}=f\left(x\right)+G^{d}d+G^{w}w$

其中，x＝[δL δλ δh δV_E δV_N δV_U δq₀ δq₁ δq₂ δq₃ ε_x ε_y ε_z]^T为状态变量；δL、δλ、δh为位置误差，分别代表纬度误差、经度误差、高度误差；δV_E、δV_N、δV_U为速度误差，分别代表东向速度误差、北向速度误差、天向速度误差；δq₀、δq₁、δq₂、δq₃为姿态误差，即计算四元数

与真实四元数(q₀、q₁、q₂、q₃)各分量的差；ε_x、ε_y、ε_z分别为载体坐标系中x、y、z轴上陀螺常值漂移；d为模型误差；G^d为模型误差分布矩阵；G^w为系统噪声输入矩阵；w为系统噪声，Q为系统噪声方差阵；f(x)的具体表达式如下

f(x)＝[f₁ f₂ f₃ f₄]^T

$f_1=\left[\begin{array}{ccc}0& 0& -\stackrel{\·}{L}/(R_M+h)\\ \stackrel{\·}{\mathrm{\λ}}\tan L& 0& -\stackrel{\·}{\mathrm{\λ}}/(R_N+h)\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\δL}\\ \mathrm{\δ\λ}\\ \mathrm{\δh}\end{array}\right]+\left[\begin{array}{ccc}0& 1/(R_M+h)& 0\\ \sec L/(R_N+h)& 0& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}{V}_E\\ \mathrm{\δ}{V}_N\\ \mathrm{\δ}{V}_U\end{array}\right]$

$f_2=-2\left[{\hat{C}}_b^n{\hat{f}}^b\right]Y^T\left({\hat{Q}}_b^n\right)\mathrm{\δQ}+2{\hat{C}}_b^n{\hat{f}}^b{\left({\hat{Q}}_b^n\right)}^T\mathrm{\δQ}-Y^T\left(\mathrm{\δQ}\right)N\left(\mathrm{\δQ}\right){\hat{f}}^b$

$-(2\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ie}}^n+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{en}}^n)V^n-(2{\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\mathrm{\δV}$

$f_3=(\frac{1}{2}<{\mathrm{\&omega;}}_{\mathrm{ib}}^b>-\frac{1}{2}[{\mathrm{\&omega;}}_{\mathrm{in}}^n\left]\right)\mathrm{\&delta;Q}+(\frac{1}{2}N\left({\hat{Q}}_b^n\right)\&CenterDot;\mathrm{\&delta;}{\mathrm{\&omega;}}_{\mathrm{ib}}^b-\frac{1}{2}Y\left({\hat{Q}}_b^n\right)\&CenterDot;\mathrm{\&delta;}{\mathrm{\&omega;}}_{\mathrm{in}}^n)$

f₄＝[0 0 0]^T

G^d、G^w和d的具体表达式如下

$d={\left[\begin{array}{c}{d}_x\\ d_y\\ d_z\end{array}\right]}_{3\&times;1},$ $G^d={\left[\begin{array}{c}{0}_{3\&times;3}\\ {\hat{C}}_b^n\\ 0_{7\&times;3}\end{array}\right]}_{13\&times;3},$ $G^w={\left[\begin{array}{c}{0}_{6\&times;3}\\ \frac{1}{2}N\left({\hat{Q}}_b^n\right)\\ 0_{3\&times;3}\end{array}\right]}_{13\&times;3}$

其中，R_M和R_N分别为沿子午圈和卯酉圈的主曲率半径；L、λ、h分别为捷联惯性导航系统捷联解算出的纬度、经度和高度；Vⁿ＝[V_E V_N V_U]^T为捷联惯性导航系统捷联解算出的载体相对地球的速度在导航坐标系中的投影，V_E、V_N、V_U分别代表东向、北向和天向速度； $\stackrel{\&CenterDot;}{L}=V_N/(R_M+h),$ $\stackrel{\&CenterDot;}{\mathrm{\&lambda;}}=V_E/(R_N+h)/\cos L;$ δ_V＝[δV_E δV_N δV_U]^T；

为SINS捷联解算出的载体坐标系到计算导航坐标系的方向余弦矩阵，也称为计算姿态矩阵； $\mathrm{\&delta;Q}={\hat{Q}}_b^n-Q_b^n={\left[\begin{array}{cccc}\mathrm{\&delta;}{q}_0& \mathrm{\&delta;}{q}_1& \mathrm{\&delta;}{q}_2& \mathrm{\&delta;}{q}_3\end{array}\right]}^T,$

和

分别为载体坐标系到导航坐标系的真实四元数和计算四元数；

为加速度计的测量值；

为地球坐标系相对惯性系的旋转角速度在导航坐标系的投影；

为导航坐标系相对地球坐标系的旋转角速度在导航坐标系的投影； ${\widehat{\mathrm{\&omega;}}}_{\mathrm{ie}}^n={\mathrm{\&omega;}}_{\mathrm{ie}}^n+\mathrm{\&delta;}{\mathrm{\&epsiv;}}_{\mathrm{ie}}^n$ 和 ${\widehat{\mathrm{\&omega;}}}_{\mathrm{en}}^n={\mathrm{\&omega;}}_{\mathrm{en}}^n+\mathrm{\&delta;}{w}_{\mathrm{en}}^n$ 分别为

和

的计算值； ${\widehat{\mathrm{\&omega;}}}_{\mathrm{ib}}^b={\mathrm{\&omega;}}_{\mathrm{ib}}^b+\mathrm{\&delta;}{\mathrm{\&omega;}}_{\mathrm{ib}}^b$ 为陀螺实际测量值， ${\mathrm{\&omega;}}_{\mathrm{ib}}^b={\left[\begin{array}{ccc}{\mathrm{\&omega;}}_x& {\mathrm{\&omega;}}_y& {\mathrm{\&omega;}}_z\end{array}\right]}^T$ 为陀螺的理论输出，ω_x、ω_y、ω_z分别代表载体坐标系中x、y、z轴上陀螺的理论输出，

表示陀螺的测量误差； ${\widehat{\mathrm{\&omega;}}}_{\mathrm{in}}^n={\mathrm{\&omega;}}_{\mathrm{in}}^n+\mathrm{\&delta;}{\mathrm{\&omega;}}_{\mathrm{in}}^n$ 为计算导航坐标系相对惯性坐标系的旋转角速度， ${\mathrm{\&omega;}}_{\mathrm{in}}^n={\left[\begin{array}{ccc}{\mathrm{\&omega;}}_E& {\mathrm{\&omega;}}_N& {\mathrm{\&omega;}}_U\end{array}\right]}^T$ 为真实导航坐标系相对惯性坐标系的旋转角速度在导航坐标系的投影，ω_E、ω_N、ω_U分别为在真实导航坐标系中东向、北向、天向轴上的投影；

为

的计算误差；d_x、d_y、d_z代表模型误差分量，包含载体坐标系中x、y、z轴上的加速度计误差等未知模型误差量；N(δQ)和Y(δQ)是由δq₀、δq₁、δq₂、δq₃组成的矩阵，

和

是由

组成的矩阵，

和

是分别由ω_x、ω_y、ω_z和ω_E、ω_N、ω_U组成的矩阵，N(δQ)、

Y(δQ)、

和

的具体定义如下

$N\left(\mathrm{\&delta;Q}\right)\&equiv;\left[\begin{array}{ccc}-\mathrm{\&delta;}{q}_1& -\mathrm{\&delta;}{q}_2& -\mathrm{\&delta;}{q}_3\\ \mathrm{\&delta;}{q}_0& -\mathrm{\&delta;}{q}_3& \mathrm{\&delta;}{q}_2\\ \mathrm{\&delta;}{q}_3& \mathrm{\&delta;}{q}_0& -\mathrm{\&delta;}{q}_1\\ -\mathrm{\&delta;}{q}_2& \mathrm{\&delta;}{q}_1& \mathrm{\&delta;}{q}_0\end{array}\right];$ $N\left({\hat{Q}}_b^n\right)\&equiv;\left[\begin{array}{ccc}-{\hat{q}}_1& -{\hat{q}}_2& -{\hat{q}}_3\\ {\hat{q}}_0& -{\hat{q}}_3& {\hat{q}}_2\\ {\hat{q}}_3& {\hat{q}}_0& -{\hat{q}}_1\\ -{\hat{q}}_2& {\hat{q}}_1& {\hat{q}}_0\end{array}\right];$

$Y\left(\mathrm{\&delta;Q}\right)\&equiv;\left[\begin{array}{ccc}-\mathrm{\&delta;}{q}_1& -\mathrm{\&delta;}{q}_2& -\mathrm{\&delta;}{q}_3\\ \mathrm{\&delta;}{q}_0& \mathrm{\&delta;}{q}_3& -\mathrm{\&delta;}{q}_2\\ -\mathrm{\&delta;}{q}_3& \mathrm{\&delta;}{q}_0& \mathrm{\&delta;}{q}_1\\ \mathrm{\&delta;}{q}_2& -\mathrm{\&delta;}{\hat{q}}_1& \mathrm{\&delta;}{q}_0\end{array}\right];$ $Y\left({\hat{Q}}_b^n\right)\&equiv;\left[\begin{array}{ccc}-{\hat{q}}_1& -{\hat{q}}_2& -{\hat{q}}_3\\ {\hat{q}}_0& {\hat{q}}_3& {-\hat{q}}_2\\ {-\hat{q}}_3& {\hat{q}}_0& {\hat{q}}_1\\ {\hat{q}}_2& {-\hat{q}}_1& {\hat{q}}_0\end{array}\right];$

$<{\mathrm{\&omega;}}_{\mathrm{ib}}^b>\&equiv;\left[\begin{array}{cccc}0& -{\mathrm{\&omega;}}_x& -{\mathrm{\&omega;}}_y& -{\mathrm{\&omega;}}_z\\ {\mathrm{\&omega;}}_x& 0& {\mathrm{\&omega;}}_z& -{\mathrm{\&omega;}}_y\\ {\mathrm{\&omega;}}_y& -{\mathrm{\&omega;}}_z& 0& {\mathrm{\&omega;}}_x\\ {\mathrm{\&omega;}}_z& {\mathrm{\&omega;}}_y& -{\mathrm{\&omega;}}_x& 0\end{array}\right];$ $\left[{\mathrm{\&omega;}}_{\mathrm{in}}^n\right]\&equiv;\left[\begin{array}{cccc}0& -{\mathrm{\&omega;}}_E& -{\mathrm{\&omega;}}_N& -{\mathrm{\&omega;}}_U\\ {\mathrm{\&omega;}}_E& 0& -{\mathrm{\&omega;}}_U& {\mathrm{\&omega;}}_N\\ {\mathrm{\&omega;}}_N& {\mathrm{\&omega;}}_U& 0& -{\mathrm{\&omega;}}_E\\ {\mathrm{\&omega;}}_U& -{\mathrm{\&omega;}}_N& {\mathrm{\&omega;}}_E& 0\end{array}\right];$

取捷联惯性导航系统捷联解算与全球导航卫星系统的速度和位置之差作为量测值，量测方程为：

y＝H(x)+v

其中y＝[δL′δλ′δh′δV′_E δV′_N δV′_U]^T为量测变量，δL′，δλ′和δh′分别表示捷联惯性导航系统与全球导航卫星系统输出的纬度、经度和高度之差；δV′_E，δV′_N和δV′_U分别表示捷联惯性导航系统与全球导航卫星系统输出的东向、北向和天向速度之差；量测噪声 $v={\left[\begin{array}{cccccc}{v}_{\mathrm{\&delta;}{L}^{\&prime;}}& v_{\mathrm{\&delta;}{\mathrm{\&lambda;}}^{\&prime;}}& v_{{\mathrm{\&delta;H}}^{\&prime;}}& v_{{\mathrm{\&delta;V}}_E^{\&prime;}}& v_{\mathrm{\&delta;}{V}_N^{\&prime;}}& v_{{\mathrm{\&delta;V}}_U^{\&prime;}}\end{array}\right]}^T,$ v_δL′、v_δλ′、v_δH′、

分别代表全球导航卫星系统的纬度、经度、高度、东向速度、北向速度、天向速度的量测噪声；量测噪声方差阵R根据全球导航卫星系统位置、速度噪声水平选取；H(x)的具体表达式为H(x)＝[(R_M+h)δL(R_N+h)cosL·δλ δh δV_E δV_N δV_U]^T；

所述的步骤(2)中对传统二阶插值滤波进行改进，并基于步骤(1)建立的数学模型对状态进行估计的步骤为：

a)计算出t_k-1时刻协方差阵

和系统噪声方差阵Q_k-1均方根的各列向量

$s_{x,k-1}^i=\sqrt{3}{\left(\sqrt{P_{k-1}^x}\right)}_i$ i＝1，2，…，13

$s_{w,k-1}^i=\sqrt{3}{\left(\sqrt{Q_{k-1}}\right)}_i$ i＝1，2，3

其中，下标k-1表示t_k-1时刻；

表示矩阵

的平方根矩阵中的第i列；

表示矩阵Q_k-1的平方根矩阵中的第i列；

b)计算系统状态一步预测

和状态一步预测协方差阵

${\hat{x}}_{k/k-1}=\frac{-13}{3}f({\hat{x}}_{k-1},w_{k-1})+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})]$

$+\frac{1}{6}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)]$

$P_{k/k-1}^x=\frac{1}{12}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}{[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})-f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})]}^2$

$+\frac{1}{12}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)-f({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}{[f({\hat{x}}_{k-1}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1}-s_{x,k-1}^i,w_{k-1})-2f({\hat{x}}_{k-1},w_{k-1})]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{[f({\hat{x}}_{k-1},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1},w_{k-1}-s_{w,k-1}^i)-2f({\hat{x}}_{k-1},w_{k-1})]}^2$

其中，下标k/k-1表示从t_k-1到t_k时刻的预测；为t_k-1时刻的状态估计，w_k-1为系统噪声w在t_k-1时刻的值。

c)计算出t_k时刻状态协方差阵和量测噪声方差阵R_k均方根的各列向量

$s_{x,k}^i=\sqrt{3}{\left(\sqrt{P_k^x}\right)}_{i}i=\mathrm{1,2},...,13$

$s_{v,k}^i=\sqrt{3}{\left(\sqrt{R_k}\right)}_{i}i=\mathrm{1,2},...,6$

其中，下标k表示t_k时刻；

表示矩阵的平方根矩阵中的第i列；表示矩阵R_k的平方根矩阵中的第i列；

d)量测修正更新，即计算量测一步预测

量测协方差阵协方差阵

滤波增益K_k、状态估计和状态协方差阵

${\hat{y}}_{k/k-1}=\frac{16}{3}H({\hat{x}}_{k/k-1},v_k)+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)+H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]$

$+\frac{1}{6}\underset{i=1}{\overset{6}{\mathrm{\&Sigma;}}}[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)+H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)]$

$P_k^y=\frac{1}{12}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)-H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]}^2$

$+\frac{1}{12}\underset{i=1}{\overset{6}{\mathrm{\&Sigma;}}}{[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)-H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)+H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)-2H({\hat{x}}_{k/k-1},v_k)]}^2$

$+\frac{1}{18}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{[H({\hat{x}}_{k/k-1},v_k+s_{v,k}^i)+H({\hat{x}}_{k/k-1},v_k-s_{v,k}^i)-2H({\hat{x}}_{k/k-1},v_k)]}^2$

$P_k^{\mathrm{xy}}=\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}{s}_{x,k}^i{[H({\hat{x}}_{k/k-1}+s_{x,k}^i,v_k)-H({\hat{x}}_{k/k-1}-s_{x,k}^i,v_k)]}^T$

$K_k=P_k^{\mathrm{xy}}{\left[P_k^y\right]}^{-1}$

${\hat{x}}_k={\hat{x}}_{k/k-1}+K_k[y_k-{\hat{y}}_{k/k-1}]$

$P_k^x=P_{k/k-1}^x-K_k{P}_k^y{K}_k^T$

其中，y_k和v_k分别为量测变量y和量测噪声v在t_k时刻的值；

e)状态修正更新，首先计算平滑值

然后用

替换

重新计算

和

${\hat{x}}_{k-1/k}={\hat{x}}_{k-1}+P_{k-1}^x{\left[f({\hat{x}}_k,w_k)\right]}^2{\left[P_{k/k-1}^x\right]}^{-1}[{\hat{x}}_k-{\hat{x}}_{k/k-1}]$

${\hat{x}}_{k/k-1}=\frac{-13}{3}f({\hat{x}}_{k-1/k},w_{k-1})+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1/k}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1/k}-s_{x,k-1}^i,w_{k-1})]$

$+\frac{1}{6}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1/k},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1/k},w_{k-1}-s_{w,k-1}^i)]$

${\hat{x}}_k={\hat{x}}_{k/k-1}+K_k[y_k-{\hat{y}}_{k/k-1}]$

综合以上a)至e)步骤，即为改进后的二阶插值滤波算法。

2.根据权利要求1所述的捷联惯性/卫星组合导航滤波方法，其特征在于：所述的步骤(1)中预测滤波估计模型误差d的步骤为：

a)由t_k-1时刻系统的状态估值

计算系统的输出估值

计算公式为

${\hat{y}}_{k-1}=H({\hat{x}}_{k-1},t_{k-1})$

b)估计出t_k-1时刻的模型误差

计算公式为

${\hat{d}}_{k-1}=-{\left\{{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)U\left({\hat{x}}_{k-1}\right)\right]}^T{R}^{-1}\right[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)U\left({\hat{x}}_{k-1}\right)]+W\}}^{-1}$

$\&CenterDot;{\left[\mathrm{\&Lambda;}\left(\mathrm{\&Delta;t}\right)U\left({\hat{x}}_{k-1}\right)\right]}^T{R}^{-1}[Z({\hat{x}}_{k-1},\mathrm{\&Delta;t})-y_k+{\hat{y}}_{k-1}]$

其中，

为6维列向量，其各分量为

$Z_i({\hat{x}}_{k-1},\mathrm{\&Delta;t})=\underset{j=1}{\overset{p_i}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&Delta;}{t}^j}{j!}{L}_f^j\left(H_i\right)$ i＝1，...6， $L_f^0\left(H_i\right)=H_i,$ $L_f^j\left(H_i\right)=\frac{\&PartialD;L_f^{j-1}\left(H_i\right)}{\&PartialD;\hat{x}}f({\hat{x}}_{k-1},t_{k-1})$ j≥1

H_i为H(x)的第i个分量；

Λ(Δt)∈R^6×6为对角阵，其对角元素为

${\mathrm{\&lambda;}}_{\mathrm{ii}}=\frac{\mathrm{\&Delta;}{t}^{p_i}}{p_i!},$ i＝1，...，6

$U\left({\hat{x}}_{k-1}\right)\&Element;R^{6\&times;3}$ 的各行元素可表示为

$U_i\left({\hat{x}}_{k-1}\right)=\left\{\begin{array}{ccc}{L}_{g_1}\left[L_f^{p_i-1}\left(H_i\right)\right]& ,...,& L_{g_q}\left[L_f^{p_i-1}\left(H_i\right)\right]\end{array}\right\}$ q＝1，2，3；i＝1，2...，6

其中，g_j为G^d的第j列，

表达式为

$L_{g_j}\left[L_f^{p_i-1}\left(H_i\right)\right]=\frac{\&PartialD;L_f^{p_i-1}\left(H_i\right)}{\&PartialD;\hat{x}}{g}_j,$ j＝1，2，3。

3.根据权利要求1所述的捷联惯性/卫星组合导航滤波方法，其特征在于：所述的步骤(2)中利用预测滤波估计出的模型误差修正改进二阶插值滤波算法的状态估计，得到修正后的姿态误差、位置误差和速度误差的步骤为：

a)在改进二阶插值滤波

的计算式中增加修正项

以修正模型误差引起的状态估计误差

${\hat{x}}_{k-1/k}={\hat{x}}_{k-1}+P_{k-1}^x{\left[f({\hat{x}}_k,w_k)\right]}^2{\left[P_{k/k-1}^x\right]}^{-1}[{\hat{x}}_k-{\hat{x}}_{k/k-1}]+G_{k-1}^d{\hat{d}}_{k-1}\mathrm{\&Delta;t}$

其中，为模型误差分布矩阵G^d在t_k-1时刻的值；

b)利用步骤a)中计算出的

重新计算

和

${\hat{x}}_{k/k-1}=\frac{-13}{3}f({\hat{x}}_{k-1/k},w_{k-1})+\frac{1}{6}\underset{i=1}{\overset{13}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1/k}+s_{x,k-1}^i,w_{k-1})+f({\hat{x}}_{k-1/k}-s_{x,k-1}^i,w_{k-1})]$

$+\frac{1}{6}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}[f({\hat{x}}_{k-1/k},w_{k-1}+s_{w,k-1}^i)+f({\hat{x}}_{k-1/k},w_{k-1}-s_{x,k-1}^i)]$

${\hat{x}}_k={\hat{x}}_{k/k-1}+K_k[y_k-{\hat{y}}_{k/k-1}]$

以上获得的即为最终修正后的状态变量，包括姿态误差、位置误差、速度误差。

4.根据权利要求1所述的捷联惯性/卫星组合导航滤波方法，其特征在于：所述的步骤(3)修正SINS捷联解算出的姿态、位置和速度的公式为：

a)姿态修正

计算数学平台失准角φ_E、φ_N和φ_U：

${\left[\begin{array}{ccc}{\mathrm{\&phi;}}_E& {\mathrm{\&phi;}}_N& {\mathrm{\&phi;}}_U\end{array}\right]}^T=-2{\left[Y\left({\hat{Q}}_b^n\right)\right]}^T\mathrm{\&delta;Q}$

式中φ_E、φ_N和φ_U为导航坐标系东向、北向和天向轴上的数学平台失准角， $\mathrm{\&delta;Q}={\hat{Q}}_b^n-Q_b^n={\left[\begin{array}{cccc}{\mathrm{\&delta;q}}_0& \mathrm{\&delta;}{q}_1& \mathrm{\&delta;}{q}_2& \mathrm{\&delta;}{q}_3\end{array}\right]}^T$ 为姿态误差， $Q_b^n={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T$ 和 ${\hat{Q}}_b^n={\left[\begin{array}{cccc}{\hat{q}}_0& {\hat{q}}_1& {\hat{q}}_2& {\hat{q}}_3\end{array}\right]}^T$ 分别为载体坐标系到导航坐标系的真实四元数和计算四元数，

是由

组成的矩阵， $Y\left({\hat{Q}}_b^n\right)\&equiv;\left[\begin{array}{ccc}-{\hat{q}}_1& -{\hat{q}}_2& -{\hat{q}}_3\\ {\hat{q}}_0& {\hat{q}}_3& -{\hat{q}}_2\\ -{\hat{q}}_3& {\hat{q}}_0& {\hat{q}}_1\\ {\hat{q}}_2& -{\hat{q}}_1& {\hat{q}}_0\end{array}\right];$

利用φ_E、φ_N和φ_U组成姿态修正矩阵

$C_{n^{\&prime;}}^n=\left[\begin{array}{ccc}1& -{\mathrm{\&phi;}}_U& {\mathrm{\&phi;}}_N\\ {\mathrm{\&phi;}}_U& 1& -{\mathrm{\&phi;}}_E\\ -{\mathrm{\&phi;}}_N& {\mathrm{\&phi;}}_E& 1\end{array}\right]$

利用

修正捷联惯性导航系统捷联解算出的姿态矩阵

$C_b^n=C_{n^{\&prime;}}^n{\hat{C}}_b^n$

其中，

即为修正后的姿态矩阵，根据

即可计算出捷联惯性导航系统/全球导航卫星系统组合导航的姿态角，所述姿态角包括航向角、俯仰角和横滚角；

b)速度修正

$\left[\begin{array}{c}{V_E}^{\&prime;}\\ {V_N}^{\&prime;}\\ {V_U}^{\&prime;}\end{array}\right]=\left[\begin{array}{c}{V}_E-\mathrm{\&delta;}{V}_E\\ V_N-\mathrm{\&delta;}{V}_N\\ V_U-\mathrm{\&delta;}{V}_U\end{array}\right]$

其中，V′_E、V′_N和V′_U为修正后的东向速度、北向速度和天向速度；

c)位置修正

$\left[\begin{array}{c}{L}^{\&prime;}\\ {\mathrm{\&lambda;}}^{\&prime;}\\ h^{\&prime;}\end{array}\right]=\left[\begin{array}{c}L-\mathrm{\&delta;L}\\ \mathrm{\&lambda;}-\mathrm{\&delta;\&lambda;}\\ h-\mathrm{\&delta;h}\end{array}\right]$

其中，L′、λ′和h′为修正后的纬度、经度和高度。

Images