優(yōu)化VMD與NLM結(jié)合的信號去噪

尹佳璠，陳小奇，李世林，邵 磊，李 季，劉宏利+

(1.天津理工大學(xué) 電氣電子工程學(xué)院 復(fù)雜系統(tǒng)控制理論及應(yīng)用重點(diǎn)實(shí)驗(yàn)室，天津 300384；2.天津喜諾生物醫(yī)藥有限公司 儀器研發(fā)部，天津 300457)

0 引 言

人體心電信號是判斷人類心臟疾病的重要依據(jù)。且心電信號具有不平穩(wěn)、隨機(jī)性、非線性等特點(diǎn)。在采集過程中又易被噪聲干擾，例如周圍環(huán)境和機(jī)器內(nèi)部電場磁場導(dǎo)致的工頻干擾，因此，需要抑制ECG信號中噪聲干擾。

針對ECG信號微弱、不平穩(wěn)的特性，經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)無疑是最合適且應(yīng)用最廣泛的，但是，經(jīng)驗(yàn)?zāi)B(tài)分解在進(jìn)行信號重構(gòu)時(shí)會存在模態(tài)混疊的缺陷。后續(xù)也有研究人員通過添加白噪聲的方法優(yōu)化經(jīng)驗(yàn)?zāi)B(tài)分解。文獻(xiàn)[1]利用改進(jìn)的小波閾值去噪方法與EEMD結(jié)合，并利用過零率選取合適的固有模態(tài)分量。文獻(xiàn)[2]提出了一種變分模態(tài)分解的參數(shù)選擇的方法，文獻(xiàn)[3]運(yùn)用CEEMDAN與小波閾值結(jié)合的去噪方法對ECG信號進(jìn)行去噪處理。但這些方法的不足之處在于噪聲分量和信號分量分離的不完全以及去噪效果不夠理想等缺點(diǎn)。研究人員提出了一種全新的變分模態(tài)分解(variational mode decomposition，VMD)，解決了模態(tài)混疊的問題。但是，這種分解方法需要通過確定合適的參數(shù)才能達(dá)到較好的分解效果。傳統(tǒng)的參數(shù)選取方法主要通過經(jīng)驗(yàn)選取，參數(shù)選取非常不準(zhǔn)確。

因此，本文設(shè)計(jì)了一種參數(shù)優(yōu)化的變分模態(tài)分解算法，將原始信號分解，引入樣本熵的概念計(jì)算分解后信號分量中含噪聲的比例，將噪聲分量與信號分量區(qū)分開，對噪聲分量進(jìn)行非局部均值去噪(non-local means，NLM)，得到的去噪后的信號分量與低階分量重構(gòu)，得到去噪后的ECG信號。通過對實(shí)際信號進(jìn)行多種方法分解重構(gòu)并與原信號對比，驗(yàn)證了參數(shù)優(yōu)化的變分模態(tài)分解后的效果最好。

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

變分模態(tài)分解(variational mode decomposition，VMD)通過迭代法搜尋變分模型中的最優(yōu)解，確定模態(tài)分量及其中心頻率和帶寬，變分模態(tài)分解本質(zhì)上是基于維納濾波、Hilbert變換以及頻率混合推廣得來的[4]。運(yùn)用這種方法可以將信號分解為K個(gè)模態(tài)分量，每個(gè)分量的中心頻率為ωK，將固有模態(tài)分量的帶寬問題變成具有約束的變分問題。采用這種不同于經(jīng)驗(yàn)?zāi)B(tài)分解的處理方式，不直接從已有的信號中分離，因此有效避免了類EMD算法中不可避免的模態(tài)混疊和端點(diǎn)效應(yīng)。變分模態(tài)分解算法步驟基本可以總結(jié)為構(gòu)造變分模型和模型的求解兩步，理論推導(dǎo)如下。

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

對原始信號x(t)分解得到K個(gè)模態(tài)分量uK(t)，進(jìn)行Hilbert變換，得到解析信號

(1)

對每個(gè)解析后的信號與預(yù)估中心頻率e-jωKt混合，然后將頻譜調(diào)制到頻帶上

(2)

計(jì)算以上信號梯度的平方L2范數(shù)得到約束條件

(3)

(4)

其中：{uK}={u1,…,uK}是分解后得到的K個(gè)IMF分量；{ωK}={ω1,…,ωK}為每個(gè)分量的中心頻率。δ(t)為單位脈沖函數(shù)；j為虛數(shù)單位；*表示卷積運(yùn)算；?K為計(jì)算偏導(dǎo)數(shù)。

1.2 變分模型求解

對于式(3)、式(4)得到的約束條件，引入懲罰因子α和Lagrange算子λ(t)，然后利用交替方向的乘子算法求出增廣函數(shù)的鞍點(diǎn)，然后通過迭代的方式求出各個(gè)參數(shù)的最優(yōu)解，最后逆變換得到各個(gè)模態(tài)分量。懲罰因子可以保證信號收到噪聲干擾時(shí)進(jìn)行分解重構(gòu)后的信號具有較高的精度，Lagrange算子λ(t)可以保證約束的嚴(yán)格性[5]，將兩個(gè)參數(shù)有帶入后，得到

(5)

運(yùn)用交替方向乘子算法(alternate direction method of multiplers, ADMM)，交替更新{uK}，{ωK}，λ

(6)

運(yùn)用傅里葉變換，假設(shè)在已知其它變量的情況下，一次求解一個(gè)變量。因此，在n+1次循環(huán)后可得

(7)

(8)

(9)

2 變分模態(tài)參數(shù)優(yōu)化方法

2.1 K-L散度

使 K-L散度又稱相對熵，可以判斷兩個(gè)概率分布的相似程度，等價(jià)于兩個(gè)概率分布的信息熵差值，公式如下

(10)

式中：p(xi)為真實(shí)數(shù)據(jù)的概率分布，q(xi)為理論概率分布。兩組數(shù)據(jù)越相近則相對熵越小，兩組數(shù)據(jù)分布差別越大，則相對熵越大[6]。

2.2 變分模態(tài)分解的參數(shù)優(yōu)化方法

根據(jù)變分模態(tài)分解方法和理論推導(dǎo)，變分模態(tài)分解算法雖然克服了傳統(tǒng)經(jīng)驗(yàn)?zāi)B(tài)分解及其改進(jìn)方法的缺點(diǎn)，但分解前需要設(shè)定幾個(gè)參數(shù)值才能得到最好的分解效果，其中分解層數(shù)K和懲罰因子α對分解結(jié)果的影響最大。通常情況下，進(jìn)行變分模態(tài)分解前需要人為選取K，α的值，但是人為選取無法使變分模態(tài)分解達(dá)到最好的效果。因此，簡單直接選取參數(shù)是急需解決的問題。文獻(xiàn)[5]提出了一種運(yùn)用遺傳算法優(yōu)化的變分模態(tài)分解參數(shù)。文獻(xiàn)[7]運(yùn)用粒子群算法優(yōu)化變分模態(tài)分解中參數(shù)的選取。文獻(xiàn)[4]提出了改進(jìn)變分模型并通過粒子群算法選取參數(shù)的優(yōu)化變分模態(tài)分解算法。但無論是基于遺傳算法還是粒子群算法的優(yōu)化，除了這兩種算法都容易陷入局部最優(yōu)的缺陷以外，其運(yùn)算所需要的時(shí)間也極長，這和算法本身的十分復(fù)雜有關(guān)。因此，這幾種方法無論是從計(jì)算復(fù)雜度運(yùn)算時(shí)間，還是可操作性上都不夠理想。

為了尋找一種思路簡單又便于運(yùn)行的算法，本文根據(jù)文獻(xiàn)[8]提出的思想，結(jié)合稀疏字典學(xué)習(xí)的思路。改進(jìn)了停止條件，通過固定分解層數(shù)，優(yōu)化懲罰因子，取得最優(yōu)參數(shù)。

2.3 參數(shù)優(yōu)化流程

參數(shù)優(yōu)化流程如圖1所示，具體操作方法主要分為兩步：首先確定懲罰因子α，設(shè)定模態(tài)K的取值范圍，K的取值范圍設(shè)定為[2,20]，初始值為2，對原始信號進(jìn)行變分模態(tài)分解，再通過計(jì)算確定每一次分解所得的本征模態(tài)分量(IMF)的相對熵，獲得相對熵最小值時(shí)所對應(yīng)的K的最優(yōu)解。再根據(jù)得到最優(yōu)的K值，確定懲罰因子α的最優(yōu)解。設(shè)定α的范圍為[100,2000]，步長為50。再尋找最小相對熵所對應(yīng)的α的值，是為最優(yōu)的α值。此時(shí)，得到的K，α值都是最優(yōu)的。變分模態(tài)分解的參數(shù)便確定了下來。這個(gè)方法簡單且易于執(zhí)行，用時(shí)少，運(yùn)行邏輯簡單。在代碼編寫和數(shù)據(jù)處理上更加方便。

圖1 參數(shù)優(yōu)化流程

3 樣本熵

由于經(jīng)過分解后的信號分量包含了信號主導(dǎo)分量和噪聲主導(dǎo)分量，本文通過對信號分量求樣本熵，判斷信號分量的中噪聲含量的大小，樣本熵的計(jì)算原理如下：

(1)對一個(gè)信號序列{x(i),i=1,2,…,N},假定一個(gè)m維向量序列

Xm(i)={x(i),x(i+1),…，x(i+m-1)},1≤i≤N-m+1

(11)

(2)將向量Xm(i),Xm(j)中對應(yīng)向量的絕對差值設(shè)為d[Xm(i),Xm(j)]，即

d[Xm(i),Xm(j)]=
maxk=0,……,m-1(|x(i+k)-x(j+k)|)

(12)

(3)給定一個(gè)參數(shù)r，對于每個(gè)Xm(i)，確定d[Xm(i),Xm(j)]

(13)

(4)定義B(m)(r)為

(14)

(5)令m=m+1，重復(fù)以上述4個(gè)步驟，得到B(m+1)(r);

(6)樣本熵的定義為，在同一參數(shù)r上，有

(15)

當(dāng)N為有限值時(shí)

(16)

4 非局部均值去噪

非局部均值去噪算法(non-local mean)能夠有效去除噪聲并保留信號的真實(shí)性。該算法是利用了圖像的自相似性，即設(shè)置一個(gè)搜索窗口，在這個(gè)搜索窗口內(nèi)尋找盡可能多的相似值，然后運(yùn)用加權(quán)平均值得到真實(shí)的值并同時(shí)去除了噪聲干擾[9]。通過樣本熵判斷出噪聲主導(dǎo)的信號分量并對這些分量進(jìn)行非局部均值(non-local mean)去噪處理得到去噪分量，在去除噪聲干擾的同時(shí)，保證了信號的細(xì)節(jié)。運(yùn)算原理如圖2所示。

圖2 一維NLM原理

(17)

(18)

式中：B是以i點(diǎn)為中心的相似區(qū)域，L為領(lǐng)域B中包含點(diǎn)的個(gè)數(shù)，λ是濾波器參數(shù)。

5 仿真結(jié)果

為了驗(yàn)證本文方法的去噪效果，本文首先對仿真信號進(jìn)行去噪實(shí)驗(yàn)。仿真信號采用MATLAB中自帶的Bumps信號進(jìn)行實(shí)驗(yàn)，采樣點(diǎn)數(shù)為1000。Bumps信號與ECG信號相似，都是非平穩(wěn)非線性的連續(xù)信號，具有代表性，可以驗(yàn)證本文算法的優(yōu)越性。經(jīng)過前文中提出的方法，確定了仿真信號的最優(yōu)分解參數(shù)為：K=7，α=250。向仿真信號中添加信噪比為20 dB的高斯白噪聲，分別對比傳統(tǒng)的小波閾值去噪、文獻(xiàn)[1]提出的EEMD-小波閾值去噪、文獻(xiàn)[3]提出的CEEMDAN小波閾值去噪以及文獻(xiàn)[9]的傳統(tǒng)NLM去噪效果。由于篇幅限制，本文只截取了最終的去噪結(jié)果作為對比，去噪效果如圖3所示。

圖3 Bumps信號去噪

圖3提供了Bumps信號在加入白噪聲信號后的去噪圖像，對比圖3中各個(gè)方法去噪后的結(jié)果能夠發(fā)現(xiàn)，本文提出的方法在信號保真度和去噪效果上都較好，有效抑制了高頻部分的噪聲和毛刺，并且完全還原了峰值部分的原始信號狀態(tài)，基線部分的噪聲也去除了大部分，但還殘留有少量的噪聲。相比之下，圖3(b)中小波閾值去噪效果最差，曲線峰值上的噪聲雖然處理干凈但是也損失了大量的有效信號，基線部分的噪聲去噪去除效果更差，甚至改變了原始信號。圖3(c)中EEMD與小波閾值結(jié)合去噪方法同樣在基線部分的去噪表現(xiàn)不好，有大量的毛刺起伏，基線部分不夠平穩(wěn)和圓滑。圖3(d)中CEEMDAN與小波閾值結(jié)合的去噪方法相比較EEMD-小波閾值去噪效果好了很多，但仔細(xì)觀察會發(fā)現(xiàn)毛刺現(xiàn)象比較嚴(yán)重，而波峰頂部去噪后也有明顯的信號失真，真實(shí)性不如VMD-NLM去噪。圖3(e)中傳統(tǒng)的NLM去噪方法雖然抑制基線部分毛刺比較出色，但是明顯嚴(yán)重?fù)p害了信號的真實(shí)性，幅值出現(xiàn)了較大的改變。因此，其它文獻(xiàn)中改進(jìn)的幾種方法去噪效果都不如本文算法。為了量化去噪效果，使得去噪效果更加的直觀。本文通過計(jì)算信號的信噪比(SNR)和均方誤差(MSE)兩個(gè)指標(biāo)判斷去噪情況。

為了驗(yàn)證去噪效果的準(zhǔn)確性以及合理性，降低偶然性，本文進(jìn)行了重復(fù)性實(shí)驗(yàn)，通過向仿真信號中分別加入信噪比為5 dB、10 dB、15 dB、20 dB的噪聲信號，將各種方法去噪后的信噪比和均方誤差列出成表。根據(jù)表1可以直觀得出，本文算法在Bumps仿真信號的去噪過程中，均能夠得到最高的信噪比數(shù)值和最低的均方誤差數(shù)值。在不同的噪聲干擾下，去噪情況都優(yōu)于其它算法。本文實(shí)驗(yàn)環(huán)境為matlab2017a，win10-64位系統(tǒng)，8 G運(yùn)行內(nèi)存。詳細(xì)結(jié)果見表1。

表1 去噪效果比較

6 實(shí)際心電信號去噪

為了驗(yàn)證本文去噪算法在實(shí)際采集的信號中的效果，驗(yàn)證算法的可行性。本文從麻省理工大學(xué)MIT-BIH心律失常數(shù)據(jù)庫中提取了噪聲較多的112 m實(shí)際ECG信號為實(shí)驗(yàn)信號，該數(shù)據(jù)庫共享了醫(yī)院采集的眾多心電信號，在國際上具有極高的權(quán)威性。為了便于觀察，本文截取了1000個(gè)采樣點(diǎn)，采樣頻率為360 Hz進(jìn)行實(shí)驗(yàn)驗(yàn)證。

為了驗(yàn)證本文提出的優(yōu)化參數(shù)的變分模態(tài)分解算法與其它分解算法相比所具有的優(yōu)越性，對ECG信號分別使用參數(shù)優(yōu)化的變分模態(tài)分解算法、經(jīng)驗(yàn)?zāi)B(tài)分解算法(EMD)、EEMD和CEEMDAN算法分解重構(gòu)后比較，對比重構(gòu)后信號與原始信號的相似程度。根據(jù)前文提出的方法，尋得最優(yōu)參數(shù)為K=17,α=100，由于篇幅限制，只給出了各個(gè)方法重構(gòu)后得信號與原始信號的相關(guān)系數(shù)對比表見表2。

表2 重構(gòu)效果

據(jù)表2提供的重構(gòu)信號和原始信號的相關(guān)系數(shù)，可直觀得知，經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)算法重構(gòu)后的效果相對而言最差，而經(jīng)過改進(jìn)的EEMD和CEEMDAN算法效果逐漸變好，但是經(jīng)過參數(shù)優(yōu)化后的變分模態(tài)分解算法效果最好，本文算法的重構(gòu)信號與原始信號的相關(guān)系數(shù)能達(dá)到100%，各處的細(xì)節(jié)都有較高的還原度。經(jīng)過以上步驟能夠驗(yàn)證經(jīng)過本文優(yōu)化過的分解算法能夠完美還原原始信號。而其它的幾種常用的分解方法雖然效果也不錯(cuò)，但是仍然比本文優(yōu)化的算法稍遜一籌。

然后對含有噪聲的ECG信號進(jìn)行優(yōu)化參數(shù)的變分模態(tài)分解(VMD)去噪，如圖4所示，再對每個(gè)分量求取樣本熵如圖5所示，本文選取的樣本熵參數(shù)為：嵌入維數(shù)m=1，相似度閾值r=0.15*A，A為每個(gè)模態(tài)分量的標(biāo)準(zhǔn)差。根據(jù)樣本熵與信號噪聲的關(guān)系，樣本熵?cái)?shù)值越大，則信號中噪聲含量越多。樣本熵?cái)?shù)值越小，有效信號含量越多的原理。有效信號主要集中在第一、二，十一，十二，十三、十四分量，余下的分量均為噪聲主導(dǎo)的分量。

圖4 VMD分量

圖5 樣本熵曲線

對噪聲主導(dǎo)分量進(jìn)行非局部均值去噪(non-local means)，根據(jù)文獻(xiàn)中提供的NLM相關(guān)參數(shù)的確定方法，濾波參數(shù)h=0.5*σ，σ為噪聲的標(biāo)準(zhǔn)方差。經(jīng)過本文的多次實(shí)驗(yàn)，確定了合適的非局部均值去噪?yún)?shù)。NLM參數(shù)設(shè)定為：搜索窗口t=5，鄰域窗口f=5，濾波參數(shù)h=200。將去噪后的分量與有效信號主導(dǎo)分量重構(gòu)得到去噪信號如圖6所示。

圖6 去噪結(jié)果

為了直觀展示本文算法與第5節(jié)中文獻(xiàn)中提出的算法在去噪效果上的差異。根據(jù)傳統(tǒng)的判斷信號噪聲含量的指標(biāo)主要時(shí)信噪比(SNR)和均方誤差(MSE)

(19)

(20)

其中，si為原信號，fi為去噪信號，N為信號長度。

通過圖7與表3所給出的圖形和數(shù)據(jù)，可以看出本文提出的去噪方法無論在圖像表現(xiàn)上還是量化的數(shù)值比較上都比其它算法更為優(yōu)秀，信號保真度較高，在去除噪聲的同時(shí)保證了有效信號不流失。與本文算法相比，圖7(a)中小波閾值去噪在實(shí)際心電信號上的去噪表現(xiàn)一般，噪聲去除不夠完整，且改變了信號的完整性，曲線也顯得十分扭曲。圖7(b)中EEMD與小波閾值結(jié)合的去噪方法在低頻信號處去噪效果不夠完美，在S波后的第一個(gè)波峰上的處理不好，毛刺較多，這一部分基本沒有去噪，但在波谷的噪聲信號去除的不錯(cuò)曲線比較圓滑。圖7(c)中CEEMDAN與小波閾值結(jié)合的去噪方法改變了信號的幅值，使信號出現(xiàn)了一些震蕩，同時(shí)在Q波波谷的去噪效果一般，保留了些許毛刺噪聲。圖7(d)中非局部均值去噪效果不錯(cuò)，但是信號失真嚴(yán)重，有效信號流失較多，并且出現(xiàn)了漂移現(xiàn)象。因此，通過對比幾種不同的算法在編號102 m心電信號上的去噪效果，根據(jù)表3所示的信噪比與均方誤差的數(shù)值，可以明顯發(fā)現(xiàn)，本文提出的改進(jìn)參數(shù)的變分模態(tài)分解算法，無論是在信噪比還是均方誤差上都要優(yōu)于其余算法。在去除了毛刺和噪聲的同時(shí)，保證了曲線的平滑性而心電信號中最關(guān)鍵的Q、R、S波收到的損害最小，保證了心電信號的真實(shí)性和可靠性。

表3 去噪結(jié)果

圖7 其它算法去噪結(jié)果

為了進(jìn)一步驗(yàn)證本文算法的普適性并且減少偶然性，提高準(zhǔn)確性和可信度，本文從MIT-BIH心率失常數(shù)據(jù)庫中分別選取了100 m、102 m、103 m、109 m、116 m、113 m等6種不同的心電信號進(jìn)行重復(fù)性去噪實(shí)驗(yàn)，分別列出幾種算法的去噪結(jié)果。由于篇幅限制去噪后的信噪比和均方誤差見表4(由于表格長度限制，信噪比精確到小數(shù)點(diǎn)后3位，均方誤差精確到小數(shù)點(diǎn)后4位)。

表4 不同信號去噪效果對比

7 結(jié)束語

本文提出了一種參數(shù)優(yōu)化的變分模態(tài)分解，通過確定一個(gè)變量搜索另一個(gè)變量的方式，利用相對熵作為停止條件，尋找最優(yōu)的參數(shù)。根據(jù)優(yōu)化過的變分模態(tài)分解將實(shí)際ECG信號分解得到多個(gè)模態(tài)分量，再通過計(jì)算各個(gè)模態(tài)分量的樣本熵，確定噪聲主導(dǎo)分量和信號主導(dǎo)分量，對噪聲主導(dǎo)的分量進(jìn)行非局部均值去噪，最后將信號分量和去噪分量重構(gòu)，得到去噪后的信號。對比分析EEMD與小波閾值結(jié)合的去噪方法、CEEMDAN與小波閾值結(jié)合的去噪方法以及普通非局部均值濾波等幾種先進(jìn)算法，本文提出的變分模態(tài)分解和非局部均值結(jié)合的去噪方法能夠自適應(yīng)不同的信號，能得到最好的信噪比和均方誤差，在最大程度上保存信號的真實(shí)性。