周志高，黃 俊，劉志勤，黎茂鋒，李光偉

(西南科技大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽 621010)

0 引言

近年來高超聲飛行器的優(yōu)越性能使其成為各科技強(qiáng)國研究的熱點(diǎn)技術(shù)，而建立準(zhǔn)確的氣動力模型是飛行器工程中進(jìn)行飛行控制、穩(wěn)定設(shè)計以及飛行仿真的重要基礎(chǔ)和前提[1]。目前投入應(yīng)用的氣動模型主要來自于風(fēng)洞試驗(yàn)數(shù)據(jù)，通過將氣動力狀態(tài)量和控制量以氣動數(shù)據(jù)庫表格的形式存儲，這樣的數(shù)據(jù)庫容量非常大，包含范圍很廣的飛行包線，可以比較精確地模擬飛行器的整體非線性空氣動力學(xué)。這種方法通過數(shù)值計算和風(fēng)洞試驗(yàn)手段來獲得氣動數(shù)據(jù)將會耗費(fèi)大量的時間、金錢、計算資源和人力，成本大[2]。因此，對于研制成本和研制周期受限的飛行器來說，盡量減少風(fēng)洞試驗(yàn)狀態(tài)，提高數(shù)據(jù)拓展使用有效性，在最少的試驗(yàn)狀態(tài)下，建立起建立準(zhǔn)確的氣動模型顯得非常必要[3]。

本文采用具有局部擬合效果的移動最小二乘法(Moving least square method,MLSM)來完成非線性氣動力數(shù)據(jù)擬合。1981年，移動最小二乘法由Lancaster和Salkauskas提出[4]，用于數(shù)據(jù)擬合和曲面構(gòu)造，而后又有許多學(xué)者對其進(jìn)行研究改進(jìn)，其中陳美娟等[5-6]在改進(jìn)的移動最小二乘法中選取帶權(quán)重的正交函數(shù)作為基函數(shù)以及提出復(fù)變量移動最小二乘法。崔小曼等[7]在改進(jìn)的移動最小二乘法中引入Tikhonov正則化，對系數(shù)矩陣施加約束項(xiàng)從而得到精確解，避免病態(tài)方程組的形成。于成龍等[8]利用遺傳算法來求取樣本點(diǎn)最優(yōu)支撐域半徑，倪洪杰[9]則利用粒子群優(yōu)化算法來得到最優(yōu)支撐域半徑，冷亞洪[10]利用對支撐域內(nèi)抽樣點(diǎn)數(shù)尋優(yōu)來獲取最佳半徑。

經(jīng)實(shí)驗(yàn)研究發(fā)現(xiàn)，除了支撐域半徑，權(quán)函數(shù)的選取對模型的精度也會產(chǎn)生很大的影響。因此，本文在移動最小二乘模型的基礎(chǔ)上，以模型誤差為優(yōu)化目標(biāo)，以支撐域半徑和權(quán)函數(shù)影響因子為優(yōu)化設(shè)計變量，采用遺傳算法對其進(jìn)行優(yōu)化，得到最佳支撐域半徑值和權(quán)函數(shù)影響因子，最后代入氣動力模型中驗(yàn)證模型結(jié)果的準(zhǔn)確度。

1 基于MLS的氣動力建模方法

1.1 氣動力建模簡介

氣動力數(shù)據(jù)建模的目的就是利用CFD模擬數(shù)值計算的數(shù)據(jù)、動態(tài)風(fēng)洞試驗(yàn)數(shù)據(jù)建立空氣動力學(xué)數(shù)學(xué)模型和參數(shù)，飛行器氣動力和力矩系數(shù)一般是攻角數(shù)、馬赫數(shù)、飛行高度等狀態(tài)參數(shù)的函數(shù)。即Ci=f(α,Ma,h,δz,…),式中，Ci表示飛行器各氣動系數(shù)，α為攻角數(shù)，Ma為馬赫數(shù)，h為飛行高度，δz為俯仰舵偏角。

目前，飛行器的氣動力建模方法主要有兩個方面：1)基于傳統(tǒng)的物理特性建模，如多項(xiàng)式模型、三角函數(shù)模型等。2)基于人工智能等新型技術(shù)的建模方法，常見的有BP神經(jīng)網(wǎng)絡(luò)模型、Kriging模型等。

1.2 移動最小二乘近似原理

移動最小二乘法相比傳統(tǒng)的最小二乘法引入了支撐域的概念。如圖1所示，球形區(qū)域表示支撐域，評估點(diǎn)只受支撐域內(nèi)的樣本點(diǎn)的影響,并且距離越近的點(diǎn)對其影響程度越高，這種方法有效解決了分段擬合不能局部化處理的缺點(diǎn)。

圖1 支撐域散點(diǎn)示意圖

具體原理如下，假設(shè)已知一系列離散樣本點(diǎn)xi∈Ω(i=1,2,3,…,N)的函數(shù)值為f(xi)，移動最小二乘法擬合函數(shù)為φ(x)可表示為(本文原理性解釋都以一維為例)：

(1)

線性基：PT(x) = [1,x],m=2

二次基：PT(x) = [1,x,x2],m=3

三次基：PT(x) = [1,x,x2,x3],m=4

對于多維問題可以依次類推，本文將會依次對二維問題和四維問題進(jìn)行研究，隨著問題個數(shù)增加，計算量也會增加。

系數(shù)ai(x)的選取是使近似函數(shù)φ(x)在計算點(diǎn)x的鄰域內(nèi)是待求函數(shù)f(x)的最佳近似,而擬合函數(shù)φ(x)在所有節(jié)點(diǎn)的誤差加權(quán)平方和為：

(2)

對J 求最小值，有：

(3)

由矩陣形式表示可得：

A(x)a(x)=B(x)f(x)

(4)

因此系數(shù)向量a(x)為：

a(x)=A-1(x)B(x)f(x)

(5)

式中：

A(x)=PT(x)W(x)P(x)

(6)

B(x)=PT(x)W(x)

(7)

f(x)=[f(x1),f(x2),…,f(xN)]T

(8)

將式(5)代入式(1)后，我們可以得出擬合函數(shù)為：

φ(x)=PT(x)a(x)=PT(x)A-1(x)B(x)f(x)

(9)

在MLS擬合函數(shù)求解中，為了保證矩陣A可逆，要求支撐域內(nèi)至少有m個樣本點(diǎn)，可以設(shè)置合適的支撐域半徑。式(2)中，w(x-xi)稱為權(quán)函數(shù)，其只與樣本點(diǎn)和預(yù)測點(diǎn)的距離有關(guān)，即只有在預(yù)測點(diǎn)周圍某鄰域內(nèi)有值。一般常用的權(quán)函數(shù)有：高斯函數(shù)如式(10)所示，五次樣條函數(shù)如式(11)所示，指數(shù)函數(shù)如式(12)所示。

(10)

(11)

(12)

權(quán)函數(shù)的選取對擬合精度的影響很大[11]，如圖2所示,曲線越平緩，權(quán)函數(shù)的全局性越強(qiáng)，而曲線越陡峭，權(quán)函數(shù)的緊支撐性越強(qiáng)。支撐域內(nèi)各樣本點(diǎn)對預(yù)測點(diǎn)的貢獻(xiàn)程度也會隨權(quán)函數(shù)的不同而不同，從而會影響模型準(zhǔn)確度。

圖2 權(quán)函數(shù)的全局性

2 移動最小二乘模型優(yōu)化算法

在MLS近似方法中，支撐域半徑的選取以及權(quán)函數(shù)選擇直接影響其對樣本點(diǎn)的擬合精度[11-12],如果支撐域半徑太大，不能體現(xiàn)局部性，而支撐域半徑太小，則會使矩陣A不可逆或者病態(tài)。目前主要有兩種方法來確定支撐域半徑，其中一種是根據(jù)經(jīng)驗(yàn)公式來確定一個相對較優(yōu)的值；另外一種是在經(jīng)驗(yàn)公式的基礎(chǔ)上對支撐域半徑進(jìn)行優(yōu)化，進(jìn)而得到最佳支撐域半徑值。

由以上分析可知，支撐域半徑影響擬合精度的本質(zhì)是改變了支撐域內(nèi)所包含的樣本點(diǎn)數(shù)，而權(quán)函數(shù)的選擇則改變了不同樣本點(diǎn)對評估點(diǎn)的影響程度。因此，為了得到較高精度的氣動模型，本文選取高斯函數(shù)作為權(quán)函數(shù)，以高斯權(quán)函數(shù)影響因子β和支撐域半徑為優(yōu)化設(shè)計變量，通過留一法計算出模型誤差，并將其作為優(yōu)化目標(biāo)，再采用遺傳算法對優(yōu)化模型中設(shè)計變量進(jìn)行優(yōu)化搜索，直至達(dá)到結(jié)束條件。

具體算法步驟如下：

1)將氣動力數(shù)據(jù)集劃分成均勻訓(xùn)練集和均勻測試集，確定訓(xùn)練集中的節(jié)點(diǎn)數(shù)N、測試集中的節(jié)點(diǎn)數(shù)M。

2) 選取高斯權(quán)函數(shù)為權(quán)函數(shù)，并確定影響因子β的優(yōu)化范圍為[1,9],在這個范圍中可以尋求合適的β來保證局部擬合效果。

3)根據(jù)經(jīng)驗(yàn)公式確定支撐域半徑d的優(yōu)化范圍[dmin,dmax]，經(jīng)驗(yàn)公式如下：

d=sacle×c

(13)

式中，scale是大于1的常量，可根據(jù)實(shí)驗(yàn)效果進(jìn)行調(diào)整，c是均勻樣本點(diǎn)間的距離。

4)將訓(xùn)練集中的初始樣本點(diǎn)依次代入MLS模型中，通過留一法循環(huán)計算出訓(xùn)練集中各節(jié)點(diǎn)的預(yù)測值φ(xi)，同時確定優(yōu)化模型如下：

(14)

式中，e為移動最小二乘模型計算值與真實(shí)值的殘差平方和,其值越小表示MLS模型精度越高。

5)以式(14)為優(yōu)化模型，設(shè)置種群大小、迭代次數(shù)、編碼方式等參數(shù)，采用遺傳算法獲取最優(yōu)支撐域半徑d和最優(yōu)權(quán)函數(shù)影響因子β。

6)基于優(yōu)化時使用的訓(xùn)練樣本集合，將最優(yōu)支撐域半徑d和最優(yōu)權(quán)函數(shù)影響因子β代入模型中,得到所求解的氣動力模型。

算法流程圖如圖3。

圖3 算法流程圖

3 算例及其結(jié)果分析

本文忽略了雷諾數(shù)、側(cè)滑角以及升降舵的影響，采用數(shù)值計算得到的數(shù)據(jù)進(jìn)行建模。對應(yīng)實(shí)驗(yàn)數(shù)據(jù)有：俯仰舵偏角δz= -20°、-15°、-10°、-5°、0°、5°、10°、15°、20°，飛行高度h=0、10、20、30、40、50(單位：km)，Ma=2、4、6、8、10、12、14、16、18、20，α=-15°～15° (共31組)。因而共有9×5×10×31=13 950個狀態(tài)數(shù)據(jù)樣本點(diǎn)。

為了驗(yàn)證MLS模型經(jīng)優(yōu)化后模型精度的提升(記為改進(jìn)模型)，同時為了減少計算時間，首先考慮在俯仰舵偏角δz= 0°、飛行高度h=0 km情況下，對攻角數(shù)α、馬赫數(shù)Ma和氣動系數(shù)之間的關(guān)系建模。此時包含31×10=310個狀態(tài)數(shù)據(jù)點(diǎn)，將樣本數(shù)據(jù)集均勻劃分成訓(xùn)練集和測試集，保證節(jié)點(diǎn)之間的距離相等，其中訓(xùn)練集樣本用于移動最小二乘建模計算。

通過MLS模型預(yù)測一個評估點(diǎn)主要經(jīng)過兩次計算，一是依次計算評估點(diǎn)與訓(xùn)練樣本點(diǎn)的距離d，通過權(quán)函數(shù)公式可以得到權(quán)值矩陣W；二是基于訓(xùn)練樣本計算出基函數(shù)PT。最后利用式(6)～(9)可以計算出對應(yīng)評估點(diǎn)的預(yù)測值。其中由經(jīng)驗(yàn)公式選取支撐域半徑值的方式得到的MLS模型稱為經(jīng)驗(yàn)?zāi)Ｐ?，而改進(jìn)模型則是經(jīng)過尋優(yōu)算法求得最佳設(shè)計變量值。

將改進(jìn)模型分別用來預(yù)測訓(xùn)練集和測試集中的數(shù)據(jù)點(diǎn)，并對比經(jīng)驗(yàn)?zāi)Ｐ偷玫降念A(yù)測結(jié)果對比結(jié)果如表1所示，定義當(dāng)前模型預(yù)測誤差計算公式如下：

(15)

表1 模型優(yōu)化前后預(yù)測誤差對比

表1中，a、b分別表示模型在測試集和訓(xùn)練集中的預(yù)測誤差，從表中可以看出，改進(jìn)模型在軸向力系數(shù)CA、法向力系數(shù)CN和俯仰力矩系數(shù)Cmz上有著更好的預(yù)測效果，并且MLS模型在測試集上也有較好的預(yù)測效果，表明了移動最小二乘法建立的模型具有一定的泛化能力。

參考文獻(xiàn)[13]中基于偏最小二乘法(PLS)的建模方法。一般來說，多項(xiàng)式模型的項(xiàng)數(shù)越多，模型項(xiàng)中信息保留度也就越大，進(jìn)而更能反應(yīng)出真實(shí)模型。為了形成對比，本文采用與移動最小二乘法中基函數(shù)相同次數(shù)的多項(xiàng)式進(jìn)行主成分分析(PCA)，通過留一法交叉驗(yàn)證確定主成分個數(shù)，得到了對比結(jié)果。

圖4和圖5給出了兩種建模方法關(guān)于3種氣動系數(shù)的建模結(jié)果對比。其中“ErrCA”表示軸向力系數(shù)在對應(yīng)樣本點(diǎn)的預(yù)測誤差值，可以看出偏最小二乘法的建模結(jié)果出現(xiàn)明顯偏差，并且移動最小二乘法的建模結(jié)果整體上優(yōu)于偏最小二乘法。圖 6和圖7是移動最小二乘法建模結(jié)果和偏最小二乘法建模結(jié)果對測試集上樣本點(diǎn)的預(yù)測，同樣可以看出偏最小二乘模型中軸向力系數(shù)的預(yù)測結(jié)果有明顯誤差。

圖4 移動最小二乘法建模結(jié)果

圖5 偏最小二乘法建模結(jié)果

圖6 移動最小二乘法預(yù)測建模結(jié)果及誤差

圖7 偏最小二乘法預(yù)測建模結(jié)果及誤差

從圖中可以看出，對訓(xùn)練集和測試集上的樣本點(diǎn)進(jìn)行預(yù)測時，移動最小二乘法建立的模型都表現(xiàn)的更準(zhǔn)確。如式(16)，移動最小二乘法的系數(shù)為節(jié)點(diǎn)x的函數(shù)，并且由式(6)可知a(x)融入了權(quán)函數(shù)帶來的局部擬合特性，因而在相同次數(shù)的多項(xiàng)式擬合中，移動最小二乘法的建模效果要優(yōu)于偏最小二乘法。偏最小二乘法擬合方程式的系數(shù)ai為常數(shù)，同時會根據(jù)主成分個數(shù)刪掉一些信息保留度小的項(xiàng)，因此在偏最小二乘法在計算量上要明顯小于移動最小二乘法。

(16)

為了得到適用于大空域、寬速域的氣動模型，以飛行高度h、俯仰舵偏角δz、攻角數(shù)α、馬赫數(shù)Ma為自變量，將所有實(shí)驗(yàn)數(shù)據(jù)分為訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集，由于參與計算的樣本點(diǎn)個數(shù)以及問題維數(shù)增加，使得優(yōu)化需要的時間過長。因此在有4個輸入情況下的氣動建模算例中略去了模型優(yōu)化步驟，直接利用經(jīng)驗(yàn)公式確定支撐域半徑及影響因子β，得到適用于飛行高度為0～50 km、飛行馬赫數(shù)在2～20 mach的氣動力模型。圖8給出了4個輸入情況下的的氣動系數(shù)建模預(yù)測結(jié)果。為了方便展示，以攻角數(shù)α為x坐標(biāo)，以馬赫數(shù)Ma為坐標(biāo)，在同樣攻角和馬赫上會存在不同飛行高度和俯仰舵偏角的響應(yīng)值，并以平面圖展示各個樣本點(diǎn)的預(yù)測誤差。從圖中可以看出，建立的氣動模型預(yù)測誤差較小，效果較優(yōu)，滿足現(xiàn)代飛行器對寬速域和大空域的需求。

圖8 四變量輸入情況下的建模預(yù)測結(jié)果

4 結(jié)束語

本文利用某飛行器多次飛行仿真的氣動力數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù)，將移動最小二乘法應(yīng)用于氣動力數(shù)據(jù)建模。考慮到移動最小二乘法的擬合精度對支撐域半徑以及權(quán)函數(shù)比較敏感，提出了一種優(yōu)化算法。首先以攻角數(shù)和馬赫數(shù)為自變量進(jìn)行建模，經(jīng)過實(shí)驗(yàn)研究分析，確定支撐域半徑和影響因子β的優(yōu)化范圍，接著通過遺傳算法獲取全局最優(yōu)設(shè)計變量。經(jīng)算例驗(yàn)證，優(yōu)化后得到的移動最小二乘模型精度有明顯提高，且該方法與偏最小二乘方法相比預(yù)測數(shù)值也更準(zhǔn)確。最后以攻角數(shù)、馬赫數(shù)、飛行高度、俯仰舵偏角為自變量進(jìn)行建模，由于計算量大，跳過優(yōu)化步驟，其預(yù)測誤差也較小，得到的模型能滿足“大空域、寬速域”的氣動模型要求。移動最小二乘法的建模效果整體上優(yōu)于偏最小二乘法，但移動最小二乘法的計算時間卻明顯長于偏最小二乘法，尚有待研究。