基于人耳聽覺特性的信號顯著圖計算方法研究

李國萌， 李允公， 王 波， 吳文壽， 安 超

(東北大學 機械工程及自動化學院，沈陽 110819)

基于人耳聽覺特性的信號顯著圖計算方法研究

李國萌， 李允公， 王 波， 吳文壽， 安 超

(東北大學 機械工程及自動化學院，沈陽 110819)

信號的幅值突變往往蘊含著設備的故障信息，為明確產生幅值突變的時頻區域，基于聽覺系統的信息處理機制，提出了一種信號顯著圖計算方法。首先,對信號進行帶通濾波、相位調整、半波整流等聽覺外周處理，然后，提取處理結果的一次或多次包絡信息，并對包絡信息進行多尺度二維濾波。繼而,利用中心-周邊差算子得到信號在不同尺度下的局部顯著度，最后，對局部顯著度進行跨尺度整合和線性合并得到信號的時頻顯著圖和全局顯著圖。仿真和實驗結果表明，該方法適用于振動信號和語音信號,可以清晰表征幅值突變的時頻區域，并可抵抗一定的噪聲干擾，具有一定的有效性和實用性。

聽覺模型；顯著圖；故障診斷；二維濾波；包絡分析

在機械故障診斷和工程信號分析領域中，突變信號的檢測與識別一直是一個重要的研究課題，信號的突變點和不規則變化部分包含許多重要信息，有助于被測對象的狀態識別。例如，齒輪斷齒、軸承松動、動靜件碰摩等，都會產生沖擊響應振動信號成分，而語音信號的每一次發聲也相當于信號的一次突變。對于成分復雜的信號，突變往往只產生于某些局部時頻區域，甚至在原始波形中難以觀察到明顯的幅值突變。因此，有必要分析各時頻區域中的信號分量是否存在幅值的顯著性變化，所得結果即為語音分析領域中常用的顯著圖。顯著圖的實質是對信號帶通濾波結果的二次處理，其主要特點是時頻局部性與尺度多樣性，相對于傳統的歪度或峭度等統計特征而言，顯著圖能夠準確描述產生突變的時頻位置，相對于小波變換等分解方法來說，顯著圖可以明確說明某時頻區域中是否發生了幅值突變，從而便于計算機自動識別。

目前，在顯著圖計算方面，夏虹等[1]提出了基于小波變換的突變信號識別方法，通過單子帶算法和過渡函數對信號的突變點進行定位；張欣等[2]提出了一種聽覺顯著圖計算模型，首先對聲音信號劃分頻帶，再在各頻帶上采用一維高斯濾波器進行時域濾波提取顯著度；王雪君等[3]提出了一種有效的自底而上的聽覺注意顯著圖計算模型。夏虹的方法雖然能找到信號的突變點，但是不能顯示信號顯著變化的整個區域，計算過程中容易丟失突變信號的某些特征；而張欣和王雪君提出的聽覺顯著圖計算模型，只適用于聲學信號。

聽覺系統具有識別、分析、處理信號的功能，并且對于突變信號具有自底而上的外源性注意[4]機制，系統各部件可自動調節自身參數以適應信號的突變，從而使聽覺皮層感知到突變信號所在的時頻區域。同時，聽覺系統的數學模型已在語音、聲吶、噪聲信號分析中獲得成功應用，且已被證明也適用于機械振動信號分析[5]。因此，本文將顯著圖概念引入到振動信號分析當中，并結合振動信號的特點，提出了一種基于聽覺信息處理機制的顯著圖計算方法，并用不同種類的實際信號進行了仿真和實驗驗證。

1 方法的基本原理

本文所提方法流程圖如圖1所示。首先，對信號作聽覺外周處理[6]，用Gammatone濾波器組[7]模擬耳蝸基底膜對輸入信號作帶通濾波、逆序濾波[8]調整相位、半波整流[9]等計算，并對經聽覺外周處理的信號作一次或多次包絡分析[10]。并用二維高斯濾波器[11]模擬初級聽覺皮層應答野[12]的感知特征，對包絡曲面做多尺度二維濾波，得到信號在不同尺度下的時頻域表示。繼而，用中心-周邊差算子[13]計算不同尺度下的局部顯著度，最后對局部顯著度進行跨尺度整合和線性合并得到信號的時頻顯著圖和全局顯著圖。

圖1 計算方法流程圖Fig.1 Calculating method flow sheet

2 方法的具體實現

2.1 帶通濾波

聽覺模型中常用Gammatone濾波器組對信號進行帶通濾波，以模擬耳蝸基底膜的頻率分解功能。設濾波器數目為N，被分析的信號為x(t),則第n個濾波器的表達式為

h(n,t)=Bαtα-1e-2πBtcos(2πfn+φn)u(t)

(1)

式中:fn為第n個濾波器的中心頻率；α為濾波器的階數，一般取為4；u(t)為階躍函數 ；φn為第n個濾波器的初相位，通常取為0；參數B的計算公式如下

B=1.109ERB(fn)=1.109(24.7+0.108fn)

(2)

式中:ERB(fn)為第n個濾波器的矩形寬度。需要說明的是Gammatone濾波器組的各個中心頻率在頻率軸上呈對數均勻分布，濾波器的等階矩形帶寬也呈對數分布。

輸入信號x(t)，通過第n個濾波器后時域輸出為y(n,t)，則有：

y(n,t)=x(t)*h(n,t)

(3)

式中:*代表卷積。

2.2 相位調整

由于同一信號在不同濾波通道中濾波會出現相位差異，需對濾波結果進行相位調整，聽覺模型中常采用逆序濾波的方法，即將y(n,t)逆序后重新輸入濾波器，再將濾波輸出反序。由于逆序濾波過程繁瑣、計算實質和互相關函數相同，因此，本文借助互相關函數實現逆序濾波，即

(4)

2.3 半波整流

在聽覺系統中，耳蝸內毛細胞具有半波整流的特性，即只對振動刺激信號波形的正半部分產生響應，實現方法如下

(5)

半波整流不僅簡化了計算過程，也為后續的分析提取提供了方便。

2.4 包絡分析

由于聽覺中樞只對傳入信號的局部極值產生響應，因此，對半波整流后的信號做包絡分析。

設y1(n,t)的離散形式為y1(n,m)(m=1,2,…,M)，尋找各個通道內信號的局部極大值點max(tp,fp)，其中tp和fp分別表示該通道內第p個極大值點的時間和幅值。利用三次樣條插值方法擬合波形，得到一次包絡e1(n,m)，可以顯示信號幅值波動的局部細微變化。

為了突出信號的整體變化趨勢，對信號進行二次包絡，尋找每個通道內一次包絡e1(n,m)的局部極大值點，利用三次樣條插值方法擬合波形，得到較為平滑的二次包絡e2(n,m)。對于波形復雜、干擾噪聲較多的信號，甚至可取到三次或四次包絡。

2.5 多尺度二維濾波

聽覺皮層中不同的神經元會對不同的波形產生最大響應，因此，可對包絡線ei(n,m)(i=1,2,3,4;)進行多尺度濾波計算，從而模擬聽覺皮層的信息處理機制。本文使用二維高斯濾波器實現這一功能。二維零均值高斯濾波器gj(n,m,u,v)(j=1,2…6)的函數表達式如下

gj(n,m,u,v)=

(6)

式中:σj和εj決定了二維高斯函數的尺度,σ1=3,ε1=7,σj=j×σ1,εj=j×ε1。

圖2 二維高斯濾波器Fig.2 Two dimensional Gaussian filter

(7)

表達方便，令

(8)

多尺度二維濾波的主要目的是全面反映空間曲面ei(n,m)中投影面積不同的局部凸起的分布情況。

2.6 局部顯著度

(9)

式中：中心尺度c∈(1,2,3)，周邊尺度s=c+φ，φ∈(2,3)，c

2.7 全局顯著圖

(10)

最后，對Qi(n,m)再進行線性合并，得到信號的全局顯著圖，即

(11)

3 仿真和實驗

為了驗證本文所提方法的有效性和實用性，采用不同種類的信號分別進行仿真和實驗驗證。

3.1 仿真驗證

設仿真信號為

x(t)=sin(2π×2 000t)+α·s(t)+no(t)

(12)

式中:no(t)干擾噪聲，α為沖擊成分增益系數，s(t)為沖擊成分，表達式為

s(t)=exp(-100t)[3sin(2π×1 000t+0.3)+

2sin(2π×2 000t+0.5)+2sin(2π×

(13)

式中:δ(t)為脈沖函數?？梢姡瑂(t)中存在若干個以0.11 s為間隔的間歇性瞬態信號，且信號中包括三個主要頻率成分，這兩個特點均與實際的沖擊響應相接近。仿真中，Gammatone濾波器個數N為200，仿真信號的采樣頻率為10 000 Hz，包絡分析環節采用一次包絡，信號長度為5 000個點。當α=1，no(t)=0時，x(t)的方差為1.528 8，其時域波形如圖3所示。

圖3 時域波形Fig.3 Time domain waveform

將一次包絡e1(n,m)通過6個不同尺度的二維高斯濾波器，對多尺度二維濾波結果進行跨尺度相減和整合后得到時頻顯著圖Q1(n,m)，如圖4所示。

圖4 時頻顯著圖Q1(n,m)Fig.4 Time-frequency saliency map Q1(n,m)

在Q1(n,m)中，可明顯看到仿真信號中存在4個沖擊成分，對Q1(n,m)線性合并，得到了仿真信號的全局顯著圖T1(m)，如圖5所示。

圖5 全局顯著圖T1(m)Fig.5 Global saliency map T1(m)

由圖5知，顯著圖中4個沖擊成分極大值對應點分別為1 121、2 219、3 320、4 422，和時域波形中沖擊極大值對應點數不完全一致，這與帶通濾波、包絡分析、多尺度二維濾波等計算環節產生延時有關，但是點數誤差在可以接受的范圍內，所以，后面參數的比較可以以顯著圖中沖擊極大值對應的點數為依據，而且，在沖擊響應的起止時刻，顯著度發生明顯的突變，顯著圖反映出沖擊響應的整體顯著變化區域，說明本文所提方法可以準確判定沖擊成份的位置和起止時間。

當α=1，no(t)為方差為6.124 7的高斯白噪聲時，仿真信號時域波形和全局顯著圖如圖6所示。從圖6(a)可知，加入強烈的高斯白噪聲后，時域波形中無法觀察到4個明顯的沖擊成分，但圖6(b)準確的表達出信號中存在4個顯著變化時段。和圖5比較，雖然兩個顯著圖在波形上存在差別，但整體上沖擊響應的顯著度突變情況幾乎一致。

(a) 時域波形 (b) T1(m) 圖6 時域波形及其全局顯著圖T1(m)Fig.6 Time domain waveform and global saliency map

(14)

式中：Δmq表示的是第q個沖擊極大值偏離點數，最后統計結果如表1和圖7所示。

表1 噪聲對全局顯著圖T1(m)的影響

圖7 噪聲方差不同時全局顯著圖峭度值Fig.7 Kurtosises of global saliency map corresponding to different noise variances

從以上數據可知，加入的高斯白噪聲方差越大，顯著圖峭度值越小，但噪聲在一定范圍內時，顯著圖依然具有良好的識別效果，而且4個沖擊成分極大值對應的點數和無噪聲時比較，偏離的誤差都在可以接受的范圍內，說明本文所提方法具有一定的抗噪聲干擾能力。而當仿真信號中加入方差為20.560 9的高斯白噪聲時，仿真信號時域波形和全局顯著圖如圖6所示，此時沖擊成分完全淹沒在噪聲中,全局顯著圖已經無法準確識別出沖擊成分。另外，之所以只用一次包絡，是因為仿真信號成分比較簡單，得到的顯著圖已經能夠準確識別出信號里面的沖擊成分。

(a) 時域波形 (b) T1(m)圖8 時域波形及其全局顯著圖T1(m)Fig.8 Time domain waveform and global saliency map T1(m)

3.2 實驗驗證

3.2.1 語音信號驗證

為了驗證本文所提方法對語音信號具有良好的適用性，取長為100 000個采樣點的中文詞組“振動”進行實驗驗證。語音信號采樣頻率為44 100 Hz，濾波器個數N為100。語音信號的時域波形如圖9所示。

圖9 時域波形Fig.9 Time domain waveform

將二次包絡e2(n,m),通過6個尺度的二維高斯濾波器,對濾波結果進行跨尺度相減和整合，得到時頻顯著圖Q2(n,m)，如圖10所示。

圖10 時頻顯著圖Q2(n,m)Fig.10 Time-frequency saliency map Q2(n,m)

圖11(a)和圖11(b)分別是在一次包絡、二次包絡分析的基礎上得到的語音信號的全局顯著圖，均能很好地識別出語音信號的顯著變化區域，全局顯著圖的波峰和語音中音節的位置一一對應，且語音信號音節開始和停頓時刻，顯著度突變。圖11(b)和圖11(a)相比較，二次包絡分析基礎上得到的顯著變化區域輪廓更為清晰簡潔，整體顯著性更好。

(a) T1(m) (b) T2(m)圖11 全局顯著圖Fig.11 Global saliency map

為觀察語音信號和全局顯著圖的細微情況，對語音信號和全局顯著圖的第25 000到第40 000個點進行局部放大，時域波形局部放大結果如圖12所示。

圖12 時域波形的局部放大Fig.12 Local enlarged time domain waveform

基于一次包絡、二次包絡得到的全局顯著圖局部放大結果如圖13所示。從圖中可以明顯看出，顯著度的變化趨勢和語音信號幅值變化趨勢相同，且變化位置一一對應。圖13(a)能較好的反映出語音信號的細微變化，圖13(b)相當于是對一次包絡做平滑處理，更能突出語音信號的整體顯著性。

(a) T1(m) (b) T2(m)圖13 全局顯著圖的局部放大Fig.13 Global enlarged global saliency map

實驗結果表明，本文所提方法對語音信號具有良好的適用性，可以對語音信號的突變點進行準確定位，利用全局顯著圖可以準確檢測和判定語音信號的起止時刻和每個音節的具體位置。

3.2.2 振動信號驗證

為了驗證本文所提方法對實際振動信號同樣具有良好的有效性和實用價值，搭建如圖14所示的二級斜齒圓柱齒輪減速器實驗臺，并在中間軸齒輪上人為制造了斷齒故障。其主要零部件參數如下:電動機額定功率750 W，輸入軸轉速937 r/min，減速器高低速級齒數分別為13，86，14，85。理論計算，中間軸的轉速為141.636 r/min, 即斷齒齒輪每轉一圈耗時為0.423 s。

(a) 試驗臺 (b) 斷齒故障圖14 減速器試驗臺和人為斷齒故障Fig.14 Retarder testing machine and human-made broken fault

用壓電式加速度傳感器采集振動信號，采樣頻率設為10 240 Hz，信號長度為15 000個點，斷齒故障加速度振動信號的時域波形如圖15所示。依據斷齒齒輪轉速，理論計算可知, 15 000個采樣點中應該包含4個與斷齒有關的沖擊響應信號。

圖15 時域波形Fig.15 Time domain waveform

Gammatone濾波器個數N為200，濾波結果如圖16所示，無法清楚的識別出沖擊成分，需要對信號作進一步處理。

圖16 Gammatone濾波結果Fig.16 Gammatone filtering result

圖17 多尺度二維濾波結果(n,m)Fig.17 Multi-scaled two dimensional filtering results

圖18 時頻顯著圖Q1(n,m)Fig.18 Result of multi-scaled integration Q1(n,m)

對Q1(n,m)線性合并得到的結果如圖19所示，全局顯著圖中可以看到4個明顯的沖擊成分，相互間的平均間隔為4 327個點，即0.42 s，與中間軸旋轉一周的時間相同，從而說明產生的4個周期性的沖擊成分是由斷齒故障造成的，說明本文所提方法同樣適用于實際振動沖擊信號，可以準確識別出沖擊響應的顯著變化區域。

圖19 全局顯著圖T1(m)Fig.19 Global saliency map

為了觀察斷齒故障振動信號的整體顯著性變化，在包絡分析環節采用二次包絡和三次包絡重新計算全局顯著圖，結果如圖20所示。對比圖19和圖20可知，一次包絡和二次包絡基礎上得到的顯著圖結果差別微小，但三次包絡得到的結果清晰整潔，4個周期性出現的沖擊響應波峰突出，顯著變化區域明顯。

(a) T2(m) (b) T3(m)圖20 全局顯著圖Fig.20 Global saliency map

為了觀察斷齒故障振動信號和計算結果的細節情況，截取振動信號和不同包絡次數下得到的全局顯著圖的第3 500到第4 100個點進行局部放大，結果如圖21所示。從圖18可知，該段信號幅值變化微小，但是圖21(b)和圖21(c)顯著度變化依然較大，顯著變化區域明顯。由此可知，本文所提方法對微弱的沖擊和幅值變化較小的信號依然具有較好的識別效果，具有一定的靈敏特性。

(a) 時域波形 (b) T1(m)

(c) T2(m) (d) T3(m)圖21 時域波形及全局顯著圖的局部放大Fig.21 Local enlarged time domain waveform and global saliency map

之所以出現這樣的結果，是因為實際采集到的振動信號成分復雜，除了因斷齒故障出現的明顯的周期性沖擊成分，還包含很多微弱沖擊(如齒側拍擊)。因此，基于一次包絡得到的顯著圖除了反映能量相對較大的由斷齒誘發的沖擊響應外，還會表達出其他的微弱的沖擊。一次包絡顯著性識別度較高，而二次和三次包絡均相當于是對前一次包絡的平滑處理，從而可使沖擊響應的顯著變化區域更為突出，并略去幅值較小的沖擊成分和干擾噪聲的影響。這也再次說明，多次包絡分析的結果整體顯著性較好，更貼近于人耳對突變信號的實際感知情況，對于不同的信息提取要求，包絡分析環節可以選擇不同包絡次數。

4 結 論

(1)基于人耳聽覺特性,提出了一種信號顯著圖計算方法，該方法既適用于語音信號，也適用于實際振動信號，可以準確識別出信號的顯著變化區域，具有一定的抗噪聲干擾能力。

(2)包絡分析環節可以選擇一次包絡或多次包絡，使用不同次數的包絡可實現不同的識別效果，一次包絡得到的顯著圖除了能識別出信號顯著變化的區域，還能表達出信號一些細微的幅值變化，而二次或多次包絡得到的結果便于觀察信號整體顯著性。

(3)對聽覺外周模型輸出進行多尺度二維濾波，可全面揭示信號中時程和頻寬不同的突變信號成分，這也說明對信號的帶通濾波結果進行較為深入的二次處理有助于提取信號中更多的特征信息，這也是聽覺系統信息處理機制與傳統信號處理方法的重要區別，具有較為樂觀的研究前景。

[ 1 ] 夏虹，陳志輝，鄧禮平,等. 基于小波的突變信號識別方法與應用研究[J]. 原子能科學技術，2012, 46(12): 1451-1456. XIA Hong, CHEN Zhihui, DENG Liping, et al .Mutation signal recognition method based on wavelet and application [J]. Atomic Energy Science and Technology, 2012, 46(12): 1451-1456.

[ 2 ] 張欣，夏秀渝，王雪君,等. 一種聽覺顯著圖提取模型[J].四川大學學報(自然科學報)，2014，51(2):292-298. ZHANG Xin, XIA Xiuyu, WANG Xuejun,et al. A saliency map extraction model for auditory attention[J]. Journal of Sichuan Univeisity(Natural Science Edition), 2013, 29(9): 1142-1147.

[ 3 ] 王雪君，夏秀渝，張欣,等. 新的聽覺注意顯著圖計算模型研究[J]. 信號處理，2013， 29(9): 1142-1147. WANG Xuejun, XIA Xiuyu, ZHANG Xin,et al. Research on a novel saliency map computational model of auditory attention [J]. Journal of SignalI Processing, 2013, 29(9): 1142-1147.

[ 4 ] 吳燕，隋光遠,曹曉華. 內源性注意和外源性注意的ERP研究[J].心理科學進展，2007, 15(1): 71-77. WU Yan, SUI Guangyuan, CAO Xiaohua . ERP studies of endogenous and exogenous attention[J]. Advances in Psychological Science, 2007, 15(1): 71-77.

[ 5 ] 李允公,張金萍，戴麗,等. 基于聽覺模型ZCPA的故障診斷特征提取方法研究[J].中國機械工程，2009，20(24): 2988-2992. LI Yungong, ZHANG Jinping, DAI Li, et al . Study on the ZCPA-auditory-model-based method of feature extraction for mechanical faults diagnosis [J].Chinese Journal of Mechanical Engineering, 2009, 20(24): 2988-2992.

[ 6 ] SHAMMA S.On the role of space and time in auditory processing[J].Trends in Cognitive Science，2001,5: 340-348．

[ 7 ] 陳世雄, 宮琴, 金慧君. 用Gammatone濾波器組仿真人耳基底膜的特性[J] . 清華大學學報(自然科學版), 2008, 46(6): 33-38. CHEN Shixiong, GONG Qin, JIN Huijun, Gammatone filter bank to simulate the characteristics of the human basilar membrane [J]. J. Tsinghua Univ. (Sci. & Tech.),2008, 46(6): 33-38.

[ 8 ] LIU H Y, HUANG W G, WANG S B, et al . Adaptive spectral kurtosis filter based on Morlet wavelet and its application for signal transients detection [J].Signal Processing, 2007,21(6): 2616-2633.

[ 9 ] 趙鶴鳴,王永琦,陳雪勤.聽覺模型反演方法及其應用[J].聲學學報,2005, 30(6): 530-535. ZHAO Heming, WANG Yongqi, CHEN Xueqin . Auditory model inversion and its application [J].Acta Acustica, 2005, 30(6):530-535.

[10] 李智,陳祥初,劉政波. 包絡分析及其在設備故障診斷中的應用[J].測試技術學報，2006, 12(6): 92-95. LI Zhi, CHEN Xiangchu, LIU Zhengbo. Envelopment analysis and its application in the fault diagnosis[J]. Journal of Test and Measurement Technology, 2006,12(6): 92-95.

[11] 馮新星,張麗艷,葉南,等.二維高斯分布光斑中心快速提取算法研究[J]. 光學學報，2012，32(5):78-85. FENG Xinxing, ZHANG Liyan, YE Nan, et al. Fast aigorithms on center location of two dimensional gaussian distribution spot [J] . Acta Optica Sinica, 2012, 32(5):78-85.

[12] 劉揚,張苗輝,鄭逢斌.聽覺選擇性注意的認知神經機制與顯著性計算模型[J]. 計算機科學, 2013, 40(6): 283-287. LIU Yang, ZHANG Miaohui, ZHENG Fengbin. Cognitive neural mechanisma and saliency computational model of auditory selective attention [J]. Compututer Science, 2013,40(6):283-287.

[13] ITTI L,KOCH C．Computational modelling of visual attention[J]. Nature Reviews Neuroscience，2001，2(3): 194-203．

Calculation method for signal saliency map based on human auditory characteristics

LI Guomeng, LI Yungong, WANG Bo, WU Wenshou, AN Chao

(Institute of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China)

A signal amplitude saltation often implies fault information of an equipment. In order to recognize the time-frequency domain where a signal amplitude saltation appears, a calculation method of signal saliency map based on the information processing mechanism of an auditory system was proposed. Firstly, auditory peripheral processings, such as, band-pass filtering, phase adjustment and half wave rectification were performed for a signal. Secondly, extracting the envelope information from the results of processings one time or multiple times, they were then processed with the multi-scaled two dimensional filtering. Nextly, the local saliency of the signal under different scales was obtained using the centre-rim difference operator. Finally, the global saliency map and time-frequency saliency map of the signal were gained by means of the multi-scaled integration and the linear combination of local saliencies. The results of tests and simulations showed that the proposed method is applicable to both vibration signals and voice signals, it is effective and practical, it can clearly represent amplitude saltation time-frequency domain of a signal and resist a certain noise interference.

auditory model; saliency map; fault diagnosis; two-dimensional filtering; envelope analysis

國家自然科學基金(51275080)

2016-01-15 修改稿收到日期：2016-06-12

李國萌 女，碩士生，1990年生

李允公 男，博士，副教授，1976年生

TN911.7

A

10.13465/j.cnki.jvs.2017.03.010