组合导航原理剖析（六）：光纤级别惯性测量单元（IMU）的误差测试与标定

内容主要参考任乾在知乎的博客《多传感器融合定位理论基础（四）：IMU内参模型及标定》，以及严恭敏的《惯性仪器测试与数据分析》、经典论文《A Robust and Easy to implement method for imu calibration without External Equipments》。本文仅针对于精度指标为0.5°/h光纤陀螺仪。

在获得厂家给到的高精度惯组之后，常需要对其进行测试性能和标定的工作，在无转台的条件下，想要准确的进行标定并不容易。

IMU的噪声来源多样，对其进行误差分析可以使用Allan方差的方式，在静止状态下对IMU进行采样，在理想情况下可以标定量化噪声、随机游走噪声、零偏不稳定性等5项主要参数。实际测试下发现往往并不能对多项参数进行估计，认为大概看一下随机游走噪声的幅值范围，对后续导航解算有所参考即可，不必拘泥于相关的理论。本人认为Allan方差的最大功能就是评估一下IMU在静态下的各项随机噪声的特性，也算是一种测试手段。

而对于IMU内参的标定，不管是高精度的光学陀螺还是低精度的MEMS，通常指的是零偏、标度因数和交轴耦合偏差（安装偏差）。目前主流的思路主要分为基于“非线性优化”离线处理的思路和基于“卡尔曼滤波”在线估计的思路2大类，前者简单的处理方式也可以采用“12位置标定法”，而后者见于严恭敏的书中。

《A Robust and Easy to implement method for imu calibration without External Equipments》中就采用非线性优化的方式对参数进行标定，思路比较简单，论文针对MEMS陀螺和加表进行了全参数的标定，本人并没有消费级IMU用于测试，因此不讨论其标定精度。

在严的书中，给出了一种借助于转台的标定方案。在与光纤陀螺厂家的相关研发人员沟通过程中，发现各厂家的标定方案都不太一样，详细标定方法不得而知，但往往生产环节的标定方案很简单直接。按照这一思路，认为严给出的经典标定方法是一种可用于生产环节的标定方案。详细内容见书本，大概可以概括为：使用12位置标定法对加速度计进行标定，然后使用速率转台对陀螺仪的标度因数进行标定，最后使用两个角度静置陀螺30min的方式标定陀螺的零偏。可以看出这种方式需要借助于6自由度转台，实际测试条件下往往缺少高精度6自由度转台进行标定，也可以有替代的标定方案。

首先介绍一下使用最小二乘法和十二位置标定法标定加速度计，详细见严的书，这里根据书中的数据给出实现方法，如下。

data = [1	-32.604	21.762	9990.421
2	-32.616	21.272	9990.464
3	-25.436	0.502	-10031.447
4	-25.431	1.004	-10031.454
5	-33.86	10017.538	-28.156
6	-34.354	10017.509	-28.151
7	-23.175	-9994.914	-12.691
8	-22.723	-9994.933	-12.685
9	9969.802	14.938	-19.102
10	9969.787	14.96	-19.614
11	-10027.294	6.902	-21.903
12	-10027.271	6.892	-21.411];

g0 = [0,0,1
    0,0,1;
    0,0,-1
    0,0,-1
    0,1,0
    0,1,0
    0,-1,0
    0,-1,0
    1,0,0
    1,0,0
    -1,0,0
    -1,0,0];
yk = reshape(g0', 12*3, []);
H = {};
for i = 1:12
    H{i} = zeros(3, 12);
    Nak = [data(i, 2:4), -1];
    H{i}(1, 1:4) = Nak;
    H{i}(2, 5:8) = Nak;
    H{i}(3, 9:12) = Nak;
end
Hk = [];
for i = 1:12
    Hk = [Hk; H{i}];
end
xk = inv(Hk'*Hk)*Hk'*yk;
bias = xk([4,8,12])*1e3;
Ka = [xk(1:3)';xk(5:7)';xk(9:11)'];
Ka_1 = inv(Ka);
K1_x = norm(Ka_1(1, :));
K1_y = norm(Ka_1(2, :));
K1_z = norm(Ka_1(3, :));

严的书中提到：在加速度计的标定过程中，当地的重力加速度矢量准确已知，可作为精确的比力激励基准使用。虽然地球自转角速度${\omega }_{ie}$是准确已知的，但是该角速度量级太小，不宜当做标定陀螺仪的角速度激励源来使用，必须借助于转台的转动提供更大的加速度激励。加表的标定如无需标定安装参数，完全可以对单轴垂向放置进行标定，只需要对标定台调平即可。忽略在无转台辅助情况下难以标定的多轴交轴耦合偏差，以z轴为例，加表标定误差的推导可以表示为（式中1为正放，即正轴方向朝向天）

$$\begin{array}{l}{a}_m=K_a{a}^b+b_a\\ a_{m,z1}=K_{a,z}\left[\begin{array}{l}0\\ 0\\ g_0\end{array}\right]+b_{a,z}\\ a_{m,z2}=K_{a,z}\left[\begin{array}{l}0\\ 0\\ -g_0\end{array}\right]+b_{a,z}\end{array}$$

式中的g0需要根据当地经纬高进行严格计算。

调平偏差对非垂直轴向的加表标定并不友好，因此采用单轴分解标定的方式是比较合理的。用下面的代码评估精度，发现和12面标定法精度差不多，仅零偏存在4ug左右的偏差。

% 通过正反置的方式进行测试惯组
up = (9990.421 + 9990.464 )/2;
down = (-10031.447  -10031.454)/2;
k2_z = (up - down)/2;
ba_z = (up + down)/k2_z/2 * 1e3;

当陀螺仪放置水平时，可以同样通过正反放置陀螺仪的方式，利用摇头陀螺敏感地球自转不变形的原理，去测量陀螺仪的零偏。但是由于地球自转的量过小，并不能辅助识别陀螺仪的标度因数。此外，还可以用当地地球自转来来校验标定结果的准确性。以z轴为例，陀螺仪标定过程的推导可以表示为

$$\begin{array}{l}{\omega }_{ib}^b={\omega }_{nb}^b+C_n^b{\omega }_{in}^n\\ {\omega }_{ib.m}^b={\omega }_{nb}^b+C_n^b{\omega }_{in}^n+b_g\\ {\omega }_{ib.m,z1}^b={\omega }_{ie}\sin L+b_{g,z}\\ {\omega }_{ib.m,z2}^b=-{\omega }_{ie}\sin L+b_{g,z}\end{array}$$

通过以上过程，就可以简易测试光纤陀螺仪与加速度计，一组结果表示如下。

天向加速度计 bias is 791.235 ug, scale is 186.287 ppm
摇头陀螺仪 bias is -0.137 deg/h, err is 0.042 deg/h
右向加速度计 bias is -724.513 ug, scale is -113.828 ppm
俯仰陀螺仪 bias is 0.227 deg/h, err is -0.029 deg/h
前向加速度计 bias is -383.168 ug, scale is -79.433 ppm
测滚陀螺仪 bias is -0.214 deg/h, err is -0.091 deg/h

其中提到的陀螺仪的err，实际上值得是通过正反向标定的方式，可以测得当地的地球自转的分量，表示为${\omega }_{ie}sinL$，与理论值之间的偏差，该项可以反应标定的误差是否在可接受的范围。比如上面的误差只有0.09deg/h，对于当前标定的陀螺而言就是可以接受的。

最后，再来讨论一下IMU的技术指标参数，仅讨论指标对应的测试方法。标定过程通常情况下虽然能够测量准确的零偏值，但是零偏值会一直变化，长期来看基本围绕着一个均值在发生变化。因此，通常统计IMU的零偏稳定性变化。通常用10s或者100s平滑方法计算零偏稳定性，下面是静置16个小时下，100s平滑得到的零偏值变化，扣除了包含地球自转、重力矢量等带来的常值偏差，图中曲线的波动范围即为常说的零偏稳定性指标。图中加速度计的零偏稳定性为200ug，陀螺仪的零偏稳定性为0.5°/h。

另一个指标是上电零偏重复性，测试方式是逐次上电统计零偏值，观察其噪声范围。

另一个指标是常值零偏，这个可以用过上述标定得到。