重要说明：本文从网上资料整理而来，仅记录博主学习相关知识点的过程，侵删。

一、参考资料

新手入门系列3——Allan方差分析方法的直观理解
惯性测量单元Allan方差分析详解
IMU误差&测量模型
IMU标定之—Allan方差
IMU误差模型简介及VINS使用说明

参考文献

[1] An overview of the Allan variance method of IFOG noise analysis
[2] Handbook of Frequency Stability Analysis
[3] Allan Variance: Noise Analysis for Gyroscopes

i2Nav

武汉大学多源智能导航实验室
i2Nav b站
b站视频：惯性导航原理
b站视频：组合导航原理

二、相关介绍

1. 白噪声

1.1 白噪声的概念

白噪声（white noise）是指功率谱密度在整个频域内是常数的噪声。 所有频率具有相同能量密度的随机噪声称为白噪声。

白噪声，是一种功率谱密度为常数的随机信号。换句话说，此信号在各个频段上的功率谱密度是一样的，由于白光是由各种频率（颜色）的单色光混合而成，因而此信号的这种具有平坦功率谱的性质被称作是“白色的”，此信号也因此被称作白噪声。

理想的白噪声具有无限带宽，因而其能量是无限大，但这在现实世界是不可能存在的。实际上，我们常常将有限带宽的平整信号视为白噪声，因为这让我们在数学分析上更加方便。白噪声在数学处理上比较方便，因此它是系统分析的有力工具。一般，只要一个噪声过程所具有的频谱宽度远远大于它所作用系统的带宽，并且在该带宽中其频谱密度基本上可以作为常数来考虑，就可以把它作为白噪声来处理。例如，热噪声和散弹噪声在很宽的频率范围内具有均匀的功率谱密度，通常可以认为它们是白噪声。

1.2 高斯白噪声

如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。即白噪声未必是高斯白噪声。

热噪声和散粒噪声是高斯白噪声。

高斯白噪声是指信号中包含从负无穷到正无穷之间的所有频率分量，且各频率分量在信号中的权值相同。白光包含各个频率成分的光，白噪声这个名称是由此而来的。它在任意时刻的幅度是随机的，但在整体上满足高斯分布函数。

2. MEMS技术

MEMS惯性传感器的研究背景与发展现状

MEMS惯性传感器

MEMS（Micro-electromechanical Systems），简称“微机电系统”，MEMS集机械元件、微传感器、信号处理与控制电路、接口电路、通信和电源为一体，一般认为MEMS是由微型机械传感器、执行机构和微电子电路组成的微机电系统。

2.1 MEMS惯性传感器简介

MEMS惯性传感器是利用MEMS微加工技术和半导体集成电路制造技术制作而成的惯性传感器，由质量块、弹性元件、敏感元件等机械结构件（传感器芯片）和信号处理电路（ASIC）等部分组成。

MEMS惯性传感器行业分类可以根据根据被测物理量，将MEMS惯性传感器分成三类：微机械加速度计（MEMS加速度计）、微机械陀螺仪（MEMS陀螺仪）、微惯性测量组合（MIMU）。

MEMS惯性传感器工作原理：MEMS惯性传感器的工作原理是经典力学中的牛顿定律。它的作用是测量运动物体(如车辆、飞机、导弹、舰船、卫星等)的质心运动和姿态运动，进而对运动物体进行控制和导航。

由于MEMS惯性传感器基本的测量功能，如加速度测量，倾斜测量，振动测量甚至旋转测量，消费电子领域有待探索的应用将不断涌现。MEMS惯性传感器可以组成低成本的GNSS/INS组合导航系统。

2.2 MEMS加速度计

MEMS加速度计是利用感测质量的惯性力进行测量的传感器，一般由标准质量块(感测元件)和检测电路组成。根据传感原理的不同，MEMS加速度计可以分为压阻式、热流式、谐振式和电容式等。压阻式MEMS 加速度计容易受到压阻材料影响温度效应严重、灵敏度低，横向灵敏度大，精度不高。热流式加速度计受传热介质本身的特性限制，器件频率响应慢、线性度差、容易受外界温度影响。因此，热流式和压阻式加速度计主要用于对精度要求不高的民用领域或军事领域中的高 g 值测量。谐振式微加速度计理论上可以达到导航级的精度，但目前技术状态还达不到实用化。而电容式硅微加速度计由于精度较高、技术成熟、且环境适应性强，是目前技术最为成熟、应用最为广泛的 MEMS 加速度计。

2.3 MEMS微机械陀螺仪

MEMS微机械陀螺仪主要有线振动型陀螺仪和谐振环型陀螺仪，线振动型陀螺仪工艺简单，利于大批量、低成本生产；谐振环型陀螺仪具有更高的理论精度，但结构及原理更为复杂。

2.4 MIMU微惯性测量组合

微惯性测量组合（Micro Inertial Messurement Unit，MIMU）是基于 MEMS技术的新型惯性测量器件，由三个MEMS加速度计和三个MEMS陀螺仪组成。MIMU用来测量物体的三轴角速度和三轴加速度信息，是实现微小型无人机、交通工具等导航制导的核心部件。

MIMU作为一种运动测量传感器，IMU的数据采集直接决定了INS的精度和动态响应性能上限。

目前，国内研制的微型惯性测量单元，主要技术路线还是采用三个单轴加速度计和三个陀螺仪立体组装方式集成，虽然产品指标能够满足一些现有领域的使用要求，但还存在体积偏大、安装精度差、成本较高等问题。

2.5 惯性微系统

惯性微系统是利用3D异构集成技术，将MEMS 陀螺、MEMS 加速度计、压力传感器、磁传感器和信号处理电路等在硅基片上进行集成，并内置导航定位算法，实现芯片级精确制导、导航、定位等功能，同时能够与卫星导航共同组成组合导航系统，是装备制导导航和定位的核心部件。

3. IMU组成部分

惯性测量单元（Inertial Measuremnt Unit，IMU）是MEMS技术结合的微惯性测量组合，因此在很多地方被称为MIMU。它主要由三个MEMS加速度传感器、三个陀螺仪和求解电路组成。

通常，IMU由以下组件组成：

• 3个加速度计；
• 3个陀螺仪；
• 数字信号处理硬件/软件；
• 功率调节器；
• 通信硬件/软件；
• 外壳。

三个加速度计相互成直角安装，这样可以在三个轴上独立测量加速度：X、Y和Z。三个陀螺仪也彼此成直角，因此可以围绕每个加速度轴测量角速率。

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

4. INS惯性导航系统

惯性导航系统(Inertial Navigation System，INS)，简称惯导，利用惯性测量单元 (IMU)的输出，将加速度和旋转信息相结合，带有位置、速度和姿态的初始信息。

5. IMU测量模型

5.1 陀螺仪测量模型

5.2 加速度计测量模型

6. IMU误差

IMU Errors and Their Effects

电子器件输出都会有误差，又叫噪声。在IMU采集数据时，会产生两种IMU误差：确定性误差和随机性误差。为获得精确的数据，需要对确定性误差和随机性误差这两种误差进行标定。Allan方差主要用于标定随机误差。

• 确定性误差：包括交轴耦合误差（Axis-misalignment），比例因子误差（Scale Factor），零偏（Bias）。其中，零偏为加性误差，交轴耦合和比例因子误差为乘性误差。
• 随机误差：主要是高斯白噪声和bias随机游走。

IMU误差主要有以下几种：

6.1 确定性误差

确定性误差，又称为系统性误差。

6.1.1 交轴耦合误差（Axis-misalignment）

交轴耦合误差是由于IMU的三个轴非正交导致的，其衡量了非正交的程度，图中绿色的是理想的正交系（定义为B系），而红色的为传感器的非正交系（定义为S系），因此，需要求解出S系到B系的映射关系，才能得到每个轴较为准确的输出。

6.1.2 比例因子误差（Scale Factor）

i2Nav

比例因子误差指的是，传感器输出值和实际值之间的比例误差。

6.1.3 零偏（Bias）

零偏指的是传感器输出与实际值之间的常值误差。

6.2 随机误差

Gyro Noise and Allan Deviation + IMU Example

6.2.1 角度随机游走（ARW）

角度随机游走（Angle Random Walk, ARW）是由陀螺仪产生的噪声，其反映了角速度信号中白噪声的特性，积分后表现为角度随机游走。

在实际情况下，当陀螺仪静止时，由于热噪声等影响，其角速度并不是一成不变的，假设其服从白噪声模型，那么其功率谱密度在整个频域内服从均匀分布。因此，对角速度积分以后，可以得到角度随时间的变化关系，也就是角度随机游走。

6.2.2 速度随机游走（VRW）

速度随机游走（Velocity Random Walk, VRM）是由加速度计产生的噪声，速度随机游走为是加速度白噪声的积分结果。

7. 陀螺零偏

新手入门系列1——如何区分惯性器件的零偏误差？

如果只用一个指标来代表一款IMU的精度的话，那毫无疑问是陀螺零偏。这是因为：

1. 惯导系统的精度主要取决于IMU中的陀螺器件精度，而不是加速度计精度；
2. 陀螺的精度指标中最重要的又是零偏误差，它基本上决定了该惯导长时间独立工作时的误差发散速度。

但是，这里需要特别注意的是，陀螺零偏有好几种，看产品指标时一定要弄清楚是哪一种陀螺零偏指标。

7.1 常值零偏

常值零偏就是这只陀螺生产出来后就一直固定不变的零偏值。对于传统的高性能惯性器件来说，该误差在出厂标定时往往就被补偿干净了，因此不会标注这个指标；但对于低端MEMS IMU芯片来说（例如用在手机中的价格不到10块钱的芯片），则不可能做逐个的标定和补偿，因此常会存在deg/s（也就是几千deg/h）量级的常值零偏。这看上去巨大无比，但我们在实际使用中很容易对付，例如在初始启动过程中利用几秒钟的静态数据求平均即可扣掉大部分。

7.2 全温零偏误差

全温零偏误差（bias error over temperature），又称为零偏的温度敏感系数（或称温漂系数），反映器件参数的温度敏感性。就是指陀螺零偏在其额定工作温度范围内相对于室温零偏值的变化量。对于传统惯性器件，一般是逐个做了温漂补偿的，因此这个全温零偏误差就是温补后的残差；而对于低端MEMS芯片，不可能逐个做温漂标定和补偿，如果给出这个全温零偏误差可能会很吓人，因此厂家往往给个零偏的温度敏感系数（例如0.01 deg/s/℃）。

零偏重复性

零偏重复性全称零偏逐次上电重复性（run-to-run repeatability），这是传统惯性器件的经典指标，是指惯性器件不同次上电运行时的零偏的不重复程度。具体测量方法是在常温下将器件多次上电，测量和记录每次上电的零偏数值，然后统计其差异。很多新手会好奇，一个传感器的零偏值每次上电的变化能有多大？这么细微的误差因素都需要考虑？这是因为传统惯性器件在出厂环节就把常值误差和温漂误差等主要误差都已仔细补偿掉了，因此逐次上电重复性这个次要误差才成为了一个不可忽视的主要误差。

7.3 零偏稳定性

严格来说应该称为零偏不稳定性（in-run instability），反映器件上电稳定后其零偏随时间变化的情况。根据具体测算方法又分为两种：

1. 我国的国军标定义的零偏不稳定性：采集几个小时的静态数据，每10秒或100秒求平均（以便抑制器件白噪声的影响），然后统计这些平均值的标准差。

2. Allan方差给出的零偏不稳定性：采集足够长时间的静态数据（一般大于10小时，越高等级的器件所需时间越长），画Allan方差曲线，取其谷底值。

前者对惯导的实际表现有比较直接的影响，有现实指导意义；而后者则只是反映器件在极端理想条件下的性能极限，缺乏现实意义。从具体数值来看，前者也比后者大几倍甚至高一个量级。

7.4 零偏的加速度敏感性

零偏的加速度敏感性英文是g-sensitivity或linear acceleration effect。陀螺的输出本来应该对加速度完全不敏感，但由于其敏感结构的加工误差等因素，多少还是会受到线加速度影响的，我们就用零偏的加速度敏感性来描述。显然，这种零偏误差只有在强动态载体上才会造成显著影响，而对于常见车载、船载低动态载体往往可以忽略。需要注意的是，越是高灵敏度的MEMS陀螺，其微机械结构的敏感质量越大，因此其加速度敏感性往往会比较大，例如某高端MEMS陀螺的加速度敏感性为18 deg/h/g。

7.5 总结

基于上述这么多种零偏误差类型，如果厂家只向用户提供一个零偏指标，那么就应该是所有零偏误差的总和，或者是所有零偏指标中最大的那个。但毫不意外的，厂家最喜欢提Allan方差零偏（尤其是MEMS厂商），因为它的数值最小，听上去最棒。而某些不负责任的商家甚至都不明确标注这是Allan方差给出的零偏，以此误导用户。因此大家一定要擦亮眼睛，认真理解具体指标及其测算方法，尤其是对MEMS器件。

**我一般会重点关注MEMS零偏的温度敏感系数这个指标，因为MEMS器件的温漂是其内在的关键误差！对温度变化不敏感的话，其它指标通常都不会太差。**然后再看零偏逐次上电重复性和零偏不稳定性（国军标）这两个指标。但MEMS器件手册往往不给这两个指标，只给Allan方差零偏不稳定性，此时无奈之下，我会把它放大5~10倍来当作国军标零偏不稳定性来用。

还有一个最容易理解但很容易被忽视的陷阱，就是给出的器件误差的具体数值含义，到底是峰峰值（p-p）、最大值（max）、RMS值或是1σ？假设某个随机误差是符合零均值正态分布的，那么会有：

峰峰值 = 2倍最大值= 6倍RMS或σ值

所以，使用不同的统计值来给出误差指标，其具体数值会有天壤之别。大家在读指标时要特别注意！毫不意外的，厂家都喜欢用1σ值，因为它最小。

​ 因此，我们调研选型惯性器件时一定要确保是把相同零偏指标的相同统计值来进行对比！千万不要因为自己的马虎而被商家忽悠。

最后想跟大家说的是，不要因为MEMS器件巨大的常值零偏和全温零偏误差而轻易否定它、拒绝它。反正MEMS惯导不会长时间单独工作，一般都是与GNSS构成组合导航系统，那么常值零偏和缓慢零偏变化（包括温漂）以及零偏重复性一般都能被组合导航算法（例如增广Kalman滤波）有效地进行在线估计和补偿，此时我们真正在乎和需要关心的是零偏不稳定性中的中短期变化成份（例如，10s~1000s时间尺度上的波动变化），而MEMS器件的这种中短期零偏不稳定性与传统惯性器件相比并不会差得太远。千万别因为我们的偏执而错过了一个物美价廉的好器件。

8. 时间序列分析方法

8.1 功率谱密度分析方法（PSD）

功率谱密度分析方法（PSD）描述误差序列在不同频率尺度上的成份，属于频域分析。PSD擅长于分析中短期误差（即中高频成份）。

8.2 自相关分析方法

自相关分析擅长于分析中短期误差（即中高频成份）。

8.3 通用的样本方差分析方法

[2]

通用的样本方差分析方法与Allan方差不同之处是在分块、取平均之后，不是采取相邻块求差，而是截取连续M块求其标准方差，然后统计这些标准方差的样本均值，作为反映1~M倍块长度这个时间尺度区间的误差指标。

容易推导，Allan方差只是上述通用样本方差分析在M=2时的一个特例，如下式：

由于Allan方差取M=2这一最小值，因此它具有最高的时间尺度分辨率，对不微附件同时间尺度上的误差成份的刻画最细腻，因此被大家用得最多。

三、GNSS/INS组合导航系统

新手入门系列4——MEMS IMU原始数据采集和时间同步的那些坑

在进行GNSS（Global Navigation Statellite System）和INS数据融合时，两者的时间同步的精度也直接影响着组合导航系统的最终导航定位精度。

GNSS/INS组合导航系统中，有两个关键点：数据采集和时间同步。

1. IMU数据采集

1.1 工业级MEMS IMU模块

一般来说，消费级MEMS芯片仅能输出速率形式的测量值，对应陀螺输出的是角速度（常用单位为deg/s），加速度计输出的是加速度（严格说应该是比力，常用单位为m/s2或g）。

工业级MEMS IMU模块，本文以ADI的ADIS16465为例。

1.2 消费级MEMS IMU芯片

对于工业级MEMS模块而言，除了能够输出速率形式的测量值以外，有些还能够输出增量形式的测量值，对应陀螺输出的是角度增量（deg），加速度计输出的是速度增量（m/s）。

事实上，离散化后的INS机械编排解算，需要的正是增量形式的测量值。因此，如果IMU提供增量形式的测量值，我们推荐直接使用；否则，我们需要人为地将速率形式的测量值转为增量形式。一般地，在采样率足够的情况下，可以使用前后两个IMU历元的速率测量值的平均值乘上时间间隔（梯形积分），即可得到对应区间的增量测量值。

消费级MEMS IMU芯片，本文以InvenSense的ICM20602为例。

2. IMU时间同步

IMU时间同步，是指不同传感器的数据都打上共同的时标（即采用统一的时间系统），而不是指所有传感器都同步采样（即采样时刻对齐）。传感器同步采样一方面很难落实，另一方面也没有必要，因为多传感器的不同步采样可以通过组合导航算法中惯导的时间传递（即状态转移）作用来解决。

IMU时间同步，对GNSS接收机和IMU之间进行时间同步。GNSS接收机本身就具备高精度授时的功能，GNSS数据都天生具有精准的GNSS时间，因此最理想的时间同步是将IMU数据也赋予GNSS时间。一般地，通过GNSS接收机输出的1PPS信号（1Hz的秒脉冲信号，也可以是其它频率），结合本地时钟或者IMU内部时钟，即可实现高精度的硬件时间同步（更确切地说是“时标统一”，精度在微秒量级）。然而，很多应用系统不具备硬件同步的条件（即无法连接1PPS脉冲信号或IMU采样信号），那么我们只能被动地采用软件同步的方式（精度只有几十毫秒量级）。

2.1 硬件同步

一般地，为了保证GNSS/INS组合导航系统的精度，都采用硬件时间同步的方式，实现GNSS信号采样时标和IMU数据采集时标的严格统一。对于GNSS接收机而言，其1PPS信号严格对应其在GNSS整秒的GNSS信号采样时刻；而对于MEMS IMU，一般都具备“数据采集状态”信号接口（具体形式为脉冲输出信号）。ICM20602有中断输出引脚（INT），可配置为采样完成（Data Ready）中断信号输出；ADIS16465则直接有一个“采样完成”信号的输出接口。

下图给出IMU硬件时间同步的示意图，嵌入式MCU（单片机）同时接收IMU的“采样完成”脉冲信号和GNSS接收机的1PPS信号，各自触发中断并记录当前的本地晶振时间。利用1PPS信号所对应的本地晶振时间（例如间隔1us的计数器）和GNSS整秒时间，即可计算出本地晶振时间相对于GNSS绝对时间的偏移值；再将这个偏移值补偿给每个IMU“采样完成”信号所对应的本地晶振时间，即可将IMU的本地时标转化为GNSS绝对时标。IMU硬件时间同步的过程如公式(1)~(3)所示。最终实现IMU数据和GNSS数据的时标统一。1PPS信号的误差一般为几十纳秒，可忽略不计；本地晶振时间一般采用微秒级别的计数器，那么最终GNSS数据和IMU数据的时间同步精度即为本地时间的精度，即微秒级。

$$\Delta t_2=(C_3-C_1)\times 1.0e-6+b_t=1.0(1)$$
可以推出：
$$b_t=1.0-(C_3-C_1)\times 1.0{\mathrm{e}}^{-6}(2)$$
那么，IMU采样时刻的绝对时标：
$$\begin{array}{rl}{t}_2& =t_1+(C_2-C_1)\times 1.0e^{-6}\\ & +b_t\times (C_2-C_1)/(C_3-C_1)(3)\end{array}$$
其中，我们假设本地时钟的额定频率为1MHz（即计数间隔为1.0e-6秒）；

$\Delta t_2$ ：两个相邻1PPS之间的物理时间长度（即1秒）；

$C_1$：某1PPS信号对应时刻的本地时钟计数值；

$C_2$：某IMU采样完成信号对应时刻的本地时钟计数值；

$C_3$：下一个1PPS信号对应时刻的本地时钟计数值；

$b_t$：本地时钟在1秒内（即两个相邻1PPS之间）的漂移值；

$t_1$：对应于C1（即1PPS时刻）的GNSS绝对时间，是由GNSS接收机提供的已知量；

$t_2$：对应于C2（即IMU采样完成时刻）的GNSS绝对时间，是待求量。

注：如果本地晶振的品质较高（例如优于10个ppm的TCXO），那么其在1秒内的漂移值bt可忽略不计。

2.2 软件同步

许多GNSS/INS系统在开发过程中，没有条件连接GNSS接收机的1PPS信号或者IMU根本不具备数据采集相关的硬件信号（诸如“触发采样”、“采样完成”等信号管脚），此时就不得不退而求其次，采用软件同步。所谓的软件同步，就是把数据处理器（例如单片机、工控机等）接收到IMU数据和GNSS数据的时刻“当作”这些数据的采样时刻，记录下本地晶振时间作为其统一的时标。当然也可以进一步通过将GNSS数据接收时刻的本地时间与该GNSS数据中的GNSS绝对时标对比来获得两者的偏差量，进而将IMU的本地时标都转化为GNSS绝对时标。

软件时间同步与硬件时间同步的道理是完全一样的，只是将精准的数据采集时刻（通过1PPS信号和“采样完成”脉冲来体现）替换为数据接收时刻（有几毫秒到几十毫秒的不确定延迟）。其中GNSS数据接收时刻会受到若干因素的影响而造成延迟，包括GNSS接收机内部的信号处理时间、接收机到嵌入式处理器的接口传输时间和嵌入式处理器的解码时间，尤其是当处理器处于忙碌状态时就更加不可控；而IMU数据接收时刻的延迟也类似，只是因为每条IMU数据量很少而使得延迟较小。总之，通过软件的方式实现的时间同步具有诸多不可控因素，导致实际上通过软件时间同步的精度低且稳定性差。

软件同步是在确实无法实现硬件时间同步的情况下的无奈选择，对于GNSS精密定位模式下的中高速载体的组合导航（例如RTK/INS车辆导航），几十毫秒的同步误差会深刻地伤害其组合导航精度，因而是无法接受的。

2.3 多源导航传感器的时间同步

在多源导航定位系统中，GNSS接收机已经成为必不可少的一部分，无论是低至几块钱的导航型GNSS接收机模块，还是价值上千元的测量型板卡，都能够提供高精度的绝对授时（误差一般为几十纳秒）。因而，GNSS接收机自然就成为了多源导航定位系统中的高精度时间基准。利用GNSS接收机的1PPS信号，还可以对精度较低的MCU本地时间系统进行校正，从而得到与GNSS时间同步的高精度MCU本地时间系统（唐海亮等，2019），进而对多源导航传感器进行时间同步（即给各传感器都打上GNSS绝对时标）。

2.3.1 相机的时间同步

对于相机而言，其光学信号采样对应着相机的曝光，因此记录相机的曝光时间就获得了相机图像的采样时刻。一般的相机都具有硬件触发信号接口，即通过外部信号触发相机曝光采样；相机同时输出一个反映曝光过程的事件信号。那么，我们可以利用经过GNSS时间校准后的本地时间系统，实现相机的时间同步。例如，每过一个固定的时间间隔，MCU输出一个相机触发信号；相机收到触发信号后，经过内部很小的处理延时后，开始曝光图像，同时输出曝光信号，曝光结束，则曝光信号结束；MCU则记录曝光开始和曝光结束的两个时间，并取均值，作为最终的相机图像采样时间，这就实现了相机图像数据的时间同步。

2.3.2 LiDAR激光雷达的时间同步

对于激光雷达而言，其时间同步方式相对简单，因为激光雷达厂商已经将大部分工作集成到激光雷达内部，为用户提供了现成的时间同步接口。一般地，我们只要给激光雷达提供GNSS接收机的1PPS信号及其对应的GNSS绝对时间消息，即可实现激光雷达点云数据的时间同步。例如，Velodyne的VLP-16激光雷达对外提供1PPS信号的输入接口和一路消息输入串口，分别连接GNSS接收机的1PPS信号和定位输出接口（输出GPRMC或GPGGA格式的NMEA消息）。而部分Livox激光雷达硬件上仅有1PPS信号输入接口，需要在上位机将额外的GPRMC消息发送给激光雷达。激光雷达通过接收的1PPS信号及其对应的GNSS定位消息中的GNSS绝对时标，再配合其内部的时间系统，即可实现激光点云数据的时间同步。

四、Allan方差

Allan方差适用于分析惯导（和时钟）看重长期稳定性的传感器，Allan方差曲线对中长期时间尺度上的误差特性有很强的表现力。Allan方差相比于普通的方差具有更好的描述长时间误差的优势，因此，Allan方差被用来当作IMU噪声标定的标准方法。

Allan方差是一个相对而言直观形象、朴实无华的概念，它当初是在实践中针对高精度时钟稳定性分析的需求而提出的。

1. 问题引入

传统的误差指标往往是采用均值误差（反映整个误差序列有无宏观偏置）、标准差（反映整个误差序列的波动情况），以及均方根（RMS，可以认为是宏观偏置与波动情况的综合）。如下图所示，它们都是反映误差序列的整体情况的指标，其中含有短时间快速变化和长时间缓变的各个成份，无法细分出不同时间尺度上的误差波动情况。

2. Allan方差的概念

Allan方差，又称为阿伦方差或阿兰方差，Allan方差法是20世纪60年代由美国国家标准局的David Allan提出的基于时域的分析方法。

Allan方差是时频分析和惯性导航领域常用的一种误差分析方法，它有效地刻画了待研究的误差时间序列在不同时间尺度上的波动水平（不稳定性），并可根据不同时间尺度上的Allan方差值所构成的曲线的形状特征来辨识其中包含的随机过程模型。Allan方差分析方法对中长期的随机波动具有很强的表现力，它完全可以作为一个通用的时间序列分析工具来推广到其它应用领域。

在我们对惯性导航器件（陀螺和加速度计）进行误差分析时，常采用Allan方差分析方法。这是当初从时频领域（高精度时钟的频率稳定性）借鉴来的一种时间序列分析方法 ，非常适合于对中长期波动（不稳定性）进行定量描述和分析，因此正是我们惯性器件误差分析所需要的。

3. Allan方差计算过程

Allan方差是将误差序列在某个指定的时间尺度上的波动情况进行了精确提取，其具体计算步骤如下：

1. 保持IMU静止， 采集N个数据点。将采样数据划分为包含不同数量采样点的子集 : $T(n)=n{T}_s，n=1,2\dots N/2$，$T(n)$ 为该子集的平均时间。
   

2. 将N个点定长分组，常用的交叠式Allen方差如下(n=3)：
   

3. 对于每个子集，计算均值后，计算相邻两组的方差。这样就会得到$n$个在每一平均时间$T(n)$对应的方差$\sigma (n)$。根据这$n$个点作图，可得到Allan标准差$\sigma$随平均时间$T$变化的双对数曲线。Allan方差有效反映了相邻两个采样段内平均频率差的起伏。
   

4. Allan方差图读取误差系数

Allan方差法可用于5种随机误差的标定：

• 量化噪声：误差系数为 $Q$，Allan方差双对数曲线上斜率为-1的线段延长线与t=1的交点的纵坐标读数为 $\frac{\sqrt{3}}{T}Q$；
• 角度随机游走：其误差系数 $N$，Allan方差双对数曲线上斜率为$-1/2$的线段延长线与$t=1$交点的纵坐标读数即为$\frac{N}{\sqrt{T}}$；
• 零偏不稳定性：其误差系数 $B$，Allan方差双对数曲线上斜率为0的线段延长线与t=1交点的纵坐标读数为$\sqrt{\frac{2\ln 2}{\pi }}B$ ，一般常取底部平坦区的最小值或取$t=10^1$或 $t=10^2$ 处的值；
• 角速率随机游走：其误差系数 $K$，斜率为1/2的线段延长线与$t=1$交点的纵坐标读数为 $\frac{K}{\sqrt{T/3}}$
• 角速率斜坡：其误差系数 $R$，斜率为1的线段延长线与$t=1$交点的纵坐标读数为$\frac{RT}{\sqrt{2}}$。
  

假设各种误差源统计独立，那总的艾伦方差为各种误差源之和，即将量化噪声的平方${\sigma }_Q$、角度随机游走的平方${\sigma }_{RAW}$、零偏不稳定性的平方${\sigma }_{bias}$、角速率随机游走的平方${\sigma }_{RRW}$、角速率斜坡的平方${\sigma }_{RR}$的总和。

5. Allan方差分析

如果我们只对某一时间尺度上的误差（即误差波动）感兴趣的话，那么比这个时间尺度更小的细节变化（短时间快速跳动）和比这个时间尺度更大的宏观变化（长时间缓慢漂移）就都不关心了，希望在我们的误差指标中都被消除掉。而Allan方差是这样做到的：

1. 通过分块确定所要考察的时间块长度；

2. 利用块内求平均的办法把短于块长度的那些快速变化成份（细节）都抹掉；

3. 再利用相邻块求差的办法把长于两块长度的那些缓慢变化成份（宏观）都抹掉；

4. 最后对差值序列统计其方值（这是处理任何随机样本的标准操作），这样统计出来的就是介于1倍块长度和2倍块长度这样一个很窄的时间尺度范围内的误差波动情况。

6. Allan方差曲线

如果我们对误差序列的各个时间尺度上的成份都感兴趣的话，可以将块长度由短到长，“扫描”一遍，得出一组Allan方差值，然后画个“Allan方差 vs. 块长度”的曲线，这样就能全面反映被研究的误差序列的特性了。具体的，实际上是“Allan方差的开方（Allan Deviation） vs. 块长度”的双对数曲线，以便在更大的范围上有更强的表现力。以下是一张经典的高精度陀螺的Allan方差示意图。

从Allan方差曲线上，我们可以根据曲线的形状特征来识别出不同类型的随机误差（即随机过程模型）并提取其参数。例如：斜率为-1/2的直线段代表白噪声。具体可参见IEEE制定的一个关于陀螺测试的技术标准中的附录C内容[1]。这个附录中还给出了一个从Allan方差曲线中分析陀螺的多种随机误差的案例（如上图）。但我想提醒大家，实际的Allan方差曲线往往只能表现出少数两三个主要误差类型，因为其它误差都被这几个主要误差给淹没了。而由于我们在工程上关心的恰恰也就是主要误差源，因此我们并不受此困扰。更多的使用Allan方差分析过程中的注意事项可参考西工大严恭敏老师的博文《Allan方差分析的使用要点》，写得非常精准务实。

7. Allan方差的应用

Allan方差分析方法并不局限于对惯性器件和时钟误差的分析，也可推广到其它误差序列分析；而且它也不局限于是对误差序列进行分析，也可应用于任何时间序列的分析（例如建筑形变、地质演化等）。

8. Allan方差与ROS

下面是读取bag包，并获取Allen方差的图片，会生成对应的Allen方差图。这里主要参考了allan_variance_ros 对应的C++文件。

sudo pip install allantools

#!/usr/bin/env python
import rospy
import sys
import allantools
import rosbag
import numpy as np
import csv
import rospkg
import os
import matplotlib.pyplot as plt  # only for plotting, not required for calculations
import math

def getRandomWalkSegment(tau,sigma):

    m = -0.5 # slope of random walk
    """""""""""""""""""""""""""""""""
    " Find point where slope = -0.5 "
    """""""""""""""""""""""""""""""""
    randomWalk = None
    i = 1
    idx = 1
    mindiff = 999
    logTau = -999
    while (logTau<0):
        logTau = math.log(tau[i],10) 
        logSigma = math.log(sigma[i],10)
        prevLogTau = math.log(tau[i-1],10)
        prevLogSigma = math.log(sigma[i-1],10)
        slope = (logSigma-prevLogSigma)/(logTau-prevLogTau)# 随机漫步的斜率
        diff = abs(slope-m)# 当前斜率与目标斜率的差值
        if (diff<mindiff):
            mindiff = diff# 更新最小差值
            idx = i
        i = i + 1

    """"""""""""""""""""""""""""""
    " Project line to tau = 10^0 "
    """"""""""""""""""""""""""""""
    x1 = math.log(tau[idx],10)# 当前点的横坐标
    y1 = math.log(sigma[idx],10)# 当前点的纵坐标
    x2 = 0
    y2 = m*(x2-x1)+y1

    return (pow(10,x1),pow(10,y1),pow(10,x2),pow(10,y2))

def getBiasInstabilityPoint(tau,sigma):
    i = 1
    while (i<tau.size):
        if (tau[i]>1) and ((sigma[i]-sigma[i-1])>0): # only check for tau > 10^0
            break
        i = i + 1
    return (tau[i],sigma[i])

def main(args):

    rospy.init_node('allan_variance_node')

    t0 = rospy.get_time()

    """"""""""""""
    " Parameters "
    """"""""""""""
    bagfile = rospy.get_param('~bagfile_path','~/data/static.bag')# 输入的bag文件路径
    topic = rospy.get_param('~imu_topic_name','/imu')# 输入的imu topic名称
    axis = rospy.get_param('~axis',0)# 输入的轴，0为x轴，1为y轴，2为z轴
    sampleRate = rospy.get_param('~sample_rate',100)# 输入的采样频率
    isDeltaType = rospy.get_param('~delta_measurement',False)# 是否是delta measurement
    numTau = rospy.get_param('~number_of_lags',1000)# 输入的tau数目
    resultsPath = rospy.get_param('~results_directory_path',None)# 输出的结果路径

    """"""""""""""""""""""""""
    " Results Directory Path "
    """"""""""""""""""""""""""
    if resultsPath is None:
        paths = rospkg.get_ros_paths()
        path = paths[1] # path to workspace's devel
        idx = path.find("ws/")# 找到路径
        workspacePath = path[0:(idx+3)]# 取到workspace的路径
        resultsPath = workspacePath + 'av_results/'# 结果输出路径

        if not os.path.isdir(resultsPath):
            os.mkdir(resultsPath)

    print("\nResults will be save in the following directory: \n\n\t %s\n"%resultsPath)

    """"""""""""""""""
    " Form Tau Array "
    """"""""""""""""""
    taus = [None]*numTau# 初始化tau数组

    cnt = 0;
    for i in np.linspace(-2.0, 5.0, num=numTau): # lags will span from 10^-2 to 10^5, log spaced
        taus[cnt] = pow(10,i)# 将tau数组赋值,维度在10^-2 到 10^5
        cnt = cnt + 1

    """""""""""""""""
    " Parse Bagfile "
    """""""""""""""""
    bag = rosbag.Bag(bagfile)

    N = bag.get_message_count(topic) # 在bag文件中找到该topic的消息数量

    data = np.zeros( (6,N) ) # 初始化数据矩阵，维度为6*N

    if isDeltaType:
        scale = sampleRate
    else:
        scale = 1.0

    cnt = 0
    for topic, msg, t in bag.read_messages(topics=[topic]):# 遍历bag文件中的消息
        data[0,cnt] = msg.linear_acceleration.x * scale
        data[1,cnt] = msg.linear_acceleration.y * scale
        data[2,cnt] = msg.linear_acceleration.z * scale
        data[3,cnt] = msg.angular_velocity.x * scale
        data[4,cnt] = msg.angular_velocity.y * scale
        data[5,cnt] = msg.angular_velocity.z * scale
        cnt = cnt + 1

    bag.close()

    print ("[%0.2f seconds] Bagfile parsed\n"%(rospy.get_time()-t0))

    """"""""""""""""""
    " Allan Variance "
    """"""""""""""""""
    if axis is 0:
        currentAxis = 1 # 循环通过所有轴1-6
    else:
        currentAxis = axis # 只需循环一次，然后中断

    while (currentAxis <= 6):
        # taus_used 对应了是否可以使用taus的数组，adev是allan deviation degree的缩写，为allan偏差；adev_err是allan偏差的误差；adev_norm是allan偏差的标准化;
        (taus_used, adev, adev_err, adev_n) = allantools.oadev(data[currentAxis-1], data_type='freq', rate=float(sampleRate), taus=np.array(taus) )# 计算allan variance

        randomWalkSegment = getRandomWalkSegment(taus_used,adev)# 计算random walk segment
        biasInstabilityPoint = getBiasInstabilityPoint(taus_used,adev)# 计算bias instability point

        randomWalk = randomWalkSegment[3]# 获取random walk segment的纵坐标
        biasInstability = biasInstabilityPoint[1]# 获取bias instability point的纵坐标

        """""""""""""""
        " Save as CSV "
        """""""""""""""
        if (currentAxis==1):
            fname = 'allan_accel_x'
            title = 'Allan Deviation: Accelerometer X'
        elif (currentAxis==2):
            fname = 'allan_accel_y'
            title = 'Allan Deviation: Accelerometer Y'
        elif (currentAxis==3):
            fname = 'allan_accel_z'
            title = 'Allan Deviation: Accelerometer Z'
        elif (currentAxis==4):
            fname = 'allan_gyro_x'
            title = 'Allan Deviation: Gyroscope X'
        elif (currentAxis==5):
            fname = 'allan_gyro_y'
            title = 'Allan Deviation: Gyroscope Y'
        elif (currentAxis==6):
            fname = 'allan_gyro_z'
            title = 'Allan Deviation: Gyroscope Z'

        print ("[%0.2f seconds] Finished calculating allan variance - writing results to %s"%(rospy.get_time()-t0,fname))

        f = open(resultsPath + fname + '.csv', 'wt')

        try:
            writer = csv.writer(f)
            writer.writerow( ('Random Walk', 'Bias Instability') )
            writer.writerow( (randomWalk, biasInstability) )
            writer.writerow( ('Tau', 'AllanDev', 'AllanDevError', 'AllanDevN') )
            for i in range(taus_used.size):
                writer.writerow( (taus_used[i],adev[i],adev_err[i],adev_n[i])  )
        finally:
            f.close()

        """""""""""""""
        " Plot Result "
        """""""""""""""
        plt.figure(figsize=(12,8))
        ax = plt.gca()
        ax.set_yscale('log')
        ax.set_xscale('log')

        plt.plot(taus_used,adev)
        plt.plot([randomWalkSegment[0],randomWalkSegment[2]],
                 [randomWalkSegment[1],randomWalkSegment[3]],'k--')
        plt.plot(1,randomWalk,'rx',markeredgewidth=2.5,markersize=14.0)
        plt.plot(biasInstabilityPoint[0],biasInstabilityPoint[1],'ro')

        plt.grid(True, which="both")
        plt.title(title)
        plt.xlabel('Tau (s)')
        plt.ylabel('ADEV')

        for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                    ax.get_xticklabels() + ax.get_yticklabels()):
            item.set_fontsize(20)

        plt.show(block=False)

        plt.savefig(resultsPath + fname)

        currentAxis = currentAxis + 1 + axis*6 # increment currentAxis also break if axis is not =0

    inp=input("Press Enter key to close figures and end program\n")

if __name__ == '__main__':
  main(sys.argv)

五、IMU标定工具

手写VIO-使用Allen方差工具标定IMU

多传感器融合定位（3-惯性导航原理及误差分析）6-Allan方差 进行随机误差分析

1. imu_utils工具

imu_utils

使用imu_utils工具生成IMU的Allan方差标定曲线

VIO 中 IMU 的标定流程 (1/3) - imu_utils 使用备忘

imu_utils是用于求取随机误差的开源工具。

2. kalibr工具

kalibr

VIO 中 IMU 的标定流程 (2/3) - kalibr_allan 使用备忘

Kalibr标定工具箱使用详细过程

Kalibr由ETHZ视觉组开源。

Kalibr可以解决以下校准问题：

1. 多摄像机校准：具有非全局共享的重叠视场的摄像机系统的内部和外部校准。
2. 视觉惯性校准（Camera-IMU）：带相机系统的IMU的时空校准。
3. 卷帘相机校准：卷帘相机的完整固有校准（投影，失真和快门参数）。

Kalibr工程有两种，一种是已经编译好的包，叫CDE package,另一种是未经编译的源文件。

前者：安装简单，不需要依赖ROS，但是功能不全。

后者：安装稍麻烦，但功能全，建议安装这种，基于ROS。

3. imu_tk工具

『imu_tk』工具标定IMU的基本原理

4. gnss-ins-sim

GNSS-INS-SIM

gnss-ins-sim是一个由Aceinna新纳公司开发的GNSS/INS仿真项目，可生成参考轨迹、IMU传感器输出、GPS输出、里程表输出和磁强计输出用户选择/设置传感器模型，定义路径点并提供算法，以及gnss-ins- sim可以生成算法所需的数据，运行算法，绘制
仿真结果，保存仿真结果，并将结果可视化。