禹志龍，李穎暉，鄭無計，徐文豐，董澤洪

空軍工程大學 航空工程學院，西安 710038

飛機結冰是飛行實踐中廣泛存在的一種物理現(xiàn)象[1]，通常會破壞飛機的動力學特性，影響飛機的操穩(wěn)性特性，降低飛行品質，嚴重威脅飛機飛行安全[2]。自20世紀20年代以來，飛機結冰一直是威脅航空安全的重要問題[3]，美國Safety Advisor對1990—2000年間因氣象原因造成的飛行事故進行統(tǒng)計，共有3 230起，其中結冰事故占12%[4]，隨后在2003—2008年間，又發(fā)生了380起因為結冰導致的飛行事故。

針對結冰對飛行安全造成的威脅，國外學者從20世紀30年代開始對結冰現(xiàn)象進行研究。Lynch和Khodadoust深入系統(tǒng)地研究了不同結冰位置、結冰形狀對飛機升、阻力系數(shù)的影響，并量化分析了結冰引起的氣動特性變化[5]。Bragg和Basar在國際結冰會議上提出了進行結冰程度量化的模型[6]，為結冰的量化分析奠定了基礎。國內以研究飛機結冰、防除冰機理[7], 操縱保護，風險量化、評估[8]為重點。

為了降低飛機在結冰條件下的飛行風險，必須建立飛機的飛行安全包線，而傳統(tǒng)飛行安全包線保護通常只用于防止飛行員操縱飛機進入超出飛機結構和氣動極限的狀態(tài)。1994年，美國鷹航公司一架ATR72-212飛機墜毀于印第安納州，當時飛機的飛行迎角僅為5°，遠低于飛控系統(tǒng)提供的迎角限制值18.1°，可是在機翼結冰的影響下，飛機仍然因失速而導致墜毀[9]。因此，必須獲得飛機在復雜結冰環(huán)境下的飛行安全包線。目前，國內外確定飛行安全包線的方法很多，其中獲得廣大研究者公認的方法之一就是可達性分析[10]?？蛇_性分析是指針對微分動力系統(tǒng)，是否可以通過有效控制策略、在特定的時間范圍內使初始集范圍內的狀態(tài)達到指定目標集合的過程。該方法應用范圍較廣，在航空領域安全分析中表現(xiàn)突出。文獻[11]利用可達集方法計算了飛機著陸階段的飛行安全包線，并分析了襟翼切換對于飛機著陸時飛行安全包線的影響，文獻[12]以最大可控不變集作為飛機的安全集，并基于安全集研究了機動邊界的保護控制，文獻[13]基于可達集方法求解飛行器筋斗動作各階段的操縱邊界，為飛行員完成標準動作提供決策支撐。由于可達集理論的實用性，一些學者將可達集理論應于對結冰飛行安全的研究，文獻[14-15]將可達集用于結冰飛機著陸階段的飛行安全研究，并分別得出不同結冰程度下的飛行風險與操縱策略，但是都是在結冰完全已知的情況下得出的結論，具有一定的局限性。

為了分析復雜結冰環(huán)境對飛機飛行安全包線的影響，首先利用時間尺度分離的方法將飛機運動方程解耦成兩個系統(tǒng)，為方便研究選取快子系統(tǒng)為研究對象，以快子系統(tǒng)輸出為虛擬輸入，建立考慮飛機橫向運動(β、φ)的縱向運動方程。然后根據(jù)結冰對飛機氣動特性的影響，建立考慮結冰位置、冰型及分布等不確定性的復雜結冰環(huán)境下的魯棒氣動導數(shù)模型。運用可達集理論，得到飛機在不同結冰程度下的飛行安全包線，對結冰條件下飛機滾轉運動進行分析。最后進一步分析考慮不同結冰程度下結冰位置、冰型及分布等不確定性結冰情況時的魯棒飛行安全包線，由此得出結冰不確定性對飛機飛行安全帶來的影響。

1 結冰飛機動力學建模

為簡化分析過程，以RCAM(Research Civil Aircraft Model)簡化模型為研究對象，飛機所有氣動參數(shù)基于文獻[18]。

圖1 飛機動力學時間尺度分離Fig.1 Time scale separation of aircraft dynamics

1.1 飛機非線性動力學模型

將運動飛機視為一個質點，其受力如圖2所示。圖中：Xb、Yb分別表示機體坐標系下的縱軸及豎軸；Xω、Yω分別表示速度坐標系下的阻力軸及升力軸；α為迎角；θ為俯仰角。

建立飛機非線性動力學方程:

(1)

式中：m為飛機質量；V為飛行速度；γ為航跡傾角；φ為滾轉角；χ為航跡方位角；FX、FY、FZ分別表示沿X、Y、Z軸方向的合力，可表示為

(2)

式中：升力L、阻力D和側力Y可表示為

(3)

圖2 飛機受力質點圖Fig.2 Point mass force diagram for aircraft

其中：ρ為空氣密度；S為機翼參考面積；CL、CD、CY分別為升力、阻力和側力系數(shù)。

(4)

1.2 飛機結冰模型

飛機結冰是指飛機機體表面的某些部位聚集冰層的現(xiàn)象，主要會導致飛機空氣動力學特性變差，升力系數(shù)減小、阻力系數(shù)增大，同時使得飛機最大飛行迎角減小，如圖3所示。

圖3 飛機結冰前后升力系數(shù)變化示意圖Fig.3 Diagram of lift coefficient of aircraft changes in clean and iced

根據(jù)Bragg和Basar提出的量化結冰模型，構建結冰程度對飛機氣動參數(shù)影響的模型[6]：

C(A)iced=(1+ηkC(A))C(A)

(5)

式中：C(A)和C(A)iced分別為飛機結冰前后的某個氣動導數(shù)值；kC(A)為飛機結冰因子，反映氣動導數(shù)受結冰影響程度，通常由飛機本身結構決定,對于特定飛機為一個常值；η為表征飛機結冰程度的參數(shù)，未結冰時η=0，隨著結冰程度的增加，η的值不斷增大。

Bragg的結冰模型雖然能夠定量反映飛機結冰的影響，但是一方面，飛機結冰程度難以準確測量，另外，結冰是一個非常復雜的過程，不同結冰程度、不同冰型、不同結冰位置等，都會對飛機的氣動特性造成不同影響，并不能簡單用η來描述。本文針對結冰難以準確量化問題，根據(jù)結冰對飛機氣動特性的影響，提出魯棒氣動導數(shù)模型,可表示為

(6)

2 可達集及飛行安全包線保護

基于可達集理論，利用可達性分析獲得的反向可達集作為結冰飛機最大飛行安全包線，一旦飛行狀態(tài)超出反向可達集范圍，飛機就有失控的危險。

2.1 可達集理論

飛機的動態(tài)特性可以表示為一個連續(xù)時間控制：

(7)

式中：x∈Rn表示n維狀態(tài)空間；u∈U?Rm表示m維控制量且f(·,·):Rn×U→Rn為有界、連續(xù)的Lipschitz函數(shù)。

定義狀態(tài)軌跡為

ξ(x0,t0,u(·),t):t→x∈Rn

(8)

?tf∈[t,T]|ξ(x,t0,u(·)，tf)∈K}

(9)

式中：K∈Rn為目標集；tf為達到目標集所需的時間。

本文研究飛機在結冰條件下的飛行安全包線，從理論上確定飛機能夠安全飛行的最大可操縱范圍，而通過可達性分析得到的反向可達集剛好能夠求得飛機可以恢復到安全飛行狀態(tài)的最大范圍，因此，選取反向可達集作為復雜結冰條件下的飛行安全包線。

2.2 可達集計算

可達集的求解基本上可以分為兩類，包括拉格朗日法和歐拉法[19]。拉格朗日方法利用緊集表示向量場的流向，因此，計算復雜度較高，一般用于計算較高維的可達集；歐拉法主要是以水平集方法為代表，水平集方法用HJB-PDE(Hamilton Jacobi Bellman-Partial Differential Equation)零水平集方程的黏性解表示可達集[20]，該方法將可達性分析與最優(yōu)控制理論緊密結合起來，能處理復雜的非線性和最優(yōu)控制策略問題[21]，在航空界應用較為廣泛。

2.2.1 水平集方法

利用水平集方程求解時間依賴的Hamilton-Jacobi方程零水平集的黏性解，并將Hamilton-Jacobi方程的黏性解作為可達集的隱式表達式。

水平集方程可以表示為

(10)

式中：φ(x,t)為水平集函數(shù)。

目標集K可由水平集函數(shù)表示為

K={x∈Rn|φ(x,0)≤0}

(11)

通過計算HJB-PDE方程的黏性解，得到目標集K對應的可達集為

(12)

式中：

(13)

最優(yōu)控制輸入為

u*(x,p)=arg maxpTf(x,t,u)

(14)

表示在狀態(tài)x下能使哈密爾頓方程H(x,p)達到最大的控制輸入。

2.2.2 最優(yōu)控制輸入計算

最優(yōu)控制輸入u*(x,p)在飛控系統(tǒng)中的作用是使飛行狀態(tài)x的軌跡始終保持在反向可達集內部，最終在t=τ時刻達到目標集。

由式(13)可得本文動力學模型的Hamilton函數(shù)為

(15)

2.2.3 考慮不確定性的可達集

?tf∈[t,T]|ξ(x,t0,u(·),tf,Δ)∈K}

(16)

HJB-PDE方程可以表示為

(17)

3 結冰飛機的飛行安全包線

以RCAM簡化模型為研究對象，對結冰飛機低空平飛階段的飛行安全包線變化及運動特性進行分析。正常飛行時飛機及飛行相關參數(shù)為：飛機質量m=120 000 kg,重力加速度g=9.81 m/s2,飛行高度H=2 000 m,大氣密度ρ=1.006 5 kg/m3，最大推力Tmax=410 920 N，最小推力Tmin=20 546 N，機翼參考面積S=260 m2,升力系數(shù)CL0=1.065 6、CLα=6.072 5，阻力系數(shù)CD0=0.159 9、CDα=0.503 5、CDα2=2.117 5，側力系數(shù)CY=-1.6。

3.1 結冰對飛行安全包線的影響

首先對確定性結冰進行分析，假設飛機機翼發(fā)生結冰，參照文獻[22]中結冰程度的區(qū)分標準，分別定義3種結冰程度η，即η=0，0.1，0.3，分別對應未結冰、輕度結冰、嚴重結冰。

計算飛機正常情況飛行時在給定目標集下的反向可達集，并以反向可達集作為此條件下飛機的飛行安全包線，如圖4所示，在可達集范圍內，所有飛行狀態(tài)都可在最優(yōu)控制操縱下達到目標集。一旦狀態(tài)超出可達集邊界，飛機將無法恢復到正常飛行狀態(tài)。

圖4 正常條件下飛行安全包線Fig.4 Flight safe envelope in clean condition

計算飛機未結冰(η=0)發(fā)生輕度結冰(η=0.1)、嚴重結冰(η=0.3)時的飛行安全包線，并進行對比，如圖5所示。圖中綠色、粉色和紅色包線分別代表未結冰、輕度結冰及嚴重結冰時飛機飛行安全包線，由仿真圖可以看出，隨著結冰程度的增加，飛機的可操縱性下降，飛行安全包線不斷收縮，在速度軸V方向上收縮最為嚴重。當飛機發(fā)生嚴重結冰時，飛機的氣動參數(shù)發(fā)生改變，最大失速迎角減小，最小失速速度增大，飛行安全包線向右收縮，飛機發(fā)生失速的危險增大。

考慮在結冰情形下飛機滾轉對于飛行安全包線的影響，分別選滾轉角φ=0°,30°,60°的情形進行研究，如圖6所示。圖中綠色、粉色和紅色包線分別代表滾轉角為φ=0°，30°，60°時的飛行安全包線，通過對不同滾轉角下飛行安全包線的比較可以看出，隨著滾轉角φ的不斷增大，飛行安全包線持續(xù)減小，且主要表現(xiàn)在γ軸方向上，飛機的爬升性能降低，這主要因為飛機爬升的力由Lcosφ提供，當滾轉角增大時，爬升力減小，爬升性能降低。因此當飛機發(fā)生結冰時，應當減少滾轉機動，以增大飛機爬升與俯沖的安全操縱范圍。

圖5 不同結冰程度下的飛行安全包線Fig.5 Flight safe envelope under different icing degrees

圖6 不同滾轉角下的飛行安全包線Fig.6 Flight safe envelope at different roll angles

3.2 不確定結冰條件下的飛行安全包線

在3.1節(jié)中討論了在確定性結冰條件下飛機飛行安全包線的變化，但在實際飛行中無法完全獲取飛機結冰信息，因此，也無法獲得飛機準確的飛行安全包線。在此對結冰程度已知，但結冰位置、冰型及分布等存在不確定性情況下的飛機飛行安全包線進行魯棒性分析。

根據(jù)第1節(jié)所建模型，選取高斯分布標準差為3σ(3σ置信度為99%)，不確定量Δ分別為0、20%和30%，對輕度結冰(η=0.1)、嚴重結冰(η=0.3)時的飛行安全包線進行分析，如圖7、圖8所示。

圖7 輕度結冰條件下飛行安全包線的變化Fig.7 Changes of flight safe envelope under mild icing condition

圖8 嚴重結冰條件下飛行安全包線的變化Fig.8 Changes of flight safe envelope under severe icing condition

圖7中藍色曲面和黑色曲面分別表示飛機在輕度結冰(η=0.1)條件下，不確定量Δ分別為20%和30%時的飛行安全包線。由于結冰位置、冰型及分布等存在不確定性，使得氣動導數(shù)變化不確定，當考慮Δ=20%時，魯棒飛行安全包線相較確定性安全包線在φ=0°收縮了30%，當考慮Δ=30%時，魯棒飛行安全包線相較確定性安全包線在φ=0°收縮了45%。

圖8中藍色曲面和黑色曲面分別表示飛機在重度結冰(η=0.3)條件下，不確定量Δ分別為20%和30%時的飛行安全包線。當考慮Δ=20%時，魯棒飛行安全包線相較確定性安全包線在φ=0°收縮了15%，當考慮Δ=30%時，魯棒飛行安全包線相較確定性安全包線在φ=0°收縮了35%。

通過對比可知，結冰位置、冰型及分布等不確定性因素對輕度結冰時的飛行安全包線的影響大于重度結冰時的影響，這主要是由于確定性輕度結冰時飛機的氣動導數(shù)變化相對較小，飛機所處于的狀態(tài)仍然為平衡狀態(tài)，但是當考慮不確定性時，氣動導數(shù)變化增大，原平衡點的性質發(fā)生了改變，由穩(wěn)定的平衡點變?yōu)椴环€(wěn)定平衡點。因此，當考慮不確定性時，輕度結冰條件下飛行安全包線收縮更加嚴重。

4 結 論

建立了考慮橫向運動(β、φ)的飛機縱向運動模型，基于可達集理論利用可達性分析研究了飛機在確定性結冰條件下和考慮結冰位置、冰型及分布等不確定性結冰條件下飛機的飛行安全包線，分析了不同結冰程度(η=0,0.1,0.3)對于飛行安全的影響，得出以下結論：

1) 以反向可達集作為飛行安全包線相較于傳統(tǒng)包線具有明顯優(yōu)勢，并且在外部環(huán)境發(fā)生變化時，仍然能夠指導飛行員操縱。

2) 隨著結冰程度的增加，飛行安全包線的范圍不斷收縮，主要體現(xiàn)在速度軸方向上，飛機的氣動參數(shù)發(fā)生改變，最大失速迎角減小，最小失速速度增大,包線向右收縮。

3) 當考慮滾轉角φ對結冰飛機飛行安全包線影響時，隨著φ的增大，飛行安全包線主要沿著γ軸方向收縮，飛行安全包線減小。因此，當飛機發(fā)生結冰時，應盡量減小滾轉角φ同時避免俯沖機動，以提高飛行安全。

4) 當考慮結冰位置、冰型及分布等不確定性結冰情況下的魯棒飛行安全包線時，不確定量Δ對于飛行安全包線影響較大，且不確定性因素對輕度結冰時的飛行安全包線的影響大于重度結冰時的影響。因此，即使飛機只是發(fā)生輕度結冰但是由于結冰位置、冰型及分布等的不確定性，飛行風險仍然較大，飛行員應該保持警惕，適當縮小操縱范圍。