基于EMD的自適應(yīng)ECG信號(hào)基線漂移消除方法研究

郭樹(shù)言，王黎明，牛曉東，張巖，李集照，趙越

中北大學(xué)信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室，山西太原030051

前言

醫(yī)學(xué)領(lǐng)域在采集心電（ECG）信號(hào)時(shí)，往往會(huì)由于各種原因疊加不同類型的噪聲，如高頻電噪聲、肌電噪聲、由于病人呼吸和活動(dòng)產(chǎn)生的基線漂移（BW）噪聲，這些噪聲的存在有時(shí)會(huì)嚴(yán)重干擾后續(xù)對(duì)ECG信號(hào)QRS 復(fù)合波的檢測(cè)，而QRS 復(fù)合波形狀的一些具體指標(biāo)是醫(yī)生診斷某些疾病的關(guān)鍵，因而在QRS復(fù)合波檢測(cè)前對(duì)采樣得到的原始ECG 信號(hào)進(jìn)行濾波是必不可少的［1-3］。然而，ECG 信號(hào)是一種典型的非線性非平穩(wěn)的信號(hào)，使用傳統(tǒng)時(shí)頻域信號(hào)處理方法（如短時(shí)傅里葉變換、小波變換）都有頻譜展寬的問(wèn)題，且原始ECG 信號(hào)疊加的各種噪聲的頻域范圍很寬，這樣就導(dǎo)致傳統(tǒng)時(shí)頻域信號(hào)處理方法在消除ECG信號(hào)各種噪聲時(shí)失真嚴(yán)重或效果不理想。

近年來(lái)，EMD 方法由于對(duì)非線性非平穩(wěn)信號(hào)良好的處理效果逐漸被國(guó)內(nèi)外學(xué)者應(yīng)用在ECG 信號(hào)去噪領(lǐng)域，尤其是在BW 去除方面提出了許多優(yōu)秀的方法，如EMD 方法與低通濾波器相結(jié)合［4］，EMD 方法與形態(tài)學(xué)濾波相結(jié)合［5］，EMD 方法與BW 特點(diǎn)相結(jié)合進(jìn)行經(jīng)驗(yàn)濾波［1，6］等。上述方法都是通過(guò)EMD 分解產(chǎn)生ⅠMF 分量，再對(duì)單獨(dú)的ⅠMF 分量進(jìn)行處理或檢驗(yàn)。在Weng 等［4］研究中，通過(guò)對(duì)高階ⅠMF 分量進(jìn)行低通濾波，得到對(duì)應(yīng)ⅠMF的低頻成分作為BW 的部分信號(hào)。但低通濾波器截止頻率的選擇沒(méi)有自適應(yīng)性，而且ⅠMF 分量本質(zhì)上也是非線性非平穩(wěn)信號(hào)，直接進(jìn)行頻域低通濾波而不考慮EMD方法中最重要的瞬時(shí)頻率的觀點(diǎn)，必然造成對(duì)單一的ⅠMF 分量，在有的時(shí)間范圍內(nèi)對(duì)信號(hào)過(guò)處理，而在有的時(shí)間范圍內(nèi)濾波不徹底的問(wèn)題；而文獻(xiàn)［5］提出的方法中為判斷每個(gè)ⅠMF 分量是否是零均值的，引入固定的閾值，從而判斷ⅠMF 分量中是否含有BW 噪聲的成分；這樣，由于固定閾值造成算法處理數(shù)據(jù)時(shí)自適應(yīng)性的缺失不可避免。文獻(xiàn)［6］采用與文獻(xiàn)［5］相同的BW 去除算法，自適應(yīng)性不強(qiáng)。本文在詳細(xì)研究EMD 分解產(chǎn)生的ⅠMF 分量物理意義的基礎(chǔ)上，提出一種自適應(yīng)的BW 去除算法，利用EMD 自適應(yīng)性和ⅠMF 頻率篩選性質(zhì)篩選出BW 噪聲，并對(duì)算法性能進(jìn)行定性和定量分析。

1 經(jīng)驗(yàn)?zāi)Ｊ椒纸?/h2>

1.1 EMD方法的物理意義

EMD 方法是1998年黃鍔教授等人提出的，它本來(lái)是為了將非平穩(wěn)非線性信號(hào)分解為一系列的固有模式函數(shù)（ⅠMF）。這樣對(duì)ⅠMF 分量進(jìn)行希爾伯特變換之后得到的瞬時(shí)頻率才是有意義的。因?yàn)檎J(rèn)為每個(gè)ⅠMF分量在任意時(shí)刻只有1個(gè)瞬時(shí)頻率，而原始信號(hào)是所有ⅠMF 分量的疊加（將剩余分量算為1 個(gè)ⅠMF），在原始信號(hào)的任意時(shí)間點(diǎn)上有多個(gè)瞬時(shí)頻率，所以對(duì)原始信號(hào)直接求希爾伯特變換得到的瞬時(shí)頻率是沒(méi)有物理意義的［7-14］。從這一意義上說(shuō)，如果將原始信號(hào)理解為多個(gè)信號(hào)源產(chǎn)生信號(hào)的疊加，那每一個(gè)ⅠMF分量就代表了1個(gè)信號(hào)源產(chǎn)生的信號(hào)，雖然每一個(gè)ⅠMF 分量也是非線性非平穩(wěn)的，但其瞬時(shí)頻率和幅度的變化一定有一個(gè)過(guò)程，即瞬時(shí)頻率和幅度的變化是連續(xù)的，而不是突變的。這樣，每一個(gè)ⅠMF分量的瞬時(shí)頻率總在一個(gè)范圍中，而高階ⅠMF的瞬時(shí)頻率所處的范圍總體上總是低于比它低階的ⅠMF 分量的。如果組成一個(gè)信號(hào)的信號(hào)源發(fā)射信號(hào)的頻率總體上的高低已知，我們可以通過(guò)ⅠMF 分量的頻率范圍來(lái)確定這個(gè)ⅠMF對(duì)應(yīng)的信號(hào)源。

1.2 EMD過(guò)程描述

EMD 方法的核心概念是時(shí)間尺度，它被定義為信號(hào)的一個(gè)極值點(diǎn)到與它相鄰的極值點(diǎn)的時(shí)間范圍（如果是數(shù)字信號(hào)，定義為數(shù)字點(diǎn)數(shù)）。如果是原始信號(hào)，這個(gè)時(shí)間尺度就被視為信號(hào)中最高頻分量（第一個(gè)ⅠMF 分量）的頻率指標(biāo)，而其他ⅠMF 分量合起來(lái)被視為是第一個(gè)ⅠMF 分量的趨勢(shì)，因?yàn)棰馦F 分量是零均值的，EMD 方法中直接將原始信號(hào)上下包絡(luò)求均值作為趨勢(shì)量，然后原始信號(hào)減去趨勢(shì)量得到第一個(gè)ⅠMF 分量［7，12，15-18］。因?yàn)榈玫降蘑馦F 分量一般不滿足零均值條件，需要將上面的過(guò)程迭代多次，具體過(guò)程如下：

（1）對(duì)一個(gè)信號(hào)x(t)，先找出他的所有極大值點(diǎn)和極小值點(diǎn)，然后利用三次樣條插值函數(shù)求得信號(hào)極大值和極小值的包絡(luò)，并求兩個(gè)包絡(luò)的均值，得到一條均值曲線m1(t)，令h1(t)=x(t)-m1(t)，得到h1(t)；

（2）對(duì)得到的h1(t)進(jìn)行步驟（1）操作，得到h2(t)，并重復(fù)10 次，得到第一個(gè)ⅠMF 分量e1(t)，并令r1=x(t)-e1(t)，r1為對(duì)應(yīng)于第一個(gè)ⅠMF分量的剩余分量；

（3）對(duì)r1(t)進(jìn)行步驟（1）、（2），得到e2(t)，e3(t),…,em(t)和r(t)，r(t)為剩余分量；

結(jié)果，EMD將信號(hào)x(t)分解成如下形式：

其中，k為ⅠMF 分量的階數(shù)，而m為總的ⅠMF 分量個(gè)數(shù)，r(t)為剩余分量［19-20］。

2 基于EMD的消除ECG信號(hào)BW的算法

2.1 ⅠMF分量頻率指標(biāo)的確定

由于ECG 信號(hào)的BW 噪聲主要由患者的呼吸或運(yùn)動(dòng)引起，相對(duì)ECG 信號(hào)來(lái)說(shuō)是一種低頻噪聲，所以由以上理論可知，通過(guò)EMD 分解，低頻的BW 噪聲被分解到高階ⅠMF分量中。

前面提到每一個(gè)ⅠMF 分量都有一個(gè)連續(xù)的頻率范圍，而且高階的ⅠMF分量頻率范圍低于低階的ⅠMF分量（可能會(huì)有重疊）。如果將每一個(gè)ⅠMF 分量相鄰兩個(gè)極大值點(diǎn)或極小值點(diǎn)中間的時(shí)間尺度視為其局部的頻率指標(biāo)，定義為Sk，k為第k個(gè)時(shí)間尺度，求出這些時(shí)間尺度的平均值，定義為每一個(gè)ⅠMF 分量的“周期”，并將其作為ⅠMF 分量的頻率指標(biāo)，用Tn來(lái)表示，n為ⅠMF階數(shù)。有：

其中N為第n個(gè)ⅠMF分量中Sk的個(gè)數(shù)。

這里必須指出，理論上高階ⅠMF 分量的Tn小于低階ⅠMF 分量，但有必要進(jìn)行驗(yàn)證，分別使用GSTA（一組表征年平均全球表面溫度異常的數(shù)據(jù)）［8］和MⅠT-BⅠH 心率失常數(shù)據(jù)庫(kù)的100 信號(hào)（MLⅠⅠ導(dǎo)聯(lián)，取3 600 個(gè)數(shù)據(jù)點(diǎn)）進(jìn)行驗(yàn)證。由表1 可以看出無(wú)論GSTA數(shù)據(jù)還是MⅠT-BⅠH 100數(shù)據(jù)的ⅠMF的“周期”都是隨著ⅠMF 分量階數(shù)的升高而增大的，說(shuō)明這里定義的“周期”在實(shí)際使用中與理論一致，可以反映ⅠMF分量的頻率指標(biāo)。

2.2 基于EMD的消除ECG信號(hào)BW的算法

基于第2.1 節(jié)給出的ⅠMF 分量的“周期”Tn，得到一種自適應(yīng)的ECG 信號(hào)的BW 去除算法。假設(shè)一種情況，1 個(gè)ECG 信號(hào)中只有1 個(gè)心率周期的數(shù)據(jù)，通過(guò)EMD 分解之后最高階的ⅠMF 分量的“周期”理應(yīng)是1個(gè)心率周期的長(zhǎng)度，拓展到1個(gè)ECG 信號(hào)中有多個(gè)心率周期的數(shù)據(jù)，當(dāng)通過(guò)EMD 分解之后產(chǎn)生的ⅠMF 分量中，有理由認(rèn)為包含ECG 信號(hào)的ⅠMF 分量的“周期”不大于心率周期，即EMD 分解產(chǎn)生的“周期”大于ECG 信號(hào)周期的ⅠMF 分量都認(rèn)為是BW 噪聲的分解。

這里值得一提的是，通過(guò)對(duì)理想ECG 信號(hào)模型進(jìn)行EMD 分解，發(fā)現(xiàn)ECG 本身存在直流分量，即理想ECG 信號(hào)模型最后的剩余分量（最后一個(gè)ⅠMF 分量）（圖1）。我們有理由認(rèn)為，醫(yī)學(xué)上測(cè)量的ECG 信號(hào)存在直流分量，同時(shí)通過(guò)EMD 分解，ECG 信號(hào)的直流分量被分解到最后一個(gè)ⅠMF 分量中。所以在基于EMD 的ECG 信號(hào)BW 去除算法中，直接將最后一個(gè)ⅠMF 分量去除是不合適的。圖1 中ⅠMF4 和ⅠMF5中明顯存在由于端點(diǎn)效應(yīng)造成的擾動(dòng)，在ⅠMF6 中邊界效應(yīng)造成的擾動(dòng)向內(nèi)移動(dòng)，在之后的ⅠMF 分量都是接近零的一條直線，直到最后的ⅠMF12，明顯偏離零值，在0.12 附近，可以視為模擬標(biāo)準(zhǔn)ECG 信號(hào)的“直流部分”。

鑒于上述分析，我們不能直接將最后一個(gè)ⅠMF分量舍棄，由于ECG信號(hào)是類周期信號(hào)，將一段ECG信號(hào)分成一個(gè)一個(gè)的周期，每個(gè)周期單獨(dú)進(jìn)行EMD分解，每個(gè)周期分解得到的最后一個(gè)ⅠMF 分量是幾乎一樣的，這樣理論上一段ECG 信號(hào)經(jīng)過(guò)EMD 分解得到的最后一個(gè)ⅠMF 分量應(yīng)該近似是一條直線。而實(shí)際上，EMD 分解ECG 信號(hào)得到的最后一個(gè)ⅠMF 分量本身就近似一條直線，同時(shí)由于BW 來(lái)源于患者呼吸或在測(cè)量信號(hào)過(guò)程中的一些動(dòng)作，EMD 分解后不可能存在近乎直流的分量，所以EMD 分解ECG 信號(hào)的最后一個(gè)ⅠMF分量就是來(lái)源于ECG信號(hào)本身的。

于是，本文基于EMD 的消除ECG 信號(hào)BW 的算法具體步驟如下：

（1）使用EMD 方法對(duì)信號(hào)x(t)進(jìn)行處理，得到n個(gè)ⅠMF分量；

（2）計(jì)算每個(gè)ⅠMF 分量的“周期”Tn，如式（2）所示；

（3）計(jì)算x(t)的平均R-R 周期，作為x(t)的心率周期t；

（4）將大于t的Tn對(duì)應(yīng)的ⅠMF 分量去除，但保留最后一個(gè)ⅠMF 分量，然后將剩余的ⅠMF 分量疊加，得到去除基線漂移噪聲的心電信號(hào)xdelBW(t)。

3 仿真結(jié)果分析與應(yīng)用

3.1 定性分析

為直觀看到本文算法的濾波效果，選用美國(guó)麻省理工學(xué)院MⅠT-BⅠH 心率異常數(shù)據(jù)庫(kù)中的100 信號(hào)和103 信號(hào)（都采用MLⅠⅠ導(dǎo)聯(lián)，取65 536 個(gè)數(shù)據(jù)點(diǎn)），用本文算法進(jìn)行處理效果如圖2 和圖3 所示。為方便觀察濾波效果，這里沒(méi)有將最后一個(gè)ⅠMF 分量加入濾波后的信號(hào)，因?yàn)樽詈蟮蘑馦F 分量是原始ECG信號(hào)的直流成分，去除之后只是將濾波后的信號(hào)向上平移，與原始信號(hào)交錯(cuò)開(kāi)，方便觀察。

由圖2 和圖3 可以看出經(jīng)過(guò)本算法處理的ECG信號(hào)明顯平穩(wěn)了許多，那些明顯的非ECG 的趨勢(shì)信號(hào)被濾除，而ECG信號(hào)的細(xì)節(jié)得以保留。

3.2 定量分析

為進(jìn)一步驗(yàn)證本文算法的性能，需要進(jìn)一步的定量分析。因?yàn)镸ⅠT-BⅠH 心率異常數(shù)據(jù)庫(kù)中的數(shù)據(jù)本身含有的BW 噪聲未知，許多研究都會(huì)對(duì)數(shù)據(jù)庫(kù)中的數(shù)據(jù)進(jìn)行一些處理，以得到近似沒(méi)有BW 噪聲的信號(hào)，但這些處理顯然無(wú)法完全得到干凈的信號(hào)，致使后續(xù)的實(shí)驗(yàn)和研究不能讓人完全信服。所以本文采用ECGSYN［9］生成的模擬理想ECG 信號(hào)x(t)（圖4a），其中加入的BW為：bw= 0.1sin(( 2π/4096)t)），并加入低頻的正弦波作為模擬的BW 噪聲bw(t)，得到帶有BW 噪聲的ECG 信號(hào)X(t（)圖4b），這樣就提前已知了信號(hào)和噪聲的全部信息。

將合成的信號(hào)用EMD方法和本文算法進(jìn)行處理后，得到信號(hào)X′(t)和噪聲信號(hào)bw′(t)。通過(guò)x(t)與X′(t)相關(guān)系數(shù)，bw(t)與bw′(t)的相關(guān)系數(shù)，基線矯正率3 個(gè)指標(biāo)來(lái)判斷本文去除BW 算法的性能，并與文獻(xiàn)［1］中的相應(yīng)算法進(jìn)行比較。其中，相關(guān)系數(shù)越高，說(shuō)明信號(hào)還原越精確；基線矯正率越低，濾出的BW越完整。

圖1 EMD分解模擬標(biāo)準(zhǔn)ECG信號(hào)Fig.1 Decomposition of analog standard electrocardiogram(ECG)signals by empirical mode decomposition(EMD)

圖2 MIT-BIH 100原始數(shù)據(jù)及去除基線漂移后的數(shù)據(jù)Fig.2 MIT-BIH 100 raw data and the data after baseline wander removal

相關(guān)系數(shù)：

基線矯正率：

其中，分別用ρnn′表示n(t)與n′(t)的相關(guān)系數(shù)，用ρxX′表示x(t)與X′(t)相關(guān)系數(shù)，BR表示BW率。

由表2 可知，本文算法分離出的ECG 信號(hào)X′(t)和BW 信號(hào)n′(t)與原始的理想ECG 信號(hào)x(t)和BW信號(hào)n(t)幾乎完全一致，相較于文獻(xiàn)［1］的算法有更佳的濾除BW的效果。

4 結(jié)束語(yǔ)

本文在前人研究基礎(chǔ)上，針對(duì)基于EMD 方法濾除ECG 信號(hào)BW 算法存在的不足，提出一種更加貼合EMD方法物理意義的自適應(yīng)的濾除ECG信號(hào)BW的算法。在提出ⅠMF 分量平均周期的概念，以及闡述它與ECG 信號(hào)關(guān)系的基礎(chǔ)上，通過(guò)排除不包含ECG 信號(hào)信息的ⅠMF 分量來(lái)濾除ECG 信號(hào)中的BW噪聲。本文首先使用ECG 數(shù)據(jù)驗(yàn)證了提出的ⅠMF 分量平均周期概念的實(shí)用性，之后通過(guò)本算法處理大量ECG 數(shù)據(jù)直觀體現(xiàn)濾除BW 噪聲的效果，最后，用理想ECG 信號(hào)疊加正弦信號(hào)合成的信號(hào)進(jìn)行定量分析。經(jīng)實(shí)驗(yàn)驗(yàn)證，本文提出的方法由于貼合EMD 的物理意義和其自適應(yīng)性，相較于其他基于EMD 算法的ECG 信號(hào)BW 去除算法有很大的優(yōu)勢(shì)，不論是信號(hào)的相關(guān)系數(shù)還是BW 噪聲的相關(guān)系數(shù)都有所提高。不過(guò)，EMD的端點(diǎn)效應(yīng)的影響依然較大，是算法得到信號(hào)的誤差的主要來(lái)源。

圖3 MIT-BIH 103原始數(shù)據(jù)及去除基線漂移后的數(shù)據(jù)Fig.3 MIT-BIH 103 raw data and the data after baseline drift removal

圖4 理想心電及疊加的模擬基線漂移信號(hào)Fig.4 Ideal ECG and superimposed analog baseline drift signal

表2 文獻(xiàn)［1］算法與本文算法性能指標(biāo)的比較Tab.2 Comparison of performance indexes between the algorithm in literature［1］and the proposed algorithm