基于移動最小二乘法的氣動力數據建模方法

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

(西南科技大學 計算機科學與技術學院，四川 綿陽 621010)

0 引言

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

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

經實驗研究發現，除了支撐域半徑，權函數的選取對模型的精度也會產生很大的影響。因此，本文在移動最小二乘模型的基礎上，以模型誤差為優化目標，以支撐域半徑和權函數影響因子為優化設計變量，采用遺傳算法對其進行優化，得到最佳支撐域半徑值和權函數影響因子，最后代入氣動力模型中驗證模型結果的準確度。

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

1.1 氣動力建模簡介

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

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

1.2 移動最小二乘近似原理

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

圖1 支撐域散點示意圖

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

(1)

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

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

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

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

系數ai(x)的選取是使近似函數φ(x)在計算點x的鄰域內是待求函數f(x)的最佳近似,而擬合函數φ(x)在所有節點的誤差加權平方和為：

(2)

對J 求最小值，有：

(3)

由矩陣形式表示可得：

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

(4)

因此系數向量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)后，我們可以得出擬合函數為：

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

(9)

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

(10)

(11)

(12)

權函數的選取對擬合精度的影響很大[11]，如圖2所示,曲線越平緩，權函數的全局性越強，而曲線越陡峭，權函數的緊支撐性越強。支撐域內各樣本點對預測點的貢獻程度也會隨權函數的不同而不同，從而會影響模型準確度。

圖2 權函數的全局性

2 移動最小二乘模型優化算法

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

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

具體算法步驟如下：

1)將氣動力數據集劃分成均勻訓練集和均勻測試集，確定訓練集中的節點數N、測試集中的節點數M。

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

3)根據經驗公式確定支撐域半徑d的優化范圍[dmin,dmax]，經驗公式如下：

d=sacle×c

(13)

式中，scale是大于1的常量，可根據實驗效果進行調整，c是均勻樣本點間的距離。

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

(14)

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

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

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

算法流程圖如圖3。

圖3 算法流程圖

3 算例及其結果分析

本文忽略了雷諾數、側滑角以及升降舵的影響，采用數值計算得到的數據進行建模。對應實驗數據有：俯仰舵偏角δ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個狀態數據樣本點。

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

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

將改進模型分別用來預測訓練集和測試集中的數據點，并對比經驗模型得到的預測結果對比結果如表1所示，定義當前模型預測誤差計算公式如下：

(15)

表1 模型優化前后預測誤差對比

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

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

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

圖4 移動最小二乘法建模結果

圖5 偏最小二乘法建模結果

圖6 移動最小二乘法預測建模結果及誤差

圖7 偏最小二乘法預測建模結果及誤差

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

(16)

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

圖8 四變量輸入情況下的建模預測結果

4 結束語

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