基于峭度和自適應滑動窗的陀螺動態特性分析方法

汪立新，朱戰輝，黃松濤

（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