基于非參數(shù)相關(guān)系數(shù)的心肌病自動(dòng)診斷

林卓琛，張晉昕

前言

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

臨床醫(yī)生通常使用12 導(dǎo)聯(lián)心電圖來評(píng)估各種心血管疾病的發(fā)病風(fēng)險(xiǎn)，心電圖機(jī)是一種將電極放在皮膚上的工具，用于收集一段時(shí)間內(nèi)心臟電活動(dòng)的信息［5］。然而，人工研判這些心電圖是單調(diào)、耗時(shí)的，對(duì)于那些關(guān)鍵的、同時(shí)也是細(xì)微的診斷細(xì)節(jié)，單純通過人工評(píng)閱可能會(huì)有遺漏，甚至導(dǎo)致心電圖研判的錯(cuò)誤［6］。因此，為了克服人工評(píng)估的局限性，越來越多的研究者開發(fā)了相應(yīng)的計(jì)算機(jī)輔助自動(dòng)診斷技術(shù)。

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

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

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

1 特征提取

1.1 數(shù)據(jù)來源

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

1.2 數(shù)據(jù)預(yù)處理與心拍分割

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

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

1.3 小波分析

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

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

1.4 非參數(shù)相關(guān)系數(shù)的計(jì)算

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

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

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

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

1.5 特征選擇

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

根據(jù)變量的F值對(duì)變量進(jìn)行排序，挑選出較重要的特征。選擇F值在3 600 以上的特征，共挑選出12個(gè)特征，特征的情況見表1，Hr3,65表示展開的第5層第3 導(dǎo)聯(lián)和第6 導(dǎo)聯(lián)的HoeffdingD測度，其余符號(hào)按相同規(guī)則標(biāo)記。

2 建模與結(jié)果分析

2.1 建模與評(píng)價(jià)介紹

分類任務(wù)涉及訓(xùn)練集和測試集的劃分，通過將訓(xùn)練集輸入到支持向量機(jī)中建立預(yù)測模型，并使用測試集評(píng)價(jià)模型的優(yōu)劣。本研究使用10折交叉驗(yàn)證評(píng)價(jià)模型效果，通過隨機(jī)分組，將計(jì)算得到的9 186個(gè)心拍片段平均分成10 份，每次取出其中1 份作為測試集，余下9份合并在一起作為訓(xùn)練集。使用訓(xùn)練集進(jìn)行建模后，對(duì)預(yù)測集進(jìn)行預(yù)測，得到1 次預(yù)測結(jié)果。這個(gè)過程重復(fù)10次，計(jì)算出其均數(shù)和標(biāo)準(zhǔn)差，可對(duì)模型預(yù)測能力做出準(zhǔn)確而穩(wěn)健的評(píng)價(jià)。

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

表1 前12個(gè)重要特征在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擴(kuò)張型心肌病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 支持向量機(jī)

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

2.3 分類結(jié)果

以真實(shí)疾病結(jié)局為橫行，預(yù)測疾病結(jié)局為縱列，對(duì)每個(gè)人預(yù)測情況進(jìn)行計(jì)數(shù)，統(tǒng)計(jì)出混淆矩陣，結(jié)果見表2。采用總體準(zhǔn)確率（Acc）以及各類結(jié)局的靈敏度（SenHC,SenDCM,SenHCM）和特異度（SpecHC,SpecDCM,SpecHCM）評(píng)價(jià)模型結(jié)果。指標(biāo)的計(jì)算公式如下：

其中，HC、DCM、HCM 分別代表健康人、擴(kuò)張型心肌病以及肥厚型心肌病。TPi表示屬于第i類結(jié)局且被正確預(yù)測的心拍個(gè)數(shù)；FPi表示不屬于第i類結(jié)局但被錯(cuò)誤地預(yù)測為第i類結(jié)局的心拍個(gè)數(shù)；TNi表示不屬于第i類結(jié)局的心拍且被預(yù)測為非第i類結(jié)局的心拍個(gè)數(shù)；FNi表示屬于第i類結(jié)局的心拍但沒有被預(yù)測為第i類結(jié)局的心拍個(gè)數(shù)?？倻?zhǔn)確率公式中分子為3類結(jié)局預(yù)測正確的心拍個(gè)數(shù)，分母為總心拍個(gè)數(shù)。靈敏度公式中分子為第i類正確預(yù)測的心拍個(gè)數(shù)，分母為真是患有第i類結(jié)局疾病的心拍個(gè)數(shù)。特異度公式中分子為不屬于第i類結(jié)局的心拍且被預(yù)測為非第i類結(jié)局的心拍個(gè)數(shù)，分母為不屬于第i類的心拍個(gè)數(shù)。由于采用10 折交叉驗(yàn)證，結(jié)果得到了10 個(gè)準(zhǔn)確率、靈敏度和特異度，分別計(jì)算其均數(shù)與標(biāo)準(zhǔn)差，結(jié)果見表3。

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

表3 模型分類的評(píng)價(jià)指標(biāo)（%）Tab.3 Evaluation measurements to the classifications(%)

當(dāng)前對(duì)于心肌病自動(dòng)診斷的研究仍較少，且多為健康人與單個(gè)心肌疾病之間的鑒別診斷。近幾年對(duì)心肌病的研究結(jié)果見表4，由表可以看出，本研究提出的特征計(jì)算簡單、方便，同時(shí)對(duì)兩種心肌病進(jìn)行預(yù)測且總準(zhǔn)確率均高于同類研究。

表4 與近期其他文獻(xiàn)結(jié)果比較Tab.4 Comparison with the results from other recent studies

3 結(jié)論

基于心電圖導(dǎo)聯(lián)間的相關(guān)性，本研究提出使用HoeffdingD測度對(duì)兩種心肌病進(jìn)行自動(dòng)診斷的算法。首先使用MODWT 對(duì)預(yù)處理后的心電圖信號(hào)進(jìn)行分解，隨后計(jì)算4 層12 個(gè)導(dǎo)聯(lián)間小波系數(shù)的HoeffdingD測度，通過F檢驗(yàn)提取最重要的12 個(gè)特征，使用支持向量機(jī)分類器實(shí)現(xiàn)對(duì)兩種心肌病與健康心拍的自動(dòng)分類。本文最終識(shí)別正確率達(dá)到99.88%，優(yōu)于同類型研究，可用于輔助臨床心肌病的預(yù)測與診斷。

本研究在當(dāng)前階段使用國際公開數(shù)據(jù)庫進(jìn)行建模分析，尚未在當(dāng)?shù)蒯t(yī)院進(jìn)行數(shù)據(jù)收集以及相關(guān)的模型檢驗(yàn)。為結(jié)合臨床的實(shí)際應(yīng)用需求，本研究未來計(jì)劃將與醫(yī)療機(jī)構(gòu)合作，結(jié)合中國人病例數(shù)據(jù)分析兩者篩選所得特征的一致程度，進(jìn)一步分析模型的泛化能力與實(shí)際應(yīng)用效果。