， ， ，志勤，

(1.西南科技大學 計算機科學與技術(shù)學院,四川 綿陽 621010；2.中國空氣動力研究與發(fā)展中心 吸氣式高超聲速技術(shù)研究中心,四川 綿陽 621010)

0 引言

高超聲速進氣道作為飛行器發(fā)動機關(guān)鍵氣動部件[1]，若發(fā)生進氣道不起動現(xiàn)象直接影響著推進系統(tǒng)整體效能的發(fā)揮，發(fā)動機不能產(chǎn)生推力。若能對進氣道工作狀態(tài)進行及時的判斷與預(yù)測，對進氣道的再起動控制和發(fā)動機安全運行起關(guān)鍵作用，因此在進氣道的研究控制中，進氣道起動/不起動判斷與預(yù)警已經(jīng)成為一個重要的研究工作[2-3]。通常使用模式分類來解決進氣道起動/不起動判斷與預(yù)警[4]，數(shù)值計算產(chǎn)生進氣道起動/不起動樣本數(shù)據(jù)作為訓練集，采用計算結(jié)果中多個壓力點和壓力值作為特征屬性和特征屬性值。但在實際應(yīng)用中對壓力點的要求是最少最靈敏的壓力點，因此在對進氣道起動/不起動判別預(yù)警中，特征選擇算法尤為重要，一般要求選擇出2個最優(yōu)壓力點，再在二維空間中訓練分類準則，對進氣道起動/不起動現(xiàn)象進行準確判別。目前可用的特征選擇算法僅有SVMRFE，分類算法有FLD和SVM[4-5]，而SVMRFE算法只能每次刪除一個特征，這對于維度較高的數(shù)據(jù)集，大大增加了計算負擔。

特征選擇算法在機器學習中常用到的一類算法，可以分為wrapper和filter兩大類[6-7]，wrapper類算法是利用后續(xù)學習算法的訓練準確率來評估特征子集，因此評估結(jié)果和后續(xù)算法性能偏差較小，但由于評估每個特征子集時都要運行學習算法，計算量很大，大的數(shù)據(jù)集耗時會很長；filter類算法和后續(xù)學習算法無關(guān)，一般直接利用所有訓練數(shù)據(jù)的統(tǒng)計性能評估特征，如相關(guān)性、一致性、信息增益等，因此filter類算法具有速度快的優(yōu)點，但評估結(jié)果和后續(xù)算法的性能有較大的偏差。而這兩類算法中最具代表性的算法是SVMRFE和Relief，SVMRFE算法準確性能有一定的保證，但是運行時間較長，而Relief耗時少，但是Relief的缺點是不能去掉冗余特征。

目前基于Relief的改進算法也較為常見，使用Relief算法快速去掉一大部分無關(guān)特征，再進一步有效地去掉冗余特征[8-10]。本文在數(shù)值計算結(jié)果中，內(nèi)流道上表面會多達幾百個壓力點，數(shù)據(jù)維度較高，所以SVMRFE算法選擇特征耗時較長，本文結(jié)合Relief算法和SVMRFE算法進行特征選擇，尋找新的處理進氣道壓力數(shù)據(jù)的方法，提高運行效率，并對結(jié)果準確率進行分析。

1 數(shù)據(jù)準備

本文的數(shù)據(jù)集是通過數(shù)值計算而獲得的，包括進氣道起動/不起動兩類樣本，壓力點為屬性，對應(yīng)壓力值為屬性值。

1.1 數(shù)值計算

1.1.1 物理模型和計算網(wǎng)格

給出高超聲速進氣道/隔離段模型見圖1，是參考文獻[11]進行建模。圖中點A為鼻錐前緣，點B為二級楔面起點，點C為進氣道外罩唇口，點D為進氣道喉道部位，點E為隔離段出口位置。L1,L2,L3為點A到點B,C,D的位置，L4為隔離段長度，H1,H2分別為喉道和出口的高度。圖中具體數(shù)據(jù)如表1所示。

圖1 進氣道/隔離段幾何模型

L1/mL2/mL3/mL4/m0.63001.30151.84750.3620H1/mH2/mσ1/(°)σ2/(°)0.0440.06002.76.2

計算網(wǎng)格如圖2，網(wǎng)格劃分軟件為Gridgen，網(wǎng)格類型為結(jié)構(gòu)網(wǎng)格，其中壁面法向第一層網(wǎng)格間距取值為0.0005 mm，y+<1。二維網(wǎng)格總數(shù)約為12.7萬。

圖2 計算網(wǎng)格

1.1.2 數(shù)值計算方法

采用fluent14.5軟件[12]進行數(shù)值模擬，其中使用的控制方程為RANS方程，本文采用的湍流模型為SST k-ω，無粘通量離散格式為AUSM，熱完全氣體模型，采用定常模擬方法，數(shù)值計算的收斂準則為連續(xù)方程的標準化殘差小于10-3。

參考飛行試驗中X-51A飛行器飛行條件，設(shè)置邊界條件如表2，設(shè)置了3組不同的動壓條件，每組動壓條件下分別對應(yīng)的馬赫數(shù)分別為4.5,5.0,5.5，攻角大小分別為[0,2°,4°,6°,8°]。計算自由流壓力條件下的流場，此時喉道質(zhì)量加權(quán)平均壓力為P，背壓比定義為隔離段出口條件中的壓力Pb與壓力P之比，即背壓比r=Pb/P。按[2,3,4,5…]逐步增大背壓比，計算流場直到不起動，隨后在相鄰的起動/不起動背壓比之間，采用二分取值法細分背壓比來計算流場。

表2 自由來流條件

1.1.3 數(shù)值計算結(jié)果數(shù)據(jù)分析

變化的背壓引起進氣道是否起動主要體現(xiàn)在隔離段激波串前沿位置的變化，變化的背壓，使得激波串幾乎完整地進行簡單滑動和平移。如果背壓過大，會出現(xiàn)激波串前沿越過喉道位置處的情況，本文將此現(xiàn)象定義為進氣道不起動。

在激波串移動的過程中，造成隔離段內(nèi)上表面壓力增大，當進氣道不起動時，隔離段內(nèi)整個上表面壓力值增大。由于楔面轉(zhuǎn)角處出現(xiàn)膨脹波，導致壓力下降，所以圖3中進氣道起動狀態(tài)下，喉道處壓力會突然減小。

圖3 條件Q= 0.5, M= 4.5, A=0下上表面壓力分布

由此，本文通過上表面壓力值變化趨勢來判別進氣道起動/不起動。最終從數(shù)值計算的結(jié)果中共獲得321條樣本數(shù)據(jù)，進氣道起動樣本數(shù)據(jù)共161條、不起動樣本數(shù)據(jù)160條，每條數(shù)據(jù)共包含576個壓力數(shù)據(jù)值。

1.2 數(shù)據(jù)集描述

定義該數(shù)據(jù)集為某二元進氣道/隔離段上表面壓力數(shù)據(jù)集，數(shù)據(jù)集中包括321條進氣道起動/不起動樣本，576個屬性(即指CFD網(wǎng)格點上576個壓力分布點)，樣本類別為起動/不起動，屬性值為數(shù)值形式。每個樣本表示成一個序列向量Xi，i的取值為1～321。

Xi=(label:x1x2…xj…x576)

其中： label是樣本類別標簽，label∈{+1,-1}，+1表示進氣道起動樣本，-1表示進氣道不起動樣本。

2 特征選擇算法

2.1 基本算法簡介

Relief算法[13]設(shè)計了一個“相關(guān)統(tǒng)計量”來度量特征的重要性，主要是通過特征間的距離來計算相關(guān)統(tǒng)計量，認為特征的分類能力越強，越能聚集同類樣本，疏遠異類樣本。有樣本集{(x1,y1), (x1,y1),…, (xn,yn)}，對于一個樣本xi，在同類中最鄰近的樣本xi,nh，異類中的鄰近樣本xi,nm，相關(guān)統(tǒng)計量對應(yīng)于屬性j的分類為：

(1)

多次隨機抽取樣本xi，對基于這些樣本得到的估計結(jié)果進行平均，就得到各個特征的相關(guān)統(tǒng)計分量，分量值越大，則對應(yīng)特征越好。

SVMREF算法[13-15]使用遞歸特征去除(Recursive Feature Elimination, REF)算法來保證在特征排序的過程中優(yōu)先保留優(yōu)化特征子集。在對特征進行排序時，一般可以使用支持向量機(SVM)判別函數(shù)中的權(quán)值來計算特征排序準則。每次迭代消去排序準則最小的特征，直到所有特征被刪除，越是重要的特征越是晚被刪除。

SVM分類器的權(quán)值w計算方法為公式(1)，m是支持向量的個數(shù)，αi是非零拉格朗日乘子，xi為訓練樣本，yi為對應(yīng)類別標簽。

(2)

常見的排序準則[16]有如下兩種，公式(2)為線性核函數(shù)SVM分類器的排序準則，wh為消去第h個特征后的權(quán)重；非線性核函數(shù)SVM分類器對應(yīng)的排序為公式(3), 其中α是所有非零拉格朗日乘子，H是一個矩陣，其元素為yiyjK(xi,xj)，H(-h)為消去第h個特征后的矩陣，K為核函數(shù)。

DJ(h)=(wh)2

(3)

(4)

粒子群優(yōu)化算法(PSO:Particle Swarm Optimization)是通過群體中個體之間的協(xié)作和信息共享來尋找最優(yōu)解。初始化一群隨機粒子(隨機解)，然后再每一次的迭代中，粒子通過跟蹤兩個極值(個體最優(yōu)和種群最優(yōu))來更新自己。粒子速度與位置更新公式[17]如下：

vi(t+1)=w·vi(t)+c1rand()·(pj(t)-xi(t))+

c2rand()·(pg(t)-xi(t))

(5)

xi(t+1)=xi(t)+vi(t+1)

(6)

其中：vi是粒子速度，t是迭代次數(shù)，w是慣性權(quán)重，rand()是(0,1)間的隨機數(shù)，Pj是個體最優(yōu)，Pg是種群最優(yōu)，c1和c2是加速系數(shù)。

2.2 基于Relief組合式特征選擇算法

2.2.1 Relief-Corre方法

Relief-Corre方法首先運行Relief算法，得到每個特征的權(quán)重wi，將取值小于設(shè)定無關(guān)閾值的特征刪除，在采用相關(guān)系數(shù)(Correlation coefficient)評估兩兩特征間的冗余度。若相關(guān)度小于設(shè)定的閾值，就刪除權(quán)重較小的一個屬性。具體算法流程如下[8]：

Step1：輸入數(shù)據(jù)集D，特征集S={si,i=1,2,…,n}，Relief無關(guān)閾值σ，相關(guān)度閾值α。

Step2：運行Relief，得到降序排列的特征權(quán)重W={wi,i=1,2,…,n}。

Step3：刪除wi>σ的特征，得到特征集S`={si,i=1,2,…,m}和數(shù)據(jù)集D`。

Step4：fori=1:m-1

Forj=i+1:m

計算兩個特征間的相關(guān)度CORR(i,j)

If(CORR(i,j)>=α)，刪除特征si

Step5：得到最優(yōu)的特征子集。

2.2.2 Relief-SVMRFE方法

Relief-SVMRFE方法仍然在第一步采用Relief算法去掉無關(guān)特征，在第二步中就采用序列后向搜索(Sequential Backward Search,SBS)方法SVMRFE算法去掉冗余特征。具體流程如下：

Step1：輸入數(shù)據(jù)集D，特征集S={si,i=1,2,…,n}，relief無關(guān)閾值σ。

Step2：運行Relief，得到降序排列的特征權(quán)重W={wi,i=1,2,…,n}，刪除wi>σ的特征，得到特征集S`={si,i=1,2,…,m}和數(shù)據(jù)集D`。

Step3：運行SVMRFE算法。

Step4：得到最優(yōu)的特征子集。

2.2.3 特征選擇算法Relief-PSO-SVMRFE

第一階段使用Relief過濾掉權(quán)重小于無關(guān)閾值的特征，在第二階段一分類器的準確率作為自己的評估標準，采用PSO逐步尋找最優(yōu)特征子集，刪除冗余特征。算法描述如下：

Step1：輸入數(shù)據(jù)集D，特征集S={si,i=1,2,…,n}，relief無關(guān)閾值σ。

Step2：Relief剔除無關(guān)特征。

Step3：隨機產(chǎn)生N個長度為P的初始粒子種群，即特征子集。

Step4：利用支持向量機的精度來評估當前粒子的優(yōu)劣性。

Step5：更新局部最優(yōu)和全局最優(yōu)，更新每個粒子的位置。

Step6：轉(zhuǎn)到Step4，直到達到最大迭代數(shù)。

Step7：輸出最優(yōu)特征子集。

3 實驗與分析

3.1 算法結(jié)果分析

基于MATLAB環(huán)境下調(diào)用fitcsvm函數(shù)完成上述原始SVMRFE算法，其中核函數(shù)選擇線性核函數(shù)，懲罰系數(shù)C=∞，優(yōu)化方法選擇序列最小優(yōu)化SMO，并運行SVMRFE算法及第二節(jié)中基于Relief算法進行特征選擇，并對分類準確率和耗時記錄分析。

在Relief階段每次隨機取值300次，共計算3次，統(tǒng)計權(quán)重W的平均值，無關(guān)閾值設(shè)置為0.5，該階段能刪除約69.5%的特征。Relief-Corre方法的相關(guān)閾值設(shè)置為0.8和0.9，Relief-PSO-SVMRFE中c1=c2=2,w=0.9，種群大小為50，迭代次數(shù)為200。首先在整個數(shù)據(jù)集上測試，重復(fù)5次統(tǒng)計最優(yōu)特征子集和平均耗時(由于PSO隨機產(chǎn)生種群的原因，5次結(jié)果中出現(xiàn)3個最優(yōu)特征子集)。對特征子集進行5次5折交叉驗證，給出SVM和樸素貝葉斯的分類錯誤率。

表3 各算法特征選擇結(jié)果耗時與SVM和NBC分類誤差

表3中顯示，所有算法都能選擇出2個或3個特征，Relief-Corre方法耗時最少，Relief-PSO-SVMRFE耗時相當長，而Relief-SVMRFE方法對于SVMRFE在效率上提高了約3.1倍。樸素貝葉斯分類時，每個特征組對應(yīng)的誤差不大，而在SVM分類時，Relief-Corre方法對應(yīng)分類誤差明顯較大，其他3個算法對應(yīng)分類誤差區(qū)別不明顯，但Relief-PSO-SVMRFE方法不能一次就選擇出唯一準確的結(jié)果。通過特征選擇過程耗時與特征子集分類誤差判斷，Relief-SVMRFE方法的性能明顯優(yōu)于其他方法。

3.2 結(jié)果應(yīng)用分析

在進氣道起動預(yù)測中，一般會選擇2個最優(yōu)、最敏感的壓力傳感器放置位置，即選擇兩個特征。但3.1節(jié)中通過Relief-SVMRFE選擇出的[406,505,506]特征并不能直接運用于實際應(yīng)用中，因此對[406,505]和[406,506]兩組特征進行討論，刪除冗余特征；由于實際應(yīng)用中壓力傳感器測到的數(shù)據(jù)存在誤差，再對得到的分類面的魯棒性進行分析。

通過五次五折對[406,505]和[406,506]特征SVM分類面判別誤差進行計算，分別為0.0031和0.00124，特征505,506對應(yīng)Relief權(quán)重分別為16.1077和15.9996，所以506特征是優(yōu)于505特征的，因此刪除505特征。令第406與506個壓力測試點分別為P1,P2，其中P1是進氣道內(nèi)某個位置上的壓力測試點，P2為靠近喉道位置某個壓力測試點?；赑1和P2兩點訓練分類面，分類面為p2-4.7058*p1+3388.6156=0，兩條邊界分別為p2-4.7058*p1+9763.7793=0(起動)和p2-4.7058*p1+2986.5481=0(不起動)。數(shù)據(jù)分類情況如圖4所示。

圖4 SVM訓練的分類結(jié)果

由于壓力測試點位置的進度已經(jīng)達到0.1 mm，在實際應(yīng)用中，不能保證壓力傳感器測試到的壓力值就是P1,P2點的值，因此應(yīng)檢驗分類面的魯棒性能。結(jié)合點的位置坐標與Relief權(quán)重因素，選擇相似的點進行考察，這些點對應(yīng)的壓力值代表在實際應(yīng)用中壓力傳感器可能會測到的壓力值，獲得兩組考察壓力點405～407和503～509，Relief權(quán)重和壓力點坐標位置如表4。一共可以獲得21個組合，用已得到的分類面對這些數(shù)據(jù)進行分類，表5中所有的分類誤差均不大于0.001。盡管壓力傳感器沒有測到P1和P2點的壓力，實際測到的壓力值具有一定的偏差，該分類面依然可以較準確的將進氣道起動/不起動分類開。

表4 兩組壓力點Relief權(quán)重和坐標位置

表5 特征子集對應(yīng)分類誤差

4 結(jié)果

本文主要解決高超聲速進氣道問題，使用Relief算法和SVMRFE算法來保證特征選擇的效率與準確性。試驗中，特征選擇使用到了Relief-Corre方法、Relief-SVMRFE方法、算法Relief-PSO-SVMRFE和SVMRFE算法。Relief-SVMRFE方法和常用的分析進氣道起動方法(SVMRFE算法)在獲得同樣結(jié)果的基礎(chǔ)上，效率上有了很大的提高，而其他的基于relief算法在該方面應(yīng)用并沒有明顯的優(yōu)勢；基于文中最優(yōu)特征子集[406,506]得到分類面具有較高分類性能和魯棒性能，有實際應(yīng)用的意義。最后特征選擇結(jié)果中特征數(shù)量是否可以進一步減少，將在接下來的研究中繼續(xù)展開討論。

參考文獻：

[1] 吳穎川, 賀元元, 賀 偉,等. 吸氣式高超聲速飛行器機體推進一體化技術(shù)研究進展[J]. 航空學報, 2015, 36(1):245-260.

[2] Heiser W, Pratt D, Daley D, et al. Hypersonic airbreathing propulsion[M]. American Institute of Aeronautics and Astronautics, 2010.

[3] Cox C, Lewis C, Pap R, et al. Prediction of unstart phenomena in hypersonic aircraft[M]. International Aerospace Planes and Hypersonics Technologies,1995.

[4] Yu D, Chang J, Bao W, et al. Optimal classification criterions of hypersonic inlet start/unstart[J]. Journal of Propulsion & Power, 2007, 23(2):310-316.

[5] 常軍濤, 于達仁, 鮑 文. 基于支持向量機的高超聲速進氣道不起動模式細化分類研究[A].中國工程熱物理學會2008年熱機氣動熱力學與流體機械學術(shù)會議[C]. 2008.

[6] Kononenko I. Estimating attributes: Analysis and extensions of RELIEF[J]. Machine Learning:ECML-94, 1994, 784:171-182.

[7] Kohavi R, John G H. Wrappers for feature subset selection[J].Artifical Intelligence,1997,97(1):173-324.

[8] 張麗新, 王家廞, 趙雁南,等. 基于Relief的組合式特征選擇[J]. 復(fù)旦學報(自然科學版), 2004, 43(5):893-898.

[9] 吳紅霞, 吳 悅, 劉宗田,等. 基于Relief和SVM-RFE的組合式SNP特征選擇[J]. 計算機應(yīng)用研究, 2012, 29(6):2074-2077.

[10] 杜洪波, 董文娟. Relief-PSO混合算法在基因微陣列特征選擇中的應(yīng)用[J]. 沈陽工程學院學報(自然科學版), 2016, 12(3):267-271.

[11] Hank J, Murphy J, Mutzman R. The X-51A Scramjet Engine Flight Demonstration Program[A]. AIAA International Space Planes and Hypersonic Systems and Technologies Conference[C]. 2008.

[12] 余華兵，唐士廷，胡仁喜. FLUENT14.5流場分析從入門到精通[M]. 北京:機械工業(yè)出版社, 2014.

[13] 周志華. 機器學習 : = Machine learning[M]. 北京:清華大學出版社, 2016.

[14] Cortes C, Vapnik V. Support-Vector Networks[J]. Machine Learning, 1995, 20(3):273-297.

[15] 鄧乃揚, 田英杰. 支持向量機:理論、算法與拓展[M]. 北京:科學出版社, 2009.

[16] 盧運梅. SVM-RFE算法在數(shù)據(jù)分析中的應(yīng)用[D]. 吉林:吉林大學, 2009.

[17] 齊立磊, 陸 劍. 基于粒子群算法改進SVM的滑?？刂蒲芯縖J]. 計算機測量與控制, 2014(10):3230-3232.