譜估計(jì)理論在彈道數(shù)據(jù)參數(shù)化建模中的應(yīng)用

張 召，荊武興，李君龍，高長生

(1. 哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 150001；2. 中國航天科工集團(tuán)有限公司第二研究院，北京 100854)

0 引 言

憑借高速度和高機(jī)動(dòng)性能，高超聲速飛行器具備極強(qiáng)的突防能力，對現(xiàn)代防御體系構(gòu)成了巨大威脅。為對其進(jìn)行有效的跟蹤、預(yù)報(bào)和攔截，需要對其彈道特性進(jìn)行深入研究，了解其運(yùn)動(dòng)規(guī)律并進(jìn)行參數(shù)化建模。

目前針對彈道特性的研究，大都基于飛行器的動(dòng)力學(xué)模型，采用基于標(biāo)稱彈道的線性化[1-2]、分岔理論[3-4]和多尺度理論[5-6]等方法，研究制導(dǎo)參數(shù)對彈道數(shù)據(jù)的影響規(guī)律。然而，對于防御方來說，進(jìn)攻方飛行器的制導(dǎo)規(guī)律和動(dòng)力學(xué)模型是不可觀測的，而且在線識別的技術(shù)難度較大[7]。因此，本文將飛行過程視作黑箱模型，針對可觀測的彈道數(shù)據(jù)開展研究，分析其變化規(guī)律并給出參數(shù)化描述。相較于傳統(tǒng)研究方法，本文方法著重于挖掘彈道數(shù)據(jù)自身的內(nèi)在規(guī)律，可以為不同飛行器提供統(tǒng)一的彈道描述方式，而且易于利用跟蹤數(shù)據(jù)進(jìn)行在線建模，為實(shí)現(xiàn)彈道層面的匹配、識別和預(yù)報(bào)奠定基礎(chǔ)[8-9]。

針對數(shù)據(jù)內(nèi)在規(guī)律的挖掘，時(shí)間序列分析以及信號處理等領(lǐng)域都有豐富的研究成果。信號處理領(lǐng)域的譜估計(jì)理論，可以在輸入未知的情況下，對平穩(wěn)輸出序列進(jìn)行分析和建模。尤其是功率譜估計(jì)，是經(jīng)常被采用的一種重要方法，反映了信號功率隨頻率的分布。具體可分為經(jīng)典譜估計(jì)和現(xiàn)代譜估計(jì)[10]，其中現(xiàn)代譜估計(jì)的參數(shù)模型法可以給出信號的參數(shù)化描述。以參數(shù)模型法描述目標(biāo)運(yùn)動(dòng)已有諸多研究，文獻(xiàn)[11-13]使用自回歸(Auto-regressive，AR)模型描述人和機(jī)器人的運(yùn)動(dòng)規(guī)律并進(jìn)行預(yù)報(bào)；文獻(xiàn)[14-16]采用AR模型提取彈道目標(biāo)的進(jìn)動(dòng)特性；文獻(xiàn)[17-18]則在模型中引入了輸入量，研究競技體育中目標(biāo)的軌跡特性。數(shù)據(jù)平穩(wěn)是應(yīng)用譜估計(jì)理論的前提，但是高超聲速飛行彈道數(shù)據(jù)帶有明顯的趨勢，不符合該前提。本文采用線性消勢法消除彈道數(shù)據(jù)的線性趨勢，將其轉(zhuǎn)變?yōu)槠椒€(wěn)信號，進(jìn)而采用譜估計(jì)理論進(jìn)行分析和建模。

應(yīng)用參數(shù)模型的過程中，最重要的步驟之一就是模型階數(shù)的選擇。作為一般規(guī)律，如果模型階數(shù)選擇太低，將會得到一個(gè)高度平滑譜；如果選擇得太高，則可能在譜中引入虛假低峰[10]。為此，相關(guān)學(xué)者提出了不同的選擇準(zhǔn)則：F-檢驗(yàn)、最終預(yù)報(bào)誤差(Final prediction error，F(xiàn)PE)準(zhǔn)則、Akaike信息準(zhǔn)則(Akaike information criterion，AIC)、貝葉斯信息準(zhǔn)則(Bayesian information criterion，BIC)等。然而，一些試驗(yàn)結(jié)果表明，模型階數(shù)選擇準(zhǔn)則不能生成確定的結(jié)果[10]。為克服該問題，本文將經(jīng)典譜估計(jì)引入模型階數(shù)的選擇過程。以經(jīng)典譜估計(jì)生成的彈道數(shù)據(jù)譜圖為參考，對參數(shù)模型階數(shù)進(jìn)行調(diào)整，以完成模型階數(shù)的確定。

本文以高超聲速飛行彈道為對象，研究參數(shù)化建模方法。以典型彈道數(shù)據(jù)為基礎(chǔ)，使用線性消勢法將其轉(zhuǎn)變?yōu)槠椒€(wěn)信號；綜合使用經(jīng)典譜估計(jì)中的Welch直接法和現(xiàn)代譜估計(jì)中的AR參數(shù)模型法，對彈道數(shù)據(jù)進(jìn)行參數(shù)化分析和建模。文章最后以彈道地心距為例進(jìn)行了仿真校驗(yàn)，驗(yàn)證了本文算法所確定模型與動(dòng)力學(xué)模型的一致性。

1 高超聲速飛行彈道

進(jìn)行彈道數(shù)據(jù)譜估計(jì)的基礎(chǔ)是充足的數(shù)據(jù)支持，既可以是先驗(yàn)的彈道數(shù)據(jù)，也可以是實(shí)時(shí)的跟蹤數(shù)據(jù)。本節(jié)以動(dòng)力學(xué)模型和典型飛行模式為基礎(chǔ)，通過數(shù)值仿真方法獲得典型高超聲速飛行器的先驗(yàn)彈道，為后文分析其變化規(guī)律并實(shí)現(xiàn)參數(shù)化建模提供數(shù)據(jù)支持。高超聲速飛行在彈道系下表示為[19]：

(1)

(2)

(3)

(4)

(5)

(6)

典型的飛行模式包括常攻角飛行、常升阻比飛行、最大升阻比飛行、平衡滑翔飛行以及指標(biāo)最優(yōu)飛行[20]。其中指標(biāo)最優(yōu)包括：射程最優(yōu)、末速最大、氣動(dòng)加熱最少以及突防性能最優(yōu)等[21-24]。

2 基于Welch算法的彈道數(shù)據(jù)分析

將彈道數(shù)據(jù)視作信號，使用信號處理領(lǐng)域的譜估計(jì)理論對其進(jìn)行研究。在實(shí)際工作中，往往假定信號是平穩(wěn)的和各態(tài)遍歷的[10]。然而，高超聲速飛行彈道除了小范圍的波動(dòng)外，具有明顯的趨勢，直接進(jìn)行譜估計(jì)會導(dǎo)致較大的偏差。本節(jié)采用線性消勢法提前消除信號的趨勢，將其轉(zhuǎn)變?yōu)槠椒€(wěn)信號，再使用經(jīng)典譜估計(jì)理論生成信號的譜圖。

2.1 消除信號趨勢

消除信號趨勢的主要方法有線性消勢和差分消勢?？紤]到差分消勢會引入運(yùn)動(dòng)的高階信息而且高超聲速飛行彈道具有波動(dòng)特性，采用差分消勢會導(dǎo)致較高的模型階數(shù)，增大計(jì)算量甚至引入虛假低峰。本文采用線性消勢法，將數(shù)據(jù)信號的趨勢表述為線性形式：

d(n)=a+bn

(7)

結(jié)合生成的彈道數(shù)據(jù)，以擬合方式確定式(7)系數(shù)；然后從原始信號中消除相應(yīng)的線性部分，即可完成消勢。消除趨勢后的數(shù)據(jù)信號可視作平穩(wěn)的和各態(tài)遍歷的，可以開展進(jìn)一步的研究。

2.2 Welch算法

(8)

Welch算法很好地改善了周期圖法譜估計(jì)的方差特性，而且概念清晰、方法簡便。本文采用Welch算法對彈道數(shù)據(jù)進(jìn)行初步分析，給出大致的功率譜特性，為后文進(jìn)行更加精確的譜估計(jì)和參數(shù)化建模提供參考。

3 基于AR模型的彈道數(shù)據(jù)參數(shù)化建模

經(jīng)典譜估計(jì)可以給出大致準(zhǔn)確的譜圖，但是無法給出數(shù)據(jù)的參數(shù)化描述，不利于開展進(jìn)一步的研究?，F(xiàn)代譜估計(jì)的AR參數(shù)模型法不僅可以提高功率譜估計(jì)的分辨率，而且可以建立數(shù)據(jù)的參數(shù)化模型。然而，仍存在模型階數(shù)選擇準(zhǔn)則引起的不確定性問題。本節(jié)采用AR模型描述彈道數(shù)據(jù)的內(nèi)在規(guī)律，以改進(jìn)的協(xié)方差法為基礎(chǔ)對其進(jìn)行求解，并結(jié)合AIC、AR模型譜圖和Welch算法譜圖確定具體的模型階數(shù)。

3.1 平穩(wěn)信號的AR模型

假定所研究信號x(n)是由一個(gè)輸入信號u(n)激勵(lì)一個(gè)線性系統(tǒng)H(z)的輸出，則x(n)和u(n)之間有如下關(guān)系：

(9)

式中：若x(n)是確定性的，那么u(n)是一個(gè)沖激序列；若x(n)是隨機(jī)的，那么u(n)是一個(gè)白噪聲序列。

如果b1,b2,…,bq全為0，則有：

(10)

在Z域上，轉(zhuǎn)移函數(shù)為：

(11)

可見，上述模型僅有極點(diǎn)，這種全極點(diǎn)模型即為AR模型。其含義是，模型當(dāng)前輸出由當(dāng)前輸入和過去p個(gè)時(shí)刻的輸出決定。

本文中輸入序列是不可觀測的，但是如果輸出序列表現(xiàn)為平穩(wěn)隨機(jī)過程，那么輸入序列也被假定為平穩(wěn)隨機(jī)過程[10]。因此，假定u(n)是一個(gè)白噪聲序列，方差為σ2；由隨機(jī)信號通過線性系統(tǒng)的理論可知，x(n)的功率譜為：

(12)

確定式(12)中的白噪聲方差σ2及模型的系數(shù)a1,a2,…,ap，即可求出x(n)的功率譜。

3.2 AR模型系數(shù)確定

目前提出的有關(guān)AR模型系數(shù)的求解及AR模型性能的討論大都是建立在線性預(yù)測理論上的，而且這些算法的性能一般要優(yōu)于自相關(guān)法[10]。而且AR模型和線性預(yù)測器是等價(jià)的，AR模型的白噪聲能量σ2等于線性預(yù)測器的最小預(yù)測誤差功率ρmin。不失一般性，令前后預(yù)測誤差功率之和：

(13)

為最小。上標(biāo)f表示前向預(yù)測，上標(biāo)b表示后向預(yù)測。式(13)中：

(14)

(15)

在令ρfb為最小時(shí)，不是僅令ρfb相對am(m)=km為最小，而是令ρfb相對am(1),am(2),…,am(m)都為最小，m=1,…,p。

由ab(k)=af*(k)，令

(16)

(17)

即

(18)

又令

(19)

寫成矩陣形式：

(20)

最小預(yù)測誤差功率可由式(21)～(22)求出：

(21)

或者

(22)

式(22)即為改進(jìn)的協(xié)方差方法的正則方程，許多學(xué)者給出了相應(yīng)的求解算法[10]。

3.3 AR模型階數(shù)選擇

為解決上述問題，本文將經(jīng)典譜估計(jì)引入AR模型階數(shù)的選擇過程：首先，模型階數(shù)由1逐步增加，計(jì)算AIC：

(23)

本節(jié)所設(shè)計(jì)的模型階數(shù)選擇算法，綜合考慮了時(shí)域指標(biāo)和頻域指標(biāo)，分別從時(shí)域和頻域兩個(gè)角度分析所建模型的精度。其中，時(shí)域指標(biāo)由改進(jìn)的協(xié)方差方法求解得到，描述了模型誤差的統(tǒng)計(jì)特性；頻域指標(biāo)由Welch算法的估計(jì)結(jié)果提供參考，用以評估AR模型譜圖的質(zhì)量。兩個(gè)指標(biāo)的綜合運(yùn)用，可以有效改善僅有時(shí)域指標(biāo)造成的不確定性問題，提高建模精度和可靠性。

4 仿真校驗(yàn)

以HTV-2最大升阻比飛行彈道的地心距為例進(jìn)行仿真分析，初始高度45 km，速度6 km/s，當(dāng)?shù)厮俣葍A角0°，飛行時(shí)長1178 s，采樣周期0.1 s，地心距如圖1所示?？梢姷匦木喑尸F(xiàn)波動(dòng)下降，并且振蕩幅值逐步衰減。下面將地心距數(shù)據(jù)視作信號，分析其功率譜并建立參數(shù)化模型，最后對所建立的模型進(jìn)行驗(yàn)證。

圖1 地心距信號時(shí)域圖Fig.1 Time domain diagram of geocentric distance signal

使用Welch算法分析地心距信號的功率譜，如圖2所示?？梢娫?.02667 rad/s處存在一個(gè)轉(zhuǎn)折頻率，而且在低頻部分存在更大的能量。對照地心距信號時(shí)域圖可知，地心距為非平穩(wěn)信號，呈現(xiàn)波動(dòng)減小趨勢：小范圍波動(dòng)即為轉(zhuǎn)折頻率所代表的部分，減小趨勢即是低頻部分。該非平穩(wěn)信號具有較大的頻域跨度，不利于對信號進(jìn)行精確分析和建模，需要消除信號的趨勢，也就是消除其低頻部分。

圖2 Welch算法確定地心距信號的譜圖Fig.2 Spectrum of geocentric distance signal via Welch algorithm

采用線性消勢法消除地心距信號的趨勢，經(jīng)擬合確定的趨勢表述為(單位：m)：

d(n)=6.4235×106-5.8370n

(24)

消除地心距信號中相應(yīng)的線性項(xiàng)，得到線性消勢地心距信號，如圖3所示?？梢娦盘栐?值附近振蕩，且為衰減振蕩，整體仍存在小幅度的非線性趨勢。

圖3 線性消勢地心距信號時(shí)域圖Fig.3 Time domain diagram of linearly detrended geocentric distance signal

使用Welch算法分析線性消勢地心距信號的功率譜，如圖4所示。可見第一峰值頻率為0.032 rad/s、第二峰值頻率為0.02667 rad/s，該頻段集中了大部分能量，在低頻部分存在一個(gè)小峰值。對照信號時(shí)域圖可知，衰減振蕩即為第一和第二峰值所代表的部分，整體小幅度非線性趨勢即是低頻峰值所代表的部分。

圖4 Welch算法確定線性消勢地心距信號的譜圖Fig.4 Spectrum of linearly detrended geocentric distance signal via Welch algorithm

使用AR模型對線性消勢地心距信號進(jìn)行功率譜估計(jì)。首先，模型階數(shù)由1逐步增加并計(jì)算AIC，確定AIC最優(yōu)階數(shù)為3，則階數(shù)選擇范圍定為2，3，4。下面分別使用這三種AR模型進(jìn)行功率譜估計(jì)，如圖5～7所示，可知線性消勢地心距信號的峰值頻率分別為0.02769 rad/s、0.02894 rad/s、0.01001 rad/s和0.02944 rad/s。

圖5 2階AR模型確定線性消勢地心距信號的譜圖Fig.5 Spectrum of linearly detrended geocentric distance signal via second-order AR model

圖6 3階AR模型確定線性消勢地心距信號的譜圖Fig.6 Spectrum of linearly detrended geocentric distance signal via third-order AR model

圖7 4階AR模型確定線性消勢地心距信號的譜圖Fig.7 Spectrum of linearly detrended geocentric distance signal via fourth-order AR model

結(jié)合時(shí)域圖以及Welch算法的譜估計(jì)結(jié)果可知，線性消勢地心距信號的峰值頻率在0.02667～0.032 rad/s之間，而且在低頻部分存在一個(gè)小峰值。2階和3階AR模型均識別出該范圍內(nèi)的一個(gè)峰值，4階AR模型則識別出該范圍內(nèi)一個(gè)峰值和范圍外一個(gè)低頻峰值。3階AR模型雖然是AIC最優(yōu)的，但其低頻部分仍存在較大的能量，因此其對信號的描述是不準(zhǔn)確的，其輸出存在低頻漂移。2階和4階模型識別出窄尖峰譜線，對信號的刻畫較為準(zhǔn)確。識別得到的線性消勢地心距AR模型分別為(單位：m)：

A(z)x(n)=u(n)

(25)

(26)

則完備的地心距模型為經(jīng)線性補(bǔ)償?shù)腁R模型(Linearly compensated AR model，LAR)：

r(n)=x(n)+d(n)=x(n)+

6.4235×106-5.8370n

(27)

下面對上述模型進(jìn)行校驗(yàn)，基本思想為：將動(dòng)力學(xué)模型作為參考，考察在相同初始條件下LAR模型輸出與動(dòng)力學(xué)模型輸出的一致性程度；其中LAR模型的輸出稱為LAR數(shù)據(jù)，動(dòng)力學(xué)模型的輸出稱為參考數(shù)據(jù)。驗(yàn)證指標(biāo)為：以均值和標(biāo)準(zhǔn)差檢驗(yàn)兩樣本的基本性能參數(shù)是否一致；采用秩和檢驗(yàn)法校驗(yàn)兩樣本概率分布的一致性[25]。

首先，由動(dòng)力學(xué)模型生成地心距的真實(shí)數(shù)據(jù)，并在相同初始條件下使用上述LAR模型對地心距進(jìn)行預(yù)報(bào)，結(jié)果如圖8～9所示?？梢?，1階模型給出的預(yù)報(bào)結(jié)果是線性的，預(yù)報(bào)誤差隨著彈道的振蕩而波動(dòng)；2階模型、3階模型和4階模型預(yù)報(bào)結(jié)果都是波動(dòng)形式，頻率和振幅與動(dòng)力學(xué)模型給出的真實(shí)值一致。

然后，考察LAR數(shù)據(jù)與參考數(shù)據(jù)的一致性程度，如表1所示。秩和檢驗(yàn)中，0假設(shè)為兩樣本概率分布一致，1假設(shè)為不一致?？梢?階和3階模型的均值偏離參考值較大，未通過秩和檢驗(yàn)，這是輸出的低頻漂移造成的；2階和4階模型通過了秩和檢驗(yàn)，與動(dòng)力學(xué)模型具有較高的一致性，驗(yàn)證了本文方案的有效性。

圖8 地心距預(yù)報(bào)結(jié)果Fig.8 Prediction results of geocentric distance

圖9 地心距預(yù)報(bào)誤差Fig.9 Prediction error of geocentric distance

表1 模型校驗(yàn)結(jié)果Table 1 Model validation results

5 結(jié) 論

本文綜合運(yùn)用信號處理領(lǐng)域的經(jīng)典和現(xiàn)代譜估計(jì)理論，對彈道數(shù)據(jù)進(jìn)行分析和建模。由Welch算法確定大致的譜圖，然后由AR模型進(jìn)行精確的功率譜估計(jì)，并給出彈道數(shù)據(jù)的參數(shù)化模型。針對彈道數(shù)據(jù)非平穩(wěn)特性，本文采用線性消勢法對信號進(jìn)行消除趨勢操作；針對模型階數(shù)選擇準(zhǔn)則無法給出確定結(jié)果的問題，本文提出綜合經(jīng)典和現(xiàn)代譜估計(jì)的譜圖，確定AR模型的階數(shù)。仿真結(jié)果表明，針對以線性趨勢為主要趨勢并帶有小幅波動(dòng)的高超聲速跳躍彈道，以本文方法確定的LAR模型與動(dòng)力學(xué)模型具有較高的一致性。