基于小波包能熵譜的爆震特征頻帶及強度分析

薛勁梓，胡春明，劉娜，趙英博，張建飛

(1.天津大學內燃機研究所，天津，300072；2.天津大學機械工程學院，天津，300350)

點燃式發(fā)動機的爆震燃燒屬于不正常的燃燒現(xiàn)象，爆震會使燃燒室內的壓力曲線在上止點附近呈現(xiàn)鋸齒狀，爆震程度越大，這種形狀越發(fā)明顯、振蕩幅度也越大。為了監(jiān)測發(fā)動機的運行狀態(tài)，需要準確評估爆震強度。相較于缸體振動信號的檢測方式，基于缸壓的爆震檢測指標是一種較為準確的判別方法。

爆震會激發(fā)燃燒室諧振，使缸壓中的高頻成分升高，所以，對于缸壓信號，使用4～20 kHz 高通濾波器提取其中的高頻分量，將濾波后高頻振蕩壓力的最大幅值(the maximum amplitude pressure oscillation，MAPO)作為爆震強度的評價指標[1-2]，對爆震有一定的指示效果。但是，當缸壓信號因外部電信號干擾出現(xiàn)毛刺、噪點時，提取出的MAPO 將會因為干擾存在而不可靠，導致錯誤的檢測結果。為了克服這一缺點，Siemens VDO爆震算法將低通濾波后的最大壓力點作為爆震始點，在爆震始點前后部分(參考窗口和爆震窗口)信號進行積分，取爆震窗口與參考窗口高頻信號的積分之比作為爆震因子(knock factors，KF)[3-4]，若KF大于一定的閥值則認為是爆震循環(huán)。Siemens VDO算法雖然在一定程度上能消除信號出現(xiàn)毛刺的影響，但是其結果依賴于爆震始點、爆震窗口寬度等關鍵參數(shù)的人為主觀設置。

近年來，小波包分析作為一種精細的時頻分析工具應用在特征信號的提取上，其具有削弱干擾，突出信號特征的優(yōu)點。史歷程等[5]將小波變換和熵原理結合，利用小波分解計算出小波能量熵，得到在不同尺度內的能量分布情況，從而提取出傳感器故障信號的特征。與小波分解相比，小波包變換能夠提供更高的分辨率，能夠對信號的低頻和高頻部分進行分解，從而得出多級小波包分解樹的結構。楊興國[6]利用小波包變換與雙譜估計相結合的方法提取發(fā)動機聲信號特征向量作為輸入，通過機器學習中支持向量機的算法對模型訓練來識別發(fā)動機的異響。劉成材等[7]在基于加速度傳感器檢測發(fā)動機爆震的研究中，對缸體振動信號進行了小波包分解并重構信號，從重構后的振動信號提取出爆震特征，尤其增強了發(fā)動機輕微爆震檢測性能。

本文作者使用小波包變換的方法處理缸內燃燒壓力的信號，提取其中高頻壓力分量，計算有效頻帶的范數(shù)能量熵并得到小波包能熵譜。然后，從分析大量爆震試驗所獲得的樣本統(tǒng)計值出發(fā)，探究不同爆震強度樣本下的小波包能熵譜分布情況。最后基于小波包能熵譜的分析提出一種新的爆震強度指標，驗證其可行性和有效性，并利用ROC曲線估計出選定工況下的爆震閾值。

1 小波包變換分析方法

1.1 小波包分解基本原理

在小波變換的基礎上，小波包變換可以進一步分解高頻區(qū)域中信號的詳細信息[8]。在每一級信號分解時，除了對低頻子帶進行分解，也對高頻子帶進行分解。隨著小波包分解的級數(shù)增加，其在頻域上的分辨率就越高，因此，可以獲取信號較高頻率的信息分量。

小波函數(shù)通過尺度變量和平移變量平移、伸縮獲得不同的頻域和時間的位置，而小波包Ψji，k是小波函數(shù)φj的線性組合[9]：

式中：i為調制參數(shù)；j為尺度參數(shù)；k為轉換參數(shù)；t為時間參數(shù)。小波函數(shù)φj由式(2)和式(3)遞歸計算而得。

其中：h(k)和g(k)分別為與尺度函數(shù)和母小波函數(shù)關聯(lián)的正交鏡像低通濾波器和高通濾波器。在母小波函數(shù)中，Daubechies小波族函數(shù)具有較好的可逆性和正交性[10]，令φ1(t)=φ(t)建立母小波函數(shù)，得到滿足可逆性和正交性的小波包函數(shù)Ψji，k。根據(jù)式(4)獲得在分解等級z上的第u個節(jié)點處的小波包系數(shù)。由小波包的重構系數(shù)cuz，k(t)和小波包函數(shù)Ψji，k的線性組合得到小波包分量信號yuz(t)由式(5)給出[11]。

1.2 小波包能熵譜

當系統(tǒng)偏離正常狀態(tài)時，信號的頻率和幅值都將產(chǎn)生較大變化，熵值也會產(chǎn)生相應的變化[12]，并且將小波包和能量熵結合能更好地便于信號特征的提取。范數(shù)熵(norm entropy)是小波包能量熵的一種，設時間序列長度為n，對于單變量時間序列xt的范數(shù)熵可由式(6)計算可得：

當p=2時，定義的范數(shù)熵可以用來代表信號的能量[13]。據(jù)此計算小波包分解和重構后各頻帶分量的時間離散序列yij的范數(shù)熵，選取高頻分量中能量熵較高的前m個有效分量，忽略能量熵較低的分量，組成小波包能熵譜。小波包能熵譜中各頻帶的總能量熵E及各頻帶的能量占比νm由式(7)和(8)可得。

2 直噴煤油發(fā)動機爆震試驗設計及數(shù)據(jù)處理

2.1 臺架試驗

爆震試驗在1臺自主設計開發(fā)的四沖程單缸航空煤油發(fā)動機上進行，煤油發(fā)動機的各項參數(shù)如表1所示。試驗臺架包括單缸試驗煤油發(fā)動機、電力測功機、燃油供給系統(tǒng)、冷卻系統(tǒng)、燃燒數(shù)據(jù)采集與分析系統(tǒng)、發(fā)動機電子控制系統(tǒng)以及上位機標定與控制系統(tǒng)組成。

表1 單缸煤油發(fā)動機參數(shù)Table 1 Parameters of single-cylinder kerosene engine

試驗中用于采集和分析缸內燃燒壓力的設備有Kistler 6113C 型火花塞式壓力傳感器、Kistler 5011B型電荷放大器和德維創(chuàng)公司的DEWESoft X3燃燒分析儀。外時鐘觸發(fā)采樣，外時鐘頻率與轉速有關，為獲取更加準確的爆震高頻信號，缸壓信號采集分辨率設置為0.1°。為了使爆震強度能夠隨著點火提前角的變化而發(fā)生顯著變化，選取發(fā)動機中高速小負荷工況進行試驗，試驗工況參數(shù)如表2所示。在爆震試驗過程中，通過上位機ECU標定軟件增加點火提前角，以獲得不同強度的爆震工況。

表2 爆震試驗工況參數(shù)Table 2 Parameters of knock condition for test

2.2 小波包提取爆震壓力高頻分量及其能熵譜

爆震發(fā)生時，燃燒室的末端混合氣自燃并釋放出大量的能量，能量的快速釋放將激發(fā)一系列特定的共振模式，即存在多個爆震頻率中心。燃燒室中的共振頻率取決于燃燒室的形狀、尺寸和局部聲速，可通過求解如下理論公式來估算[14]：

一天，我和表姐帶曉玥出門逛街。路過的天橋上有一個滿身泥巴的乞丐。曉玥老遠就抓著我們的衣角，躲在我們身后，盡可能地遠離乞丐。我們以為她嫌棄這個臟兮兮的乞丐，特意停下來給她2元硬幣，要她送到乞丐的碗里去。她怯怯地走到乞丐跟前，小心地把硬幣放進去，然后小跑著離開。我不禁感嘆：孩子的童真和善良都被這個功利的社會吞噬了！

式中：ρη，δ為振動模態(tài)因子，下標η為切向模態(tài)數(shù)，下標δ為徑向模態(tài)數(shù)；B為氣缸直徑；C為氣體聲速。本文使用的發(fā)動機內徑B=100 mm，氣體聲速C=1 000 m/s。表3所示為爆震發(fā)生時不同共振模態(tài)下的預測頻率。

表3 圓柱形燃燒室不同振動模式下的理論頻率Table 3 Theoretical frequency under different vibration modes of cylindrical combustion chamber

缸壓信號的采樣頻率為210 kHz，在高頻分量的提取過程中使用具有正交性和緊支撐性的db5小波函數(shù)[9]，本文依據(jù)表3中的爆震頻帶范圍(5.86～20.4 kHz)使用6級小波包分解。經(jīng)過每層分解，高頻濾波部分在下采樣后頻譜發(fā)生翻轉，導致頻率順序與節(jié)點順序不同[15]，小波包分解樹如圖1所示。

圖1 六級小波包分解樹Fig.1 Six-level wavelet packet decomposition tree

由于能量熵較小的小波包分量成分容易受到信號噪聲的干擾，忽略26.3 kHz 以上的頻率成分。根據(jù)表3中的理論頻率選取6級小波包分解后的有效高頻分量：d1(4.6～6.6 kHz)，d2(9.9～13.2 kHz)，d3(6.6～9.9 kHz)，d4(19.8～23.1 kHz)，d5(13.2～16.5 kHz)和d6(16.5～19.8 kHz)。取這6 個高頻頻帶信號的范數(shù)能量熵構成小波包能熵譜。

圖2所示為某一爆震循環(huán)的缸壓信號分解得到的不同頻帶信號的時域圖和頻域圖。圖2(a)中每個頻帶信號是通過對原始缸壓信號執(zhí)行6級小波包分解和重構而得，圖2(b)中頻譜由快速傅里葉變換(FFT)獲得。由圖2(b)可知，小波包分解能夠有效提取相應頻帶的信號。圖3所示為由此爆震循環(huán)的高頻頻帶d1～d6能量熵構建的能熵譜，揭示出爆震循環(huán)中各高頻頻帶的能量分布情況。

圖2 六級小波包分解信號及其頻率譜Fig.2 Six-level wavelet packet decomposition signal and its frequency spectrum

圖3 缸壓信號高頻頻帶能量熵譜圖Fig.3 Energy entropy spectrum diagram of high frequency band of cylinder pressure signal

3 基于小波包能熵譜的爆震分析

從試驗獲得的4組樣本中，利用所采集到的循環(huán)缸壓曲線基于Siemens VDO 爆震算法計算出爆震因子KF，然后，根據(jù)爆震因子KF的大小將爆震強度分為4個等級：嚴重爆震、中等爆震、輕微爆震、未爆震，如表4所示。由于爆震強度在循環(huán)間表現(xiàn)出隨機特性，本文將從所得樣本的統(tǒng)計值出發(fā)，探究爆震現(xiàn)象對小波包能熵譜中高頻頻帶能量分布的影響，并提出一種新的表征爆震強度的指標。

表4 爆震試驗不同樣本分類情況Table 4 Classification of different samples in knock test個

3.1 爆震強度對能熵譜特征的影響

圖4所示為樣本中所有不同爆震強度分類下的各頻帶能量占比的平均統(tǒng)計值。從能量成分上看，d1，d2，d3和d5是爆震頻率的主要成分，即爆震能量主要集中在4.6～16.5 kHz頻帶內。從能量變化幅度上看，隨著爆震強度的增加，頻帶d2(9.9～13.2 kHz)和d5(13.2～16.5 kHz)的能量占比不斷增大，d3(6.6～9.9 kHz)和d6(16.5～19.8 kHz)的能量占比卻不斷下降，d1(4.6～6.6 kHz)和d4(19.8～23.1 kHz)的能量占比基本不變，不受爆震強度的影響。這說明不同頻帶能量對爆震的敏感性不同，d2和d5頻帶的能量增長幅度較大，受爆震影響最為顯著，d1和d4次之。

圖4 不同爆震強度樣本下各頻帶能量占比的平均統(tǒng)計值Fig.4 Average statistical value of energy ratio of each frequency band under different knock intensity samples

圖5所示為樣本1～3中各頻帶能量熵值的變化情況。從圖5可知：無論何種爆震水平，爆震能量均主要分布在9.9～16.5 kHz的爆震頻帶。在同一點火提前角下，隨著爆震等級的提升，各爆震頻帶的能量熵不斷增加。點火提前角的不同，并不會改變頻帶能量的分布情況，這說明爆震特征頻率與點火提前角無關。隨著點火提前角增大，各頻帶具有更大的能量熵，這是因為點火角增加使得末端混合氣在壓縮行程末期溫度和壓力不斷增加，自燃后瞬間釋放出更多的熱量，缸內燃燒壓力急劇升高，在燃燒室內產(chǎn)生壓力振蕩沖擊波導致出現(xiàn)更多的高頻壓力能量[16]。

圖5 不同點火提前角下爆震頻帶能量熵的樣本統(tǒng)計值Fig.5 Sample statistical value of energy entropy of knock frequency band under different ignition advance angles

各頻帶的能量反映出爆震強度的等級，這主要是因為爆震發(fā)生的瞬間，可燃混合氣自燃釋放的熱量越多，會產(chǎn)生更多的高頻振蕩壓力能量，從而出現(xiàn)不同強度等級的爆震現(xiàn)象。值得注意的是，在d2(9.9～13.2 kHz)和d5(13.2～16.5 kHz)頻帶內的能量在發(fā)生嚴重爆震強度時明顯增加，此頻帶與燃燒室一階徑向，二階、三階切向的理論共振頻率相對應。

圖6 不同爆震強度循環(huán)的缸壓信號及其爆震頻帶能熵譜Fig.6 Cylinder pressure signal of different knock intensity cycles and its knock frequency band energy entropy spectrum

多個循環(huán)的樣本統(tǒng)計值和單個循環(huán)樣本的小波包能熵譜均能反映出此款煤油發(fā)動機的爆震特征頻率為9.9～16.5 kHz；沒有爆震發(fā)生時能熵譜中的能量分布較為均勻，爆震發(fā)生時能量分布發(fā)生明顯變化，爆震特征頻帶9.9～16.5 kHz所具有的能量熵顯著增加。

3.2 基于譜分析的爆震強度表征

不同的循環(huán)因為爆震強度不同，各個頻帶中能量熵最大值也發(fā)生相應變化?；诖爽F(xiàn)象，本文提出將小波包能熵譜中最大能量熵(the maximum of entropy,MOE)作為爆震發(fā)生的強度評價指標，即MOE 表征的是爆震時壓力信號分布在各個高頻頻帶中能量的最大值。為了檢驗爆震指標評價的效果，計算樣本1～3的各個循環(huán)下的爆震指標。與使用Siemens VDO算法得到的爆震指標KF進行相關性檢驗[17]，并用最小二乘法線性擬合得到分布圖，如圖7所示。從圖7可見：在不同爆震工況循環(huán)下，KF 和本文提出的MOE 指標的Pearson 相關系數(shù)r均在0.91之上，說明爆震指標KF和MOE指標具有高度的線性正相關關系，所以，MOE 同樣能對爆震強度進行合理評價。

圖7 樣本1～3中各循環(huán)的MOE和KF相關性檢驗圖Fig.7 Correlation test chart of MOE and KF for each cycle in samples 1-3

通過對爆震樣本數(shù)據(jù)分析，爆震頻帶能量熵通常在頻帶d2(9.9～13.2 kHz)和d5(13.2～16.5 kHz)處取得最大值，即所提出的爆震指標反映的是與爆震相關特征頻帶的能量。相比于Siemens VDO 算法的準確度由爆震窗口寬度和起點確定，本文基于缸壓信號提出的爆震表征指標MOE可以克服這一不足，且計算簡便，對于點燃式發(fā)動機爆震強度評價具有一定的參考價值。

圖8所示為點火提前角從27°增加到36°的4個樣本的爆震強度評價指標MOE的分布圖。從圖8可以看出：從無爆震工況到爆震工況，評價指標MOE 能清晰地將不同爆震強度等級的工況區(qū)分出來，點火提前角增加至36°的過程中，爆震循環(huán)的比例逐漸提高。樣本標準差(σ)逐漸增大，MOE分布更加離散，即爆震發(fā)生的頻率增加；樣本均值(μ)也不斷增大，即高頻壓力分量的能量水平逐漸提高。這符合點燃式發(fā)動機爆震指標分布規(guī)律[18]。處于爆震臨界狀態(tài)的循環(huán)位于圖8中紅色方框內，此區(qū)域內的MOE 在0.3～0.4 之間，既有無爆震循環(huán)，又有輕微爆震循環(huán)。

圖8 無爆震至爆震工況中各循環(huán)的MOE分布圖Fig.8 MOE distribution of each cycle from knock-free to knock conditions

劃定合理的閾值會減少爆震檢測的錯誤發(fā)生率，本文通過ROC 曲線評估不同閾值對爆震檢測的準確率。從0.3開始每隔0.05設置爆震閾值，檢測的負正類率(false positive rate，F(xiàn)PR)和真正類率(true positive rate，TPR)會隨之改變，觀測結果如圖9所示。從圖9可見：爆震閾值設置在0.355 時曲線出現(xiàn)轉折，識別出真正類的比例達到94.22%，設置的爆震閾值減小會使負正類率FPR陡然上升，正確檢測概率TPR 趨緩，故將MOE 設置為0.355時，在當前工況下，爆震檢測有較高的正確率。需要指出的是，在不同工況下，爆震閾值有所不同，需要經(jīng)過試驗標定確定合適的閾值。

圖9 不同爆震閾值對應的正誤率ROC曲線Fig.9 ROC curve of correct error rate correspond-ing to different knock thresholds

4 結論

1)通過小波包分解和重構算法可以提取爆震發(fā)生時不同頻帶下的高頻壓力信號，進而獲得小波包能熵譜。該方法為分析爆震各頻帶的能量變化提供了有效的途徑。

2)爆震能量主要集中在4.6～16.5 kHz高頻壓力內。當爆震強度增加時，小波包各頻帶分量的能量熵值及其分布會受到影響，9.9～13.2 kHz 和13.2～16.5 kHz 頻帶能量熵增長迅速，隨著爆震的發(fā)生而出現(xiàn)明顯的變化。

3)本文提出的爆震評價指標MOE 與Siemens VDO 算法計算得到的爆震因子(KF)有較好的相關性，可以作為爆震強度的評價指標。在此款煤油發(fā)動機上，在試驗工況下，MOE處于0.3～0.4為爆震臨界區(qū)間，利用ROC 曲線估計時，將爆震閾值設置為0.355作為爆震監(jiān)測條件更合適。