改進Kriging模型及在寬馬赫數變幾何進氣道的應用*

徐珊珊，金玉華，張慶兵

(1.北京電子工程總體研究所，北京 100854；2.中國航天科工集團有限公司 第二研究院，北京 100854)

0 引言

寬馬赫數范圍的變幾何進氣道擁有不同的調節形式，如跳躍式唇口、轉動式唇口、楔板角度可調和楔板位置可調[1-3]。為選擇較為適合的調節方法，研究人員需要經歷仿真、建模和優化過程。整個研究周期需得到不同馬赫數下、不同攻角以及不同調節參數對不同進氣道性能參數的影響模型。盧燕[4]的建模方法采用馬赫數每增加0.1，仿真一次進氣道性能，最后使用最小二乘法擬合得到馬赫數和二級楔板角度變化的數學模型。而隨著飛行包線的進一步拓展以及組合調節方式的使用，一次建模需要大量實驗數據或仿真數據。這將消耗過高的經濟成本和計算成本。雖然文獻[5]提出了寬馬赫數變幾何進氣道的快速計算方法，但采用預測精度足夠高和計算速度足夠快的建模方法仍然十分重要。

最先被用于預測礦產儲量的Kriging代理模型方法不僅能提供最優無偏預測，還能對預測誤差進行定量估計。因此，該方法被大量運用到計算機輔助設計和工程建模領域[6-7]。相比于多項式響應面、BP (back propagation)神經網絡和徑向基神經網絡，Kriging方法對于非線性復雜問題的較好預測效果和較高精度，適宜于氣動特性預測[8-9]。通用的Kriging模型多采用模式搜索法進行參數優化，不再對相關參數求偏導，但由于較為依賴初始點位置，易收斂于局部極小值，從而嚴重影響預測精度[10]。陳志英等引入粒子群優化算法(particle swarm optimization，PSO)替換模式搜索方法，依靠PSO算法的群體搜索能力克服模式搜索法單點序列搜索方式產生的局限性和嚴重依賴初猜解的缺點，有效保證了在任意初始條件下都能獲取最佳預測精度的相關參數[11]。

雖然PSO具有程序簡單，多維并行優化，優化速度較快的優點，但其也具有一定局限，比如容易局部收斂、收斂速度較慢等，不適合運算量較大的建模過程。而SUN等提出的量子粒子群優化算法(quantum-behaved particle swarm optimization，QPSO)[12]，計算更為簡便，但隨機性較大，在多維優化問題中，更易收斂到局部最優解。這都會極大影響Kriging建模的精確度。

本文將帶全局判據的改進量子粒子群優化算法(improved quantum-particle swarm optimization with global criterion, GCIQPSO)與Kriging建模方法結合，提出了GCIQPSO-Kriging建模方法。GCIQPSO算法不僅降低了計算成本，而且可以保證每次優化結果均為全局最優解。為驗證GCIQPSO-Kriging建模方法的精確度，本文與標準Kriging方法和PSO-Kriging進行了比對。結果表明，本文方法計算速度快，精確度高，誤差較小，具有一定的工程應用價值。

針對進氣道調節問題，本文采用文獻[5]的快速計算方法和GCIQPSO-Kriging代理模型建模方法，對馬赫數范圍為2.5～4的二楔板反折式變幾何進氣道進行了仿真計算和近似模型建立。結果表明，本文方法對變幾何進氣道性能參數的非線性適應性較強，建模精度較高，可適用于變幾何進氣道調節方法的優化研究。

1 方法研究

1.1 Kriging方法[11]

Kriging方法假設所需要建立模型具有如下形式：

y(x)=fT(x)β+z(x)

,

(1)

式中:f(x)=(f1(x),f2(x),…,fp(x))T為回歸基函數；β=(β1,β2,…,βp)T為回歸系數；z(x)為高斯隨機過程，具有下列統計特性：

E[z(x)]=0，Var[z(x)]=σ2，
E[z(xi)z(xj)]=σ2R(xi,xj)

，

(2)

式中:xi與xj為2個任意輸入向量；R(*,*)為相關模型；σ2為過程方差。相關模型為各個維度上相關函數的乘積

(3)

設生成的樣本點為S={s1,s2,…,sm}，對應的仿真程序輸出為Y=(y(s1),y(s2),…,y(sm))T。對任意給定的θ=(θ1,θ2,…,θm)T，β與σ2的估值分別為

β=(FTR-1F)-1FTR-1Y

，

(4)

(5)

式中：基函數矩陣F∈Rm×p，第i行第j列元素為fj(si)；相關矩陣R∈Rm×m，第i行第j列的元素為R(si,sj)。

極大似然意義下的最優相關參數θ*應使似然函數ψ(θ)取最大值

(6)

完成對θ*的求解后即可得出建立在m次數值仿真結果基礎上的Kriging近似模型：

(7)

r(x)=(R(x,s1),R(x,s2),…,R(x,sm))T.

(8)

可以證明，在確保式(6)取極值的情況下，式(7)即為待測響應的最優無偏預測。此外，Kriging還能對預測結果的誤差進行定量估計

u(x)=FTR-1r(x)-f(x)

.

(9)

1.2 GCIQPSO方法

為保證每次優化結果均為全局最優解，并提高高維優化問題的建模精度，在QPSO算法和慣性權重自適應調整的量子粒子群算法(dynamically changing weight’s quantum-behaved particle swarm optimization，DCWQPSO)[13]基礎上，本文擴大粒子搜索范圍，并添加全局判據以判斷優化結果的全局性，建立了GCIQPSO優化算法。

QPSO算法引入δ勢阱模型，通過求解薛定諤方程得到相關波函數Ψ(x,t)，進而計算出粒子在設計空間內某一處的概率密度函數，最后由蒙特卡羅方法確定粒子的新位置。粒子的位置方程為

(10)

式中：xiop(t)=φixip(t)+(1-φi)xg(t) ；t為迭代次數；φi和ui為[0,1]上服從均勻分布的隨機數，粒子每一維取不同隨機數；xip(t)為第t次迭代時粒子的第i維局部最優位置，i=1，2，…，N；N為粒子群中粒子維數；xg(t)為第t次迭代時，粒子群的全局最優位置；L為勢阱的特征長度，代表了粒子更新的搜索范圍大小，Li(t+1)=2β|xpm(t)-x(t)|，其中xpm(t)表示為第t次迭代時種群粒子每一維度的局部最優位置平均值，表達式為

DCWQPSO算法把慣性權重β的表達式改進為與進化速度因子sd和粒子聚集度因子jd有關的自適應權重。表達式為

β=f(sd,jd)=β0-sdβ1+jdβ2

，

(11)

為增大粒子多樣性，避免粒子群陷入局部最優，GCIQPSO算法加入β值和粒子位置周期性變異。位置周期性變異每迭代10次變異一次：

xi(t)=xi(t)[1+(0.5-randi(0,1))δ]，

n=1,2,…，

randi(0,1)表示0和1之間的隨機數。

在50步以后，若β值在5步時均保持一個值，則對β值進行參數變異，具體方法如下：

βi(t)=βi(t)[1+(0.5-randi(0,1))δ]，

randi(0,1)表示0和1之間的隨機數。

為了保證優化結果均為全局最優解，GCIQPSO算法加入全局判據：①優化迭代中：粒子聚集度jd變化越劇烈，粒子的全局搜索能力越強，若其連續10步(粒子位置變異周期為10步)均保持不變，則跳出此次優化，重新初始化種群；②當優化參數小于等于5時，收斂速度較快，在迭代步數最大之前，已經收斂。因此迭代結束時的收斂判據主要考慮連續5次出現提前跳出，且數值近似，取最優值；在計算多維優化問題(優化參數大于5)時，迭代結束時：若優化過程沒有中斷，迭代步數達到最大迭代值，則判斷粒子聚集度jd均值和標準差是否滿足判斷條件。當迭代步驟為150～200步之間，則判斷條件為jd均值小于0.5，標準差小于0.3。若以上2個步驟均滿足，則本次優化結果收斂于全局最優解。

為驗證本文算法，本文選擇了非線性和局部極值個數依次遞增的4個不同函數：

Sphere函數：

Rosenbrock函數：

Rastrigin函數：

Griewank函數：

min(f(x*))=f(0,0,…,0)=0.

本文采用了改進粒子群優化算法[14](improvedparticleswarmoptimization，IPSO),QPSO算法和DCWQPSO算法與本文算法對比。由于QPSO和DCWQPSO算法主要針對優化參數不大于5的優化問題，為驗證本文算法精度，首先對4個函數進行了4維30個粒子組成的粒子群150步30次獨立優化。然后，針對多維優化問題，本文采用4種優化算法進行了30維30個粒子組成的粒子群迭代200步的30次獨立優化。所有優化結果的平均值列于表1。

從表1可以看出，當優化參數為4時，4種優化算法的Sphere,Rastrigin和Griewank函數優化結果均較好，但前3種算法隨著函數的非線性增加，優化精度均有一定下降。而QPSO算法的粒子量子性能使粒子位置更新隨機性更大，加快了算法收斂速度，降低了優化時間，但也易導致優化失敗。這使得Rosenbrock函數的優化結果極易不穩定，其30次獨立優化出現數值極大的失敗結果，致使平均優化結果非常大。而DCWQPSO算法有效地消除了一定隨機性，使得Rosenbrock函數優化結果有較大改進，但仍低于IPSO算法。而在4種算法中，本文算法的Rosenbrock函數優化結果最好，且速度也較快。其原因是全局判據可以有效消除粒子隨機性導致的優化失敗，而且雙重變異方法也可有效擴大全局搜索能力，增大收斂速度。由此可見，在低維優化問題中，雖然GCIQPSO算法基于QPSO算法，但有效避免了算法的不穩定性，并且在提高優化精度的同時，計算時間并未有明顯增加。

而當優化參數增加到30時，前3種優化方法的優化精度明顯下降，而GCIQPSO算法優化精度一直較高，雖然隨著函數非線性度的增加，優化時間有所提高，但仍短于IPSO算法。這是由于增大粒子全局搜索能力的雙重變異方法增大了高維粒子同時收斂到最優位置的概率，而全局判據有效地篩除了不良的優化結果。由此可見，在高維優化問題中，GCIQPSO算法的精度較高，適用性更強。現有優化算法均采用多次優化，人工篩選最優解的方式，但GCIQPSO算法可省略此步驟。

表1 4種函數的4種優化算法結果

在低維和高維優化問題中，GCIQPSO算法性能具有較短的優化時間和較高的優化精度，且擁有較好的函數非線性度適應性。

1.3 GCIQPSO-Kriging方法

本文參考文獻[11]的PSO-Kriging建模流程，建立了GCIQPSO-Kriging的建模方法，除2個優化算法的區別外，步驟均一致。具體步驟如下：

步驟1：根據待重構函數(即基于仿真程序的隱式函數)確定：輸入維數n；各維輸入的取值范圍；用于構造Kriging模型的樣本數m。對于實際問題，應該在主要關鍵點進行加密。

步驟2：生成樣本輸入矩陣S={s1,s2,…,sm}，將仿真程序重復運行m次，得出與S相對應的輸出向量Y=(y(s1),y(s2),…,y(sm))T。

步驟3：使用GCIQPSO算法搜索極大似然意義下的最優相關參數θ*：似然函數ψ(θ)為適應度函數；粒子維度為樣本點的維度n，粒子取值范圍采用對數分布，搜索范圍取為[-2,2]，運算形式為[lgθ1,lgθ2,…,lgθn]；種群規模按N≈5n。當建模參數較少時，優化收斂速度較快，jd保持不變的步驟提前，易跳出優化程序，因此，下一次優化初始粒子群加入本次Gbest粒子，若5次優化均跳出程序，采用5次優化的最優解。

步驟4：根據得到的最優相關參數θ*確定出回歸系數β及相關矩陣R，任意向量x的相應依據公式(7)預測得到。

步驟5：若精度不足，采用式(9)確定出的誤差極值區域補充樣本并返回步驟3，重新優化。

2 方法驗證

為驗證不同Kriging代理模型的精確度，本文采用以下兩種測試函數進行了重構：

f1(x)函數非線性程度較高，用以測試代理模型對非線性問題的預測精確度。式f2(x)用以驗證不同參數個數對建模精確度的影響。模型精確度采用3個參數衡量代理模型的預測精度，包括：代理模型在測試樣本點的百分比最大相對誤差(maximum relative error，MAE)、百分比平均相對誤差(average relative error，ARE)和百分比均方根差(root mean square error，RMSE)，其表達式為

為消除因額外補充的樣本點增加的模型精度，本節函數重構不添加步驟4。f1(x)函數的初始樣本點個數選擇自變量的200倍，f2(x)函數的初始樣本點個數選擇自變量的20倍，隨機選擇441個預測點建立模型精度比對標準。為消除建模計算的偶然性，本文采用10次建模的平均數據具體數據，如表2所示。其中，式f2(x)采用不同參數個數n=2,3,5,10，對應于a,b,c和d，IPSO-Kriging方法的優化算法采用XIA[14]的改進粒子群算法(IPSO)，使得全局搜索能力更強，收斂更快。

由表2可見，在消除了對優化初始值的依賴后，后2種模型的平均精度較Kriging方法均有所提高。IPSO優化算法收斂速度較慢，在迭代步數下降的情況下，更易陷入局部最優點，這大大影響了IPSO-Kriging模型的精度。而由于樣本點選擇個數較少，Kriging和IPSO-Kriging模型精度較文獻[11]有所降低。但GCIQPSO-Kriging代理模型方法的3個參考值明顯最優，且對非線性更高的函數適應性較好。而隨著模型參數個數的增加，GCIQPSO-Kriging優勢也有所增加，這與GCIQPSO優化方法在多參數優化的更好性能有關。若添加步驟5，進一步補充樣本點個數，GCIQPSO-Kriging模型精度可進一步提高。

表2 3種建模方法精確度

3 變幾何進氣道建模

寬馬赫數變幾何進氣道建模的參數較多、非線性度較高，因此，導致計算量大且精度不足。本文采用以上3種方法進行建模，以檢驗本文算法的精度。

對于混壓進氣道而言，變幾何需要在保持較大流量系數的同時，盡量提高臨界總壓恢復系數。而臨界總壓恢復系數關鍵是喉道高度改變[15]，在正激波產生之前，盡量提高斜激波壓縮效率，可有效提高臨界總壓恢復。因此，本文選用第2級楔板角度和位置均可調的變幾何方式提高臨界總壓恢復系數，跳躍唇口提高流量系數(如圖1所示)。以2.0為接力馬赫數，2.8為設計馬赫數Mad(即激波封口馬赫數)。根據空氣動力學理論，接力狀態下(Ma=2.0)氣流的最大折角為26°，考慮5°的攻角裕度，確定進氣道的外壓楔板總折角δ最大為10°，一級楔板折角δ1設計為5°。在Ma=2.0時，二級楔板折角δ2最大為5°，其與喉道高度可隨飛行條件變化。在一級楔板和二級楔板聯接處，裝有一滑動鉸，使二級楔板在轉動的同時，沿一級楔板滑動，亞聲速擴壓段楔板可伸縮，以補償喉道平直段牽連運動所造成的水平位移。本文運用文獻[5]的二級楔板反壓式進氣道快速計算方法得到進氣道性能參數，采用GCIQPSO-Kriging代理模型方法進行建模，并在非線性較高的區域，預先進行樣本點加密。

圖1 變幾何進氣道示意圖Fig.1 Schematic of geometry-variable inlet

模型的5個自變量為：馬赫數Ma(2～4)、攻角α(-5°～5°)、二級楔板角度δ2(下界為5°，上界為來流Ma下、5°攻角裕度時最大角度)、一級楔板長度OT1(上界以不造成喉道堵塞為準，下界為接力狀態起始位置)和唇口位置(上下界分別以接力馬赫數和最大馬赫數封口位置為準)。所求的2個預測值為進氣道臨界總壓恢復系數σ和流量系數φ。由文獻[5]可知，唇口位置的變化對σ沒有影響，故在σ的建模過程中，去除唇口位置這一變量。本文仍然采用400個預測比對點，以驗證模型的精確度，3種方法的精確度數據對比列于表3。

表3 3種建模方法精確度

從文獻[5]可以看出，作為變幾何進氣道性能參數的流量系數φ與臨界總壓恢復系數σ較第2節的2個測試函數均較弱，因此，表3顯示的3種代理模型方法精度均較高。而臨界總壓恢復系數σ非線性程度稍高，與前2個代理模型方法相比，本文方法誤差更小。

可見，對于寬馬赫數變幾何進氣道的建模問題，本文方法精度較高，非線性適應性更佳，若用于該問題的優化，將大大降低計算時間，并可有效提高優化結果可信度。

4 結束語

針對寬馬赫數變幾何進氣道優化過程中，建模參數多、精度不夠、計算時間長和非線性度高的問題，本文改進了Kriging代理模型方法，采用帶全局判據的改進量子粒子群算法替代原有的優化部分，以提高優化搜索時間，擴大搜索范圍，進一步降低對初猜解的依賴，得到最大似然函數的全局最優解。本文采用驗證函數和變幾何進氣道算例，對該方法進行了驗證。結果表明，本文方法有效提高了代理模型對多參數建模和非線性問題的適應力，對于寬馬赫數變幾何進氣道建模問題精度較高，可適用于該問題的優化計算。

[1] BOUCHEZ M,FALEMPIN F.Air breathing Space Launcher Interest of a Fully Variable Geometry Propulsion System Status 1999[R].AIAA-99-2376,1999.

[2] HUEBNER L D,ROCK K E,RUF E G,et al.Hyper-X Flight Engine Ground Testing for X-43 Flight Risk Reduction[R].AIAA-2001-1809,2001.

[3] FALEMPIN F,WENDLING E,GOLDFEILD M,et al.Experimental Investigation of Starting Process for a Variable Geometry Air Inlet Operating from Mach 2 to Mach 8[R].AIAA-2006-4513,2006.

[4] 盧燕,樊思齊,馬會民.超聲速進氣道數學模型研究[J].推進技術,2002,23(6):468-471. LU Yan,FAN Si-qi,MA Hui-min.Study for the Mathematical Model of Supersonic Inlet[J].Journal of Propulsion Technology,2002,23(6):468-471.

[5] 徐珊珊,金玉華,張慶兵,等.寬馬赫數變幾何進氣道性能快速計算方[J].現代防御技術,2017,45(2):74-79+92. XU Shan-shan,JIN Yu-hua,ZHANG Qing-bing,et al.A Mathematical Model for Fast Design of the Variable Geometry Supersonic Inlets with Large Mach Number Range[J].Modern Defence Technology,2017,45(2):74-79+92.

[6] KAYMAZ J.Application of Kriging Method to Structural Reliability Problems[J].Structural Safety,2005,27(2):133-151.

[7] DELLINO G,LINO P,MELONI C,et al.Kriging Metamodel Management in the Design Optimization of a CNC Injection System[J].Mathematics and Computers in Simulation,2009,29(10):2345-2360.

[8] 趙軻.基于CFD的復雜氣動優化與穩健設計方法研究[D].西安：西北工業大學,2014. ZHAO Ke.Complex Aerodynamic Optimization and Robust Design Method Based on Computational Fluid Dynamics[D].Xi’an:Northwestern Polytechnical University,2014.

[9] 吳先宇,羅世彬,陳小前.基于替代模型的高超聲速進氣道優化[J].彈箭與制導學報,2008,28(1):148-152. WU Xiang-yu,LUO Shi-bin,CHEN Xiao-qian.Optimization of a Hypersonic Inlet Based on Surrogate Models.Journal of Projectiles,Rockets,Missiles and Guidance,2008,28(1):148-152.

[10] MARTIN J D,SIMPSON T W.Use of Kriging Models to Approximate Deterministic Computer Models[J].AIAA Journal,2005,43(4):854-863.

[11] 陳志英,任遠,白廣忱,等.粒子群優化的Kriging近似模型及其在可靠性分析中的應用[J].航空動力學報,2011,26(7):1522-1530. CHEN Zhi-ying,REN Yuan,BAI Guang-chen,et al.Particle Swarm Optimized Kriging Approximate Model and its Application to Reliability Analysis[J].Journal of Aerospace Power，2011,26(7):1522-1530.

[12] SUN Jun,XU Wen-bo,FENG Bin.A Global Search Strategy of Quantum-Behaved Particle Swarm Optimization[C]∥Proceedings of the IEEE Congress on Cybernetics and Intelligent System，Singapore:IEEE Press,2004:111-116.

[13] 黃澤霞,俞攸紅,黃德才,等.慣性權自適應調整的量子粒子群優化算法[J].上海交通大學學報，2012，46(2):228-232. HUANG Ze-xia,YU You-hong,HUANG De-cai,Quantum-behaved Particle Swarm Alogorithm with Self-Adapting Adjustment of Inertia Weight[J].Journal of Shanghai Jiaotong University,2012，46(2):228-232.

[14] XIA Chen-chao,JIANG Ting-ting,CHEN Wei-fang.Particle Swarm Optimization of Aerodynamic Shapes with Nonuniform Shape Parameter-Based Radial Basis Function[J].Journal of Aerospace Engineering,2017,30(3):04016089 1-12.

[15] 蔡飛超,陳鳳明,徐東來,等.寬Ma數范圍固定幾何進氣道設計問題研究[J].固體火箭技術,2010,33(2):163-166. CAI Fei-chao,CHEN Feng-ming,XU Dong-lai,et al.Study on Fixed-Geometry Supersonic Inlet Design for Wide Mach Number Range Application[J].Journal of Solid Rocket Technology,2010,33(2):163-166.