CN103090867B - 相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法 - Google Patents

Abstract

本发明公开一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法。所述方法具体步骤为：首先，确定旋转方案，其次根据载体初始位置参数建立惯导系统初始捷联转换矩阵，再次计算得到惯导系统初始转换矩阵；然后在旋转机构转动过程中根据捷联惯导系统实时解算载体位置信息以及转换矩阵，更新旋转机构下一时刻旋转角速度，直到完成调制型捷联惯导系统相对惯性系旋转的调制过程。本发明能够完全消除惯性器件各常值误差项对系统导航精度影响，特别是消除陀螺刻度因数误差和安装误差与地球转速耦合项误差，提高了惯导系统导航精度。

Description

相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法

技术领域

本发明设计的是一种基于旋转调制的惯性导航系统误差抑制方法，更确切地说，是通过相对地心惯性系的三轴旋转调制方式对光纤陀螺捷联惯性导航系统误差的一种抑制方法。

背景技术

捷联惯导系统SINS是一种全自主导航系统，它利用惯性敏感元件（陀螺仪和加速度计）测量载体相对惯性空间的线运动和角运动，并在已知初始条件下，用计算机计算出载体的速度、位置和姿态等导航参数。由于捷联惯导系统具有体积小、重量轻、隐蔽性好、制造和维护成本低、可靠性高等特点，在航海、航空、航天等许多领域都得到广泛应用。然而，根据SINS基本原理可知，SINS在导航过程中由于惯性组件常值偏差的存在而导致系统定位误差随时间增长而发散是影响系统导航精度的主要因素之一。

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

旋转调制技术是一种惯性器件偏差自补偿方法，该方法通过旋转机构带动惯性组件有规律的转动，对惯性器件常值偏差的调制来抵消该误差项对系统的影响，进而提高系统导航精度。基于旋转调制技术的捷联惯导系统根据转轴数目划分为单轴系统、双轴系统、三轴系统。目前应用较为广泛的旋转型捷联惯导系统，主要是绕惯性组件水平轴和天向轴旋转的单轴旋转和双轴旋转系统。

但是，由于地球转速的存在，无论对于单轴系统还是多轴系统，只要惯性测量单元是相对于地理坐标系旋转，总是存在光纤陀螺刻度因数误差和安装误差与地球转速的耦合项，并且利用前面所述的旋转调制方式无法完全消除该耦合项对系统的影响。因此，如何完全消除光纤陀螺刻度因数误差和安装误差与地球转速的耦合项对导航系统的精度影响具有重要意义。

CNKI库中目前已有部分与旋转调制相关的研究报道：（1）申请号为200910071733.3，名称为“基于单轴旋转的光纤陀螺捷联惯性导航系统误差抑制方法”的申请专利。该专利主要论述了相对地理系单轴连续旋转的调制型捷联惯导系统。（2）申请号为200910073241.8，名称为“基于单轴四位置转停方案的捷联系统误差抑制方法”的申请专利。该专利主要论述了一种相对地理系旋转的单轴四位置转停方案。上述两个专利申请的不足之处是无法抑制天向轴IMU常值误差对系统导航精度的影响。（3）国防科技大学博士学位论文《四频激光陀螺旋转式惯导系统研究》，该文章主要以静电陀螺的旋转调制方式为例讨论了相对地理系旋转的单轴、双轴旋转调制方法，但仍无法消除陀螺仪常值误差与地球自转角速度的耦合误差项对系统导航精度影响。

发明内容

本发明的目的在于提供一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，该方法不仅能够抑制惯导系统的发散式定位误差，还可以完全抵消惯性组件（光纤陀螺仪和加速度计）常值误差，特别是光纤陀螺刻度因数误差和安装误差与地球转速耦合项对导航精度的影响。

本发明提供的相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，旋转调制过程中，旋转机构通过三转轴同时变速转动来达到IMU相对地心惯性系调制转动的目的；根据载体位置以及转换矩阵实时更新转换矩阵和旋转机构旋转角速度ω^s，并施加给旋转机构带动IMU转动，即可完全调制IMU常值偏差对导航误差的影响。其具体步骤如下：

(1)系统导航前，设计捷联惯导系统相对地心惯性系的四位置八次序旋转方案。导航过程中，根据设计旋转方案确定旋转机构相对惯性坐标系的旋转角速度其中，b代表载体坐标系，s表示IMU坐标系（与旋转机构坐标系重合），i代表地心惯性系，表示s系相对b系旋转角速度在i系投影，即旋转方案设计的旋转机构旋转角速度。

(2)利用全球定位系统GPS确定载体初始位置参数，将它们装订至导航计算机中，建立初始时刻转换矩阵其中，n代表导航坐标系，代表i系到n系转换矩阵。

(3)将旋转机构转动至IMU系与载体系重合的位置，有其中b表示载体坐标系，表示b系到s系转换矩阵，I表示单位阵。光纤陀螺捷联惯导系统进行预热后采集光纤陀螺仪和石英加速度计（简称加速度计）输出的数据；根据加速度计输出与重力加速度关系以及陀螺仪输出与地球自转角速度关系确定载体姿态角，完成系统初始对准，建立惯导系统初始捷联转换矩阵其中，代表n系到b系的转换矩阵。

(4)根据计算得到惯导系统初始转换矩阵其中，表示i系到s系的转换矩阵；

(5)依据步骤（1）中设计的旋转方案，以及步骤（4）中得到的转换矩阵利用角速度投影计算方法得到旋转机构转动角速度旋转机构带动IMU开始转动，其中，表示s系相对b系旋转角速度在s系投影；

(6)转动过程中，根据捷联惯导系统实时解算载体位置信息以及转换矩阵通过更新步骤（4）中的转换矩阵

(7)依据步骤（1）中的旋转方案已设定的旋转角速度以及步骤（6）中得到的更新转换矩阵通过角速度投影计算方法更新旋转机构下一时刻旋转角速度；

(8)将步骤（7）中得到旋转角速度施加给旋转机构，旋转机构带动惯性组件转动；

(9)旋转过程中不断重复步骤（6）～步骤（8），更新旋转机构的旋转角速度，进而完成调制型捷联惯导系统相对惯性系旋转的调制过程。

所述的IMU采用三轴四位置八次序转动次序（定义逆时针旋转为正）为一个旋转周期的转位方案为：

次序1，IMU从A位置出发，绕oz_i轴正向旋转360°，到达位置A；次序2，IMU从A位置出发，绕oy_i轴正向旋转180°，到达位置B；次序3，IMU从B位置出发，绕oz_i轴反向旋转360°，到达位置B；次序4，IMU从B位置出发，绕oz_i轴反向旋转180°，到达位置C；次序5，IMU从C位置出发，绕oz_i轴反向旋转360°，到达位置C；次序6，IMU从C位置出发，绕oy_i轴反向旋转180°，到达位置D；次序7，IMU从D位置出发，绕oz_i轴正向旋转360°，到达位置D；次序8，IMU从D位置出发，绕oz_i轴正向旋转180°，到达位置A。IMU按照此转动顺序循环进行。

所述的相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，旋转机构相对地心惯性系坐标轴转动调制。旋转机构转速由两部分组成，

${\mathrm{\ω}}_{\mathrm{bs}}^s={\mathrm{\ω}}^s+{\mathrm{\ω}}_0^s---\left(1\right)$

其中，ω^s是为使IMU整周旋转、翻转等运动而设定的角速度值；与地球转速大小相同、方向相反的恒定角速度在s系投影，具体形式为：

${\mathrm{\ω}}_0^s=-{\mathrm{\ω}}_{\mathrm{ie}}^s=-C_i^s{\mathrm{\ω}}_{\mathrm{ie}}^i---\left(2\right)$

其中，e表示地球坐标系； 分别表示地球自转角速度在s系和i系投影；表示i系到s系的转换矩阵。

地球自转角速度在i系投影为常值形式：

${\mathrm{\ω}}_{\mathrm{ie}}^i={\left[\begin{array}{ccc}0& 0& \mathrm{\Ω}\end{array}\right]}^T---\left(3\right)$

其中，Ω表示常值地球自转角速度。

本发明与现有技术相比的优点在于：本发明打破了传统的相对地理坐标系旋转的调制方式，提出了一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其优点是不需要任何外界信息，完全自主的抵消了惯性器件常值偏差、刻度因数误差、安装误差对系统导航精度的影响，并且能够完全抑制相对地理系旋转无法消除的陀螺仪刻度因数误差和安装误差分别与地球自转角速度的耦合误差项，进而降低了系统发散式定位误差，使系统导航精度不再受器件常值偏差影响，提高导航精度。

附图说明

图1为本发明相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法流程图；

图2为附表1中相对地理系旋转的双轴四位置旋转方案四个转位示意图；

图3为附表2中相对地心惯性系的三轴四位置旋转方案四个转位示意图；

图4（a）为IMU静止状态下，由加速度计误差引起的系统定位误差曲线；

图4（b）为采用相对地理系的双轴四位置旋转方案时，由加速度计误差引起的系统定位误差曲线；

图4（c）为采用相对地心惯性系的三轴四位置旋转方案时，由加速度计误差引起的系统定位误差曲线；

图5（a）为IMU静止状态下，由光纤陀螺仪误差引起的系统定位误差曲线；

图5（b）为采用相对地理系的双轴四位置旋转方案时，由光纤陀螺仪误差引起的系统定位误差曲线；

图5（c）为采用相对地心惯性系的三轴四位置旋转方案时，由光纤陀螺仪误差引起的系统定位误差曲线；

图6(a)、6(b)分别为结合图4(a)～4(c)、图5(a)～5(c)得到IMU在三种状态下72小时引起系统加速度误差和脱落误差的比较直方图。

具体实施方式

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

本发明提供一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，流程图如图1所示，所述的误差抑制方法通过如下步骤实现：

第一步，系统导航前，设计捷联惯导系统相对地心惯性系的四位置八次序旋转方案。

IMU采用四位置八次序转动次序为一个旋转周期的旋转方案为：

次序1，IMU从位置A出发，绕oz_i轴正向旋转360°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置A；次序2，IMU从位置A出发，绕oy_i轴正向旋转180°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& \mathrm{\ω}& -\mathrm{\Ω}\end{array}\right]}^T,$ 到达位置B；次序3，IMU从位置B出发，绕oz_i轴反向旋转360°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置B；次序4，IMU从位置B出发，绕oz_i轴反向旋转180°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置C；次序5，IMU从位置C出发，绕oz_i轴反向旋转360°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置C；次序6，IMU从位置C出发，绕oy_i轴反向旋转180°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& -\mathrm{\ω}& -\mathrm{\Ω}\end{array}\right]}^T,$ 到达位置D；次序7，IMU从位置D出发，绕oz_i轴正向旋转360°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置D；次序8，IMU从位置D出发，绕oz_i轴正向旋转180°旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置A。IMU按照此转动顺序循环进行。其中，Ω＝0.00414°/s，ω＝6°/s。

上述旋转调制方案中旋转机构转速由两部分组成，

${\mathrm{\ω}}_{\mathrm{bs}}^s=C_i^s{\mathrm{\ω}}_{\mathrm{bs}}^i={\mathrm{\ω}}^s+{\mathrm{\ω}}_0^s---\left(4\right)$

其中，表示i系到s系的转换矩阵；表示s系相对于b系的旋转角速度在s系投影，即旋转机构的旋转角速度在s系投影；ω^s是为使IMU整周旋转、翻转等运动而设定的角速度值；与地球转速大小相同、方向相反的恒定角速度在s系投影，具体形式为：

${\mathrm{\ω}}_0^s=-{\mathrm{\ω}}_{\mathrm{ie}}^s=-C_i^s{\mathrm{\ω}}_{\mathrm{ie}}^i---\left(5\right)$

其中，i表示地心惯性系，e表示地球坐标系； 分别表示地球自转角速度在s系和i系投影。

地球自转角速度在i系投影为常值形式：

${\mathrm{\ω}}_{\mathrm{ie}}^i={\left[\begin{array}{ccc}0& 0& \mathrm{\Ω}\end{array}\right]}^T---\left(6\right)$

其中，Ω表示常值地球自转角速度。

对光纤陀螺仪输出误差加速度计输出误差δf_n在一个旋转周期内进行积分，

${\∫}_0^T{\mathrm{\δ\ω}}_{\mathrm{is}}^n\mathrm{dt}={\∫}_0^T{C}_s^b\mathrm{dt}\·{\mathrm{\ϵ}}^s+{\∫}_0^T{C}_s^n\mathrm{dt}\·({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g){\mathrm{\ω}}_{\mathrm{bs}}^s+{\∫}_0^T\left[C_s^n({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g)C_i^s\right]\mathrm{dt}\·{\mathrm{\ω}}_{\mathrm{ie}}^i---\left(7\right)$

${\∫}_0^T{\mathrm{\δf}}^n\mathrm{dt}={\∫}_0^T{C}_s^n\mathrm{dt}\·[{\▿}^s+({\mathrm{\δK}}_a+{\mathrm{\Δ}}_a)f^s]---\left(8\right)$

式中，T为旋转周期，δfⁿ表示光纤陀螺仪和加速度计沿导航坐标系测量误差；ε^s、δK_g、Δ_g分别为陀螺仪的常值漂移、刻度因数误差、安装误差；δK_a、Δ_a分别为加速度计的零偏、刻度因数误差、安装误差，f^s为加速度计敏感轴输入信息，为s系到n系转换矩阵。

若采用（4）式的旋转角速度，即绕惯性系旋转方式，变为

${\∫}_0^T{\mathrm{\δ\ω}}_{\mathrm{is}}^n={\∫}_0^T{C}_s^n\mathrm{dt}\·{\mathrm{\ϵ}}^s+{\∫}_0^T{C}_s^n\mathrm{dt}\·({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g){\mathrm{\ω}}^s---\left(9\right)$

将惯性测量单元IMU旋转后光纤陀螺仪生成的数据转换到导航坐标系下，得到惯性器件常值偏差的调制形式。结合（5）式、（6）式以及旋转方案，IMU正反转动八个转动次序为一个旋转周期，该旋转周期内转换矩阵在每一个旋转过程的积分结果为：

其中，s·表示sin·，c·表示cos·；λ表示任意时刻载体所在位置纬度、经度信息，T_i(i＝1,...,8)分别表示每个转动次序对应的转动时间，且T＝T₁+…+T₈。

将一个完整转动周期T内的IMU转动过程的转换矩阵积分结果相加，得到导航系下转换矩阵的累积作用结果。

${\∫}_0^{T_1}{C}_{s(A\→A)}^n\mathrm{dt}+{\∫}_0^{T_2}{C}_{s(A\→B)}^n\mathrm{dt}+{\∫}_0^{T_3}{C}_{s(B\→B)}^n\mathrm{dt}+{\∫}_0^{T_4}{C}_{s(B\→C)}^n\mathrm{dt}+---\left(11\right)$

${\∫}_0^{T_5}{C}_{s(C\→C)}^n\mathrm{dt}+{\∫}_0^{T_6}{C}_{s(C\→D)}^n\mathrm{dt}+{\∫}_0^{T_7}{C}_{s(D\→D)}^n\mathrm{dt}+{\∫}_0^{T_8}{C}_{s(D\→A)}^n\mathrm{dt}=0$

将（11）式代入（8）式和（9）式中，得到

$\left\{\begin{array}{c}{\∫}_0^T{\mathrm{\δf}}^n\mathrm{dt}=0\\ {\∫}_0^T{\mathrm{\δ\ω}}_{\mathrm{is}}^n\mathrm{dt}=0\end{array}\right.---\left(12\right)$

至此，导航坐标系上各轴陀螺仪和加速度计常值偏差得到调制，其在一个旋转调制周期内沿导航系的投影积分结果为零，即IMU常值偏差不影响捷联惯导系统的导航精度。

第二步，利用全球定位系统GPS确定载体初始位置参数，将它们装订至导航计算机中，并建立初始时刻转换矩阵

所述的载体初始位置参数包括初始时刻载体所在位置的纬度和经度信息。其中，λ₀表示初始时刻载体所在位置纬度、经度信息；，i代表地心惯性系，n代表导航坐标系。

第三步，将旋转机构转动至IMU坐标系与载体坐标系重合位置，有其中b代表载体坐标系，s表示惯性组件坐标系，表示b系到s系转换矩阵，I表示单位阵。光纤陀螺捷联惯导系统进行预热后采集光纤陀螺仪和石英加速度计输出的数据；根据加速度计输出与重力加速度关系以及光纤陀螺仪输出与地球自转角速度关系确定载体姿态角，完成系统初始对准，建立惯导系统初始捷联矩阵

$C_n^b=\left[\begin{array}{ccc}\cos {\mathrm{\φ}}_{y0}\cos {\mathrm{\φ}}_{z0}-\sin {\mathrm{\φ}}_{x0}\sin {\mathrm{\φ}}_{y0}\sin {\mathrm{\φ}}_{z0}& {\sin \mathrm{\φ}}_{z0}{\cos \mathrm{\φ}}_{y0}+{\sin \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{y0}{\cos \mathrm{\φ}}_{z0}& -{\cos \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{y0}\\ -{\cos \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{z0}& {\cos \mathrm{\φ}}_{x0}\cos {\mathrm{\φ}}_{z0}& {\sin \mathrm{\φ}}_{x0}\\ {\sin \mathrm{\φ}}_{y0}{\cos \mathrm{\φ}}_{z0}+\mathrm{si}{\mathrm{n\φ}}_{z0}{\sin \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}& {\sin \mathrm{\φ}}_{z0}{\sin \mathrm{\φ}}_{y0}+{\cos \mathrm{\φ}}_{z0}{\sin \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}& {\cos \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}\end{array}\right]---\left(14\right)$

其中，φ_x0、φ_y0、φ_z0分别表示初始时刻载体俯仰角、横滚角、航向角，n代表导航坐标系，b代表载体坐标系。

第四步，建立惯导系统初始转换矩阵

$C_i^s=C_b^s{C}_n^b{C}_i^n---\left(15\right)$

其中，n表示导航坐标系，这里采用当地地理坐标系； 分别表示b系到s系、n系到b系、i系到n系转换矩阵。

第五步，依据第一步中设计的旋转方案，以及第四步中得到的转换矩阵计算初始时刻旋转机构旋转角速度，旋转机构带动IMU开始转动。具体计算方法为：

${\mathrm{\ω}}_{\mathrm{bs}}^s=C_i^s{\mathrm{\ω}}_{\mathrm{bs}}^i---\left(16\right)$

第六步：转动过程中，根据捷联惯导系统实时解算载体位置信息以及转换矩阵通过更新步骤（4）中的转换矩阵

实时采集光纤陀螺仪和石英加速度计测量载体运动的线速度和角速度信息，导航解算得到导航信息；

通过旋转调制状态下采集的陀螺仪数据，更新转换矩阵具体为：

角速度更新：

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

其中，e表示地球坐标系，n表示导航坐标系，这里采用当地地理坐标系；表示s系到n系转换矩阵；（m＝n,i,e，p＝s，e,n，q＝s,n）表示p系相对m系旋转角速度在q系投影。

四元数姿态矩阵更新：

设任意时刻载体坐标系相对平台坐标系的转动四元数为：

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

其中，Q为四元数；q₀、q₁、q₂、q₃为四元数的四个实数；i_b、j_b、k_b分别表示IMU坐标系ox_s轴、oy_s轴、oz_s轴上的单位方向向量。

根据k时刻载体坐标系相对平台坐标系的转动四元数q₀(k)、q₁(k)、q₂(k)、q₃(k)，求取k时刻转动四元数的变化率为：

$\left[\begin{array}{c}{\stackrel{\·}{q}}_0\left(k\right)\\ {\stackrel{\·}{q}}_1\left(k\right)\\ {\stackrel{\·}{q}}_2\left(k\right)\\ {\stackrel{\·}{q}}_3\left(k\right)\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{nsx}}^s\\ {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsz}}^s& 0& {\mathrm{\ω}}_{\mathrm{msx}}^s\\ {\mathrm{\ω}}_{\mathrm{nsz}}^s& {\mathrm{\ω}}_{\mathrm{nsy}}^s& -{\mathrm{\ω}}_{\mathrm{nsx}}^s& 0\end{array}\right]\left[\begin{array}{c}{q}_0\left(k\right)\\ q_1\left(k\right)\\ q_2\left(k\right)\\ q_3\left(k\right)\end{array}\right]---\left(19\right)$

其中， 分别表示旋转机构相对导航系的运动角速度在IMU坐标系ox_s轴、oy_s轴、oz_s轴上的分量。 分别表示q₀(k)、q₁(k)、q₂(k)、q₃(k)的变化率。

在k+1时刻载体的转动四元数具体为：

$\left\{\begin{array}{c}{q}_0(k+1)=q_0\left(k\right)+{\stackrel{\·}{q}}_0\left(k\right)\\ q_1(k+1)=q_1\left(k\right)+{\stackrel{\·}{q}}_1\left(k\right)\\ q_2(k+1)=q_2\left(k\right)+{\stackrel{\·}{q}}_2\left(k\right)\\ q_3(k+1)=q_3\left(k\right)+{\stackrel{\·}{q}}_3\left(k\right)\end{array}\right.---\left(20\right)$

利用得到的q₀(k+1)、q₁(k+1)、q₂(k+1)、q₃(k+1)更新转换矩阵

$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(21\right)$

其中，（21）式中省略了四元数中的(k+1)部分。

更新载体姿态信息，具体为：

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

其中，c_ij（i＝1,2,3，j＝1,2,3）表示中第i行第j列矩阵元素；φ_x、φ、φ_z表示载体纵摇角、横摇角、航向角。

将加速度计沿IMU坐标系测量的比力信息，通过转换矩阵进行投影转换：

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

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

根据k时刻的载体东向水平速度v_x(k)、北向水平速度v_y(k)和天向速度v_z(k)，求取k时刻载体速度变化率为：

$\left[\begin{array}{c}{\stackrel{\·}{v}}_x\left(k\right)\\ {\stackrel{\·}{v}}_y\left(k\right)\\ {\stackrel{\·}{v}}_z\left(k\right)\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\left(k\right)\\ v_y\left(k\right)\\ v_z\left(k\right)\end{array}\right]---\left(24\right)$

其中，下角标x、y、z分别表示在导航系oxn轴、oy_n轴、oz_n轴上的分量；v_x(k)、v_y(k)、v_z(k)表示解算载体速度； 表示v_x(k)、v_y(k)、v_z(k)的变化率，即载体运动加速度； 表示加速度计测量比力；g为重力加速度。 表示地球自转角速度在导航系的投影； 表示由于载体运动而导致导航系相对地球系变化的旋转角速度在导航系投影。

在k+1时刻载体速度和位置分别为：

$\left\{\begin{array}{c}{v}_x(k+1)=v_x\left(k\right)+{\stackrel{\·}{v}}_x\left(k\right)\\ v_y(k+1)=v_y\left(k\right)+{\stackrel{\·}{v}}_y\left(k\right)\\ v_z(k+1)=v_z\left(k\right)+{\stackrel{\·}{v}}_z\left(k\right)\end{array}\right.---\left(25\right)$

其中，R表示地球半径；λ分别表示计算载体所在地理位置的纬度和经度信息，当k＝1时，v_x(1)、v_y(1)、v_z(1)为第二步中利用GPS获得的载体初始速度，λ(1)为第二步中利用GPS获得的载体初始位置。

转动过程中，根据捷联惯导系统实时解算载体位置以及转换矩阵更新转换矩阵

1）利用实时解算载体位置信息更新转换矩阵建立任意时刻i系到n系转换矩阵

2）更新转换矩阵

$C_i^s={\left(C_s^n\right)}^T{C}_i^n---\left(28\right)$

第七步，计算旋转机构下一时刻旋转角速度：

依据旋转方案已设定的旋转角速度以及前面更新的转换矩阵利用角速度投影计算方法确定旋转机构下一时刻旋转角速度；

第八步，将第九步中得到旋转角速度施加给旋转机构，旋转机构带动惯性组件转动；

第九步，旋转过程中不断重复第六步～第八步，实时更新旋转机构的旋转角速率，完成调制型捷联惯导系统相对惯性系旋转的调制过程。

下面通过实施例对本发明的有益效果进行说明验证：

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

旋转机构带动惯性组件分别处于以下三种状态进行导航解算：（1）静止状态，（2）按照相对地理系旋转的双轴四位置旋转方案转动，如图2。（3）按照相对地心惯性系的三轴四位置旋转方案转动，如图3。其中，x_i、y_i、z_i表示i系的三个坐标轴。其中，相对地理系旋转的双轴四位置旋转方案如附表1，相对地心惯性系的三轴四位置旋转方案如附表2。

表1相对地理系旋转的双轴四位置旋转方案

表2相对地心惯性系的三轴四位置旋转方案

载体初始位置：北纬45.7796°，东经126.6705°；

载体真实姿态角：φ_x＝0°，φ_y＝0°，φ_z＝30°；

赤道半径：R＝6378393.0m；

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

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

常数：π＝3.1415926535；

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

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

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

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

加速度计零偏：10^-4g₀；

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

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

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

仿真时间：t＝72h；

采样频率：Hn＝0.01s；

不同旋转方案下，由加速度计误差和光纤陀螺误差引起的系统定位误差曲线分别如图4(a)~4(c)、图5(a)~5(c)所示。其中图4(a)和图5（a）为IMU静止状态，图4(b)和图5（b）为相对地理系的双轴四位置旋转方案，图4(c)和图5（c）为相对地心惯性系的三轴四位置旋转方案。图6是结合图4(a)～4(c)、图5(a)~5(c)得到的IMU在三种状态下72小时引起系统定位误差比较直方图。

通过对图4(a)～4(c)、图5(a)~5(c)、图6对比可以看出，本发明中提出的相对地心惯性系的旋转方案不仅可以有效的减小定位误差，相比传统相对地理系的旋转方式可以进一步提高捷联惯导系统的定位精度。

Claims (7)

1.一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于，包括如下步骤：

第一步，系统导航前，设计捷联惯导系统相对地心惯性系的四位置八次序旋转方案；导航过程中，根据设计旋转方案确定旋转机构相对惯性坐标系的旋转角速度其中，b代表载体坐标系，s表示IMU坐标系，i代表地心惯性系，表示s系相对b系旋转角速度在i系投影，即旋转方案设计的旋转机构旋转角速度；

第二步，利用全球定位系统GPS确定载体初始位置参数，将它们装订至导航计算机中，建立初始时刻转换矩阵其中，n代表导航坐标系，代表i系到n系转换矩阵：

所述的载体初始位置参数包括初始时刻载体所在位置的纬度和经度信息，其中，λ₀表示初始时刻载体所在位置纬度、经度信息；i代表地心惯性系，n代表导航坐标系；

第三步，将旋转机构转动至IMU系与载体系重合的位置，有其中b表示载体坐标系，表示b系到s系转换矩阵，I表示单位阵；光纤陀螺捷联惯导系统进行预热后采集光纤陀螺仪和石英加速度计输出的数据；根据加速度计输出与重力加速度关系以及陀螺仪输出与地球自转角速度关系确定载体姿态角，完成系统初始对准，建立惯导系统初始捷联转换矩阵其中，代表n系到b系的转换矩阵：

$C_n^b=\left[\begin{array}{ccc}{\cos \mathrm{\φ}}_{y0}\cos {\mathrm{\φ}}_{z0}-{\sin \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{y0}{\sin \mathrm{\φ}}_{z0}& {\sin \mathrm{\φ}}_{z0}{\cos \mathrm{\φ}}_{y0}+{\sin \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{y0}{\cos \mathrm{\φ}}_{z0}& -{\cos \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{y0}\\ -{\cos \mathrm{\φ}}_{x0}{\sin \mathrm{\φ}}_{z0}& {\cos \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{z0}& {\sin \mathrm{\φ}}_{x0}\\ {\sin \mathrm{\φ}}_{y0}{\cos \mathrm{\φ}}_{z0}+{\sin \mathrm{\φ}}_{z0}{\sin \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}& {\sin \mathrm{\φ}}_{z0}{\sin \mathrm{\φ}}_{y0}+{\cos \mathrm{\φ}}_{z0}{\sin \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}& {\cos \mathrm{\φ}}_{x0}{\cos \mathrm{\φ}}_{y0}\end{array}\right]---\left(2\right)$

其中，φ_x0、φ_y0、φ_z0分别表示初始时刻载体俯仰角、横滚角、航向角，n代表导航坐标系，b代表载体坐标系；

第四步，根据计算得到惯导系统初始转换矩阵其中，表示i系到s系的转换矩阵；

第五步，依据第一步中设计的旋转方案，以及第四步中得到的转换矩阵利用角速度投影计算方法得到旋转机构旋转角速度旋转机构带动IMU开始转动，其中，表示s系相对b系旋转角速度在s系投影；

第六步，转动过程中，根据捷联惯导系统实时解算载体位置信息以及转换矩阵通过更新第四步中的转换矩阵

第七步，依据第一步中的旋转方案已设定的旋转角速度以及第六步中得到的更新转换矩阵通过角速度投影计算方法更新旋转机构下一时刻旋转角速度；

第八步，将第七步中得到的旋转角速度施加给旋转机构，旋转机构带动惯性组件转动；

第九步，旋转过程中不断重复第六步～第八步，更新旋转机构的旋转角速度，进而完成调制型捷联惯导系统相对惯性系旋转的调制过程。

2.根据权利要求1所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：所述的四位置八次序旋转方案为：

次序1，IMU从位置A出发，绕oz_i轴正向旋转360°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置A；次序2，IMU从位置A出发，绕oy_i轴正向旋转180°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& \mathrm{\ω}& -\mathrm{\Ω}\end{array}\right]}^T,$ 到达位置B；次序3，IMU从位置B出发，绕oz_i轴反向旋转360°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置B；次序4，IMU从位置B出发，绕oz_i轴反向旋转180°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置C；次序5，IMU从位置C出发，绕oz_i轴反向旋转360°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}-\mathrm{\ω}\end{array}\right]}^T,$ 到达位置C；次序6，IMU从位置C出发，绕oy_i轴反向旋转180°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& -\mathrm{\ω}& -\mathrm{\Ω}\end{array}\right]}^T,$ 到达位置D；次序7，IMU从位置D出发，绕oz_i轴正向旋转360°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置D；次序8，IMU从位置D出发，绕oz_i轴正向旋转180°，旋转角速度 ${\mathrm{\ω}}_{\mathrm{bs}}^i={\left[\begin{array}{ccc}0& 0& -\mathrm{\Ω}+\mathrm{\ω}\end{array}\right]}^T,$ 到达位置A；IMU按照此转动顺序循环进行；其中，Ω＝0.00414°/s，ω＝6°/s。

3.根据权利要求1所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：所述的旋转机构旋转角速度由两部分组成，

${\mathrm{\ω}}_{\mathrm{bs}}^s={\mathrm{\ω}}^s+{\mathrm{\ω}}_0^s---\left(3\right)$

其中，ω^s是为使IMU整周旋转、翻转运动而设定的角速度值；与地球转速大小相同、方向相反的恒定角速度在s系投影，具体形式为：

${\mathrm{\ω}}_0^s=-{\mathrm{\ω}}_{\mathrm{ie}}^s=-C_i^s{\mathrm{\ω}}_{\mathrm{ie}}^i---\left(4\right)$

其中，e表示地球坐标系；分别表示地球自转角速度在s系和i系投影；表示i系到s系的转换矩阵；

地球自转角速度在i系投影为常值形式：

${\mathrm{\ω}}_{\mathrm{ie}}^i={\left[\begin{array}{ccc}0& 0& \mathrm{\Ω}\end{array}\right]}^T---\left(5\right)$

其中，Ω表示地球自转角速度常值。

4.根据权利要求3所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：所述的调制过程，具体为：

对光纤陀螺仪输出误差加速度计输出误差δfⁿ在一个旋转周期内进行积分，

${\∫}_0^T\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{is}}^n\mathrm{dt}={\∫}_0^T{C}_s^n\mathrm{dt}\·{\mathrm{\ϵ}}^s+{\∫}_0^T{C}_s^n\mathrm{dt}\·({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g){\mathrm{\ω}}_{\mathrm{bs}}^s+{\∫}_0^T\left[C_s^n({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g)C_i^s\right]\mathrm{dt}\·{\mathrm{\ω}}_{\mathrm{ie}}^i---\left(6\right)$

${\∫}_0^T\mathrm{\δ}{f}^n\mathrm{dt}={\∫}_0^T{C}_s^n\mathrm{dt}\·[{\▿}^s+({\mathrm{\δK}}_a+{\mathrm{\Δ}}_a)f^s]---\left(7\right)$

式中，T为旋转周期，δfⁿ表示光纤陀螺仪和加速度计沿导航坐标系测量误差；ε^s、δK_g、Δ_g分别为陀螺仪的常值漂移、刻度因数误差、安装误差；▽^s、δK_a、Δ_a分别为加速度计的零偏、刻度因数误差、安装误差，f^s为加速度计敏感轴输入信息，为s系到n系转换矩阵；

若采用(3)式的旋转角速度，即绕惯性系旋转方式，变为

${\∫}_0^T\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{is}}^n\mathrm{dt}={\∫}_0^T{C}_s^n\mathrm{dt}\·{\mathrm{\ϵ}}^s+{\∫}_0^T{C}_s^n\mathrm{dt}\·({\mathrm{\δK}}_g+{\mathrm{\Δ}}_g){\mathrm{\ω}}^s---\left(8\right)$

将惯性测量单元IMU旋转后光纤陀螺仪生成的数据转换到导航坐标系下，得到惯性器件常值偏差的调制形式；结合(4)式、(5)式以及旋转方案，IMU正反转动八个转动次序为一个旋转周期，该旋转周期内转换矩阵在每一个旋转过程的积分结果为：

其中，s·表示sin·，c·表示cos·；λ表示任意时刻载体所在位置纬度、经度信息，T_i分别表示每个转动次序对应的转动时间，i＝1,…,8，且T＝T₁+…+T₈；

将一个完整转动周期T内的IMU转动过程的转换矩阵积分结果相加，得到导航系下转换矩阵的累积作用结果：

$\begin{array}{c}{\∫}_0^{T_1}{C}_{s(A\→A)}^n\mathrm{dt}+{\∫}_0^{T_2}{C}_{s(A\→B)}^n\mathrm{dt}+{\∫}_0^{T_s}{C}_{s(B\→B)}^n\mathrm{dt}+{\∫}_0^{T_4}{C}_{s(B\→C)}^n\mathrm{dt}+\\ {\∫}_0^{T_5}{C}_{s(C\→C)}^n\mathrm{dt}+{\∫}_0^{T_6}{C}_{s(C\→D)}^n\mathrm{dt}+{\∫}_0^{T_7}{C}_{s(D\→D)}^n\mathrm{dt}+{\∫}_0^{T_8}{C}_{s(D\→A)}^n\mathrm{dt}=0\end{array}---\left(10\right)$

将(10)式代入(7)式和(8)式中，得到：

$\left\{\begin{array}{c}{\∫}_0^T{\mathrm{\δf}}^n\mathrm{dt}=0\\ {\∫}_0^T\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{is}}^n\mathrm{dt}=0\end{array}\right.---\left(11\right)$

至此，导航坐标系上各轴陀螺仪和加速度计常值偏差得到调制。

5.根据权利要求1所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：第六步中更新转换矩阵是指在转动过程中，根据捷联惯导系统实时解算载体位置信息以及转换矩阵更新转换矩阵具体为：

1)利用实时解算得到载体位置信息更新转换矩阵建立任意时刻i系到n系转换矩阵

λ表示任意时刻载体所在位置纬度、经度信息；

2)更新转换矩阵

$C_i^s={\left(C_s^n\right)}^T{C}_i^n.---\left(13\right)$

6.根据权利要求5所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：更新转换矩阵具体为：

角速度更新：

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

其中，e表示地球坐标系，n表示导航坐标系，这里采用当地地理坐标系；表示s系到n系转换矩阵；表示p系相对m系旋转角速度在q系投影，m＝n,i,e，p＝s,e,n，q＝s,n；

四元数姿态矩阵更新：

设任意时刻载体坐标系相对平台坐标系的转动四元数为：

Q＝q₀+q₁i_b+q₂j_b+q₃k_b   (15)

其中，Q为四元数；q₀、q₁、q₂、q₃为四元数的四个实数；i_b、j_b、k_b分别表示IMU坐标系ox_s轴、oy_s轴、oz_s轴上的单位方向向量；

根据k时刻载体坐标系相对平台坐标系的转动四元数q₀(k)、q₁(k)、q₂(k)、q₃(k)，求取k时刻转动四元数的变化率为：

$\left[\begin{array}{c}{\stackrel{\·}{q}}_0\left(k\right)\\ {\stackrel{\·}{q}}_1\left(k\right)\\ {\stackrel{\·}{q}}_2\left(k\right)\\ {\stackrel{\·}{q}}_3\left(k\right)\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{nsx}}^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\left(k\right)\\ q_1\left(k\right)\\ q_2\left(k\right)\\ q_3\left(k\right)\end{array}\right]---\left(16\right)$

其中，分别表示旋转机构相对导航系的运动角速度在IMU坐标系ox_s轴、oy_s轴、oz_s轴上的分量；分别表示q₀(k)、q₁(k)、q₂(k)、q₃(k)的变化率；

在k+1时刻载体的转动四元数具体为：

$\left\{\begin{array}{c}{q}_0(k+1)=q_0\left(k\right)+{\stackrel{\·}{q}}_0\left(k\right)\\ q_1(k+1)=q_1\left(k\right)+{\stackrel{\·}{q}}_1\left(k\right)\\ q_2(k+1)=q_2\left(k\right)+{\stackrel{\·}{q}}_2\left(k\right)\\ q_3(k+1)=q_3\left(k\right)+{\stackrel{\·}{q}}_3\left(k\right)\end{array}\right.---\left(17\right)$

利用得到的q₀(k+1)、q₁(k+1)、q₂(k+1)、q₃(k+1)更新转换矩阵

$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(18\right)$

其中，(18)式中省略了四元数中的(k+1)部分。

7.根据权利要求1所述的一种相对地心惯性系旋转的光纤陀螺捷联惯性导航系统误差抑制方法，其特征在于：解算载体位置信息，具体为：

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

其中，c_ij表示中第i行第j列矩阵元素，i＝1,2,3，j＝1,2,3；φ_x、φ_y、φ_z表示载体纵摇角、横摇角、航向角；

将加速度计沿IMU坐标系测量的比力信息，通过转换矩阵进行投影转换：

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

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

根据k时刻的载体东向水平速度v_x(k)、北向水平速度v_y(k)和天向速度v_z(k)，求取k时刻载体速度变化率为：

$\left[\begin{array}{c}{\stackrel{\·}{v}}_x\left(k\right)\\ {\stackrel{\·}{v}}_y\left(k\right)\\ {\stackrel{\·}{v}}_z\left(k\right)\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\left(k\right)\\ v_y\left(k\right)\\ v_z\left(k\right)\end{array}\right]---\left(21\right)$

其中，下角标x、y、z分别表示在导航系ox_n轴、oy_n轴、oz_n轴上的分量；v_x(k)、v_y(k)、v_z(k)表示解算载体速度；表示v_x(k)、v_y(k)、v_z(k)的变化率，即载体运动加速度； 表示加速度计测量比力；g为重力加速度；表示地球自转角速度在导航系的投影； 表示由于载体运动而导致导航系相对地球系变化的旋转角速度在导航系投影；

在k+1时刻载体速度和位置分别为：

$\left\{\begin{array}{c}{v}_x(k+1)=v_x\left(k\right)+{\stackrel{\·}{v}}_x\left(k\right)\\ v_y(k+1)=v_y\left(k\right)+{\stackrel{\·}{v}}_y\left(k\right)\\ v_z(k+1)=v_z\left(k\right)+{\stackrel{\·}{v}}_z\left(k\right)\end{array}\right.---\left(22\right)$

其中，R表示地球半径；λ分别表示计算载体所在地理位置的纬度和经度信息，当k＝1时，v_x(1)、v_y(1)、v_z(1)为第二步中利用GPS获得的载体初始速度，λ(1)为第二步中利用GPS获得的载体初始位置。