基于雙高斯衰減模型的超聲回波處理方法*

王大為 王召巴 陳友興 李海洋 王浩坤

1) (中北大學信息與通信工程學院,太原 030051)

2) (山西師范大學物理與信息工程學院,臨汾 041000)

信號降噪與特征提取是超聲檢測數據處理的關鍵技術.基于超聲信號有特定結構而噪聲和超聲信號的結構無關,本文提出一種旨在解決強噪聲背景下超聲回波的參數估計和降噪問題的方法.該方法將超聲回波的參數估計和降噪問題轉換為函數優化問題,首先根據工程經驗建立超聲信號的雙高斯衰減數學模型,然后根據觀測回波和建立的超聲信號模型確定目標函數,接著選擇人工蜂群算法對目標函數進行優化從而得到參數的最優估計值,最后由估計出的參數根據建立的超聲信號數學模型重構出無噪的超聲估計信號.通過仿真和實驗表明本文方法可以準確估計出信噪比大于—10 dB的含噪超聲回波中的無噪信號,且效果優于基于自適應閾值的小波降噪方法和經驗模態分解方法; 此外相比常用的指數模型和高斯模型,本文提出的雙高斯衰減超聲信號模型與實測超聲信號更接近,其均方誤差為9.4 × 10—5,波形相似系數為0.98.

1 引言

超聲檢測是重要的無損檢測技術之一,也是目前國內外使用頻率最高且發展較快的一種無損檢測技術,在工業生產等實踐中得到了廣泛應用[1,2].信號降噪與特征提取是超聲檢測數據處理的關鍵技術,同時也是超聲無損檢測的核心環節,其性能優劣直接影響著無損檢測結果的表征,因此對其進行研究有重要的學術意義和實用價值[3].Wang等[4]為評價點焊接頭強度,首先利用快速傅里葉變換和小波包變換對檢測到的超聲信號進行處理,然后利用數理統計方法提取超聲檢測數據的時域、頻域和小波域特征,再用支持向量機的分類器對提取的特征進行分類.Meng等[5]為分類碳纖維增強復合材料中的孔洞和分層,首先對采集到的超聲信號進行小波分解,然后利用深度卷積神經網絡從小波系數中學習每個信號緊湊而有效的表示,最后利用支持向量機的方法對學習結果進行分類.孫靈芳等[6]針對換熱管薄層污垢超聲檢測信號降噪問題提出了一種基于改進完備總體經驗模態分解的小波收縮閾值信號處理方法.該方法首先引入夾角余弦計算原始信號和固有模態函數的相似程度,判斷信號和噪聲主導模態的分界點,然后利用小波收縮閾值方法提取噪聲主導模態中的細節信息,最后重構得到降噪后的超聲信號.這些方法的共同特點是先對檢測到的超聲回波降噪,然后再進行特征提取和表征.對于信噪比相對較高的超聲回波,這些方法可以取得期望的處理效果,但對于信噪比較低的超聲回波,例如回波中信號被噪聲淹沒時,這些方法難以對超聲回波中的特征進行提取和表征.為解決該問題,本文作者在前期工作中提出了一種強噪聲背景下微弱超聲信號提取方法[7],然而該方法對實測超聲信號處理結果的優劣依賴于所建立的超完備字典中是否包含能夠和實測超聲信號相匹配的原子.為增強對實測超聲信號的處理能力,本文提出雙高斯衰減超聲信號模型,并基于該模型提出一種超聲回波處理方法.該方法將超聲信號處理問題轉換為對含參函數的優化問題,通過利用人工智能優化算法對目標函數求解而實現從觀測到的含噪回波中估計無噪信號參數.

2 問題描述

在超聲脈沖回波檢測中,觀測到的超聲回波y(t)是測試樣本中各界面反射回波sk(t)的線性組合[8],即

其中sk(t)為發射脈沖遇到第k個界面后的反射回波信號,ak為對應的反射系數,ξ(t)表示來自于測試系統和被測材料中的噪聲.由于實際測量到的觀測回波y(t)中通常包含噪聲,本文主要解決如何從包含噪聲的觀測回波中提取超聲信號特征參數及重構無噪超聲信號.根據工程經驗每一個超聲回波sk(t)都是由超聲傳感器調制發射的,故其有特定的結構并可以表示為一組含參數的非線性函數;而噪聲ξ(t)通常是一個和超聲信號結構無關的高斯隨機過程[9-11]因此,在超聲無損檢測中由觀測回波估計出無噪原始信號的問題可以轉化成函數在參數集上的優化問題.本文通過四個環節來解決該問題,即建立超聲信號模型、確定目標函數、優化目標函數并估計參數和重構信號.流程圖如圖1所示.

3 本文方法

3.1 超聲回波模型

圖1 超聲信號處理流程圖Fig.1.Flow chart of ultrasonic signal processing.

在超聲檢測中,超聲回波通常是一個被超聲換能器中心頻率調制的寬帶信號.當前,超聲回波數學模型有指數模型[10]和高斯模型[11],其包絡分別按照指數規律和高斯規律變化.指數模型和高斯模型在超聲換能器特性模擬和超聲信號數值分析方面有重要作用,并得到了廣泛應用.但實際中超聲換能器發射的超聲脈沖包絡并非嚴格的指數包絡或高斯包絡,而是產生了非線性失真.為提高本文方法對實測超聲信號參數估計的準確性,希望建立的超聲信號數學模型包絡盡可能地接近超聲換能器發射脈沖的包絡.基于通過兩個不同中心位置的高斯窗函數疊加來模擬高斯模型和發射脈沖包絡之間存在的非線性失真,本文提出超聲信號的雙高斯衰減模型,其數學表達式如下：

式中參數為 θ=[A,a,b,T,τ,fc,φ].其中 A 是超聲回波幅度; a,b和T與超聲換能器特性有關,決定超聲回波的包絡,因此定義為包絡參數; 包絡參數有直觀的物理意義,a,b分別是第一個和第二個高斯包絡的中心位置,T是高斯包絡的尺度參數,其決定著高斯包絡的時寬; τ為渡越時間; fc是超聲發射脈沖的中心頻率,φ是初相位.雙高斯衰減模型是高斯模型的拓展,當a=b時,該模型退化為高斯模型.不同包絡參數[a,b,T ]對應的雙高斯包絡如圖2所示.

從圖2中可以看出,本文提出的雙高斯衰減模型的包絡結構具有可塑性,本文正是利用該特性模擬超聲換能器發射脈沖相對于高斯模型的非線性失真.

3.2 目標函數

目標函數的選擇直接影響模型參數估計的準確性.文獻[8]中采用內積作為目標函數評價超聲回波和字典中原子的相似度; 文獻[10]中以均方誤差函數作為目標函數采用指數模型估計超聲回波渡越時間.在本文的應用背景中以內積或均方誤差函數作為目標函數存在目標函數最優收斂值隨被處理回波信噪比變化的問題,不利于判斷優化算法是否正確收斂.余弦相似度是一種有效衡量數據序列相似性的方法[12],本質上是對內積的歸一化; 此外,以余弦相似度作為目標函數,則目標函數值是否為1可成為判斷優化算法是否正確收斂的標準.設f,g分別為兩個長度相同的離散時間序列,則余弦相似度定義為

圖2 超聲回波模型包絡 (a) [a,b,T ]=[1.5,2.5,1]; (b) [a,b,T ]=[3,3,1]; (c) [a,b,T ]=[1.7,3.2,1]; (d) [a,b,T ]=[1,3,1]Fig.2.Envelope of ultrasonic echo model： (a) [a,b,T ]=[1.5,2.5,1]; (b) [a,b,T ]=[3,3,1]; (c) [a,b,T ]=[1.7,3.2,1]; (d) [a,b,T ]=[1,3,1].

當兩個序列有相同的形狀時余弦相似度為1,當兩個序列無關時余弦相似度為0.在[0,1]區間上,余弦相似度值越大表明兩個序列形狀越接近.設由測試系統觀測到的回波為y(t),待估計的超聲信號為s(θ; t),則觀測回波和待估計超聲信號的余弦相似度可表示為

為分析問題簡單起見,假設y(t),s(θ; t)是經過范數歸一化處理的,同時利用內積運算的線性性質得

由于超聲信號有特定的結構,而噪聲和超聲信號的結構無關,所以理論上噪聲和超聲信號不相關,故等式右邊第二項為零,即

其中bk是范數歸一化系數.因此由觀測回波y(t)求超聲信號s(θ; t)可以轉化為在參數空間θ=[A,a,b,T,τ,fc,φ]上對(8)式所示目標函數的最大值優化,

綜上,采用余弦相似度作目標函數時,從包含多個回波的觀測回波中提取超聲信號特征參數及重構無噪超聲信號的問題可以轉化為對多個單個回波參數估計問題的線性組合.而對從單個回波中提取超聲信號特征參數的問題即是在θ張成參數空間中尋找一組參數使得目標函數(9)式取最大值,

顯然,該目標函數的理論最大值為1.

3.3 目標函數優化及信號重構

很多優化算法都可以用來求解本文的目標函數,例如模擬退火算法[13-15]、粒子群優化算法[16-18]、人工蜂群算法[19-21]等.不失一般性,選擇人工蜂群(artificial bee colony,ABC)算法優化本文目標函數并以目標函數作為評價每個解的適應度函數,其基本過程如下.

首先在由θ張成的參數空間中隨機生成初始解xi(i=1,2,···,SN ),SN是解的個數,每個解xi是一個與θ結構相同的d維向量.根據(10)式計算每個解的適應度值.初始化之后整個種群將進行引領蜂、跟隨蜂和偵查蜂搜尋過程的重復循環,直到達到設定的最大迭代次數,

1) 引領蜂階段

引領蜂對當前解進行鄰域搜索,產生新解并貪婪選擇較優解.每個引領蜂由(11)式產生一個新解,

其中,i∈{ 1,2,···,SN },j∈{ 1,2,···,d },k 是隨機產生的整數,k∈{ 1,2,···,SN }且 k≠i; φij為[—1,1]之間的隨機數.檢驗新解的適應度值并根據貪心原則保留優質解,即若新解的適應度值比舊解大,則引領蜂保留新解舍棄舊解; 反之,保留舊解.在所有引領蜂完成搜索過程之后會把解的信息及適應度值信息分享給跟隨蜂.

2) 跟隨蜂階段

跟隨蜂根據(12)式計算每個解的選擇概率,

然后在[0,1]區間內產生一個隨機數,如果解的概率大于該隨機數,則跟隨蜂由(11)式產生一個新解,并根據貪心原則保留適應度值大的解.

3) 偵查蜂階段

在所有跟隨蜂完成搜索過程之后,如果一個解通過有限次循環不能被進一步改良,則該解被舍棄.設解xi被舍棄,則此解對應的引領蜂將變成一個偵查蜂.此偵查蜂由(13)式產生一個新解代替被舍棄解,

其中lj和uj分別是變量xij的下界和上界.然后返回引領蜂搜索過程,開始重復循環.

利用ABC算法對本文目標函數進行優化的基本步驟如下：

Step 1隨機初始化種群xi,i=1,2,···,SN ;

Step 2根據(10)式計算種群中各初始解xi的適應度值fit(xi);

Step 3引領蜂根據(11)式產生新解vi,并計算適應度值fit(vi),根據貪心原則保留適應度值大的解;

Step 4跟隨蜂根據(12)式計算解xi的選擇概率p(xi),依輪盤賭方式選擇解xi;

Step 5若解xi被選擇,根據(11)式產生新解vi,并計算新解的適應度; 跟隨蜂根據貪心原則保留優質解;

Step 6記錄最優解,即適應度值最大的解θbest;

Step 7偵查蜂決定是否存在需要放棄的解,如果存在,根據(13)式隨機產生一個新解替代舊解;

Step 8重復Step 3至Step 7,直到達到最大循環次數.

經人工蜂群算法優化后得到最優解θbest=[A,a,b,T,τ,fc,φ],不僅能直接獲取超聲回波的特征參數,還可以利用特征參數根據(3)式重構出超聲信號s(θbest; t).

4 數值仿真與分析

目前,廣泛用于評價仿真信號的指標有均方誤差(mean square error,MSE)、估計信噪比(estimated signal-to-noise ratio,ESNR)和波形相似系數(normalized correlation coefficient,NCC)[22,23].其定義如下：

4.1 超聲信號仿真

為定量分析本文算法對超聲回波的處理能力,以(3)式中的超聲信號與3個不同強度的脈沖序列進行卷積模擬發射的超聲脈沖通過三層不同介質界面后產生的超聲回波.無噪超聲信號及加入噪聲后的超聲回波如圖3所示.

圖3 仿真超聲回波 (a)理想超聲信號; (b)含噪超聲回波Fig.3.Simulation of ultrasonic echo： (a) Ideal ultrasonic signal; (b)noisy ultrasonic echo.

圖3(a)所示超聲信號由三個回波組成,第一個回波渡越時間為 1 00 μs,第二個回波渡越時間為500 μs,第三個回波渡越時間為 8 00 μs.但在工程實際中獲取的觀測回波通常是如圖3(b)所示的被噪聲污染的信號,經計算其信噪比為5 dB,從圖中已經無法準確分辨出第二個回波和第三個回波的渡越時間.

4.2 渡越時間估計

利用本文方法可以直接提取超聲回波渡越時間.對圖3(a)所示的信號,給其加入不同大小的噪聲,利用本文方法分別對不同信噪比的超聲回波渡越時間進行估計,結果列于表1.

表1 渡越時間τ估計誤差Table 1. Estimation error of transit time τ.

從表1可以看出,在相同信噪比下估計第一個回波渡越時間的絕對誤差小于第二個回波,第二個回波小于第三個回波.這是因為三個回波幅值依次減小,第三個回波最先被噪聲淹沒(見圖3(b)).隨信噪比的降低估計誤差增大,但在信噪比低至—10 dB時,估計第一個回波渡越時間的絕對誤差為1.3809 μs,相對誤差為1.3809%,估計第二個回波渡越時間的絕對誤差為 4 .4061 μs,相對誤差為0.8812%,由于噪聲太大,無法估計出第三個回波的渡越時間.

4.3 降噪性能分析

為了驗證本文方法的降噪能力,將本文方法和當前超聲信號處理中常用的小波閾值降噪方法、經驗模態分解(empirical mode decomposition,EMD)方法進行對比.為方便研究,不失一般性地選擇圖3(a)中的第二個回波作為研究對象對其加入噪聲,當信噪比SNR=—10 dB時,本文方法、小波閾值降噪方法以及EMD方法的降噪結果如圖4所示.

如圖4所示,當SNR=—10 dB時,超聲信號已完全被噪聲淹沒,采用本文方法可以不失真地重構出超聲回波,該回波的渡越時間為 5 00 μs ; 而小波方法和EMD方法已產生失真,從圖中無法準確分辨出回波的渡越時間.為進一步驗證本文方法對超聲信號降噪處理的有效性,對上述原始信號加入不同信噪比的噪聲后,分別采用本文方法、小波方法和EMD方法進行降噪處理,得到的降噪性能評價指標結果如表2所列.

圖4 不同方法降噪結果對比 (a)黑色實線為含噪回波,紅色虛線為原始信號; (b)本文方法處理結果; (c)小波方法處理結果;(d) EMD方法處理結果Fig.4.Comparison of de-noising by different methods： (a) Noisy echo plotted in a black solid line,original signal plotted in a red dotted line; (b) signal de-noised by our proposed method; (c) signal de-noised by wavelet method; (d) signal de-noised by EMD method.

均方誤差MSE反映了原始信號和去噪后信號之間的平均誤差,波形相似系數NCC衡量原始信號和去噪后信號之間的波形相似度,估計信噪比ESNR反映了降噪后信號的凸顯程度.MSE越小,ESNR越大,NCC越接近與1,說明降噪效果越好[23].從表2中可以看出,信噪比較高時本文方法、小波方法和EMD方法都可以取得較好的降噪結果,在信噪比較低時本文方法降噪結果明顯優于小波方法和EMD方法.

5 實測信號驗證

在實驗室中用RITEC公司生產的RITECRAM-5000-SNAP超聲檢測系統采集了用透射法測試金屬構件拉伸疲勞實驗的超聲回波進行處理,從而驗證本文方法是否可以有效處理實測超聲信號.如圖5所示,超聲回波采集裝置由RAMSNAP系統、匹配電阻、衰減器、低通濾波器、超聲換能器、示波器及計算機等組成.實驗中設置RAM-SNAP系統激發頻率為5 MHz、周期數為8的脈沖信號,經過匹配電阻和低通濾波器傳輸到超聲換能器上,再由換能器將電信號轉換成超聲波發射到被測金屬試件上表面,超聲波在被測試件內部傳播到達試件下表面后發生透射,透射的超聲波被接收換能器接收并轉換為電信號,最后將接收到的電信號在示波器上顯示.

表2 不同方法降噪結果對比Table 2. Comparison of denoising results using different methods.

圖5 超聲回波采集裝置Fig.5.Ultrasonic echo acquisition device.

5.1 超聲回波模型性能對比

為驗證本文提出的雙高斯衰減模型的有效性,將目前常用的指數模型、高斯模型及雙高斯衰減模型分別和采集到的實測超聲回波進行匹配,結果如圖6所示.

圖6 不同超聲信號模型對比Fig.6.Comparison of different ultrasonic signal models.

從圖6可以看出,相比指數模型和高斯模型,本文提出的雙高斯衰減模型與實測超聲回波匹配度更高; 此外,經計算本文模型和實測回波的均方誤差 MSE=9.4 × 10—5,波形相似系數 NCC=0.98; 指數模型和實測回波的均方誤差MSE=3.5 × 10—4,波形相似系數 NCC=0.93; 高斯模型和實測回波的均方誤差 MSE=2.9 × 10—4、波形相似系數NCC=0.94.因此,本文提出的雙高斯衰減模型和實驗室中采集到的實際超聲回波信號更接近.

5.2 實測回波處理

由于和工程實際相比實驗室環境中的噪聲很小,為驗證本文方法對實測超聲回波的降噪處理能力,將無激勵時測試到的系統噪聲放大后疊加到實測回波中以增加實測回波中的噪聲.實測超聲回波、加入噪聲的實測超聲回波及本文方法、小波方法和EMD方法處理結果對比如圖7所示.

圖7 實測超聲回波處理結果對比 (a)實測超聲回波;(b)混入系統噪聲后的實測超聲回波,信噪比為—10 dB;(c)本文方法結果; (d)小波方法結果; (e) EMD方法結果Fig.7.Comparison of measured ultrasonic echo processing results： (a) Measured ultrasonic echo; (b) measured ultrasonic echo contaminated by system noise; (c) echo processed by our proposed method; (d) echo processed by wavelet method; (e) echo processed by EMD method.

經計算,小波處理結果的MSE=1.6 × 10—3,NCC=0.62,ESNR=1.64 dB; EMD方法的MSE=2.6 × 10—3,NCC=0.65,ESNR=—0.49 dB;本文方法的 MSE=0.2 × 10—3,NCC=0.97,ESNR=11.58 dB.此外,從圖7(c)中可以分辨出本文處理結果回波渡越時間為 7 .29 μs,和圖7(a)中結果一致.因此在信噪比低至—10 dB時本文方法仍可以提取出實測含噪超聲回波中的超聲信號及對應的渡越時間,而此時小波方法和EMD方法已經失真.

6 結 論

基于超聲信號有特定結構而噪聲和超聲信號的結構無關,本文提出了一種超聲回波處理方法,并通過對仿真信號和實測回波的處理驗證了該方法能有效解決強噪聲中的超聲回波參數估計和降噪問題且處理效果優于小波方法和經驗模態分解方法.本文主要創新如下：

1) 將超聲回波的參數估計和降噪問題轉換為函數優化問題,通過建立超聲信號模型、確定目標函數、優化目標函數和根據優化結果重構無噪信號四步來實現信號參數估計和降噪;

2) 提出了雙高斯衰減超聲信號模型,相比于目前廣泛應用的指數模型和高斯模型,該模型產生的仿真超聲信號與實際測量到的超聲信號更接近;

3) 基于超聲信號有特定結構而噪聲和超聲信號的結構無關這一事實,在理論上證明了以余弦相似度作為目標函數可以將多回波信號估計問題轉化為多個單回波信號估計問題的線性組合,對于單回波估計問題可通過目標函數值是否為1判斷優化算法是否正確收斂.

同時還存在一些問題,例如在實驗中發現信噪比低于—10 dB時本文方法的穩定性和準確性都變差.這可能是由于采用的優化算法不能正確收斂導致的.開發適合解決本文目標函數的智能優化算法以增強本文方法處理效果將是我們下一步要重點研究解決的問題.