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

摘要

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

主权项

1.一种基于双解算程序的光纤陀螺常值误差标定方法，具体包括如下步骤：步骤1：将捷联惯导系统中的惯性组件，陀螺仪和加速度计安装在三轴转台上；步骤2：三轴转台带动惯性组件旋转至与当地东北天地理坐标系一致的位置，设该位置为位置A，停留时间大于两小时；在位置A处，陀螺仪沿导航系z轴的输出误差为：${ϵ}_{\mathrm{zA}}^n={ϵ}_{z0}+{\mathrm{\delta K}}_{\mathrm{gz}}{\omega }_z---\left(1\right)$其中，表示在位置A处陀螺仪沿导航系z轴的输出误差；ε_z0、δK_gz分别表示方位轴陀螺仪的常值漂移和刻度因数误差；ω_z表示方位轴陀螺仪敏感角速度信息，Ω表示地球自转角速度，表示载体所在地理纬度；步骤3：在位置A的停留时间内，通过导航解算程序1进行导航解算，利用惯性组件测量信息解算载体姿态信息，具体过程如下：首先更新角速度：${\omega }_{\mathrm{ns}}^s={\omega }_{\mathrm{is}}^s-{\left(C_{s\left(A\right)}^n\right)}^T({\omega }_{\mathrm{ie}}^n+{\omega }_{\mathrm{en}}^n)-{\left(C_{s\left(A\right)}^n\right)}^T{\omega }_c^n---\left(2\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             （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& {-\omega }_{\mathrm{nsx}}^s& {-\omega }_{\mathrm{nsy}}^s& {-\omega }_{\mathrm{nsz}}^s\\ {\omega }_{\mathrm{nsx}}^s& 0& {\omega }_{\mathrm{nsz}}^s& {-\omega }_{\mathrm{nsy}}^s\\ {\omega }_{\mathrm{nsy}}^s& {-\omega }_{\mathrm{nsz}}^s& 0& {\omega }_{\mathrm{nsx}}^s\\ {\omega }_{\mathrm{nsz}}^s& {\omega }_{\mathrm{nsy}}^s& -{\omega }_{\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(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(8\right)$更新载体姿态信息：$\left\{\begin{array}{c}{\theta }_{A1}=\arcsin \left(c_{33}\right)\\ {\phi }_{A1}=\arctan (c_{32}/c_{31})\\ {\psi }_{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\omega }_{\mathrm{iez}}^n& -({2\omega }_{\mathrm{iey}}^n+{\omega }_{\mathrm{eny}}^n)\\ {-\omega }_{\mathrm{iez}}^n& 0& {2\omega }_{\mathrm{iex}}^n+{\omega }_{\mathrm{enx}}^n\\ {2\omega }_{\mathrm{iey}}^n+{\omega }_{\mathrm{eny}}^n& -({2\omega }_{\mathrm{iex}}^n+{\omega }_{\mathrm{enx}}^n)& 0\end{array}\right]\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]---\left(11\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进行导航解算，其中导航参数设置不同，其中：角速度更新为：${\omega }_{\mathrm{ns}}^{\prime s}={\omega }_{\mathrm{is}}^s-{\left(C_{s\left(A\right)}^{\prime n}\right)}^T({\omega }_{\mathrm{ie}}^n+{\omega }_{\mathrm{en}}^n)-{\left(C_{s\left(A\right)}^{\prime n}\right)}^T{\omega }_c^{\prime n}---\left(11\right)$其中，表示位置A处导航解算程序2中解算s系到n系转换矩阵；为导航解算程序2中解算控制角速率在n系上的投影；为导航解算程序2中IMU相对导航系旋转角速度在IMU系投影的计算值；（11）式中，更新过程与步骤3中相同；控制角速率更新过程与步骤3中有所区别，更新过程为：其中，k₁′＝k₂′＝2ξ₂ω_n2其中，ξ₂、ω_n2为导航解算程序2中的导航参数；k′₁、k₂′、k′_E、k′_N、k′_U为在导航解算程序2的罗经参数；更新捷联矩阵修正四元数：$\left[\begin{array}{c}{\stackrel{·}{q}}_0^{\prime }\\ {\stackrel{·}{q}}_1^{\prime }\\ {\stackrel{·}{q}}_2^{\prime }\\ {\stackrel{·}{q}}_3^{\prime }\end{array}\right]=\frac{1}{2}\left[\begin{array}{cccc}0& {-\omega }_{\mathrm{nsx}}^{\prime s}& {-\omega }_{\mathrm{nsy}}^{\prime s}& {-\omega }_{\mathrm{nsz}}^{\prime s}\\ {\omega }_{\mathrm{nsx}}^{\prime s}& 0& {\omega }_{\mathrm{nsz}}^{\prime s}& {-\omega }_{\mathrm{nsy}}^{\prime s}\\ {\omega }_{\mathrm{nsy}}^{\prime s}& {-\omega }_{\mathrm{nsz}}^{\prime s}& 0& {\omega }_{\mathrm{nsx}}^{\prime s}\\ {\omega }_{\mathrm{nsz}}^{\prime s}& {\omega }_{\mathrm{nsy}}^{\prime s}& {-\omega }_{\mathrm{nsx}}^{\prime s}& 0\end{array}\right]\left[\begin{array}{c}{q}_0^{\prime }\\ q_1^{\prime }\\ q_2^{\prime }\\ q_3^{\prime }\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)}^{\prime n}=\left[\begin{array}{ccc}{q}_0^{\prime 2}+q_1^{\prime 2}-q_2^{\prime 2}-q_3^{\prime 2}& 2(q_1^{\prime }{q}_2^{\prime }-q_0^{\prime }{q}_3^{\prime })& 2(q_1^{\prime }{q}_3^{\prime }+q_0^{\prime }{q}_2^{\prime })\\ 2(q_1^{\prime }{q}_2^{\prime }+q_0^{\prime }{q}_3^{\prime })& q_0^{\prime 2}-q_1^{\prime 2}+q_2^{\prime 2}-q_3^{\prime 2}& 2(q_2^{\prime }{q}_3^{\prime }-q_0^{\prime }{q}_1^{\prime })\\ 2(q_1^{\prime }{q}_3^{\prime }-q_0^{\prime }{q}_2^{\prime })& 2(q_2^{\prime }{q}_3^{\prime }+q_0^{\prime }{q}_1^{\prime })& q_0^{\prime 2}-q_1^{\prime 2}-q_2^{\prime 2}+q_3^{\prime 2}\end{array}\right]---\left(15\right)$更新载体姿态信息：$\left\{\begin{array}{c}{\theta }_{A2}=\arcsin \left(c_{33}^{\prime }\right)\\ {\phi }_{A2}=\arctan (c_{32}^{\prime }/c_{31}^{\prime })\\ {\psi }_{A2}=\arctan (c_{13}^{\prime }/c_{23}^{\prime })\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={-ϵ}_{z0}+\delta K_{\mathrm{gz}}{\omega }_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={ϵ}_{x0}+\delta K_{\mathrm{gx}}{\omega }_z---\left(21\right)$其中，表示在位置C处陀螺仪沿导航系z轴的输出误差；ε_xs0、δ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={-ϵ}_{x0}+\delta K_{\mathrm{gx}}{\omega }_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={ϵ}_{y0}+{\mathrm{\delta K}}_{\mathrm{gy}}{\omega }_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={-ϵ}_{y0}+{\mathrm{\delta K}}_{\mathrm{gy}}{\omega }_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轴的常值漂移和刻度因数误差。