一、概念及意义

拿到一个IMU，首要任务是对其器件误差进行分析，包括陀螺仪的误差和加速度计的误差。由于陀螺仪与加速度计的误差特性类似，仅是测量的物理量含义不同，因此我们就以陀螺仪为准来介绍随机误差成分、原理及分析方法。加速度计的误差分析直接照葫芦画瓢就行。

陀螺仪的误差主要包括这样几项：

1. 量化噪声

量化噪声是数字传感器必然出现的噪声，我们通过AD采集把连续时间信号采集成离散信号，在这个过程中，精度就会损失，损失的精度大小和AD采样的精度有关（这里具体指的是模数转换时，AD器件的位数，位数越高采样越精确），精度越高，量化噪声越小。

2. 角度随机游走

陀螺敏感角速率并输出时是有噪声的，这个噪声里面的白噪声成分叫宽带角速率白噪声，我们计算姿态时，本质上是对角速率做积分，这必然会对噪声也做了积分。白噪声的积分并不是白噪声，而是一个马尔可夫过程，即这一次的误差是在上一次误差的基础上累加一个随机白噪声得到的。角度误差所包含的这种马尔可夫性质的误差就叫做角度随机游走。

3. 角速率随机游走

从理解上和角度随机游走一样，角速率里面并不全是白噪声，它也有马尔可夫性质的误差成分，而这个误差是由宽带角加速率白噪声累积的结果。

4. 零偏不稳定性噪声

这应该是大家再熟悉不过的一个误差项了，如果一个陀螺只让你用一个指标来体现精度，那必然就是它了。但是这个指标的理解上却不像前几个参数那样直白。

我们可以先把它理解为零偏随时间的缓慢变化，假设在刚开始时零偏大小是某个值，那么过一段时间之后，零偏便发生了变化，具体变化成了多少，无法预估，所以就要给他一个概率区间，来描述它有多大的可能性落在这个区间内，时间越长，区间越大。

实际上，如果你真的测的时间足够长，会发现它也不会无限制增长下去，所以，这个对概率区间的描述只是近似有效，或者一定时间内有效，由于这个有效时间比较长，所以我们一般仍然使用这种方式来描述，只是在理解上要知道这一点的存在。

5. 速率斜坡

看到斜坡这种描述词，我们一般会想它是不是一种趋势项。实际上，它确实是趋势性误差，而不是随机误差。所谓随机误差，是指你无法用确定性模型去拟合并消除它，最多只能用概率模型去描述它，这样得到的预测结果也是概率性质的。而趋势性误差是可以直接拟合消除的，在陀螺里，这种误差最常见的原因是温度引起零位变化，可以通过温补来消除。

6. 零偏重复性

这项误差很好理解，IMU断电之后重新上电时，它当前的bias并不与断电之前的bias一致，即使断电与上电两个动作之间的间隔时间很短。这个和器件本身的结构特性相关，我们不必深入到这个层面去追究它的原因，理解这个现象即可。

不同原理的器件，零偏不稳定性差异非常大，对于激光陀螺仪这种，它的零偏重复性要远小于其零偏不稳定性，而对于MEMS这种器件，零偏重复性甚至会比零偏不稳定性大两个数量级，因此MEMS每次上电之后，就要去重新估计它的bias。

二、常用误差分析方法

1. Allan方差

这应该是当前最常用的误差分析方法了，我们得到的数据都是各种误差混合在一起的，这种方法能够把各种混合误差很好地分离出来，使之不会互相影响。比如，有的陀螺温度特性比较好，有的比较差，用别的方法可能分析出的误差会因此产生很大的区别，而用Allan方差，则能够把温变误差隔离开，并准确地识别零偏不稳定性等随机误差。

需要注意的是，该方法可以分析上面介绍的前5项误差，而不能分析第6项，即零偏不稳定性。其实该项误差很好统计，找一个器件，反复上电、断电，计算输出的变化就可以了，不必使用很复杂的误差分析方法。

2. 直接方差统计

具体原理就是，把数据按一定时间区间分成多段，常见的分段区间长度是100s一段，每段计算出一个均值，然后所有段的均值计算出一个方差，就作为陀螺的精度指标了。

这种方法得到的统计值，在理论上和上面介绍的6项误差中的任何一项都不对等，它是在Allan方差方法出现之前所使用的，而且之所以出现Allan方差，就是因为这种方法统计的结果不准，才被迫另寻高招。

按照上面介绍的各种随机误差原理，我们能看到，这种方法统计出来的结果，必然是统计时间越长，指标越差。而且冷启动和热启动结果也不一样，因为有温度漂移误差在里面。所以这种方法本应被淘汰的，之所以还有些在用，多数是以前的流程里沿用下来的，变成历史遗留问题了。

由于不同方法统计出的结果不一样，所以我们在比较器件精度的时候，一定要确认好他们是否是同一种方法统计出来的。另外需要注意的是，在直接方差统计方法中，区间分段时间（也就是上面说的100s）不同时，统计的结果也会不同，区间段越长，精度越好。所以不仅要确认方法相同，还要确认参数相同，有的告诉别人说自己是直接方差统计方法，却不说自己使用的是1000s分段，这就是典型的耍流氓了。

三、Allan方差原理介绍

我们都知道，计算Allan方差，是通过拟合方差随时间间隔变化的双对数曲线来实现的。所以，我们所计算的Allan方差，其实是一个系数，它体现了方差和时间间隔之间的关系。

上面我们介绍了五种误差，所以就存在五种关系，这里我们直接给出结果公式，推导过程可参考西工大严恭敏老师的《惯性仪器测试与数据分析》。

以上五种误差，和时间间隔在双对数坐标轴上分别对应一条直线，斜率依次分别是-1，-1/2，1/2，0，1，如果我们把这五条曲线画在一张图上，会看到下面的图形。（图太丑，凑合看吧）

五种误差对应了五条虚线，实际数据是五条虚线的叠加。而我们用实际数据测出来的方差曲线往往对应的都是图中实线的形状，之所以会这样是因为这是双对数坐标，在这样的坐标系下，纵轴上半部分的值是比下半部分的值大很多的，所以累加之后就只剩下上半部分对应的形状了。

需要注意的是，实际数据分析曲线中，不一定能完整体现出这全部的五种误差，如果有某一项比较小，就可能被其他误差所淹没。反之，如果某一项特别大，则可能把其他几项全部淹没。理想形状出现的概率并不是很高，不必过于纠结。

四、Allan方差代码实现

1. 实现方式

Allan方差识别误差的方式是画“方差-时间间隔”双对数曲线，这包含三方面：方差、时间间隔、双对数(好像是一句废话，哈哈)。

在测试中，我们按照一定的周期T采集了一段时间数据，根据这些离散的数据，便能计算其方差。这个周期T和方差在二维坐标系中便构成了一个点的横纵坐标。为了得到曲线，我们需要得到更多这样的点，把周期缩短是不可能了，因为不能凭空提高已有数据频率，所以只能把周期拉长。

如果我们对每两个相邻的数取平均，合并成一个，便构成了周期为2T的离散数据，同样也可以得到它的方差，也即又多得到了一个点。按照这样的方法，循环下去，可以得到周期为4T、8T、16T对应的点，直到数据少到无法再合并为止，这样就得到系列的点，也就是曲线了。

这样做出来的曲线，其横纵坐标都会很长，各个点之间的间隔也极其不均匀，这也就是为什么用双对数曲线的原因了。

2. 参数拟合

我们先给出一张Allan方差曲线图，对着图来说吧。

在理论上，图中的曲线是由五个直线组成的，五个直线的参数分别包含了我们想要提取的五项误差。所以，只要把直线公式拟合出来，那么参数自然就知道了。从公式上可以看出，我们只要找到直线在任一横坐标处的值，便可提取参数，方便起见，一般取横坐标为1处的值，因为这样，公式右边和时间相关的项就是零。

3. 代码实践

我们使用开源系统imu_utils为例，来看提取参数的具体实现方式。

开源代码地址：github.com/gaowenliang/

1) 代码思路

以上我们拟合直线，再求交点的方式，适合人工操作，按这种方式，给出一个曲线图，很快就可以计算出各误差参数。而对于机器来讲，则似乎不太适合，首先让它把曲线自适应地分成五段就要费一番功夫，此时就需要换一种思路。

既然曲线是五条直线组成，而五条直线的公式已知，那么五条直线加一起不就是曲线的公式了吗，有了曲线公式，又有了曲线实际数据，那么这就变成了一个曲线拟合问题，用优化的方式就可以直接解决。

比较来讲，这两种方式无好坏之分，只能说一个更适合人工，一个更适合机器，为各自量身打造。

2) 代码细节

首先是把原始数据划分成不同的时间间隔，思路就是上面提到的倍增法。

int mode = numData / 2;
unsigned int maxStride = 1;
int shft  =0;
while(mode){
    mode = mode >> 1;
    maxStride = 1 << shft;
    shft++;
}

然后每一个间隔下都计算一个方差，并把方差和真实曲线建立残差模型。

for ( int i = 0; i < num_samples; ++i ) {
    ceres::CostFunction* f = new ceres::AutoDiffCostFunction< AllanSigmaError, 1, 5 > (new AllanSigmaError( sigma2s_tmp[i], m_taus[i] ) );
    problem.AddResidualBlock( f, NULL, param );
}

其中残差的计算如下:

// 计算残差
template< typename T >
bool operator( )( const T* const _paramt, T* residuals ) const {
    T _Q   = T( _paramt[0] );
    T _N   = T( _paramt[1] );
    T _B   = T( _paramt[2] );
    T _K   = T( _paramt[3] );
    T _R   = T( _paramt[4] );
    T _tau = T( tau );

    T _sigma2    = calcSigma2( _Q, _N, _B, _K, _R, _tau );
    T _dsigma2   = T( calcLog10( _sigma2 ) ) - T( calcLog10( sigma2 ) );
    residuals[0] = _dsigma2;
        
    return true;
}

// 叠加五个直线，得到曲线方程
template< typename T >
T calcSigma2( T _Q, T _N, T _B, T _K, T _R, T _tau ) const {
    return  _Q * _Q / ( _tau * _tau )
          + _N * _N / _tau
          + _B * _B
          + _K * _K * _tau
          + _R * _R * _tau * _tau;
}

最后，通过迭代优化，就可以直接计算出那五个参数了。

五、一些实践经验

至此为止，我们应该已经弄清了器件误差提取方法，然而，实践总是会给人惊喜，出现一些意料之外，而且我们也不必完全拘泥于理论，要看我们的目的是什么。

1) 有的时候不一定五个斜率全都能体现出来，我们知道，每个斜率对应一项误差，如果某项误差相对于其他项的量级比较小，那么曲线中就没有他什么事了，直接被掩盖，极端情况下，一个曲线只有一个斜率也是有可能的。

2) Allan方差对误差项提取的比较多，但是很多误差项的分析都是为了提高器件生产工艺用的，每种误差成分产生的原理不同，通过分析哪种误差成分为主，我们就知道为了提高精度，应该重点改进哪项工艺。如果只是为传感器融合找个方差的参数，大可不必计算那么细致，在图上看一下大致量级，就可以直接用了，因为参数早晚都是要调的。