基于非參數相關系數的心肌病自動診斷

林卓琛，張晉昕

前言

2015年，由心血管疾病所致的死亡人數占全球死亡人數的31%。2016年，世界衛生組織將其確定為全球居民主要死因，2001年至2015年間每年約有1 750 萬人死于該病［1］。這個數字隨著世界人口的年齡增長，在2030年預計將增加到2 220 萬［2］。調查顯示，2010年由心血管疾病引起的全球直接醫療費用約為8 630 億美元，并且在2030年累計將達到20 萬億美元［3］。而在中國，心臟病的患病率和死亡率一直處于不斷上升狀態，及早發現心肌病在內的心臟疾病，顯得尤為重要［4］。

臨床醫生通常使用12 導聯心電圖來評估各種心血管疾病的發病風險，心電圖機是一種將電極放在皮膚上的工具，用于收集一段時間內心臟電活動的信息［5］。然而，人工研判這些心電圖是單調、耗時的，對于那些關鍵的、同時也是細微的診斷細節，單純通過人工評閱可能會有遺漏，甚至導致心電圖研判的錯誤［6］。因此，為了克服人工評估的局限性，越來越多的研究者開發了相應的計算機輔助自動診斷技術。

當前使用心電圖數據進行自動分類的關鍵在于信號特征的提取。有的研究者直接對原始ECG 數據進行特征提取，也有的研究者首先進行小波分析等處理得到小波系數后再進行分析，雖然小波分析存在小波基選擇等問題，其多分辨率分析的特點使其在時頻域上有效提取信號的局部特征，因此仍是當前研究者選擇的主要心電數據處理方法。

人們對于心電數據特征的提取主要是以統計特征為主，主要包含描述波形信息的特征如偏度、峰度［7］，描述能量信息的方差［8］，描述導聯信息量的熵［9］等。除統計特征外，有研究者嘗試使用多項式擬合系數［10］、分型系數［9,11］等作為信號的特征。然而，人們對特征的研究多集中于單個導聯提取的特征，對于導聯間的相關性開展深入挖掘的力度不足。Maharaj 等［12］曾在2014年研究了小波分析系數的導聯間相似性對心肌梗死的影響，但是沒有考慮到小波系數在各個分解層級中可能正態性欠佳，也沒有考慮到導聯的系數間不是線性關系。本研究提出一種基于非參數的相關系數特征，穩健地刻畫了導聯間的相關性，對分解后系數不作正態性假設。本研究的分析流程圖見圖1，本研究采用Matlab2018b 進行數據分析。

圖1 本研究分析流程圖Fig.1 Flow diagram of research design

1 特征提取

1.1 數據來源

德國國家計量署的PTB數據庫來自柏林本杰明富蘭克林（Benjamin Franklin）醫院心臟內科，常用于心臟疾病相關算法的研究，是本研究的合適選擇。該數據庫包含148名心肌梗死患者、18名冠心病及心衰患者、15名束支傳導阻滯患者、14名心率失常患者、52名健康人以及21名其余心臟相關疾病患者。每個記錄采用并列儲存方式，包含15列，分別為標準12導聯（I,II,III,avR,avL,avF,V1,V2,V3,V4,V5,V6）和3富蘭克林導聯（vx,vy,vz），值域范圍為-1 000～3 000 mV，表示該時間點下心電信號的強度。每個記錄時長約2 min（120 000個采樣點/行），采樣率為1 000 Hz，分辨率為16 bit（±16.384 mV），儲存格式為Format 16。本研究使用其中的標準12導聯信號進行分析。本文旨在考察鑒別診斷方法的診斷能力，與其他報道進行對比。所以，也按照已有報道的做法［13-14］，納入數據庫中所有的健康人、肥厚型心肌病患者和擴張型心肌病患者。納入患者以臨床醫生的診斷結果為準劃分組別，排除掉合并其他心臟系統疾病的患者，最后選擇了來自PTB數據庫的52名健康人、7名肥厚型心肌病患者和9名擴張型心肌病患者進行分析研究。

1.2 數據預處理與心拍分割

為了消除心電信號帶有的各類型噪聲干擾（工頻干擾、基線漂移和高頻干擾等）［15］，本研究使用Daubechies 6（db6）小波基函數對心電信號進行去噪和消除基線漂移［16］。與主流的分析方法相一致，本研究用0.351～125.000 Hz 帶通濾波進行過濾，隨后使用5 到11 層的細節系數重構信號。預處理后的信號高頻噪聲被過濾，局部變得光滑，信號整體處于同一水平。

去噪后的整體信號需要根據獨立的心動周期分割為單獨的心電片段，本研究使用Pan-Tompkin 算法進行R 峰的識別來達到心拍分割的目的［17］。Pan-Tompkin 算法是一個穩健的R 峰識別算法，計算簡單，易于實現。將R 峰前的255 個采樣點、R 峰后的256 個采樣點以及R 峰本身作為一個單獨的片段進行后續分析（共512個采樣點），最終每個人的心電信號將被分割為80～160個完整的心拍信號。為了確保片段的長度符合標準，我們刪除了每個人的第一個片段和最后一個片段，使得所有的片段具有相同的長度。最終從7名肥厚型心肌病患者中提取片段891段，9名擴張型心肌病患者中提取片段1 196段，52名健康人中提取片段7 099段。

1.3 小波分析

獲取心電圖片段后，對該片段512 個采樣點進行歸一化處理，使其總方差等于1。處理后不同碎片之間總能量的差異被消除，有利于后續挖掘有意義的特征。我們使用極大重疊離散小波變換（MODWT）來分解原始信號［18］。MODWT 是對離散小波變換的一種改進，它不使用下采樣過程，因此所有尺度上的小波系數都具有與原始時間序列長度相等的長度，在處理相關問題上MODWT 的這個特性會有一定的優勢。每個片段經過以dmey 為小波基的MODWT后被分解為7個層級［19］。我們只使用了5到7層的細節系數以及第7 層的逼近系數。這4 個系數序列包含0.351～62.500 Hz 頻帶范圍，涵蓋了原始信號幾乎所有的信息。每一個小波系數序列包含512 個值。關于ECG 經過MODWT 產生4 個小波系數序列的一個示例見圖2。

圖2 基于dmey小波的原始信號和分解信號Fig.2 Raw signal and decomposed signals based on dmey wavelet

1.4 非參數相關系數的計算

經過小波分解后，每個導聯可計算出4個小波系數。對于同一種小波系數，我們計算12 個導聯兩兩間的非參數相關系數。每一種小波系數計算出66個特征，共計264個特征。HoeffdingD測度的直觀思想是通過計算零假設（兩變量獨立）下邊際分布乘積與經驗二元分布之間的距離來檢驗數據集的獨立性［20］。這兩個分布之間的差異越大，D就越大，導聯間相關性就越大。HoeffdingD測度有助于識別非線性、非單調關聯。其計算公式如下：

其中，記原始信號X 在MODWT 分解后的第j條系數序列為Wx,j,·，原始信號Y 在MODWT 分解后的第j條系數序列為Wy,j,·，本研究每個導聯使用4 條系數序列。Ri是Wx,j,i在系數序列Wx,j,·由小到大排序后得到的秩。而Si是Wy,j,i在系數序列Wy,j,·由小到大排序后得到的秩。將序列Wx,j,·與Wy,j,·結成的n對點對(Wx,j,t,Wy,j,t)(t= 1,…,n)，Qi是點對中的兩者同時小于(Wx,j,t,Wy,j,t)的點對個數。

通過計算、整理，最終每個心拍片段計算得到264 個HoeffdingD測度，組成9 186×264 的矩陣。部分計算結果見圖3。其中結局標簽的含義為1：健康人，2：擴張型心肌病，3：肥厚型心肌病；受試者id 為每個受試者的編號，片段id 為同一個受試者內心拍片段的編號；緊接的264 列為HoeffdingD測度，編碼形式為x_y_z，其中x 的含義為5：第5 層細節系數，6：第6層細節系數，7：第7層細節系數，8：第7層的逼近系數；y 與z 表示第y 導聯與第z 導聯，第1 至12 導聯依次對應導聯I,II,III,avR,avL,avF,V1,V2,V3,V4,V5,V6。

圖3 Hoeffding D測度計算結果（部分）Fig.3 Calculated results of Hoeffding D measure(part of contents)

1.5 特征選擇

得到的264個特征中含有對疾病預測有價值的特征，需要通過變量選擇的方式挑選出來。對于每個特征，本研究通過方差分析（ANOVA）計算其關于3個結局的F值。F值越大，說明特征越重要。使用10折交叉驗證分別計算F閾值為1 100～4 900時的總準確率均值以及標準差，具體見圖4。隨著F值增大，在抵達輔助線（縱向虛線F=3 600處）之前，被保留下來特征數在減少，總正確率并未減低；當F值進一步增大，越過輔助線后，總準確率均值開始下降，且標準差開始增大。可見3 600是比較合適的F閾值，此處對應的預測正確率高，且穩健。

根據變量的F值對變量進行排序，挑選出較重要的特征。選擇F值在3 600 以上的特征，共挑選出12個特征，特征的情況見表1，Hr3,65表示展開的第5層第3 導聯和第6 導聯的HoeffdingD測度，其余符號按相同規則標記。

2 建模與結果分析

2.1 建模與評價介紹

分類任務涉及訓練集和測試集的劃分，通過將訓練集輸入到支持向量機中建立預測模型，并使用測試集評價模型的優劣。本研究使用10折交叉驗證評價模型效果，通過隨機分組，將計算得到的9 186個心拍片段平均分成10 份，每次取出其中1 份作為測試集，余下9份合并在一起作為訓練集。使用訓練集進行建模后，對預測集進行預測，得到1 次預測結果。這個過程重復10次，計算出其均數和標準差，可對模型預測能力做出準確而穩健的評價。

表1 前12個重要特征在3組間的情況（± s）Tab.1 The result of the first 12 important features among 3 groups(Mean±SD)

表1 前12個重要特征在3組間的情況（± s）Tab.1 The result of the first 12 important features among 3 groups(Mean±SD)

特征Hr3,6 5 Hr3,8 5 Hr3,9 5 Hr6,10 5 Hr1,3 7 Hr3,8 8 Hr4,10 8 Hr4,12 8 Hr5,7 8 Hr5,8 8 Hr5,9 8 Hr7,10 8肥厚型心肌病0.613±0.268 0.406±0.215 0.378±0.169 0.475±0.182 0.344±0.217 0.658±0.195 0.399±0.161 0.543±0.242 0.440±0.241 0.441±0.242 0.564±0.267 0.362±0.277擴張型心肌病0.583±0.236 0.240±0.148 0.411±0.180 0.242±0.223 0.211±0.132 0.619±0.194 0.132±0.169 0.281±0.223 0.231±0.212 0.383±0.209 0.566±0.273 0.573±0.216健康人0.217±0.151 0.066±0.080 0.085±0.108 0.705±0.130 0.052±0.060 0.252±0.154 0.056±0.072 0.169±0.083 0.051±0.057 0.079±0.090 0.158±0.070 0.166±0.100 F值3 666.829 4 514.751 4 947.168 5 322.891 4 361.455 3 856.695 3 905.770 3 751.385 4 970.278 4 990.929 6 571.701 4 328.031

2.2 支持向量機

支持向量機是定義在特征空間上使得類間間隔最大化的線性分類器，通過引入選定的核函數，將輸入空間的特征映射到一個高維的特征向量空間中，并在該空間中構造最優分類平面，實現健康信號與兩種心肌病的分類。本研究采用徑向基函數作為核函數，對于需要選定的超參數c和γ，使用網格搜索法得到理想值，c等于2.83，γ等于4。

2.3 分類結果

以真實疾病結局為橫行，預測疾病結局為縱列，對每個人預測情況進行計數，統計出混淆矩陣，結果見表2。采用總體準確率（Acc）以及各類結局的靈敏度（SenHC,SenDCM,SenHCM）和特異度（SpecHC,SpecDCM,SpecHCM）評價模型結果。指標的計算公式如下：

其中，HC、DCM、HCM 分別代表健康人、擴張型心肌病以及肥厚型心肌病。TPi表示屬于第i類結局且被正確預測的心拍個數；FPi表示不屬于第i類結局但被錯誤地預測為第i類結局的心拍個數；TNi表示不屬于第i類結局的心拍且被預測為非第i類結局的心拍個數；FNi表示屬于第i類結局的心拍但沒有被預測為第i類結局的心拍個數。總準確率公式中分子為3類結局預測正確的心拍個數，分母為總心拍個數。靈敏度公式中分子為第i類正確預測的心拍個數，分母為真是患有第i類結局疾病的心拍個數。特異度公式中分子為不屬于第i類結局的心拍且被預測為非第i類結局的心拍個數，分母為不屬于第i類的心拍個數。由于采用10 折交叉驗證，結果得到了10 個準確率、靈敏度和特異度，分別計算其均數與標準差，結果見表3。

表2 分類結果的混淆矩陣Tab.2 Confusion matrix of classification results

表3 模型分類的評價指標（%）Tab.3 Evaluation measurements to the classifications(%)

當前對于心肌病自動診斷的研究仍較少，且多為健康人與單個心肌疾病之間的鑒別診斷。近幾年對心肌病的研究結果見表4，由表可以看出，本研究提出的特征計算簡單、方便，同時對兩種心肌病進行預測且總準確率均高于同類研究。

表4 與近期其他文獻結果比較Tab.4 Comparison with the results from other recent studies

3 結論

基于心電圖導聯間的相關性，本研究提出使用HoeffdingD測度對兩種心肌病進行自動診斷的算法。首先使用MODWT 對預處理后的心電圖信號進行分解，隨后計算4 層12 個導聯間小波系數的HoeffdingD測度，通過F檢驗提取最重要的12 個特征，使用支持向量機分類器實現對兩種心肌病與健康心拍的自動分類。本文最終識別正確率達到99.88%，優于同類型研究，可用于輔助臨床心肌病的預測與診斷。

本研究在當前階段使用國際公開數據庫進行建模分析，尚未在當地醫院進行數據收集以及相關的模型檢驗。為結合臨床的實際應用需求，本研究未來計劃將與醫療機構合作，結合中國人病例數據分析兩者篩選所得特征的一致程度，進一步分析模型的泛化能力與實際應用效果。