應(yīng) 銘，馮國(guó)勝 ，賈素梅 ，霍肖楠 ，馬春庭

(1. 北京理工大學(xué) 機(jī)械與車(chē)輛學(xué)院，北京 100081；2. 石家莊鐵道大學(xué) 機(jī)械工程學(xué)院，河北 石家莊 050043； 3. 石家莊鐵道大學(xué) 河北省工程機(jī)械動(dòng)力與傳動(dòng)控制重點(diǎn)實(shí)驗(yàn)室，河北 石家莊 050043； 4. 河北軍濤科技有限公司，河北 石家莊 050000)

柴油機(jī)已廣泛應(yīng)用于工業(yè)、農(nóng)業(yè)、交通運(yùn)輸業(yè)和國(guó)防建設(shè)事業(yè)，是為多種機(jī)械提供動(dòng)力的核心構(gòu)件．柴油機(jī)故障會(huì)造成經(jīng)濟(jì)損失和人員傷亡．傳統(tǒng)故障診斷的過(guò)程一般為：采集數(shù)據(jù)、分析數(shù)據(jù)和檢出故障．在實(shí)際應(yīng)用中，分析數(shù)據(jù)的速度往往要慢于采集數(shù)據(jù)的速度．這導(dǎo)致了故障的發(fā)現(xiàn)存在一定滯后性，發(fā)動(dòng)機(jī)在一定時(shí)間內(nèi)帶故障繼續(xù)運(yùn)行，且兩次采集的數(shù)據(jù)之間存在間斷，可能導(dǎo)致間歇性故障漏診斷．這種現(xiàn)象在需要高采樣頻率的診斷信號(hào)中尤為顯著，采樣頻率越高，短時(shí)間內(nèi)獲得的數(shù)據(jù)量就越大，處理大量數(shù)據(jù)所需時(shí)間也越多．

柴油機(jī)缸蓋振動(dòng)信號(hào)包含豐富的狀態(tài)信息[1-2]，且傳感器安裝容易、成本較低，是理想的監(jiān)測(cè)信號(hào)來(lái)源．但柴油機(jī)缸蓋振動(dòng)的激勵(lì)源較多，要獲取其完整頻譜，需在較高采樣頻率下采集[3-4]，且振動(dòng)信號(hào)包含燃燒噪聲、活塞敲擊噪聲、進(jìn)/排氣門(mén)落座噪聲、噴油泵噪聲和正時(shí)齒輪噪聲等多種噪聲[5]．噪聲不僅會(huì)干擾缸內(nèi)壓力、燃燒狀態(tài)等工作參數(shù)的提取[6-7]，也會(huì)影響故障的檢出率[8]．

經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)和小波變換(WT)是處理非平穩(wěn)信號(hào)的常用手段[9-10]．EMD缺乏嚴(yán)格的數(shù)學(xué)理論模型支持，各個(gè)固有模態(tài)函數(shù)(IMF)分量不正交，信號(hào)重構(gòu)時(shí)易產(chǎn)生混疊模態(tài)[11]，且其在運(yùn)算時(shí)間上不穩(wěn)定，不適合在嵌入式平臺(tái)上移植．小波變換已被廣泛運(yùn)用于分析具有瞬態(tài)、非穩(wěn)態(tài)或時(shí)變特性的信號(hào)[11-12]．但傳統(tǒng)小波變換實(shí)現(xiàn)算法的運(yùn)算速度較慢，難以滿(mǎn)足高采樣頻率下的實(shí)時(shí)降噪需求．提升小波也稱(chēng)為第二代小波，同時(shí)具有實(shí)時(shí)性強(qiáng)和原位運(yùn)算的優(yōu)點(diǎn)[13-14]，但以往對(duì)小波降噪實(shí)時(shí)性的評(píng)價(jià)依賴(lài)統(tǒng)計(jì)小波變換運(yùn)算過(guò)程中乘法和加法的次數(shù)[15]，并未考慮工程應(yīng)用中算法運(yùn)行載體—數(shù)字信號(hào)處理器(DSP)的占用時(shí)間．影響了對(duì)小波降噪實(shí)時(shí)性評(píng)價(jià)的

準(zhǔn)確性，使得評(píng)價(jià)結(jié)果不能真實(shí)反映工程應(yīng)用情況．且小波降噪?yún)?shù)大多僅根據(jù)降噪效果進(jìn)行選擇，并未考慮實(shí)時(shí)性的影響．目前，從實(shí)時(shí)性角度對(duì)小波進(jìn)行的研究較少，有必要構(gòu)造高實(shí)時(shí)性的小波降噪系統(tǒng)進(jìn)而提高檢測(cè)儀器的診斷實(shí)時(shí)性．

筆者以工程應(yīng)用中的算法載體DSP為研究對(duì)象，分析了小波算法的實(shí)時(shí)性評(píng)價(jià)方法及考慮實(shí)時(shí)性的小波參數(shù)優(yōu)選方法，以期解決工程實(shí)踐中實(shí)時(shí)小波降噪系統(tǒng)的構(gòu)建問(wèn)題．

1 小波變換概述

對(duì)一能量有限信號(hào)[16]f(t)，滿(mǎn)足

式中L2(R)為實(shí)數(shù)域中的平方可積空間．

記f(t)的傅里葉變換為Ψ(ω)，若Ψ(ω)滿(mǎn)足容許性條件CΨ，即

則稱(chēng)Ψ(t)為一個(gè)基小波或母小波，將基小波Ψ(t )進(jìn)行平移、伸縮可得到一個(gè)小波函數(shù)Ψa,b(t )，即

式中：a為尺度因子；b為平移因子．

則能量有限信號(hào)f(t)的連續(xù)小波變換(CWT)為

式中Ψ*為對(duì)Ψ的共軛運(yùn)算．其逆變換定義為

在實(shí)際應(yīng)用中，所采集的信號(hào)都是離散化的，故需對(duì)連續(xù)小波進(jìn)行離散化．對(duì)于小波函數(shù)(式(3))，假設(shè)a,b∈R，a＞0且a≠1，Ψ(t )滿(mǎn)足容許性條件．將尺度因子a和平移因子b離散化，a以?xún)缂?jí)數(shù)的形式離散化，b以均勻抽樣的方式離散化．

離散化的小波函數(shù)為

式中：j、k為整數(shù)；a0、b0為常數(shù)，且a0≠1．

則離散小波變換(DWT)定義為

其逆變換定義為

式中c為與原始信號(hào)無(wú)關(guān)的常數(shù)．

2 小波變換指令周期模型

指令周期是指DSP取出并執(zhí)行一條指令所需的時(shí)間，指令周期模型即算法執(zhí)行總計(jì)所需的指令周期個(gè)數(shù)．小波變換包括取值、賦值、乘法和加法等指令，在DSP中這些指令耗時(shí)相等，則傳統(tǒng)的忽略賦值和取值指令的統(tǒng)計(jì)方法不適合在DSP中使用．乘積累加運(yùn)算(MAC)的定義為

式中：A、B和C為常數(shù)．

MAC包括取值、賦值、乘法和加法，基于硬件設(shè)計(jì)，DSP可在一個(gè)指令周期內(nèi)完成一次MAC．

2.1 經(jīng)典小波指令周期模型

經(jīng)典小波的實(shí)現(xiàn)依賴(lài)于Mallat算法，該算法可分為延拓、卷積和下采樣3個(gè)步驟[17]，有

式中：n為小波變換序列點(diǎn)數(shù)；m為當(dāng)前分解層數(shù)；cm和dm為第m層分解的近似系數(shù)和細(xì)節(jié)系數(shù)；g(j)和h(j)為低通和高通濾波器的第j個(gè)系數(shù)；l為濾波器長(zhǎng)度．將式(10)改寫(xiě)為乘積累加形式，有

式中：i =1，2，3,…,n ；j=1,2,3,…,l．可知，使用Mallat算法對(duì)點(diǎn)數(shù)為n的序列進(jìn)行小波變換需要(l·n)次MAC和(l2-l)次延拓賦值，故Mallat算法的指令周期模型為

2.2 提升小波指令周期模型

提升小波也稱(chēng)第二代小波，其保留了經(jīng)典小波多分辨率的特性，提升小波的構(gòu)建經(jīng)過(guò)奇偶分裂、預(yù)測(cè)和更新3個(gè)步驟[18]，有

式中：splits(sn)表示將序列sn按位置奇偶分裂；dn-1為奇數(shù)序列；sn-1為偶數(shù)序列；P為預(yù)測(cè)算子；U為更新算子．將式(13)改寫(xiě)為乘積累加形式，有

式中：i=1，2，3…，n/2；j=1，2，3…，p；k=1，2，3…，lj；Pj(k)和Uj(k)分別為預(yù)測(cè)、更新算子第j次預(yù)測(cè)、更新步驟的第k個(gè)系數(shù)；p為預(yù)測(cè)、更新步驟的次數(shù)；lj為預(yù)測(cè)、更新算子第j次預(yù)測(cè)、更新步驟的長(zhǎng)度．可知，一次提升小波變換需要次的MAC．

2.2.1 原位運(yùn)算提升小波指令周期模型

“搭腳手架”實(shí)踐其實(shí)并不難，我們大部分家長(zhǎng)只是沒(méi)有這個(gè)意識(shí)。下面是幾個(gè)我自己在教學(xué)中總結(jié)出來(lái)的主要核心點(diǎn)。

提升小波可以實(shí)現(xiàn)原位運(yùn)算，原位運(yùn)算是指僅使用原序列空間，不需其他輔助空間的運(yùn)算過(guò)程．但原位運(yùn)算可能導(dǎo)致運(yùn)算速度降低．文獻(xiàn)[19]給出了奇偶分裂的原位運(yùn)算方法，該方法所需的最小指令周期為3n2/16．

則原位運(yùn)算提升小波指令周期模型為

2.2.2 非原位運(yùn)算提升小波指令周期模型

奇偶分裂非原位運(yùn)算的快速實(shí)現(xiàn)方法是將原始信號(hào)的奇數(shù)、偶數(shù)序號(hào)放入輔助數(shù)組的前半部分和后半部分，非原位運(yùn)算的奇偶分裂需要n條賦值指令．

則非原位運(yùn)算提升小波指令周期模型為

2.3 多層小波變換指令周期模型

對(duì)于分解層數(shù)為q的多層小波變換，可近似認(rèn)為每層變換的序列點(diǎn)數(shù)為上一層的一半，則q層小波變換的指令周期模型為

式中：q為小波分解總層數(shù)；T(21-k·n)為第k層分解時(shí)所用小波算法的指令周期模型．

將式(12)、(15)和(16)代入式(17)計(jì)算，可知，Mallat算法和非原位運(yùn)算提升小波所需的指令周期與小波變換序列點(diǎn)數(shù)呈正相關(guān)；原位運(yùn)算提升小波所需的指令周期與小波變換序列點(diǎn)數(shù)呈二次相關(guān)．在相同預(yù)測(cè)、更新算子下，當(dāng)小波變換序列點(diǎn)數(shù)超過(guò)8時(shí)，原位運(yùn)算提升小波所需的運(yùn)算時(shí)間將超過(guò)非原位運(yùn)算提升小波；若記為預(yù)測(cè)、更新算子的長(zhǎng)度，則Mallat算法與非原位運(yùn)算提升小波所需運(yùn)算時(shí)間的長(zhǎng)短取決于該小波基函數(shù)對(duì)應(yīng)的濾波器與預(yù)測(cè)、更新算子的長(zhǎng)度．

對(duì)于多層小波變換，Mallat算法和非原位運(yùn)算提升小波第k層分解所需的指令周期為第(k-1)層的一半；原位運(yùn)算提升小波第k層分解所需的指令周期為第(k-1)層的1/4．

3 柴油機(jī)缸蓋振動(dòng)信號(hào)的小波降噪

小波變換的分解層數(shù)、預(yù)測(cè)、更新算子和濾波器的長(zhǎng)度對(duì)小波實(shí)時(shí)性有顯著影響．筆者以CA4DF3-13E3型柴油機(jī)為對(duì)象研究柴油機(jī)實(shí)時(shí)監(jiān)測(cè)系統(tǒng)中小波降噪的合適參數(shù)，表1為柴油機(jī)的主要技術(shù)參數(shù)．

表1 發(fā)動(dòng)機(jī)主要技術(shù)參數(shù) Tab.1 Engine specifications

圖1為CA4DF3-13E3型柴油機(jī)缸蓋振動(dòng)加速度的測(cè)點(diǎn)位置．試驗(yàn)設(shè)置采樣頻率為25kHz，分別采集800、1000、1200、1400、1600、1800、2000和2300r/min空載工況下第1缸缸蓋振動(dòng)加速度，每次 采集5000點(diǎn)數(shù)據(jù)．圖2為低、中和高速3種典型工況下缸蓋振動(dòng)的功率譜密度．

圖1 測(cè)點(diǎn)位置Fig.1 Location of measuring point

小波分解層數(shù)會(huì)影響降噪效果與實(shí)時(shí)性．分解層數(shù)過(guò)多，會(huì)導(dǎo)致原始信號(hào)被過(guò)度降噪失去原有特征，且實(shí)時(shí)性較差；分解層數(shù)過(guò)少，會(huì)導(dǎo)致原始信號(hào)降噪不完全．柴油機(jī)缸蓋振動(dòng)信號(hào)能量主要集中在頻段0～1.6kHz(圖2)．一次小波分解會(huì)將原始信號(hào)等分為高、低兩個(gè)頻帶，圖3為根據(jù)功率譜密度計(jì)算 不同層數(shù)小波分解的低頻能量相對(duì)總能量的占比．

圖2 功率譜密度Fig.2 Power spectral density

圖3 頻帶能量比例Fig.3 Band energy ratio

低頻能量占比在第4層小波分解時(shí)出現(xiàn)顯著下降，說(shuō)明振動(dòng)信號(hào)的主體特征在第4層小波分解時(shí)被劃歸到需降噪的高頻部分．為保留信號(hào)主要特征，防止過(guò)度降噪，選擇小波分解層數(shù)為3層．

小波基函數(shù)決定了Mallat算法中濾波器的長(zhǎng)度和提升小波中的預(yù)測(cè)、更新算子，合適的小波基函數(shù)應(yīng)同時(shí)具有較好的降噪效果和較強(qiáng)的實(shí)時(shí)性．為對(duì)降噪效果進(jìn)行評(píng)價(jià)，通常將信噪比作為評(píng)價(jià)降噪效果的主要指標(biāo)，但信噪比的計(jì)算需以不含噪信號(hào)作為參考，而在缸蓋振動(dòng)測(cè)量過(guò)程中，不含噪信號(hào)是不可測(cè)的．因而可通過(guò)組合均方根誤差和平滑度給出一種小波降噪復(fù)合評(píng)價(jià)指標(biāo)[20]，均方根誤差與平滑度呈負(fù)相關(guān)，在組合的過(guò)程中會(huì)出現(xiàn)一個(gè)極值，該極值表示去噪后信號(hào)保留的細(xì)節(jié)信息和逼近信息達(dá)到了最 佳比例，其對(duì)應(yīng)的為尋找的最優(yōu)小波基函數(shù)為

式中：PRMSE為歸一化后的均方根誤差；PR為歸一化后的平滑度；WPRMSE和為通過(guò)變異系數(shù)定權(quán)法所確定的權(quán)值．T最小時(shí)表示降噪效果最優(yōu)．變異系數(shù)定權(quán)法是指在評(píng)價(jià)指標(biāo)體系中數(shù)值差異越大的指標(biāo)，即越難實(shí)現(xiàn)的指標(biāo)，因而要賦以更大的權(quán)值．權(quán)值的計(jì)算方法為

式中：σ和μ分別為求標(biāo)準(zhǔn)差操作與均值操作．

式(18)的評(píng)價(jià)指標(biāo)僅對(duì)降噪效果進(jìn)行評(píng)價(jià)，為綜合評(píng)價(jià)小波基函數(shù)的降噪時(shí)間和降噪效果，利用變異系數(shù)定權(quán)法將降噪效果和降噪時(shí)間線性組合，進(jìn)一步提出綜合評(píng)價(jià)指標(biāo)S，有

式中：PT和Pt為多工況下歸一化的平均降噪效果和平均降噪時(shí)間；和為通過(guò)變異系數(shù)定權(quán)法確定的權(quán)值．由于小波基函數(shù)的降噪時(shí)間一般與降噪效果呈負(fù)相關(guān)，且兩者均為最小時(shí)最優(yōu)，故S存在極小值，該極值表示降噪效果和降噪時(shí)間達(dá)到了最佳比例，此時(shí)該小波基函數(shù)為綜合最優(yōu)．

若函數(shù)僅在有限的定義域范圍內(nèi)存在非零函數(shù)值，在該區(qū)間外的函數(shù)值全為0，稱(chēng)該函數(shù)在這個(gè)區(qū)間上緊支，具有該性質(zhì)的小波稱(chēng)為緊支撐小波，只有緊支撐小波才可用于離散小波變換．正交小波具有去除各子帶數(shù)據(jù)相關(guān)性的特性，且計(jì)算量相對(duì)非正交小波要小，更適合用于實(shí)時(shí)降噪．筆者篩選出工程應(yīng)用中針對(duì)振動(dòng)信號(hào)應(yīng)用較為廣泛的緊支撐正交小波基函數(shù)db1～db25、sym1～sym25、coif1～coif5與緊支撐雙正交小波基函數(shù)bior1.1～bior6.8，圖4為各小波基函數(shù)在多種工況下的平均降噪時(shí)間及平均降噪效果．

圖4 不同小波基函數(shù)降噪時(shí)間與降噪效果Fig.4 Noise reduction time and noise reduction effect of different wavelet basis

圖5為進(jìn)一步計(jì)算得到小波基函數(shù)的綜合評(píng)價(jià)指標(biāo)S．小波基函數(shù)序號(hào)從左至右分別為db1～ db25、sym1～sym25、bior1.1～bior6.8和coif1～coif5．可知，S存在4處較小點(diǎn)，其中最小值出現(xiàn)在2號(hào)處，坐標(biāo)為(53，0.029061)，其對(duì)應(yīng)的小波基函數(shù)為bior2.2，故bior2.2為綜合最優(yōu)小波基函數(shù)．

圖5 不同小波基函數(shù)綜合評(píng)價(jià)指標(biāo)Fig.5 Comprehensive evaluation index of different wavelet basis

4 基于DSP的監(jiān)測(cè)系統(tǒng)實(shí)時(shí)性試驗(yàn)

試驗(yàn)選用可提供150MIPS指令速度的32位數(shù)字信號(hào)處理器TMS320C28x作為信號(hào)處理的核心芯片．為了實(shí)時(shí)評(píng)價(jià)所選參數(shù)，定義實(shí)時(shí)性特征δ為

式中：ts為采樣時(shí)間；ta為分析時(shí)間；δ越小，實(shí)時(shí)性越強(qiáng)，δ≤1時(shí)可實(shí)現(xiàn)實(shí)時(shí)處理．

試驗(yàn)設(shè)置采樣頻率為25kHz，選用3層小波分解和bior2.2小波基函數(shù)．通過(guò)算法相應(yīng)的指令周期模型計(jì)算其實(shí)時(shí)性特征，并在TMS320C28x上通過(guò)試驗(yàn)獲取了不同算法在不同采樣點(diǎn)數(shù)下的實(shí)時(shí)性特征，如圖6所示．

圖6 不同小波實(shí)現(xiàn)算法的實(shí)時(shí)性對(duì)比Fig.6 Real-time comparison of different wavelet algorithms

可知，試驗(yàn)結(jié)果與建立的小波變換指令周期模型吻合較好．但試驗(yàn)所得的實(shí)時(shí)性特征值略高于指令周期模型預(yù)測(cè)值，這是由于在函數(shù)運(yùn)行過(guò)程中存在臨時(shí)變量建立、循環(huán)變量累加等原因，導(dǎo)致實(shí)時(shí)性變差．

不同的小波變換序列點(diǎn)數(shù)對(duì)實(shí)時(shí)性有較大影響．為獲得較好的降噪效果，需合理確定小波變換序列點(diǎn)數(shù)，小波變換序列點(diǎn)數(shù)的統(tǒng)計(jì)特征應(yīng)與被測(cè)信號(hào)整體的統(tǒng)計(jì)特征基本一致．柴油機(jī)缸蓋振動(dòng)信號(hào)隨曲軸旋轉(zhuǎn)呈周期性變化，一個(gè)完整旋轉(zhuǎn)周期能較好地代表柴油機(jī)短時(shí)間內(nèi)的運(yùn)行狀態(tài)．取怠速時(shí)曲軸旋轉(zhuǎn)一周時(shí)間(0.075s)，則小波變換序列點(diǎn)數(shù)應(yīng)大于0.075×25000＝1875．此時(shí)使用非原位運(yùn)算提升小波的實(shí)時(shí)性最優(yōu)且滿(mǎn)足δ≤1，可以滿(mǎn)足在25kHz采樣頻率下的實(shí)時(shí)降噪需求．

5 結(jié) 論

(1) 原位運(yùn)算提升小波運(yùn)行所需的時(shí)間隨變換序列增加呈二次增長(zhǎng)關(guān)系，非原位運(yùn)算提升小波和Mallat算法運(yùn)行所需時(shí)間隨變換序列增加呈線性增長(zhǎng)關(guān)系；當(dāng)序列點(diǎn)數(shù)超過(guò)8時(shí)，非原位運(yùn)算提升小波的實(shí)時(shí)性將優(yōu)于原位運(yùn)算提升小波．

(2) 指令周期模型能較好地?cái)M合DSP中小波算法的實(shí)時(shí)性，但根據(jù)模型計(jì)算所得的實(shí)時(shí)性會(huì)稍強(qiáng)于實(shí)際運(yùn)行時(shí)的實(shí)時(shí)性．

(3) 提供150MIPS指令速度時(shí)，選用3層小波分解和bior2.2小波基函數(shù)，小波實(shí)現(xiàn)算法選用非原位運(yùn)算提升小波快于Mallat算法和原位運(yùn)算提升小波，且使用非原位運(yùn)算提升小波可滿(mǎn)足25kHz采樣頻率下的實(shí)時(shí)降噪需求．