基于多項式混沌展開方法的翼傘飛行不確定性

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

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

0 引言

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

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

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

1 翼傘系統多體動力學建模

1.1 模型描述

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

圖1 翼傘系統9自由度模型

1.2 受力分析

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

(1)

(2)

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

(3)

(4)

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

(5)

式中：vL為荷載質心速度；SL為荷載阻力面積；CDL為荷載阻力系數。

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

1.3 動力學方程

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

引入9自由度模型下的系統廣義質量矩陣Am及廣義力矩陣Bf：

Am=

(6)

(7)

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

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

(8)

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

(9)

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

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

2 PCE方法

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

2.1 模型描述

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

(10)

(11)

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

2.2 PCE系數的確定

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

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

Ab=G，

(12)

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

根據線性最小二次回歸求PCE系數：

(13)

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

b=(ATA)-1ATG.

(14)

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

3 翼傘飛行不確定性分析

3.1 問題描述

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

圖2 某小型載彈翼傘

表1 翼傘基本參數

3.2 模型階數與樣本數量分析

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

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

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

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

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

圖4 不同樣本點結果的相對變化值

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

3.3 與Monte Carlo方法對比

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

如圖5所示，比較兩種方法飛行過程中位移、速度以及偏航角的均值和標準差，結果表明兩種方法的結果一致。圖6給出了PCE方法計算結果相對Monte Carlo方法的偏差。偏航角的均值相對偏差在一段時間后達到穩定且保持在在10-4量級，標準差穩定時也低于2.5%；速度均值和標準差變化有周期性規律，其相對偏差也呈周期變化，其均值相對偏差不超過1.5%，標準差相對偏差也小于4.2%；位移的均值相對誤差在0.025%以內且標準差相對偏差也小于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方法建立了翼傘飛行過程不確定性傳播的代理模型，該模型基于簡單多項式的線性組合，避免了對原響應復雜動力學過程的計算，顯著提高了計算效率。利用已經建立的PCE模型可以快速計算翼傘飛行過程中的狀態變量分布。分別利用Monte Carlo方法和PCE方法預測受氣動力系數不確定性影響的翼傘系統落點，Monte Carlo方法歷時3 h，PCE方法基于已有模型僅耗時0.44 s.圖7為兩種方法計算落點分布密度圖，PCE方法得到的落點分布與Monte Carlo方法結果一致，落點大部分集中在確定條件下的著陸位置附近，偏離距離越遠分布點越少，Monte Carlo方法結果的最大誤差距離為125.33 m，PCE方法最大誤差距離為126.83 m.

圖7 兩種方法的落點分布

4 結論

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

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

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