常戰(zhàn)國,蒲寶明,李相澤,王 帥,楊 朔

1(中國科學(xué)院大學(xué),北京 100049)

2(中國科學(xué)院 沈陽計算技術(shù)研究所,沈陽 110168)

3(東北大學(xué) 計算機(jī)科學(xué)與工程學(xué)院,沈陽 110004)

心電圖是臨床上醫(yī)生診斷心臟疾病的一個常規(guī)而 且高效的技術(shù)手段,心電圖各波與心肌動作電位有著密切關(guān)系.心肌梗死根據(jù)病變程度不同,在心電圖上可以體現(xiàn)出不同的形態(tài).目前在心肌梗死的檢測上,多采用傳統(tǒng)的信號處理方法,如傅里葉變換,小波變換等,該傳統(tǒng)方法對于波形質(zhì)量要求高,對于有噪聲干擾和形態(tài)不好的波形效果不佳.也有一些學(xué)者提出使用SVM[1],KNN,決策樹等機(jī)器學(xué)習(xí)方法[2],這一方法前期依賴于傳統(tǒng)方法,需要提取相關(guān)特征,對于非典型性心肌梗死的檢測[3],由于特征不明顯,故準(zhǔn)確率不高.目前也有學(xué)者單純使用卷積神網(wǎng)絡(luò)診斷心肌梗死,但該方法對于典型心肌梗死的檢測準(zhǔn)確率沒有直接提取特征的效果好,并且時間復(fù)雜度較大.由于非典型心肌梗死波形復(fù)雜多變,特征模糊不易提取,所以BP 神經(jīng)網(wǎng)絡(luò)在非典型心肌梗死的檢測上表現(xiàn)一般.

本文針對上述問題,綜合考慮準(zhǔn)確率和算法時間復(fù)雜度因素,提出BP 和CNN 結(jié)合的方法.該方法在非典型性和典型性心肌梗死的檢測上都具有較好效果.首先對ECG 信號做前期處理,主要包括濾波,去基線,差分等.之后針對典型心肌梗死,對各個波做準(zhǔn)確定位和識別,再根據(jù)心肌梗死的典型特征,提取出相應(yīng)的特征值,采用BP 神經(jīng)網(wǎng)絡(luò)訓(xùn)練；對于非典型性心肌梗死,由于特征模糊,在對ECG 信號前期處理后,采用卷積神經(jīng)網(wǎng)絡(luò)訓(xùn)練[4].算法整體流程如圖1所示.

圖1 算法流程圖

1 研究方法

本文主要從數(shù)據(jù)獲取,數(shù)據(jù)前期處理,波形檢測,特征提取,構(gòu)建神經(jīng)網(wǎng)絡(luò)及識別等幾個方面敘述.

1.1 數(shù)據(jù)獲取

數(shù)據(jù)主要來源于MIT-BIH 數(shù)據(jù)庫和中國醫(yī)科大學(xué)附屬第一醫(yī)院.

1)MIT-BIH 是國際上公認(rèn)的心電數(shù)據(jù)庫,由美國麻省理工大學(xué)提供.在其ECG Database 中包含PTB Diagnostic ECG Database 數(shù)據(jù)庫,其中有268 個病人數(shù)據(jù),心肌梗死的共有148 人.該數(shù)據(jù)庫以其固定格式編碼存儲,需要特定工具處理.每個病人數(shù)據(jù)中包含.dat(ECG 數(shù)據(jù)),.hea(病人詳細(xì)信息)和.xyz(Frank 導(dǎo)聯(lián)數(shù)據(jù))三種格式文件,采樣率為1000.除了常規(guī)的12 導(dǎo)聯(lián)外,還有vx,vy,vz三個Frank 導(dǎo)聯(lián),共15 個導(dǎo)聯(lián).

2)中國醫(yī)科大學(xué)附屬第一醫(yī)院提供了從2008年到2015年間心電科的病人數(shù)據(jù).從中提取出所有的心電數(shù)據(jù),并對心電數(shù)據(jù)按醫(yī)生診斷結(jié)果進(jìn)行詳細(xì)分類,提取出所有患有心肌梗死的病人的心電數(shù)據(jù),并進(jìn)一步核對數(shù)據(jù).其中,診斷為心肌梗死(包含疑似心肌梗死)的病人481 人.該數(shù)據(jù)為xml 格式,需要從中提取信息,文件包含病人姓名,性別,年齡,診斷結(jié)果等信息.采樣率為500,包含常規(guī)12 導(dǎo)聯(lián)數(shù)據(jù).

1.2 ECG 數(shù)據(jù)前期處理

對于MIT-BIH 需要使用WFDB 工具包根據(jù)不同數(shù)據(jù)格式提取相應(yīng)數(shù)據(jù).對于醫(yī)大一院數(shù)據(jù),需要從xml 中根據(jù)關(guān)鍵字段提取ECG 數(shù)據(jù)和診斷結(jié)果.之后根據(jù)二者不同采樣率和幅值設(shè)置不同的濾波系數(shù)等.數(shù)據(jù)先經(jīng)過帶通濾波器去除一定的噪聲干擾,之后經(jīng)過改進(jìn)的均值濾波平滑后,保證數(shù)據(jù)形態(tài)特征不發(fā)生較大改變.然后,將數(shù)據(jù)做去基線處理,防止基線漂移.最后將數(shù)據(jù)做放大處理,以便后期提取特征值.

1.2.1 濾波

本文根據(jù)濾波效果,波形特點(diǎn)和對時間復(fù)雜度的要求嘗試了多種濾波方法,最終比較了均值濾波和butter worth 濾波.

(1)均值濾波如下：

假如濾波器移動平均窗寬為5,經(jīng)過平滑處理返回值yy,則yy為：

(2)Butterworth 濾波如下：

其中,n是濾波器的階數(shù),wc是截止頻率,wp是通頻帶邊緣頻率,是|H(w)|2在通頻帶邊緣的數(shù)值.

如果令s=jw,拉普拉斯變換在虛軸s=jω上的性質(zhì)可以得到：

由極點(diǎn)在單位圓上可以得到：

因此傳遞函數(shù)為：

1到10 階的Butterworth 多項(xiàng)式因子可以通過查表得到.

圖2-圖4為兩種方法的對比結(jié).

圖2 原始ECG 信號

圖3 Smooth 平滑后ECG

圖4 Butterworth 濾波后ECG

從圖中可以看出,經(jīng)過均值濾波后,ECG 波形基本沒有發(fā)生形變,只是隨著平滑點(diǎn)數(shù)的增加,幅值有所降低,為此采用一種補(bǔ)償方法,將其補(bǔ)償?shù)皆瓉淼臄?shù)值.但是經(jīng)過Butterworth 濾波后,除了波形幅值有所降低外,在波形拐點(diǎn)處會有一定的形變,尤其在S 波會有一個明顯向下的尖小的波,因此導(dǎo)致波形在一定程度上有所失真.

1.2.2 去基線

由于存在基線漂移,對提取特征值造成很大干擾,尤其是心肌梗死判斷中,準(zhǔn)確識別出ST 段是否抬高或壓低是及其關(guān)鍵的,如果存在基線漂移,就不能準(zhǔn)確識別這一重要特征.在一個周期中采用插值方法找到基線,然后將所有數(shù)據(jù)減去基線,得到去完基線的ECG信號.

1.3 QRS 波群檢測

理論上每一個完整節(jié)拍ECG 信號都有P 波,QRS 波群,T 波組成.由于在一些特殊情況下,可能出現(xiàn)P 波消失以及T 波不明顯的情況,但都會有QRS 波群,所以,我們先檢測QRS 波群,然后,由QRS 波群確定P 波和T 波.文獻(xiàn)[5]提出一種改進(jìn)的DWT 算法來做波形匹配.首先,對波形求一階導(dǎo)數(shù),得到QRS 波群的大概位置序列,之后在對原始波形求二階導(dǎo),找到更為精確的位置,然后根據(jù)位置序列的索引范圍從原始波形中找到QRS 波群.之后,在經(jīng)過濾波等處理后的波形中按照索引范圍,以及原始波形高度的閾值等多維特征共同確定QRS 波群位置,并對其修正.圖5-圖7分別為濾波后原始數(shù)據(jù),一階導(dǎo)數(shù)據(jù)和二階導(dǎo)數(shù)據(jù).

圖5 濾波后ECG

圖6 一階導(dǎo)ECG

由圖可以看到,原始波形中QRS 波群波峰的位置索引靠近一階導(dǎo)峰值位置,理論上在一階導(dǎo)數(shù)值為0 對應(yīng)的位置上,而一階導(dǎo)為0 的位置變化率最快,因此與二階導(dǎo)谷值對應(yīng)索引更加接近,故選擇在二階導(dǎo)谷值附近查找原始信號的QRS 波峰.由QRS 波群進(jìn)而定位P 波和T 波.

圖7 二階導(dǎo)ECG

1.4 ECG 特征提取

在準(zhǔn)確識別了QRS 波群的前提下,需要對ECG波形做進(jìn)一步的特征提取.國際上公認(rèn)的常規(guī)12 導(dǎo)聯(lián)是標(biāo)準(zhǔn)I,II,III,加壓單極肢體導(dǎo)聯(lián)aVR,aVL,aVF 和單極胸腔壁導(dǎo)聯(lián)V1-V6[6].12 導(dǎo)聯(lián)之間存在如下向量運(yùn)算關(guān)系：

每一個單極導(dǎo)聯(lián)都由周期性的P,QRS,T 波等組成.個別還有U 波,但目前沒有發(fā)現(xiàn)明確的生理學(xué)意義(U 波異常可能與缺鉀有關(guān),但不是絕對的),我們暫不研究.一個完整心動周期如圖8所示.

圖8 心動周期示意圖

其中,P 波代表心房肌除極的電位變化,PR 間期代表心房開始除極的時間,QRS 波群代表心室肌除極的電位變化,ST 段代表心室緩慢復(fù)極的過程,T 波代表心室快速復(fù)極的電位變化.結(jié)合醫(yī)學(xué)知識和專家意見,為了能通過心電圖準(zhǔn)確判斷出心肌梗死或心律失常等疾病,我們需要準(zhǔn)確提取的特征有P 波,QRS 波群,T 波的振幅,寬度,斜率等,以及PQ 間期,PR 間期,ST 段等.對于心肌梗死,我們著重提取Q 波振幅和寬度,S-T 段抬高或壓低幅度,并將其轉(zhuǎn)換為正弦值,T 波是否倒置.

1.5 心肌梗死特征

心肌梗死又叫心肌梗塞,是冠狀動脈閉塞,血流中斷,使部分心肌因嚴(yán)重的持久性缺血而發(fā)生局部壞死.在ECG 上的表現(xiàn)大概可以分兩類：急性心肌梗死和非典型心肌梗死.

1.5.1 急性心肌梗死在ECG 的表現(xiàn)

1)異常Q 波(Q 波的深度大于等于R 波高度的四分之一).

2)S-T 段抬高或壓低或T 波倒置.(ST 段抬高超過0.1 mv 或壓低超過0.05 mv)[7].

1.5.2 非典型心肌梗死在ECG 的表現(xiàn)

1)II III avF 導(dǎo)聯(lián)QRS 波振幅降低且QRS 波起始出現(xiàn)頓挫或切跡[8].

2)R 波振幅對比下降明顯,出現(xiàn)頓挫或切跡.(主要為V1-V3 導(dǎo)聯(lián)).

3)P 波增寬,粗鈍,切跡成W 或M 型.

4)T 波高聳,T/R＞1.

1.6 神經(jīng)網(wǎng)絡(luò)構(gòu)建

針對心肌梗死程度不同,我們將其分為4 類.采用one-hot 編碼方式,如表1所示.

表1 訓(xùn)練數(shù)據(jù)統(tǒng)計表

對于非數(shù)值型特征值,例如T 波是否倒置,若T 波倒置,該特征值置1,否則,置-1.最后綜合所有特征值來確定心肌梗死的程度.初步提取的特征有：異常Q 波(Q 波振幅和寬度),ST 段抬高,ST 段壓低,T 波倒置,R 波振幅等.對于急性心肌梗死,采用BP 神經(jīng)網(wǎng)絡(luò),其中激活函數(shù)使用logsig,其函數(shù)原型為：

訓(xùn)練函數(shù)采用trainlm,其函數(shù)為L-M 算法,它同時具有梯度法和牛頓法的優(yōu)點(diǎn).當(dāng)λ很小時,步長等于牛頓法步長,當(dāng)λ很大時,步長約等于梯度下降法的步長.算法如下：

考慮如下信賴域模型：

其中,hk為信賴域半徑.該方程的解可有如下方程得到：

BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖9.

圖9 BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

對于非典型性心肌梗死,由于判斷條件復(fù)雜多變,在ECG 表現(xiàn)的特征不一,醫(yī)學(xué)上至今也未能給出明確的指標(biāo),為了更充分的提取特征,采用卷積神經(jīng)網(wǎng)絡(luò)：

微積分中卷積的表達(dá)式為：

離散形式：

用矩陣表示為(*表示卷積運(yùn)算)：

如果是二維的卷積,則表示式為：

其中,W為卷積核,X為輸入.如果X是一個二維輸入的矩陣,而W也是一個二維的矩陣.如果X是多維張量,那么W也是一個多維的張量[9].

本文卷積神經(jīng)網(wǎng)絡(luò)部分結(jié)構(gòu)如圖10所示[10,11].

圖10 卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

將波形按10 s 一次采樣,經(jīng)過平滑,濾波,去基線等前期處理后,送入卷積神經(jīng)網(wǎng)絡(luò)[12],分別經(jīng)過卷積層,dropout層,池化層,BN 層,如此循環(huán),最后經(jīng)過Flatten 層,全連接層輸出結(jié)果,本文采用1-D 卷積網(wǎng)絡(luò)[13].

2 實(shí)驗(yàn)結(jié)果

為了最小化模型結(jié)構(gòu)風(fēng)險,本文采用K折交叉驗(yàn)證的方式,隨機(jī)地將數(shù)據(jù)切分為K個大小相同但互不相交的子集,然后利用K-1 個子集訓(xùn)練數(shù)據(jù),剩余的子集測試模型,如此重復(fù)選擇,選出K次測試中誤差最小的模型[14].為了防止過擬合或欠擬合,在訓(xùn)練過程中加入L1,L2 混合使用的正則化.并且本文采用訓(xùn)練誤差和預(yù)測誤差隨著訓(xùn)練次數(shù)變化的動態(tài)曲線,當(dāng)訓(xùn)練誤差和預(yù)測誤差都在降低時,說明還處于欠擬合,如果預(yù)測誤差開始上升,說明處于過擬合階段,應(yīng)該停止訓(xùn)練.圖11為誤差隨訓(xùn)練次數(shù)變化的曲線圖.

圖11 訓(xùn)練效果圖

由圖可以看到,當(dāng)訓(xùn)練次數(shù)為700 次左右時,雖然訓(xùn)練誤差還在降低,但預(yù)測誤差已有上升趨勢,此時網(wǎng)絡(luò)泛化誤差最小,網(wǎng)絡(luò)訓(xùn)練效果最好.

本文分別使用MIT-BIH 和醫(yī)大一院數(shù)據(jù)進(jìn)行預(yù)測.其中,MIT-BIH 數(shù)據(jù)中沒用明確區(qū)分心肌梗死類型,所以我們只預(yù)測兩類.醫(yī)大一院根據(jù)醫(yī)生診斷結(jié)果分為四類.網(wǎng)絡(luò)預(yù)測的統(tǒng)計結(jié)果如表2、表3所示.

表2 MIT-BIH 數(shù)據(jù)預(yù)測結(jié)果統(tǒng)計表

表3 醫(yī)大一數(shù)據(jù)預(yù)測結(jié)果統(tǒng)計表

由預(yù)測結(jié)果可以看出,對于急性心梗和正常人群預(yù)測效果理想,因?yàn)槠涮卣鞑顒e較大,易于區(qū)分.但對于非典型性心肌梗死,由于其特征不明顯,變化較多,所以正確率有所降低.

測試相同PTB 數(shù)據(jù)集,與其他方法的效果對比見表4.

表4 不同算法對比結(jié)果統(tǒng)計表

由表3可以看出,本文算法(CNN＆BP)相對于其他方法,正確率有所提高,主要原因是在非典型性心肌梗死的檢測上,明顯優(yōu)于其他方法.

為了更好的評估模型,本文采用精準(zhǔn)率和召回率,F值衡量模型預(yù)測結(jié)果,并做出RoC 曲線圖,也有采用MAE 評價模型.定義如下：

TP：預(yù)測為正,實(shí)際為正；

FP：預(yù)測為正,實(shí)際為負(fù)；

TN：預(yù)測為負(fù),實(shí)際為負(fù)；

FN：預(yù)測為負(fù),實(shí)際為正.

精準(zhǔn)率P為：

召回率R為：

F值為：

由上述公式可知,精準(zhǔn)率描述預(yù)測為真,并且預(yù)測結(jié)果對的概率.召回率描述預(yù)測為真占所有真集的概率.而F值是精準(zhǔn)率和召回率的調(diào)和平均,只有當(dāng)二者都高時,F才會大.以TPR為y軸,以FPR為x軸,我們就直接得到了RoC 曲線,TPR和FPR公式如下：

從FPR和TPR的定義可以理解,TPR越高,FPR越小,模型和算法就越高效,或者說在RoC 曲線下所圍面積越大,算法效果越好,圖12為幾種算法的RoC 曲線：

圖12 RoC 曲線

由上圖RoC 曲線可知,CNN＆BP 算法在精準(zhǔn)率和召回率上都優(yōu)于其他算法.

3 總結(jié)

本文首先對ECG 信號做濾波,去基線等處理,根據(jù)波形特點(diǎn)在準(zhǔn)確識別QRS 波群后,提取相應(yīng)特征.針對心肌梗死不同的類型,采用BP 神經(jīng)網(wǎng)絡(luò)和卷積神經(jīng)網(wǎng)絡(luò)[15]相結(jié)合的方法,并對12 導(dǎo)聯(lián)做分析.實(shí)驗(yàn)結(jié)果表明,針對不同心肌梗死類型采用不同神經(jīng)網(wǎng)絡(luò)預(yù)測在準(zhǔn)確率上提高2%,但由于非典型心肌梗死特征多樣且不明顯,所以對此仍有待改進(jìn)和完善,文獻(xiàn)[16]提出可以根據(jù)12 導(dǎo)聯(lián)確定心肌梗死具體起源位置.雖然心電圖自動分析程序已經(jīng)達(dá)到了較高的準(zhǔn)確率,但是由于疾病種類繁多,臨床表現(xiàn)不一,給診斷帶來一定困難,仍不能完全取代醫(yī)生,只是起到一個輔助診斷的作用.對于一些復(fù)雜情況,仍需要結(jié)合臨床診斷.