CN103090865A - 一种调制型捷联惯性导航系统姿态误差抑制方法 - Google Patents

Abstract

本发明提供了一种调制型捷联惯性导航系统姿态误差抑制方法，旋转机构带动惯性组件旋转至IMU坐标系与载体坐标系重合的位置，惯性组件进行单轴正反转停运动，依据旋转调制下导航参数设计原则，设计导航参数，实时采集光纤陀螺仪和石英加速度计测量载体运动的角速度和线加速度，修正控制角速率，利用IMU测量信息以及修正的控制角速率进行导航解算，得到载体姿态信息，更新载体运动速度，得到解算的速度误差，将得到的载体姿态信息和速度误差作为系统最终输出的导航信息。本发明方法消除了解算姿态信息中与旋转调制有关的振荡误差，提高了导航精度，增强了系统解算姿态信息的适用性。

Description

一种调制型捷联惯性导航系统姿态误差抑制方法

技术领域

本发明涉及一种调制型光纤陀螺捷联惯性导航系统姿态误差抑制方法，属于惯性技术领域中减小导航信息误差的抑制方法。

背景技术

捷联惯导系统SINS作为能够连续输出载体速度、姿态、位置信息的全自主导航系统被广泛应用于航空、航天、航海等领域。其主要是利用陀螺仪和加速度计测量载体的角运动和线运动信息，经导航解算后得到导航信息。然而，由于惯性组件输出值与测量值总是存在一定偏差，即存在惯性组件常值偏差，导致系统解算定位误差随导航时间的增长而逐渐发散，影响系统导航精度，制约了捷联惯导系统的长时间导航能力。

为了提高系统定位精度，一方面可以提高惯性元件精度，但是由于受加工技术水平的限制，无限制地提高元件精度是很难实现的；另一方面就是采取捷联惯性导航系统的误差抑制技术，自动抵消惯性器件的误差对系统精度的影响。这样就可以应用现有精度的惯性元件构成较高精度的捷联惯性导航系统。

旋转调制技术作为一种误差抑制技术，主要是通过旋转机构带动惯性组件按照已经设定好的旋转方案旋转，使得惯性组件常值偏差沿导航系投影呈周期振荡形式，一个旋转周期内积分结果为零，进而抵消惯性组件常值输出误差对定位误差影响，提高导航精度。虽然旋转调制能够抑制系统发散式定位误差，但又为系统解算姿态信息带来了与旋转周期有关的新振荡误差，使姿态精度降低，制约了该信息的可用性。并且旋转周期越长，该姿态误差振荡幅值越大。

《中国惯性技术学报》2009年17卷第1期由翁海娜等人撰写的《旋转式光学陀螺捷联惯导系统的旋转方案设计》，针对惯性测量组件输出误差的特性，设计单轴正反转停和双轴转位的系统旋转方案；《北京航空航天大学学报》2012年38卷第4期由杨国梁等人撰写的《基于双轴旋转的惯导系统误差自补偿技术》，提出了一种通过单元体的连续正反旋转，双单元体结构的误差自补偿方案，可以实现对所有方向上陀螺常值漂移的调制；公开号为101514899的中国发明专利在2009年8月26日公开的《基于单轴旋转的光纤陀螺捷联惯性导航系统误差抑制方法》，主要针对旋转调制这一环境，设计了单轴正反转停旋转方案，该方案能够有效地抑制与水平轴惯性组件常值偏差。以上文献都是提出了能够抑制定位误差的旋转方案，但并没有分析研究旋转状态下，系统解算姿态信息的精度和姿态误差的形式，也没有给出相应的误差抑制方案。

发明内容

本发明针对旋转状态下的姿态误差的抑制，提出了一种调制型捷联惯性导航系统姿态误差抑制方法。本发明方法依据旋转调制下导航参数设计原则，在制定旋转方案的基础上设计导航参数，通过调整修正控制角速率值来消除系统解算姿态信息中与旋转调制有关的姿态误差，提高姿态精度，增强系统解算姿态信息的适用性。

一种调制型捷联惯性导航系统姿态误差抑制方法，包括如下步骤：

步骤1：旋转机构带动惯性组件旋转至惯性组件（Inertial Measurement Unit，简称IMU）坐标系与载体坐标系重合的位置，有其中b表示载体坐标系，s表示IMU坐标系，

表示s系到b系转换矩阵，I表示单位阵；

步骤2：将光纤陀螺捷联惯导系统进行充分预热后，旋转机构带动惯性组件以角速度ω进行单轴正反转停运动；惯性组件采用四个转停次序为一个旋转周期的旋转方案；

步骤3：依据导航参数设计原则，根据惯性组件的旋转调制周期T设计导航参数；具体是：首先，确定导航解算过程中固有频率ω₀的取值范围：

ω₁表示惯性组件的旋转调制周期的调制频率ω₁＝2π/T；ξ为光纤陀螺捷联惯导系统阻尼系数，取值范围0～1；k表示光纤陀螺捷联惯导系统解算姿态误差角受IMU常值偏差的影响降低倍数；然后，确定根据下式确定导航参数K₁、K₂、K_E、K_N和K_U：g表示重力加速度；Ω＝0004167°/s，表示地球自转角速度；

表示载体所在位置的地理纬度；

步骤4：实时采集光纤陀螺仪和石英加速度计测量载体运动的角速度和线加速度，依据下式修正控制角速率：

其中， ${\mathrm{\ω}}_{\mathrm{cx}}^n,{\mathrm{\ω}}_{\mathrm{cy}}^n,{\mathrm{\ω}}_{\mathrm{cz}}^n$ 分别表示控制角速率

在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量，δv_j（j＝x,y）表示光纤陀螺捷联惯导系统解算速度误差；p表示复数域参变量；

步骤5：利用IMU测量信息以及修正的控制角速率进行导航解算，得到载体姿态信息；

步骤6：更新载体运动速度信息，得到光纤陀螺捷联惯导系统解算速度误差；

步骤7：重复步骤4至步骤6，直至光纤陀螺捷联惯导系统导航结束，不断更新修正控制角速率、速度误差和载体姿态信息，将抑制后的载体姿态信息和速度误差作为光纤陀螺捷联惯导系统最终输出导航信息。

本发明提供了一种针对调制型捷联惯导系统的姿态误差抑制方法，相比现有技术，具有如下优点和积极效果：本发明方法依据导航参数设计原则，根据一个旋转调制周期设计导航参数，利用该导航参数实时更新修正控制角速率，进而消除解算姿态信息中与旋转调制有关的振荡误差；此外，根据惯导系统基本原理可知，系统解算各导航信息包含与舒勒和地球周期有关的振荡误差，按本发明提出的导航参数设计方法可一并消除这两项振荡误差，使得系统解算姿态稳定，提高导航精度。

附图说明

图1为本发明的调制型捷联惯性导航系统姿态误差抑制方法的流程图；

图2为利用本发明的步骤2中IMU四位置转停示意图；

图3为本发明实施例1中利用Visual C++仿真得到未使用本发明方法解算的姿态误差曲线；

图4为本发明实施例1中利用Visual C++仿真得到使用本发明方法解算的姿态误差曲线；

图5为本发明实施例2中利用单轴转台试验，未使用本发明方法解算的姿态误差曲线；

图6为本发明实施例2中利用单轴转台试验，使用本发明方法解算的姿态误差曲线。

具体实施方式

下面结合附图对本发明的具体实施方式进行详细描述。

如图1所示，本发明提供一种针对调制型捷联惯导系统的姿态误差抑制方法，具体包括如下步骤：

步骤1：旋转机构带动惯性组件旋转至IMU坐标系与载体坐标系重合的位置，有

其中b表示载体坐标系，s表示IMU坐标系，

表示s系到b系转换矩阵，I表示单位阵。

步骤2：将光纤陀螺捷联惯导系统进行充分预热后，旋转机构带动惯性组件以ω进行单轴正反转停运动。采用四个转停次序为一个旋转周期的旋转方案。

所述惯性组件转动过程采用四个转停次序为一个旋转周期的转位方案，如图2所示，图中，①～④表示前4个旋转过程，A、B、C、D表示四个停留位置，x_b、y_b表示载体坐标系的水平轴，并要求旋转初始时刻IMU与载体坐标系完全重合。具体为：

次序1，IMU从A点出发顺时针转动90，到达位置C，停位时间T_r；

次序2，IMU从C点出发顺时针转动180°，到达位置B，停位时间T_r；

次序3，IMU从B点出发逆时针转动90，到达位置D，停位时间T_r；

次序4，IMU从D点出发逆时针转动180，到达位置A，停位时间T_r。

IMU按照此转动顺序循环进行。

由此得到惯性组件的旋转调制周期T为

式中，T表示一个旋转调制的调制周期，ω表示转动过程中的转动角速度，T_r表示停位时间。

旋转机构提供其实时转过的角度，得到s系到b系转换矩阵

为

$C_s^b=\left[\begin{array}{ccc}\cos \mathrm{\α}& -\sin \mathrm{\α}& 0\\ \sin \mathrm{\α}& \cos \mathrm{\α}& 0\\ 0& 0& 1\end{array}\right]---\left(2\right)$

其中，α表示旋转机构绕IMU的z轴实时转过的角度。

步骤3：依据导航参数设计原则，根据旋转周期时间设计导航参数。设计方法如下：

${\mathrm{\ω}}_0\≤{\mathrm{\ω}}_1\·\sqrt[4]{\frac{1}{k{\mathrm{\ξ}}^2}}---\left(3\right)$

式中，ω₀表示导航解算过程中固有频率，是参数设计过程中的主要设计对象；ω₁＝2π/T，表示一个旋转调制周期的调制频率；ξ为光纤陀螺捷联惯性导航系统的阻尼系数，取值范围0～1；k表示光纤陀螺捷联惯性导航系统解算姿态误差角受IMU常值偏差的影响降低倍数，可采用任意值。

依据（3）式可以得到ω_n的取值范围。在该范围内ω_n取值，按下式计算得到导航参数的最终设计结果，

其中，K₁、K₂、K_E、K_N、K_U都为导航参数，主要用于下一步的修正控制角速率计算。Ω＝0004167°/s，表示地球自转角速度；

表示载体所在位置的地理纬度；g表示重力加速度。

步骤4：实时采集光纤陀螺仪和石英加速度计测量载体运动的角速度和线加速度信息，修正控制角速率。修正控制角速度的公式为：

其中，

分别表示控制角速率

在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量，δv_j（j＝x,y）表示光纤陀螺捷联惯导系统解算速度误差；p表示复数域参变量。

步骤5：利用IMU测量信息以及修正控制角速率导航解算，得到载体姿态信息。

首先更新角速度：

${\mathrm{\ω}}_{\mathrm{ns}}^s={\mathrm{\ω}}_{\mathrm{is}}^s-{\left(C_s^n\right)}^T({\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)-{\left(C_s^n\right)}^T{\mathrm{\ω}}_c^n---\left(6\right)$

其中，i表示地心惯性系，e表示地球坐标系，s表示IMU坐标系，n表示导航坐标系，

表示s系到n系转换矩阵；·^T表示矩阵转置；（m＝n,i,e，p₁＝s,e，n，q＝s,n）表示p₁系相对m系的旋转角速度在q系投影；

为地球自转角速度在n系投影。

地球自转角速度在导航坐标系投影

更新为：

表示载体运动导致导航坐标系相对地球坐标系变化的旋转角速度在导航坐标系投影,更新过程为

其中，v_j（j＝x,y）表示光纤陀螺捷联惯导系统解算的载体沿导航坐标系ox_n轴和oy_n轴的运动速度；R表示地球半径。

然后采用更新四元数法更新转换矩阵

设载体坐标系相对导航坐标系的转动四元数Q为：

Q＝q₀+q₁i_b+q₂j_b+q₃k_b            （9）

其中，q₀、q₁、q₂和q₃为四元数的四个实数；i_b、j_b和k_b分别为载体坐标系ox_b轴、oy_b轴和oz_b轴上的单位方向向量。

四元数Q的及时修正：

$\left[\begin{array}{c}{\stackrel{\·}{q}}_0\\ {\stackrel{\·}{q}}_1\\ {\stackrel{\·}{q}}_2\\ {\stackrel{\·}{q}}_3\end{array}\right]=\frac{1}{2}\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_{\mathrm{nsx}}^s& -{\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsz}}^s\\ {\mathrm{\ω}}_{\mathrm{nsx}}^s& 0& {\mathrm{\ω}}_{\mathrm{nsz}}^s& -{\mathrm{\ω}}_{\mathrm{nsy}}^s\\ {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsz}}^s& 0& {\mathrm{\ω}}_{\mathrm{nsx}}^s\\ {\mathrm{\ω}}_{\mathrm{nsz}}^s& {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsx}}^s& 0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]---\left(10\right)$

其中，

表示IMU系相对导航坐标系的旋转角速度在IMU坐标系ox_s轴、oy_s轴、oz_s轴上的分量。

分别表示q₀、q₁、q₂、q₃的微分量。

利用求出的q₀、q₁、q₂、q₃更新转换矩阵

$C_s^n=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]---\left(11\right)$

转换矩阵

更新：

$C_b^n=C_s^n{\left(C_s^b\right)}^T---\left(12\right)$

更新载体姿态信息：

$\left\{\begin{array}{c}\mathrm{\θ}=\arcsin \left(c_{33}\right)\\ \mathrm{\φ}=\arctan (c_{32}/c_{31})\\ \mathrm{\ψ}=\arctan (c_{13}/c_{23})\end{array}\right.---\left(13\right)$

其中，θ、φ、ψ为载体纵摇角、横摇角和航向角；c_ij（i,j＝1,2,3）表示

中第i行第j列矩阵元素。

步骤6：更新载体运动速度信息，得到系统解算速度误差。具体过程为，

将加速度计测量比力通过转换矩阵

转换，得到：

$f^n=C_s^n{f}^s---\left(14\right)$

其中，fⁿ、f^s分别表示加速度计测量比力在n系和s系投影。

利用下列微分方程求解载体运动速度：

$\left[\begin{array}{c}{\stackrel{\·}{v}}_x\\ {\stackrel{\·}{v}}_y\\ {\stackrel{\·}{v}}_z\end{array}\right]=\left[\begin{array}{c}{f}_x^n\\ f_y^n\\ f_z^n\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ g\end{array}\right]+\left[\begin{array}{ccc}0& 2{\mathrm{\ω}}_{\mathrm{iez}}^n& -(2{\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n)\\ -{\mathrm{\ω}}_{\mathrm{iez}}^n& 0& 2{\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n\\ 2{\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n& -(2{\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n)& 0\end{array}\right]\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]---\left(15\right)$

其中，

分别表示加速度计测量比力fⁿ在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量。

和

分别表示地球自转角速度

在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量。分别表示由于载体运动导致导航坐标系相对地球坐标系变化的旋转角速度在导航坐标系ox_n轴、oy_n轴上的投影。v_x、v_y和v_z分别表示光纤陀螺捷联惯导系统解算的载体沿导航坐标系ox_n轴、oy_n轴和oz_n轴的运动速度，

分别表示v_x、v_y、v_z的微分量。

光纤陀螺捷联惯导系统解算的速度误差：

$\left\{\begin{array}{c}\mathrm{\δ}{v}_x=v_x-v_{\mathrm{rx}}\\ \mathrm{\δ}{v}_y=v_y-v_{\mathrm{ry}}\end{array}\right.---\left(16\right)$

其中，v_rx、v_ry分别表示外界测量载体运动速度，这里作为速度基准。

步骤7：重复步骤4到步骤6，直至光纤陀螺捷联惯导系统导航结束。利用步骤6得到的系统解算速度误差δv_x和δv_y更新修正控制角速率再利用

求取转换矩阵

最终得到载体姿态信息，得到的载体姿态信息和速度误差作为光纤陀螺捷联惯导系统的输出。

对本发明的有益效果进行验证如下：

（1）在Visual C++仿真条件下，对该方法进行仿真实验：

载体初始位置：北纬457796°，东经1266705°；

载体匀速直航运动，运动速度为v＝15m/s；

赤道半径：R＝6378393.0m；

由万有引力可得的地球表面重力加速度：g＝9.78049m/s²；

地球自转角速度：Ω＝72921158×10^-5rad/s；

常数：π＝3.1415926535；

光纤陀螺常值漂移：001°/h；

光纤陀螺白噪声误差：0005°/h；

光纤陀螺刻度因数误差：10ppm；

光纤陀螺安装误差：1×10^-3rad；

加速度计零偏：10^-4g₀；g₀表示赤道重力加速度；

加速度计白噪声误差：5×10^-5g₀；

加速度计刻度因数误差：10ppm；

加速度计安装误差：1×10^-3rad；

仿真时间：t＝48h；

采样频率：Hn＝0.01s；

IMU四位置转停方案的参数：

四个位置的停顿时间：T_r＝800s；

转动180°和90°的转动角速度：ω＝6°/s；

转动180°和90°的过程中，每一个转位中的角加（减）速度：α＝3°/s²；

导航参数：k＝50、ξ＝0707、ω₀＝00005。

利用发明所述方法，得到姿态误差比较曲线如图3和图4：图3为未使用本发明方法的捷联惯导系统解算的姿态误差曲线，图4为使用本发明方法的捷联惯导系统解算的姿态误差曲线。结果表明采用本发明方法能够较好地抑制调制状态下捷联惯导系统解算的姿态误差，提高导航精度，增强姿态信息可用性。

（2）光纤陀螺捷联惯导系统单轴转台试验

采用920E型单轴测试转台和自行研制的光纤陀螺惯导系统构建试验系统。

a.所用光纤陀螺惯导系统主要技术指标如下：

动态范围：±100°/s；

零偏稳定性：≤0005°/h；

随机游走：

标度因数非线性度：≤5ppm。

b.920E型单轴转台台主要技术指标如下：

面直径：450mm；

负载要求：重量50kg；

台体回转精度：±2″；

台体转角范围：连续无限；

位置精度：±3″；

位置分辨力：00001°；

速率范围：0005-200°/s；

速率精度：5×10^-5（360°平均）、5×10^-4（10°平均）、1×10^-2（1°平均）。

c.其他实验条件

试验时间：t＝48h；

IMU四位置转停方案的参数：

四个位置的停顿时间：T_r＝800s；

转动180°和90°的转动角速度：ω＝6°/s；

导航参数：k＝50、ξ＝0707、ω₀＝00005。

本发明实施例2得到的姿态误差比较曲线如图5和图6所示：图5为未使用本发明方法解算的姿态误差曲线，图6为使用本发明方法解算的姿态误差曲线。结果表明采用本发明方法抑制姿态误差能力较好，可以满足实际需求。

Claims (4)

1.一种调制型捷联惯性导航系统姿态误差抑制方法，其特征在于，包括以下步骤：

步骤1：旋转机构带动惯性组件旋转至惯性组件(IMU)坐标系与载体坐标系重合的位置，有

其中b表示载体坐标系，s表示IMU坐标系，表示s系到b系转换矩阵，I表示单位阵；

步骤2：将光纤陀螺捷联惯导系统进行充分预热后，旋转机构带动惯性组件以角速度ω进行单轴正反转停运动，惯性组件采用四个转停次序为一个旋转周期的旋转方案；

步骤3：根据惯性组件的旋转调制周期T设计导航参数，具体是：

首先，确定导航解算过程中固有频率ω₀的取值范围：

ω₁表示惯性组件的旋转调制周期的调制频率ω₁＝2π/T；ξ为光纤陀螺捷联惯导系统阻尼系数，取值范围0～1；k表示光纤陀螺捷联惯导系统解算姿态误差角受IMU常值偏差的影响降低倍数；

然后，确定根据下式确定导航参数K₁、K₂、K_E、K_N和K_U：

表示载体所在位置的地理纬度；g表示重力加速度；Ω＝0004167°/s，表示地球自转角速度；

步骤4：实时采集光纤陀螺仪和石英加速度计测量载体运动的角速度和线加速度，修正控制角速率，具体根据下式来修正：

其中，n表示导航坐标系，分别表示控制角速率

在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量，δv_j（j＝x,y）表示光纤陀螺捷联惯导系统解算速度误差；p表示复数域参变量；

步骤5：利用IMU测量信息以及修正的控制角速率进行导航解算，得到载体姿态信息；

步骤6：更新载体运动速度，得到光纤陀螺捷联惯导系统解算速度误差；

步骤7：重复步骤4至步骤6，直至导航结束，不断更新修正控制角速率、速度误差和载体姿态信息。

2.根据权利要求1所述的一种调制型捷联惯性导航系统姿态误差抑制方法，其特征在于，步骤2中所述的惯性组件采用四个转停次序为一个旋转周期的旋转方案，具体是：

次序1，IMU从A点出发顺时针转动90°，到达位置C，停位时间T_r；

次序2，IMU从C点出发顺时针转动180°，到达位置B，停位时间T_r；

次序3，IMU从B点出发逆时针转动90°，到达位置D，停位时间T_r；

次序4，IMU从D点出发逆时针转动180°，到达位置A，停位时间T_r；

惯性组件的旋转调制周期

3.根据权利要求1所述的一种调制型捷联惯性导航系统姿态误差抑制方法，其特征在于，所述的步骤5具体为：

首先更新角速度：

其中，i表示地心惯性系，e表示地球坐标系；

表示s系到n系转换矩阵；·^T表示矩阵转置；（m＝n,i,e，p₁＝s,e,n，q＝s,n）表示p₁系相对m系的旋转角速度在q系投影；

为地球自转角速度在n系投影，的更新公式为：

的更新公式为：其中，v_j（j＝x,y）表示光纤陀螺捷联惯导系统解算速度；R表示地球半径；

然后采用更新四元数法更新转换矩阵

设载体坐标系相对导航坐标系的转动四元数Q为：Q＝q₀+q₁i_b+q₂j_b+q₃k_b，q₀、q₁、q₂和q₃为四元数的四个实数；i_b、j_b和k_b分别为载体坐标系ox_b轴、oy_b轴和oz_b轴上的单位方向向量；

对四元数Q的及时修正： $\left[\begin{array}{c}{\stackrel{\·}{q}}_0\\ {\stackrel{\·}{q}}_1\\ {\stackrel{\·}{q}}_2\\ {\stackrel{\·}{q}}_3\end{array}\right]=\frac{1}{2}\left[\begin{array}{cccc}0& -{\mathrm{\ω}}_{\mathrm{nsx}}^s& -{\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsz}}^s\\ {\mathrm{\ω}}_{\mathrm{nsx}}^s& 0& {\mathrm{\ω}}_{\mathrm{nsz}}^s& -{\mathrm{\ω}}_{\mathrm{nsy}}^s\\ {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsz}}^s& 0& {\mathrm{\ω}}_{\mathrm{nsx}}^s\\ {\mathrm{\ω}}_{\mathrm{nsz}}^s& {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsx}}^s& 0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right],$ ${\mathrm{\ω}}_{\mathrm{nsx}}^s,{\mathrm{\ω}}_{\mathrm{nsy}}^s,{\mathrm{\ω}}_{\mathrm{nsz}}^s$ 表示IMU系相对导航坐标系的旋转角速度在IMU坐标系ox_s轴、oy_s轴、oz_s轴上的分量，

分别表示q₀、q₁、q₂、q₃的微分量；

利用求出的q₀、q₁、q₂、q₃更新转换矩阵

$C_s^n=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2(q_1{q}_2-q_0{q}_3)& 2(q_1{q}_3+q_0{q}_2)\\ 2(q_1{q}_2+q_0{q}_3)& q_0^2-q_1^2+q_2^2-q_3^2& 2(q_2{q}_3-q_0{q}_1)\\ 2(q_1{q}_3-q_0{q}_2)& 2(q_2{q}_3+q_0{q}_1)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]$

转换矩阵

更新为：

进一步得到更新的载体姿态信息： $\left\{\begin{array}{c}\mathrm{\θ}=\arcsin \left(c_{33}\right)\\ \mathrm{\φ}=\arctan (c_{32}/c_{31})\\ \mathrm{\ψ}=\arctan (c_{13}/c_{23})\end{array}\right.;$ 其中，θ、φ、ψ为载体纵摇角、横摇角和航向角；c_ij（i,j＝1,2,3）表示转换矩阵

中第i行第j列矩阵元素。

4.根据权利要求1或3所述的一种调制型捷联惯性导航系统姿态误差抑制方法，其特征在于，所述的步骤6具体是：

将加速度计测量比力通过转换矩阵

转换：

其中，fⁿ、f^s分别表示加速度计测量比力在n系和s系投影；

利用微分方程求解载体运动速度：

$\left[\begin{array}{c}{\stackrel{\·}{v}}_x\\ {\stackrel{\·}{v}}_y\\ {\stackrel{\·}{v}}_z\end{array}\right]=\left[\begin{array}{c}{f}_x^n\\ f_y^n\\ f_z^n\end{array}\right]-\left[\begin{array}{c}0\\ 0\\ g\end{array}\right]+\left[\begin{array}{ccc}0& 2{\mathrm{\ω}}_{\mathrm{iez}}^n& -(2{\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n)\\ -{\mathrm{\ω}}_{\mathrm{iez}}^n& 0& 2{\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n\\ 2{\mathrm{\ω}}_{\mathrm{iey}}^n+{\mathrm{\ω}}_{\mathrm{eny}}^n& -(2{\mathrm{\ω}}_{\mathrm{iex}}^n+{\mathrm{\ω}}_{\mathrm{enx}}^n)& 0\end{array}\right]\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]$

其中，

分别表示加速度计测量比力fⁿ在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量；

和

分别表示地球自转角速度

在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量；

分别表示导航坐标系相对地球坐标系变化的旋转角速度在导航坐标系ox_n轴、oy_n轴上的投影；v_x、v_y和v_z分别表示光纤陀螺捷联惯导系统解算的载体沿导航坐标系ox_n轴、oy_n轴和oz_n轴的运动速度，

分别表示v_x、v_y、v_z的微分量；

最后得到系统解算速度误差： $\left\{\begin{array}{c}\mathrm{\δ}{v}_x=v_x-v_{\mathrm{rx}}\\ \mathrm{\δ}{v}_y=v_y-v_{\mathrm{ry}}\end{array}\right.,$ 其中，v_rx、v_ry分别表示外界测量的载体运动速度，为速度基准。

Images