胡偉杰，黃增輝，劉學(xué)軍,*，呂宏強(qiáng)

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

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

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

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

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

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

1 高斯過程回歸

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

1.1 模型推導(dǎo)

高斯過程是隨機(jī)變量的集合，任意有限維隨機(jī)變量的聯(lián)合分布是高斯分布[17]。因此，直接從函數(shù)空間角度推理，指定目標(biāo)函數(shù)f(x)的先驗(yàn)為一個(gè)高斯過程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)

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

(5)

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

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

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

(6)

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

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

(7)

后驗(yàn)分布為高斯分布，均值和方差公式具體為

(8)

(9)

1.2 模型訓(xùn)練

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

(10)

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

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

(11)

(12)

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

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

2 自動(dòng)核構(gòu)造算法

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

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

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

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

2.1 簡(jiǎn)單核及其先驗(yàn)假設(shè)

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

(13)

(14)

(15)

(16)

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

(17)

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

每一類簡(jiǎn)單核都編碼了對(duì)目標(biāo)函數(shù)不同的假設(shè)。Lin核假設(shè)函數(shù)具有線性變化結(jié)構(gòu)；SE核假設(shè)函數(shù)具有局部變化結(jié)構(gòu)；Per核假設(shè)函數(shù)具有重復(fù)變化結(jié)構(gòu)；RQ核假設(shè)函數(shù)具有多尺度變化結(jié)構(gòu)；WN核假設(shè)函數(shù)具有獨(dú)立噪聲結(jié)構(gòu)。

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

2.2 構(gòu)造規(guī)則

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

①S→S+B

②S→S×B

③S→B

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

2.3 AKC算法描述

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

圖1 構(gòu)造過程示例Fig.1 Example of construction procedure

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

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

算法1 自動(dòng)核構(gòu)造Algorithm 1 Automatic kernel construction

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

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

3.1 幾何建模

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

圖2 導(dǎo)彈幾何構(gòu)型Fig.2 Geometry of missile

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

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

3.2 網(wǎng)格劃分

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

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

3.3 數(shù)值求解

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

3.4 后處理

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

表2 部分?jǐn)?shù)據(jù)集Table 2 Part of dataset

4 核函數(shù)的選擇

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

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

5 實(shí)驗(yàn)結(jié)果及分析

5.1 GPR模型非線性擬合能力測(cè)試

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

(18)

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

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

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

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

5.2 不同模型對(duì)比

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

(19)

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

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

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

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

圖8 預(yù)測(cè)結(jié)果穩(wěn)定性對(duì)比Fig.8 Comparison of stability of predicted results

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

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

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

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

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

圖9 不同劃分比例下的實(shí)驗(yàn)結(jié)果Fig.9 Experimental results at different partition ratios

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

5.4 殘差分析

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

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

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

5.5 氣動(dòng)性能快速預(yù)測(cè)應(yīng)用

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

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

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

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

6 結(jié) 論

針對(duì)傳統(tǒng)導(dǎo)彈氣動(dòng)性能預(yù)測(cè)方法的局限，本文使用高斯過程回歸代理模型對(duì)定常流場(chǎng)條件下的典型導(dǎo)彈氣動(dòng)性能進(jìn)行預(yù)測(cè)，得到結(jié)論如下：

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

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

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