基于深度神經網絡的高階非線性激波判別式

楊 華，陳樹生，2，*，李 猛，龐 超，高正紅，向星皓

(1.西北工業大學 航空學院，西安 710072；2.中航工業第一飛機設計研究院，西安 710089；3.空氣動力學國家重點實驗室，綿陽 621000)

0 引言

激波是廣泛存在于可壓縮流動中的重要現象之一，對飛行器的氣動特性影響巨大。計算流體力學（Computational Fluid Dynamics，CFD）作為可壓縮流動問題研究的主要手段之一[1]，近幾十年來，其數值方法主要是圍繞如何準確捕捉激波以及如何能夠在光滑區域滿足高精度的要求中不斷發展的[2-3]。

對存在激波、接觸間斷等現象的流場進行數值計算時，為了保證數值格式的非線性、穩定性，同時保持較高分辨率，可以采用一種高分辨率混合格式[4]，該混合方法將非線性激波捕捉格式與針對光滑流動區域的高分辨率線性格式結合，能準確且高效地識別激波位置，在包含間斷的流場中清晰地捕捉激波，最終得到穩定的計算結果。為了提高混合格式的計算效率和求解精度，有必要發展更加精準的激波識別方法。傳統的激波識別方法均為人工構造的判別式且多含有經驗參數[5]，因此在實際應用中存在一定局限性。近年來，機器學習在偏微分方程（partial differential equation,PDE）、反問題等領域的應用日漸增多。借助機器學習方法構建新型間斷識別器可以有效避免對人工經驗參數的依賴，具有巨大潛力。

Ray 等[6]構建了基于多層感知機（multilayer perceptron,MLP）的間斷識別器并將其作為黑箱，在采用Runge-Kutta 間斷 Galerkin（Runge-Kutta discontinuous Galerkin,RKDG）格式框架下，將其與minmod 限制器和基于minmod 的TVB 限制器進行比較，數值結果表明MLP 識別器具有更好的魯棒性，并且不依賴于問題相關的參數。隨后Ray 等[7]在此基礎上，將該方法擴展到二維非結構網格中檢測間斷位置，并縮減了人工神經網絡的結構，使其具有更高的計算效率。Abgrall 等[8]利用RKDG 方法和高階限制器的模擬數據訓練神經網絡，并將其作為激波識別器；數值結果表明其在笛卡爾網格以及非結構網格上中均有較高的精度。Sun 等[9]構建了基于卷積神經網絡的通用間斷檢測器，在間斷處采用高階WENO 格式，光滑處采用高階中心差分格式，求解一維和二維問題的數值結果表明該檢測器能夠更加清晰地識別出間斷位置。

在前人的基礎上，本文基于深度前饋神經網絡模型構建了高階非線性激波判別式，與Huynh[10]提出的二階無振蕩格式（centered non-oscillatory scheme of second order,Cn2）中的激波判別方法進行了對比，并將其應用于二維超聲速無黏圓柱繞流算例中。

1 訓練及測試數據集的生成

本文構建的基于神經網絡的高階非線性激波判別式流程如圖1 所示，分別是生成訓練集和驗證集、訓練模型以及驗證。

圖1 基于神經網絡的激波判別式建模流程Fig.1 Process of shock discriminant modeling based on Neural Network

生成訓練數據集與測試數據集的步驟如下：

Step1：選擇一個函數u(x)，在區間 [a,b]內包含其規律性的完整信息，如選擇正弦函數、階躍函數等。

Step2：選取任意點xi∈[a,b]及 網格大小h，確保a≤xi-5/2·h＜xi＜xi+5/2·h≤b。以xi為中心，在區間 [a,b]內創建5 個網格單元，如圖2 所示。

圖2 在區間[a,b]內隨機抽取網格點示意圖Fig.2 Random extraction of grid points within an interval [a,b]

Step4：對應于輸入向量X的實際輸出取決于選取的函數性質及規律。理想情況下，如果網格單元Ii處的解不連續，則該網格單元會被標記為間斷點，輸出為Y=(1,0)T，否則輸出為Y=(0,1)T。

Step5：重復2～4 步，每次重復需改變點xi的位置，網格大小h。

Step6：重復1～5 步，每次重復使用不同的函數u(x)。

需要注意的是，為了使得所生成樣本更符合物理現象，在生成訓練集與測試集中的間斷函數數據時，默認參數a、b保持相同的正負號。

根據上述所列步驟，分別選取三角函數、一次函數、二次函數和間斷函數等生成訓練數據集T 與測試數據集V，如表1 所示。

表1 用于模型訓練和測試的數據集Table 1 Data sets for model training and testing

2 神經網絡模型

本節采用深度前饋神經網絡（feedforward neural network,FNN）[12]構建激波識別器并開展神經網絡模型的訓練及驗證，最后選取預測精度最高的模型作為最終實際使用的模型。

2.1 深度前饋神經網絡模型

深度前饋神經網絡也叫做前饋神經網絡或者多層感知機，是典型的深度學習模型。前饋網絡的目標是近似某個函數f*，通過定義一個映射y=f(x,θ)，并且學習參數 θ的值，使它能夠得到最佳的函數近似。

在前饋神經網絡中，各神經元分別屬于不同的層，每一層的神經元可以接收前一層神經元的信號，并產生信號輸出到下一層。第0 層稱為輸入層，最后一層為輸出層，其他中間層稱為隱藏層。整個網絡中沒有反饋，信號從輸入層向輸出層單向傳播。其中網絡的層數和各隱藏層的單元個數均為超參數，通常超參數是按照人為經驗設定，并不斷試錯調整最終確定[12]。

本文所構建的深度前饋神經網絡模型，選擇7 層神經網絡，輸入層為7 個神經元，輸出層為2 個神經元，中間有5 個隱藏層，其神經元個數分別為256、128、64、32、16。此網絡結構中的參數還包 括W={W1,W2,W3,W4,W5,W6} 和b={b1,b2,b3,b4,b5,b6}。其中W1∈R256×7、W2∈R128×256、W3∈R64×128、W4∈R32×64、W5∈R16×32、W6∈R2×16，b1∈R256、b2∈R128、b3∈R64、b4∈R32、b5∈R16、b6∈R2。激活函數選 擇LeakyReLU 函數。網絡結構如圖3 所示。

圖3 FNN 模型網絡結構示意圖Fig.3 FNN model network structure diagram

2.2 訓練神經網絡

本文選擇交叉熵損失作為模型訓練的損失函數：

式（1）中，S是一次訓練所使用的樣本 數；為所給樣本第k個網格單元是否為間斷點的實際標簽，其值只有0 和1 之分；分別為模型預測該樣本第k個網格單元為間斷點或光滑點的概率，其中。通過交叉熵損失可以計算出兩個概率分布與實際值之間的誤差，最小化訓練樣本中預測值與真實值之間的誤差，從而達到提高模型預測精度的目的。

模型訓練的優化算法選擇Adam 算法，此外為了防止模型過擬合，提高泛化能力，同時應用正則化與提前終止方法。訓練的殘差變化如圖4 所示。

圖4 模型訓練殘差圖Fig.4 Model training residuals

在實際訓練中，為了避免標準梯度下降中存在的無法充分利用計算機并行計算能力的問題，使用小批量梯度下降法減少參數更新的波動并節省計算資源，每一小批量選擇200 個訓練樣本以實現效果更好更穩定的收斂。其中經過每一小批訓練樣本后的測試精度如圖5 所示，選取其中測試精度最高的模型參數作為此次訓練的結果。

圖5 模型訓練中的測試精度變化圖Fig.5 Test accuracy variation in model training

模型訓練中選擇重復訓練參數R=10，即按照上述流程重復訓練10 次模型，選取其中測試數據集預測精度最高的參數作為最終模型參數。

3 模型驗證

本節采用多個算例（Lax 激波問題、Sod 激波管問題及二維超聲速無黏圓柱繞流算例）對FNN 模型進行驗證，檢驗所訓練模型的預測精度及其適用性。

3.1 一維算例

本小節首先將深度前饋神經網絡模型（FNN）和原Cn2 格式[10,13]中的激波識別方法（此處簡稱為Cn2）進行對比，通過上述生成的訓練和測試數據集檢驗兩種激波判別式的識別效果。接著將FNN 模型嵌入Cn2 格式中以識別激波位置，并形成新的數值格式（此處簡稱為FNN 格式），將該格式與原Cn2 格式進行對比。通過求解一維的Lax 激波問題以及Sod 激波管問題[14]以驗證FNN 模型在數值格式中的有效性。

原Cn2 格式中的激波識別方法[10]如下：

以密度場為例，令：

定義該網格單元Ij為光滑區域，否則標記為間斷點。

3.1.1 訓練和測試數據集

本文生成的訓練數據集包含20 000 個光滑樣本和20 000 個間斷樣本，測試數據集包含12 000 個光滑樣本和12 000 個間斷樣本。應用FNN 和Cn2 激波識別方法對訓練數據集和測試數據集中的光滑樣本和間斷樣本分別進行預測，并對比兩種激波識別方法的預測精度，如表2 所示。

表2 FNN 和Cn2 激波識別方法的預測精度對比Table 2 Comparison of predictive accuracy between FNN and Cn2 shock wave identification methods

一般神經網絡模型在訓練過程中，各類樣本的權重是相等的，從而使得模型在對各類樣本進行預測時保持預測精度基本一致。但針對激波識別問題，為保證數值求解的穩定性，通常要求在所有可能出現大梯度或者間斷的流場處引入數值耗散來抑制數值非物理振蕩；同時為保證流場模擬的精度，通常希望在光滑區域采用高分辨率低耗散格式。因此在FNN 模型訓練過程中，經驗性地設置光滑樣本和間斷樣本的權重分別為[0.04,1.0]，在保證神經網絡模型將所有間斷樣本都能識別出來的情況下，盡可能提高其對光滑樣本的預測精度。

從表1 可以看到，通過在模型訓練過程中給定各個類別樣本不同的權重，放寬對光滑樣本的預測精度，從而使得模型在預測間斷樣本時的預測精度更高。相較于Cn2 激波識別方法，FNN 模型提高了在光滑樣本處的預測精度，并保持了正確識別間斷樣本的特性。

3.1.2 Lax 激波問題

計算域為[0,1]，計算網格數為100，初始條件為：

計算終止時間t=1.3。

圖6 中給出了FNN 與Cn2 格式求解Lax 激波問題時的無量綱密度在計算域[0,1]的分布曲線，其中紅色圓點為1 表示識別的激波位置。可以看到，FNN 格式分別將兩個間斷附近的7 個網格點和6 個網格點識別為激波區域，Cn2 格式分別將兩個間斷附近的10 個網格點和7 個網格點識別為激波區域。相較而言，FNN 格式識別出的激波位置更加精準，Cn2 格式在光滑區域會產生更多誤判，這和3.1.1小節從訓練及測試數據集中得到的結論是一致的。

圖6 Lax 問題密度分布曲線Fig.6 Density distributions of Lax’s problem

3.1.3 Sod 激波管問題

計算域為[0,1]，計算網格數為100，初始條件為：

計算終止時間為t=0.2。

圖7 中給出了FNN 與Cn2 格式求解Sod 激波管問題時的無量綱密度在計算域[0,1]的分布曲線。可以看到，FNN 格式分別將兩個間斷附近的4 個網格點識別為激波區域，Cn2 格式分別將兩個間斷附近的8 個網格點和6 個網格點識別為激波區域。兩種激波識別方法均能很好地識別激波位置，但FNN 格式相較而言更加精準。

圖7 Sod 問題密度分布曲線Fig.7 Density distributions of Sod’s problem

3.1.4 一維算例結果討論

表3 為FNN 和Cn2 格式在求解上述Lax 激波問題和Sod 激波管問題時的耗時情況，選擇CPU_TIME函數估計兩種格式在同一計算平臺上的計算時間。

表3 運用FNN 和Cn2 格式的計算時間對比Table 3 Comparison of computing time between FNN and Cn2 schemes

從表3 的對比可以看到，相對于Cn2 格式，FNN格式在求解Lax 激波問題時，消耗時間增加了22%；在求解Sod 激波管問題時，消耗時間增加了25%。由于FNN 模型在每個網格點處判別其是否為間斷時，均需要重新計算498 個神經元的值才能給出判別結果，而Cn2 格式中激波識別方法為一次代數表達式，相較而言其消耗時間更少。

對于本文提到的Lax 和Sod 問題，FNN 模型所增加的時間消耗較少，完全可以忽略不計。對于二維或三維復雜流場數值模擬，由于流場網格數量激增以及可能出現激波反射、激波干擾等現象，采用FNN 方法相較于Cn2 方法會產生更大的時間代價；當時間代價達到一定量級時，可以使用多步判斷以節省數值模擬時間。此外，可以考慮修改神經網絡模型結構，通過適當減少網絡層數及神經元個數，以減小程序運行的時間消耗。

3.2 二維超聲速無黏圓柱繞流算例

本小節采用超聲速無黏圓柱繞流[15-16]（不同來流馬赫數、不同網格分布）來驗證深度前饋神經網絡模型（FNN）在二維真實流場中判別激波間斷的有效性。數值計算方法采用Roe 格式加剪切黏性[17]，剪切黏性激活函數采用本文的FNN 激波識別方法。

圖8 為來流馬赫數Ma分別為5、10、15、20 時二維圓柱流場示意圖。其中上半部為壓力云圖，下半部為使用FNN 模型識別的流場激波結果。由圖中可以看出，訓練得到的FNN 模型對于不同來流馬赫數的流場都有很好的激波識別能力。其中網格采用161 ×161（周向×壁面法向）均勻分布網格。

圖8 不同馬赫數下超聲速二維圓柱繞流流場Fig.8 Flow field around a supersonic two-dimensional cylinder at different Mach numbers

圖9 為不同網格分布下超聲速二維圓柱計算所得的流場圖（Ma=20）。其中網格所示的第一個參數表示周向網格數，第二個參數表示壁面法向網格數，從左至右、從上至下網格長寬比分別為8、1/8、4、1/4、4、1/4、2、1/2、1。可以看到不同網格分布下，FNN 模型都能準確地識別出激波位置。其中圖9（i）為161 ×161（周向×壁面法向）的均勻分布網格，其識別效果最好。

圖9 不同網格分布下超聲速二維圓柱繞流流場（Ma=20）Fig.9 Flow field around a supersonic two-dimensional circular cylinder with different grid distribution (Ma=20)

圖9（b、d、f、h）分別對應長寬比為1/8、1/4、1/4、1/2 的網格計算結果。使用FNN 模型在中段正激波處識別激波位置較為精確，在末段斜激波處識別范圍較大。這是由于壁面法向網格布點較多而周向網格布點較少，導致在壁面法向垂直方向的激波分辨率較高而在周向垂直方向激波分辨率較差，在末段斜激波附近即體現為所識別激波范圍較大，因此激波識別精度依賴于網格分辨率。

4 結論

本文針對超聲速流動中的激波識別問題，采用深度前饋神經網絡構建了高階非線性激波判別式，得到以下結論：

1）本文構建的高階非線性激波判別式FNN 相較于二階振蕩格式Cn2 中人工構造的激波識別器，能夠更加精準地識別間斷區域和光滑區域，但一定程度上會增加時間成本。

2）本文以激波判別式FNN 驗證了其對不同來流馬赫數、不同網格分布下二維超聲速無黏圓柱繞流算例中激波位置的預測能力。算例結果表明，在不同來流馬赫數下FNN 模型均能很好地捕捉激波位置；對于不同網格分布，FNN 模型都能識別出激波位置，但激波識別精度依賴于網格分辨率，其中均勻分布網格表現最好。

綜上，本文激波判別式在不同來流馬赫數及不同網格分布下均有良好的適用性，但在將其應用于工程三維算例計算時，存在計算耗時多的問題，因此需在保持較高激波識別精度基礎上，繼續縮減模型結構，以便該激波判別式的后續應用，進一步為精確捕捉激波以及構造高精度數值格式奠定基礎。