機器學習方法在氣動特性建模中的應用

何 磊,錢煒祺,汪 清,陳 海,楊 俊

(中國空氣動力研究與發展中心,綿陽 621000)

0 引 言

氣動數據建模是飛行器飛行性能仿真評估的基礎。對于飛機、導彈等飛行器，其常規飛行包線內的氣動特性通常是馬赫數、總迎角、滾轉角及三個等效舵偏角的六維函數，如果要通過計算和風洞試驗手段

飛行器的氣動特性建模主要有兩大類方法：一類是基于機理建模的方法，如多項式模型、基于幾何對稱關系的三角函數模型、微分/積分方程模型等，比較典型的例子是大迎角非定常氣動力的微分方程模型[2-7]。另一類方法是基于人工智能的“黑箱”建模方法，如決策樹方法、BP神經網絡、RBF神經網絡建模方法、模糊邏輯建模方法、Kriging方法及高斯過程法等[8-11]，其基本思想是建立決策樹、神經網絡等輸入-輸出之間的映射模型，通過訓練樣本數據的學習，對模型的結構、權值參數進行調整，使其對數據具有最佳的擬合及預測效果，因而該方法也通常稱為機器學習方法。相比較而言，基于機器學習的建模方法對非線性數據的擬合能力通常優于機理建模方法，但泛化能力較弱，而機理建模方法中的模型各項都有一定的物理意義，具有較好的普適性。近年來，深度學習方法研究取得了重大進展，該方法也是機器學習方法，目前已在流體力學和氣動特性建模中得到了初步應用[12]。為加以區分，BP神經網絡、RBF神經網絡建模方法、Kriging方法等通常稱為淺層學習方法。本文采用三類建模方法——分類與回歸樹方法、淺層學習方法和深度學習方法，對典型案例的氣動數據進行了建模計算，并對這幾類建模方法的優勢和不足進行了比較分析。

1 分類與回歸樹方法

1.1 基本思想

分類與回歸樹(Classification and Regression Tree Method，CART)方法是一種決策樹方法，也是一種非常有趣并且十分有效的非參數分類和回歸方法，是以實例為基礎，利用已知的多變量數據構建預測準則，再根據其他變量值進行預測的方法[13]。該方法采用了與傳統統計學完全不同的二叉決策樹的形式來構建預測準則，其決策樹的主要結構如圖1所示，其基本思想是根據樣本各變量值Xi與閾值Ti的比較進行分化，構建起二叉樹，樹的葉子結點Li代表的就是預測值，其中，分化的優劣及閾值Ti的選取由它產生的子結點值的離散程度來確定。決策樹的最終構建包含兩個過程：建樹和剪枝。建樹就是通過部分的訓練樣本數據建立起決策樹，直到每一個葉子結點都對應樣本點為止。剪枝的過程則是通過剩余的訓練樣本數據，對已經構建的決策樹進行結點的增減，使得決策樹在分類性能上能夠達到最優。該方法目前在航空科學研究領域已得到了應用[14-15]。

圖1 CART方法中的決策樹示例Fig.1 Decision tree in CART method

1.2 火箭氣動特性建模算例

考慮如圖2所示的軸對稱火箭，采用DATCOM軟件[16]，計算馬赫數M=1.5、總迎角α=6°情況下三個舵面都偏轉時的軸向力，作為實測數據進行建模。數據共有216組舵面組合，對應滾轉舵偏δx= 0°,5°,10°；偏航舵偏δy=-15°,-10°,-5°,0° ,5°,10°,15°,20°；俯仰舵偏δz=-20°,-15°,-10°,-5°,0° ,5°,10°,15°,20°的全組合。同時由于每組舵面組合取了8個不同的滾轉角φ=-135°,-45°,0°,22.5°,45°,67.5°,90°,180°。因而共計216×8=1728個狀態數據點。文獻[17]中基于均勻性原則從1728組數據中選出50個狀態點作為訓練樣本點，采用正交最小二乘法(OLS)進行建模，對其余1678個狀態點進行預測，得到了比較好的結果。

對這一算例考慮CART建模方法，首先分析了訓練樣本數對建模效果的影響。圖3中給出了不同訓練樣本數情況下CART模型對所有1728組樣本的預測結果，可以看到，隨著訓練樣本點的增加，預測結果的均方差逐步減小。圖4中的“CART”給出了訓練樣本數為250時建立的CART模型對一組舵偏組合狀態的預測結果，及其與正交最小二乘建模結果的比較。從圖4中可以看到該建模方法是可行的。

圖2 典型軸對稱火箭Fig.2 Typical axisymmetric rocket

圖3 訓練樣本數對預測均方差影響Fig.3 Influence of training samples number on modelling result

圖4 建模預測結果對比Fig.4 Comparison of model prediction results

2 淺層學習建模方法

2.1 RBF神經網絡建模方法

在函數逼近領域，早先最普遍采用的神經網絡是BP網絡，它是一種是多層前饋神經網絡，因網絡訓練采用誤差反向傳播(Back-Propagation，BP)學習算法而得名。近年來，隨著徑向基函數(Radial Basis Function，RBF)神經網絡的提出和發展，為函數逼近提供了一種更有效的方法，它在逼近能力、分類能力、學習速度等方面均優于BP網絡[18]。

RBF網絡是一種三層前向神經網絡：輸入層由信號源節點組成；第二層為隱層，隱單元的個數由所描述的問題而定，隱層采用徑向基函數作為網絡的激勵函數；第三層為輸出層，采用線性激勵函數。

設一個多輸入多輸出的RBF網絡有m個輸入節點、n個隱層節點、p個輸出節點。RBF網絡的輸出為：

(1)

(2)

式中cj為第j個徑向基函數的中心；σj為與第j個徑向基函數寬度相關的變量；‖·‖為L2范數，表示x與cj的歐拉距離。徑向基函數的輸出位于區間(0 ,1)內，輸入與中心向量的距離越近，則基函數的輸出越大。對于給定的輸入x，只有靠近中心部分才被激活。

2.2 Kriging建模方法

Kriging模型是從采礦和地理統計學應用中發展起來的一種響應面函數模型，是在已知近似函數的基礎上添加了局部的“擾動”來擬合樣本點[19]。該方法實際上是高斯過程建模方法(Gaussian Process，GP)的特例。高斯過程定義為一系列隨機變量的集合，且其中任意兩個隨機變量間都滿足聯合高斯分布。高斯過程建模方法核心思想是將訓練樣本視為高斯隨機過程，當系統有新輸入時，高斯過程方法對輸出的預測值是一個符合高斯分布的隨機變量。高斯過程建模需要確定均值函數、協方差函數中的未知參數，這些未知參數被統稱為超參數，通常基于貝葉斯原理來確定[20]。

Kriging建模的具體方法如下[19]：設訓練樣本的自變量x為n維向量，ns個樣本輸出構成的輸出向量為y，則其近似函數模型為：

(3)

rT(x)=[R(x,x1),R(x,x2),...,R(x,xns)]T；

(4)

這里，x的上標表示樣本點序號，下標表示矢量的分量；R(xi,xj)是任意兩個樣本點xi和xj的互相關函數；參數θk即為超參數。

(5)

(6)

而參數θk的最大似然估計則由下式給出：

(7)

(8)

式中rT(xpre)利用式(4)計算。

為提高Kriging模型的泛化能力，可通過交叉校驗的方法來對模型中的超參數進行優化。交叉校驗的基本思想是從所有樣本點中依次選取不同的樣本點進行建模后，對剩余樣本點進行預測，計算其均方預測誤差。然后基于均方預測誤差最小化的原則對超參數θk進行優化，以舍一交叉校驗為例，其均方預測誤差如下式(9)。

(9)

其中e(ns-i)表示建模數據不包含(xi,yi)時，利用其余ns-1個樣本點建模結果對xi的預測誤差。

2.3 支持向量機SVM方法

支持向量機(Support Vector Machines,SVM)是在統計學習理論的基礎上發展起來的新型算法，由V.Vapnik等在20世紀90年代提出，建立在嚴密的統計學理論基礎上，引入了結構風險最小化原則，具有較好的泛化能力。該模型的基本思想是在線性可分情況下構建最優分類超平面，獲得支持向量,然后以支持向量的核展開表示模型，從而達到模型的稀疏表示。對于非線性問題，則是將原問題映射到一個高維特征空間，然后在特征空問中進行線性運算[21-22]。最小二乘支持向量機(Least Square SVM)是標準SVM的一個變形，是將標準支持向量機算法中的不等式約束化為等式約束得到的。由此，求解二次優化問題轉化成解一次線性方程組，使收斂速度更快，計算的復雜性降低。具體算法為：

給定N個訓練樣本{xi,yi}，其中xi∈Rn為n維的訓練樣本輸入，yi∈R為訓練樣本輸出。其目標函數為：

s.t.yi(wTg(xi)+b)=1-ξi，i=1,...,N

(10)

式中，w為l維權重矢量，g(x)是將x從輸入空間映射到特征空間的函數，ξi是xi的松弛系數，b為偏置項，γ是邊際系數。

將上述目標函數增廣為拉格朗日函數形式后，對各變量求偏導并令其偏導數為0，可將求解的優化問題轉化為求解如下線性方程組，即：

(11)

其中，l=[1,1,…,1]T；y=[y1,y2,…,yN]T；α=[α1,α2,…,αN]T；I為單位矩陣；Ω為核函數矩陣，矩陣元素Ωij=f(xi,xj)，f(·)為滿足Mercer條件的對稱正定核函數，通常取為:

f(xi,xj)=exp(-‖xi-xj‖2/δ)，δ為核半徑。

由式(11)解出b和α后，LS-SVM的回歸函數為：

(12)

上述模型中有兩個超參數——正則化參數γ和核半徑δ，需要通過交叉校驗等方法來確定[23]。

2.4 算例

2.4.1 火箭氣動特性建模

沿用1.2節中的算例，將基于均勻性原則選取的50個狀態點作為訓練樣本點，采用Kriging建模方法進行建模，對其余1678個狀態點進行預測，預測的均方誤差為0.1924。進一步針對50個建模樣本采用舍一交叉校驗方法來優化超參數，對其余1678個狀態點進行預測，可得出均方誤差為0.1219。圖5中給出了兩組典型舵偏狀態下的建模預測結果對比。其中“Kriging”表示標準Kriging模型建模結果，“Kriging(with CV)”表示通過交叉校驗方法優化超參數后的建模結果。可以看到，對這一氣動建模問題能采用Kriging模型進行建模，且通過交叉校驗可提高模型的泛化能力。

圖5 Kriging模型建模預測結果Fig.5 Prediction results of Kriging model

2.4.2 三角翼大迎角非定常氣動特性建模

以NASA阿姆斯研究中心7 ft×10 ft低速風洞中開展的展弦比λ=2尖前緣三角翼俯仰振蕩動態試驗[24]為研究對象。迎角變化歷程為：

α(τ)=45°+45°sin(kτ)

(13)

首先用靜態試驗數據和減縮頻率k=0.01、0.03、0.04的大振幅振動試驗數據訓練Kriging模型、RBF神經網絡模型和SVM支持向量機模型。不同減縮頻率下的不同氣動系數訓練結果與試驗數據的比較示于圖6，圖中“Kriging”對應Kriging模型，“RBF”對應RBF神經網絡模型，“SVM”對應SVM支持向量機模型，“Exp.”對應風洞試驗數據。由圖6可見，訓練結果與試驗數據擬合得較好。

再用訓練后的模型參數預測縮減頻率k=0.02的工況。不同模型的預測結果與試驗數據比較示于圖7，具體的預測誤差值示于表1。從圖表中看出，預測結果與試驗數據也符合較好。其中，由于SVM進行了交叉校驗，一定程度上提升了模型的泛化能力，使得SVM的預測結果略優于Kriging模型和RBF神經網絡模型。

(a)k=0.01,阻力系數建模結果

(b)k=0.03,升力系數建模結果

(c)k=0.04,俯仰力矩系數建模結果

圖6 建模結果對比
Fig.6 Comparison of modeling results

(a)k=0.02,阻力系數預測結果

(b)k=0.02,升力系數預測結果

(c)k=0.02,俯仰力矩系數預測結果

圖7 預測結果對比
Fig.7 Comparison of model prediction results

2.4.3 氣動熱試驗數據融合

以某飛行器的表面熱流數據為例，對數值計算獲得的連續分布熱環境和風洞試驗點熱流值進行融合。要求融合后的結果在風洞試驗點位置處同測量結果一致，對CFD數值仿真網格點處的值進行適當的修正，融合后的結果能夠在整體上體現試驗結果的數據特性。圖8是所有氣動熱數據樣本點的位置，其中32 922個藍色點為CFD網格點，計算數據如圖9所示，118紅色點位置為風洞試驗測量點，試驗數據如圖10示。

表1 k=0.02工況的不同模型預測誤差Table 1 Comparison of prediction error for different model for k=0.02 case

圖8 數據采集點Fig.8 Data sampling points

圖9 數值計算數據Fig.9 Data of CFD

圖10 風洞試驗數據Fig.10 Wind tunnel test data

數據融合過程中需要建立風洞試驗數據和CFD計算數據的差量模型，并使用差量模型對CFD網格點處的值進行修正。由于Kriging模型具有準確擬合建模數據點的特性，因此被用于建立差量模型。另外，由于風洞試驗點位置與CFD網格點并非完全對應，因此對每個風洞試驗點，選取附近的CFD數據點利用RBF方法進行插值。之所以采用RBF方法是由于實驗表明在數據較少時，RBF建模效率比Kriging更高，且建模精度相當。

由此可得氣動熱數據融合方法：(1)針對每一個風洞試驗點，若有與之對應的CFD數據點，則直接使用，否則，選取附近少量CFD數據，使用RBF方法建立插值模型，并計算相應的值;(2)計算風洞試驗數據和CFD計算數據的差量，并使用Kriging方法建立差量模型;(3)使用Kriging模型對所有的CFD數據點進行修正，得到融合結果。

為了驗證本文所提方法的有效性，隨機選取10%的風洞試驗數據作為預測驗證數據，其余90%的數據作為樣本數據與CFD計算數據一起構建融合模型。實驗進行了兩次，圖11(a)和圖11(b)的紅色點是兩次實驗隨機選取的10%風洞試驗數據在空間的分布情況。圖12(a)和圖12(b)是驗證結果對比，其中橫坐標表示用于驗證的風洞試驗數據的真實值，縱坐標表示在這些驗證位置通過融合模型計算得到的預測值。第一次和第二次實驗驗證的均方根誤差分別為0.04782和0.02762。由圖12和均方根誤差可知預測數據與驗證數據基本一致，說明模型達到了數據融合的目的。

(a)第一次實驗

(b)第二次實驗圖11 兩次實驗驗證數據分布圖Fig.11 The positions of validation data for two experimental cases

(a)第一次實驗

(b)第二次實驗圖12 預測數據與驗證數據對比Fig.12 Comparison of validation data and prediction

2.5 建模方法對比分析

通過上述機器學習模型算例分析，可以看到，幾種方法都能對氣動數據進行較好的建模。但幾種模型也有各自的特點和優缺點，如表2所示。

(1)Kriging建模方法可對訓練樣本完全擬合，RBF神經網絡方法和SVM支持向量機方法則不要求訓練樣本完全擬合。

(2)Kriging建模方法、RBF神經網絡方法和SVM支持向量機方法都可以通過交叉校驗優化超參數來提高模型的泛化能力。對Kriging模型而言，超參數是θk；對RBF模型而言，超參數是隱層節點數；對SVM模型而言，超參數是γ和δ。CART方法則對訓練樣本量有較高要求，建模結果的精度和模型的泛化能力隨訓練樣本數增加而提高。

(3)Kriging方法、RBF神經網絡方法和SVM支持向量機方法都需要對樣本數據進行大量的矩陣運算，但RBF神經網絡方法的矩陣維數取決于隱層節點數，而Kriging方法和SVM支持向量機方法的矩陣維數取決于樣本點數，當樣本點數較多時，矩陣運算的計算資源要求會增加，計算效率和計算精度會下降。因此，Kriging方法和SVM支持向量機方法中的建模樣本點數不宜過多，在普通個人微機上一般不超過2000個樣本點，具體數量視微機性能而定，而RBF神經網絡方法則對樣本點數無特殊限制。CART建模方法中沒有矩陣運算，且樣本點越多，建模效果越好，但隨著樣本點的增加，邏輯判斷運算量較大。

表2 幾種建模方法比較Table 2 Comparison of modeling methods

3 氣動特性建模的深度學習方法

前面的淺層學習方法存在兩方面問題，一是需要人為提取輸入特征參數，二是泛化能力有所不足。深度神經網絡通過增加隱層層數，一方面使得模型能更有效地表征非線性函數的豐富特征；另一方面則通過層與層之間的無監督學習，實現特征的逐級提取，可較好地克服淺層學習方法存在的不足。但是，隨著網絡層數的增加，如果采用傳統的神經網絡訓練方法會出現三個問題，一是“梯度消失”現象；二是優化易于陷入局部最優值；三是如果網絡全連接的話，會帶來網絡權值參數數目的幾何增長。因此，在深度神經網絡訓練方法中，使用自下而上非監督學習與自上而下監督學習結合的方法，對訓練數據特征提取后再進行網絡權值調整，解決了“梯度消失”和陷入局部最優的問題；同時基于局部感知和權值共享，則大幅減少了待優化網絡權值數[25-26]。

近年來，隨著計算資源的增強和大數據的興起，以深度神經網絡為代表的深度學習方法在圖像與語音識別、語義理解、智能控制等方面得到了廣泛應用，并逐步與空氣動力學研究結合[27-33]。文獻[33]將深度神經網絡模型用于氣動特性建模，基于CNN模型對翼型氣動系數進行了預測，但該方法只能處理給定馬赫數和迎角的情況，沒有將流動條件作為輸入，限制了其使用范圍。為了在輸入中引入流動條件，文獻[34]中利用迎角度數對翼型圖像進行旋轉，利用馬赫數對圖像顏色進行改變，構造了一幅人工圖像，作為預測模型的輸入。本節提出另一種方法，將流動條件參數組成一個向量，認為該向量是一幅分辨率為1×nfc的圖像，其中nfc為流動條件的參數個數。通過雙三次插值，將該向量的分辨率調整為同翼型圖像一致，與翼型圖像一起構成“合成圖像”，即2個二維矩陣，作為預測模型的輸入。同時，由于翼型圖像與流動參數向量對應圖像的特征豐富程度不同，選用如圖13所示的有向無環圖(Directed Acyclic Graph，DAG)深度神經網絡[35]，該網絡結構在卷積層1后分為兩個支路，能較好地處理圖像特征相差較大的情況，提高預測精度[11]。其訓練和測試方法與常規CNN預測模型類似。具體的網絡參數設置：圖像分辨率為85×85；卷積層1選取卷積核大小為8×8，卷積核數量為40，步長為1；卷積層2選取卷積核大小為5×5，卷積核數量為80，步長為1；卷積層3選取卷積核大小為5×5，卷積核數量為120，步長為2；卷積層4選取卷積核大小為2×2，卷積核數量為120，步長為2；所有卷積層所用到的非線性激勵函數均選取為ReLU函數；池化區域大小為2×2，步長為2，即相鄰池化區域不重疊，選擇的池化方法為最大池化法。

圖13 DAG模型結構Fig.13 Model structure of DAG network

算例中用到的翼型為NACA0012翼型疊加擾動函數后的變形體，共300個；流動條件中迎角取值為2°、6°、9°、12°和15°，馬赫數取值為0.1、0.3和0.6，雷諾數固定為6.5×106；翼型形狀和流動條件任意組合，共構成翼型數據樣本4200組。圖14中給出了隨機選取的8個翼型圖像。圖15給出了3個迎角和馬赫數組成的向量[6,0.3]、[9,0.6]、[12,0.1]乘以放大倍數后變為[102,127.5]、[153,255]、[204,42.5]，視為三張圖片將分辨率調整為85×85后的圖像。

圖14 任意8個翼型圖像(翼型厚度放大10倍)Fig.14 Eight random selected airfoils (thickness amplified 10 times)

圖15 任意3個由迎角和馬赫數構成的圖像Fig.15 Pictures corresponding to three sets of angle of attack and Mach number

將這4200組樣本的順序打亂，隨機選擇3360組翼型圖像和對應的流動條件作為訓練樣本，將氣動系數作為期望值，對預測模型進行訓練；剩下的840組翼型數據作為測試樣本，對預測模型進行測試。結果顯示，840組翼型升力系數、阻力系數、俯仰力矩系數預測結果的標準差絕對值分別為0.0418、0.0049和0.0042，相對值分別為4.38%、7.23%和12.73%。隨機給出一個翼型在不同迎角和馬赫數下CFD計算氣動系數和預測氣動系數的對比如圖16所示。圖16中，“CFD”表示通過CFD計算得到的氣動系數，“DAG”表示通過模型預測得到的氣動系數，α表示迎角。從圖16中可以看出，本文建立的預測模型能夠比較準確地對3個氣動系數進行預測。

(a)升力系數對比

(b)阻力系數對比

在此需要指出的是，深度學習模型雖然從建模便捷性和泛化能力等方面看優于淺層學習模型，但也存在不足，還需要進一步的改進完善：一是對訓練數據量的要求較高；二是模型中超參數，如網絡層數、隱層節點數的選取，較為困難，這些超參數的選取對建模結果有直接影響；三是所提取特征的可解釋性較弱。

4 小 結

本文對氣動特性“黑箱”建模的三類機器學習方法——分類與回歸樹方法、淺層學習方法和深度學習方法,進行了算法說明與分析應用。將分類與回歸樹方法、Kriging建模方法、RBF神經網絡方法及SVM支持向量機方法分別應用于火箭氣動特性建模、三角翼大迎角非定常氣動特性建模、氣動熱試驗數據融合。結果顯示淺層學習方法可通過交叉校驗對模型中的超參數進行優化來提高泛化能力；Kriging建模方法和SVM支持向量機方法的訓練樣本數受矩陣運算限制較大。CART方法對訓練樣本量有較高要求，但其建模結果的精度和模型的泛化能力隨訓練樣本數增加而提高。

同時，將流動條件參數組成向量，再映射為圖像，與翼型圖像構成“合成圖像”，建立了以翼型幾何圖像、來流馬赫數、迎角為輸入，翼型氣動特性為輸出的DAG深度神經網絡模型，得到了比較好的預測效果，拓展了氣動特性深度學習建模方法的使用范圍。

從上述對淺層學習和深度學習模型的分析應用看，淺層學習在模型結構簡潔性、計算效率和建模數據量方面有優勢，而深度學習模型在泛化能力、特征自動提取方面有優勢。因此，后續可考慮將兩種建模方法結合起來開展研究。例如，提高深度學習模型的可解釋性，分析其提取特征的物理含義并將這些特征作為淺層學習的輸入，使淺層學習模型在較少數據量的情況下也能具有較強的泛化能力。