高效粒子群算法研究及飛翼無人機氣動隱身優化設計

樊華羽，詹浩，程詩信，米百剛，姚會勤

(1.西北工業大學 航空學院，西安 710072) (2.清華大學 航天航空學院，北京 100084)

0 引 言

飛行器外形多目標優化設計是一個跨學科多目標的優化任務，包含了總體、氣動、結構、隱身等多個領域和專業的關鍵技術。這些因素之間相互交叉干擾，進一步增加了現代飛行器多目標優化設計的復雜性。隨著電子計算機技術的蓬勃發展和計算流體力學(CFD)、電磁計算精度的提高，數值優化算法結合先進的CFD與電磁計算，為氣動隱身一體化優化設計提供了優秀的平臺，進而針對各類型作戰飛機的氣動隱身一體化研究成為國內外關注的熱點[1-4]。

相比于常規布局的飛行器氣動隱身優化設計，飛翼布局自身的結構特性使得其在隱身性能上優于常規布局飛行器。但也因此在缺失尾翼、垂尾等結構的情況下，其氣動及結構載荷需要通過機翼來補足。國內外針對飛翼布局無人機做了大量研究。K.Hyoungjin等[5]針對一體化翼身融合體進行了詳細地氣動分析；王豪杰等[6]通過風洞實驗對比，結合CFD完成對某型無尾飛翼布局無人機的氣動力設計；在飛翼構型氣動優化方面，王榮等[7]采用高精度粘性氣動計算模型建立了飛翼布局無人機外形精度下的氣動隱身多目標優化；張樂等[8]針對雙發布局下的飛翼無人機大鼓包機身，提出了一種減小翼型前緣半徑的機身前緣類“鷹嘴”形飛翼布局優化構型；陳曦等[9]開展了飛翼布局無人機在考慮隱身迎角下的氣動隱身綜合優化。

在針對飛翼布局無人機的多目標優化中，由于需要頻繁調用目標函數，優化過程需要投入較大的時間成本和計算機資源。如何在平衡計算精度與計算效率的情況下建立一套高效的多目標優化算法，是目前亟待解決的問題。本文針對典型飛翼布局，結合基于EHVI(Expected Hyper-Volume Improvement)加點準則的高精度多目標優化算法，建立一種新的能夠在較少調用目標函數的情況下完成優化任務的高效多目標粒子群優化算法，并對某型飛翼布局無人機進行隱身氣動多目標優化設計研究。

1 高效多目標粒子群算法

1.1 多目標優化問題

在n維搜索空間S∈Rn內，定義一個目標數量為m的多目標優化向量fi(x),i=1,…,m，那么帶約束的最小化多目標優化問題的通用數學模型可以表示為

minf(x)=[f1(x),f2(x),…,fm(x)]

(1)

s.t.gj(x)≤0,j=1,…,l

hk(x)=0,k=1,…,p

式中：x=[x1,x2,…,xn]，為設計變量向量；gj(x)為第j個不等式約束；hk(x)為第k個等式約束；l和p分別為不等式約束和等式約束的數量。

1.2 多目標粒子群優化算法

粒子群優化(Particle Swarm Optimization，簡稱PSO)算法最初是為了模擬鳥群的飛行覓食行為而提出的一種群智能算法。其初始化為一群隨機粒子，通過迭代尋找最優解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自己：第一個極值是粒子個體所找到的歷史最優解，稱為個體最優值Pb；另一個極值是整個種群目前找到的最優解，稱為全局最優值Gb。在找到上述兩個最優值時，粒子根據式(2)～式(3)來更新自己的速度和位置。

(2)

(3)

式中：c1和c2分別為認知和社會加速系數，取值范圍為[0,4]，一般取c1=c2=2；r1和r2為兩個在[0,1]內服從均勻分布的隨機變量。

經過眾多研究人員的深入研究，現在粒子群算法已可以方便地用于處理多目標優化問題。

CDMOPSO(Crowding Distance Multi-objective PSO)算法是C.R.Raquel等在C.A.C.Coello等提出的多目標粒子群算法(MOPSO)[10]的基礎上，使用擁擠度算子代替超網格來維護外部檔案而改進(如圖1所示)的一種多目標優化算法[11]。擁擠度算子[12]能夠通過選擇引導粒子運動的最優解以及維護外部檔案來增加種群的多樣性，避免算法過早陷入局部最優。這種算法的優點是結構簡單，易于實現，缺陷是基于均勻分布的變異操作使其易于陷入局部最優，解決多模態優化問題的能力較差。氣動隱身多目標優化設計是一種高維強非線性問題，因此本文選擇α-stable變異函數來對其進行改進，通過α-stable分布[13-14]產生的隨機數，對PSO種群中的個體進行變異。在變異的過程中，通過動態調整函數的穩定系數α來調整不同優化階段變異的范圍和幅度。穩定系數α描述了分布的尾跡大小，決定了隨機數的范圍。α的變化范圍是[1,2]。在開始階段使用較小的α值，增強全局搜索能力，隨著優化循環的進程，α值越來越大，全局搜索能力減弱，局部搜索能力增強，為尋找高精度解提供幫助。其具體變異操作過程以及優化結果對比詳見文獻[15]，流程圖如圖2所示。本文將其命名為ASMOPSO(Alpha-stable Multi-objective PSO)。

圖1 擁擠距離計算示意圖Fig.1 Schematic diagram of congestion distance calculation

圖2 ASMOPSO流程圖Fig.2 Flow chart of ASMOPSO

1.3 基于Kriging模型的EHVI加點準則

為了增加局部地區的代理模型精度，減少對初始樣本空間精度的依賴，采用基于自適應加點的動態Kriging代理模型對粒子群中的未觀測點進行近似評估。在動態更新代理模型樣本點的過程中，使用EHVI準則選擇粒子群多目標優化過程中若干個最可靠的未觀測點進行真實函數評估，更新Pareto解集，同時提高Kriging代理模型局部區域精度。

EHVI是M.Emmerich等[16-17]在超體積理論[18]的基礎上結合單目標EI(Expected Improvement)加點準則[19]提出的一種處理高耗時優化問題的新型多目標加點準則。

給定一個新的解y，假設它不被Pareto解集P中的任意個體所支配，那么解集P的超體積改善為

H(y,P)=H(P∪{y})-H(P)

(4)

進而就可以定義超體積改善函數為

(5)

那么基于Kriging代理模型的多目標響應可視為互相獨立且服從高斯分布的隨機變量Yj(x)，有：

(6)

式中：m為目標數量。

在此基礎上，可以得到多目標期望改善函數：

(7)

式中：Vnd為由Pareto解集確定的非支配區域，即如圖3所示的細實線與坐標軸封閉的區域。

圖3 超體積改善示意圖Fig.3 Diagram of hypervolume improvement

從式(6)可以看出：如果要精確計算EHVI值，需要在非支配區域進行多維積分。由于非支配區域的不規則性，將超體積區域分割成多個矩形單元是必不可少的步驟。當Pareto解較多或者目標維度較高時，精確地識別這些矩形單元并對其進行EI積分是一件非常耗時的工作。

為此，Cheng Shixin等[20]提出了一種動態的EHVI值計算方式，通過對多個標準函數的對比測試結果分析，相比其基本型CDMOPSO優化算法，結合動態超體積期望改善的優化算法能在大幅度減少調用真實函數次數的情況下保持計算精度，極大地提高了優化效率。本文所使用的加點準則就是基于此動態計算方法的EHVI加點準則。

1.4 約束的處理

為了克服罰函數系數設定難的問題，優化算法中采用不可行度(Infeasibility Degree，簡稱IFD)法來處理優化問題中的約束[21]，當一個粒子處在不可行域內時，IFD值可以表示為粒子與可行域邊界的接近程度，則定義第i個粒子的不可行度值為

(8)

式中：γ為等式約束的容限區間，γ≥0。

如果一個粒子位于可行域內，則它的IFD值為0，這種粒子稱為可行粒子。基于不可行度和Pareto優勝比較的粒子選擇機制是：①如果一個粒子可行，而另一個粒子不可行，那么選擇可行粒子；②如果兩個粒子都不可行，則選擇IFD值較小的那個粒子；③如果兩個粒子都可行，執行Pareto優勝比較，選擇非支配粒子。

基于式(8)將每個粒子的近似不可行度改寫為

(9)

1.5 基于EHVI的多目標粒子群算法

區別于傳統的尋找當前Kriging代理模型下的最大EHVI值的子優化過程，本文計算當前種群中所有個體的EHVI值并對它們進行降序排列，選取前幾個個體進行真實函數評估并加入樣本集，用于更新外部非支配解檔案，引導ASMOPSO中個體的移動。在這個過程中Kriging代理模型只是作為提供未觀測點的均值和方差的工具，用于計算EHVI值，代理模型的精度不是首要考慮的因素。

綜上，結合多目標粒子群算法、Kriging代理模型、EHVI加點形成的高效多目標優化算法(EHVIMOPSO)的流程如圖4所示。

圖4 EHVIMOPSO流程圖Fig.4 Flow chart of EHVIMOPSO

2 氣動隱身多目標優化設計

2.1 優化模型的建立

使用飛翼無人機作為初始外形，如圖5所示，以翼面為研究對象，以阻力和RCS(Radar Cross Section)為優化目標，約束為升力系數不小于原始翼型和俯仰力矩系數絕對值不減小。設計點氣動計算狀態為：Ma=0.8,α=2.0°,Re=2.78×107。

圖5 飛翼布局無人機幾何外形Fig.5 Geometry of flywing UAV

由此可得優化的數學模型為

(10)

式中：X為設計變量向量；CD為阻力系數；ARCS為RCS均值；CL為升力系數；CM為力矩系數；上標*表示初始翼型的計算值。

2.2 計算網格及仿真計算方法

流體計算方面，流場求解基于三維非定常雷諾平均N-S方程，湍流模型采用k-ωSST模型，物面邊界采用無滑移邊界條件，遠場采用壓力遠場邊界條件。

本文使用半模計算，采用由ANSYS ICEM CFD軟件生成的空間六面體網格對流動區域進行離散，網格數量為165萬，計算網格如圖6所示。在優化設計循環中，使用該軟件自帶的腳本功能實現網格的運動[21]。

電磁隱身特性求解以RCS為衡量標準，它描述了物體因被電磁波照射而向各個方向散射，被雷達捕捉的雷達回波強度及其電磁特性。在計算目標RCS的過程中，需要平衡計算效率與計算精度。由于不同的計算方法對于不同尺寸的目標具有不同的適應性，在選擇計算方法時需要考慮目標的電尺寸大小。

圖6 計算網格Fig.6 Computational grid

綜合考慮，本文選用大面元物理光學法LEPO(Large Element PO)配合一致性幾何繞射理論 (UTD)計算邊緣繞射場。計算狀態為單站，極化方式為水平極化，雷達頻率選擇10 GHz。采用半模計算，計算范圍(θ)為0°～180°，步長為1°。RCS計算示意圖如圖7所示。

圖7 RCS計算示意Fig.7 Diagram of RCS calculation

2.3 參數化方法

在優化設計中，采用空間變形能力較強的FFD(Free Form Deformation)自由曲面變形法[22]作為參數化方法。該方法最早由T.W.Sederberg和S.R.Parry提出，是一種針對三維可變形物體的有效建模工具。其基本方法是通過構建并操縱空間三維框架與被操縱面的映射關系，通過改變空間三維框架來改變被操縱面的形狀。由于此方法能夠適用于非常復雜的外形控制且模擬精度很高，在飛行器三維參數化中有較好的應用。

飛翼優化模型中的FFD參數化[23]空間控制點為翼面上下各25個控制點(如圖8(a)所示)，在坐標系中，x、y、z三個方向上布置的控制點數為(5,5,2)，由于控制變量明顯偏多，考慮在外形優化中將整個變形框外框固定不變，選擇上下翼面中心位置的9個點(如圖8(b)所示)，一共18個點，坐標系中表示在x、y、z方向上的點數為(3,3,2)，變形方向限定為僅在z方向上。為了不讓變形過于劇烈，變形范圍設定在±0.05之內。最后控制點一共為3×3×2=18個。

(b) FFD參數化翼面控制點示意圖圖8 翼面FFD參數化控制點Fig.8 Parameterized control point of airfoil FFD

在FFD參數化后，通過曲面差值生成新的飛翼翼面。參數化結果如表1所示，可以看出：FFD參數化后的代理模型與原始外形的翼面力學特性幾乎完全吻合，隱身性能存在0.3%的誤差，符合計算中的誤差范圍，故認為該參數化方法能夠精確表達原始模型的氣動外形。

表1 參數化結果比較Table 1 Comparison of parametric results

2.4 優化過程及優化結果分析

使用最優拉丁超立方抽樣(Optimal Latin Hypercube Sampling，簡稱OLHS)[24]生成40個初始樣本點，建立初始代理模型并通過EHVIMOPSO多目標優化算法搜索Pareto前沿。其中粒子群種群規模為200，外部檔案大小為50；慣性系數從0.75逐漸減小到0.25；認知加速系數和社會加速系數c1=c2=2.0；穩定系數α的變化范圍為1.0～2.0。迭代步數為60步，每代使用EHVI值選擇3個最穩定的個體加入到樣本空間中，并更新代理模型。優化結果如圖9和表2所示。

圖9 Pareto前沿Fig.9 Forward position of Pareto 表2 Pareto前沿數值統計Table 2 Frontier point of Pareto

序號CDARCS/m210.008 5151.065 95120.008 5471.053 66030.008 5221.061 07740.008 5251.059 75450.008 5361.053 8496(A點)0.008 5321.055 57170.008 5271.057 25880.008 5301.056 98090.008 5421.053 797100.008 5451.053 757

從圖9和表2可以看出：Pareto解分布較為均勻，但是范圍還不夠寬廣。

優化前后數據比較如表3所示。由于此飛翼布局無人機為低可探測低阻外形，其原始外形的氣動及隱身性能非常優秀(如表1所示)，因此在此優化算法下的氣動優化效果沒有特別明顯的提升。

表3 優化前后數據比較Table 3 Data comparison before and after optimization

在平衡減阻與降低RCS兩者情況下選擇如圖9所示的ParetoA與原始構型進行比較分析。

從表2可以看出：在ParetoA的構型下阻力減少了0.34%，RCS縮減了4.41%。

在氣動計算結果對比方面，為了更好地說明，選取Y分別為1.0、2.5、4.5三個展向站位的翼型截面進行對比分析，如圖10所示。優化后的翼型剖面形狀、壓力系數與原始外形對比結果如圖11～圖13所示。

圖10 翼面截取示意Fig.10 Sketch of wing interception

(a) 翼型對比

(b) 壓力系數對比圖11 對比結果(Y=1.0)Fig.11 Contrast result(Y=1.0)

(a) 翼型對比

(b) 壓力系數對比圖12 對比結果(Y=2.5)Fig.12 Contrast result(Y=2.5)

(a) 翼型對比

(b) 壓力系數對比圖13 對比結果(Y=4.5)Fig.13 Contrast result(Y=4.5)

從圖11～圖13可以看出：優化構型三個站位上的翼型剖面與原始翼型相比，在上翼面前部與下翼面后部都有小量縮減，翼面最大厚度向后移動且最大厚度值都減小、變薄，其中Y=2.5、Y=4.5兩個站位翼型變薄的情況非常明顯；翼型的前后緣位置壓力分布較陡峭，中段壓力系數曲線分布較為平緩，優化后構型仍能在此基礎上有一定程度的減緩翼面激波強度；Y=1.0站位，后緣附近有一個向下的加載力區域，使的整個優化構型的低頭力矩減小，提高了無人機的可操控性。

RCS優化結果對比如圖14所示，可以看出：在入射角為57°和140°處有峰值存在，57°時的峰值為機翼前緣的電磁散射相干疊加而形成的，140°處的峰值為機翼端面造成的鏡面反射；RCS縮減部位主要集中在60°～90°范圍內，這是由于本文選擇的翼面控制點在翼面中部位置，對于前緣和端面形成的電磁散射無法做到有效縮減。

圖14 優化結果A與原始外形的RCS計算結果對比Fig.14 Comparison of RCS calculation results between Pareto A and original profile

3 結 論

本文使用計算流體力學和計算電磁學等數值模擬手段計算飛翼布局無人機的氣動和隱身性能，結合一種基于EHVI加點的高效多目標粒子群優化算法使用FFD法進行參數化表達，對一種低阻、低可探測飛翼布局無人機進行氣動隱身多目標優化設計。在200次調用目標函數的情況下，降低了無人機的阻力和雷達反射截面積，表明本文提出的高效優化設計方法在解決類似飛翼式布局無人機氣動隱身多目標優化設計等昂貴優化問題時具有較大的應用潛力。