張秋昭，張書畢* ，侯東陽，劉萬利

(1.中國礦業(yè)大學(xué)國土環(huán)境與災(zāi)害監(jiān)測國家測繪局重點(diǎn)實(shí)驗(yàn)室，江蘇徐州 221116;2.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221116)

在慣性技術(shù)中，常把慣性測量單元(Inertial Measurement Unit，IMU)誤差分為系統(tǒng)性誤差和隨機(jī)性誤差，其中系統(tǒng)性誤差可以通過測試和標(biāo)定進(jìn)行補(bǔ)償;而隨機(jī)性誤差隨時(shí)間變化，成為影響IMU精度的關(guān)鍵因素。常用的隨機(jī)誤差辨識(shí)方法主要包括功率譜密度(PSD)［1］、Allan 方差［2－3］、時(shí)間序列法［5－7］和自相關(guān)函數(shù)估計(jì)［6］法。

譜密度刻畫了一個(gè)平穩(wěn)時(shí)間序列的功率譜，并特征化其二階矩的性質(zhì)。通過研究估計(jì)后的譜密度，可以識(shí)別引起數(shù)據(jù)最大變化的頻率范圍，因此可以用信號(hào)的功率譜密度來表征信號(hào)的頻域特性，可以應(yīng)用于辨識(shí)慣性傳感器的隨機(jī)誤差［1，7］。Allan方差是20世紀(jì)60年代提出的一種方法，最初被用于研究高精度振蕩器的頻率穩(wěn)定性，后來用于辨識(shí)慣性傳感器隨機(jī)誤差，由于其可以辨識(shí)出較多的誤差項(xiàng)，得到了大家的廣泛重視［2－3］。自相關(guān)函數(shù)法由于其辨識(shí)能力和精度有限［6］，應(yīng)用不多。

PSD和Allan方差法都需要用近似擬合曲線得到隨機(jī)模型，但這些估計(jì)是兩次統(tǒng)計(jì)建模，含有較大統(tǒng)計(jì)誤差［8］，而時(shí)間序列分析方法(AR模型等)只需做一次統(tǒng)計(jì)分析，建模精度高，并可以得到等效的連續(xù)狀態(tài)方程表達(dá)式，常用來作為激光陀螺等高精度慣導(dǎo)系統(tǒng)的建模方法［4］。但AR模型常要求序列為穩(wěn)態(tài)隨機(jī)過程，而中低精度的光纖陀螺(Fiber Optic Gyro，F(xiàn)OG)常含有趨勢(shì)項(xiàng)或周期性誤差，必須對(duì)隨機(jī)序列進(jìn)行去平穩(wěn)化處理［8］。

目前常見的平穩(wěn)化處理的主要方法［9］:多項(xiàng)式回歸、小波分析等。多項(xiàng)式回歸由于原理簡單易于實(shí)現(xiàn)，在平穩(wěn)化處理中得到了廣泛應(yīng)用。常見的多項(xiàng)式擬合多假設(shè)序列含有現(xiàn)行趨勢(shì)項(xiàng)或二次項(xiàng)，多項(xiàng)式的階數(shù)如何取，取的的階數(shù)太少，擬合精度不高，太高會(huì)產(chǎn)生過擬合現(xiàn)象［9］;小波分析可以提取序列的趨勢(shì)項(xiàng)［10］，但小波基、閾值、分解的層數(shù)的選擇受主觀性影響比較嚴(yán)重。

經(jīng)驗(yàn)?zāi)B(tài)分解方法(Empirical Mode Decomposition，EMD)是一種自適應(yīng)的局部時(shí)頻分析方法［11］，它根據(jù)信號(hào)自身的特性將信號(hào)自適應(yīng)地分解成有限個(gè)經(jīng)驗(yàn)?zāi)B(tài)函數(shù)，特別適合非平穩(wěn)信號(hào)分析，在信號(hào)去噪、趨勢(shì)項(xiàng)剔除等方面得到了廣泛應(yīng)用［11－12］。黨淑雯等引入基于經(jīng)驗(yàn)?zāi)B(tài)分解的提升小波消噪技術(shù)，去除光纖陀螺中的1/fγ類型分形噪聲［13］。曲從善等采用一種基于經(jīng)驗(yàn)?zāi)B(tài)分解的濾波方法對(duì)激光陀螺隨機(jī)信號(hào)進(jìn)行消噪［14］。論文引入經(jīng)驗(yàn)?zāi)B(tài)分解，以低精度光纖陀螺KVH-CPT為實(shí)驗(yàn)對(duì)象，對(duì)陀螺輸出的序列進(jìn)行平穩(wěn)化處理，采用逆序法驗(yàn)證平穩(wěn)化效果，與小波分析等方法的對(duì)比實(shí)驗(yàn)結(jié)果表明，EMD方法能有效提前陀螺隨機(jī)序列的趨勢(shì)項(xiàng)，實(shí)現(xiàn)陀螺序列的平穩(wěn)化處理。

1 AR模型

描述隨機(jī)序列{xt}某一時(shí)刻t和前p個(gè)時(shí)刻序列值之間的相互關(guān)系表示為

若隨機(jī)系列{εt}是白噪聲且和前一時(shí)刻序列xk(k＜t)不相關(guān)，則這樣的模型稱為p階自回歸模型，記為AR(p)。

AR(p)模型描述的是這樣的一種平穩(wěn)隨機(jī)過程:該隨機(jī)過程在t時(shí)刻的觀測值xt與t時(shí)刻之前的觀測 εt－1，εt－2，…εt－p存在相關(guān)性，即與過去p個(gè)時(shí)刻的觀測值存在相關(guān)性。因此，AR模型常用來描述陀螺的隨機(jī)噪聲［7，9］。

2 平穩(wěn)化檢驗(yàn)

AR模型描述的是隨機(jī)平穩(wěn)序列的方法，因此在用AR對(duì)隨機(jī)系列進(jìn)行建模時(shí)必須首先對(duì)隨行序列進(jìn)行平穩(wěn)性檢驗(yàn)。

平穩(wěn)性檢驗(yàn)的主要目的是檢驗(yàn)陀螺隨機(jī)誤差時(shí)間系列是否具有不隨時(shí)間原點(diǎn)的推移而變化的統(tǒng)計(jì)特性。如果隨機(jī)誤差系列是平穩(wěn)的，再假設(shè)序列服從正態(tài)分布，則對(duì)陀螺隨機(jī)誤差的研究，就可以采用單個(gè)樣本的時(shí)間序列代替總體平均，這就給數(shù)據(jù)處理帶來很大方便。如果不平穩(wěn)，則需要對(duì)數(shù)據(jù)進(jìn)行平穩(wěn)后處理［8－9］。

造成隨機(jī)序列不平穩(wěn)的主要原因是其中包含有隨時(shí)間緩慢變化的趨勢(shì)項(xiàng)。檢驗(yàn)這種非平穩(wěn)趨勢(shì)項(xiàng)的一個(gè)很有效的方法是逆序法［9］，簡述如下:

(1)對(duì)于測試數(shù)據(jù)x1，x2，…，xn，將其分為M小段，然后求各小段的平均值，得到一個(gè)大致不相關(guān)的均值序列y1，y2，…，yM。

(2)對(duì)于下標(biāo)為i的yi，每當(dāng)出現(xiàn)yj＞yi(j＞i，i=1，2，…，M－1)時(shí)就將其定義為yi的一個(gè)逆序，與yi對(duì)應(yīng)的逆序的個(gè)數(shù)Ai稱為yi的逆序數(shù)。

(3)記序列y1，y2，…，yM的逆序總數(shù)定義為:

(4)記隨機(jī)整數(shù)序列A的均值和方差分別為:

如果λ值處于±2之內(nèi)，則可接受序列無趨勢(shì)的假設(shè)(在0.05顯著水平上);否則拒絕該假設(shè)，認(rèn)為序列為非平穩(wěn)序列［9］。

3 經(jīng)驗(yàn)?zāi)B(tài)分解原理

經(jīng)驗(yàn)?zāi)B(tài)分解將信號(hào)分解為滿足以下條件的模態(tài)函數(shù)(Intrinsic Mode Function，IMF):零點(diǎn)數(shù)目相同與極值點(diǎn)數(shù)目相同或至多相差1;函數(shù)關(guān)于局部平均對(duì)稱;最終將信號(hào)表示成若干IMF與單調(diào)殘差函數(shù)之和:

其中，imfi(t)表示i個(gè)IMF，rn(t)是單調(diào)殘差函數(shù)。EMD算法采用所謂的“篩”完成，其基本步驟為:

①求取信號(hào)x(t)的極值點(diǎn);

②分別由極小值點(diǎn)與極大值點(diǎn)求得包絡(luò)線emin(t)和emax(t);

③計(jì)算上下包絡(luò)的平均值m(t)=(emax(t)+emin(t))/2;

④提取細(xì)節(jié)d(t)=x(t)－m(t);

⑤對(duì)細(xì)節(jié)d(t)重復(fù)進(jìn)行①至④步驟，直至d(t)為一個(gè)零均值，此時(shí)d(t)即為一個(gè)IMF;

⑥計(jì)算殘差r(t)=x(t)－imfi(t);

⑦重復(fù)運(yùn)算直至殘差不含有IMF。

應(yīng)用EMD進(jìn)行多尺度信號(hào)分解時(shí)，定義如下按尺度標(biāo)準(zhǔn)化模量的累積均值［12］(Mean of the Accumulated Standardized Modes，MASM):

其中m≤nimfi(t)是第i尺度的模量。如果明顯偏離零均值，則認(rèn)為從尺度m開始是信號(hào)的趨勢(shì)變化所致。去趨勢(shì)項(xiàng)后的信號(hào)可以表示為:

4 實(shí)驗(yàn)分析

為了驗(yàn)證論文的方法，以KVH－CPT型號(hào)開環(huán)光纖陀螺為實(shí)驗(yàn)對(duì)象。靜態(tài)采集了該光纖陀螺的原始觀測值，處理后得到實(shí)測陀螺隨機(jī)信號(hào)序列，如圖1所示。對(duì)其采用逆序法進(jìn)行平穩(wěn)性檢驗(yàn)，求得統(tǒng)計(jì)參數(shù) λ=－5.2613，由前所述 λ?［－2，2］，拒絕假設(shè)，認(rèn)為該序列為非平穩(wěn)序列。采用EMD提取趨勢(shì)項(xiàng)，結(jié)果如圖2所示，當(dāng)分層至13層時(shí)，其累積均值明顯偏離0值，如圖3所示?？梢哉J(rèn)為剩余的信號(hào)為趨勢(shì)項(xiàng)見圖2(b)。對(duì)去除趨勢(shì)項(xiàng)重構(gòu)后的序列采用逆序法求取統(tǒng)計(jì)量 λ=0.4935，λ∈［－2，2］，則接受檢驗(yàn)，認(rèn)為該序列為平穩(wěn)序列。

圖1 陀螺信號(hào)去偏后的序列

圖2 數(shù)據(jù)一EMD分解圖及分解尺度與累加均值的關(guān)系

采用偏態(tài)系數(shù)和峰態(tài)系數(shù)對(duì)其進(jìn)行正態(tài)性檢驗(yàn)。偏態(tài)系數(shù)(Skew Value)ξ和峰態(tài)系數(shù)(Kurtosis Value)υ的定義為:

其中μxσx分別為序列xt的均值和中誤差。偏態(tài)系數(shù)ξ和峰態(tài)系數(shù)υ的理想取值為0和3［9］。

同時(shí)采用Matlab中函數(shù)庫中的去趨勢(shì)項(xiàng)函數(shù)dtrend和小波分析方法對(duì)該序列進(jìn)行處理。經(jīng)過多次實(shí)驗(yàn)，發(fā)現(xiàn)當(dāng)小波分析采用‘db4’小波基和分解層數(shù)為20層時(shí)效果較好。比較得到的結(jié)果如表1所示。

采用EMD對(duì)隨機(jī)序列處理前后結(jié)果的數(shù)值比較見表1。

表1 數(shù)據(jù)一幾種方法處理前后結(jié)果比較

為了使論文采用的方法更具說服力，采用上述方法對(duì)另一種數(shù)據(jù)進(jìn)行同樣的處理，結(jié)果如圖3和表2所示。其中小波分析方法中采用‘db4’小波基，分解層數(shù)為18層。

圖3 數(shù)據(jù)二EMD分解圖及分解尺度與累加均值的關(guān)系

表2 數(shù)據(jù)二幾種方法處理前后結(jié)果比較

由表1和表2可以看出，dtrend函數(shù)對(duì)隨機(jī)序列的趨勢(shì)項(xiàng)提取效果較差，小波分析方法和EMD方法，均能分辨出原始隨機(jī)序列的趨勢(shì)項(xiàng)，均能明顯改善原隨機(jī)序列的平穩(wěn)性，正態(tài)性也略有改善。但小波分析方法需要多次實(shí)驗(yàn)才能確定最佳的小波基和分解層數(shù)，而EMD方法能自動(dòng)判斷分解層數(shù)，效果較好。

5 結(jié)論

AR模型等時(shí)間序列分析方法常用來對(duì)光纖陀螺隨機(jī)誤差建模，但其只能對(duì)平穩(wěn)序列進(jìn)行建模和處理。論文以低精度光纖陀螺為實(shí)驗(yàn)對(duì)象，采用經(jīng)驗(yàn)?zāi)B(tài)分解對(duì)非平穩(wěn)信號(hào)序列進(jìn)行平穩(wěn)化處理，并與小波分析等方法進(jìn)行了比較實(shí)驗(yàn)。結(jié)果表明:該方法能有效辨識(shí)非平穩(wěn)信號(hào)序列中的趨勢(shì)項(xiàng)，平穩(wěn)化后的信號(hào)可以采用時(shí)間序列分析方法進(jìn)行建模。

［1］Yi Yudan.On Improving the Accuracy and Reliability of GPS/INSBased DirectSensorGeoreferencing.［D］.TheOhioState University Ph.D Dissertation.2007.

［2］Naser El-Sheimy，Hou Haiying，Niu Xiaoji.Analysis and Modeling of Inertial Sensors Using Allan Variance［J］.IEEE Transactions on Instrumentation and Measurement，2008，57(1).:140－149.

［3］熊凱，雷擁軍，曾海波.基于Allan方差法的光纖陀螺建模與仿真［J］.空間控制技術(shù)與應(yīng)用，2010，36(3):7－13.

［4］Tang Jianghe，F(xiàn)u Zhenxian，Deng Zhenglong.Indntification Method for RLG Random Errors Based on Allan Variance and Equivalent Theorem［J］.Chinese Journal of Aeronautics.2009，(22):273－278.

［5］方靖，尚捷，顧啟泰.光纖陀螺隨機(jī)誤差建模的實(shí)驗(yàn)研究［J］.傳感技術(shù)學(xué)報(bào)，2008，21(9):1514－1518.

［6］Nassar 2003 Modeling Inertial Sensor Errors Using Autorgressive Models［C］//Proceedings of the Institute of Navigation(ION)National Technical Meeting(NTM 2003)，2003:116－225.

［7］吳富梅.GNSS/INS組合導(dǎo)航誤差補(bǔ)償與自適應(yīng)濾波理論的拓展［D］.解放軍信息工程大學(xué)博士論文.2010.

［8］朱奎寶，張春熹，張小躍.光纖陀螺隨機(jī)漂移ARMA模型研究［J］.宇航學(xué)報(bào)，2006，27(5):1118－1121.

［9］馮慶奇.激光陀螺捷聯(lián)慣性導(dǎo)航系統(tǒng)組合導(dǎo)航及零速修正技術(shù)研究［D］.國防科學(xué)技術(shù)大學(xué)碩士論文.2009.

［10］周哲，莊良杰，熊正南，等.基于小波分析的陀螺漂移趨勢(shì)項(xiàng)提?。跩］.中國慣性技術(shù)學(xué)報(bào)，1999，7(4):58－61.

［11］Huang N E，Shen Z，Long S R.The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series［J］.Proceedings of the Royal Society A，1998，(454):903－995.

［12］王堅(jiān)，高井祥，王金玲.基于經(jīng)驗(yàn)?zāi)B(tài)分解的GPS基線解算模型［J］.測繪學(xué)報(bào)，2008，37(1):10－14.

［13］曲從善，于鴻，許化龍，等.基于經(jīng)驗(yàn)?zāi)B(tài)分解的激光陀螺隨機(jī)信號(hào)消噪［J］.紅外與激光工程，2009，38(5):859－863.

［14］黨淑雯，田蔚風(fēng)，金志華.基于經(jīng)驗(yàn)?zāi)B(tài)分解的光纖陀螺1/fγ類型分形噪聲濾除方法［J］.紅外，2009，30(10):12－17.