龐 宇，汪立宇，陳亞軍

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

0 引 言

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

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

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

1 R波檢測(cè)定位算法

1.1 VMD處理

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

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

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

(1)

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

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

(2)

(3)

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

(4)

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

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

(5)

(6)

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

其中τ為雙上升時(shí)間步長(zhǎng)，ε為收斂性判別標(biāo)準(zhǔn)的公差。

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

1.2 K值選取

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

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

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

(7)

式中：η為幅值標(biāo)準(zhǔn)差0.1的加性高斯白噪聲。不同K值下，對(duì)信號(hào)f進(jìn)行分解，各BIMF分量對(duì)應(yīng)中心頻率ω分布見表1。

表1 信號(hào)f在不同K值下的中心頻率分布

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

表2 信號(hào)103在不同K值下的中心頻率

1.3 R波定位

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

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

(8)

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

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

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

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

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

(5)尋找到的最大幅值點(diǎn)對(duì)應(yīng)的即是R波在該模態(tài)分量中的索引位置，在原始信號(hào)中標(biāo)記R波。

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

圖1 R波檢測(cè)定位算法框架

2 信號(hào)仿真分析

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

圖2 103信號(hào)

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

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

圖4 各模態(tài)分量對(duì)應(yīng)頻譜

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

圖5 BIMF3歸一化再平方處理

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

圖6 邏輯R波對(duì)應(yīng)的索引位置標(biāo)記

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

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

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

圖8 103信號(hào)的R波定位

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

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

3 結(jié)果與討論

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

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

(9)

(10)

(11)

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

(12)

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

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

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

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

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

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