做IMU标定的人应该都有过这种经历陀螺和加速度计静态采集出来的数据看起来就是一条带毛刺的直线零偏也扣过了六面静态标定也做了转到导航解算里却还是能看到姿态和位置慢漂。这堆“毛刺”并不是同一种东西有的积分后变成随机游走有的带有低频相关性直接表现为零偏不稳定。想要把这些成分区分开最常用的手段就是Allan方差Allan Variance。Allan方差是1966年David Allan研究原子钟频率稳定性时提出的方法后来被惯性导航领域大规模拿来分析IMU的随机噪声。它最大的优势在于不需要转台、不需要额外激励只要有一段静止放置时采集的数据通过画一条双对数曲线就能把角度随机游走、零偏不稳定性、速率随机游走这些参数一个个读出来。这篇文章我会按从头到尾的实操逻辑讲清楚包括数据怎么采、代码怎么写、曲线怎么判读、参数怎么换算以及我实际踩过的一些坑。不管你是做组合导航、LIO/VIO还是相机IMU联合标定这套方法都能用得上。1. 为什么IMU标定流程里总会出现Allan方差1.1 IMU噪声不是同一回事很多人刚开始接触IMU时以为零偏就是零偏扣掉一个固定值就结束了。实际测试下来会发现陀螺仪静止时的输出不是一条绝对平稳的直线它会在某个均值附近随机波动。这里面既有测量电路带来的白噪声也有低频漂移还有温度变化引起的慢变分量。如果只用一个方差去描述等于把所有噪声搅成一锅粥没法指导后续的滤波设计。工程上更合理的做法是把随机噪声拆成几个独立的成分白噪声积分后形成角度随机游走低频的零偏漂移形成零偏不稳定性还有更慢的速率斜坡和量化噪声。不同的噪声成分在时间尺度上的行为完全不同。比如白噪声在短时间尺度内占主导而零偏不稳定性在长时段采样里会越来越明显。正因如此单靠一个静态方差值根本说明不了问题必须有一个能把“时间尺度”考虑进去的工具Allan方差就是干这个的。1.2 Allan方差到底在算什么理解Allan方差的关键在于“簇时间”这个概念。假设你以200Hz采样持续采集了2个小时数据量大概有144万个点。现在把这串数据按时间窗口切成一段段每一段的长度叫做簇时间τ。比如τ1秒就是把数据切成1秒一个簇然后对每个簇求平均值。求完平均值之后你会得到一串“簇平均值”序列。Allan方差做的事情很简单把相邻两个簇平均值的差求出来平方后取平均再除以2。这个值会随着τ的变化而变化。当τ从0.01秒一直增加到几千秒时不同噪声成分会依次在曲线不同位置显露出来画成双对数图后每段斜率对应的就是一类噪声。为什么这里非要除以2因为相邻两个簇平均值之差同时混入了两个簇各自的噪声波动除以2之后得到的才是“单个簇平均值”的波动方差。这个细节初学者很容易忽略但它恰恰是Allan方差和普通方差的本质区别。普通方差描述的是数据偏离均值的程度Allan方差描述的是“不同时间尺度下均值随时间变化的剧烈程度”。1.3 为什么标定流程绕不开它在做IMU误差建模时确定性误差和随机误差是分开处理的。确定性误差包括零偏、刻度因子误差、安装失准角这些大多可以用六面静态标定、多位置标定、转台标定等手段求出来。随机误差则不一样它带有统计特性需要从长时间静态数据里估计。Allan方差就是目前业界接收度最高的随机误差参数提取方法。更关键的是后面接滤波器时需要IMU噪声参数作为先验。松组合导航的EKF、紧组合导航、LIO/VIO中的预积分都要用到角度随机游走和零偏不稳定性作为噪声协方差矩阵的初值。如果这几个参数拍脑袋乱填滤波器收敛速度和稳态精度都会受影响。所以实际项目里拿到新IMU的第一件事往往是“先静态采几个晚上跑一遍Allan方差”拿到参数之后再谈其他标定。2. 原始数据怎么采直接决定分析上限2.1 采集环境与硬件安置Allan方差分析的假设是IMU在静止环境中输出只包含自身噪声。这个“静止”的要求听起来简单实际上有很多细节。首先设备必须放在稳定的平面上不要放在桌上因为桌子的微弱震动会被高精度IMU敏锐捕捉到。我一般会找一块大理石平台或者厚实的铸铁平台上面垫一层泡沫棉再把IMU放上去。如果是测试板卡还要注意风扇之类旋转部件带来的振动。其次采集过程中周围不要有人走动不要开关门空调风口不要正对设备。温度变化对MEMS IMU的影响非常明显风吹一下可能造成局部温漂反映在Allan曲线上就是长τ段明显上翘看起来像速率随机游走实际上是环境干扰。所以有条件的话最好把设备放在封闭柜子里或者等实验室人员下班之后开始采集。2.2 采样率、时长和放置位置的工程取舍采样率方面一般MEMS IMU用100Hz到200Hz足够了光纤陀螺或更高精度的器件也不需要刻意提高。Allan方差曲线的横轴是簇时间采样率主要影响最短可分析的τ。比如200Hz采样最短簇时间能到0.005秒但实际上我们更关心τ从0.1秒到几千秒这一段。过高的采样率只会增加数据量对分析结果帮助不大。采集时长则是决定长τ段可靠性的关键。为了画出完整的Allan曲线理论上采集时长最好是最长簇时间的10倍以上。如果希望看到τ1000秒处的零偏不稳定性平台至少需要采集2~3小时。我做工程分析时有个经验值MEMS IMU至少采4小时高精度光纤IMU最好采12小时以上。实际操作中经常是晚上睡觉前开始采第二天早上起来收数据一晚上正好8小时左右足够覆盖绝大多数情况。放置方向也有讲究。为了同时分析陀螺和加速度计一般让IMU的z轴朝上或朝下因为MEMS加速度计通常对重力方向敏感不同姿态下的温漂特性会有差异。如果是双轴或三轴器件建议采集一组正放的数据再采一组侧放或倒放的数据做交叉验证。不过Allan方差分析本身不依赖姿态正放侧放的结果差别不会太大真正影响大的是振动和温变。2.3 数据预处理与温漂处理拿到原始数据后不要直接丢给脚本算Allan方差。我的习惯是先看一眼完整波形的趋势确认没有明显的掉线、跳变和异常尖峰。如果中间有插拔线材或者重启设备最好直接把这段数据抠掉不要拼接。现场采集时也一定要记录下时间戳方便后面定位异常段。温漂处理是个容易纠结的问题。严格来说Allan方差分析要求信号是平稳随机过程温度斜坡会导致均值缓慢变化这种趋势项会把长τ端的曲线抬高造成“假速率斜坡”的误导。有两种处理方式一种是在数据预处理阶段对整段序列做去趋势只保留零均值附近的波动另一种是把数据按温度变化区间切段只取温度相对平稳的那一段。我更推荐后者因为去趋势可能把真实的低频零偏漂移也一并滤掉后面读取零偏不稳定性时就不准了。3. 从原始序列算出Allan方差的完整过程3.1 通俗手算版把“簇平均”讲明白先拿手工推导来理解计算过程。设原始数据为x₁, x₂, …, x_N采样周期为τ₀。选择簇时间τ m·τ₀也就是把每m个连续数据点组成一个簇。第k个簇的平均值为ȳ_k (1/m)·Σ_{i1}^{m} x_{(k-1)mi}。得到一串簇平均值序列之后相邻两个平均值的差为d_k ȳ_{k1} - ȳ_kAllan方差就定义为σ²(τ) (1/2) · E[(ȳ_{k1} - ȳ_k)²]实际计算时用样本均值代替期望值对d_k的平方求和再除以2。如果数据量足够同时为了充分利用数据实践中更常使用“重叠式Allan方差”也就是让簇窗口逐个滑动而不是不重叠地切分。重叠式公式为σ²(τ) 1 / [2(N - 2m)] · Σ_{k1}^{N-2m} (x_{k2m} - 2x_{km} x_k)²这个公式利用二次差分直接计算比先求簇均值再两两相减的做法更高效而且数值更加平滑。从数学上看x_{k2m} - 2x_{km} x_k 就是三个相距m的点的二阶差分它等价于相邻两个簇平均值之差的某种扩展形式。这也是很多开源工具内部实际采用的实现方式。3.2 Python实现的脚本与注释用Python实现重叠式Allan方差非常简洁。我通常只依赖numpy和matplotlib不需要任何第三方专用库。下面这个脚本可以直接用来处理陀螺或加速度计的单轴静态数据import numpy as np import matplotlib.pyplot as plt def overlapping_allan_deviation(data, fs200.0, max_pointsNone): data np.asarray(data, dtypenp.float64) n data.size if max_points is None: max_points int(n / 4) tau_list [] adev_list [] m 1 while 2 * m n and m max_points: # 重叠式Allan方差二次差分法 diff data[2*m:] - 2.0 * data[m:-m] data[:-2*m] avar np.mean(diff**2) / (2.0 * (n - 2*m)) tau_list.append(m / fs) adev_list.append(np.sqrt(avar)) m m * 2 # 簇时间按2倍递增画对数图足够 return np.array(tau_list), np.array(adev_list) # 加载静态数据假设只有一列单位是 deg/s gyro_x np.loadtxt(gyro_x_static.txt) fs 200.0 tau, adev overlapping_allan_deviation(gyro_x, fsfs) plt.loglog(tau, adev, markero, linestyle-) plt.xlabel(Cluster Time τ (s)) plt.ylabel(Allan Deviation σ(τ)) plt.grid(True, whichboth, linestyle--, linewidth0.5) plt.show()这段代码把簇时间按2的幂次取点画出来的曲线在对数坐标下均匀分布不会在某些区段挤成一团。如果你希望更精细的曲线可以把 m m * 2 改成 m m 1但计算量会大很多对于几小时的数据来说没必要。实际我用下来按2倍递增已经能看清所有关键斜率段了。3.3 曲线拟合与参数换算Allan方差画出来之后横轴是簇时间τ纵轴是Allan标准差σ(τ)双对数坐标。在典型IMU噪声模型下曲线会呈现几个不同斜率的区段。工程上最常读取的参数有三个角度随机游走、零偏不稳定性和速率随机游走。读取方式并不需要做复杂拟合直接在对应的直线段上读值即可。角度随机游走对应斜率-1/2的区段。在这个区段内满足σ(τ) N / √τ其中N就是角度随机游走系数。所以直接看τ1秒处曲线的纵坐标就是角随机游走数值。需要注意的是单位换算如果原始数据单位是deg/s那么读出来的ARW单位是deg/s/√Hz习惯上再乘以3600的平方根即60变成deg/√h。加速度计对应的参数叫速度随机游走VRW单位是m/s/√h。零偏不稳定性对应曲线斜率接近0的“平台段”也就是曲线最低点的位置。这个平台段的纵坐标值近似代表了陀螺零偏的波动幅度单位是deg/h。严格意义上零偏不稳定性的精确值是平台值除以0.664但很多工程实现图省事直接读平台值。对于一致性要求不高的场景这个近似够用如果要做严谨误差模型建议还是按0.664换算。4. 判读曲线时最容易踩的坑4.1 采集时间不够长后半段曲线都是幻觉Allan方差曲线的长τ端依赖长时间数据。如果你只采了20分钟数据却想看到τ200秒以后的零偏不稳定性平台那基本是不可能的。因为最长可用簇时间通常限制在总时长的四分之一到三分之一超过这个范围统计样本数量急剧减少曲线会剧烈抖动读出来的数没有意义。我见过有人用10分钟数据算出来的曲线到后面直线飙升还以为是速率斜坡很大实际上只是采样段太短导致的方差估计偏差。所以判断一条Allan曲线能否采信先看横轴最右端对应的τ值是否达到了你关心的尺度。如果关心零偏不稳定性至少要保证曲线在最低平台附近有若干个连续的采样点。否则就老老实实加长采集时间。我在实际操作中宁可采12个小时也不愿意用外推或拟合去“脑补”长τ段行为。4.2 温度、安装与数据截取干扰温度波动是MEMS传感器长τ段曲线走样最常见的原因。设备刚上电时内部温度还没稳定零点在缓慢漂移这段数据会引入明显的低频趋势。我的做法是上电后先预热15~30分钟再开始采集让温度场基本稳定。如果项目时间紧也可以直接采集完整过程但在预处理时把前10分钟的数据切掉只分析后半段。对比过几次预热后的零偏平台明显更平。安装松动也会造成类似问题。IMU和支架之间如果有微小的机械间隙热胀冷缩或者外界轻微振动都会导致传感器姿态发生细微变化这种变化在Allan曲线里会被“误判”为零偏不稳定性或速率随机游走。固定设备时最好用螺钉加胶固定而不是双面胶或直接放在平面上。机械上的微小问题最终都会在统计曲线上以噪声的形式暴露出来。4.3 别把确定性误差和随机噪声混为一谈Allan方差分析的基本前提是已经完成了扣除零偏的预处理。如果陀螺数据里包含一个大的固定零偏直接做Allan方差曲线形态不会有本质区别因为Allan方差考查的是均值的变化而不是均值的大小。但如果数据里存在明显的标定残差、刻度因子误差或者安装失准角这些误差在整段采集过程中不怎么变化只会轻微影响曲线绝对值而不会改变斜率结构。反过来说Allan方差分析并不能替代传统多位置标定或转台标定。多位置标定解决的是比例因子、安装误差、零偏残差这些确定性项Allan方差解决的是随机噪声建模。两者是互补关系我在实践中通常先做多位置标定再做长时间静态数据分析和温漂测试。如果只跑Allan方差而不做确定性标定那只是在噪声层面做了工作IMU系统精度还是上不去。4.4 与外参标定、多位置标定的分工这几年相机IMU联合标定、Lidar-IMU外参标定很流行很多人拿到新IMU就把精力全放在外参标定工具上。实际上外参标定解决的是传感器坐标系之间的旋转平移关系跟IMU自身噪声参数完全是两码事。开源标定工具通常会要求你提供IMU噪声参数甚至有些工具内部会默认给定一个初值如果直接用默认值状态估计里的协方差矩阵就是错的精调外参也很难解决。我的做法是把流程拆开第一步做温度试验和Allan方差分析把角度随机游走和零偏不稳定性这些底噪参数记录下来第二步做多位置静态标定求解零偏、刻度因子和安装误差第三步再做相机或激光雷达与IMU的外参标定。顺序不能反因为外参标定过程中的轨迹估计依赖IMU内参和噪声模型内参不准外参结果也会被带偏。5. 拿到Allan方差参数之后怎么用5.1 给状态估计的Q矩阵一个靠谱初值组合导航和LIO系统里IMU噪声参数直接进入过程噪声协方差矩阵Q。角度随机游走影响姿态预测中的角速度白噪声强度零偏不稳定性影响零偏随机游走模型的强度。这两个值设置太大会让滤波器过于信任观测、对IMU预测不信任导致高频抖动设置太小则会让滤波器过度信任IMU观测更新滞后甚至出现发散。用Allan方差得到参数之后还需要做一个简单的换算填充到系统模型中。以EKF为例陀螺角度随机游走N_g通常转化成连续噪声谱密度然后乘上离散化时间步长再填入姿态估计对应的协方差块。这个换算公式很多开源库已经封装好了但输入参数的来源往往写得很模糊。我自己会写一个配置模板把所有IMU器件的Allan参数、采样率、离散化方式都备注清楚免得几个月后回头看代码时一头雾水。5.2 器件选型与标定质量的横向对比Allan方差分析还有一个很实用的场景横向对比不同IMU器件的底噪水平。选型阶段如果只看数据手册的零偏稳定性指标往往会发现手册写的是“典型值”实际器件之间离散度很大。通过同一套静态采集和Allan方差分析流程把几款候选IMU的ARW和零偏不稳定性排在一起谁好谁差一目了然。用来验证标定质量也很直观。比如同一颗IMU在温度补偿前后分别跑Allan方差如果补偿后的零偏不稳定性平台明显下降说明温补算法确实把低频漂移压下去了。如果补偿前后曲线几乎没变化那就要检查温补模型是否真的起作用或者标定数据本身有没有问题。这个方法不需要额外设备完全是数据分析手段但比单纯看RMS误差更能定位问题。5.3 与各类标定工具配合的使用心得最后聊一下实际跟标定工具链配合的经验。做Lidar-IMU标定时很多流程要求先初始化IMU噪声模型。过去我偷懒直接用默认参数后来发现点云配准位姿里的高频抖动明显重新跑一遍Allan方差更新参数后才好转。原因很简单预积分协方差给错了后端优化就会对IMU预积分产生错误信任漂移就被分配到了外参或者位姿里。相机IMU联合标定也是一样标定包里的优化过程中需要IMU积分约束Allan参数对VIO偏移量的可观测性有影响。严格来说如果你在一个系统里同时估计陀螺零偏和加速度计零偏噪声参数给得太乐观会导致零偏收敛过慢给得太悲观又会把真实运动误识别为零偏。所以我建议每个新项目开始前都重新针对当前硬件做一次Allan方差分析和内参标定不要跨项目复用旧参数。我个人的习惯是把Allan方差分析脚本固化成一个独立的小工具每次拿到新IMU或者固件版本有更新时花一晚上采静态数据第二天花十几分钟出曲线、读参数。这套流程成本很低但能帮你省掉后期排查“系统精度为什么上不去”的几天时间。Allan方差不是万能的它针对的是传感器自身的随机噪声特性解决不了安装误差、外参错误和算法缺陷。但它确实给了你一个把噪声水平讲清楚的依据至少在跟同事讨论IMU选型、标定质量和滤波参数时它是个能说服人的数字。