王培臣，張睿軒，閆循良

(1.西北工業(yè)大學(xué) 航天學(xué)院·西安·710072；2.中國航空工業(yè)空氣動力研究院·沈陽·110034)

0 引 言

高超聲速滑翔飛行器通常采用大升阻比氣動外形，可進(jìn)行長距離機動飛行，具有飛行速度快、機動靈活、突防效率高等特點，因此得到了各軍事強國的廣泛關(guān)注。彈道優(yōu)化設(shè)計作為高超聲速滑翔飛行器研究的難點與熱點問題之一，吸引了大量國內(nèi)外學(xué)者對此開展研究[1-5]。俯沖攻擊段作為高超聲速滑翔飛行器的最后飛行段[6]，對于飛行器能否精確命中目標(biāo)至關(guān)重要，因此也得到了相應(yīng)的關(guān)注。徐明亮等[7]利用偽譜法對臨近空間飛行器鉸鏈力矩最小的俯沖彈道進(jìn)行了優(yōu)化設(shè)計，得到近似垂直的攻擊彈道，但未考慮實際飛行過程中的隨機干擾。喬浩等[8]利用高斯偽譜法對高超聲速滑翔飛行器俯沖段翻身下壓問題進(jìn)行了優(yōu)化，得到了高安全性的快速下壓軌跡，但同樣未考慮不確定因素影響，設(shè)計結(jié)果過于理想。

雖然針對俯沖段軌跡優(yōu)化問題已有大量研究，但是傳統(tǒng)的優(yōu)化設(shè)計方法尚未考慮不確定性因素的影響[7-8]。然而，俯沖段飛行速度大，飛行狀態(tài)及環(huán)境變化劇烈、復(fù)雜，所受氣動力熱環(huán)境惡劣，且高超聲速滑翔飛行器的控制能力有限，因此，不能忽視不確定性因素對彈道的影響和制導(dǎo)控制系統(tǒng)的干擾。為保證俯沖攻擊落點精度，提高彈道抗干擾能力，有必要對俯沖攻擊段軌跡進(jìn)行魯棒優(yōu)化。

目前，分析不確定性因素對軌跡影響的主要方法有：蒙特卡羅打靶分析法[9]、無跡變換法[10]、混沌多項式法[11-12]以及線性協(xié)方差分析法[13-15]等，而針對如何提高彈道抗干擾能力的相關(guān)研究較少。文獻(xiàn)[14]將魯棒性指標(biāo)引入到滑翔炮彈軌跡優(yōu)化中，得到的計算結(jié)果較為理想。因此，本文借鑒其思路，在傳統(tǒng)高超聲速滑翔飛行器彈道優(yōu)化方法的基礎(chǔ)上，引入魯棒性設(shè)計指標(biāo),以降低不確定性因素干擾的影響，進(jìn)而提升彈道抗干擾能力，同時降低制導(dǎo)控制系統(tǒng)負(fù)擔(dān)。

1 俯沖運動數(shù)學(xué)模型

1.1 滑翔飛行器俯沖段動力學(xué)建模

由于目前高超聲速滑翔飛行器普遍采用腹部防熱設(shè)計，考慮到飛行器在俯沖段的熱載荷較大，可在此過程中采用翻身下壓策略[8]。即飛行器采用正攻角翻身飛行，利用有防熱設(shè)計的腹部迎接來流，同時也可滿足彈道快速下壓需求。此外，考慮到俯沖段飛行時間較短，因此可忽略地球自轉(zhuǎn)及扁率影響。因此，俯沖段三自由度質(zhì)心運動模型可描述為

(1)

式中，h、λ、φ為飛行器位置參數(shù)，分別為飛行高度、地心經(jīng)度、地心緯度；V、θ、σ為飛行器速度參數(shù)，分別為飛行速度、速度傾角、速度偏角；m為飛行器質(zhì)量；α、υ為飛行器控制變量，分別為攻角和傾側(cè)角；L、D為飛行器的升、阻力，計算公式如下

(2)

式中，ρ為大氣密度；S為飛行器參考面積；CL、CD為飛行器氣動升、阻力系數(shù)。

1.2 不確定性因素建模

在實際飛行過程中，高超聲速滑翔飛行器會受到多種不確定性因素干擾，從而使得實際飛行彈道偏離理論軌跡。結(jié)合滑翔俯沖段飛行環(huán)境特點和飛行特性，本文考慮的不確定性因素主要有氣動力模型偏差、大氣模型偏差、陣風(fēng)干擾偏差以及初始點狀態(tài)偏差。

1.2.1 氣動力模型偏差

由于地面實驗無法準(zhǔn)確模擬實際飛行環(huán)境，以及制造加工工藝和飛行加熱燒蝕帶來的誤差，使得彈體氣動力系數(shù)存在較大的偏差和不確定性，同時大氣參數(shù)偏差也會導(dǎo)致飛行器實際氣動力與理論設(shè)計值存在偏差，但氣動力系數(shù)隨馬赫數(shù)與攻角的變化規(guī)律不變，所以實際氣動力系數(shù)可表示為

(3)

式中，基準(zhǔn)值CD0、CL0為計算或?qū)嶒炈蔑w行器氣動力系數(shù)；Ncd、Ncl為相互獨立，且服從標(biāo)準(zhǔn)正態(tài)分布的隨機變量；σcd、σcl代表實際氣動力系數(shù)偏離基準(zhǔn)值的程度，假設(shè)偏離值置信度(3σ)為基準(zhǔn)值的15%，則有

σcd=σcl=15%/3

(4)

1.2.2 大氣模型偏差

通常彈道仿真計算中應(yīng)用標(biāo)準(zhǔn)大氣模型，但實際大氣參數(shù)與標(biāo)準(zhǔn)參數(shù)會存在一定偏差，主要為溫度偏差和大氣密度偏差。由于當(dāng)?shù)芈曀僦饕蓽囟葲Q定，因此溫度偏差會對聲速的計算產(chǎn)生影響，進(jìn)而影響馬赫數(shù)的計算。而氣動力系數(shù)通過馬赫數(shù)與攻角插值得到，所以溫度偏差帶來的干擾可在氣動力偏差中加以考慮計算，故此處只考慮大氣密度偏差。參考文獻(xiàn)[14]給出了不同海拔下實際大氣密度的標(biāo)準(zhǔn)差與標(biāo)準(zhǔn)大氣密度的比值關(guān)系，基于此關(guān)系，實際大氣密度可計算如下

(5)

式中，ρ0為標(biāo)準(zhǔn)大氣密度；Nρ為服從標(biāo)準(zhǔn)正態(tài)分布的隨機變量。

1.2.3 陣風(fēng)干擾偏差

在飛行器翻身下壓飛行過程中，陣風(fēng)干擾會使得質(zhì)心運動產(chǎn)生隨機擾動，主要體現(xiàn)在對速度大小和方向的影響，因此要考慮陣風(fēng)對速度參數(shù)產(chǎn)生的影響，即有以下關(guān)系

(6)

1.2.4 初始點狀態(tài)偏差

(7)

2 基于線性協(xié)方差分析的誤差傳播建模

(8)

E(w(t)w(t)T)=R(w)

(9)

式中，R(w)為譜密度矩陣，為

(10)

故根據(jù)不確定環(huán)境下的系統(tǒng)模型，可得狀態(tài)量的協(xié)方差矩陣為

(11)

(12)

式中，對角線各項為實際狀態(tài)偏離基準(zhǔn)值的方差。當(dāng)參考軌跡給定后，協(xié)方差矩陣P可通過如下李雅普諾夫方程求解

(13)

式中，A為雅可比矩陣，由于數(shù)值求解雅可比矩陣計算量大，為提高計算速度，推導(dǎo)如下解析表達(dá)式進(jìn)行計算

(14)

(15)

3 彈道優(yōu)化模型

落點散布大小是衡量彈道抗干擾能力的重要指標(biāo)之一。為增強俯沖彈道對不確定性因素的抗干擾能力，可以將落點散布大小作為代價函數(shù)。在滿足過程約束的前提下，當(dāng)優(yōu)化得到代價函數(shù)的值最小時，即可認(rèn)為得到的優(yōu)化彈道抗干擾能力最強，因此構(gòu)建如下代價函數(shù)

(16)

式中，右端第一項代表落點散布，ω為權(quán)重系數(shù)，其取值可表征軌跡抗干擾能力的大??；第二項保證飛行器實際工作中，控制量相對光滑且不會發(fā)生突變。

分析代價函數(shù)式可知，ω=0時，即代價函數(shù)中不考慮軌跡抗干擾能力，問題將退化為傳統(tǒng)的俯沖攻擊段最優(yōu)彈道求解問題；ω>0時，則該問題為考慮軌跡抗干擾能力的魯棒軌跡優(yōu)化問題。

同時，優(yōu)化過程中應(yīng)滿足下列約束條件：

(17)

式中，K為常值，其取值與飛行器結(jié)構(gòu)相關(guān)。

2)控制約束

(18)

(19)

3)終端約束

(20)

通過以上建模過程，即可將提高彈道抗干擾能力問題轉(zhuǎn)化為多約束條件下的最優(yōu)控制問題。目前求解最優(yōu)控制問題主要有間接法與直接法兩種。前者很難處理復(fù)雜約束條件下的最優(yōu)控制問題，相反地，直接法更適合求解這類問題。在眾多直接法中，高斯偽譜法以較少的參數(shù)和較高的精度優(yōu)勢，已廣泛應(yīng)用于復(fù)雜約束條件下的軌跡優(yōu)化問題[5,16]。因此，本文選用高斯偽譜法作為求解最優(yōu)控制問題的方法。

4 算例仿真與分析

表1 初始狀態(tài)

表2 陣風(fēng)干擾高斯白噪聲功率譜密度

為驗證LinCov方法的可行性，取ω=0進(jìn)行彈道優(yōu)化仿真，并基于優(yōu)化結(jié)果利用LinCov方法和蒙特卡羅打靶法(Monte Carlo，MC)計算得到終端誤差，如表3所示?？梢园l(fā)現(xiàn)，由LinCov方法得到的落點經(jīng)度方差與蒙特卡羅打靶法得到的結(jié)果相對誤差為1.77%，得到的緯度方差相對誤差為6.99%。圖1分別給出了兩種方法得到的落點分布3σ誤差橢圓，也可以看出兩種方法的落點分布范圍相差不大，驗證了線性協(xié)方差分析法估計終端經(jīng)緯度散布的可行性。

表3 線性協(xié)方差分析法與蒙特卡羅打靶法結(jié)果對比

圖1 線性協(xié)方差方法和蒙特卡羅打靶法落點3σ誤差圓對比Fig.1 Comparison of 3σ position error ellipses between linear covariance analysis and Monte Carlo shooting methods at target without considering the influence of uncertain factor

為驗證魯棒軌跡優(yōu)化算法的可行性，取ω=0.1和1，結(jié)合線性協(xié)方差分析法和高斯偽譜法進(jìn)行彈道優(yōu)化，所得最優(yōu)彈道對應(yīng)的攻角、傾側(cè)角控制量曲線分別如圖2、圖3所示。對比控制量曲線可以發(fā)現(xiàn)，隨著權(quán)重系數(shù)的增大，控制量的變化更劇烈，說明為了提高軌跡對不確定性因素的抗干擾能力，需要消耗一定的控制量裕度，且隨著權(quán)重系數(shù)ω的增大，控制裕度消耗的程度增加。

圖2 攻角隨時間變化曲線Fig.2 Time histories of the attack angle variables

圖3 傾側(cè)角隨時間變化曲線Fig.3 Time histories of the tilt angle variables

下面對不同權(quán)重下的彈道優(yōu)化結(jié)果進(jìn)行蒙特卡羅打靶仿真，并統(tǒng)計落點分布情況。表4、表5分別給出了落點分布的均值和方差。當(dāng)ω=0.1時，落點分布均值與目標(biāo)值在經(jīng)度方向偏差為0.009°，緯度方向偏差為0.008°；其落點分布方差與ω=0的結(jié)果相比，經(jīng)度分布方差減小55.12%，緯度分布方差減小52.57%，落點分布范圍明顯減小，表明本文方法可有效地提高彈道抗干擾能力，降低彈道對不確定性因素的敏感度。而對比ω=0.1和ω=1時的落點分布均值和方差可知，后者的數(shù)據(jù)較前者都有所提升，但差異并不大，說明隨著權(quán)重系數(shù)ω的進(jìn)一步增加，落點散布減小但效果有限。因此，在實際應(yīng)用時，需要設(shè)置權(quán)重系數(shù)ω以權(quán)衡控制量裕度與落點密集度之間的關(guān)系，在控制允許范圍內(nèi)可適當(dāng)增大權(quán)重系數(shù)ω以盡可能降低落點散布。

表4 不同權(quán)重情況下優(yōu)化結(jié)果的蒙特卡羅打靶落點均值對比

表5 不同權(quán)重情況下優(yōu)化結(jié)果的蒙特卡羅打靶落點方差對比

為更直觀地表示上述蒙特卡羅打靶落點分布情況，圖4與圖5分別給出了ω=0和ω≠0時的落點分布3σ誤差橢圓及落點散布，與表4、表5中的結(jié)果相對應(yīng)。由圖4、圖5可知，ω≠0時落點分布范圍明顯減小，說明彈道抗干擾能力明顯增加。但隨著權(quán)重系數(shù)的進(jìn)一步增加，落點散布橢圓幾乎不變。這是由于打靶過程使用開環(huán)控制，沒有加入反饋信號，不確定因素的存在使得落點散布的減小有一個極限。

圖4 不同權(quán)重情況下優(yōu)化結(jié)果蒙特卡羅打靶落點3σ誤差圓對比Fig.4 Comparison of 3σ position error ellipses based on optimization results under different weights

(a) ω=0

5 結(jié) 論

本文將線性協(xié)方差分析方法和高斯偽譜法相結(jié)合，建立了一種考慮不確定性因素的彈道魯棒優(yōu)化模型及算法。對魯棒優(yōu)化彈道結(jié)果進(jìn)行蒙特卡羅打靶驗證，得到了以下結(jié)論：

1)當(dāng)優(yōu)化目標(biāo)中考慮落點精度時，蒙特卡羅打靶的落點分布范圍明顯減小，說明本文優(yōu)化方法能明顯提高彈道的抗干擾能力；

2)由于打靶過程使用開環(huán)控制，沒有加入反饋信號，不確定因素的存在使得誤差減小有一個極限。因此，進(jìn)一步增加目標(biāo)函數(shù)中的權(quán)重系數(shù)對于打靶結(jié)果改善有限。