王 迪,李 強,孫亞楠,韋杰英,張 瑞
(西北大學 醫(yī)學大數據研究中心, 陜西 西安 710127)
心房顫動(atrial fibrillation, AF),簡稱房顫,是臨床上最為常見的心律失常之一,往往會導致嚴重的并發(fā)癥(如心力衰竭和中風等)[1]。心電圖(electrocardiogram,ECG)是由心電圖儀記錄反映心臟興奮的電活動過程,對房顫的診斷與治療具有重要的參考價值[2]。圖1顯示一個正常心動周期內的心電圖波形,其中P波表示心房去極化過程,QRS復合波表示心室去極化過程,T波表示心室復極化過程。房顫發(fā)作時,心電圖的改變主要表現在兩個方面:①RR間期絕對不規(guī)則;②P波缺失,代之以大小不等、間距不均、形態(tài)各異的f波。圖2顯示了正常心律和心房顫動的心電圖片段。
圖1 一個正常心動周期的心電圖波形Fig.1 ECG of a normal cycle
圖2 竇性心律與房顫心電Fig.2 Sinus rhythm and atrial fibrillation
根據房顫發(fā)作持續(xù)時長,可分為陣發(fā)性房顫、持續(xù)性房顫和永久性房顫。其中陣發(fā)性房顫發(fā)作突然且持續(xù)時間短暫,臨床上易被漏診,因此往往需要長時程的心電監(jiān)測。傳統的房顫檢測主要是醫(yī)生根據經驗進行視覺診斷,但長時程的心電數據使得這一過程較為耗時。此外,僅依靠視覺觀察往往難以捕捉陣發(fā)性房顫在心電圖上所表現出的短暫且微弱的波形變化,從而造成較多漏診。因此,研究陣發(fā)性房顫的自動檢測具有十分重要的臨床意義。
關于陣發(fā)性房顫自動檢測的研究可大體分為兩類:基于RR間期不規(guī)則的方法和基于P波缺失的方法。1983年Moody G.B.采用馬爾可夫過程模型對RR間期序列進行分解并得到RR間期預測數組,當該數組的平均誤差超過某一給定的閾值時,則認為對應于一個房顫片段[3]。Duverney D. 利用小波多尺度分析技術將RR間期時間序列分解成不同的頻域水平,進而突出高頻段的房顫起點與終點,最后利用分形法來分類竇性心律與房顫[4]。Tateno K.等提取了RR間期和RR間期差的房顫標準密度直方圖,并采用柯爾莫哥洛夫史密諾夫檢驗評估待測密度直方圖與標準直方圖的差異度,如果差異度小于某一給定的閾值,則認為該段信號為房顫心電[5]。然而,值得注意的是,RR間期不規(guī)則并非唯一對應于房顫,這一現象也會出現在其他心律失常之中(如房撲、多源性房速等)。此外,當房顫發(fā)作時,如果出現房室傳導阻滯等心臟異常情況時,其心電圖中RR間期也可能是規(guī)則的。因而,僅依賴房顫不規(guī)則這一特征顯然不能全面刻畫房顫?;诖?結合心房活動(即P波特征)進行研究分析可大大提高房顫檢測的準確率[6]。已有的基于P波特征的房顫檢測方法主要分為兩類:時域法和頻域法。時域法主要是對P波間期、P波離散度、P波變異性等概念進行定量分析。此類方法對P波形態(tài)的依賴性較小,而且重復性較好。其中P波離散度最早由Dilavefis P.E.等人提出,是指在常規(guī)的12導聯心電圖中最大與最小P波間期之差,其基準值一般設定為4s[7]。P波變異性[8]由Andrikopoulos G.K.等人中提出,將其定義為P波間期標準差的平方。該參數類似于P波離散度,能夠反映心房活動在傳導過程中的變異性。頻域法則是通過對心電圖中的P波或f波作快速傅里葉變換,比較分析房顫發(fā)作時各個頻域段的特征進行房顫檢測。文獻[9]和[10]利用傅里葉變換對信號時域和頻域之間進行全域變換,進而對信號進行相干分析后得到時間段內的相干系數。Ortigosa等人對心電信號進行廣義傅里葉變換,提取心房活動特征并結合支撐向量機完成陣發(fā)性房顫的自動檢測[11]。白鵬飛等人利用小波包變換處理濾波后的心電信號,同時重構出P波主頻帶波形并確定P波的始末位置,最后將所提取P波的形態(tài)學特征與人工神經網絡結合完成陣發(fā)性房顫的自動檢測[12]。
為了充分運用心電信號的時頻局部特性及瞬時特性,本文提出一種基于小波相干分析的陣發(fā)性房顫自動檢測方法。首先,對所有的心電信號進行預處理,包括去除噪聲、分段處理、選取模板等;其次,對模板與待測心電信號進行小波相干分析得到小波相干圖;進而,計算小波相干值的均值、比率以及交叉小波相位角方差,構成房顫心電特征;最后,將所提取特征結合超限學習機完成陣發(fā)性房顫的自動檢測。
陣發(fā)性房顫自動檢測的本質就是一個兩類分類問題。本小節(jié)將系統介紹相關基礎理論知識、心電信號的預處理、房顫心電特征提取方法以及超限學習機。
1.1.1 連續(xù)小波變換 連續(xù)小波變換(continuous wavelet transform, CWT)是一種時間窗和頻率窗都可改變的局部化分析方法[13],可以在預先定義的尺度上對原始信號進行多尺度濾波,能夠對非平穩(wěn)的心電信號由粗及細地進行局部化分析[14-15]。
給定一段心電信號x(t)∈L2(R),則該心電信號的連續(xù)小波變換可以定義為
(1)
式(1)中:a為尺度因子,b為平移因子,Ψa,b(t)為小波基函數。通過對信號進行連續(xù)小波變換,可以得到信號在不同尺度下對應帶通的信息。
1.1.2 小波相干分析 小波相干分析[16]是將小波變換和相干分析相結合的一種分析方法,可以分析兩列信號的相互依賴關系尤其可用于刻畫模板與待測心電信號在不同頻率點上的相關程度。
任意給定兩段長度為N的信號x和y,分別實施連續(xù)小波變換得到其小波系數為CWTx(a,b)與CWTy(a,b)。進一步,采用交叉小波變換得到交叉小波譜WCSx,y(a,b)為
(2)
其中CWTy*(a,b)表示CWTy(a,b)的復共軛。基于此,可計算交叉小波相位角為
(3)
其中I和R分別代表交叉小波譜的虛部和實部,φxy∈[-π,π]。
為了更好地得到信號之間的同步信息,在計算小波相干值之前,我們需要對小波譜做平滑處理。定義平滑函數為
S(W)=Sa[St(W)]
(4)
其中Sa表示對尺度軸做平滑,St表示對時間軸作平滑,分別表示為
(5)
Sa(CWTx(a,b))=CWTx(a,b)*c2Π(0.6a)。
(6)
在式(5)中,c1是標準化系數,*代表卷積運算。在式(6)中,c2是標準化系數,Π是矩形函數,參數0.6是根據經驗確定的尺度[17]。
在上述工作的基礎上,信號x和y的小波相干值(wavelet coherence,WTC)定義為
WTCxy(a,b)=
(7)
1.2.1 預處理 心電信號是一種非平穩(wěn)信號,在采集過程中容易受到各種背景噪聲的影響(主要包括工頻干擾、肌電干擾和基線漂移等)。本文首先采用FIR數字濾波器[18]去除心電信號中的基線漂移;其次,采用小波變換[19]對原始心電信號進行分解并對高頻系數進行軟閾值處理;最后通過小波重構得到去噪后干凈的心電信號。圖3顯示了一段心電信號(時長10s)去噪前后的效果。將去噪之后的心電信號進行無重疊的分段處理,分段后每個心電片段長度為10s(2500個采樣點),并從數據集中隨機選取一段10s正常心電片段作為模板。
圖3 心電信號去噪效果Fig.3 Thedenoising performance of ECG
1.2.2 特征提取 給定模板心電片段Y={y(1),y(2),…,y(m)},并設X={x(1),x(2),…,x(n)}為待測心電信號,其中m,n表示信號長度且m 房顫發(fā)作時P波缺失而代之以f波(房顫波),且P波與f波重疊的頻帶范圍一致。因此,可以提取并分析P波所在尺度上的小波相干圖,如圖5所示。其中,圖5(a)是模板與正常心電信號的小波相干圖,圖5(b)是模板與異常心電信號的小波相干圖。圖中色標代表小波相干值,黑色箭頭表示交叉小波相位角,這里僅標出了WTC≥0.5時的交叉小波相位角。從圖5可以看出,相比于正常心電與模板匹配的小波相干圖,房顫心電與模板匹配的小波相干值整體接近于0,且交叉小波相位角分布較為分散?;诖?本文采用3個度量指標來刻畫上述差異,并將其作為所提取的房顫心電特征。具體步驟可總結為如下算法。 算法1(房顫心電特征提取方法) 給定心電信號X。 步驟2采用小波多尺度分析技術分別對Y和Xi進行連續(xù)小波變換,得到各自r到t層的低頻小波系數,并根據公式(2)和(3)分別計算Y和Xi之間的交叉小波譜WCSXi,Y(a,b)和交叉小波相位角φXiY。記交叉小波相位角為 步驟3利用平滑函數S對交叉小波譜做平滑處理,進而得到小波相干值WTCXiY(a,b),記為 步驟4分別以樣本點與周期為x軸與y軸,小波相干值為坐標(x,y)處的值,繪制模板與待測心電信號的小波相干圖。 為心電信號X的房顫心電特征。 1.2.3 超限學習機 超限學習機(extreme learning machine,ELM)是2006年由Huang等人提出的一種新的學習算法[20]。原始ELM算法針對單隱層前饋神經網絡,輸入層與隱層的連接權重和閾值隨機選取,采用可微或不可微的激活函數。與其他學習算法相比,ELM學習速度快、泛化精度高,且不易陷入局部最小值。 圖4 待測信號與模板匹配對應的小波相干圖Fig.4 The wavelet coherence map sbetween the template EEG and the present EEG 圖5 特定尺度下待測信號與模板匹配對應的小波相干圖Fig.5 The wavelet coherence maps between the template EEG and the present EEGat specific scales (8) 這里wj=[w1j,w2j,…,wnj]T是連接輸入層與第j個隱單元的權值向量,bj是第j個隱單元的閾值,βj=[βj1,βj2,…,βjm]T是連接第j個隱單元與輸出層的權值向量;wj·xi表示wj與xi的內積。假設網絡中的實際輸出與期望輸出相同,則網絡的學習過程即為求解到如下線性方程組: Hβ=T。 (9) 這里,H,β,T的定義如下: H= (10) β=[β1,β2,…,βL]T,T=[t1,t2,…,tN]T。 (11) 本文采用MIT-BIH心房顫動數據集[6],該數據集包括25個時長為1h的心電圖記錄,采樣率250Hz,分辨率12-bit。本文從中隨機選取8個心電記錄進行數值實驗,所有實驗均在Matlab R2016中運行。 數值實驗中,所有數據被隨機均分為訓練集和測試集。實驗執(zhí)行50次并取其平均結果作為最終分類性能的度量。 在特征提取過程中,首先對心電信號進行等長分段處理,每段心電數據長10s(2500個點),并做連續(xù)小波變換。正常與房顫心電圖中最大區(qū)別在于P波缺失代之以f波(房顫波),且P波、f波重疊的頻帶范圍為3~12Hz,故在小波相干圖對應的P波尺度范圍內進行特征提取。由于存在個體差異性,且模板的選擇具有一定隨機性,不適用于所有心電信號,因此根據臨床驗證對α進行取值。 表1 參數選取表Tab.1 Parameter selection table 首先,驗證小波相干值的均值、比率、交叉小波相位角方差這3個單一特征的有效性與可行性,圖6展示上述3個特征的箱線圖。由圖6可見,房顫心電的小波相干值均值和比率低于正常心電的對應值,房顫心電的交叉小波相位角方差高于正常心電的對應值,且均無異常值。這說明本文所提取特征可以有效區(qū)分正常心電與房顫心電模式。 其次,采用交叉驗證來驗證所提自動檢測算法的可行性與有效性。在所選取的8個心電記錄中隨機選擇7個作為訓練集,剩余1個作為測試集。以此類推,最終取8組實驗結果的平均值作為最終分類結果。表2列出了每組實驗的準確率、特異性和敏感性以及平均值。從表中可以看出,本文所提算法具有優(yōu)良的分類性能。 圖6 本文所提取特征的箱線圖Fig.6 The boxplot of the wavelet coherence analysis 編號準確率/%敏感性/%特異性/% 0404398.5699.1998.730404898.7397.8594.830401597.2899.1599.96 0493697.6998.1098.76 0642697.1299.3899.89 0645398.8399.241000791097.5498.3298.860837896.7497.0897.88 平均值9 78198.5498.61 表3 本文所提方法與已有方法的檢測性能比較Tab.3 Performance comparison between the proposed method and some existing methods 最后,將本文所提算法與已有算法進行比較(見表3)。從表中可以看出,本文所提算法的平均敏感性及特異性達到SE=98.54%及SP=98.61%,均高于其他算法。 本文提出了一種基于小波相干分析的陣發(fā)性房顫自動檢測方法。首先,對所有的心電信號進行預處理,包括去除噪聲、分段處理、選取模板等;其次,對模板與待測心電信號進行小波相干分析得到小波相干圖;進而,計算小波相干值的均值、比率以及相位角方差,構成房顫心電特征;最后,將所提取特征結合超限學習機完成陣發(fā)性房顫的自動檢測。本文采用MIT-BIH心房顫動數據庫進行數值實驗,實驗結果表明,本文所提算法具有良好的檢測性能,可以為臨床診治提供一定的輔助作用。2 數值實驗結果與分析
2.1 心電數據
2.2 結果與分析
3 結 語