基于隱馬爾科夫模型的滑動窗口投票策略的QRS波群形態識別

宋鑫海,韓京宇,郎 杭,毛 毅

(南京郵電大學計算機學院,江蘇 南京 210023)

1 引言

據世界衛生組織統計,心腦血管疾病已成為威脅人類健康的頭號殺手,預計到2030年,全球心臟病致死人數將達2 330萬[1,2]。心電圖ECG(ElectroCardioGraphy)作為一種無創檢測手段,在臨床心臟病預警和診斷中得到了廣泛應用[3]。人工識別心電異常需要專業醫生,存在效率低和成本高的問題,因此自動心電異常診斷受到廣泛關注[4]。

在自動心電異常診斷中,主要根據心拍從前到后的不同的波形特征來識別可能的癥狀。一個完整的心拍從前到后依次分為P波、QRS波和T波,如圖1所示。

Figure 1 Waveform of a complete heart beat圖1 一個完整的心拍波形

圖1中QRS波群波形高大,是判斷心拍的主要標志,能反映左右心室的除極狀況,是診斷各種心臟疾病的重要依據。如圖2所示,右束支阻滯患者通常在V1或V2導聯的QRS波群表現出rSr′形態;下壁心肌梗死患者在aVR導聯的QRS波群多呈rs或rS形態[5,6]。因此,QRS波群形態的準確識別是提高心電異常診斷精確度的關鍵。自動QRS波群形態識別面臨的主要挑戰有：(1)波群形態種類多達數十種,不同種類之間形態差異細小,容易誤判;(2)許多波群形態的區分依賴前后相繼的波形識別,為實現對不同波形的對齊,需要對波形整體進行縮放,增加了自動識別難度。

Figure 2 Different QRS waveforms of patients圖2 患者不同的QRS波形示例

針對QRS形態識別,早期采用差分法、斜率判別法和小波變換[7,8]等對信號進行時頻分析,以確定波峰和波谷等關鍵位置,完成波群形態分析。這些方法能較好地識別正常QRS形態,但對于異常或復雜波形,容易誤判。后期,一些工作從ECG波形曲線擬合的角度解決問題,如采用高階多項式曲線逼近ECG信號[9,10]、組合多種函數模擬波形[11]。這些方法簡單直觀,但前者參數較多,難以控制波形擬合效果,后者只能識別少數波形。

針對上述問題,本文提出一種基于狀態轉移受限的隱馬爾科夫模型滑動窗口投票策略SWVHMM(Sliding Window Voting based on Hidden Markov Model)識別QRS波群形態。該方法對前后相繼的波段窗口進行拼接,構建狀態轉移受限的隱馬爾科夫模型來識別每個拼接的形態,進而確定最大可能的波群類型。該方法主要包含波段劃分和窗口樣本提取、構建狀態轉移受限的隱馬爾科夫模型以及投票預測3個部分。首先,根據峰值和極值點提取4類前后相繼的波段,并對每個波段設置滑動窗口,抽取所有窗口樣本。其次,將窗口樣本作為觀測(Observation),波形作為狀態(State),構建隱馬爾科夫模型。最后,為待預測樣本構建多個拼接,對所有拼接進行投票,得票最多的即為最終波群類型。

本文方法的主要優點在于：(1)每個波段取多個滑動窗口,能夠捕獲各波段最有效的局部觀測。(2)隱馬爾科夫模型可以充分體現波段前后相繼的依賴關系,根據波群的整體特征來識別波形,不會囿于局部波段來判斷整體。(3)通過對波段的序列進行拼接,產生多個樣本代表,進而通過投票進行判斷,保證了方法的泛化性和魯棒性。

2 相關工作

2.1 去噪方法

常用心電信號去噪方法有數字濾波器、基于小波變換的閾值去噪、數學形態學法(Mathematical Morphology)和經驗模態分解EMD(Empirical Mode Decomposition)。

數字濾波器針對不同噪聲頻率,選擇合適的高通或低通濾波器去噪,實時處理,但濾波后T波容易被誤判為R波[12];小波變換對信號進行時頻分析,通過伸縮平移對信號多尺度細化,實踐證明對非平穩信號效果較好[8];數學形態學法以積分幾何和代數等數學理論為基礎,使用被稱為“探針”的結構元素來最大程度地保留信號形態特征,較好地保證了ECG完整性,但是去除低頻噪聲效果不佳[13];EMD將信號分解為多個本征模態函數IMF(Intrinsic Mode Function),再通過組合IMF達到去噪目的[14]。此外,還有許多組合改進算法,例如文獻[15]將EMD與小波變換相結合,取得了較好的去噪效果。

2.2 QRS波形識別

傳統檢測QRS波群主要有閾值檢測和模板匹配等方法。閾值檢測因為占用內存小、運行速度快,被廣泛用于可穿戴式心電設備,但該方法對信號質量要求高,容易受噪聲干擾,準確率有待提高[16]。模板匹配直觀簡單,但檢測速度慢,且嚴重依賴模板[17]。此外,還有多種新方法用于QRS檢測,例如文獻[18]結合小波變換和希爾伯特變換檢測含噪聲的心電信號;文獻[19]提出了基于一維卷積神經網絡1D CNN(one Dimensional Convolutional Neural Network)的QRS檢測算法。這些方法對于正常波形能取得較好檢測效果,但對于異常或不規則波形的檢測效果不佳。

針對傳統檢測方法的不足,面對QRS復雜形態,還有研究從擬合曲線角度進行檢測。毛玲等[10]提出了一種QRS波群關鍵點結合有限狀態機的形態分析算法FSM-QRS(Finite State Machine for QRS complex),采用二次多項式曲線分段逼近心電信號;然后分析各段曲線的陡峭性和單調程度等特征,以此檢測并分析波峰和波谷等關鍵點;最后,構建有限自動機,輸入關鍵點序列,實現QRS形態識別。上述算法能夠識別多種波形,但需要設定多種參數,如預定角度和閾值等,效果難以控制。Caldas等[11]提出基于數學函數組合建模殘差RCMF(Residual modeling with Composition of Mathematical Functions)的特征提取方法,利用Gaussian、Rayleigh和Mexican Hat函數構建數學模型擬合波群形態,對所選函數進行標準化使其峰值相等,便于合并函數,然后將所得模型與輸入QRS波群形態進行對比,計算二者的均方誤差,將均方誤差最小值所對應的模型作為波形分類的最終模型。上述方法僅能夠識別少數波形,且時間復雜度較高。

3 方法概述

3.1 方法概述和相關概念

本文方法主要包含3個階段,即波段識別和窗口樣本提取、構建狀態轉移受限的隱馬爾科夫模型和多窗口樣本投票,本文方法框架如圖3 所示。第1階段,在清除切跡后,根據極值點和R波識別4個波段,各波段利用滑動窗口提取多個窗口樣本;第2階段,對各波段窗口樣本聚類以構建狀態轉移受限的隱馬爾科夫模型;第3階段,針對待預測波群,利用滑動窗口獲得多個候選樣本,每個候選樣本通過模型得到預測波形,最后投票所有波形,輸出最終結果。

Figure 3 Framework of the proposed method圖3 本文方法框架

一個QRS波群是一個心拍中的連續點序列〈(X1,Y1),…,(Xi,Yi),…,(Xn,Yn)〉,其中Yi(1≤i≤n)表示Xi時刻的電壓。一個QRS波群中的波段定義如下：

定義1(波段) 一個QRS波群從前到后分為Q、R、S和R′ 4段。各波段呈現不同波形,記為{Q,q,Nq},{R,r,Nr},{S,s,Ns},{R′,r′,Nr′}。其中Nq、Nr、Ns和Nr′表示對應波段不存在相應波形。

在QRS波群中,每個波段的波形根據形態、方向、幅度及前后相繼關系來判定,如圖4所示。第1個向下波屬于Q波段,振幅大于0.5 mV的記為Q波,小于0.5 mV的記為q波,表現為平直線段的記為Nq;第1個向上波屬于R波段,第2個向下波屬于S波段,在S波段后出現的第1個正向波屬于R′波段,各波段波形描述和Q波段相同。當心臟傳導系統異常,QRS波群還會出現切跡和頓挫現象。其中,切跡指QRS波群中各波的上升或下降部分有方向和斜率發生改變的節段;頓挫指QRS波群中各波的上升或下降部分有斜率發生改變,而方向不變的節段[10]。根據歐洲定量心電圖通用標準CSE(European project entitled Common Standards in quantitative Electrocardiography)[20]規定,若切跡等各種波形存在,必須滿足時間不少于6 ms,振幅不小于0.1 mV 2個條件。

Figure 4 Typical QRS shape圖4 典型QRS形態

3.2 波形去噪

心電信號受基線漂移、工頻干擾和肌電干擾等噪聲影響[3],波形識別前要先濾波去噪,使信號平滑。使用小波變換去除工頻干擾和肌電干擾[8]：采用bior4.4小波基,根據噪聲頻率分布對小波8階分解后的低頻和高頻分量進行軟閾值去噪,最后小波重構得到連續信號。之后對該連續信號使用6階Butterworth 濾波器,其截止頻率為0.5 Hz,去除基線[21]。

之后采用經典的Pan-Tompkins雙閾值算法[22]進行R波定位。通常QRS波群時長不超過0.2 s,故以R波為中心,前后各取0.1 s作為形態檢測區域。整個去噪前后對比如圖5所示。

Figure 5 Waveforms comparison before and after denoising圖5 去噪前后波形對比

3.3 各波段定位和窗口樣本提取

QRS波群形態復雜多變,各波段檢測會受到切跡的干擾。如圖6所示,在R波下降部分出現方向向上的曲線ABC,會讓極小值點A識別為S波段,導致波形誤判。因此,在利用極值點定位波段前要先排除切跡干擾。

Figure 6 qR wave with notch圖6 伴隨切跡的qR波形

(1)

(2)

其中,sample指采樣率。

搜索波群內所有極值點,在排除切跡后,檢測剩余極值點位置來判斷波段類型。如圖6的Q波段,在峰值R左側存在首個極小值點D,取D前后一段時間作為Q波段[6]。關于波段長度,文獻[6]指出正常Q波段存在0.04 s,不同波段存在時間不同。

定義2(權重窗口樣本) 每個波段對應一個區域〈(X1,Y1),…,(Xn,Yn)〉,給定滑動窗口長度m,可以獲取n-m+1個窗口樣本,第i個窗口對應的子序列是〈(Xi,Yi),…,(Xi+m,Yi+m)〉。

本文根據每個樣本相對波段中心的位置來確定其權重,因為越靠近波段中心的樣本,形狀特征越強,越能夠有效識別波段形態。記波段區域長度Lc=Xn-X1,中心位置Cc=X1+Lc/2;第i個窗口中心位置Dc=Xi+(Xi+m-Xi)/2,其權重如式(3)所示：

(3)

樣本距離波段中心越近,權重越大,當距離為0時,權重為1。所有樣本權重在[0,1]變換,便于計算。

4 構建狀態轉移受限的隱馬爾科夫模型

4.1 波形聚類構建觀測代表

由于每個波段的滑動窗口會取得多個窗口樣本,所有實例對應的窗口樣本數目會非常龐大,不便于直接采用隱馬爾科夫建模。為此,提出通過聚類識別窗口樣本的類簇,用每個類簇表示一類窗口樣本,從而用有限個類簇中心來表征可能的觀測。

常見的聚類方法有：基于距離的劃分聚類方式,代表算法為K-means(K-means clustering algorithm),要求隨機選取k個對象作為聚類中心,之后計算每個對象和聚類中心的距離[23];基于密度的聚類方式,如DBSCAN(Density Based Spatial Clustering of Applications with Noise),其任意選擇一個核心點,然后尋找所有從該點密度可達的樣本,形成類簇[23];基于圖論的聚類方式,如譜聚類,利用無向帶權圖表示數據集中對象,之后根據相關準則劃分該無向圖,將聚類轉換為圖的構建與切分問題[24]。其中,K-means適用于發現特定形狀的數據集,而DBSCAN和譜聚類可以發現任意類型的數據。相比于譜聚類,DBSCAN還具有較強的抗噪性,更符合心電信號微弱、易受噪聲影響的特點。同時,由于時間序列的滯后性,如歐氏距離、曼哈頓距離等不足以描述序列相似性,為更好地度量窗口樣本的相似性,使用規范化后的互相關系數作為距離度量[25]。

給定2個長度為m的窗口樣本S和T,為搜索二者在形態上的相似性,固定T,令S左右滑動產生相位偏移,如式(4)所示計算每次偏移后二者內積Rt：

(4)

整個滑動過程共產生2m-1個內積,如式(5)所示：

C(S,T)=max(Rt)=max(Rw-m),

w=1,2,…,2m-1

(5)

其中,t<0表示S向左滑動|t|個單位長度,最大值即為最佳相位匹配,然后如式(6)進行歸一化。互相關系數Nor越大,說明樣本相似性越強。

(6)

傳統DBSCAN聚類要求設置密度半徑等參數,若參數設置不當,則會影響聚類效果[26]。本文提出基于樣本稀疏程度來進行DBSCAN聚類,使用樣本平均密度和標準差描述樣本稠密性,具體步驟如下所示：

(1)離群點判斷。每個樣本,包括隨機初始樣本在拓展形成小類簇前,都要判斷是否為離群點,如果是離群點就忽略該點,其判斷依據為樣本的稠密程度,定義如下所示：

定義3(窗口樣本密度) 給定窗口樣本xi,其密度Dxi形式化如式(7)所示：

(7)

其中,N是以xi為圓心、Rxi為密度半徑的圓內樣本個數。

定義4(平均密度) 給定窗口樣本xi,平均密度是其密度半徑內所有樣本的平均密度,如式(8)所示：

(8)

定義5(密度標準差) 給定窗口樣本xi,密度標準差用來衡量其密度半徑內樣本的離散程度,如式(9)所示：

(9)

(3)將類簇C內每個樣本作為新的簇心向外拓展,重復步驟(2)形成多個類簇,直到不能發現新的小類簇。之后合并這些小類簇組成最終大類簇Z。對于剩余樣本,重復步驟(2)和步驟(3)不斷形成新的大類簇,直至所有樣本聚類完成。

算法1 波形聚類輸入：波段窗口樣本集合Ls;簇內樣本個數k。輸出：觀測集合Observation。1：Z←?;2：while L.size>0：3： /?隨機選擇樣本P,根據式(7)～式(9)判斷是否為離群點?/4： if P is not noise：5： search k samples which are close to P;6： calculate k samples′ mean as R;7： //得到以P為圓心、R為半徑的小類簇C8： Z.add(C);9： for each Ci in C：10： repeat step 6 and step 7 to get clusters Bi;11： Z.add(Bi);12： end for13： calculate the cluster center Tc for Z;14： Observation.add(Tc);15： Ls.remove(Z);16： else：17： Ls.remove(P);18： end if 19：end while 20：return Observation;

算法的時間復雜度為O(n*m2),其中n為波段窗口樣本個數,m為窗口樣本的長度。

4.2 構建狀態轉移受限的隱馬爾科夫

本文方法所構建的隱馬爾科夫模型狀態的轉移是受限的,只能從當前波段狀態轉移到后繼波段狀態,如從Q波段轉移至R波段,從R波段轉移至S波段等。一個轉移受限的隱馬爾科夫模型可表示為五元組HMM=(H,S,O,TR,EM)。其中,H表示初始狀態(Q,q,Nq);S代表所有狀態,包括Q,q,Nq,R等;O代表觀測,即聚類獲得的類簇中心;TR是狀態轉移概率;EM代表發射概率,即在每個狀態下能夠看到觀測的概率。假設存在波群qRs和qrs各200個,qRS和qrS各100個,其狀態轉移概率如圖7所示。

Figure 7 Illustration of state transition probability圖7 狀態轉移概率示意圖

采用監督學習方法構建隱馬爾科夫模型：在所有訓練樣本上統計各個觀測值,就可以構建對應的隱馬爾科夫模型。其中,初始狀態分布H可由訓練集樣本中Q波段出現頻率估計;狀態轉移概率TRu,v和發射概率EMv,o可形式化為式(10)和式(11)所示：

(10)

(11)

其中,Au,v表示當前時刻狀態u在下一時刻轉移到狀態v的出現頻數;Bv,o表示狀態v時觀測為o的出現頻數。

5 基于多窗口波群形態預測

采用候選樣本投票方法來決定最終的波群形態。首先,構建樣本代表,然后對每個樣本代表計算預測波群類型,最后多個波群投票,具體步驟如下所示：

(1)構建候選樣本。

定義7(候選樣本) 給定一個QRS波群,設其Q波段對應X個觀測樣本,分別記為{Q1,…,Qi,…,QX};R波段對應Y個觀測樣本,記為{R1,…,Rj,…,RY};S波段對應ZS個觀測樣本,記為{S1,…,Sm,…,SZS};R′波段對應W個觀測樣本,記為{R′1,…,R′n,…,R′W}。則一個候選樣本是4段的拼接,記為O=(Qi,Rj,Sm,R′n),即一個待預測QRS波群會產生X*Y*ZS*W個候選樣本。

(2)模型預測。

對每個候選樣本,應用隱馬爾科夫模型獲取其對應的波群形態,即在已知模型HMM和觀測序列O情況下,求解最大可能波群形態T,可形式化為式(12)所示：

(12)

P(O|HMM)為常數,問題轉化為求解maxP(T,O|HMM)。利用Viterbi算法求解該最優化問題,定義μtb(u)為波段tb(tb∈{Q,R,S,R′})時,終止狀態為u的所有路徑中最大概率,則波段tb+1最大概率表示如式(13)所示：

μtb+1(u)=max(EMu,Otb+1·μtb(v)·TRv,u)

(13)

求解到R′波段預測狀態后,向前回溯求解波群類型T。

(3)投票。

算法2 QRS波群形態預測輸入：待預測波群Wy;隱馬爾科夫模型HMM。輸出：最大概率波群類型。1：dict←?; // key為預測波形,value為概率2：generate all candidate samples S for Wy;3：for each s in S：// 各候選樣本s4： calculate weight probabilities P for s;5： predict the waveform T according to the Viterbi algo-rithm;6： if T in dict：7： dict[T]←dict[T]+P;8： else：9： dict[T]←P;10： end if11：end for12：sort dict in descending order by value;13：return dict.keys()[0];

時間復雜度為O(q*m),其中O(m)為單次運行維特比算法的時間復雜度,q為候選樣本個數。

6 實驗與結果分析

6.1 實驗設置

為檢驗本文方法性能,采用12導聯心電數據DS-COM[27]進行測試。該數據集共包含6 000個樣本,采樣率為500 Hz,時長為10 s。為測試結果正確與否,請心電專家對QRS波群進行人工識別。實驗共抽取11種形態的QRS波群,分別為：RS波350個,qRs波300個,qRsr′和rSr′波各200個,Qr、QS、R和rS波各150個,qR,Rs和qrs波各100個,總計1 950個。

本文實驗采取6輪交叉驗證,即將數據集均分為6份,依次分6輪進行模型構建和測試,最后取平均實驗結果。每輪中,5份數據用于構建馬爾科夫模型,另外1份用做測試集。為驗證本文方法的有效性,采用精確度PR(PRecision)、召回率RL(RecaLl)和FS(F1-Score)度量指標,其定義分別如式(14)～式(16)所示：

(14)

(15)

(16)

其中,TP表示正樣本被模型預測為正類的個數,FP表示負樣本被預測為正類的個數,FN表示正樣本被預測為負類的個數。

6.2 方法性能評估

(1)簇內樣本個數k對波段聚類效果的影響。

使用輪廓系數表征k對各波段聚類效果。圖8展示了k對Q波段聚類的影響,隨著k增大,輪廓系數先上升后下降,并且幅度較大,說明聚類對k較敏感。

Figure 8 Effect of k on the clustering silhouette coefficient圖8 k對聚類輪廓系數的影響

(2)滑動窗口尺寸效果分析。

本文方法通過滑動窗口提取各波段樣本,其窗口尺度會影響最終預測效果。圖9展示了窗口尺度對11種波形平均預測結果的影響。當窗口尺度在[4,6]時,波形識別性能沒有較大波動;但當窗口過大或過小時波形預測效果不佳。

Figure 9 Effect of window size on waveform prediction圖9 窗口尺度對波形預測的影響

(3)聚類算法的比較。

對于K-means,利用手肘法找出最佳k值并聚類[23];對于譜聚類,采用KNN(K-Nearest Neighbor)算法遍歷所有樣本構建無向圖,并使用NCut(Normalized Cut)方式切圖聚類[24]。使用輪廓系數來表征3種算法對Q波段聚類的效果。如表1所示,K-means效果最差,其輪廓系數為0.389 5;譜聚類次之,輪廓系數為0.501 4;本文算法最佳,輪廓系數為0.520 7。這是因為相比于K-means,譜聚類和DBSCAN對數據集不敏感,同時DBSCAN抗噪聲能力更強。

Table 1 Results comparison of clustering algorithms

6.3 對比實驗

與FSM-QRS[10]、RCMF[11]和1DLBPAT(1-Dimension Local Binary Pattern with Adaptive Threshold)[28]3種方法在DS-COM數據集上進行對比測試,相關參數與實驗結果如下所示：

(1)FSM-QRS：連續擬合點的個數k=5,誤差閾值ε=0.01,預定閾值thr=0.004,方向角度α=0.03,β=0.02。

(2)RCMF：Rayleigh和Mexican Hat函數的參數λ∈{0.1,0.5,…,10};Gaussian函數的參數μ=0,σ2∈{0.1,0.5,…,10}。

(3)SWVHMM：簇內樣本個數k為4,窗口尺度為5。

(4)1DLBPAT：對經過EMD分解后的IMF使用LBPAT算法檢測R波,其閾值AT由截止頻率為3 Hz的低通濾波器濾波后的信號和該信號的平均值線性疊加構成。

如圖10～圖14所示,對比FSM-QRS,SWVHMM在RS波形上,PR、RL和FS分別提升了11.07%,6.37%和8.54%;在qRs波形上,PR、RL和FS各自提升了9.50%,2.60%和5.72%;Qr波形的PR、RL和FS分別提升了4.89%,4.67%和4.78%;8種波形PR、RL和FS平均提升了6.46%,3.42%和4.85%。

Figure 10 Experimental results of RS waveform comparison圖10 RS波形對比實驗結果

相比RCMF,SWVHMM在RS波形上,PR、RL和FS分別提升了7.13%,6.28%和6.68%;在qRs波形上,PR、RL和FS各自提升7.00%,6.81%和6.94%;Qr波形的RL下降了1.96%,PR和FS上升了7.52%和2.86%。對比1DLBPAT,SWVHMM在R波形上,PR、RL和FS分別提升了1.33%,3.19%和2.27%。

結果表明,相比于FSM-QRS、RCMF和1DLBPAT,SWVHMM使得指標FS的平均值提高了5.97%,5.49%和2.27%。該方法能夠取得更好效果的原因在于：相比于RCMF和1DLBPAT,SWVHMM利用滑動窗口提取波段局部特征,再利用隱馬爾科夫模型描述波段前后相繼的整體關系,從而能夠識別更多形態;相比FSM-QRS控制參數較少,降低了參數敏感性。

7 結束語

QRS波群形態是許多心臟疾病診斷的主要依據。早期的差分法和斜率判別法等只能處理簡單形態;后期從波形曲線入手,利用各種函數擬合心電信號,但存在只能識別少數波形或參數過多、難以控制的缺點。本文提出基于隱馬爾科夫的滑動窗口投票策略,包含波段劃分及窗口樣本的提取、構建狀態轉移受限的隱馬爾科夫模型和多窗口預測3個部分,可識別出多種QRS波群形態,并且準確率更高、魯棒性更好。后續將考慮如何更好地描述窗口樣本相似性和序列聚類,以及如何更好地定位R波。