CN103852085B - 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法 - Google Patents

Abstract

本发明公开了一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，属于惯性技术领域。本发明采用9次翻转路径设计，对系统输出惯性数据应用最小二乘拟合的方法求解得到光纤捷联惯导系统18项误差系数；再通过一次翻转得到零偏误差，进而得到全部21项器件误差参数；本发明提供的技术方案利用六面体或其它相似的可翻转装置即可完成现场标定试验，克服了传统实验室标定的不足，提高了系统实际使用精度。

Description

一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法

技术领域

本发明属于惯性技术领域，涉及一种光纤捷联惯导系统现场标定方法，具体地说，是指一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法。

背景技术

光纤陀螺具有精度高、启动快、动态范围大、抗振动冲击能力强及成本低等优点，是惯性仪表领域的发展趋势。近年来，光纤陀螺技术的迅猛发展推动了光纤捷联惯导系统在陆、海、空、天领域的应用。光纤捷联惯导在使用前必须通过实验室转台标定试验确定出其核心部件即光纤陀螺和加速度计的标度因数和各项误差系数，在后续的导航计算中进行补偿。

但是，通过实验室转台试验标定出的光纤惯导系统各项误差系数并不是固定不变的，包括光纤陀螺零位漂移误差、标度因数、安装误差和加速度计常值误差、标度困数、安装误差等。这些误差参数随着系统使用或存放时间的推移而变化，尤其是光纤陀螺零位漂移和加速度计偏值，每次上电启动都不相同，且时间间隔越长变化越大。光纤陀螺和加速度计误差参数的改变直接导致光纤惯导系统的精度降低，使惯导系统无法使用要求。

因此，通常需要对光纤惯导系统进行半年或三个月一次的定期标定，而且传统的基于精密转台的标定方法复杂、耗时长，这为使用单位增添了巨大的工作量及人力、物力、财力的消耗。因此，在使用现场对光纤惯导的各项误差系数进行标定，不仅可以减少甚至取消定期标定，还可以提高光纤捷联惯导的实际使用精度。但是，在现场没有精密的转台作为测试基准，不能对充纤陀螺捷联惯导进行精确定向，所以传统的基于精密转台的静态多位置标定方法和速率标定方法都无法实施，为光纤捷联惯导的现场标定带来了极大难度。

在没有精确的基准设备和复杂的转位机构的情况下，现场标定就不能采用传统的速率位置法，必须依靠系统级方法来克服基准信息缺乏的困难。系统级方法以光纤惯导的导航误差为观测量，通过惯导系统误差参数与导航误差之间的关系，建立惯导系统误差参数与导航误差之间的方程式，进而求解出惯导系统的误差参数。因此，系统级标定方法可以克服使用现场设备条件不足的缺点，并获得高精度的系统误差标定接过。

参考文献[1]（光电工程，刘百奇,房建成.光纤陀螺IMU的六位置旋转现场标定新方法[J].35(1),2008:60-65）公开了一种光纤陀螺IMU的六位置旋转现场标定新方法，本文采用光纤陀螺IMU在六个位置上进行12次旋转，然后根据光纤陀螺IMU的误差模型建立42个非线性输入输出方程，通过旋转积分和对称位置误差相消，消除方程中的非线性项，最终求出陀螺标度因数、陀螺常值漂移、陀螺安装误差和加速度计常值偏置等15个误差系数。但是该方法不能够标定出加速度计通道的标度因数和安装误差。

参考文献[2]（测控技术，颜开思,李岁劳,龚柏春,贾继超.[J].30(5),2011:106-109）公开了一种基于平台和正六面体的惯导系统现场标定技术，该文献通过翻转正六面体使对称位置误差相消，并且在对准中获取姿态信息，同时精确标定出陀螺漂移和加速度计零偏。最后对理论分析结果进行了仿真验证，仿真结果表明该方案可以实现外场条件下的陀螺漂移和加速度计零偏的精确标定。但是该方法不能够标定出陀螺、加速度的标度因数和安装误差。

参考文献[3]（吴赛成,秦石乔,王省书,胡春生.[J].中国惯性技术学报,19(2),2011:185-189）公开了一种激光陀螺惯性测量单元系统级标定方法，该文献建立了附加约束条件的陀螺和加速度计安装坐标系数学模型，根据陀螺和加速度计的输出误差方程，从惯性导航基本误差方程出发推导了惯性测量单元的系统级误差参数标定Kalman滤波模型，该模型包含了陀螺和加速度计零偏、比例因子、安装误差在内共21维标定误差状态变量，且仅以速度解算误差为观测量。但是该方法标定步骤较多，标定时间过长，缺少实验数据验证。

参考文献[4]（专利文献号CN102607594A）公开了一种捷联惯导光纤陀螺系统误差参数现场标定方法，所述现场标定方法通过姿态测量仪器给出载体姿态角，选取姿态作为观测量，标定出光纤捷联惯导系统光纤陀螺各项误差系数。但是该方法需要现场提供姿态测量辅助器件，实时精确测量载体姿态角，并且要与光纤陀螺输出值保持一致。

发明内容

本发明的目的在于提供一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，减少或者甚至取消光纤捷联惯导系统周期性地返实验室校标，提高系统实际使用精度。

本发明采用9次翻转路径设计，包括转动轴、转动顺序和转动角度等；应用最小二乘拟合的方法得到光纤捷联惯导系统全部21项器件误差参数；利用六面体或其它相似的可翻转装置即可完成现场标定试验。具体方法步骤如下：

第一步：将光纤捷联惯导系统通过工装安装在正六面体装置上，锁紧。连接系统、电源和采集计算机之间的线缆，并检查正确。

第二步：将正六面体装置置于水平面上，上电预热使光纤捷联惯导系统达到热平衡状态，并装订光纤捷联惯导系统的初始位置参数，包括初始的经度、纬度和高度。

第三步：采用“静止-转动-静止”进行手动翻转正六面体装置，按照转动路径序列完成9次翻转。转动前后每个位置静止3～5min，并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据。

第四步：对惯性器件误差进行建模。

第五步：对9组惯性器件数据进行处理，采用最小二乘拟合的方法求解，得到除陀螺零偏以外的其他18项误差系数；

第六步：根据标定得到18项误差系数，对保存的9组惯性器件数据进行补偿，得到新的9组惯性器件数据。

第七步：对于第六步中得到的新的9组惯性器件数据，重复第五步和第六步3～5次，即迭代器件18项误差系数值进行重复最小二乘拟合计算，直至收敛。累加每次迭代估计得到器件18项误差系数值，即为最终标定参数值。

第八步：根据上述标定的18项误差系数，对光纤捷联惯导系统进行补偿。

第九步：重复第一步到第二步，按照序列1进行手动手动翻转六面体依次，转动前后分别静止时间30min，保存陀螺仪输出数据，基于该输出数据，计算陀螺零偏误差。

第十步：根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿，采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可，完成了光纤捷联惯导系统21项误差参数的现场标定。

本发明的有益效果在于：

本发明所提出的方法可以在现场完成光纤捷联惯导系统21项误差参数的标定，克服了传统实验室标定的不足，提高了系统实际使用精度。

附图说明

图1为本发明提供的基于最小二乘拟合的光纤捷联惯导系统现场标定方法流程图；

图2A和图2B分别为本发明实施例中静态和摇摆情况下现场标定补偿前后20min导航定位误差对比曲线。

具体实施方式

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

本发明提供一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，如图1所示，所述方法包括如下步骤：

第一步：将光纤捷联惯导系统通过工装安装在正六面体装置上，锁紧。连接系统、电源和采集计算机之间的线缆，并检查正确。

第二步：将正六面体装置置于水平面上，上电预热使光纤捷联惯导系统达到热平衡状态，并装订光纤捷联惯导系统的初始位置参数，包括初始的经度、纬度和高度。

第三步：按照表1转动路径序列，表1中转动轴X、Y、Z，光纤捷联惯导系统初始姿态为0时，光纤捷联惯导系统XYZ轴与导航坐标系东北天位置重合。

采用“静止-转动-静止”进行手动翻转正六面体装置，完成9次翻转，转动角允许存在±10°误差。转动前后每个位置静止3～5min，并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据。

表1最小二乘拟合法转动路径序列（转动角单位：度）

第四步：对惯性器件误差进行建模，包括光纤陀螺误差模型和加速度计误差模型，分别如下：

${{\mathrm{\δ\ω}}_{\mathrm{ib}}}^b=\left[\begin{array}{c}{{\mathrm{\δ\ω}}_{\mathrm{ibx}}}^b\\ {{\mathrm{\δ\ω}}_{\mathrm{iby}}}^b\\ {{\mathrm{\δ\ω}}_{\mathrm{ibz}}}^b\end{array}\right]=\left[\begin{array}{c}{\mathrm{gB}}_x\\ {\mathrm{gB}}_y\\ {\mathrm{gB}}_z\end{array}\right]+\left[\begin{array}{ccc}{\mathrm{gSF}}_x& {\mathrm{gMA}}_{\mathrm{xy}}& {\mathrm{gMA}}_{\mathrm{xz}}\\ {\mathrm{gMA}}_{\mathrm{yx}}& {\mathrm{gSF}}_y& {\mathrm{gMA}}_{\mathrm{yz}}\\ {\mathrm{gMA}}_{\mathrm{zx}}& {\mathrm{gMA}}_{\mathrm{zy}}& {\mathrm{gSF}}_z\end{array}\right]\left[\begin{array}{c}{{\mathrm{\ω}}_{\mathrm{ibx}}}^b\\ {{\mathrm{\ω}}_{\mathrm{iby}}}^b\\ {{\mathrm{\ω}}_{\mathrm{ibz}}}^b\end{array}\right]$

${\mathrm{\δf}}_{\mathrm{ib}}^b=\left[\begin{array}{c}{\mathrm{\δf}}_{\mathrm{ibx}}^b\\ {\mathrm{\δf}}_{\mathrm{iby}}^b\\ {\mathrm{\δf}}_{\mathrm{ibz}}^b\end{array}\right]=\left[\begin{array}{c}{\mathrm{aB}}_x\\ {\mathrm{aB}}_y\\ a{B}_z\end{array}\right]+\left[\begin{array}{ccc}{\mathrm{aSF}}_x& 0& 0\\ {\mathrm{aMA}}_{\mathrm{yx}}& {\mathrm{aSF}}_y& 0\\ {\mathrm{aMA}}_{\mathrm{zx}}& {\mathrm{aMA}}_{\mathrm{zy}}& {\mathrm{aSF}}_z\end{array}\right]\left[\begin{array}{c}{f}_{\mathrm{ibx}}^b\\ f_{\mathrm{iby}}^b\\ f_{\mathrm{ibz}}^b\end{array}\right]$

式中δω_ib ^b表示陀螺仪的误差输出矢量；δω_ibx ^b、δω_iby ^b、δω_ibz ^b分别表示由三轴陀螺误差引起的误差角速度。ω_ibx ^b、ω_iby ^b、ω_ibz ^b分别表示三轴陀螺测量值；gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间安装误差角；gB_x、gB_y、gB_z分别表示三轴陀螺仪零偏误差；

式中表示加速度计的误差输出；δf_ibx ^b、δf_iby ^b、δf_ibz ^b分别表示由三轴加速度计误差引起的误差加速度；f_ibx ^b、f_iby ^b、f_ibz ^b分别表示三轴加速度计测量值；aSF_x、aSF_y、aSF_z分别为三轴加速度计的标度因数误差；aB_x、aB_y、aB_z分别为三轴加速度计的零偏；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角；

第五步：对保存的9组惯性器件数据进行处理，采用最小二乘拟合的方法求解，得到除陀螺零偏以外的其他18项误差系数。

所述最小二乘拟合法包含以下几个步骤：

步骤1：建立器件误差与系统速度误差一阶导数变化量的数学模型。

K·X＝A

其中K表示转动系数矩阵，X表示误差向量，A表示观测矩阵。误差向量X＝[X_a;X_g]；

X_a＝[aB_x aB_y aB_z aSF_x aMA_yx aSF_y aMA_zx aMA_zy aSF_z]^T

X_g＝[gSF_x gMA_xy gMA_xz gMA_yx gSF_y gMA_yz gMA_zx gMA_zy gSF_z]^T

式中gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间的安装误差角；aSF_x、aSF_y、aSF_z分别为三轴加速度计标度因数误差；aB_x、aB_y、aB_z分别为三轴加速度计零偏；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角；

步骤2：求解观测矩阵A。

观测矩阵表示T₂时刻速度误差一阶导数矢量，表示T₁时刻速度误差一阶导数矢量。

观测矩阵由三向速度误差一阶导数组成，光纤捷联惯导系统静止时，速度误差即为光纤捷联惯导系统导航输出的速度值。采用标准卡尔曼滤波计算速度误差一阶导数值。

建立状态方程： $\frac{d}{\mathrm{dt}}\left[\begin{array}{c}{\mathrm{\δV}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·}{V}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·\·}{V}}_i\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\δV}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·}{V}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·\·}{V}}_i\left(t\right)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ w_i\end{array}\right]$

建立量测方程： $Z_i=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\δV}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·}{V}}_i\left(t\right)\\ \mathrm{\δ}{\stackrel{\·\·}{V}}_i\left(t\right)\end{array}\right]+{\mathrm{\μ}}_i$

其中，Z_i表示观测矢量；δV_i(t)＝V_i(t),i＝x,y,z，δV_i(t)表示速度误差矢量，V_i(t)表示导航解算速度矢量，w_i和μ_i分别表示状态噪声和观测噪声。采用标准卡尔曼滤波方程进行迭代，滤波完毕即可得到速度误差一阶导数的估计值求取转动前后两次速度误差一阶导数并作差，即可得到系统观测阵A。

步骤3：求解转动系数矩阵K。

转动系数矩阵 $K=\left[\begin{array}{cc}{\mathrm{df}}_x& -g\·{\mathrm{d\φ}}_y\\ {\mathrm{df}}_y& g\·{\mathrm{d\φ}}_x\\ {\mathrm{df}}_z& 0_{1\×9}\end{array}\right]$

式中，df＝[df_x df_y df_z]^T表示加速度计误差系数阵，dφ＝[dφ_x dφ_y dφ_z]^T表示陀螺误差系数阵。

$d_f=C_b^n\left(T_2\right)M_2-C_b^n\left(T_1\right)M_1,$ 其中，

$M_i=\left[\begin{array}{ccccccccc}1& 0& 0& f_{\mathrm{ibx}}^b\left(T_i\right)& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& f_{\mathrm{ibx}}^b\left(T_i\right)& f_{\mathrm{iby}}^b\left(T_i\right)& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& f_{\mathrm{ibx}}^b\left(T_i\right)& f_{\mathrm{iby}}^b\left(T_i\right)& f_{\mathrm{ibz}}^b\left(T_i\right)\end{array}\right],$ 其中i＝1,2；

$d_{\mathrm{\φ}}=-{\∫}_{T_1}^{T_2}{C}_{b\left(T_1\right)}^n\·C_{b\left(t\right)}^{b\left(T_1\right)}\·\left[\begin{array}{ccc}{\mathrm{\ω}}_x\left(t\right)& 0& 0\\ 0& {\mathrm{\ω}}_y\left(t\right)& 0\\ 0& 0& {\mathrm{\ω}}_z\left(t\right)\end{array}\right]M_g\·\mathrm{dt}$

其中，f_ibx ^b(T_i)、f_iby ^b(T_i)、f_ibz ^b(T_i)分别表示T_i时刻三轴加速度计测量值；为惯组初始时刻T₁的姿态矩阵，ω_x(t)、ω_y(t)、ω_z(t)为t时刻的转动角速度，而为t时刻b系到T₁时刻b系的转换矩阵，M_g为系数矩阵，这两者与转动轴向有关；所述b系是指载体坐标系，n代表n系，是指导航坐标系，本发明中定义导航坐标系为东北天，即当地地理坐标系。

若光纤捷联惯导系统绕X轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left({\mathrm{\ω}}_x\left(t\right)\right)& \sin \left({\mathrm{\ω}}_x\left(t\right)\right)\\ 0& -\sin \left({\mathrm{\ω}}_x\left(t\right)\right)& \cos \left({\mathrm{\ω}}_x\left(t\right)\right)\end{array}\right],$ $M_g=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right]$

若光纤捷联惯导系统绕Y轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}|=\left[\begin{array}{ccc}\cos \left({\mathrm{\ω}}_y\left(t\right)\right)& 0& -\sin \left({\mathrm{\ω}}_y\left(t\right)\right)\\ 0& 1& 0\\ \sin \left({\mathrm{\ω}}_y\left(t\right)\right)& 0& \cos \left({\mathrm{\ω}}_y\left(t\right)\right)\end{array}\right],$ $M_g=\left[\begin{array}{ccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\end{array}\right]$

若光纤捷联惯导系统绕Z轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}|=\left[\begin{array}{ccc}\cos \left({\mathrm{\ω}}_z\left(t\right)\right)& \sin \left({\mathrm{\ω}}_z\left(t\right)\right)& 0\\ -\sin \left({\mathrm{\ω}}_z\left(t\right)\right)& \cos \left({\mathrm{\ω}}_z\left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right],$ $M_g=\left[\begin{array}{ccccccccc}0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

步骤4：采用最小二乘方法求解误差向量：X＝K^T·(KK^T)^-1·A。由于转动系数矩阵K的秩rank(K)＝15，得到的误差系数有15项，包括三轴陀螺仪标度因数误差gSF_x、gSF_y、gSF_z，各轴陀螺仪间安装误差角gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy，三轴加速度计零偏aB_x、aB_y、aB_z；各轴加速度计间安装误差角aMA_yx、aMA_zx、aMA_zy；

步骤5：求解三轴加速度计标度因数；

采用第1次转动数据： ${\mathrm{aSF}}_z=\frac{1}{2(T_1+T_2)g}{(\underset{t=0}{\overset{t=T_2}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right)+\underset{t=0}{\overset{t=T_1}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right))}_1$

采用第4次转动数据： ${\mathrm{aSF}}_y=\frac{1}{2(T_1+T_2)g}{(\underset{t=0}{\overset{t=T_2}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right)+\underset{t=0}{\overset{t=T_1}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right))}_4$

采用第7次转动数据： ${\mathrm{aSF}}_x=\frac{1}{2(T_1+T_2)g}{(\underset{t=0}{\overset{t=T_2}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right)+\underset{t=0}{\overset{t=T_1}{\mathrm{\Σ}}}\mathrm{\δ}{\stackrel{\·}{V}}_z\left(t\right))}_7$

式中，表示转动前从0时刻到T₁时刻的时间内天向速度误差一阶导数之和，表示转动后从0时刻到T₂时刻的时间内天向速度误差一阶导数之和。

第六步：根据上述步骤4和步骤5中标定得到18项误差系数，对保存的9组惯性器件数据进行补偿，得到新的9组惯性器件数据。其中补偿公式如下：

$\left[\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right]=(I+\left[\begin{array}{ccc}{\mathrm{gSF}}_x& {\mathrm{gMA}}_{\mathrm{xy}}& {\mathrm{gMA}}_{\mathrm{xz}}\\ {\mathrm{gMA}}_{\mathrm{yx}}& {\mathrm{gSF}}_y& {\mathrm{gMA}}_{\mathrm{yz}}\\ {\mathrm{gMA}}_{\mathrm{zx}}& {\mathrm{gMA}}_{\mathrm{zy}}& {\mathrm{gSF}}_z\end{array}\right])\left[\begin{array}{c}{{\mathrm{\ω}}_{\mathrm{ibx}}}^b\\ {{\mathrm{\ω}}_{\mathrm{iby}}}^b\\ {{\mathrm{\ω}}_{\mathrm{ibz}}}^b\end{array}\right]$

$\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]=(I+\left[\begin{array}{ccc}{\mathrm{aSF}}_x& 0& 0\\ {\mathrm{aMA}}_{\mathrm{yx}}& {\mathrm{aSF}}_y& 0\\ {\mathrm{aMA}}_{\mathrm{zx}}& {\mathrm{aMA}}_{\mathrm{zy}}& {\mathrm{aSF}}_z\end{array}\right])\left[\begin{array}{c}{f}_{\mathrm{ibx}}^b-{\mathrm{aB}}_x\\ f_{\mathrm{iby}}^b-{\mathrm{aB}}_y\\ f_{\mathrm{ibz}}^b-{\mathrm{aB}}_z\end{array}\right]$

式中gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差标定结果；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间安装误差角标定结果；aSF_x、aSF_y、aSF_z分别为三轴加速度计标度因数误差标定结果；aB_x、aB_y、aB_z分别为三轴加速度计零偏标定结果；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角标定结果；ω_ibx ^b、ω_iby ^b、ω_ibz ^b分别表示三轴陀螺原测量值；f_ibx ^b、f_iby ^b、f_ibz ^b分别表示三轴加速度计原测量值；I表示三阶单位阵；w_x、w_y、w_z分别表示三轴陀螺补偿后测量值；f_x、f_y、f_z分别表示三轴加速度计补偿后测量值。

第七步：重复第五步和第六步3～5次，即迭代器件18项误差系数值进行重复最小二乘拟合计算，直至收敛。累加每次迭代估计得到器件18项误差系数值，即为最终标定参数值。

第八步：根据上述最终标定的18项误差系数值，按照第六步中补偿公式对光纤捷联惯导系统进行补偿。

第九步：重复第一步到第二步，按照表1的序列1进行手动翻转六面体，转动前后分别静止时间30min，保存陀螺仪输出数据。基于输出数据，计算陀螺零偏误差gB＝[gB_x gB_ygB_z]^T。计算公式如下：

$\mathrm{gB}=\frac{1}{2}\{({\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_1\right)+{\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_2\right))-{[C_b^n\left(T_1\right)+C_b^n\left(T_1\right)]}^T(I-\mathrm{\φ}\×){\mathrm{\ω}}_{\mathrm{ie}}^n\}$

其中，φ＝[0 0 φ_z]^T，φ_z表示航向角误差， ${\mathrm{\φ}}_z=\frac{{\mathrm{\Δ\ω}}_{\mathrm{ibx}}^b-{\mathrm{\ΔC}}_b^n\left(\mathrm{1,2}\right){\mathrm{\ω}}_{\mathrm{iey}}^n-{\mathrm{\ΔC}}_b^n\left(\mathrm{1,3}\right){\mathrm{\ω}}_{\mathrm{iez}}^n}{{\mathrm{\ΔC}}_b^n\left(\mathrm{1,1}\right){\mathrm{\ω}}_{\mathrm{iey}}^n},$

${\mathrm{\Δ\ω}}_{\mathrm{ib}}^b={\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_2\right)-{\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_1\right),$ ${\mathrm{\ΔC}}_b^n=C_b^n\left(T_2\right)-C_b^n\left(T_1\right)$

式中分别表示T₁、T₂时刻三轴陀螺测量矢量，表示转动前后两时刻陀螺测量误差矢量、表示第一行分量，分别表示T₁、T₂时刻载体系b系到导航系n系的方向余弦矩阵，表示转动前后两时刻方向余弦矩阵误差，表示第i行第j列元素值，表示导航系下地球速度矢量，分别表示导航系下地球速度矢量的北向分量和天向分量。

第十步：根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿，采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可，完成了光纤捷联惯导系统21项误差参数的现场标定。

实施例

第一步：将光纤捷联惯导系统通过工装安装在正六面体装置上，锁紧。连接系统、电源、采集计算机之间的线缆，并检查正确。

第二步：将正六面体装置置于水平面上，上电预热使光纤捷联惯导系统达到热平衡状态，并装订光纤捷联惯导系统的初始位置参数，包括初始的经度、纬度和高度。

第三步：按照表1转动路径，采用“静止-转动-静止”进行手动翻转正六面体装置，完成9次翻转，转动角允许存在±10°误差。转动前后每个位置静止3～5min，并保存9次转动数据。

第四步：对保存的9组数据进行处理，采用最小二乘拟合的方法求解，并迭代5次，得到除陀螺零偏以外的其他18项误差参数。

第五步：先补偿18项误差参数，重复第一步到第二步，按照表1的序列1进行手动翻转六面体，转动前后分别静止时间30min，保存陀螺输出数据。

第六步：基于第五步数据，计算三轴陀螺零偏误差值，并对光纤捷联惯导系统陀螺输出进行补偿，完成光纤捷联惯导系统21项误差参数的现场标定。

第七步：将系统断电，1天后重新启动光纤捷联惯导系统，首先静态采集23min惯性器件数据，然后再将系统断电3-5h，将系统安装于双轴转台上，先静态3min再摇摆20min，一共采集23min惯性器件数据。

第八步：离线处理两组惯性器件数据，将第五步得到的现场标定结果对数据进行补偿，采用解析式粗对准3min和纯惯性导航，对比补偿前后两组数据的纯惯性导航结果。

结果分析：

（1）采用最小二乘拟合法迭代5次得到的18项误差参数值，和最后补偿得到的三轴陀螺零偏值如表2所示。从表2中可以看出，经过5次迭代，参数估计值渐近收敛，其中三轴加速度计零偏误差值收敛至60μg以内，加速度计标度因数收敛值在5ppm以内，陀螺标度因数收敛值在35ppm以内，各轴间安装误差收敛至14″以内。

表2最小二乘拟合法现场标定实验结果

（2）对比现场标定补偿前后的数据导航结果如图2A和图2B所示。图2A是20min静态水平定位误差对比曲线，图2B是摇摆情况下水平定位误差对比曲线。从图2A和图2B中可以看出，不管是静态还是动态情况下，系统的导航定位误差减小了1倍以上，因此采用本发明提供的现场标定方法补偿后的数据精度更高。

Claims (5)

1.一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，其特征在于通过如下步骤实现：

第一步：将光纤捷联惯导系统通过工装安装在正六面体装置上，锁紧；连接系统、电源和采集计算机之间的线缆，并检查正确；

第二步：将正六面体装置置于水平面上，上电预热使光纤捷联惯导系统达到热平衡状态，并装订光纤捷联惯导系统的初始位置参数，包括初始的经度、纬度和高度；

第三步：采用“静止-转动-静止”进行手动翻转正六面体装置，按照转动路径序列完成9次翻转；转动前后每个位置静止3～5min，并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据；所述转动路径序列如下：

第四步：对惯性器件误差进行建模；

第五步：对保存的9组惯性器件数据进行处理，采用最小二乘拟合法求解，得到除陀螺零偏以外的其他18项误差系数；所述最小二乘拟合法包含以下几个步骤：

步骤1：建立器件误差与系统速度误差一阶导数变化量的数学模型：

K·X＝A

其中K表示转动系数矩阵，X表示误差向量，A表示观测矩阵；误差向量X＝[X_a；X_g]；

X_a＝[aB_x aB_y aB_z aSF_x aMA_yx aSF_y aMA_zx aMA_zy aSF_z]^T

X_g＝[gSF_x gMA_xy gMA_xz gMA_yx gSF_y gMA_yz gMA_zx gMA_zy gSF_z]^T

式中gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间的安装误差角；aSF_x、aSF_y、aSF_z分别为三轴加速度计标度因数误差；aB_x、aB_y、aB_z分别为三轴加速度计零偏；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角；

步骤2：求解观测矩阵A；

观测矩阵 表示T₂时刻速度误差一阶导数矢量，表示T₁时刻速度误差一阶导数矢量；

观测矩阵由三向速度误差一阶导数组成，光纤捷联惯导系统静止时，速度误差即为光纤捷联惯导系统导航输出的速度值；采用标准卡尔曼滤波计算速度误差一阶导数值；

建立状态方程：

建立量测方程：

其中，Z_i表示观测矢量；δV_i(t)＝V_i(t),i＝x,y,z，δV_i(t)表示速度误差矢量，V_i(t)表示导航解算速度矢量，w_i和μ_i分别表示状态噪声和观测噪声；采用标准卡尔曼滤波方程进行迭代，滤波完毕即得到速度误差一阶导数的估计值求取转动前后两次速度误差一阶导数并作差，即得到系统观测阵A；

步骤3：求解转动系数矩阵K；

转动系数矩阵

式中，df＝[df_x df_y df_z]^T表示加速度计误差系数阵，dφ＝[dφ_x dφ_y dφ_z]^T表示陀螺误差系数阵；

其中，分别表示T₁、T₂时刻载体系b系到导航系n系的方向余弦矩阵，

其中i＝1,2；

$d_{\mathrm{\φ}}=-{\∫}_{T_1}^{T_2}{C}_{b\left(T_1\right)}^n\·C_{b\left(t\right)}^{b\left(T_1\right)}\·\left[\begin{array}{ccc}{\mathrm{\ω}}_x\left(t\right)& 0& 0\\ 0& {\mathrm{\ω}}_y\left(t\right)& 0\\ 0& 0& {\mathrm{\ω}}_z\left(t\right)\end{array}\right]M_g\·dt$

其中，f_ibx ^b(T_i)、f_iby ^b(T_i)、f_ibz ^b(T_i)分别表示T_i时刻三轴加速度计测量值；为惯组初始时刻T₁的姿态矩阵，ω_x(t)、ω_y(t)、ω_z(t)为t时刻的转动角速度，而为t时刻b系到T₁时刻b系的转换矩阵，M_g为系数矩阵，这两者与转动轴向有关；

若光纤捷联惯导系统绕X轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\left({\mathrm{\ω}}_x\right(t\left)\right)& sin\left({\mathrm{\ω}}_x\right(t\left)\right)\\ 0& -sin\left({\mathrm{\ω}}_x\right(t\left)\right)& \cos \left({\mathrm{\ω}}_x\right(t\left)\right)\end{array}\right],M_g=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right]$

若光纤捷联惯导系统绕Y轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}|=\left[\begin{array}{ccc}cos\left({\mathrm{\ω}}_y\right(t\left)\right)& 0& -sin\left({\mathrm{\ω}}_y\right(t\left)\right)\\ 0& 1& 0\\ sin\left({\mathrm{\ω}}_y\right(t\left)\right)& 0& \cos \left({\mathrm{\ω}}_y\right(t\left)\right)\end{array}\right],M_g=\left[\begin{array}{ccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\end{array}\right]$

若光纤捷联惯导系统绕Z轴转动，则：

$C_{b\left(t\right)}^{b\left(T_1\right)}|=\left[\begin{array}{ccc}cos\left({\mathrm{\ω}}_z\right(t\left)\right)& sin\left({\mathrm{\ω}}_z\right(t\left)\right)& 0\\ -sin\left({\mathrm{\ω}}_z\right(t\left)\right)& cos\left({\mathrm{\ω}}_z\right(t\left)\right)& 0\\ 0& 0& 1\end{array}\right],M_g=\left[\begin{array}{ccccccccc}0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$

步骤4：采用最小二乘方法求解误差向量：X＝K^T·(KK^T)^-1·A；由于转动系数矩阵K的秩rank(K)＝15，得到的误差系数有15项，包括三轴陀螺仪标度因数误差gSF_x、gSF_y、gSF_z，各轴陀螺仪间安装误差角gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy，三轴加速度计零偏aB_x、aB_y、aB_z；各轴加速度计间安装误差角aMA_yx、aMA_zx、aMA_zy；

步骤5：求解三轴加速度计标度因数；

采用第1次转动数据：

采用第4次转动数据：

采用第7次转动数据：

式中，表示转动前从0时刻到T₁时刻的时间内天向速度误差一阶导数之和，表示转动后从0时刻到T₂时刻的时间内天向速度误差一阶导数之和；

第六步：根据标定得到18项误差系数，对保存的9组惯性器件数据进行补偿，得到新的9组惯性器件数据；

第七步：对于第六步中得到的新的9组惯性器件数据，重复第五步和第六步3～5次，即迭代器件18项误差系数值进行重复最小二乘拟合计算，直至收敛；累加每次迭代估计得到器件18项误差系数值，即为最终标定参数值；

第八步：根据上述标定的18项误差系数，对光纤捷联惯导系统进行补偿；

第九步：重复第一步到第二步，按照第三步中的序列1进行手动翻转正六面体装置，转动前后分别静止时间30min，保存陀螺仪输出数据，基于该输出数据，计算三轴陀螺零偏误差；

第十步：根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿，采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可，完成了光纤捷联惯导系统21项误差参数的现场标定。

2.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，其特征在于：第三步中所述转动角允许存在±10°误差。

3.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，其特征在于：第四步中所述的对惯性器件误差进行建模，包括光纤陀螺误差模型和加速度计误差模型，分别如下：

${{\mathrm{\δ\ω}}_{ib}}^b=\left[\begin{array}{c}{{\mathrm{\δ\ω}}_{ibx}}^b\\ {{\mathrm{\δ\ω}}_{iby}}^b\\ {{\mathrm{\δ\ω}}_{ibz}}^b\end{array}\right]=\left[\begin{array}{c}{\mathrm{gB}}_x\\ {\mathrm{gB}}_y\\ {\mathrm{gB}}_z\end{array}\right]+\left[\begin{array}{ccc}{\mathrm{gSF}}_x& {\mathrm{gMA}}_{xy}& {\mathrm{gMA}}_{xz}\\ {\mathrm{gMA}}_{yx}& {\mathrm{gSF}}_y& {\mathrm{gMA}}_{yz}\\ {\mathrm{gMA}}_{zx}& {\mathrm{gMA}}_{zy}& {\mathrm{gSF}}_z\end{array}\right]\left[\begin{array}{c}{{\mathrm{\ω}}_{ibx}}^b\\ {{\mathrm{\ω}}_{iby}}^b\\ {{\mathrm{\ω}}_{ibz}}^b\end{array}\right]$

${\mathrm{\δf}}_{ib}^b=\left[\begin{array}{c}{\mathrm{\δf}}_{ibx}^b\\ {\mathrm{\δf}}_{iby}^b\\ {\mathrm{\δf}}_{ibz}^b\end{array}\right]=\left[\begin{array}{c}{\mathrm{aB}}_x\\ {\mathrm{aB}}_y\\ {\mathrm{aB}}_z\end{array}\right]+\left[\begin{array}{ccc}{\mathrm{aSF}}_x& 0& 0\\ {\mathrm{aMA}}_{yx}& {\mathrm{aSF}}_y& 0\\ {\mathrm{aMA}}_{zx}& {\mathrm{aMA}}_{zy}& {\mathrm{aSF}}_z\end{array}\right]\left[\begin{array}{c}{f}_{ibx}^b\\ f_{iby}^b\\ f_{ibz}^b\end{array}\right]$

式中δω_ib ^b表示陀螺仪的误差输出矢量；δω_ibx ^b、δω_iby ^b、δω_ibz ^b分别表示由三轴陀螺误差引起的误差角速度；ω_ibx ^b、ω_iby ^b、ω_ibz ^b分别表示三轴陀螺原测量值；gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间安装误差角；gB_x、gB_y、gB_z分别表示三轴陀螺仪零偏误差；

式中表示加速度计的误差输出；δf_ibx ^b、δf_iby ^b、δf_ibz ^b分别表示由三轴加速度计误差引起的误差加速度；f_ibx ^b、f_iby ^b、f_ibz ^b分别表示三轴加速度计原测量值；aSF_x、aSF_y、aSF_z分别为三轴加速度计的标度因数误差；aB_x、aB_y、aB_z分别为三轴加速度计的零偏；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角。

4.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，其特征在于：第六步中所述的对保存的9组惯性器件数据进行补偿，得到新的9组惯性器件数据，其中补偿公式如下：

$\left[\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right]=\left(I+\left[\begin{array}{ccc}{\mathrm{gSF}}_x& {\mathrm{gMA}}_{xy}& {\mathrm{gMA}}_{xz}\\ {\mathrm{gMA}}_{yx}& {\mathrm{gSF}}_y& {\mathrm{gMA}}_{yz}\\ {\mathrm{gMA}}_{zx}& {\mathrm{gMA}}_{zy}& {\mathrm{gSF}}_z\end{array}\right]\right)\left(\begin{array}{c}{{\mathrm{\ω}}_{ibx}}^b\\ {{\mathrm{\ω}}_{iby}}^b\\ {{\mathrm{\ω}}_{ibz}}^b\end{array}\right)$

$\left[\begin{array}{c}{f}_x\\ f_y\\ f_z\end{array}\right]=\left(I+\left[\begin{array}{ccc}{\mathrm{aSF}}_x& 0& 0\\ {\mathrm{aMA}}_{yx}& {\mathrm{aSF}}_y& 0\\ {\mathrm{aMA}}_{zx}& {\mathrm{aMA}}_{zy}& {\mathrm{aSF}}_z\end{array}\right]\right)\left[\begin{array}{c}{f}_{ibx}^b-{\mathrm{aB}}_x\\ f_{iby}^b-{\mathrm{aB}}_y\\ f_{ibz}^b-{\mathrm{aB}}_z\end{array}\right]$

式中gSF_x、gSF_y、gSF_z分别表示三轴陀螺仪标度因数误差；gMA_xy、gMA_xz、gMA_yx、gMA_yz、gMA_zx、gMA_zy分别表示各轴陀螺仪间安装误差角；aSF_x、aSF_y、aSF_z分别为三轴加速度计标度因数误差；aB_x、aB_y、aB_z分别为三轴加速度计零偏；aMA_yx、aMA_zx、aMA_zy分别表示各轴加速度计间安装误差角；ω_ibx ^b、ω_iby ^b、ω_ibz ^b分别表示三轴陀螺原测量值；f_ibx ^b、f_iby ^b、f_ibz ^b分别表示三轴加速度计原测量值；I表示三阶单位阵；w_x、w_y、w_z分别表示三轴陀螺补偿后测量值；f_x、f_y、f_z分别表示三轴加速度计补偿后测量值。

5.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法，其特征在于：计算陀螺零偏误差计算公式如下：

$gB=\frac{1}{2}\{\left({\mathrm{\ω}}_{ib}^b\right(T_1)+{\mathrm{\ω}}_{ib}^b(T_2\left)\right)-{\[C_b^n\left(T_1\right)+C_b^n\left(T_1\right)\]}^T(I-\mathrm{\φ}\×){\mathrm{\ω}}_{ie}^n\}$

其中，φ＝[0 0 φ_z]^T，φ_z表示航向角误差，

$\mathrm{\Δ}{\mathrm{\ω}}_{\mathrm{ib}}^b={\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_2\right)-{\mathrm{\ω}}_{\mathrm{ib}}^b\left(T_1\right),\mathrm{\Δ}{C}_b^n=C_b^n\left(T_2\right)-C_b^n\left(T_1\right);$

式中分别表示T₁、T₂时刻三轴陀螺测量矢量，表示转动前后两时刻陀螺测量误差矢量、表示第一行分量，分别表示T₁、T₂时刻载体系b系到导航系n系的方向余弦矩阵，表示转动前后两时刻方向余弦矩阵误差，表示第i行第j列元素值，表示导航系下地球速度矢量，分别表示导航系下地球速度矢量的北向分量和天向分量。