王哲，李國輝，趙善祿

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

0 引言

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

1 最大熵譜估計(jì)與一致性驗(yàn)證

1.1 最大熵譜估計(jì)

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

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

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

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

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

1.2 一致性驗(yàn)證

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

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

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

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

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

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

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

在檢驗(yàn)過程中，原則上對(duì)每個(gè)頻率點(diǎn)都要進(jìn)行檢驗(yàn)，如果每個(gè)頻率點(diǎn)都接受該檢驗(yàn)，則說明兩序列相容。檢驗(yàn)的公式為：

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

2 實(shí)例驗(yàn)證

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

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

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

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

采用Marple算法解yule-walker方程，利用MATLAB軟件計(jì)算出相關(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式對(duì)仿真數(shù)據(jù)與試飛數(shù)據(jù)進(jìn)行最大熵譜估計(jì)，得到如圖1所示功率譜對(duì)比曲線。

圖1nx隨Vi變化最大熵譜估計(jì)對(duì)比曲線

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

其中，M1,M2為AR模型階次，N1,N2為采樣點(diǎn)數(shù)量，α為顯著水平，統(tǒng)計(jì)量D對(duì)于該假設(shè)檢驗(yàn)的接受域?yàn)閇 ]

-Zα/2,Zα/2，按公式 11進(jìn)行檢驗(yàn)：

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

圖2nx隨Vi變化最大熵譜估計(jì)驗(yàn)證曲線

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

3 結(jié)語

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

[1]李鵬波.時(shí)間序列樣本的總體一致性檢驗(yàn)——頻域方法[J].飛行器測(cè)控學(xué)報(bào)，1999，18（4）.

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

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

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

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

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

[7]王建華,符文星,董敏周.最大墑譜估計(jì)在空空導(dǎo)彈仿真模型驗(yàn)證中的應(yīng)用[J].彈箭與制導(dǎo)學(xué)報(bào),2005,25（4）.