基于高斯過程回歸的空天飛行器多精度氣動建模方法

季廷煒，查旭，謝芳芳，吳雨思，張鑫帥，蔣逸陽，杜昌平，鄭耀

(浙江大學 航空航天學院，浙江 杭州 310027)

氣動性能評估是空天飛行器[1-6]設計階段的重要環節，其中氣動力數據的獲取是空天飛行器氣動布局、結構設計、軌跡仿真等相關工作的基礎.傳統上，空天飛行器氣動力數據可以通過地面風洞實驗、飛行實驗、計算流體動力學（computational fluid dynamics，CFD）數值模擬、工程估算等多種途徑獲得.其中，工程估算計算結果置信度不高，依賴于風洞試驗結果修正，對不同外形的適應性差，一般只用于簡單外形或飛行器概念設計的氣動估算.地面試驗存在靜溫低、雷諾數低、有洞壁干擾、支架干擾等問題.飛行試驗產生的數據是客觀的，但數據量通常偏少，且來流條件存在擾動以及傳感器精度問題也會導致辨識的數據存在偏差.雖然準確可靠的氣動力數據可以通過建造更先進的風洞、采用更高性能的數值模擬系統以及進行更多更精細的飛行試驗來獲取，但在短時間內采用單一方式提高氣動力數據精準度的程度有限，不能完全滿足未來飛行器研制的需要，而且將付出高昂的成本[7].研究者希望可以將這些復雜來源的氣動力數據進行融合、建模，降低氣動數據獲取成本，最大限度的提升氣動模型建模精度.

對于氣動力多源數據融合方法的研究，目前主要圍繞以下2 類問題展開：1）氣動力數據源之間存在明顯精度差異的數據融合問題（有相對真解的數據融合問題）.受到計算成本、變量多樣、魯棒性等條件限制，高精度的氣動數據（比如飛行試驗）難以充分獲得，而低精度氣動數據的獲取較便利.考慮到低精度氣動數據與高精度氣動數據的匹配特性，低精度和高精度模型具有相近的或者有關聯的內在特征，因此可以通過混合高精度與低精度數據（或模型）建立數據融合模型以逼近更高的數據精度[8-13].2）氣動數據精度區分不明確、狀態不匹配的數據融合問題.在采用數值計算氣動數據、風洞試驗數據與飛行試驗數據對飛行器氣動力或模型進行確定時，須考慮不同狀態、不同數據來源的影響.不同氣動數據間的狀態與參數很難達到統一，這時就須針對不同狀態、不同來源的數據進行融合，提升多源氣動力數據庫的一致性和精度.近年來，機器學習在多源數據融合工作中逐步得到了嘗試和應用[14-18].它不僅能通過融合不同精度的數據實現高精確度預測，而且可以對模型的不確定性進行量化估計.因此，將傳統方法與機器學習技術緊密結合，為空天飛行器氣動問題求解，甚至流體力學的發展都提供了嶄新的思路.

傳統的工程估算方法計算速度快，但是無法考慮飛行中的諸多特殊復雜的物理現象，預測精度低；CFD 數值模擬可以高分辨率解析飛行中的物理現象，但是計算成本大、周期長.本研究基于工程估算和CFD 數值模擬方法，引入高斯過程回歸理論，提出面向空天飛行器的多精度氣動建模方法，以求在較低計算成本下實現高精度氣動特性預測.隨后，將該方法應用于FTB 空天飛行器的氣動建模問題中，構建了該飛行器的多精度氣動模型，并將氣動模型與再入軌跡仿真模型相耦合.通過對比分析仿真結果，凸顯多精度氣動建模方法的工程價值.

1 空天飛行器氣動分析方法

1.1 CFD 數值模擬

基于物理學三大守恒定律可以對流體運動進行嚴格的數學描述，從而推及由質量守恒方程、動量守恒方程、能量守恒方程構成的流體力學基本方程形式即N-S 方程組.由于空天飛行器的高速運動特性，忽略其中黏性項，進而得到三維非定常可壓縮流體的Euler 方程組，其微分形式如下：

式中：W為守恒變量，F為二階對流通量張量.

采用由美國斯坦福大學開發的開源CFD 求解器SU2，它是一款基于C++編程語言和非結構網格的高精度偏微分方程求解器，能夠解決非結構網格拓撲下的多物理場分析和優化問題.SU2 求解器的流動分析能力在RAM-C II 高超音速飛行器的試驗案例中[19]已經過嚴格驗證和確認，可確保對復雜外型的氣動特性實現準確評估.

1.2 基于工程估算的快速預測方法

針對空天飛行器的復雜外形，在應用物面壓力估算方法進行氣動力計算前，須將其表面劃分為很多小的網格面元.一般劃分為三角形或四邊形網格面元，本研究采用典型的三角形非結構網格面元.

在高超聲速條件下，基于工程估算理論求解飛行器的氣動力系數的基本流程如圖1 所示.該流程主要包括4 個部分，通過幾何生成軟件定義氣動外形，再利用網格生成軟件作非結構網格面元的劃分.將飛行器分解為機身、機翼區域，再結合迎風/背風情況，在不同區域的網格面元上采用適合的方法計算壓力系數并求和，最后進行數據轉化，將氣動力系數由機體系映射為慣性系，轉換完畢后存儲或輸出.其中，機身迎風面采用修正牛頓法：

圖1 快速預測方法的基本流程Fig.1 Basic process of fast predicting method

式中：C p 為壓力系數，δ 為物面傾角，C pmax為正激波后滯止點的壓力系數.機身背風面采用Prandtl-Meyer 法：

式中：Ma∞為來流馬赫數，γ0為比熱比.

機翼背風面采用激波膨脹波法，迎風面采用切楔/切錐法.在存在攻角 α 的情況下，為了修正物面近似后得到的錐體的物面傾角 δ，采用近似后等價錐的半頂角 δTC取代物面傾角δ：

式中：φ為等價錐的徑向角.

采用自主開發的高超聲速飛行器氣動力快速預測平臺，它遵循工程估算中氣動力的一般求解流程.該平臺是基于非結構網格文件進行計算求解的，它主要由輸入模塊、幾何處理模塊、求解模塊、數據轉換模塊、輸出模塊等模塊構成.

2 多精度高斯過程回歸

2.1 高斯過程回歸

高斯過程回歸（Gaussian process regression，GPR）是通過對樣本數據構建高斯過程，再基于貝葉斯理論求解后驗概率分布實現的，本質上是利用高斯過程先驗知識對樣本數據進行回歸分析的非參數回歸模型[20].訓練集為觀測值輸入為x，輸出為y.輸入與輸出之間存在一個未知的映射關系：

假設該映射關系服從均值為0 的高斯分布，則

式中：x'為核函數中心，k為由超參數 θ 確定的核函數.

由核函數k確定其協方差矩陣：

式中：k(xi,xj)為樣本xi、xj之間的協方差.

協方差矩陣由核函數k確定，一般選用徑向基函數作為高斯過程的核函數，定義如下：

式中：l0和 σ0為該核函數的超參數.該超參數可以通過最大化邊緣對數似然（marginal log-likelihood）來找到最優值.該對數似然函數表達式如下：

式中：n為樣本數量，p(y|x,θ) 為條件x、θ 下出現y的概率.

本研究采用L-BFGS 算法來求解該優化問題的最優參數.

2.2 基于高斯過程的多精度回歸實現

針對不同信息源或不同模型提供的數據，根據模型本身在計算精度和計算成本方面的特征，可以將數據劃分為多個精度的數據集.其中，多精度高斯過程回歸方法最適合處理由具備如圖2所示特征的模型生成的數據，它能夠利用不同精度數據實現多精度建模，獲得的代理模型精度不僅高于低精度模型，而且求解速度遠快于高精度計算模型，這是多精度高斯過程回歸模型的優勢所在.

圖2 高低精度模型的計算精度與成本分布Fig.2 Calculation fidelity versus cost distribution for high fidelity and low fidelity model

多精度高斯過程回歸是在高斯過程回歸的基礎上，基于不同精度的數據集實現的.假設不同精度數據集為Dt={(xt,yt)},t=1,···,s.其中，ys表示精度最好的數據，y1表示精度最低的數據.采用自回歸的形式來融合不同精度的數據，表達式如下：

式中：ft(x) 和ft-1(x) 分別表示精度t和精度t-1 的高斯過程回歸模型；ρ 為參數，將不同精度數據集間聯系起來；δt為隨機選取的高斯分布，與高斯分布ft-1相互獨立.假設 δt服從高斯分布，即

式中：μ 表示高斯分布對應的均值.基于以上假設，此時，精度t的ft(x) 只與下一級精度t-1相關，而與其他精度的模型無關.因此，該多精度模型在新的輸入點x*的后驗分布表達式如下：

式中：X為觀測值輸入集合.

根據貝葉斯定理，該后驗分布在精度t中的均值和方差表達式如下：

式中：nt表示精度t的訓練數據點個數，k**=k(x*,x*),.同樣的，該模型的超參數也可以由最大化式(9)的邊緣對數似然函數得到.

3 空天飛行器再入段運動模型

3.1 運動模型

在再入飛行過程中，飛行器可以視作質量不變的剛體，對其運動模型作簡化處理，僅考慮它的平移運動.所采用的算例初始再入高度為70 km，飛行器發動機推力降低為0，由于科氏加速度較小忽略不計，此時地球自轉影響可忽略不計，飛行器受到的外力僅為氣動力和重力.在再入過程中假設側滑角為0°，速度矢量始終保持在飛行器的縱向平面內.因此，可以建立飛行器在航跡坐標系下的三自由度運動模型如下：

式中：m為飛行器質量，V為速度，r為飛行器與地心之間的距離，L為升力，D為阻力，λ 為經度，Φ為緯度，σ 為傾斜角，μ 為地球引力常量，γ 為航跡角，ψ 為航向角.其中，升力和阻力表達式如下：

式中：Cl 為飛行器的升力系數，Cd 為阻力系數，Sref為參考面積.

3.2 大氣模型

為了便于軌跡仿真中大氣參數的獲取，采用基于1976 年標準大氣表發展出的表達公式[21]，該數學模型以海平面為基準，將91 km 下的大氣層分為8 層，在每層以高度h為自變量定義了大氣參數如密度、溫度、氣壓、聲速等.其中，部分高度內的大氣密度模型如圖3 所示.隨著高度的下降，大氣密度呈指數增長，其顯著變化發生在高度約為30 km 處.

圖3 0～70 km 高度的大氣密度模型Fig.3 Density model of atmosphere between 0 and 70 km

4 仿真分析

4.1 FTB 模型描述

FTB 是由意大利開發的可重復使用運載飛行器（reusable launch vehicle，RLV），如圖4 所示為FTB 模型的三視圖.機體總長為7.15 m，翼展為3.6 m，機翼面積為5.18 m2，鼻錐半徑為0.12 m.機翼前緣后掠角65°且與機身腹部呈5°二面角，該設計可以增強機體的橫向穩定性.機身的上下側具有簡單的錐形、球體幾何形狀和光滑的流線型表面，采用三角形機翼和V 形尾翼，機身和機翼均有尖銳的前緣[21-24].

圖4 FTB 模型三視圖Fig.4 Three views of FTB model

4.2 氣動分析方法比較

4.2.1 計算精度 以FTB 模型為例，比較2 種氣動分析方法的求解精度，在馬赫數等于16 的條件下CFD 數值模擬、工程估算結果與文獻數據[22]的對比如圖5 所示.2 種方法的結果隨攻角的變化趨勢基本與文獻數據吻合，其中CFD 數值模擬結果與文獻數據更吻合，工程估算結果則與文獻數據之間存在較大偏差.

圖5 升、阻力系數的計算結果對比Fig.5 Comparison of calculation results of lift and drag coefficients

式中：nt為測試集點數；為第i個測試點通過CFD 數值模擬或工程估算獲得的氣動系數，即升、阻力系數；為第i個測試點對應的測試集中的數據；分別為CFD 數值模擬和工程估算的平均誤差.2 種方法求解FTB 模型氣動系數的平均誤差如表1 所示.可以看出，升、阻力系數的計算中CFD 數值模擬的精度均比工程估算的高.在升力系數的計算中兩者差距更小，而在阻力系數方面CFD 數值模擬相對于工程估算的精度提升比率為9.56%.

表1 2 種方法求解氣動系數的平均誤差對比Tab.1 Comparison of average errors of solving aerodynamic coefficients by two methods

4.2.2 計算成本 2 種方法對于計算成本的要求如表2 所示.表中，Ns為網格數量，NC為CPU 使用數，t為計算時間.CFD 數值模擬中的網格數量更多，且在時間、空間維度上迭代求解，故在計算設備配置和計算時間方面的要求均遠高于工程估算的，尤其在計算時間方面，CFD 數值模擬所需時間是工程估算的104數量級倍.綜合計算方法的平均誤差和計算成本的對比結果，CFD 數值模擬方法相較于基于工程估算的快速預測方法可以作為高精度計算方法使用.

表2 2 種方法求解氣動系數的計算成本對比Tab.2 Comparison of calculation cost of solving aerodynamic coefficients by two method

4.3 多精度氣動建模

根據前面的結論，將通過CFD 數值模擬和基于工程估算的快速預測方法獲得的數據分別用作高精度數據和低精度數據.多精度氣動建模如圖6所示，將2 種信息源數據通過GPR 完成多精度建模.由496 個低精度樣本點和42 個高精度樣本點共同組成多精度氣動建模的數據集，從數據集中隨機選取42 個高精度樣本構建高精度訓練集，選取256 個低精度樣本構建低精度訓練集.在-5°≤α ≤25°,5 ≤Ma≤20的參數空間內構建二維氣動力系數的多精度預測模型.

4.3.1 擬合效果 選取高精度訓練集中的數據為參考值，通過繪制高精度樣本與多精度預測值間的相關性圖展現多精度高斯過程回歸對數據的擬合效果，如圖7 所示.圖中，C lHF表示升力系數的高精度樣本，ClLF、ClMF分別為升力系數的低精度樣本、多精度預測值，實線表示完全相關曲線y=x，散點表征高精度樣本與多精度預測值的聯合位置.當散點遠離完全相關曲線時，表明在該處的擬合精度較低，反之，擬合精度較高.其中，高精度樣本數N=0 時實際上為僅對256 個低精度樣本的單精度回歸結果.在多精度GPR 的結果圖中，散點與完全相關曲線幾乎完全重合，表明相較于僅用低精度數據的單精度GPR 結果，多精度GPR 對升力系數擬合效果更好.

將預測結果與訓練集中高精度樣本作差得到升力系數的擬合誤差eCl，并以云圖形式呈現.如圖8所示為2 種高精度樣本數下，預測模型給出的升力系數擬合誤差.當N=0 時，低精度代理模型在參數空間內的最大擬合誤差穩定在約0.12；當N=42 時，多精度代理模型在參數空間內的最大誤差僅為0.002 5.多精度升力系數模型的擬合精度得到了極大的提升，將升力系數的擬合誤差降低至原來的1/48，相對誤差小于1%.

4.3.2 預測效果 在參數空間內隨機選取12 個有別于訓練集的位置點，基于CFD 數值模擬構建升力系數多精度預測模型的測試集，將該測試集作為參照用以說明多精度高斯過程回歸對參數空間的預測能力.如圖9 所示為N=0，42 時升力系數測試集數據與多精度預測值之間的相關性結果.當N=42 時升力系數預測結果幾乎與完全相關曲線重合，多精度代理模型對升力系數表現出良好的預測能力.

圖9 不同數量高精度樣本下升力系數的測試集數據與多精度預測值的相關性Fig.9 Correlation between data of test set and multi-fidelity predicted values for lift coefficient under different numbers of high-fidelity samples

進一步地，以測試集中高精度樣本數據為參照，對多精度氣動模型的誤差進行定量的分析.其中，多精度氣動模型的誤差（平方誤差的平均值）定義如下：

基于測試集數據計算得到的多精度氣動模型的平方誤差的平均值及平方誤差的方差隨訓練集使用的高精度樣本數變化，如表3、4 所示.同時，僅使用訓練集中高精度樣本通過單精度GPR 構建代理模型，再基于測試集數據計算該模型的與多精度結果形成對比.根據對比結果，在FTB 氣動外形的多精度升力系數建模過程中，多精度高斯過程回歸方法基本上遵循隨著使用的高精度樣本數增加，模型精度也得到提升的規律，更重要的是它比僅用高精度數據構建的單精度模型更具精確性和穩定性.當訓練集高精度樣本數為42、低精度樣本數為256 時，升力系數代理模型平方誤差的平均值降到10-5數量級水平.此時，多精度氣動模型的預測誤差波動很小，預測穩定性高.

表3 多精度升力系數模型的平方誤差的平均值Tab.3 Mean of squared error for multi-fidelity model of lift coefficient

表4 多精度升力系數模型的平方誤差的方差Tab.4 Variance of squared error for multi-fidelity model of lift coefficient

4.3.3 氣動建模結果 類似地，對FTB 模型的阻力系數也進行多精度建模，當訓練集高精度樣本數為42、低精度樣本數為256 時，升、阻力系數的多精度建模結果如圖10 所示.此時阻力系數與升力系數模型預測效果相仿，在高精度樣本數小于等于42 范圍內達到最優.

圖10 升、阻力系數的多精度建模結果云圖Fig.10 Cloud figure of multi-fidelity modeling results for lift and drag coefficients

4.4 基于多精度氣動模型的軌跡仿真

針對空天飛行器的再入問題，將已經構建的多精度氣動模型應用在軌跡仿真中，模擬出FTB 模型在控制攻角條件下的飛行狀態變化.根據飛行動力學方程求解流程，構建軌跡仿真模型，再入軌跡仿真流程如圖11 所示.其中，初始參數包括飛行器質量m、初始高度h、飛行速度V、參考面積S、航跡角γ、航向角ψ、經度λ、緯度Φ、傾斜角σ 等.

圖11 再入軌跡仿真流程Fig.11 Process of reentry trajectory simulation

在再入過程中，存在嚴重的氣動加熱現象、巨大的動壓以及過載情況，超過一定限制會對機體本身產生破壞，一般氣動加熱取駐點處熱流率密度作為參考，熱流率密度、動壓 q、過載n定義如下：式中：Rn為駐點曲率半徑，g0為重力加速度，K為常數.

4.4.1 問題描述 在再入過程中，通常會預先定義飛行器制導策略，而定義制導策略常見的方式就是使用預先定義的攻角曲線.這里給出攻角控制方案，假設它是隨速度變化的函數，表達式如下：

在再入過程初期飛行速度和初始高度都較大，初始仿真條件如下：機體駐點曲率半徑為0.12 m，質量為2 000 kg，參考面積為5.18 m2，初始高度為70 km，速度為6 km/s，航跡角為-2°，航向角為0°，起始經緯度均為0°，傾斜角σ 為0°.根據Detra-Hidalgo 模型[21]，計算駐點熱流率密度的常數取K=5.16×10-5.由于多精度氣動模型基于參數空間-5°≤α ≤25°，5 ≤Ma≤20 中的高低精度數據構建，為了保證氣動力系數預測的有效性，仿真結束條件設置為飛行馬赫數小于5.

4.4.2 軌跡仿真結果 將僅用42 個高精度樣本、僅用256 個低精度樣本以及使用256 個低精度樣本和42 個高精度樣本構建的3 種氣動模型分別作為再入軌跡仿真的氣動數據輸入，并根據仿真條件和初始參數設置再入軌跡模型.在多精度氣動建模中，本研究已經從定性和定量的角度證明了多精度相對于單精度模型具有更好的預測效果，因此這里以基于多精度模型的軌跡仿真結果為參照，將單精度的仿真結果分別與多精度的仿真結果作差獲得再入特征變量的偏差，仿真結果及偏差對比如圖12 所示.

圖12 基于單、多精度氣動模型的再入特征曲線及偏差比較Fig.12 Curve of reentry characteristic and comparison of errors based on single-fidelity and multi-fidelity aerodynamic models

根據仿真結果，高度隨時間的變化呈震蕩曲線下降，這是空天飛行器再入過程中的典型路徑特征，它是下降過程中升阻力系數、速度、空氣密度等綜合影響的結果.根據熱流率密度曲線、過載曲線、動壓曲線，熱流率密度和過載為再入初期影響飛行的主要因素，隨著飛行速度下降和高度下降，空氣密度急速上升，動壓對飛行器狀態影響逐漸顯著.偏差對比結果顯示，通過多精度GPR、高精度GPR 構建的氣動模型分別作為輸入獲得的再入曲線相差較小，通過低精度GPR 構建的氣動模型作為輸入獲得的再入曲線與前面兩者相差較大，采用低精度GPR 的最大仿真偏差約為采用高精度GPR 的10～100 倍.具體地，初始高度、速度較大，熱流率密度、動壓為速度的函數，因此低精度GPR 構建的氣動模型對再入高度、速度以及熱流率密度、動壓偏差均帶來顯著影響.其中，熱流率密度偏差絕對值最大至105數量級，多精度氣動建模的工程意義也因此得到體現.由于不能獲取軌跡仿真所需要的全部高精度氣動數據，無法對多精度氣動模型獲得的再入仿真結果給出定量的誤差分析.

5 結論

(1)以FTB 氣動外形為研究對象，分析得出相較于工程估算，CFD 數值模擬方法求解阻力系數的精度提升9.56%，而升力系數的結果差距略小.但是，CFD 數值模擬計算成本要遠高于工程估算的.

(2)采用多精度高斯過程回歸方法，能以低計算成本完成該空天飛行器的多精度氣動建模，且模型精度較高.當訓練集使用42 個高精度樣本、256 個低精度樣本時，所提出的多精度氣動模型均方誤差低至10-5數量級，且比僅用相同數量高精度數據構建的單精度模型精確度和穩定性更高.

(3)將多精度氣動模型與再入軌跡仿真耦合，模擬再入過程中的重要狀態參數變化，對比分析單、多精度氣動模型的仿真結果，發現氣動模型的微小差異對再入高度、速度、熱流密度、動壓影響顯著，進一步驗證了所提出的多精度氣動建模方法對于空天飛行器再入過程軌跡高精度仿真的工程意義.

(4)本研究在多精度氣動建模中，并未對訓練集和測試集的選用策略作深入考量，最終采用模型是否為全局最優模型也未分析論證.在未來研究中可將優化過程和多精度氣動建模方法相結合，選用不同采樣策略，通過迭代尋優給出目標誤差下的最優氣動模型.