CN103983274A - 一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法 - Google Patents

Abstract

本发明公开了一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，属于惯性技术领域。本标定方法使用低精度无方位基准双轴转位设备，整个标定旋转编排需19个位置，然后以各个位置上的速度误差和天向姿态误差拟合出一阶中间参数Δ_g以及二阶中间参数最后依据中间参数与误差参数的关系，由最小二乘法计算出各个器件误差参数，为消除由转台引起的定位误差，将前一次迭代计算得到的误差参数和原有的惯性测量单元原始输出数据代入导航方程，再进行一次观测量、中间参数和误差参数残差的解算并进行残差补偿，依此类推，直至迭代计算得到的误差参数残差小于阈值。该标定方法可以大幅降低标定对转台精度的依赖性，具有很好的工程实用性。

Description

一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法

技术领域

本发明属于航空航天捷联惯性导航技术中的惯性测量组合测试技术领域，具体涉及一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法。与传统标定方法相比，该标定方法适用于双轴转位设备，对转位设备的基准精度要求低且不需要对北。

背景技术

捷联惯性导航系统具有反应时间短、可靠性高、体积小、重量轻等优点，广泛应用于飞机、舰船、导弹等军用和民用导航领域，具有重要的国防意义和巨大的经济效益。

惯性测量组合是捷联惯性导航系统的核心部件，主要由3个加速度计和3个陀螺组成。

标定技术是惯性导航领域的核心技术之一，是一种对误差的辨识技术，即建立惯性器件和惯导系统的误差模型，通过一系列的试验求解出误差模型中的误差项，进而通过软件算法对误差进行补偿。惯性测量组合的标定结果好坏直接影响捷联惯性导航系统的精度。

惯性测量组合标定方法按层次可分为分立式标定和系统级标定两种。当前分立式标定方法的研究已经非常成熟，而系统级标定方法是由20世纪80年代发展起来的，目前正成为标定技术研究的热点。

分立标定方法是根据陀螺和加速度计的误差模型，利用三轴转台提供的精确速率、姿态和位置，采集惯性测量组合的输出，然后利用最小二乘法辨识误差模型系数。然而分立式标定过分依赖转台的精度，当转台精度不高时，标定结果不理想。

系统级标定是建立捷联惯导系统导航输出误差与惯性器件误差参数之间的关系，充分考虑惯性器件误差系数的可辨识性，合理安排实验位置，进而辨识出惯性器件的各项误差系数。该方法可以大幅减小甚至克服标定对转台精度的依赖，适合现场标定使用。

早在上世纪80-90年代，系统级标定方法就已经在国外的工程中得到了推广应用。国内的相关研究起步较晚，近年随着捷联惯导技术的成熟度不断提高，国内也出现了很多介绍系统级标定的文献和资料，但大多数停留在理论研究和仿真验证的阶段。在公开的文献和资料中，国内一般采用低精度的三轴或双轴转台，在引北的条件下在实验室内进行系统级标定。尚未有发现免对北双轴系统级标定算法的相关资料。

发明内容

本发明提供一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，与国内外其他系统级标定方法相比，该标定方法不需对北，同时可以大幅降低标定对转台精度的依赖性，具有很好的工程实用性。

本发明适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，包含如下步骤：

步骤一：将惯性测量单元（Inertial Measurement Unit，IMU）安装在双轴转位设备上，所述惯性测量单元初始位置保证X轴朝上或朝下，惯性测量单元通电预热后开始采集输出的原始数据，所述惯性测量单元先在第0个位置上静止3-5分钟后，再转动到第1个位置静止3-5分钟，随后转动到第2个位置，依此类推，直至在第18个位置上静止3-5分钟后停止采集惯性测量单元输出的原始数据；

步骤二：利用步骤一采集的惯性测量单元原始数据，在第0位置上使用双矢量定姿法进行初始对准，进而得到第0位置上导航起始时刻的天向转角然后利用对准结果和第0位置上的采集数据进行导航解算，进而得到第0位置上导航过程中的实时速度以及实时天向转角设第0位置上导航起始时刻的速度 均为0，以速度和天向转角为观测结果拟合出第0位置上的和一阶中间参数所述包含和所述和分别为第0位置上的参数在x轴、y轴和z轴上投影的标量，所述包含和所述和分别为第0位置上的一阶中间参数在x轴、y轴和z轴上投影的标量；

步骤三：利用步骤一采集的惯性测量单元原始数据，在第i位置上使用双矢量定姿法进行初始对准，所述i=0,1,2……17，然后在第i个位置到第i+1个位置的转动过程中以及第i+1个位置上的静止过程中进行连续导航，通过导航获取转动到达第i+1个位置瞬间的速度和以及瞬间的天向转角以及转动完成后在第i+1个位置静止过程中的实时速度和以及实时天向转角θⁿ⁽ⁱ⁺¹⁾，

$v_x^{n(i+1)}=v_{x0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gx}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vx}}\frac{{\mathrm{gT}}^2}{2}$

$v_y^{n(i+1)}=v_{y0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gy}}\mathrm{gT}$

$v_z^{n(i+1)}=v_{z0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gz}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vz}}\frac{{\mathrm{gT}}^2}{2}$

${\mathrm{\θ}}^{n(i+1)}={\mathrm{\θ}}_0^{n(i+1)}+{\mathrm{\ω}}_{\mathrm{vy}}T$

式中：g是重力加速度，T是实时时间，

ω_vx、ω_vy和ω_vz分别为系数ω_v在x轴、y轴和z轴上的分量，

以速度和天向转角为观测，拟合出第i+1位置上的和一阶中间参数其中，i=0,1,2……17，所述包含和 包含和所述和分别为第i+1位置上的参数在x轴、y轴和z轴上投影的标量，所述和分别为第i+1位置上的一阶中间参数在x轴、y轴和z轴上投影的标量；

步骤四：在惯性测量单元坐标系内，加速度计的误差模型为：

$\left[\begin{array}{c}\mathrm{\δ}{f}_x\\ \mathrm{\δ}{f}_y\\ \mathrm{\δ}{f}_z\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{ax}}\\ B_{\mathrm{ay}}\\ B_{\mathrm{az}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{axx}}& 0& 0\\ K_{\mathrm{ayx}}& K_{\mathrm{ayy}}& 0\\ K_{\mathrm{azx}}& K_{\mathrm{azy}}& K_{\mathrm{azz}}\end{array}\right]\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{ax}2}& 0& 0\\ 0& K_{\mathrm{ay}2}& 0\\ 0& 0& K_{\mathrm{az}2}\end{array}\right]\left[\begin{array}{c}{f_x}^2\\ {f_y}^2\\ {f_z}^2\end{array}\right]$

上述误差模型的向量形式为：

$\mathrm{\δ}{f}^b=B_a^b+K_a{f}^b+K_{a2}{\left(f^b\right)}^2$

其中，

f^b为载体坐标系下加速度计测得的比力，

f_x ^b、f_y ^b和f_z ^b分别为f^b在x轴、y轴和z轴上的投影，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

δf^b为载体坐标系下加速度计测得的比力误差；

陀螺的误差模型为：

$\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{gx}}\\ B_{\mathrm{gy}}\\ B_{\mathrm{gz}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{gxx}}& K_{\mathrm{gxy}}& K_{\mathrm{gxz}}\\ K_{\mathrm{gyx}}& K_{\mathrm{gyy}}& K_{\mathrm{gyz}}\\ K_{\mathrm{gzx}}& K_{\mathrm{gzy}}& K_{\mathrm{gzz}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]$

上述误差模型的向量形式为：

${\mathrm{\ϵ}}^b={\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b$

其中，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

ε^b为载体坐标系下陀螺测得的角速度误差；

然后将加速度计零偏B_ax、B_ay、B_az，加速度计标度因数K_axx、K_ayy、K_azz，加速度计失准角K_ayx、K_azx、K_azy，加速度计二次项系数K_ax2、K_ay2、K_az2，陀螺标度因数K_gxx、K_gyy、K_gzz，陀螺失准角K_gxy、K_gxz、K_gyx、K_gyz、K_gzx、K_gzy共计21个误差参数记为一阶误差参数K，其中，B_ax、B_ay、B_az分别为加速度计零偏B_a在x轴、y轴和z轴上投影的标量；

在每个位置上，依据一阶中间参数Δ_g与一阶误差参数K的关系，利用步骤二得到的构建方程即：

${\mathrm{\Δ}}_g^{\left(0\right)}=\frac{{\left[\mathrm{\δ}{f}^{n\left(0\right)}\right]}_V}{g}=\frac{{\left[C_{b\left(0\right)}^n(B_a^b+K_a{C}_n^{b\left(0\right)}{f}^n+K_{a2}{\left(C_n^{b\left(0\right)}{f}^n\right)}^2)\right]}_V}{g}$

其中，

δfⁿ⁽⁰⁾为导航坐标系（n系，navigation）下第0个位置的比力误差，

g是重力加速度，

是在第0个位置上从载体坐标系（b系，body）到导航坐标系的方向余弦矩阵，

是在第0个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

fⁿ为导航坐标系下的比力，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_V表示垂直分量与X相同，水平分量为0的矢量；

用步骤三中的构建方程其中，中i=0,1,2…17，中i=1,2…18，即：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_g^{\left(i\right)}=\frac{-{\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H+f^n\×\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}+\mathrm{\δ}{f}^{n\left(i\right)}}{g}\\ -{\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H=-{\left[C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\right]}_H\\ \mathrm{\δ}{f}^{n\left(i\right)}=C_{b\left(i\right)}^n(B_a^b+K_a{C}_n^{b\left(i\right)}{f}^n+K_{a2}{\left(C_n^{b\left(i\right)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

δf^n(i-1)和δfⁿ⁽ⁱ⁾分别为导航坐标系下第i-1和第i个位置下的比力误差，

fⁿ为导航坐标系下的比力，

ΔΦ^n(i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差，

g是重力加速度，

和分别是在第i-1个和第i个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

和分别是在第i-1个和第i个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

将以上方程联立得到如下方程：

${\mathrm{\Δ}}_g={\left[\begin{array}{c}{\mathrm{\Δ}}_g^{\left(0\right)}\\ {\mathrm{\Δ}}_g^{\left(1\right)}\\ {\mathrm{\Δ}}_g^{\left(2\right)}\\ \·\\ \·\\ \·\\ {\mathrm{\Δ}}_g^{\left(18\right)}\end{array}\right]}_{57\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ \·\\ \·\\ \·\\ A^{\left(18\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程：

Δ_g=AK

步骤五：利用步骤四中的联立方程计算A的最小二乘逆矩阵进而通过计算出一阶误差参数K；

步骤六：利用步骤五计算的K求解每个静止位置i对应的其中i=0,1,2…18，计算方法如下：

当i=0时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(0\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(0\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f_z}^{n\left(0\right)}}{g}& \frac{\tan L\·\mathrm{\δ}{f_z}^{n\left(0\right)}}{g}& -\frac{\mathrm{\δ}{f_x}^{n\left(0\right)}}{g}\end{array}\right]}^T\end{array}\right.$

其中，

为在第0个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第0个位置的姿态误差，包括第0个位置上初始对准引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

[X]_V代表垂直分量与X相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

δf_x ⁿ⁽⁰⁾和δf_z ⁿ⁽⁰⁾分别为导航坐标系下第0个位置计算的比力误差在X轴和Z轴的分量，

L为纬度；

当i=1,2……18时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(i\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f_z}^{n(i-1)}}{g}& \frac{\tan L\·\mathrm{\δ}{f_z}^{n(i-1)}}{g}& -\frac{\mathrm{\δ}{f_x}^{n(i-1)}}{g}\end{array}\right]}^T+\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}\\ \mathrm{\δ}{f}^{n(i-1)}=C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

为在第i个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与x轴相同，垂直分量为0的矢量，

[X]_V代表垂直分量与x轴相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

δf_x ^n(i-1)和δf_z ^n(i-1)分别为导航坐标系下第i-1个位置计算的比力误差在X轴和Z轴的分量，

L为纬度，

ΔΦ^n(i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差，

δf^n(i-1)为导航坐标系下第i-1个位置下的比力误差，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

是在第i-1个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

然后通过步骤二中的及步骤三中的计算出每个静止位置上的二阶中间参数其中，中i=0,1,2……17，中i=0,1,2,……,18；

步骤七：将二阶误差参数即陀螺零偏B_gx、B_gy、B_gz记为列向量ω，依据二阶中间参数和二阶误差参数ω之间的关系，利用步骤六中构建方程其中i=0,1,2……18，即：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{*\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)})]}_H}{g}+{[{\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)}]}_V\\ {\mathrm{\Φ}}_0^{*n\left(i\right)}={\left[\begin{array}{ccc}0& -\frac{{\mathrm{\ω}}_{0z}^{n(i-1)}}{{\mathrm{\ω}}_{\mathrm{ie}}\cos L}& 0\end{array}\right]}^T\\ {\mathrm{\ω}}_0^{n(i-1)}=C_{b^{(i-1)}}^n{\mathrm{\ω}}_0^b\end{array}\right.$

其中，

表示第i个位置的二阶中间参数，

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与二阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

和分别为导航坐标系下第i-1个和第i个位置的等效陀螺零偏，

导航坐标系下第i-1个位置的等效陀螺零偏在Z轴上的投影，

ω_ie为地球自转角速率，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

为载体坐标系下的陀螺零偏；

将以上方程联立得到如下方程：

${\mathrm{\ω}}_v^{*}=\left[\begin{array}{c}{{\mathrm{\ω}}_v^{*}}^{\left(0\right)}\\ {{\mathrm{\ω}}_v^{*}}^{\left(1\right)}\\ {{\mathrm{\ω}}_v^{*}}^{\left(2\right)}\\ \·\\ \·\\ \·\\ {{\mathrm{\ω}}_v^{*}}^{\left(18\right)}\end{array}\right]=\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ \·\\ \·\\ \·\\ B^{\left(18\right)}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}$

步骤八：利用步骤七中的联立方程计算B的最小二乘逆矩阵进而通过计算出二阶误差参数ω；

步骤九：当一阶误差参数K和二阶误差参数ω的残差大于阈值时，用一阶误差参数K和二阶误差参数ω残差补偿前次标定的误差参数，然后将得到的一阶误差参数K和二阶误差参数ω以及步骤一中采集的惯性测量单元原始数据代入到导航方程中，再进行一次一阶中间参数Δ_g、二阶中间参数一阶误差参数K和二阶误差参数ω残差的解算，然后对一阶误差参数K和二阶误差参数ω进行残差补偿，依此类推，经过多次迭代直至某一次迭代计算得到的一阶误差参数K和二阶误差参数ω残差小于阈值。

在上述技术方案中，在步骤一中，导航坐标系选取北天东（N-U-E，North-Up-East）坐标系。

在上述技术方案中，在步骤一中，标定旋转顺序如表1所示：

表1标定旋转顺序

位置序号	旋转过程
		0	+90Z
1	+90Z

2	+90Z
		3	-90Z
4	-90Z
		5	-90Z
6	+90X
		7	+90Y
8	+90Y
		9	+90Y
10	-90Y
		11	-90Y
12	-90Y
		13	+90X
14	+90X
		15	-90X
16	-90X
		17	-90X
18

在上述技术方案中，惯性测量单元坐标系是：X轴与X加速度计输入轴方向相同，Y轴位于X加速度计和Y加速度计输入轴构成的平面内，接近Y加速度计输入轴方向，Z轴方向由右手定则确定。

在上述技术方案中，在步骤一中，所述的惯性测量单元通电预热时间为30分钟，原始数据的采样周期为0.01s。

在上述技术方案中，在步骤一中，停止采集惯性测量单元后关闭惯性测量单元。

本方法原理描述如下：

标定方法利用采集的原始数据，在第i（i=0,1,2……17）位置上进行初始对准，然后在第i个位置到第i+1个位置的转动过程中以及第i+1个位置上的静止过程中进行连续导航。在每个静止位置上，天向的速度误差和姿态误差呈线性增长，水平方向的速度误差呈二次曲线增长。又在静止位置上，真实的速度和绕天向的转角为0，因此导航计算得到的速度增量即为速度误差增量，绕天向的转角增量即为天向的姿态误差增量。由此，可以将其作为观测量，按下式对第i个静止位置上导航计算得到的速度和天向转角进行拟合，即：

$v_x^{n\left(i\right)}=v_{x0}^{n\left(i\right)}+{\mathrm{\Δ}}_{\mathrm{gx}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vx}}\frac{{\mathrm{gT}}^2}{2}$

$v_y^{n\left(i\right)}=v_{y0}^{n\left(i\right)}+{\mathrm{\Δ}}_{\mathrm{gy}}\mathrm{gT}$

$v_z^{n\left(i\right)}=v_{z0}^{n\left(i\right)}+{\mathrm{\Δ}}_{\mathrm{gz}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vz}}\frac{{\mathrm{gT}}^2}{2}$

${\mathrm{\θ}}^{n\left(i\right)}={\mathrm{\θ}}_0^{n\left(i\right)}+{\mathrm{\ω}}_{\mathrm{vy}}T$

式中：

为到达静止位置瞬间的速度（带圆括号的上标i表示第i个位置，下同）；

为到达静止位置瞬间的天向转角。

上式中的各个系数与误差参数相关，将Δ_g（包含Δ_gx、Δ_gy、Δ_gz）等系数称为一阶中间参数，它们与加速度计误差参数、陀螺的标度因数误差和失准角相关，加速度计误差参数、陀螺的标度因数误差和失准角又称为一阶误差参数。为降低系数ω_v（包含ω_vx、ω_vy、ω_vz）与误差参数的耦合度，将其分解为和两个分量，称为二阶中间参数，它与陀螺的零偏误差相关，后者称为二阶误差参数。

具体来说，以各个位置上的速度误差和天向姿态误差拟合出一阶中间参数构成的列向量Δ_g以及二阶中间参数构成的列向量然后依据中间参数与误差参数的关系，由最小二乘法计算出各个器件误差参数。设各一阶误差参数构成列向量K，各二阶误差参数构成列向量ω。可以用矩阵形式将其关系表示为：

Δ_g=AK

${\mathrm{\ω}}_v^{*}=\mathrm{B\ω}$

为了有效消除由转台引起的定位误差，可将计算得到的误差参数K、ω以及采集的惯性测量单元原始数据代入到导航方程中，再进行一次观测量、中间参数和误差参数残差的解算，然后对误差参数进行残差补偿。依此类推，经过多次迭代直至某一次迭代计算得到的误差参数残差小于一定的阈值为止。

本发明适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法的有益效果在于：

（1）该标定方法使用双矢量定姿法进行初始对准（采用g及g×ω_ie作为参考矢量），所以不用对北，适用于在低精度双轴转位设备上进行标定；

（2）该标定方法采用迭代算法，同时标定旋转通过编排19个位置，可以大幅降低标定对转台精度的依赖性。

附图说明

图1为本发明适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法的流程图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。

本发明提供一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，具体标定步骤如下：

步骤一：将惯性测量单元安装在双轴转位设备上，所述惯性测量单元初始位置保证X轴朝上或朝下，优选的，导航坐标系选取北天东坐标系；惯性测量单元通电预热30分钟开始采集输出的原始数据，采样周期为0.01s，标定旋转顺序如表一所示：所述惯性测量单元先在第0个位置上静止3-5分钟后，再转动到第1个位置静止3-5分钟，随后转动到第2个位置，依此类推，直至在第18个位置上静止3-5分钟后停止采集惯性测量单元输出的原始数据，然后关闭惯性测量单元。

表1标定旋转顺序

位置序号	旋转过程
		0	+90Z
1	+90Z
		2	+90Z
3	-90Z
		4	-90Z
5	-90Z
		6	+90X
7	+90Y
		8	+90Y
9	+90Y
		10	-90Y
11	-90Y
		12	-90Y
13	+90X
		14	+90X
15	-90X
		16	-90X
17	-90X
		18

步骤二：利用步骤一采集的惯性测量单元原始数据，在第0位置上使用双矢量定姿法进行初始对准，进而得到第0位置上导航起始时刻的天向转角然后利用对准结果和第0位置上的采集数据进行导航解算，进而得到第0位置上导航过程中的实时速度以及实时天向转角θⁿ⁽⁰⁾，设第0位置上导航起始时刻的速度 均为0，以速度和天向转角为观测结果拟合出第0位置上的和一阶中间参数所述包含所述和分别为矢量在x轴、y轴和z轴上投影的标量，所述包含所述和分别为矢量在x轴、y轴和z轴上投影的标量。

步骤三：利用步骤一采集的惯性测量单元数据，在第i位置上使用双矢量定姿法进行初始对准，所述i=0,1,2……17，然后在第i个位置到第i+1个位置的转动过程中以及第i+1个位置上的静止过程中进行连续导航，通过导航获取转动到达第i+1个位置瞬间的速度 和以及瞬间的天向转角以及转动完成后在第i+1个位置静止过程中的实时速度和以及实时天向转角θⁿ⁽ⁱ⁺¹⁾，

$v_x^{n(i+1)}=v_{x0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gx}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vx}}\frac{{\mathrm{gT}}^2}{2}$

$v_y^{n(i+1)}=v_{y0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gy}}\mathrm{gT}$

$v_z^{n(i+1)}=v_{z0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gz}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vz}}\frac{{\mathrm{gT}}^2}{2}$

${\mathrm{\θ}}^{n(i+1)}={\mathrm{\θ}}_0^{n(i+1)}+{\mathrm{\ω}}_{\mathrm{vy}}T$

式中：g是重力加速度，T是实时时间，

ω_vx、ω_vy和ω_vz分别为系数ω_v在x轴、y轴和z轴上的分量，

以速度和天向转角为观测，拟合出第i+1位置上的和一阶中间参数其中，i=0,1,2……17，所述包含和 包含和所述和分别为第i+1位置上的参数在x轴、y轴和z轴上投影的标量，所述和分别为第i+1位置上的一阶中间参数在x轴、y轴和z轴上投影的标量。

步骤四：惯性测量单元坐标系是：X轴与X加速度计输入轴方向相同，Y轴位于X加速度计和Y加速度计输入轴构成的平面内，接近Y加速度计输入轴方向，Z轴方向由右手定则确定；

在此坐标系内，加速度计的误差模型为：

$\left[\begin{array}{c}\mathrm{\δ}{f}_x\\ \mathrm{\δ}{f}_y\\ \mathrm{\δ}{f}_z\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{ax}}\\ B_{\mathrm{ay}}\\ B_{\mathrm{az}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{axx}}& 0& 0\\ K_{\mathrm{ayx}}& K_{\mathrm{ayy}}& 0\\ K_{\mathrm{azx}}& K_{\mathrm{azy}}& K_{\mathrm{azz}}\end{array}\right]\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{ax}2}& 0& 0\\ 0& K_{\mathrm{ay}2}& 0\\ 0& 0& K_{\mathrm{az}2}\end{array}\right]\left[\begin{array}{c}{f_x}^2\\ {f_y}^2\\ {f_z}^2\end{array}\right]$

上述误差模型的向量形式为：

$\mathrm{\δ}{f}^b=B_a^b+K_a{f}^b+K_{a2}{\left(f^b\right)}^2$

其中，

f^b为载体坐标系下加速度计测得的比力，

f_x ^b、f_y ^b和f_z ^b分别为f^b在x轴、y轴和z轴上的投影，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

δf^b为载体坐标系下加速度计测得的比力误差；

陀螺的误差模型为：

$\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{gx}}\\ B_{\mathrm{gy}}\\ B_{\mathrm{gz}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{gxx}}& K_{\mathrm{gxy}}& K_{\mathrm{gxz}}\\ K_{\mathrm{gyx}}& K_{\mathrm{gyy}}& K_{\mathrm{gyz}}\\ K_{\mathrm{gzx}}& K_{\mathrm{gzy}}& K_{\mathrm{gzz}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]$

上述误差模型的向量形式为：

${\mathrm{\ϵ}}^b={\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b$

其中，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

ε^b为载体坐标系下陀螺测得的角速度误差；

然后将各一阶误差参数：加速度计零偏B_ax、B_ay、B_az，加速度计标度因数K_axx、K_ayy、K_azz，加速度计失准角K_ayx、K_azx、K_azy，加速度计二次项系数K_ax2、K_ay2、K_az2，陀螺标度因数K_gxx、K_gyy、K_gzz，陀螺失准角K_gxy、K_gxz、K_gyx、K_gyz、K_gzx、K_gzy共计21个误差参数记为列向量K，其中，B_ax、B_ay、B_az分别为加速度计零偏B_a在x轴、y轴和z轴上投影的标量；

在每个位置上，依据一阶中间参数Δ_g和一阶误差参数K之间的关系，利用步骤二得到的构建方程即：

${\mathrm{\Δ}}_g^{\left(0\right)}=\frac{{\left[\mathrm{\δ}{f}^{n\left(0\right)}\right]}_V}{g}=\frac{{\left[C_{b\left(0\right)}^n(B_a^b+K_a{C}_n^{b\left(0\right)}{f}^n+K_{a2}{\left(C_n^{b\left(0\right)}{f}^n\right)}^2)\right]}_V}{g}$

其中，

δfⁿ⁽⁰⁾为导航坐标系（n系，navigation）下第0个位置的比力误差，

g是重力加速度，

是在第0个位置上从载体坐标系（b系，body）到导航坐标系的方向余弦矩阵，

是在第0个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

fⁿ为导航坐标系下的比力，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_V表示垂直分量与X相同，水平分量为0的矢量；

用步骤三中的构建方程其中，中i=0,1,2…17，中i=1,2…18，即：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_g^{\left(i\right)}=\frac{-{\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H+f^n\×\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}+\mathrm{\δ}{f}^{n\left(i\right)}}{g}\\ -{\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H=-{\left[C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\right]}_H\\ \mathrm{\δ}{f}^{n\left(i\right)}=C_{b\left(i\right)}^n(B_a^b+K_a{C}_n^{b\left(i\right)}{f}^n+K_{a2}{\left(C_n^{b\left(i\right)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

δf^n(i-1)和δfⁿ⁽ⁱ⁾分别为导航坐标系下第i-1和第i个位置下的比力误差；

fⁿ为导航坐标系下的比力

ΔΦ^n(i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差；

g是重力加速度，

和分别是在第i-1个和第i个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

和分别是在第i-1个和第i个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

将以上方程联立得到如下方程：

${\mathrm{\Δ}}_g={\left[\begin{array}{c}{\mathrm{\Δ}}_g^{\left(0\right)}\\ {\mathrm{\Δ}}_g^{\left(1\right)}\\ {\mathrm{\Δ}}_g^{\left(2\right)}\\ \·\\ \·\\ \·\\ {\mathrm{\Δ}}_g^{\left(18\right)}\end{array}\right]}_{57\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ \·\\ \·\\ \·\\ A^{\left(18\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程：

Δ_g=AK

步骤五：利用步骤四中的联立方程计算A的最小二乘逆矩阵进而通过计算出一阶误差参数K。

步骤六：利用步骤五中的K计算出每个静止位置i对应的其中i=0,1,2…18，计算方法如下：

当i=0时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(0\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(0\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f_z}^{n\left(0\right)}}{g}& \frac{\tan L\·\mathrm{\δ}{f_z}^{n\left(0\right)}}{g}& -\frac{\mathrm{\δ}{f_x}^{n\left(0\right)}}{g}\end{array}\right]}^T\end{array}\right.$

其中，

为在第0个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第0个位置的姿态误差，包括第0个位置上初始对准引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

[X]_V代表垂直分量与X相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

δf_x ⁿ⁽⁰⁾和δf_z ⁿ⁽⁰⁾分别为导航坐标系下第0个位置计算的比力误差在X轴和Z轴的分量，

L为纬度；

当i=1,2……18时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(i\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f_z}^{n(i-1)}}{g}& \frac{\tan L\·\mathrm{\δ}{f_z}^{n(i-1)}}{g}& -\frac{\mathrm{\δ}{f_x}^{n(i-1)}}{g}\end{array}\right]}^T+\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}\\ \mathrm{\δ}{f}^{n(i-1)}=C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

为在第i个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与x轴相同，垂直分量为0的矢量，

[X]_V代表垂直分量与x轴相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

δf_x ^n(i-1)和δf_z ^n(i-1)分别为导航坐标系下第i-1个位置计算的比力误差在X轴和Z轴的分量，

L为纬度，

ΔΦ^n(i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差，

δf^n(i-1)为导航坐标系下第i-1个位置下的比力误差，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

是在第i-1个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

然后通过步骤二中的及步骤三中的计算出每个静止位置上的二阶中间参数其中，中i=0,1,2……17，中i=0,1,2,……,18。

步骤七：将二阶误差参数即陀螺零偏B_gx、B_gy、B_gz记为列向量ω，依据二阶中间参数和二阶误差参数ω之间的关系，利用步骤六中构建方程其中i=0,1,2……18，即：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{*\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)})]}_H}{g}+{[{\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)}]}_V\\ {\mathrm{\Φ}}_0^{*n\left(i\right)}={\left[\begin{array}{ccc}0& -\frac{{\mathrm{\ω}}_{0z}^{n(i-1)}}{{\mathrm{\ω}}_{\mathrm{ie}}\cos L}& 0\end{array}\right]}^T\\ {\mathrm{\ω}}_0^{n(i-1)}=C_{b^{(i-1)}}^n{\mathrm{\ω}}_0^b\end{array}\right.$

其中，

表示第i个位置的二阶中间参数，

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与二阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

和分别为导航坐标系下第i-1个和第i个位置的等效陀螺零偏，

导航坐标系下第i-1个位置的等效陀螺零偏在Z轴上的投影，

ω_ie为地球自转角速率，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

为载体坐标系下的陀螺零偏；

将以上方程联立得到如下方程：

${\mathrm{\ω}}_v^{*}=\left[\begin{array}{c}{{\mathrm{\ω}}_v^{*}}^{\left(0\right)}\\ {{\mathrm{\ω}}_v^{*}}^{\left(1\right)}\\ {{\mathrm{\ω}}_v^{*}}^{\left(2\right)}\\ \·\\ \·\\ \·\\ {{\mathrm{\ω}}_v^{*}}^{\left(18\right)}\end{array}\right]=\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ \·\\ \·\\ \·\\ B^{\left(18\right)}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}$

步骤八：利用步骤七中的联立方程计算B的最小二乘逆矩阵进而通过计算出二阶误差参数ω。

步骤九：当一阶误差参数K和二阶误差参数ω的残差大于阈值时，用一阶误差参数K和二阶误差参数ω残差补偿前次标定的误差参数。然后将得到的一阶误差参数K和二阶误差参数ω的以及步骤一中采集的惯性测量单元原始数据代入到导航方程中，再进行一次一阶中间参数Δ_g、二阶中间参数一阶误差参数K和二阶误差参数ω残差的解算，然后对一阶误差参数K和二阶误差参数ω进行残差补偿，依此类推，经过多次迭代直至某一次迭代计算得到的一阶误差参数K和二阶误差参数ω残差小于阈值。

Claims (6)

1.一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：包含如下步骤：

步骤一：将惯性测量单元安装在双轴转位设备上，所述惯性测量单元初始位置保证X轴朝上或朝下，惯性测量单元通电预热后开始采集输出的原始数据，所述惯性测量单元先在第0个位置上静止3-5分钟后，再转动到第1个位置静止3-5分钟，随后转动到第2个位置，依此类推，直至在第18个位置上静止3-5分钟后停止采集惯性测量单元输出的原始数据；

步骤二：利用步骤一采集的惯性测量单元原始数据，在第0位置上使用双矢量定姿法进行初始对准，进而得到第0位置上导航起始时刻的天向转角然后利用对准结果和第0位置上的采集数据进行导航解算，进而得到第0位置上导航过程中的实时速度以及实时天向转角θⁿ⁽⁰⁾，设第0位置上导航起始时刻的速度 均为0，以速度和天向转角为观测结果拟合出第0位置上的和一阶中间参数所述包含和所述和分别为第0位置上的参数在x轴、y轴和z轴上投影的标量，所述包含和所述和分别为第0位置上的一阶中间参数在x轴、y轴和z轴上投影的标量；

步骤三：利用步骤一采集的惯性测量单元原始数据，在第i位置上使用双矢量定姿法进行初始对准，所述i＝0，1，2……17，然后在第i个位置到第i+1个位置的转动过程中以及第i+1个位置上的静止过程中进行连续导航，通过导航获取转动到达第i+1个位置瞬间的速度和以及瞬间的天向转角以及转动完成后在第i+1个位置静止过程中的实时速度和以及实时天向转角θⁿ⁽ⁱ⁺¹⁾，

$v_x^{n(i+1)}=v_{x0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gx}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vx}}\frac{g{T}^2}{2}$

$v_y^{n(i+1)}=v_{y0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gx}}\mathrm{gT}$

$v_z^{n(i+1)}=v_{z0}^{n(i+1)}+{\mathrm{\Δ}}_{\mathrm{gz}}\mathrm{gT}+{\mathrm{\ω}}_{\mathrm{vz}}\frac{g{T}^2}{2}$

${\mathrm{\θ}}^{n(i+1)}={\mathrm{\θ}}_0^{n(i+1)}+{\mathrm{\ω}}_{\mathrm{vy}}T$

式中：g是重力加速度，T是实时时间，

ω_vx、ω_vy和ω_vz分别为系数ω_v在x轴、y轴和z轴上的分量，

以速度和天向转角为观测，拟合出第i+1位置上的和一阶中间参数其中，i＝0，1，2……17，所述包含和 包含和所述和分别为第i+1位置上的参数在x轴、y轴和z轴上投影的标量，所述和分别为第i+1位置上的一阶中间参数在x轴、y轴和z轴上投影的标量；

步骤四：在惯性测量单元坐标系内，加速度计的误差模型为：：

$\left[\begin{array}{c}\mathrm{\δ}{f}_x\\ \mathrm{\δ}{f}_y\\ \mathrm{\δ}{f}_z\end{array}\right]=\begin{array}{c}\left[\begin{array}{c}{B}_{\mathrm{ax}}\\ B_{\mathrm{ay}}\\ B_{\mathrm{az}}\end{array}\right]+\begin{array}{c}\left[\begin{array}{ccc}{K}_{\mathrm{axx}}& 0& 0\\ K_{\mathrm{ayx}}& K_{\mathrm{ayy}}& 0\\ K_{\mathrm{azx}}& K_{\mathrm{azy}}& K_{\mathrm{azz}}\end{array}\right]\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]+\end{array}\end{array}\left[\begin{array}{ccc}{K}_{\mathrm{ax}2}& 0& 0\\ 0& K_{\mathrm{ay}2}& 0\\ 0& 0& K_{\mathrm{az}2}\end{array}\right]\left[\begin{array}{c}{f_x}^2\\ {f_y}^2\\ {f_z}^2\end{array}\right]$

上述误差模型的向量形式为：

$\mathrm{\δ}{f}^b=B_a^b+K_a{f}^b+K_{a2}{\left(f^b\right)}^2$

其中，

f^b为载体坐标系下加速度计测得的比力，

f_x ^b、f_y ^b和f_z ^b分别为f^b在x轴、y轴和z轴上的投影，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

δf^b为载体坐标系下加速度计测得的比力误差；

陀螺的误差模型为：

$\left[\begin{array}{c}{\mathrm{\ϵ}}_x\\ {\mathrm{\ϵ}}_y\\ {\mathrm{\ϵ}}_z\end{array}\right]=\left[\begin{array}{c}{B}_{\mathrm{gx}}\\ B_{\mathrm{gy}}\\ B_{\mathrm{gz}}\end{array}\right]+\left[\begin{array}{ccc}{K}_{\mathrm{gxx}}& K_{\mathrm{gxy}}& K_{\mathrm{gxz}}\\ K_{\mathrm{gyx}}& K_{\mathrm{gyy}}& K_{\mathrm{gyz}}\\ K_{\mathrm{gzx}}& K_{\mathrm{gzy}}& K_{\mathrm{gzz}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ω}}_x\\ {\mathrm{\ω}}_y\\ {\mathrm{\ω}}_z\end{array}\right]$

上述误差模型的向量形式为：

${\mathrm{\ϵ}}^b={\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b$

其中，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

ε^b为载体坐标系下陀螺测得的角速度误差；

然后将加速度计零偏B_ax、B_ay、B_az，加速度计标度因数K_axx、K_ayy、K_azz，加速度计失准角K_ayx、K_azx、K_azy，加速度计二次项系数K_ax2、K_ay2、K_az2，陀螺标度因数K_gxx、K_gyy、K_gzz，陀螺失准角K_gxy、K_gxz、K_gyx、K_gyz、K_gzx、K_gzy共计21个误差参数记为一阶误差参数K，其中，B_ax、B_ay、B_az分别为加速度计零偏B_a在x轴、y轴和z轴上投影的标量；

在每个位置上，依据一阶中间参数Δ_g与一阶误差参数K的关系，利用步骤二得到的构建方程即：

${\mathrm{\Δ}}_g^{\left(0\right)}=\frac{[\mathrm{\δ}{f}^{n\left(0\right)}{]}_V}{g}=\frac{{\left[C_{b\left(0\right)}^n(B_a^b+K_a{C}_n^{b\left(0\right)}{f}^n+K_{a2}{\left(C_n^{b\left(0\right)}{f}^n\right)}^2)\right]}_V}{g}$

其中，

δfⁿ⁽⁰⁾为导航坐标系(n系，naｖigation)下第0个位置的比力误差，

g是重力加速度，

是在第0个位置上从载体坐标系(b系，body)到导航坐标系的方向余弦矩阵，

是在第0个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

fⁿ为导航坐标系下的比力，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_V表示垂直分量与X相同，水平分量为0的矢量；

用步骤三中的构建方程其中，中i＝0，1，2…17，中i＝1，2…18，即：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_g^{\left(i\right)}=\frac{{-\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H+f^n\×\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}+\mathrm{\δ}{f}^{n\left(i\right)}}{g}\\ -{\left[\mathrm{\δ}{f}^{n(i-1)}\right]}_H=-{\left[C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\right]}_H\\ \mathrm{\δ}{f}^{n\left(i\right)}=C_{b\left(i\right)}^n(B_a^b+K_a{C}_n^{b\left(i\right)}{f}^n+K_{a2}{\left(C_n^{b\left(i\right)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

δf^n(i-1)和δfⁿ⁽ⁱ⁾分别为导航坐标系下第i-1和第i个位置下的比力误差，

fⁿ为导航坐标系下的比力，

ΔΦ^n(i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差，

g是重力加速度，

和分别是在第i-1个和第i个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

和分别是在第i-1个和第i个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

将以上方程联立得到如下方程：

${\mathrm{\Δ}}_g={\left[\begin{array}{c}{\mathrm{\Δ}}_g^{\left(0\right)}\\ {\mathrm{\Δ}}_g^{\left(1\right)}\\ {\mathrm{\Δ}}_g^{\left(2\right)}\\ \·\\ \·\\ \·\\ {\mathrm{\Δ}}_g^{\left(18\right)}\end{array}\right]}_{57\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ \·\\ \·\\ \·\\ A^{\left(18\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程：

Δ_g＝AK

步骤五：利用步骤四中的联立方程计算A的最小二乘逆矩阵进而通过计算出一阶误差参数K；

步骤六：利用步骤五计算的K求解每个静止位置i对应的其中i＝0，1，2…18，计算方法如下：

当i＝0时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(0\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(0\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(0\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f}_z^{n\left(0\right)}}{g}& \frac{\tan L\·\mathrm{\δ}{f}_z^{n\left(0\right)}}{g}& -\frac{\mathrm{\δ}{f}_x^{n\left(0\right)}}{g}\end{array}\right]}^T\end{array}\right.$

其中，

为在第0个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第0个位置的姿态误差，包括第0个位置上初始对准引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与X相同，垂直分量为0的矢量，

[X]_V代表垂直分量与X相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

和分别为导航坐标系下第0个位置计算的比力误差在X轴和Z轴的分量，

L为纬度；

当i＝1，2……18时，的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{0\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n)]}_H}{g}+{[{\mathrm{\Φ}}_0^{0n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n]}_V\\ {\mathrm{\Φ}}_0^{0n\left(i\right)}={\left[\begin{array}{ccc}\frac{\mathrm{\δ}{f}_z^{n(i-1)}}{g}& \frac{\tan L\·\mathrm{\δ}{f}_z^{n(i-1)}}{g}& -\frac{\mathrm{\δ}{f}_x^{n(i-1)}}{g}\end{array}\right]}^T+\mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}\\ \mathrm{\δ}{f}^{n(i-1)}=C_{b(i-1)}^n(B_a^b+K_a{C}_n^{b(i-1)}{f}^n+K_{a2}{\left(C_n^{b(i-1)}{f}^n\right)}^2)\\ \mathrm{\Δ}{\mathrm{\Φ}}^{n(i-1~i)}=-\underset{t=0}{\overset{T}{\∫}}{\mathrm{\ϵ}}^n\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n{\mathrm{\ϵ}}^b\mathrm{dt}=-\underset{t=0}{\overset{T}{\∫}}{C}_{b(i-1~i)}^n({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}\end{array}\right.$

其中，

为在第i个位置上的

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与一阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

g是重力加速度，

[X]_H代表水平分量与x轴相同，垂直分量为0的矢量，

[X]_V代表垂直分量与x轴相同，水平分量为0的矢量，

[X]^T为矩阵X或向量X的转置，

和分别为导航坐标系下第i-1个位置计算的比力误差在X轴和Z轴的分量，

L为纬度，

ΔΦⁿ(^i-1～i)为在导航坐标系下从第i-1到第i个位置旋转过程中引入的姿态误差，

δf^n(i-1)为导航坐标系下第i-1个位置下的比力误差，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

是在第i-1个位置上从导航坐标系到载体坐标系的方向余弦矩阵，

为载体坐标系下的加速度计零偏，

K_a包括加速度计标度因数误差和加速度计失准角，

K_a2为加速度计二次项系数，

εⁿ为陀螺测得的角速度误差在导航坐标系下的投影，

ε^b为载体坐标系下陀螺测得的角速度误差，

为从第i-1个位置到第i个位置的旋转过程中，载体坐标系到导航坐标系的方向余弦矩阵，

ω^b为载体坐标系下陀螺测得的角速度，

为载体坐标系下陀螺零偏，

K_g包括陀螺标度因数误差和陀螺失准角，

T是实时时间；

然后通过步骤二中的及步骤三中的计算出每个静止位置上的二阶中间参数其中，中i＝0，1，2……17，中i＝0，1，2，……，18；

步骤七：将二阶误差参数即陀螺零偏B_gx、B_gy、B_gz记为列向量ω，依据二阶中间参数和二阶误差参数ω之间的关系，利用步骤六中构建方程其中i＝0，1，2……18，即：

$\left\{\begin{array}{c}{\mathrm{\ω}}_v^{*\left(i\right)}=\frac{{[f^n\×({\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)})]}_H}{g}+{[{\mathrm{\Φ}}_0^{*n\left(i\right)}\×{\mathrm{\ω}}_{\mathrm{ie}}^n-{\mathrm{\ω}}_0^{n\left(i\right)}]}_V\\ {\mathrm{\Φ}}_0^{*n\left(i\right)}={\left[\begin{array}{ccc}0& -\frac{{\mathrm{\ω}}_{0z}^{n(i-1)}}{{\mathrm{\ω}}_{\mathrm{ie}}\cos L}& 0\end{array}\right]}^T\\ {\mathrm{\ω}}_0^{n(i-1)}=C_{b^{(i-1)}}^{n^8}{\mathrm{\ω}}_0^b\end{array}\right.$

其中，

表示第i个位置的二阶中间参数，

fⁿ为导航坐标系下的比力，

为导航坐标系下第i个位置的姿态误差，包括第i-1个位置上初始对准引入的姿态误差及第i-1个位置到第i个位置旋转过程中引入的姿态误差，

为中与二阶误差参数相关的项，

为地球自转角速度在导航坐标系下的投影，

和分别为导航坐标系下第i-1个和第i个位置的等效陀螺零偏，

导航坐标系下第i-1个位置的等效陀螺零偏在Z轴上的投影，

ω_ie为地球自转角速率，

是在第i-1个位置上从载体坐标系到导航坐标系的方向余弦矩阵，

为载体坐标系下的陀螺零偏；

将以上方程联立得到如下方程：

${\mathrm{\ω}}_v^{*}=\left[\begin{array}{c}{\mathrm{\ω}}_v^{*\left(0\right)}\\ {\mathrm{\ω}}_v^{*\left(1\right)}\\ {\mathrm{\ω}}_v^{*\left(2\right)}\\ .\\ .\\ .\\ {\mathrm{\ω}}_v^{*\left(18\right)}\end{array}\right]=\begin{array}{c}\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ .\\ .\\ .\\ B^{\left(18\right)}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}\end{array}$

步骤八：利用步骤七中的联立方程计算B的最小二乘逆矩阵进而通过计算出二阶误差参数ω；

步骤九：当一阶误差参数K和二阶误差参数ω的残差大于阈值时，用一阶误差参数K和二阶误差参数ω残差补偿前次标定的误差参数，然后将得到的一阶误差参数K和二阶误差参数ω以及步骤一中采集的惯性测量单元原始数据代入到导航方程中，再进行一次一阶中间参数Δ_g、二阶中间参数一阶误差参数K和二阶误差参数ω残差的解算，然后对一阶误差参数K和二阶误差参数ω进行残差补偿，依此类推，经过多次迭代直至某一次迭代计算得到的一阶误差参数K和二阶误差参数ω残差小于阈值。

2.根据权利要求1所述的一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：在步骤一中，导航坐标系选取北天东坐标系。

3.根据权利要求1或2所述的一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：在步骤一中，标定旋转顺序如下表所示：

标定旋转顺序

4.根据权利要求1或2所述的一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：惯性测量单元坐标系是：X轴与X加速度计输入轴方向相同，Y轴位于X加速度计和Y加速度计输入轴构成的平面内，接近Y加速度计输入轴方向，Z轴方向由右手定则确定。

5.根据权利要求1或2所述的一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：在步骤一中，所述的惯性测量单元通电预热时间为30分钟，原始数据的采样周期为0.01s。

6.根据权利要求1或2所述的一种适用于低精度无方位基准双轴转位设备的惯性测量单元标定方法，其特征在于：在步骤一中，停止采集惯性测量单元后关闭惯性测量单元。