王大為，王召巴，李 鵬，陳友興，李海洋

(1.中北大學(xué)信息與通信工程學(xué)院，山西 太原 030051; 2.山西師范大學(xué)物理與信息工程學(xué)院，山西 臨汾 041000;3.西安近代化學(xué)研究所，陜西 西安 710065)

0 引 言

變分模態(tài)分解[1]（variational mode decomposition，VMD）是Dragomiretskiy于2014年提出的一種信號處理理論，該方法認(rèn)為任意信號都可分解為一系列帶寬有限的模態(tài)函數(shù)的線性組合，將信號處理問題轉(zhuǎn)化為變分模型求解問題，通過迭代搜索變分模型的極值確定每個模態(tài)的中心頻率和帶寬，從而實(shí)現(xiàn)對原信號的分解。變分模態(tài)分解因其有嚴(yán)密的數(shù)學(xué)理論推導(dǎo)，克服了經(jīng)驗(yàn)?zāi)B(tài)分解[2]（empirical mode decomposition，EMD）缺乏理論支撐、出現(xiàn)模態(tài)混疊和端點(diǎn)效應(yīng)等不足，一經(jīng)提出就引起了國內(nèi)外學(xué)者廣泛關(guān)注，已在生物組織健康診斷[3-4]、信號增強(qiáng)和辨識[5-6]、軸承故障診斷[7-8]等領(lǐng)域進(jìn)行了試驗(yàn)性研究和成功的應(yīng)用。超聲檢測是目前國內(nèi)外使用頻率最高且發(fā)展較快的一種無損檢測技術(shù)，在工業(yè)生產(chǎn)等實(shí)踐中得到了廣泛應(yīng)用[9-10]。信號降噪與特征提取是超聲檢測數(shù)據(jù)處理的關(guān)鍵技術(shù)，同時也是超聲無損檢測的核心環(huán)節(jié)，其性能優(yōu)劣直接影響著無損檢測結(jié)果的表征，因此對其進(jìn)行研究有重要的學(xué)術(shù)意義和實(shí)用價值[11]。

本文將變分模態(tài)分解方法引入到超聲檢測信號處理中，針對變分模態(tài)分解可以將待處理信號分解為一系列有限帶寬的模態(tài)函數(shù)這一特征，利用變分模態(tài)分解原理對超聲檢測信號和噪聲分解得到的各模態(tài)頻域特性的差異，即超聲信號被分解為信號頻帶內(nèi)一系列窄帶信號，而高斯噪聲被分解到[-π,π]整個頻帶內(nèi)的一系列寬帶信號，以此為依據(jù)提出了基于變分模態(tài)分解的超聲信號降噪方法。該方法通過變分模態(tài)分解，計算各模態(tài)帶寬和中心頻率，選擇信號模態(tài)，由信號模態(tài)重構(gòu)無噪信號4步實(shí)現(xiàn)對超聲信號的降噪。此外，為提高超聲仿真信號和實(shí)測超聲信號的匹配度，本文提出了一種改進(jìn)的高斯超聲信號模型。

1 變分模態(tài)分解原理

變分模態(tài)分解以每個本征模態(tài)函數(shù)是其中心頻率附近的窄帶信號為準(zhǔn)則構(gòu)造變分模型，通過交替方向乘子算法（alternate direction method of multipliers,ADMM）[12]在其頻域進(jìn)行迭代求解，進(jìn)而將輸入信號分解為一組本征模態(tài)函數(shù)的線性組合。

1.1 變分模型構(gòu)造

在變分模型中本征模態(tài)函數(shù)（intrinsic mode functions，IMF）定義為調(diào)幅-調(diào)頻信號，其表達(dá)式為：

其中瞬時幅值A(chǔ)k(t)為非負(fù)函數(shù)，相位?k(t)為單調(diào)非減函數(shù)。

瞬時頻率為：

瞬時頻率和瞬時幅值相對于相位應(yīng)是緩慢變化的，即在足夠長的時間間隔[t-δ,t+δ]（δ≈2π/ωk(t)）上，本征模態(tài)函數(shù)可以看作幅值為Ak(t)，頻率為ωk(t)的余弦波。

變分模態(tài)分解的目標(biāo)是將輸入的實(shí)信號f分解為一組本征模態(tài)函數(shù)uk(t)的線性組合。分解過程中要求本征模態(tài)函數(shù)能表征原信號的同時應(yīng)具備足夠的稀疏性。通過以下3步來構(gòu)造變分模型：

1）為得到信號的單邊頻譜，對每個模式uk(t)進(jìn)行希爾伯特變換并求其對應(yīng)的解析信號：

2）估計每個模式的中心頻率ωk(t)，然后通過混頻將其解調(diào)到基帶：

3）利用H1高斯平滑度估計解調(diào)信號的帶寬[1]。

因此，基于以上3點(diǎn)可得IMF的變分約束模型為：

1.2 變分模型求解

為求解上述約束變分模型，通過引入二次懲罰項(xiàng)和拉格朗日乘子使構(gòu)造的約束變分模型變?yōu)闊o約束變分模型。

通過交替方向乘子算法在頻域求解該方程[1]，其主要步驟如下：

輸出：分解得到的K個模式。

2） 令n=n+1，并根據(jù)下式分別更新每個模式及其中心頻率，其中k=1，2，???，K。

4）根據(jù)下式判斷是否滿足收斂條件，如果收斂則輸出分解得到的K個模式，否則返回步驟2）繼續(xù)迭代分解：

以上就是變分模態(tài)分解的基本過程。需要特別說明的是上述過程是在頻域進(jìn)行求解，其中=FT(x(t))表示對x(t)做傅里葉變換。

2 VMD在超聲信號降噪中的應(yīng)用

2.1 超聲信號模型

在超聲脈沖回波檢測中超聲回波信號通常是一個被超聲探頭中心頻率調(diào)制的寬帶信號[13]，其表達(dá)式為

其中A為反射回波幅度；T為帶寬因子；τ為回波到達(dá)時間；fc為發(fā)射超聲脈沖的中心頻率；φ為初相位。

該模型一經(jīng)提出就得到廣泛應(yīng)用，其在超聲信號數(shù)值分析和仿真方面有著重要作用[14-15]。但是該模型存在兩點(diǎn)不足之處：1）根據(jù)參數(shù)的物理意義，在t＞τ的范圍內(nèi)該模型才有波形，而在t＜τ范圍內(nèi)該模型無波形或波形無意義；2）該模型僅能模擬時域半邊帶的高斯包絡(luò)，這與實(shí)測中的超聲信號有很大差別。為提高在信號起始時刻仿真超聲信號和實(shí)測超聲信號的相似度，本文提出一種改進(jìn)的超聲信號模型，其表達(dá)式為：

其中，a為包絡(luò)函數(shù)B(t)的起點(diǎn)，U(t)為單位階躍函數(shù)，δ(t)為Dirac函數(shù)，*表示卷積運(yùn)算，其他參數(shù)物理意義同式（11）。

2.2 VMD降噪原理

由于超聲信號具有特定的結(jié)構(gòu)，而高斯噪聲沒有固定的結(jié)構(gòu)；表現(xiàn)在頻域超聲信號是一有限帶寬信號，而高斯噪聲對應(yīng)的頻譜分布在數(shù)字頻率ω∈[-π,π]整個頻帶內(nèi)。如圖1和圖2所示分別為采用VMD對超聲信號和高斯噪聲進(jìn)行6層分解得到的頻譜圖。

圖1 無噪信號VMD分解譜

圖2 噪聲VMD分解譜圖

圖1中s為無噪超聲信號的頻譜圖，u1~u6分別為VMD分解后的6個模態(tài)的頻譜圖，從圖中可以看出VMD將原信號s的頻譜在頻域內(nèi)分解為6個中心頻率不同的窄帶信號，每個模態(tài)對應(yīng)于原始信號s頻譜的一部分。

圖2中n為噪聲的頻譜圖，u1~u6為對噪聲進(jìn)行VMD分解后得到的各模態(tài)的頻譜圖。從圖中可以看出VMD將噪聲分解為頻帶內(nèi)不同中心頻率的寬帶信號。對比圖1和圖2可得出結(jié)論：對于K層VMD，信號和噪聲都被等間隔地分解為各自頻帶內(nèi)K個模式。由于超聲信號帶寬相比于高斯噪聲帶寬較窄，對于相同的分解層數(shù)K，VMD分解后噪聲模態(tài)對應(yīng)的帶寬遠(yuǎn)大于信號模態(tài)的帶寬。根據(jù)這一原則可在頻域分辨超聲信號模式和噪聲模式，基于VMD的超聲信號降噪方法具體步驟如下：

1）初始化，設(shè)定待分解層數(shù)K。

2）對給定的含噪超聲信號s(t)，求其傅里葉變換，并令=FT(s(t))，根據(jù)1.2中變分模態(tài)分解過程對其進(jìn)行VMD分解。

3） 計算每個模式的帶寬和中心頻率，依據(jù)信號對應(yīng)窄帶模式，噪聲對應(yīng)寬帶模式原則，對分解后的模式進(jìn)行分類，選擇具有窄帶特性的模式。

4） 根據(jù)選擇的窄帶模式重構(gòu)超聲信號。

3 實(shí)驗(yàn)與仿真分析

當(dāng)前廣泛用于信號降噪結(jié)果評價的指標(biāo)有均方誤差MSE、波形相似系數(shù)NCC和重構(gòu)信噪比ESNR[16]。均方誤差越小，重構(gòu)信噪比越大，波形相似系數(shù)越接近于1表明降噪效果越好。

3.1 超聲信號模型

當(dāng)本文提出的超聲信號模型參數(shù)[A,T,fc,φ]=[0.8, 100, 0.1, 0]時，參數(shù)a取不同值時對應(yīng)的超聲仿真信號如圖3所示。

圖3 超聲仿真信號

當(dāng)a=0時，本文提出的超聲信號模型退化為高斯超聲信號模型。對比圖3中各波形可以看出，改進(jìn)的超聲信號模型相比高斯模型包含的超聲信號包絡(luò)信息更豐富，隨參數(shù)a變化可以模擬起始時刻不同包絡(luò)的超聲信號；而高斯模型起始時刻包絡(luò)是固定不變的。因此，改進(jìn)的超聲信號模型是高斯模型的拓展，比高斯模型包含的包絡(luò)信息更豐富。

3.2 不同模型比較

為驗(yàn)證本文提出的超聲信號模型性能，將本文模型、指數(shù)模型[15]分別和在實(shí)驗(yàn)室中用RITEC RAM-5000-SNAP超聲檢測系統(tǒng)采集的超聲回波信號[16]進(jìn)行匹配，其結(jié)果如圖4所示。經(jīng)計算，指數(shù)模型和實(shí)測超聲信號的均方誤差為0.022 3，波形相似度為0.916 0；本文改進(jìn)的高斯模型和實(shí)測超聲信號的均方誤差為0.008 3，波形相似度為0.973 3。因此，無論直觀對比效果還是客觀評價指標(biāo)，本文提出的超聲信號模型都優(yōu)于指數(shù)信號模型。

圖4 不同模型效果對比

3.3 降噪性能分析

當(dāng)K=5時，對含噪超聲信號（SNR=10 dB）進(jìn)行VMD分解后得到的各模態(tài)頻譜圖，如圖5所示。

圖5 含噪信號分解頻譜

從圖中可以看出，u1~u3模式的頻譜帶寬窄且幅值大，由前述分析知這是信號模式的典型特征；而u4、u5中頻譜帶寬較寬且幅值相對小，為噪聲模式的特征；由u1~u3模式對應(yīng)的信號在時域相加后得到的重構(gòu)超聲信號如圖6所示。

經(jīng)計算，相比原始無噪超聲信號，降噪后信號的均方誤差MSE=0.000 7，波形相似系數(shù)NCC=0.994 0，重構(gòu)信噪比ESNR=19.237 4 dB。

圖6 降噪效果對比

3.4 不同方法性能比較

為驗(yàn)證本文方法的降噪性能，分別采用本文提出的VMD降噪方法、小波降噪方法及EMD降噪方法對圖7（a）中無噪超聲信號加入噪聲后進(jìn)行降噪處理，當(dāng)信噪比為0 dB時降噪結(jié)果如圖7所示。

圖7 不同方法比較

不同方法對應(yīng)的評價指標(biāo)如表1所示。

表1 降噪效果評價指標(biāo)

由圖7和表1知，在對如圖7(a)中所示的超聲檢測信號進(jìn)行降噪處理時，本文提出的方法效果要好于小波方法和EMD方法。

4 結(jié)束語

本文將變分模態(tài)分解理論應(yīng)用到超聲信號處理中，提出了一種基于變分模態(tài)分解的超聲信號降噪方法。通過變分模態(tài)分解，計算各模態(tài)帶寬和中心頻率，選擇信號模態(tài)，根據(jù)信號模態(tài)重構(gòu)無噪信號4步實(shí)現(xiàn)對超聲檢測信號的降噪。此外，為提高超聲仿真信號和實(shí)測超聲信號的匹配度還提出了一種改進(jìn)的超聲信號仿真模型，相比指數(shù)模型和高斯模型，本文提出的仿真模型與實(shí)測超聲信號更匹配。