高超聲速飛行器平穩滑翔彈道擾動運動伴隨分析

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

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

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

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

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

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

1 伴隨仿真方法

為了闡述伴隨仿真方法的原理，本節從線性時變系統的伴隨系統的定義出發，利用線性系統與其伴隨系統之間的本質關系(定義式)，給出伴隨仿真方法的統一解釋。

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

(1)

式中：〈·,·〉為空間U或Y上的內積。

1.1 伴隨系統狀態空間表示

考慮如下狀態空間描述的線性時變系統：

(2)

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

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

(3)

式中：Φ(t,τ)為狀態轉移矩陣。

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

(4)

式中：⊕為直和。

定義內積為

(5)

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

(6)

存在如下狀態空間實現[21]：

(7)

為了驗證系統式(7)是式(2)的伴隨的狀態空間實現，進行如下推導運算：

(-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)左右兩邊進行簡單移項，并結合內積定義式(5)，可以看出式(8)正是伴隨系統定義中式(1)的狀態空間變量表示，這說明系統式(7)和式(2)互為各自伴隨的狀態空間實現。式(8)通常稱為Bliss公式[19]或伴隨定理[17]。

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

(9)

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

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

(10)

式(10)的解為

ν(τ)dτ

(11)

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

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

(12)

1.2 伴隨仿真統一解釋

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

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

(13)

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

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

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

(14)

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

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

(15)

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

(16)

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

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

y(tf)=w(tgo)

(17)

式(17)對所有tgo∈[0,tf]成立，這說明原系統在仿真時刻為tf-tgo，或者說在剩余仿真時間為tgo時，加入脈沖輸入，得到終端輸出y(tf)，等于伴隨仿真在tgo時刻輸出；為了得到原系統在不同剩余仿真時間時脈沖輸入響應，原系統需要按照tgo循環，進行多次仿真，而伴隨仿真，因為初值和輸入選擇給定，只需要一次連續仿真就可以得到所有tgo時刻的輸出。

原系統為MISO時，結論類似，還可以得到各輸入導致終端輸出的誤差預算。假設原系統x(0)=0，可以通過增加階躍輸入或脈沖輸入轉化得到[6]。

1.3 協方差分析的伴隨

(18)

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

(19)

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

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

(20)

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

(21)

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

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

(22)

容易證明，系統式(18)與式(21)之間存在類似于式(12)的如下關系：

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

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

(23)

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

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

(24)

同時考慮相應的確定性伴隨系統式(10)，即如果選取伴隨系統式(21)中Z(0)=0，V(t)=δ(t)，同樣選取伴隨系統式(10)中z(0)=0，ν(t)=δ(t)；如果選取伴隨系統式(21)中Z(0)=CT(tf)C(tf)，V(t)=0，同樣選取伴隨系統式(10)中z(0)=CT(tf)，ν(t)=0。那么，結合式(11)、式(22)，同時利用跡的循環性質，式(24)可寫為

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

(25)

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

(26)

至此，利用伴隨系統的數學定義式，對確定性線性系統和隨機線性系統的伴隨仿真給出了統一的解釋；單從仿真結果來考慮，這里給出的解釋本質上不需要借助于狀態轉移矩陣或脈沖響應矩陣，形式上也更簡單。

2 滑翔動力學模型

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

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

滑翔動力學微分方程可表示為

(27)

式中：

(28)

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

CL=0.046 75α+0.105 68

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

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

3 平穩滑翔彈道

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

(29)

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

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

表1 不同時間常數tc的優化結果Table 1 Optimization results with different time constants tc

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

4 擾動運動及線性化

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

(30)

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

分析干擾對標準平穩滑翔彈道的影響，最直觀的方法就是在實際干擾條件下對方程式(30)求解，再與標準平穩滑翔彈道作差，但是該方法的缺點是計算量大，需要數值求解非線性微分方程，并不是直接分析干擾與彈道偏差的關系，不方便應用于制導[20]。另一種方法是基于小擾動假設的攝動法，將帶有誤差干擾模型的動力學方程式(30)在標準平穩滑翔彈道附近線性化，得到如下狀態空間線性系統：

(31)

式中：

其中：

a12=vscosγs

a13=sinγs

a32=-gcosγs

式中：下標為s的量(包括L、D、g和所有偏導數)指標準平穩滑翔彈道上的數據。所以，系數矩陣中所有分量元素都是時間t的已知函數。線性系統式(31)近似描述了真實彈道偏差的動態特性，即

δh≈h-hs

δγ≈γ-γs

δv≈v-vs

δs≈s-ss

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

(32)

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

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

5 伴隨仿真算例分析

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

5.1 確定性常值小擾動

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

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

圖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增加振蕩衰減；合擾動結果在tgo較小時主要由δh0決定，在tgo較大時，主要取決于δv0、εL、εD等，δγ0影響很小。當需要另計算一組新的小擾動：

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

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

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

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

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

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

5.2 隨機擾動

(33)

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

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

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

表3 單位各擾動導致的終端速度偏差

表4 單位各擾動導致的終端射程偏差

圖7 各隨機擾動輸入下終端高度標準差、終端速度標準差、終端射程標準差伴隨仿真結果Fig.7 Adjoint simulation results of standard deviation of terminal height, velocity and range for each random disturbance

圖8 合隨機擾動下終端高度標準差和終端射程標準差Fig.8 Total standard deviation of terminal height and range with all random disturbances

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

6 結 論

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

1) 利用伴隨系統的數學定義式，從新的角度給出了伴隨仿真方法的統一解釋，包括誤差預算性質和伴隨一次仿真結果一般意義。

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

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

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