何 磊,錢煒祺,易 賢,王 強,張顯才
(1. 中國空氣動力研究與發展中心 計算空氣動力研究所, 四川 綿陽 621000;2. 中國空氣動力研究與發展中心 低速空氣動力研究所, 四川 綿陽 621000)
結冰問題一直以來是影響航空飛行器飛行安全的重要隱患之一[1]。飛機結冰后,機翼、尾翼和舵面上的積冰破壞了物面附近流場,嚴重影響飛機氣動性能和飛行性能,增加飛行風險,危害飛行安全。我國幅員遼闊,地勢和氣象條件復雜多變,先后發生了多起因結冰導致的墜機事故[2]。由于飛機結冰與氣象條件密切相關,若能根據氣象條件實時預測可能發生的結冰外形情況,不僅可以為飛行安全包線確定提供依據,為飛行員提供結冰信息和安全預警信息,也能支撐飛機的防除冰系統設計和智能結冰系統研制。因此,國內外進行了大量飛機結冰研究。
長期以來,開展結冰研究主要依靠結冰風洞試驗[3]、飛行試驗[4]和CFD模擬[5-6]。歐美一些發達國家紛紛建造結冰風洞并開發結冰數值模擬軟件[7-8],獲取了大量飛機結冰數據[9]。這些數據被很多結冰研究作為驗證或建模樣本。國內結冰實驗設備和數值模擬研究方面起步較晚,但中國空氣動力研究與發展中心、北京航空航天大學、西北工業大學、南京航空航天大學等單位都在積極開展相關研究[10]。
隨著機器學習算法的發展,國內外許多專家學者將機器學習方法應用到飛機結冰研究領域,包括自組織特性映射,徑向基函數、BP神經網絡、概率神經網絡等,用于解決結冰嚴重性探測[11]、結冰位置和體積探測[12-13]、結冰冰形預測[14]、結冰后飛機氣動特性影響[15-16]等問題。在結冰冰形預測方面,人工神經網絡類算法應用較為廣泛和成功,取得了許多有價值的成果。
結冰冰形預測研究內容主要包括對冰形曲線進行描述和建立預測模型兩個方面。Ogretim等[17]最早提出一種基于神經網絡的翼型結冰預測方法:首先需要進行坐標變換得到新的冰形曲線,將坐標變換后的冰形曲線看作一個非周期的復雜信號,展開為傅里葉級數形式,將冰形與飛行狀態和氣象參數之間映射關系的建模轉化為了傅里葉系數與飛行狀態和氣象參數之間的建模;然后利用神經網絡進行建模和預測,模型輸入為來流速度、溫度、液態水含量、平均水滴等效直徑和結冰時間5個參數,輸出為冰形曲線傅里葉級數的正余弦系數;該方法主要針對NACA00XX 系列翼型,在坐標變換算法方面有一定局限性。潘環等[18]也用類似方法對冰形預測的建模與方法進行了研究,并增加了相對濕度和攻角2項參數。Chang等[14]和李珺[19]改進了Ogretim等的方法,分別提出用小波包變換方法和多值變量擬合函數替代傅里葉變換對冰形曲線的描述。2種方法的冰形預測效果不僅取決于神經網絡模型的設計和訓練效果,也取決于對冰形的描述能力。另外,李珺僅僅采用了NASA格林結冰研究中心的80組風洞結冰數據,樣本偏少,限制了神經網絡的訓練效果。鮑雨晨等[20]基于BP神經網絡研究了通用飛機機翼的冰形預測方法,該方法對翼型特征點上的結冰厚度進行預測,再通過計算公式轉換為冰形特征點坐標,對坐標點連線可得冰形輪廓;該研究為冰形預測提供了不錯的實踐途徑,但研究中僅標記了14個機翼特征點,對冰形輪廓的描述能力稍顯不足;另外研究中僅包含43組訓練樣本,而輸入的實驗變量包含空速、姿態角、總溫等10個參數,訓練樣本數量明顯偏少,訓練獲得的模型可能無法應對實驗變量值變化稍大的情況。
近年來,深度學習作為機器學習的重要分支發展迅速,并繼承了神經網絡對非線性關系的描述能力,在眾多領域得到了廣泛應用。本文的主要研究工作就是基于深度學習技術構建翼型結冰冰形的預測模型,實現對冰形的圖像化預測能力。
對翼型結冰冰形進行智能預測需要建立預測模型,其核心就是描述影響飛機結冰因素和冰形之間的映射關系,如圖1所示。

圖1 翼型結冰預測問題Fig.1 Problem of airfoil ice accretion prediction
影響飛機結冰冰形的參數主要有2類:大氣環境參數和飛行狀態參數。對于飛行狀態參數,主要考慮飛行攻角α、飛行速度v;而大氣環境參數,主要考慮液態水含量(Liquid Water Content, LWC)、水滴平均直徑(Median Volumetric Diameter, MVD) 、環境溫度T、結冰時長t。
以往采用的描述映射關系的方法,如Kriging、BP神經網絡、RBF神經網絡等,輸出冰形曲線都是數值類型的,需要通過翼面特征點法向的冰形值去描述冰形曲線,或通過類似傅里葉變換的方法對冰形曲線做進一步處理,這些方法都要求翼面同一位置處法線方向冰形厚度只能存在單值,從而無法解決復雜冰形在翼面同一位置處法線方向冰形厚度存在多值的問題,如圖2所示[15]。

圖2 復雜冰形多值問題Fig.2 Multivalued problem of complex ice shape
為克服這一問題,提出采用圖像方式對冰形進行描述,圖像方式比傳統數值方式更直觀,描述能力更強,因此建模的關鍵就在于設計直接將圖像作為輸出的模型架構。
根據上述對結冰預測問題的分析,提出如圖3所示的建模和預測框架。主要包括數據獲取、數據預處理、模型設計、模型訓練、模型預測等環節,框架的特點是將圖像化的冰形作為模型訓練標簽和預測的輸出。

圖3 建模與預測框架Fig.3 Framework of modeling and prediction
翼型結冰冰形預測模型的輸入輸出數據格式是模型結構和接口設計的關鍵。對于輸入參數,即影響飛機結冰的飛行狀態和大氣環境因素,其值存在符號和數量級的差異,因此需要通過歸一化處理將其數值范圍調整到一致的范圍之內。對于輸出的冰形圖像,考慮將其灰度化,既能滿足對冰形幾何特征描述的需求,又可避免多通道帶來的額外計算量。冰形圖像尺寸設置為512×256,每個像素點的取值在0~255之間,實際運算中將其歸一化到0~1之間。
如圖4所示,為避免由于圖像中冰形的坐標軸取值范圍不一致對模型帶來的計算誤差,冰形灰度圖像寬度對應的x坐標軸取值范圍固定為[-0.2, 1.0],圖像高度對應的y坐標軸取值范圍固定為[-0.2, 0.2],坐標軸取值范圍可根據實際情況調整。

圖4 結冰翼型圖像坐標軸范圍固定Fig.4 Fixed the coordinate range of iced airfoil image
卷積神經網絡(Convolutional Neural Network,CNN)是深度神經網絡的典型結構,卷積操作為卷積神經網絡提供了強大的圖像特征提取能力。轉置卷積(又名反卷積)是卷積操作的相反過程,可以對卷積操作提取的編碼器中的特征進行解碼。利用多個轉置卷積操作可以實現生成圖片的目的,深度卷積對抗生成網絡的生成器就是使用多個轉置卷積層對隨機噪聲值進行操作以實現生成完整的圖片功能。
因此基于轉置卷積神經網絡設計了如圖5所示冰形預測模型結構。輸入的結冰條件參數通過2個全連接網絡層映射到更大尺寸的數據結構,2個全連接網絡層均使用ReLU激活函數。接著通過5個轉置卷積網絡層將低維特征向量向高維特征空間映射。除了預測模型的輸出層外,其他的轉置卷積網絡層后都連接了批標準化(Batch Normalization, BN)網絡層[21]和ReLU激活函數層,批標準化技術通過規范化手段,使得每一層神經網絡的輸入保持相同分布,以保證梯度傳播到每一層,避免出現梯度消失現象。從結冰條件參數到冰形圖像是一個典型的非線性映射問題,ReLU激活函數層的功能就是為神經網絡提供非線性映射功能,從而提升冰形預測模型的預測能力,并使用Dropout網絡層隨機丟棄神經網絡單元,提高泛化能力;輸出層在轉置卷積層后連接一個sigmoid激活函數層。

圖5 翼型冰形預測模型網絡結構Fig.5 Network structure of airfoil ice shape prediction model
模型采用二元交叉熵(binary cross entropy)損失函數,為了防止模型出現過擬合,進一步提升泛化性能,損失函數中加入L2正則化懲罰項,損失函數如式(1)所示。式中,右邊第一項為交叉熵損失項,右邊第二項為L2正則化項,它表示了模型的復雜度。
(1)
式中:yi表示類別i的真實標簽;pi表示模型計算出類別i的概率值;N表示訓練樣本總數;k表示類別數,二元交叉熵中k=2;w表示神經網絡的權重;λ表示L2正則化率;n表示整個網絡的神經元總數。
考慮以某運輸機尾翼翼型(NACA0012)為實驗對象,通過數值模擬方法生成實驗所需冰形數據。
對于飛行參數取值范圍,考慮固定飛行攻角條件下飛行速度的取值,所研究的運輸機最小飛行速度約為268 km/h(74.44 m/s)、巡航飛行速度為550 km/h(152.77 m/s)。
對于氣象參數取值范圍,有研究表明:飛機最容易發生結冰的溫度范圍為0~-20 ℃,特別是在-2~-10 ℃ 范圍內遭遇結冰的次數通常最多,而-2~-8 ℃是飛機發生強烈結冰的主要溫度范圍[22]。結冰中常見的過冷水滴平均直徑在20~40 μm之間[19],如《中國民用航空規章第 25 部運輸類飛機適航標準》(CCAR 25部)附錄C中通常考慮的過冷水滴尺寸就在15~40 μm范圍。液態水含量是影響結冰最重要的因素之一,影響云中液態水含量的因素較多,隨水汽凝結和降水而變化,云層中不同位置液態水含量也不相同。根據CCAR 25部附錄C中給出的大氣約束條件,通常液態水含量的考慮范圍為0.2~0.8 g/m3。
綜上分析并結合具體研究需求,確定的結冰參數取值范圍如表1所示。

表1 結冰參數取值范圍
本節對典型翼型結冰進行數值模擬,并將冰形計算結果與風洞試驗結果進行對比,以驗證數值模擬方法的可靠性。風洞試驗數據取自美國NASA格林結冰研究中心和 FAA 威廉·J.休斯技術中心聯合發起的現代翼型項目的翼型結冰資料[15],該資料記錄了3種典型飛機翼型在多種結冰條件下的冰形和氣動特性風洞試驗數據。
3.2.1 數值計算
翼型采用格林結冰研究中心的現代翼型項目中記錄的商用噴氣飛機(business jet)翼型[24-25],如圖6所示,該翼型是現代商用飛機的典型翼型。

圖6 商用飛機翼型Fig.6 Airfoil of business jet
共安排2組實驗,結冰條件參數包括液態水含量、水滴平均直徑、環境溫度、結冰時長、飛行速度、飛行攻角。2組實驗結冰條件參數取值如表2所示,分別對應結冰資料中在結冰研究風洞的試驗序號“run213”和“run214”[15],除結冰時長不同外,其他條件一致。

表2 驗證算例計算條件
翼型結冰數值計算采用中國空氣動力研究與發展中心的計算方法和計算軟件irc2d,主要包括流場計算、過冷水滴運動計算、結冰計算和物面外形更新4個步驟[16, 26]。如圖7所示,計算采用C型網格,在網格生成階段已經加入攻角。

圖7 結冰計算網格(α=6.2°)Fig.7 Computational grid of ice accretion(α=6.2°)
3.2.2 結果驗證
表3是風洞試驗的冰形幾何特征參數值,圖8分別給出了2組實驗數值計算冰形和風洞試驗冰形的對比情況。對于冰體輪廓,兩組實驗的數值計算結果與風洞試驗結果基本一致,冰角厚度、角度也符合較好;但對于冰體在翼型物面的結冰上極限和下極限位置,計算結果和試驗結果對比稍有差異,兩組實驗中計算冰形上極限稍大,下極限稍小;實驗1的駐點結冰厚度符合較好,實驗2的駐點結冰厚度略有差異。總體而言,計算冰形和試驗冰形在冰體輪廓、冰形體積和主要特征方面符合較好,說明所采用的翼型結冰數值模擬方法是可靠的。

表3 冰形幾何特征參數值

(a) 實驗1(a) Experiment 1

(b) 實驗2(b) Experiment 2圖8 計算結果與風洞試驗結果對比Fig.8 Comparison of computational and wind tunnel test results
使用上節所述irc2d翼型結冰數值仿真軟件計算NACA0012翼型在α=2°的情況下的結冰外形。計算采用C型網格,如圖9所示。

圖9 結冰計算網格(α=2°)Fig.9 Computational grid of ice accretion(α=2°)
計算共獲得11 200組訓練樣本數據,對應飛行速度為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;液態水含量為0.2 g/m3, 0.35 g/m3, 0.5 g/m3, 0.65 g/m3, 0.8 g/m3;結冰時長為6 min, 9 min, 12 min, 15 min, 18 min, 20 min, 22.5 min的全組合。
共獲得768組驗證樣本數據,對應飛行速度為85 m/s, 98 m/s, 120 m/s, 145 m/s;溫度為-18 ℃, -10 ℃, -5 ℃, -3 ℃;水滴直徑為23 μm, 32 μm, 52 μm;液態水含量為0.3 g/m3, 0.4 g/m3, 0.6 g/m3, 0.7g/m3;結冰時長為7 min, 10 min, 12.5 min, 18.5 min的全組合。
獲得冰形樣本數據后,按照建模框架所述數據規范將訓練樣本集轉化為如圖10所示的灰度圖像。

圖10 翼型結冰冰形灰度圖像Fig.10 Gray scaled images of iced airfoil shape
根據圖5所示神經網絡模型建立冰形預測模型,并使用上節生成的11 200組訓練樣本對網絡模型進行訓練,768組驗證樣本對訓練好的網絡模型進行測試。由于實驗數據的飛行攻角固定值為2°,可以只考慮其他5個影響因素,因此網絡輸入層尺寸為5。
論文基于開源深度學習架構tensorflow的keras高級接口實現所提的翼型結冰預測模型的構建、訓練和預測。預測模型訓練參數設置如下:訓練優化方法選擇Adam[27],其中學習率(learning rate)的初始值為0.001,參數β1設置為0.9,β2設置為0.999,為穩定訓練過程使用學習率衰減,設置值為1.0×10-8;分批大小為20,即每批輸入20組結冰條件參數數據和冰形圖像標簽進行訓練;迭代次數為200。
模型采用CPU訓練模式,用于建模和訓練的計算機配置為:Intel Core i7-7700、3.6 GHz、4核8線程CPU、16 GB內存。基于以上訓練參數,對11 200個樣本的訓練耗時約為88.2 h。訓練過程中,樣本的損失函數值隨迭代次數的變化曲線如圖11所示。從圖11中可以看出,模型在約前10次迭代中,損失函數值迅速趨于收斂,之后隨著迭代次數的增加,有繼續小幅下降的趨勢,最終穩定在0.007 8附近。

圖11 訓練過程Fig.11 Training history
模型訓練完成后,將訓練樣本集和驗證樣本集的結冰條件參數輸入預測模型,預測其對應的翼型結冰冰形圖像,測試模型的預測能力,檢驗預測模型的訓練和泛化效果。訓練樣本集預測耗時576 s(單個樣本耗時約0.048 s),驗證樣本集預測耗時38 s(單個樣本預測耗時約0.049 s)。
模型預測輸出結果為類似圖10的灰度圖像,為便于觀察和對比分析,利用圖像技術對結果進行了加入原始翼型物面、顏色替換、剪裁等后處理操作。
圖12展示了3組典型訓練冰形的模型預測結果,以及與CFD計算冰形進行對比的情況。從圖12中可見,在冰體輪廓、結冰上下極限位置、冰角厚度和角度等主要幾何特征方面,3組實驗預測結果,以及與CFD計算結果都基本重合,符合較好,說明模型在訓練樣本集上表現良好。

(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 典型訓練冰形預測結果Fig.12 Prediction of typical training ice shape
圖13展示了4組典型測試冰形的模型預測結果,以及與CFD計算冰形進行對比的情況。對于冰體輪廓,4組實驗的模型預測結果與CFD計算結果基本重合,符合較好。冰形冰角角度、結冰上下極限位置等幾何特征參數也都符合較好。對于結冰厚度這一特征,部分實驗模型預測結果與CFD計算結果略有差異,且這種差異隨著冰形厚度的增加而變大。但總體而言,模型預測冰形和CFD計算冰形在冰體輪廓、冰形體積和主要特征方面符合較好。可見,預測模型的整體泛化性能較好,但對結冰較厚情況的泛化性能仍有進一步提升空間。

(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 典型測試冰形預測結果Fig.13 Prediction of typical testing ice shape
本文針對飛機結冰問題,開展了翼型結冰冰形預測方法研究,提出了建模和預測框架,設計了基于深度神經網絡的圖像化冰形預測模型,用于建立冰形與飛行狀態參數、氣象參數之間的映射關系,主要考慮了飛行速度、攻角、大氣中液態水含量、水滴平均直徑、結冰時的溫度、結冰時長等多物理參數對冰形的影響。
以某運輸機水平尾翼(NACA0012翼型)為對象,利用CFD數值模擬生成的冰形作為訓練和驗證樣本,對所建立的冰形預測模型進行了訓練和測試。結果表明:
1)提出的翼型結冰冰形圖像化預測方法是可行的,預測冰形與CFD數值計算的冰形在冰形輪廓、結冰上下極限、上下冰角位置、結冰厚度等主要幾何特征參數方面都符合較好,但在結冰較厚情況下,模型泛化性能還可進一步提高。
2)雖然模型訓練耗時較多,但預測模型訓練完成后,便能夠快速預測一定范圍內的冰形,且計算速度快(單工況計算耗時約50 ms),能適應機載要求;另外在模型訓練過程中,若使用圖形處理器計算能大幅減少訓練時間。
總之,深度學習方法在飛機結冰研究領域具有很強的應用前景,下一步考慮將風洞試驗數據加入樣本集對預測模型進行訓練,增強模型的工程實用性。同時,在現有研究基礎上,繼續基于深度學習方法研究飛機結冰后對氣動特性的影響。