三、Allan方差分析

1. 概述

在上一篇文章，我们对器件误差的成分进行了分析，并且介绍了最主要的五种误差成分在“方差-时间间隔”的双对数曲线上的

直观表现及产生这种表现的原理。同时也提到，在两种误差分析方法中，Allan方差方法比直接方差统计要更有优势，

所以本篇文章我们就一起讨论下Allan方差的具体实现方式。

2. 实现方式

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

在测试中，我们按照一定的周期T采集了一段时间数据，根据这些离散的数据，便能计算其方差。这个周期T和方差在二维坐标系中便构成了一个点的横纵坐标。

为了得到曲线，我们需要得到更多这样的点，把周期缩短是不可能了，因为不能凭空提高已有数据频率，所以只能把周期拉长。

如果我们对每两个相邻的数取平均，合并成一个，便构成了周期为2T的离散数据，同样也可以得到它的方差，也即又多得到了一个点。按照这样的方法，循环下去，

可以得到周期为4T、8T、16T对应的点，直到数据少到无法再合并为止，这样就得到系列的点，也就是曲线了。

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

3. 参数拟合

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

在理论上，图中的曲线是由五个直线组成的，五个直线的参数分别包含了我们想要提取的五项误差。

所以，只要把直线公式拟合出来，那么参数自然就知道了。我们把上一篇文章的公式再放过来看一下。

从公式上可以看出，我们只要找到直线在任一横坐标处的值，便可提取参数，方便起见，一般取横坐标为1处的值，因为这样，公式右边和时间相关的项就是零。

4. 代码实践

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

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

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;
}

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

 

5. 一些实践经验

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

1)  有的时候不一定五个斜率全都能体现出来，我们知道，每个斜率对应一项误差，如果某项误差相对于其他项的量级比较小，那么曲线中就没有他什么事了，

直接被掩盖，极端情况下，一个曲线只有一个斜率也是有可能的。

2)  Allan方差对误差项提取的比较多，但是很多误差项的分析都是为了提高器件生产工艺用的，每种误差成分产生的原理不同，通过分析哪种误差成分为主，

我们就知道为了提高精度，应该重点改进哪项工艺。如果只是为传感器融合找个方差的参数，大可不必计算那么细致，在图上看一下大致量级，

就可以直接用了，因为参数早晚都是要调的。

严格来讲，滤波的方差其实应该用白噪声强度，它不属于这五项误差的任何一项，但是很多时候手册上不给出这个强度值，

而方差参数在用的时候都会结合系统再一次进行调整，所以有些时候直接用零偏不稳定性做方差初值也不是不可以。

3)  随机游走这种东西理论上是可以建模，但是实际滤波模型中却很少使用，甚至不建议使用。一是、因为它并不是构成零偏趋势项的主要成分，

往往都被温变等的趋势掩盖，即使做了温补也是这样。二是、因为Allan方差识别出的参数并不完全准确，这会带来模型误差，

不准确的模型还不如没有模型，因为它会使滤波变得不稳定，效果反而变差。

 

四、IMU误差标定之基于转台的标定