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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

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

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

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

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

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

发明内容

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

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

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

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

步骤三：利用步骤一采集的惯性测量单元原始数据，在第i位置上使用双矢量定姿法进行初始对准，i＝0，1，2，3，然后在第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，3，所述包含和 包含和所述和分别为第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，3，中i＝1，2，3，4，即： 

$\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=-[C_{\text{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){]}_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(3\right)}\\ {\mathrm{\Δ}}_g^{\left(4\right)}\end{array}\right]}_{15\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ A^{\left(3\right)}\\ A^{\left(4\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程： 

Δ_g＝AK 

步骤五：根据一阶中间参数与一阶误差参数的关系式，设A＝[a₁ a₂ … a_n]、K＝[k₁ k₂ … k_n]^T，其中，a_i为矩阵A的列向量，i＝1，2…n，k_i为向量K的各元素，i＝1，2…n，[k₁ k₂ … k_n]^T为行向量[k₁ k₂ … k_n]的转置， 

K中包含B_ay、K_ayy、K_azy、K_ax2、K_ay2、K_az2、K_gxx、K_gzz、K_gxy、K_gxz、K_gyx、K_gyz、K_gzx13个误差参数直接由最近一次的完整标定结果给出，这些由外部给出的误差参数在向量K中的序号为e_l，l＝1，2，...，n_e，其余参与标定计算的误差参数在向量K中的序号为c_j，j＝1，2，...，n_c，且n_c+n_e＝n；然后将矩阵A中所有序号非c_j的列置为零后形成的矩阵设为A_cal，j＝1，2，...，n_c，将矩阵A中所有序号非e_l的列置为零后形成的矩阵设为A_ext，l＝1，2，...，n_e，将向量K中所有序号非c_j的元素置为零后形成的向量设为K_cal，j＝1，2，...，n_c，将向量K中所有序号非e_l的元素置为零后形成的向量设为K_ext，l＝1，2，...，n_e，则Δ_g＝AK可以记为： 

${\mathrm{\Δ}}_g=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{a}_i=\underset{j=1}{\overset{n_c}{\mathrm{\Σ}}}{k}_{c_j}{a}_{c_j}+\underset{l=1}{\overset{n_e}{\mathrm{\Σ}}}{k}_{e_l}{a}_{e_l}=A_{\mathrm{cal}}{K}_{\mathrm{cal}}+A_{\mathrm{ext}}{K}_{\mathrm{ext}}$

然后将A_cal去掉元素全为零的行和列得到A_ca ln z，并记下这些行和列的序号，接着求A_ca ln z的最小二乘逆矩阵并对扩充元素全为零的行和列得到其中，扩充的行序号与A_cal中全为零的列序号相同，列序号与A_cal中全为零的行序号相同， 

由下式求解参与标定的误差参数： 

计算得到K_cal，其中，参与标定计算的误差参数序号上元素为相应的标定计算值，外部给出的误差参数序号上元素全为零，将K_cal与外部给出的K_ext相加得到一阶误差参数K； 

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

当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，3，4时，的计算方法如下： 

$\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为纬度， 

ΔΦ^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，3，由i＝0，1，2，3，4； 

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

$\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(3\right)}\\ {\mathrm{\ω}}_v^{*\left(4\right)}\end{array}\right]=\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ B^{\left(3\right)}\\ {B^{\left(4\right)}}^{}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}$

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

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

在上述技术方案中，在步骤一中，所述惯性测量单元初始位置保证Y轴朝向与东西向夹角为15～90°。 

在上述技术方案中，当步骤一中惯性测量单元初始位置的Y轴朝向与东西向夹角为0～15°时，步骤八中的B接近奇异，此时二阶误差参数ω中的B_gy项可由外部给出，然后按照步骤五的方法求取B的最小二乘逆矩阵，进而得到ω中的另外两个需标定的误差参数B_gx、B_gz。 

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

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

标定旋转顺序 

旋转序号	旋转过程
		1	-90Y
2	-90Y
		3	270Y
4	-90Y

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

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

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

本方法原理描述如下： 

标定方法利用采集的原始数据，在第i(i＝0，1，2，3)位置上进行初始对准，然后在第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{g{T}^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{g{T}^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作为参考矢量)，所以不用对北，适用于在低精度单轴转位设 备上进行标定，同时标定旋转只需通过编排5个位置，并可在给定历史标定参数的条件下经过多次迭代运算求解出部分所需标定的参数(这是考虑到部分惯性测量单元误差参数的长期稳定性较好，或者该部分惯性测量单元误差参数对惯性导航精度影响较小，采用单轴转位设备标定惯性测量单元主要的误差参数是可取的)，可以大幅度降低标定成本，缩短标定时间； 

(2)该标定方法采用迭代算法，可以大幅降低标定对转台精度的依赖性。 

与国内外其他系统级标定方法相比，免对北单轴系统级标定方法不需对北，且只需在单轴转位设备上标定，甚至在不依赖单轴转位设备的情况下，直接在大理石平台上进行手动翻转标定即可，这样惯性测量组合就可以直接在外场标定，标定场所不必局限于实验室。这种外场标定将给工程人员带来极大的便利，同时彻底摆脱了标定对转台的依赖。另外，相对于双轴转位设备而言，本发明不仅对器材和采集的原始数据要求低，而且单轴转位设备费用较低，大幅度降低了标定成本，意义重大。 

附图说明

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

具体实施方式

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

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

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

表一 标定旋转顺序 

旋转序号	旋转过程
		1	-90Y
2	-90Y
		3	270Y
4	-90Y

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

步骤三：利用步骤一采集的惯性测量单元数据，在第i位置上使用双矢量定姿法进行初始对准，i＝0，1，2，3，然后在第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，3，所述包含和包含和所述和分别为第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^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}}\\ g_{\mathrm{gy}}\\ g_{\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，3，中i＝1，2，3，4，即： 

$\left\{\begin{array}{c}{\mathrm{\Δ}}_g^{\left(i\right)}=\frac{-{\left[{\mathrm{\δf}}^{n(i-1)}\right]}_H+f^n\×{\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{\Δ\Φ}}^{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(3\right)}\\ {\mathrm{\Δ}}_g^{\left(4\right)}\end{array}\right]}_{15\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ A^{\left(3\right)}\\ A^{\left(4\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程： 

Δ_g＝AK 

步骤五：根据一阶中间参数与一阶误差参数的关系式，设A＝[a₁ a₂ … a_n]、K＝[k₁ k₂ … k_n]^T，其中，a_i为矩阵A的列向量，i＝1，2…n，k_i为向量K的各元素，i＝1，2…n，[k₁ k₂ … k_n]^T为行向量[k₁ k₂ … k_n]的转置， 

K中包含B_ay、K_ayy、K_azy、K_ax2、K_ay2、K_az2、K_gxx、K_gzz、K_gxy、K_gxz、K_gyx、K_gyz、K_gzx的13个误差参数直接由最近一次的完整标定结果给出，这些由外部给出的误差参数在向量K中的序号为e_l，l＝1，2，...，n_e，其余参与标定计算的误差参数在向量K中的序号为c_j，j＝1，2，...，n_c，且n_c+n_e＝n；然后将矩阵A中所有序号非c_j的列置为零后形成的矩阵设为A_cdl，j＝1，2，...，n_c，将矩阵A中所有序号非e_l的列置为零后形成的矩阵设为A_ext，l＝1，2，...，n_e，将向量K中所有序号非c_j的元素置为零后形成的向量设为K_cal，j＝1，2，...，n_c，将向量K中所有序号非e_l的元素置为零后形成的向量设为K_ext，l＝1，2，...，n_e，则Δ_g＝AK可以记为： 

${\mathrm{\Δ}}_g=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{a}_i=\underset{j=1}{\overset{n_c}{\mathrm{\Σ}}}{k}_{c_j}{a}_{c_j}+\underset{l=1}{\overset{n_e}{\mathrm{\Σ}}}{k}_{e_l}{a}_{e_l}=A_{\mathrm{cal}}{K}_{\mathrm{cal}}+A_{\mathrm{ext}}{K}_{\mathrm{ext}}$

然后将A_cal去掉元素全为零的行和列得到A_ca ln z，并记下这些行和列的序号，接着求A_ca ln z的最小二乘逆矩阵并对扩充元素全为零的行和列得到其中，扩充的行序号与A_cal中全为零的列序号相同，列序号与A_cal中全为零的行序号相同， 

由下式求解参与标定的误差参数： 

计算得到K_cal，其中，参与标定计算的误差参数序号上元素为相应的标定计算值，外部给出的误差参数序号上元素全为零，将K_cal与外部给出的K_ext相加得到一阶误差参数K。 

步骤六：利用步骤五中的K计算出每个静止位置i对应的补偿分量其中i＝0，1，2，3，4，计算方法如下： 

当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，3，4时，的计算方法如下： 

$\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{\Δ}\text{\Φ}}^{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为纬度， 

ΔΦ^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，3，  ${\mathrm{\ω}}_v^{*\left(i\right)}={\mathrm{\ω}}_v^{\left(i\right)}-{\mathrm{\ω}}_v^{0\left(i\right)}$ 中i＝0，1，2，3，4。 

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

$\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(3\right)}\\ {\mathrm{\ω}}_v^{*\left(4\right)}\end{array}\right]=\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ B^{\left(3\right)}\\ B^{\left(4\right)}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}$

步骤八：利用步骤七中的联立方程计算B的最小二乘逆矩阵进而通过计算出二阶误差参数ω。另外，当步骤一中惯性测量单元初始位置的Y轴朝向与东西向夹角为0～15°时，B接近奇异，此时二阶误差参数ω中的B_gy项可由外部给出，然后按照步骤五的方法求取B的最小二乘逆矩阵进而得到ω中的另外两个需标定的误差参数B_gx、B_gz。 

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

Claims (8)

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

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

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

步骤三：利用步骤一采集的惯性测量单元原始数据，在第i位置上使用双矢量定姿法进行初始对准，i＝0，1，2，3，然后在第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，3，所述包含和包含和所述和分别为第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]+\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}\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{{[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){]}_V}{g}$

其中，

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

g是重力加速度，

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

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

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

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

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

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

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

用步骤三得到的构建方程其中中i＝0，1，2，3，中i＝1，2，3，4，即：

$\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}_{0(i-1~i)({\mathrm{\ω}}_0^b+K_g{\mathrm{\ω}}^b)\mathrm{dt}}^b\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(3\right)}\\ {\mathrm{\Δ}}_g^{\left(4\right)}\end{array}\right]}_{15\×1}=\left[\begin{array}{c}{A}^{\left(0\right)}\\ A^{\left(1\right)}\\ A^{\left(2\right)}\\ A^{\left(3\right)}\\ A^{\left(4\right)}\end{array}\right]K=\mathrm{AK}$

最终构建如下方程：

Δ_g＝AK

步骤五：根据一阶中间参数与一阶误差参数的关系式，设A＝[a₁ a₂ … a_n]、K＝[k₁ k₂ … k_n]^T，其中，a_i为矩阵A的列向量，i＝1，2…n，k_i为向量K的各元素，i＝1，2…n，[k₁ k₂ … k_n]^T为行向量[k₁ k₂ … k_n]的转置，

K中包含B_ay、K_ayy、K_azy、K_ax2、K_ay2、K_az2、K_gxx、K_gzz、K_gxy、K_gxz、K_gyx、K_gyz、K_gzx的13个误差参数直接由最近一次的完整标定结果给出，这些由外部给出的误差参数在向量K中的序号为e_l，l＝1，2，...，n_e，其余参与标定计算的误差参数在向量K中的序号为c_j，j＝1，2，...，n_c，且n_c+n_e＝n；然后将矩阵A中所有序号非c_j的列置为零后形成的矩阵设为A_cal，j＝1，2，...，n_c，将矩阵A中所有序号非e_l的列置为零后形成的矩阵设为A_ext，l＝1，2，...，n_e，将向量K中所有序号非c_j的元素置为零后形成的向量设为K_cal，j＝1，2，...，n_c，将向量K中所有序号非e_l的元素置为零后形成的向量设为K_ext，l＝1，2，...，n_e，则Δ_g＝AK可以记为：

${\mathrm{\Δ}}_g=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{a}_i=\underset{j=1}{\overset{n_c}{\mathrm{\Σ}}}{k}_{c_j}{a}_{c_j}+\underset{l=1}{\overset{n_e}{\mathrm{\Σ}}}{k}_{e_l}{a}_{e_l}=A_{\mathrm{cal}}{K}_{\mathrm{cal}}+A_{\mathrm{ext}}{K}_{\mathrm{ext}}$

然后将A_cal去掉元素全为零的行和列得到A_ca ln z，并记下这些行和列的序号，接着求A_ca ln z的最小二乘逆矩阵并对扩充元素全为零的行和列得到其中，扩充的行序号与A_cal中全为零的列序号相同，列序号与A_cal中全为零的行序号相同，

由下式求解参与标定的误差参数：

计算得到K_cal，其中，参与标定计算的误差参数序号上元素为相应的标定计算值，外部给出的误差参数序号上元素全为零，将K_cal与外部给出的K_ext相加得到一阶误差参数K；

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

当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，3，4时，的计算方法如下：

$\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{\Δ\Φ}}^{n(i-1~i)}\begin{array}{ccc}& & \end{array}\\ {\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{\Δ\Φ}}^{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为纬度，

ΔФ^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，3，中i＝0，1，2，3，4；

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

$\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(3\right)}\\ {\mathrm{\ω}}_v^{*\left(4\right)}\end{array}\right]=\left[\begin{array}{c}{B}^{\left(0\right)}\\ B^{\left(1\right)}\\ B^{\left(2\right)}\\ B^{\left(3\right)}\\ B^{\left(4\right)}\end{array}\right]\mathrm{\ω}=\mathrm{B\ω}$

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

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

2.根据权利要求1所述的一种适用于低精度无方位基准单轴转位设备的惯性测量单元标定方法，其特征在于：在步骤一中，所述惯性测量单元初始位置保证Y轴朝向与东西向夹角为15～90°。

3.根据权利要求1所述的一种适用于低精度无方位基准单轴转位设备的惯性测量单元标定方法，其特征在于：当步骤一中惯性测量单元初始位置的Y轴朝向与东西向夹角为0～15°时，步骤八中的B接近奇异，此时二阶误差参数ω中的B_gy项可由外部给出，然后按照步骤五的方法求取B的最小二乘逆矩阵进而得到ω中的另外两个需标定的误差参数B_gx、B_gz。

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

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

标定旋转顺序

旋转序号 旋转过程 1 -90Y 2 -90Y 3 270Y 4 -90Y

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

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

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