曾 祥， 吳學雷， 周曉軍， 楊辰龍， 陳越超

(1. 浙江大學 流體動力與機電系統(tǒng)國家重點實驗室，杭州 310027; 2. 北京航天發(fā)射技術研究所，北京 100076)

基于S變換和時頻譜空間濾波的超聲缺陷回波檢測方法

曾 祥1， 吳學雷2， 周曉軍1， 楊辰龍1， 陳越超1

(1. 浙江大學 流體動力與機電系統(tǒng)國家重點實驗室，杭州 310027; 2. 北京航天發(fā)射技術研究所，北京 100076)

缺陷回波的檢測是超聲探傷的一項重要內容，為減弱噪聲的影響準確檢測缺陷回波，提出基于S變換時頻分析和時頻譜空間濾波的信號處理方法。討論了高斯回波模型下的到達時間和中心頻率與S變換時頻譜的關系，說明了利用S變換時頻譜幅值矩陣的極值提取回波到達時間和中心頻率的合理性；為檢測回波，首先對原始信號作S變換，然后對得到的時頻譜幅值矩陣應用最大熵法自適應選擇去噪閾值，對S變換時頻譜作空間濾波完成降噪；從降噪后的區(qū)域中提取反映缺陷的到達時間和中心頻率；對降噪后的時頻譜作S逆變換，獲得缺陷回波明顯的時域信號。仿真研究表明，基于S變換和時頻譜空間濾波的方法能夠有效去除噪聲，檢測回波。棒材試塊的實驗結果同樣表明了該方法在缺陷檢測上的有效性。

超聲；S變換；最大熵閾值；空間濾波

超聲脈沖反射法是一種重要的無損檢測方法。缺陷回波提供了材料內部的缺陷信息，到達時間(Time of Arrival, TOA)和中心頻率(Center Freguency, CF)是反映回波的兩個重要參數。由于噪聲的存在，缺陷回波在時域上容易被掩蓋，在時頻域上噪聲容易形成多個峰值。為此，國內外學者對含噪超聲信號中的缺陷回波檢測進行了研究[1]，采用了小波變換[2]、自適應濾波[3]、經驗模態(tài)分解[4]、盲源分離[5]等方法。

本文提出一種基于廣義S變換時頻分析和最大熵閾值分割算法的超聲缺陷回波檢測方法。首先對回波信號作廣義S變換，獲得時頻譜系數矩陣；然后對時頻譜系數矩陣作空間濾波降噪，基于最大熵法確定雙閾值并處理；最后從降噪后的時頻譜系數矩陣中提取回波參數，作S逆變換得到降噪后的缺陷回波信號。

1 基于S變換的回波參數提取

1.1 廣義S變換

S變換是短時傅里葉變換和小波變換的繼承和發(fā)展，由Stockwell等[6]在研究地球物理數據時提出。相比短時傅里葉變換，S變換具有多分辨率特性，即在信號低頻段具有高的頻域分辨率，在信號高頻段具有高的時間分辨率。相比連續(xù)小波變換，S變換保持了原始信號的相位信息。因此，S變換被廣泛應用于心電信號(Electrocardiogram, ECG)[7]和心音信號處理[8]、滾動軸承故障信號分析[9]、地震波分析[10]等領域。

(1)

當p=1時，即為S變換。當plt;1時，相對S變換窗寬變寬，時間分辨率下降，頻率分辨率上升；pgt;1時相反。調節(jié)因子p，可以實現時頻分辨率的靈活調整、能量優(yōu)化集中。為描述簡便，不論p取值，統(tǒng)稱為S變換。

1.2 回波參數提取原理

單個高斯回波為

s0(τ)=βe-α(τ-b)2cos[2πfc(τ-b)+φ]

(2)

用復數表示為

s0(τ)=βe-α(τ-b)2ej[2πfc(τ-b)+φ]

(3)

信號能量Es0為

(4)

式中：β為信號強度；α為帶寬因子；b為到達時間；fc為中心頻率；φ為初相位。回波為窄帶信號，99%以上的能量集中在非指數項，滿足f2c≥0.24α[12]。

對s0(τ)作S變換，得到幅值矩陣元素

(5)

其中,

(6)

(7)

(8)

一般認為，回波的到達時間b和中心頻率fc可以通過求|S(t,f)|的極值位置提取。令?|S(t,f)|/?t=0，有

(9)

t=b與回波的其他參數和S變換的窗寬因子p無關。因此，S變換非常適合對超聲回波信號的TOA提取。

(10)

容易驗證t=b，f=fc時A+B=0且?(A+B)/?f= 0，有

(11)

可見當t=b，f=fc時，若p=0(Gabor變換)，?W/?f=0，極值點坐標滿足f=fc，與現有結論一致[13]。若pgt;0，有?W/?fgt;0，可以推測極值點坐標fgt;fc。因此|S(t,f)|取極值時的頻率高于信號的實際CF。盡管S變換對CF的提取有偏移，但它的多分辨率分析能力使得它在實際信號分析中應用廣泛。

為定量分析頻率的偏移量與其他參數的依賴關系，令?W/?f=0，并將t=b代入，可以得到

pqw2+πf(α+q)w+pα(α+q)/2=0

(12)

式中：w=π(fc-f)=πΔf；q=f2p/2。求解得到

(13)

Δf為α、f、p的函數且Δflt;0，如圖1(a)。從上到下的4個曲面依次對應p=0.25、p=0.5、p=0.75和p=1。圖1(b)為p=1時情形，從左到右10條實線依次對應f為1 MHz, 2 MHz, …, 10 MHz的Δf(虛線表示不滿足f2c≥0.24α條件)。

(a) Δf隨α、f、p的值

(b) p=1時的Δf隨α、f的值圖1 CF偏差分析Fig.1 Bias analysis of CF

偏差Δf的變化趨勢為隨p的增大和α的增大而增大，隨f增大而減小。當p、f均確定時，Δf隨α近似成正比例下降。

可見由極值條件得到的f是CF的近似值。取最大相對誤差ef=Δf/(f-Δf)×100%，由圖1(b)可知，p=1時ef約為15%。實際中高頻回波由于α有限，ef將降低。如α∈[0,80]，f為10 MHz的最大偏差約為0.4 MHz，ef=-0.4/(10-0.4)×100%=-4.2%。因此以f近似地作為CF能夠滿足實際檢測要求。如果需要CF的精確值，可以考慮Gabor變換，求解優(yōu)化問題得到TOA、CF，精度較高，但需要額外計算；減小p可實現增大頻率分辨率、降低Δf和ef，但時間分辨率將下降。

若對S變換的窗函數作能量歸一化，此時需在式(5)的前面乘以能量歸一化因子η

(14)

同樣可以得到極值條件下TOA沒有偏差，CF存在偏差

(15)

兩種情形下Δf相差很小。讀者可自行分析，與圖1(b)對比。

考慮實際檢測中，本文不作能量歸一化，直接按極值條件確定TOA、CF，且p=1。值得注意的是，本文忽略了實際信號由于采樣頻率、信號長度有限造成的偏差和不同p值下時頻分辨率不同造成的偏差。

2 時頻譜空間濾波

2.1 問題轉換

對x(t)作S變換，得到的時頻譜系數矩陣S是復數矩陣，它的幅值矩陣反映了信號能量在不同時間段、頻率段的分布情況，可以視作一幀像素灰度值可為任意非負值的灰度圖像I，可以采用數字圖像的空間濾波方法進行處理。這樣，由于時頻圖的“可視化”效果，相比小波閾值去噪和模極大值去噪，過程逐步直觀。

閾值法是常用的降噪方法，通過選擇合理的閾值T，將S分成回波時頻譜和噪聲時頻譜兩部分。該過程相當于對灰度圖像I用閾值T作分割，結果為I與J的乘積，其中J為乘積模板，滿足

(16)

若J(i,j)=1，則稱為前景；反之為背景。初始時J為零矩陣。采用單閾值時，在閾值附近的值處于“模糊”狀態(tài)，為保留較多信息采用雙閾值。設閾值TH、TL(THgt;TL)，分別得到JH、JL。JH中的前景屬于回波，JL的背景全為噪聲，而JH的背景和JL的前景的交集C可以利用灰度在像素空間的分布信息判定。將C中的所有判定得到的前景點加入到JH即為J。強噪聲下J中可能出現小面積前景，屬于偽回波區(qū)域，需要去除。

將時頻譜系數矩陣S與J逐點相乘，即可得到降噪后的時頻譜系數矩陣Sn。

2.2 雙閾值選擇

灰度圖像的自適應全局閾值分割算法有大津法(OTSU)[14]、最大熵法等，其中最大熵法對不同信噪比的圖像均能產生較好的分割效果。本文采用最大Kapur熵法[15]確定分割閾值。

設灰度圖像I有L個灰度級，閾值T將其分割為背景A和前景B兩部分。A區(qū)域像素的灰度值k1,k2, …,km均不高于T，B區(qū)域像素的灰度值km+1,km+2,…,kL均高于T。再記I的總像素數為N，第i個灰度級的像素數為Ni。那么A、B區(qū)域像素灰度的概率分布為

式中：pi=Ni/N為第i個灰度級的概率，而

A、B部分的熵計算為

(18)

總的熵H為

H=H(A)+H(B)

(19)

H取得最大值意味著獲取的前景區(qū)域A和背景區(qū)域B的總信息量最大，對應的閾值即為最佳全局閾值。將此閾值作為高閾值TH。

認為I中小于TH的其他元素基本由噪聲形成，計算它們的均值TM，將低閾值TL取作TH和TM的均值，即TL=(TH+TM)/2。

由此TH、TL均被自適應確定。應用TH、TL對I按式(16)進行二值化分別得到JH、JL。

2.3 交集處理和偽回波區(qū)域去除

對C中像素，若為回波，應該“環(huán)繞”著JH的前景且距離不能過遠。對JH作距離變換，計算各點到最近前景點的最短歐氏距離，結果記為DJH。給定最大允許距離D≥ 0，任意像素C(i,j)，若滿足DJH(i,j)lt;D，則認為對應回波取作前景，即令JH(i,j)=1；反之則認為對應噪聲。顯然，D=0時C中的點均被認為是噪聲，相當于用高閾值TH進行單一閾值去噪；D足夠大時，C中的點均被認為是回波，相當于用低閾值TL進行單一閾值去噪。

若JH存在偽回波小區(qū)域，可給定系數κ，對所有N個區(qū)域面積從大到小排序，記A(n)為區(qū)域n的面積，記面積比qn=A(n)/A(n-1)，若qm≤κ均成立(n=N,N-1,…,2;m=n,n+1,…,N)，則將區(qū)域n設為背景?；蛘呓涍^有限次腐蝕后再做有限次膨脹。實際中小區(qū)域面積和回波區(qū)域往往相差明顯，可以直接指定系數κ，本文采用這種方法。最終得到的二值矩陣即為乘積模板J。

3 基于S變換和時頻譜空間濾波的回波檢測

從受到噪聲干擾的信號x(τ)中檢測缺陷回波的步驟如下：

步驟1 S變換。對x(τ)作S變換得到時頻譜系數矩陣S，由S計算譜系數幅值得到幅值矩陣I。記I中元素的最小值、最大值分別為u、v。

步驟2 低通濾波和雙閾值選擇。估計回波的上限頻率fu，取I中頻率低于fu的區(qū)域作為感興趣區(qū)域(Region of Interest, ROI)。選擇迭代步長Δ，計算閾值為T0=u+ (i-1)Δ時的熵Vi，其中正整數i為計算次數。遞增i直到T0≥v結束迭代。得到閾值序列{T0}和對應的熵序列{H}。

步驟3 閾值化。計算雙閾值TH、TL，根據TH、TL對I執(zhí)行閾值化，得到JH、JL；處理交集C，選擇D將C

中距離小于D的區(qū)域加入JH的前景；如有必要，指定系數κ，去除JH的小區(qū)域。最終得到乘積模板J。

步驟4 空間濾波。將時頻譜系數矩陣S與J逐點相乘，得到降噪后的時頻譜系數矩陣Sn。

步驟5 回波參數提取和信號重建。從Sn中提取TOA、CF。對Sn作S逆變換，得到降噪后的信號x1(τ)。

4 仿真研究

為驗證提出的超聲缺陷回波檢測算法的有效性，構造仿真信號x(t)進行研究。設具有K重回波的超聲信號x(t) =s1(t) +…+sK(t)+n(t)。n(t)為零均值，方差為σ2的高斯白噪聲。記E0為回波信號的能量，信噪比SNR單位dB，SNR=10lg(E0/σ2)。

考慮到本文方法φ對回波檢測沒有影響，所以仿真中忽略φ。對單個缺陷回波si(t)，設λi=[βiαibifci]T，取K=4，λ1=[1 15 3 7.5]T，λ2=[0.95 15 4 6.5]T，λ3=[0.90 12 6 5.5]T，λ4=[0.90 10 8 4.5]T構造s(t)，并向s(t)中添加SNR=-2 dB的噪聲信號n(t)，獲得含噪信號x(t)。仿真信號采樣頻率fs=50 MHz，選擇低通濾波頻率為20 MHz。

圖2 仿真信號處理結果Fig.2 Processing result of simulation signal

可見對低信噪比混合信號x(t)(見圖2(b))，S變換后缺陷回波信號明顯(見圖2(c))。經雙閾值THTL處理和交集處理和去除偽回波小區(qū)域后得到乘積模板J(見圖2(d))。J與S逐點相乘得到降噪后的時頻譜Sn(見圖2(e))，經S逆變換后重建時域波形(見圖2(f))。可見，重建的信號波形完備清晰?；赟tein無偏似然估計和極大極小方法選擇的閾值將部分噪聲信號誤認為缺陷回波信號予以了保留；而固定閾值的降噪方法去除了較多有效信息，得到的缺陷回波近似呈現脈沖狀(見圖2(g))。另外，由于本文方法僅保留了信號能量較高的時間頻率段成分，因此降噪后的信號幅值降低；但對于強噪聲下的缺陷檢測和定位并無影響。可以預見，弱噪聲下，由于閾值的降低，降噪后的信號幅值降低較小。

為定量評價本文回波特征提取方法的有效性，對原始信號和處理后的信噪比(Signal-to-Noise Ratio, SNR)、均方根誤差(Root-Mean-Square-Error, RMSE)列入表1。TOA、CF的提取結果列入表2。

表1 信噪比和均方差

表2 參數提取結果

由表1可知，相比給出的3種小波閾值去噪方法，本文方法的信噪比提升明顯，均方差較小。其他3種方法，采用固定閾值和極大極小閾值的信噪比較高，均方差較小。由表2可知，經S變換可以方便地獲取回波信號的TOA和CF，且偏差較小。偏差的來源主要有3方面，包括時頻分辨率的限制、噪聲與回波在時頻域上的混疊，以及S變換極值處頻率對CF的偏移。但總體而言，在混疊較小的情形下，TOA和CF的提取偏差較小。

為進一步比較本文方法和3種小波閾值方法在不同SNR下的降噪效果，改變仿真信號中加入噪聲的強度，取SNR = 4 dB，2 dB，0 dB，-2 dB，-4 dB。為減弱處理的隨機性，對每個SNR重復3次，計算3次處理后的SNR和MSE的平均值作為最終結果，如圖3所示。

可見，相比使用的3種小波閾值去噪方法，本文方法去噪的效果更好。因此，數字圖像處理技術應用于回波信號去噪是可行的。事實上，本文方法與小波閾值去噪法是相通的，都是在時頻域上去除或減弱幅值低于閾值的局部成分。借助于數字圖像處理技術，本文方法過程更加直觀。

(a) 處理前后的SNR

(b) 處理前后的MSE圖3 不同SNR下降噪效果比較Fig. 3 Comparison of denoising effects with different SNR

5 實際回波檢測

現在將S變換和時頻譜系數矩陣空間濾波方法應用于實際缺陷回波檢測。如圖4所示，采用20#鋼棒材試塊，公稱直徑為Φ120 mm，實測數據為Φ119.94 mm。3個平底孔直徑依次為Φ0.8 mm、Φ1.2 mm、Φ2.0 mm作為人工缺陷。本文選取Φ1.2 mm孔作為代表進行檢測，其公稱尺寸為Φ1.2×埋深90 mm，實測數據為Φ1.22×埋深89.9 mm。材料聲速事先已經過測定為c=5 934 m/s。超聲探頭中心頻率為7.5 MHz，信號采樣頻率為100 MHz。始波和底波之間為背散射信號s(t)，由平底孔缺陷回波s1(t)、結構噪聲和微小缺陷散射反射回波s2(t)、噪聲n0(t)組成，三者幅度依次減小?，F在需檢測s1(t)。選擇Turkey-Hanning窗函數對s(t)進行截取，截取的采樣點區(qū)間為1 501～5 000。處理結果如圖5所示。

圖4 20#棒材試塊Fig. 4 20# bar specimen

(a) Turkey-Hanning窗函數

(b) 實測信號和截取后的背散射信號

(c) S變換時頻圖

(d) 還原的信號圖5 20#棒材試塊人工缺陷回波檢測Fig. 5 Echo detection of artificial flaw in 20# bar specimen

從S變換時頻圖中可以看到明顯的缺陷區(qū)域(見圖5(c))，對應平底孔缺陷，回波在圖5(c)的坐標為(2 650, 8.086)。在圖5(b)實測信號中時間點為4 201。始波幅值近似達到5V的首末時間點分別為1 107和1 167，取始波位于均值1 137處。計算得到平底孔的位置為：d=cΔt/2=cN/(2fs)=1 000×5 934×(4 201-1 137)/(2×100×106)=90.9 mm，與實測埋深89.9 mm偏差很小，表明TOA的提取是較準確的。還原的信號中缺陷回波更加明顯，幾乎無噪聲，這是因為超聲在金屬材料中衰減較小，平底孔反射強，從而缺陷回波在時頻圖上能量集中度很高 (見圖5(c))，CF相比探頭中心頻率上升了約0.5 MHz。

如果需要進一步降噪研究s2(t)，可以用s(t)的S變換時頻譜減去s1(t)的S變換時頻譜，差值近似作為s2(t)和n0(t)的疊加信號的S變換時頻譜(因為s1(t)的檢測存在誤差)，再進行時頻譜系數空間濾波。最終結果如圖6所示。

圖6 信號成分s2(t)檢測Fig. 6 Detection of signal component s2(t)

6 結 論

(1) S變換是一種優(yōu)秀的時頻分析方法，得到的時頻譜系數矩陣能夠很好地揭示信號的時頻特性，可用來提取TOA和CF，適合超聲缺陷回波檢測。

(2) 基于最大熵法的閾值選擇具有自適應、分割優(yōu)良的優(yōu)點，能夠較好地去除噪聲，保留回波信號。

(3) 將數字圖像處理技術應用于S變換的時頻譜系數矩陣的閾值降噪是有效的，且過程直觀可見?？梢哉J為，如短時傅里葉變換、連續(xù)小波變換得到的時頻圖也可以視作灰度圖像進行降噪處理。

[ 1 ] ZHANG G M，HARVEY D M．Contemporary ultrasonic signal processing approaches for nondestructive evaluation of multilayered structures[J]．Nondestructive Testing and Evaluation，2012，27(1)：1-27．

[ 2 ] MATZ V，SMID R，STARMAN S，et al．Signal-to-noise ratio enhancement based on wavelet filtering in ultrasonic testing[J]．Ultrasonics，2009，49(8)：752-759．

[ 3 ] 楊克己. 基于神經網絡的自適應濾波技術及其在超聲檢測中的應用[J]．儀器儀表學報，2005，26(8)：813-817．

YANG Keji．An adaptive filter based on artificial neural net and its application in ultrasonic testing[J]．Chinese Journal of Scientific Instrument，2005，26(8)：813-817．

[ 4 ] LU Y F, ORUKLU E, SANIIE J．Chirplet signal and empirical mode decompositions of ultrasonic signals for echo detection and estimation [J]．Journal of Signal and Information Processing，2013，4(2)：149-157．

[ 5 ] LIU Q K，QUE P W，GUO H W，et al．Noise cancellation of ultrasonic NDE signals using blind source separation[J]．Russian Journal of Nondestructive Testing, 2006, 42(1)：

63-68．

[ 6 ] STOCKWELL R G，MANSINHA L，LOWE R P．Localization of the complex spectrum：the S transform[J]．IEEE Transaction on Signal Process，1996，44(4)：998-1001．

[ 7 ] ARI S，DAS M K, CHACKO A．ECG signal enhancement using S-Transform[J]．Computers in Biology and Medicine，2013，43(6)：649-660.

[ 8 ] 李戰(zhàn)明，韓陽，韋哲，等. 基于S變換的心音信號特征提取[J]．振動與沖擊，2012，31(21)：179-183．

LI Zhanming, HAN Yang, WEI Zhe, et al．Heart sound feature extraction based on S tranformation[J]．Journal of Vibration and Shock，2012，31(21)：179-183．[ 9 ] 郭遠晶，魏燕定，周曉軍，等. 基于S變換譜閾值去噪的沖擊特征提取方法[J].振動與沖擊，2014，33(21): 44-50.

GUO Yuanjing, WEI Yanding, ZHOU Xiaojun, et al.An impact feature extracting method based on S transformation spectrum threshold denoising[J].Journal of Vibration and Shock，2014，33(21): 44-50.

[10] 樊劍，呂超，張輝. 基于S變換的地震波時頻分析及人工調整[J].振動工程學報，2008，21(4): 381-386.

FAN Jian, Lü Chao, ZHANG Hui.Time-frequency analysis and artificial simulation of earthquake groud motions via S-transform[J].Journal of Vibration Engineering，2008，21(4): 381-386.

[11] DJUROVIC I, SEJDIC E, JIANG J.Frequency-based window width optimization for S-transform[J].AEU-International Journal of Electronics and Comunications，2008，62(4): 245-250.

[12] DEMIRLI R, SANIIE J.Model-based estimation of ultrasonic echoes Part I: analysis and algrithms[J].IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2001，48(3): 787-802.

[13] LU Z K, YANG C, QIN D H.Estimating the parameters of ultrasonic echo signal in the Gabor transform domain and its resolution analysis[J]. Signal Processing, 2016, 120: 607-619.

[14] OTSU N.A threshold selection method from gray-level histograms[J].IEEE Trans on Systems，Man，and Cybernetics，1979，9(1): 62-66.

[15] KAPUR J N，SAHOO P K，WONG A K C.A new method for gray-level picture thresholding using the entropy of the histogram[J].Computer Vision，Graphics and Image Processing，1985，29(3): 273-285.

AflawechodetectionmethodbasedonS-transformationandtime-frequencyspectrumspatialfiltering

ZENG Xiang1, WU Xuelei2, ZHOU Xiaojun1, YANG Chenlong1, CHEN Yuechao1

(1. The State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou 310027, China;2. Beijing Institute of Space Launch, Beijing 100076, China)

Detection of flaw echoes is an important task in ultrasonic testing. In order to eliminate the effects of noise and detect echoes correctly, a signal processing method based on the S-transformation (ST) and the time-frequency spectrum (TFS) spatial filtering was proposed．Based on the Gaussian echo model, the relationships between the time of arrival (TOA) together with the center frequency (CF) and ST TFS were discussed. The reasonableness of the method using the extremum of ST TFS amplitude matrix for TOA and CF extraction were illustrated. To detect the echoes, firstly the ST was performed on the original signal．Then the dual-threshold were determined by the Maximum Entropy Thresholding method adaptively, and the spatial filtering was performed on ST TFS for denoising. The TOA and CF could be extracted from the denoised region. Finally, the inverse ST was performed on the denoised TFS and the signal with clear flaw echoes were gained. Simulation results show that using ST and TFS spatial filtering can remove noise and detect echoes effectively. And the experimental results of bar specimen also show the effectiveness of the method in flaw detection.

ultrasonic; S-transformation; maximum entropy thresholding; spatial filtering

浙江省自然科學基金(LY14E050013)；浙江省公益技術研究工業(yè)項目(2015C31052)；重慶齒輪箱有限責任公司"海面平臺洋流發(fā)電裝備研制與開發(fā)"

2016-03-29 修改稿收到日期： 2016-07-14

曾祥 男，博士生，1989年生

周曉軍 男，教授，博士生導師，1958年生

TB553

A

10.13465/j.cnki.jvs.2017.22.006