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

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

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

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

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

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

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

1 稀疏PC展開

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

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

(1)

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

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

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

(2)

式中：0范數(shù)表示α中非零元素的個數(shù)。然而P0優(yōu)化問題中的目標(biāo)函數(shù)是非凸函數(shù)，這是一個NP-hard問題，在實(shí)際中很難求解[33]。為此將P0優(yōu)化問題轉(zhuǎn)化為凸松弛的P1優(yōu)化問題，即

(3)

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

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

(4)

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

(5)

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

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

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

1) 定義殘差向量與每個基函數(shù)向量的內(nèi)積，即

并找到與當(dāng)前殘差最相關(guān)的特征：λ=argmaxjλj。

2) 更新迭代次數(shù)指標(biāo)和特征索引集：

k←k+1

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

2 不確定性量化分析

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

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

第1個算例是Spalart-Allmaras(SA)模型系數(shù)的不確定性對跨聲速翼型模擬的影響。在該算例中，假定模型系數(shù)不再為常數(shù)，而是在一定區(qū)間內(nèi)變化的隨機(jī)變量。關(guān)注重點(diǎn)是不確定性的傳播研究，即模擬結(jié)果的統(tǒng)計特性分析。

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

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

文獻(xiàn)[37]通過完整PC展開得到了系統(tǒng)輸出的響應(yīng)。這里提到的完整PC展開是指使用完全多項(xiàng)式展開，需要的樣本點(diǎn)數(shù)目為

(6)

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

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

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

表1 不同采樣點(diǎn)下的截斷誤差設(shè)置值和LOO誤差

Table 1 Truncation error settings and LOO error for different samples

采樣點(diǎn)N截斷誤差設(shè)定值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

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

稀疏PC展開的精度可以進(jìn)一步通過預(yù)測值和CFD計算值對比驗(yàn)證，表2給出了模型參數(shù)為標(biāo)準(zhǔn)取值時預(yù)測值和計算值的對比。不同采樣點(diǎn)下的預(yù)測值都跟CFD計算非常接近，說明稀疏PC展開也能夠準(zhǔn)確預(yù)測系統(tǒng)輸出的響應(yīng)。

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

表2 PC展開得到的升力、阻力系數(shù)與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

在不確定性量化中，比較關(guān)注的是輸出的平均值和標(biāo)準(zhǔn)差等統(tǒng)計信息。圖3和圖4給出了不同采樣點(diǎn)下，升、阻力系數(shù)的統(tǒng)計信息，其中N=110代表完整PC展開的結(jié)果。從圖可以看出，升、阻力系數(shù)的平均值隨著樣本點(diǎn)數(shù)目的增加變化不大，標(biāo)準(zhǔn)差的變化稍大，但仍在可以接受的范圍內(nèi)。

圖3 升力系數(shù)的平均值和標(biāo)準(zhǔn)差Fig.3 Mean value and standard deviation of lift coefficients

圖4 阻力系數(shù)的平均值和標(biāo)準(zhǔn)差Fig.4 Mean value and standard deviation of drag coefficients

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

圖5 升力系數(shù)預(yù)測中各輸入變量的Sobol指標(biāo)Fig.5 Sobol indices for input parameters in lift coefficient prediction

2.2 炭化材料燒蝕熱響應(yīng)不確定性傳播分析

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

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

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

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

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

本文使用拉丁超立方抽樣方法，在輸入?yún)?shù)構(gòu)成的9維隨機(jī)空間里抽樣獲得30個樣本點(diǎn)，通過稀疏PC方法構(gòu)建了目標(biāo)變量與隨機(jī)輸入?yún)?shù)的響應(yīng)關(guān)系，借此來分析目標(biāo)變量的統(tǒng)計特性。同時為了驗(yàn)證稀疏PC方法的精度，使用蒙特卡洛采樣方法采樣500次，其統(tǒng)計結(jié)果作為參考值來進(jìn)行對比分析。需要指出的是，蒙特卡洛方法非常依賴于樣本空間的大小，這里使用500個樣本進(jìn)行直接統(tǒng)計分析可能會有偏差，因此也同時使用這些樣本點(diǎn)進(jìn)行完整PC展開來進(jìn)行統(tǒng)計分析。

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

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

表4 均值和標(biāo)準(zhǔn)差的比較

表5 輸入?yún)?shù)的敏感性分析Table 5 Sensitivity information of input properties

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

(7)

為進(jìn)一步分析該問題的不確定性傳播規(guī)律，對目標(biāo)變量的概率密度函數(shù)進(jìn)行了研究，結(jié)果如圖8所示。其中藍(lán)色分布是通過對目標(biāo)變量的稀疏PC展開進(jìn)行重采樣得到的?？梢钥闯鱿∈鑀C方法得到的目標(biāo)變量概率密度分布與500個采樣點(diǎn)基本一致，都呈現(xiàn)出近似的高斯分布形態(tài)。

圖7 目標(biāo)變量CFD模擬與PC預(yù)測的比較Fig.7 Comparison of CFD simulation and PC prediction for quantity of interest

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

3 結(jié) 論

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

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

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