何 磊,錢煒祺,易 賢,王 強,張顯才
(1. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000;2. 中國空氣動力研究與發(fā)展中心 低速空氣動力研究所, 四川 綿陽 621000)
結(jié)冰問題一直以來是影響航空飛行器飛行安全的重要隱患之一[1]。飛機結(jié)冰后,機翼、尾翼和舵面上的積冰破壞了物面附近流場,嚴重影響飛機氣動性能和飛行性能,增加飛行風險,危害飛行安全。我國幅員遼闊,地勢和氣象條件復雜多變,先后發(fā)生了多起因結(jié)冰導致的墜機事故[2]。由于飛機結(jié)冰與氣象條件密切相關(guān),若能根據(jù)氣象條件實時預測可能發(fā)生的結(jié)冰外形情況,不僅可以為飛行安全包線確定提供依據(jù),為飛行員提供結(jié)冰信息和安全預警信息,也能支撐飛機的防除冰系統(tǒng)設計和智能結(jié)冰系統(tǒng)研制。因此,國內(nèi)外進行了大量飛機結(jié)冰研究。
長期以來,開展結(jié)冰研究主要依靠結(jié)冰風洞試驗[3]、飛行試驗[4]和CFD模擬[5-6]。歐美一些發(fā)達國家紛紛建造結(jié)冰風洞并開發(fā)結(jié)冰數(shù)值模擬軟件[7-8],獲取了大量飛機結(jié)冰數(shù)據(jù)[9]。這些數(shù)據(jù)被很多結(jié)冰研究作為驗證或建模樣本。國內(nèi)結(jié)冰實驗設備和數(shù)值模擬研究方面起步較晚,但中國空氣動力研究與發(fā)展中心、北京航空航天大學、西北工業(yè)大學、南京航空航天大學等單位都在積極開展相關(guān)研究[10]。
隨著機器學習算法的發(fā)展,國內(nèi)外許多專家學者將機器學習方法應用到飛機結(jié)冰研究領(lǐng)域,包括自組織特性映射,徑向基函數(shù)、BP神經(jīng)網(wǎng)絡、概率神經(jīng)網(wǎng)絡等,用于解決結(jié)冰嚴重性探測[11]、結(jié)冰位置和體積探測[12-13]、結(jié)冰冰形預測[14]、結(jié)冰后飛機氣動特性影響[15-16]等問題。在結(jié)冰冰形預測方面,人工神經(jīng)網(wǎng)絡類算法應用較為廣泛和成功,取得了許多有價值的成果。
結(jié)冰冰形預測研究內(nèi)容主要包括對冰形曲線進行描述和建立預測模型兩個方面。Ogretim等[17]最早提出一種基于神經(jīng)網(wǎng)絡的翼型結(jié)冰預測方法:首先需要進行坐標變換得到新的冰形曲線,將坐標變換后的冰形曲線看作一個非周期的復雜信號,展開為傅里葉級數(shù)形式,將冰形與飛行狀態(tài)和氣象參數(shù)之間映射關(guān)系的建模轉(zhuǎn)化為了傅里葉系數(shù)與飛行狀態(tài)和氣象參數(shù)之間的建模;然后利用神經(jīng)網(wǎng)絡進行建模和預測,模型輸入為來流速度、溫度、液態(tài)水含量、平均水滴等效直徑和結(jié)冰時間5個參數(shù),輸出為冰形曲線傅里葉級數(shù)的正余弦系數(shù);該方法主要針對NACA00XX 系列翼型,在坐標變換算法方面有一定局限性。潘環(huán)等[18]也用類似方法對冰形預測的建模與方法進行了研究,并增加了相對濕度和攻角2項參數(shù)。Chang等[14]和李珺[19]改進了Ogretim等的方法,分別提出用小波包變換方法和多值變量擬合函數(shù)替代傅里葉變換對冰形曲線的描述。2種方法的冰形預測效果不僅取決于神經(jīng)網(wǎng)絡模型的設計和訓練效果,也取決于對冰形的描述能力。另外,李珺僅僅采用了NASA格林結(jié)冰研究中心的80組風洞結(jié)冰數(shù)據(jù),樣本偏少,限制了神經(jīng)網(wǎng)絡的訓練效果。鮑雨晨等[20]基于BP神經(jīng)網(wǎng)絡研究了通用飛機機翼的冰形預測方法,該方法對翼型特征點上的結(jié)冰厚度進行預測,再通過計算公式轉(zhuǎn)換為冰形特征點坐標,對坐標點連線可得冰形輪廓;該研究為冰形預測提供了不錯的實踐途徑,但研究中僅標記了14個機翼特征點,對冰形輪廓的描述能力稍顯不足;另外研究中僅包含43組訓練樣本,而輸入的實驗變量包含空速、姿態(tài)角、總溫等10個參數(shù),訓練樣本數(shù)量明顯偏少,訓練獲得的模型可能無法應對實驗變量值變化稍大的情況。
近年來,深度學習作為機器學習的重要分支發(fā)展迅速,并繼承了神經(jīng)網(wǎng)絡對非線性關(guān)系的描述能力,在眾多領(lǐng)域得到了廣泛應用。本文的主要研究工作就是基于深度學習技術(shù)構(gòu)建翼型結(jié)冰冰形的預測模型,實現(xiàn)對冰形的圖像化預測能力。
對翼型結(jié)冰冰形進行智能預測需要建立預測模型,其核心就是描述影響飛機結(jié)冰因素和冰形之間的映射關(guān)系,如圖1所示。
圖1 翼型結(jié)冰預測問題Fig.1 Problem of airfoil ice accretion prediction
影響飛機結(jié)冰冰形的參數(shù)主要有2類:大氣環(huán)境參數(shù)和飛行狀態(tài)參數(shù)。對于飛行狀態(tài)參數(shù),主要考慮飛行攻角α、飛行速度v;而大氣環(huán)境參數(shù),主要考慮液態(tài)水含量(Liquid Water Content, LWC)、水滴平均直徑(Median Volumetric Diameter, MVD) 、環(huán)境溫度T、結(jié)冰時長t。
以往采用的描述映射關(guān)系的方法,如Kriging、BP神經(jīng)網(wǎng)絡、RBF神經(jīng)網(wǎng)絡等,輸出冰形曲線都是數(shù)值類型的,需要通過翼面特征點法向的冰形值去描述冰形曲線,或通過類似傅里葉變換的方法對冰形曲線做進一步處理,這些方法都要求翼面同一位置處法線方向冰形厚度只能存在單值,從而無法解決復雜冰形在翼面同一位置處法線方向冰形厚度存在多值的問題,如圖2所示[15]。
圖2 復雜冰形多值問題Fig.2 Multivalued problem of complex ice shape
為克服這一問題,提出采用圖像方式對冰形進行描述,圖像方式比傳統(tǒng)數(shù)值方式更直觀,描述能力更強,因此建模的關(guān)鍵就在于設計直接將圖像作為輸出的模型架構(gòu)。
根據(jù)上述對結(jié)冰預測問題的分析,提出如圖3所示的建模和預測框架。主要包括數(shù)據(jù)獲取、數(shù)據(jù)預處理、模型設計、模型訓練、模型預測等環(huán)節(jié),框架的特點是將圖像化的冰形作為模型訓練標簽和預測的輸出。
圖3 建模與預測框架Fig.3 Framework of modeling and prediction
翼型結(jié)冰冰形預測模型的輸入輸出數(shù)據(jù)格式是模型結(jié)構(gòu)和接口設計的關(guān)鍵。對于輸入?yún)?shù),即影響飛機結(jié)冰的飛行狀態(tài)和大氣環(huán)境因素,其值存在符號和數(shù)量級的差異,因此需要通過歸一化處理將其數(shù)值范圍調(diào)整到一致的范圍之內(nèi)。對于輸出的冰形圖像,考慮將其灰度化,既能滿足對冰形幾何特征描述的需求,又可避免多通道帶來的額外計算量。冰形圖像尺寸設置為512×256,每個像素點的取值在0~255之間,實際運算中將其歸一化到0~1之間。
如圖4所示,為避免由于圖像中冰形的坐標軸取值范圍不一致對模型帶來的計算誤差,冰形灰度圖像寬度對應的x坐標軸取值范圍固定為[-0.2, 1.0],圖像高度對應的y坐標軸取值范圍固定為[-0.2, 0.2],坐標軸取值范圍可根據(jù)實際情況調(diào)整。
圖4 結(jié)冰翼型圖像坐標軸范圍固定Fig.4 Fixed the coordinate range of iced airfoil image
卷積神經(jīng)網(wǎng)絡(Convolutional Neural Network,CNN)是深度神經(jīng)網(wǎng)絡的典型結(jié)構(gòu),卷積操作為卷積神經(jīng)網(wǎng)絡提供了強大的圖像特征提取能力。轉(zhuǎn)置卷積(又名反卷積)是卷積操作的相反過程,可以對卷積操作提取的編碼器中的特征進行解碼。利用多個轉(zhuǎn)置卷積操作可以實現(xiàn)生成圖片的目的,深度卷積對抗生成網(wǎng)絡的生成器就是使用多個轉(zhuǎn)置卷積層對隨機噪聲值進行操作以實現(xiàn)生成完整的圖片功能。
因此基于轉(zhuǎn)置卷積神經(jīng)網(wǎng)絡設計了如圖5所示冰形預測模型結(jié)構(gòu)。輸入的結(jié)冰條件參數(shù)通過2個全連接網(wǎng)絡層映射到更大尺寸的數(shù)據(jù)結(jié)構(gòu),2個全連接網(wǎng)絡層均使用ReLU激活函數(shù)。接著通過5個轉(zhuǎn)置卷積網(wǎng)絡層將低維特征向量向高維特征空間映射。除了預測模型的輸出層外,其他的轉(zhuǎn)置卷積網(wǎng)絡層后都連接了批標準化(Batch Normalization, BN)網(wǎng)絡層[21]和ReLU激活函數(shù)層,批標準化技術(shù)通過規(guī)范化手段,使得每一層神經(jīng)網(wǎng)絡的輸入保持相同分布,以保證梯度傳播到每一層,避免出現(xiàn)梯度消失現(xiàn)象。從結(jié)冰條件參數(shù)到冰形圖像是一個典型的非線性映射問題,ReLU激活函數(shù)層的功能就是為神經(jīng)網(wǎng)絡提供非線性映射功能,從而提升冰形預測模型的預測能力,并使用Dropout網(wǎng)絡層隨機丟棄神經(jīng)網(wǎng)絡單元,提高泛化能力;輸出層在轉(zhuǎn)置卷積層后連接一個sigmoid激活函數(shù)層。
圖5 翼型冰形預測模型網(wǎng)絡結(jié)構(gòu)Fig.5 Network structure of airfoil ice shape prediction model
模型采用二元交叉熵(binary cross entropy)損失函數(shù),為了防止模型出現(xiàn)過擬合,進一步提升泛化性能,損失函數(shù)中加入L2正則化懲罰項,損失函數(shù)如式(1)所示。式中,右邊第一項為交叉熵損失項,右邊第二項為L2正則化項,它表示了模型的復雜度。
(1)
式中:yi表示類別i的真實標簽;pi表示模型計算出類別i的概率值;N表示訓練樣本總數(shù);k表示類別數(shù),二元交叉熵中k=2;w表示神經(jīng)網(wǎng)絡的權(quán)重;λ表示L2正則化率;n表示整個網(wǎng)絡的神經(jīng)元總數(shù)。
考慮以某運輸機尾翼翼型(NACA0012)為實驗對象,通過數(shù)值模擬方法生成實驗所需冰形數(shù)據(jù)。
對于飛行參數(shù)取值范圍,考慮固定飛行攻角條件下飛行速度的取值,所研究的運輸機最小飛行速度約為268 km/h(74.44 m/s)、巡航飛行速度為550 km/h(152.77 m/s)。
對于氣象參數(shù)取值范圍,有研究表明:飛機最容易發(fā)生結(jié)冰的溫度范圍為0~-20 ℃,特別是在-2~-10 ℃ 范圍內(nèi)遭遇結(jié)冰的次數(shù)通常最多,而-2~-8 ℃是飛機發(fā)生強烈結(jié)冰的主要溫度范圍[22]。結(jié)冰中常見的過冷水滴平均直徑在20~40 μm之間[19],如《中國民用航空規(guī)章第 25 部運輸類飛機適航標準》(CCAR 25部)附錄C中通常考慮的過冷水滴尺寸就在15~40 μm范圍。液態(tài)水含量是影響結(jié)冰最重要的因素之一,影響云中液態(tài)水含量的因素較多,隨水汽凝結(jié)和降水而變化,云層中不同位置液態(tài)水含量也不相同。根據(jù)CCAR 25部附錄C中給出的大氣約束條件,通常液態(tài)水含量的考慮范圍為0.2~0.8 g/m3。
綜上分析并結(jié)合具體研究需求,確定的結(jié)冰參數(shù)取值范圍如表1所示。
表1 結(jié)冰參數(shù)取值范圍
本節(jié)對典型翼型結(jié)冰進行數(shù)值模擬,并將冰形計算結(jié)果與風洞試驗結(jié)果進行對比,以驗證數(shù)值模擬方法的可靠性。風洞試驗數(shù)據(jù)取自美國NASA格林結(jié)冰研究中心和 FAA 威廉·J.休斯技術(shù)中心聯(lián)合發(fā)起的現(xiàn)代翼型項目的翼型結(jié)冰資料[15],該資料記錄了3種典型飛機翼型在多種結(jié)冰條件下的冰形和氣動特性風洞試驗數(shù)據(jù)。
3.2.1 數(shù)值計算
翼型采用格林結(jié)冰研究中心的現(xiàn)代翼型項目中記錄的商用噴氣飛機(business jet)翼型[24-25],如圖6所示,該翼型是現(xiàn)代商用飛機的典型翼型。
圖6 商用飛機翼型Fig.6 Airfoil of business jet
共安排2組實驗,結(jié)冰條件參數(shù)包括液態(tài)水含量、水滴平均直徑、環(huán)境溫度、結(jié)冰時長、飛行速度、飛行攻角。2組實驗結(jié)冰條件參數(shù)取值如表2所示,分別對應結(jié)冰資料中在結(jié)冰研究風洞的試驗序號“run213”和“run214”[15],除結(jié)冰時長不同外,其他條件一致。
表2 驗證算例計算條件
翼型結(jié)冰數(shù)值計算采用中國空氣動力研究與發(fā)展中心的計算方法和計算軟件irc2d,主要包括流場計算、過冷水滴運動計算、結(jié)冰計算和物面外形更新4個步驟[16, 26]。如圖7所示,計算采用C型網(wǎng)格,在網(wǎng)格生成階段已經(jīng)加入攻角。
圖7 結(jié)冰計算網(wǎng)格(α=6.2°)Fig.7 Computational grid of ice accretion(α=6.2°)
3.2.2 結(jié)果驗證
表3是風洞試驗的冰形幾何特征參數(shù)值,圖8分別給出了2組實驗數(shù)值計算冰形和風洞試驗冰形的對比情況。對于冰體輪廓,兩組實驗的數(shù)值計算結(jié)果與風洞試驗結(jié)果基本一致,冰角厚度、角度也符合較好;但對于冰體在翼型物面的結(jié)冰上極限和下極限位置,計算結(jié)果和試驗結(jié)果對比稍有差異,兩組實驗中計算冰形上極限稍大,下極限稍小;實驗1的駐點結(jié)冰厚度符合較好,實驗2的駐點結(jié)冰厚度略有差異??傮w而言,計算冰形和試驗冰形在冰體輪廓、冰形體積和主要特征方面符合較好,說明所采用的翼型結(jié)冰數(shù)值模擬方法是可靠的。
表3 冰形幾何特征參數(shù)值
(a) 實驗1(a) Experiment 1
(b) 實驗2(b) Experiment 2圖8 計算結(jié)果與風洞試驗結(jié)果對比Fig.8 Comparison of computational and wind tunnel test results
使用上節(jié)所述irc2d翼型結(jié)冰數(shù)值仿真軟件計算NACA0012翼型在α=2°的情況下的結(jié)冰外形。計算采用C型網(wǎng)格,如圖9所示。
圖9 結(jié)冰計算網(wǎng)格(α=2°)Fig.9 Computational grid of ice accretion(α=2°)
計算共獲得11 200組訓練樣本數(shù)據(jù),對應飛行速度為70 m/s, 80 m/s, 90 m/s, 100 m/s, 110 m/s, 125 m/s, 140 m/s, 150 m/s;溫度為-20 ℃, -17 ℃, -14 ℃, -11 ℃, -8 ℃, -6 ℃, -4 ℃, -2 ℃;水滴直徑為20 μm, 30 μm, 40 μm, 50 μm, 60 μm;液態(tài)水含量為0.2 g/m3, 0.35 g/m3, 0.5 g/m3, 0.65 g/m3, 0.8 g/m3;結(jié)冰時長為6 min, 9 min, 12 min, 15 min, 18 min, 20 min, 22.5 min的全組合。
共獲得768組驗證樣本數(shù)據(jù),對應飛行速度為85 m/s, 98 m/s, 120 m/s, 145 m/s;溫度為-18 ℃, -10 ℃, -5 ℃, -3 ℃;水滴直徑為23 μm, 32 μm, 52 μm;液態(tài)水含量為0.3 g/m3, 0.4 g/m3, 0.6 g/m3, 0.7g/m3;結(jié)冰時長為7 min, 10 min, 12.5 min, 18.5 min的全組合。
獲得冰形樣本數(shù)據(jù)后,按照建??蚣芩鰯?shù)據(jù)規(guī)范將訓練樣本集轉(zhuǎn)化為如圖10所示的灰度圖像。
圖10 翼型結(jié)冰冰形灰度圖像Fig.10 Gray scaled images of iced airfoil shape
根據(jù)圖5所示神經(jīng)網(wǎng)絡模型建立冰形預測模型,并使用上節(jié)生成的11 200組訓練樣本對網(wǎng)絡模型進行訓練,768組驗證樣本對訓練好的網(wǎng)絡模型進行測試。由于實驗數(shù)據(jù)的飛行攻角固定值為2°,可以只考慮其他5個影響因素,因此網(wǎng)絡輸入層尺寸為5。
論文基于開源深度學習架構(gòu)tensorflow的keras高級接口實現(xiàn)所提的翼型結(jié)冰預測模型的構(gòu)建、訓練和預測。預測模型訓練參數(shù)設置如下:訓練優(yōu)化方法選擇Adam[27],其中學習率(learning rate)的初始值為0.001,參數(shù)β1設置為0.9,β2設置為0.999,為穩(wěn)定訓練過程使用學習率衰減,設置值為1.0×10-8;分批大小為20,即每批輸入20組結(jié)冰條件參數(shù)數(shù)據(jù)和冰形圖像標簽進行訓練;迭代次數(shù)為200。
模型采用CPU訓練模式,用于建模和訓練的計算機配置為:Intel Core i7-7700、3.6 GHz、4核8線程CPU、16 GB內(nèi)存?;谝陨嫌柧殔?shù),對11 200個樣本的訓練耗時約為88.2 h。訓練過程中,樣本的損失函數(shù)值隨迭代次數(shù)的變化曲線如圖11所示。從圖11中可以看出,模型在約前10次迭代中,損失函數(shù)值迅速趨于收斂,之后隨著迭代次數(shù)的增加,有繼續(xù)小幅下降的趨勢,最終穩(wěn)定在0.007 8附近。
圖11 訓練過程Fig.11 Training history
模型訓練完成后,將訓練樣本集和驗證樣本集的結(jié)冰條件參數(shù)輸入預測模型,預測其對應的翼型結(jié)冰冰形圖像,測試模型的預測能力,檢驗預測模型的訓練和泛化效果。訓練樣本集預測耗時576 s(單個樣本耗時約0.048 s),驗證樣本集預測耗時38 s(單個樣本預測耗時約0.049 s)。
模型預測輸出結(jié)果為類似圖10的灰度圖像,為便于觀察和對比分析,利用圖像技術(shù)對結(jié)果進行了加入原始翼型物面、顏色替換、剪裁等后處理操作。
圖12展示了3組典型訓練冰形的模型預測結(jié)果,以及與CFD計算冰形進行對比的情況。從圖12中可見,在冰體輪廓、結(jié)冰上下極限位置、冰角厚度和角度等主要幾何特征方面,3組實驗預測結(jié)果,以及與CFD計算結(jié)果都基本重合,符合較好,說明模型在訓練樣本集上表現(xiàn)良好。
(a) α=2°, v=70 m/s, T=-20 ℃, MVD=30 μm,LWC=0.5 g/m3, t=20 min
(b) α=2°, v=80 m/s, T=-20 ℃, MVD=20 μm,LWC=0.8 g/m3, t=18 min
(c) α=2°, v=140 m/s, T=-17 ℃, MVD=60 μm,LWC=0.8 g/m3, t=15 min圖12 典型訓練冰形預測結(jié)果Fig.12 Prediction of typical training ice shape
圖13展示了4組典型測試冰形的模型預測結(jié)果,以及與CFD計算冰形進行對比的情況。對于冰體輪廓,4組實驗的模型預測結(jié)果與CFD計算結(jié)果基本重合,符合較好。冰形冰角角度、結(jié)冰上下極限位置等幾何特征參數(shù)也都符合較好。對于結(jié)冰厚度這一特征,部分實驗模型預測結(jié)果與CFD計算結(jié)果略有差異,且這種差異隨著冰形厚度的增加而變大。但總體而言,模型預測冰形和CFD計算冰形在冰體輪廓、冰形體積和主要特征方面符合較好??梢?,預測模型的整體泛化性能較好,但對結(jié)冰較厚情況的泛化性能仍有進一步提升空間。
(a) α=2°, v=85 m/s, T=-10 ℃, MVD=35 μm,LWC=0.7 g/m3, t=12.5 min
(b) α=2°, v=98 m/s, T=-18 ℃, MVD=35 μm,LWC=0.4 g/m3, t=18.5 min
(c) α=2°, v=145 m/s, T=-18 ℃, MVD=52 μm,LWC=0.6 g/m3, t=10 min
(d) α=2°, v=145 m/s, T=-18 ℃, MVD=23 μm,LWC=0.7 g/m3, t=12.5 min 圖13 典型測試冰形預測結(jié)果Fig.13 Prediction of typical testing ice shape
本文針對飛機結(jié)冰問題,開展了翼型結(jié)冰冰形預測方法研究,提出了建模和預測框架,設計了基于深度神經(jīng)網(wǎng)絡的圖像化冰形預測模型,用于建立冰形與飛行狀態(tài)參數(shù)、氣象參數(shù)之間的映射關(guān)系,主要考慮了飛行速度、攻角、大氣中液態(tài)水含量、水滴平均直徑、結(jié)冰時的溫度、結(jié)冰時長等多物理參數(shù)對冰形的影響。
以某運輸機水平尾翼(NACA0012翼型)為對象,利用CFD數(shù)值模擬生成的冰形作為訓練和驗證樣本,對所建立的冰形預測模型進行了訓練和測試。結(jié)果表明:
1)提出的翼型結(jié)冰冰形圖像化預測方法是可行的,預測冰形與CFD數(shù)值計算的冰形在冰形輪廓、結(jié)冰上下極限、上下冰角位置、結(jié)冰厚度等主要幾何特征參數(shù)方面都符合較好,但在結(jié)冰較厚情況下,模型泛化性能還可進一步提高。
2)雖然模型訓練耗時較多,但預測模型訓練完成后,便能夠快速預測一定范圍內(nèi)的冰形,且計算速度快(單工況計算耗時約50 ms),能適應機載要求;另外在模型訓練過程中,若使用圖形處理器計算能大幅減少訓練時間。
總之,深度學習方法在飛機結(jié)冰研究領(lǐng)域具有很強的應用前景,下一步考慮將風洞試驗數(shù)據(jù)加入樣本集對預測模型進行訓練,增強模型的工程實用性。同時,在現(xiàn)有研究基礎上,繼續(xù)基于深度學習方法研究飛機結(jié)冰后對氣動特性的影響。