劉安民，高峰，張青斌，袁天保，張國斌

(1.國防科技大學 空天科學學院，湖南 長沙 410073；2.96796部隊，寧夏 銀川 750000；3.96901部隊，北京 100095)

0 引言

翼傘精確空投系統(tǒng)不僅具有普通降落傘彈道系數(shù)小、可折疊的優(yōu)異特性，還具備更好的滑翔能力和可操控性，并且通過一定的制導控制，還能夠?qū)崿F(xiàn)精準著陸。因此，翼傘系統(tǒng)在部隊突擊、補給物資定點投送等領域都被廣泛運用[1]，如加拿大MMIST公司研發(fā)的Snow Goose空投系統(tǒng)，美國Atair Aerospace公司研制的自主導航空投系統(tǒng)、美國Strong Enterprises公司研制的Screamer空投系統(tǒng)，以及加拿大米斯特機動綜合系統(tǒng)技術(shù)公司研制的Sherpa系統(tǒng)等[2]。但是，飛行環(huán)境和系統(tǒng)參數(shù)等不確定性因素給翼傘空投系統(tǒng)的精確歸航帶來挑戰(zhàn)，阻礙翼傘精確空投系統(tǒng)著陸精度的進一步提升，甚至會導致空投任務的失敗[3]。在軌跡規(guī)劃時進行不確定性優(yōu)化可以增強翼傘系統(tǒng)歸航策略的魯棒性、提升歸航性能[4]，而對飛行過程進行不確定性分析并研究翼傘系統(tǒng)對諸多不確定性因素的響應是整個優(yōu)化過程的基礎。因此，對翼傘系統(tǒng)飛行的不確定性分析具有重要意義。

工程應用中已經(jīng)發(fā)展了大量系統(tǒng)不確定性分析的方法，其中基于大樣本概率統(tǒng)計的Monte Carlo方法是應用最廣泛的方法，其可靠性在工程實踐中得到充分的驗證，但是該方法的可靠性依賴于龐大的樣本數(shù)量，當模型參數(shù)較為復雜時，存在計算量過大、效率低下的缺點。多項式混沌展開(PCE)是不確定性分析中較新的一種方法，該方法主要基于多項式混沌理論，具有堅實的數(shù)學基礎和良好性能，在實踐中也逐漸得到了廣泛的運用。PCE方法最初在研究湍流問題中被Wiener[5]提出，這也被認為是PCE方法的起源，之后該方法也在計算流體力學領域被廣泛運用。Wu等[6]利用非干涉式多項式混沌的方法對跨聲速空氣動力學的不確定性和靈敏度進行分析，研究了馬赫數(shù)和翼型的不確定性對氣動特性的影響，并且還利用稀疏網(wǎng)格法提高樣本利用率，發(fā)現(xiàn)該方法比干涉式多項式混沌方法和Monte Carlo方法具有更高的效率。Prabhakar等[7]將PCE用于分析彈道系數(shù)、大氣密度和升阻比等隨機因素對高超聲速飛行器軌跡的影響，證明了PCE方法可用于高超聲速飛行中的不確定性研究，并且也得到該方法比Monte Carlo方法具有更高效率的結(jié)論。此外，PCE方法在翼型繞流[8]、火星著陸設計[9]和飛行器軌跡優(yōu)化[10]等領域都得到運用，充分證明了其可靠性和計算效率。目前對翼傘空投系統(tǒng)的不確定性分析研究較少。

本文針對某小型載彈翼傘飛行中的氣動參數(shù)不確定性問題，利用翼傘系統(tǒng)9自由度動力學模型，研究了PCE方法在翼傘系統(tǒng)飛行不確定性分析中的應用?？紤]計算精度和效率，分析了多項式階數(shù)和樣本點數(shù)量對PCE模型計算結(jié)果的影響，得到適用于翼傘飛行不確定性分析的PCE模型。此外，通過與傳統(tǒng)Monte Carlo法對比，驗證了PCE模型的計算精度和可靠性。

1 翼傘系統(tǒng)多體動力學建模

1.1 模型描述

翼傘系統(tǒng)主要由傘衣、傘繩、吊帶和荷載4部分組成。在穩(wěn)定飛行階段，傘衣充氣后具有一定剛度而保持特定形狀，并且傘繩和吊帶材料強度大、變形小，故可以將傘衣和傘繩記為剛體P，質(zhì)量為mP,將荷載和吊帶記為剛體L,質(zhì)量為mL.將地面記為零剛體，且零剛體質(zhì)心坐標系與慣性系固連。剛體P具有質(zhì)心平動和繞質(zhì)心轉(zhuǎn)動6個自由度,剛體L通過球鉸Oc與剛體P連接,相對剛體P有3個轉(zhuǎn)動自由度,建立如圖1所示的翼傘系統(tǒng)9自由度多體動力學模型。

圖1 翼傘系統(tǒng)9自由度模型

1.2 受力分析

假設傘衣氣動壓心與其幾何中心重合，位于弦向距離前緣1/4處[11]，則傘衣氣動力FP和氣動力矩MP計算分別為

(1)

(2)

式中：ρ為氣體密度；va是翼傘相對大氣的運動速度；SP為翼傘的參考面積；CD、CL和CY分別為傘衣受到的氣動阻力、升力和側(cè)力系數(shù)，Cl、Cm和Cn為傘衣受到的滾轉(zhuǎn)、俯仰和偏航氣動力矩系數(shù)，用多項式形式[12]表示為

(3)

(4)

荷載為帶螺旋槳的剛體，為簡化計算，不考慮其氣動力矩和螺旋槳對荷載氣動力的影響，荷載所受氣動力記為FL，其值為

(5)

式中：vL為荷載質(zhì)心速度；SL為荷載阻力面積；CDL為荷載阻力系數(shù)。

除氣動力，荷載還受螺旋槳推力Ft，兩剛體重力GP和GL分別作用于各自質(zhì)心位置。

1.3 動力學方程

傘衣充氣后平均密度與空氣密度相近，所以必須考慮附加質(zhì)量效應。本文采用Barrows計算附加質(zhì)量的方法[13]：在平直翼附加質(zhì)量計算方式的基礎上，考慮翼傘沿展向圓弧型的彎曲影響，計算得到翼傘的平動附加質(zhì)量矩陣與繞各主慣性軸的轉(zhuǎn)動附加質(zhì)量矩陣，并將附加質(zhì)量平動分量和轉(zhuǎn)動分量寫成矩陣形式Ma和Ja.在傘體坐標系下，相對于坐標原點Oc的附加質(zhì)量慣量張量矩陣為Ja,Oc.

引入9自由度模型下的系統(tǒng)廣義質(zhì)量矩陣Am及廣義力矩陣Bf：

Am=

(6)

(7)

式中：ΩP和ΩL分別為剛體P和剛體L的角速度；IP和IL分別為兩剛體相對質(zhì)心的慣量張量；vOc為鉸點速度；E3×3為3階單位矩陣。

將傘體坐標系原點建立在鉸點Oc處，避免了翼傘和荷載兩體之間相對約束力的出現(xiàn)。在慣性坐標系下，以鉸點Oc的位置坐標XOc=[xOc,yOc,zOc]T、剛體P姿態(tài)角ψP=[φP,θP,ζP]T和剛體L的姿態(tài)角ψL=[φL,θL,ζL]T為廣義坐標，φ、θ、ζ分別為對應剛體的滾轉(zhuǎn)角、俯仰角和偏航角。根據(jù)Newton-Euler方程，建立翼傘系統(tǒng)的動力學方程為

(8)

由于翼傘和荷載兩體是通過鉸點Oc聯(lián)結(jié)，補充Oc點在飛行過程中位移的運動學關系:

(9)

式中：P=[xP,yP,zP]T為Oc點在傘體坐標系下的位置坐標。

建立的9自由度翼傘系統(tǒng)多體動力學模型已得到研究驗證，證明了該模型的有效性[11-12]。但是，風洞試驗或計算流體力學方法得到的氣動參數(shù)與真實飛行數(shù)據(jù)存在一定誤差，翼傘系統(tǒng)動力學模型中的氣動力系數(shù)存在不確定性[14]。因此，本文將翼傘系統(tǒng)的氣動參數(shù)偏差作為隨機輸入，假設真實氣動力系數(shù)和氣動力矩系數(shù)均為服從均勻分布的獨立隨機變量，標準差都為0.10[3].

2 PCE方法

本文采用非干涉式的PCE方法，該方法可以分為建立PCE模型、求解PCE系數(shù)和概率特性分析。

2.1 模型描述

PCE的原理是將受各分量相互獨立的多維隨機變量ξ=[ξ1ξ2…ξn]影響的系統(tǒng)輸出表示為所有隨機變量的多項式加權(quán)線性組合的形式[15]，將隨機變量的個數(shù)n稱為PCE的維數(shù)，對于隨機響應Y(ξ)可以展開為

(10)

(11)

根據(jù)隨機變量分布類型，從Wiener-Askey系列超幾何正交多項式中選取相應多項式作為基函數(shù)，從而確保所得PCE方法的收斂速度，如均勻分布基函數(shù)為Legendre多項式、正態(tài)分布采用Hermite多項式作為基函數(shù)等。

2.2 PCE系數(shù)的確定

由2.1節(jié)確定的PCE模型中，基函數(shù)根據(jù)隨機變量的分布對應選取，正交多項式由基函數(shù)組合得到，方程中更重要的是確定PCE系數(shù)。目前已經(jīng)發(fā)展多種求解PCE系數(shù)的方法，Galerkin投影法和回歸法是最常用的兩種?；貧w法利用樣本真實響應與PCE模型預測值的相對誤差最小為條件求解系數(shù)[16]，本文采用回歸方法求解PCE系數(shù)。

在標準隨機空間中選取N個樣本點{ξ1,ξ2,…,ξN}，并將標準隨機空間的樣本變換到原隨機空間得到樣本點{X1,X2,…,XN}，在原隨機空間的樣本點上調(diào)用真實響應函數(shù)，得到各樣本點的真實響應值G[17]，則PCE方程在樣本點上可以表示為

Ab=G，

(12)

通常樣本點數(shù)目遠大于待求系數(shù)個數(shù)[18]，則(12)式為超定方程，可以通過線性最小二次回歸方法解得系數(shù)矩陣。對于高維隨機空間，需要的樣本數(shù)量也隨之增加，為避免“維數(shù)災難”則需要通過特殊的采樣方法抽取樣本點進行回歸分析，例如拉丁超立方抽樣、加權(quán)的高斯積分點抽樣以及單項求容積法則抽樣等。本文采用各維均勻分布的拉丁超立方抽樣方法，將各維隨機變量均勻劃分為N個等概率區(qū)間，并在各區(qū)間內(nèi)隨機生成1個1維樣本點，最后將各維度的樣本點隨機組合得到隨機空間中的樣本向量。該方法在降低所需樣本點數(shù)量的同時保證了其均勻性。

根據(jù)線性最小二次回歸求PCE系數(shù)：

(13)

式中：g(Xj)為真實響應函數(shù)值；(ξj)為基于PCE模型的各樣本點的隨機響應函數(shù)值，則最小二次回歸得到系數(shù)為

b=(ATA)-1ATG.

(14)

求得PCE系數(shù)后就可以構(gòu)建完整的PCE模型，相當于構(gòu)建了真實響應的代理模型，輸出量的概率特性(統(tǒng)計矩、概率密度函數(shù)等)可以通過在隨機空間調(diào)用代理模型求得。

3 翼傘飛行不確定性分析

3.1 問題描述

翼傘模型基于某小型載彈翼傘，如圖2所示，翼傘系統(tǒng)基本參數(shù)如表1所示。翼傘系統(tǒng)初始高度為1 000 m，相對氣流具有7.5 m/s的初始水平速度，下落初始速度3.8 m/s，無初始角速度。剛體L相對剛體P初始時刻無相對速度和角速度，初始偏航角為30°，其余姿態(tài)角為0°.翼傘系統(tǒng)飛行總時間為400 s，在滑翔20 s后施加轉(zhuǎn)彎控制。

圖2 某小型載彈翼傘

表1 翼傘基本參數(shù)

3.2 模型階數(shù)與樣本數(shù)量分析

由2.2節(jié)可知，建立考慮6個氣動力系數(shù)隨機性的6維h階PCE模型，一共有s=(6+h)!/(6!h!)項待求PCE系數(shù)，而求解精度依賴于樣本點的選取。模型階數(shù)和樣本點數(shù)量直接影響計算精度，因此針對以上因素展開分析。

為對比各階PCE模型的精度，設置足夠大的樣本點數(shù)量(500個樣本點)，分別計算各階模型結(jié)果的均值和標準差，并統(tǒng)計h階模型計算結(jié)果相對h-1階模型的變化值(h=3,4,5)。圖3給出了各階模型相對變化值的結(jié)果，可以看出均值和標準差的相對變化值隨著階數(shù)升高而逐漸收斂。4階PCE模型相對于3階PCE模型的結(jié)果相對變化值達到收斂：位移均值相對變化值在0.042 m以內(nèi)，標準差的相對變化值也小于0.74 m；速度均值相對變化值低于0.003 6 m/s，標準差相對變化值在0.055 m/s量級；偏航角的動力學方程非線性度較低，低階模型的相對變化值就達到10-4量級。

圖3 不同多項式階數(shù)的相對變化值比較

實際運用中，在達到收斂后提高模型階數(shù)所帶來的精度提升有限，但是計算量增大，因此需要綜合計算精度和效率來選擇合適的模型階數(shù)。對于本文翼傘系統(tǒng)的不確定性分析可以選擇3階PCE模型。

為研究樣本點數(shù)量對PCE模型精度的影響，設置3階PCE模型的樣本點數(shù)量N分別為s的整數(shù)倍(N=s,2s,3s,4s,5s)，并分別計算翼傘系統(tǒng)飛行過程中結(jié)果均值與標準差的相對變化值。計算結(jié)果表明，如果達到精度要求，則樣本點數(shù)量至少為2s，2s個樣本點相對于s個樣本點的位移、速度均值相對變化值分別高達3.82 m和0.21 m/s，標準差相對變化量分別為0.69 m和0.056 m/s.樣本點數(shù)量大于2s的相對變化值如圖4所示，從圖中可以看出過程中樣本點數(shù)量在3s時，相對變化值就達到收斂，繼續(xù)增加樣本點數(shù)量對于精度提升作用不明顯。

圖4 不同樣本點結(jié)果的相對變化值

通過以上分析，對于翼傘飛行的不確定性研究選擇3階PCE模型、樣本點數(shù)量為3s可以滿足計算精度。

3.3 與Monte Carlo方法對比

基于3.2節(jié)對PCE模型階數(shù)和樣本點數(shù)量的分析，PCE方法可選擇多項式截斷于第3階，均勻選取252個樣本點。Monte Carlo方法則直接在三維隨機空間中生成5 000個樣本點。分別采用PCE方法和Monte Carlo方法計算翼傘系統(tǒng)飛行過程中位移的均值和標準差，并將Monte Carlo方法計算結(jié)果作為參考值，計算PCE方法的相對偏差。

如圖5所示，比較兩種方法飛行過程中位移、速度以及偏航角的均值和標準差，結(jié)果表明兩種方法的結(jié)果一致。圖6給出了PCE方法計算結(jié)果相對Monte Carlo方法的偏差。偏航角的均值相對偏差在一段時間后達到穩(wěn)定且保持在在10-4量級，標準差穩(wěn)定時也低于2.5%；速度均值和標準差變化有周期性規(guī)律，其相對偏差也呈周期變化，其均值相對偏差不超過1.5%，標準差相對偏差也小于4.2%；位移的均值相對誤差在0.025%以內(nèi)且標準差相對偏差也小于5.6%.通過與Monte Carlo方法的對比，證明了PCE方法在翼傘飛行過程中不確定性分析的適用性。

圖5 PCE方法與Monte Carlo方法的均值和標準差比較

圖6 PCE方法與Monte Carlo方法的相對誤差比較

更重要的是，Monte Carlo方法計算了5 000個樣本點，耗時近12.5 h，而PCE方法僅利用252個樣本點，耗時僅0.5 h，二者達到幾乎相同的精度。說明PCE方法在相同精度下的效率遠高于Monte Carlo方法，在翼傘飛行或歸航的不確定性分析和靈敏度分析中，該種方法具有更大的價值。

通過PCE方法建立了翼傘飛行過程不確定性傳播的代理模型，該模型基于簡單多項式的線性組合，避免了對原響應復雜動力學過程的計算，顯著提高了計算效率。利用已經(jīng)建立的PCE模型可以快速計算翼傘飛行過程中的狀態(tài)變量分布。分別利用Monte Carlo方法和PCE方法預測受氣動力系數(shù)不確定性影響的翼傘系統(tǒng)落點，Monte Carlo方法歷時3 h，PCE方法基于已有模型僅耗時0.44 s.圖7為兩種方法計算落點分布密度圖，PCE方法得到的落點分布與Monte Carlo方法結(jié)果一致，落點大部分集中在確定條件下的著陸位置附近，偏離距離越遠分布點越少，Monte Carlo方法結(jié)果的最大誤差距離為125.33 m，PCE方法最大誤差距離為126.83 m.

圖7 兩種方法的落點分布

4 結(jié)論

本文針對翼傘飛行的不確定性問題，利用翼傘系統(tǒng)9自由度動力學模型，建立了考慮氣動參數(shù)不確定性的PCE模型，并分析了PCE模型的階數(shù)和樣本點數(shù)量對計算精度的影響。得出以下主要結(jié)論：

1)多項式截止于3階、樣本點數(shù)量為3s時，既能滿足計算精度要求又具有較高計算效率。

2)利用Monte Carlo方法對PCE方法結(jié)果進行對比驗證，結(jié)果顯示均值相對誤差不超過1.5%，標準差相對誤差也在5.6%以內(nèi)，兩種方法可以達到相同的精度。但是，PCE方法僅調(diào)用252次原響應函數(shù)就可以得到Monte Carlo方法計算5 000次的結(jié)果，前者具有更高的效率。分別用Monte Carlo方法和PCE方法計算翼傘系統(tǒng)落點分布，兩方法結(jié)果一致，最大誤差可以達到120 m以上。PCE方法可以為翼傘精確空投任務的性能評估與優(yōu)化提供一種新思路和理論參考。