王哲，李國輝，趙善祿

（空軍航空大學軍事仿真研究所，長春130022）

0 引言

經(jīng)典譜估計雖然計算效率高，但是由于窗函數(shù)的局限性，在計算過程中不可避免地出現(xiàn)頻譜泄露、譜峰轉(zhuǎn)移、分辨率低等缺點。針對經(jīng)典譜估計的局限，現(xiàn)代譜估計方法應運而生，其中最大熵譜估計法是一種常用的現(xiàn)代譜估計方法。

1 最大熵譜估計與一致性驗證

1.1 最大熵譜估計

1967年，Burg最早將“熵”這一概念從信息論中引入到譜估計中，提出最大熵譜估計[1]。最大熵譜估計具有分辨率高，使用范圍廣的特點。

“熵”是對于隨機性進行度量的量，熵值越大，則隨機性越大。假設p(x)為連續(xù)隨機變量的概率函數(shù)密度，則“熵”可以定義為：

如果平穩(wěn)時間序列{xt,t∈Z}是正態(tài)隨即過程，其功率譜為S(ω)，那么功率譜熵為：

與熵的原始定義相似，I（S）表示隨即序列不確定性的度量，I（S）值越大，表明隨機性越強。根據(jù)譜估計知識可知，樣本序列自相關(guān)函數(shù)和功率譜密度構(gòu)成傅里葉變換對，在已知相關(guān)函數(shù)r(m),m=0,1,2,…,p情況下，對m>p的自相關(guān)函數(shù)進行外推，當功率譜S(ω)滿足如下約束條件：

當在上數(shù)條件下使得H(S(ω))最大時，這樣的過程稱為最大熵譜估計，估計出的S?()ω為最大熵譜?？梢宰C明，正態(tài)隨機過程的最大熵譜估計與AR模型譜估計是等效的[2]。AR模型的功率譜為：要想求信號功率譜密度pn(ω)，需要求出上式中的參數(shù)φ和σ。可以通過Yule-Walker法或Burg法求得相關(guān)參數(shù)，具體過程見參考文獻[3]。當確定相關(guān)參數(shù)后，就可通過4式計算功率譜。

1.2 一致性驗證

對最大熵譜估計進行一致性檢驗，有常規(guī)的非參數(shù)檢驗、窗譜估計一致性檢驗和最大熵譜估計一致性檢驗等多種方法[4]。本文采用最大熵譜估計一致性檢驗法與模型相結(jié)合進行檢驗。

最大熵譜估計的統(tǒng)計結(jié)果是有限的近似結(jié)果，是漸進無偏與漸進正態(tài)的[5]。當采樣數(shù)目N與AR模型的階數(shù)M足夠大時，最大熵譜的均值與方差為：

對 S?(ω )作對數(shù)變換得到logS?(ω )，由于 S?(ω)是正態(tài)的，所以 logS?(ω )仍為正態(tài)的。logS?(ω)的均值和方差為：

該對數(shù)譜 logS?(ω)的100(1 -α)%置信區(qū)間近似為

對于兩個獨立時間序列的功率譜密度均漸進正態(tài)情況下，第i個帶寬中對數(shù)的抽樣分布近似為：

假設兩個時間序列具有相同功率譜密度函數(shù)，即S1(ω1)=S2(ω2)，根據(jù)式8可以得出：

該假設檢驗的接受域是[ ]-Zα/2≤D≤Zα/2，α 為該檢驗的顯著水平，N1、N2是兩次試驗的樣本容量，M1、M2是AR模型的階次。

在檢驗過程中，原則上對每個頻率點都要進行檢驗，如果每個頻率點都接受該檢驗，則說明兩序列相容。檢驗的公式為：

但是，與窗譜估計一致性檢驗類似，在實際驗證過程中，由于頻率點過多并且有的點功率譜較弱，所以工程中常常關(guān)注于頻率譜較強的點。

2 實例驗證

由前面章節(jié)的知識可知，對于正態(tài)隨機過程，可以通過建立與之等效的AR模型進行最大熵譜估計。本節(jié)仍采用飛機水平加速的縱向過載nx隨飛行速度ViV變化的仿真數(shù)據(jù)與試飛數(shù)據(jù)進行一致性檢驗。故需先建立AR模型。

在對仿真數(shù)據(jù)與試飛數(shù)據(jù)進行最大熵譜估計前，首先需要建立AR模型。由于經(jīng)過數(shù)據(jù)預處理，仿真數(shù)據(jù)與試飛數(shù)據(jù)均可視為正態(tài)隨機過程，可以參考文獻[6]中時序建模的方法建立AR模型。對于模型階數(shù)的確定，目前有多種準則與方法，主要有赤池FPE準則、AIC準則和Parzen的自回歸傳遞函數(shù)準則[7]。本文對于模型階數(shù)的確定，采用工程上常用的經(jīng)驗公式由于采樣數(shù)目為N1,N2=3000，則M≈18，采樣頻率為

階數(shù)為 P，模型參數(shù)為 φ1,φ2,…,φP的 AR（p）可以表示為

其中，{εt}是白噪聲序列，滿足 EXsεt=0 對一切s

采用Marple算法解yule-walker方程，利用MATLAB軟件計算出相關(guān)參數(shù)，得到仿真數(shù)據(jù)的AR模型為：

Xt=1.5761Xt-1-0.7562Xt-2-0.5123Xt-3-0.3165Xt-4+0.2119Xt-5+0.0352Xt-6-0.0158Xt-7-0.0264Xt-8+0.0322Xt-9+0.0286Xt-10-0.0346Xt-11+0.0153Xt-12-0.0234Xt-13+0.0315Xt-14+0.0446Xt-15-0.0061Xt-16+0.0019Xt-17+0.0211Xt-18

試飛數(shù)據(jù)的AR模型為：

Xt=1.3213Xt-1-0.8613Xt-2-0.7449Xt-3-0.4881Xt-4-0.3591Xt-5-0.1011Xt-6-0.0695Xt-7+0.0987Xt-8+0.0413Xt-9-0.0363Xt-10-0.0303Xt-11-0.0134Xt-12+0.0006Xt-13+0.0201Xt-14+0.0362Xt-15-0.0061Xt-16+0.0109Xt-17+0.0211Xt-18

根據(jù)4式對仿真數(shù)據(jù)與試飛數(shù)據(jù)進行最大熵譜估計，得到如圖1所示功率譜對比曲線。

圖1nx隨Vi變化最大熵譜估計對比曲線

根據(jù)歸一化處理后的功率譜密度，確定主要頻率區(qū)間為[0,0.025]，因此可以對該區(qū)間內(nèi)的頻率點進行一致性檢驗。要驗證nx隨Vi變化的試飛數(shù)據(jù)功率譜S1(ω1)和仿真數(shù)據(jù)功率譜S2(ω2)的相容性，需要檢驗S1(ω1)=S2(ω2)，構(gòu)造統(tǒng)計量D：

其中，M1,M2為AR模型階次，N1,N2為采樣點數(shù)量，α為顯著水平，統(tǒng)計量D對于該假設檢驗的接受域為[ ]

-Zα/2,Zα/2，按公式 11進行檢驗：

當顯著水平 α取 0.05時，查表得 Zα/2=1.96，N1=N2=3000，M1=M2=18，計算得B≈0.30。因此，只需要驗證在主要頻率區(qū)間內(nèi)每個頻率點處A的值小于0.30，就可以說明試飛數(shù)據(jù)與仿真數(shù)據(jù)的功率譜密度具有一致性。圖2為A在[0,0.025]區(qū)間內(nèi)的曲線圖。

圖2nx隨Vi變化最大熵譜估計驗證曲線

由圖2可知，在主要頻率區(qū)間內(nèi)，有A≤B，所以在置信水平為95%情況下，S1(ω1)=S2(ω2)成立，即 nx隨Vi變化的試飛數(shù)據(jù)與仿真數(shù)據(jù)具有一致性。

3 結(jié)語

本文以飛機在某高度平飛加速過程中過載隨速度的變化為例，采用最大熵譜估計法對仿真數(shù)據(jù)和試飛數(shù)據(jù)進行一致性檢驗，取譜峰附近[0,0.025dB]作為主要頻率區(qū)間，檢驗結(jié)果顯示，功率譜的主要頻率點在95%的置信水平下，nx隨Vi變化的仿真數(shù)據(jù)與試飛數(shù)據(jù)在統(tǒng)計意義下一致，該方法適用于其他動態(tài)飛行性能的一致性檢驗。

[1]李鵬波.時間序列樣本的總體一致性檢驗——頻域方法[J].飛行器測控學報，1999，18（4）.

[2]曾鳴,李雪青,謝保川,等.基于飛參數(shù)據(jù)的飛行仿真模型驗證[J].指揮控制與仿真,2011,6（12）.

[3]徐慧娟.自回歸AR模型的整體最小二乘分析研究[D].上海:東華理工大學,2012.

[4]李鶴.基于試飛數(shù)據(jù)的模型飛行模擬器飛行性能驗證研究[D].空軍工程大學,2009.

[5]李鵬波,高霞.應用最大熵譜估計進行導彈系統(tǒng)的仿真模型驗證[J].國防科技大學學報，1999,21（2）.

[6]李鶴,呂巖,李國輝.飛行仿真動態(tài)數(shù)據(jù)驗證的時序建模方法[J].兵工自動化,2008,27（12）.

[7]王建華,符文星,董敏周.最大墑譜估計在空空導彈仿真模型驗證中的應用[J].彈箭與制導學報,2005,25（4）.