姚晰童，代 煜?，張建勛，葛錦濤，陳 通，楊 灝

1) 南開大學(xué)機(jī)器人與信息自動化研究所，天津 300310 2) 南開大學(xué)電子信息與光學(xué)工程學(xué)院，天津 300310

陡脈沖是指短時高壓電脈沖.通過向細(xì)胞施加陡脈沖可以使得細(xì)胞出現(xiàn)瞬時的電穿孔.基于這種特性，陡脈沖廣泛的應(yīng)用于靶向藥物傳遞[1]、癌癥治療[2]等多個領(lǐng)域，其可行性已經(jīng)在人體實驗中得到了證實[3].除此之外，陡脈沖在材料領(lǐng)域[4]、食品工業(yè)[5]也有著廣泛的應(yīng)用.

心電是心臟活動的重要表現(xiàn)，是醫(yī)生對患者生理狀態(tài)進(jìn)行判斷的重要佐證.然而，心電信號（Electrocardiogram, ECG）受到陡脈沖的干擾，波形可能發(fā)生變異，這導(dǎo)致病患的生理狀況無法被掌握，影響治療效果.以不可逆電消融療法[6]為例，心電信號的作用除監(jiān)測患者生理狀態(tài)外，還包括同步陡脈沖的施放，如果不將陡脈沖的施放時間與心電周期同步，手術(shù)范圍將被嚴(yán)格的限制在遠(yuǎn)離心臟的20 cm的范圍外.

目前，對心電信號的濾波算法的研究，主要針對高頻噪聲、基線漂移、工頻干擾，對于陡脈沖干擾的濾波算法研究較為缺乏.由于陡脈沖產(chǎn)生的相關(guān)噪聲與心電信號主成分頻帶重合，經(jīng)典的數(shù)字濾波器[7]雖然簡單易實現(xiàn)，但是效果較差；小波算法[8]對心電信號濾波是基于信號與噪聲頻譜分離特性，運(yùn)用于此場景效果不理想，且小波基的選擇影響最終結(jié)果；基于EMD的濾波算法[9]改進(jìn)了小波算法的缺陷，但是對分解信號的進(jìn)一步處理需要依賴經(jīng)驗，且對采樣頻率和噪聲敏感.因此，本文針對陡脈沖干擾的特性，基于VMD[10]算法提出了解決方案，實現(xiàn)了對陡脈沖干擾下心電信號的濾波并滿足了心電周期同步需求.

1 陡脈沖干擾下心電信號的分析

陡脈沖具有如下特征：脈沖幅值要求達(dá)到kV·cm-1，脈沖上升沿時間為ns級別，脈沖寬度為μs級別[11].可用表達(dá)式（1）表示：

式中，δpluse(t)表示高壓陡脈沖的時域曲線，S為陡脈沖幅值，Δ為陡脈沖脈寬，h為陡脈沖施放時間點(diǎn).陡脈沖對ECG的干擾如圖1所示.時域上可將陡脈沖干擾分為兩個部分，即：前支和后支.其中，前支頻率高，時間短，位于心電信號中的不應(yīng)期，相當(dāng)于陡脈沖對心電采集設(shè)備等效RC電路的充電過程；后支以前支的終點(diǎn)為起點(diǎn)，頻率低，時間長，對心電信號的波形影響較大，其數(shù)學(xué)模型如式（2）：

式中，Ufall為陡脈沖沖激下降沿的幅值，s為陡脈沖沖激的峰值電壓，R、C分別表示人體等效電阻和等效電容，t0為陡脈沖沖激后支起點(diǎn)時間.此外，在醫(yī)療環(huán)境中，隨機(jī)噪聲伴隨陡脈沖干擾心電信號.隨機(jī)噪聲的能量均勻的分布在整個頻帶上，在此以nw表示.經(jīng)過以上分析可知，受陡脈沖相關(guān)噪聲干擾的心電信號的數(shù)學(xué)模型如式（3）：

其中，xreal與xideal分別表示實際心電信號和理想心電信號，npluse表示陡脈沖干擾.

圖1 受陡脈沖干擾的ECG.（a）受到陡脈沖干擾的心電信號；（b）未受陡脈沖干擾的心電信號Fig.1 ECG disturbed by steep pulses: (a) ECG disturbed by steep pulses; (b) ECG not disturbed by steep pulses

2 基于VMD的陡脈沖干擾抑制

VMD算法最早由Dragomiretskiy等于2014年提出[10]，可以將任意信號x(t)分解成離散數(shù)量K的子信號uk（Band-limited intrinsic mode functions, BLIMF）.每一個子信號的能量都集中在相應(yīng)的中心頻率wk周圍，具有特殊的稀疏特性[10].其基本表達(dá)式如式（4）.這是一種后驗的、非遞歸算法，并具有帶通濾波器屬性[12].目前，這種算法被廣泛的應(yīng)用在軸承故障檢測[13]，地震信號分析[14]，經(jīng)濟(jì)走勢預(yù)測[15]等領(lǐng)域.

式中，K表示分割層數(shù)，uk表示第k個子信號，wk表示uk的中心頻率，x表示原信號，δ(t)為狄拉克函數(shù)，?t代表梯度運(yùn)算，*代表卷積運(yùn)算，j為虛數(shù)符號.

2.1 心電信號的預(yù)處理

根據(jù)香農(nóng)采樣定律，信號的頻帶寬度是采樣頻率的0.5倍.以500 Hz的采樣頻率為例，實際信號的寬度為0～250 Hz.而心電信號的主要成分集中于0.05～49 Hz[16].因此，不對50～250 Hz的信號進(jìn)行分析，以提高實時性.在此，采用了一個4階IIR濾波器對這部分信號進(jìn)行衰減，截止頻率為50 Hz.

2.2 陡脈沖干擾的抑制

對式（2）進(jìn)行傅里葉變換，可得其幅頻曲線A(f)，如式（5）.

式中，f表示頻率.由于心電采樣設(shè)備通常要求大輸入電阻[17]，因此1/RC為常數(shù)且較小，在此為方便分析將其忽略.幅頻特性如圖2.

圖2 衰減模型單邊幅頻圖Fig.2 Diagram of a single-sided amplitude-frequency attenuationmodel

可以看到，陡脈沖干擾的后支能量主要分布在10 Hz前.VMD算法基于對混合信號的主成分估計，結(jié)合其帶通濾波器屬性，包含陡脈沖干擾后支的子信號對應(yīng)中心頻率wk應(yīng)小于10 Hz.由于此部分干擾對心電波形干擾較大，是抑制的重點(diǎn).圖3展示了當(dāng)K=6時，各子信號時域圖.其中，子信號u1的中心頻率w1接近于0，包含直流分量、基線漂移、陡脈沖干擾的低頻干擾分量.將u1拋棄，可以有效地消除基線漂移和直流分量.u2的中心頻率滿足w2<10 Hz，包含陡脈沖干擾分量、P波、T波分量[16].完全拋棄u2不利于恢復(fù)心電信號.

圖3 VMD 6階分解子信號.（a）原信號；（b）u1；（c）u2；（d）u3；（e）u4；（f）u5；（g）u6Fig.3 Use of VMD to divide the signal into six layers: (a) original signal; (b) u1; (c) u2; (d) u3; (e) u4; (f) u5; (g) u6

設(shè)心電采集系統(tǒng)共模抑制比為-100 dB，陡脈沖的峰值大于1000 V，理論上產(chǎn)生噪聲的幅值約大于10 mV.實際測得幅值由于損耗而略小.以肢體導(dǎo)聯(lián)為例，P波幅值不超過0.25 mV，T波幅值不超過0.5 mV[16]，遠(yuǎn)低于同頻段陡脈沖干擾幅值.基于以上分析，本文設(shè)計了一種自適應(yīng)閾值，對陡脈沖干擾在u2中的陡脈沖干擾分量進(jìn)行消除.閾值λsp如式（6）.

上式表示對子信號u2中絕對幅值最大的10個采樣點(diǎn)sort(|u2|)[1:10]取 均值.sort(||)表示對輸入信號序列的絕對值降序排列，mean(·)表示均值運(yùn)算.

圖4 u2中干擾分量消除效果.（a） 閾值處理前；（b） 閾值處理后Fig.4 Interference component elimination effect in u2: (a) before using thresholds; (b) after using thresholds

對于陡脈沖干擾前支，由于其在時域上與后支存在緊密關(guān)聯(lián)，在u2之后的子信號中，以[zon-Δt,zon+Δt]為范圍，搜索極大值（由陡脈沖特性，陡脈沖干擾前支時域上表現(xiàn)為時長短的特點(diǎn)，因此，Δt取500 Hz采樣頻率下的5個采樣點(diǎn)，即10 ms）.將極大值所屬脈沖置零以完全消除陡脈沖沖激.

2.3 消除隨機(jī)噪聲

隨機(jī)噪聲由電子器件產(chǎn)生，在頻域分布均勻，根據(jù)VMD算法的帶通濾波器組特性[12]，uk中的隨機(jī)噪聲分量的幅值遠(yuǎn)遠(yuǎn)小于原信號隨機(jī)噪聲幅值.受基于wavelet的閾值濾波算法[8]啟發(fā)，本文提出了一種基于VMD的閾值法去隨機(jī)噪聲算法.由于P、T波幅值小，采用閾值消除相應(yīng)子信號中的隨機(jī)噪聲分量可能得不償失.因此需uk對應(yīng)的wk滿足式（8）：

式中，won和woff分別對應(yīng)P、T波頻帶的上限與心電信號主要頻帶的上限，分別取12 Hz和49 Hz[16].對滿足條件的子信號uk施加閾值，如式（9）.

其中，i表示采樣點(diǎn)序號，T為采樣周期，sgn(·)代表階躍函數(shù)，λk代表uk的閾值.如下式（10）所示.

式中，N代表uk中采樣點(diǎn)的數(shù)量，a和b均為常數(shù)[18].τ為比例因子.施加閾值后效果如圖5.

最終，抑制噪聲后的心電信號用xfilter表示.如式（11）所示：

圖5 消除隨機(jī)噪聲分量之后的子信號.（a） u4；（b） u5Fig.5 Sub-signal after eliminating random noise components: (a) u4;(b) u5

式中，w是對心電信號中隨機(jī)噪聲的估計，pluse是算法對陡脈沖沖激的估計.表示消除隨機(jī)噪聲后的uk.

3 基于VMD算法的QRS提取

在心電信號中，QRS波群的能量集中在8～16 Hz之間[16].如式（12）所示，則可以提取QRS波群主要成分.

在重組信號xQRS的基礎(chǔ)上，對QRS波群進(jìn)行檢測，可避免P、T波等信號對檢測過程的干擾.為了突出R波的位置，對重組信號進(jìn)行預(yù)處理.預(yù)處理算法如式（13）、（14）、（15）所示，依次對重組信號序列表示差分、平方、積分運(yùn)算.xdiff、xsq、xmvw分別是差分、平方、積分運(yùn)算的輸出.式中，T表示采樣周期，W表示移動滑窗寬度，在此取80 ms效果較為理想，這與R波的實際寬度有關(guān)[19].效果如圖6（c）所示.

圖6 R波檢測的算法驗證.（a） 原信號；（b） 重組的QRS；（c） R波增強(qiáng)Fig.6 Verification of the R-wave detection algorithm: (a) original signal; (b) restructured QRS wave; (c) enhanced R wave

處理之后的心電信號，采用閾值法對R波峰值點(diǎn)進(jìn)行提取.初始閾值λ0由最初10個1 s內(nèi)預(yù)處理后序列極大值的均值決定.隨后，每檢測到一個R波即對閾值進(jìn)行更新.閾值更新算法如式（16）.

式中，λR為當(dāng)前R波的檢測閾值，Pnoise和PQRS分別是最近n個非R峰值點(diǎn)幅值的均值和最近n個R峰值點(diǎn)幅值的均值.γ為比例常數(shù)，當(dāng)n取8時，γ取0.475效果最佳[19].算法流程圖如圖7.圖中，rm-1和rm分別表示第m-1個和第m個已經(jīng)被定位的R波所在位置，rrmean表示連續(xù)n個RR間期長度；tRP為不應(yīng)期，此期間無R波，時長為200 ms[20]；max[.]表示求括號中的最大值.

4 實驗和數(shù)據(jù)分析

實驗部分采用MIT-BIH數(shù)據(jù)庫數(shù)據(jù)與實際采集的心電信號共同驗證本文中提出的算法.圖8（a）所示，是天津市機(jī)器人與智能重點(diǎn)實驗室研制的不可逆電脈沖消融設(shè)備樣機(jī)，可以釋放峰峰值范圍1000～3000 V的陡脈沖.為采集受陡脈沖相關(guān)噪聲干擾的心電信號，本文采用TI公司ADS1298芯片，自行設(shè)計制作了心電信號采集設(shè)備，如圖8（b）所示.本文的實驗部分，陡脈沖的峰峰值為1000 V.

圖7 R波檢測的算法流程圖Fig.7 Flowchart of the R-wave detection algorithm

圖8 實驗設(shè)備展示.（a） 不可逆電消融樣機(jī)；（b） 心電信號采集板Fig.8 Experimental equipment: (a) prototype of irreversible electrical pulse ablation surgery; (b) ECG signal sampling circuit

4.1 ECG濾波算法實驗

4.1.1 預(yù)處理對實時性的提升

實驗平臺基本參數(shù)：主頻2.9 GHz，i5處理器，運(yùn)行內(nèi)存8 G，Windows10操作系統(tǒng).本次實驗選用500 Hz采樣頻率下16.384 s的包含陡脈沖干擾的ECG信號為實驗對象.

圖9對比了心電信號是否預(yù)處理對w2的影響.根據(jù)VMD相關(guān)理論，分解層數(shù)K應(yīng)當(dāng)盡可能的少以節(jié)約計算時間；同時u1表示基線信息因而w1約等于0，w2的值意味著分解的有效性.當(dāng)w2<10 Hz時有效分離高壓陡脈沖干擾.若不進(jìn)行預(yù)處理，分解層數(shù)K≥14時可將陡脈沖干擾與QRS波群分離.當(dāng)K=14時，耗時30.1939 s.如對ECG進(jìn)行預(yù)處理，當(dāng)K≥5時即滿足要求.當(dāng)K=5時，僅需2.511 s.顯然，對帶噪心電信號預(yù)處理可以有效提高實時性且并較少內(nèi)存耗費(fèi).

圖9 信號預(yù)處理效果對比Fig.9 Effect of signal preprocessing interference

4.1.2 懲罰因子α和分解層數(shù)K的確定

懲罰因子α和分解層數(shù)K是VMD算法兩個重要的參數(shù)，決定了濾波器通頻帶寬度[10].由圖9，對心電信號進(jìn)行預(yù)處理的條件下，分解層數(shù)K應(yīng)滿足{K|K∈[5,9]}.為進(jìn)一步確定參數(shù)，設(shè)計以下實驗：選取MIT-BIH數(shù)據(jù)庫中第100號數(shù)據(jù)，取7200個采樣點(diǎn)，在隨機(jī)噪聲級別均為5 dB的條件下，添加一定數(shù)量的模擬高壓陡脈沖干擾.在此，以信噪比（Signal to noise ratio, SNR）作為評價指標(biāo)，如式（17）.式中，xclean表示未添加噪聲的心電信號，xnoise表示添加了噪聲后的心電信號，Nal表示全部參加測試的采樣點(diǎn)數(shù).

實驗結(jié)果如圖10所示.實驗結(jié)果表明，當(dāng)K=6且α=2500時，濾波效果較為理想.

圖10 對分解層數(shù)和懲罰因子進(jìn)行驗證Fig.10 Values of K and quadratic penalty were verified through experiments

4.1.3 陡脈沖干擾抑制

圖11為對心電信號中陡脈沖干擾抑制效果.其中，圖11（a）為原信號圖，圖11（b）和圖11（c）分別是選擇保留和不保留陡脈沖前支的心電信號圖.可以看到，運(yùn)用本算法，對陡脈沖沖激有著明顯的消除效果，且最大程度上保持了心電信號應(yīng)有波形.圖12對比了本文提出算法與小波包閾值去燥[8]、帶阻濾波算法（采用FIR帶阻濾波器濾除0～10 Hz信號）在消除陡脈沖干擾時的效果.其中圖12（a）為本文提出算法對陡脈沖干擾的消除效果，圖12（b）為小波包閾值去燥算法對陡脈沖干擾的消除效果，雖然這種算法相對較好的保留了心電信號的特征，但是對陡脈沖干擾消除并不徹底；圖12（c）為采用帶阻濾波器算法對陡脈沖干擾的消除效果，此算法則完全不能保留心電信號特征，且抑噪效果差.

為進(jìn)一步驗證算法效果，選擇MIT-BIH數(shù)據(jù)庫中10段含噪聲較少的心電數(shù)據(jù)，每段數(shù)據(jù)長度為16.384 s；由式（5）的數(shù)學(xué)模型，分別向心電信號中加入不同數(shù)量的模擬高壓陡脈沖干擾.分別采用本文提出的算法、小波包閾值法、帶阻濾波器對噪聲進(jìn)行處理.實驗結(jié)果如表1所示.

圖11 施放陡脈沖期間相關(guān)噪聲抑制.（a） 陡脈沖干擾消除前；（b）消除陡脈沖干擾后支；（c） 完全消除陡脈沖干擾Fig.11 Result of eliminating steep pulse interference: (a) ECG disturbed by steep pulses; (b) forehead is retained; (c) fore branches are eliminated

表中，第一列表示加入不同數(shù)量模擬陡脈沖干擾后的原始信號信噪比，其后三列分別表示采用不同濾波算法處理后的濾波效果，以信噪比表示.可以從以上實驗結(jié)果中得出結(jié)論，本文設(shè)計算法可有效去除心電信號中的高壓陡脈沖干擾，且相比于其他算法具有明顯的優(yōu)勢.

4.1.4 隨機(jī)噪聲消除效果

閾值公式（10）由文獻(xiàn)[18]提出，式中參數(shù)a，b分別取0.083和1.39.再此對參數(shù)正確性進(jìn)行驗證.采用一段來自MIT-BIH數(shù)據(jù)庫第103條心電數(shù)據(jù)，加入一定數(shù)量的模擬陡脈沖干擾以及5 dB的隨機(jī)噪聲模擬真實環(huán)境.實驗結(jié)果如圖13所示：當(dāng)a和b改變時，信噪比不同程度降低.由此，可以認(rèn)為文獻(xiàn)[18]中a，b的取值是無誤的.

圖12 陡脈沖干擾抑制效果對比.（a） 本文提出算法效果；（b） 小波算法效果；（c） 帶阻濾波器算法效果Fig.12 Comparison of the results of steep pulse noise filtering: (a)algorithm in this article; (b) wavelet algorithm; (c) bandpass filter

再分別將信噪比為1～10 dB的隨機(jī)噪聲和相同數(shù)量的模擬陡脈沖加入到等長的心電數(shù)據(jù)中，以此數(shù)據(jù)對算法進(jìn)行對比驗證.表2為將本文提出的算法與wavelet去噪法[8]、EMD去燥法[9]的效果對比.信噪比越大說明信號包含隨機(jī)噪聲的強(qiáng)度越小，也就意味著對噪聲的消除效果越好.

可以得出結(jié)論，在消除消融設(shè)備工作帶來的隨機(jī)噪聲方面，本算法可以明顯的提高信號的信噪比；相較于另外兩種算法，提升效果更佳明顯.

4.2 QRS波群檢測實驗和數(shù)據(jù)分析

對于R波檢測準(zhǔn)確度，采用兩個指標(biāo)進(jìn)行評價：

（1）誤檢率（Error detection, ED）：表達(dá)式如式（18）所示；

（2）敏感度（Sensitivity, SE）：表達(dá)式如式（19）所示.

表1 陡脈沖干擾消除效果對比Table 1 Results obtained after a variety of filtering processes

圖13 a，b最優(yōu)值驗證.（a） b恒定a改變；（b） a恒定b改變Fig.13 Verification of the optimal values of a and b: (a) when b is constant and a is changed; (b) when a is constant and b is changed

式中，TP表示正確檢測，F(xiàn)P表示漏檢的情況，F(xiàn)N表示誤檢的情況，Ntotal表示R波總數(shù)量.

為驗證本算法的實際表現(xiàn)，采用8段來自不同患者的心電信號數(shù)據(jù).每組數(shù)據(jù)集的時間長度均為30 min.考慮相關(guān)陡脈沖療法的治療位置可能位于軀干，因此，采用肢體導(dǎo)聯(lián)作為心電信號來源，避免干擾正常治療過程.經(jīng)過人工標(biāo)記，共標(biāo)記16218個R峰值點(diǎn).此期間，共記錄釋放陡脈沖522次.將本文提出算法與差分閾值法[20]，帶通濾波器法[21]的對比結(jié)果，實驗數(shù)據(jù)如表3所示.經(jīng)過對比發(fā)現(xiàn)，相比于本文提出的算法，其他兩種算法的誤檢率遠(yuǎn)高于本算法，這與未對陡脈沖干擾進(jìn)行抑制由直接關(guān)系.可以通過數(shù)據(jù)得出結(jié)論，本文提出的算法，在應(yīng)對陡脈沖干擾下的心電信號時，有著更好的準(zhǔn)確度和敏感度.

表2 去噪效果對比Table 2 Comparison of the results of the filtering processes

表3 與常用算法準(zhǔn)確度對比Table 3 Comparison of the accuracy of the algorithms

5 結(jié)論

本文以陡脈沖的應(yīng)用為背景，對這種特殊環(huán)境下采集的心電信號進(jìn)行了研究.創(chuàng)新點(diǎn)在于：針對陡脈沖干擾特性，本文提出了一種基于VMD的心電信號陡脈沖干擾抑制算法；基于VMD特性，提出了一種閾值算法，有效的消除了電子設(shè)備帶來的隨機(jī)噪聲；并根據(jù)分割后子信號的頻域特性，提出了一種基于VMD的QRS波群提取算法.最終，以在天津市機(jī)器人與智能重點(diǎn)實驗室研制的不可逆電消融設(shè)備干擾下采集的心電信號，與來自MIT-BIH數(shù)據(jù)庫中的心電數(shù)據(jù)共同對算法性能進(jìn)行了驗證.結(jié)果證明，算法可以實現(xiàn)預(yù)期目標(biāo)，且效果較為理想.其中，對陡脈沖相關(guān)噪聲干擾下的ECG，提取QRS波群的準(zhǔn)確度達(dá)到了99%以上.