翼型不確定性量化中正交匹配追蹤的應(yīng)用

胡漢鐸，宋彥萍，俞建陽，劉瑤，陳浮，高文秀

哈爾濱工業(yè)大學 能源科學與工程學院，哈爾濱 150001

隨著航空技術(shù)的不斷發(fā)展，對飛行器的飛行性能要求不斷提高［1］。機翼作為飛行器升力的主要提供部件，其性能水平直接決定了飛行器的氣動性能和安全裕度［2］。然而，在實際加工、裝配和飛行過程中，由于制造誤差［3］、結(jié)冰［4-5］和腐蝕［6］等因素的存在，機翼的外形不可避免地與設(shè)計的理想形狀存在偏差，翼型幾何在一定的范圍內(nèi)隨機變化，即存在幾何不確定性。不確定性因素的存在會使機翼性能發(fā)生波動，進而導(dǎo)致飛行器性能惡化，甚至帶來安全隱患［7］，因此開展不確定性量化、獲得機翼性能參數(shù)在考慮不確定性因素下的隨機響應(yīng)特性，對于評估飛行器安全性能、指導(dǎo)機翼魯棒性設(shè)計具有重要研究意義［8］。

在不確定性量化方法方面，蒙特卡洛采樣（Monte Carlo Sampling，MCS）和混沌多項式展開（Polynomial Chaos Expansion，PCE）已經(jīng)被國內(nèi)外學者廣泛應(yīng)用于翼型的不確定性研究中。Loeven等［9］于2007年提出用概率配點法（Probabilistic Collocation Method，PCM）采樣求解混沌多項式系數(shù)，并將其用于NACA0012翼型的不確定性量化，只考慮來流馬赫數(shù)的變化，根據(jù)概率配點法，只需要5個配點即可量化翼型升力和阻力的概率密度分布；不久后，Loeven等［10］將概率配點法應(yīng)用于考慮3個幾何不確定性因素的NACA5412翼型不確定性量化，配點數(shù)增加到27個。高遠等［11］使用3階非嵌入式混沌多項式展開，針對考慮來流馬赫數(shù)和攻角不確定性的RAE2822翼型開展不確定性量化，使用25個樣本可以在效率與精度之間取得平衡。于佳鑫等［12］針對S809翼型，考慮了12個變量的幾何不確定性，對比了蒙特卡洛采樣和稀疏網(wǎng)格配點法獲得收斂量化結(jié)果所需的樣本數(shù)和精度，對于蒙特卡洛采樣方法，當采樣數(shù)達到2 000時才收斂，而3階Hermite多項式需要313個樣本就可以得到與蒙特卡洛采樣一致的結(jié)果?？紤]到蒙特卡洛采樣需要大量樣本，Liu等［13］利用代理模型技術(shù)，首先用100個訓練樣本構(gòu)建Kriging代理模型，該模型可以代替耗時的CFD仿真預(yù)測RAE5243翼型上不同壓敏涂層厚度下的升力系數(shù)，誤差不超過0.6%，然后使用容易計算的代理模型進行了107次蒙特卡洛采樣，獲得了幾何不確定性對翼型表面流動的影響。相比而言，蒙特卡洛采樣直接模擬實際隨機過程，獲得的統(tǒng)計量精度高，可信度好，但樣本量需求十分龐大；混沌多項式展開將翼型性能參數(shù)的隨機響應(yīng)在混沌多項式基底上展開，利用多項式系數(shù)計算均值和標準差，樣本數(shù)需求有所降低，而精度與蒙特卡洛采樣相近。盡管如此，通過概率配點法求解多項式系數(shù)要求樣本數(shù)量不少于多項式基底數(shù)，限制了其在考慮高維不確定性因素的復(fù)雜系統(tǒng)中的實際應(yīng)用。

近年來，壓縮感知技術(shù)在信號處理領(lǐng)域取得迅速發(fā)展［14］，其利用原始信號的稀疏性質(zhì)，只需要很少的觀測樣本就可以用特定的重構(gòu)算法精確重構(gòu)原始信號［15］。針對不確定性量化中樣本數(shù)需求較多的問題，研究人員開始探索將壓縮感知技術(shù)應(yīng)用于不確定性量化中。Liang等［16］詳細探討了壓縮感知理論對研究對象的基本要求，包括系數(shù)的稀疏性和觀測矩陣的受限等距性質(zhì)（Restricted Isometry Property，RIP），指出混沌多項式展開在實際不確定性量化研究中滿足上述條件，并通過測試函數(shù)驗證。Yang等［17］利用壓縮感知重構(gòu)混沌多項式系數(shù)，并在若干測試函數(shù)和偏微分方程組上測試，輸入變量維度涵蓋20～500范圍，所提出的方法即使在小樣本集上也能夠獲得良好的精度。陳江濤等［18］利用稀疏混沌多項式展開方法，在翼型繞流湍流模型稀疏不確定性和燒蝕材料物性參數(shù)不確定性2個算例上進行研究，對于輸出量的統(tǒng)計特性稀疏混沌多項式展開能夠進行準確預(yù)測。Lu等［19］將壓縮感知中的正交匹配追蹤算法（Orthogonal Matching Pursuit，OMP）用于重構(gòu)混沌多項式系數(shù)，分別在考慮幾何不確定性的NLR7301翼型流動和考慮壁面熱流不確定性的強迫對流換熱上對比驗證，結(jié)果表明，OMP算法的優(yōu)勢在于達到相同的精度要求時收斂更快，觀測樣本量需求更少，適合復(fù)雜問題的不確定性分析。趙歡等［20］采用稀疏混沌多項式重構(gòu)算法，對高速自然層流翼型LRN1015開展考慮來流馬赫數(shù)和升力系數(shù)不確定性的穩(wěn)健性優(yōu)化，得益于自適應(yīng)前向-后向選擇的稀疏重構(gòu)算法，不確定性分析效率得到改進，獲得了更好的穩(wěn)健性設(shè)計結(jié)果。盡管在壓縮感知技術(shù)與混沌多項式展開相結(jié)合方面國內(nèi)外學者已經(jīng)開展了初步探索，但在翼型幾何不確定性量化領(lǐng)域的研究仍相對較少。

本文考慮RAE2822翼型的幾何不確定性，將壓縮感知技術(shù)用于氣動性能參數(shù)隨機響應(yīng)的不確定性量化，利用正交匹配追蹤算法重構(gòu)混沌多項式系數(shù)，并且將獲得的統(tǒng)計特性與傳統(tǒng)的蒙特卡洛采樣和滿秩概率配點法進行對比分析，探討正交匹配追蹤重構(gòu)算法在收斂性能和樣本需求等方面的優(yōu)勢，為壓縮感知技術(shù)應(yīng)用于翼型不確定性量化領(lǐng)域的進一步研究提供參考。

1 不確定性量化方法

1.1 混沌多項式展開

不確定性量化中的混沌多項式展開（Polynomial Chaos Expansion，PCE）方法將系統(tǒng)不確定性輸入和隨機響應(yīng)之間的映射關(guān)系在一系列正交多項式基底上展開。盡管現(xiàn)實中大多數(shù)信號都不具備稀疏性質(zhì)，但經(jīng)過混沌多項式展開后的系數(shù)將逐漸衰減，因此可以認為具有一定的稀疏性。

混沌多項式展開可以表示為［21］

式中：y為系統(tǒng)隨機響應(yīng)輸出；x為系統(tǒng)不確定性輸入；Hi為與不確定性分布相對應(yīng)的混沌多項式基底。在本文研究范圍內(nèi)，假定幾何不確定性因素滿足正態(tài)分布，對應(yīng)的多項式基底為Hermite基底，其形式可以按照遞推公式［22］推導(dǎo)；ai為第i項基底對應(yīng)的多項式系數(shù)，需要通過采樣獲得若干組系統(tǒng)輸入-輸出響應(yīng)值后，利用最小二乘法求解：

其中：H可以寫為矩陣形式，行數(shù)為采樣樣本數(shù)N，列數(shù)（Npc+1）為總共的多項式項數(shù)：

Npc與不確定性輸入變量的維數(shù)m和多項式基底階數(shù)p有關(guān)：

從式（4）可見，隨著輸入維數(shù)的多項式基底階數(shù)的增加，基底總數(shù)呈階乘增長，當不確定性因素較多、系統(tǒng)響應(yīng)非線性較強時，多項式基底個數(shù)劇增，求解系數(shù)所需的樣本數(shù)也隨之增多，導(dǎo)致計算成本過高，即面臨“維數(shù)災(zāi)難”的困境。

1.2 壓縮感知理論

在求解多項式系數(shù)a時，研究人員希望需要的樣本量盡可能少，以提高不確定性量化效率。當樣本數(shù)N小于多項式基底項數(shù)時，式（3）中的H矩陣行數(shù)小于列數(shù)，式（1）的求解轉(zhuǎn)變?yōu)榍范ǚ匠探M問題，壓縮感知理論正是關(guān)注此類問題的求解。從線性方程組的角度考慮，欠定方程組沒有唯一解，但壓縮感知理論證明，如果待求解向量滿足一定的條件，例如具有稀疏性，那么通過特定的重構(gòu)算法可以獲得唯一解。對于不確定性量化問題而言，如果系統(tǒng)響應(yīng)在混沌多項式基底上展開的系數(shù)足夠稀疏，那么系數(shù)向量可以通過如下l0范數(shù)最小化問題來求解：

式中：H矩陣起到觀測矩陣的作用；||a||0表示系數(shù)向量a中非零元素的個數(shù)，反映了系數(shù)向量的稀疏程度；ε為允許的擬合容差。

然而，直接求解上述l0范數(shù)最小化問題是一個NP-hard難題。為此，開發(fā)了一系列重構(gòu)算法，以可接受的計算成本獲得近似的最稀疏解，包括凸優(yōu)化算法、貪婪類算法等。其中，正交匹配追蹤（Orthogonal Matching Pursuit，OMP）是貪婪類重構(gòu)算法的代表，其主要優(yōu)勢是計算成本低、收斂速度快，適合交叉驗證。

1.3 正交匹配追蹤算法

重構(gòu)算法是壓縮感知能否精確重構(gòu)原始信號的核心，根據(jù)應(yīng)用場景選擇合適的重構(gòu)算法至關(guān)重要。本文采用貪婪重構(gòu)算法中的正交匹配追蹤重構(gòu)混沌多項式系數(shù)。

由于遍歷H矩陣中所有列向量基底組合過于復(fù)雜，正交匹配追蹤算法采用貪婪思想，每次從列向量基底中選擇與殘差最相關(guān)的基底，添加到已有的支持集中，在選出的基底上嘗試重構(gòu)多項式系數(shù)，系數(shù)向量中只有選出的基底對應(yīng)的系數(shù)不為零，而其余系數(shù)均為零，以此實現(xiàn)重構(gòu)系數(shù)的稀疏性。上述步驟迭代若干次，直至殘差降至指定的容差范圍內(nèi)，算法流程如圖1所示。

圖1 正交匹配追蹤算法流程Fig.1 Flowchart of OMP

由于在每一個迭代步，OMP算法都是尋找能使殘差下降最快的基向量，因此能夠以較快的速度收斂到唯一的稀疏解。

通過重構(gòu)算法可以獲得系數(shù)向量a，進而獲得混沌多項式展開完整的表達式，可以根據(jù)式（1）計算系統(tǒng)輸出響應(yīng)。對于不確定性量化問題，研究人員更加關(guān)注隨機響應(yīng)的統(tǒng)計特性，包括均值μ和標準差σ，兩者可以直接通過混沌多項式系數(shù)分量來估計：

1.4 標準測試算例驗證

綜合上述方法，建立了一個可以減少樣本數(shù)需求的不確定性量化框架，如圖2所示。首先確定不確定性的輸入維度和概率分布，然后根據(jù)概率分布和多項式階數(shù)生成混沌多項式基底，所需的樣本從高斯積分點中隨機選取，確定觀測矩陣H。對于翼型流動問題，通過計算流體力學（Computational Fluid Dynamics，CFD）仿真獲得各個樣本下的氣動力系數(shù)和其他關(guān)注的流場量，最后在建立的樣本集上使用正交匹配追蹤重構(gòu)算法，求解混沌多項式系數(shù)，并獲得所關(guān)注流動參數(shù)的統(tǒng)計特性。

圖2 不確定性量化流程Fig.2 Uncertainty quantification framework

為了驗證上述不確定性量化框架的精度，本文首先在Ishigami測試函數(shù)上進行算例驗證。該函數(shù)具有很強的非線性，在不確定性量化研究中廣泛用于精度校驗［23］。

該函數(shù)的定義為

其均值μ和標準差σ可以通過理論分析精確計算，作為評估不確定度量化誤差的基值。將上述不確定性量化框架用于該測試算例，并與理論值和傳統(tǒng)的不確定性量化方法，包括蒙特卡洛采樣和滿秩概率配點法（Full Rank Probabilistic Collocation，F(xiàn)RPC）進行對比。

圖3給出了均值和標準差相對誤差隨樣本數(shù)增加的變化過程。相比之下，3種方法中蒙特卡洛采樣收斂最慢，重構(gòu)誤差隨著樣本數(shù)的增加緩慢降低。混沌多項式展開基底輸入維數(shù)為3，多項式階數(shù)取13，根據(jù)式（4）可知共有560個基底，當變量數(shù)小于基底數(shù)時，方程組欠定，難以獲得精確解，故均值和標準差的重構(gòu)誤差都較大；而當樣本數(shù)等于基底數(shù)時，方程組滿秩，可以用最小二乘法求得精確解，因而重構(gòu)誤差曲線陡降。而采用同樣混沌多項式基底的正交匹配追蹤算法，僅需140個樣本即可收斂到與滿秩概率配點法相同的精度水平，這是因為560個混沌多項式基底對應(yīng)的系數(shù)向量具有稀疏性，即大部分系數(shù)幅值接近于0，表明這些系數(shù)對應(yīng)的多項式基底影響作用較小，而正交匹配追蹤算法可以優(yōu)先篩選占主要影響的基底并求解對應(yīng)的系數(shù)，從而實現(xiàn)用較少的樣本獲得精確的不確定性量化結(jié)果。

圖3 蒙特卡洛采樣、滿秩概率配點和正交匹配追蹤算法對Ishigami測試算例的預(yù)測誤差Fig.3 Prediction errors of Ishigami function by MCS，F(xiàn)RPC and OMP

2 RAE2822不確定性量化分析

在通過測試算例驗證所用方法的精度后，將上述框架用于RAE2822翼型的幾何不確定性量化分析。首先對所采用的數(shù)值模型進行驗證，將仿真結(jié)果與試驗結(jié)果對比。為了描述翼型幾何不確定性變化，采用類函數(shù)/形函數(shù)變換（Class-Shape Transformation，CST）將 原 始 翼 型 參 數(shù)化，將翼型幾何的不確定性變化轉(zhuǎn)變?yōu)镃ST參數(shù)的隨機波動，并確定各參數(shù)概率分布規(guī)律和變化范圍。最后，在本文所研究的翼型幾何不確定性因素下，用混沌多項式展開和正交匹配追蹤算法開展不確定性量化，獲得翼型升力系數(shù)、阻力系數(shù)和其他主要流場量的統(tǒng)計特性，并與蒙特卡洛采樣和概率配點法在收斂性、量化精度和樣本數(shù)需求等方面進行對比，以驗證所用方法的有效性。

2.1 數(shù)值模型驗證

RAE2822翼型繞流是典型的二維跨聲速可壓縮流動，被歐洲EUROVAL項目和AGARD選為經(jīng)典驗證算例。Cook等［24］在RAE 8 ft×6 ft（1 ft=0.304 8 m）跨聲速風洞中進行流場實驗測量，共選擇10個工況，通過壓力探針測量壓力分布，而升力、阻力和力矩系數(shù)通過對翼型表面壓力進行積分獲得。此外，還用邊界層探針測量了邊界層速度剖面。

RAE2822翼型的幾何形狀以上下表面離散點的形式表示，為了描述翼型幾何形狀的不確定性變化，需要用參數(shù)化方法將離散的翼型幾何點轉(zhuǎn)化為控制參數(shù)。本文采用Kulfan［25］提出的類形狀/形函數(shù)變換方法來擬合原始翼型。CST方法使用類函數(shù)C（φ）和形函數(shù)S（φ）來定義翼型的幾何形狀：

式中：φ和ζ是基于弦長歸一化的翼型坐標；ζT是尾緣厚度，對于RAE2822類型的圓頭尖尾翼型，尾緣厚度為0。

C（φ）是用于定義翼型整體類型的類函數(shù)，其表達式為

通過調(diào)整N1和N2，可以生成不同類型的翼型，如圖4所示。

圖4 類函數(shù)C（φ）決定的不同N1、N2取值對應(yīng)的基本翼型Fig.4 Basic airfoil determined by class function C（φ）with different N1， N2

對于前緣圓形、尾緣尖形的RAE2822翼型，N1取0.5，N2取1。形函數(shù)S（φ）通過多階Bernstein多項式的疊加來描述翼型的具體形狀：

式中：Ai是第i項的權(quán)重系數(shù)；n是Bernstein多項式 階 數(shù)（Bernstein Polynomial Order，BPO）。CST參數(shù)化方法即是通過控制形函數(shù)中的Ai來擬合成不同形狀的翼型，擬合誤差隨著多項式階數(shù)的增加而減小，但與此同時變量Ai的個數(shù)也隨之增加。

圖5給出了不同階數(shù)下的對原始翼型的擬合誤差，可見4階CST已經(jīng)能夠在風洞模型公差范圍內(nèi)擬合RAE2822翼型，可以滿足實際應(yīng)用需要。

圖5 不同Bernstein基函數(shù)階數(shù)的CST方法對RAE2822翼型的擬合誤差Fig.5 Fitting error of RAE2822 airfoil using CST with different BPO

采用4階CST參數(shù)化，上下表面各需要5個參數(shù)（A0～A9）來描述無量綱翼型形狀；此外，為了保持與實驗?zāi)Ｐ鸵恢乱则炞C數(shù)值模型精度，同時考慮弦長不確定性的變化，翼型的弦長c也作為幾何參數(shù)。因此，總共確定了11個幾何參數(shù)。

為了驗證數(shù)值模型，首先進行確定性仿真。在不確定性量化研究中，通常采用來流馬赫數(shù)Ma為0.734，攻角為2.79°，對應(yīng)雷諾數(shù)Re為6.5×106的工況［26］。在數(shù)值仿真驗證時，選取與其最為接近的實驗測量工況進行仿真，對應(yīng)的Ma為0.730，攻角為3.19°。為了同時驗證經(jīng)過CST擬合的翼型氣動性能，對原始翼型和擬合翼型均進行雷諾平均納維-斯托克斯方程（Reynold-Averaged Navier-Stokes equations，RANS）模擬，計算域邊界距離翼型表面25倍弦長，以避免計算域邊界對翼型附近流動的干擾。對計算域劃分C型拓撲網(wǎng)格，邊界層第1層網(wǎng)格高度由y+值控制，本研究采用SSTk-ω湍流模型，要求壁面上y+＜1。網(wǎng)格節(jié)點數(shù)通過網(wǎng)格無關(guān)性驗證確定，進行了一系列網(wǎng)格單元總數(shù)逐漸遞增的仿真，表1列出了升力系數(shù)（Lift coefficient，CL）和阻力系數(shù)（Drag coefficient，CD）與實驗結(jié)果的對比。表中結(jié)果表明，當網(wǎng)格單元數(shù)達到120 000后，計算結(jié)果隨網(wǎng)格單元的進一步增加變化不大。因此，在后續(xù)仿真將基于129 258單元數(shù)的網(wǎng)格開展，以平衡仿真精度和時間成本。該套網(wǎng)格在翼型附近區(qū)域的分布情況如圖6所示，在翼型上表面中間弦長后的激波存在區(qū)域進行網(wǎng)格節(jié)點加密，以提高對激波的辨識度。

表1 網(wǎng)格無關(guān)性驗證Table 1 Grid independence test

圖6 RAE2822翼型計算網(wǎng)格Fig.6 Computation grid around RAE2822 airfoil

在確定了網(wǎng)格單元數(shù)后，用相同的網(wǎng)格拓撲對擬合翼型也劃分了計算域網(wǎng)格，圖7中給出了原始翼型和擬合翼型表面壓力系數(shù)（Pressure coefficient，Cp）分布與實驗結(jié)果的對比，結(jié)果表明，經(jīng)過CST擬合的翼型仿真結(jié)果與原始翼型取得一致，從仿真結(jié)果層面驗證了擬合精度。

圖7 原始翼型和擬合翼型壓力系數(shù)的對比Fig.7 Comparison of original and fitting airfoils on pressure coefficient

2.2 幾何不確定性范圍和概率分布

在驗證數(shù)值模型的基礎(chǔ)上，可以開展考慮幾何不確定性的RAE2822翼型氣動性能不確定性量化。由于翼型幾何已經(jīng)通過CST方法進行參數(shù)化，因此翼型幾何不確定性的變化可以轉(zhuǎn)化為CST參數(shù)的隨機波動。參考其他有關(guān)翼型幾何不確定性的研究［27］，針對加工、裝配和實際飛行工況等因素導(dǎo)致的實際翼型與原始翼型差異，將翼型弦長的不確定性變化也考慮在內(nèi)，確定了本研究中RAE2822的幾何厚度變化范圍為±1 mm。圖8中給出翼型前緣和尾緣的不確定性變化邊界。

圖8 RAE2822在前緣和尾緣的幾何不確定性變化范圍Fig.8 Geometrical uncertainty range of RAE2822 airfoil at leading edge and trailing edge

11個CST幾何參數(shù)的范圍由圖中最薄和最厚的翼型確定，這些參數(shù)假設(shè)均服從高斯分布N（μ，σ）。為了實現(xiàn)CST參數(shù)隨機變化時翼型幾何能夠保持在最薄和最厚翼型之間，采用最小二乘法分別反求加厚翼型和減薄翼型對應(yīng)的CST參數(shù)值，由此確定了11個CST參數(shù)能夠變化的上下限范圍。根據(jù)3σ原則，假定此范圍為±3σ的邊界，由此確定了各變量標準差σ的取值，均值μ則由原始翼型反求CST參數(shù)獲得。各CST參數(shù)的均值μ和標準差σ在表2中列出。

表2 11個CST參數(shù)的均值μ和標準差σTable 2 Mean μ and standard deviation σ of 11 CST parameters

圖9給出了按照上述CST參數(shù)化方法生成的一系列隨機翼型，這些翼型分布在加厚翼型和減薄翼型之間，其分布規(guī)律表明所采用的方法生成的隨機翼型可以覆蓋±1 mm的幾何不確定性變化區(qū)間，可以用于描述翼型幾何不確定性變化。

圖9 隨機翼型幾何分布圖Fig.9 Distribution of random airfoils geometry

幾何不確定性變化范圍確定后，采用圖2所示的不確定性量化框架，用正交匹配追蹤算法求解混沌多項式系數(shù)并獲得RAE2822翼型氣動性能的統(tǒng)計特性，與蒙特卡洛采樣（MCS）和滿秩概率配點法（FRPC）對比。

2.3 氣動力系數(shù)對比

表3中列出了3種方法求解的升力系數(shù)和阻力系數(shù)的統(tǒng)計特性。由于對于實際翼型的不確定性量化問題尚無理論解，考慮到蒙特卡洛采樣法模擬真實概率分布采樣，樣本信息豐富，因此以蒙特卡洛采樣獲得的結(jié)果為參考計算相對誤差。

表3 蒙特卡洛采樣、滿秩概率配點法和正交匹配追蹤預(yù)測的氣動力系數(shù)均值μ和標準差σTable 3 Mean μ and standard deviation σ of aerodynamic coefficients predicted by MCS， FRPC and OMP

由表3可知，在幾何不確定性作用下，翼型升力系數(shù)均值相比于原型下降了約1%，表明幾何不確定性的存在導(dǎo)致翼型升力降低；另一方面，幾何不確定性導(dǎo)致升力系數(shù)和阻力系數(shù)出現(xiàn)一定程度的波動，其中阻力系數(shù)變化的標準差可以達到均值的4%。比較而言，滿秩概率配點法所獲得的升力系數(shù)和阻力系數(shù)與蒙特卡洛采樣較為接近，但需要364個樣本才能使H矩陣滿秩；而正交匹配追蹤算法只需要120個樣本，即可取得與滿秩概率配點法基本一致的結(jié)果，表明壓縮感知方法可以用相對較少的樣本獲得較高的精度。

上述結(jié)果體現(xiàn)了樣本量充足時各方法的預(yù)測精度，而圖10中則給出不同樣本數(shù)下的相對誤差來比較各方法的收斂性，其中相對誤差是以最多樣本數(shù)下蒙特卡洛采樣的結(jié)果為基準計算的。3種方法中，蒙特卡洛采樣收斂速度最慢，幾乎在所有樣本數(shù)水平下相對誤差都是最高的。基于概率配點法的PCE方法直到樣本數(shù)達到364，即滿足滿秩條件時才快速收斂，364個樣本的結(jié)果與1 200個樣本量的蒙特卡洛采樣方法相比，相對誤差減少了一個數(shù)量級；而在樣本數(shù)小于多項式系數(shù)個數(shù)時，由于方程組欠定，求解誤差較大。

圖10 蒙特卡洛采樣、滿秩概率配點和正交匹配追蹤算法收斂情況對比Fig.10 Convergence comparison of MCS， FRPC and OMP

比較而言，正交匹配追蹤算法可以用更少的樣本獲得與概率配點法相接近的精度，特別是當樣本數(shù)低于滿秩要求，相比于概率配點法用最小二乘擬合多項式系數(shù)，正交匹配追蹤的誤差更低或持平。不過，在樣本數(shù)增加到一定規(guī)模之前，正交匹配追蹤算法和概率配點法之間的差距并不明顯。具體對于本算例而言，均值μ的預(yù)測出現(xiàn)顯著差異是在樣本數(shù)達到30以上，而對于標準差σ則是60個樣本，表明標準差的準確預(yù)測比均值難度更大。

正交匹配追蹤算法能夠表現(xiàn)出上述優(yōu)勢得益于多項式系數(shù)的稀疏性。盡管在實際應(yīng)用中假設(shè)了稀疏性的前提，但重構(gòu)獲得的多項式系數(shù)是否足夠稀疏仍然需要校驗。圖11對比了2種PCE方法中多項式基底對應(yīng)系數(shù)的幅值，按基底序號排序。圖中的2組系數(shù)分別是364個樣本的滿秩概率配點法PCE和120個樣本的正交匹配追蹤獲得的，為了直觀顯示多項式系數(shù)的相對大小，橫縱坐標軸都采用對數(shù)刻度。圖中系數(shù)幅值分布表明，大部分PCE系數(shù)都接近于0，如果假定小于10-6的系數(shù)可被視為0，那么系數(shù)向量確實滿足稀疏性條件。從另一角度分析，所有通過正交匹配追蹤求解出的系數(shù)幅值都大于10-6，這表明該算法能夠從所有的多項式基底中篩選出起主導(dǎo)作用的子集，從而利用系數(shù)向量的稀疏性精確重構(gòu)對應(yīng)的非零系數(shù)。

圖11 滿秩概率配點和正交匹配追蹤重構(gòu)的混沌多項式系數(shù)對比Fig.11 Comparison of polynomial chaos coefficients reconstructed by FRPC and OMP

除了對比均值和標準差等統(tǒng)計特性之外，還需要進一步驗證概率密度分布。圖12中給出了采用蒙特卡洛采樣、滿秩法PCE和正交匹配追蹤預(yù)測的升力系數(shù)和阻力系數(shù)的概率密度分布函數(shù)（Probability Density Function，PDF）。圖中結(jié)果表明，滿秩概率配點法PCE和正交匹配追蹤能夠精確重構(gòu)出與蒙特卡洛采樣結(jié)果一致的氣動力系數(shù)概率分布，且正交匹配追蹤只需要120個觀測樣本，遠少于需要364個樣本才能滿秩的PCE和3 000個樣本才能收斂的蒙特卡洛采樣。

圖12 蒙特卡洛采樣、滿秩概率配點和正交匹配追蹤算法預(yù)測的升力系數(shù)、阻力系數(shù)概率密度函數(shù)對比Fig.12 Comparison of probability density function of CL and CD predicted by MCS， FRPC and OMP

綜上，在考慮幾何不確定性的翼型氣動力系數(shù)的預(yù)測中，壓縮感知技術(shù)在樣本需求量少方面展現(xiàn)出顯著優(yōu)勢，能夠高效靈活地進行不確定性量化分析。

2.4 流場參數(shù)對比

氣動力系數(shù)是將翼型表面所受壓力積分獲得的，為了進一步對比流場重構(gòu)結(jié)果，提取翼型上下表面壓力系數(shù)沿弦長的分布，如圖13所示。

圖13 翼型表面壓力系數(shù)沿弦長的均值和標準差分布Fig.13 Mean and standard deviation of Cp on airfoil surfaces along chord

RAE2822作為典型的超臨界翼型，氣流在翼型上表面加速降壓，從前緣起有一個較為穩(wěn)定的低壓區(qū)，直到出現(xiàn)激波后壓力驟升，而在下表面上壓力系數(shù)沒有突變，且在翼型尾緣附近壓力有所升高，以提供更大升力。在均值方面，滿秩概率配點法和正交匹配追蹤都可以獲得與蒙特卡洛采樣相同的壓力系數(shù)分布規(guī)律，結(jié)果一致性較好；而在標準差方面，對于壓力波動不明顯的區(qū)域，包括翼型上表面沒有激波的區(qū)域和整個翼型下表面，標準差σ的預(yù)測結(jié)果也很接近，3種方法的主要差異體現(xiàn)在60%弦長附近的激波所在區(qū)域。幾何不確定性的存在使得翼型偏離設(shè)計形狀，不但會使激波的位置發(fā)生改變，激波強度也會因波前氣流加速程度不同而變化。由于氣流經(jīng)過激波后壓力驟升，激波強度和位置的變化將導(dǎo)致激波所在區(qū)域壓力系數(shù)變化劇烈，即標準差增大，同時不同不確定性量化方法之間標準差預(yù)測的差異也被放大。相比而言，在標準差最大的位置，正交匹配追蹤的結(jié)果更加接近蒙特卡洛采樣的參考值，而滿秩法PCE的結(jié)果則高估了壓力系數(shù)的波動程度。除了激波所在區(qū)域之外，翼型前緣壓力波動也較為顯著，這是因為翼型下表面駐點位置因幾何形狀發(fā)生微小改變而移動。

幾何不確定性的存在不僅影響著翼型表面壓力，同時也改變了翼型附近的其他流場參數(shù)。在圖14中，給出了蒙特卡洛采樣、滿秩法和正交匹配追蹤算法重構(gòu)的馬赫數(shù)均值和標準差分布，與前述的氣動力系數(shù)不同，流場馬赫數(shù)的重構(gòu)是在每個網(wǎng)格節(jié)點上進行的。正交匹配追蹤算法對馬赫數(shù)均值的預(yù)測在大多數(shù)位置上與蒙特卡洛采樣非常接近，上表面空氣從前緣開始加速至超聲速，經(jīng)過激波后變?yōu)閬喡曀?，馬赫數(shù)驟變。標準差的差異主要體現(xiàn)在60%弦長附近的激波區(qū)域，這也與壓力系數(shù)分布呈現(xiàn)的規(guī)律一致，幾何不確定性影響了激波出現(xiàn)的位置以及激波前后馬赫數(shù)差值，重構(gòu)激波區(qū)域內(nèi)馬赫數(shù)的概率特性也更具挑戰(zhàn)。正交匹配追蹤的結(jié)果表明，翼型附近區(qū)域的流場參數(shù)變化同樣可以通過稀疏重構(gòu)的方式以一定的精度準確預(yù)測，且所需樣本量更少，從而在不確定性量化中展現(xiàn)出顯著的優(yōu)勢。

3 結(jié) 論

本文將壓縮感知技術(shù)與不確定性量化相結(jié)合，對正交匹配追蹤算法用于RAE2822翼型幾何不確定性量化進行了研究，可以得到以下結(jié)論：

1）飛行器機翼在加工、裝配和實際飛行過程中，會因為各種因素存在幾何外形的不確定性，將RAE2822翼型經(jīng)過CST方法參數(shù)化后，翼型幾何的不確定性可以轉(zhuǎn)換為CST參數(shù)的隨機波動。

2）經(jīng)過不確定性量化研究，幾何外形不確定性的存在會導(dǎo)致翼型升力系數(shù)降低，且翼型上表面激波的強度和位置發(fā)生變化，激波所在區(qū)域的流場參數(shù)出現(xiàn)波動，影響飛行器的正常飛行。

3）通過測試函數(shù)和RAE2822翼型驗證，與傳統(tǒng)的蒙特卡洛采樣、滿秩概率配點法等不確定性量化方法相比，正交匹配追蹤可以用相對較少的樣本獲得與相同精度水平的量化結(jié)果。

4）正交匹配追蹤算法用稀疏重構(gòu)的方式求解混沌多項式系數(shù)，不僅能夠預(yù)測氣動力系數(shù)的統(tǒng)計特性和概率分布，還可以重構(gòu)翼型附近的流場參數(shù)變化規(guī)律，有助于進一步研究不確定性對翼型流動的影響。

通過測試函數(shù)和RAE2822翼型幾何不確定性的研究，驗證了壓縮感知技術(shù)用于不確定性量化在效率、收斂性和精度等方面的優(yōu)勢。后續(xù)的工作將在此基礎(chǔ)上開展參數(shù)敏感度分析，進一步研究不同幾何不確定性參數(shù)對翼型氣動性能的影響程度，為開展翼型穩(wěn)健性設(shè)計提供參考。