基于峭度和自适应滑动窗的陀螺动态特性分析方法

汪立新，朱战辉，黄松涛

（1. 第二炮兵工程大学，西安 710025；2. 中国人民解放军96401部队，宝鸡 721006）

基于峭度和自适应滑动窗的陀螺动态特性分析方法

汪立新1，朱战辉1，黄松涛2

（1. 第二炮兵工程大学，西安 710025；2. 中国人民解放军96401部队，宝鸡 721006）

针对运用动态Allan方差法进行陀螺随机误差分析时由于采用固定窗函数截取信号，导致信号跟踪效果与方差估计置信度不能同时兼顾的问题，提出了一种根据信号短时非平稳度自动调节窗宽的改进算法。首先运用截断窗内数据的峭度值表征陀螺输出的短时非平稳性，以随时间变化的峭度值为变量构造窗宽截取函数，再将陀螺量测信号中的随机误差分离出来，应用窗宽函数智能选取合适的窗长，用其来截取信号并按照一定的时间顺序分段计算Allan方差，最终将其绘制在时间、相关时间和Allan方差三维一体的图中。仿真和陀螺实测数据表明：新算法在保证动态跟踪效果不变的前提下，将中、低频噪声系数的辨识准确度提高了50%以上。

陀螺；峭度；自适应滑动窗；动态特性；动态Allan方差

惯性系统静动态性能的优劣直接影响到导弹、火箭及其它飞行器的工作精度，而陀螺又是惯性系统的核心部件。由于陀螺载体实际工作时会经受阵风、海浪的影响，内部发动机的振动、冲击（如头体分离、尾罩分离）、温度变化及电磁辐射等因素也会造成陀螺的输出表现出非平稳性，因而对其动态环境下随机误差进行辨识和补偿可以有效提高导航精度。Allan方差是一种时域分析技术，本质是将惯性传感器输出的随机误差信号输入到带通参数为τ的Allan方差滤波器，得到一组滤波输出，进而辨识出较多的误差项，从而对各种误差源及其噪声统计特性进行细致的表征和辨识[1]。但该方法只能用来分析理想的平稳信号，不能表征随机误差及噪声系数随时间的变化过程。动态Allan方差（DAVAR，dynamic Allan variance）是Allan方差的扩展，同样可以用来分析仪器的振荡稳定性[2-3]。近年来，该方法被引入到惯性器件的非平稳随机过程辨识中来，以期更好地分析动态条件下陀螺性能变化规律。 文献[4-5]分别将该方法应用到光学陀螺的动态性能评价上，对温度、振动及其它动态干扰下的陀螺性能变化进行了分析，有效地表征了惯性器件输出信号受环境影响所表现出的时变特性。文献[6]用按指数变化的相关时间τ替代线性变化的τ，提高了该算法的运算速度。文献[7]用动态Allan方差法追踪MEMS惯性测量单元随机误差变化，并用获取的零偏不稳定性误差值实时更新卡尔曼滤波算法中的过程协方差矩阵，通过实际数据证明新算法比传统参数固定的卡尔曼滤波算法具有更为优良的处理性能。然而由于动态Allan方差计算时采用固定长度的窗函数截取信号，导致参与方差估计的样本数据量大幅减小，长相关时间下方差估计置信度降低，使动态跟踪能力和方差估计值置信度提高二者之间形成了对立的关系。

DAVAR方法提出者Gallean指出：DAVAR方法中截断窗的长度和类型可以自由选择，尤其是窗口长度的选取对准确分析信号的动态特性起到了至关重要的作用。长窗宽截断窗可使方差计算样本包含更多数据，从而增加方差估计置信度，但同时长窗宽又降低了动态跟踪效果；短窗宽可以及时追踪信号突变，但方差估计的置信度又因此降低，通常情况下两者很难兼顾，需要经过多次试验，才能找到一个相对合适的窗长[8]。然而截断窗函数一旦选定，其长度和类型就无法再改变，算法缺乏灵活性。本文提出了一种用峭度表征截断窗内数据短时平稳性，实时监测信号平稳程度并自动调节截断窗长度的算法，既改善了信号的跟踪效果，又提高了方差估计值的置信度。通过仿真和对陀螺非平稳信号进行动态特征分析，验证了本文算法的有效性。

1 动态Allan方差和峭度

Allan方差法是一种基于时域的分析方法，是IEEE公认的陀螺参数分析的标准方法。不过该算法的缺点是计算获得的噪声系数是一组恒定的值，无法表征信号的非平稳特征。动态Allan方差法以一种更直观的方式将信号的时域稳定性和频域稳定性以三维的形式表征出来。DAVAR的具体实现步骤如下：

① 确定分析随机噪声信号x(t)时间起始点t1；

② 用中心点为t1，宽度为L(t1)的窗函数截断随机信号x(t)，获得窗口截断信号yT(t1)；

⑤ 选择另外一个窗口截断信号，即将窗口滑动到t2，重复步骤②③④，得到σy(t2,τ)，以此类推，可以得到Allan方差集合σy(tN,τ)，tn+1截断窗按照一定的步长与以tn为中心的截断窗相重叠；

峭度(Kurtosis)是反映信号分布特性的数值统计量，对信号标准差和幅值的变化格外敏感，不少学者已将其运用到动态信号特征提取中[9]。其表达式为

2 基于峭度的K-DAVAR算法设计

DAVAR法的本质就是带有滑动窗的Allan方差，目前的分析方法都是在窗函数类型和长度都固定的情况下进行，不能做到针对不同信号特征及时转变窗函数的长度及类型，造成了不必要的能量泄漏、置信度低以及不能及时追踪突变信号等缺点。

2.1 窗长函数的构造

Allan方差的估计是基于有限长度的数据，估计的可信度依赖于数据的独立组数。对于给定的随机序列，窗宽越短，窗宽内样本数据越少，方差估计的置信度就会越低。因而在理想的平稳随机过程中，各个噪声系数相对稳定，在较长时间内也不会有大的变化，显然窗长是越长越好，可以有效提高方差估计的置信度。然而在发生非平稳变化时，用较大的窗长计算方差又不可避免地把信号突变成份包含进去并平均掉，导致跟踪信号超前或延迟，甚至动态变化完全不能体现，因而此时选用短窗宽更为合适。对于非平稳随机序列来说，窗宽越短，跟踪动态变化的能力也越强，但数据划分的独立子集个数会变少，长相关时间下的方差估计置信度会降低。本文提出在信号平坦的时段选择长窗宽，在发生动态变化的时段自动调节为短窗宽的思想，兼顾动态信号跟踪能力和方差估计值的置信度。

改进算法的核心就是通过计算上一个截断窗内数据的峭度来预测下一个截断窗的长度。在峭度波动小，也就是相对平稳的数据段，采用长窗宽截断数据计算Allan方差，增加置信度；在峭度大，平坦度差的数据段，采用短窗宽截取数据，跟踪动态变化。结合惯性传感器非平稳信号随机误差的特点，本文提出了下面的窗长计算函数：

式中：λ1和λ2是常数，且λ2＞λ1，为自适应窗口长度的上下限，可以根据待分析信号的长度及信号平稳状况的先验信息来取值；阈值k1和k2为峭度值稳定区间的上下界（在k1和k2之间的峭度值被认为是稳定的）；ΔL为窗口调节的步长，可以决定窗口每次增加或者缩小的范围。从公式中可以看出，如果t时刻发生动态变化，其峭度值K(t)就会偏离k1和k2所包含的区间，t+1时刻窗口会自动变窄，偏离越多窗长缩小的越快；相反，如果t时刻计算得到的峭度值K(t)位于稳态时的峭度值k1和k2之间，窗口会自动变宽，逐步增加参与计算的样本个数。这样就可以保证窗长能够随着平稳性的变化自动伸长或缩短。

2.2 K-DAVAR改进算法实现

本文提出了以滑动窗内数据的峭度值为变量建立窗宽函数，智能截取数据进行Allan方差计算的K-DAVAR（Kurtosis-Dynamic Allan Variance）方法。改进算法设计如下：

① 确定随机信号x(t)分析时间起始点t1；

② 用中心点为t1，起始宽度为L(t1)的窗函数截断随机信号x(t)，获得窗口截断信号yT(t1)，支撑变量t′描述窗内渐逝的时间；

PL(t′)为长度L(t1)的矩形窗函数，其定义为

截取得到的信号为

③ 依据公式(1)计算截断信号yT(t1,t′)的峭度K(t1)，按照窗宽函数方程(2)确定t2段信号的截取长度L(t2)；

④ 将yT(t1,t′)同Allan窗口hτ(t′)做卷积建立增量过程Δ(t1,t′,τ)：

这时变量t′的范围变为

0≤τ≤τmax，通常τmax=L/3，则t1时刻估计值为

Allan方差可以被定义为上式的总体期望值

⑤ 此外，t1时刻陀螺的主要噪声系数可以通过最小二乘拟合分离出来，公式如下：

⑥ 将窗口滑动到t2，按照步骤③计算所得的窗宽L(t2)截断信号x(t)，获得窗口截断信号yT(t2)，即重复步骤②③④⑤，得到σy(t2,τ)，以此类推，可以得到时间域的方差序列σy(tN,τ)和噪声拟合系数A(tN)1,…,A(tN)5，其中以tn+1为中心的截断窗最好与以tn为中心的截断窗相重叠。公式(11)为离散信号下的DAVAR计算公式，其中τ0为采样周期，m为平均因子，1≤m≤N-1。

3 仿真试验与分析

为验证K-DAVAR方法对非平稳信号分析的有效性，用图1所示的有幅值突变的分段高斯白噪声x(t)来模拟陀螺随机误差信号，采样时间为1 s，共采样2.3小时。仿真信号3个平稳过程方差为1，两个突变段（1000~2000 s，3000~4000 s）的方差分别设为2和3，三个平稳段长度分别为1000 s和4000 s。

图1 陀螺输出仿真信号Fig.1 Simulate data of gyro output

对该信号分别用窗长为401、801的经典DAVAR法和K-DAVAR法进行分析，其中K-DAVAR法上下限1λ和2λ分别取401和801，ΔL取2，k取平稳条件下的峭度值3.15。图2是用K-DAVAR方法分析的三维图。

图2 仿真信号的K-DAVAR分析Fig.2 K-DAVAR analysis of simulate data

图3 峭度和窗宽变化过程Fig.3 Change process of kurtosis and window length

图3是在K-DAVAR算法分析中，峭度和窗长随时间的变化过程，可以看到峭度计算准确地捕获到了1000 s、2000 s、3000 s和4000 s的突变，算法在突变处逐渐缩小窗宽，离开突变点后窗宽又逐步变大，随着峭度的变化自适应的进行调节。

通过对Allan方差曲线进行最小二乘拟合，可以得到各项陀螺仪噪声系数，分别代表不同的误差来源[10]。以惯性传感器性能指标中常用到的角度随机游走（Rate random walk，简写为N）及零偏不稳定性（Bias instability，简写为B）为研究对象，经DAVAR和K-DAVAR二维计算，这两个噪声系数在不同窗长条件下量值随时间变化的趋势如图4所示。

图4 角度随机游走在不同窗长下的分析结果Fig.4 Analysis results of N with different length windows

图5 零偏不稳定性不同窗长下的分析结果Fig.5 Analysis results of B with different length windows

可以清晰地看到：DAVAR短窗宽401情况下由于参与计算的样本数据少，造成方差计算值波动很大，置信度较低，但对突变起始点和结束点的定位相对比较准确。长窗宽801情况下，样本数据较大，置信度高，方差计算值相对稳定平滑，但也因为窗口长的原因，截断窗较早的把突变信号包含了进去，又较晚地退了出来，以致跟踪时间上出现了较大的超前和滞后。而K-DAVAR方法对突变点的定位基本与DAVAR401方法相同，而噪声系数在平稳过程的波动程度又与DAVAR801方法基本相同，做到了两者兼顾。零偏不稳定性的动态跟踪能力比较如表1所示。

表1 K-DAVAR和DAVAR动态跟踪能力比较表Tab.1 Comparison on dynamic tracking capabilities

表2 稳定条件下噪声系数的估计值Tab.2 HRG noise coefficients under stationary condition

可以看到，在动态跟踪能力上，K-DAVAR方法与DAVAR401窗长基本一致，远远好于DAVAR801窗长。

理想的平稳随机过程下，动态算法拟合出的噪声系数时间变化曲线应该在Allan方差值（用全部数据估计出的数值，置信度很高）的附近波动[11-12]。但实际上因动态算法截断窗内数据量大幅减少造成功率泄漏以及长相关时间τ下估计值置信度降低等原因，噪声曲线的均值肯定要偏离标准值，我们要做的就是比较哪种算法拟合出的噪声曲线均值和标准值更为接近。本文将最后一个平稳段，也就是采样点4800到7800作为研究对象，综合比较两种算法效能。首先对该3000个数据计算Allan方差值，并拟合各个噪声系数值，此时样本数据量大，置信度相对较高，方差估计值可以作为标准值。再分别计算出DAVAR401窗长、801窗长和K-DAVAR方法下的噪声系数波动曲线，并对噪声序列取均值。

图6是零偏不稳定性在不同分析方法下4800~ 7800 s的变化过程，其中数值0.0612对应的红色虚线是用该3000个点计算的Allan方差零偏不稳定性标准值。可以看到DAVAR401窗长波动剧烈，且偏离标准值最多，而K-DAVAR方法和DAVAR801窗长基本保持在相同的估计水平。

图6 平稳条件下下噪声系数变化过程Fig.6 Noise coefficients under stationary condition

表2是4800 s到7800 s各个噪声系数时间序列均值同Allan方差标准值的比较表，可以看到K-DAVAR与标准值最为接近。长相关时间τ下拟合得到的系数误差特别大，是DAVAR方法的缺陷，因为长相关时间下拟合的噪声（如速率斜坡R）是低频噪声，必须要较长的数据才能得到准确结果，而DAVAR方法用截断窗把原始信号截成小段，从而造成能量泄漏，低频噪声系数拟合精度降低。但该方法对高频噪声的辨识还是比较准确的，而高频噪声系数如量化噪声Q、零偏不稳定性N和角度随机游走B又有更为广泛的应用价值。仿真结果充分验证了该算法既具有DAVAR短窗宽的动态跟踪能力，又具有DAVAR方法长窗宽的置信度，使两者之间的矛盾得到了完美的解决。

4 动态试验验证

对动态测试试验中采集到的某型号光学陀螺的量测信号进行预处理，主要包括奇点的剔除、趋势项和周期项的去除，得到随机误差信号如图7所示。信号采样周期为20 ms，采集数据约90000个，采样时间1800 s。运用自适应窗宽的K-DAVAR进行三维分析的结果如图8所示，可以看出，两个动态峰值和初始阶段较高的随机噪声都在三维图中得到了表征。

图7 陀螺随机误差输出Fig.7 Output of gyro’s random error

图8 陀螺随机误差的K-DAVAR分析Fig.8 K-DAVAR analysis of gyro’s random errors

图9 峭度和窗宽变化过程图Fig.9 Changing processes of Kurtosis and windows length

图9是陀螺峭度和滑动窗长随时间的变化过程。在峭度偏离稳定分布的地方，窗宽自动收缩，随着峭度回复平稳，窗宽也逐渐回到它的上限（即1601个采样点），实现了窗长伴随信号平稳程度自适应调节。

图10 角度随机游走在不同窗长下的分析结果Fig.10 Analysis results of rate random walk with different length windows

图11 零偏不稳定性在不同窗长下的分析结果Fig.11 Analysis results of bias instability with different length windows

零偏不稳定性（Bias instability，本文简称B）是光学陀螺的一个非常重要的随机误差分量，主要是由环境扰动和光学陀螺的残余非互异性引起，表现为输出数据中零偏值的波动。角度随机游走（Rate random walk，简写为N）也是影响陀螺性能的主要误差之一。图10和图11以这两个噪声系数为例来对比研究K-DAVAR方法的动态跟踪能力。特别值得注意的是：K-DAVAR方法甚至比DAVAR401窗长具有更为优秀的动态跟踪能力，因为其方差上下波动小，更容易定位动态发生点；而DAVAR1601窗长因其窗宽过长，使得第二个突变信号几乎没有跟踪到，信号细节也失去了很多，动态跟踪效果极差。

图12 平稳条件下噪声系数变化过程Fig.12 Noise coefficients under stationary condition

图12是专门截取零偏不稳定性B在静态平稳过程（信号1000~1600 s时间段），分别用自适应窗长的K-DAVAR方法和经典DAVAR方法窗长为401采样点和1601采样点辨识得到的噪声系数时间序列曲线。其中数值0.013是该段信号30000个数据的Allan方差计算值，理论上零偏不稳定性的时间变化曲线应在0.013附近上下波动。通过图12的比较可以看到，对于估计值的准确度来说，在稳定随机过程，K-DAVAR方法具有和1601窗长相同的估计置信度。陀螺五个随机噪声系数在该稳定段的曲线均值和Allan方差标准值的对比结果如表3。

表3 稳定条件下噪声系数的估计值Tab.3 Gyro noise coefficients under stationary condition

从表3中看到，K-DAVAR方法与DAVAR1601窗长方法在平稳过程中的估计值准确度基本一致，远远好于DAVAR401窗长，同时又具有了与DAVAR401窗长相同的跟踪能力。

5 结 论

K-DAVAR算法能有效地依据信号平稳程度不同自动调节滑动窗的窗长，解决了经典DAVAR方法动态跟踪能力和方差估计值置信度提高不能兼顾的问题。通过仿真和光学陀螺实验数据验证了算法的有效性。在动态跟踪能力不下降的前提下，对于用较短长度相关时间τ下就能准确辨识的噪声，也就是对于中、低频随机噪声（如角度随机游走和零偏不稳定性）来说，K-DAVAR方法的辨识准确度提高了50%以上。对于长相关时间τ下才能准确辨识的噪声系数，也就是部分低频噪声（如速率随机游走和速率斜坡）来说，虽然K-DAVAR算法的性能优于传统DAVAR方法，但基于截断窗的动态算法确实很难给出准确的辨识结果来，这也是动态算法不能完全代替Allan方差的重要原因。但是对于高、中频随机噪声也就是量化噪声、角度随机游走和零偏不稳定性，K-DAVAR方法具有比传统DAVAR方法准确更高，实时性更好的特点。用K-DAVAR方法对惯性传感器动态特性进行分析，对下一步的传感器结构改进，导航补偿算法和滤波算法精度提高都具有重要的应用价值。

（References）：

[1] 李杨, 胡柏青, 覃方君, 等. 光纤陀螺信号的解耦自适应Kalman滤波降噪方法[J]. 中国惯性技术学报, 2014, 22(2): 260-264. Li Y, Hu B Q, Qin F J, et al. De-noising method of decoupling adaptive Kalman filter for FOG signal[J]. Journal of Chinese Inertial Technology, 2014, 22(2): 260-264.

[2] Galleani L. The dynamic Allan variance III: Confidence and detection surfaces[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2011, 58(8): 1550-1558.

[3] Galleani L. The dynamic Allan variance IV: Characterization of atomic clock anomalies[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2015, 62(5): 791-801.

[4] 李冀辰, 高凤岐, 王广龙, 等. 光纤陀螺振动和变温条件下的DAVAR分析[J]. 中国激光, 2013, 40(9): 09080041-09080047 Li J C, Gao F Q, Wang G L, et al. Analysis of dynamic Allan variance for fiber optic gyro under vibration and variable temperature conditions[J]. Chinese Journal of Laser, 2013, 40(9): 09080041-09080047.

[5] 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.

[6] 张谦, 王玮, 王蕾, 等. 基于动态Allan 方差的光纤陀螺随机误差分析及算法改进[J]. 光学学报, 2015, 35(4): 0406003. Zhang Q, Wang W, Wang L, et al. Research on random errors of fiber optic gyro based on dynamic Allan variance and algorithm improvement[J]. Acta Optica Sinica, 2015, 35(4): 0406003.

[7] 韦官余, 徐伯健, 丁阳, 等. 动态阿伦方差辅助的卡尔曼滤波算法在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.

[8] Galleani L, Tavella P. The dynamic Allan variance[J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2009, 56(3): 450-460.

[9] Kohei H, Masato N, Takanobu N, et al. Close/distant talker discrimination based on kurtosis of linear prediction residual signals[C]//2014 IEEE International Conference on Acoustic, Speech and Signal Processing (ICASSP), 2014: 2327-2331.

[10] 白俊卿, 张科, 卫育新. 光纤陀螺随机漂移建模与分析[J]. 中国惯性技术学报, 2012, 20(5): 621-624. Bai J Q, Zhang K, Wei Y X. Modeling and analysis of fiber optic gyroscope random drifts[J]. Journal of Chinese Inertial Technology, 2012, 20(5): 621-624.

[11] 徐定杰, 苗志勇, 沈峰, 等. MEMS陀螺随机漂移误差系数的动态提取[J]. 宇航学报, 2015, 36(2): 217-223. Xu D J, Miao Z Y, Sheng F, et al. Dynamic extraction MEMS gyro random error coefficients[J]. Journal of Astronautics, 2015, 36(2): 217-223.

[12] Galleani L. The dynamic Allan variance II: A fast computation algorithm[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2010, 57(1): 182-187.

Analysis method for gyroscope dynamic characteristics based on kurtosis and adaptive sliding window

WANG Li-Xin1, ZHU Zhan-Hui1, HUANG Song-tao2
(1. The Second Artillery Engineering University, Xi’an 710025, China; 2. Unit 96401 of the Chinese People’s Liberation Army, Baoji 721006, China)

In view that the traditional dynamic Allan variance method in processing gyroscope’s random errors is difficult to make a good tradeoff between dynamic tracking capabilities and have a good confidence on the estimate due to the fixed length analyzing windows, an improved method was presented which can change length of truncation window automatically according to the instantaneous non-stationary of signal. Firstly, the kurtosis was introduced to describe the non-stationary process of signal, and the window length function which would be used to truncate signal was built by taking kurtosis as variables. Secondly, the random errors were separated from the outputs of gyroscopes. Then the length of the moving window was determined through the application of the widow length function, and the estimate of Allan variance was obtained with the data captured by the window. Simulation and dynamic experimental results show that the identification accuracy of low and intermediate frequency noise coefficients can be increased by 50% under the same dynamic tracking effect.

gyroscope; kurtosis; adaptive window; dynamic characteristics; dynamic Allan variance

V241.5

A

1005-6734(2015)04-0533-07

10.13695/j.cnki.12-1222/o3.2015.04.021

2015-04-12；

2015-07-29

二炮装备技术基础项目（EP114054）

汪立新（1966—），男，教授，博士生导师，从事惯性技术及测试方面的研究。E-mail：wanglixin066@sina.cn