葉文斌,郝 佳,2+,龍 輝,何仕培,王國新,滿佳寧

(1.北京理工大學 機械與車輛學院,北京 100086;2.北京理工大學長三角研究院,浙江 嘉興 314019;3.航天科工集團貴州航天技術研究院總體設計部, 貴陽 貴州 550081)

0 引言

臨近空間高超聲速飛行器(Near Space Hypersonic Vehicle,NSHV)指飛行高度在20 km～100 km,飛行速度在5 Ma～25 Ma的一類飛行器,主要包括空天飛機、臨空無人飛行器、超高速巡航導彈等[1-2]。為了應對臨近空間不同武器的威脅,發(fā)展和建立新型武器防御系統(tǒng)是必要之需,其中地基攔截是一種有效的遏制手段[3]。然而臨近空間空氣稀薄,攔截器的氣動力不足,傳統(tǒng)舵面控制方法無法滿足攔截器機動過載要求。為了克服傳統(tǒng)舵面控制方法的不足、實現(xiàn)控制指令的快速響應,直接力/氣動力復合控制是目前較理想的方法[4]。該方法使攔截器控制系統(tǒng)具有更快的響應速度和過載能力,有效提高了打擊精度,大幅提升攔截器的機動性能[5-6]。

側向直接力噴流與高超聲速來流相互干擾的流動特征非常復雜,存在激波與激波干擾、激波與附面層干擾等復雜流動,對攔截彈產(chǎn)生嚴重的力和力矩干擾,呈現(xiàn)極強非線性特性[7]。計算流體力學(Computational Fluid Dynamics,CFD)是預測側向噴流干擾效應的主要手段之一,但是噴流干擾下流動特征評估方法對CFD計算模型與網(wǎng)格有著極高要求,導致CFD計算時間急劇提升,大大延長臨近空間高超聲速攔截器研制周期。因此,如何快速、精確地構建臨近空間高超聲速攔截器側向噴流干擾氣動模型,具有十分重要的學術與工程價值[8]。

臨近空間高超聲速攔截器側向噴流干擾氣動模型具有以下特征:高維、強非線性、昂貴[9]。這些特性導致了巨大的設計空間需要被探索,設計人員必須花費大量時間和成本進行計算機仿真,才能找到最佳的設計方案。為了快速得到方案,有效方法之一是構建代理模型(也被稱為元模型、響應曲面模型或模擬器)來代替計算機仿真,然后應用優(yōu)化技術來獲得最佳設計方案[11-12]。使用代理模型可以避免大量調(diào)用耗時、昂貴的數(shù)值分析模型,從而降低計算成本。

目前,國內(nèi)外關于代理模型研究較多,在高超聲速攔截器氣動性能方面也有相關研究。早期研究人員對攔截器氣動性能評估問題進行簡化處理(降低馬赫數(shù)、無噴流干擾、降低問題維度等),構建該條件下的氣動性能代理模型。DOWELL等[15]和SILVA等[16]提出了基于CFD數(shù)值模擬構造攔截器流場降階代理模型(Reduced-Order Model,ROM)的思想,實驗證明該模型可以快速、準確地預測氣動彈性問題中非定常流場的氣動力參數(shù);AHMED等[17]以氣動外形和氣動性能(阻力和氣動熱)之間的影響關系為研究對象,利用Kriging方法構建映射模型、NSGA-II算法構建優(yōu)化模型,完成鈍頭攔截器頭部優(yōu)化設計;CARPENTER等[18]使用神經(jīng)網(wǎng)絡對攔截器氣動性能參數(shù)進行快速預測,針對不同馬赫數(shù)下的性能參數(shù),訓練不同的神經(jīng)網(wǎng)絡。實驗證明在低馬赫數(shù)情況下,實驗方法可以保證攔截器氣動特性的快速準確評估。但是隨著建模技術的發(fā)展,不少學者開始研究完整工況下,受噴流干擾影響的氣動性能代理模型構建方法。歐敏[19]對逆向噴流下的再入攔截器進行參數(shù)化分析,對主要影響因素構建二次響應面模型,在其基礎上進行減阻桿的優(yōu)化設計;王毅等[20]使用Kriging插值模型進行了帶噴流構型的氣動阻力優(yōu)化,采用拉丁超立方采樣方法來獲得訓練數(shù)據(jù),并在此基礎上進行氣動方案優(yōu)化;LEE等[21]結合普通克里金與協(xié)同克里金方法,研究側噴干擾效應下氣動性能代理模型構建方法,通過實際氣動數(shù)據(jù)集進行驗證,證明了方法的有效性。

從以上應用可以看出,基于代理模型的氣動性能快速預測已經(jīng)成為攔截器設計領域快速發(fā)展的前沿方向。國內(nèi)外對氣動性能代理模型的構建方法進行了一定的研究和實驗,但基本都是在樣本數(shù)據(jù)充足的條件下開展的。針對噴流干擾下的攔截器設計問題來說,由于計算機仿真的耗時性/昂貴性,可用數(shù)據(jù)量顯著減少。在小樣本情況下,樣本缺失了部分工況和氣動力之間的映射規(guī)律,無論如何選擇超參數(shù),得到的預測誤差始終較高[22]。因此,找到一種小樣本下代理模型構建方法,對高超聲速攔截器性能的精準預測尤為重要。

考慮到設計人員在設計攔截器時,往往可以總結出氣動參數(shù)與氣動性能之間的規(guī)律,本文提出一種融合工程經(jīng)驗與小樣本數(shù)據(jù)來構建代理模型的方法(Experience-informed Bayesian Neural Network,EBNN)。針對攔截器氣動性能參數(shù)之間非線性關系較強且工程經(jīng)驗存在主觀性、不確定性的特點,本文采用貝葉斯神經(jīng)網(wǎng)絡作為代理模型構建方法,既對非線性關系有較好的擬合能力,又能保留模型不確定性的表達能力。通過某高超聲速攔截器具體案例驗證了所提方法的有效性,分別在有/無工程經(jīng)驗條件下進行代理模型構建,研究工程經(jīng)驗對代理模型預測精度的提升效果。

1 高超聲速攔截器氣動性能評估分析

1.1 攔截器氣動性能評估模型

氣動參數(shù)是高超聲速攔截器最重要的設計參數(shù)。不同外形在不同工況下具有不同的氣動性能,如圖1所示。高超聲速攔截器設計過程中需要確定攔截器噴流控制方式。根據(jù)側向發(fā)動機在高超聲速攔截器上安裝位置的不同,可以將噴流的控制方式分為姿態(tài)控制、軌道控制以及結合前兩者的姿態(tài)軌道復合控制。在確定噴流控制方法后,還需要構建噴流干擾的數(shù)學參數(shù)模型。噴流干擾效應有兩種數(shù)學描述方法:①將噴流干擾效應折算到對發(fā)動機直接力、力矩的影響,按照干擾因子形式描述[25],干擾因子物理意義清晰,但該方法針對姿控、軌控、姿軌控耦合3種情況需要分別進行建模,使得氣動建模以及攔截器動力學建模復雜程度增加;②將噴流干擾效應折算到對氣動力、力矩的影響,該方法可以針對不同噴流干擾情況采用統(tǒng)一的數(shù)學描述,且有利于控制系統(tǒng)設計。因此,本文采用第②種方法來構建高超聲速攔截器側向噴流干擾效應數(shù)學模型。

高超聲速攔截器氣動性能評估涉及到的參數(shù)可以分為3類:第1類:外形參數(shù)。由于本文針對某一具體高超聲速攔截器作為研究對象,因此在性能評估實驗中將外形參數(shù)設定為不變量。第2類:工況狀態(tài)參數(shù),主要有馬赫數(shù)Ma、攻角α,高度h,滾轉角φ、舵偏角DEL、姿控噴流壓力比pzj、姿軌控噴流壓力比pzgj等,如表1所示。第3類:氣動特性數(shù)據(jù),主要參數(shù)有全彈法向力系數(shù)cy1、全彈俯仰力矩系數(shù)mz1、舵偏引起的法向力系數(shù)Δcy1(δ)、舵偏引起的俯仰力矩系數(shù)Δmz1(δ)四個氣動基礎量,以及噴流引起的全彈法向力系數(shù)增量Δcy1(j)、噴流引起的全彈俯仰力矩增量Δmz1(j)、噴流對舵偏引起的法向力系數(shù)影響量Δ(Δcy1(δ))j、噴流對舵偏引起的俯仰力矩系數(shù)影響量Δ(Δmz1(δ))j等四個氣動干擾量,如表2所示。

表1 工況狀態(tài)變量表

表2 氣動特性變量表

氣動干擾量需要根據(jù)氣動性能進行計算,計算方式具體如下:

噴流引起的全彈法向力系數(shù)增量

Δcy1(j)=cy1(j)-cy1。

(1)

式中:cy1(j)表示噴流存在時全彈法向力系數(shù);cy1表示無噴流時全彈法向力系數(shù)。

噴流引起的全彈俯仰力矩增量

Δmz1(j)=mz1(j)-mz1。

(2)

式中:mz1(j)表示噴流存在時全彈俯仰力矩系數(shù);mz1表示無噴流時全彈俯仰力矩系數(shù)。

噴流對舵偏引起的法向力系數(shù)影響量

Δ(Δcy1(δ))j=Δcy1(δ)j-Δcy1(δ)。

(3)

式中:Δcy1(δ)j表示噴流存在時舵偏引起的法向力系數(shù);Δcy1(δ)表示無噴流時舵偏引起的法向力系數(shù)。

噴流對舵偏引起的俯仰力矩系數(shù)影響量

Δ(Δmz1(δ))j=Δmz1(δ)j-Δmz1(δ)。

(4)

式中:Δmz1(δ)j表示噴流存在時舵偏引起的俯仰力矩系數(shù);Δmz1(δ)表示無噴流時舵偏引起的俯仰力矩系數(shù)。

1.2 攔截器氣動性能評估問題

氣動特性快速評估是高超聲速攔截器外形設計的關鍵一環(huán),但是由于噴流氣動干擾量受較多因素的影響(Δcy1(j)、Δmz1(j)是馬赫數(shù)、攻角、滾轉角、高度、姿控開度、軌控開度等6個變量的函數(shù);Δ(Δcy1(δ))j、Δ(Δmz1(δ))j是馬赫數(shù)、攻角、滾轉角、開度、姿控壓比、軌控壓比、舵偏角等7個變量的函數(shù)),建模維度高,仿真耗時長[26]。攔截器外形設計過程中,往往只能進行有限次的仿真(例如性能評估過程中,姿/軌控開度只有5水平:0-25-50-75-100),需要耗時270天以上(20臺128核高性能系統(tǒng)并行計算)。在有限次仿真條件下,設計空間存在大量未知區(qū)域。針對仿真耗時的問題,部分研究人員提出用代理模型代替仿真軟件,以達到快速評估的效果。但就超高聲速攔截器而言,數(shù)據(jù)少、成本高,已有數(shù)據(jù)無法滿足高精度代理模型的需求。

針對上述問題,本文提出基于工程經(jīng)驗與數(shù)據(jù)融合的高超聲速攔截器氣動性能代理模型構建方法,該方法將工況狀態(tài)作為輸入?yún)?shù),氣動性能作為輸出參數(shù),通過挖掘工程經(jīng)驗中蘊含的規(guī)律信息,補充小數(shù)據(jù)情況下缺失的映射規(guī)律,構建高精度代理模型,進而幫助設計人員快速探索設計空間,構建優(yōu)化方案,實現(xiàn)高超聲速攔截器快速設計的目標。

2 基于仿真數(shù)據(jù)和工程經(jīng)驗的代理模型構建

本文建立了一種融合工程經(jīng)驗的貝葉斯神經(jīng)網(wǎng)絡,該方法利用設計人員積累的設計經(jīng)驗,補充小樣本數(shù)據(jù)缺失的映射規(guī)律。考慮到高超聲速攔截器問題維度高、非線性強的特點,提出用貝葉斯神經(jīng)網(wǎng)絡作為模型框架,一方面利用神經(jīng)網(wǎng)絡對于非線性問題表征能力強的優(yōu)勢,另一方面結合貝葉斯方法增強其不確定性表征能力。

2.1 貝葉斯神經(jīng)網(wǎng)絡模型

傳統(tǒng)神經(jīng)網(wǎng)絡是點估計模型,通過單點輸入,預測對應的單點輸出。如式(5)所示,其訓練方式可以視為通過極大似然估計(Maximum Likelihood Estimation, MLE)來實現(xiàn)網(wǎng)絡參數(shù)(權重/偏置)修正。

(5)

式中:D表示訓練數(shù)據(jù)集;xi表示訓練數(shù)據(jù)的輸入部分;yi表示訓練數(shù)據(jù)的輸出部分;w表示神經(jīng)網(wǎng)絡的權重。

貝葉斯神經(jīng)網(wǎng)絡通過貝葉斯方法以及KL(Kullback-Leibler)散度來構造損失函數(shù),進而對神經(jīng)網(wǎng)絡參數(shù)進行訓練(如圖2)。與傳統(tǒng)點估計的神經(jīng)網(wǎng)絡不同,貝葉斯神經(jīng)網(wǎng)絡是一種概率估計模型,其權重和偏置是一種概率分布。通過這種形式,貝葉斯神經(jīng)網(wǎng)絡不僅可以給出預測結果,還能給出結果的置信度與置信區(qū)間。圖2上半部分所示為貝葉斯神經(jīng)網(wǎng)絡神經(jīng)元節(jié)點的計算方式,先引入先驗假設,即假設模型參數(shù)(權重/偏置)服從某種分布(一般設為標準正態(tài)分布),然后通過采樣得到具體的權重/偏置數(shù)值,最后根據(jù)神經(jīng)元節(jié)點的輸入、權重/偏置與激活函數(shù)fact,計算得到該神經(jīng)元節(jié)點的輸出。

如圖2下半部分所示,貝葉斯神經(jīng)網(wǎng)絡的前向過程與傳統(tǒng)神經(jīng)網(wǎng)絡類似,但是引入先驗概率后,神經(jīng)網(wǎng)絡參數(shù)的估計方法從極大似然估計變成了最大后驗估計(Maximum a Posteriori Probability estimate,MAP),如式(6)所示。

wMAP=argmaxlogP(D|w)=

argmaxlogP(D|w)+logP(w)。

(6)

同時貝葉斯神經(jīng)網(wǎng)絡的參數(shù)是概率分布的形式,因此不局限于argmax值。根據(jù)貝葉斯法則,神經(jīng)網(wǎng)絡模型參數(shù)θ的后驗分布為:

(7)

式中:p(D|θ)為參數(shù)θ的似然函數(shù);p(θ)為參數(shù)θ的先驗概率;p(D)為一個歸一化常數(shù);D為訓練數(shù)據(jù)集。

2.2 基于貝葉斯神經(jīng)網(wǎng)絡的攔截器代理模型構建過程

本節(jié)以高超聲速攔截器為建模對象,結合仿真數(shù)據(jù)和工程經(jīng)驗,通過貝葉斯神經(jīng)網(wǎng)絡構建高精度代理模型,構建流程如圖3所示。

建模流程主要由兩部分組成,下面對兩部分內(nèi)容進行詳細介紹。

2.3 仿真數(shù)據(jù)/經(jīng)驗獲取

仿真數(shù)據(jù)的獲取主要依靠仿真軟件,其獲取過程包括網(wǎng)格劃分、數(shù)值求解和后處理3部分。首先構建攔截器模型文件,將其導入到網(wǎng)格生成軟件中,對計算域進行網(wǎng)格劃分;然后采用國家數(shù)值風洞(National Numerical Wind tunnel,NNW)工程套裝軟件FlowStar獲取臨近空間攔截器側向噴流干擾氣動特性;最后根據(jù)CFD仿真云圖計算無噴流攔截器氣動特性與有噴流攔截器氣動特性

工程經(jīng)驗的獲取則主要依靠設計人員。本文將工程經(jīng)驗定義為:工況變量與氣動性能之間的關聯(lián)關系?？紤]到貝塞爾曲線能夠準確量化描述變量之間的變化規(guī)律,本文提出基于貝塞爾曲線的工程經(jīng)驗的獲取方法[27],為了方便計算,以梯度的形式表征工程經(jīng)驗。

2.3.1 基于貝塞爾曲線的工程經(jīng)驗獲取方法

貝塞爾曲線的原理是伯恩斯坦多項式(由伯恩斯坦多項式可證明:在[a,b] 區(qū)間上所有的連續(xù)函數(shù)都可以用多項式來逼近),因此貝塞爾曲線可在指定區(qū)間內(nèi)表征任意的連續(xù)函數(shù)。以三階貝塞爾曲線為例說明其基本原理。

(8)

(9)

(10)

(11)

進一步推導:

(12)

(13)

B(t)=P0(1-t)3+3P1(1-t)2t+3P2(1-t)t2+P3t3,
t∈[0,1]。

(14)

結合貝塞爾曲線的表達式可以發(fā)現(xiàn),擬合曲線由一系列控制點Pi所定義,曲線上點的數(shù)量由t所決定。本文利用曲線由控制點定義的特點,開發(fā)了由專家操作的交互式知識獲取工具。利用知識獲取工具,專家通過增減控制點的數(shù)量、拖動控制點的位置可得到擬合曲線。專家完成知識繪制后,知識獲取工具保存控制點坐標,完成知識獲取整個過程。如圖5所示為知識獲取工具界面,初始曲線如圖5a所示,依次拖動控制點2至控制點7,擬合曲線由圖5a改變成圖5b。

2.3.2 工程經(jīng)驗梯度化表征

由于獲取到的工程經(jīng)驗只表示不同參數(shù)間的變化趨勢,其變化區(qū)間沒有實際物理意義,本文采用梯度的形式表征其變化趨勢,并基于梯度構建貝葉斯神經(jīng)網(wǎng)絡損失函數(shù)的經(jīng)驗項。根據(jù)貝塞爾曲線解析式,貝塞爾曲線的梯度信息就是擬合點縱坐標關于橫坐標的導數(shù)信息,因此本文根據(jù)鏈式求導法則,將t作為中間變量,計算擬合曲線上點B(t)的梯度gexp,如式(15)所示:

(15)

梯度化表征后的工程經(jīng)驗具體如表3所示,其余氣動性能參數(shù)梯度化工程經(jīng)驗請見附錄。

表3 全彈法向力工程經(jīng)驗表

2.4 代理模型訓練

本文以貝葉斯神經(jīng)網(wǎng)絡作為攔截器代理模型構建方法,首先對貝葉斯神經(jīng)網(wǎng)絡參數(shù)進行初始化,包括神經(jīng)網(wǎng)絡層數(shù)、每層節(jié)點數(shù)、激活函數(shù)、權重、偏置的概率分布等。然后根據(jù)權重/偏置的概率分布進行采樣,獲取實際權重/偏置值。之后利用仿真數(shù)據(jù)與工程經(jīng)驗分別構建數(shù)據(jù)項損失函數(shù)和經(jīng)驗項損失函數(shù),結合二者構建最終損失函數(shù)。最后基于損失函數(shù)計算誤差,通過誤差反向傳播訓練貝葉斯神經(jīng)網(wǎng)絡的模型參數(shù)。

2.4.1 參數(shù)初始化

構建貝葉斯神經(jīng)網(wǎng)絡第一步是模型參數(shù)初始化,其中神經(jīng)網(wǎng)絡層數(shù)、每層節(jié)點數(shù)、激活函數(shù)設置方式與傳統(tǒng)神經(jīng)網(wǎng)絡相同。然而貝葉斯神經(jīng)網(wǎng)絡還是一種概率估計模型,初始化時需要設定網(wǎng)絡權重及偏置的概率模型,一般將其設置為高斯分布w/b～N(μ,σ2)。但如果貝葉斯神經(jīng)網(wǎng)絡的權重、偏置直接從高斯分布中采樣獲得,會使得μ和σ變得不可微。目前解決該問題的方法是對其進行重參數(shù)化操作,以權重為例,將w重寫為:

w=μ+σε。

(16)

式中:w表示神經(jīng)網(wǎng)絡的權重;μ表示該權重概率分布的均值;σ表示該權重概率分布的方差;ε服從標準正態(tài)分布,ε～N(0,I)。偏置b獲取公式同理。

在重參數(shù)化中,ε到w只涉及到線性操作(平移縮放),采樣操作在神經(jīng)網(wǎng)絡計算圖之外,對于神經(jīng)網(wǎng)絡,ε只是一個常數(shù),不會影響損失函數(shù)對μ和σ的求導。這樣便可以先從標準高斯分布采樣出隨機量,然后可導地引μ,和σ,通過這種方法,完成貝葉斯神經(jīng)網(wǎng)絡的初始化操作。

2.4.2 參數(shù)設置

由于本文采用誤差反向傳播對參數(shù)(μ,σ)進行優(yōu)化,而反向傳播有可能使得方差σ小于0,對σ進行特殊處理,將σ重寫為:

σ=log(1+eρ)。

(17)

優(yōu)化對象由(μ,σ)變成(μ,ρ),同時保證方差σ恒大于0。

2.4.3 構建損失函數(shù)

為了挖掘工程經(jīng)驗中蘊含的規(guī)律信息,并將其融入到貝葉斯神經(jīng)網(wǎng)絡的訓練過程,本文定義了一種新的損失函數(shù)公式,其損失函數(shù)值由數(shù)據(jù)項和經(jīng)驗項兩部分組成。損失函數(shù)值越低,表示模型訓練效果更好。損失函數(shù)公式如下:

LossFunction=lossdata+lossexp。

(18)

式中:LossFunction表示貝葉斯神經(jīng)網(wǎng)絡的損失函數(shù);lossdata表示基于仿真數(shù)據(jù)構建的損失函數(shù),稱為數(shù)據(jù)項;lossexp表示基于工程經(jīng)驗構建的損失函數(shù),稱為經(jīng)驗項。

(1)數(shù)據(jù)項

由2.1節(jié)可知,貝葉斯神經(jīng)網(wǎng)絡模型訓練的訓練過程就是尋找權重參數(shù)w的最大后驗概率MAP。直接求后驗概率p(W|D)比較困難,本文采樣的是變分推斷的方式,利用一個由一組參數(shù)θ(μ,σ)控制的高斯分布q(W|θ)來逼近后驗p(W|D),利用KL散度度量q(W|D)、p(W|D)兩個分分布之間的相似性。KL散度越小,則兩個分布越相似。通過這種方式將求后驗概率的問題轉化為求最優(yōu)θ的優(yōu)化問題。KL散度公式如下:

(19)

經(jīng)過分解轉為,原等式變?yōu)?

(20)

此時最小化式(20)就可以近似求得后驗p(W|D)。同時,因為貝葉斯神經(jīng)網(wǎng)絡每一個權重/偏置都是一個獨立的概率分布,式(20)可以進一步簡化為:

(21)

式中:q(wi|θi)表示給定正態(tài)分布的參數(shù)后權重參數(shù)的分布;p(wi)表示權重的先驗;p(D|wi)表示給定網(wǎng)絡參數(shù)后,觀測數(shù)據(jù)的似然。

(2)經(jīng)驗項

經(jīng)驗項的計算是通過自動微分方法,計算貝葉斯神經(jīng)網(wǎng)絡輸入?yún)?shù)節(jié)點(例如馬赫數(shù)、攻角)與輸出參數(shù)節(jié)點(例如全彈法向力)之間的梯度關系,并將該梯度關系與實際工程經(jīng)驗中的梯度關系(2.2節(jié))進行比較,計算均方誤差,進而得到基于經(jīng)驗的損失函數(shù)項。

損失函數(shù)經(jīng)驗項用到的自動微分是將一個復雜的數(shù)學運算過程分解為一系列簡單的基本運算,每一項基本運算都可以通過查表得出。自動微分有兩種形式前向模式 (forward mode)和反向模式 (reverse mode),本文用到的就是自動微分的反向模式。該方法首先正向遍歷整個圖,計算出每個節(jié)點的值;然后逆向(從上到下)遍歷整個圖,計算出節(jié)點的偏導值,步驟如圖6所示;節(jié)點內(nèi)藍色圓圈的數(shù)值表示正向計算出的結果,為了方便表達,從下到上、從左到右依次標注為n1到n7,可以看到最后的值n7(頂部節(jié)點)為42。

在逆向求導過程中使用鏈式求導方法:

(22)

(23)

通過自動微分可以快速獲得輸出變量和輸入變量之間的微分。將該方法應用于神經(jīng)網(wǎng)絡,可以快速獲得輸出層任意節(jié)點與輸入層任意節(jié)點之間的微分。因此,可以得到知識項的損失函數(shù)如式(25)所示。

(24)

lossexp=(gmodel-gexp)2。

(25)

其中:gmodel表示貝葉斯神經(jīng)網(wǎng)絡的輸出節(jié)點y到輸入節(jié)點x的偏導數(shù);gexp表示輸出節(jié)點y與輸入節(jié)點x之間的工程經(jīng)驗(梯度關系)。

2.4.4 參數(shù)優(yōu)化

貝葉斯神經(jīng)網(wǎng)絡訓練方法與普通神經(jīng)網(wǎng)絡總體相同,首先計算損失函數(shù),然后計算權重/偏置關于損失函數(shù)的導數(shù)。與普通神經(jīng)網(wǎng)絡不同的是,其訓練對象是權重和偏置的概率分布參數(shù)N(μ,σ2),而非權重和偏置的數(shù)值。計算過程如下。

首先計算概率分布的期望μ關于損失函數(shù)的梯度:

(26)

然后計算概率分布的方差σ關于損失函數(shù)的梯度:

(27)

最后對概率分布參數(shù)進行更新:

μ←μ-αΔμ,

(28)

ρ←ρ-αΔρ。

(29)

其中:α表示學習率;μ表示對應權重/偏置的期望;ρ表示對應權重/偏置的方差,∈表示標準正太分布。

重復計算過程,直到模型訓練誤差或者訓練代數(shù)達到指定要求。

4 實驗設置與結果分析

本文針對某高超聲速攔截器氣動性能預測問題,以融合工程經(jīng)驗的貝葉斯神經(jīng)網(wǎng)絡(Experience-informed Bayesian Neural Network,EBNN)為代理模型構建方法,根據(jù)仿真結果對代理方法的氣動性能預測結果進行誤差分析。實驗的輸入?yún)?shù)為攔截器工況狀態(tài),如表1所示,輸出參數(shù)為攔截器氣動性能,如表2所示。

本文分別對有噴流影響與無噴流影響情況開展實驗,研究工程經(jīng)驗融入對模型預測精度的影響。本實驗采用的模型參數(shù)設置如表4所示。

表4 算法參數(shù)設置表

4.1 實驗指標

(1)均方根誤差

(30)

式中:m為預測值的數(shù)量;yi為第i個預測點的真實值;y′i為第i個預測點的預測值。

均方根誤差(RMSE)又叫標準誤差,它表示預測值和觀測值之間差異(稱為殘差)的樣本標準差。一般來說,RMSE越小表示模型精度越高。

(2)平均絕對百分比誤差

(31)

式中:m為預測值的數(shù)量,yi表示第i個預測點的真實值,y′i表示第i個預測點的預測值。

平均絕對誤差表示(MAPE)是相對誤差度量值。一般來說,MAPE越小,表示模型精度越高。

(3)模型擬合度

(32)

式中:m為預測值的數(shù)量;yi為第i個預測點的真實值;y′i為第i個預測點的預測值。

模型擬合度(R-Squared)表示模型擬合效果。若結果為0,說明模型擬合效果很差;若結果為 1,說明模型無誤差。一般來說,R-Squared 越大,表示模型擬合效果越好。

4.2 實驗結果分析

4.2.1 工程經(jīng)驗對代理模型的影響分析

為了驗證本文所提方法的有效性,針對氣動性能參數(shù)展開對比實驗,對有/無噴流干擾情況下進行測試。首先通過仿真軟件獲取訓練數(shù)據(jù),共300組數(shù)據(jù)。本文使用留出法來對代理模型進行訓練,將數(shù)據(jù)分為兩組互斥數(shù)據(jù),其中200組用來訓練,100組用做測試。圖7展示了不同參數(shù)的訓練過程。

在圖7中,曲線表示30 000次訓練過程中模型預測誤差的變化情況。圖7a～圖7d四個子圖表示無噴流干擾下的EBNN訓練情況,圖7e～圖7h四個子圖表示有噴流干擾下的EBNN訓練情況。其中,紫色曲線表示沒有融合工程經(jīng)驗的訓練過程,黃色表示融入工程經(jīng)驗的訓練過程。從圖7可以看到,對于不同氣動性能參數(shù)來說,黃色曲線大多低于紫色曲線,這說明融入工程經(jīng)驗的模型誤差始終低于無工程經(jīng)驗的模型誤差。在同樣迭代次數(shù)下,融合工程經(jīng)驗的代理模型具有更高的模型精度。同時,可以看到在引入噴流干擾因素后,EBNN訓練出現(xiàn)明顯的震蕩,但是融入工程經(jīng)驗的預測誤差仍然低于無工程經(jīng)驗融入的預測誤差。綜上,可以得出以下結論:融入工程經(jīng)驗能降低EBNN的訓練誤差。

除了訓練誤差,本文也對模型泛化性能展開了實驗,實驗方式是利用300組數(shù)據(jù)中的100組測試數(shù)據(jù)測試訓練完成的代理模型,進而得到測試誤差,并以絕對百分比誤差作為測試誤差指標,如圖8所示。圖8中黃色部分表示在訓練過程中沒有加入經(jīng)驗的代理模型的測試誤差,綠色部分表示訓練過程中加入經(jīng)驗的代理模型的測試誤差。從圖8中可以看出在加入了工程經(jīng)驗后,不同氣動性能的測試誤差都有一定幅度的降低,例如cy1(全彈法向力系數(shù))由16.58%降低到11.23%。由此,可以得出結論:融入工程經(jīng)驗也可以降低EBNN的泛化誤差,提高其泛化性能。

為了進一步研究融入工程經(jīng)驗對代理模型的影響,本文以有/無噴流干擾和有/無工程經(jīng)驗融入作為標準進行分類,以均方根誤差(RSME)為指標,通過100組數(shù)據(jù)對訓練好的代理模型進行誤差測試。測試結果如圖9所示。

由圖9可知,在有噴流干擾情況下,無論是有/無經(jīng)驗融入,代理模型整體誤差都大于無噴流干擾情況,均方誤差分別從7.06提升到21.67,5.03提升到14.09。這說明噴流干擾會給攔截器氣動性能預測模型訓練帶來較大的困難,因為在噴流干擾情況下,設計參數(shù)增加了姿控開度和軌控開度參數(shù),數(shù)據(jù)維度上升,進一步加重數(shù)據(jù)稀缺問題,提高了代理模型的測試誤差。同時,在融入工程經(jīng)驗后,無噴流干擾下,代理模型預測誤差從7.06降到5.03,降幅28.75%;有噴流干擾下,代理模型的預測誤差從21.67降到14.09,降幅34.98%。這說明樣本數(shù)據(jù)稀缺問題越嚴重,融入工程經(jīng)驗帶來的誤差降幅效果越好。

綜上,可以得到如下結論:①噴流干擾影響比較大,引入噴流干擾因素會增加數(shù)據(jù)維度,使得數(shù)據(jù)缺失問題更嚴重,進而增加代理模型擬合難度;②融入工程經(jīng)驗可以降低代理模型的測試誤差,提高模型訓練效果;③工程經(jīng)驗對信息缺失更嚴重的代理模型提升效果更好。在有噴流干擾情況,經(jīng)驗的融入對模型精度有著更好的效果,預測誤差降低幅度大于無噴流干擾的誤差降低幅度。

4.2.2 不同代理模型的預測結果對比分析

為了進一步說明方案的優(yōu)越性,本文還將EBNN與常用的代理模型方法進行比較,包括人工神經(jīng)網(wǎng)絡方法(Artificial Neural Network, ANN)、線性插值方法(Linear Interpolation, LI)、克里金插值方法(Kriging, KG),以RMSE與R2為指標。為了降低隨機性的影響,代理模型訓練方式采用五折交叉驗證(5-fold cross validation),并重復10次取平均結果,比較結果如圖10和表5所示。

表5 不同算法模型擬合度(R2)

從圖10中可以看出,針對cy1、mz1、Δcy1(δ)、Δmz1(δ),最低模型擬合度為0.965(LI)。這說明在無噴流干擾情況下,4種代理模型方法都能保證較高的模型擬合度。而受到噴流干擾時,代理模型的模型擬合度(R2)普遍有降低趨勢,例如全彈法向力系數(shù)cy1到噴流引起的法向力系數(shù)增量cy1(j),R2從0.99分別降低到了0.923(ANN),0.862(LI)、KG(0.887)和0.956(EBNN)。同時在表5中還可以看到,在有噴流干擾時,EBNN的R2始終能保證較高的水準,最低為0.775。其他算法則有較大的起伏波動,最低分別為0.692(ANN),0.561(LI)和0.681(KG)。這說明EBNN相比于其他算法,在有/無噴流干擾情況中都能有較高的模型擬合度。

除了使用交叉驗證法預測不同代理模型的預測精度和模型擬合度,本文還使用留出法來對不同代理模型進行訓練,一部分作為訓練集,另一部分用作測試集。圖11展示了不同代理模型在不同劃分比下的預測誤差。橫坐標表示訓練集的數(shù)據(jù)比例,縱坐標表示所有氣動性能預測誤差的均值,對每種劃分比例的實驗重復進行10次,取平均值。從圖11中可以看到,隨著訓練集數(shù)據(jù)占比的增加,4種代理模型的預測誤差都呈下降趨勢,且在所有數(shù)據(jù)劃分比下,EBNN都有著最低的預測誤差。由此可見,EBNN相較于其他代理模型方法,確實能夠更準確地預測攔截器氣動性能。

除了預測精度,本文還以誤差均值和方差為依據(jù)進行了誤差分析,以驗證方法的魯棒性,結果如圖12所示。同樣是五折交叉驗證法,先是記錄每種代理模型對所有氣動性能參數(shù)的預測誤差,然后計算誤差的均值和方差。橫軸表示用到的代理模型,左邊縱軸表示代理模型的8種氣動性能預測誤差均值,右邊縱軸表示預測誤差的方差。從圖12可以看出,面對不同氣動性能參數(shù),在4種代理模型中,EBNN的方差最小。由此可見EBNN的魯棒性優(yōu)于其他代理模型算法,其性能最穩(wěn)定?？紤]到攔截器氣動力復雜多變,因此筆者認為:在噴流與非噴流條件都表現(xiàn)較好的EBNN算法更適合攔截器氣動外形設計問題。

5 結束語

本文針對噴流干擾下的攔截器氣動性能快速評估問題,提出一種基于貝葉斯神經(jīng)網(wǎng)絡的代理模型構建方法。首先分析攔截器飛行過程,從中總結出重要設計參數(shù)。然后針對傳統(tǒng)設計過程中的不足,提出基于代理模型的快速設計方法,并用融合工程經(jīng)驗的代理模型代替原本耗時昂貴的仿真軟件。以某高超聲速攔截器為對象,基于設計過程中常用的代理模型構建方法為對比,以均方誤差和模型擬合度為驗證指標展開試驗。得到結論如下:

(1)本文所提出的代理模型構建方法在某高超聲速攔截器上有著良好表現(xiàn)。通過融合工程經(jīng)驗,可以有效降低代理模型的預測誤差。無論是訓練過程,還是測試過程,有經(jīng)驗融合的代理模型都表現(xiàn)出更為優(yōu)秀的預測精度。目前只用到了梯度形式的工程經(jīng)驗,在后續(xù)工作中,可以考慮挖掘分析更多類型的工程經(jīng)驗,進一步提供代理模型的預測精度。

(2)從實驗結果中可以看出,噴流干擾因素將會影響代理模型的預測精度。在考慮舵偏影響情況下,代理模型預測誤差會有小幅提升。而在考慮噴流干擾因素時,代理模型的預測誤差會進一步提升。但是本文所提出的融合工程經(jīng)驗的代理模型構建方法對噴流干擾影響下的性能預測有著更加優(yōu)秀的效果。相較于無噴流干擾情況,在融合工程經(jīng)驗后,最終訓練得到的代理模型的測試誤差在5%以內(nèi),表明本文所建立的代理模型能夠獲得較高精度的氣動性能參數(shù)預測結果,并具有良好的泛化能力。

(3)本文所提出的代理模型構建方法在魯棒性上也有更為優(yōu)秀的表現(xiàn)。相較于傳統(tǒng)代理模型構建方法,本文提出的融合工程經(jīng)驗的貝葉斯神經(jīng)網(wǎng)絡方法(EBNN)有著更高的魯棒性。通過殘差分析,我們發(fā)現(xiàn)對于不同的氣動性能參數(shù),EBNN的異常誤差更少。并且在不考慮異常誤差時,EBNN依舊有著更低的測試誤差。

除了以上結論,作者在實驗過程中還遇到了一些研究問題,計劃在后續(xù)工作中解決。首先是貝葉斯神經(jīng)網(wǎng)絡的冷后驗問題。貝葉斯神經(jīng)網(wǎng)絡對于先驗分布有一定的要求,不同的先驗分布會影響訓練結果。標準高斯先驗通常不是最優(yōu)的,如何找到最優(yōu)的先驗分布是一個需要解決的問題。其次是經(jīng)驗類型問題。工程設計問題中存在多種多樣的設計經(jīng)驗,目前本文只用到了梯度形式的經(jīng)驗,如何統(tǒng)籌不同類型的工程經(jīng)驗,將其都融入到代理模型訓練過程中也是一個需要解決的問題。

附錄

表2 舵偏引起的法向力系數(shù)工程經(jīng)驗表

表3 舵偏引起的俯仰力矩系數(shù)工程經(jīng)驗表

表4 噴流引起的法向力系數(shù)工程增量經(jīng)驗表

表5 噴流引起的俯仰力矩系數(shù)增量經(jīng)驗表

表6 噴流對舵偏引起的法向力系數(shù)增量經(jīng)驗表

表7 噴流對舵偏引起的俯仰力矩系數(shù)增量經(jīng)驗表