赫泰龍， 陳萬春,*， 劉芳

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

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

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

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

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

1 伴隨仿真方法

為了闡述伴隨仿真方法的原理，本節(jié)從線性時變系統(tǒng)的伴隨系統(tǒng)的定義出發(fā)，利用線性系統(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ài)空間描述的線性時變系統(tǒng)：

(2)

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

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

(3)

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

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

(4)

式中：⊕為直和。

定義內(nèi)積為

(5)

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

(6)

存在如下狀態(tài)空間實現(xiàn)[21]：

(7)

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

(-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)

對上式從t0到tf積分可得到

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

(8)

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

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

(9)

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

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

(10)

式(10)的解為

ν(τ)dτ

(11)

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

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

(12)

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

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

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

(15)

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

(16)

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

性質(zhì)2(伴隨一次仿真的通用解釋) 為了簡化敘述，不失一般性，這里假設(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)對所有tgo∈[0,tf]成立，這說明原系統(tǒng)在仿真時刻為tf-tgo，或者說在剩余仿真時間為tgo時，加入脈沖輸入，得到終端輸出y(tf)，等于伴隨仿真在tgo時刻輸出；為了得到原系統(tǒng)在不同剩余仿真時間時脈沖輸入響應(yīng)，原系統(tǒng)需要按照tgo循環(huán)，進行多次仿真，而伴隨仿真，因為初值和輸入選擇給定，只需要一次連續(xù)仿真就可以得到所有tgo時刻的輸出。

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

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

(18)

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

(19)

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

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

(20)

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

(21)

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

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

(22)

容易證明，系統(tǒng)式(18)與式(21)之間存在類似于式(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)可以稱為協(xié)方差分析式(18)的伴隨，利用式(23)能夠得到類似于確定性線性系統(tǒng)的性質(zhì)1、性質(zhì)2的伴隨解釋，只不過輸入輸出換成功率譜密度和協(xié)方差矩陣，不但可以分析均方響應(yīng)(對角線元素)，還可以分析各輸入分量之間相關(guān)性(非對角線元素)對輸出的影響。作為例子，考慮MISO線性系統(tǒng)式(2)的協(xié)方差傳播微分方程式(18)，給定初值X(0)，輸入U(t)；選取伴隨系統(tǒng)式(21)的輸入V(t)=δ(t)，初值Z(0)=0，或者V(t)=0，Z(0)=CT(tf)C(tf)，代入式(23)可得

(24)

同時考慮相應(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)，同時利用跡的循環(huán)性質(zhì)，式(24)可寫為

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

(25)

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

(26)

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

2 滑翔動力學(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 變量符號及受力分析示意圖Fig.1 Schematic of variable symbol and force analysis

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

(27)

式中：

(28)

其中：ρ為大氣密度，通常采用指數(shù)大氣模型；ρ0=1.225 kg/m3；β=1.389×10-4m-1；g0為海平面重力加速度，取g0=9.806 65 m/s2；S為氣動參考面積；CL為升力系數(shù)；CD為阻力系數(shù)。本文在仿真計算時采用通用航空飛行器CAV-H[22]的相關(guān)模型數(shù)據(jù)，該飛行器質(zhì)量m=907 kg，氣動參考面積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)滑翔彈道

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

(29)

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

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

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

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

4 擾動運動及線性化

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

(30)

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

分析干擾對標(biāo)準(zhǔn)平穩(wěn)滑翔彈道的影響，最直觀的方法就是在實際干擾條件下對方程式(30)求解，再與標(biāo)準(zhǔn)平穩(wěn)滑翔彈道作差，但是該方法的缺點是計算量大，需要數(shù)值求解非線性微分方程，并不是直接分析干擾與彈道偏差的關(guān)系，不方便應(yīng)用于制導(dǎo)[20]。另一種方法是基于小擾動假設(shè)的攝動法，將帶有誤差干擾模型的動力學(xué)方程式(30)在標(biāo)準(zhǔn)平穩(wěn)滑翔彈道附近線性化，得到如下狀態(tài)空間線性系統(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ù)矩陣中所有分量元素都是時間t的已知函數(shù)。線性系統(tǒng)式(31)近似描述了真實彈道偏差的動態(tài)特性，即

δh≈h-hs

δγ≈γ-γs

δv≈v-vs

δs≈s-ss

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

(32)

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

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

5 伴隨仿真算例分析

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

5.1 確定性常值小擾動

仍然分析第3節(jié)中得到的v0=7 000 m/s的平穩(wěn)滑翔彈道，研究確定性常值小擾動式(32)對終端高度hf=h(tf)的影響，tf=3 044 s，此時滑翔高度大約為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(0)=CT(tf)，仿真結(jié)果如圖4所示。

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

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

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

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

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

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

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

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

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

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

5.2 隨機擾動

(33)

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

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

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

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

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

圖7 各隨機擾動輸入下終端高度標(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 合隨機擾動下終端高度標(biāo)準(zhǔn)差和終端射程標(biāo)準(zhǔn)差Fig.8 Total standard deviation of terminal height and range with all random disturbances

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

6 結(jié) 論

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

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

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

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

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