翼傘系統5段歸航軌跡優化研究

高峰, 郭銳, 豐志偉, 揭錦亮, 張青斌

(1.國防科技大學 空天科學學院， 湖南 長沙 410073; 2. 96796部隊， 寧夏 銀川 750001)

0 引言

翼傘空投系統在軍事以及民用上都具有很大的應用價值，可用于跳傘救生、災區物資投放定點著陸、航天器回收精確歸航等。作為近年來空投的新型工具，翼傘折疊體積小、質量輕，并且具有優良的滑翔性能和穩定性。結合導航、制導與控制技術，翼傘系統可在一定初始位置偏差和風場等干擾條件下實現定點無損著陸。

翼傘系統的工作過程主要有拉直、充氣、歸航控制和雀降，歸航方法主要有分段歸航法、錐形可行域控制方法和最優控制歸航法[1]，在歸航過程中施加控制操作改變翼傘飛行方向，以保證雀降著陸的精確性和穩定性。翼傘空投著陸點環境復雜，著陸范圍廣，著陸的不確定性大，歸航方法對著陸的成功與否影響較大。針對該問題，Li等[2]對翼傘歸航軌跡進行分段處理，采用定線目標尋的方法有效解決了著陸點范圍問題。滕海山等[3]針對運載火箭助推器一子級無控墜落地面落點散布較大的情況，提出了一種翼傘系統的線目標歸航方法。Gao等[4]采用分段歸航方法設計了翼傘歸航軌跡，并利用基于種群的量子差分進化算法對軌跡優化。Sun等[5]設計了飛行復合尋的法，先對翼傘空軌跡采用經典徑向歸航方法分段，再考慮風場作用，對軌跡跟蹤控制。蒲志剛等[6]提出了翼傘分段歸航控制的方向控制方法，探討了位置控制和方向控制在歸航過程中的優先權問題。胡文治等[7]以耗能最少和落地誤差最小為優化目標，對翼傘采用左側下拉和右側下拉耗能情況進行了比較。Suresh Babu等[8]以耗能最小和落地點誤差最小為目標，利用直接多次射擊法求得最優路徑。Yang等[9]采用遺傳算法求解了翼傘分段歸航的最優軌跡，并以貝塞爾曲線為基礎對翼傘末制導進行路徑規劃，以處理變滑翔比的情況。Zou等[10]針對翼傘精確空投問題，為提高翼傘軌跡跟蹤控制精度，提出一種區間型模糊邏輯自整定參數方法，并通過試驗驗證了方法的有效性和實用性。分段歸航使得實際控制量函數變得簡單，利于工程實現。經典的分段歸航一般分為向心飛行段、能量保持段和著陸段[11]。周靚[12]根據歸航任務需求不同，進一步細分為5段歸航法，在3段歸航基礎上增加了初始準備階段和方向調整階段，使得歸航軌跡更加接近真實情況。

盡管翼傘系統具有成本低、易于工程實現的優點，但其也具有可控能力弱、受環境因素影響大的缺點。因此尋求適合于工程實現的歸航策略一直是翼傘系統的關鍵技術。為解決翼傘系統的實時歸航問題， 本文采用擬坐標法建立翼傘系統高保真度的兩體9自由度動力學模型，并通過空投試驗驗證模型的有效性；在此動力學模型基礎上，利用數值實驗結果獲得翼傘系統穩態飛行特性，給出了翼傘系統飛行的滑翔比和轉彎速率限值，得到了近似4自由度模型；采用5段歸航策略，以第1、第5階段飛行時間和第2、第4階段偏航角速度為優化參數，以5段歸航軌跡為初始條件，以偏航角加速度不超過極限值為過程約束條件，以歸航耗能最少為優化目標，利用偽譜法給出了翼傘最優歸航方案；對直接歸航和5段歸航兩種方法的耗能情況以及控制量函數的難易實現程度做了比較，為空投過程提供理論參考。

1 翼傘系統建模

本文研究對象為某型彈載荷翼傘系統，結構如圖1所示。翼傘系統主要分為傘衣、傘繩、吊帶以及載荷4部分。將傘衣、傘繩整體記為剛體A，將吊帶、載荷整體記為剛體B. 剛體A具有3個平動自由度和3個轉動自由度，剛體B通過球鉸和剛體A連接，系統共計9個自由度。為得到翼傘系統飛行過程的穩定速度、滑翔比、轉彎半徑等特征參數，本文首先建立較精確的9自由度模型。

1.1 坐標系定義

如圖2所示，將發射坐標系(慣性系)記為OIxIyIzI，傘衣坐標系記為OPxPyPzP，剛體A體坐標系記為OAxAyAzA，載荷坐標系記為OSxSySzS，剛體B體坐標系記為OBxByBzB，剛體A和剛體B連接處為鉸點O，速度坐標系記為OVxVyVzV.

圖2 翼傘系統坐標系示意圖Fig.2 Schematic diagram of coordinate system of parafoil

記φ、θ和ψ分別為剛體A在慣性系下的滾轉角、俯仰角和偏航角，p、q、r為對應在慣性系下的角速度；φr、θr和ψr分別為剛體B體坐標系相對剛體A在慣性系下的滾轉角、俯仰角和偏航角，pr、qr、rr為對應的相對角速度；α、β分別為翼傘攻角和側滑角；ΓP為傘衣安裝角，ΓS為載荷安裝角。則從慣性系到剛體A體坐標系轉換矩陣為EAI(φ,θ,ψ)；從剛體A體坐標系到剛體B體坐標系轉換矩陣為EBA(φr、θr、ψr)；從剛體A體坐標系到傘衣坐標系轉換矩陣為EPA(ΓP)；從剛體B體坐標系到載荷坐標系轉換矩陣ESB(ΓS)；從傘衣坐標系到速度坐標系轉換矩陣EVP(β,-α)。為便于分析，本文假設傘衣安裝角ΓP和載荷安裝角ΓS均為0°，則EPA=ESB=I，其中I為3階單位矩陣。

1.2 9自由度多體動力學建模

翼傘在飛行過程中受力主要包括自身重力G、氣動力Fa以及發動機推力Ft. 圖2所示剛體A、剛體B兩體模型中，系統所受合外力在剛體A體坐標系下表示為

(1)

式中：GA、GB分別為剛體A、剛體B所受的重力；FaA、FaB分別為剛體A、剛體B受到的氣動力；FPM為傘衣的附加質量對傘衣質心的作用力；Ft為作用于載荷的螺旋槳推力。

傘衣氣動力系數CL、CD、CY和力矩系數Cl、Cm、Cn采用多項式模型描述，即

(2)

(3)

式中：δs為對稱舵偏；δa為反對稱舵偏；vPV為傘衣質心相對氣流的速度；b、c為傘衣展長和弦長，其他參數k、CL0、CLα、CLδS、CLδa、CD0、CDδS、CDδa、CYβ、CYp、CYr、CYδa、Clβ、Clp、Clr、Clδa、Cm0、Cmα、Cmq、CmδS、Cmδa、Cnβ、Cnp、Cnr、Cnδa均為氣動力常系數。

作用在傘衣質心的氣動力為

(4)

式中：S為傘衣參考面積。

剛體A所受合力矩在其體坐標系下表示為

(5)

式中：rAM為在剛體A體坐標系下剛體A的質心到鉸接點之間的矢徑；rMB為在剛體B體坐標系下鉸接點到剛體B質心的矢徑；rBS為在剛體B體坐標系下剛體B質心到載荷質心矢徑；MaA、MPM分別為剛體A受到的氣動力矩、傘衣附加質量在剛體A質心處產生的力矩。氣動力矩表示形式為

(6)

式中：ρ為空氣密度。

剛體B所受合力矩在其體坐標系下表示為

(7)

載荷氣動力形式為

(8)

式中：v為載荷質心速度；Sd為載荷參考面積；Cd為載荷阻力系數。

為便于分析建模，選取翼傘系統兩體模型的擬坐標為

(9)

式中：RA、θA、θBA分別為剛體A在慣性系下的質心位移、姿態角以及剛體B相對于剛體A的姿態角，具體形式為

(10)

相應的，選取擬速度為

(11)

式中：vA為剛體A在其體坐標系下的質心速度，ωA為剛體A在其體坐標系下的角速度，ωBA為剛體B相對于剛體A的角速度，具體形式為

(12)

則翼傘系統擬坐標形式拉格朗日動力學方程為

(13)

式中：L為拉格朗日函數，對擬速度ω偏導為

(14)

(15)

(16)

C為慣性系到剛體A體坐標系的坐標轉換矩陣，DA、DBA分別為剛體A和剛體B的歐拉角時間導數到其體角速度的坐標轉換矩陣；廣義力矩陣Q*表示為

(17)

系統的運動學方程為

(18)

建立9自由度動力學模型的目的是為了開展精確的仿真，以便確定更接近于實際情況的翼傘穩定飛行水平速度、穩定滑翔比以及偏航角速度限值等參數。為保證建立模型的準確性，對系統進行空投試驗，選擇試驗的當天天氣狀況良好。空投環境風場在發射坐標系下水平方向分解為xI軸和yI軸兩個方向，分別記為wxI=-0.5 m/s，wyI=0.25 m/s. 繪制試驗與仿真的軌跡坐標隨時間變化曲線，如圖3所示。

圖3 試驗與仿真坐標位置隨時間變化曲線比較Fig.3 Experimental and simulaed position coordinates over time

試驗與仿真的xI軸、yI軸、zI軸方向坐標隨時間變化曲線趨勢一致，在xI軸和yI軸方向曲線偏差一部分是仿真風場與實際風場差別引起，一部分是所建立模型與真實模型的差異導致。3條飛行曲線與仿真曲線總體趨勢一致說明了該模型參數的有效性。

利用上述翼傘參數，采用表1數據作為初始條件，對9自由度翼傘系統仿真，得到軌跡三維曲線如圖4所示，得到偏航角速度隨時間變化如圖5所示。考慮到實際情況，翼傘下拉偏量須有一定余量，故本文取翼傘轉彎過程下拉偏量最大為總偏量的60%.

表1 9自由度模型仿真初始條件

圖4 翼傘9自由度系統模型仿真三維軌跡Fig.4 3-dimensional trajectory of 9-degree-of-freedom parafoil system model simulation

圖5 60%下拉偏量下偏航角速度隨時間變化曲線Fig.5 Variation curve of yaw angle velocity with time under 60% pull-down deviation

翼傘在t為30～380 s時段轉彎，對稱舵偏為0，反對稱舵偏取最大60%，由圖4可知翼傘系統的最小轉彎半徑Rmin=55 m，由圖5可知穩定飛行時的最大偏航角速度ωmax=0.106 2 rad/s. 仿真得到穩定飛行水平速度大小為5.894 m/s，穩態滑翔比為2.3.

1.3 4自由度質點模型

為便于工程中進行實時歸航設計，需要將翼傘系統簡化為4自由度質點動力學模型[8,14-15]，該模型包括3個平動自由度和1個轉動自由度，如(19)式所示：

(19)

式中：vh為翼傘飛行水平速度大小，vh=5.894 m/s；ψs等價于偏航角；ω等價于偏航角速度大小，-ωmax≤ω≤ωmax；u等價于(偏航角加速度)控制量，-umax≤u≤umax，umax為偏航最大角加速度，由控制系統電機驅動，umax=0.1 rad/s2；垂直速度vzI與水平速度vh滿足

(20)

LD為滑翔比，LD=2.3.

控制量u和轉彎半徑R以及飛行水平速度vh之間滿足關系式：

(21)

2 偽譜法軌跡優化

2.1 5段歸航法簡介

為了能夠進行實時歸航計算，本文采用基于4自由度質點動力學模型的翼傘系統分段歸航策略，該方法所需要考慮的參數較少，且每一段軌跡對應的控制量為常數。如圖6所示為5段歸航示意圖，已知初始點OI(xIOI,yIOI,zIOI)，初始偏航角大小ψ1，水平速度大小vh，豎直速度大小vzI，目標點F(xIF,yIF,zIF)，著陸偏航角大小為ψF.

在第1階段，翼傘從歸航軌跡的起始點O滑翔飛行到盤旋區域上空P點，該過程對翼傘不施加控制；在第2階段，翼傘消耗多余的高度，從P點開始盤旋下降直到高度適中再從C點飛出盤旋區，該過程對翼傘施加單側下拉控制；在第3階段，翼傘從盤旋區C點飛出，經過一段時間滑翔，到達目標點附近區域D點，該過程對翼傘不施加控制；在第4階段，外部風場干擾使得翼傘飛行方向偏離目標點，需要進行180°轉彎以校正方向，之后到達E點，進入最終的逆風雀降段，該過程對翼傘施加單側下拉控制；在第5階段，翼傘逆風雀降著陸，該過程對翼傘施加雙側下拉控制。在本文中，對翼傘在轉彎過程的控制，只采用單一的左側下拉方式。相比于左右兩側交替下拉，單側下拉方法對能量的損耗較小，同時也可以降低翼傘系統的不穩定性[16]。

2.2 5段歸航解析解

如果不考慮最優軌跡優化問題，給定轉彎半徑，就可以獲得5段歸航法的解析解。由于歸航過程zI軸方向速度為常數，因此zI坐標可由時間t唯一確定，故而此處不作求解，進而可僅對二維平面軌跡求解，如圖6所示， 從2.1節5段歸航介紹知，第2階段和第4階段轉彎，設轉彎半徑分別為R2、R4，第1、第3、第5階段為直飛。設第i階段飛行時間為Δti，i=1,2,3,4,5. 根據歸航軌跡特性，只需要確定第1段飛行時間Δt1、第5段飛行時間Δt5、R2和R4，整條軌跡即可解析求解得到。圖6中，vw為風速，由幾何關系得ψ′2=ψ3=ψ4，ψ1=ψ2.

圖6 翼傘系統5段歸航示意圖Fig.6 Schematic diagram of five segments homing of parafoil system

求解歸航過程特征點P、C、D、E、C2、C4以及特征角度ψ2、ψ′2、ψ3、ψ4如下：

特征點P坐標為

(22)

特征點C2點坐標為

(23)

特征點E坐標為

(24)

特征點C4點坐標為

(25)

特征點C坐標表示為

(26)

特征點D坐標表示為

(27)

向量CD與xI軸正方向的夾角為

(28)

第2階段盤旋圈數n可由(29)式確定：

(29)

給出整個歸航軌跡的初始點和終點的坐標及偏航角、水平速度以及滑翔比，第1、第5階段的飛行時間，如表2所示，其中Rmin=55 m. 從初始點開始經過Δt1=20 s后開始單側下拉以調整前進方向和消耗高度，同時為保證雀降能順利實施，要求最終雀降時間Δt5不小于10 s.

表2 初始條件

由表2所給條件，結合圖6，依據(22)式～(28)式，轉彎半徑取R2=R4=100 m，所有角度大小限制在0～180°，求解得到5段歸航軌跡的特征點及特征角度如表3所示。

需要說明的是，上述計算中沒有考慮能量消耗情況，并非最優解。下文將在考慮耗能最小情況下對軌跡進行優化求解。

表3 軌跡解析解

2.3 軌跡最優化求解

高斯偽譜法將狀態量(xI,yI,zI,ψ,ω)和控制量u進行離散，然后在離散點上構造拉格朗日插值多項式，采用適當的非線性規劃得到最優解，在處理含初始和終端約束的問題上具有較大優勢[17-18]。

采用高斯偽譜法進行軌跡優化，設定優化目標為耗能最小，即

(30)

式中：t0和tf分別為優化的初始和結束時間點。

設控制約束為

(31)

5段歸航代數約束方程如(19)式～(21)式。

為了降低偽譜法中非線性規劃的迭代次數，提高計算效率，將表3中的軌跡特征數據作為偽譜法優化的初始解。此外，將5段歸航法和直接歸航法進行對比。直接歸航法以起始點坐標和偏航角為初始條件，以落地目標點坐標和偏航角為終端約束，同樣以偏航角速度、角加速度不超過最大值為過程約束，通過優化得到一條滿足條件的參考軌跡。相比于分段歸航法，直接歸航法不考慮中間過程，只要軌跡滿足約束條件即可。

仿真得到優化后的軌跡三維曲線如圖7所示。在5段歸航模式下，翼傘第1階段飛行20 s下降到948.7 m開始盤旋下降，經過255.6 s盤旋兩圈下降到293.80 m的高度，之后朝著目標點飛行64.60 s到達方向修正起始點D，再經過39.8 s完成方向修正，然后逆風向目標點飛行，再經過10 s后雀降著陸到達目標點F. 在直接歸航模式下，翼傘從初始點開始經過一系列轉彎過程到達目標點。三維軌跡地面投影如圖8所示，翼傘最終著陸方向沿著風場反方向。考慮實際偏航角速度不能突變，本文設其按照一定速率線性增大或減小，如圖9和圖10所示。第2階段偏航角速度從0 rad/s開始以0.006 rad/s2的加速度增大到0.060 46 rad/s后穩定，此時轉彎半徑為97.5 m，而后又以-0.001 7 rad/s2的加速度減小到0 rad/s；第4階段先以0.006 5 rad/s2的加速度增加到0.077 63 rad/s穩定，此時轉彎半徑為76 m，而后以-0.006 4 rad/s2的加速度減小到0 rad/s. 整個過程中偏航角速度連續線性變化，偏航角加速度通過電機操控在一定小范圍內可突變，在工程上切實可行。

圖7 歸航三維軌跡曲線Fig.7 Homing 3D trajectory curves

圖8 歸航軌跡在水平面投影Fig.8 Homing trajectory projected on the horizontal plane

圖9 角速度隨時間變化曲線Fig.9 Changing curves of angular speed with time

圖10 控制量隨時間變化曲線Fig.10 Changing curves of control quantity with time

將5段歸航軌跡作為偽譜法優化的初始解，利用MATLAB軟件運算耗時約4.5 s左右，而僅給出初始點和目標點然后進行5段優化，偽譜法需要在很大范圍內尋找最優解，耗時約在30 s附近，且耗時、能否得到最優解均與設定范圍大小有關。因此，對翼傘系統分段歸航而言，預先求解軌跡的特征點，可在一定程度上較大地提高運算效率。

直接歸航法和給出特征點的5段歸航法兩種方法在優化過程中耗能如表4所示。以偏航角加速度在時間上的積分為耗能參考，則5段歸航法耗能約為直接歸航法的50倍。對圖9和圖10進行比較分析可得：直接歸航法過程簡單，耗能少，但控制量u隨時間變化復雜，實現難度大；5段歸航過程總耗能較大，但每段的控制量均為常數，便于實際操控。在實際應用中，可根據翼傘系統能達到的最小控制精度以及耗能要求合理選擇歸航方案，以求在精度能達到的情況下耗能最少。

表4 兩種歸航方法耗能比較

3 結論

本文建立了翼傘系統擬坐標形式的兩體9自由度動力學模型方程，通過飛行試驗驗證了模型的有效性；采用數值方法仿真求解了系統滑翔比及最大轉彎速率，進而得到翼傘系統簡化的4自由度模型特征參數；考慮5段歸航方法，并將5段歸航軌跡特征點作為偽譜法優化的初始解，以耗能最小為目標求得最優軌跡，最后將直接歸航軌跡及耗能情況與之作比較。得出如下主要結論：

1)所建立的9自由度模型真實有效，可作為翼傘系統數值分析的模型參照。

2)翼傘4自由度模型可由復雜高自由度模型適當簡化得到。

3)與直接歸航法的軌跡優化比較，將5段歸航法軌跡的特征點作為偽譜法優化初始解可顯著提高計算效率。

4)5段歸航偽譜法優化軌跡相比于直接歸航方法耗能較多，但其每段的控制量均為常數，在能耗允許條件下5段歸航法更具有工程實用價值。