基于VMD與包絡(luò)峰值提取的R波檢測算法

龐 宇，汪立宇，陳亞軍

(重慶郵電大學 光電工程學院，重慶 400065)

0 引 言

心電(electrocardiogram，ECG)信號的波形是診斷心血管疾病的重要參考。傳統(tǒng)上心電圖的分析是由QRS波群開始的，它的能量高、振幅大，因此最容易被檢測。而QRS波群的定位檢測的第一步就是先檢測R峰的位置，完成該處理后才能根據(jù)規(guī)律檢測Q波和S波[1]。

對于ECG信號的R波檢測，有差分閾值[2]、小波變換[3]等算法，均具有良好的檢測定位效果。但是，傳統(tǒng)的差分閾值法檢測R波原理簡單，對低頻干擾抑制作用好，但容易受高頻干擾影響從而降低識別率；小波變換檢測計算量略大且對硬件平臺要求高。此外還有香農(nóng)能量變換法[4]、粒子群優(yōu)化算法[5]、神經(jīng)網(wǎng)絡(luò)[6]的模式識別算法等，各種方法都在特定的檢測環(huán)境下表現(xiàn)出不同的優(yōu)勢。經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)及相關(guān)理論的提出，為非平穩(wěn)信號的處理找到了一條新途徑，也適用于ECG的波形檢測[7]，但算法本身存在模態(tài)混疊的缺點。Dragomiretskiy K等提出了變分模態(tài)分解算法(variational mode decomposition，VMD)[8]，該算法可有效消除EMD存在的模態(tài)混疊問題，在生物體征信號的分析處理領(lǐng)域具有很大的應用價值。

VMD是在EMD基礎(chǔ)上，結(jié)合數(shù)學變分法的理論提出的一種信號分析方法，適用于非平穩(wěn)信號的分析，在心電信號R波檢測中具有較強的抗干擾能力與較高的識別率[9]。本文通過對預設(shè)尺度K值選取設(shè)定最優(yōu)化算法，將ECG信號分解提取出特征化模態(tài)分量，然后對該分量進行包絡(luò)峰值提取且設(shè)定復檢規(guī)則，能有效排除背景噪聲的影響，提高R波定位的準確率。

1 R波檢測定位算法

1.1 VMD處理

VMD是一種處理非線性、非平穩(wěn)信號的有效方法，可以將原始ECG信號分解成有限數(shù)量的二維本征模態(tài)函數(shù)(bi-dimensional intrinsic mode function，BIMF)[10]，當再現(xiàn)輸入信號時，經(jīng)分解得到的子信號具有特定的稀疏特性[11]。

VMD將EMD方法中本征模態(tài)函數(shù)定義為窄帶BIMF，即

xk(t)=Ak(t)cos[φk(t)]

(1)

其中，φk(t)非遞減，φ′k(t)≥0；包絡(luò)線非負，即Ak(t)≥0；且Ak(t)與瞬時頻率ωk(t)=φ′k(t)相對于相位φk(t)的變換是緩變的。

VMD方法將信號放入約束變分模型中求取最優(yōu)解，確定限帶本征模函數(shù)的中心頻率及帶寬，模型表達式為

(2)

(3)

其中，xk(t)、ωk分別表示信號的第k個模態(tài)分量及其對應的中心頻率，{xk}={x1,x2,…,xK}、{ωk}={ω1,ω2,…,ωK}則分別表示K個模態(tài)分量及中心頻率的集合。上述約束優(yōu)化問題引用增廣拉格朗日修正方程求解，表達式如下

(4)

式中：α——懲罰參數(shù)；λ——拉格朗日因子。

根據(jù)文獻[10，11]中的方法，從譜域解得到的模態(tài)及對應的更新中心頻率表示為

(5)

(6)

然后利用交替方向乘子算法[12]求該模型最優(yōu)解，得到的結(jié)果即是信號分解后的K個限帶模態(tài)分量。算法步驟如下：

其中τ為雙上升時間步長，ε為收斂性判別標準的公差。

原始ECG信號包含有PQRST各個波段有效信息，以及工頻、肌電、基線漂移、運動偽跡干擾等各種噪聲，經(jīng)過以上過程即可分解為從低頻區(qū)到高頻區(qū)的各個模態(tài)分量，對于R波的檢測定位只需分析特定模態(tài)分量以排除噪聲與其它無效信息的干擾。

1.2 K值選取

VMD技術(shù)中，K值需要預先設(shè)定，其大小直接影響了信號分解的結(jié)果。文獻[13]指出，當K值最優(yōu)時，BIMF末分量的中心頻率將達到最大值。因此，可以通過觀察中心頻率隨不同K值的變化，確定最優(yōu)預設(shè)尺度。

這里設(shè)定確定K值的算法：初設(shè)一個較大k值，令K=k和K=k-1，若兩次分解得到的BIMF末分量的中心頻率相近，則說明設(shè)定的k值偏大；再取K=k-1和K=k-2，依次循環(huán)，當兩次分解得到的最后一個BIMF分量中心頻率差異較大時，可給出最優(yōu)K值。通過設(shè)定一個閾值a(視信號頻域而定)來判斷中心頻率是否相近。在此通過一個仿真信號進行分析。仿真信號如下

f=0.5cos(2*pi*2*t)+0.25cos(2*pi*24*t)+
0.2cos(2*pi*288*t)+η

(7)

式中：η為幅值標準差0.1的加性高斯白噪聲。不同K值下，對信號f進行分解，各BIMF分量對應中心頻率ω分布見表1。

表1 信號f在不同K值下的中心頻率分布

表1指出，當K=4時，BIMF4的中心頻率與K=5時BIMF5的基本相等。當K=5時，出現(xiàn)過分解，其中BIMF3和BIMF4的中心頻率接近，所以K值取4。該實驗為K值選取算法提供了依據(jù)，使用MIT-BIH數(shù)據(jù)庫中的103心電信號仿真分析，如表2所示，K=4時，中心頻率ω4幾乎達到最大值，考慮到ECG信號90%的能量集中在0.05 Hz～45 Hz且P波(0 Hz～17 Hz)、QRS波群(0 Hz～33 Hz)與T波(0 Hz～9 Hz)等特征波形的頻率分布，所以選擇K=4。

表2 信號103在不同K值下的中心頻率

1.3 R波定位

對于R波的定位，核心思想是：首先以變分模態(tài)分解為基礎(chǔ)，確定預設(shè)尺度K值的自適應準則；K個不同的模態(tài)分量是原始信號在不同頻域(從低頻到高頻)的表達，然后根據(jù)頻域相關(guān)性提取出R波所在的有效分量，進行歸一化并平方處理；最后通過平滑濾波獲得包絡(luò)峰值，標記為邏輯R波所在的索引位置。

邏輯上R波波峰僅僅可以初步縮小R波波峰所在位置的區(qū)域，不能做到準確定位。平滑濾波在峰值點的定位過程中存在時移偏差，需要附加峰值搜尋算法定位真正的R波波峰[14]，該算法需要在上一步歸一化平方后、平滑處理前的第三模態(tài)信號區(qū)域峰值群(每個平滑包絡(luò)峰值點的附近存在幾個很陡峭的尖峰)中搜索最大幅值點，設(shè)定該復檢規(guī)則可確保算法的準確性。可在提取到的平滑包絡(luò)峰值前后特定窗口區(qū)選取N個采樣點，最大的幅值點對應原始信號的索引位置即為真正R波峰值位置。故對應的R波定位表達式為

(8)

所以R波的后續(xù)定位算法如下：

(1)信號有效分量歸一化平方處理；

(2)滑動平均濾波器對{N(BIMF)}2進行包絡(luò)平滑處理；

(3)提取平滑包絡(luò)的峰值得到R波的索引位置，標記為邏輯上的R波；

(4)在{N(BIMF)}2搜索包絡(luò)峰值區(qū)域最大的幅值點，在平滑包絡(luò)峰值的窗口Wt分別取N個采樣點，取最大幅值點；

(5)尋找到的最大幅值點對應的即是R波在該模態(tài)分量中的索引位置，在原始信號中標記R波。

綜上，本文的算法框架如圖1所示。

圖1 R波檢測定位算法框架

2 信號仿真分析

為了驗證所提出算法的有效性，本文首先選取MIT-BIH數(shù)據(jù)庫中103信號做過程模擬。圖2為103信號時域與頻域圖。

圖2 103信號

1.2節(jié)有關(guān)K值選取的實驗給出103信號的最優(yōu)模態(tài)分解數(shù)為4，圖3以及組圖4為VMD分解給出的每個子信號與對應頻譜。

圖3 103信號VMD分解后的模態(tài)分量

圖4 各模態(tài)分量對應頻譜

分解后的BIMF1模態(tài)主要由輸入信號的最低頻率組成，BIMF3由信號的高頻QRS區(qū)組成，T波和P波完全減弱，終模態(tài)BIMF4由心電其余高頻成分或是50 Hz工頻噪聲組成，結(jié)合1.2節(jié)中所提到的ECG信號頻率分布情況，提取BIMF3作為有效分量做R波定位分析。圖5所示為BIMF3經(jīng)過歸一化再平方處理效果圖。

圖5 BIMF3歸一化再平方處理

歸一化平方處理能夠突出信號幅值特征，放大大的值，縮小小的值。接著，將歸一化平方后的BIMF3取平滑包絡(luò)的峰值，定義為邏輯上的R波在模態(tài)分量上的索引位置，標記如圖6所示。

圖6 邏輯R波對應的索引位置標記

用滑動平均濾波器對歸一化平方的BIMF3處理得到的平滑包絡(luò)，其標記的邏輯峰值是將波形的每一處峰值群的幾個尖峰平均處理得到的一個峰值集合，必然存在對于真正R波索引位置的微小時移，每一個標記的邏輯峰值附近存在不等的尖峰如圖7所示。需要根據(jù)R波定位算法的復檢規(guī)則，搜尋真正的最大幅值點標記。

圖7 包絡(luò)峰值的多個實際尖峰

一般人的心率絕對小于300 次/分鐘，即在一定間隔內(nèi)，只能出現(xiàn)一段QRS波，醫(yī)學上稱之為“不應期”。心電圖學中，R峰的不應期在200 ms。MIT-BIH庫中信號的采樣率為360 Hz，以及QRS波群的時間寬度為100 ms，可設(shè)定在初始定位位置前后取50 ms窗口，依次取25個采樣點[15]，標記最大幅值點，對應真正R波的定位位置標記如圖8所示。

圖8 103信號的R波定位

以上仿真實驗復刻了算法的基本流程。信號103屬于數(shù)據(jù)庫中較為純凈的ECG信號，為了檢驗算法的魯棒性，再選取庫中幾個典型的對于R波檢測干擾較大的信號仿真驗證，信號來源于108、109、113、209，分別對應特征為R波倒置、R波倒置加基線漂移、T波尖聳、強噪聲。其定位效果如圖9所示。

圖9 具有不同特征的R波定位效果

3 結(jié)果與討論

計算算法的靈敏度與準確率，需要多樣本參考，本實驗繼續(xù)采用MIT-BIH庫中樣本仿真，主要針對5類典型ECG信號：正常竇性心律(NSR)、左束支傳導阻滯(LBBB)、右束支傳導阻滯(RBBB)、室性早搏(PVC)與房性早搏(APB)。

本文選用靈敏度(SEN)、陽性準確率(+P)以及準確率(Acc)評估算方法指標，引入真陽性(TP)、假陽性(FP)、真陰性(TN)、假陰性(FN)4個參數(shù)，分別代表正確檢測出的R峰數(shù)量、誤將尖峰噪聲視為R峰的數(shù)量、正確判斷為尖峰噪聲數(shù)量、誤將R峰視為尖峰噪聲的數(shù)量，其計算式如下

(9)

(10)

(11)

考慮到真陰性(TN)指標，正確判斷尖峰噪聲數(shù)量的工作難以進行并且不是本文研究重點，所以可忽略該指標影響，然后結(jié)合文獻[1，8，10，14]對評價指標的規(guī)定重新定義式(11)

(12)

樣本信號采用30 min時程數(shù)據(jù)，對應的R波識別統(tǒng)計結(jié)果見表3。

表3 本文算法的R波識別統(tǒng)計結(jié)果

為進一步表明算法優(yōu)勢，將本文算法與參考文獻中出現(xiàn)的算法比對靈敏度、陽性準確率及準確率指標(部分文獻數(shù)據(jù)不足)，得到對照結(jié)果見表4。

4 結(jié)束語

對于原始ECG信號噪聲強R波不便識別的問題，本文通過采用VMD技術(shù)對信號實現(xiàn)各特征頻域的分解，排除噪聲的干擾從而提取到R波所在的有效模態(tài)分量，對該分量進行分析能大大簡化問題思路。接下來通過對分量信號歸一化平方處理進一步濾除干擾，平滑濾波提取包絡(luò)峰值初步確定R波位置，然后根據(jù)窗口采樣方法校正R波索引位置從而使定位更加準確。借助MIT-BIH庫中原始信號的支持，采用多樣本仿真獲取統(tǒng)計結(jié)果驗證了算法的有效性。本文算法在實驗仿真中表現(xiàn)出了高靈敏度高準確率的優(yōu)勢，但相對于傳統(tǒng)算法計算量較大且過程相對繁瑣，適合應用于對定位要求相對較高的領(lǐng)域。

表4 本文算法與相關(guān)文獻算法結(jié)果對比