基于新型高維代理模型的氣動外形設計方法

趙歡，高正紅，夏露

西北工業大學 航空學院，西安 710072

飛行器氣動外形優化對象越來越復雜，需要考慮飛行器各部件甚至各學科之間的影響，這使得氣動設計不得不轉向全局精細優化，也導致優化建模所需的設計變量個數急劇增加，遠超代理模型和優化算法的發展速度，使代理優化算法面臨嚴重的“維度災難”難題，其主要來源于高維代理模型的困境［1-2］。大規模設計變量一方面使基于代理模型的全局優化算法收斂性降低，計算復雜度和計算時間劇烈增加；同時使得代理模型預測精度急劇下降，達到相同樣本密度所需的樣本數以指數形式增加，從而導致代理優化算法難以尋找到全局最優解［2-3］。目前對于大多數復雜氣動外形優化問題（包含多目標和復雜約束要求，同時減小激波阻力和黏性阻力），其設計呈現復雜的多峰或多極值、強非線性等特征［4］，如跨聲速機翼優化［5-6］、層流超臨界翼型優化［7-9］、多目標旋翼翼型優化［10-11］等，獲得全局最優解而不是局部可行解是目前追求的目標。一些研究者試圖通過多點開始的梯度優化方法尋找全局最優解，然而初始點的選擇及選樣方式對最優結果有明顯影響。因此，面對越來越嚴苛的氣動性能要求，提高代理優化算法對更高維問題的優化效率和能力勢在必行。

目前，在工程問題中廣泛使用具有不同特點的代理模型［12-14］，如Kriging、徑向基函數（Radial Basis Function，RBF）、支持向量回歸（Support Vector Regression，SVR）及多項式混沌展開（Polynomial Chaos Expansion，PCE）等，均隨著設計變量維數的增加，達到相同預測精度所需樣本數不同程度地急劇增加，能達到的最好精度（突破樣本數限制）甚至也在不同程度地嚴重變差，直至無法滿足優化設計精度的要求。一些設計者提出在建立代理模型過程中引入梯度信息以增強高維空間預測精度，如直接梯度增強Kriging（Gradient-Enhanced Kriging，GEK）和間接GEK 模型［15］，權重梯度增強Kriging （Weight Gradient-Enhanced Kriging，WGEK）模型［16］，梯度增強PCE［17］等。然而，目前這些梯度增強代理模型應用到氣動優化中，其有效維數最高僅在60 ～ 100 維，更高的維數將使這些“改進的”代理模型也會面臨嚴峻的“維度災難”挑戰［2］。因此，代理模型的發展速度與日益提高的精細化氣動優化設計要求難以相匹配，甚至成為阻礙代理優化算法向高維氣動優化及多學科優化發展的瓶頸。

高維代理模型技術目前發展主要包括3 個方向，高維模型代表（High-Dimensional Model Representation，HDMR）、降維代理建模（Dimension-Reduction Surrogate Model，DRSM）、深度神經網 絡（Deep Neural Network，DNN）。 其 中，HDMR 最 早 由Sobol［18］提 出，被 引 入 到 高 維 問題［19］。Cut-HDMR 在中心點附近展開為一系列階次和變量數逐漸增加的無窮項級數之和，隨著階次的增加但貢獻迅速減小。HDMR 與代理模型相結合，如Kriging-HDMR，RBF-HDMR，SVRHDMR 等，獲得了對高維輸出響應的近似。HDMR 方法雖然已得到了一些應用，如渦輪葉片優化和機翼優化等［20-21］，但隨著設計變量個數的增多，用于近似HDMR 的低維代理模型個數也在快速增加。例如，對于Cut-HDMR 保留到l階展開，其時間花費為其中s代表沿X每個方向所取樣本點數。對于更加復雜的工程問題（包含高階項交互），保留高階項則要求更多的計算花費，樣本點數量呈幾何級數增加。因此，HDMR 仍是一種昂貴的高維代理建模途徑。

與HDMR 相比，DRSM 技術在控制計算成本前提下擁有較好的近似效果。DRSM 通常包含兩步，即變量降維和代理建模，其重點在第1步。一些研究者提出將輸入空間樣本投影（映射）到特征空間，特征空間能保留原始空間最多的信息。根據映射方式的不同可以分為線性和非線性降維、監督和非監督學習。經典的線性非監督降維方法如Karhunen-Loève（K-L）展開［22］、主 成 分 分 析（Principal Component Analysis，PCA）、本征正交分解（Proper Orthogonal Decomposition，POD）［23］等。其 中，K-L 與Kriging［22］、POD 與PCE［24］、POD 與Kriging［25］等相結合獲得了對高維數據良好的近似能力。然而，對于更多變量高階相關性等復雜問題，線性降維方法通常很難找到有效特征空間。因此非線性降維技術得到了青睞［26］，包括活躍子空間方法（Active Subspace Method，ASM）、核主成分分析 （Kernel Principal Component Analysis，KPCA）、最大方差展開（Maximum Variance Unfolding，MVU）等。其中監督DRSM 技術如ASM 在降維過程中充分考慮了輸入對輸出的影響，相對于非監督式降維技術，通常能獲得更好的表現［27-28］。基于ASM 的DRSM 是一種監督式建模途徑，其利用了輸出對所有輸入方向的梯度信息（?f(X)∈Rn）來選擇和構成活躍子空間。但由于ASM 方法引入了多個高維積分和求導過程，導致了其計算時間花費長和適用性降低的問題［29-30］。盡管目前可利用的非線性降維方法非常多，但為獲得穩定和準確的結果，仍要求大量的訓練樣本（與維數密切相關）和較多的訓練時間，如ASM-Kriging，難以滿足工程需要。

DNN［31］是近些年來得到關注最多的一種建模技術。相對于早期的人工神經網絡（Artificial Neural Network，ANN），它通過更多的隱藏層、反饋方式及學習算法，使得其能近似更加復雜的高維應用。例如，DNN 已應用在代理和分析包含高維不確定性的偏微分方程和雷諾平均湍流模型［32］以及控制方程［33］等。但DNN 需要大量的樣本完成對成千上萬參數的訓練，這對于昂貴的高可信度CFD 分析難以承擔。這幾種方法雖然在一些高維代理問題中得到了應用，但是其均不同程度地面臨適用范圍窄、訓練花費非常高，以及在復雜氣動優化問題中缺乏驗證等難題。

因此，針對以上問題，為建立高效準確并對復雜氣動問題適應性強的工程應用型高維代理模型，本文提出了監督式非線性降維代理建模（Supervised Nonlinear Dimension-Reduction Surrogate Modeling，SN-DRSM）方法。該方法將KPCA 非線性降維模型與代理模型統一訓練，建立了從高維輸入到輸出的準確映射，有效解決了HDMR、DNN 以及非監督式DRSM 等高維代理建模方法訓練時間花費高、適應性差等問題，提高了代理優化方法在高維問題中的優化效率和能力。

1 新型監督式非線性降維代理模型

1.1 基于線性降維高維代理模型的問題

目前代理優化算法使用的代理模型如RBF、Kriging 和GEK 等面臨嚴重的“維度災難”問題［2，34］，其主要體現在：①設計變量個數的增加導致設計空間以指數增加，達到相同樣本密度所需樣本數劇烈增加，使計算花費面臨“維度災難”難題；②隨著變量維數增加，一些代理模型本身出現了病態問題，比如隨著維數增加，Kriging 和GEK等相關矩陣規模過大，訓練時間變長，同時相關矩陣條件數變大，預測精度下降。因此，代理優化方法遭遇“維度災難”的挑戰已經愈加嚴重，目前針對高維代理優化問題已經開展的研究卻非常少。

基于PCA 的非監督DRSM 方法的主要問題包括2 個方面：①線性降維方法對輸入數據建立了最優的線性空間近似，在一些線性或近似線性數據分布中得到了廣泛應用，一般復雜輸入分布均是線性不可分的；②非監督式降維途徑將降維與代理建模分開執行，降維過程沒有考慮輸入對輸出的最有效影響，并且主成分的選擇也是人為截斷，不合適的主成分使代理模型對原始輸入數據預測精度下降。無監督式降維代理模型包括2 個步驟，首先利用降維方法對輸入數據降維；然后在降維后的低維空間修建代理模型。無監督目標包括輸入重構誤差最小、樣本方差最大化、映射前后距離不變等。其具體步驟為先利用PCA 方法建立輸入空間的高維數據到特征空間的線性映射，通過保留特征空間中有效維度（即主成分）以達到降維目的，然后利用經典代理模型建立所有樣本的主成分到輸出的代理，具體過程如圖1 所示。因為該過程僅關注輸入空間數據的特征，找到能近似原始空間的最優特征空間，而忽略了降維有效空間對輸入-輸出復雜映射關系的影響，也影響了它與代理模型結合后的表現。

圖1 基于非監督式PCA 降維的代理建模過程Fig.1 Flowchart of unsupervised PCA-based surrogate model

以跨聲速RAE2822 翼型流場為例，測試狀態計算狀態如表1 所示。測試使用類函數/形函數（Class/Shape Transformation，CST）方法對翼型參數化及拉丁超立方抽樣（Latin Hypercube Sampling，LHS）并在設計范圍內抽取滿足要求的訓練和驗證樣本。結果表明，隨著變量維度增加，使用相同設計變量倍數的樣本數，Kriging 模型預測相對均方根誤差（Relative Root Mean Square Error，RRMSE）增加劇烈。對于較高維變量（n≥32），Kriging 模型預測誤差RRMSE 隨著樣本數（線性）的增加，下降非常緩慢，預測精度難以滿足高效地優化設計需求。

表1 RAE2822 翼型測試狀態Table 1 Test conditions of RAE2822 airfoil

而對于基于PCA 的非監督式降維代理模型（結合Kriging 代理模型，即PCA-Kriging）的預測精度如圖2 所示，PCA 降維均按照99%的能量準則進行截斷。結果表明：相對于原始Kriging 代理模型，PCA-Kriging 模型在較低變量維度時預測誤差明顯更大；隨著輸入變量維數增加，PCAKriging 代理模型預測誤差增加緩慢，但其預測精度并沒有比原始Kriging 模型顯著提升。在高維變量（n≥48）時，隨著樣本數成倍增加，預測誤差有所減少，直至32 倍（32n）設計變量時，預測誤差雖然較原始Kriging 模型有所減少，但預測平均誤差仍然較高（遠高于5counts 的平均阻力系數預測誤差），導致其仍然難以滿足高效的代理優化 設 計 需 求。為 此，一 些 研 究 者［11，22］首 先 基 于Kriging 代理模型優化高維變量問題，在獲得充足數量接近優化解的樣本后（數量遠大于32n），然后使用PCA 進行降維分析后建立DRSM 模型，最后基于DRSM 模型完成優化設計直至收斂。這個過程雖然相對于直接基于Kriging 模型優化提高了優化效率，但其要求充足的數量優化樣本，這對于使用昂貴的高可信度CFD 分析依然造成了巨大的計算花費難題。

圖2 Kriging 和非監督式PCA-Kriging 預測的相對均方根誤差（RRMSE）和相關系數（Corrcoef）隨設計變量維數變化曲線Fig.2 RRMSE and Corrcoef predicted by Kriging and unsupervised PCA-Kriging varying with number of design variables

1.2 新型基于KPCA 的監督式非線性降維代理模型

無監督降維途徑沒有考慮輸入-輸出間復雜的映射關系及輸出對輸入的影響等。因此其降維空間也許并不是最有利于輸出預測的，甚至可能損失重要的信息，導致其很難適用于復雜工程問題代理的建模。針對這個問題，為建立有效的高維代理優化方法，本文引入基于KPCA 的非線性空間降維方法。KPCA 方法建立了輸入空間到特征空間的非線性映射，可滿足廣泛工程問題中數據降維的需求。KPCA 通過選擇合適的核函數和參數即可獲得對輸入空間的有效近似。然而，核函數、核函數參數及特征空間有效維數的選擇對DRSM 的精度影響非常明顯。因此，本文提出了監督式非線性降維代理建模（SNDRSM）方法，即在代理模型訓練中，引入對非線性降維模型KPCA 的優化，使降維后的有效空間輸入提供給代理模型最有用的信息，以獲得最優的DRSM 近似準確率。

1.2.1 基于KPCA 的非線性降維方法

PCA 方法能對原始數據進行線性變換，投影到正交坐標軸系中，以找到方差最大的特征方向，減小數據分析維度。在實際問題中，大多數數據線性不可分（相關性大于二階），即難以通過線性變換找到問題的有效維度。

KPCA 相對于PCA 引入了非線性核函數，將當前空間不可分數據映射到更高維度的特征空間至線性可分，以進行降維。KPCA 降維思路如圖3 所示。對原始空間D中N個n維向量樣本即S={X(i)∈Rn，i=1，2，…，N}， 其 中X(i)=代表第i個樣本，通過非線性核函數?(X)映射到更高維特征空間H中。然后在H中計算協方差矩陣，以及協方差矩陣的特征值和相應的特征向量，得到降維后的有效空間。

圖3 基于KPCA 非線性降維的思路圖Fig.3 Illustration of KPCA-based nonlinear dimensionreduction method

特征空間中的所有運算僅涉及內積操作，通過定義合適的核函數，可在原始空間中完成相應運算，即

核函數K應該滿足Mercer 條件，對原始數據進行非線性映射后，可得到特征空間中的協方差矩陣，即

式中：?與ψj均為列向量。

假設n維輸入樣本X()i∈Rn映射到高維特征空 間 后 的 維 數 為nh，則 有?∈Rnh，ψi∈Rnh，CH∈Rnh×nh，其中

然后有

其中：

假設vi是/N的特征向量（列向量模式），λi是對應vi的特征值，則有

式中：vi和λi可通過對矩陣/N進行特征分解得到，即

可求得N個正交的特征向量和非負特征值，N個特征值按從大到小排列為λ1≥λ2≥…λN≥0，對應于N個特征向量v1，v2，…，vN。然后利用求得的特征值和特征向量有

將式（2）代入式（8）可得：

得到協方差矩陣CH的特征向量為，以及相應的特征值為λi。得到特征空間主成分后，對每個非線性映射?(X(j))投影即可得到特征坐標，即可得到第j個樣本X(j)降維后對應變量的第i個分量y(j)i。

式中：是矩陣的列向量。通過選取前q個最大特征值對應的特征向量（主成分），即可得到降維后的數據為

1.2.2 SN-DRSM 框架

使用KPCA 和Kriging 代理模型一起訓練，構成基于KPCA 的監督式非線性降維代理模型，其具體流程如圖4 所示。KPCA-Kriging 方法的核心部分是執行一個嵌套優化循環：外循環優化KPCA 降維模型參數（如KPCA 核函數參數及特征空間有效維數等），使預測誤差盡可能小，直至收斂，其中在每次降維分析后向Kriging 模型輸入降維后的數據；內循環接收KPCA 降維后的樣本輸入優化代理模型參數，直至輸出最優代理模型并輸出預測值。預測誤差使用留一交叉驗證誤 差（Leave-One-Out Cross-Validation error，LOOCV）εLOO評定標準，以減少高可信度樣本CFD 分析時間。在內循環優化中，代理模型的訓練通常包含許多待優化參數，包括KPCA 核函數參數、特征空間有效維數及Kriging 模型超參數等。最常見的Kriging模型超參數訓練方法包括最小化LOOCV 誤差和最大化似然估計等。本文主要使用最大似然估計方法訓練超參數，以獲得對一般問題的適應性［35］。如此的嵌套循環使代理模型自動選擇和訓練降維模型，直至尋找到滿足精度要求的最佳降維模型。同時在每個子優化中，優化參數變少，優化時間變短。KPCA-Kriging 也是完全非嵌入式的，在不同的問題中可方便使用。

圖4 基于KPCA 的監督式非線性降維代理建模流程Fig.4 Flowchart of KPCA-based supervised nonlinear dimension-reduction surrogate modeling

仍使用跨聲速RAE2822 翼型流場（表1 所示的測試狀態）檢測KPCA-Kriging 代理模型的實際表現，結果如圖5 所示。其中KPCA-Kriging模型使用Gaussian 核函數。結果表明，KPCAKriging 模型在較高維度時（n≥32），其預測精度比原始Kriging 模型顯著提高，并隨著樣本數成倍增加（≥16n）時，其平均阻力系數預測誤差小于10 counts 以及相關系數大于95%，可為高維變量氣動優化設計提供比原始Kriging 模型更為準確的預測。從結果可知，監督式方法在構建降維代理模型過程中充分考慮了輸入對輸出的影響，使輸入空間降維模型的訓練與代理模型輸出空間緊密聯系，從而使降維模型提供給代理模型最有用的信息，以獲得最優的降維代理模型近似準確率。而非監督式方法將降維模型的訓練和代理模型的訓練分別單獨執行，缺乏對復雜問題的適應性，導致預測精度也難以保證。

圖5 Kriging 和監督式KPCA-Kriging 預測的相對均方根誤差（RRMSE）和相關系數（Corrcoef）隨設計變量維數變化曲線Fig.5 RRMSE and Corrcoef predicted by Kriging and supervised KPCA-Kriging varying with number of design variables

1.3 高維氣動算例測試

本節使用美國國家航空航天局（National Aeronautics Space and Administration，NASA）CRM（Common Research Model）機翼［36］作為測試算例以驗證所提出的代理建模方法的全局預測精度。CRM 機翼測試狀態Ma∞=0.85，CL=0.50，Re=5.0×106，計算網格量為193 382，湍流模型使用SST。設計變量數為n=40～248，及相應的自由變形（Free Form Deformation，FFD）框如圖6 所示。每套控制框的變形包括控制點沿z方向平移及各個剖面繞前緣點的扭轉變形。在代理模型構建過程中，使用拉丁超立方抽樣得到數量為Nr(n)=4n的訓練樣本以及Nv=1 000 的驗證樣本。使用Kriging、無監督PCAKriging 以及監督式KPCA-Kriging 這3 種代理建模方法來驗證阻力系數預測精度。其中KPCAKriging 使 用Gaussian 核 函 數 訓 練，而PCA 降 維均按照99%能量準則截斷。

圖6 不同設計變量數的FFD 框Fig.6 FFD frames corresponding to different numbers of design variables

降維后得到的有效維數與輸入維數的對比如圖7 所示，可以看出在高維時KPCA-Kriging訓練后得到的有效維數小于PCA-Kriging 模型選擇的有效維數以及遠小于實際輸入變量的維數。本文針對118 維設計變量輸入，KPCAKriging 訓練后有效維數僅為24 維，而此時PCAKriging 訓練得到的有效維數為51。各個設計維度下，3 種代理模型得到RRMSE 誤差對比曲線如圖8 所示，結果表明在高維各個輸入變量維度下KPCA-Kriging 方法交叉驗證誤差明顯小于Kriging 方法和PCA-Kriging 得到的預測誤差。而PCA-Kriging 模型相對于Kriging 模型并沒有顯著優勢，導致其在優化設計中難以被使用。尤其Kriging 方法在設計變量維數大于118 后預測精度就急劇變差；而KPCA-Kriging 方法此時預測精度仍然保持在一個可接受的范圍內，并明顯小于Kriging 方法和PCA-Kriging 的預測誤差。比如在118 維時，Kriging 方法得到的RRMSE 為0.030 223，PCA-Kriging 模型得到的RRMSE 為0.028 9，而KPCA-Kriging 方 法 預 測RRMSE 為0.016 520，其精度較前兩者分別提升了45.34%和42.84%。隨著變量維數的增加，KPCAKriging 模型優勢更加顯著。原因之一為要維持相同樣本密度，所需樣本數指數增加，而KPCAKriging 降維后有效維數遠少于輸入變量維數，相同樣本數下樣本密度顯著增加，有效提升了預測精度。

圖7 KPCA-Kriging 訓練后選擇的有效維數隨輸入維數變化曲線Fig.7 Effective dimensions of feature space selected by tuned KPCA-Kriging model varying with input dimensions

圖8 KPCA-Kriging 與Kriging 模型驗證誤差隨設計維數變化曲線Fig.8 Validation errors of KPCA-Kriging model and Kriging model varying with number of design variables

118 維設計變量輸入時，KPCA-Kriging 方法訓練得到的核矩陣特征值隨特征空間中有效維數衰減曲線如圖9 所示，可以看出在特征空間（更高維的空間）中，幾乎所有能量（特征值的貢獻）集中在前面很少的主特征方向里，前24 維占據了所有能量的90% 更多。KPCA-Kriging 預測誤差隨截斷的有效維數變化曲線如圖10 所示。隨著選擇的有效維數（主特征方向）逐漸增加，預測誤差快速下降，直至到24 維左右，預測誤差開始小幅波動，并不再顯著減少。而繼續增加截斷維數到50 維以后時，預測誤差開始略有增加。因此，通過監督式式KPCA-Kriging 訓練后，代理模型自適應地選擇了最優模型參數和最佳有效維數，使得選擇的特征空間中有效維數（主特征方向）遠小于實際輸入變量維數，進而有效緩解了代理模型的維度災難問題。

圖9 核矩陣特征值隨有效維數衰減曲線（118 維設計變量）Fig.9 Eigenvalues of kernel matrix varying with effective dimensions （118 design variables）

圖10 KPCA-Kriging 模型驗證誤差隨選擇的有效維數變化曲線（118 維設計變量）Fig.10 Validation errors of KPCA-Kriging model varying with selected numbers of effective dimensions （118 design variables）

2 基于代理模型的自適應加點優化框架

基于代理模型的優化設計方法首先需要建立可用的代理模型，以便在代理空間內搜索。隨著設計變量維數的增加，滿足代理模型全局精度要求的樣本數也急劇升高，這使得優化代價大大增加。并且當使用的樣本數較少時，難以保證代理模型的精度，進而導致優化解容易陷入局部最優。因此，自適應加點方法通過每次利用已有樣本點、已有代理模型信息和已有優化結果信息選擇新的樣本點增加到樣本庫，這使得其需要平衡兩方面內容：一是增加樣本點到未探索到但很有可能包含最優解的設計空間，即全局探索能力（Exploration）；二是選擇點在已經探索到并感興趣的區域（如可能包含最優值的區域，或者精度很差的高非線性區域等），即增加局部空間代理預測精度（Exploitation）。

圖11給出了自適應空間擴展的代理優化流程圖，每次優化結束后，通過自適應加點技術選擇加點位置，評估新樣本點并添加到樣本庫，構建代理模型開始再次優化，這個過程為自適應加點給代理模型和優化搜索算法提供信息，直到收斂為止。目前針對代理優化的自適應加點方法包括最小化代理預測（Minimizing Surrogate Prediction，MSP），最大化期望改進（Expected Improvement，EI）策略，置信下界（Lower Confidence Bounds，LCB）策略，以及多種加點策略混合的加點方式等。筆者前期工作［1］給出了單一加點和混合加點的詳細對比，MSP+EI 和MSP+LCB 的混合加點策略相對于單一加點方法能獲得更好的優化結果及更好的魯棒性，如圖12 所示。而單一的加點技術優缺點非常明顯，如MSP具有良好的收斂特性但容易陷入局部最優，而EI具有全局優化能力但隨著加點過程到最后其收斂性變差。而組合加點準則通常兼具了兩種加點策略的優點并提高了優化搜索的魯棒性。本文主要使用EI+MSP 混合加點策略進行優化設計，即一步EI 加點后一步MSP 加點，如此重復直至最大EI 值小于閥值后，使用MSP 加點直至收斂。

圖11 基于自適應空間擴展的代理優化流程圖Fig.11 Flowchart of surrogate-based adaptive spaceextension optimization method

圖12 不同加點方法優化結果對比Fig.12 Comparison of converged curves of different infilling strategies in surrogate-based optimization

3 復雜氣動外形優化設計應用

3.1 跨聲速無黏NACA0012 翼型

該氣動算例定義為跨聲速無黏狀態下進行單點減阻優化，設計狀態為Ma=0.85，α=0o。優化模型為

式（12）要求優化翼型在各個坐標點處厚度不小于初始NACA0012 翼型，該狀態下，翼型的阻力來源于翼型的激波阻力，優化翼型為對稱分布，采用半模外形和半模計算網格進行氣動分析和優化。網格如圖13 所示［37］，初始翼型阻力為471 counts。采用CST 參數化方法變形，設計變量數為24，初始設計空間范圍如圖14 所示。采用圖11 所示的基于自適應空間擴展（Adaptive Design space Extension，ADE）的代理優化框架進行優化設計，并分別使用KPCA-Kriging 和Kriging 代理模型，設計結果分別記為ADE-KPCAKriging 和ADE-Kriging。采用給定的設計空間和使用Kriging 代理模型進行對比優化設計，記為Kriging。3 種方法均使用相同的初始樣本，樣本數為Nini=4d=96，并由LHS 抽樣方法產生。同時均利用EI+MSP 加點優化選取了200 個新樣本點添加到樣本庫中，優化歷程對比如圖15 所示，可以看出基于自適應設計空間擴展的2 種優化方法阻力保持了持續的下降，直到150 步以后接近收斂，而基于固定設計空間的方法在50 次迭代后就已經收斂，阻力值不再下降。3 種優化方法得到的阻力系數對比如表2 所示，基于自適應設計空間擴展方法獲得了更低的阻力值，其中基于KPCA-Kriging 模型優化方法獲得了最低的阻力系數42 counts，顯著低于目前國內基于Kriging方法和基于多水平分層Kriging 方法得到的收斂阻 力 值（259.7 counts）［38］。結 果 顯 示3 種 方 法 在總花費相同的情況下（均使用246 個樣本），本文發展的方法得到了更低的阻力值。外形對比如圖16 所示，3 種方法得到的外形后緣有顯著的差異。基于自適應空間擴展的優化方法在每次子迭代結束后，自動判斷每個設計變量是否到達邊界，如果到達邊界則自動對設計變量邊界進行適當拓展后再迭代優化，否則繼續優化。在本算例優化過程中，基于自適應空間擴展的優化方法能自動選擇并拓展翼型后緣的設計變量邊界，避免了因邊界限制導致收斂過程變慢的問題，從而顯著提高了優化效率。

圖13 翼型計算網格［37］Fig.13 Computational grid of airfoil ［37］

圖14 設計空間展示Fig.14 Legend of design space

圖15 阻力系數收斂曲線對比Fig.15 Converged curves of drag coefficients

表2 優化翼型與設計翼型阻力系數對比（全模）Table 2 Comparison of drag coefficients of baseline and optimized airfoils （whole model）

圖16 優化設計翼型外形對比Fig.16 Comparison of baseline and optimized airfoils

圖18 NACA0012 壓力系數云圖Fig.18 Contour of pressure coefficients of NACA0012 airfoil

優化翼型與NACA0012 翼型壓力分布對比如圖17～圖19 所示，結果顯示，優化翼型前緣壓縮和后緣恢復，后緣激波變弱，壓力分布更趨近“對稱”。從上述優化過程可以看出，最小阻力對應的外形應該是趨同的，極有可能為單極值優化問題。KPCA-Kriging 通過監督式訓練過程能更好地找到對阻力影響最敏感的一些變量或方向，在加入越多的訓練樣本點后能準確地預測阻力值，提高了代理優化效率。

圖17 優化翼型上表面壓力分布對比Fig.17 Comparison of pressure coefficients distributions on upper surfaces of optimized airfoils

圖19 使用ADE-KPCA-Kriging 得到的壓力系數云圖Fig.19 Contour of pressure coefficients using ADEKPCA-Kriging method

3.2 跨聲速CRM 機翼

CRM 機 翼（NASA Common Research Model）減阻優化算例為ADODG 定義的第4 個優化問題。 CRM 機翼外形參數為Cref=1.0，Sref=191.845，Xref=33.678，其 外 形 和 計算網格如圖20 所示，計算網格量為2 050 560。設計狀態為Ma∞=0.85，CL=0.50，Re=5.0×106。本算例使用SST 湍流模型進行計算，初始機翼阻力為234.1 counts。設計問題定義為在每個翼剖面最大厚度、力矩系數以及升力系數不減小的情況下最小化阻力系數，其優化數學模型為

圖20 CRM 機翼計算空間網格Fig.20 Computational grid of CRM wing model

采用FFD方法對機翼參數化，如圖21所示，選擇7個剖面對機翼變形，每個剖面定義了16個控制點，通過控制每個控制點沿z方向移動從而產生新外形，機翼剖面不扭轉。為對比基于2 種代理模型方法優化效率的差異，均采用相同的初始樣本，并優化迭代100 次，采用MSP+EI 混合加點策略。初始樣本均通過LHS 抽樣產生，樣本數為Nini=3×112=336。2 種方法迭代加點100 次，在100 次迭代后2 種方法優化得到的機翼計算結果如表3 所示，結果表明，基于KPCA-Kriging 代理優化得到的機翼阻力為205.2 counts，較初始機翼降低了28.9 counts。而基于Kriging 代理優化的機翼阻力為219.4 counts，較初始機翼僅僅降低了14.7 counts。并且基于Kriging 代理模型的優化方法收斂很慢，直到100 次迭代加點得到的優化機翼阻力也高于基于KPCA-Kriging 代理建模方法的結果。在總花費相同的情況下（均使用436 個樣本），本文提出的方法得到了更好的結果，表明了其效率相對于傳統方法有顯著的改進。

圖21 CRM 機 翼FFD 框Fig.21 FFD frame of CRM wing model

表3 初始機翼與優化機翼氣動力系數對比Table 3 Comparison of aerodynamic force coefficients of baseline and optimized wings

2 種方法優化得到機翼上表面壓力云圖和初始外形壓力云圖對比如圖22 和圖23 所示。由圖可知，初始CRM 機翼上表面存在較強激波，基于Kriging 代理優化得到的機翼上表面激波強度有所減弱，但仍然較強；而基于KPCA-Kriging 代理優化得到的機翼上表面僅僅存在非常弱的激波。2 種方法收斂曲線對比如圖24 所示，結果顯示基于KPCA-Kriging 代理模型方法的優化過程較基于Kriging 代理優化過程收斂顯著更快，最終收斂阻力也顯著更低。進一步可以發現，在80 次加點后繼續迭代，阻力值下降非常緩慢，隨機搜索算法效率變低。基于2 種代理模型優化過程中加點位置處代理模型預測誤差對比曲線如圖25 所示。顯然KPCA-Kriging 代理模型預測誤差始終小于Kriging 代理模型預測誤差，并且隨著加點次數接近100 次時，Kriging 模型預測誤差仍然在20 counts 左右波動，而KPCA-Kriging 模型預測誤差下降到5 counts 以內。

圖22 CRM 機翼與優化機翼上表面壓力云圖對比Fig.22 Comparison of pressure coefficients contours of CRM and optimized wings

圖23 2 種方法優化機翼上表面壓力云圖對比Fig.23 Comparison of pressure coefficients contours of optimized wings using different optimization methods

圖24 基于兩種代理模型的阻力系數優化收斂歷程對比Fig.24 Comparison of converged curves of drag coefficient of two optimization methods

圖25 優化過程中阻力系數預測誤差對比Fig.25 Comparison of prediction errors of drag coefficient during two optimization processes

機翼上要對比的8 個翼剖面位置如圖26 所示，不同方法優化的機翼與初始機翼在這8 個剖面處外形對比如圖27 所示，2 種方法優化得到的機翼與原始機翼在關鍵剖面位置處的壓力分布與翼型對比如圖28 所示。從結果可以看出各個剖面位置處基于KPCA-Kriging 代理優化得到的翼型激波強度明顯弱于基于Kriging 方法得到的剖面翼型激波強度和原始機翼剖面翼型激波強度。比如D、F和H處的壓力分布對比顯示，基于KPCA-Kriging 代理優化得到的翼型激波強度非常弱，幾乎被抹平。設計結果體現了基于KPCA-Kriging 代理優化方法在高維問題中相對于原始基于Kriging 代理優化方法的巨大優勢和效率的顯著改進。本文基于新型高維代理模型優化結果也明顯好于文獻［39］給出的基于梯度加強Kriging 代理模型優化結果。因此，針對此典型算例，基于KPCA-Kriging 代理模型優化方法相對于基于傳統Kriging 代理模型優化方法效率顯著提升，其優化得到的機翼激波更弱，阻力系數也明顯更低。

圖26 機翼8 個剖面位置示意圖Fig.26 Sketch map of eight constraint sections of wings

圖27 設計機翼與CRM 機翼8 個剖面翼型對比Fig.27 Comparison of eight section airfoils of CRM and optimized wings

圖28 優化機翼與原始CRM 機翼各個剖面位置翼型與壓力分布對比Fig.28 Comparison of airfoil and pressure coefficients of section airfoils of CRM and optimized wings

此外，基于KPCA-Kriging 代理模型優化方法得到的機翼壓力分布圖也顯示機翼上表面的局部激波沒有完全消除，全機的氣動特性并沒有達到全局最優。這是因為該算例氣動優化空間可能呈現一個大的“碗狀”，最優阻力值對應的外形位于“碗底”某位置，從初始外形開始搜索時，剛開始阻力值下降迅速，然后進入碗底附近阻力值就下降的非常慢，這使得隨機搜索算法每次子迭代阻力值只能下降一點（甚至小于0.1 counts），可能導致需要數以千次迭代和添加高可信度樣本點才能搜索到全局最優值，這個計算花費對于昂貴的機翼優化是難以承擔的［40］。在相同的總計算花費下，發展的方法相對于傳統的方法能找到更低阻力值的機翼外形；而傳統方法確在更高阻力值下接近收斂，使之更難以找到最優氣動外形。

4 結 論

本文提出了基于KPCA 的監督式非線性降維代理建模方法。該方法將降維模型與代理建模統一訓練，建立了從高維輸入到輸出的準確映射，解決了高維代理模型在復雜工程應用中精度差和穩定性低的難題，并發展了適用于復雜氣動外形設計的高維全局優化平臺，通過AIAA 氣動優化小組發布的典型算例驗證。本文代理優化平臺對2 個典型算例均能取得良好的優化效果，滿足工程使用要求。

通過對比發現基于Kriging 代理優化方法隨著設計變量維度從24 維增加到112 維時，收斂速率明顯降低，要求更多的新樣本點和迭代次數才能獲得收斂解，甚至遠遠更多的迭代次數才可能獲得和基于KPCA-Kriging 代理優化接近的結果。而基于KPCA-Kriging 代理優化方法對24維翼型優化算例能獲得相對于基于Kriging 代理優化方法更快的收斂速率。在對112 維的CRM機翼優化中，基于KPCA-Kriging 的代理優化相對于Kriging 代理優化效率提升顯著，在相同的加點次數下能獲得明顯更好的優化解。當加點到100 次時，基于KPCA-Kriging 代理優化方法幾乎收斂，而Kriging 優化方法收斂速度非常慢，需要更多的樣本才可能收斂到相近解。這是因為KPCA-Kriging 降維后有效維數為24 維，在同樣的樣本數下，降維空間樣本密度遠高于初始設計空間，有利于提高代理預測精度。同時通過監督式訓練過程，降維后的有效空間選擇是最有利于輸出預測的，提供了對輸出最有用的信息，提高了代理優化的收斂效率。因此，本文提出的基于監督式非線性降維代理建模方法，有效改進了代理優化方法在100 維左右氣動優化問題中的優化效率，拓展了代理優化方法在高維問題中的適應能力，有效緩解了“維度災難”難題對代理優化方法的限制。