IMU噪声参数辨识-艾伦方差
IMU噪声参数辨识-艾伦方差
- 前言
- 推导过程
- 1. 量化噪声（Quantization Noise，QN）
- 2．角度随机游走（Angular Random Walk，ARW）
- 3．零偏不稳定性噪声（Bias Instability，BI）
- 4．角速率随机游走（Rate Random Walk，RRW）
- 5．速率斜坡（Rate Ramp，RR）
- 5种艾伦方差汇总介绍与单位转换
- 艾伦方差辨识参数的误差问题
- 艾伦方差曲线的绘制方法
- 艾伦方差相关代码
- 艾伦方差应用实例
- 1.ADIS16465
- 1.ADIS16470
- 3.某型号光纤惯导
- 10s平滑计算零偏不稳定性
- 参考文献
前言
在统计学中描述随机变量的两个经典参数是均值和方差，早期在定量表征原子钟的频率稳定度时采用的就是经典方差方法。1996年，学者D.W.Allan在分析铯原子钟频标的频率稳定度时发现经典方差随着时间的增长而发散，为了解决该问题，提出了一种新的评定方法，后来称为艾伦方差。由于惯性器件也具有振荡器的特征，Allan方差分析也被广泛应用于惯性器件的随机误差建模，IEEE标准中就将Allan方差方法引入到了激光陀螺的建模分析。
Allan方差的物理意义以及应用本质来源于它与功率谱之间的关系，本文首先从功率谱入手推导艾伦方差，然后教大家如何从艾伦方差log曲线中辨识IMU的各种噪声，最后会给出Allan方差的实例应用，包括光纤惯导和Mems惯导，并且给出了和10s平滑方法的差异。
推导过程
因为内容较多，这里只会给出一个大概的推导流程，更多细节的推导可以参考严恭敏老师《惯性仪器测试与数据分析》一书。
从时域测量方面来分析，Allan方差定义如下:
σy2(τ)=12<[yi+1(τ)−yi(τ)]2>\sigma_{y}^{2}(\tau)=\frac{1}{2}<\left[y_{i+1}(\tau)-y_{i}(\tau)\right]^{2}>σy2(τ)=21<[yi+1(τ)−yi(τ)]2>
然后是艾伦方差和功率谱密度之间的关系，如下所示：
σy2(τ)=2∫−∞∞Sy(f)sin4(πfτ)(πfτ)2df=4∫0∞Sy(f)sin4(πfτ)(πfτ)2df\sigma_{y}^{2}(\tau)=2 \int_{-\infty}^{\infty} S_{y}(f) \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} f=4 \int_{0}^{\infty} S_{y}(f) \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} fσy2(τ)=2∫−∞∞Sy(f)(πfτ)2sin4(πfτ)df=4∫0∞Sy(f)(πfτ)2sin4(πfτ)df
上式中Sy(f)S_{y}(f)Sy(f)为功率谱密度。
从艾伦方差曲线中可以辨识出IMU的五种噪声，分别为：量化噪声、角度随机游走、零偏不稳定性噪声，角速率随机游走，速率斜坡，一般在IMU噪声辨识中用的比较多的是中间3种。
下面分别介绍5种噪声和Allan方差之间的关系。
1. 量化噪声（Quantization Noise，QN）
量化噪声产生的来源是传感器在经过A/D转换环节，即存在模拟量向数字量的转化的量化过程，然后这一过程是必然存在信息损失的，所以量化噪声代表了传感器检测的最小分辨率水平。
角速率的量化噪声功率谱为：
SΩ(f)=τ0Q2(2πf)2S_{\Omega}(f)=\tau_{0} Q^{2}(2 \pi f)^{2}SΩ(f)=τ0Q2(2πf)2
代入艾伦方差和功率谱密度关系公式中得到：
σQN2(τ)=4∫0∞τ0Q2(2πf)2sin4(πfτ)(πfτ)2df=16τ0Q2πτ3∫0∞sin4(πfτ)d(πfτ)\sigma_{Q N}^{2}(\tau)=4 \int_{0}^{\infty} \tau_{0} Q^{2}(2 \pi f)^{2} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} f=\frac{16 \tau_{0} Q^{2}}{\pi \tau^{3}} \int_{0}^{\infty} \sin ^{4}(\pi f \tau) \mathrm{d}(\pi f \tau)σQN2(τ)=4∫0∞τ0Q2(2πf)2(πfτ)2sin4(πfτ)df=πτ316τ0Q2∫0∞sin4(πfτ)d(πfτ)
上公式化简整理可得：
log10σQN(τ)=log103Q−log10τ\log _{10} \sigma_{Q N}(\tau)=\log _{10} \sqrt{3} Q-\log _{10} \taulog10σQN(τ)=log103Q−log10τ
由此可见，在Allan方差τ−σQN\tau-\sigma_{Q N}τ−σQN双对数图上，量化噪声对应的斜率为-1，它（或延长线）与τ=1\tau=1τ=1交点的纵坐标读数为3Q\sqrt{3} Q3Q，如下图所示。
这里需要注意一下，严格意义上应该称σ(τ)\sigma(\tau)σ(τ)为Allan标准差，而称σ(τ)2\sigma(\tau)^{2}σ(τ)2为Allan方差，但习惯上常常直接称σ(τ)\sigma(\tau)σ(τ)为Allan方差，以后在不引起歧义的情况下所说的Allan方差亦指σ(τ)\sigma(\tau)σ(τ)，在概念上不要造成混乱。
量化噪声具有很宽的带宽，属于高频噪声，在实际应用中可进行低通滤波处理或大部分被导航姿态更新（积分）环节所滤除，因而一般对系统精度的影响不大。
2．角度随机游走（Angular Random Walk，ARW）
如果陀螺仪角速度为高斯白噪声，那么积分得到的角度就会表现为角度随机游走现象。根据随机过程理论，随机游走是一种独立增量过程，含义是：角速率白噪声在两相邻采样时刻进行积分，不同时段的积分值之间互不相关。
从角速率来看，其功率谱为常值，公式如下：
SΩ(f)=N2S_{\Omega}(f)=N^{2}SΩ(f)=N2
公式中的N为角度随机游走系数，同样的，代入艾伦方差和功率谱密度关系公式中可得：
σARW2(τ)=4∫0∞N2sin4(πfτ)(πfτ)2df=4N2πτ∫0∞sin4(πfτ)(πfτ)2d(πfτ)=N2τ\sigma_{A R W}^{2}(\tau)=4 \int_{0}^{\infty} N^{2} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} f=\frac{4 N^{2}}{\pi \tau} \int_{0}^{\infty} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d}(\pi f \tau)=\frac{N^{2}}{\tau}σARW2(τ)=4∫0∞N2(πfτ)2sin4(πfτ)df=πτ4N2∫0∞(πfτ)2sin4(πfτ)d(πfτ)=τN2
可见，在τ−σARW(τ)\tau-\sigma_{A R W}(\tau)τ−σARW(τ)，角度随机游走的斜率为−1/2-1 / 2−1/2，它（或延长线）与 τ=1\tau=1τ=1的交点纵坐标读数即为角度随机游走系数NNN，如下图所示。
角度随机游走是陀螺仪非常重要的一个噪声参数，在卡尔曼滤波中设置P阵中需要辨识这一参数。如果使用的是陀螺仪角增量信号，因为增量输出是一个积分过程，所以Allan方差的角度随机游走系数和采样频率无关，但是如果直接采样角速率数据，建议尽量提高采样频率来减少信息的损失，比如取两倍的带宽频率，能达到4到6倍或者以上更佳。
3．零偏不稳定性噪声（Bias Instability，BI）
根据《光纤陀螺仪指标（国军标）》的定义，零偏不稳定性主要表征了光纤陀螺输出量围绕其均值变化的离散程度。其功率谱密度与频率成反比，零偏不稳定性噪声的角速率功率谱为：
SΩ(f)=B2/(2πf)S_{\Omega}(f)=B^{2} /(2 \pi f)SΩ(f)=B2/(2πf)
上式中BBB为零偏不稳定性系数，对应的方差计算公式如下：
σBI2(τ)=4∫0∞B22πfsin4(πfτ)(πfτ)2df=2B2π∫0∞sin4(πfτ)(πfτ)3d(πfτ)≈4B29\sigma_{B I}^{2}(\tau)=4 \int_{0}^{\infty} \frac{B^{2}}{2 \pi f} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} f=\frac{2 B^{2}}{\pi} \int_{0}^{\infty} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{3}} \mathrm{d}(\pi f \tau) \approx \frac{4 B^{2}}{9}σBI2(τ)=4∫0∞2πfB2(πfτ)2sin4(πfτ)df=π2B2∫0∞(πfτ)3sin4(πfτ)d(πfτ)≈94B2
所以，在τ−σBI(τ)\tau-\sigma_{B I}(\tau)τ−σBI(τ)双对数图上，零偏不稳定性的斜率为0，它（或延长线）与 τ=1\tau=1τ=1的交点纵坐标读数为2B/32 B / 32B/3，如下图所示。
零偏不稳定性噪声具有低频特性，在陀螺输出中表现为零偏随时间的缓慢波动。
4．角速率随机游走（Rate Random Walk，RRW）
和角速度随机游走的概念类似，如果陀螺角加速率建模为高斯白噪声，则角速率表现为随机游走现象。
角加速率的功率谱为：
SΩ′(f)=K′S_{\Omega^{\prime}}(f)=K^{\prime}SΩ′(f)=K′
根据功率谱相关性质，得到角速率的功率谱公式如下：
SΩ′(f)=K′S_{\Omega^{\prime}}(f)=K^{\prime}SΩ′(f)=K′
所以：
σRRW2(τ)=4∫0∞K2(2πf)2⋅sin4(πfτ)(πfτ)2df=K2τπ∫0∞sin4(πfτ)(πfτ)4d(πfτ)=K2τ3\sigma_{R R W}^{2}(\tau)=4 \int_{0}^{\infty} \frac{K^{2}}{(2 \pi f)^{2}} \cdot \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{2}} \mathrm{d} f=\frac{K^{2} \tau}{\pi} \int_{0}^{\infty} \frac{\sin ^{4}(\pi f \tau)}{(\pi f \tau)^{4}} \mathrm{d}(\pi f \tau)=\frac{K^{2} \tau}{3}σRRW2(τ)=4∫0∞(2πf)2K2⋅(πfτ)2sin4(πfτ)df=πK2τ∫0∞(πfτ)4sin4(πfτ)d(πfτ)=3K2τ
因此，在τ−σRRW(τ)\tau-\sigma_{R R W}(\tau)τ−σRRW(τ)双对数图上，角速率随机游走的斜率为1/21 / 21/2，它（或延长线）与 τ=1\tau=1τ=1交点的纵坐标读数为K/3K / \sqrt{3}K/3，如下图所示。
5．速率斜坡（Rate Ramp，RR）
如果陀螺仪的角速率输出随时间缓慢变化，比如由于环境温度变化，角速率Ω\OmegaΩ与测试时间ttt之间呈现线性关系，即：
Ω(t)=Ω(0)+Rt\Omega(t)=\Omega(0)+R tΩ(t)=Ω(0)+Rt
上公式中RRR为速率斜坡系数，根据艾伦方差的定义公式：
yi(τ)=[Ω(0)+Rti]+[Ω(0)+R(ti+τ)]2=Ω(0)+Rti+Rτ2y_{i}(\tau)=\frac{\left[\Omega(0)+R t_{i}\right]+\left[\Omega(0)+R\left(t_{i}+\tau\right)\right]}{2}=\Omega(0)+R t_{i}+\frac{R \tau}{2}yi(τ)=2[Ω(0)+Rti]+[Ω(0)+R(ti+τ)]=Ω(0)+Rti+2Rτ
yi+1(τ)=[Ω(0)+R(ti+τ)]+[Ω(0)+R(ti+2τ)]2=Ω(0)+Rti+3Rτ2y_{i+1}(\tau)=\frac{\left[\Omega(0)+R\left(t_{i}+\tau\right)\right]+\left[\Omega(0)+R\left(t_{i}+2 \tau\right)\right]}{2}=\Omega(0)+R t_{i}+\frac{3 R \tau}{2}yi+1(τ)=2[Ω(0)+R(ti+τ)]+[Ω(0)+R(ti+2τ)]=Ω(0)+Rti+23Rτ
所以有：σRR2(τ)=12⟨[yi+1(τ)−yi(τ)]2⟩=12⟨(Rτ)2⟩=R2τ22\sigma_{R R}^{2}(\tau)=\frac{1}{2}\left\langle\left[y_{i+1}(\tau)-y_{i}(\tau)\right]^{2}\right\rangle=\frac{1}{2}\left\langle(R \tau)^{2}\right\rangle=\frac{R^{2} \tau^{2}}{2}σRR2(τ)=21⟨[yi+1(τ)−yi(τ)]2⟩=21⟨(Rτ)2⟩=2R2τ2
可见，当角速率误差随时间线性变化时，将在τ−σRR(τ)\tau-\sigma_{R R}(\tau)τ−σRR(τ)双对数图上得到斜率为1的直线，它（或延长线）与τ=1\tau=1τ=1的交点纵坐标读数为R/2\mathrm{R} / \sqrt{2}R/2，如下图所示。
实际上，角速率斜坡更像是一种确定性的误差，而不是随机误差。角速率斜坡常常由系统误差引起的，比如环境温度的缓慢变化，通过严格的环境温度控制或者引入补偿机制可以降低此类误差。
5种艾伦方差汇总介绍与单位转换
上面介绍了5种噪声参数和Allan方差之间的关系，实际系统输出信号的功率谱不一定是单一斜率的，而可能由几种统计独立的线性功率谱叠加组合而成，若将该总功率谱作为Allan方差滤波器的输入，则陀螺误差的Allan方差分析结果可写为：
σA2(τ)=σQN2(τ)+σARW2(τ)+σBI2(τ)+σRRW2(τ)+σRR2(τ)\sigma_{A}^{2}(\tau)=\sigma_{Q N}^{2}(\tau)+\sigma_{A R W}^{2}(\tau)+\sigma_{B I}^{2}(\tau)+\sigma_{R R W}^{2}(\tau)+\sigma_{R R}^{2}(\tau)σA2(τ)=σQN2(τ)+σARW2(τ)+σBI2(τ)+σRRW2(τ)+σRR2(τ)
整理可得到：
σA2(τ)=σQN2(τ)+σARW2(τ)+σBI2(τ)+σRRW2(τ)+σRR2(τ)\sigma_{A}^{2}(\tau)=\sigma_{Q N}^{2}(\tau)+\sigma_{A R W}^{2}(\tau)+\sigma_{B I}^{2}(\tau)+\sigma_{R R W}^{2}(\tau)+\sigma_{R R}^{2}(\tau)σA2(τ)=σQN2(τ)+σARW2(τ)+σBI2(τ)+σRRW2(τ)+σRR2(τ)
需要注意的是上公式是在国际单位制上推导出来的，上式中的 σA(τ)\sigma_{A}(\tau)σA(τ)和τ\tauτ的单位分为是 rad/sr a d / srad/s和sss，由此可得到各种噪声系数的单位，即：
Q−rad,N−rad/s1/2,B−rad/s,K−rad/s3/2,R−rad/s2Q-\mathrm{rad}, \quad N-\mathrm{rad} / \mathrm{s}^{1 / 2}, \quad B-\mathrm{rad} / \mathrm{s}, \quad K-\mathrm{rad} / \mathrm{s}^{3 / 2}, \quad R-\mathrm{rad} / \mathrm{s}^{2}Q−rad,N−rad/s1/2,B−rad/s,K−rad/s3/2,R−rad/s2
但是我们习惯上σA(τ)\sigma_{A}(\tau)σA(τ)常以(∘)/h\left(^{\circ}\right) / \mathrm{h}(∘)/h为单位，并且并且各项误差系数也常使用°°°和hhh的组合作为单位，为了达到该目的，根据换算关系lrad/s=180/π1/3600(∘)/h\operatorname{lrad} / \mathrm{s}=\frac{180 / \pi}{1 / 3600}\left(^{\circ}\right) / \mathrm{h}lrad/s=1/3600180/π(∘)/h} ，进行如下转换：
σA2(τ)(180/π1/3600)2=3(Q×180/π×3600)2τ2+[N180/π(1/3600)1/2×60]2τ+4(B180/π1/3600)29+[K180/π(1/3600)3/2×160]2τ3+[R180/π(1/3600)2×13600]2τ22\begin{array}{l} \sigma_{A}^{2}(\tau)\left(\frac{180 / \pi}{1 / 3600}\right)^{2}=\frac{3(Q \times 180 / \pi \times 3600)^{2}}{\tau^{2}}+\frac{\left[N \frac{180 / \pi}{(1 / 3600)^{1 / 2}} \times 60\right]^{2}}{\tau} \\ +\frac{4\left(B \frac{180 / \pi}{1 / 3600}\right)^{2}}{9}+\frac{\left[K \frac{180 / \pi}{(1 / 3600)^{3 / 2}} \times \frac{1}{60}\right]^{2} \tau}{3}+\frac{\left[R \frac{180 / \pi}{(1 / 3600)^{2}} \times \frac{1}{3600}\right]^{2} \tau^{2}}{2} \end{array}σA2(τ)(1/3600180/π)2=τ23(Q×180/π×3600)2+τ[N(1/3600)1/2180/π×60]2+94(B1/3600180/π)2+3[K(1/3600)3/2180/π×601]2τ+2[R(1/3600)2180/π×36001]2τ2
通过相应的变量替换：
σA′(τ)=σA(τ)180/π1/3600((∘)/h),Q′=Q×180/π×3600(")N′=N180/π(1/3600)1/2((∘)/h1/2),B′=B180/π1/3600K′=K180/π(1/3600)3/2((∘)/h3/2)R′=R180/π(1/3600)2((∘)/h2)\begin{aligned} &\sigma_{A}^{\prime}(\tau)=\sigma_{A}(\tau) \frac{180 / \pi}{1 / 3600}\left(\left(^{\circ}\right) / \mathrm{h}\right), Q^{\prime}=Q \times 180 / \pi \times 3600(")\\ &N^{\prime}=N \frac{180 / \pi}{(1 / 3600)^{1 / 2}}\left(\left(^{\circ}\right) / \mathrm{h}^{1 / 2}\right), \quad B^{\prime}=B \frac{180 / \pi}{1 / 3600}\\ &K^{\prime}=K \frac{180 / \pi}{(1 / 3600)^{3 / 2}}\left(\left(^{\circ}\right) / \mathrm{h}^{3 / 2}\right) \quad R^{\prime}=R \frac{180 / \pi}{(1 / 3600)^{2}}\left(\left({ }^{\circ}\right) / \mathrm{h}^{2}\right) \end{aligned}σA′(τ)=σA(τ)1/3600180/π((∘)/h),Q′=Q×180/π×3600(")N′=N(1/3600)1/2180/π((∘)/h1/2),B′=B1/3600180/πK′=K(1/3600)3/2180/π((∘)/h3/2)R′=R(1/3600)2180/π((∘)/h2)
整理可得：
σA′2(τ)=3Q′2τ2+(60N′)2τ+4B′29+(K′/60)2τ3+(R′/3600)2τ22=∑k=−22Ak′τk\sigma_{A}^{\prime 2}(\tau)=\frac{3 Q^{\prime 2}}{\tau^{2}}+\frac{\left(60 N^{\prime}\right)^{2}}{\tau}+\frac{4 B^{\prime 2}}{9}+\frac{\left(K^{\prime} / 60\right)^{2} \tau}{3}+\frac{\left(R^{\prime} / 3600\right)^{2} \tau^{2}}{2}=\sum_{k=-2}^{2} A_{k}^{\prime} \tau^{k}σA′2(τ)=τ23Q′2+τ(60N′)2+94B′2+3(K′/60)2τ+2(R′/3600)2τ2=k=−2∑2Ak′τk
上式非常重要，有了它，当我们画出Allan方差曲线后，就可以直接通过读图的方式，辨识出各种噪声参数。
艾伦方差辨识参数的误差问题
理论上，Allan方差是随机过程一个现实的一种时间平均，但是由于实际应用时总是基于有限长度样本数据计算的，因而与经典方差估计一样，Allan方差也可看作是一种估计或统计量，服从一定的概率分布。研究表明，Allan方差估计σ^A(τ)\hat{\sigma}_{A}(\tau)σ^A(τ)的1σ1 \sigma1σ均方根误差可近似为（以相对百分比表示）
EA(τ)=∣σ^A(τ)−σA(τ)∣σA(τ)≈12(N/M−1)×100%E_{A}(\tau)=\frac{\left|\hat{\sigma}_{A}(\tau)-\sigma_{A}(\tau)\right|}{\sigma_{A}(\tau)} \approx \frac{1}{\sqrt{2(N / M-1)}} \times 100 \%EA(τ)=σA(τ)∣σ^A(τ)−σA(τ)∣≈2(N/M−1)1×100%
上式中 MMM为样本长度，NNN为分组数据长度，N/MN / MN/M为分组数，上式表明表明分组数越少Allan方差的估计误差越大。
艾伦方差曲线的绘制方法
假设我们获得了一组平均角速率样本序列：
Ωˉ1(τ0),Ωˉ2(τ0),Ωˉ3(τ0),⋯,ΩˉN(τ0)\bar{\Omega}_{1}\left(\tau_{0}\right), \bar{\Omega}_{2}\left(\tau_{0}\right), \bar{\Omega}_{3}\left(\tau_{0}\right), \cdots, \bar{\Omega}_{N}\left(\tau_{0}\right)Ωˉ1(τ0),Ωˉ2(τ0),Ωˉ3(τ0),⋯,ΩˉN(τ0)
公式中τ0\tau_{0}τ0为采样间隔，一般IMU的采样频率一般在几百HzHzHz甚至上千HzHzHz,一般在计算艾伦方差时会静止放置几个小时，但是由于样本序列的长度总是有限的，因此Allan方差计算只能给出理论值的一个估计。
下面给出实现Allan方差计算的具体步骤：
- 首先计算取样时间为τ0\tau_{0}τ0时的Allan方差：
σ^A2(τ0)=12(N−1)∑k=1N−1[Ωˉk+1(τ0)−Ωˉk(τ0)]2\hat{\sigma}_{A}^{2}\left(\tau_{0}\right)=\frac{1}{2(N-1)} \sum_{k=1}^{N-1}\left[\bar{\Omega}_{k+1}\left(\tau_{0}\right)-\bar{\Omega}_{k}\left(\tau_{0}\right)\right]^{2}σ^A2(τ0)=2(N−1)1k=1∑N−1[Ωˉk+1(τ0)−Ωˉk(τ0)]2 - 其次将取样时间间隔加倍，记τ1=21τ0\tau_{1}=2^{1} \tau_{0}τ1=21τ0和N1=[N/2]N_{1}=[N / 2]N1=[N/2]，（[⋅][\cdot][⋅]表示取整），在相继奇偶序号角速率之间作算术平均，即
Ωˉk(τ1)=Ωˉ2k−1(τ0)+Ωˉ2k(τ0)2k=1,2,3,⋯,N1\bar{\Omega}_{k}\left(\tau_{1}\right)=\frac{\bar{\Omega}_{2 k-1}\left(\tau_{0}\right)+\bar{\Omega}_{2 k}\left(\tau_{0}\right)}{2} \quad k=1,2,3, \cdots, N_{1}Ωˉk(τ1)=2Ωˉ2k−1(τ0)+Ωˉ2k(τ0)k=1,2,3,⋯,N1
组成新的取样时间间隔为τ1\tau_{1}τ1的平均角速率序列，即：
Ωˉ1(τ1),Ωˉ2(τ1),Ωˉ3(τ1),⋯,ΩˉN1(τ1)\bar{\Omega}_{1}\left(\tau_{1}\right), \bar{\Omega}_{2}\left(\tau_{1}\right), \bar{\Omega}_{3}\left(\tau_{1}\right), \cdots, \bar{\Omega}_{N_{1}}\left(\tau_{1}\right)Ωˉ1(τ1),Ωˉ2(τ1),Ωˉ3(τ1),⋯,ΩˉN1(τ1)
显然新序列的长度减半（但可能相差1个数据，下同），计算取样时间为τ1\tau_{1}τ1时的Allan方差，即：
σ^A2(τ1)=12(N1−1)∑k=1N1−1[Ωˉk+1(τ1)−Ωˉk(τ1)]2\hat{\sigma}_{A}^{2}\left(\tau_{1}\right)=\frac{1}{2\left(N_{1}-1\right)} \sum_{k=1}^{N_{1}-1}\left[\bar{\Omega}_{k+1}\left(\tau_{1}\right)-\bar{\Omega}_{k}\left(\tau_{1}\right)\right]^{2}σ^A2(τ1)=2(N1−1)1k=1∑N1−1[Ωˉk+1(τ1)−Ωˉk(τ1)]2 - 再次将取样时间间隔加倍，记τ2=2τ1=22τ0\tau_{2}=2 \tau_{1}=2^{2} \tau_{0}τ2=2τ1=22τ0和N2=[N1/2]N_{2}=\left[N_{1} / 2\right]N2=[N1/2] ,计算平均角速率，即：
Ωˉk(τ2)=Ωˉ2k−1(τ1)+Ωˉ2k(τ1)2k=1,2,3,⋯,N2\bar{\Omega}_{k}\left(\tau_{2}\right)=\frac{\bar{\Omega}_{2 k-1}\left(\tau_{1}\right)+\bar{\Omega}_{2 k}\left(\tau_{1}\right)}{2} \quad k=1,2,3, \cdots, N_{2}Ωˉk(τ2)=2Ωˉ2k−1(τ1)+Ωˉ2k(τ1)k=1,2,3,⋯,N2
组成新的取样时间间隔为τ2\tau_{2}τ2的平均角速率序列，即：
Ωˉ1(τ2),Ωˉ2(τ2),Ωˉ3(τ2),⋯,ΩˉN2(τ2)\bar{\Omega}_{1}\left(\tau_{2}\right), \bar{\Omega}_{2}\left(\tau_{2}\right), \bar{\Omega}_{3}\left(\tau_{2}\right), \cdots, \bar{\Omega}_{N_{2}}\left(\tau_{2}\right)Ωˉ1(τ2),Ωˉ2(τ2),Ωˉ3(τ2),⋯,ΩˉN2(τ2)
新序列的长度再次减半，计算取样时间为KaTeX parse error: Can't use function '$' in math mode at position 9: \tau_{2}$̲，时的Allan方差，即： \hat{\sigma}{A}^{2}\left(\tau{2}\right)=\frac{1}{2\left(N_{2}-1\right)} \sum_{k=1}{N_{2}-1}\left[\bar{\Omega}_{k+1}\left(\tau_{2}\right)-\bar{\Omega}_{k}\left(\tau_{2}\right)\right]{2}$$ - 如此反复将取样时间间隔不断加倍，记τL=2τL−1=2Lτ0\tau_{L}=2 \tau_{L-1}=2^{L} \tau_{0}τL=2τL−1=2Lτ0和NL=[NL−1/2]N_{L}=\left[N_{L-1} / 2\right]NL=[NL−1/2],直至最终序列的长度不小于2，得平均角速率序列为:
Ωˉ1(τL),⋯,ΩˉNL(τL)\bar{\Omega}_{1}\left(\tau_{L}\right), \cdots, \bar{\Omega}_{N_{L}}\left(\tau_{L}\right)Ωˉ1(τL),⋯,ΩˉNL(τL)
并计算取样时间为τL\tau_{L}τL时的Allan方差，即:
σ^A2(τL)=12(NL−1)∑k=1NL−1[Ωˉk−1(τL)−Ωˉk(τL)]2\hat{\sigma}_{A}^{2}\left(\tau_{L}\right)=\frac{1}{2\left(N_{L}-1\right)} \sum_{k=1}^{N_{L}-1}\left[\bar{\Omega}_{k-1}\left(\tau_{L}\right)-\bar{\Omega}_{k}\left(\tau_{L}\right)\right]^{2}σ^A2(τL)=2(NL−1)1k=1∑NL−1[Ωˉk−1(τL)−Ωˉk(τL)]2
至此获得一系列的点对，τi∼σ^A2(τi)\tau_{i} \sim \hat{\sigma}_{A}^{2}\left(\tau_{i}\right)τi∼σ^A2(τi)或τi∼σ^A(τi)\tau_{i} \sim \hat{\sigma}_{A}\left(\tau_{i}\right)τi∼σ^A(τi)，(i=1,2,3,⋯,L)(i=1,2,3, \cdots, L)(i=1,2,3,⋯,L)，完成Allan方差估计，并将结果绘制成双对数图。
从上述计算过程可以看出，在τ0\tau_{0}τ0基础上取样间隔时间是以2的倍数递增的，即取样时间为2的幂次方倍，而在一般应用中不需计算其他时间间隔上的Allan方差值。还容易看出，Allan方差的计算过程比较简单，并且计算量也不大。
艾伦方差相关代码
下面给出艾伦方差曲线绘制的matlab代码，可以对应上面的步骤。
function [sigma, tau, Err] = avar(y0, tau0)
% 计算Allan方差
% 输入：y -- 数据(一行或列向量), tau0 -- 采样周期
% 输出：sigma -- Allan方差(量纲单位与输入y保持一致), tau -- 取样时间, Err -- 百分比估计误差
% 作者: Yan Gong-min, 2012-08-22
% example:
% y = randn(100000,1) + 0.00001*[1:100000]';
% [sigma, tau, Err] = avar(y, 0.1);N = length(y0);y = y0; NL = N;for k = 1:infsigma(k,1) = sqrt(1/(2*(NL-1))*sum([y(2:NL)-y(1:NL-1)].^2));tau(k,1) = 2^(k-1)*tau0;Err(k,1) = 1/sqrt(2*(NL-1));NL = floor(NL/2);if NL<3break;endy = 1/2*(y(1:2:2*NL) + y(2:2:2*NL)); % 分组长度加倍(数据长度减半)endsubplot(211), plot(tau0*[1:N], y0); gridxlabel('\itt \rm/ s'); ylabel('\ity');subplot(212), loglog(tau, sigma, '-+', tau, [sigma.*(1+Err),sigma.*(1-Err)], 'r--'); gridxlabel('\itt \rm/ s'); ylabel('\it\sigma_A\rm( \tau )');
从艾伦方差曲线中辨识噪声参数可以通过手动读图的方式，当然也可以通过自动化的形式，自动辨识出参数，下面代码是一个自动辨识角度随机游走参数的代码，主要步骤是：通过辨识的噪声参数对应到log曲线频率，通过查找最为接近的频率对应的纵坐标数据。
需要注意的是下面代码输入的陀螺仪单位是国际单位，为rad/sr a d / srad/s，如果使用deg/h\operatorname{deg} / hdeg/h，需要注意单位的转换。
%% Noise Parameter Identification
% Find the index where the slope of the log-scaled Allan deviation is equal to the slope specified.% Angle Random Walk
% The above equation is a line with a slope of -1/2 when plotted on a log-log plot of σ(τ)versusτ.
% The value of N can be read directly off of this line at τ=1.
slope = -0.5;
logtau = log10(tau);
logadev = log10(adev);
dlogadev = diff(logadev) ./ diff(logtau);
[~, i] = min(abs(dlogadev - slope)); %找到slope=-0.5 最接近的点；% Find the y-intercept of the line.
b = logadev(i) - slope*logtau(i);
% Determine the angle random walk coefficient from the line.
logN = slope*log(1) + b;
% disp('random walks');
N =[10^logN]; %角度随机游走值，单位为：(rad/s/root-Hz)% Plot the results.
tauN = 1;
lineN = N ./ sqrt(tau);
figure
loglog(tau, adev, tau, lineN, '--', tauN, N, 'o')
title('Allan Deviation with Angle Random Walk')
xlabel('\tau')
ylabel('\sigma(\tau)')
legend('\sigma', '\sigma_N')
text(tauN, N, 'N')
grid on
axis equal
艾伦方差应用实例
1.ADIS16465
下图是是ADI的一款6轴MEMS IMU的数据手册，型号是ADIS16465，手册中给出的噪声参数如下，注意其单位：
- In-Run Bias Stability: 2.5∘/h2.5^{\circ} / h2.5∘/h
- Angular Random Walk: 0.15∘/h0.15^{\circ} / \sqrt{h}0.15∘/h
另外还可以看到其零偏重复性是0.4∘/sec0.4^{\circ} / \mathrm{sec}0.4∘/sec，，还是挺大的，所以一般需要在开始上电后静置数秒得到零偏重复性。根据国军标的定义，零偏重复性为：在同样条件下及规定间隔时间内，多次通电过程中，陀螺仪零偏相对其均值的离散程度，以多次测试所得零偏的标准偏差表示。
下图是ADIS16465静置2.5小时ZZZ轴的数据，采样率为200Hz200Hz200Hz，从图中可以看出：曲线大致可以分成三段，最左边的斜率为 −1/2-1/2−1/2，即角度随机游走噪声，因为这里的纵坐标单位是°/h°/h°/h，横坐标τ=1\tau=1τ=1对应的纵坐标为 7.2∘/h7.2^{\circ} / h7.2∘/h，根据前面介绍的如下公式：
σA′2(τ)=3Q′2τ2+(60N′)2τ+4B′29+(K′/60)2τ3+(R′/3600)2τ22=∑k=−22Ak′τk\sigma_{A}^{\prime 2}(\tau)=\frac{3 Q^{\prime 2}}{\tau^{2}}+\frac{\left(60 N^{\prime}\right)^{2}}{\tau}+\frac{4 B^{\prime 2}}{9}+\frac{\left(K^{\prime} / 60\right)^{2} \tau}{3}+\frac{\left(R^{\prime} / 3600\right)^{2} \tau^{2}}{2}=\sum_{k=-2}^{2} A_{k}^{\prime} \tau^{k}σA′2(τ)=τ23Q′2+τ(60N′)2+94B′2+3(K′/60)2τ+2(R′/3600)2τ2=k=−2∑2Ak′τk
从上式可知：
N=σARW⋅τ60=σARW13600h=7.2∘h60=0.12∘hN=\sigma_{A R W} \cdot \frac{\sqrt{\tau}}{60}=\sigma_{A R W} \sqrt{\frac{1}{3600} h}=\frac{\frac{7.2^{\circ}}{\sqrt{h}}}{60}=\frac{0.12^{\circ}}{\sqrt{h}}N=σARW⋅60τ=σARW36001h=60h7.2∘=h0.12∘
所以和数据手册中的Angular Random Walk相差不大。
下面看零偏不稳定性，斜率为0对应的横坐标为 41s41s41s，纵坐标为2.4°/h2.4°/h2.4°/h ，所以
B=2.4∘/(2/3)/h=3.6∘/hB=2.4^{\circ} /(2 / 3) / h=3.6^{\circ} / hB=2.4∘/(2/3)/h=3.6∘/h
可以看到和数据手册也相差不大，因为艾伦方差只能当成一种粗略的参数辨识手段，所以为了方便识图，这里当成2.4°/h2.4°/h2.4°/h 更简单点。
下面看角速率随机游走参数，这里看横坐标为100s100s100s时的数据，对应的纵坐标为3.0°/h3.0°/h3.0°/h，所以：
K=603σRRWτ=603⋅3.0100/h3/2=31∘/h3/2K=\frac{60 \sqrt{3} \sigma_{R R W}}{\sqrt{\tau}}=\frac{60 \sqrt{3} \cdot 3.0}{\sqrt{100}} / h^{3 / 2}=31^{\circ} / h^{3 / 2}K=τ603σRRW=100603⋅3.0/h3/2=31∘/h3/2
另外图中显示的数据平均值为329°/h329°/h329°/h，因为16465的零偏稳定性很小，只有3.0°/h3.0°/h3.0°/h，基本可以忽略不计，所以零偏重复性大概为0.09∘/sec0.09^{\circ} / \mathrm{sec}0.09∘/sec，小于数据手册中的0.4∘/sec0.4^{\circ} / \mathrm{sec}0.4∘/sec，但是仍然说明零偏重复性对系统影响巨大，在实际应用中需要很好的处理。
下图是ADIS16465的加速度计参数
下面是加速度三个轴的艾伦方差曲线，纵坐标单位为ggg，从图中可以看出，横坐标0.1到10之间的曲线的斜率大约为−1/2-1/2−1/2，所以该段曲线为速度随机游走误差， 横坐标10010^{0}100所对应的的纵坐标大约为20ug20 u g20ug，即0.012msh0.012 \frac{m}{s \sqrt{h}}0.012shm， 和数据手册给出的值相当。
横坐标为10 的地方的斜率大约为0，所以该点对应的是零偏不稳定性噪声，读图可知，零偏不稳定性噪声 X轴略大，为34ug34 u g34ug，Y轴和Z轴轴的数据大致相同，为20ug20 u g20ug，比数据手册中给出的略大。
1.ADIS16470
下图是ADIS16470静置2小时Z轴的采样的数据，采样率为200Hz200 H z200Hz，可以看到大致可以分成两段，角速率随机游走参数已经很小了。
数据手册参数如下：
- In-Run Bias Stability:8∘/h8^{\circ} / h8∘/h
- In-Run Bias Stability:0.34∘/h0.34^{\circ} / \sqrt{h}0.34∘/h
N=σARW⋅τ60=σARW13600h=12∘h60=0.2∘hB=5∘/(2/3)/h=7.5∘/h\begin{array}{c} N=\sigma_{A R W} \cdot \frac{\sqrt{\tau}}{60}=\sigma_{A R W} \sqrt{\frac{1}{3600} h}=\frac{\frac{12^{\circ}}{\sqrt{h}}}{60}=\frac{0.2^{\circ}}{\sqrt{h}} \\ B=5^{\circ} /(2 / 3) / h=7.5^{\circ} / h \end{array}N=σARW⋅60τ=σARW36001h=60h12∘=h0.2∘B=5∘/(2/3)/h=7.5∘/h
下图是ADIS16470的加速度计噪声参数。
下图是ADIS16470的加速度艾伦方差曲线，可以看到横坐标10之前的曲线的斜率大概为 −1/2-1 / 2−1/2，即为速度随机游走部分曲线，当横坐标为1时，纵坐标为60ug60 u g60ug，可得速度随机游走噪声为：0.036msh0.036 \frac{m}{s \sqrt{h}}0.036shm，和数据手册中的 噪声相当。
当横坐标为10时，斜率为0，对应的是零偏不稳定性，可得零偏不稳定性噪声大约为45ug45 u g45ug，比数据手册中的参数略大。
3.某型号光纤惯导
下图是光纤陀螺静置1.5个小时的采样数据，采样率为100Hz100 H z100Hz，可以看出曲线的后半段为角度随机游走噪声，零偏不稳定性在图中没有显示出来，可能采样时间还是短了些。
N=σARW⋅τ60=0.14⋅100∘60h=0.0023∘hN=\sigma_{A R W} \cdot \frac{\sqrt{\tau}}{60}=0.14 \cdot \frac{\sqrt{100^{\circ}}}{60 \sqrt{h}}=\frac{0.0023^{\circ}}{\sqrt{h}}N=σARW⋅60τ=0.14⋅60h100∘=h0.0023∘
另外从光纤陀螺的的数据均值可以看出，相比Mems陀螺仪的零偏重复性已经非常小了，图中的均值并不能反映零偏，因为没有去除地球自转角速度在Z轴的分量。
10s平滑计算零偏不稳定性
国军标中常用10s平滑方法计算零偏不稳定性，这里比较下它和艾伦方差在计算零偏稳定性之间的数值差异。10s平滑主要思路是每10s计算下均值，再计算这些均值的方差，例如有1000s的数据，需要计算100个均值，并求这100个均值的方差。主要参考代码如下：
for i=1:3%Bias stability and BiasSmo=10;%N 10秒平滑num=Smo/t0;%%m每组内数据个数m=floor(length(Wb(:,i))/num);%%共可以分成m组数据gx=zeros(m,1);Xbais_all = mean(Wb(:,i));for j=1:m-1gx(j)=mean(Wb(1+(j-1)*num:1+j*num,i)); %取平均endgx= gx(2:end-2,:);Xbais_stability = std(gx); % 计算标准差fprintf('Bias Stability 10s:%0.2e [deg/h]\n',Xbais_stability);
end
- 光纤陀螺的10s平滑结果为：4.62∘/h4.62^{\circ} / h4.62∘/h
- ADIS16470的10s平滑结果为：13.3∘/h13.3^{\circ} / h13.3∘/h；艾伦方差：7.5∘/h7.5^{\circ} / h7.5∘/h
- ADIS16465的10s平滑结果为：29.1∘/h29.1^{\circ} / h29.1∘/h；艾伦方差：3.6∘/h3.6^{\circ} / h3.6∘/h
可看到10平滑方法计算出来的零偏稳定性相对艾伦方差较大，不知道有没有理论证明。
下面是加速度计的10s平滑结果：
- ADIS16465：XXX轴99ug99 u g99ug，YYY轴26ug26u g26ug，ZZZ轴29ug29 u g29ug;
- ADIS16470：XXX轴4.24mg4.24 m g4.24mg，YYY轴6.17mg6.17m g6.17mg，ZZZ轴0.46mg0.46 m g0.46mg;
- 光纤惯导：XXX轴97ug97 u g97ug，YYY轴53ug53u g53ug，ZZZ轴150ug150 u g150ug;
参考文献
- 惯性仪器测试与数据分析 严恭敏
- Analysis and Modeling of Inertial Sensors Using Allan Variance
IMU噪声参数辨识-艾伦方差相关推荐
- 【学习记录】IMU内参标定：Allan方差与代码
本文仅用于记录自己学习IMU内参标定过程中的一些总结. 参考 关于IMU参数: 死磕陀螺仪之(一)陀螺仪参数意义以及工程转换 关于Allan方差: 多传感器融合定位理论基础(三):惯性器件误差分析 I ...
- 基于Matlab使用艾伦方差来确定MEMS陀螺仪的噪声参数（附源码）
目录 一.背景 二.艾伦方差计算 二.噪声参数识别 2.1 角度随机游走 三.速率随机游走 四.偏置不稳定性 五. 陀螺仪模拟 六.程序 本示例展示了如何使用艾伦方差来确定MEMS陀螺仪的噪声参数.这 ...
- 在Hammerstein非线性模型中，基于PSO的参数辨识系统
Hammerstein非线性模型的基于PSO的参数辨识系统的本质就是将参数的辨识问题转换为参数空间优化问题,对整个参数域进行搜索并最终获得最优的参数估计.我们需要的参数辨识模型具体描述如下所示: 将H ...
- IMU特性参数、误差模型及卡尔曼滤波参数设置
文章目录 IMU指标参数与卡尔曼滤波参数设置 两种IMU误差模型 PSINS代码中的IMU误差生成 参考文献 本文主要总结了IMU的两种常用误差模型,以及如何通过IMU的指标对卡尔曼滤波中的参数进行设 ...
- matlab传递函数参数辨识,基于matlab/Simulink的参数辨识
基于Simulink的辨识 Simulink自带Parameter Estimation功能.可以对Simulink模型中的参数进行估计.MATLAB的Parameter Estimation官方说明 ...
- 第5章 Python 数字图像处理(DIP) - 图像复原与重建8 - 估计噪声参数
标题 估计噪声参数 估计噪声参数 周期噪声的参数通常是通过检测图像的傅里叶谱来估计的. 只能使用由传感器生成的图像时,可由一小片恒定的背景灰度来估计PDF的参数. 来自图像条带的数据的最简单用途是,计 ...
- 【参数辨识】永磁同步电机的参数辨识
辨识的参数为: 定子电阻Rs, 永磁磁链φf, dq轴电感Lq 主要步骤可如下: 1.建立矢量控制的电机模型 2.确定辨识方法 3.将定子电压方程转化为辨识方程 4.搭建辨识仿真模型,开始仿真 第一步 ...
- 锂离子电池离线参数辨识（基于二阶RC电池模型）
1基本参数概念及特性分析 1.1电池电压主要包括:电动势.开路电压和端电压等. (1)电动势(或理论电压),是指电池内部的电化学反应完全停止,电池达到稳定状态时,两极间的电位差.电动势是以电池内部化学 ...
- 电机仿真系列-基于最小二乘法的永磁同步电机参数辨识
基于最小二乘法的永磁同步电机参数辨识 永磁同步电机具有体积小.转动惯量低.结构简单等优点,被广泛应用于控制系统中.然而在实际应用过程中,控制系统会受到高温.负载等外界因素的影响,永磁同步电机的电感 ...
- 永磁同步电机在线参数辨识综述
优点 实时性:在线参数辨识能够在电机实际运行时进行参数估计,可以实时地获取电机的参数信息.这使得在线参数辨识更适用于需要实时调节和优化电机控制策略的应用场景. 动态适应性:在线参数辨识可以根据电机的实 ...
最新文章
- latex中常使用符号大全
- 关于我喜欢计算机的作文600字,关于我喜欢的字作文600字5篇
- 具体knn算法概念参考knn代码python实现
- 怎样自己写一个MVC框架
- apache+webdav的安装配置
- ios 替换数组中元素_leetcode169 数组中的主要元素
- 二叉树的前序中序后序 递归与非递归解法
- 华硕笔记本电脑＜举例：华硕玩家国度G752VS＞启动时/重装系统后开机会自动进入BIOS界面？
- 计算机作业动画flash,计算机动画设计:Flash
- dbv连接mysql_mysql数据库版本控制dbv使用_MySQL
- 51单片机波特率计算c语言,8051单片机波特率计算公式(配套C语言例程
- python 输入一个整数，将该整数逆向输出
- MIS课设 JavaFX考试管理系统
- iconv()和mb_conver_encoding()字符编码转换函数
- Mac系统-重置MySQL登陆密码
- vue脚手架项目使用element-ui
- 频繁模式挖掘Frequent Pattern Mining
- 一个详细测试用例的文档
- 联想Thinkpad重装系统的详细操作指南
- Proxy error: Could not proxy request /students from localhost:8080 to http://localhost:5000/.See ht