基于LMD 降噪和互相關的聲波飛行時間測量法

明靖川，顏 華，魏元焜

（沈陽工業(yè)大學 信息科學與工程學院，沈陽 110870）

基于聲波飛行時間測量的溫度場重建技術具有非接觸、維護方便等優(yōu)點，在工業(yè)爐等高溫領域有突出優(yōu)勢。該項技術還可用于海底熱液口、倉儲糧食、集成電路晶圓溫度監(jiān)測等領域[1-4]。飛行時間測量屬時延估計，常采用互相關法。為提高估值的抗噪性能，文獻[5]將三次相關與廣義互相關結(jié)合；文獻[6]將變分模態(tài)分解與廣義二次互相關結(jié)合；文獻[7]將譜減法與互功率譜頻域加權法結(jié)合。但這些研究中未考慮相關噪聲的影響，算法驗證時添加的是互不相關的白噪聲?，F(xiàn)實世界持續(xù)性或間歇性的噪音無處不在。同一噪音源，如汽車鳴笛聲，在各聲波接收器中導致的觀測噪聲彼此相關，而互相關法估值對相關噪聲的抵抗能力很弱。

基于局部均值分解（LMD）[8]具有自適應性好、端點效應小、信號分解后噪聲集中于高頻分量的特點，本文提出一種將信號的LMD 重構(gòu)和互相關（CC）結(jié)合的聲波飛行時間測量法，簡稱LMD-CC 法。首先依據(jù)環(huán)境噪音決定是否短截信號以降低模態(tài)混疊的概率；而后對信號進行LMD 分解，基于本文提出的有用分量篩選法重構(gòu)降噪信號；最后通過對降噪信號的互相關運算獲得聲波飛行時間估值，并采用模擬糧倉中的實測數(shù)據(jù)，驗證了所提方法的有效性。

1 LMD-CC 算法原理

1.1 互相關時延估計

互相關法是時延估計最基本、最常用的方法。假設揚聲器發(fā)出的有用聲波信號s（n）先后被近端傳聲器1 和遠端傳聲器2 接收，如混響可以忽略，則傳聲器1、2 接收的信號可表示為[3]

式中：G 為聲波傳輸過程中的衰減因子；D 為兩傳聲器接收信號間的時間延遲；n1（n）、n2（n）分別為傳聲器1、2 接收信號中的噪聲。

n1（n）、n2（n）的互相關函數(shù)可以表示為

式中：Rss（·）為s（n）的自相關函數(shù)，當m=D 時Rss（m-D）達到最大值；Rn1s（·）、Rsn2（·）、Rn1n2（·）分別為n1（n）和s（n）、s（n）與n2（n）、n1（n）和n2（n）的互相關函數(shù)，當s（n）、n1（n）、n2（n）彼此互不相關時，這些互相關函數(shù)的值皆為0，此時有：

互相關法通過Rx1x2（m）峰值點的橫坐標，確定2個接收信號間的時間延遲D[3]。實際測量中有各種噪音干擾，同一噪音源在各傳聲器接收信號中對應的噪聲彼此相關，因此有Rn1n2（m）≠0，導致互相關法性能下降甚至失效。

1.2 局部均值分解法

LMD 法[8-9]可以根據(jù)信號自身特點自適應地將一個由多分量構(gòu)成的復雜非平穩(wěn)信號分解為一系列瞬時頻率由高到低且具有物理意義的PF（product function）分量。每一個PF 分量可由一個包絡信號與一個純調(diào)頻信號之積得出。包絡信號就是該PF 分量的瞬時幅值，而PF 分量的瞬時頻率則可以由純調(diào)頻信號直接求出。存在某一PF 分量，其前、其后的分量分別以噪聲和有用信號為主要成分，因此可通過剔除噪聲為主的PF 分量后的信號重構(gòu)，實現(xiàn)降噪。

假設原始信號為x（t），令c（t）=x（t），p=q=1，角標p、q 代表第p 層PF 分量的第q 次迭代，本文的LMD 分解降噪步驟可簡述如下[9]：

步驟1初始化包絡函數(shù)fp（t）=1。

步驟2尋找c（t）的局部極大、極小值點，用三次樣條插值法獲得極大值和極小值所構(gòu)向量的插值函數(shù)mpq（t）和npq（t），進而獲得局部均值函數(shù)lpq（t）、局部包絡函數(shù)fpq（t）、包絡函數(shù)fp（t）。與常用的滑動平均法相比，三次樣條插值法獲得的局部均值函數(shù)和局部包絡函數(shù)對應著更好的LMD 分析精度[10]。

步驟3計算spq（t），若spq（t）為局部包絡函數(shù)為1 的純調(diào)頻信號，跳到步驟4；若不是則令c（t）=spq（t），q=q+1 回到步驟2。

步驟4計算第p 層的PF 分量PFp（t）和對應的殘余分量rp（t），若rp（t）為單調(diào)函數(shù)或分解層數(shù)p達到預設閾值（本文閾值為15），則停止迭代，結(jié)束LMD 分解過程，并令n=p；若不滿足迭代停止條件則令c（t）=rp（t），p=p+1 回到步驟1。

LMD 分解過程結(jié)束后，去掉前k-1 個PF 分量的重構(gòu)信號，如式（10）所示，顯然k=1 時x（t）。

1.3 有用PF分量的篩選

LMD 降噪的關鍵是尋找有用信號分量與噪聲分量的分界點，確定參與信號重構(gòu)的有用PF 分量。文獻[11]提出一種基于CMSE 的信噪分界點確定法。若CSMEk隨k 變化的曲線在全局極小值之前有局部極小值，則分界點取第1 個局部極小值對應位置加1；若沒有則取全局極小值對應位置加1。但實踐發(fā)現(xiàn)原始信號是在安靜環(huán)境下采集的高信噪比信號時，該方法確定的分界點會將較多的有用分量當做噪聲剔除，導致時延估值誤差很大。

本文提出一種以PF 分量的累加能量（累加CMSE）判據(jù)為主，以PF 分量與原始信號相關性為輔的有用PF 分量篩選法，具體步驟簡述如下：

步驟1計算前j 個PF 分量的累加能量并記為ACMSEj：

步驟2前j 個PF 分量的累加能量達到所有PF 分量累加能量的二分之一時，令dp=j。

步驟3先以dp 為信噪分界點篩除高頻噪聲分量，而后計算各剩余PF 分量與原始信號x（t）的相關系數(shù)，剔除相關系數(shù)小于剩余PF 分量中最大相關系數(shù)的十分之一的低相似性分量，最后將剩余的PF 分量視為有用PF 分量，與殘余分量一起重構(gòu)出降噪后的信號。

1.4 減少模態(tài)混疊的影響

汽車鳴笛聲、飛機引擎聲以及嘈雜人聲等都以間歇形式出現(xiàn)，描述間歇事件的模態(tài)也須具備間歇特性，但LMD 算法不具備間歇特性，這也是間歇現(xiàn)象引起模態(tài)混疊的內(nèi)在原因[12]。模態(tài)混疊會使PF 模態(tài)分量中出現(xiàn)虛假的模態(tài)分量，甚至會導致其失去物理意義。因此環(huán)境噪音較大時只截取采集信號中最具特征的部分參與LMD 分解，可降低模態(tài)混疊出現(xiàn)的概率。但環(huán)境噪音較小、模態(tài)混疊概率很低時，使用完整的采集信號可獲得更穩(wěn)定的時延估值?；谀M糧倉中實測數(shù)據(jù)的時延估計，也驗證了這一推斷。傳聲器接收到信號被人為地分為有用信號和噪聲信號，但在能量角度二者并無差異。為此本文在近端傳聲器接收的信號能量大于安靜環(huán)境下的α倍時，只截取采集信號中最具特征的部分參與LMD分解。由于聲學法測溫中傳聲器位置固定，所以信號中最具特征的部分可事先確定。

1.5 結(jié)合LMD 降噪和互相關的的測量法（LMDCC 法）

基于以上研究，本文提出一種結(jié)合LMD 降噪和基本互相關的聲波飛行時間測量法，簡稱LMD-CC法。首先，計算近端傳聲器所接收信號的能量，若大于安靜環(huán)境下近端信號能量J0的α 倍，則判定環(huán)境噪音比較大，分別對近、遠端采集到的信號做短截處理；而后，分別對近、遠端信號（若滿足短截條件則是對短截的近、遠端信號）進行LMD 分解、篩選有用PF 分量進行信號重構(gòu)；最后對降噪后的近、遠端信號進行互相關運算，得到聲波飛行時間。圖1 所示為LMD-CC 算法流程。

圖1 LMD-CC 法流程Fig.1 LMD-CC method flow chart

2 LMD-CC 法的實驗驗證

2.1 實驗數(shù)據(jù)的獲取

圖2 給出了揚聲器1、2 和近端傳聲器1、遠端傳聲器2 布局示意圖。揚聲器1 用來產(chǎn)生有用的聲波信號，和傳聲器1、2 位于同一水平線上且皆置于大豆中。揚聲器2 用來產(chǎn)生噪音信號，置于糧槽外，糧槽中裝滿大豆。傳聲器1、2 的輸出經(jīng)信號調(diào)理器后由PCI-6123（采樣率設置為500 KS/s）同步采集到計算機。

圖2 揚聲器、傳聲器布局示意圖Fig.2 Layout diagram of speakers and microphones

實驗數(shù)據(jù)的采集均在安靜環(huán)境下進行，每個信號的采集長度均為30000 點。首先，令揚聲器2 靜默，揚聲器1 發(fā)出時長為0.3 s 的chirp 信號，通過傳聲器1、2 采集一組近、遠端有用聲波信號s1（n）和s2（n），重復該過程獲得100 組有用聲波信號。而后，令揚聲器1 靜默，揚聲器2 發(fā)出汽車鳴笛聲、飛機引擎聲以及人聲的混合噪音，通過傳聲器1、2 采集一組近、遠端噪音信號n1（n）和n2（n），重復該過程獲得100 組噪聲信號。最后，由相同組號的安靜環(huán)境下采集的有用聲波信號和噪音信號合成出噪音環(huán)境下傳聲器i（i=1，2）采集到聲波信號，即xi（n）=si（n）+gni（n）。該實驗數(shù)據(jù)獲取方法的優(yōu)點是方便控制帶噪信號的信噪比，減少了實驗中不可控的環(huán)境變量干擾，而且安靜環(huán)境下采集到的有用聲波信號間的時間延遲也可為各降噪算法在不同信噪比下的有效性提供比較基準。

以其中一組信號為例，圖3（a）～圖3（c）給出了安靜環(huán)境下的近、遠端聲波信號以及它們的互相關圖；圖3（d）～圖3（f）給出了噪音環(huán)境下的近、遠端聲波信號以及它們的互相關圖，圖3（d）、圖3（e）中標注了信號滿足短截條件時被短截的區(qū)域，遠端信號的信噪比為-5.076 dB；圖3（g）～圖3（i）給出了LMD降噪后的近、遠端聲波信號以及它們的互相關圖，由于滿足信號短截條件，所以圖3（g）、圖3（h）中聲波信號被短截?；ハ嚓P圖中給出了峰值點坐標，峰值點的橫坐標即為近、遠端信號之間的時延估值。由圖3 可以看出：①LMD 降噪較好地還原了原信號的特征；②直接對噪音環(huán)境下采集的信號做互相關得到的時延估值與安靜環(huán)境下的時延估值有很大偏差，為23.01 ms，但LMD 降噪后的時延估值與安靜環(huán)境下的時延估值一致，皆為5.79 ms。

圖3 聲波信號和互相關波形圖Fig.3 Acoustic signals and cross correlation waveforms

2.2 算法對比

為對本文算法的有效性做更充分的檢驗，采用本文提出的LMD-CC 法，以及CC（互相關）法、EMD（經(jīng)驗模態(tài)分解）降噪與互相關結(jié)合的EMD-CC 法、VMD（變分模態(tài)分解）降噪與互相關結(jié)合的VMDCC 法，計算了安靜環(huán)境和噪音環(huán)境下的聲波飛行時間估值。EMD 降噪時采用了文獻[10]所述CMSE準則進行重構(gòu)層數(shù)的選取，VMD 降噪時則采用了利用相關系數(shù)確認分解層數(shù)、用粒子群算法獲得懲罰因子的自適應VMD。

表1 給出了聲波飛行時間估值的平均值和標準差，均值和標準差對應的樣本數(shù)皆是100。估值的標準差越小，估值的穩(wěn)定性越好；估值的平均值與安靜環(huán)境下的平均值越接近，估值的準確性越好。由表1 可以看出，遠端信噪比為10 dB、5 dB 即環(huán)境噪音較小時，4 種算法均具有較好的穩(wěn)定性和準確性；當遠端信噪比為0.08 dB、-5.08 dB 即環(huán)境噪音較大時，LMD-CC 的結(jié)果顯著優(yōu)于其它3 種算法。

表1 相關噪聲干擾下不同方法對應的聲波飛行時間估值評價Tab.1 Evaluation of acoustic flight time of different methods under related noise interferenc

各方法運算1 組數(shù)據(jù)的平均運行時間如下：CC法為0.024 s；EMD-CC 法為1.92 s；VMD-CC 為78.86 s；LMD-CC 法在高信噪比環(huán)境下為0.28 s，低信噪比環(huán)境為0.10 s。

3 結(jié)語

為提高聲學法溫度場重建系統(tǒng)對環(huán)境噪音的抵抗能力，本文提出了一種LMD 降噪和互相關結(jié)合的聲波飛行時間測量法，并通過模擬糧倉中獲取的實驗數(shù)據(jù)，驗證了該方法對環(huán)境噪音導致的相關噪聲有良好的抵抗能力，且具有較快的運行速度。

本文所提方法可提高聲學法溫度場重建系統(tǒng)的環(huán)境適應能力，因此有較大的實用價值。