楊國濤，岳振江，2，，劉莉，2

1．北京理工大學 宇航學院，北京 100081

2．北京理工大學 飛行器動力學與控制教育部重點實驗室，北京 100081

隨著航空航天技術(shù)的發(fā)展，高超聲速飛行器由于速度快、突防能力強、覆蓋范圍廣、隱身效果顯著等優(yōu)點，具有廣泛的應(yīng)用前景[1］。由于高超聲速飛行器在大氣層內(nèi)高速飛行，氣動熱防護是需要妥善解決的關(guān)鍵問題之一[2］。傳統(tǒng)氣動熱獲取主要包括試驗研究、工程計算、數(shù)值模擬3 種方法[3］。飛行試驗和地面風洞試驗能夠獲取飛行器氣動熱的真實狀況，但長周期、高成本限制了其在設(shè)計階段初期的應(yīng)用[3］。氣動熱工程計算基于試驗數(shù)據(jù)或邊界層方程相似解，通過合理假設(shè)，采用經(jīng)驗公式、半經(jīng)驗公式獲得，具有快速高效的優(yōu)點，但只能處理簡單外形的飛行器[4］。數(shù)值模擬是利用計算流體力學（Computational Fluid Dynamic，CFD）方法求解Navier-Strokes 方程及其近似，對飛行器外形適應(yīng)性強，能充分考慮氣流黏性、真實氣體效應(yīng)等，但較低的計算效率在設(shè)計迭代過程中應(yīng)用受限[3］。針對上述傳統(tǒng)方法存在的局限性，為有效解決高耗時高成本氣動熱獲取的工程問題，出現(xiàn)了利用相對簡單的數(shù)學模型來表達復雜全階空氣動力學模型輸入輸出關(guān)系的方法[5］，以達到快速預示氣動熱的目的，并逐漸成為近年來氣動研究的熱點[6］。陳海昕等[7］總結(jié)了機器學習在氣動優(yōu)化設(shè)計中的應(yīng)用，并進一步探討了利用深度代理模型和深度降階模型進行氣動數(shù)據(jù)的歸納分析和流場重建。Yan 等[8］基于POD 和修正的切比雪夫多項式提出一種POD-Chebyshev 方法，應(yīng)用于飛行器控制面將平均誤差降到了2%以內(nèi)。張智超等[9］基于徑向基神經(jīng)網(wǎng)絡(luò)提出一種氣動熱預測代理模型方法，應(yīng)用在NASA 火星實驗室的橢圓鈍化高超聲速飛行器中，關(guān)鍵區(qū)域預示誤差小于10%。

氣動熱數(shù)學近似模型主要通過采樣獲取物理系統(tǒng)的主要特征，因此基于采樣技術(shù)獲取的樣本質(zhì)量很大程度決定樣本所包含的系統(tǒng)信息，進而影響模型的預示精度?，F(xiàn)有的采樣方法包括一 次 性 采 樣（One-Shot）和 序 列 采 樣[6，10］。One-Shot 一次性獲取具有良好填充均勻性和投影均勻性的樣本，建模次數(shù)少、效率高，常用的方法有全析因設(shè)計、正交設(shè)計、均勻設(shè)計、拉丁超立方設(shè)計（Latin Hypercube Design， LHD）和最優(yōu)LHD（OLHD）等[6，10-11］，但是存在欠采樣或過采樣的局限[12］。序列采樣利用少量初始樣本構(gòu)建近似模型，并通過迭代增加樣本更新近似模型，以達到有效控制采樣規(guī)模和兼顧模型精度的目的[12］。自適應(yīng)采樣是一類重要的序列采樣方法，在利用已有樣本空間信息的同時，還充分考慮系統(tǒng)的物理特征在具有較大誤差或高度非線性等局部欠采樣區(qū)域重點新增樣本，構(gòu)建更加準確高效的近似模型，并達到節(jié)約計算成本的目的，相較于一次性采樣，自適應(yīng)采樣具有更加廣泛的應(yīng)用前景[13］。

根據(jù)預示模型的用途，自適應(yīng)采樣方法主要包括兩類，一類是面向優(yōu)化以提高模型在其關(guān)注的局部區(qū)域的預示精度[14-16］；另一類是面向全局近似以提高整體預示精度[17］。針對全局近似任務(wù)，自適應(yīng)采樣通過順序迭代獲取新增樣本信息，以盡可能小的采樣規(guī)模獲取整個參數(shù)空間內(nèi)充分的系統(tǒng)信息，支持響應(yīng)的準確預示。全局探索與局部開發(fā)的平衡對全局近似效果至關(guān)重要，有效的平衡不僅可以在已識別的區(qū)域中進行充分的局部開發(fā)，還可以及時地提供全局探索以避免錯過一些未檢測到的區(qū)域[18］。Xiong 等[19］利用基于非平穩(wěn)協(xié)方差的 Kriging 模型的預測方差輔助采樣以減小模型的不確定性，從而提高預測質(zhì)量。Yao 等[20］提出一種基于梯度的徑向基神經(jīng)網(wǎng)絡(luò)方法，利用現(xiàn)有模型的梯度信息依次擴展樣本集并細化模型，從而有效提高模型近似精度。Farhang-Mehr 等[21］提出了最大熵自適應(yīng)采樣的方法，旨在將樣本放置在響應(yīng)的不規(guī)則區(qū)域。Turner 等[22］采用模擬退火來平衡局部開發(fā)和全局探索，以達到多準則自適應(yīng)采樣的目的。Li 等[23］提出一種基于預測誤差的Kriging 近似自適應(yīng)方法，以改善響應(yīng)中近似不佳的部分。Aute 等[24］提出了基于交叉驗證誤差的序列采樣方法，根據(jù)預測誤差在高度非線性區(qū)域采集更多的點。Chen 等[25］提出了一種基于模糊聚類的自適應(yīng)加點策略（Adding-Point Strategy based on Fuzzy Clustering， APSFC），通過對具有較大誤差的樣本集進行聚類，并在聚類內(nèi)部產(chǎn)生新增樣本，一定程度提高了構(gòu)建高超聲速升力面氣動熱預示模型的效率和精度。Crombecq 等[12］提出了一種通用可靠的LOLA-Voronoi 混合序列采樣方法，通過局部線性近似衡量樣本周圍非線性程度和通過Voronoi 圖衡量樣本周圍密度相結(jié)合的方法獲取新增樣本，所提方法針對不同模型類型都能取得很好的效果，但是在估計樣本梯度時需要進行復雜的鄰域計算。Xu 等[26］將LOO（Leave-One-Out）交叉驗證和Voronoi 圖相結(jié)合提出一種通過誤差追蹤逐步提高預示模型準確性的CV-Voronoi 方法，與LOLA-Voronoi 相比，該方法效果更優(yōu)。Kaminsky 等[27］還 結(jié) 合k-fold 交 叉 驗 證 和 批 量 采樣技術(shù)提出了一種改善的KFCV-Voronoi 方法，降低了采樣過程中的計算成本。

為加快高超聲速飛行器氣動熱環(huán)境分析，本文提出了一種基于模糊聚類超球體影響域和當?shù)卣`差評分系數(shù)加權(quán)拒絕域的自適應(yīng)采樣方法（adaptive sampling method based on Fuzzy Clustering HyperSphere influence domain and local error scoring coefficient Weighted refused domain， FCHSW），在此基礎(chǔ)上建立氣動熱全局快速預示模型（Rapid Predicting Model，RPM）。首先，對一次采樣獲取的少量高精度工況構(gòu)建氣動熱RPM；然后，評估模型的誤差分布特征，建立基于模糊聚類的超球影響域和基于局部誤差評分系數(shù)加權(quán)的拒絕域，進而結(jié)合maxmin 準則在綜合確定的重點采樣域新增高精度樣本，提高采樣質(zhì)量，實現(xiàn)RPM 的自適應(yīng)更新，高效提升全局預示精度。

1 氣動熱全局快速預示框架

高超聲速飛行器氣動熱快速準確預示是高超聲速熱防護和防隔熱設(shè)計分析的重要前提。面向高精度高效率氣動熱設(shè)計的需求，本文采用基于數(shù)據(jù)驅(qū)動的氣動熱全局快速預示框架，根據(jù)有限的高精度氣動熱數(shù)值計算結(jié)果構(gòu)造簡單數(shù)學模型來替代復雜CFD 模型的評估，實現(xiàn)參數(shù)空間內(nèi)飛行器氣動熱環(huán)境快速準確預示。

氣動熱全局快速預示框架按順序分為CFD氣動熱計算與RPM 構(gòu)建、氣動熱RPM 測試與更新2 個關(guān)鍵步驟，基本流程如圖1 所示。

圖1 氣動熱RPM 框架示意Fig.1 Schematic of aerothermodynamics RPM framework

1.1 CFD 氣動熱計算與RPM 構(gòu)建

為有效控制飛行器高精度數(shù)值計算成本，需要從眾多流場控制參數(shù)中確定對物理系統(tǒng)有重要影響的最小維度飛行狀態(tài)參數(shù)設(shè)計空間，并利用采樣技術(shù)獲取少量具有良好空間填充性和投影均勻性的初始訓練樣本，以代表整個飛行狀態(tài)空間。

根據(jù)選定的飛行器氣動外形，基于CFD 計算策略獲取所有初始飛行工況的高精度氣動熱數(shù)值計算結(jié)果。然后基于給定狀態(tài)下流場的輸入輸出特性確定數(shù)學近似模型，并采用參數(shù)估計法對數(shù)學模型進行訓練，完成飛行器表面氣動熱流RPM 的構(gòu)造。

1.2 氣動熱RPM 測試與更新

初始構(gòu)建的訓練模型往往精度較低，為了使RPM 能充分實現(xiàn)對大規(guī)模高耗時復雜CFD 系統(tǒng)的近似，根據(jù)訓練樣本的空間分布信息和訓練模型的誤差分布特性，基于自適應(yīng)采樣方法添加新的訓練樣本，并調(diào)用CFD 模型獲取新增樣本響應(yīng)，據(jù)此迭代更新訓練模型以逐步提高其預示精度。模型的誤差分布特性可以通過預示方差、交叉驗證誤差以及梯度等信息來衡量，其中交叉驗證誤差由于易于獲取且不依賴模型，更廣泛的用于自適應(yīng)指導模型更新[13，24］。

為了提供系統(tǒng)的主要特征，要求RPM 與CFD 模型誤差應(yīng)當足夠小，因此獲取額外的測試樣本對每一輪更新的模型進行精度校驗，以驗證模型的準確性，并判斷迭代終止。

由于訓練樣本規(guī)?；蚋呔扔嬎愠杀镜南拗?，為了提升高精度樣本采樣質(zhì)量，進一步改善RPM 的全局預示精度，更好的適用于可能存在的多極值、強非線性問題，自適應(yīng)采樣是指導RPM高效更新的關(guān)鍵技術(shù)。

2 FCHSW 基本方法介紹

針對有限計算成本下提升RPM 全局預示精度的問題，本文提出的FCHSW 方法基本流程如圖2 所示，影響域和拒絕域分別為新增樣本的潛在區(qū)域和拒絕區(qū)域，同時結(jié)合maxmin 準則，綜合權(quán)衡參數(shù)空間的全局探索和局部開發(fā)，實現(xiàn)RPM精度的高效提升。

圖2 FCHSW 基本流程Fig.2 Basic flow of FCHSW

FCHSW 具體步驟如下：

步驟 1確定nv維設(shè)計變量、參數(shù)空間Ω、測試樣本數(shù)nt及方法相關(guān)參數(shù)。方法相關(guān)參數(shù)主要包括：初始訓練樣本數(shù)ns0、每輪新增樣本數(shù)nadd、波 動 系 數(shù)μ以 及 終 止 條 件（預 期 誤 差 閾 值e?est和終止樣本數(shù)閾值nest）。

步驟 2基于OLHD 在參數(shù)空間生成相互獨立的初始訓練樣本集x和測試樣本集xt，計算真實響應(yīng)y和yt，并初始化樣本數(shù)據(jù)庫。

步驟 3利用初始訓練樣本集構(gòu)建RPM，利用測試樣本集進行模型精度校驗。判斷是否滿足終止條件，即當預期誤差e?＜e?est或訓練樣本數(shù)ns＞nest，則迭代終止，否則進行步驟4。本文將平均相對誤差作為模型精度評估的停機準則：

其中：上標U 和L 分別表示真實響應(yīng)yt的最大值和最小值；下標Nt 表示響應(yīng)在區(qū)間[1，2]的歸一化值，表示為

步驟 4利用k-fold 或LOO 獲取所有樣本的交叉驗證誤差CV，挑選出大于平均誤差的樣本做模糊聚類，進而得到每個聚類的影響域DEi(1 ≤i≤nc)，其 中nc為聚 類個數(shù)。

步驟 5利用超球體對參數(shù)空間進行分割，獲取每個樣本的拒絕域DRi(1 ≤i≤ns)。

步驟 6根據(jù)獲取的影響域和拒絕域，新的設(shè)計域為初始參數(shù)空間對拒絕域的補集與影響域的交集?？紤]到新的設(shè)計域為不規(guī)則區(qū)域，采用Monte Carlo 仿真在整個參數(shù)空間生成隨機分布的若干樣本xrand，并利用接受-拒絕采樣將設(shè)計域Dc內(nèi)的隨機樣本作為候選樣本xcandidate；進而根據(jù)maxmin 準則[28］從xcandidate中選出每輪迭代的新增樣本xadd并計算響應(yīng)yadd；最后，更新樣本數(shù)據(jù)庫并返回步驟3 以實現(xiàn)RPM 的更新。

隨機樣本個數(shù)定義為N(xrand)=ns0×nv×ω，ω為控制參數(shù)，為保證方法有效性的同時盡可能降低計算成本，參考文獻[26］的取值，本文取ω=100。

圖3 展示了二維情況下FCHSW 的一輪迭代加點示意，實心圓為訓練樣本，圓半徑越大表示誤差越大，將具有較大誤差的樣本集分為2 個聚類，三角形為聚類中心，實線圓表示樣本的拒絕域邊界，填充部分為新的設(shè)計域，五角星代表新增樣本。

圖3 FCHSW 加點過程示意Fig.3 Adding-point process of FCHSW

2.1 模糊聚類

基于模糊聚類[29］按照相似度對數(shù)據(jù)集進行分類，建立樣本與類的不確定性描述，以實現(xiàn)對參數(shù)空間的有效劃分和縮減。模糊c-均值聚類（Fuzzy C-Means clustering method，F(xiàn)CM）將問題歸結(jié)為一個非線性規(guī)劃問題，通過極小化各數(shù)據(jù)點到聚類中心的加權(quán)平均和得到每個數(shù)據(jù)點對所有類中心的隸屬度，實現(xiàn)對樣本自動分類[29］。

考 慮 樣 本 集X=[xi](i=1，2，…，ns)，將樣本集中元素根據(jù)模糊聚類的方法劃分為nc(2 ≤nc≤ns)個模糊子集，每個數(shù)據(jù)點對所有類中心的整體差異可表示為

式中：U=[uij] (i=1，2，…，nc；j=1，2，…，ns)為隸屬矩陣，uij表示第j個樣本xj在第i個聚類的隸屬度，滿足，表達式為

其中：d為影響隸屬度矩陣模糊化程度的指數(shù)權(quán)重（通常取d=2）。求解式（3）中的極小值問題，即可得到聚類中心和隸屬度矩陣U*。

通過最大輪廓系數(shù)法來預先確定FCM 的最佳聚類個數(shù)，以實現(xiàn)方法的無監(jiān)督學習。輪廓系數(shù)定義如下[30］：

式中：a(i)為第i個樣本xi與其所在簇內(nèi)其他樣本的平均距離，代表簇內(nèi)凝聚度；b(i)為樣本xi與其他每個簇樣本平均距離的最小值，代表簇間分離度。輪廓系數(shù)滿足-1 ≤S≤1，越接近1 說明聚類效果越好，反之越接近-1 聚類效果越差。

2.2 影響域

影響域是由模糊聚類所屬較大誤差樣本集拓展得到的一個參數(shù)空間，定義為超球體?？紤]第i 個影響域DEi(1 ≤i ≤nc)包含的樣本集為xci，如果xci有多個樣本，定義影響半徑rci為聚類中心到該聚類所有樣本距離的最大值，超球體影響域表示為

如果xci只有一個樣本，定義超球體影響域為以該樣本的拒絕半徑r（拒絕半徑在2.3 節(jié)中給出了介紹）作為半?yún)^(qū)間寬度的超立方體的外接超球體，表示為

2.3 拒絕域與誤差評分系數(shù)

為充分考慮樣本的空間分布特性，為每個樣本定義一個拒絕域，拒絕域為新增樣本需要規(guī)避的區(qū)域，從而避免其與初始樣本空間分布過近造成樣本冗余。

誤差評分系數(shù)與拒絕域和maxmin 準則聯(lián)合使用，一方面，誤差越大的樣本其拒絕域應(yīng)當越小，以允許候選樣本向其適當靠近；另一方面，利用maxmin 準則從若干候選點中挑選新增樣本時，通過誤差評分系數(shù)對樣本間的幾何距離進行乘法加權(quán)，更容易將靠近較大誤差樣本的候選點挑選出來，從而在其附近增加更多的樣本，以達到有效權(quán)衡全局探索和局部開發(fā)的目的。

樣本所屬的拒絕域是以該樣本點為球心的超球體，拒絕半徑為該點與剩余樣本點最小距離的一半乘以誤差評分系數(shù)，樣本xj(1 ≤j ≤ns)的拒絕半徑和拒絕域分別表示為

式中：λj為樣本xj的誤差評分系數(shù)，根據(jù)當?shù)卣`差對每個樣本進行評分。λj表示為

式中：μ 為波動系數(shù)，當μ=1 時不考慮誤差評分，此時誤差評分系數(shù)λj=1。

根據(jù)每個聚類的影響域和所有樣本的拒絕域，第i 個聚類的設(shè)計域即為

添加樣本時，根據(jù)誤差評分系數(shù)分配各設(shè)計域的新增樣本數(shù)，實現(xiàn)在較大誤差區(qū)域增添更多的樣本。第i 個設(shè)計域的誤差評分系數(shù)表示為

其 中：e?CVij表示 第i 個 設(shè) 計 域 第j 個 樣 本 的 相 對 誤差。第i 個設(shè)計域新增樣本數(shù)為niadd=λic·nadd。

此外，基于maxmin 準則能夠獲取較好空間均勻性的樣本，但是直觀上非線性程度較大的區(qū)域?qū)颖镜拿芗纫蟾?，因此計算候選樣本xi∈xcandidate與初始樣本xj∈x 的距離時，引入誤差評分系數(shù)對距離進行加權(quán)，表示為

3 算例驗證及分析

本文所提自適應(yīng)采樣方法是一種非侵入方法，可以適配不同的數(shù)學近似模型，為驗證方法的效率和可靠性，以Kriging[31-32］為例，應(yīng)用于數(shù)值測試算例，并討論相關(guān)參數(shù)的影響。

3.1 數(shù)值測試算例

附錄A 給出了測試函數(shù)用例，采用LOO 獲取交叉驗證誤差。設(shè)置每輪新增樣本數(shù)nadd=nv和波動系數(shù)μ=2，收斂準則為歸一化預期誤差閾值e?est=1% 或終止樣本數(shù)閾值nest=1 000。選取1 000 個測試樣本進行精度校驗，將所提方法與One-Shot、CV-Voronoi 和APSFC 方 法 進 行 對比。為減小隨機性影響，更客觀地評價方法的效果，每個算例運行50 次，分析其結(jié)果的統(tǒng)計特性。文獻[13］總結(jié)了初始樣本數(shù)的合適取值，考慮到初始樣本數(shù)不是本文方法的主要關(guān)注點，故不再討論其對采樣結(jié)果的影響，后續(xù)不做特殊說明處，均選取一個較小的初始樣本數(shù)ns0=5nv計算。

圖4 展示了初始樣本為6 終止樣本為16 的具有大范圍平坦區(qū)域以及多個局部極值的一維算例結(jié)果，其中圓為One-Shot 采樣，叉為初始樣本，下三角、上三角、五角星分別表示CV-Voronoi、APSFC 和FCHSW 對 應(yīng) 新 增 樣 本。One-Shot 在整個參數(shù)空間均勻采樣，APSFC 有多個新增樣本落在平坦區(qū)域且新增樣本間存在空間分布過近造成一定程度的資源浪費，F(xiàn)CHSW 和CV-Voronoi重點探索非線性區(qū)域，在平坦區(qū)域上的點較少，且FCHSW 預示響應(yīng)與真實響應(yīng)曲線基本重合，相比One-Shot、CV-Voronoi 和APSFC，F(xiàn)CHSW具有更好的全局近似效果。

圖4 一維測試算例結(jié)果對比Fig.4 Results comparison of one-dimensional numerical examples

表1 展示了9 個二維～六維測試算例在One-Shot、CV-Voronoi、APSFC 和FCHSW 這4 種方法下終止樣本數(shù)的收斂均值和均方差對比結(jié)果，對應(yīng)的缺口箱線圖如圖5 所示。本文所提FCHSW 方法在9 個測試算例中均取得最優(yōu)效果。FCHSW 與基于One-Shot 的采樣相比，采樣規(guī)模平均減少14.2%，其中最大減少幅度相比ES 達到35.7%；與基于CV-Voronoi 的采樣相比，采樣規(guī)模平均減少6.8%，其中最大減少幅度相比MZ 達到21.5%；與基于APSFC 的采樣相比，采樣規(guī)模平均減少26.1%，其中最大減少幅度相比MZ 達到45.5%。此外，根據(jù)均方差對比結(jié)果，APSFC魯棒性最差，F(xiàn)CHSW 則具有很好的魯棒性。

表1 FCHSW 與One-Shot、CV-Voronoi、APSFC 測 試結(jié)果對比Table 1 Results comparison of One-Shot， CV-Voronoi， APSFC and FCHSW

圖5 FCHSW 與One-Shot、CV-Voronoi、APSFC 測試結(jié)果缺口箱線圖對比Fig.5 Notched boxplots of calculations by One-Shot， CV-Voronoi， APSFC and FCHSW

為進一步驗證FCHSW 方法的有效性，選取具有局部高度非線性和大面積平坦分布特征的ES、PK 和MZ 測試函數(shù)為例展示了自適應(yīng)采樣詳細結(jié)果，其中ES 和PK 分別為單峰和多峰函數(shù)，MZ 除了擁有多個局部極值還具有狹長波谷帶。圖6 展示了ES、PK、MZ 計算收斂迭代圖，與One-Shot、CV-Voronoi 和APSFC 相 比，F(xiàn)CHSW能夠迅速收斂到期望誤差閾值。圖7 展示了FCHSW 分 別 針 對ES、PK 和MZ 算 例 的 新 增 樣本結(jié)果，其中正方形表示初始樣本，圓表示新增樣本。圖7（a）和圖7（b）展示了樣本在中部非線性區(qū)域采樣密集在邊緣平坦區(qū)域采樣稀疏，但是在整個平坦區(qū)域樣本分布較為均勻。圖7（c）展示了FCHSW 在非線性波谷帶采樣密集，尤其是在2 個局部極值附近更為明顯，在平坦波峰帶區(qū)域采樣相對均勻稀疏。結(jié)果表明，F(xiàn)CHSW 方法對形態(tài)良好的區(qū)域?qū)崿F(xiàn)稀疏、均勻且全覆蓋的采樣；對多極值、強非線性區(qū)域?qū)崿F(xiàn)重點密集采樣，實現(xiàn)了全局探索和局部開發(fā)的高效平衡。此外，針對不同的模型特征以及高維問題，F(xiàn)CHSW 均表現(xiàn)出良好的性能，從而驗證了方法的有效性和廣泛適用性。

圖6 典型算例下FCHSW 與One-Shot、CV-Voronoi、APSFC 的迭代收斂圖Fig.6 Convergence histories of One-Shot， CV-Voronoi， APSFC and FCHSW for typical numerical cases

圖7 FCHSW 所得典型算例的樣本分布圖Fig.7 Sample distribution of typical test cases using FCHSW

針對計算效率，通過代理模型建模次數(shù)與計算耗時2 個方面對方法進行評價。One-Shot 方法作為一次性采樣方法理論上只需要進行一次模型訓練，但是面對未知問題，通常事先難以估計達到期望精度所需樣本數(shù)，故實際操作中需要多次嘗試，且由于此類采樣方法通常追求分布均勻性與投影正交性，不能針對相關(guān)物理問題進行重點區(qū)域的針對性采樣，樣本利用效率無法保證，下文討論計算效率時不再考慮。

表2 展示了3 種自適應(yīng)采樣方法在測試算例所需的建模次數(shù)及訓練時間的統(tǒng)計結(jié)果。在建模次數(shù)方面，本文FCHSW 方法在所有測試算例中得到最小的均值與方差。在計算時間方面，在多數(shù)算例中取得了最佳或相當?shù)谋憩F(xiàn)，特別是在存在多平臺死區(qū)的復雜算例MZ 中，CV-Voronoi 和APSFC 平均耗時分別是本文所提FCHSW 方法的4.7 倍和11.7 倍。以上結(jié)果表明本文所提方法在獲得指定精度的代理模型過程中，不僅可以減少對樣本數(shù)量的需求，同時也具有良好的計算效率。

表2 計算過程中產(chǎn)生的建模次數(shù)及所需時間Table 2 Number of modeling times and time required in calculation process

3.2 進一步討論

為了驗證FCHSW 方法的有效性，本節(jié)對影響域的選擇及每輪新增樣本數(shù)nadd和波動系數(shù)μ對結(jié)果的影響進行進一步討論。

作為對比，在其他參數(shù)設(shè)置相同的情況下，圖8 展示了FCHSW 方法的2 種簡化形式。其中圖8（a）展示了不考慮誤差評分系數(shù)的拒絕域，此時各樣本的誤差評分系數(shù)λi=1(i=1，2，…，ns)，相應(yīng)的波動系數(shù)μ=1。圖8（b）展示了影響域為超立方體的情況，超立方體影響域為聚類所屬樣本在所有維度空間的跨度，定義為

圖8 FCHSW 方法的2 種簡化形式Fig.8 Two simplified forms of FCHSW method

3.2.1 關(guān)于影響域的討論

針對超立方體采樣區(qū)域使得新增樣本趨于在已有樣本內(nèi)空間，不利于全局探索的局限性，所提FCHSW 方法采用超球域，可在給定誤差分布特性下拓展不同維度的采樣范圍，提高新增樣本的潛在可探索性，進而改善RPM 的全局預示精度。

分別針對超球體和超立方體兩種影響域計算方式，對比測試算例集的自適應(yīng)采樣結(jié)果，如表3 所示。超球體影響域在9 個算例中均表現(xiàn)出不同程度的優(yōu)異性，其采樣規(guī)模相比超立方體影響域平均減少6.7%，其中最大性能提升在PK 算例中取得，達到了13.6%。

表3 影響域?qū)Ρ冉Y(jié)果Table 3 Results comparison of influential domain

3.2.2 每輪新增樣本數(shù)對結(jié)果的影響

為考察每輪新增樣本數(shù)對方法有效性影響，選取nadd=nv，2nv，3nv，4nv，5nv，結(jié)果對比如 圖9所示，具體計算結(jié)果如表4 所示。

表4 每輪新增樣本數(shù)的影響Table 4 Summary of test results about number of new samples per round

圖9 每輪新增樣本數(shù)對比結(jié)果Fig.9 Results comparison of number of new samples per round

FCHSW 整體上呈現(xiàn)隨著每輪新增樣本數(shù)增加，終止樣本數(shù)緩慢增加的趨勢，nadd=nv在9個測試算例中均取得最佳性能，且當nadd分別取2nv、3nv、4nv、5nv與nadd=nv相比，采樣規(guī)模平均增加幅度分別達到5.0%、6.9%、10.1%、13.6%，F(xiàn)CHSW 針對每輪新增樣本數(shù)具有較強魯棒性。另外，隨著nadd增大，所需迭代次數(shù)和交叉驗證計算成本成倍降低。為有效兼顧計算效率和采樣質(zhì)量，本文建議nadd取nv～3nv進行計算。

3.2.3 波動系數(shù)對結(jié)果的影響

為分析波動系數(shù)對方法有效性的影響，依次選取μ=1，1.2，1.5，2，2.5，3，5，計算結(jié)果如表5所示，結(jié)果對比如圖10 所示。

圖10 波動系數(shù)對比結(jié)果Fig.10 Results comparison of fluctuation coefficient

表5 波動系數(shù)的影響Table 5 Summary of test results about fluctuation coefficient

FCHSW 整體上隨著波動系數(shù)增大，終止樣本數(shù)先減小后增大。當μ=1.2 時，F(xiàn)CHSW 針對TG、AK、GP 算例具有最佳性能，相較于不考慮誤差評分系數(shù)分別有1.1%、0.6%、1.4%的性 能 提 升；當μ=2 時，F(xiàn)CHSW 針 對ES、PK、MZ、HT、DP 算例具有最佳性能，相較于不考慮誤差評分系數(shù)分別有31.7%、6.7%、2.2%、10.2%、2.8% 的 性 能 提 升；當μ=2.5 時，F(xiàn)CHSW 針對CH 算例具有最佳性能，相較于不考慮誤差評分系數(shù)有8.2%的性能提升。事實上，增大波動系數(shù)，意味著允許更大的樣本聚集度，進而減小全局探索同時增大了局部開發(fā)能力，應(yīng)選取一個合適的波動系數(shù)值來有效權(quán)衡二者以達到縮減采樣規(guī)模的目的。綜合考慮，本文建議μ取1.2～2.5 進行計算。

4 高超聲速飛行器氣動熱全局快速預示

本節(jié)基于FCHSW 方法進行高超聲速飛行器氣動熱RPM 的高效構(gòu)建，實現(xiàn)高精度高效率低成本獲取氣動熱的目的，對所提方法的工程有效性與實用性進行驗證。

4.1 氣動熱物理建模與高精度計算

選取類HTV-2 高超聲速飛行器作為氣動外形。HTV-2 氣動布局具有乘波體與升力體相結(jié)合的特征，外形扁平，具有尖前緣、大后掠[33］，幾何尺寸參考文獻[34］，如圖11 所示。

圖11 類HTV-2 飛行器幾何模型Fig.11 Geometry of HTV-2 type shape

飛行狀態(tài)參數(shù)為馬赫數(shù)、攻角和飛行高度，參數(shù)空間如表6 所示。給定側(cè)滑角0°，設(shè)置壁面為等溫壁面300 K 計算冷壁熱流。遠場邊界由飛行參數(shù)控制，為節(jié)省計算資源，根據(jù)流場對稱特性，采用對稱邊界進行求解，計算域及邊界條件如圖12 所示。采用344.6 萬網(wǎng)格作為最終計算網(wǎng)格，網(wǎng)格劃分如圖13 所示。

表6 飛行狀態(tài)參數(shù)空間Table 6 Design space of flight parameters

圖12 計算域及邊界條件Fig.12 Computational domain and boundary conditions

圖13 飛行器網(wǎng)格劃分示意Fig.13 Schematic of vehicle meshing

選擇定常求解，黏性控制為層流，本文不考慮真實氣體效應(yīng)等。利用有限體積法空間離散求解質(zhì)量守恒方程、Navier-Stokes 方程和能量方程，采用Ausmdv 空間離散格式和minmod 限制器格式，MUSCL 插值精度為2 階迎風。時間離散由庫朗數(shù)控制，起始庫朗數(shù)、結(jié)束庫朗數(shù)、變庫朗數(shù)步數(shù)分別設(shè)為0.01、10、500。

4.2 氣動熱RPM 構(gòu)建與更新

利用OLHD 選取15 個樣本作為初始訓練集，選取10 個測試樣本用于模型精度校驗。在飛行器表面每個離散化的網(wǎng)格節(jié)點單獨構(gòu)造一個RPM，然后所有網(wǎng)格節(jié)點的氣動熱RPM 協(xié)同預示飛行器表面熱流。

根據(jù)預示誤差分布特性，利用所提FCHSW方法獲取新增訓練樣本并更新RPM，直至模型預示精度滿足要求，其中FCHSW 相關(guān)參數(shù)設(shè)置為：nadd=2nv和μ=2，迭代終止條件為測試集平均相對誤差e?≤5%。

4.3 結(jié)果分析

圖14 為相對誤差迭代收斂過程，由圖可知每輪加點后模型精度得到顯著改善。7 次迭代后，F(xiàn)CHSW 方法通過51 個訓練樣本將RPM 的平均誤差從103.3%降低到4.8%，作為對比，相同樣本數(shù)下APSFC 方法的誤差為8.3%；而APSFC需要通過10 次迭代后的69 個訓練樣本達到4.4% 的精度。 為滿足RPM 的精度要求，F(xiàn)CHSW 所需樣本是APSFC 的73.9%，使得高精度模型計算成本明顯減小，提高了構(gòu)建RPM的效率。

此外，使用最終得到的RPM 對飛行器表面熱流進行預示，單工況耗時小于3 s，而采用CFD進行高精度計算，單個工況約需240 核時，表明本文基于FCHSW 建立的氣動熱RPM 與高精度CFD 計算相比體現(xiàn)出明顯的計算效率優(yōu)勢。

表7 給出了RPM 在10 個測試工況下精度對比，駐點熱流相對誤差在5%以內(nèi)，整體平均相對誤差在10%以內(nèi)，結(jié)果表明基于FCHSW 方法構(gòu)建的氣動熱RPM 在有限高精度計算成本下通過縮減采樣規(guī)模達到了高精高效預示氣動熱的目的。

表7 RPM 在10 個測試集上誤差對比Table 7 Results comparison of 10 test cases using RPM

考慮到駐點預示誤差和整體平均誤差的最大值均在工況Test3 取得，不失一般性，對工況Test3 進行分析。圖15 給出了工況Test3 的高精度CFD 計算熱流與預示熱流qw及誤差分布對比，可以看到在飛行器的頭部出現(xiàn)高熱流區(qū)，在側(cè)翼（前緣）出現(xiàn)狹長熱流分布帶。RPM 能準確預示熱流分布，且與高精度計算結(jié)果保持一致。圖16（a）和圖16（b）分別為RPM 與CFD在工況Test3 的迎風中心線和前緣中心線熱流對比，RPM 在迎風中心線和前緣中心線平均預示誤差分別為12.6% 和7.8%，能準確實現(xiàn)高熱流區(qū)重建，提供系統(tǒng)的主要物理特征。相對誤差較大的部分則主要集中在飛行器后半段靠近對稱軸附近的低熱流區(qū)。事實上，對于工況Test3，在所有相對誤差大于15%的區(qū)域，最大熱流為75.3 kW/m2；進一步對于相對誤差大于20%的區(qū)域，最大熱流僅為8.4 kW/m2，比駐點熱流小兩個數(shù)量級，因而可對認識飛行器氣動熱變化規(guī)律、支持熱防護系統(tǒng)設(shè)計提供有意義的結(jié)果。

圖15 Test3 工況下RPM 和CFD 表面熱流分布及誤差對比Fig.15 Comparison of surface heat flux and error contours between RPM and CFD （Test3）

圖16 Test3 工況下RPM 和CFD 迎風中心線和前緣中心線熱流對比Fig.16 Comparison of heat flux along center 1ine on windward wall and leading edge between RPM and CFD （Test3）

5 結(jié) 論

1）針對現(xiàn)有自適應(yīng)采樣方法更新RPM 難以兼顧計算效率和計算精度問題，提出一種基于模糊聚類超球體影響域和誤差評分系數(shù)加權(quán)拒絕域的批量自適應(yīng)采樣方法。數(shù)值測試算例表明，所提方法與One-Shot、APSFC、CV-Voronoi相比能有效減小RPM 構(gòu)建所需樣本規(guī)模，高效提升預示精度。

2）基于模糊聚類的超球體影響域可有效改善全局探索能力，提升采樣質(zhì)量；對FCHSW 的2個主要參數(shù)（每輪新增樣本數(shù)和波動系數(shù)）進行了對比分析，結(jié)果表明采用本文建議值可兼顧計算效率和采樣質(zhì)量，有效應(yīng)對不同維度、不同類型問題，具有較強的適應(yīng)性。

3）類HTV-2 高超聲速飛行器氣動熱快速預示實例表明，基于FCHSW 構(gòu)建的氣動熱RPM 在有限采樣規(guī)模下能夠獲得良好的全局預示精度，其中駐點預示誤差和整體平均預示誤差分別在5%和10%以內(nèi)。且所提FCHSW 方法相比APSFC 方法節(jié)約了26.1%的高精度數(shù)值計算成本，有效提升了氣動熱RPM 的構(gòu)建效率。

附錄A：數(shù)值測試函數(shù)

1）一維測試算例（多峰，高度非線性），如圖 A1 所示。

圖 A1 一維測試算例Fig.A1 1D test function

2）Trigonometric function（TG，二維，多峰，均勻分布），如圖 A2 所示。

3）Easom function （ES，二維，單峰，非線性），如圖 A3 所示。

圖 A2 TG 測試算例Fig.A2 Test function TG

圖 A3 ES 測試算例Fig.A3 Test function ES

4）Peaks function（PK，二維，多峰，高度非線性）， 如圖 A4 所示。

圖 A4 PK 測試算例Fig.A4 Test function PK

5）Ackley’s path function（AK，二維，多峰，高度非線性），如圖 A5 所示。

圖 A5 AK 測試算例Fig.A5 Test function AK

式中：a=20，b=0.2，c=2π。

6）Goldstein-Price’s function（GP，二維，多峰，尖峰）， 如圖 A6 所示。

圖 A6 GP 測試算例Fig.A6 Test function GP

7）Michalewicz’s function（MZ，二維，狹長峰谷帶），如圖 A7 所示。

圖 A7 MZ 測試算例Fig.A7 Test function MZ

8）Hartman function（HT，三維，非線性）

式中：

9） Colville-Himmelblau function （CH，四維）

10）Dixon-Price function（DP，六 維，均 勻分布）