基于自動核構造高斯過程的導彈氣動性能預測

胡偉杰，黃增輝，劉學軍,*，呂宏強

1. 南京航空航天大學 計算機科學與技術學院/人工智能學院，模式分析與機器智能工業(yè)和信息化部重點實驗室，南京 211106 2.軟件新技術與產(chǎn)業(yè)化協(xié)同創(chuàng)新中心，南京 210023 3.南京航空航天大學 航空學院，南京 210016

在飛行器的初步設計階段，探索一種快速且精確的導彈氣動性能預測方案對導彈氣動設計具有重要價值。合理有效的氣動性能預測方法對導彈外形設計、彈道設計、彈體機動性和控制系統(tǒng)的穩(wěn)定性等具有重要意義。傳統(tǒng)的氣動性能預測方案大致分為3類，即飛行試驗、風洞試驗和理論計算(包括工程計算和數(shù)值計算)[1]。飛行試驗和風洞試驗可以準確地模擬導彈在真實飛行環(huán)境中的流場條件，能夠得到精確的氣動數(shù)據(jù)。但這種方法過度依賴設計人員的經(jīng)驗，另外試驗周期長、成本高，限制了導彈外形設計空間的推廣。

為了彌補這一缺陷，以美國空軍力學實驗室開發(fā)的Missile Datcom為代表的氣動分析和快速估算軟件，利用美國空軍幾十年的風洞試驗數(shù)據(jù)建立典型氣動模型，基于部件組合法并采用經(jīng)驗與理論相結合的方式快速預測氣動性能，對常規(guī)導彈配置具有較高的氣動表征能力，能夠基本滿足導彈設計初期階段的精度要求，因而在過去幾十年間得到廣泛的應用[2]。但對于氣動外形復雜、氣動干擾嚴重或伴有大攻角飛行姿態(tài)條件下的導彈氣動性能預測，仍具有很大的局限性，需要結合風洞數(shù)據(jù)進行矯正[3-6]。文獻[3]指出，超聲速下對于與尾翼緊密耦合的低展弦比邊條翼，Missile Datcom嚴重高估邊條翼渦對尾翼的下洗效應，導致導彈法向力預測值偏低，壓力中心前移。文獻[6]指出，跨聲速下的激波使得流場具有強烈的非線性效應，Missile Datcom基于線性理論模型的方法難以對高度非線性流體流動產(chǎn)生解析結果。與風洞試驗數(shù)據(jù)對比，尾翼法向力系數(shù)曲線斜率預測誤差達到50%。因此，Missile Datcom盡管具有快速預測的能力，但在很多配置和飛行條件下的預測精度很差，在工程應用上有待改進。

隨著流體力學和計算機科學相互融合，計算流體力學(Computational Fluid Dynamics， CFD)得到快速發(fā)展，借助高性能計算機的計算能力，在流動基本方程的控制下對復雜物理流動進行數(shù)值模擬，可以得到滿足工程要求的高精度流場解，彌補了工程估算軟件預測精度不足的缺陷。近年來，CFD技術大量應用于導彈氣動性能分析任務，如文獻[1]使用以Fluent為代表的CFD流場計算軟件求解某型導彈的氣動性能，實驗結果表明，比工程估算方法具有更高的計算精度；文獻[7]利用CFD計算155 mm自旋穩(wěn)定射彈在偏航角為40°時的氣動性能，減少了風洞試驗的成本。CFD數(shù)值計算方法雖然可以得到高精度的氣動解，但計算時間太長，對于網(wǎng)格量較大的模型，求解一個算例通常需要數(shù)小時甚至幾天。如果解的收斂性不好，需重新調(diào)整網(wǎng)格或設置求解器參數(shù)，將耗費更長的時間，這在導彈設計初期階段是不能接受的。因而，迫切需要一種既快速又精確且能處理高度非線性問題的氣動性能預測方案。

隨著數(shù)據(jù)科學、機器學習技術的發(fā)展和計算機計算能力的提高，基于數(shù)據(jù)驅動的氣動建模方案在航空航天各個領域得到應用[8-9]。為彌補傳統(tǒng)CFD計算時間長的缺陷，代理模型(Surrogate Model)技術成為新的氣動性能預測選擇[10]。代理模型是基于一組可信數(shù)據(jù)，運用相適應的機器學習方法對問題抽象建模，從數(shù)據(jù)中提取輸入?yún)?shù)與輸出響應之間近似關系的模型。常用的代理模型有多項式回歸(Polynomial Regression)[10]、支持向量回歸(Support Vector Regression，SVR)[11]、神經(jīng)網(wǎng)絡(Neural Network，NN)[12-13]、Kriging模型[10,14-16]、高斯過程回歸(Gaussian Process Regression， GPR)[17]和深度網(wǎng)絡[18-19]等。近年來，代理模型在工業(yè)設計中得到了大量應用。文獻[20]基于代理模型的優(yōu)化算法提出一種面向工程應用的多輪次高效全局氣動優(yōu)化設計方法，該方法應用在大型民機超臨界機翼氣動設計中，驗證了方法有效性和工程實用性。文獻[21]基于梯度增強型Kriging模型的代理優(yōu)化算法開展了翼型設計變量的氣動反設計，得到很好的優(yōu)化效果。文獻[22]基于高斯過程代理模型提出一種混合差分進化和雜草算法的翼型優(yōu)化框架，在升力系數(shù)和橫截面積不會小很多的約束下，顯著提高了翼型的升阻比。可見，代理模型在氣動優(yōu)化設計領域的研究已經(jīng)很深入，但應用于導彈氣動性能預測的研究相對較少。文獻[11]基于支持向量回歸的方法構建火箭氣動學科代理模型，氣動性能預測的效率遠高于CFD模型，預測精度也能夠滿足精度要求，但無法提供對預測結果的置信度。文獻[12-13]基于人工神經(jīng)網(wǎng)絡構建導彈氣動性能快速預測代理模型，獲得了較高的預測精度，但神經(jīng)網(wǎng)絡參數(shù)多，模型訓練容易過擬合，且神經(jīng)元初始參數(shù)設置的隨機性導致預測結果不穩(wěn)定。文獻[16]對比了多種代理模型用于導彈氣動性能預測的結果，得出各向異性的Kriging模型建模效果最優(yōu)，預測結果穩(wěn)定且能夠提供對預測結果的置信度。本質(zhì)上，具有高斯隨機場假設的Kriging模型等價于正態(tài)似然下的GPR模型[23]。

除文獻[16]之外，還未見有GPR代理模型在導彈氣動性能預測上的應用，鑒于GPR模型相對其他模型的優(yōu)勢，本文基于GPR模型完成對典型導彈幾何外形與對應氣動性能關系的建模。當對問題本身缺乏足夠的先驗(Prior)知識，傳統(tǒng)的GPR方法很難構建適合數(shù)據(jù)內(nèi)在結構特性的模型，導致數(shù)據(jù)擬合能力不強或者模型的泛化能力差等各種問題。針對這一問題，借鑒文獻[24]的方法：在傳統(tǒng)的GPR模型框架下，嵌入一種自動核構造(Automatic Kernel Construction，AKC)算法，該算法能夠自動捕獲數(shù)據(jù)內(nèi)在結構特性，充分挖掘數(shù)據(jù)中隱含的模式信息，從而建立更精準的代理模型。目前，該方法已在很多領域得到廣泛應用。例如，文獻[25]將該方法用于預測大氣環(huán)境中超細粒子數(shù)濃度，根據(jù)特定數(shù)據(jù)自動構造基于線性核和有理二次核相加的非平穩(wěn)核函數(shù)，能夠準確地描述粒子數(shù)濃度的動態(tài)變化規(guī)律。文獻[26]基于該方法預測鋰電池的剩余使用壽命，強調(diào)了使用自動核構造函數(shù)來捕捉復雜行為的優(yōu)勢。文獻[27]基于此方法自動構造出準周期核函數(shù)，用于建模和預測全球水平輻射，結果表明該模型優(yōu)于簡單核的高斯過程模型。本文針對導彈特定應用的問題，將這種改進的GPR模型(簡稱AKC-GPR)與傳統(tǒng)的GPR模型在特定導彈氣動建模任務上進行氣動性能預測對比試驗。試驗結果驗證了使用AKC-GPR模型在導彈氣動性能快速預測應用上的良好性能。

1 高斯過程回歸

高斯過程(Gaussian Process，GP)模型是一種融合機器學習和統(tǒng)計學于一體的應用貝葉斯推理的概率建模框架[17]。GPR是在GP先驗假設下對數(shù)據(jù)進行回歸分析的非參數(shù)化模型，能夠給出對預測樣本的完整后驗分布，提供預測結果的不確定性，使得模型具有良好的可解釋性。當似然為正態(tài)分布時，預測樣本的后驗分布具有高斯形式的閉式解。與線性回歸(Linear Regression，LR)相比，GPR具有更強的非線性表達能力；與SVR相比，GPR模型預測結果具有更好的解釋性；與NN相比，GPR需要估計的參數(shù)較少，從而減少了對復雜優(yōu)化或正則化方案的需求。此外，GPR具有容易訓練，預測結果穩(wěn)定等優(yōu)點。但是GPR模型訓練涉及到協(xié)方差矩陣求逆運算，計算復雜度為O(n3)，其中n表示訓練樣本數(shù)量。因此，GPR適用于小樣本學習問題，目前已經(jīng)在眾多領域得到成功應用[28-29]。

1.1 模型推導

高斯過程是隨機變量的集合，任意有限維隨機變量的聯(lián)合分布是高斯分布[17]。因此，直接從函數(shù)空間角度推理，指定目標函數(shù)f(x)的先驗為一個高斯過程GP，f(x)的分布為

f(x)～GP(m(x),k(x,x′))

(1)

式中：m(x)為GP的均值函數(shù)(Mean Function)；k(x,x′)為GP的協(xié)方差函數(shù)(Covariance Function)或核函數(shù)(Kernel Function)，定義為

m(x)=Ε[f(x)]

(2)

k(x,x′)=cov(f(x),f(x′))

(3)

考慮一般的回歸模型：

yi=f(xi)+εi

(4)

為了推理結果表達的簡潔性，一般地，將均值函數(shù)m(x)設為0均值函數(shù)。因此，對于給定的訓練樣本輸入X，其對應輸出y的多元高斯分布為

(5)

式中：K(X,X)表示n×n的對稱正定協(xié)方差矩陣，矩陣元素為k(x,x′)，利用輸入變量x、x′之間的關系模擬隨機變量f(x)與f(x′)間的相關性；I是n×n的單位矩陣。

根據(jù)訓練樣本集D，回歸模型的目標是給定新的樣本輸入向量x*求解輸出y*，其中x*∈Rd，y*∈R。

則訓練樣本的觀測值y和待測樣本輸出y*的聯(lián)合先驗高斯分布為

(6)

式中：K(X,x*)=K(x*,X)Τ,度量訓練樣本輸入X與待測樣本輸入x*間的協(xié)方差；k(x*,x*)表示x*自身的方差。

求式(6)關于觀測向量y的條件概率分布，得到待測樣本輸出y*的后驗分布為

(7)

后驗分布為高斯分布，均值和方差公式具體為

(8)

(9)

1.2 模型訓練

根據(jù)式(8)和式(9)可知，0均值高斯過程先驗下的GPR模型在給定任務上的性能完全由核函數(shù)k(x,x′)決定。因而，GPR模型的訓練過程就是基于一組訓練數(shù)據(jù)，對核函數(shù)的形式和超參數(shù)(Hyper-parameter)進行推理，最終得到能夠描述數(shù)據(jù)內(nèi)在隱含模式的模型[17]。假定選擇的核函數(shù)形式為平方指數(shù)(Squared Exponential，SE)核函數(shù)，其表達式為

(10)

當協(xié)方差函數(shù)確定之后，GPR訓練的目標就變成優(yōu)化合適的超參數(shù)θ。常用最大化邊際似然(Maximum Marginal Likelihood)方法求解超參數(shù)的值。在GP假設下，利用高斯過程的邊際化屬性(Marginal Property)，通過先驗與似然的積分可以得到邊際似然p(y|X,θ)。為了求導方便，常取對數(shù)形式，表達式為

L(θ)=lgp(y|X,θ)=

(11)

(12)

式(12)是一個非凸優(yōu)化(Non-convex Optimization)問題，會存在多個局部最優(yōu)解。為了跳出局部最優(yōu)，最終獲得全局最優(yōu)解，可采用諸如遺傳算法、模擬退火等隨機搜索算法，但同時伴隨搜索速度慢，優(yōu)化效率低等問題。基于梯度的優(yōu)化算法容易陷入局部最優(yōu)，但其優(yōu)化效率比上述隨機搜索算法高。因此從效率角度考慮，本文采用基于梯度的優(yōu)化算法優(yōu)化超參數(shù)θ的值。

為降低超參數(shù)θ陷入局部最優(yōu)的可能性，需要設置適當?shù)某瑓?shù)初始值。一種有效的方法是假設超參數(shù)初始值服從某一先驗分布p(θ)，從此分布中采樣若干候選的超參數(shù)初始值，并采用最大化邊際似然的方法選擇使得似然最大的候選值作為最終超參數(shù)的初始值，文獻[28]驗證了這種方法的可行性。

2 自動核構造算法

傳統(tǒng)的GPR學習通常需要給定核函數(shù)的參數(shù)化形式，而選擇合適的核函數(shù)是GPR建模的關鍵。核函數(shù)編碼了對目標函數(shù)的先驗假設，表達了對數(shù)據(jù)內(nèi)在結構特性的描述[30]。如果對建模問題的先驗知識掌握不夠，很難在廣闊的核結構空間(Kernel Structure Space)內(nèi)找到適合特定問題的最優(yōu)核函數(shù)。

有理論表示，如果給定充足的數(shù)據(jù)，在某些條件下，選擇SE核的GPR能夠學習任意的連續(xù)函數(shù)[31]。但是，SE核是一種簡單核(Simple Kernel)，只能捕捉到數(shù)據(jù)中的一種結構特性。當面臨復雜的建模問題，SE核將無法編碼內(nèi)含多種屬性或結構特性的數(shù)據(jù)，所以需要選擇更復雜的核函數(shù)來描述問題。

有效的核函數(shù)必須是半正定(Positive Semidefinite)核，半正定核在“加”(Addition)和“乘”(Multiplication)運算符下是封閉的[17]。因此可以通過“加”“乘”運算將簡單核組合成復合核(Composite Kernel)，進而捕捉數(shù)據(jù)中更復雜的結構特性。文獻[17]總結了基于“加”“乘”運算符組合新核的方法，這種組合方式需要人為指定復合核中每一組件所使用的簡單核，因此依然需要先驗知識。文獻[32]提出一種稱為加性核(Additive Kernel，Add)的復合核，它是基于“加”“乘”運算符，為每一維特征指定一個簡單核，由所有可能的一維簡單核乘積的加權和組成。它考慮了一維簡單核在不同階數(shù)下的交互作用，能夠靈活地捕捉到數(shù)據(jù)中隱含的結構信息。文獻[28]將該方法應用于風洞馬赫數(shù)控制，驗證了Add核相較于傳統(tǒng)的SE核具有更強的學習能力。但是這種核的組合形式相對固定，并不能保證在各種數(shù)據(jù)上都能學到最優(yōu)的核函數(shù)形式。

針對上述簡單核和復合核的建模能力的局限性，本文使用AKC算法，在特定數(shù)據(jù)上自動構造出核結構空間中最優(yōu)的核函數(shù)。AKC算法的核心是定義簡單核的類型和構造規(guī)則(Construct Rules)，核的構造過程不需要依賴先驗知識。

2.1 簡單核及其先驗假設

常用的簡單核有SE核、線性核(Linear，Lin)、周期核(Periodic，Per)、有理二次核(Rational Quadratic，RQ)和白噪聲核(White Noise，WN)。AKC算法中的簡單核是應用在x的某一特征維度上，因此只給出簡單核的一維表示形式：

(13)

(14)

(15)

(16)

kWN(x,x′)=σ2δx,x′

(17)

式中：{σ,,σb,σv,p,α}表示核函數(shù)的超參數(shù)；δx,x′表示克羅內(nèi)克函數(shù)(Kronecker Delta Function)。

每一類簡單核都編碼了對目標函數(shù)不同的假設。Lin核假設函數(shù)具有線性變化結構；SE核假設函數(shù)具有局部變化結構；Per核假設函數(shù)具有重復變化結構；RQ核假設函數(shù)具有多尺度變化結構；WN核假設函數(shù)具有獨立噪聲結構。

使用“加”“乘”運算符對具有不同假設的簡單核按照構造規(guī)則自動組合，形成具有不同復雜結構特性的核函數(shù)。例如，Lin核“加”Lin核的組合假設函數(shù)是二次的；SE核“加”Per核的組合假設函數(shù)是局部周期變化的；Lin核“乘”SE核的組合假設函數(shù)具有加速變化的結構特征。

2.2 構造規(guī)則

構造規(guī)則定義了核函數(shù)的組合方式或構造方式，不同的構造規(guī)則能夠產(chǎn)生具有不同形式的核函數(shù)，進而能夠描述不同的復雜結構特性。一般地，由“加”“乘”和“替換”運算符定義的構造規(guī)則，即可構造出具有各種結構的核函數(shù)。

①S→S+B

②S→S×B

③S→B

式中：S表示當前核函數(shù)中的任意組件(Component)，B表示任意類型的簡單核。規(guī)則①代表“加”規(guī)則(Addition Rule)，表示任意組件S可以加任意的簡單核B；規(guī)則②代表“乘”規(guī)則(Multiplication Rule)，表示任意的組件S可以乘任意的簡單核B；規(guī)則③代表“替換”規(guī)則(Replacement Rule)，表示任意組件S可以被任意簡單核替換。

2.3 AKC算法描述

基于簡單核和構造規(guī)則的描述，圖1展示核構造的大致過程。其中假定簡單核選擇SE核和Lin核，構造規(guī)則為2.2節(jié)描述的3條規(guī)則，初始核選擇WN核。

圖1 構造過程示例Fig.1 Example of construction procedure

AKC算法中涉及到的變量和函數(shù)，具體含義如表1所示。算法流程大致分為初始化、搜索和訓練3個階段。在初始化階段，設定最大迭代次數(shù)等參數(shù)；在搜索階段，某次迭代中，基于簡單核集合和構造規(guī)則集合進行搜索，形成核結構空間；在訓練階段，使用某種評分準則如貝葉斯信息準則(Bayesian Information Criterion，BIC)，對核結構空間內(nèi)的每一核函數(shù)進行評分[33]，選擇評分最高的核函數(shù)，進入下一個迭代過程或達到停止條件退出。詳細的算法描述如算法1所示。

表1 變量和函數(shù)描述Table 1 Description of variables and functions

算法1 自動核構造Algorithm 1 Automatic kernel construction

3 數(shù)據(jù)生成

本文使用CFD方法生成導彈氣動性能模擬數(shù)據(jù)集，在該數(shù)據(jù)集上對傳統(tǒng)的GPR和AKC-GPR進行對比實驗。數(shù)據(jù)生成過程包括幾何建模、網(wǎng)格劃分、數(shù)值求解和后處理4部分。

3.1 幾何建模

使用CAD幾何建模軟件進行建模，繪制導彈的幾何形體，生成幾何模型文件。采用典型的旋成體加尾翼的構型，尾翼沿軸線呈“×”形分布，彈體總長8 m保持不變。整體構型如圖2所示。

圖2 導彈幾何構型Fig.2 Geometry of missile

基于上述幾何構型，定義3個參數(shù)描述導彈外形,分別為：導彈首部長度L、彈體半徑R和尾翼縱向相對長度F(如圖3所示)。L取1.6 m、2.0 m 和2.4 m；R取0.4 m、0.5 m和0.6 m；F為以初始設計翼型長度為基準進行的等比例放縮值，分別取0.8、1和1.25，由此構建27組不同尺寸構型的導彈模型。

圖3 導彈幾何參數(shù)Fig.3 Geometric parameters of missile

3.2 網(wǎng)格劃分

將模型文件導入網(wǎng)格生成軟件中，對計算域進行網(wǎng)格劃分，設置附面層第1層高度為0.02 mm，整個流場計算域總長約為10倍的彈體長度。根據(jù)導彈幾何外形，生成結構化面網(wǎng)格如圖4所示。

圖4 導彈面網(wǎng)格Fig.4 Surface mesh of missile

3.3 數(shù)值求解

導入網(wǎng)格文件到流體計算軟件中，選取基于密度的隱式求解器，采用implicit格式以保證迭代穩(wěn)定性，同時打開能量方程，采用SSTk-ε湍流模型。密度和黏性設為隨溫度而變化，馬赫數(shù)設置為5始終不變，攻角變化范圍在6°～16.5°。設定收斂控制參數(shù)，取Courant數(shù)為1。欠松弛因子的設置：湍動能為0.8，湍流耗散率為0.8，湍流黏度為0.8。通過迭代法對離散化代數(shù)方程組求解，得到滿足收斂性要求的數(shù)值解。

3.4 后處理

根據(jù)理論計算公式，求解出相應的氣動性能數(shù)值，本文選取升力系數(shù)CL、阻力系數(shù)CD和力矩系數(shù)Cm作為氣動性能預測對象。由此生成以外形參數(shù)R、L、F和攻角α作為輸入，3種氣動系數(shù)為輸出的模擬數(shù)據(jù)集，共182個算例。部分數(shù)據(jù)如表2所示。

表2 部分數(shù)據(jù)集Table 2 Part of dataset

4 核函數(shù)的選擇

由于對CFD生成的導彈氣動數(shù)據(jù)內(nèi)在結構特性沒有足夠的先驗信息，很難指定合適的核函數(shù)，因此，本文使用AKC算法，基于給定數(shù)據(jù)集自動構造核函數(shù)，并在GPR模型框架下，與其他通用的核函數(shù)進行對比。本文使用的核函數(shù)如表3所示，各核函數(shù)具體表達式見文獻[34]。

表3 各種類型的核函數(shù)Table 3 Varieties of kernel functions

5 實驗結果及分析

5.1 GPR模型非線性擬合能力測試

GPR模型在小樣本情形下表現(xiàn)出比其他回歸模型更強的非線性擬合能力，尤其是核函數(shù)選擇AKC時，將進一步增強模型的非線性擬合能力。為了驗證該結論，本文設計了一個一維非線性函數(shù)f(x)=sin(x)+cos(2x)+0.5，x取值范圍設定為[0,2π]，在此區(qū)間上隨機采樣10個點作為訓練樣本。采樣時，為模擬受真實世界噪聲的影響，假設樣本被高斯噪聲侵蝕。使用這10個訓練樣本進行模型訓練并預測非線性函數(shù)曲線。使用根均方誤差(Root Mean Square Error，RMSE)指標度量曲線擬合精度。

(18)

圖5展示了SVR、NN和GPR 3種模型對函數(shù)曲線的擬合結果，其中GPR模型的核函數(shù)選擇常規(guī)的SE核函數(shù)，RMSE分別為0.761 1、0.387 1和0.104 6。可知，GPR模型的擬合能力最強。在采樣點附近，3種模型尚且能夠給出較好的預測，但在沒有采樣到樣本的區(qū)間，只有GPR模型仍然能夠較好地表征函數(shù)曲線的非線性變化規(guī)律。由此可見，GPR模型具有在小樣本情形下非線性擬合的獨到優(yōu)勢。

圖5 不同模型的曲線擬合結果Fig.5 Curve fitting results of different models

圖6展示了使用SE核函數(shù)和AKC核函數(shù)的GPR模型對非線性函數(shù)曲線的擬合結果。從圖中看出，與真實曲線偏差較大的區(qū)域，使用AKC核函數(shù)擬合的曲線偏差更小一些，RMSE分別為0.104 6和0.089 2。可見，AKC核函數(shù)的GPR模型曲線擬合能力進一步增強。

圖6 不同核函數(shù)的曲線擬合結果Fig.6 Curve fitting results using different kernel functions

5.2 不同模型對比

為了驗證GPR模型在特定導彈數(shù)據(jù)集上的預測性能優(yōu)于其他常用的回歸模型，實驗對比了SVR、LR、NN和GPR 4種模型的預測結果。分別使用4種模型預測給定導彈外形參數(shù)和攻角下的升力系數(shù)、阻力系數(shù)和力矩系數(shù)，使用平均絕對百分比誤差(Mean Absolute Percentage Error，MAPE)指標度量預測結果的精度。

(19)

由于數(shù)據(jù)量較小，為了使得預測結果具有可信度，本文使用5折交叉驗證(Cross Validation)對各模型進行訓練，GPR模型的核函數(shù)選擇最常用的SE-ARD核函數(shù)。

圖7給出4種不同代理模型的預測精度對比。SVR模型對3種氣動系數(shù)的預測MAPE在1.10%～2.39%之間；LR模型的預測MAPE在1.03%～1.91%之間；NN模型的預測MAPE在0.52%～0.80%之間；GPR模型的預測MAPE在0.24%～0.36%之間。由此可見，使用SE-ARD核函數(shù)的GPR模型的預測精度高于LR和SVR模型，略高于NN模型。從而驗證了在導彈氣動性能預測精度上，GPR模型是最優(yōu)的。

圖7 不同模型性能度量Fig.7 Performance measures for different models

為驗證GPR模型在預測結果穩(wěn)定性方面的優(yōu)勢，對NN和GPR兩種代理模型進行對比實驗。圖8給出在重復10次5折交叉驗證實驗下，GPR和NN模型對3種系數(shù)預測精度的對比。可以看出，GPR模型在重復10次實驗下MAPE值保持不變，說明GPR模型預測結果非常穩(wěn)定；而NN模型的10次預測結果差異很大，MAPE最小值為0.29%，最大值達到1.04%，說明NN的預測結果很不穩(wěn)定。由于NN屬于參數(shù)化模型，網(wǎng)絡層數(shù)、每層神經(jīng)元數(shù)量、參數(shù)初始化等都沒有一個標準可循，尤其是在小樣本情況下，NN很難學得一個穩(wěn)定可靠的模型。

圖8 預測結果穩(wěn)定性對比Fig.8 Comparison of stability of predicted results

由此可以得出結論：無論從預測精度還是預測結果穩(wěn)定性角度來看，作為非參數(shù)化模型的GPR模型是最優(yōu)的，可以作為導彈氣動性能快速預測的一個合理選擇。

5.3 不同核函數(shù)對比

為了準確挖掘導彈氣動數(shù)據(jù)中潛在的模式信息，彌補傳統(tǒng)GPR模型不能準確捕獲數(shù)據(jù)內(nèi)在結構特性的缺陷，本文基于AKC算法從數(shù)據(jù)中預先構造具有特定結構的核函數(shù)，將此核函數(shù)和表3中的其他核函數(shù)同時嵌入GPR模型框架下，對比不同核函數(shù)的預測性能。其中，AKC算法中簡單核選擇SE核和WN核，構造規(guī)則為2.2節(jié)的3條規(guī)則，最大迭代次數(shù)設置為10，評分準則為BIC，初始核為WN核。

表4展示了應用各種核函數(shù)的GPR模型在5折交叉驗證下的預測結果對比。由表中數(shù)據(jù)可知，多項式核函數(shù)(Poly-3)和SE乘SE核函數(shù)(Prod)的預測誤差超過0.5%，預測結果較好，但不能準確捕捉數(shù)據(jù)內(nèi)在的結構特性，且預測結果的不確定性較大。各向異性Per核(Per-ARD)的預測誤差在0.3%左右，預測精度很高，但其預測結果的不確定性較大，說明氣動數(shù)據(jù)中不存在周期結構特性。其他類型的核函數(shù)預測結果都能達到很高的精度要求，標準差都基本控制在2×10-4左右，預測結果不確定性很小。特別地，基于構造核(AKC)的GPR模型的氣動性能預測結果在3種系數(shù)的預測上都達到最高的精度，并且對應的預測不確定性也最小。由此可見，基于自動構造核的方法得到的核函數(shù)比依靠人工先驗選擇的核函數(shù)更能夠捕捉到數(shù)據(jù)內(nèi)在的結構特性，模型表達能力更強且預測結果更穩(wěn)定可靠。

除了使用交叉驗證法預測不同核函數(shù)的GPR模型的預測效果，本文還使用留出法(Hold-out)對不同GPR模型進行訓練，即按不同比例將數(shù)據(jù)集劃分為互不相交的兩部分，分別用于訓練和測試。圖9分別展示了3種氣動性能在不同的數(shù)據(jù)劃分比例下，各核函數(shù)的GPR模型預測結果比較。橫坐標表示訓練集所占比例，縱坐標表示MAPE。為了降低隨機性的影響，對于每種比例的訓練數(shù)據(jù)集重復采樣10次，取平均結果。由圖9可知，隨著訓練集比例的增加，各核函數(shù)的GPR

圖9 不同劃分比例下的實驗結果Fig.9 Experimental results at different partition ratios

模型的預測誤差基本呈下降趨勢，且AKC-GPR模型對3種系數(shù)的預測結果在所有訓練集比例下都能得到最高的預測精度。由此可見，無先驗的AKC-GPR模型相對于傳統(tǒng)的GPR模型確實能夠更加準確地捕獲特定數(shù)據(jù)中的內(nèi)在結構特性，使得GPR模型的預測性能進一步得到提升。

5.4 殘差分析

為了檢驗AKC-GPR建模方法的合理性，需要對預測結果的殘差作正態(tài)性檢驗。正態(tài)性檢驗就是檢驗模型產(chǎn)生的殘差是否滿足正態(tài)分布假設，即是否服從以0為均值、某一常數(shù)為方差的正態(tài)分布。如果滿足該假設，則說明模型產(chǎn)生的殘差是隨機的，不是模型自身原因造成的，從而說明建模方法是合理的。本文對AKC-GPR的預測殘差進行分析。

圖10顯示了3種氣動系數(shù)的殘差直方圖，從圖中明顯看出，3種系數(shù)的殘差均服從均值為0、方差為某一常數(shù)的正態(tài)分布，即滿足正態(tài)假設，進而說明基于AKC-GPR的建模方法對于導彈氣動性能預測任務是可行且合理的。

圖10 殘差直方圖Fig.10 Residual histogram

5.5 氣動性能快速預測應用

為了展示AKC-GPR模型在實際應用中的預測效果，隨機選擇其中26組外形數(shù)據(jù)對AKC-GPR模型進行訓練，用訓練好的模型對剩下的1組 外形完全未見的導彈氣動性能進行預測。本文選擇R=0.6 m、L=2.4 m、F=1.0作為未見外形，從其余26組訓練外形中選擇一個未見外形所在樣本空間范圍內(nèi)附近的已見外形作對比，例如R=0.6 m、L=2.4 m、F=0.8。圖11分別展示了使用AKC-GPR模型對已見外形進行預測得到的3種氣動系數(shù)曲線。從圖中可以看出，預測曲線與真實氣動特性曲線完全重合，由此說明AKC-GPR模型能極其精確地描述隱含于訓練數(shù)據(jù)中的特性，充分挖掘導彈外形與氣動性能之間的對應關系。

圖11 氣動系數(shù)曲線(已見外形)Fig.11 Aerodynamic coefficient curves (shape seen)

圖12分別展示了使用AKC-GPR模型對未見外形進行預測得到的3種氣動系數(shù)曲線。從圖中可以看出，3種系數(shù)預測曲線與真實曲線基本吻合，且真實值基本落在95%的置信區(qū)間內(nèi)。由此可以說明通過對已有外形的氣動數(shù)據(jù)進行學習，AKC-GPR模型可以獲得導彈外形和氣動性能之間的對應關系，實現(xiàn)在給定工況范圍內(nèi)對新外形的導彈氣動性能的準確預測。

圖12 氣動系數(shù)曲線(未見外形)Fig.12 Aerodynamic coefficient curves (shape unseen)

6 結 論

針對傳統(tǒng)導彈氣動性能預測方法的局限，本文使用高斯過程回歸代理模型對定常流場條件下的典型導彈氣動性能進行預測，得到結論如下：

1) 對比線性回歸(LR)、支持向量回歸(SVR)、神經(jīng)網(wǎng)絡(NN)和SE核的高斯過程回歸(GPR)4種代理模型，對升力系數(shù)、阻力系數(shù)和力矩系數(shù)進行預測。實驗結果表明高斯過程回歸預測精度最高，誤差控制在0.24%～0.36%之間，同時高斯過程回歸預測結果比神經(jīng)網(wǎng)絡更加穩(wěn)定。

2) 將基于自動核構造(AKC)算法學習得到的核函數(shù)與其他常用的核函數(shù)進行對比，通過交叉驗證法和留出法兩種數(shù)據(jù)集劃分方式展開實驗，對3種氣動性能參數(shù)進行預測。實驗結果顯示AKC-GPR模型的預測精度最高，說明AKC算法能夠不依賴人工先驗而自動地構造出最“恰當”的核函數(shù)，用于準確捕獲隱含于數(shù)據(jù)中的內(nèi)在結構特性，從而提升GPR模型的性能。

3) 通過對一個新外形導彈氣動性能進行快速預測的實例，給出了基于AKC核的高斯過程回歸模型在實際中的應用方案，相較于傳統(tǒng)的氣動性能預測方法，AKC-GPR代理模型方法可為導彈設計初期階段的氣動性能快速預測提供方案參考。