CN103076025B - 一种基于双解算程序的光纤陀螺常值误差标定方法 - Google Patents

Abstract

本发明公开了一种基于双解算程序的光纤陀螺常值误差标定方法。该方法利用旋转机构带动惯性组件分别旋转至各轴陀螺沿导航系z轴正向和反向的六个位置，每个位置停留过程中，将一组惯性组件的测量值同时作为两组导航解算程序的输入值，其中两组程序中导航参数设定值不同。再将两组解算姿态信息进一步耦合运算，估算各轴陀螺常值漂移和刻度因数误差。本发明则是提出了在导航计算机中以一组惯性组件的测量值作为输入，同时进行两组导航程序解算，利用两组解算结果进一步估算器件误差，不需要任何外界基准信息，计算量小、简单易操作，并且两组导航解算程序具有相同的输入信息，不存在任何安装偏差和时间延迟，输出信息具有相关性。

Description

一种基于双解算程序的光纤陀螺常值误差标定方法

技术领域

本发明涉及捷联惯性导航系统的惯性组件误差测试领域，具体是一种基于双解算程序的光纤陀螺仪常值误差确定方法。

背景技术

捷联惯性导航系统SINS作为一种不需要任何外界信息，能够连续输出载体速度、位置、姿态信息的全自主导航系统被广泛应用于航空、航天、航海等军事和民用领域。其将惯性组件（InertialMeasurementUnit,简称IMU，包括陀螺仪和加速度计）直接安装在运载体上，利用陀螺仪敏感的载体角运动信息解算得到载体系到导航系的转换矩阵，再将加速度计测量沿载体系的加速度信息投影到导航系上，投影后的加速度一次积分得到载体运动速度，二次积分得到载体位置信息。此外，利用转换矩阵与载体姿态关系可以得到载体姿态信息。

陀螺仪作为捷联惯导系统的核心器件，由于各种误差的存在，导致其输出信息与敏感角速度总是存在一定偏差，导致转换矩阵不准确，影响系统导航精度。标定技术作为一种惯性组件误差项测量技术，主要是通过器件误差与系统解算信息的数学关系，利用包含器件误差的导航信息来推算惯性器件各误差项。标定技术主要分为分立式标定和系统级标定。分立式标定是直接以光纤陀螺仪输出为观测量，用最小二乘法标定其系数，该方法数据处理复杂，标定精度低；系统级标定则是利用导航误差作为观测量，通过滤波等手段来估计陀螺仪误差系数，但该方法需要外界参考基准信息以得到导航误差，若基准信息不准确则会导致观测量存在误差，进而影响标定精度。因此，如何快速、有效、高精度的估算陀螺仪误差系数是一项非常重要的课题。

《光电工程》2008年第35卷第1期中由刘百奇等人撰写的《光纤陀螺IMU的六位置旋转现场标定新方法》，该文章将IMU在六个位置上进行十二次旋转，根据光纤陀螺IMU的误差模型建立42个非线性输入输出方程求解出陀螺标度因数、陀螺常值漂移等误差系数。该方法不仅建立数学模型多，而且计算量大；《中国惯性技术学报》2005年第13卷第4期中由尚捷等人撰写的《捷联惯导现场最优标定方法研究》，该文章研究了一种基于虚拟噪声的现场最优标定方法，虽然该方法结构简单、易于实现，但文中并没有提及虚拟噪声设定方法，不同惯性组件中引入同一组虚拟噪声是否适用有待考虑；申请号为201010523436.0的中国发明专利《一种光纤陀螺常值漂移标定方法》中，针对旋转调制型捷联惯导系统，提出了一种利用单轴旋转机构变换IMU姿态，在建立惯性器件偏差与失准角估值的二元一次方程组基础上，求解陀螺水平常值漂移的估算方法。该方法只能估算常值漂移，并且若系统中存在未知误差源，导致建立惯性器件偏差与失准角估值模型不准确，则会影响估算精度。

发明内容

本发明的目的是为了解决上述问题，提供一种基于双解算程序的光纤陀螺仪常值误差确定方法。为了估算三轴陀螺仪的常值漂移和刻度因数误差，利用旋转机构带动惯性组件分别停留在三轴陀螺仪处于朝上和朝下的六个位置。每个位置以一组惯性组件的测量值为导航计算机的输入信息，导航计算机中同时进行两组参数设定值不同的导航程序并行解算，解算姿态值做差，计算得到陀螺仪常值漂移和刻度因数误差。

本发明提供的一种基于双解算程序的光纤陀螺常值误差标定方法，包括如下步骤：

步骤1：将捷联惯导系统中的惯性组件，陀螺仪和加速度计安装在三轴转台（旋转机构）上；

步骤2：三轴转台带动惯性组件旋转至与当地东北天地理坐标系一致的位置，此时方位轴陀螺仪沿地理系的oz_n轴方向，定义该位置为位置A，停留时间大于两小时；

步骤3：在位置A的停留时间内，以陀螺仪和加速度计敏感的角速度和线加速度信息作为导航计算中导航解算程序1的输入信息进行导航解算，其中导航参数设定为ξ₁、ω_n1，得到载体姿态θ_A1、φ_A1、ψ_A1。其中，θ、φ、ψ分别表示纵摇角、横摇角、方位角，角标A表示在位置A解算结果，角标1表示导航解算程序1解算结果；

步骤4：在位置A的停留时间内，与步骤3同步进行。以步骤3中同一组惯性组件测量值作为导航计算机中导航解算程序2的输入信息进行导航解算，其中导航参数设定为ξ₂、ω_n2，得到载体姿态θ_A2、φ_A2、ψ_A2。角标2表示导航解算程序2解算结果；

步骤5：利用步骤3与步骤4得到的两组姿态信息，耦合运算，得到陀螺仪沿oz_n轴输出误差。其中，n表示导航系，oz_n表示沿n系oz轴方向。计算方式如下，

式中，表示位置A处陀螺仪沿oz_n轴输出误差，Ω表示地球自转角速度，表示当地纬度。

步骤6：相对于位置A，旋转机构带动惯性组件绕ox_s轴旋转180°（定义顺时针旋转为正），使oz_s轴陀螺仪沿导航系-oz_n方向，定义该位置为位置B，停留时间大于两小时。其中，s表示IMU坐标系，ox_s表示沿s系ox轴方向；

步骤7：在位置B的停留时间内，重复步骤3-步骤5。即在该位置以一组惯性组件的输出同时进行两组导航程序解算，解算姿态耦合运算，得到捷联惯导旋转至位置B处陀螺仪沿oz_n轴输出误差。计算方式为

其中，角标B表示在位置B处两组导航程序解算姿态信息。

步骤8：将步骤5与步骤7解算结果按下式计算，得到oz_s陀螺仪的常值漂移和刻度因数误差；

其中，ε_z0表示oz_s轴陀螺常值漂移，δK_gz表示oz_s陀螺仪刻度因数误差。

至此，通过步骤2-步骤8，利用旋转机构带动IMU分别旋转至方位轴陀螺仪沿oz_n和-oz_n方向，进而估算出方位轴陀螺仪的常值漂移和刻度因数误差项。下面利用相同的旋转方式，估算其余两轴陀螺仪的常值漂移和刻度因数误差。

步骤9：相对于位置B，旋转机构带动惯性组件绕oy_s轴旋转-90°，使ox_s轴陀螺仪沿导航系oz_n方向，定义该位置为位置C，停留时间大于两小时。停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置C处陀螺仪沿oz_n轴输出误差。计算方式为

其中，角标C表示在位置C处两组导航程序解算姿态信息。

步骤10：相对于位置C，旋转机构带动惯性组件绕oz_s轴旋转180°，使ox_s轴陀螺仪沿导航系-oz_n方向，定义该位置为位置D，停留时间大于两小时。停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置D处陀螺仪沿oz_n轴输出误差。计算方式为

其中，下角标D表示在位置D处两组导航程序解算姿态信息。

步骤11：将步骤9与步骤10估算结果按下式计算得到ox_s轴陀螺仪的常值漂移和刻度因数误差；

其中，ε_x0表示ox_s轴陀螺常值漂移，δK_gx表示ox_s轴陀螺刻度因数误差。至此，通过步骤9-步骤11估算出ox_s轴陀螺仪的常值漂移和刻度因数误差。

步骤12：相对于位置D，旋转机构带动惯性组件绕oz_s轴旋转-90°，使oy_s轴陀螺仪沿导航系oz_n方向，定义该位置为位置E，停留时间大于两小时。停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置E处陀螺仪沿oz_n轴输出误差。计算方式为

其中，角标E表示在位置E处两组导航程序解算姿态信息。

步骤13：相对于位置E，旋转机构带动惯性组件绕oz_s轴旋转180°，使oy_s轴陀螺仪沿导航系-oz_n方向，定义该位置为位置F，停留时间大于两小时。停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置F处陀螺仪沿oz_n轴输出误差。计算方式为

其中，角标F表示在位置F处两组导航程序解算姿态信息。

步骤14：将步骤13与步骤14估算结果按下式计算，得到oy_s轴陀螺仪的常值漂移和刻度因数误差；

其中，ε_y0表示oy_s轴陀螺常值漂移，δK_gy表示oy_s轴陀螺刻度因数误差。

至此，通过步骤12-步骤14估算出oy_s轴陀螺仪的常值漂移和刻度因数误差。且通过以上所有步骤估算出三轴陀螺仪的误差项。

本发明与现有技术相比的优点在于：

本发明提出了一种基于双解算程序的光纤陀螺常值误差标定方法。该方法利用旋转机构带动惯性组件分别旋转至各轴陀螺沿导航系z轴正向和反向的六个位置，每个位置停留过程中，将一组惯性组件的测量值同时作为两组导航解算程序的输入值，其中两组程序中导航参数设定值不同。再将两组解算姿态信息进一步耦合运算，估算各轴陀螺常值漂移和刻度因数误差。本发明属于系统级标定方法，传统的系统级标定方法是利用惯性器件测量值作为导航计算机的输入，在导航计算机中通过一组导航解算程序得到的导航误差作为观测量，通过滤波等手段来估计陀螺仪误差系数。该方法需要外界参考基准信息以得到导航误差，若基准信息不准确则会导致观测量存在误差，进而影响标定精度。本发明则是提出了在导航计算机中以一组惯性组件的测量值作为输入，同时进行两组导航程序解算，利用两组解算结果进一步估算器件误差。其优点在于：（1）不需要任何外界基准信息；（2）计算量小、简单易操作；（3）两组导航解算程序具有相同的输入信息，不存在任何安装偏差和时间延迟，输出信息具有相关性。

附图说明

图1为本发明确定陀螺仪误差参数的标定流程图；

图2为本发明确定陀螺仪误差参数的标定原理图；

图3为利用本发明进行室内转台试验，将标定结果补偿前后系统解算姿态误差对比曲线；

图4为利用本发明进行室内转台试验，将标定结果补偿前后系统解算速度误差对比曲线；

图5为利用本发明进行室内转台试验，将标定结果补偿前后系统解算定位误差对比曲线。

具体实施方式

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

如图1所示，本发明提供的一种基于双解算程序的光纤陀螺常值误差标定方法，具体包括如下步骤：

步骤1：将捷联惯导系统中的惯性组件，陀螺仪和加速度计安装在三轴转台（旋转机构）上；

步骤2：三轴转台带动惯性组件旋转至与当地东北天地理坐标系一致的位置，设该位置为位置A，停留时间大于两小时。在位置A处，陀螺仪沿导航系z轴的输出误差为：

${\mathrm{\ϵ}}_{\mathrm{zA}}^n={\mathrm{\ϵ}}_{z0}+{\mathrm{\δK}}_{\mathrm{gz}}{\mathrm{\ω}}_z---\left(10\right)$

其中，表示在位置A处陀螺仪沿导航系z轴的输出误差；ε_z0、δK_gz分别表示方位轴陀螺仪的常值漂移和刻度因数误差；ω_z表示方位轴陀螺仪敏感角速度信息，Ω表示地球自转角速度，表示载体所在地理纬度。

步骤3：在位置A的停留时间内，通过导航解算程序1进行导航解算，利用惯性组件测量信息解算载体姿态信息，如图2所示，具体过程如下：

首先更新角速度：

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

其中，i表示地心惯性系，e表示地球坐标系，s表示IMU坐标系，n表示导航坐标系，表示位置A处s系到n系转换矩阵；·^T表示矩阵转置；为控制角速率在n系上的投影；表示IMU相对导航系旋转角速度在IMU系投影；表示由于载体运动导致导航系变化角速度在导航系投影；为地球自转角速度在n系投影。

地球自转角速度在导航系投影更新：

由于标定过程中，载体没有线运动，因此惯导系统解算速度就是系统解算速度误差，即v_x＝δv_x，v_y＝δv_y。其中，v_j、δv_j（j＝x,y）分别表示惯导系统解算速度和速度误差。

更新过程为：

其中，R表示地球半径。

控制角速率在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量和更新为：

其中，k₁＝k₂＝2ξ₁ω_n1、 k₁、k₂、k_E、k_N、k_U为在导航解算程序1的罗经参数，ξ₁、ω_n1表示导航参数，人为设定，取值范围分别为ξ₁∈(0,1)、ω_n1∈(0,1)。s表示复数域参变量。g表示重力加速度

采用更新四元数法更新捷联矩阵

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

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

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

四元数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(16\right)$

其中，表示IMU系相对导航坐标系的旋转角速度在载体坐标系ox_s轴、oy_s轴、oz_s轴上的分量。分别表示q₀、q₁、q₂、q₃的微分量。

通过（16）式，利用陀螺仪测量值间接计算得到的求解微分方程，得到四元数q₀、q₁、q₂、q₃更新结果。更新捷联矩阵

$C_{s\left(A\right)}^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(17\right)$

更新载体姿态信息：

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

其中，θ_A1、φ_A1、ψ_A1为导航程序1解算载体的纵摇角、横滚角、航向角；c_ij（i，j＝1,2,3）表示中第i行第j列矩阵元素。

利用加速度计测量比力通过矩阵转换：

$f^n=C_{s\left(A\right)}^n{f}^s---\left(19\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(20\right)$

其中，分别表示加速度计测量比力在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量；g为重力加速度。和分别表示地球自转角速度在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量。分别表示由于载体运动导致导航坐标系相对地球坐标系变化的旋转角速度在导航坐标系ox_n轴、oy_n轴上的投影。分别表示v_x、v_y、v_z的微分量。

步骤4：同样在位置A的停留时间内，利用与步骤3相同的计算方法，通过导航解算程序2进行导航解算，其中导航参数设置不同，具体过程如下：

首先更新角速度：

${\mathrm{\ω}}_{\mathrm{ns}}^{\′s}={\mathrm{\ω}}_{\mathrm{is}}^s-{\left(C_{s\left(A\right)}^{\′n}\right)}^T({\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)-{\left(C_{s\left(A\right)}^{\′n}\right)}^T{\mathrm{\ω}}_c^{\′n}---\left(21\right)$

其中，表示位置A处导航解算程序2中解算s系到n系转换矩阵；为导航解算程序2中解算控制角速率在n系上的投影；为导航解算程序2中IMU相对导航系旋转角速度在IMU系投影的计算值。

（21）式中，更新过程与步骤3中相同。控制角速率更新过程与步骤3中有所区别，更新过程为：

其中，

k₁′＝k₂′＝2ξ₂ω_n2

其中，ξ₂、ω_n2为导航解算程序2中的导航参数；k′₁、k₂′、k′_E、k′_N、k′_U为在导航解算程序2的罗经参数。

更新捷联矩阵修正四元数：

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

其中，表示导航解算程序2中，IMU系相对导航坐标系的旋转角速度在载体坐标系投影的计算值。q′_k为导航解算程序2中解算四元数，为q′_k的微分量，k＝1,2,3,4。

通过（24）式，利用陀螺仪测量值间接计算得到的求解微分方程，得到四元数q₀′、q₁′、q₂′、q₃′更新结果。

更新捷联矩阵

$C_{s\left(A\right)}^{\′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(25\right)$

更新载体姿态信息：

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

其中，θ_A2、φ_A2、ψ_A2为导航程序2解算载体的纵摇角、横滚角、航向角；c′_ij（i，j＝1,2,3）表示中第i行第j列矩阵元素。

步骤5：利用步骤3与步骤4得到的两组姿态信息，进行耦合运算，计算方式如下：

式中，表示位置A处陀螺仪沿oz_n轴输出误差，Ω表示地球自转角速度，表示位置A的地理纬度。

结合（10）式可知，利用（27）式计算得到的项包含ε_z0和δK_gz。

步骤6：相对于位置A，旋转机构带动惯性组件绕ox_s轴旋转180°，其中，顺时针旋转为正，使oz_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置B，停留时间大于两小时。其中，s表示IMU坐标系，ox_s表示沿s系ox轴方向。

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{\mathrm{zB}}^n=-{\mathrm{\ϵ}}_{z0}+\mathrm{\δ}{K}_{\mathrm{gz}}{\mathrm{\ω}}_z---\left(28\right)$

其中，表示在位置B处陀螺仪沿导航系z轴的输出误差。

步骤7：在位置B的停留时间内，重复步骤3-步骤5。

IMU停留在位置B时，由于oz_s轴陀螺仪沿导航系-oz_n方向，因此相对于位置A，方位轴陀螺仪的测量输出值不同，即方位轴陀螺常值误差与刻度因数误差沿导航系输出的体现形式也不同。因此在位置B处利用重复步骤3-步骤5的过程中，虽然计算过程完全相同，计算结果也不相同。

IMU停留在位置B处以一组惯性组件的输出同时进行两组导航程序解算，解算姿态耦合运算，得到捷联惯导旋转至位置B处陀螺仪沿oz_n轴输出误差为：

其中，下角标B1和B2分别表示在位置B处导航计算机中导航解算程序1、导航解算程序2解算姿态信息。

结合（28）式可知，利用（29）式计算得到的包含-ε_z0和δK_gz项。

步骤8：根据步骤5与步骤7得到两组陀螺仪沿oz_n轴输出误差和结合（10）式、（28）式，得到方位陀螺仪的常值漂移和刻度因数误差；

至此，通过步骤2-步骤8，估算出方位轴陀螺仪的常值漂移和刻度因数误差。

下面利用相同的方式，旋转机构带动IMU转动，将纵摇轴陀螺仪，即ox_s轴陀螺仪旋转至沿oz_n轴方向的位置C和-oz_n轴方向的位置D，估算该陀螺仪的刻度因数误差和常值误差。

步骤9：相对于位置B，旋转机构带动惯性组件绕oy_s轴旋转-90°，使ox_s轴陀螺仪沿导航系oz_n方向，设该位置为位置C，停留时间大于两小时。停留过程中，重复步骤3-步骤5。

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{\mathrm{zC}}^n={\mathrm{\ϵ}}_{x0}+{\mathrm{\δK}}_{\mathrm{gx}}{\mathrm{\ω}}_z---\left(31\right)$

其中，表示在位置C处陀螺仪沿导航系z轴的输出误差；ε_xs0、δK_gx分别表示纵摇轴陀螺仪的常值漂移和刻度因数误差。

由于IMU停留在位置C时，ox_s轴陀螺仪沿导航系oz_n方向，因此利用步骤3-步骤5中的计算方式得到陀螺仪沿导航系oz_n轴的输出误差只包含由ox_s轴陀螺仪误差项ε_xs0和δK_gx引起的输出误差。

由此得到停留过程中，IMU旋转至位置C处陀螺仪沿oz_n轴输出误差为

其中，下角标C1和C2分别表示在位置C处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息。

步骤10：相对于位置C，旋转机构带动惯性组件绕oz_s轴旋转180°，使ox_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置D，停留时间大于两小时。

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{\mathrm{zD}}^n={-\mathrm{\ϵ}}_{x0}+{\mathrm{\δK}}_{\mathrm{gx}}{\mathrm{\ω}}_z---\left(33\right)$

其中，表示在位置D处陀螺仪沿导航系z轴的输出误差。

IMU停留在位置D时，ox_s轴陀螺仪沿-oz_n方向，利用步骤3-步骤5的计算方式得到IMU沿导航系z轴的输出误差只包含由ox_s轴陀螺仪引起的误差项-ε_x0和δK_gxω_z，且相对于位置C，ox_s轴陀螺仪常值漂移沿oz_n方向投影相反。

停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置D处陀螺仪沿oz_n轴输出误差计算方式为

其中，下角标D1和D2分别表示在位置D处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息。

步骤11：根据步骤9与步骤10得到的两组陀螺仪沿oz_n轴输出误差和结合（31）和（33）式，得到ox_s陀螺仪的常值漂移和刻度因数误差；

其中，ε_x0表示ox_s轴陀螺常值漂移，δK_gx表示ox_s轴陀螺刻度因数误差。

至此，通过步骤9-步骤11，估算出ox_s轴陀螺仪的常值漂移和刻度因数误差。下面利用相同的方法，通过将横摇轴陀螺仪旋转至沿oz_n轴方向的位置E和-oz_n轴方向的位置F，估算该陀螺仪的刻度因数误差和常值误差；

步骤12：相对于位置D，旋转机构带动惯性组件绕oz_s轴旋转-90°，使oy_s轴陀螺仪沿导航系oz_n方向，设该位置为位置E，停留时间大于两小时。

IMU置于位置E时，oy_s轴陀螺仪，即横摇轴陀螺仪沿导航系oz_n方向，因此陀螺仪沿oz_n方向的输出误差形式为

${\mathrm{\ϵ}}_{\mathrm{zE}}^n={\mathrm{\ϵ}}_{y0}+{\mathrm{\δK}}_{\mathrm{gy}}{\mathrm{\ω}}_z---\left(36\right)$

其中，表示在位置E处陀螺仪沿导航系z轴的输出误差；ε_y0、δK_gv分别表示横摇轴陀螺仪的常值漂移和刻度因数误差。

停留过程中，重复步骤3-步骤5得到捷联惯导旋转至位置E处陀螺仪沿oz_n轴输出误差。计算方式为

其中，下角标E1和E2分别表示在位置E处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息。

（36）、（37）式中可以看出，IMU置于位置E处，利用一组惯性组件测量值、两组导航程序解算值耦合计算得到沿导航系z轴陀螺仪输出误差包含ε_y0和δK_gy项。

步骤13：相对于位置E，旋转机构带动惯性组件绕oz_s轴旋转180°，使oy_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置F，停留时间大于两小时。

IMU置于位置F时，oy_s轴陀螺仪沿导航系-oz_n方向，陀螺仪沿oz_n方向的输出误差为

${\mathrm{\ϵ}}_{\mathrm{zF}}^n=-{\mathrm{\ϵ}}_{y0}+{\mathrm{\δK}}_{\mathrm{gy}}{\mathrm{\ω}}_z---\left(38\right)$

其中，表示在位置F处陀螺仪沿导航系z轴的输出误差。

停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置F处陀螺仪沿oz_n轴输出误差为

其中，下角标F1和F2分别表示在位置F处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息。（36）、（37）式中可以看出，IMU置于位置F处，利用一组惯性组件测量值、两组导航程序解算值耦合计算得到沿导航系z轴陀螺仪输出误差包含-ε_y0和δK_gy项。

步骤14：根据步骤12与步骤13得到的两组陀螺仪沿oz_n轴输出误差和结合（36）式和（38）式得到oy_s陀螺仪的常值漂移和刻度因数误差；

其中，ε_y0表示oy_s轴陀螺常值漂移，δK_gy表示o_ys轴陀螺刻度因数误差。

至此，通过步骤12-步骤14，估算出oy_s轴陀螺仪的常值漂移和刻度因数误差。

通过上述步骤，得到陀螺仪方位轴的常值漂移和刻度因数误差、ox_s轴的常值漂移和刻度因数误差、oy_s轴的常值漂移和刻度因数误差，在完成本法发明的标定过程基础上，在捷联惯导系统的导航过程中，利用估算的三轴陀螺仪常值漂移和刻度因数误差对其测量输出信息进行补偿，以提高陀螺仪测量载体角运动速度的精度，进而抑制捷联惯导系统中由陀螺误差为各导航信息带来的影响。

实施例：

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

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

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

赤道半径：R＝6378393.0m；

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

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

常数：π＝3.1415926535；

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

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

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

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

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

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

导航参数：ξ₁＝0.707，ω_n1＝0.02，ξ₂＝0.707，ω_n2＝0.002；

采样频率：Hn＝0.01s；

利用发明所述方法，得到光纤陀螺仪误差的标定结果如表1，结果表明本发明能够较好估计陀螺常值漂移和刻度因数误差。

表1光纤陀螺仪误差的标定结果

（2）光纤陀螺捷联惯导系统室内转台标定试验

采用光纤陀螺捷联惯导系统进行的室内转台标定试验，三轴转台带动惯性组件旋转到本发明中设计的六个位置进行标定试验，每个位置停留时间2小时。标定结束后将标定结果补偿给器件输出值，进行72小时无阻尼静态试验，比较补偿前后系统解算各导航信息误差来验证标定结果的准确性。表2为试验标定结果。图3、图4、图5分别为陀螺仪误差项标定结果补偿前后系统解算姿态误差、速度误差、定位误差曲线。

表2光纤陀螺仪误差的标定结果

根据图3-图5，光纤陀螺仪的标定结果补偿后能够有效的降低系统解算各导航信息误差，即本发明提出的基于双程序解算的陀螺仪六位置标定方法能够有效的估算陀螺各误差参数项。

Claims (1)

1.一种基于双解算程序的光纤陀螺常值误差标定方法，具体包括如下步骤：

步骤1：将捷联惯导系统中的惯性组件安装在三轴转台上，所述惯性组件包括陀螺仪和加速度计；

步骤2：三轴转台带动惯性组件旋转至与当地东北天地理坐标系一致的位置，设该位置为位置A，停留时间大于两小时；在位置A处，陀螺仪沿导航系z轴的输出误差为：

${\mathrm{\ϵ}}_{zA}^n={\mathrm{\ϵ}}_{z0}+{\mathrm{\δK}}_{gz}{\mathrm{\ω}}_z---\left(1\right)$

其中，表示在位置A处陀螺仪沿导航系z轴的输出误差；ε_z0、δK_gz分别表示方位轴陀螺仪的常值漂移和刻度因数误差；ω_z表示方位轴陀螺仪敏感角速度信息，Ω表示地球自转角速度，表示载体所在地理纬度；

步骤3：在位置A的停留时间内，通过导航解算程序1进行导航解算，利用惯性组件测量信息解算载体姿态信息，具体过程如下：

首先更新角速度：

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

其中，i表示地心惯性系，e表示地球坐标系，s表示IMU坐标系，n表示导航坐标系，表示位置A处s系到n系转换矩阵；ω_is ^s表示陀螺仪测量输出；T表示矩阵转置；为控制角速率在n系上的投影；表示IMU相对导航系旋转角速度在IMU系投影；表示由于载体运动导致导航系变化角速度在导航系投影；为地球自转角速度在n系投影；

地球自转角速度在导航系投影更新：

其中：v_x＝δv_x，v_y＝δv_y，v_j、δv_j分别表示惯导系统解算速度和速度误差，j＝x,y；

更新过程为：

其中，R表示地球半径；

控制角速率在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量和更新为：

${\mathrm{\ω}}_{cx}^n=\frac{k_N}{s+k_1}s\·{\mathrm{\δv}}_x$

${\mathrm{\ω}}_{cz}^n=\frac{k_U}{(s+k_1)(s+k_2)}s\·{\mathrm{\δv}}_x$

其中，k₁＝k₂＝2ξ₁ω_n1、 k₁、k₂、k_E、k_N、k_U为在导航解算程序1的罗经参数，ξ₁、ω_n1表示导航参数，取值范围分别为ξ₁∈(0,1)、ω_n1∈(0,1)；s表示复数域参变量；g表示重力加速度；

采用更新四元数法更新捷联矩阵

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

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

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

四元数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{\ω}}_{nsx}^s& -{\mathrm{\ω}}_{nsy}^s& -{\mathrm{\ω}}_{nsz}^s\\ {\mathrm{\ω}}_{nsx}^s& 0& {\mathrm{\ω}}_{nsz}^s& -{\mathrm{\ω}}_{nsy}^s\\ {\mathrm{\ω}}_{nsy}^s& -{\mathrm{\ω}}_{nsz}^s& 0& {\mathrm{\ω}}_{nsx}^s\\ {\mathrm{\ω}}_{nsz}^s& {\mathrm{\ω}}_{nsy}^s& -{\mathrm{\ω}}_{nsx}^s& 0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]---\left(7\right)$

其中，表示IMU系相对导航坐标系的旋转角速度在载体坐标系ox_s轴、oy_s轴、oz_s轴上的分量；分别表示q₀、q₁、q₂、q₃的微分量；

通过(7)式，利用陀螺仪测量值间接计算得到的求解微分方程，得到四元数q₀、q₁、q₂、q₃更新结果；更新捷联矩阵

$C_{s\left(A\right)}^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\left(q_1{q}_2+q_0{q}_3\right)& 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(8\right)$

更新载体姿态信息：

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

其中，θ_A1、φ_A1、ψ_A1为导航程序1解算载体的纵摇角、横滚角、航向角；c_ij表示中第i行第j列矩阵元素，i,j＝1,2,3；

利用加速度计测量比力通过矩阵转换：

$f^n=C_{s\left(A\right)}^n{f}^s---\left(10\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{\ω}}_{iez}^n& -(2{\mathrm{\ω}}_{iey}^n+{\mathrm{\ω}}_{eny}^n)\\ -{\mathrm{\ω}}_{iez}^n& 0& 2{\mathrm{\ω}}_{iex}^n+{\mathrm{\ω}}_{enx}^n\\ 2{\mathrm{\ω}}_{iey}^n+{\mathrm{\ω}}_{eny}^n& -(2{\mathrm{\ω}}_{iex}^n+{\mathrm{\ω}}_{enx}^n)& 0\end{array}\right]\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]$

其中，分别表示加速度计测量比力在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量；g为重力加速度；和分别表示地球自转角速度在导航坐标系ox_n轴、oy_n轴、oz_n轴上的分量；分别表示由于载体运动导致导航坐标系相对地球坐标系变化的旋转角速度在导航坐标系ox_n轴、oy_n轴上的投影；分别表示v_x、v_y、v_z的微分量；

步骤4：同样在位置A的停留时间内，利用与步骤3相同的计算方法，通过导航解算程序2进行导航解算，其中导航参数设置不同，其中：

角速度更新为：

${\mathrm{\ω}}_{ns}^{\′s}={\mathrm{\ω}}_{is}^s-{\left(C_{s\left(A\right)}^{\′n}\right)}^T({\mathrm{\ω}}_{ie}^n+{\mathrm{\ω}}_{en}^n)-{\left(C_{s\left(A\right)}^{\′n}\right)}^T{\mathrm{\ω}}_c^{\′n}---\left(11\right)$

其中，表示位置A处导航解算程序2中解算s系到n系转换矩阵；为导航解算程序2中解算控制角速率在n系上的投影；为导航解算程序2中IMU相对导航系旋转角速度在IMU系投影的计算值；

(11)式中，更新过程与步骤3中相同；控制角速率更新过程与步骤3中有所区别，更新过程为：

${\mathrm{\ω}}_{cx}^{\′n}=\frac{k_N^{\′}}{s+k_1^{\′}}s\·{\mathrm{\δv}}_x$

${\mathrm{\ω}}_{cz}^{\′n}=\frac{k_U^{\′}}{(s+k_1^{\′})(s+k_2^{\′})}s\·{\mathrm{\δv}}_x$

其中，

k′₁＝k′₂＝2ξ₂ω_n2

$k_E^{\′}=k_N^{\′}={\mathrm{\ω}}_{n2}^2(1+{\mathrm{\ξ}}_2^2)/g---\left(13\right)$

其中，ξ₂、ω_n2为导航解算程序2中的导航参数；k′₁、k′₂、k′_E、k′_N、k′_U为在导航解算程序2的罗经参数；

更新捷联矩阵修正四元数：

$\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{\ω}}_{nsx}^{\′s}& -{\mathrm{\ω}}_{nsy}^{\′s}& -{\mathrm{\ω}}_{nsz}^{\′s}\\ {\mathrm{\ω}}_{nsx}^{\′s}& 0& {\mathrm{\ω}}_{nsz}^{\′s}& -{\mathrm{\ω}}_{nsy}^{\′s}\\ {\mathrm{\ω}}_{nsy}^{\′s}& -{\mathrm{\ω}}_{nsz}^{\′s}& 0& {\mathrm{\ω}}_{nsx}^{\′s}\\ {\mathrm{\ω}}_{nsz}^{\′s}& {\mathrm{\ω}}_{nsy}^{\′s}& -{\mathrm{\ω}}_{nsx}^{\′s}& 0\end{array}\right]\left[\begin{array}{c}{q}_0^{\′}\\ q_1^{\′}\\ q_2^{\′}\\ q_3^{\′}\end{array}\right]---\left(14\right)$

其中，表示导航解算程序2中，IMU系相对导航坐标系的旋转角速度在载体坐标系投影的计算值；q′_k为导航解算程序2中解算四元数，为q′_k的微分量，k＝1,2,3,4；

通过(14)式，利用陀螺仪测量值间接计算得到的求解微分方程，得到四元数q′₀、q′₁、q′₂、q′₃更新结果；

更新捷联矩阵

$C_{s\left(A\right)}^{\′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\left(q_1^{\′}{q}_2^{\′}+q_0^{\′}{q}_3^{\′}\right)& 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\left(q_2^{\′}{q}_3^{\′}+q_0^{\′}{q}_1^{\′}\right)& q_0^{\′2}-q_1^{\′2}-q_2^{\′2}+q_3^{\′2}\end{array}\right]---\left(15\right)$

更新载体姿态信息：

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

其中，θ_A2、φ_A2、ψ_A2为导航程序2解算载体的纵摇角、横滚角、航向角；c′_ij表示中第i行第j列矩阵元素；

步骤5：利用步骤3与步骤4得到的两组姿态信息，进行耦合运算，计算方式如下：

式中，表示位置A处陀螺仪沿oz_n轴输出误差，Ω表示地球自转角速度，表示位置A的地理纬度；

结合(1)式可知，利用(17)式计算得到的项包含ε_z0和δK_gz；

步骤6：相对于位置A，旋转机构带动惯性组件绕ox_s轴旋转180°，其中，顺时针旋转为正，使oz_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置B，停留时间大于两小时；其中，s表示IMU坐标系，ox_s表示沿s系ox轴方向；

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{zB}^n=-{\mathrm{\ϵ}}_{z0}+{\mathrm{\δK}}_{gz}{\mathrm{\ω}}_z---\left(18\right)$

其中，表示在位置B处陀螺仪沿导航系z轴的输出误差；

步骤7：在位置B的停留时间内，重复步骤3-步骤5；

IMU停留在位置B处以一组惯性组件的输出同时进行两组导航程序解算，解算姿态耦合运算，得到捷联惯导旋转至位置B处陀螺仪沿oz_n轴输出误差为：

其中，下角标B1和B2分别表示在位置B处导航计算机中导航解算程序1、导航解算程序2解算姿态信息；

结合(18)式可知，利用(19)式计算得到的包含-ε_z0和δK_gz项；

步骤8：根据步骤5与步骤7得到两组陀螺仪沿oz_n轴输出误差和结合(1)式、(18)式，得到方位陀螺仪的常值漂移和刻度因数误差；

至此，通过步骤2-步骤8，估算出方位轴陀螺仪的常值漂移和刻度因数误差；

步骤9：相对于位置B，旋转机构带动惯性组件绕oy_s轴旋转-90°，使ox_s轴陀螺仪沿导航系oz_n方向，设该位置为位置C，停留时间大于两小时；停留过程中，重复步骤3-步骤5；

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{zC}^n={\mathrm{\ϵ}}_{x0}+{\mathrm{\δK}}_{gx}{\mathrm{\ω}}_z---\left(21\right)$

其中，表示在位置C处陀螺仪沿导航系z轴的输出误差；ε_x0、δK_gx分别表示纵摇轴陀螺仪的常值漂移和刻度因数误差；

由此得到停留过程中，IMU旋转至位置C处陀螺仪沿oz_n轴输出误差为

其中，下角标C1和C2分别表示在位置C处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息；

步骤10：相对于位置C，旋转机构带动惯性组件绕oz_s轴旋转180°，使ox_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置D，停留时间大于两小时；

此时陀螺仪沿导航系z轴的输出误差为

${\mathrm{\ϵ}}_{zD}^n=-{\mathrm{\ϵ}}_{x0}+{\mathrm{\δK}}_{gx}{\mathrm{\ω}}_z---\left(23\right)$

其中，表示在位置D处陀螺仪沿导航系z轴的输出误差；

停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置D处陀螺仪沿oz_n轴输出误差计算方式为

其中，下角标D1和D2分别表示在位置D处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息；

步骤11：根据步骤9与步骤10得到的两组陀螺仪沿oz_n轴输出误差和结合(21)和(23)式，得到ox_s陀螺仪的常值漂移和刻度因数误差；

其中，ε_x0表示ox_s轴陀螺常值漂移，δK_gx表示ox_s轴陀螺刻度因数误差；

通过步骤9-步骤11，估算出ox_s轴陀螺仪的常值漂移和刻度因数误差；

步骤12：相对于位置D，旋转机构带动惯性组件绕oz_s轴旋转-90°，使oy_s轴陀螺仪沿导航系oz_n方向，设该位置为位置E，停留时间大于两小时；

IMU置于位置E时，oy_s轴陀螺仪，即横摇轴陀螺仪沿导航系oz_n方向，因此陀螺仪沿oz_n方向的输出误差形式为

${\mathrm{\ϵ}}_{zE}^n={\mathrm{\ϵ}}_{y0}+{\mathrm{\δK}}_{gy}{\mathrm{\ω}}_z---\left(26\right)$

其中，表示在位置E处陀螺仪沿导航系z轴的输出误差；ε_y0、δK_gy分别表示横摇轴陀螺仪的常值漂移和刻度因数误差；

停留过程中，重复步骤3-步骤5得到捷联惯导旋转至位置E处陀螺仪沿oz_n轴输出误差；计算方式为

其中，下角标E1和E2分别表示在位置E处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息；

(26)、(27)式中可以看出，IMU置于位置E处，利用一组惯性组件测量值、两组导航程序解算值耦合计算得到沿导航系z轴陀螺仪输出误差包含ε_y0和δK_gy项；

步骤13：相对于位置E，旋转机构带动惯性组件绕oz_s轴旋转180°，使oy_s轴陀螺仪沿导航系-oz_n方向，设该位置为位置F，停留时间大于两小时；

IMU置于位置F时，oy_s轴陀螺仪沿导航系-oz_n方向，陀螺仪沿oz_n方向的输出误差为

${\mathrm{\ϵ}}_{zF}^n=-{\mathrm{\ϵ}}_{y0}+{\mathrm{\δK}}_{gy}{\mathrm{\ω}}_z---\left(28\right)$

其中，表示在位置F处陀螺仪沿导航系z轴的输出误差；

停留过程中，重复步骤3-步骤5，得到捷联惯导旋转至位置F处陀螺仪沿oz_n轴输出误差为

其中，下角标F1和F2分别表示在位置F处利用导航计算机中导航解算程序1、导航解算程序2解算的姿态信息；(26)、(27)式中可以看出，IMU置于位置F处，利用一组惯性组件测量值、两组导航程序解算值耦合计算得到沿导航系z轴陀螺仪输出误差包含-ε_y0和δK_gy项；

步骤14：根据步骤12与步骤13得到的两组陀螺仪沿oz_n轴输出误差和结合(26)式和(28)式得到oy_s陀螺仪的常值漂移和刻度因数误差；

其中，ε_y0表示oy_s轴陀螺常值漂移，δK_gy表示oy_s轴陀螺刻度因数误差；

通过步骤12-步骤14，估算出oy_s轴陀螺仪的常值漂移和刻度因数误差；

通过上述步骤，得到陀螺仪方位轴的常值漂移和刻度因数误差、ox_s轴的常值漂移和刻度因数误差、oy_s轴的常值漂移和刻度因数误差。