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

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

0 引 言

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

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

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

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

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

1.1 構造變分模型

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

(1)

對每個解析后的信號與預估中心頻率e-jωKt混合，然后將頻譜調制到頻帶上

(2)

計算以上信號梯度的平方L2范數得到約束條件

(3)

(4)

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

1.2 變分模型求解

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

(5)

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

(6)

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

(7)

(8)

(9)

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

2.1 K-L散度

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

(10)

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

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

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

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

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

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

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

3 樣本熵

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

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

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

(11)

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

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

(12)

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

(13)

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

(14)

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

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

(15)

當N為有限值時

(16)

4 非局部均值去噪

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

圖2 一維NLM原理

(17)

(18)

式中：B是以i點為中心的相似區(qū)域，L為領域B中包含點的個數，λ是濾波器參數。

5 仿真結果

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

圖3 Bumps信號去噪

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

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

表1 去噪效果比較

6 實際心電信號去噪

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

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

表2 重構效果

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

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

圖4 VMD分量

圖5 樣本熵曲線

對噪聲主導分量進行非局部均值去噪(non-local means)，根據文獻中提供的NLM相關參數的確定方法，濾波參數h=0.5*σ，σ為噪聲的標準方差。經過本文的多次實驗，確定了合適的非局部均值去噪參數。NLM參數設定為：搜索窗口t=5，鄰域窗口f=5，濾波參數h=200。將去噪后的分量與有效信號主導分量重構得到去噪信號如圖6所示。

圖6 去噪結果

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

(19)

(20)

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

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

表3 去噪結果

圖7 其它算法去噪結果

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

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

7 結束語

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