2. 中国人民解放军 96401部队, 宝鸡 721006
2. Unit 96401 of the Chinese People's Liberation Army, Baoji 721006, China
惯性导航系统的准确度很大程度上取决于其核心部件陀螺的精度[1]。陀螺误差通常可以分为确定性误差和随机误差。确定性误差可以通过试验标定其误差系数以进行补偿。随机误差是对陀螺精度影响较大又难以处理的误差项，通过对其分析，一方面可以辨识出其噪声来源，优化内部构造；另一方面可以检测其对外界随机激励的响应特征，改进滤波算法，进而达到提高导航精度的目的。目前，IEEE推荐使用Allan方差方法对陀螺随机误差进行分析，但该方法只能对静态条件下随机误差进行分析，不具有时变特性，难以表征动态条件下陀螺输出的非平稳变化特性。近年来，有研究人员运用动态Allan方差(Dynamic Allan Variance，DAVAR)对光学陀螺、MEMS陀螺及半球谐振陀螺进行动态特征分析，提取出动态条件下随时间变化的陀螺随机误差系数，观测随机误差在外界干扰下的变化规律[2-4]。
经典的Allan方差方法存在一定的不足，对于固定长度的数据，其可划分的独立组数越多，方差估计的置信度越高，估计误差会随着独立分组数量的减小而增大；此外其用于计算的相关时间只能达到样本总时间的二分之一，在长相关时间下随机误差的估计准确度相对较低[5]。而DAVAR方法用固定长度的窗函数去截断量测信号，使参与计算的数据大大减少，再用Allan方差方法对截断窗内数据进行噪声分离提取，必将进一步放大Allan方差的缺点，使估计值的置信度进一步下降。增加截断窗长又会引起动态跟踪能力的下降，与DAVAR方法的设计原则相背离[6-7]。因此，如何在样本数据量不变的情况下提高方差估计的置信度，就成为亟待解决的问题。
实际上，针对样本数据量固定不变情况下的方差估计置信度提高方法已经展开了不少研究，主要有以下几种方法：①交叠Allan方差，该方法利用重叠采样获得相关时间τ下尽可能多的子序列，提高了样本内现有数据的利用率[8]。②Allan总方差，该方法同Allan方差的采样方式相同，通过对数据两端进行延拓的方法来增加参与计算的数据量，从而提高长相关时间下随机误差的估计精度[9]。③Theo1方差和TheoH方差，采用了二次采样的方式来增加方差估计置信度[10-11]。TheoH方差是偏差修正后的Theo1方差与Allan方差合成的一种方差估计方法，因此又称混合理论方差。该方法在噪声分量成分复杂以及相关时间长的情况下，能够从同一样本中提取出更多的特征信息[12]。文献[5]用TheoH方差对静态条件下光学陀螺的误差统计特性进行了分析，同Allan方差相比，TheoH方差对长相关时间误差系数具有更高的辨识能力[5]。
本文运用TheoH方差替代Allan方差对截断窗内样本数据进行分析，改善了DAVAR方法对陀螺随机误差的估计准确度。对仿真信号和陀螺实测信号的动态特征提取结果证明，相对于DAVAR方法来说本文方法有更好的误差系数提取准确度。
1 Allan方差和TheoH方差 1.1 Allan方差Allan 方差是一种经典的数据统计方法，主要用于对信号的时域稳定性分析。利用Allan方差同功率谱密度之间的定量关系，可以分离出陀螺随机误差的类型和大小，即将Allan方差表示为各类噪声的平方和：
|
(1) |
式中： σQ2为量化噪声(简称Q)的方差；σN2为角度随机游走(简称N)的方差；σB2为零偏不稳定性(简称B)的方差；σK2为角速率随机游走(简称K)的方差；σR2为速率斜坡(简称R)的方差；τ为相关时间。其中，量化噪声和角度随机游走为高频噪声系数，属于在短相关时间下就可以辨识出的噪声系数，用Allan方差或交叠Allan方差可以取得很好的效果。而零偏不稳定性、角速率随机游走和速率斜坡属于中、低频噪声系数，需要较长的相关时间才能辨识出来。DAVAR用截断窗将数据截成更短的样本，必然会大大降低后者的估计精度。
1.2 TheoH方差Theo1方差可以更为准确地对长相关时间下随机误差的特征进行辨识[13]。但由于该方法与Allan方差的采样方式不完全相同，造成Theo1方差的估计值偏离了Allan方差，必须对其进行修正。Theo1方差相关时间从10τ0到0.75 Nxτ0 (其中，τ0为采样周期，Nx为总采样个数)，比Allan方差要长。对于时间序列x(t)(t=1,2,…,Nx)，采样间隔为τ0，总时间为T=Nxτ0，Theo1方差采取二次采样的方式来增加参与计算的独立数组，也就是在相关时间τ=mτ0内，10≤m≤Nx-1，m为偶数，采样的步幅变为τs=0.75 mτ0，进行第2次采样。Theo1方差有2个采样时间：τ0和τs，τs与Allan方差中的相关时间意义相同。图 1为Theo1方差变步长采样过程。
|
| 图 1 Theo1方差的采样原理 Fig. 1 Sampling principle of Theo1 variance |
Theo1方差的计算公式如下：
|
(2) |
当m较小时，Theol方差的估计置信度与Allan方差基本一致；但在提取长相关时间随机误差系数时，为得到较大的m，在对量测信号的采样过程中，尽可能地使τ0更小以得到较大的数据量，从而分离出更多噪声特征信息。
TheoBR(m,τ0,Nx)=
|
(3) |
TheoH(m,τ0,Nx)=
|
(4) |
式中：TheoBR(Bias-Removed version of Theo1)为自动补偿Theo1方差与交叠Allan方差之间偏差的方差计算函数，TheoBR方差与交叠Allan方差合成后的曲线称为TheoH(Hybrid-TheoBR)；k为τ≤10%T条件下的最大值(T为数据的总时间长度)。选择10%T为切换点是因为Allan方差在相关时间≤10%T时能提供足够的置信度[12]。
方差的开平方是标准差，用标准差曲线图可以更形象地表征信号的变化特征。交叠Allan标准差和TheoBR标准差组合过程如图 2所示。
|
| 图 2 TheoH标准差的合成过程 Fig. 2 Forming process of TheoH deviation |
与文献[5]不同的是，在相关时间小于10%T时，我们采用交叠Allan方差替代经典Allan方差，在大于10%T的相关时间下，用TheoBR方差替代Allan方差，这样既能使短相关时间下辨识的随机误差系数具有比经典Allan方差更高的准确度，又通过引入TheoBR方差提高了长相关时间的方差估计置信度。
2 动态TheoH方差方法设计 2.1 动态Allan方差DAVAR实际上是带有滑动窗的Allan方差，本质是用一个固定窗宽的截断窗在时间序列信号上滑动并分段截取数据，再分别计算截断窗内数据的Allan方差，然后按照时间、相关时间和方差σ2(t,τ)的顺序将Allan方差绘制在一幅三维分析图上，以观测方差随时间和相关时间的动态变化规律，其流程见图 3。
|
| 图 3 DAVAR方法设计流程图 Fig. 3 Flowchart of DAVAR method design |
本文提出用TheoH方差替代图 3灰色框中的Allan方差进行方差的估计，简称动态TheoH方差(Dynamic TheoH Variance，DHVAR)。TheoH方差的统计特性和Allan方差相似，也可以绘制双对数曲线图，并用最小二乘法来辨识随机误差系数。其优势在于对中、长相关时间下噪声系数的辨识上具有更高的置信度，可以更为准确地提取噪声系数，它是相对较新也较好的随机误差统计方法，可以有效地改善DAVAR截断数据造成的随机误差系数辨识准确度下降的问题，可以更为准确地描述出随机误差系数在时间域的变化规律，其步骤如下：
1) 选择待分析动态信号x(t)的时域起始点t1。
2) 用中心点为t1、窗口长度为L(t1)的窗函数去截取信号x(t)，得到截断信号yT(t1)，支撑变量代表截断窗内的时间。
|
(5) |
PW(t′)为长度L(t1)的矩形窗函数，其定义为
|
(6) |
截取得到的信号为
|
(7) |
3) 将yT(t1,t′)同替代方法(TheoH方差)hτ(t′)做卷积建立增量过程Δ(t1,t′,τ)：
|
(8) |
这时变量t′的范围变为
|
(9) |
0≤τ≤τmax，通常τmax=L/3，则t1时刻估计值为
|
(10) |
TheoH方差可以被定义为式(10)的总体期望值：
|
(11) |
4) 通过最小二乘法可以将t1时刻陀螺的主要噪声系数辨识出来：
|
(12) |
|
(13) |
5) 将截断窗平移到t2，得到截断信号yT(t2)，再次计算步骤2)~步骤5)，获得σy2(t2,τ)，不断平移截断窗，能够得到时间域的方差序列σy2(t,τ)和噪声系数A(t)1,A(t)2,…,A(t)5。原则上，为保持数据连贯，通常使以tn+1为中点的滑动截断窗与以tn为中点的滑动窗相交叠。
6) 将随时间变化的方差序列σy2(t,τ)或标准差序列(动态TheoH标准差(Dynamic TheoH Deviation,DHDEV))σy(t,τ)绘制在同一副三维图上(三维为标准差、相关时间和时间)，还可以将各个随机误差系数Q(t1,t2，…,tn)、B(t1,t2，…,tn)等绘制在随时间变化的二维图中，使随机误差系数具备时变特征。
3 试验结果及讨论 3.1 仿真试验通常Allan方差或DAVAR的仿真试验都是用高斯白噪声来模拟，本文借鉴文献[14]，用一阶马尔可夫过程来对噪声系数零偏不稳定性进行模拟，用高斯白噪声驱动的数学模型来仿真其余噪声系数，构建陀螺输出的数学仿真模型。重点模拟零偏不稳定性有3个方面原因：①噪声系数的提取是在假设其为统计独立的前提下进行，只模拟单一噪声系数更有利于验证及比较不同方法的提取效果；②在5个主要噪声源中，零偏不稳定性不只是陀螺相对重要的随机误差，而且经常用于补偿算法或滤波算法的改进[15]；③零偏不稳定性属于中等长度的相关时间，最容易与相邻相关时间的误差项发生耦合，对其辨识效果进行比较，更能客观评判方法性能。
零偏不稳定性仿真信号的方差在0～400 s和400～800 s均为0.001，模拟静态下较小的漂移；而在400～800 s之间方差为0.01，模拟受到振动、冲击、温度变化或其他干扰而发生的较大的漂移。
图 4为模拟的陀螺输出信号，其中浅色部分为仿真的零偏不稳定性时变信号，深色部分为与零偏不稳定性统计不相关的其他随机误差信号，主要由高斯白噪声模拟。
|
| 图 4 陀螺输出随机误差仿真模型 Fig. 4 Stochastic error simulation model of gyro output |
对该陀螺输出仿真信号进行Allan标准差、TheoH标准差和Theo1标准差分析，其方差的双对数曲线如图 5所示。可以看到，在900 s左右时，Theo1方差和TheoH方差的估计值相对稳定。但在不到100 s时，经典Allan方差的估计值就开始波动。Theo1方差在相关时间10 s以后逐渐偏离了Allan方差变化的趋势。只有TheoH方差既将相关时间提高到了3T/4，又可以自动对偏离进行修正，更接近方差的真实变化趋势。然而，这3种方法辨识得到的噪声系数都为固定值，无法表征动态条件下随机误差系数的时变特性。因此，下面分别用DAVAR方法、DHVAR-A方法和DHVAR方法对该仿真数据进行分析，根据辨识结果评价方法的性能。其中，DHVAR-A为经典Allan方差与TheoBR方差混合的方差估计方法，可以用于验证在短相关时间用交叠Allan方差比经典Allan方差能获得更好的辨识准确度的推论。
|
| 图 5 仿真信号的Allan标准差、Theo1标准差和TheoH标准差分析 Fig. 5 Analysis by Allan deviantion，Theo1 deviantion and TheoH deviantion for simulated signal |
由图 6可以看到，DAVAR方法绘制的三维图波动范围为10-1~10-4，而DHVAR方法和DHVAR-A方法绘制的三维图波动范围为10-1~10-2，大大减小了长相关时间下双对数曲线的振荡，也就是提高了方差估计的置信度。
|
| 图 6 实测信号的DAVAR、DHVAR-A和DHVAR分析 Fig. 6 Analysis by DAVAR,DHVAR-A andDHVAR for measured signal |
理论上，DAVAR方法、DHVAR-A方法和DHVAR方法都应该可以把统计独立的噪声源从随机误差里分离出来。也就是说，我们仿真的动态变化的零偏不稳定性噪声系数值应该出现在设定的时间段(400~800 s)及相关时间区域(斜率为0)。然而，DAVAR方法由于可辨识的相关时间短及方差估计置信度低等原因，使得零偏不稳定性动态变化的发生区域与其他随机误差项发生了交叠，其动态变化范围在相关时间1～10 s内都有分布，必将影响到下一步零偏不稳定性噪声系数的准确提取。同时可以看到，DHVAR-A方法和DHVAR方法提取到的动态变化更加独立，与周边分量重合少，准确地把噪声定位在相关时间3～6 s之间，也就是双对数曲线中斜率为0的地方，有效避免了能量泄漏及误差项的相互干扰，但DHVAR-A方法的波动要大于DHVAR方法。
根据文献[2]的结论，DAVAR方法辨识出的时变随机误差系数曲线应围绕经典Allan方差计算得到的误差系数值(用所有数据计算获取的固定值)上下波动，且理想情况下，其时变曲线的均值应该等于Allan方差方法得到的固定值。我们把大数据下Allan方差计算所得的值当作期望值，通过比较不同方法绘制的时变曲线偏离期望值的程度来评价方法的性能[16]。
图 7为分别用DAVAR方法、DHVAR-A方法和DHVAR方法以401窗宽截取仿真信号，进而提取到的零偏不稳定性时间变化曲线。0.026 5为用经典Allan方差对全部时间段仿真信号辨识得到的零偏不稳定性的期望值。通过图 7和表 1可以看到，对于量化噪声和角度随机游走来说，3种方法辨识的结果基本一致，对于中、长相关下可以辨识的噪声系数来说，DHVAR方法的动态曲线均值更贴近期望值。
|
| 图 7 零偏不稳定性的动态提取 Fig. 7 Dynamic extraction of bias instability |
表 1为分别用DAVAR方法、DHVAR-A方法和DHVAR方法对仿真信号的随机误差系数进行辨识所得到的结果，各个随机误差系数值取的是动态变化曲线的均值。可见，用交叠Allan方差替代经典Allan方差可以进一步提高辨识准确度。
| 系数 | 单位 | 期望值 | DAVAR | DHVAR-A | DHVAR |
| Q | μrad | 0.006 3 | 0.006 7 | 0.006 8 | 0.006 7 |
| N | ( )/h1/2 | 0.001 0 | 0.001 1 | 0.001 1 | 0.001 1 |
| B | ( )/h | 0.026 5 | 0.035 4 | 0.030 8 | 0.029 9 |
| K | (°)/h3/2 | 0.451 9 | 1.052 9 | 0.735 7 | 0.706 2 |
| R | (°)/h2 | 1.119 3 | 9.797 3 | 5.806 5 | 5.622 4 |
3.2 动态试验验证
为验证对陀螺信号的分析效果，本课题组对动态试验中采集的某型号光学陀螺的激励响应信号进行了分析。对该信号进行预处理，得到随机误差信号如图 8所示。该实测信号采样周期为40 ms，共采集了1 600 s，约40 000个数据，用401窗长的DHVAR方法展开三维分析，结果如图 9所示。可以看到，振动对陀螺高频噪声系数的影响比较大，对中、低频噪声系数影响较小。
|
| 图 8 陀螺随机误差输出 Fig. 8 Output of stochastic error of gyro |
|
| 图 9 陀螺实测信号随机误差的DHVAR分析 Fig. 9 DHVAR analysis of gyro measued signal stochastic error |
从图 10中分析可得，在原始信号中出现的第1个波动主要是由量化噪声引起，主要原因是信号采集系统分辨力不能满足高频干扰的采样要求造成的误差。角度随机游走波动不大，说明该陀螺仪采用的光学器件比较稳定，抗干扰能力比较强。第2个高量级振动对除角度随机游走外的其他噪声分量都带来了影响，说明随着振动量级的增加，陀螺的各个噪声系数都受到了影响。
|
| 图 10 各随机误差系数随时间变化曲线 Fig. 10 Time-variant curves of different stochasticerror coefficients |
与对仿真信号的处理方法一样，本文对达到平稳状态后的陀螺随机误差用不同方法进行辨识。也就是以1 000～1 600 s之间这15 000个数据为研究对象，以该段所有Allan方差误差系数估计值为期望值，比较DAVAR方法和DHVAR方法的性能。图 11为零偏不稳定性的时间变化曲线，401窗宽，0.52是用该段所有数据计算得到的期望值。
|
| 图 11 零偏不稳定性的DAVAR和DHVAR动态提取 Fig. 11 Dynamic extraction of bias instabilitywith DAVAR and DHVAR |
表 2为DAVAR方法和DHVAR方法提取出的随机误差系数时间变化曲线均值与期望值的比较。可以看到，2种方法对高频随机误差系数(量化噪声和角度随机游走)的提取效果基本一致，而对低频随机误差系数(零偏不稳定性、角速率随机游走和速率斜坡)来说，DHVAR方法和DHVAR-A方法提取出的估计值更加接近于期望值，结合陀螺实测信号进一步验证了TheoH方差确实是一种高置信度的方差估计方法。
| 系数 | 单位 | 量纲 | 期望值 | DAVAR | DHVAR-A | DHVAR |
| Q | μrad | ×10-3 | 0.400 | 0.047 | 0.043 | 0.041 |
| N | (°)/h1/2 | ×10-4 | 0.011 | 0.011 | 0.011 | 0.011 |
| B | ( )/h | ×10-2 | 0.520 | 0.732 | 0.545 | 0.533 |
| K | ( )/h3/2 | 14.1 | 38.9 | 20.7 | 18.3 | |
| R | (°)/h2 | 106 | 735 | 256 | 219 |
4 结 论
TheoH方差是一种混合了Allan方差和去除偏移的Theo1方差的新的方差估计方法。把该方法用到动态信号随机误差分析中，可以改善DAVAR方法由于加窗截取数据，以致提取出的误差系数估计置信度降低的问题。
1) 在短相关时间上，DHVAR提取出的随机误差系数与Allan方差基本一致。中、长相关时间下，可以将随机误差系数的提取准确度提高30%以上。
2) DHVAR绘制的三维分析图更加清晰准确，可以观测到更为准确的动态信号的细节变化。
3) 对于速率斜坡这样需要较长相关时间才能辨识的噪声系数，动态方法的辨识结果置信度依然不高，需要进一步深入研究。
| [1] | 吕琳, 全伟. 基于GP+GA的陀螺仪随机误差建模分析[J]. 北京航空航天大学学报, 2015, 41 (6) : 1135 –1140. LV L, QUAN W. Modeling and analysis of gyroscope's random drift based on GA+GP method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41 (6) : 1135 –1140. (in Chinese) |
| [2] | 徐定杰, 苗志勇, 沈锋, 等. MEMS陀螺随机漂移误差系数的动态提取[J]. 宇航学报, 2015, 36 (2) : 217 –223. XU D J, MIAO Z Y, SHEN F, et al. Dynamic extraction MEMS gyro random error coefficients[J]. Journal of Astronautics, 2015, 36 (2) : 217 –223. (in Chinese) |
| [3] | ZHANG C X, WANG L, GAO S, et al. Dynamic Allan variance analysis for stochastic errors of fiber optic gyroscope[J]. Infrared and Laser Engineering, 2014, 43 (9) : 3081 –3088. |
| [4] | 王新龙, 李娜. MEMS陀螺随机误差的建模与分析[J]. 北京航空航天大学学报, 2012, 38 (2) : 170 –174. WANG X L, LI N. Error modeling and analysis for random drift of MEMS gyroscopes[J]. Journal of Beijing University of Aeronautics and Astronautics, 2012, 38 (2) : 170 –174. (in Chinese) |
| [5] | 汤霞清, 程旭维, 高军强. 光学陀螺随机误差特性的混合理论方差方法分析[J]. 兵工学报, 2015, 36 (9) : 1688 –1695. TANG X Q, CHENG X W, GAO J Q. Hybrid theoretical variance analysis for random error properties of optic gyroscope[J]. Acta Armamentarii, 2015, 36 (9) : 1688 –1695. (in Chinese) |
| [6] | GALLEANI L. Dynamic Allan variance III:Confidence and detection surfaces[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2011, 58 (8) : 1550 –1558. DOI:10.1109/TUFFC.2011.1982 |
| [7] | GALLEANI L, TAVELLA P. The dynamic Allan variance[J]. IEEE Transaction on Ultrasonics, Ferroelectrics,and Frequency Control, 2009, 56 (3) : 450 –460. DOI:10.1109/TUFFC.2009.1064 |
| [8] | LI J T, FANG J C. Not fully overlapping Allan variance and total variance for inertial sensor stochastic error analysis[J]. IEEE Transactions on Instrumentation and Measurement, 2013, 62 (10) : 2659 –2672. DOI:10.1109/TIM.2013.2258769 |
| [9] | HOWE D A. The total deviation approach to long-term characterization of frequency stability[J]. IEEE Transactions on Ultrasonic,Ferroelectrics,and Frequency Control, 2000, 47 (5) : 1102 –1110. DOI:10.1109/58.869040 |
| [10] | HOWE D A,TASSET T N.Theol:Characterization of very long-term frequency stability[C]//Proceedings of the 18th European Frequency and Time Forum (EFTF 04),Guildford,2004:581-587. |
| [11] | 程旭维, 汤霞清, 黄湘远. 基于#1理论方差的光学陀螺长期随机误差分析[J]. 中国激光, 2014, 41 (10) : 146 –153. CHENG X W, TANG X Q, HUANG X Y. Investigation on random error properties of optic gyroscope based on theoretical variance #1[J]. Chinese Journal of Lasers, 2014, 41 (10) : 146 –153. (in Chinese) |
| [12] | HOWE D A. ThêoH:A hybrid,high-confidence statistic that improves on the Allan deviation[J]. Metrologia, 2006, 43 (4) : 322 –331. DOI:10.1088/0026-1394/43/4/S20 |
| [13] | LEVINE J,PARKER T E.The algorithm used to realize UTC(NIST)[C]//2002 IEEE International Frequency Control Symposium. Piscataway,NJ:IEEE Press,2002:537-542. |
| [14] | 金毅, 吴训忠, 谢聂. 基于Allan方差的光纤陀螺随机漂移建模与仿真[J]. 应用光学, 2014, 35 (3) : 547 –551. JIN Y, WU X Z, XIE N. Modeling and simulation of FOG random drift based on Allan variance[J]. Journal of Applied Optics, 2014, 35 (3) : 547 –551. (in Chinese) |
| [15] | 韦官余,徐伯健,丁阳,等.动态阿伦方差辅助的卡尔曼滤波算法在GPS/INS组合导航中的应用[C]//第三届中国卫星导航学术年会电子文集,2012:330-334. WEI G Y,XU B J,DING Y,et al.Dynamic Allan variance aided Kalman filter in GPS/INS integrated navigation[C]//The Third China Satellite Navigation Conference,2012:330-334(in Chinese). |
| [16] | 汪立新, 朱战辉, 黄松涛. 基于峭度和自适应滑动窗的陀螺动态特性分析方法[J]. 中国惯性技术学报, 2015, 23 (4) : 533 –539. WANG L X, ZHU Z H, HUANG S T. Dynamic characteristics analysis method for gyroscope based on kurtosis and adaptive sliding window[J]. Journal of Chinese Inertial Technology, 2015, 23 (4) : 533 –539. (in Chinese) |