基于代理模型全局優化的自適應參數化方法

張偉，高正紅，周琳，夏露

西北工業大學 航空學院，西安 710072

飛翼布局[1-4]因其具有良好的氣動效率和隱身特性，被新一代軍用飛機采用。相較于常規布局，由于尾翼等部件的缺失以及翼身融合的特征，翼型的作用對于飛翼布局更加重要。在氣動特性方面，相較于傳統翼型設計主要考慮升阻特性要求，飛翼氣動布局則要求在翼型設計中增加量化的俯仰力矩以滿足飛機縱向力矩配平的需要。在隱身特性方面，翼型的隱身設計是飛翼布局雷達散射截面(Radar Cross Section,RCS)減縮的重要手段之一[5]。因此，飛翼布局翼型設計需要在滿足氣動、隱身及控制等多學科性能要求下進行。

對于翼型設計而言，設計變量(Design Variable，DV)的選擇直接影響設計結果[6]。在開展滿足氣動、隱身與控制等多學科影響的翼型設計時，需要在更寬范圍內對設計變量的影響進行研究。本文以基于代理模型的全局優化設計(Surrogate Base Optimization，SBO)[5-12]方法為基礎，開展了設計空間對翼型氣動與隱身設計結果的影響分析，并提出了基于自適應參數化翼型氣動隱身一體化優化設計方法。

影響氣動外形優化結果的設計空間包括尺度和維度兩個方面。為了解決空間尺度不足導致無法找到理想設計結果的問題，傳統大多采用多輪優化方法[13]，即以上一輪優化的結果作為新的初始外形，擴大設計空間尺度再次進行優化，但這種“人在回路”的方法無疑會造成大量的計算資源消耗以及設計者工作量的增加。為此，王超等[6]在研究了設計變量每一維的尺度對設計空間及代理模型影響的基礎上，提出了自適應設計空間尺度擴展方法，即在設計過程中選擇性擴張設計空間尺度，有效提升了優化設計質量和效率。

然而在氣動設計中，設計變量及其數量的選擇對設計結果影響很大[14-15]。以往設計者多是根據經驗或直接布置大量設計變量，但對于新型翼型而言，設計者一般缺乏必要的經驗，而大幅增加設計變量，將對基于代理模型的全局優化造成“維度災難”現象。

為了緩解設計空間維度擴張對優化算法的壓力，設計者們開展了自適應參數化方法研究。Desideri[16]、Masters[17]等采用逐步加密設計變量的方法增加設計空間維度，但這樣增加設計變量缺乏針對性，造成計算資源的浪費；Anderson[18]、Zingg等[19]以設計變量梯度信息作為敏感性信息，選擇增加梯度較大的設計變量，但對于進化類算法而言，設計變量梯度難以直接獲取，并且梯度信息只能反映設計變量在局部設計空間的敏感性。此外，設計空間維度擴張后，已有的低維空間產生的樣本難以在高維空間內直接利用，重新在高維空間取樣則會帶來額外的計算負擔，這更加限制了自適應參數化方法在全局優化算法中的應用。

針對現有的自適應參數化方法難以在代理模型全局優化中應用的問題，本文首先研究不同維度設計空間對翼型設計結果影響，對比分析了不同設計變量對翼型氣動與隱身特性影響的敏感性，在此基礎上提出了適用于代理模型全局優化方法的自適應參數加點方法：為了避免設計變量分布過度集中，引入密度閾值控制設計變量分布；針對維度擴展后低維樣本無法利用的問題，應用節點插入算法將已有的低維樣本轉化為高維樣本，避免了重新取樣的工作量，實現了樣本的高效配置。最后通過飛翼布局翼型氣動隱身設計，并與不同設計空間維度設計結果進行對比，驗證了該方法的可靠性。

1 基本方法

1.1 基本效應法

基本效應法能夠憑借較少的計算量識別計算模型Y=f(X)中對輸出量Y影響較大的輸入變量X，該方法尤其適用于計算量較大的物理模型，是應用最為廣泛的全局敏感性分析方法之一[20]。基本效應法通過OAT(One at A Time)[21]的方式計算任一維輸入量Xi對輸出量Y的影響，即任一維輸入量在取值范圍內每次變化Δi，其對輸出量的敏感性信息計算方式為

EEi=(f(X1,X2,…，Xi+Δi，…，Xn)-

f(X1,X2,…，Xi，…，Xn))/Δi

(1)

文獻[20]提出徑向取樣(Radial Sample)法,通過在整體設計空間內進行r次隨機取樣，可以得到每一維設計變量的全局敏感度信息SAi：

(2)

由式(2)可知，SAi的計算結果存在正負問題，可能導致計算結果不能充分反映設計變量的敏感性信息。因此，Campolongo等將式(2)轉化為式(3),以消除正負號對結果的影響。

(3)

對于氣動優化問題而言，通過優化過程獲得樣本中的設計變量對目標函數的敏感度分析，可以獲得不同設計變量對目標特性影響的敏感度，并根據敏感度選擇擴展設計空間維度。

1.2 Bspline 曲線參數化方法

Bspline[22]方法憑借靈敏的局部擾動特性，廣泛應用于工業設計中。通過一組基函數及其系數實現對外形的擬合和擾動，基函數的系數稱為控制點，其位置由橫、縱坐標共同決定。在翼型參數化中通常采用三次樣條曲線。

Bspline通常由一組u∈[0,1]的標量來組成。

(4)

(5)

(6)

式中：遞增節點序列r={r0,r2，…,rl},ri≤ri+1，u={r0,…,r0,rk+1,…,rn-1,rl,…rl}，重節點r0、rl的數量均為k+1，其余節點一般均勻分布。

相較于其他經典參數化方法，如CST[23]、Hicks-Henne[24]等,Bspline曲線具備良好的局部擾動特性[25]。

由式(6)可以看出，Bspline曲線每段曲線的特性由相應節點矢量決定，因此Bspline曲線控制點的增加是從節點開始的。在節點區間[ri,ri+1]插入新的節點r，從而得到新的節點標量u1={r0,…,ri,r,ri+1,…，rl}，并由新的節點矢量確定新的基函數和控制點矢量[26]：

Pnew,j=

(7)

式中：αj=(t-tj)/(tj+k-1-tj)，這樣一來，原來的樣條曲線就可以用新的基函數與控制點矢量表示，并與初始曲線完全一致。控制點的增加僅影響其左右相鄰控制點的位置。因此，Bspline曲線節點插入算法的特性能夠實現低維樣本在高維設計空間的重構，避免了重新取樣，提高了樣本利用率[26]。

1.3 二維矩量法

本文采用基于電場積分方程的二維矩量法[27]進行翼型RCS特性計算。

(8)

矩量法計入了各部分感應電流相互之間的影響，是麥克斯韋方程的精確求解方法。由于矩量法離散的面元數與照射波長成正比，對于飛行器這種大尺寸物體的高頻計算，三維矩量法需要的計算量和存儲量往往難以承受。而二維情況下計算量較小，雖然電磁散射的二維情況在實際中很少遇到，但是其對分析相應的三維問題具有很好的參考意義。

二維情況下，散射體為無限長柱體，雷達散射截面σ變為雷達散射寬度σ′：

(9)

式中：r為柱體到散射觀察點的距離；ρ為表面電荷密度；Es和Hs均為目標散射電磁波在雷達處的電磁場強度。

三維情況下長度為有限值l的柱形結構的雷達散射截面σ與無限長相同結構的二維雷達散射寬度σ′之間有如下簡單關系[27]：

(10)

式中：λ是照射波長。

在隱身數值計算中，EMSS公司發布的三維全波電磁仿真商用軟件FEKO得到了較為廣泛的認可。為了驗證本文應用的二維矩量法的可靠性，下面以NACA65，3-018翼型為例，分別采用商用計算軟件FEKO(三維矩量法)和二維矩量法程序，在入射頻率分別為3 GHz、5 GHz、7 GHz，入射角度為±90°的條件下進行計算，其中，二維矩量法的計算結果經式(10)轉化為三維結果，并與FEKO軟件計算結果進行對比。FEKO計算模型為三維模型，如圖1所示，模型展長為1 m。計算環境為Inter(R)Core(TM) i7-8700 CPU @3.2 GHz RAM 32 G。

圖1 FEKO三維計算模型Fig.1 3D computation model of FEKO

圖2為計算結果對比，其中橫坐標為入射方位角θ。可以看出，二維矩量法程序與FEKO軟件的計算結果吻合良好，計算精度滿足要求。因此，在針對翼型的氣動隱身協同設計中，本文應用二維矩量法程序進行翼型RCS特性的數值模擬。

圖2 FEKO與二維矩量法計算結果對比Fig.2 Results comparison of FEKO and 2D moment method

2 設計空間維度對翼型氣動隱身特性影響

2.1 不同設計空間維度設計結果對比

本文以相對厚度為18%的對稱翼型NACA65，3-018為例[28]，選取典型跨聲速設計狀態馬赫數Ma=0.70,升力系數CL=0.25。翼型氣動特性要求在保持量化的抬頭力矩時提高升阻特性；隱身特性要求在入射頻率9 GHz下，前向±30°均值降低，入射角θ如圖3所示。

圖3 入射角示意圖Fig.3 Sketch map of incident angle

設計目標為減小阻力系數CD和提高前向RCS特性，氣動約束為力矩系數Cm不小于0.03，幾何約束為翼型最大厚度thickmax不減小，設計模型為

objCD,RCS

(11)

在取樣過程中，氣動特性計算采用基于RANS方程的Cfl3d程序，隱身特性計算采用二維矩量法[27]。針對減阻和降低前向RCS設計目標，本文將兩者進行正則化，然后采用線性加權的方式轉換為一個目標進行優化搜索，權重比例為1∶1。

為了對比分析設計變量對目標特性的影響，本文首先采用設計變量均勻分布的方式，設計變量數(DV)即空間維度分別為12、18、24、30。利用拉丁超立方(LHS)方法在各設計空間內分別取樣100個，加點500次。表1為不同設計空間維度下的設計結果，可以看出在滿足力矩約束下，各組設計外形阻力特性和隱身特性都有很大的提升，圖4為前向RCS對比。從圖5對比中可以看出，低維設計空間(DV_12)收斂效率高，但低維空間內難以搜索到理想解，因此目標值很快陷入停滯；隨著設計空間維度增加(DV_12→DV_24)，較之低維空間此時得到了更好的設計結果；當設計空間維度均勻增加30維(DV_30)時，由于設計變量分布缺乏針對性，且空間維度擴張導致代理模型精度降低，因此在相同的計算量下，得到的結果反而不如低維設計空間的結果。

圖4 設計翼型前向RCS特性對比Fig.4 Forward RCS comparison of optimized airfoils

圖5 目標值收斂歷程Fig.5 Convergence history of fitness

表1 翼型計算結果Table 1 Computational results of airfoils

設計翼型外形和壓力分布對比分別如圖6、圖7所示。各組設計結果外形接近：前緣半徑都有所變小，具備明顯的前緣正加載特征和中后部反加載，同時最大厚度位置后移。從結果對比中可知各組設計外形基本相似，說明優化搜索的方向是正確的，但設計變量的選擇直接影響設計結果，簡單地增加設計變量并不能確保改善設計結果。

圖6 設計翼型外形對比Fig.6 Shape comparison of optimized airfoils

圖7 設計翼型壓力分布對比Fig.7 Pressure distribution comparison of optimized airfoils

2.2 設計變量關于氣動隱身特性敏感性分析

設計空間由設計變量及其數量決定，為了分析設計空間內關于氣動、隱身特性敏感的區域，本節以DV_12算例為例，對設計變量進行敏感性分析。

圖8、圖9分別為設計變量關于翼型阻力系數和前向RCS的敏感性評估結果，圖中前6個為上表面設計變量，后6個為下表面變量。對于阻力特性而言，翼型上表面外形主要影響著激波的強弱和發展，因此這部分設計變量敏感性更強；對于隱身特性而言，上、下表面第1個設計變量的影響遠大于其他變量，而這兩個設計變量主要影響翼型前緣20%左右的區域，因此這部分區域對前向RCS特性影響更大。

圖8 設計變量關于阻力系數敏感性Fig.8 Sensitivity to CD of design variables

圖9 設計變量關于前向RCS敏感性Fig.9 Sensitivity to forward RCS of design variables

結合減阻和降低前向RCS的設計目標，為了得到更好的設計結果，需要在設計空間敏感區域擴展維度。而均勻分布的設計空間內，設計變量的分布缺乏針對性，為了得到理想的結果往往需要更多的設計變量，即更高的維度，而維度的擴展意味著設計空間呈指數增加，導致代理模型效率急劇降低，需要更多的樣本才可能得到理想的設計結果。因此，在設計空間敏感的區域進行維度擴展，使得每個設計變量充分發揮作用，避免作用不大的設計變量的增加導致空間的過度增長，提高代理模型的效率的同時，得到更加理想的設計結果。

3 自適應參數化方法

對基于代理模型的全局氣動優化而言，設計空間的構造對優化設計結果和效率有很大的影響，而設計者一般難以直接給出合適的設計空間。本文提出的自適應參數化方法，設計初始階段只需憑借較少的設計變量定義低維設計空間，在優化過程中逐步在設計空間敏感區域增加維度，自適應構造設計空間，提升優化設計質量。

本文構建的自適應參數化方法基本框架如圖10所示，具體步驟為

圖10 自適應參數化方法流程Fig.10 Flowchart of adaptive parameterization method

步驟1初始設計空間確定，利用較少的設計變量確定低維設計空間。

步驟2初始樣本庫構建，利用LHS方法在初始設計空間內生成樣本，并加入樣本庫中。

步驟3訓練代理模型并優化，利用樣本庫內的樣本訓練代理模型，在當前設計空間內進行優化搜索，將搜索到的最優值進行評估，更新樣本庫和代理模型。

步驟4設計空間維度擴展判斷，維度擴展需滿足如下條件：① 目標值長期陷入停滯；② 代理模型精度滿足要求。

本文定義目標函數值的更新效率Δt和代理模型精度f*，作為維度擴展的判斷：

Δt=fgi-t+1-fgi

(12)

(13)

Tr=f*/t

(14)

本文設置目標函數值更新的閾值objtr，代理模型精度閾值c。

當Δt≤objtr，即目標函數優化效率不滿足要求，以及Tr≤c，即代理模型精度滿足要求的情況下，進行設計空間維度的擴展。

步驟5生成候補設計變量，進行維度擴展時，首先根據已有節點[ri,ri+1]生成新的節點r：

(15)

然后利用節點插入算法，在設計變量[Pi,Pi+1]間生成候補設計變量Pnew，在此過程中根據設計變量密度分布閾值dρ對候補設計變量進行第一次遴選：

Pnew=

(16)

式中：|Pnew-Pi|為設計變量橫坐標距離。

將不滿足密度分布閾值要求的候補變量去掉(即為空集[])，留下滿足要求的變量進行下一步。

步驟6設計空間維度擴展，通過對候補設計變量敏感性分析進行第2次遴選，敏感性分析過程如第2節所示，將滿足敏感度閾值的候補變量加入已有設計變量中，實現設計空間維度擴展。

步驟7確定新增設計變量取值范圍，新增候補變量的取值范圍VRnew由其相鄰變量的取值范圍VRi，以及設計變量橫坐標距離|Pnew-Pi|共同決定，計算公式為

(17)

步驟8樣本維度變換，結合新增候補設計變量的位置，同樣利用節點插入算法將樣本庫內的低維樣本轉換為相應的高維樣本，避免了在新的高維設計空間內重新取樣的工作量。

步驟9繼續進行優化，利用變換后的高維樣本重新訓練代理模型，在高維設計空間繼續進行優化。

步驟10直至收斂，程序收斂條件為預設迭代數，滿足上述條件即終止程序，否則轉到步驟3直至滿足停止條件。

4 NACA 65,3-018翼型氣動隱身設計

采用本文的自適應參數化方法對NACA 65,3-018翼型重新進行優化設計，算例初始設置與DV_12保持一致，設計目標依舊為減阻和提高前向RCS特性，氣動約束為力矩系數不小于0.03，幾何約束為翼型最大厚度不減小，設計模型為

objCD，RCS

(18)

將設計結果與固定設計空間維度方法進行對比。

4.1 自適應參數化方法設計結果

表2為設計翼型與初始翼型的特性評估對比，相較于初始翼型，設計翼型阻力系數降低26 counts，力矩特性滿足約束，隱身特性提升92%，氣動與隱身特性均有明顯提升。應用自適應參數化方法的優化收斂歷程如圖11所示，設計空間維度經歷12-18-24-30的變化過程，在此過程中，目標函數最優值保持了較高的收斂效率，在120代左右達到收斂。

圖11 自適應參數化方法目標值收斂歷程Fig.11 Fitness convergence history of adaptive parameterization method

表2 自適應參數化方法設計翼型計算結果Table 2 Computational results of optimized airfoil of adaptive parameterization method

設計翼型與初始翼型的外形和壓力分布對比分別如圖12、圖13所示。應用自適應參數化方法設計的翼型外形與第2節設計結果相似，前緣半徑呈現了較為明顯的“鷹嘴”特征，上表面趨于平坦，最大厚度位置后移，具備明顯的前加載和中后部反加載；從壓力分布對比可知，由于設計翼型的前緣外形呈現從鈍頭體到“鷹嘴”的變化，其頭部尤其是最前緣的外形曲率變化較為劇烈，使設計翼型壓力分布出現了較小的吸力峰，為了滿足抬頭力矩的約束，設計翼型壓力分布前緣正升力區域擴大，中后部出現負升力區域，初始翼型的強激波削弱為兩道弱激波。設計翼型的RCS特性如圖14所示。

圖12 設計翼型外形Fig.12 Shape of optimized airfoil

圖13 設計翼型壓力分布Fig.13 Pressure distribution of optimized airfoil

圖14 設計翼型前向RCS特性Fig.14 Forward RCS of optimized airfoil

圖15和圖16展示了優化過程中節點和設計變量位置的變化。相較于初始，設計變量數量明顯增加，并且設計變量的分布更加具有針對性。翼型前緣設計變量明顯增加，有利于提升設計翼型的前向RCS特性；為了減弱激波，翼型上表面中部區域的設計變量分布更加密集；此外，翼型后緣下表面的變量也有增加，從最終外形可知，翼型后緣的設計變量使得其具有明顯的后緣反加載特征，配合前緣的正加載，使設計翼型在滿足力矩特性的約束下，獲得更好的阻力和隱身特性。

圖15 自適應參數化過程增加的節點Fig.15 Added knots in adaptive parameterization method

圖16 初始與最終階段設計變量橫坐標位置對比Fig.16 Design variable position comparison between initial and final stages

4.2 設計結果對比

為了驗證自適應參數化方法的可靠性，本文將其與固定設計空間維度的設計結果進行比較。圖17和圖18分別為各組設計翼型的外形和收斂歷程對比，可以看出，各組設計的翼型外形特征相似。各組設計翼型的特性評估對比如表3所示，在滿足力矩系數的約束下，應用自適應參數化方法調用的CFD計算次數更少，設計的翼型其阻力特性和隱身特性更優于固定設計空間維度的設計結果。

圖17 不同方法設計翼型外形對比Fig.17 Shape comparison of optimized airfoils of different methods

圖18 不同方法優化收斂歷程對比Fig.18 Comparison of fitness convergence histories in different methods

表3 優化設計翼型計算結果Table 3 Computational results of optimized airfoils

通過以上對比結果可以發現，自適應參數化方法能夠在設計空間的敏感區域逐漸擴展設計空間維度，相較于均勻分布的設計空間，自適應構造的設計空間能夠更加精準地描述目標外形，反映目標函數的變化趨勢，充分發揮設計變量的作用，提高設計質量。

5 結 論

本文針對飛翼布局翼型的氣動隱身設計問題，探討了設計空間維度對設計結果的影響。針對設計空間維度構造問題，提出一種適用于全局優化的自適應參數化方法，并進行了算例驗證，得到結論如下：

1) 基于設計空間維度對優化結果和效率的影響分析，提出一種適用于代理模型全局優化的自適應參數化方法。該方法在設計過程中從低維設計空間出發，結合設計變量的全局敏感性信息，在設計空間敏感性強的區域擴展維度。擴展后的高維設計空間能夠更加精準地描述目標外形，反映目標變化的趨勢，滿足高效精細化設計需求。

2) 在設計空間維度擴張后，利用節點插入算法，將樣本庫內已有的低維樣本在高維空間內重構，以重新訓練代理模型，避免了在高維空間重新取樣，實現了樣本的高效配置。

3) 通過NACA 65,3-018翼型氣動隱身一體化設計算例，驗證了本文提出的自適應參數化方法能夠在優化過程中自適應配置設計變量，構造合理的設計空間，得到更加理想的設計結果。并與固定設計空間維度方法進行了對比，從設計質量和設計效率的角度分別驗證了該方法的優越性。