段效聰, 徐浩軍, 張久星, 易 賢, 宣遠(yuǎn)勃

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

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

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

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

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

1 機(jī)翼結(jié)冰后人-機(jī)-環(huán)動力學(xué)建模

機(jī)翼結(jié)冰使得飛機(jī)氣動特性發(fā)生變化，失速迎角降低，駕駛員的操縱策略、飛機(jī)力學(xué)特性與環(huán)境因素相互影響，考慮到三者之間的相互作用，本文基于飛機(jī)運(yùn)動模型、結(jié)冰條件下典型駕駛員模型、結(jié)冰影響模型構(gòu)建了人-機(jī)-環(huán)復(fù)雜系統(tǒng)動力學(xué)模型。其中飛機(jī)運(yùn)動模型采用基于四元數(shù)法構(gòu)建的六自由度模型，具體可參考文獻(xiàn)[16]。

1.1 結(jié)冰條件下典型駕駛員模型

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

圖1 駕駛員補(bǔ)償操縱模式下的人-機(jī)閉環(huán)系統(tǒng)

(1)

1.2 結(jié)冰影響模型

Bragg教授基于經(jīng)驗估算和風(fēng)洞試驗，提出了結(jié)冰影響參量模型。結(jié)冰前后的氣動導(dǎo)數(shù)關(guān)系如式(2)所示：

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

(2)

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

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

2.1 蒙特卡羅仿真流程

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

1)設(shè)定仿真初始狀態(tài)的參數(shù)。

2)在機(jī)翼結(jié)冰數(shù)據(jù)庫中根據(jù)飛行階段和任務(wù)提取結(jié)冰模型關(guān)鍵參數(shù)，并計算對飛機(jī)氣動導(dǎo)數(shù)和氣動力的影響，同時飛機(jī)運(yùn)動狀態(tài)的改變亦反饋到結(jié)冰模型，提取下一時間的結(jié)冰參數(shù)。

3)利用蒙特卡羅方法將其他內(nèi)外部影響條件變量按照其出現(xiàn)頻率進(jìn)行隨機(jī)抽樣。

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

5)進(jìn)行第i次飛行仿真計算。

6)通過系統(tǒng)中總線實時提取多維飛行參數(shù)，辨識極值點(diǎn)并存入數(shù)據(jù)庫。

7)本次計算結(jié)束后，返回步驟2)開啟下一次計算過程，直至迭代循環(huán)完成特定點(diǎn)上全部計算次數(shù)。

2.2 飛行安全關(guān)鍵參數(shù)確定方法

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

圖2 第21次飛行仿真中飛行狀態(tài)參數(shù)變化情況

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

2.3 基于一維飛參極值的機(jī)翼結(jié)冰飛行風(fēng)險判定條件與量化評估

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

(3)

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

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

(4)

(5)

(6)

3 結(jié)冰遭遇情形下三維極值參數(shù)Copula模型

結(jié)冰導(dǎo)致飛行事故的發(fā)生是一個多因素耦合作用下的動態(tài)過程，僅僅依靠一維極值參數(shù)開展飛行風(fēng)險定量評估是不全面的，有可能因為參數(shù)的選取而導(dǎo)致風(fēng)險判定的不準(zhǔn)確，單一參數(shù)的超限不一定導(dǎo)致飛行事故的發(fā)生，因此需要綜合考慮多個關(guān)鍵參數(shù)。

依靠單參數(shù)開展研究無法分析參數(shù)之間的相依性關(guān)系對飛行風(fēng)險的影響，而多參數(shù)條件下的聯(lián)合分布形式考慮了參數(shù)間的相依性結(jié)構(gòu)，能夠分析表征參數(shù)耦合情形下的分布特征和飛行風(fēng)險，本文采用的Copula模型便是利用了相同的思路。

3.1 三維極值參數(shù)辨識

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

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

F(Vmin,φmax,αmax)

(7)

根據(jù)Archimedean族Copular函數(shù)性質(zhì)，構(gòu)建三維非對稱Copular函數(shù)結(jié)構(gòu)：

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

(8)

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

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

(9)

C(u,v,w)=

(10)

(11)

(12)

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

表1 參數(shù)辨識結(jié)果

3.2 擬合優(yōu)度檢驗

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

表2 Copula函數(shù)擬合優(yōu)度檢驗結(jié)果

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

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

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

3.3 三維參數(shù)極值Copula模型飛行風(fēng)險概率計算

根據(jù)構(gòu)建的飛行風(fēng)險判定條件，依據(jù)Joe Copula模型，求解飛行風(fēng)險概率如式(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模型。

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

4 結(jié)論

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

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