基于Copula模型的機翼結冰飛行風險定量評估方法研究

段效聰, 徐浩軍, 張久星, 易 賢, 宣遠勃

(1.空軍工程大學航空工程學院，西安，710038；2.93756部隊，天津，300131;3.中國空氣動力研究寫發展中心結冰與防除冰重點實驗室，四川綿陽，621000)

飛機遭遇結冰時容易引發諸多不利因素(氣動特性惡化、飛行性能下降、空速管等外置傳感器讀數失真等)的耦合連鎖反應，導致人-機-環系統失去穩定性，飛行員操縱不當易引起飛行事故。波音公司統計結果表明，飛機結冰是誘發飛行失控(LOC-I)嚴重事故的三大因素之一。

飛機結冰后的飛行安全保障是當前研究的熱點問題。目前最有效的方式之一是通過開展風險評估，提前預測潛在的結冰區域[1-2]對飛行風險影響的嚴重程度，進而達到規避結冰高風險區域，保障飛行安全的作用[3]。當前學術界有2種典型的結冰條件下風險預測方法，一種是通過冰形積聚情況來對結冰區潛在飛行風險進行預測[4-6]。另一種是通過結冰后的典型空氣動力學特性的改變來判斷結冰條件下的風險[7-8]。然而外界條件、積冰冰形、氣動參數的改變、飛機氣動布局的差異、駕駛員操縱等諸多影響飛行安全的因素最終都要體現在人-機-環復雜系統動力學特性中。因此本文基于人-機-環復雜系統的動力學仿真研究機翼結冰條件下飛行風險的定量評估問題。

飛機遭遇結冰后的飛行動力學特性[9-10]計算是研究飛機結冰致災風險的前提。如：Bragg教授基于經驗估算提出了結冰參數模型[11]，并開發了智能防冰系統[12-13]，有效提升了結冰自動邊界保護和結冰自適應飛行控制能力。NASA開展了平尾結冰項目(TIP)，研究了平尾結冰對空氣動力學及飛行動力學特性的影響[14-15]。

針對機翼結冰后飛行風險定量評估問題，本文結合蒙特卡羅飛行仿真實驗與Copula理論，對結冰條件下人-機-環復雜系統風險定量評估方法進行了研究。對機翼非對稱結冰情形下的飛行風險進行了定量計算。所提出的方法完善了現有的飛行安全分析手段，并為結冰條件下飛行風險的指示與規避提供參考與借鑒。

1 機翼結冰后人-機-環動力學建模

機翼結冰使得飛機氣動特性發生變化，失速迎角降低，駕駛員的操縱策略、飛機力學特性與環境因素相互影響，考慮到三者之間的相互作用，本文基于飛機運動模型、結冰條件下典型駕駛員模型、結冰影響模型構建了人-機-環復雜系統動力學模型。其中飛機運動模型采用基于四元數法構建的六自由度模型，具體可參考文獻[16]。

1.1 結冰條件下典型駕駛員模型

飛機遭遇機翼結冰，影響飛行狀態時，駕駛員為恢復姿態進行一定的精準補償操縱，如圖1所示，根據飛機響應輸出Y(s)與操縱指令R(s)之間的系統誤差E(s)進行修正。駕駛員操縱過程可以分為3個階段：①經過時間t1發現飛機結冰和飛行姿態改變；②操縱駕駛桿至狀態參數反向變化，操縱過程時間為t2；③保持飛機姿態穩定。結冰條件下駕駛員輸入輸出模型可表示為式(1)，其中F(s)與G(s)為不同的McRuer模型[17]。

圖1 駕駛員補償操縱模式下的人-機閉環系統

(1)

1.2 結冰影響模型

Bragg教授基于經驗估算和風洞試驗，提出了結冰影響參量模型。結冰前后的氣動導數關系如式(2)所示：

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

(2)

式中：C(A)為任一氣動導數，C(A)iced為相應的結冰后氣動導數；kCA為飛機結冰因子，對于特定飛機，其為定值；ηice為氣象因子，通常取值為0～0.3[11]。

2 基于蒙特卡羅仿真實驗提取飛參極值

2.1 蒙特卡羅仿真流程

研究結冰條件下內外部影響因素的復雜隨機性需要的數據量較大，進行真機自然結冰條件的飛行試驗難度高，通過蒙特卡羅仿真實驗[18]提取相應的飛行參數極值成為有效的解決途徑，仿真實驗流程如圖2所示。具體步驟如下：

1)設定仿真初始狀態的參數。

2)在機翼結冰數據庫中根據飛行階段和任務提取結冰模型關鍵參數，并計算對飛機氣動導數和氣動力的影響，同時飛機運動狀態的改變亦反饋到結冰模型，提取下一時間的結冰參數。

3)利用蒙特卡羅方法將其他內外部影響條件變量按照其出現頻率進行隨機抽樣。

4)將抽樣的變量數值作為檢索條件提取飛行員操縱行為特性參數及其他影響飛行狀況的條件數據。

5)進行第i次飛行仿真計算。

6)通過系統中總線實時提取多維飛行參數，辨識極值點并存入數據庫。

7)本次計算結束后，返回步驟2)開啟下一次計算過程，直至迭代循環完成特定點上全部計算次數。

2.2 飛行安全關鍵參數確定方法

通過選取一次典型的飛機遭遇結冰后駕駛員進行操縱的仿真案例來說明飛行安全關鍵參數的確定方法。背景飛機選取與A320、波音737和C919為同一量級的典型大型客機。仿真背景設定為：在3 000 m高度，120 m/s速度平飛，在仿真開始5 s后遭遇機翼結冰，設定右側機翼除冰系統故障。圖2所示的是第21次仿真過程中，飛行狀態參數的變化情況。

圖2 第21次飛行仿真中飛行狀態參數變化情況

從圖2可以看出，在機翼結冰發生告警之前，由于升力降低、阻力增大，飛機空速略有降低；右側機翼除冰系統故障無法除冰，飛機在不對稱升力和阻力作用下，滾轉力矩致使飛機向右側滾轉產生側滑，并開始急劇掉高度。駕駛員采取增加油門，反向壓桿的操縱策略，飛機逐漸恢復到初始狀態。在這一過程中，飛機迎角、滾轉角、飛行空速和高度變化較大。由于失速迎角降低、橫向鉸鏈力矩增大，飛行員難以控制航向穩定。基于此文中選擇變化劇烈且對飛行安全較為敏感的迎角、滾轉角和飛行空速作為飛行安全關鍵參數。

2.3 基于一維飛參極值的機翼結冰飛行風險判定條件與量化評估

通過對得到的3個飛行參數極值樣本進行統計學分析，發現均具有厚尾分布的特征。采用常見的極值分布模型對提取的飛行參數極值的概率分布和數據序列的邊際概率分布尾部進行了建模，通過擬合優度檢驗驗證了3個飛行參數極值樣本均服從式(3)所示的GEV分布[18]模型。

(3)

在仿真實驗中通常利用飛行關鍵參數超出安全許用范圍來表征飛行風險。由于3個參數的取值范圍相差較大，通過歸一化處理便于對比分析。飛機迎角與飛行馬赫數和襟翼偏度狀態有關，考慮結冰對失速迎角的影響，通過αmax/αc(δf，Ma，ice)對迎角極值進行歸一化處理。氣動手冊規定背景飛機滾轉角臨界值為Φmax=85°；飛機空速左邊界值與襟翼偏度和結冰影響有關。

因此結冰條件下飛行風險發生概率可通過式(4)～(6)計算。

(4)

(5)

(6)

3 結冰遭遇情形下三維極值參數Copula模型

結冰導致飛行事故的發生是一個多因素耦合作用下的動態過程，僅僅依靠一維極值參數開展飛行風險定量評估是不全面的，有可能因為參數的選取而導致風險判定的不準確，單一參數的超限不一定導致飛行事故的發生，因此需要綜合考慮多個關鍵參數。

依靠單參數開展研究無法分析參數之間的相依性關系對飛行風險的影響，而多參數條件下的聯合分布形式考慮了參數間的相依性結構，能夠分析表征參數耦合情形下的分布特征和飛行風險，本文采用的Copula模型便是利用了相同的思路。

3.1 三維極值參數辨識

由Sklar定量和Copula函數的基本性質可知，一定存在一個Copula函數C，滿足：

C(F1(Vmin),F2(φmax),F3(αmax))=

F(Vmin,φmax,αmax)

(7)

根據Archimedean族Copular函數性質，構建三維非對稱Copular函數結構：

C(u,v,w)=C1(w,C2(u,v))=

(8)

式中φ(·)為阿基米德Copular函數，u,v,w為最小空速極值、滾轉角極值、迎角極值的邊緣分布，即u=F1(Vmin)，v=F2(φmax)，w=F3(αmax)。φ1和φ2分別對應Copular函數C1和C2，函數φ2由u和v生成，函數φ1由w和φ2生成，假設其參數分別為θ1和θ2，根據式(8)將C2(u,v,θ1)代入C1(w,C2,θ2)即可得到三維Copular函數模型。

對于本文三維極值參數的Copula模型選擇，主要有Gumbel Copula模型、Clayton Copula 模型、GS Copula模型、Joe Copula模型。其母函數和Copula三維模型如式(9)～(12)所示。

(9)

C(u,v,w)=

(10)

(11)

(12)

通過辨識得到上述4種Copula模型中的未知參數如表1所示。

表1 參數辨識結果

3.2 擬合優度檢驗

表2為Copula函數擬合優度檢驗結果。分析表2可得，Clayton模型和Joe模型的P(K-S)值均大于0.05，在置信水平0.01、0.02、0.05的情況下均能有效通過檢驗，而其他模型的辨識精度則無法通過檢驗。進一步比較這2個模型的AIC、BIC和χ2檢驗值，Joe模型的計算值更低，描述精度更高。

表2 Copula函數擬合優度檢驗結果

因四維圖像不易直觀地表征極值參數樣本的尾部特性，選擇w=0.8，進一步給出Joe Copula模型和Clayton模型的概率密度圖，如圖3和圖4所示。Joe模型和Clayton模型具有很強的耦合性，在高尾部分的風險值較高，而其他區域風險值較低。相較于Clayton模型，Joe模型的風險概率分布更加密集，高尾部分的分布密度梯度更大，對后尾特性描述更好，因此選擇Joe模型作為機翼非對稱結冰情形下最小空速極值、滾轉角極值和迎角極值的描述模型。

圖3 Joe模型分布概率密度和累積概率圖

圖4 Clayton模型分布概率密度和累積概率圖

3.3 三維參數極值Copula模型飛行風險概率計算

根據構建的飛行風險判定條件，依據Joe Copula模型，求解飛行風險概率如式(13)所示：

Pr=1-CJoe(F1(Vmin/Vc(δf,ice)>1，

F2(φmax/85>1),F3(αmax/αc(δf,Ma,ice)>1)

(13)

式中：Cjoe代表基于式(12)求解的Joe Copula模型。

將辨識得到的參數帶入式(12)和(13)，求得飛行風險概率為0.062 4，參考MIL-STD-882D標準，飛行風險已經達到B極，此時需要飛行員密切關注飛行狀態的變化。運用Joe Copula函數計算出的風險概率值較高，一定程度上反映了參數耦合特性對飛行風險事件演化的影響，考慮多個飛行參數極值構建聯合概率分布模型，能夠更加全面地分析樣本數據所蘊含的風險信息。

4 結論

1)基于仿真案例確定了機翼非對稱結冰條件下迎角、滾轉角和飛行空速對飛行安全影響最大。利用蒙特卡羅仿真實驗提取多組飛參極值，基于安全許用范圍建立風險判定依據，利用辨識得到的三維極值模型參數進行擬合優度檢驗，并計算出符合該三維極值模型的飛行風險概率。

2)所提出的機翼結冰飛行風險定量評估方法可以定量評估機翼非對稱結冰對飛行安全的影響。由于飛行安全受到多種因素影響，在評估飛行風險時不可能將所有因素均考慮在內。文中所計算的風險值在某種程度上只是一個參考值，但可以作為一種客觀的指標，為分析飛機結冰科目試飛風險提供理論依據，科學地衡量機翼結冰情形對飛行安全影響的嚴重程度。