惯性测量单元Allan方差分析详解

• ⭐ Zeal's Blog
• 🛠 Zeal's 知乎专栏
• 🌀 Zeal's Github

前言

原本以为Allan方差应该是一个简单的方差计算方法，但是经过两天的各种探索和资料学习，发现这玩意的原理远比想象中的复杂。详细的公式推导需要设计到随机信号的频谱分析、IMU各种误差模型等内容。并且，Allan方差一开始用来分析原子钟晶体振荡频率的误差分析，后续这种分析方法被用到IMU的误差分析，并且IEEE也制订了相应标准（见参考资料1），这一系列的发展过程其实很长，同时也有不同的定义。以上种种，让读者放弃了摸清Allan方差的详细来龙去脉和公式推导，毕竟，工程师和数学家还是有区别的……但是关于Allan方差的解释和应用网上的博文也都不够全面，单看某一篇文章很难获取对Allan方差的直观理解（没有找到一篇很好的文章能够说明Allan方差的定义和双对数曲线不同斜率段对应的噪声项是怎么来的）。

本文对以下内容做了梳理：

1. IMU的噪声项及其含义
2. Allan方差的定义和推导
3. 对陀螺仪的数据计算Allan方差的方法
4. Allan方差双对数曲线的含义和分析
5. 关键概念解读

PS: 笔者水平有限，本文只是笔者对Allan方差学习的总结，并不是一篇严格的理论推导文章

惯性测量单元（IMU：Inertial Measuremnt Unit）误差项

这部分可以参考参考资料2

IMU模型

IMU的输出可以看上图，典型的六轴IMU输出三个轴的加速度和三个轴的角速度。对加速度做积分可以得到速度、再做一次积分可以得到位移；对角速度积分可以得到角度。

电子器件输出都会有误差，又叫噪声。IMU的噪声项主要有以下几种：

IMU误差示意图

1. 零偏（Bias），单位deg/hr,m/s^2 零偏这个词似乎不太准确，因为在非零输出时也会有偏差。如果把IMU的输出当作一个遵循高斯分布的随机信号，在静止时该随机信号的均值应该是零，但是由于bias，该均值并不是0。见上面的IMU误差示意图可以较好地理解。要注意的是零偏是会随着时间改变的，因此在很多卡尔曼滤波模型中会把Bias当作一个状态量。
2. 零偏可重复性（Bias Repeatability），单位deg/hr,m/s^2 每次上电零偏的稳定程度。
3. 零偏稳定性（Bias Stability），单位deg/hr/hr,m/s^2/hr 有时也叫零偏不稳定性（Bias Instability）。描述零偏的稳定程度，也可以理解为零偏的噪声。（有些文章也会叫做零偏的角度随机游走）。
4. 比例因子（Scale Factor） 比如比例因子是1.2,那么当真值是10m/s^{2的时候，IMU的输出可能是12m/s}2。
5. 角度随机游走（ARW:Angle Random Walk），单位：$deg/\sqrt{hr}$,$m/s/\sqrt{hr}$ 这是一个比较难理解的量，实际上笔者也没法完全理解。简单理解就是当信号叠加上一个白噪声之后，对该信号的积分就包含了白噪声的积分，白噪声的积分不再是一个白噪声，而是一个马尔科夫模型，也就是此时的噪声是上一时刻的噪声加上一个白噪声。那么放在这里，对于陀螺仪的输出（角速度）一般是用来做积分得到角度，那么角速度的白噪声的积分就会得到角度的偏差，因此叫做角度随机游走。有时候角度随机游走的值又叫做白噪声强度，比如在kalibr:IMU-Noise-Model中就有这样的描述： > White Noise Terms: The parameters for the "white noise" processes are often specified in the datasheet of the sensor manufacturer. A bit misleading, they are commonly denoted as angular random walk in case of the gyro, and velocity random walk for the accel
6. 速率随机游走（RRW:Rate Random Walk） 与角度随机游走类似，速率随机游走可以看作是角加速度的误差积分导致的角速度的误差。速率随机游走对应的是陀螺仪在长时间下才会发生改变的误差项，因此可以用来近似当作Bias使用。在kalibr:IMU-Noise-Model中使用速率随机游走系数当作零偏参数。

角度随机游走（ARW:Angle Random Walk）v.s. 角速率随机游走（ARRW）

这个很容易弄混淆的一个概念。对应到加速度计上应该叫速度随机游走和加速度随机游走。因为讲解Allan方差的时候一般会以陀螺仪为示例，因此在文中可能只出现角度随机游走和角速率随机游走。角速率随机游走又会被简写成速率随机游走，比如参考资料1的IEEE标准。下面这段百度百科的描述反而清晰： > 角速率传感器的性能参数中包含角度随机游走（ARW）和角速率随机游走（ARRW），前者反映了角速率信号中白噪声的特性，积分后表现为角度随机游走；后者为角速率本身的随机游走，可认为是角加速度白噪声的积分结果。

Allan方差

Allan方差又叫阿伦方差，可以用Allan方差用来描述物理器件的不同噪声参数。Allan方差相比于普通的方差具有更好的描述长时间误差的优势（具体也不太清楚）。Allan方差现在好像基本被用来当作IMU噪声标定的标准方法了。

数据方差

数据方差及标准差计算方法

（上图引用自新手入门系列3——Allan方差分析方法的直观理解）

系数$\frac{1}{N-1}$使得这样计算出来的方差是数据方差的无偏估计（概率统计原理）。

对不同时间尺度上的误差刻画

上面的数据方差计算方法对数据的整个序列计算误差水平，无法细分出不同时间尺度上的误差波动情况。为了能够刻画不同时间周期下数据的波动情况，我们可以先对数据分块，对每块数据取平均作为该块数据的代表。假设原始数据的采样时间是${\tau }_0$，每3个数据分成一块，则每块的间隔时间$\tau =3{\tau }_0$。分块之后，我们再将$M$个数据块计算方差，最后求这个方差的数学期望，作为我们最后的方差

通用样本方差计算方法

$$\begin{array}{}\text{(1)}& \begin{array}{r}\begin{array}{rcl}{\sigma }_g^2(M,\tau )& =& <{\sigma }^2(M,\tau )>\\ & =& <\frac{1}{M-1}\sum _{i=1}^M({\overline{y}}_i-\frac{1}{M}\sum _{j=1}^M{\overline{y}}_j{)}^2>\end{array}\end{array}\end{array}$$

上式中$<\cdot >$表示求统计平均（期望）

Allan方差定义

Allan方差是上述通用样本方差分析在$M=2$时的一个特例，也就是在公式（1）中令$M=2$

$$\begin{array}{}\text{(2)}& \begin{array}{r}\begin{array}{rcl}{\sigma }^2(\tau )& =& <{\sigma }^2(2,\tau )>\\ & =& <\sum _{i=1}^2({\overline{y}}_i-\frac{1}{2}\sum _{j=1}^2{\overline{y}}_j{)}^2>\\ & =& <\frac{1}{2}({\overline{y}}_2-{\overline{y}}_1{)}^2>\\ & =& \frac{1}{2}<({\overline{y}}_2-{\overline{y}}_1{)}^2>\end{array}\end{array}\end{array}$$

注意，Allan方差是$\tau$的函数

Allan方差定义

角速度的Allan方差计算

这里只推导一种方式，根据IEEE Standard

令陀螺仪的角度输出为$\mathrm{\Omega }(t^{\prime })$，则角度为

$$\theta (t)={\int }^t\mathrm{\Omega }(t^{\prime })d{t}^{\prime }$$

对于离散时间形式，时间点$t=k{\tau }_0,k=1,2,3,...,N$。其中${\tau }_0$是陀螺仪的输出时间间隔，$N$为总体样本数量。则，在时间段$t_k$到$t_k+\tau$之间，角速度的平均值是：

$$\begin{array}{}\text{(3)}& {\overline{\mathrm{\Omega }}}_k(\tau )=\frac{{\theta }_{k+m}-{\theta }_k}{\tau }\end{array}$$

其中$\tau =m{\tau }_0$

将公式（3）带入上文Allan方差的定义公式（2）中可以得到：

$$\begin{array}{}\text{(4)}& \begin{array}{r}\begin{array}{rcl}{\sigma }^2(\tau )& =& \frac{1}{2}<({\overline{\mathrm{\Omega }}}_{k+m}-{\overline{\mathrm{\Omega }}}_k{)}^2>\\ & =& \frac{1}{2{\tau }^2}<({\theta }_{k+2m}-2{\theta }_{k+m}+{\theta }_k{)}^2>\end{array}\end{array}\end{array}$$

将公式（4）的期望代替为在整个数据样本（数量为N）上求平均，可以得到如下形式

$$\begin{array}{}\text{(5)}& \begin{array}{r}\begin{array}{rcl}{\sigma }^2(\tau )& =& \frac{1}{2{\tau }^2(N-2m)}\sum _{k=1}^{N-2m}({\theta }_{k+2m}-2{\theta }_{k+m}+{\theta }_k{)}^2\end{array}\end{array}\end{array}$$

重叠分段法计算Allan方差

这部分主要参考：Allan Variance: Noise Analysis for Gyroscopes 在公式（5）中，对每$m$个原始数据作为一段。切割方式有很多种，在图Allan方差定义中，每段数据都不重叠。另一种分割方法是让每段数据都重叠，如下图所示

重叠采样计算Allan方差

上图中$m=3$，也就是每4个原始数据作为一组，组与组之前相差${\tau }_0$长度。

Allan方差的计算方式与公式（5）保持一致。

Allan方差双对数曲线

根据公式（5），Allan方差是时间周期$\tau =m{\tau }_0$的函数，我们根据不同的$\tau$，计算出对应的Allan方差，并将$lo{g}_{10}(\tau )$作为x，$lo{g}_{10}(\sigma (\tau ))$作为y，画出曲线，如下图：

Allan方差双对数曲线

注意，y轴是Allan方差的开根，又叫Allan

Allan方差双对数曲线可以用不同斜率的直线拟合，不同斜率部分的直线代表不同的含义。

这里讲述一下如何从双对数曲线分析信号的不同误差成分（以陀螺仪为例）

角度随机游走/白噪声强度

角度随机游走对应到双对数曲线中斜率为$-\frac{1}{2}$的直线，取值可以直接读取$\tau =1$处的值。

零偏不稳定性

零偏不稳定性对应双对数曲线中斜率为0的直线，也就是双对数曲线中最小值位置。

速率随机游走

速率随机游走对应双对数曲线中斜率为$\frac{1}{2}$的直线，取值时取该直线与$\tau =3$的交点。

上面三个误差项是我们用得比较多的误差项，画成图如下：

三种误差项在双对数曲线中的示意图

为什么不同的斜率对应不同的误差项？？

说实话，笔者对这个部分也是一知半解，说实话，也没有必要完全弄明白，里面涉及到时域信号分析太多基础知识了。这里根据IEEE Standard里一些内容根据自己的理解简单讲述一下。 噪声的双边功率谱密度(PSD:Power Sepctral Density)与Allan方差有如下关系：

功率谱密度与Allan方差的关系

上式中可以看作Allan方差是噪声能量经过一个$si{n}^4(x)/(x{)}^2$为转移函数的滤波器，这个函数是由Allan方差计算决定的。这个滤波器的带宽跟$\tau$有关，因此，不同的$\tau$对应的Allan方差来分析噪声的不同频段特性。

以白噪声为例，其功率谱密度为

白噪声PSD

将其带入功率谱密度与Allan方差关系的公式中，并进行公式中的积分，可以得到如下结果:

白噪声PSD转换成Allan方差

对上式 ${\sigma }^2(\tau )=\frac{N^2}{\tau }$ 开根后求对数后可以得到

$$log(\sigma (\tau ))=N-\frac{1}{2}log(\tau )$$

刚好可以得到一条斜率为$-\frac{1}{2}$的直线，且当$\tau =1$的时候可以得到$N$。这也对应了上文中从Allan方差双对数曲线中求角度随机游走的方式。

应用重点

不同的IMU建模方式需要不同的系数

强调这一点是因为imu_utils工具标定的输出结果和kalibr的模型要求是不一样的。

在GTSAM中

在GTSAM中使用IMUFactor时需要IMU的几个参数：

1. IMU零偏Bias的初始值。但是这个参数可以先设置为0,后续逐步更新。（Bias是IMUFactor中的一个变量，会随时间改变）。
2. IMU加速度（3轴）和陀螺仪（3轴）的噪声参数。
3. IMU加速度（3轴）和陀螺仪（3轴）的Bias噪声参数。

上面2、3两点的参数可以使用imu_utils的标定结果。

其中第2点对应的是Allan方差双对数曲线（后文详细讲解）中斜率为-1/2的直线与 $\tau =1s$的交点（IEEE标准）。（可以直接取双对数曲线 $\tau =1s$处的值，因为"据说"惯性单元在1Hz的时候对应的噪声强度主要由白噪声强度左右。这段话出自kalibr:IMU-Noise-Model，imu_utils也是直接取这一点的值） >This is only true since the noise power in most inertial sensors is dominated by "white noise" at a frequency of approximately 1Hz.

第3点描述的是Bias的不稳定性，有些地方也叫Bias随机游走。这个参数对应的是Bias [In]Stability。Bias Stability可以直接从Allan方差双对数曲线中读出来（后文详解）

之所以GTSAM需要这些参数是因为GTSAM中对IMU的建模方法跟大部分算法相同，都是把IMU的Bias当作状态，认为Bias会随着时间变化而缓慢改变（这也是较为正确的做法）。因此需要标定出Bias的稳定性。

在[kalibr]中 在[kalibr]中对IMU的建模是使用下面的公式：

$$\tilde{w}(t)=w(t)+b(t)+n(t)$$

其中$w(t)$代表真值，$b(t)$是Bias，$n(t)$是噪声。因此在[kalibr]中使用角度随机游走系数作为$n(t)$。同时使用角速率随机游走作为$b(t)$。一开始对这点很迷惑，为什么角速率随机游走可以用来作为$b(t)$。后面我在一个issue的回答中得到答案：issues63。[kalibr]对IMU的建模形式实际上是不够准确的，因为这里将$b(t)$当作一个恒定的量求解，实际上这是一个随着时间缓慢变化的量。而角速率随机游走则可以部分用力刻画这一偏差。事实上还是很难理解

单位

First, it is important to realize that Hertz (Hz) is defined as the inverse of seconds, which means that a noise density specification of X°/s/√Hz is exactly equivalent to an angle random walk specification of X°/√s with no conversion necessary. 出自：https://www.vectornav.com/resources/inertial-navigation-primer/specifications--and--error-budgets/specs-imuspecs

参考资料

【1】IEEE Standard Specification Format Guide and Test Procedure for Single-Axis Laser Gyros, AnnexC, IEEE Std 647-1995, 1995

【2】IMU Errors and Their Effects

【3】kalibr:IMU-Noise-Model

【4】imu_utils

【5】新手入门系列3——Allan方差分析方法的直观理解

【6】Allan Variance: Noise Analysis for Gyroscopes