董騰飛，潘 松，吳 松，郝珠珠，趙偉杰，楊 勇

(杭州電子科技大學(xué)自動化學(xué)院，杭州310018)

基于不可逆電穿孔技術(shù)(Irreversible Electroporation，IRE)的納米刀在治療肝癌、胰腺癌以及其他惡性腫瘤的場景下發(fā)揮著越來越強大的作用[1-2]。但是IRE治療中釋放的高壓脈沖對心電信號的產(chǎn)生和傳導(dǎo)會產(chǎn)生極大干擾，增大了引發(fā)心律失常等術(shù)中并發(fā)癥的機率，故需要對患者的心電進行實時監(jiān)測[3]，并使IRE治療脈沖在心臟的有效不應(yīng)期內(nèi)釋放，減小對心臟的正常工作過程的影響。故能否準(zhǔn)確識別患者心臟的有效不應(yīng)期，及時觸發(fā)納米刀釋放高壓脈沖，對降低術(shù)中并發(fā)癥有著重要的意義。

由于在術(shù)中無法直接對心臟的有效不應(yīng)期進行檢測，本文根據(jù)體表心電圖和有效不應(yīng)期的對應(yīng)關(guān)系，將心臟有效不應(yīng)期的實時檢測問題轉(zhuǎn)化為R波實時檢測問題。但由個體差異導(dǎo)致的心電波形變化，以及由高壓脈沖刺激引起的骨骼肌強直收縮或心率失常等情況，會使采集到的心電波形與正常的心電波形產(chǎn)生較大的差異，從而給R波的準(zhǔn)確識別帶來困難；此外，通過檢測心臟有效不應(yīng)期來指導(dǎo)IRE治療脈沖的釋放時機，要求輸出R波檢測結(jié)果的實時性較高，算法的計算量較小。

目前，國內(nèi)外常用的臨床心電信號R波檢測方法有硬件方法和軟件方法兩種。硬件方法響應(yīng)速度快、實時性好，但是缺乏靈活性，魯棒性較低。軟件方法有時域分析法[4]、模板匹配法[5]、小波變換法[6]和神經(jīng)網(wǎng)絡(luò)法[7-9]等。時域分析法是應(yīng)用最廣、實現(xiàn)最簡單、計算量最小的心電信號分析方法之一，但是變異的大P波、大T波，無法被濾波器濾除干凈的噪聲都有可能會降低檢測結(jié)果的準(zhǔn)確率，故該方法還有進一步提升的空間；小波變換是另一種應(yīng)用較廣、研究較深的方法，但是計算量大，在實際應(yīng)用過程中軟硬件設(shè)計復(fù)雜，實現(xiàn)困難；神經(jīng)網(wǎng)絡(luò)法多用于波形識別、病例診斷方向且計算量大。

以上幾種方法并不能滿足心臟有效不應(yīng)期檢測高準(zhǔn)確率、高實時性、低計算量的要求。本文在時域分析法的基礎(chǔ)上，從預(yù)處理步驟和R波識別步驟兩個角度出發(fā)，以提高實時性和準(zhǔn)確率、簡化計算量為目標(biāo)對算法加以改進創(chuàng)新，設(shè)計出了一套可以應(yīng)用于IRE治療的R波實時檢測算法，從而達(dá)到識別延時低于15ms，且當(dāng)識別到含有大量噪聲的心電數(shù)據(jù)段時自動停止R波檢測，提高算法魯棒性的要求，進而通過固定延時實現(xiàn)了心臟有效不應(yīng)期的實時檢測定位，指導(dǎo)IRE治療脈沖釋放時機。

1 心電信號預(yù)處理

心電信號是微弱的低頻信號，對外界干擾比較敏感[10]，故需要采用一定的預(yù)處理步驟濾除噪聲，突出目標(biāo)信號特征，簡化計算量。這是準(zhǔn)確分析心電信號的重要前提，本文采用的預(yù)處理流程如圖1所示。

圖1 預(yù)處理流程圖

1.1 頻率歸一化處理

MIT-BIH心率失常數(shù)據(jù)庫[11]的采樣頻率為360 Hz，QT數(shù)據(jù)庫[12]的采樣頻率為250 Hz，采用頻率不同會影響算法的預(yù)處理效果，故采用重采樣的方法歸一化為500 Hz。

1.2 濾波處理

由于本算法實時性要求較高，且只需要檢測R波，故在保證檢測效果的前提下，可以選用低階IIR濾波器級聯(lián)的方法進行濾波，降低有效信號頻段通過濾波器的群延時。

常見心電信號噪聲主要由肌電干擾、基線漂移、50 Hz工頻干擾組成，本算法的濾波部分主要針對這幾類噪聲進行設(shè)計，濾波器參數(shù)如表1所示。

表1 濾波器參數(shù)

已知心電信號的主要能量集中在0～40(±7)Hz，QRS復(fù)波的主要能量在 6 Hz～18 Hz，使用MATLAB的fdatool工具分析得知，心電信號經(jīng)過低通濾波器后，0～40 Hz的頻率分量群延時約6 ms；經(jīng)過高通濾波器后，5 Hz以上的頻率分量群延時低于5 ms；經(jīng)過帶阻濾波器后，40 Hz以下的頻率分量群延時低于3 ms。所以經(jīng)過濾波處理后，QRS波段的群延時理論上不超過15 ms，且波形不會發(fā)生較大的非線性失真，滿足了R波識別的實時性要求，為進一步的預(yù)處理做好了準(zhǔn)備。

1.3 微分處理

本文使用微分函數(shù)來突出QRS波群在正常心電信號中斜率最大這一特征，將QRS波群從其他心電信號中區(qū)分出來。微分函數(shù)為:

式中:sig_f為經(jīng)過濾波器預(yù)處理后的心電數(shù)據(jù)；sig_d為對sig_f進行微分處理后的數(shù)據(jù)。

1.4 低閾值歸零處理

由于本文只需要檢測R波，故可采用低閾值歸零法對其進行處理，突出R波微分值較大的特征，同時濾除部分噪聲干擾，簡化算法需要處理的數(shù)據(jù)量，提高處理速度。低閾值歸零處理的公式為:

式中:sig_z為低閾值歸零處理后的數(shù)據(jù)；Thr_z為歸零閾值。

1.5 平方處理

本文采取平方處理對信號進行非線性濾波，進一步突出R波微分值較大的特征，減弱其他波形或噪聲的干擾。

式中:sig_s為sig_z平方后的結(jié)果。

1.6 移動窗口平均處理

對于正常心電信號而言，經(jīng)過以上步驟的處理，R波分量已經(jīng)非常突出，但當(dāng)原始信號的信噪比較低或噪聲的頻段與心電信號的頻段相重疊時，會導(dǎo)致濾波效果不佳，影響R波檢測的準(zhǔn)確率。

本文利用移動窗口平均法的特性將每個心跳周期分為有效信號區(qū)段和疑似噪聲區(qū)段，并分別對其進行處理，達(dá)到提高檢測結(jié)果陽性預(yù)測度的目的。移動窗口平均的運算公式為:

式中:sig_w為移動窗口平均后的數(shù)據(jù)；N為移動窗口的窗寬。當(dāng)窗寬設(shè)計過大時，會導(dǎo)致噪聲檢測過于靈敏，將含有正常R波的有效信號區(qū)段誤判為疑似噪聲區(qū)段而造成漏檢；當(dāng)窗寬設(shè)計過小時，會使得噪聲檢測不夠靈敏，將疑似噪聲區(qū)段誤判為有效信號區(qū)段而造成誤判。正常人心電信號QRS波段最長不超過0.12 s，故窗寬時長選為0.15 s，采樣頻率為500 Hz時N的值為75。

2 R波檢測算法

本文通過對MIT-BIH心率失常數(shù)據(jù)庫和QT數(shù)據(jù)庫的數(shù)據(jù)進行分析，設(shè)計了一套完整的R波實時檢測方案，同時對部分變異心電信號和含有大量噪聲的心電信號進行特殊處理，提高檢測的陽性預(yù)測度。R波檢測算法的整體框架圖如圖2所示。

圖2 R波實時檢測算法流程圖

2.1 閾值初始化

當(dāng)sig_f的長度超過5 s時，以1 s為間隔將其分為5段，保證每段中至少包含一個QRS波群，然后分別求最大值，再取中間大小的三個值求平均，記為SPK_f，根據(jù)SPK_f計算幅度閾值Thr_f，計算公式為:

微分閾值Thr_d1、Thr_d2和移動窗口平均閾值Thr_w的計算方法與上述過程相同，計算公式如下:

式中:SPK_d為對前5 s微分?jǐn)?shù)據(jù)sig_d處理后的平均值；SPK_w為對前5 s移動窗口平均值sig_w處理后的平均值；a為比例系數(shù)，一般選取a＝0.035。

低閾值歸零處理的歸零閾值Thr_z計算公式如下:

式中:b為比例系數(shù)，一般選取b＝0.3。

2.2 噪聲及室顫檢測

本算法提出根據(jù)移動窗口平均后的數(shù)據(jù)sig_w與閾值Thr_w的關(guān)系，將每個心跳周期分為有效信號區(qū)段和疑似噪聲區(qū)段。

第一階段為有效信號區(qū)段，噪聲標(biāo)志位為0，當(dāng)sig_w持續(xù)大于Thr_w的時間不超過150 ms時，進行R波檢測，否則進入疑似噪聲區(qū)段。程序運行框圖如圖3所示。

圖3 有效信號區(qū)段程序流程圖

第二階段為疑似噪聲區(qū)段，噪聲標(biāo)志位為1，當(dāng)sig_w持續(xù)大于Thr_w的時間超過500 ms時，表明心電信號中含有較大噪聲或患者出現(xiàn)心率過快、室顫等病癥，此時為了防止誤判暫停漏檢機制，直到檢測到下一個R波；如果持續(xù)時間超過5 s，則啟動糾正機制；當(dāng)sig_w持續(xù)小于Thr_w的時間超過150 ms時，返回第一階段，開始R波檢測。程序運行框圖如圖4所示。

如圖5(a)、5(b)所示，本算法能夠有效識別MIT-BIH心率失常數(shù)據(jù)庫104號數(shù)據(jù)中含有噪聲干擾的數(shù)據(jù)段和205號數(shù)據(jù)中含有異常心電波形的數(shù)據(jù)段，并自動停止檢測，待心電信號恢復(fù)正常后，再重新開始檢測。

圖4 疑似噪聲區(qū)程序流程圖

圖5 噪聲數(shù)據(jù)

2.3 針對 QS、Qr、rS、rSr波形的處理

在一組數(shù)據(jù)中，若所有的心電波形都是如圖6所示的 QS、Qr、rS、rSr′波，可以采用更換導(dǎo)聯(lián)的方法顯示向上的qRs波波形；若僅有部分心電波形是此類波形，則表明患者出現(xiàn)早搏或其他癥狀，需要重新判斷是否適合進行IRE治療，且根據(jù)這類心電波形來判斷心臟有效不應(yīng)期位置的準(zhǔn)確率較低，故當(dāng)檢測到此類波形時，暫停R波檢測。

算法實現(xiàn)如下:若存在某一時刻的心電數(shù)據(jù)同時滿足式(10)和式(11)，且該時刻距離上一個R波大于200 ms，則判斷出現(xiàn)此類心電波形。根據(jù)心臟有效不應(yīng)期原則，同時為了避免將此類異常波形之后出現(xiàn)的大T波誤判為R波，采取延時360 ms后再繼續(xù)檢測R波，且暫停漏檢機制的策略，效果如圖7所示。

圖 6 QS、Qr、rS、rSr心電波形

圖7 233號數(shù)據(jù)

2.4 平均RR間期

取最近8個RR間期的平均值為RR_mean，計算公式為:

式中:RR(i)為兩個R波峰之間的采樣點數(shù)。

2.5 R波檢測

本算法在幅值、微分值兩個維度設(shè)立動態(tài)閾值檢測R波，方法如下:①如果連續(xù)四個點的微分值sig_d大于閾值Thr_d1，且該點距離相鄰的R波超過200 ms，則表示檢測到一個R波上升沿；②存在一點的微分值sig_d小于0，幅值sig_f大于閾值 Thr_f；③從滿足條件1的點開始到找到一個滿足條件2的點為止的這段時間內(nèi)，任何一點的幅值和微分值都在上一個R波幅值和微分值的0.4倍到2倍之間。

當(dāng)某一時刻的檢測結(jié)果滿足條件1，且50 ms內(nèi)存在一點能夠同時滿足條件2和條件3，則判斷找到一個新的R波波峰。其中，條件3是針對心電信號采集過程中突然出現(xiàn)的接觸噪聲或其他能夠引起心電采集系統(tǒng)檢測到信號突變的情況而設(shè)立的。

由于QRS波的特征參數(shù)會隨著生理情況變化而變化，故采用動態(tài)閾值法提高算法魯棒性。在檢測到R波后同步更新閾值 Thr_f、Thr_d1、Thr_d2、Thr_w，動態(tài)閾值更新公式如下。

式中:sig_f_max、sig_d_max、sig_w_max 分別為檢測到R波上升沿到檢測到R波波峰這段時間內(nèi)sig_f、sig_d和sig_w的最大值。

2.6 漏檢機制及糾正機制

漏檢機制:由于實時性要求較高，本文并未采用具有較大延時的修正和回檢算法提高檢出率，而是通過動態(tài)閾值法減少漏檢。具體算法如下:

如果持續(xù)1.66倍RR間期沒有檢測到R波，則將微分閾值Thr_d1降低一半至Thr_d2，重復(fù)1.3.5節(jié)的算法繼續(xù)進行檢測，檢測到R波后更新閾值Thr_f、Thr_d1、Thr_d2、Thr_w。 其中 Thr_f、Thr_d1、Thr_d2計算公式如下，Thr_w按照式(19)和式(20)進行更新。

糾正機制:如果持續(xù)5 s沒有檢測到R波且噪聲標(biāo)志位為1，則停止檢測并提示原始信號信噪比過低；如果持續(xù)5 s沒有檢測到R波且噪聲標(biāo)志位為0，則暫停檢測并自動重新初始化閾值，再利用新的閾值進行檢測。

3 有效不應(yīng)期檢測

心臟的有效不應(yīng)期是指使用3～5倍的閾刺激值對心肌細(xì)胞或組織進行刺激時，不會引起興奮反應(yīng)的一段時間，約200 ms～300 ms，對應(yīng)于單細(xì)胞動作電位的0～2相和3相的前部，以心室肌為例相當(dāng)于從體表心電圖的QRS波開始一直持續(xù)到T波的前支。本文所述心臟有效不應(yīng)期檢測特指根據(jù)體表心電圖檢測心室肌有效不應(yīng)期。同時，在心房肌、心室肌相對不應(yīng)期剛開始的一段的時間內(nèi)，應(yīng)用較強的閾上刺激會有較大概率引發(fā)心房或心室顫動，這段期間被稱為易顫期[13-14]。在對人體施加高于閾刺激值的電刺激時，應(yīng)當(dāng)避開心房肌和心室肌的易顫期。

我們可以根據(jù)體表心電圖和有效不應(yīng)期的對應(yīng)關(guān)系，在檢測到R波后延遲一段時間避開心房肌和心室肌的易顫期，準(zhǔn)確地將IRE治療脈沖的釋放時間點控制在心臟有效不應(yīng)期內(nèi)。通過對算法延時、個體差異等因素的綜合考慮，延時時間可以設(shè)置為50 ms～100 ms。

單細(xì)胞動作電位和心電圖與不應(yīng)期的對應(yīng)關(guān)系如圖8所示[15]，圖中兩塊灰色區(qū)域分別表示心房肌易顫期和心室肌易顫期；a點表示R波波峰的位置；b點表示經(jīng)過計算和固定延時后釋放IRE治療脈沖的位置。

圖8 單細(xì)胞動作電位和心電圖與不應(yīng)期的對應(yīng)關(guān)系

4 實驗結(jié)果分析與討論

4.1 數(shù)據(jù)庫信號仿真實驗

QT數(shù)據(jù)庫標(biāo)注了各個波形的起止點和峰值點，最適合用來進行波形定位準(zhǔn)確性的評估，但部分標(biāo)注的位置有較大誤差，本文根據(jù)標(biāo)記點位置的準(zhǔn)確性，人工挑選出其中67組數(shù)據(jù)作為數(shù)據(jù)源進行仿真。

4.1.1 評價標(biāo)準(zhǔn)

4.1.1.1 準(zhǔn)確性的評價標(biāo)準(zhǔn)

為了分析R波實時檢測算法的定位準(zhǔn)確性，本文采用如下幾個參數(shù):0257真陽性(TP):準(zhǔn)確檢測出的R波峰；0258假陽性(FP):非R波峰誤判為R波峰；0259假陰性(FN):漏檢的R波峰。

因為本算法要求檢測延時低于15 ms，故取消了具有較大延時才能實現(xiàn)的R波修正和回檢機制，且設(shè)計思路之一是通過過濾部分非標(biāo)準(zhǔn)R波心拍以及噪聲較大的數(shù)據(jù)段以提高算法的魯棒性，所以會出現(xiàn)部分R波漏檢的情況，漏檢數(shù)量由信號采樣質(zhì)量以及非標(biāo)準(zhǔn)R波心拍的數(shù)量決定，故不采用假陰性FN及其相關(guān)指標(biāo)作為定量判斷標(biāo)準(zhǔn)。

理論上本算法在R波峰后的15 ms內(nèi)輸出檢測結(jié)果，但由于數(shù)據(jù)庫官方標(biāo)記位置有略微的誤差，本文設(shè)定在檢測結(jié)果之前20 ms，之后10 ms這一區(qū)間內(nèi)存在官方R波標(biāo)記點即為真陽性TP，反之為假陽性FP。

最終，本文通過計算陽性預(yù)測度(＋P)來定量評估算法，陽性預(yù)測度計算公式如下:

4.1.1.2 實時性的評價標(biāo)準(zhǔn)

在分析實時性時，本文采用平均延時(Td)來表示輸出檢測結(jié)果的延遲時間，并采用如下兩個指標(biāo)對該參數(shù)進行分析:

①均值:

4.1.2 仿真結(jié)果及分析

采用MATLAB進行仿真的結(jié)果如表2所示。

表2 QT數(shù)據(jù)庫R波檢測結(jié)果

續(xù)表2

根據(jù)表2計算平均陽性預(yù)測度、平均延時和方差的結(jié)果如表3所示。

表3 仿真結(jié)果分析表

對以上兩表分析可知，仿真結(jié)果的平均陽性預(yù)測度為99.88%，最高為100%，最低為98.37%；平均延時為13.02 ms，最小為9.49 ms，最大為14.67 ms，達(dá)到了實時性和準(zhǔn)確性的要求。

4.2 實測信號檢測實驗

將該算法移植到自主開發(fā)的心電信號采集分析系統(tǒng)，并采用SKX-2000G＋型心電信號模擬儀和人體實測心電作為數(shù)據(jù)源進行檢測，實驗結(jié)果如下:

4.2.1 使用心電信號模擬器作為數(shù)據(jù)源的測試結(jié)果

如下圖所示，黑色實線代表心電波形；虛線代表不應(yīng)期檢測結(jié)果，虛線的上升沿代表檢測到R波的時刻，下降沿代表理論上IRE治療脈沖輸出的時刻，該時刻位于心臟的有效不應(yīng)期內(nèi)。

①正常心電和竇性心電?？梢?00%的準(zhǔn)確檢出R波，實現(xiàn)不應(yīng)期定位。

圖9 正常心電

圖10 竇性心電

②含有接觸噪聲。在含有大量噪聲的時間段內(nèi)，自動停止檢測R波，噪聲結(jié)束后再重新開始檢測R波。

③疊加50 Hz工頻干擾的心電信號。由圖可見，50 Hz工頻干擾被完美濾除，不影響R波的正常檢測效果。

④高大T波。當(dāng)出現(xiàn)高大T波時，能夠準(zhǔn)確檢測出R波的波峰。

⑤倒尖角波形，模擬 QS、Qr、rS、rSr這類波形，可以完美濾除，不進行R波識別。

圖11 含有大量噪聲的心電

圖12 疊加50 Hz工頻干擾的心電

圖13 大T波

圖14 倒尖角波形

4.2.2 使用心電信號模擬器作為數(shù)據(jù)源的測試結(jié)果

利用心電信號采集分析系統(tǒng)對人體的心電信號進行采集識別，結(jié)果表明本算法不但可以有效的對正常人體心電進行識別，而且當(dāng)心電信號含有較大的運動偽跡或其他噪聲干擾時，同樣可以起到較好的識別效果，具有較強的魯棒性。

圖15 正常心電

圖16 含有運動偽跡

圖17 含有大量噪聲的心電

如上圖所示，虛線的上升沿對應(yīng)于R波下降沿頂端的位置，表明本算法可以在極小的延時內(nèi)識別到R波，達(dá)到了實時性要求；折線的下降沿一般在S波上升沿過后的位置，表明本算法可以避開易顫期，準(zhǔn)確的將IRE治療脈沖釋放的時間控制在心臟有效不應(yīng)期內(nèi)，達(dá)到了較好的效果。

5 總結(jié)

針對IRE治療需要實時檢測患者的心臟有效不應(yīng)期這一需求，本文在時域分析法的基礎(chǔ)上進行改進和創(chuàng)新，實現(xiàn)了復(fù)雜噪聲條件下的R波實時檢測功能和心臟有效不應(yīng)期定位功能，平均陽性預(yù)測度達(dá)到99.84%，平均延遲在15 ms以內(nèi)，且能夠準(zhǔn)確識別并排除部分變異心電信號和噪聲干擾。對比于傳統(tǒng)心電信號檢測算法，在保證準(zhǔn)確率的同時，實時性有了明顯的提高，與楊衍菲等人針對除顫器設(shè)計的 QRS波斜率檢測方法相比[16]，平均延時從17.3 ms降低到了13.02 ms。

本文算法測試環(huán)境與實際IRE治療時的心電信號略有差異。首先，納米刀消融的高頻脈沖會為心電信號引入一些噪聲。關(guān)于這一問題，本團隊已經(jīng)開展基于動物實驗的相關(guān)研究，后續(xù)進展將在下一篇論文中呈現(xiàn)。其次，即使在精確識別心臟有效不應(yīng)期的情況，外加全身麻醉及深度神經(jīng)肌肉阻滯的情況下，納米刀手術(shù)過程中仍有可能由于高壓脈沖的作用而偶然引發(fā)自限性心律失常[17]。針對這一問題，算法設(shè)計檢測到非正常心電信號后，會發(fā)送指令停止一段時間的檢測，以暫停納米刀的輸出。當(dāng)專業(yè)醫(yī)療人員觀測到心電信號回歸正常可以繼續(xù)手術(shù)后，再次啟用算法進行檢測。

本文算法可用于指導(dǎo)臨床納米刀消融手術(shù)過程，降低對心臟正常生理活動的影響，提高手術(shù)安全性。下一步可以使用心電信號采集分析系統(tǒng)結(jié)合高頻脈沖發(fā)生器展開相關(guān)動物實驗，對本算法和心電信號采集分析系統(tǒng)做進一步的改進和優(yōu)化。