Allan方差方法分析环形激光陀螺仪噪声的性能评估
黎奇1, 白征东1, 赵思浩2, 戴东凯3, 邢海峰4     
1. 清华大学 土木工程系, 北京 100084;
2. 清华大学 电子工程系, 北京 100084;
3. 国防科技大学 前沿交叉学科学院, 长沙 410073;
4. 清华大学 精密仪器系, 北京 100084
摘要:为验证Allan方差方法用于不同状态下不同类型陀螺仪噪声分析的适用性,该文系统总结了Allan方差方法用于分析室温(~25℃)下静态环形激光陀螺仪噪声的方法,对该陀螺仪常见的5个主要噪声项以外的其他噪声项作了补充说明。对中国自主生产的某型号环形激光陀螺仪和目前使用较广泛的MPU 9250微系统惯性测量单元中的陀螺仪分别进行静态和动态下的多次实验与对比分析,结果表明:在对陀螺仪噪声特性充分了解的前提下,Allan方差方法可以对静态下不同类型陀螺仪的若干主要噪声项进行建模与概算,以满足全球导航卫星系统/惯性导航系统(GNSS/INS)的组合系统中Kalman滤波器参数设置的需要。若干实验分析结论对《IEEE单轴干涉式光纤陀螺仪试验程序和格式指南的标准规范》中经典Allan方差方法作了补充解释。通过对比Allan方差及相关常用数据分析方法,概述了其一般适用性。Allan方差方法可以为改进精密仪器的设计与制造、提高惯性产品实际使用精度等提供依据。
关键词环形激光陀螺仪    微系统惯性测量单元    MPU 9250    Allan方差(AVAR)    评估    
Performance evaluation of the Allan variance method for ring laser gyroscope noise analyses
LI Qi1, BAI Zhengdong1, ZHAO Sihao2, DAI Dongkai3, XING Haifeng4     
1. Department of Civil Engineering, Tsinghua University, Beijing 100084, China;
2. Department of Electronic Engineering, Tsinghua University, Beijing 100084, China;
3. College of Advanced Interdisciplinary Studies, National University of Defense Technology, Changsha 410073, China;
4. Department of Precision Instrument, Tsinghua University, Beijing 100084, China
Abstract: The Allan variance method for various types of gyroscopes in various conditions was evaluated for ring laser gyroscope (RLG) noise at room temperature (~25℃) and steady conditions with analyses of the minor RLG noise terms besides the 5 major noise terms. Tests of a Chinese RLG and the widely used MPU 9250 micro-electro-mechanical systems (MEMS) inertial measurement unit (IMU) for static and dynamic conditions show that the Allan variance method can be used to estimate the main noise terms of various types of gyroscopes for static conditions to set the Kalman filter parameters for integrated global navigation satellite system/inertial navigation system (GNSS/INS) with the parameter values depending on the degree of understanding of the gyroscope physics. Several conclusions are given to supplement the classical Allan variance method in the IEEE Standard Specification Format Guide and Test Procedure for Single Axis Interferometric Fiber Optic Gyros. The wide applicability of the Allan variance method is contrasted with some commonly used data analysis methods. The Allan variance method has been widely recognized for metrology of precise instruments to improve the design and manufacture of precise instruments and to improve the precision of inertial measurements.
Key words: ring laser gyroscope (RLG)     micro-electro-mechanical systems (MEMS) inertial measurement unit (IMU)     MPU 9250     Allan variance (AVAR)     assessment    

在概率论与数理统计中,人们常用方差(又称经典方差)来衡量随机变量或一组数据的离散程度。然而,在使用经典方差分析铯原子频标的频率稳定性时,经典方差会随着时间变长而发散。为此,J. A. Barnes和D. W. Allan共同提出了一种时域分析技术,即双样本方差,通常也称为Allan方差(Allan variance, AVAR)[1],并成功将其用于石英振荡器频率稳定性研究[2]。国际电气和电子工程师协会(IEEE)在《IEEE单轴干涉式光纤陀螺仪试验程序和格式指南的标准规范》(IEEE Standard Specification Format Guide and Test Procedure for Single Axis Interferometric Fiber Optic Gyros, 简称IEEE规范)[3]中,正式将Allan方差方法推荐为频率稳定度的时域分析方法。

相比经典方差只能从总体上分析误差,Allan方差可用来分析导致误差的潜在的随机过程的数学特征。这不仅有助于识别观测数据中存在的已知误差项的来源,还有助于确定仪器自身固有的或其内部不明的误差来源。Allan方差可单独作为一种数据分析的方法,也可作为频域分析技术的一种补充,相关技术及其优劣势分析如表 1所示。Allan方差可用于振荡器及具有相似特征的精密仪器的噪声分析,其应用范围主要取决于对仪器物理特性的了解程度,目前Allan方差已用于原子钟、陀螺仪、加速度计、全球导航卫星系统(GNSS)接收机等多种精密仪器的噪声分析[4-12]

表 1 常用数据分析方法的优劣势比较
分析方法 量化噪声 角度随机游走/角速率白噪声 零偏不稳定性 速率随机游走/加速率白噪声 速率斜坡/趋势项 周期项 指数相关噪声/Markov噪声
Allan方差[2, 5, 6] ☆☆☆ ☆☆☆ ☆☆☆ ☆☆☆ ☆☆☆
频谱/功率谱密度分析[4] ☆☆ ☆☆☆
N秒平均法[13] ☆☆ ☆☆ ☆☆
相关分析/时间序列分析[4] ☆☆☆
回归分析[4] ☆☆☆
数据时域图[4] ☆☆
注:☆的数量表示适用性强度。

本文系统总结了Allan方差方法用于分析室温(~25℃)下静态环形激光陀螺仪噪声的经典方法,对该陀螺仪常见的5个主要噪声项以外的其他噪声项作了补充说明[3]。为验证Allan方差方法用于不同状态下、不同类型陀螺仪噪声分析的适用性,对中国自主生产的某型号环形激光陀螺仪和目前使用较为广泛的MPU 9250微系统(micro-electro-mechanical systems, MEMS)惯性测量单元(inertial measurement unit, IMU)中的陀螺仪分别进行了静态和动态下的多次实验与对比分析,对IEEE规范中经典Allan方差方法作了性能分析和补充解释。本研究为深入理解Allan方差分析方法利弊、改进精密仪器的设计与制造、提高惯性产品实际使用精度等提供了参考。

1 使用Allan方差分析环形激光陀螺仪噪声方法概述 1.1 Allan方差基本原理

Allan方差与经典方差的区别在于:经典方差是各样本yi(τ)与总体均值y(τ)作差再取数学期望,即

$ \sigma _y^2\left( \tau \right) = \left\langle {{{\left[ {{y_i}\left( \tau \right) - \bar y\left( \tau \right)} \right]}^2}} \right\rangle . $ (1)

其中〈 〉为总体均值。Allan方差是用一子集均值yk(τ)与下一子集均值yk+1(τ)作差再取数学期望的1/2,即

$ \sigma _{\bar y}^2\left( \tau \right) = \frac{1}{2}\left\langle {{{\left[ {{{\bar y}_{k + 1}}\left( \tau \right) - {{\bar y}_k}\left( \tau \right)} \right]}^2}} \right\rangle . $ (2)

其中乘以1/2是为使Allan方差与经典方差在形式上相等。

使用Allan方差分析环形激光陀螺仪观测数据时,先假设数据中的不确定性是由仪器某些具有特定特征的噪声产生的,再根据数据估计每个噪声的协方差[5]。具体而言,该陀螺仪的Allan方差主要与以下5个噪声项有关:量化噪声(quantization noise, QN)、角度随机游走(angle random walk, ARW)、零偏不稳定性(bias instability, BI)、速率随机游走(rate random walk, RRW)和速率斜坡(rate ramp, RR)。

将陀螺仪采样时长为τ0N个数据样本(又称时间序列或数据流)分为时长为τ0, 2τ0, …, 0(其中m<(N-1)/2)的数据簇,并获得在该簇长度上包含在每个簇中的全部数据之和的平均值。Allan方差为簇时间函数的数据簇平均值的双样本方差。

Allan方差是关于陀螺仪速率Ω(t)或角度θ(t)(实际测量值)的函数,且$\theta(t)=\int^{t} \mathit{\Omega}\left(t^{\prime}\right) \mathrm{d} t^{\prime}$。具体而言,环形激光陀螺仪的输出通常为角度增量,此处Allan方差为关于角度增量的函数,因此未设定积分下限。角度测量是在基于t=0 (k=1, 2, 3, …, N)给出的离散时间上进行的,记为θk=θ(0)。时间tktk+τ之间的陀螺仪平均速率为${\mathit{\bar \Omega }_k}(\tau) = \frac{{{\theta _{k + m}} - {\theta _k}}}{\tau }$,其中t=0。因此,环形激光陀螺仪的Allan方差定义为

$ \begin{array}{*{20}{c}} {{\sigma ^2}\left( \tau \right) = \frac{1}{2}\left\langle {{{\left( {{{\mathit{\bar \Omega }}_{k + m}} - {{\mathit{\bar \Omega }}_k}} \right)}^2}} \right\rangle = }\\ {\frac{1}{{2{\tau ^2}}}\left\langle {{{\left( {{\theta _{k + 2m}} - 2{\theta _{k + m}} + {\theta _k}} \right)}^2}} \right\rangle .} \end{array} $ (3)

其估值为

$ \sigma _\mathit{\Omega }^2\left( \tau \right) = \frac{1}{{2{\tau ^2}\left( {N - 2m} \right)}}\sum\limits_{k = 1}^{N - 2m} {{{\left( {{\theta _{k + 2m}} - 2{\theta _{k + m}} + {\theta _k}} \right)}^2}} . $ (4)

由此可知,环形激光陀螺仪的Allan方差与其原始观测数据噪声的功率谱密度有关。Allan方差与双边功率谱密度的关系式为

$ {\sigma ^2}\left( \tau \right) = 4\int_0^\infty {{S_\mathit{\Omega }}} \left( f \right)\frac{{{{\sin }^4}\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}}{{{{\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}^2}}}{\rm{d}}f. $ (5)

其中SΩ(f)为双边速率噪声功率谱密度,即平稳随机过程的功率谱密度。对于非平稳过程,如闪变噪声,需要使用时间平均功率谱密度。

式(5)是在Allan方差整个计算过程中用于表征速率噪声功率谱密度的关键公式,其物理意义为当通过具有sin4x/x2形式的传递函数的滤波器时,Allan方差与陀螺仪速率输出的总噪声功率成比例。

由式(5)和上述解释可以看出,滤波器的带通取决于τ。这表明不同类型的随机过程都可通过调整滤波器的带通(即通过改变τ)检验。因此,Allan方差提供了一种识别和量化数据中存在的各个噪声项的方法,通常在τ与Allan方差平方根σ(τ)的log-log双对数曲线图中对其作进一步分析。

1.2 环形激光陀螺仪各噪声项及其Allan方差分析 1.2.1 量化噪声

量化噪声是由电子仪器数字量化编码采样输出的特性导致的,反映的是仪器的最小分辨率。环形激光陀螺仪示值器只有在陀螺仪相位变化量为预设值时(如2π、π/2、…)才会计数。量化噪声的角度功率谱密度为

$ {S_\theta }\left( f \right) = \tau {Q^2}\left( {\frac{{{{\sin }^2}\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}}{{{{\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}^2}}}} \right) \approx \tau {Q^2}. $ (6)

其中:f<1/(2τ);Q为量化噪声系数,其理论极值为$S/\sqrt {12} $S为陀螺仪比例因子。量化噪声的速率功率谱密度与角度功率谱密度的关系式为SΩ(2πf)=(2πf)2Sθ(2πf),则量化噪声的速率功率谱密度为

$ {S_\mathit{\Omega }}\left( f \right) = \frac{{4{Q^2}}}{\tau }{\sin ^2}\left( {{\rm{ \mathsf{ π} }}f\tau } \right) \approx {\left( {2{\rm{ \mathsf{ π} }}f} \right)^2}{Q^2}\tau . $ (7)

其中f<1/(2τ)。将式(7)代入式(5)得

$ \sigma _{{\rm{QN}}}^2\left( \tau \right) = \frac{{3{Q^2}}}{{{\tau ^2}}}. $ (8)

量化噪声在图 1中对应斜率为-1的部分,它(及其延长线)与$\tau = \sqrt 3 $交点的纵坐标为Q值。

图 1 Allan方差方法分析环形激光陀螺仪噪声示例曲线

相同τ的Allan方差中可能存在闪变角度噪声、白角度噪声等具有不同光谱特性的其他噪声。

1.2.2 角度随机游走

角度随机游走主要是由环形激光陀螺仪的随机高频抖动、光子自发辐射及其他可能的具有比采样时间短得多的时间相关的高频噪声所导致的,绝大多数可通过设计消除。这些噪声项均表现为陀螺仪速率输出中的白噪声频谱。角度随机游走的速率功率谱密度可表示为

$ {S_\mathit{\Omega }}\left( f \right) = {N^2}. $ (9)

其中N为角度随机游走系数。将式(9)代入式(5)得

$ \sigma _{{\rm{ARW}}}^2\left( \tau \right) = \frac{{{N^2}}}{\tau }. $ (10)

角度随机游走在图 1中对应曲线斜率为-1/2的部分,它(及其延长线)与τ=1交点的纵坐标为N值。

1.2.3 零偏不稳定性

零偏不稳定性又称“1/f噪声”,即其功率谱密度与频率成反比。该噪声主要是由环形激光陀螺仪的放电组件、电子设备或易受随机闪变影响的其他组件所导致的。由于零偏不稳定性具有低频特性,因此表现为数据中的偏差波动。零偏不稳定性的速率功率谱密度为

$ {S_\mathit{\Omega }}\left( f \right) = \left\{ {\begin{array}{*{20}{l}} {\left( {\frac{{{B^2}}}{{2{\rm{ \mathsf{ π} }}}}} \right)\frac{1}{f},}&{f \le {f_0};}\\ {0,}&{f > {f_0}.} \end{array}} \right. $ (11)

其中:B为零偏不稳定系数,f0为截止频率。将式(11)代入式(5)得

$ \sigma _{{\rm{BI}}}^2\left( \tau \right) = \frac{{2{B^2}}}{{\rm{ \mathsf{ π} }}}\int_0^\infty {\frac{{{{\sin }^4}\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}}{{{{\left( {{\rm{ \mathsf{ π} }}f\tau } \right)}^3}}}{\rm{d}}\left( {{\rm{ \mathsf{ π} }}f\tau } \right)} \approx \frac{{4{B^2}}}{9}. $ (12)

零偏不稳定性在图 1中对应曲线斜率为0的部分,它(及其延长线)与τ=1交点的纵坐标为0.664B。图 1曲线显示, 当τ>1/f0时,零偏不稳定性的Allan方差曲线近似水平。因此,检查曲线中近似水平区域,可估计零偏不稳定性的极限值和潜在的闪变噪声的截止频率。

1.2.4 速率随机游走

速率随机游走是一个不确定起源的随机过程,可能是具有较长相关时间的指数相关噪声的极限情况。机械陀螺仪和速率偏置激光陀螺仪都表现出该噪声。速率随机游走的速率功率谱密度为

$ {S_\mathit{\Omega }}\left( f \right) = {\left( {\frac{K}{{2{\rm{ \mathsf{ π} }}}}} \right)^2}\frac{1}{{{f^2}}}. $ (13)

其中K为速率随机游走系数。将式(13)代入式(5)得

$ \sigma _{{\rm{RRW}}}^2\left( \tau \right) = \frac{{{K^2}\tau }}{3}. $ (14)

速率随机游走在图 1中对应曲线斜率为1/2的部分,它(及其延长线)与τ=3交点的纵坐标为K值。

1.2.5 速率斜坡

速率斜坡不是随机噪声,更像是确定性噪声。观测数据表明,在持续较长时间时,环形激光陀螺仪的强度存在非常缓慢且单调变化,这可能是由于在一固定方向并较长时间段测试平台有一非常小的加速度所致。速率斜坡可能是环形激光陀螺仪的真实输入,即

$ \mathit{\Omega } = Rt. $ (15)

其中R为速率斜坡系数。通过形成和操作包含由式(15)给出的输入的数据簇可得

$ \sigma _{{\rm{RR}}}^2\left( \tau \right) = \frac{{{R^2}{\tau ^2}}}{2}. $ (16)

速率斜坡在图 1中对应曲线斜率为1的部分,它(及其延长线)与$\tau = \sqrt 2 $交点的纵坐标为R值。速率斜坡的速率功率谱密度为

$ {S_\mathit{\Omega }}\left( f \right) = \frac{{{R^2}}}{{{{\left( {2{\rm{ \mathsf{ π} }}f} \right)}^3}}}. $ (17)

相同τ的Allan方差中还可能存在具有1/f3功率谱密度的闪变加速度噪声。

1.2.6 其他噪声项

在环形激光陀螺仪模型方程中,除了包括以上5个主要噪声项之外,还可能存在指数相关噪声、正弦噪声等其他噪声项[3],如图 1所示。

1.2.7 所有噪声的综合影响

实验Allan方差示例曲线及环形激光陀螺仪主要噪声项的Allan方差分析如图 1表 2所示。

表 2 环形激光陀螺仪主要噪声项及其Allan方差分析
曲线斜率 主要噪声项 噪声系数 τ/s 交点纵坐标
-1 量化噪声 Q/(°) $\sqrt 3 $ Q
-1/2 角度随机游走 N/((°)·s-1/2) 1 N
0 零偏不稳定性 B/((°)·s-1) 0 0.664B
1/2 速率随机游走 K/((°)·s-3/2) 3 K
1 速率斜坡 R/((°)·s-2) $\sqrt 2 $ R

绝大多数情况下,不同的噪声项出现在不同τ区域,因此可以很容易地识别出观测数据中存在的各种随机过程。假设已知的各随机过程均相互独立,则任意给定τ的Allan方差即为相同τ的各独立随机过程的Allan方差之和,

$ \begin{array}{*{20}{c}} {\sigma _{{\rm{total}}}^2\left( \tau \right) = \sigma _{{\rm{QN}}}^2\left( \tau \right) + \sigma _{{\rm{ARW}}}^2\left( \tau \right) + \sigma _{{\rm{BI}}}^2\left( \tau \right) + }\\ {\sigma _{{\rm{RRW}}}^2\left( \tau \right) + \sigma _{{\rm{RR}}}^2\left( \tau \right) + \sigma _{{\rm{others}}}^2\left( \tau \right).} \end{array} $ (18)

因此,估计在任何τ区域中给定随机噪声的幅度时,需要知道相同区域中的其他随机噪声的幅度。

1.3 估计精度和测试设计

因为有限长度的数据集只能生成有限数量的数据簇,且任何噪声的Allan方差只能通过一定长度的数据簇的总量进行估计,所以给定τ的Allan方差的精度取决于数据集内独立数据簇的数量。当使用具有n个点的数据集中包含k个点的n/k个数据簇时,其σ(τ)的1σ均方根误差为

$ \sigma = \frac{1}{{\sqrt {2\left( {\frac{n}{k} - 1} \right)} }}. $ (19)

由此可见,数据簇长度越短,则独立数据簇数量越多,估计误差越小;反之,数据簇长度越长,则独立数据簇数量越少,估计误差越大。

式(19)可用于设计测试来观察在给定精度范围内具有明确特性的特殊噪声,例如验证一个随机过程的存在。假设一随机过程存在于观测数据中,且疑似特征时间为100s,误差限定在25%以内,则由式(19)得σ=0.25,kmax=n/9,创建相同长度的数据簇,则测试所需总时长为100s × 9=900s。

2 实验分析 2.1 实验方法

在采样频率设置上,根据频谱分析理论,改变采样频率会影响Allan方差分析结果,为此要在陀螺仪带宽范围内尽量提高采样频率,通常取2倍带宽频率(4~6倍或以上更佳)。环形激光陀螺仪输出为角度增量,其中隐含了积分(平均)过程,因此其采样频率并不影响Allan方差的角度随机游走系数。在采样时长设置上,可根据实际需要并结合式(19)确定,亦可根据工程经验(通常观测超过2h)确定。

本文采用一款中国自主研发的某型号环形激光陀螺仪,室温下通电静置2h后,设置采样频率1000Hz,并记录其5h观测数据,其输出数据分别为协调世界时(UTC)(单位:s)、陀螺仪三轴角度增量(单位:(°))、加速度计三轴速度增量(单位:m/s)。另一款为广泛使用的MPU 9250 MEMS IMU,室温下通电静置2h后,设置采样频率100Hz,并记录其8h观测数据。输出数据分别为GPS时(GPST)(单位:s)、陀螺仪三轴角速率(单位:(°)/s)、加速度计三轴加速度(g)。

2.2 Allan方差计算方法

使用Allan方差计算数据时,输入量为平均角速率序列Ω1(τ0), Ω2(τ0), …, ΩN(τ0)(单位:(°)/s),取样频率为f (单位:Hz)(即1/τ0τ0为时间间隔,单位为s)。环形激光陀螺仪输出量为角增量A(单位:(°)),需先将其除以采样时间间隔τ0;MEMS陀螺仪和光纤陀螺仪输出量为离散采样点处的瞬时角速率(单位:(°)/s),可近似作为τ0时间段内的平均角速率。在采样频率满足Shannon采样定理条件下误差可忽略;反之,在不满足该条件时应考虑误差的影响。

按“2k”方法计算数据并绘图分析Allan方差最为高效[14],具体计算步骤如下:

步骤1    由式(4)计算取样时间间隔为τ0时的Allan方差${\hat \sigma ^2}\left({{\tau _0}} \right) = \frac{1}{2}\frac{1}{{N - 1}}\sum\limits_{k = 1}^{N - 1} {{{\left[{{{\mathit{\bar \Omega }}_{k + 1}}\left({{\tau _0}} \right) - {{\mathit{\bar \Omega }}_k}\left({{\tau _0}} \right)} \right]}^2}} $

步骤2    将取样时间间隔加倍,即τ1=2τ0,则$N_{1}=\lfloor N / 2\rfloor$,其中$\left\lfloor {\; \;} \right\rfloor $表示向下取整,相邻两个序号的平均角速率为${\mathit{\bar \Omega }_k}\left({{\tau _1}} \right) = \frac{1}{2}\left[{{{\mathit{\bar \Omega }}_{2k - 1}}\left({{\tau _0}} \right) + {{\mathit{\bar \Omega }}_{2k}}\left({{\tau _0}} \right)} \right]$,其中k=1, 2, 3, …, N1。由此组成新的平均角速率序列为Ω1(τ1), Ω2(τ1), Ω3(τ1), …, ΩN1(τ1),取样时间间隔为τ1。显然,新的平均角速率序列长度减半(亦可能相差1个数据)。

由式(4)计算取样时间间隔为τ1时的Allan方差${\hat \sigma ^2}\left({{\tau _1}} \right) = \frac{1}{2}\frac{1}{{{N_1} - 1}}\sum\limits_{k = 1}^{{N_1} - 1} {{{\left[{{{\mathit{\bar \Omega }}_{k + 1}}\left({{\tau _1}} \right) - {{\mathit{\bar \Omega }}_k}\left({{\tau _1}} \right)} \right]}^2}} $

步骤3     同理,将取样时间间隔再次加倍,即τ2=2τ1=22τ0,则${N_2} = \left\lfloor {{N_1}/2} \right\rfloor $,新的平均角速率序列为Ω1(τ2), Ω2(τ2), Ω3(τ2), …, ΩN1(τ2),取样时间间隔为τ2,其中${{\mathit{\bar \Omega }}_k}\left({{\tau _2}} \right) = \frac{1}{2}\left[{{{\mathit{\bar \Omega }}_{2k - 1}}\left({{\tau _1}} \right) + {{\mathit{\bar \Omega }}_{2k}}\left({{\tau _1}} \right)} \right]$k=1, 2, 3, …, N2。根据式(4)计算取样时间间隔为τ2时的Allan方差${{\hat \sigma }^2}\left({{\tau _2}} \right) = \frac{1}{2}\frac{1}{{{N_2} - 1}}\sum\limits_{k = 1}^{{N_2} - 1} {{{\left[{{{\mathit{\bar \Omega }}_{k + 1}}\left({{\tau _2}} \right) - {{\mathit{\bar \Omega }}_k}\left({{\tau _2}} \right)} \right]}^2}} $

步骤4     重复上述步骤,将取样时间间隔继续加倍,即τi=2τi-1=2iτ0,则${N_i} = \left\lfloor {{N_{i - 1}}/2} \right\rfloor $,直到最终所得序列长度不小于2,则新的平均角速率序列Ω1(τi), Ω2(τi), Ω3(τi), …, ΩNi(τi)取样时间间隔为τi。由式(4)计算取样时间间隔为τi时的Allan方差${{\hat \sigma }^2}\left({{\tau _i}} \right) = \frac{1}{2}\frac{1}{{{N_i} - 1}}\sum\limits_{k = 1}^{{N_i} - 1} {{{\left[{{{\mathit{\bar \Omega }}_{k + 1}}\left({{\tau _i}} \right) - {{\mathit{\bar \Omega }}_k}\left({{\tau _i}} \right)} \right]}^2}} $

步骤5     此时已经得到一系列点对,即${\tau _\alpha } \sim {{\hat \sigma }^2}\left({{\tau _\alpha }} \right)$${\tau _\alpha } \sim \hat \sigma \left({{\tau _\alpha }} \right)$,其中α=1, 2, 3, …, i (τα单位为s,$\hat \sigma \left({{\tau _\alpha }} \right)$单位为(°)/s)。然而,$\hat \sigma \left({{\tau _\alpha }} \right)$的单位通常为(°)/h,需要用$1\; \left({^ \circ } \right)/{\rm{s = }}\left({\frac{{180}}{{\rm{ \mathit{ π} }}}} \right)/\left({\frac{1}{{3600}}} \right)\left({^ \circ } \right)/{\rm{h}}$进行单位换算。以此绘制图 2,即Allan方差${{\hat \sigma }^2}\left({{\tau _\alpha }} \right)$随子集时长τα变化的双对数曲线。

图 2 环形激光陀螺仪三轴原始角度增量随时间变化及其Allan方差分析

按照相同步骤,对比分析环形激光陀螺仪和MEMS陀螺仪在不同情况下的Allan方差曲线,具体情况分别如图 35所示。

图 3 环形激光陀螺仪与MEMS陀螺仪Allan方差对比

图 4 MEMS陀螺仪原始观测值低通滤波前后Allan方差对比

图 5 环形激光陀螺仪动态原始角度增量随时间变化及其Allan方差分析

2.3 数据分析

1) 对比实际分析曲线(见图 2)和示例分析曲线(见图 1)可知,在某一陀螺仪某一时段的观测数据中,5个主要噪声项不一定同时全部表现出来,可能只表现出其中几个。可能原因如下:(1)某些斜率的Allan方差不存在,(2)某些噪声的斜率与主要噪声的斜率相同导致实验曲线重叠,(3)某些噪声量级太小被量级较大的噪声掩盖,(4)由于观测时间不足导致某些斜率的Allan方差还未表现出来。

2) 由实际分析曲线(见图 2)可知,在τ取10-3~10时,其曲线斜率约为-1,环形激光陀螺仪主要表现为量化噪声;在τ取10~103时,其曲线斜率约为-1/2,主要表现为角度随机游走;而零偏不稳定性、速率随机游走、速率斜坡及其他噪声项未明显显示。由此可知,环形激光陀螺仪的噪声主要为量化噪声,其各噪声项系数分析如表 3所示;而MEMS陀螺仪的量化噪声并不明显,但其角度随机游走明显,如图 3所示。由此可知,一般情况下不同类型陀螺仪的Allan方差曲线表现出的噪声项不同。

表 3 环形激光陀螺仪Allan方差分析结果
仪器类型 标称精度/((°)·h-1) 采样频率/Hz 采样时间/h Allan方差分析各误差项系数 实测精度/((°)·h-1)
Q/(″) N/((°)·h-1/2) B/((°)·h-1) K/((°)·h-3/2) R/((°)·h-2)
环形激光陀螺仪 0.05 1 000 5 ~0.29 ~0.09 ~0.05

3) 在使用Allan方差评价陀螺仪性能时,某些国家学者采用曲线最低点纵坐标值作为陀螺仪精度[15]。该方法最主要的弊端在于不符合实际操作情况,通常不会将陀螺仪静置“103 s”如此长时间,因此不建议采用该方法。在此问题上,中国学者通常采用τ=10s时σ(τ)的平滑结果作为陀螺仪实际使用精度,更符合实际情况。因此,建议采用此方法。据此分析图 2,结果见表 3

4) 不同噪声项的Allan方差通常在实验曲线中显示的时间轴坐标位置不同。如实际分析曲线图 2所示,该环形激光陀螺仪的量化噪声为“s”量级、角度随机游走为“10s”量级,以此类推,零偏不稳定性、速率随机游走和速率斜坡等可能均为“100s”量级。由此可以推断,主要噪声项之间相互独立或者各噪声项之间相关性很小。同时,其他噪声项,如与量化噪声具有相同斜率的闪变角度噪声、白角度噪声,与零偏不稳定性具有相同斜率的闪变噪声,与速率斜坡具有相同斜率的闪变加速度噪声,以及指数相关噪声、正弦噪声等,均未明显显示。

5) 对比分析环形激光陀螺仪与MEMS陀螺仪的Allan方差曲线(见图 3)可知,陀螺仪精度越高,其Allan方差曲线总体上离横轴越近;陀螺仪采样频率越高,其Allan方差曲线总体上离纵轴越近。由此可以直观概略地判断陀螺仪精度和采样设置情况。

6) 某些Allan方差分析曲线中的“N”形部分(如图 4中标记处)并不是Markov噪声,而是由高频信号经低通滤波处理后所导致的,该部分反映的是滤波特性而非陀螺仪性能。为此,Allan方差通常仅用于分析未经滤波处理的陀螺仪原始观测数据。然而,根据工程经验,目前商用MEMS IMU输出的陀螺仪数据绝大多数都已经过滤波处理,为此在使用Allan方差检验时少有示例曲线(如图 1所示)出现。

7) 某些情况下,一些实际观测数据的Allan方差分析曲线可能为“W”形或“M”形曲线(如图 5所示),与示例曲线(如图 1所示)中的“U”形或“V”形曲线区别较大,此时Allan方差不能用于陀螺仪噪声分析,应查看输入数据或观测数据的正误。

8) Allan方差分析的是陀螺仪静态时的观测数据,而在全球导航卫星系统/惯性导航系统(GNSS/INS)组合系统导航定位等动态应用中的陀螺仪噪声(如图 5所示),还需根据具体情况具体分析:(1)动态情况下陀螺仪各项噪声的存在与否、变化情况、量级大小等需进一步通过具体实验验证分析;(2)在设置Kalman滤波噪声参数时,主要使用的是角度随机游走系数(用于设置Q阵)和零偏不稳定性系数(用于设置一阶Markov过程方差),过多噪声项并不能显著提高精度,甚至可能显著降低计算效率,因此其他噪声项应酌情处理。

9) 由实验方法可知,Allan方差分析的是静基座陀螺仪随机漂移误差。这是由于动基座变角速率会引起额外的功率谱分量,从而影响陀螺仪Allan方差性能参数分析的准确性。Allan方差滤波器在零频率处是截止的,这是由于计算Allan方差时式(3)中进行了一次求差运算所致,为此静基座常值角速率干扰并不影响Allan方差计算结果,但也因此无法通过Allan方差分析获得陀螺仪随机常值漂移参数。

3 结论

通过对Allan方差方法分析环形激光陀螺仪噪声的系统梳理,以及对环形激光陀螺仪、MEMS陀螺仪等不同类型不同状态陀螺仪的对比实验分析,可以得知:Allan方差可以对静态陀螺仪的主要噪声项进行识别和概算(如表 1所示),并对陀螺仪实际使用精度进行估计,从中得到的陀螺仪角度随机游走系数和零偏不稳定性系数等参数可用于GNSS/INS组合系统融合算法研究。

但同时也要看到,在使用IEEE规范经典方法上要注意“深”和“浅”的结合:

“深”是指深入理解,即:1)由陀螺仪噪声Allan方差公式推导可知,陀螺仪中不仅存在5个主要噪声项,还有与5个主要噪声项同斜率的其他噪声项,以及除5个主要噪声项以外的其他噪声项。2)根据Allan方差分析原理可知,某型号陀螺仪中具体包含的噪声项需要事先深入了解该陀螺仪的物理特性,再有针对性地进行识别和分析。3)从Allan方差分析原理可知,Allan方差分析曲线的斜率仅表示噪声过程的类型,并不一定是其原因。例如,白相位噪声可能是由量化噪声引起的,也可能是由其他原因引起的。4)国外学者采用Allan方差分析曲线最低点对应的σ(τ)值作为陀螺仪实际精度的方法并不合理,以此评定一些超高精度陀螺仪的实际使用精度仍有待进一步验证。中国学者通常采用τ=10s时σ(τ)平滑结果作为陀螺仪的实际使用精度,这种做法更符合实际情况。

“浅”是指概略分析,即使用Allan方差分析各噪声项时,通常采取图 1中经典的“五段线”法概算主要噪声项即可满足一般应用需求[16]。使用最小二乘回归、Kalman滤波及其衍生方法对各项噪声严格分析评定的方法效率不高[17]

4 展望

鉴于IEEE规范解释和已有实验验证,关于Allan方差分析陀螺仪或其他精密仪器噪声方法的使用在以下方面仍有待进一步补充完善:

1) 目前Allan方差主要用于分析室温下陀螺仪静态数据,由观测条件对噪声影响的基本原理可知,观测环境温度变化、载体运动状态等对Allan方差分析陀螺仪噪声影响情况需补充分析。

2) 陀螺仪噪声是一次启动稳定性,而不具有逐次启动重复性,逐次启动噪声无法确定。但在实际使用基于GNSS/INS组合系统时,不可避免出现逐次启动情况,该噪声确定方法有待进一步研究。

3) Allan方差方法通过实验曲线分析各噪声项时假设其相互独立,但各噪声项是否严格相互独立,或者各噪声项之间的相关程度,仍需通过严密的数学推导进一步证明。

4) 采用Allan方差方法分析精密仪器噪声的前提是对精密仪器特性已经深入了解,该方法对于高精度钟、陀螺仪、加速度计、GNSS接收机等及其他精密仪器的适用性,还有待进一步补充分析。

致谢

本研究还得到了美国国家标准与技术研究所Judah Levine研究员、西北工业大学严恭敏副教授的指导与帮助,谨致谢意。

参考文献
[1]
ALLAN D W, LEVINE J. A historical perspective on the development of the Allan variances and their strengths and weaknesses[J]. IEEE Transactions on Ultrasonics, Ferroelectrics & Frequency Control, 2016, 63(4): 513-519.
[2]
ALLAN D W. Statistics of atomic frequency standards[J]. Proceedings of the IEEE, 1966, 54(2): 221-230. DOI:10.1109/PROC.1966.4634
[3]
IEEE. IEEE standard specification format guide and test procedure for single axis interferometric fiber optic gyros: IEEE Std 647-2006[S]. Piscataway, USA: IEEE, 2006.
[4]
严恭敏, 李四海, 秦永元. 惯性仪器测试与数据分析[M]. 北京: 国防工业出版社, 2012.
YAN G M, LI S H, QIN Y Y. Instrument test and data analysis of inertial devices[M]. Beijing: National Defense Industry Press, 2012. (in Chinese)
[5]
危志英, 汪世林.基于Allan方差方法激光陀螺仪性能评价方法[C]//全国第12届空间及运动体控制技术学术年会论文集.北京: 北京自动控制设备研究所, 2006: 500-504.
WEI Z Y, WANG S L. Performance evaluation method of laser gyroscope based on Allan variance method[C]//Proceedings of the 12th Annual Conference on Space and Motion Control Technology. Beijing: Beijing Institute of Automatic Control Equipment, 2006: 500-504. (in Chinese)
[6]
刘巧光, 许辅义, 滕云鹤, 等. 环形激光陀螺仪随机误差模型的研究[J]. 清华大学学报(自然科学版), 1999, 39(2): 71-74.
LIU Q G, XU F Y, TENG Y H, et al. Investigation on random error model of ring laser gyro[J]. Journal of Tsinghua University (Science and Technology), 1999, 39(2): 71-74. DOI:10.3321/j.issn:1000-0054.1999.02.019 (in Chinese)
[7]
EL-SHEIMY N, HOU H, NIU X. Analysis and modeling of inertial sensors using Allan variance[J]. IEEE Transactions on Instrumentation and Measurement, 2008, 57(1): 140-149. DOI:10.1109/TIM.2007.908635
[8]
赵思浩, 陆明泉, 冯振明. MEMS惯性器件误差系数的Allan方差分析方法[J]. 中国科学:物理学力学天文学, 2010, 40(5): 672-675.
ZHAO S H, LU M Q, FENG Z M. Allan variance analysis on error coefficients of MEMS inertial components[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2010, 40(5): 672-675. (in Chinese)
[9]
ZHANG Q, NIU X, CHEN Q, et al. Using Allan variance to evaluate the relative accuracy on different time scales of GNSS/INS systems[J]. Measurement Science and Technology, 2013, 24(8): 085006. DOI:10.1088/0957-0233/24/8/085006
[10]
张小红, 朱锋, 薛学铭, 等. 利用Allan方差分析GPS非差随机模型特性[J]. 测绘学报, 2015, 44(2): 119-127.
ZHANG X H, ZHU F, XUE X M, et al. Using Allan variance to analyze the zero-differenced stochastic model characteristics of GPS[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(2): 119-127. (in Chinese)
[11]
刘乾坤, 隋立芬, 陈泉余, 等. 基于Allan方差的BDS随机误差分析[J]. 测绘通报, 2016(5): 18-21, 65.
LIU Q K, SUI L F, CHEN Q Y, et al. Stochastic errors analysis of BDS based on Allan variance[J]. Bulletin of Surveying and Mapping, 2016(5): 18-21, 65. (in Chinese)
[12]
吴富梅, 张晓东. 几种不同时频分析法对INS信号的分析和比较[J]. 测绘科学技术学报, 2010, 27(2): 92-96.
WU F M, ZHANG X D. Analysis and comparison of INS signal based on several methods[J]. Journal of Geomatics Science and Technology, 2010, 27(2): 92-96. DOI:10.3969/j.issn.1673-6338.2010.02.005 (in Chinese)
[13]
吕品, 刘建业, 赖际舟, 等. 光纤陀螺的随机误差性能评价方法研究[J]. 仪器仪表学报, 2014, 35(2): 412-418.
LV P, LIU J Y, LAI J Z, et al. Research on the performance evaluation methods of fiber optical gyro stochastic errors[J]. Chinese Journal of Scientific Instrument, 2014, 35(2): 412-418. (in Chinese)
[14]
DU X, ZENG C, LI H. Comparison of random error analysis methods for fiber optic gyro based on Allan variance[J]. Science Discovery, 2017, 5(5): 375-379. DOI:10.11648/j.sd.20170505.22
[15]
SENSONOR. STIM210 datasheet[EB/OL]. (2018-08-01)[2019-09-10]. https://www.sensonor.com/products/gyro-modules/stim210/.
[16]
张代兵. 一种基于Allan方差方法的激光陀螺性能评价方法[J]. 仪器仪表学报, 2004(S1): 715-717.
ZHANG D B. Laser gyro performance estimate method based on Allan variance[J]. Chinese Journal of Scientific Instrument, 2004(S1): 715-717. (in Chinese)
[17]
REHDERJ. IMU noise model[EB/OL]. (2016-10-11)[2019-09-11]. https://github.com/ethz-asl/kalibr/wiki/IMU-Noise-Model.