機器學習實現心肌梗死的自動檢測①

常戰國,蒲寶明,李相澤,王 帥,楊 朔

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

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

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

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

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

圖1 算法流程圖

1 研究方法

本文主要從數據獲取,數據前期處理,波形檢測,特征提取,構建神經網絡及識別等幾個方面敘述.

1.1 數據獲取

數據主要來源于MIT-BIH 數據庫和中國醫科大學附屬第一醫院.

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

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

1.2 ECG 數據前期處理

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

1.2.1 濾波

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

(1)均值濾波如下：

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

(2)Butterworth 濾波如下：

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

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

由極點在單位圓上可以得到：

因此傳遞函數為：

1到10 階的Butterworth 多項式因子可以通過查表得到.

圖2-圖4為兩種方法的對比結.

圖2 原始ECG 信號

圖3 Smooth 平滑后ECG

圖4 Butterworth 濾波后ECG

從圖中可以看出,經過均值濾波后,ECG 波形基本沒有發生形變,只是隨著平滑點數的增加,幅值有所降低,為此采用一種補償方法,將其補償到原來的數值.但是經過Butterworth 濾波后,除了波形幅值有所降低外,在波形拐點處會有一定的形變,尤其在S 波會有一個明顯向下的尖小的波,因此導致波形在一定程度上有所失真.

1.2.2 去基線

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

1.3 QRS 波群檢測

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

圖5 濾波后ECG

圖6 一階導ECG

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

圖7 二階導ECG

1.4 ECG 特征提取

在準確識別了QRS 波群的前提下,需要對ECG波形做進一步的特征提取.國際上公認的常規12 導聯是標準I,II,III,加壓單極肢體導聯aVR,aVL,aVF 和單極胸腔壁導聯V1-V6[6].12 導聯之間存在如下向量運算關系：

每一個單極導聯都由周期性的P,QRS,T 波等組成.個別還有U 波,但目前沒有發現明確的生理學意義(U 波異?？赡芘c缺鉀有關,但不是絕對的),我們暫不研究.一個完整心動周期如圖8所示.

圖8 心動周期示意圖

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

1.5 心肌梗死特征

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

1.5.1 急性心肌梗死在ECG 的表現

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

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

1.5.2 非典型心肌梗死在ECG 的表現

1)II III avF 導聯QRS 波振幅降低且QRS 波起始出現頓挫或切跡[8].

2)R 波振幅對比下降明顯,出現頓挫或切跡.(主要為V1-V3 導聯).

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

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

1.6 神經網絡構建

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

表1 訓練數據統計表

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

訓練函數采用trainlm,其函數為L-M 算法,它同時具有梯度法和牛頓法的優點.當λ很小時,步長等于牛頓法步長,當λ很大時,步長約等于梯度下降法的步長.算法如下：

考慮如下信賴域模型：

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

BP 神經網絡結構如圖9.

圖9 BP 神經網絡結構

對于非典型性心肌梗死,由于判斷條件復雜多變,在ECG 表現的特征不一,醫學上至今也未能給出明確的指標,為了更充分的提取特征,采用卷積神經網絡：

微積分中卷積的表達式為：

離散形式：

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

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

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

本文卷積神經網絡部分結構如圖10所示[10,11].

圖10 卷積神經網絡結構

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

2 實驗結果

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

圖11 訓練效果圖

由圖可以看到,當訓練次數為700 次左右時,雖然訓練誤差還在降低,但預測誤差已有上升趨勢,此時網絡泛化誤差最小,網絡訓練效果最好.

本文分別使用MIT-BIH 和醫大一院數據進行預測.其中,MIT-BIH 數據中沒用明確區分心肌梗死類型,所以我們只預測兩類.醫大一院根據醫生診斷結果分為四類.網絡預測的統計結果如表2、表3所示.

表2 MIT-BIH 數據預測結果統計表

表3 醫大一數據預測結果統計表

由預測結果可以看出,對于急性心梗和正常人群預測效果理想,因為其特征差別較大,易于區分.但對于非典型性心肌梗死,由于其特征不明顯,變化較多,所以正確率有所降低.

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

表4 不同算法對比結果統計表

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

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

TP：預測為正,實際為正；

FP：預測為正,實際為負；

TN：預測為負,實際為負；

FN：預測為負,實際為正.

精準率P為：

召回率R為：

F值為：

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

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

圖12 RoC 曲線

由上圖RoC 曲線可知,CNN＆BP 算法在精準率和召回率上都優于其他算法.

3 總結

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