王哲,李國輝,趙善祿
(空軍航空大學軍事仿真研究所,長春130022)
經典譜估計雖然計算效率高,但是由于窗函數的局限性,在計算過程中不可避免地出現頻譜泄露、譜峰轉移、分辨率低等缺點。針對經典譜估計的局限,現代譜估計方法應運而生,其中最大熵譜估計法是一種常用的現代譜估計方法。
1967年,Burg最早將“熵”這一概念從信息論中引入到譜估計中,提出最大熵譜估計[1]。最大熵譜估計具有分辨率高,使用范圍廣的特點。
“熵”是對于隨機性進行度量的量,熵值越大,則隨機性越大。假設p(x)為連續(xù)隨機變量的概率函數密度,則“熵”可以定義為:
如果平穩(wěn)時間序列{xt,t∈Z}是正態(tài)隨即過程,其功率譜為S(ω),那么功率譜熵為:
與熵的原始定義相似,I(S)表示隨即序列不確定性的度量,I(S)值越大,表明隨機性越強。根據譜估計知識可知,樣本序列自相關函數和功率譜密度構成傅里葉變換對,在已知相關函數r(m),m=0,1,2,…,p情況下,對m>p的自相關函數進行外推,當功率譜S(ω)滿足如下約束條件:
當在上數條件下使得H(S(ω))最大時,這樣的過程稱為最大熵譜估計,估計出的S?()ω為最大熵譜。可以證明,正態(tài)隨機過程的最大熵譜估計與AR模型譜估計是等效的[2]。AR模型的功率譜為:要想求信號功率譜密度pn(ω),需要求出上式中的參數φ和σ??梢酝ㄟ^Yule-Walker法或Burg法求得相關參數,具體過程見參考文獻[3]。當確定相關參數后,就可通過4式計算功率譜。
對最大熵譜估計進行一致性檢驗,有常規(guī)的非參數檢驗、窗譜估計一致性檢驗和最大熵譜估計一致性檢驗等多種方法[4]。本文采用最大熵譜估計一致性檢驗法與模型相結合進行檢驗。
最大熵譜估計的統(tǒng)計結果是有限的近似結果,是漸進無偏與漸進正態(tài)的[5]。當采樣數目N與AR模型的階數M足夠大時,最大熵譜的均值與方差為:
對 S?(ω )作對數變換得到logS?(ω ),由于 S?(ω)是正態(tài)的,所以 logS?(ω )仍為正態(tài)的。logS?(ω)的均值和方差為:
該對數譜 logS?(ω)的100(1 -α)%置信區(qū)間近似為
對于兩個獨立時間序列的功率譜密度均漸進正態(tài)情況下,第i個帶寬中對數的抽樣分布近似為:
假設兩個時間序列具有相同功率譜密度函數,即S1(ω1)=S2(ω2),根據式8可以得出:
該假設檢驗的接受域是[ ]-Zα/2≤D≤Zα/2,α 為該檢驗的顯著水平,N1、N2是兩次試驗的樣本容量,M1、M2是AR模型的階次。
在檢驗過程中,原則上對每個頻率點都要進行檢驗,如果每個頻率點都接受該檢驗,則說明兩序列相容。檢驗的公式為:
但是,與窗譜估計一致性檢驗類似,在實際驗證過程中,由于頻率點過多并且有的點功率譜較弱,所以工程中常常關注于頻率譜較強的點。
由前面章節(jié)的知識可知,對于正態(tài)隨機過程,可以通過建立與之等效的AR模型進行最大熵譜估計。本節(jié)仍采用飛機水平加速的縱向過載nx隨飛行速度ViV變化的仿真數據與試飛數據進行一致性檢驗。故需先建立AR模型。
在對仿真數據與試飛數據進行最大熵譜估計前,首先需要建立AR模型。由于經過數據預處理,仿真數據與試飛數據均可視為正態(tài)隨機過程,可以參考文獻[6]中時序建模的方法建立AR模型。對于模型階數的確定,目前有多種準則與方法,主要有赤池FPE準則、AIC準則和Parzen的自回歸傳遞函數準則[7]。本文對于模型階數的確定,采用工程上常用的經驗公式由于采樣數目為N1,N2=3000,則M≈18,采樣頻率為
階數為 P,模型參數為 φ1,φ2,…,φP的 AR(p)可以表示為
其中,{εt}是白噪聲序列,滿足 EXsεt=0 對一切s 采用Marple算法解yule-walker方程,利用MATLAB軟件計算出相關參數,得到仿真數據的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 試飛數據的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 根據4式對仿真數據與試飛數據進行最大熵譜估計,得到如圖1所示功率譜對比曲線。 圖1nx隨Vi變化最大熵譜估計對比曲線 根據歸一化處理后的功率譜密度,確定主要頻率區(qū)間為[0,0.025],因此可以對該區(qū)間內的頻率點進行一致性檢驗。要驗證nx隨Vi變化的試飛數據功率譜S1(ω1)和仿真數據功率譜S2(ω2)的相容性,需要檢驗S1(ω1)=S2(ω2),構造統(tǒng)計量D: 其中,M1,M2為AR模型階次,N1,N2為采樣點數量,α為顯著水平,統(tǒng)計量D對于該假設檢驗的接受域為[ ] -Zα/2,Zα/2,按公式 11進行檢驗: 當顯著水平 α取 0.05時,查表得 Zα/2=1.96,N1=N2=3000,M1=M2=18,計算得B≈0.30。因此,只需要驗證在主要頻率區(qū)間內每個頻率點處A的值小于0.30,就可以說明試飛數據與仿真數據的功率譜密度具有一致性。圖2為A在[0,0.025]區(qū)間內的曲線圖。 圖2nx隨Vi變化最大熵譜估計驗證曲線 由圖2可知,在主要頻率區(qū)間內,有A≤B,所以在置信水平為95%情況下,S1(ω1)=S2(ω2)成立,即 nx隨Vi變化的試飛數據與仿真數據具有一致性。 本文以飛機在某高度平飛加速過程中過載隨速度的變化為例,采用最大熵譜估計法對仿真數據和試飛數據進行一致性檢驗,取譜峰附近[0,0.025dB]作為主要頻率區(qū)間,檢驗結果顯示,功率譜的主要頻率點在95%的置信水平下,nx隨Vi變化的仿真數據與試飛數據在統(tǒng)計意義下一致,該方法適用于其他動態(tài)飛行性能的一致性檢驗。 [1]李鵬波.時間序列樣本的總體一致性檢驗——頻域方法[J].飛行器測控學報,1999,18(4). [2]曾鳴,李雪青,謝保川,等.基于飛參數據的飛行仿真模型驗證[J].指揮控制與仿真,2011,6(12). [3]徐慧娟.自回歸AR模型的整體最小二乘分析研究[D].上海:東華理工大學,2012. [4]李鶴.基于試飛數據的模型飛行模擬器飛行性能驗證研究[D].空軍工程大學,2009. [5]李鵬波,高霞.應用最大熵譜估計進行導彈系統(tǒng)的仿真模型驗證[J].國防科技大學學報,1999,21(2). [6]李鶴,呂巖,李國輝.飛行仿真動態(tài)數據驗證的時序建模方法[J].兵工自動化,2008,27(12). [7]王建華,符文星,董敏周.最大墑譜估計在空空導彈仿真模型驗證中的應用[J].彈箭與制導學報,2005,25(4).3 結語