Allan 方差计算方法实现

计算Allan方差的Python 实现，两种分组方式（标准分组和交叠式分组）都有考量，并且还提供了简单的调用示例，很容易上手。

流浪猪头拯救地球

4471人浏览 · 2022-10-14 11:05:53

流浪猪头拯救地球 · 2022-10-14 11:05:53 发布

文章目录

Part.I 原理
- Chap.I 计算过程
- Chap.II 注意事项
Part.II 代码实现
- Chap.I 实现函数
- Chap.II 调用示例
Reference

Part.I 原理

Chap.I 计算过程

具体原理可看 Reference 的前两个，下面仅仅列举一下计算过程：

步骤：分块，块内求平均、平均值相减、求差值的均方值（RMS）
根据分块方式不同可以分为标准分块方式和交叠式分块方式，其计算方法分别如下：

在这里插入图片描述
由随机过程理论可知，各种随机噪声都是由白噪声激励得到，一个谱密度为 $q_c$ 的白噪声输入不同的系统函数，就会得到不同的随机噪声。Allan 方差分析的目的就是要识别随机噪声类型以及提取相应的白噪声谱密度 $q_c$ 。不同性质的噪声具有不同的 Allan 方差，其双对数曲线具有不同的斜率。
在这里插入图片描述
其中

斜率 -1：Quantization noise 量化噪声
斜率 -0.5：Angle Random Walk 角度随机游走
斜率 0：Bias Instability 零偏不稳定性
斜率 0.5：Rate Random Walk 角速率随机游走
斜率 1：Rate Ramp 速率斜坡

对上述几种阿伦方差分析结果进行比较整理如下表所示。
在这里插入图片描述

Chap.II 注意事项

摘自参考 3

Allan方差分析主要用于分析陀螺量化噪声、角度随机游走噪声（角速率白噪声）、零偏不稳定性、角速率随机游走（角加速度白噪声）和速率斜坡（角速率趋势项），等五种典型误差。有些文献也提到可用于分析周期项和一阶马尔可夫误差，但这不是Allan方差分析的强项，周期项更适合用功率谱分析，马尔可夫更适合于用相关分析（时间序列分析），不要试图以己之短攻人所长！
Allan方差分析的是静态误差，陀螺必须在静基座下进行采集数据，否则要是动基座下采集的陀螺数据，分析的结果到底是陀螺误差还是基座运动特性？
Allan方差计算中采用了“差分”计算操作，因此Allan方差无法分析陀螺随机常值漂移参数，若是将几个陀螺的Allan方差曲线同时绘制在一个双对数图上，曲线在左下角的对应的陀螺性能一般较好些，而右上角的差些。
Allan方差可用于分析五种典型误差，但是某个陀螺中并不一定五种误差都有表现，不画Allan方差图/不看图而盲目的采用最小二乘回归方法强行地就想计算五个系数的人是典型的教条主义者，往往得不到好的结果。
如果画出的Allan方差图形不像典型的U型或V型图，不要感觉到诧异，也不要再试图用 Allan 方差去分析它（否则又是教条主义），没意义！此时Allan 分析失效，应试图从原始数据进行分析和解释。
Allan 方差分析的一个用途是分析陀螺的性能或对比不同陀螺的性能，相比于其它分析法Allan法还是很好用的，比较全面。另一个用途是获得噪声参数，用于组合导航的 Kalman 滤波噪声参数设置。不是所有的 Allan 方差噪声系数都有用，主要有用的是角度随机游走系数（用于设置 Q 阵）和零偏不稳定性系数（用于设置一阶马氏过程的方差），其实这两个系数量级大小差不多就行了，太精细也没用，毕竟 Allan 方差分析得出的陀螺静态性能，鬼知道实际应用动态的时候陀螺误差会变化多大，存在数量级差别都很有可能。

Part.II 代码实现

下面使用 Python 对其计算方式进行了实现，函数如下

Chap.I 实现函数

def allan_cross(data,m):
    """
    Calculate Allan standard deviation, Overlapping grouping.
    
    > @param[in] data:          the data list, one line.
    > @param[in] m:             the length of a block.
    return: 
    < @param[out] sigma:        the allan standard deviation
    """
    N=len(data)
    m=math.floor(m)
    ylist=[]
    for i in range(N-m):
        j=i+m
        ave=np.mean(data[i:j])
        ylist.append(ave)
    sigma_sum=0
    for i in range(N-2*m):
        y1=ylist[i+m]
        y2=ylist[i]
        sigma_sum+=(y1-y2)**2
    sigma2=sigma_sum/(N-2*m+1)/2
    return math.sqrt(sigma2)


def allan_standard(data,m):
    """
    Calculate Allan standard deviation, standard grouping.
    
    > @param[in] data:          the data list, one line.
    > @param[in] m:             the length of a block.
    return: 
    < @param[out] sigma:        the allan standard deviation
    """
    N=len(data)
    m=math.floor(m)
    ylist=[]
    for i in range(0,N,m):
        j=i+m
        ave=np.mean(data[i:j])
        ylist.append(ave)
    sigma_sum=0
    n=len(ylist)
    for i in range(n-1):
        y1=ylist[i]
        y2=ylist[i+1]
        sigma_sum+=(y1-y2)**2
    sigma2=sigma_sum/(n-1)/2
    return math.sqrt(sigma2)

def allan_curve(Data,mode=0):
    """
    Get the Allan standard deviation curve.
    
    > @param[in] Data:          [time,data] the raw data
    > @param[in] mode:          0:standard; 1:cross
    return: 
    < @param[out] [Tlist,allanList]:    the result for plotting
    """
    x=Data[0]
    dx=x[1]-x[0]
    data=Data[1]
    N=len(data)
    mi=math.log10(N)
    allanList=[]
    Tlist=[]
    for i in range(math.floor(10*mi)):
        m=10**(0.1*i)
        if mode==0:
            allanVal=allan_standard(data,m)
        else:
            allanVal=allan_cross(data,m)
        Tlist.append(dx*math.floor(m))
        allanList.append(allanVal)
    return [Tlist,allanList]

Chap.II 调用示例

调用示例如下：

def allan_proc():
    """
    The main function to call allan function.
    """
    sdir=r'C:\Users\OHanlon\Desktop\temp'
    f1=sdir+'\\'+'allan.png'
    N=1000
    x=list(range(1,N+1))
    data=[random.random() for i in range(N)]
    Data=[x,data]
    [x,y]=allan_curve(Data,1)
    plt.loglog(x,y)
    plt.savefig(f1,dpi=300)
    return 0