基于稀疏多項式混沌方法的不確定性量化分析

陳江濤，章超，劉驍，趙輝，胡星志，吳曉軍

中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000

在工程復雜外形的CFD模擬過程中，存在著多種來源的不確定性輸入，包括計算網格、湍流模型和數值格式等[1-2]，這使得最終的模擬結果也有不可忽視的不確定性。在工程外形的性能評估和優化設計、CFD與試驗數據的相互確認等過程中，需要考慮CFD結果的不確定性，這樣才能更為準確地評價設計外形能否滿足性能需求，優化過程是否有效可靠。因此CFD結果的不確定性量化研究受到越來越廣泛和持續的關注。

蒙特卡洛(Monte Carlo，MC)方法[3]構造簡單、易于實現，是一種被廣泛使用的不確定性量化方法。但是MC方法收斂速度慢，計算開銷巨大，需要大量的采樣點計算才能得到精度滿足要求的結果，因此限制了其在多變量不確定性量化問題中的應用。多項式混沌(Polynomial Chaos，PC)方法[4-6]是一種更加高效的不確定性量化方法。Wiener[7]在研究湍流問題中隨機變量的譜空間展開時提出了各向同性混沌理論，這也被認為是PC方法的起源。Ghanem和Spanos[8]在分析固體力學中的不確定性問題時將該方法應用到有限元計算中。Xiu和Karniadakis[9]將適合高斯隨機過程的Hermite PC推廣到Wiener-Askey PC，適用于滿足更一般概率密度分布的隨機過程。基于PC方法的不確定性量化已經有了很多成功的應用[10-16]。

PC展開的項數隨著輸入變量的維數和展開階數的增加呈幾何級數式急劇增加，導致“維數災難”現象產生，這限制了PC在多變量不確定性量化問題中的應用。近些年來，各國學者們嘗試了更加高效的方法，包括自適應方法[17-19]、稀疏重構方法[20-21]和減基法[22-23]等，以此來降低隨機空間的維數。計算表明這些方法與完整的PC方法相比確實能夠在相對較小的計算成本下產生精度大體相當的結果。本文主要研究了稀疏多項式混沌方法。在包含多輸入變量的隨機問題中，經常會遇到目標響應的展開是稀疏的情況[24-25]，展開中占主導的通常是獨立的輸入變量和變量之間的低階交叉項，某些自由度的量值可能非常小，完全可以舍去，這也是稀疏多項式混沌方法能夠成功應用的前提。Blatman和Sudret[20]使用稀疏效應準則，提出了一種雙曲算法來截斷PC展開。后來信號處理領域的稀疏感知算法被用來重構稀疏的PC展開[26-32]。該算法的效率取決于PC展開的稀疏性或者說是展開式中自由度的衰減程度。

本文發展了非嵌入式的稀疏多項式混沌方法，提出了一種截斷誤差設定方法，并在湍流模型系數不確定性傳播和炭化材料燒蝕熱響應模擬不確定性分析中開展了應用研究。文章首先介紹了基于1范數最小的稀疏重構算法，然后對兩個算例進行了分析研究，最后給出了研究結論。

1 稀疏PC展開

PC方法將整個數值模擬過程看作是一個隨機過程，根據隨機輸入參數的概率密度分布，選擇對應的正交基函數序列，然后將該隨機過程的輸出在基函數空間中展開。通過嵌入式或者非嵌入式方法，確定PC展開式中的自由度，從而得到輸出的各種統計信息，包括平均值、方差以及每個輸入變量的全局敏感性等。

假定隨機的輸入變量ξ，滿足某種概率密度分布f(ξ)。對于任意的隨機輸出變量y，可以在輸入變量ξ的正交基函數序列構成的譜空間中展開，即

(1)

式中：αj為待求的展開式自由度；ψj為確定的隨機變量的基函數。

在回歸法中，通過確定性的CFD計算得到采樣點{ξ1,ξ2,…,ξN}對應的響應值為{y1,y2,…,yN}。通常情況下，采樣點的數目大于PC展開中自由度的個數，因此形成了一個關于自由度的超定方程組：Ψα=Y，其中Ψ是測量矩陣，Ψij=ψj(ξi)，α是自由度組成的向量，Y是響應值組成的向量。PC方法更詳細的介紹可以參考文獻[6]。

壓縮感知(Compressed Sensing，CS)方法是圖像和信號處理領域興起的新方法，能夠高效地重構稀疏信號，需要的采樣點數目小于自由度的個數。CS通過選擇有限個對系統輸出影響較大的基函數，實現由不完全的觀測樣本準確或者近似準確地還原稀疏信號。如果PC展開是稀疏的，即自由度非零的個數是有限的，可以通過求解P0優化問題得到完整的展開，即

(2)

式中：0范數表示α中非零元素的個數。然而P0優化問題中的目標函數是非凸函數，這是一個NP-hard問題，在實際中很難求解[33]。為此將P0優化問題轉化為凸松弛的P1優化問題，即

(3)

結果表明基于1范數最小的算法能夠有效地處理擁有有限個非零自由度的PC展開問題。

在實際應用中，考慮到測量噪聲的情況下，不要求約束Ψα=Y精確滿足。因此得到了有噪聲信號的P1優化問題：

(4)

在求解過程中，需要指定截斷誤差ε。ε的設定對于稀疏PC重構的精度十分關鍵。如果ε設置得過大，則重構的PC不夠準確；如果ε設置得過小，則重構的PC有可能出現過擬合現象。ε的設定可以認為是模型評估和選擇問題。本文通過機器學習中常用的交叉驗證方法得到較優的截斷誤差設定。具體做法為：首先給定有限個初選的截斷誤差，分別求解有噪聲信號的P1優化問題，得到每個截斷誤差對應的留一法(Leave-One-Out，LOO)誤差。最終選擇最小的LOO誤差對應的截斷誤差為優化問題的設定值。LOO誤差是在機器學習領域廣泛使用的誤差估計。將N個樣本分為N-1個學習樣本和1個驗證樣本，通過學習樣本建立PC展開，在驗證樣本上評估該展開的泛化誤差。在PC法中，LOO誤差表達式為

(5)

大體上有兩類算法來解決上述優化問題:線性規劃算法和貪婪算法。兩種方法的優劣本文不作深入研究,這里使用貪婪算法中的正交匹配追蹤算法[33]。該算法通過計算殘差向量和每個基函數向量的相關性來得到與當前殘差最相關的基函數，然后通過求解最小二乘問題得到響應的展開。

算法的輸入是測量矩陣Ψ、觀測向量Y、截斷誤差ε。定義迭代次數指標k，殘差向量r，特征索引集，具體實施過程為

初始化：k=0,r(0)=Y,=?

1) 定義殘差向量與每個基函數向量的內積，即

并找到與當前殘差最相關的特征：λ=argmaxjλj。

2) 更新迭代次數指標和特征索引集：

k←k+1

r(k)=Y-Ψα(k)

2 不確定性量化分析

本文選擇兩個算例來驗證稀疏多項式混沌方法和其中的截斷誤差設定策略，分別研究了湍流模型系數的不確定性對RAE2822翼型跨聲速繞流模擬的影響和材料物性參數的不確定性對燒蝕熱響應預測的影響。

2.1 SA模型系數的不確定性對跨聲速翼型模擬的影響

第1個算例是Spalart-Allmaras(SA)模型系數的不確定性對跨聲速翼型模擬的影響。在該算例中，假定模型系數不再為常數，而是在一定區間內變化的隨機變量。關注重點是不確定性的傳播研究，即模擬結果的統計特性分析。

計算外形是RAE2822翼型[34]，本文的計算狀態是：馬赫數Ma=0.729，雷諾數Rec=6.5×106，迎角α=2.31°，計算使用的程序是課題組自行發展的MFlow[34]，使用的網格和算法可以參考文獻[35]。

計算使用SA一方程模型[36]，假定流動為全湍流，忽略了原始模型中的轉捩項(Trip Term)。因此模型中有9個常數，分別為：cb1、σ、cb2、κ、cw2、cw3、cv1、ct3、ct4。本文假定每個參數在各自的支撐集內為均勻分布，具體參數的取值范圍與文獻[34]保持一致。

文獻[37]通過完整PC展開得到了系統輸出的響應。這里提到的完整PC展開是指使用完全多項式展開，需要的樣本點數目為

(6)

式中：n為輸入隨機變量的維數；p為混沌多項式的階數；np定義為過采樣率。該算例中，假定展開多項式為2階，過采樣率np=2，因此完整PC展開需要的樣本點數目為N=110。

為了演示有噪聲信號的P1優化問題中截斷誤差ε的設定過程，圖1給出了采樣點N=50時，截斷誤差設定值和升力系數LOO誤差的關系。計算中設定ε從10-8逐漸增加到10-3，在此過程中LOO誤差先減小后增加。選取最小的LOO誤差對應的截斷誤差作為優化問題的設定值。在本文中分別嘗試了N=20，30，40，50，表1給出了不同采樣點下的截斷誤差設定值。

圖1 LOO誤差與截斷誤差設定值的關系 (N=50)Fig.1 LOO error variation with truncation error settings (N=50)

表1 不同采樣點下的截斷誤差設置值和LOO誤差

Table 1 Truncation error settings and LOO error for different samples

采樣點N截斷誤差設定值LOO誤差201.000×10-73.9256×10-13304.650×10-61.6669×10-10401.298×10-52.0523×10-9509.322×10-54.1732×10-9

確定好不同采樣點下優化問題的截斷誤差后，可以得到關注變量的稀疏PC展開。圖2給出了不同采樣點下得到的升力系數CL展開自由度，這里分別按自由度幅值大小排序。完整PC展開并不代表真正的精確解，只是作為精度較高的結果對比參考。首先從完整PC展開(圖中Full PC曲線)的結果可以看到，升力系數的PC展開存在明顯的稀疏特性，都只有少數若干項的自由度量值較大，其他展開項的自由度可以忽略不計。從圖可以看出，當樣本點數目分別為N=20，30，40，50時，本文發展的稀疏PC展開能夠比較準確地還原量級較大的若干個自由度，捕捉到輸出的主要特征。

稀疏PC展開的精度可以進一步通過預測值和CFD計算值對比驗證，表2給出了模型參數為標準取值時預測值和計算值的對比。不同采樣點下的預測值都跟CFD計算非常接近，說明稀疏PC展開也能夠準確預測系統輸出的響應。

圖2 不同采樣點下升力系數PC展開的自由度Fig.2 Freedoms of lift coefficient obtained with PC under different samples

表2 PC展開得到的升力、阻力系數與CFD比較

Table 2 Comparison of lift and drag coefficients obtained with PC and CFD

方法CLCDPC展開N=200.693190.013135N=300.693160.013152N=400.692950.013133N=500.693140.013126 完整PC0.693160.013130CFD0.693370.013137

在不確定性量化中，比較關注的是輸出的平均值和標準差等統計信息。圖3和圖4給出了不同采樣點下，升、阻力系數的統計信息，其中N=110代表完整PC展開的結果。從圖可以看出，升、阻力系數的平均值隨著樣本點數目的增加變化不大，標準差的變化稍大，但仍在可以接受的范圍內。

圖3 升力系數的平均值和標準差Fig.3 Mean value and standard deviation of lift coefficients

圖4 阻力系數的平均值和標準差Fig.4 Mean value and standard deviation of drag coefficients

在不確定性量化分析中，另一個關注的是每個隨機輸入變量對輸出變化的貢獻度。本文使用Sobol指標[38-39]來分析每個輸入變量對輸出方差的貢獻大小。圖5給出了9個輸入變量在升力系數分析中的Sobol指標。從圖中可以明顯看出，稀疏PC相對于完整的PC，也能夠比較準確地預測每個輸入變量的貢獻大小。隨著樣本點數目的增加，預測的Sobol指標變化不大，充分證明了稀疏PC展開的精度。cb1、σ、κ、cw2、cv1這5個參數對于升力系數的不確定性影響都相對較大，κ是其中貢獻最大的參數，ct3、ct4、cb2、cw3這4個參數的貢獻可以忽略不計。

圖5 升力系數預測中各輸入變量的Sobol指標Fig.5 Sobol indices for input parameters in lift coefficient prediction

2.2 炭化材料燒蝕熱響應不確定性傳播分析

第2個算例選取了一維燒蝕熱響應算例，以此來考核稀疏PC方法的適應性和精度。此處選取4th Abaltion Workshop官方網站提供的標準算例[40]進行分析，計算條件如圖6所示。

計算使用的控制方程包括組分的熱解過程、動量方程、質量守恒方程和能量方程，數值方法可以參考文獻[41]。分析中選擇材料25 mm厚度位置的溫度值作為目標變量，材料物性參數為輸入不確定性參數，進行不確定性傳播分析。

本算例所使用的材料為TACOT，是為方便數據對比而提供的一種假想材料，材料物性參數的不確定性也非本文的重點研究對象。因此此處對材料參數的不確定性，在參考以往相關工作的基礎上，采取如表3所示的假設值。

圖6 測試算例的計算條件Fig.6 Calculation condition of test case

表3 輸入隨機變量及其變化范圍Table 3 Input properties and their uncertainties

本文使用拉丁超立方抽樣方法，在輸入參數構成的9維隨機空間里抽樣獲得30個樣本點，通過稀疏PC方法構建了目標變量與隨機輸入參數的響應關系，借此來分析目標變量的統計特性。同時為了驗證稀疏PC方法的精度，使用蒙特卡洛采樣方法采樣500次，其統計結果作為參考值來進行對比分析。需要指出的是，蒙特卡洛方法非常依賴于樣本空間的大小，這里使用500個樣本進行直接統計分析可能會有偏差，因此也同時使用這些樣本點進行完整PC展開來進行統計分析。

表4給出了稀疏PC方法(N=30)、完整PC方法(N=500)和蒙特卡洛方法(N=500)預測的目標變量的均值和標準差。3種方法得到的統計量均相差不大，說明稀疏PC方法能夠較為準確地刻畫目標變量的統計特性。

9個物性參數中，究竟哪一些參數對于目標變量預測影響較大，需要進行敏感度分析。這里分別使用稀疏PC和完整PC得到每個輸入參數的Sobol指標，如表5所示，以此來衡量對輸出方差的貢獻大小。兩種方法得到的敏感性大小基本一致，x1和x5是最重要的兩個參數，其次是x2和x3，剩余的5個參數的貢獻很小。

表4 均值和標準差的比較

表5 輸入參數的敏感性分析Table 5 Sensitivity information of input properties

為了驗證稀疏PC方法的精度，進一步地在500個樣本點上進行評估，如圖7所示。其中橫坐標為通過CFD模擬得到的目標變量，縱坐標為通過稀疏PC方法預測的目標值，兩者呈現出很好的線性關系。定義R2為決定系數來評估預測的準確度：

(7)

為進一步分析該問題的不確定性傳播規律，對目標變量的概率密度函數進行了研究，結果如圖8所示。其中藍色分布是通過對目標變量的稀疏PC展開進行重采樣得到的。可以看出稀疏PC方法得到的目標變量概率密度分布與500個采樣點基本一致，都呈現出近似的高斯分布形態。

圖7 目標變量CFD模擬與PC預測的比較Fig.7 Comparison of CFD simulation and PC prediction for quantity of interest

圖8 目標變量的概率密度估計Fig.8 Probability density estimation for quantity of interest

3 結 論

為了解決多變量不確定性量化問題中的“維數災難”問題，本文發展了基于1范數最小的稀疏PC方法。該方法可以識別出若干個對系統輸出影響較大的基函數，從而得到近似的輸出響應。實施過程中通過交叉驗證方法確定有噪聲信號的P1規劃問題的截斷誤差，使用正交匹配追蹤算法求解該規劃問題。

通過SA模型系數的不確定性對RAE2822翼型跨聲速繞流模擬的影響和材料物性參數的不確定性對炭化材料燒蝕熱響應影響這兩個算例，證實了發展的稀疏PC方法能夠準確預測系統輸出。對于輸出的統計信息，包括均值、標準差和敏感性指標，稀疏PC都能夠取得和完整PC基本一致的結果。

稀疏PC方法證實能夠在樣本點數目小于自由度的不利條件下預測輸出的響應，為解決工程中多變量的不確定性量化問題提供了很好的解決方案。不過需要指出的是，稀疏PC方法適用的前提是輸出的PC展開存在稀疏特性。對于一個未知的復雜系統，如何使用稀疏性判斷準則預判系統的稀疏性需要進一步研究討論。對于某些各展開項自由度都有著不可忽視貢獻的問題，稀疏PC方法不再適用，需要發展其他更加高效的算法。