赫泰龍， 陳萬(wàn)春,*， 劉芳

(1. 北京航空航天大學(xué) 宇航學(xué)院, 北京 100083； 2. 北京航天長(zhǎng)征飛行器研究所， 北京 100076)

高超聲速飛行器，相比于傳統(tǒng)彈道式飛行器，具有遠(yuǎn)距離快速打擊、覆蓋范圍廣、全過(guò)程機(jī)動(dòng)能力強(qiáng)、彈道難以被預(yù)報(bào)等優(yōu)點(diǎn)，被認(rèn)為是突破目前防空反導(dǎo)系統(tǒng)的有效手段，在軍事上有廣闊的應(yīng)用前景。高超聲速飛行器無(wú)動(dòng)力滑翔段彈道設(shè)計(jì)是當(dāng)前研究熱點(diǎn)方向之一，平衡滑翔是再入的一種重要飛行模式，最早由德國(guó)科學(xué)家S?nger和Bredt[1]提出，平衡滑翔彈道高度變化平緩、熱流密度和動(dòng)壓峰值小，被廣泛應(yīng)用于再入彈道的分析與設(shè)計(jì)[1-2]。胡錦川和陳萬(wàn)春[3]在傳統(tǒng)平衡滑翔概念的基礎(chǔ)上，通過(guò)分析給定攻角曲線(xiàn)和傾側(cè)角曲線(xiàn)再入彈道族的特點(diǎn)，以縱向加速度導(dǎo)數(shù)的平方的積分來(lái)衡量彈道的平滑程度，提出了平穩(wěn)滑翔彈道概念，并對(duì)平穩(wěn)滑翔彈道的解析解和彈道動(dòng)態(tài)特性進(jìn)行了分析，獲得了快速生成平穩(wěn)滑翔彈道的方法[4]。平穩(wěn)滑翔彈道是本文主要研究對(duì)象。

對(duì)于高超聲速滑翔彈道擾動(dòng)運(yùn)動(dòng)的分析，通常采用非線(xiàn)性仿真方法[3-4]，針對(duì)不同的干擾因素和不同參數(shù)(剩余飛行時(shí)間等)，需要分別進(jìn)行多次仿真；而對(duì)于隨機(jī)擾動(dòng)，則需要大量的蒙特卡羅仿真，計(jì)算效率低，但是方法簡(jiǎn)單，適用條件廣。在小擾動(dòng)假設(shè)下，對(duì)平穩(wěn)滑翔動(dòng)態(tài)方程線(xiàn)性化，利用伴隨仿真方法分析，可以大大減少計(jì)算量。

伴隨仿真方法是一種計(jì)算機(jī)仿真分析與設(shè)計(jì)工具，主要用于分析線(xiàn)性時(shí)變系統(tǒng)[5-6]。伴隨仿真方法的理論基礎(chǔ)是線(xiàn)性系統(tǒng)與其伴隨系統(tǒng)之間的本質(zhì)關(guān)系(內(nèi)積等式)[7]，一般的線(xiàn)性時(shí)變系統(tǒng)可以采用2種方法來(lái)描述：狀態(tài)空間法和脈沖響應(yīng)矩陣。所以，關(guān)于伴隨仿真方法的解釋也可以分為基于狀態(tài)空間描述(狀態(tài)轉(zhuǎn)移矩陣)[8-9]和基于脈沖響應(yīng)矩陣[5-6]。伴隨系統(tǒng)一次仿真可以得到原線(xiàn)性系統(tǒng)各不同干擾輸入對(duì)某一輸出的影響，通常稱(chēng)為誤差預(yù)算或靈敏度分析，在電路網(wǎng)絡(luò)分析中有廣泛應(yīng)用。伴隨仿真另一個(gè)重要性質(zhì)是一次仿真可以得到原線(xiàn)性系統(tǒng)在各干擾輸入不同的加入時(shí)間下某一特定時(shí)間(如仿真終端時(shí)刻)的輸出，該性質(zhì)廣泛應(yīng)用于導(dǎo)彈制導(dǎo)回路的分析與設(shè)計(jì)中[6]。伴隨仿真方法還可以用于分析輸入為隨機(jī)過(guò)程(如白噪聲信號(hào))的線(xiàn)性系統(tǒng)的均方響應(yīng)，本質(zhì)上可以理解為線(xiàn)性系統(tǒng)協(xié)方差分析[5-6]的伴隨仿真方法。Zarchan[10]基于伴隨仿真方法，提出了用于分析非線(xiàn)性導(dǎo)彈制導(dǎo)系統(tǒng)的統(tǒng)計(jì)線(xiàn)性化伴隨仿真方法(SLAM)。Weiss和Bucco[11]用伴隨仿真方法對(duì)戰(zhàn)術(shù)導(dǎo)彈的交班進(jìn)行了分析，并對(duì)如何在非標(biāo)準(zhǔn)輸入信號(hào)的制導(dǎo)系統(tǒng)應(yīng)用伴隨仿真方法進(jìn)行了研究[12]。Bucco[13]還系統(tǒng)地介紹了伴隨仿真方法在航空航天領(lǐng)域的應(yīng)用。Sarachik和Kreindler[14]、Brogan[15]、Zarchan[6]、Weiss和Bucco[16]分別研究了離散系統(tǒng)、連續(xù)-離散(混成)系統(tǒng)的伴隨仿真方法。林曉輝、崔乃剛和劉暾[17]對(duì)伴隨理論有關(guān)定理用于隨機(jī)運(yùn)動(dòng)體仿真技術(shù)作了研究。鄒暉、陳萬(wàn)春和邢曉嵐[18]基于MATRIXx開(kāi)發(fā)出了能夠自動(dòng)生成原系統(tǒng)的伴隨系統(tǒng)的工具軟件。伴隨仿真方法在攝動(dòng)制導(dǎo)方法的分析與設(shè)計(jì)中也有重要應(yīng)用[19-20]。

本文利用伴隨仿真方法，研究狀態(tài)初值和氣動(dòng)力等存在擾動(dòng)因素(包括確定性常值小擾動(dòng)和隨機(jī)擾動(dòng))時(shí)，對(duì)平穩(wěn)滑翔彈道終端高度、速度和射程的影響。首先，從伴隨系統(tǒng)基本定義出發(fā)，給出了伴隨仿真方法統(tǒng)一的解釋。在滑翔動(dòng)力學(xué)建模和平穩(wěn)滑翔彈道定義基礎(chǔ)上，研究了平穩(wěn)滑翔彈道的基本性質(zhì)。考慮擾動(dòng)因素，建立帶有誤差干擾的非線(xiàn)性動(dòng)力學(xué)模型，并根據(jù)小偏差彈道攝動(dòng)理論，在標(biāo)準(zhǔn)平穩(wěn)滑翔彈道附近線(xiàn)性化，得到關(guān)于擾動(dòng)運(yùn)動(dòng)的線(xiàn)性化微分方程。然后，采用伴隨仿真方法，同時(shí)對(duì)比非線(xiàn)性仿真和蒙特卡羅法，分析了各擾動(dòng)因素對(duì)平穩(wěn)滑翔彈道的影響；針對(duì)平穩(wěn)滑翔彈道存在初始高度、初始彈道傾角、初始速度及升力、阻力等多個(gè)擾動(dòng)因素，伴隨仿真方法只需要一次仿真就可以得到各擾動(dòng)在不同剩余飛行時(shí)間加入時(shí)對(duì)終端高度、終端速度或終端射程的影響，而非線(xiàn)性仿真需要分別對(duì)每個(gè)擾動(dòng)進(jìn)行仿真、對(duì)加入擾動(dòng)時(shí)剩余飛行時(shí)間進(jìn)行循環(huán)，而且對(duì)于隨機(jī)擾動(dòng)則需要進(jìn)行大量蒙特卡羅打靶，所以伴隨仿真方法的計(jì)算量要遠(yuǎn)小于非線(xiàn)性仿真或蒙特卡羅仿真。研究結(jié)果有利于加深對(duì)高超聲速飛行器平穩(wěn)滑翔彈道特性的認(rèn)識(shí)，并能為基于平穩(wěn)滑翔彈道的制導(dǎo)方法設(shè)計(jì)分析及仿真評(píng)估提供參考。

1 伴隨仿真方法

為了闡述伴隨仿真方法的原理，本節(jié)從線(xiàn)性時(shí)變系統(tǒng)的伴隨系統(tǒng)的定義出發(fā)，利用線(xiàn)性系統(tǒng)與其伴隨系統(tǒng)之間的本質(zhì)關(guān)系(定義式)，給出伴隨仿真方法的統(tǒng)一解釋。

〈Gu,ν〉Y=〈u,G*ν〉U?u∈U,?ν∈Y

(1)

式中：〈·,·〉為空間U或Y上的內(nèi)積。

1.1 伴隨系統(tǒng)狀態(tài)空間表示

考慮如下?tīng)顟B(tài)空間描述的線(xiàn)性時(shí)變系統(tǒng)：

(2)

式中：u(t)∈Rm為控制輸入；x(t)∈Rn為狀態(tài)向量；y(t)∈Rp為輸出；A(t)、B(t)和C(t)為關(guān)于時(shí)間t連續(xù)的函數(shù)矩陣。

系統(tǒng)式(2)的解可以表示為

(3)

式中：Φ(t,τ)為狀態(tài)轉(zhuǎn)移矩陣。

假設(shè)u∈U，y∈Y，信號(hào)空間U和Y都是有限時(shí)域勒貝格-2空間；為了考慮初值，將系統(tǒng)式(2)寫(xiě)成如下抽象形式：

(4)

式中：⊕為直和。

定義內(nèi)積為

(5)

容易驗(yàn)證，Rn⊕U和Rn⊕Y均是希爾伯特空間。利用線(xiàn)性系統(tǒng)伴隨的定義式(1)，可以推導(dǎo)出線(xiàn)性系統(tǒng)式(2)的伴隨系統(tǒng)。

(6)

存在如下?tīng)顟B(tài)空間實(shí)現(xiàn)[21]：

(7)

為了驗(yàn)證系統(tǒng)式(7)是式(2)的伴隨的狀態(tài)空間實(shí)現(xiàn)，進(jìn)行如下推導(dǎo)運(yùn)算：

(-pT(t)A(t)-rT(t)C(t))x(t)+

pT(t)(A(t)x(t)+B(t)u(t))=

-rT(t)C(t)x(t)+pT(t)B(t)u(t)=

-rT(t)y(t)+uT(t)q(t)

對(duì)上式從t0到tf積分可得到

pT(tf)x(tf)-pT(t0)x(t0)=

(8)

式(8)對(duì)任意x(t0)∈Rn，u∈U，p(tf)∈Rn，r∈Y 都成立。對(duì)式(8)左右兩邊進(jìn)行簡(jiǎn)單移項(xiàng)，并結(jié)合內(nèi)積定義式(5)，可以看出式(8)正是伴隨系統(tǒng)定義中式(1)的狀態(tài)空間變量表示，這說(shuō)明系統(tǒng)式(7)和式(2)互為各自伴隨的狀態(tài)空間實(shí)現(xiàn)。式(8)通常稱(chēng)為Bliss公式[19]或伴隨定理[17]。

注意到，伴隨系統(tǒng)式(7)的時(shí)間進(jìn)程是反向的(從tf到t0)，為了仿真方便，作如下變量替換：

(9)

為了書(shū)寫(xiě)方便，不失一般性，后文將假設(shè)t0=0。

利用式(7)進(jìn)行簡(jiǎn)單求導(dǎo)運(yùn)算和變量代換，可得到

(10)

式(10)的解為

ν(τ)dτ

(11)

式中：Ψ(t,τ)=ΦT(tf-τ,tf-t)為系統(tǒng)式(10)的狀態(tài)轉(zhuǎn)移矩陣。如無(wú)特別說(shuō)明，線(xiàn)性系統(tǒng)式(2)的伴隨系統(tǒng)是指式(10)，相應(yīng)地，伴隨系統(tǒng)定義式(8)則寫(xiě)為

zT(0)x(tf)-zT(tf)x(0)=

(12)

1.2 伴隨仿真統(tǒng)一解釋

利用伴隨的定義式(12)，能夠以更簡(jiǎn)單、統(tǒng)一的方式來(lái)解釋伴隨系統(tǒng)式(10)的仿真結(jié)果與原系統(tǒng)式(2)之間的關(guān)系。

注意到，式(12)只顯含原系統(tǒng)和伴隨系統(tǒng)的初始狀態(tài)、終端狀態(tài)、控制輸入和輸出，且由伴隨定義可知式(12)對(duì)任意的初始狀態(tài)x(0)、z(0)，輸入u(t)、ν(t)都成立。利用式(12)來(lái)解釋伴隨仿真方法的基本思路是：對(duì)原系統(tǒng)和伴隨系統(tǒng)選取適當(dāng)或特殊的狀態(tài)初值和輸入，討論原系統(tǒng)和伴隨系統(tǒng)仿真結(jié)果(終端狀態(tài)、輸出)之間的關(guān)系。對(duì)于具體待研究的線(xiàn)性系統(tǒng)，狀態(tài)初值和輸入都給定，則只需要對(duì)伴隨系統(tǒng)選取適當(dāng)或特殊的狀態(tài)初值和輸入。

(13)

式中：f(t)為任意函數(shù)。

伴隨仿真的基本解釋可以概括為如下2條性質(zhì)。記x=(xk)∈Rn，z=(zk)∈Rn，u=(uk)∈Rm，ν=(νk)∈Rp，y=(yk)∈Rp，w=(wk)∈Rm。

性質(zhì)1(誤差預(yù)算) 假設(shè)原系統(tǒng)式(2)為多輸入-單輸出(MISO)，給定任意控制輸入u(t)和狀態(tài)初值x(0)，則可以選取伴隨系統(tǒng)式(10)的輸入ν(t)=δ(t)，初值z(mì)(0)=0；代入式(12)，可得到

(14)

還可以選取伴隨系統(tǒng)ν(t)=0，z(0)=CT(tf)，注意到原系統(tǒng)式(2)為MISO，所以z(0)和CT(tf)是同維數(shù)的，利用原系統(tǒng)輸出方程y(t)=C(t)x(t)和式(12)，同樣可以得到式(14)。

特別地，如果原系統(tǒng)式(2)的u(t)=0，式(14)可寫(xiě)為

(15)

如果原系統(tǒng)式(2)中u(t)的每個(gè)分量都是δ(t)，即uk(t)=δ(t)，k=1,2,…,m，則式(14)可記為

(16)

式(14)表明原系統(tǒng)式(2)在終端時(shí)刻的輸出y(tf)，可以利用伴隨系統(tǒng)式(10)的仿真結(jié)果得到；而且伴隨一次仿真可以得到原系統(tǒng)輸入和狀態(tài)初值各分量對(duì)終端輸出y(tf)的貢獻(xiàn)大小，這一點(diǎn)可以從式(15)、式(16)更清楚地看到。性質(zhì)1也給出伴隨仿真初值和輸入的具體選取方法，特別地，對(duì)于原系統(tǒng)輸入為典型信號(hào)的[12]，可先等價(jià)轉(zhuǎn)化為零輸入、階躍輸入或脈沖輸入的線(xiàn)性系統(tǒng)[6]，然后對(duì)此系統(tǒng)進(jìn)行伴隨仿真，更方便操作和理解。

性質(zhì)2(伴隨一次仿真的通用解釋) 為了簡(jiǎn)化敘述，不失一般性，這里假設(shè)原系統(tǒng)式(2)為單輸入-單輸出(SISO)，且狀態(tài)初值x(0)=0，控制輸入u(t)=δ(t-(tf-tgo))，tgo∈[0,tf]；選取伴隨系統(tǒng)式(10)的ν(t)=δ(t)，z(0)=0或ν(t)=0，z(0)=CT(tf)，代入式(12)可得到

y(tf)=w(tgo)

(17)

式(17)對(duì)所有tgo∈[0,tf]成立，這說(shuō)明原系統(tǒng)在仿真時(shí)刻為tf-tgo，或者說(shuō)在剩余仿真時(shí)間為tgo時(shí)，加入脈沖輸入，得到終端輸出y(tf)，等于伴隨仿真在tgo時(shí)刻輸出；為了得到原系統(tǒng)在不同剩余仿真時(shí)間時(shí)脈沖輸入響應(yīng)，原系統(tǒng)需要按照tgo循環(huán)，進(jìn)行多次仿真，而伴隨仿真，因?yàn)槌踔岛洼斎脒x擇給定，只需要一次連續(xù)仿真就可以得到所有tgo時(shí)刻的輸出。

原系統(tǒng)為MISO時(shí)，結(jié)論類(lèi)似，還可以得到各輸入導(dǎo)致終端輸出的誤差預(yù)算。假設(shè)原系統(tǒng)x(0)=0，可以通過(guò)增加階躍輸入或脈沖輸入轉(zhuǎn)化得到[6]。

1.3 協(xié)方差分析的伴隨

(18)

式中：X(t)=E(x(t)xT(t))；Y(t)=E(y(t)·yT(t))。利用該方程進(jìn)行仿真分析通常稱(chēng)為協(xié)方差分析。協(xié)方差矩陣X(t)、Y(t)對(duì)角線(xiàn)元素表示各分量的方差(均方)，非對(duì)角線(xiàn)元素表示各分量之間的協(xié)方差。式(18)是一個(gè)矩陣微分方程，如果將X(t)作為狀態(tài)，U(t)作為輸入，Y(t)作為輸出，仍是線(xiàn)性時(shí)變系統(tǒng)，這可從如下等價(jià)系統(tǒng)得到：

(19)

式中：vec(·)表示矩陣?yán)边\(yùn)算，即將矩陣按照列的順序，一列接一列地組成一個(gè)長(zhǎng)向量；?指克羅內(nèi)克積。顯然，系統(tǒng)式(19)就是常見(jiàn)的狀態(tài)空間描述形式。線(xiàn)性矩陣微分方程式(18)的解為

X(t)=Φ(t,0)X(0)ΦT(t,0)+

(20)

類(lèi)似式(2)與式(10)的關(guān)系，利用式(19)，可以直接得到系統(tǒng)式(18)的伴隨系統(tǒng)的一個(gè)矩陣微分方程實(shí)現(xiàn)為

(21)

系統(tǒng)式(21)微分方程的解為

Z(t)=Ψ(t,0)Z(0)ΨT(t,0)+

(22)

容易證明，系統(tǒng)式(18)與式(21)之間存在類(lèi)似于式(12)的如下關(guān)系：

tr(ZT(0)X(tf))-tr(ZT(tf)X(0))=

tr(UT(t)W(tf-t)))dt

(23)

式中：tr (·)表示方陣的跡。

系統(tǒng)式(21)可以稱(chēng)為協(xié)方差分析式(18)的伴隨，利用式(23)能夠得到類(lèi)似于確定性線(xiàn)性系統(tǒng)的性質(zhì)1、性質(zhì)2的伴隨解釋?zhuān)徊贿^(guò)輸入輸出換成功率譜密度和協(xié)方差矩陣，不但可以分析均方響應(yīng)(對(duì)角線(xiàn)元素)，還可以分析各輸入分量之間相關(guān)性(非對(duì)角線(xiàn)元素)對(duì)輸出的影響。作為例子，考慮MISO線(xiàn)性系統(tǒng)式(2)的協(xié)方差傳播微分方程式(18)，給定初值X(0)，輸入U(xiǎn)(t)；選取伴隨系統(tǒng)式(21)的輸入V(t)=δ(t)，初值Z(0)=0，或者V(t)=0，Z(0)=CT(tf)C(tf)，代入式(23)可得

(24)

同時(shí)考慮相應(yīng)的確定性伴隨系統(tǒng)式(10)，即如果選取伴隨系統(tǒng)式(21)中Z(0)=0，V(t)=δ(t)，同樣選取伴隨系統(tǒng)式(10)中z(0)=0，ν(t)=δ(t)；如果選取伴隨系統(tǒng)式(21)中Z(0)=CT(tf)C(tf)，V(t)=0，同樣選取伴隨系統(tǒng)式(10)中z(0)=CT(tf)，ν(t)=0。那么，結(jié)合式(11)、式(22)，同時(shí)利用跡的循環(huán)性質(zhì)，式(24)可寫(xiě)為

Y(tf)=zT(tf)X(0)z(tf)+

(25)

如果初值x(0)各分量之間不相關(guān)，即X(0)的非對(duì)角線(xiàn)元素為0；U(t)=I，即u(t)為單位功率譜密度白噪聲過(guò)程，則有

(26)

至此，利用伴隨系統(tǒng)的數(shù)學(xué)定義式，對(duì)確定性線(xiàn)性系統(tǒng)和隨機(jī)線(xiàn)性系統(tǒng)的伴隨仿真給出了統(tǒng)一的解釋?zhuān)粏螐姆抡娼Y(jié)果來(lái)考慮，這里給出的解釋本質(zhì)上不需要借助于狀態(tài)轉(zhuǎn)移矩陣或脈沖響應(yīng)矩陣，形式上也更簡(jiǎn)單。

2 滑翔動(dòng)力學(xué)模型

考慮地球模型為靜止均質(zhì)圓球，高超聲速飛行器縱平面內(nèi)再入滑翔的受力分析如圖1所示。圖中：Oexiyi為地心坐標(biāo)系；Obxgyg為地理坐標(biāo)系；Obxb為彈體坐標(biāo)系x軸；h為高度；γ為彈道傾角；s為射程；α為攻角；v為速度；L為升力；D為阻力；g為重力加速度；m為質(zhì)量；Re為地球半徑。

圖1 變量符號(hào)及受力分析示意圖Fig.1 Schematic of variable symbol and force analysis

滑翔動(dòng)力學(xué)微分方程可表示為

(27)

式中：

(28)

其中：ρ為大氣密度，通常采用指數(shù)大氣模型；ρ0=1.225 kg/m3；β=1.389×10-4m-1；g0為海平面重力加速度，取g0=9.806 65 m/s2；S為氣動(dòng)參考面積；CL為升力系數(shù)；CD為阻力系數(shù)。本文在仿真計(jì)算時(shí)采用通用航空飛行器CAV-H[22]的相關(guān)模型數(shù)據(jù)，該飛行器質(zhì)量m=907 kg，氣動(dòng)參考面積S=0.483 9 m2，升力系數(shù)和阻力系數(shù)關(guān)于攻角的擬合公式[3]為

CL=0.046 75α+0.105 68

CD=0.000 508α2+0.004 228α+0.016 1

式中：攻角單位為(°)。

3 平穩(wěn)滑翔彈道

針對(duì)常值攻角控制輸入，分析平穩(wěn)滑翔彈道的性質(zhì)。按照定義[3]，給定滑翔初始速度v0，以初始彈道傾角γ0和高度h0為優(yōu)化變量，以式(29)為性能指標(biāo)，同時(shí)滿(mǎn)足微分方程約束式(27)，優(yōu)化得到的彈道稱(chēng)為平穩(wěn)滑翔彈道。

(29)

考慮時(shí)間常數(shù)tc對(duì)優(yōu)化結(jié)果的影響，表1列出了選取不同時(shí)間常數(shù)tc優(yōu)化得到的初始彈道傾角和初始高度。本算例優(yōu)化條件是：初始速度v0=7 000 m/s，常值攻角輸入α=15°，采用序列二次規(guī)劃(SQP)方法進(jìn)行優(yōu)化。從表1可以看出，不同的時(shí)間常數(shù)tc得到的優(yōu)化結(jié)果之間相差非常小，從后面分析還將會(huì)看到，這些微小差異對(duì)整個(gè)平穩(wěn)滑翔彈道的影響也是可以忽略的。基于此，后文仿真計(jì)算中，將時(shí)間常數(shù)固定為tc=1 000 s來(lái)獲得平穩(wěn)滑翔彈道。

根據(jù)定義，在每一個(gè)初始速度都可以通過(guò)優(yōu)化得到對(duì)應(yīng)的平穩(wěn)滑翔彈道。圖2給出了不同初始速度下優(yōu)化得到的平穩(wěn)滑翔彈道對(duì)比曲線(xiàn)。可以明顯看出，平穩(wěn)滑翔彈道各狀態(tài)量變化平穩(wěn)；從圖2(e)、(f)還可以得到平穩(wěn)滑翔彈道另一重要特征，即不同初始速度下優(yōu)化得到的平穩(wěn)滑翔彈道高度-速度曲線(xiàn)、彈道傾角-速度曲線(xiàn)幾乎重合，或者說(shuō)初始速度小的平穩(wěn)滑翔彈道是初始速度大的一部分，以較小的初始速度優(yōu)化得到平穩(wěn)滑翔彈道初值是在初始速度大的彈道上，這表明平穩(wěn)滑翔彈道上任一點(diǎn)狀態(tài)(h、γ、v)作為初值都滿(mǎn)足平穩(wěn)滑翔彈道的定義，或者說(shuō)平穩(wěn)滑翔彈道上的任一段仍然是平穩(wěn)滑翔彈道，這可以稱(chēng)為平穩(wěn)滑翔彈道定義的一致性。

表1 不同時(shí)間常數(shù)tc的優(yōu)化結(jié)果Table 1 Optimization results with different time constants tc

圖2 不同初始速度下的平穩(wěn)滑翔彈道對(duì)比Fig.2 Comparison of steady glide trajectories with different initial velocities

4 擾動(dòng)運(yùn)動(dòng)及線(xiàn)性化

從第3節(jié)分析可以得出，如果飛行器再入滑翔的起始高度、彈道傾角和速度滿(mǎn)足一定關(guān)系(在某條平穩(wěn)滑翔彈道上)，那么飛行器將會(huì)一直沿著平穩(wěn)滑翔彈道飛行。然而實(shí)際再入時(shí)滑翔初始高度、彈道傾角和速度很可能與期望的平穩(wěn)滑翔初始狀態(tài)存在偏差，并且飛行過(guò)程中也會(huì)存在大氣密度、風(fēng)等多種干擾因素，為了研究這些擾動(dòng)因素對(duì)平穩(wěn)滑翔彈道的影響，建立如下帶有誤差干擾模型的動(dòng)力學(xué)方程：

(30)

式中：εL、εD分別為比例形式的升力、阻力干擾輸入。同時(shí)考慮相對(duì)于平穩(wěn)滑翔彈道存在初始高度擾動(dòng)δh0、初始彈道傾角擾動(dòng)δγ0和初始速度擾動(dòng)δv0的情況。

分析干擾對(duì)標(biāo)準(zhǔn)平穩(wěn)滑翔彈道的影響，最直觀的方法就是在實(shí)際干擾條件下對(duì)方程式(30)求解，再與標(biāo)準(zhǔn)平穩(wěn)滑翔彈道作差，但是該方法的缺點(diǎn)是計(jì)算量大，需要數(shù)值求解非線(xiàn)性微分方程，并不是直接分析干擾與彈道偏差的關(guān)系，不方便應(yīng)用于制導(dǎo)[20]。另一種方法是基于小擾動(dòng)假設(shè)的攝動(dòng)法，將帶有誤差干擾模型的動(dòng)力學(xué)方程式(30)在標(biāo)準(zhǔn)平穩(wěn)滑翔彈道附近線(xiàn)性化，得到如下?tīng)顟B(tài)空間線(xiàn)性系統(tǒng)：

(31)

式中：

其中：

a12=vscosγs

a13=sinγs

a32=-gcosγs

式中：下標(biāo)為s的量(包括L、D、g和所有偏導(dǎo)數(shù))指標(biāo)準(zhǔn)平穩(wěn)滑翔彈道上的數(shù)據(jù)。所以，系數(shù)矩陣中所有分量元素都是時(shí)間t的已知函數(shù)。線(xiàn)性系統(tǒng)式(31)近似描述了真實(shí)彈道偏差的動(dòng)態(tài)特性，即

δh≈h-hs

δγ≈γ-γs

δv≈v-vs

δs≈s-ss

利用第3節(jié)獲得的初始速度v0=7 000 m/s的平穩(wěn)滑翔彈道，飛行時(shí)間為tf=3 044 s(此時(shí)滑翔高度大約為30 km)。以此為標(biāo)準(zhǔn)平穩(wěn)滑翔彈道，加入常值擾動(dòng)，然后利用式(30)進(jìn)行非線(xiàn)性仿真，利用式(31)進(jìn)行線(xiàn)性系統(tǒng)仿真，得到各狀態(tài)量的偏差結(jié)果如圖3所示。本算例中小擾動(dòng)分別為

(32)

由圖3中可直觀看出，加入擾動(dòng)后，滑翔彈道呈現(xiàn)出明顯的圍繞平穩(wěn)滑翔彈道振蕩的特性；而且線(xiàn)性化系統(tǒng)仿真得到的彈道偏差結(jié)果與真實(shí)彈道偏差很接近，這說(shuō)明得到的線(xiàn)性化系統(tǒng)可以很好地近似非線(xiàn)性系統(tǒng)在平穩(wěn)滑翔彈道附近的動(dòng)態(tài)特性。這為利用線(xiàn)性系統(tǒng)伴隨仿真方法分析平穩(wěn)滑翔彈道奠定了基礎(chǔ)。

圖3 平穩(wěn)滑翔彈道小擾動(dòng)偏差Fig.3 Trajectory deviation from steady glide in case of small disturbance

5 伴隨仿真算例分析

利用線(xiàn)性化方程式(31)可以對(duì)平穩(wěn)滑翔彈道進(jìn)行伴隨仿真分析，下面分別對(duì)確定性常值小擾動(dòng)和隨機(jī)擾動(dòng)輸入進(jìn)行仿真算例研究。

5.1 確定性常值小擾動(dòng)

仍然分析第3節(jié)中得到的v0=7 000 m/s的平穩(wěn)滑翔彈道，研究確定性常值小擾動(dòng)式(32)對(duì)終端高度hf=h(tf)的影響，tf=3 044 s，此時(shí)滑翔高度大約為30 km。根據(jù)式(10)構(gòu)造系統(tǒng)式(31)的伴隨系統(tǒng)，由于考慮高度偏差，取式(31)的輸出為y(t)=δh，相應(yīng)的輸出矩陣為C(t)=[1 0 00]；選取伴隨系統(tǒng)式(10)的控制v(t)=0，狀態(tài)初值z(mì)(0)=CT(tf)，仿真結(jié)果如圖4所示。

圖4給出了各干擾因素分別導(dǎo)致的終端高度偏差，即(δhf|δh0、(δhf|δγ0、(δhf|δv0、(δhf|εL、(δhf|εD，這得益于伴隨仿真的誤差預(yù)算性質(zhì)(式(14)～式(16))，(δhf|T表示合擾動(dòng)結(jié)果；圖中橫軸tgo表示剩余飛行時(shí)間，曲線(xiàn)上對(duì)應(yīng)點(diǎn)可以解釋為在剩余飛行時(shí)間是tgo時(shí)，加入相應(yīng)擾動(dòng)得到終端tf時(shí)刻的高度偏差。為了與實(shí)際終端高度偏差對(duì)比，利用式(30)對(duì)各擾動(dòng)因素單獨(dú)(其他擾動(dòng)設(shè)置為0)進(jìn)行非線(xiàn)性仿真，得到實(shí)際終端高度hf與標(biāo)準(zhǔn)平穩(wěn)滑翔彈道終端高度hsf之差，結(jié)果如圖5所示?？梢钥闯觯蔷€(xiàn)性仿真與伴隨仿真結(jié)果很接近，這也再次表明線(xiàn)性化系統(tǒng)式(31)能夠很好地反映平穩(wěn)滑翔彈道在小擾動(dòng)輸入下的動(dòng)態(tài)特性。圖5中伴隨仿真所有數(shù)據(jù)是通過(guò)一次連續(xù)仿真得到，不同擾動(dòng)因素的影響可以由性質(zhì)1解釋?zhuān)煌Ｓ囡w行時(shí)間tgo的結(jié)果可以由性質(zhì)2得到；而非線(xiàn)性仿真需要針對(duì)各不同干擾輸入分別進(jìn)行，在每一個(gè)干擾輸入下，又需要對(duì)tgo循環(huán)多次(本算例取循環(huán)間隔為10s，大約需要300次)仿真得到，即開(kāi)始時(shí)按照無(wú)擾動(dòng)平穩(wěn)滑翔彈道仿真，當(dāng)飛行時(shí)間為tf-tgo，或者說(shuō)剩余飛行時(shí)間為tgo時(shí)加入相應(yīng)的小擾動(dòng)，仿真得到tf時(shí)刻的結(jié)果。顯然，伴隨仿真運(yùn)算量要遠(yuǎn)小于非線(xiàn)性仿真?？紤]到平穩(wěn)滑翔彈道定義的一致性，tgo可以對(duì)應(yīng)到不同初始速度v0優(yōu)化得到的平穩(wěn)滑翔彈道，這樣伴隨仿真的結(jié)果可以解釋為不同初始速度的平穩(wěn)滑翔彈道在加入同樣的小擾動(dòng)式(32)時(shí)，導(dǎo)致終端(標(biāo)準(zhǔn)平穩(wěn)滑翔彈道高度大約為30 km時(shí))高度的偏差；相應(yīng)的數(shù)據(jù)可以用于平穩(wěn)滑翔彈道制導(dǎo)方法的設(shè)計(jì)與分析。

圖4 伴隨一次仿真輸出各擾動(dòng)帶來(lái)的終端高度偏差Fig.4 One adjoint simulation yields terminal height deviation for each disturbance

圖5 各擾動(dòng)因素對(duì)終端高度偏差的影響(非線(xiàn)性仿真和伴隨仿真)Fig.5 Influence of each disturbance on terminal height deviation(nonlinear and adjoint simulation)

從圖4和圖5還可以得到，在加入常值小擾動(dòng)式(32)的條件下，初始速度、升力、阻力等小擾動(dòng)(δv0、εL、εD)對(duì)終端高度偏差的影響隨tgo增加而增大，而初始高度和初始彈道傾角等擾動(dòng)(δh0、δγ0)對(duì)終端高度偏差的影響隨tgo增加振蕩衰減；合擾動(dòng)結(jié)果在tgo較小時(shí)主要由δh0決定，在tgo較大時(shí)，主要取決于δv0、εL、εD等，δγ0影響很小。當(dāng)需要另計(jì)算一組新的小擾動(dòng)：

導(dǎo)致的終端高度偏差時(shí)，不需要重新仿真，利用式(14)～式(16)或線(xiàn)性疊加原理，對(duì)已經(jīng)得到的仿真結(jié)果進(jìn)行簡(jiǎn)單的線(xiàn)性組合即可。

(kv·δhf|δv0+(kL·δhf|εL+kD·(δhf|εD

類(lèi)似地，還可以分析確定性常值小擾動(dòng)式(32)對(duì)終端速度、終端射程等影響，只是系統(tǒng)式(31)相應(yīng)的輸出矩陣分別取C(t)=[0 0 10]和C(t)=[0 0 0 1]。仿真結(jié)果如圖6所示?？梢钥闯?，在小擾動(dòng)式(32)的條件下，影響平穩(wěn)滑翔終端速度、射程的主要因素為δv0、εL、εD，而且小擾動(dòng)會(huì)帶來(lái)較大的射程偏差(百公里級(jí))。

為了更清晰地表達(dá)各擾動(dòng)對(duì)終端狀態(tài)的影響，表2～表4分別列出了在不同剩余飛行時(shí)間加入單位各擾動(dòng)導(dǎo)致的終端高度偏差、終端速度偏差和終端射程偏差。其中單位初始高度擾動(dòng)設(shè)為1 km,單位初始彈道傾角擾動(dòng)設(shè)為1°,單位初始速度擾動(dòng)設(shè)為1 m/s,單位升力或阻力擾動(dòng)設(shè)為1%。對(duì)比表中數(shù)據(jù)可以得到，終端狀態(tài)的偏差關(guān)

圖6 伴隨一次仿真輸出各擾動(dòng)帶來(lái)的終端速度偏差和終端射程偏差Fig.6 One adjoint simulation yields terminal range deviation and terminal velocity deviation for each disturbance

于初始彈道傾角擾動(dòng)最敏感，即相比于其他擾動(dòng)，較小初始彈道傾角擾動(dòng)會(huì)帶來(lái)較大的終端狀態(tài)偏差。

5.2 隨機(jī)擾動(dòng)

(33)

表2 單位各擾動(dòng)導(dǎo)致的終端高度偏差Table 2 Terminal height deviation with respect to each disturbance

為便于對(duì)比，給出了在合隨機(jī)擾動(dòng)條件下終端高度和射程的非線(xiàn)性蒙特卡羅仿真結(jié)果，如圖8所示。特別地，以終端速度為例，針對(duì)每一個(gè)單獨(dú)的隨機(jī)擾動(dòng)因素，按照tgo間隔50 s循環(huán)(大約60次)，每一個(gè)tgo下進(jìn)行1 000次非線(xiàn)性蒙特卡羅仿真，得到終端速度的散布情況(標(biāo)準(zhǔn)差)，結(jié)果如圖9所示。

從圖8和圖9可以看出，伴隨仿真結(jié)果與多次非線(xiàn)性蒙特卡羅仿真十分吻合，但是蒙特卡羅仿真計(jì)算量(非線(xiàn)性仿真次數(shù)大約是6×60×1 000=360 000)要遠(yuǎn)大于伴隨一次仿真。

表3 單位各擾動(dòng)導(dǎo)致的終端速度偏差

表4 單位各擾動(dòng)導(dǎo)致的終端射程偏差

圖7 各隨機(jī)擾動(dòng)輸入下終端高度標(biāo)準(zhǔn)差、終端速度標(biāo)準(zhǔn)差、終端射程標(biāo)準(zhǔn)差伴隨仿真結(jié)果Fig.7 Adjoint simulation results of standard deviation of terminal height, velocity and range for each random disturbance

圖8 合隨機(jī)擾動(dòng)下終端高度標(biāo)準(zhǔn)差和終端射程標(biāo)準(zhǔn)差Fig.8 Total standard deviation of terminal height and range with all random disturbances

圖9 各隨機(jī)擾動(dòng)單獨(dú)作用時(shí)終端速度標(biāo)準(zhǔn)差(非線(xiàn)性蒙特卡羅仿真和伴隨仿真對(duì)比)Fig.9 Standard deviation of terminal velocity for each random disturbance (comparison between nonlinear Monte Carlo simulations and adjoint simulation)

6 結(jié) 論

本文探討了伴隨仿真方法，并對(duì)高超聲速飛行器平穩(wěn)滑翔彈道進(jìn)行了伴隨仿真分析，主要結(jié)論有：

1) 利用伴隨系統(tǒng)的數(shù)學(xué)定義式，從新的角度給出了伴隨仿真方法的統(tǒng)一解釋?zhuān)ㄕ`差預(yù)算性質(zhì)和伴隨一次仿真結(jié)果一般意義。

2) 隨機(jī)線(xiàn)性系統(tǒng)伴隨仿真方法，本質(zhì)上等同于線(xiàn)性系統(tǒng)協(xié)方差分析的伴隨；利用矩陣向量化運(yùn)算，協(xié)方差分析的伴隨與確定性線(xiàn)性系統(tǒng)的伴隨數(shù)學(xué)形式相同。

3) 對(duì)于給定常值攻角，平穩(wěn)滑翔彈道定義具有一致性，即平穩(wěn)滑翔彈道上的任一點(diǎn)作為初值仍然滿(mǎn)足平穩(wěn)滑翔彈道的定義。

4) 在小擾動(dòng)假設(shè)下，伴隨仿真分析了滑翔初始高度、彈道傾角、速度及過(guò)程中的升力、阻力，在有確定性常值小擾動(dòng)和隨機(jī)擾動(dòng)時(shí)，對(duì)平穩(wěn)滑翔彈道的終端高度、速度、射程的影響。結(jié)果表明，終端狀態(tài)的偏差關(guān)于初始彈道傾角擾動(dòng)最敏感；并且對(duì)比非線(xiàn)性仿真和蒙特卡羅仿真，結(jié)果吻合，但是伴隨仿真方法的計(jì)算效率優(yōu)勢(shì)明顯。