Allan 方差分析以及python代码实现

阿伦（Allan）方差分析广泛应用于惯性器件的随机误差建模，常见的有陀螺误差分析。我需要对光纤陀螺的随机误差进行分析，故对阿伦方差的的原理以及python代码实现进行了整理，其中理论参考严恭敏老师的《惯性仪器测试与数据分析》第八章，python代码参考Allen B. Downey的《Think DSP》¹

统计学中用均值和方差描述随机变量，但是，Allan 在1966年分析铯原子频标的频率稳定度时发现经典方差随着测量时间的增长而发散，为了解决该问题，提出了阿伦方差(Allan Variance, AVAR)分析法，后来IEEE在1971年正式推荐Allan方差作为频率频率稳定度的时域分析方法，后来被广泛应用于惯性器件的分析中。

随机信号Allan方差的物理意义及应用本质上源于它和功率谱之间的关系。故本文分为三个部分，1.功率谱的幂律模型以及python代码实现 2. 经典方差和Allan方差 3. 陀螺随机漂移误差的Allan分析以及python代码实现。

功率谱的幂律模型以及python代码实现

数学定义的连续时间白噪声$w(t)$，其在整个频率轴$f\in (-\infty ,\infty )$上功率谱都为常值$S_{wt}$ ，并且仅在$t=0$时自相关函数取值无限大，其他任何间隔时间下均不相关。在实际系统中，如果处理信号的带宽远远小于噪声源的带宽，则可以对带宽噪声源的功率谱进行合理扩展，以简化系统分析，通常数学上定义的白噪声是实际物理系统在感兴趣频段上平直功率谱外推后的理想化建模。

对连续时间信号$x(t)$及其傅里叶变换$X(jf)$,根据时域微分或者积分性质有

$$CTFT[x\prime (t)]=j2\pi f\cdot X(jf)$$

$CTFT[\underset{-\infty }{\overset{t}{\int }}x(\tau )d\tau ]=\pi X(0)\delta (f)+\frac{1}{j2\pi f}X(jf)=\frac{1}{j2\pi f}X(jf)(X(0)=0)$

由此可得，积分和微分运算均为线性变换运算，可将运算结果看作是原信号经过线性系统之后的输出。如果将白噪声$w(t)$输入微分线性系统$H(jf)=j2\pi f$(在数学上白噪声不可微，这里只是形式上)则系统输出随机微分和积分信号的功率谱为：

$S_{wt\prime }=(2\pi f{)}^2{S}_{wt}$

$S_{\int wt}(f)=(2\pi f{)}^{-2}{S}_{wt}$

容易看出，微分系统具有高通滤波器特性，白噪声随机微分后，高频分量增强。积分系统具有低通滤波特性，低频分量增强。将白噪声进行一次随机微分和积分，分别得到了与$f^2$和$f^{-2}$成比例的功率谱，如果对白噪声进行$k$次随机积分或者微分，所得功率谱将具有幂次$S_k=A_k{F}^{2k}(k=0,\pm 1,\pm 2...)$形式，如果绘制log-log图，则有$lo{g}_{10}{S}_k=lo{g}_{10}{A}_k+2klo{g}_{10}f$ ,其为斜率是$2k$的线性功率谱。在实际的物理系统中，除了斜率为$2k$的功率谱外，在大量电子设备中普遍存在着$1/f$噪声。下图给出了几种斜率的功率谱，其中$f^0(beta=0)$对应白噪声序列，$f^{-2}(beta=-2)$对应白噪声序列的累计(即随机游走序列)，$f^2(beta=2)$对应白噪声序列的差分，$f^{-1}(beta=-1)$和$f^1(beta=1)$通过对频域直接处理得到。

下面对python代码实现进行说明,需要将thinkdsp.py文件放在运行目录中，文件获取后台回复阿伦方差

import numpy as np
import matplotlib.pyplot as plt

from thinkdsp import decorate #导入作图帮助函数

from thinkdsp import BrownianNoise  # 引入随机游走噪声序列

signal = BrownianNoise() 
wave = signal.make_wave(duration=0.5, framerate=11025) # 将随机游走生成波
wave.make_audio() # 在jupyter中会生成音频文件，可以试听

# 可以在时域看看其序列
wave.plot(linewidth = 1)
decorate(xlabel='Time(s)',
         ylabel='Amplitude')

# 对随机游走序列进行频域分析，并且作图
spectrum = wave.make_spectrum()
spectrum.plot_power(linewidth=0.5)
decorate(xlabel='Frequency (Hz)',
         ylabel='Power')

# 由于其低频分量占比较多，故采用loglog图
spectrum.hs[0] = 0 # log 0 不存在，故修改
spectrum.plot_power(linewidth=0.5)

loglog = dict(xscale='log', yscale='log')
decorate(xlabel='Frequency (Hz)',
         ylabel='Power',
         **loglog)spectrum.hs[0] = 0
spectrum.plot_power(linewidth=0.5)

loglog = dict(xscale='log', yscale='log')
decorate(xlabel='Frequency (Hz)',
         ylabel='Power',
         **loglog)

result = spectrum.estimate_slope() # 对功率谱斜率估计
[out]: -1.7954213807211257 #接近-2

其他斜率的功率生成与上述随机游走类似，故不在赘述。需要注意的是

1. BlueNoise和 PinkNoise分别表示斜率区间为(0,2)和(-2,2),其中BlueNoise是我自己改的，如果引用作者源代码的化可能会报错
2. 由于代码是参考的音乐的处理方式，上述的斜率和和参数beta是相反的，如果不确定对应关系，可以通过斜率估算函数得到对应的斜率。
3. 建议对源码进行阅读。

如果若干个随机过程$X_a$之间相互独立且设总和$X=\Sigma X_a$,则有$PSD[X]=\Sigma PSD[X_a]$ 。若$X_a$为线性功率谱，即$PSD[X_a]=h_a{f}^a$,并且实际应用中$a$取$-2,-1,0,1,2$,则随机过程$f-S$为

$lo{g}_{10}S=lo{g}_{10}(\Sigma h_a{f}^a)\approx \Sigma (lo{g}_{10}{h}_a+alo{g}_{10}f),f_{aL}<f<f_{aH}$

对理想的幂律谱噪声按不同比例进行混合,得到的s-log可以大概看到不同斜率的占比。

from thinkdsp import Wave

def get_sum_wave(K, wave_red_noise=wave_red_noise, wave_pink_noise=wave_pink_noise, wave_white_noise=wave_white_noise, wave_blue_noise=wave_blue_noise, wave_purple_noise=wave_purple_noise):
    wave_sum_ys = wave_red_noise.ys*K[0] + wave_pink_noise.ys*K[1] + wave_white_noise.ys*K[2] + wave_blue_noise.ys*K[3] + wave_purple_noise.ys*K[4]
    wave_sum = Wave(wave_sum_ys, wave_white_noise.ts, wave_white_noise.framerate)
    return wave_sum

K1 = [10,0,10,0,10]
K2 = [10,0,0,0,10]

spectrum_sum1 = get_sum_wave(K1).make_spectrum()
spectrum_sum2 = get_sum_wave(K2).make_spectrum()
spectrum_sum1.plot_power(linewidth=0.5, label='K1')
spectrum_sum2.plot_power(linewidth=0.5, label='K2')
loglog = dict(xscale='log', yscale='log')
decorate(xlabel='Frequency (Hz)',
         ylabel='Power',
         **loglog,
         ylim=(1e-1,1e8))

结果如下图所示

K1相比于K2在于有更多的白噪声分量，表现出结果就是在倒角处的斜率更加接近平行（白噪声斜率为0）

经典方差和Allan方差

频标输出的准周期信号通用模型为：
$$u(t)=[A_0+A(t)]sin[2\pi v_{0}t+\varphi (t)]$$
式中$A_0,v_0$分别为信号输出的理想幅值和标称频率，$A(t),\varphi (t)$分别为幅值起伏的相位抖动，一般情况下忽略$A(t)\ll A_0$,分析时忽略$A(t)$的影响，信号的瞬时频率为相位相对时间的导数，即
$$v(t)=\frac{1}{2\pi }\cdot \frac{d[2\pi v_{0}t+\varphi (t)]}{dt}=v_0+\frac{1}{2\pi }\cdot \frac{d\varphi (t)}{dt}$$
定义瞬时相对频率偏差：
$$y(t)=\frac{v(t)-v_0}{v_0}=\frac{1}{2\pi v_0}\cdot \frac{d\varphi (t)}{dt}$$
在实际情况中，通常无法测得频率的瞬间偏差值，而只能测一段时间$[t_i,t_i+\tau ]$内的平均值 $y_i(\tau )$,即
$$y_i(\tau )=\frac{1}{\tau }\underset{t_i}{\overset{t_i+\tau }{\int }}y(t)dt$$
其中$\tau$为平均时间，也就是平均频率的取样时间，经典方差也就是$y_i^2$,当无间隔取样时有
$$<y_i^2(\tau )>=2\underset{0}{\overset{\infty }{\int }}{S}_y(f)\frac{si{n}^2(\pi f\tau )}{(\pi f\tau {)}^2}df$$
其中$S_y(f)$为功率谱函数，具体推导过程参考严老师《惯性仪器测试与数据分析》

由上分析可得，平均时间为$\tau$的时域经典方差$y_i^2(\tau )$表示频率偏差信号$y(t)$通过滤波器$H(jf)$后的平均功率输出，其中$H(jf)$为

$$H(jf)=\frac{sin(\pi f\tau )}{\pi f\tau }=Sa(\pi f\tau )$$

---

若$y(t)$非平稳，如果含有趋势项（均值时变但方差定常）或者表现为维纳过程，这是$y_i(\tau )$也将包含趋势项或表现为随机游走过程，对于非平稳信号的功率谱$S_y(f)$在零频率附近为无穷大，且不可积，这样导致经典方差$y_i^2(\tau )$也趋于无穷，因此对于非平稳频标，不能再使用经典方法来表示频率稳定性。

现定义Allan方差如下：
$${\sigma }_y^2(\tau )=\frac{1}{2}<[y_{i+1}(\tau )-y_i(\tau ){]}^2>$$
略去推理过程可以得到
$${\sigma }_y^2(\tau )=4\underset{0}{\overset{\infty }{\int }}{S}_y(f)\frac{si{n}^4(\pi f\tau )}{(\pi f\tau {)}^2}df$$
同理，采样时间为$\tau$​的时域阿伦方差${\sigma }_i^2(\tau )$​表示频率偏差信号$y(t)$​通过滤波器$H_A(jf)$​后的平均功率输出，其中$H_A(jf)$​为
$$H_A(jf)=-j\sqrt{2}\frac{si{n}^2(\pi f\tau )}{\pi f\tau }=Sa(\pi f\tau )\times [-j\sqrt{2}sin(\pi f\tau )]$$
其中$Sa(\pi f\tau )$为低通滤波器，对应均值采样过程，$[-j\sqrt{2}sin(\pi f\tau )]$为带通滤波器，对应差分$y_{i+1}(\tau )-y_i(\tau )$过程。所以Allan方差计算的实质是使信号通过一个带通滤波器，然后做某种时间统计平均$[\cdot {]}^2$

通过与上一部分的幂律谱函数结合，则可以得到Allan方差是一种对功率谱进行加权后的自相关，加权函数为$|H_A(jf){|}^2$。

陀螺随机漂移误差的Allan分析以及python代码实现

综合上述的分析，可以得到，如果要对陀螺随机漂移误差进行阿伦方差建模，则首先需要得到影响陀螺随机漂移噪声源的功率谱，带入上节公式(8)得到对应噪声的阿伦方差，最后根据幂律谱模型进行建模。

本节对影响陀螺的噪声源物理背景进行说明，最后将阿伦方差整理为表格，略去推理过程，推理过程参见教材。

1. 量化噪声（Quantization Noise，QN）

   量化噪声是一切量化操作所固有的噪声，只要进行数字量化编码采样，传感器输出的理想值与量化值之间就必然会存在微笑差别，量化噪声代表了传感器检测的最小分辨率。

2. 角度随机游走（Angular Random Walk， ARW）

   角度随机游走是宽带角速率白噪声积分的结果，即陀螺从零时刻其累积的总角增量误差表现为随机游走，而每一时刻的等效角速率误差表现为白噪声。

3. 角速率随机游走（Rate Random Walk， RRW）

   角速率随机游走是宽带叫角加速率白噪声积分的结果，即陀螺角加速率误差表现为白噪声，而角速率误差表现为随机游走。

4. 零偏不稳定性噪声（Bias Instability， BI）

   零偏不稳定性噪声又称为闪变噪声或$1/f$噪声，顾名思义，其功率谱密度与频率成反比。

5. 速率斜坡（Rate Ramp, RR）

   陀螺角速率输出随时间缓慢变化，对于光纤陀螺来说，环境温度变化会引起角速率的线性变化。

噪声源的阿伦方差分析如下

噪声类型	功率谱$S_{\Omega }(f)$	Allan方差${\sigma }_A^2(\tau )$	与$\tau =1$交点纵坐标	斜率
量化噪声	${\tau }_0{Q}^2(2\pi f{)}^2$	$3Q^2/{\tau }^2$	$\sqrt{3}Q$	$-1$
f噪声	$F(2\pi f)$	$\frac{3F\ln (e\pi f_H\tau )}{\pi {\tau }^2}$	$\sqrt{3F\ln (e\pi f_H)/\pi }$
角度随机游走	$N^2$	$N^2/\tau$	$N$	$-1/2$
零偏不稳定性	$B^2/(2\pi f)$	$4B^2/9$	$2B/3$	$0$
角速率随机游走	$K^2/(2\pi f{)}^2$	$K^2\tau /3$	$K/\sqrt{3}$	$1/2$
速率斜坡	$---$	$R^2{\tau }^2/2$	$R/\sqrt{2}$	$1$

根据上表以及幂律谱模型，可以得到陀螺的Allan 方差模型为
$${\sigma }_A^2(\tau )={\sigma }_{QN}^2(\tau )+{\sigma }_{ARW}^2(\tau )+{\sigma }_{BI}^2(\tau )+{\sigma }_{RRW}^2(\tau )+{\sigma }_{RR}^2(\tau )$$
将上表的第三列带入得：
$${\sigma }_A^2(\tau )=\frac{3Q^2}{{\tau }^2}+\frac{N^2}{\tau }+\frac{4B^2}{9}+\frac{K^2\tau }{3}+\frac{R^2{\tau }^2}{2}={\Sigma }_{k=-2}^2{A}_k{\tau }^k$$

---

对实际数据进行阿伦方差，需要获得一组采样间隔为${\tau }_0$内的平均角速率值，通常取${\tau }_0=1$。

python代码如下

def get_allen(y, tau0=1):
    N = len(y)
    NL = N 
    Tau  = [] # 保存不同的tau
    Sigma = [] # 保存不同tao下的阿伦方差值
    Err = []
    for k in np.arange(1, 1000):
        sigma_k = np.sqrt(1/(2*(NL-1)) * np.sum(np.power(y[1:NL]-y[0:NL-1], 2))) #Allan的时域表达式
        Sigma.append(sigma_k)
        tau_k = 2 ** (k-1) * tau0 #将取样时间加倍， tau2 = 2tau1
        Tau.append(tau_k)
        err = 1 / np.sqrt(2* (NL-1))
        Err.append(err)
        NL = np.floor(NL/2) 
        NL = int(NL) 
        if NL < 3: 
            break
        y = 1/2 * (y[0:2*NL:2] + y[1:2*NL:2]) # 对应的序列长度减半
    return  Sigma, Tau

对随机数进行Allan分析得到如下结果

if __name__ == '__main__':

    data_size_ex = 5
    y = np.random.standard_normal((10**data_size_ex),) + (10 ** -data_size_ex) * np.arange(1,10**data_size_ex + 1)
    fig, ax = plt.subplots()
    Sigma, Tau = get_allen(y)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_adjustable("datalim")ax.plot(Tau, Sigma)
    ax.grid()
    plt.show()

---

1. Think DSP – Green Tea Press ↩︎