針對陣發性房顫早期出現很短的初次發作較難被檢出的問題,本文提出了一種基于黎曼流形稀疏編碼的檢測算法。本文算法考慮到非線性流形幾何結構更接近真實的特征空間結構,計算協方差矩陣用于表征心率變異性(RR 間期變化),使數據處于黎曼流形空間中。在流形上應用稀疏編碼,將每個協方差矩陣表示為黎曼字典原子的稀疏線性組合,其中稀疏重建損失由仿射不變黎曼度量定義,黎曼字典由迭代的方式學習得到。本文算法與現有算法相比,使用較短心率變異性信號,計算簡單且沒有對參數的依賴,并取得了更優的預測精度。在 MIT-BIH 房顫數據庫上最終分類結果為靈敏度 99.34%、特異度 95.41%、準確率 97.45%,同時在 MIT-BIH 竇性心律數據庫中實現了 95.18% 的特異度。本文提出的高精度陣發性房顫檢測算法在可穿戴設備的長期監測中具有潛在的應用前景。
引用本文: 孟憲輝, 劉明, 熊鵬, 陳健, 楊林, 劉秀玲. 基于黎曼流形稀疏編碼的陣發性房顫檢測算法. 生物醫學工程學雜志, 2020, 37(4): 683-691. doi: 10.7507/1001-5515.201907001 復制
引言
心房顫動(atrial fibrillation,AF)簡稱房顫,為臨床上最常見的持續性心律失常,導致醫療負擔的日益加重[1]。房顫發生時,心臟會失去有效的心房收縮,雖然房顫本身并不危及生命,但會使中風風險增加五倍和心力衰竭發生率增加三倍,從而導致房顫患者的死亡風險超過同齡健康人[2]。按照發作時長房顫可分為陣發性(7 天內自發終止)、持續性(不自愈或持續 7 天以上)和永久性(未終止或終止又復發)三個階段。最初的陣發性房顫(paroxysmal atrial fibrillation,PAF)具有偶發性且持續時間多變,大部分患者在這個階段沒有任何的癥狀表現。隨著病情的發展,發作時間更長、更頻繁,直至演變為持續性房顫,與此同時,治療的效果也在減弱。到最嚴重的階段,即使通過藥物治療或電復律也無法治愈房顫[3]。在此背景下,早期檢測對改善房顫的治療和降低患者的死亡風險具有重要的臨床意義。
然而,早期的陣發性房顫檢測并不是一件容易的事情。臨床上,房顫的診斷主要是依據房顫發作時在心電圖(electrocardiogram,ECG)中的兩個表現,即 P 波消失代之以一系列不規則的纖顫 f 波和 RR 間期(相鄰 R 峰之間的時間間隔)絕對不齊。圖1 顯示了正常心律和心房顫動的 ECG 信號片段。然而,早期間歇性發作時,有些患者在記錄 ECG 數據過程中可能剛好沒有經歷房顫發作,導致 ECG 數據無法記錄到房顫片段。因此,為了及時診斷和排除疾病,需要對心臟活動進行長期監測,如 Holter、智能手環或事件記錄器。而長時程的 ECG 數據量巨大,傳統的人工診斷方法是非常耗時耗力的,并且最初的陣發性房顫持續時間短暫,所以表現在 ECG 信號上的波形變化同樣短暫且微弱,很難在視覺上做出精準的判斷,將會導致較高的漏診率。因此,開展可以提供及時診斷的陣發性房顫自動檢測技術具有重要的現實意義。

目前,已有的房顫自動檢測算法要么依賴心房活動分析,要么依賴心室反應分析,又或者將兩種分析相結合。其中基于心房活動的算法主要包括 P 波缺失[4]、f 波存在[5]和小波熵[6]等。但由于 P 波幅值極微弱,當 ECG 信號中引入噪聲時,這些算法顯示出較差的性能。然而在實際情況下使用連續監測設備時,由于鉛脫離、呼吸或運動等必然因素會產生不可避免的噪音。基于心室反應的算法主要是心率變異性(如 RR 間期)分析[7],與基于心房活動的算法相比,該算法具有更好的抗噪聲能力。因此,利用 RR 間期檢測房顫已成為近年來的首選方法。Tateno 等[8]采用 K-S 檢驗方法來比較房顫信號的 RR 間期和 ΔRR 間期(RR 間期差分值)的密度直方圖,然后采用人工神經網絡來自動地完成房顫檢測。Li 等[9]提出了一種基于圖像的房顫檢測算法,通過將 RR 間期的變化繪制成網格窗格來構建地圖,然后,利用線性支持向量機對空窗格比率、信息熵和概率密度分布進行分類。Kennedy 等[10]研究了兩種基于多變量的分類技術,即隨機森林法和 k 近鄰法,使用 RR 間期樣本熵系數、方差系數、逐次差的均方根,以及絕對中位差等統計特征作為模型的輸入。以往這些算法雖然檢測性能都不錯,但它們通常需要長段數據(40~100 s)來識別長房顫發作,在處理非常短的房顫事件時存在局限性。這是這些算法的一個主要限制,因為在陣發性房顫的第一階段,短暫的發作是非常常見的。這一缺陷也可能會增加房顫的檢出時間,并可能導致對房顫邊界的不準確判斷。
基于此,學者們開始針對早期陣發性房顫持續時間短暫較難檢出的問題展開研究。Petr?nas 等[11]提出了一種基于模糊邏輯的房顫檢測算法,輸入包括 RR 間期不規則、P 波缺失、f 波存在和噪聲水平四個參數,除 RR 不規則外,所有參數均由回聲狀態網絡產生的信號確定。該算法可以檢測到短至 5 拍的房顫發作,但結合了 P 波等多種特征,計算復雜,因此不適用于內存和處理能力有限的動態監測設備中。Shashikumar 等[12]設計了一種基于雙向循環神經網絡并增加了軟注意機制的深度學習模型,該模型將 ECG 信號轉換為頻譜圖來處理,然后采用卷積神經網絡進行時序分析。該算法無需檢測波形 R 峰,沒有依賴任何心房活動或心室反應的特征,而是基于頻譜特征進行深度學習,輸入數據窗口為 30 s 的 ECG 數據記錄。Cui 等[13]提出了一種基于信息的相似性分析和集成方案,通過將 RR 間期時間序列映射為二進制符號序列,設計一個由不同對房顫和正常心律模板組成的集成分類器,比較未知 RR 間期時間序列與模板之間基于信息的差異指數。該算法分別對窗口 50、100、150、200 心搏數做了對比試驗,當截取 150 拍心搏時效果最好。但這兩種算法需要各種參數的調整,才能保證高精度的房顫檢測。
因此,本文的目的是提出一種能夠在長時程 ECG 數據記錄中精準檢測早期較短陣發性房顫的算法。本文算法考慮到特征空間幾何結構的重要性,通過計算 RR 間期時間序列的協方差矩陣,構建非線性流形空間,引入仿射不變黎曼度量來定義稀疏編碼公式,這樣通過黎曼幾何的透鏡可以更好地捕獲理想的數據屬性,以期得到更準確的稀疏表示。
1 方法
陣發性房顫的自動檢測算法主要包含以下幾個技術步驟。首先由 ECG 信號獲取 RR 間期時間序列,采用滑動窗口的方法分割成若干個子序列,并進行注釋,挑戰盡可能小的分割窗口使算法可以檢測非常短的房顫發作;然后使用不同時間序列濾波器對每個子序列進行濾波處理,不同的濾波器提取不同的特征以全面探尋序列傳播變化的本質信息;基于濾波后的特征向量矩陣計算協方差矩陣作為子序列的特征描述子,協方差矩陣自然地將多個特征融合,既捕獲了空間信息又捕獲了統計信息,由于協方差矩陣自身具有對稱半正定屬性,使得特征空間構成黎曼流形空間,不同于歐式空間的是流形空間是彎曲的并且能夠對特征做出更好的解釋;從處理好的所有協方差矩陣中選取一部分數據進行黎曼字典學習,字典中的元素也為協方差矩陣同在黎曼流形空間上,我們使用迭代的方式學習字典,首先采用 k-means 聚類算法進行字典初始化,由初始字典得到一個稀疏表示,固定稀疏表示利用非線性共軛梯度方法逐個更新字典原子,再固定字典求得新的稀疏表示,直至重構的損失函數達到最小值時停止迭代,得到的最優字典原子庫代表全部數據的所有特征;接下來采用譜投影梯度方法將每個協方差矩陣表示為字典中元素的稀疏線性組合,核心思想為使用盡可能少的字典原子來簡潔地表達輸入樣本,以此來揭示其內在結構和本質屬性;最后用得到的稀疏向量與不同類別字典元素分別進行重構,比較不同類別重構值與輸入樣本的殘差來預測類別,從而完成陣發性房顫的自動檢測。
1.1 數據預處理
RR 間期是相鄰搏動之間的時間間隔,在檢測某些心律失常方面顯示出良好的性能。將數據庫中的原始 ECG 信號轉換為 RR 間期時間序列,只存儲關于 R 峰差的信息,而不是使用原始信號中的所有中間樣本,用以降低計算復雜度并突出房顫行為。
將 ECG 信號轉換為 RR 間期序列需要對 ECG 信號進行 R 峰檢測以獲得 R 峰位置。R 峰檢測算法不是本研究的一部分,R 峰位置是從數據庫注釋文件中得到的。在定位 R 峰后,RR 間期計算為:
![]() |
其中, 是第 n 個 R 峰在樣本中的位置,
是以 Hz 為單位的采樣頻率。
在得到數據庫中每條記錄的 RR 間期后,對每條完整的 RR 間期序列采用滑動窗口分割的方法分割成若干個子序列,分割窗口大小用 表示,重疊率為 0。分割后的子序列等于
個 RR 間期,用作算法的輸入。因此算法直接對子序列進行分類,有必要將數據庫中對節拍的注釋轉換為子序列的注釋。轉換方式為:只有當子序列中房顫心搏數占總心搏數比重超過預定義的閾值(在本研究中設定為 0.5)時,才將此子序列分類為房顫。
1.2 黎曼流形空間的構建
1.2.1 協方差矩陣描述子
近年來,區域協方差矩陣(region covariance matrices,RCM)作為黎曼流形的實體廣泛應用于圖像分類、機器學習、圖像檢索、目標識別等領域[14-15],因為黎曼流形空間可以保留特征空間的幾何結構。其次,RCM 的矩陣形式較普通向量形式的特征蘊含更多的信息,多出行和列兩個方向的數據關聯性。
設 代表長度為
的 RR 間期子序列,本文通過結合理論與實驗選擇了 5 種能夠體現 RR 間期序列傳播變化的濾波器提取濾波特征,它們分別為一階差分、二階差分、最大值濾波、最小值濾波和高斯濾波。加上子序列原始值一共得到 6 維特征向量矩陣,標記為
。基于特征向量矩陣中的每個子序列可以構建協方差矩陣
為:
![]() |
其中, 和
分別是第
維特征的均值和第
維特征的平均值。
由公式可以看出協方差矩陣是一個對稱半正定矩陣,其對角線的值表示每個特征的獨立方差,而非對角線上的值表示特征間的相關性。RCM 的構造自然地融合了多個特征,而不必對特征進行規范化或使用混合權重,且同時將空間信息和統計信息嵌入其中。RCM 的大小僅取決于選取的特征向量維數,所以其不再表征子序列的順序與個數,因此算法通常具有很好的魯棒性[16]。
1.2.2 仿射不變黎曼度量
記 是
實對稱正定矩陣構成的空間,賦予黎曼度量的
構成黎曼流形。測地距離定義為連接兩點之間最短曲線的長度,相比歐式距離能更好地反映空間點的拓撲結構。仿射不變黎曼度量(affine invariant Riemannian metric,AIRM)[17]是微分流形上測地距離的一種,該度量構造了流形
和切空間
的微分同胚,
由
上的任意一點
處全體切向量的集合構成,那么流形 M 上從點
出發的等長同向測地線則可由切空間
中的向量
映射得到。在黎曼流形上計算需要的兩個基本運算算子:① 指數映射
;② 對數映射
。前者在流形的切空間上投射對稱點,后者則反向。
黎曼流形的每個切空間有自然的線性結構,因此都有光滑變化的內積。其內積定義為:對任意的 和
,
![]() |
則對于流形上任意兩點X和Y,它們之間的測地距離為:
![]() |
其中, 表示 Frobenius 范數,
為矩陣對數運算。
1.3 黎曼流形上的字典學習與稀疏編碼方法
設 表示一組 N 個對稱正定數據矩陣,其中每個
。設
表示n個對稱正定(symmetric positive definite,SPD)矩陣流形的笛卡兒積的乘積流形,即
。我們的目標是:① 學習三階張量(字典)
,其中每個切片代表 SPD 字典原子
;② 將每個
近似為
中原子的稀疏組合,也就是說
,其中
。因此,在黎曼流形上稀疏編碼與字典學習的總目標函數為[18]:
![]() |
其中 和 Ω 分別是系數向量
和字典張量
的正則化項。
1.3.1 黎曼字典學習
對于所有的數據矩陣,如果存在可用的系數向量 ,那么字典原子更新的子問題可以與式(4)分開并寫成:
![]() |
其中 ,
為正則化參數。
當式(6)取得最小值時得到最優的字典B,從而得到最優B下的最優編碼系數。主要分為兩個步驟:稀疏編碼和字典更新,迭代的進行優化。稀疏編碼過程中求解最優系數時,字典是固定的;更新字典過程中,編碼系數是固定的,然后依次固定 ,采用黎曼共軛梯度方法[19]逐個更新字典中的元素
。采用k均值聚類算法得到字典的初始化值。
要想得到式(6)對于 的黎曼梯度
,首先需要計算它的歐幾里德梯度
,令
,
可得:
![]() |
通過黎曼梯度與歐幾里德梯度的關系式 ,即可導出黎曼梯度。由于黎曼梯度特殊的幾何特性,導致
(這里我們使用
來表示第k次迭代的字典張量)和
分別屬于兩個不同的切空間
和
,不能將它們直接進行運算。因此,需要采用矢量傳輸的方法,這里對于
和
,向量傳輸的定義由以下映射
給出:
。
黎曼共軛梯度方法在第k + 1 次迭代時使用如下遞歸來更新字典:
![]() |
其中 代表下降方向經矢量傳輸可得:(
)
![]() |
其中 。
1.3.2 黎曼稀疏編碼
對于給定 SPD 矩陣 X 和黎曼字典 B,黎曼稀疏編碼的目標是找到一組非負稀疏向量 ,
, 使得
在測地線距離 AIRM 下近似于矩陣
。因此,對于數據矩陣
,我們可以編寫稀疏編碼目標為:
![]() |
其中 是稀疏誘導函數。為簡單起見,我們使用的稀疏懲罰是
,其中
是正則化參數。由于我們規定
,所以直接用可微分的
來代替這個懲罰函數。在文獻[20]中已經證明在滿足集合條件
時,
是一個凸函數。
給定 SPD 矩陣 B、C 和 X,對函數
求導可得:(
)
![]() |
將 改寫為
,可以得到稀疏編碼優化問題的梯度
為:
![]() |
然后利用譜投影梯度(spectral projected gradient,SPG)方法[18]運行以下迭代來更新系數:
![]() |
其中步長 采用 Barzila-Borwein 步長[17],
表示簡化后的投影算子,定義為:
![]() |
1.4 分類的實現
在黎曼字典 B 中的原子被注釋的情況下,算法得到的稀疏向量可直接用于分類。假設字典 B 中有指定類別 i 的 n 個原子,則設 為指定類別的稀疏向量,其中
表示字典原子
的標簽和
表示離散狄拉克
函數。計算指定類別字典與對應指定類別稀疏向量的重構與輸入樣本的殘留誤差,殘差值最小的指定類別的標簽即為輸入樣本的預測標簽。輸入樣本 X 關于類別 i 的殘差計算公式為:
![]() |
2 結果與討論
2.1 評估指標
本文利用準確率(accuracy,Acc)、靈敏度(sensitivity,Se)和特異度(specificity,Sp)三個參數對算法的分類性能進行評估,并且采用混淆矩陣對分類結果進行分析。四個參數具體計算公式分別如式(16)(17)(18)所示:
![]() |
![]() |
![]() |
其中 TP 表示被正確檢測為房顫的信號數量;TN 表示被正確檢測為正常的信號數量;FN 表示房顫信號被錯誤檢測為正常的信號數量;FP 表示正常 ECG 信號被錯誤檢測為房顫發作的信號數量。
2.2 實驗數據
本文以麻省理工學院-貝斯以色列醫院(Massa-chusetts Institute of Technology-Beth Israel Hospital,MIT-BIH)房顫數據庫(MIT-BIH Atrial Fibrillation,MIT-BIH AF)和竇性心律數據庫(MIT-BIH Normal Sinus Rhythm,MIT-BIH NSR)中的數據作為實驗數據來源。MIT-BIH AF 數據庫共包含 23 例房顫患者(多數為陣發性房顫患者)的長期 ECG 數據記錄。單個記錄的持續時間約為 10 h,采樣頻率為 250 Hz,分辨率為 12 bit,采樣帶寬為 0.1~40 Hz。MIT-BIH NSR 數據庫包含 18 條長期 ECG 信號記錄數據,每條記錄持續時間約為 24 h,采樣頻率為 128 Hz。該數據庫中的受試者沒有明顯的心律失常,因此被認為是 NSR,在本文中用于驗證算法的特異度。與每條記錄一起的還有一個參考注釋文件,其中包含每個節拍的位置和類型。以上所有數據庫均可在網址 http:/www.physionet.org/phyobank/database/[21]上免費查閱下載。表1 歸納了數據庫的主要信息。

2.3 實驗結果與討論
對于字典學習,本文從 MIT-BIH AF 數據庫中選用前 6 個人的房顫數據樣本與正常心律數據樣本,分別進行房顫字典學習與正常心律字典學習,即字典原子庫是有注釋的。將學習到的兩類子類字典合并到一起構成一個大字典,求 MIT-BIH AF 數據庫和 MIT-BIH NSR 數據庫所有數據樣本在大字典上的稀疏表示。
2.3.1 實驗結果
本文算法檢測陣發性房顫片段取得了很好的效果。表 2 統計了滑動窗口 、字典原子數
和稀疏編碼正則化項
的變化對 MIT-BIH AF 數據庫檢測結果的影響。由表中可以看出,滑動窗口
的情況下,實驗結果準確率的變化范圍為 77.18%~96.13%,當
較小時結果隨著原子數
的增大而有所增加,當
增加到一定值時結果趨于平穩,而
較大時結果受原子數
的影響變小,最終得到的最佳實驗結果為靈敏度 97.70%、特異度 94.44%、準確率 96.13%。滑動窗口
的情況下,參數的變化對結果的影響沒有大范圍的波動,準確率保持在 93%~97% 的范圍內。在這種情況下得到了本文的最佳房顫檢測結果,其中參數為
、
,此時的靈敏度為 99.34%,特異度為 95.41%,總準確率為 97.45%。滑動窗口
的情況下,大部分結果也趨于平穩,準確率在 90%~95% 之間波動,最佳結果為靈敏度 98.45%、特異度 92.69%、總準確率 95.68%。由以上分析可得,在窗口短至 16 的情況下表現出的性能比窗口為 64 時更好,說明該算法對較短數據更加適用。由此可見,利用本文算法可以對早期陣發性房顫進行精確的檢測。

該模型還在 MIT-BIH NSR 數據庫上進行了驗證。根據表 1,此數據庫中沒有房顫發作,表明 TP=0 和 FN=0。因此,公式(16)等于公式(18),使得 Acc 等于 Sp。此數據庫用來驗證算法的特異度,由表 2 可得在不同滑動窗口下的最佳參數值,表 3 顯示了最佳參數值下的檢測結果。其中滑動窗口 時效果最好,特異度為 96.01%,其次是
的情況下特異度為 95.18%,
時效果最差為 89.96%。對于相同數據量,滑動窗口越小截取到的樣本總數越大,再加上更優的檢測性能導致滑動窗口越小檢測效果越好。

2.3.2 與其他檢測算法的比較
近年來,基于 RR 間期變異性的房顫檢測算法已經發展起來。這里,我們比較了過去 10 年中已發布的基于 RR 間期的算法,使用了相同的數據庫和相同的評估指標,對比結果如表 4[9-10,13,22-25]所示。其中大多數檢測算法使用數據窗口長度在 128 次心搏或 30 s(約 50 次心搏)以上,這種檢測算法會忽略比這更短暫的陣發性房顫發作。而本文采用的數據窗口長度為 33 次心搏,也就是 32 個 RR 間期,增加了短暫陣發性房顫發作的檢出率。這也是本文主要解決的問題,早期的精準檢測可以幫助患者在最佳治療期間及時醫治。

此外,也有與本文使用的數據窗口長度相似的算法被提出來。其中,文獻[25]繪制了 RR 間期與 RR 間期變化的網格地圖,然后通過對非空單元格計數來識別陣發性房顫,使用了與本文相同的 RR 間期數據長度,最終在 MIT-BIH AF 數據庫中達到了 94.30% 的靈敏度和 95.10% 的特異度,但在 MIT-BIH NSR 數據庫中的準確率僅為 84.10%,檢測效果遠不及本文算法。文獻[22]采用將卷積神經網絡和遞歸神經網絡相結合的端到端模型來檢測房顫,使用窗口大小為 31 次心搏,在 MIT-BIH AF 數據庫中靈敏度為 98.98%,特異度為 96.95%,在 MIT-BIH NSR 數據庫中的準確率為 95.01%。雖然文獻[22]獲得了較高的精度,但由于其采用深度學習的算法,需要大量的訓練數據才能訓練出性能較優的檢測模型,從而使得模型性能完全依賴于網絡參數的設置,并不能適用于臨床應用。而本文采用稀疏編碼的算法,將測試樣本轉化為稀疏表示形式,從而簡化了學習任務,參數設置少,因此算法性能僅依賴于實時輸入數據本身的結構。此外將稀疏編碼應用于流形空間使得本文算法能夠更好地提取到房顫患者 ECG 信號的特異性特征,因此對于 MIT-BIH AF 數據庫本文算法獲得的靈敏度最高(99.34%),提示對疾病的檢出較為敏感。
對于 MIT-BIH NSR 數據庫,本文算法在數據長度最小的情況下得到了較高的特異度(96.01%),其中結果略優于本文的算法所需數據長度較長,同理會忽略更短的早期陣發性房顫發作。綜合來看,本文所提房顫檢測算法在兩個數據庫中,不僅所需數據長度較短且取得了優秀的性能表現。因此,該算法能夠滿足陣發性房顫的檢測要求,可以實現早期房顫的精準檢測。
2.3.3 分析與討論
由結果分析我們發現在兩個數據庫上的特異度都有待提高,因此我們對正常心律數據做了進一步研究,發現在檢查原始 ECG 信號時顯示出多個有噪聲的部分,所以導致 R 峰自動檢測算法未能檢測到 R 峰。這將直接轉化為 RR 間期時間序列中的偽影,因為心室節律突然出現不規則,導致部分正常竇性心律被錯誤地歸類為房顫,如圖 2 所示。一種解決方案為在有噪聲的 ECG 信號中,采用性能更好的 R 峰檢測算法,比數據庫官方注釋的檢測算法更先進。另一種解決方案為對每段 ECG 信號中的噪聲水平進行分類,并拒絕超過估計閾值的部分。這樣,含有高噪音水平的片段就會被丟棄,而不是被錯誤地歸類為房顫片段。

a. RR 間期時間序列以及真實和預測的標簽;b. ECG 信號正常心律期間的錯誤預測
Figure2. Error of R peak auto-detectiona. RR interval time series and real and prediction labels; b. error prediction during the NSR period of ECG signal
據我們所知,本研究首次將黎曼流形學習引入陣發性房顫檢測中。流形空間作為歐氏空間中曲線、曲面等概念的推廣,通過黎曼測度來衡量點與點之間的距離,從而得到高維數據的低維流形結構,有助于我們尋找 ECG 信號中 RR 間期變化的本質,找到數據產生的內在規律。通過構建房顫信號字典與正常心律信號字典,字典也存在于黎曼流形空間中,用來稀疏地表示輸入信號具有更強的表達能力,使得算法檢測性能更高。此外,房顫和正常心律模式在不同受試者之間是不完全一樣的,即使同一受試者在不同的活動中,例如在清醒期和睡眠期,ECG 數據也會產生不同類型的波動。因此,如果能夠得到更多不同受試者和同一受試者不同活動下的 ECG 數據,以構建更完備的字典,就可以更加充分地表示房顫和正常心律模式。這種方案將得到更準確的稀疏表示,可以進一步改善本文算法的檢測精度。
3 結論
本文提出了一種基于黎曼流形稀疏編碼的陣發性房顫檢測算法,僅利用 RR 間期數據在分割窗口短至 16 的情況下就表現出良好的性能,是一種適用于早期陣發性房顫的快速、準確、有效的檢測算法。該算法通過構建協方差矩陣將 RR 間期子序列表示為黎曼流形上一點,采用黎曼稀疏編碼算法得到最優的稀疏表示,稀疏表示向量直接用來分類。相比于傳統的分類算法,降低了復雜度,且不需要很多的參數設置,可以實現長期的動態監測。我們未來的工作將專注于基于可穿戴設備陣發性房顫患者的長期記錄提出的技術進行前瞻性評估。
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
引言
心房顫動(atrial fibrillation,AF)簡稱房顫,為臨床上最常見的持續性心律失常,導致醫療負擔的日益加重[1]。房顫發生時,心臟會失去有效的心房收縮,雖然房顫本身并不危及生命,但會使中風風險增加五倍和心力衰竭發生率增加三倍,從而導致房顫患者的死亡風險超過同齡健康人[2]。按照發作時長房顫可分為陣發性(7 天內自發終止)、持續性(不自愈或持續 7 天以上)和永久性(未終止或終止又復發)三個階段。最初的陣發性房顫(paroxysmal atrial fibrillation,PAF)具有偶發性且持續時間多變,大部分患者在這個階段沒有任何的癥狀表現。隨著病情的發展,發作時間更長、更頻繁,直至演變為持續性房顫,與此同時,治療的效果也在減弱。到最嚴重的階段,即使通過藥物治療或電復律也無法治愈房顫[3]。在此背景下,早期檢測對改善房顫的治療和降低患者的死亡風險具有重要的臨床意義。
然而,早期的陣發性房顫檢測并不是一件容易的事情。臨床上,房顫的診斷主要是依據房顫發作時在心電圖(electrocardiogram,ECG)中的兩個表現,即 P 波消失代之以一系列不規則的纖顫 f 波和 RR 間期(相鄰 R 峰之間的時間間隔)絕對不齊。圖1 顯示了正常心律和心房顫動的 ECG 信號片段。然而,早期間歇性發作時,有些患者在記錄 ECG 數據過程中可能剛好沒有經歷房顫發作,導致 ECG 數據無法記錄到房顫片段。因此,為了及時診斷和排除疾病,需要對心臟活動進行長期監測,如 Holter、智能手環或事件記錄器。而長時程的 ECG 數據量巨大,傳統的人工診斷方法是非常耗時耗力的,并且最初的陣發性房顫持續時間短暫,所以表現在 ECG 信號上的波形變化同樣短暫且微弱,很難在視覺上做出精準的判斷,將會導致較高的漏診率。因此,開展可以提供及時診斷的陣發性房顫自動檢測技術具有重要的現實意義。

目前,已有的房顫自動檢測算法要么依賴心房活動分析,要么依賴心室反應分析,又或者將兩種分析相結合。其中基于心房活動的算法主要包括 P 波缺失[4]、f 波存在[5]和小波熵[6]等。但由于 P 波幅值極微弱,當 ECG 信號中引入噪聲時,這些算法顯示出較差的性能。然而在實際情況下使用連續監測設備時,由于鉛脫離、呼吸或運動等必然因素會產生不可避免的噪音。基于心室反應的算法主要是心率變異性(如 RR 間期)分析[7],與基于心房活動的算法相比,該算法具有更好的抗噪聲能力。因此,利用 RR 間期檢測房顫已成為近年來的首選方法。Tateno 等[8]采用 K-S 檢驗方法來比較房顫信號的 RR 間期和 ΔRR 間期(RR 間期差分值)的密度直方圖,然后采用人工神經網絡來自動地完成房顫檢測。Li 等[9]提出了一種基于圖像的房顫檢測算法,通過將 RR 間期的變化繪制成網格窗格來構建地圖,然后,利用線性支持向量機對空窗格比率、信息熵和概率密度分布進行分類。Kennedy 等[10]研究了兩種基于多變量的分類技術,即隨機森林法和 k 近鄰法,使用 RR 間期樣本熵系數、方差系數、逐次差的均方根,以及絕對中位差等統計特征作為模型的輸入。以往這些算法雖然檢測性能都不錯,但它們通常需要長段數據(40~100 s)來識別長房顫發作,在處理非常短的房顫事件時存在局限性。這是這些算法的一個主要限制,因為在陣發性房顫的第一階段,短暫的發作是非常常見的。這一缺陷也可能會增加房顫的檢出時間,并可能導致對房顫邊界的不準確判斷。
基于此,學者們開始針對早期陣發性房顫持續時間短暫較難檢出的問題展開研究。Petr?nas 等[11]提出了一種基于模糊邏輯的房顫檢測算法,輸入包括 RR 間期不規則、P 波缺失、f 波存在和噪聲水平四個參數,除 RR 不規則外,所有參數均由回聲狀態網絡產生的信號確定。該算法可以檢測到短至 5 拍的房顫發作,但結合了 P 波等多種特征,計算復雜,因此不適用于內存和處理能力有限的動態監測設備中。Shashikumar 等[12]設計了一種基于雙向循環神經網絡并增加了軟注意機制的深度學習模型,該模型將 ECG 信號轉換為頻譜圖來處理,然后采用卷積神經網絡進行時序分析。該算法無需檢測波形 R 峰,沒有依賴任何心房活動或心室反應的特征,而是基于頻譜特征進行深度學習,輸入數據窗口為 30 s 的 ECG 數據記錄。Cui 等[13]提出了一種基于信息的相似性分析和集成方案,通過將 RR 間期時間序列映射為二進制符號序列,設計一個由不同對房顫和正常心律模板組成的集成分類器,比較未知 RR 間期時間序列與模板之間基于信息的差異指數。該算法分別對窗口 50、100、150、200 心搏數做了對比試驗,當截取 150 拍心搏時效果最好。但這兩種算法需要各種參數的調整,才能保證高精度的房顫檢測。
因此,本文的目的是提出一種能夠在長時程 ECG 數據記錄中精準檢測早期較短陣發性房顫的算法。本文算法考慮到特征空間幾何結構的重要性,通過計算 RR 間期時間序列的協方差矩陣,構建非線性流形空間,引入仿射不變黎曼度量來定義稀疏編碼公式,這樣通過黎曼幾何的透鏡可以更好地捕獲理想的數據屬性,以期得到更準確的稀疏表示。
1 方法
陣發性房顫的自動檢測算法主要包含以下幾個技術步驟。首先由 ECG 信號獲取 RR 間期時間序列,采用滑動窗口的方法分割成若干個子序列,并進行注釋,挑戰盡可能小的分割窗口使算法可以檢測非常短的房顫發作;然后使用不同時間序列濾波器對每個子序列進行濾波處理,不同的濾波器提取不同的特征以全面探尋序列傳播變化的本質信息;基于濾波后的特征向量矩陣計算協方差矩陣作為子序列的特征描述子,協方差矩陣自然地將多個特征融合,既捕獲了空間信息又捕獲了統計信息,由于協方差矩陣自身具有對稱半正定屬性,使得特征空間構成黎曼流形空間,不同于歐式空間的是流形空間是彎曲的并且能夠對特征做出更好的解釋;從處理好的所有協方差矩陣中選取一部分數據進行黎曼字典學習,字典中的元素也為協方差矩陣同在黎曼流形空間上,我們使用迭代的方式學習字典,首先采用 k-means 聚類算法進行字典初始化,由初始字典得到一個稀疏表示,固定稀疏表示利用非線性共軛梯度方法逐個更新字典原子,再固定字典求得新的稀疏表示,直至重構的損失函數達到最小值時停止迭代,得到的最優字典原子庫代表全部數據的所有特征;接下來采用譜投影梯度方法將每個協方差矩陣表示為字典中元素的稀疏線性組合,核心思想為使用盡可能少的字典原子來簡潔地表達輸入樣本,以此來揭示其內在結構和本質屬性;最后用得到的稀疏向量與不同類別字典元素分別進行重構,比較不同類別重構值與輸入樣本的殘差來預測類別,從而完成陣發性房顫的自動檢測。
1.1 數據預處理
RR 間期是相鄰搏動之間的時間間隔,在檢測某些心律失常方面顯示出良好的性能。將數據庫中的原始 ECG 信號轉換為 RR 間期時間序列,只存儲關于 R 峰差的信息,而不是使用原始信號中的所有中間樣本,用以降低計算復雜度并突出房顫行為。
將 ECG 信號轉換為 RR 間期序列需要對 ECG 信號進行 R 峰檢測以獲得 R 峰位置。R 峰檢測算法不是本研究的一部分,R 峰位置是從數據庫注釋文件中得到的。在定位 R 峰后,RR 間期計算為:
![]() |
其中, 是第 n 個 R 峰在樣本中的位置,
是以 Hz 為單位的采樣頻率。
在得到數據庫中每條記錄的 RR 間期后,對每條完整的 RR 間期序列采用滑動窗口分割的方法分割成若干個子序列,分割窗口大小用 表示,重疊率為 0。分割后的子序列等于
個 RR 間期,用作算法的輸入。因此算法直接對子序列進行分類,有必要將數據庫中對節拍的注釋轉換為子序列的注釋。轉換方式為:只有當子序列中房顫心搏數占總心搏數比重超過預定義的閾值(在本研究中設定為 0.5)時,才將此子序列分類為房顫。
1.2 黎曼流形空間的構建
1.2.1 協方差矩陣描述子
近年來,區域協方差矩陣(region covariance matrices,RCM)作為黎曼流形的實體廣泛應用于圖像分類、機器學習、圖像檢索、目標識別等領域[14-15],因為黎曼流形空間可以保留特征空間的幾何結構。其次,RCM 的矩陣形式較普通向量形式的特征蘊含更多的信息,多出行和列兩個方向的數據關聯性。
設 代表長度為
的 RR 間期子序列,本文通過結合理論與實驗選擇了 5 種能夠體現 RR 間期序列傳播變化的濾波器提取濾波特征,它們分別為一階差分、二階差分、最大值濾波、最小值濾波和高斯濾波。加上子序列原始值一共得到 6 維特征向量矩陣,標記為
。基于特征向量矩陣中的每個子序列可以構建協方差矩陣
為:
![]() |
其中, 和
分別是第
維特征的均值和第
維特征的平均值。
由公式可以看出協方差矩陣是一個對稱半正定矩陣,其對角線的值表示每個特征的獨立方差,而非對角線上的值表示特征間的相關性。RCM 的構造自然地融合了多個特征,而不必對特征進行規范化或使用混合權重,且同時將空間信息和統計信息嵌入其中。RCM 的大小僅取決于選取的特征向量維數,所以其不再表征子序列的順序與個數,因此算法通常具有很好的魯棒性[16]。
1.2.2 仿射不變黎曼度量
記 是
實對稱正定矩陣構成的空間,賦予黎曼度量的
構成黎曼流形。測地距離定義為連接兩點之間最短曲線的長度,相比歐式距離能更好地反映空間點的拓撲結構。仿射不變黎曼度量(affine invariant Riemannian metric,AIRM)[17]是微分流形上測地距離的一種,該度量構造了流形
和切空間
的微分同胚,
由
上的任意一點
處全體切向量的集合構成,那么流形 M 上從點
出發的等長同向測地線則可由切空間
中的向量
映射得到。在黎曼流形上計算需要的兩個基本運算算子:① 指數映射
;② 對數映射
。前者在流形的切空間上投射對稱點,后者則反向。
黎曼流形的每個切空間有自然的線性結構,因此都有光滑變化的內積。其內積定義為:對任意的 和
,
![]() |
則對于流形上任意兩點X和Y,它們之間的測地距離為:
![]() |
其中, 表示 Frobenius 范數,
為矩陣對數運算。
1.3 黎曼流形上的字典學習與稀疏編碼方法
設 表示一組 N 個對稱正定數據矩陣,其中每個
。設
表示n個對稱正定(symmetric positive definite,SPD)矩陣流形的笛卡兒積的乘積流形,即
。我們的目標是:① 學習三階張量(字典)
,其中每個切片代表 SPD 字典原子
;② 將每個
近似為
中原子的稀疏組合,也就是說
,其中
。因此,在黎曼流形上稀疏編碼與字典學習的總目標函數為[18]:
![]() |
其中 和 Ω 分別是系數向量
和字典張量
的正則化項。
1.3.1 黎曼字典學習
對于所有的數據矩陣,如果存在可用的系數向量 ,那么字典原子更新的子問題可以與式(4)分開并寫成:
![]() |
其中 ,
為正則化參數。
當式(6)取得最小值時得到最優的字典B,從而得到最優B下的最優編碼系數。主要分為兩個步驟:稀疏編碼和字典更新,迭代的進行優化。稀疏編碼過程中求解最優系數時,字典是固定的;更新字典過程中,編碼系數是固定的,然后依次固定 ,采用黎曼共軛梯度方法[19]逐個更新字典中的元素
。采用k均值聚類算法得到字典的初始化值。
要想得到式(6)對于 的黎曼梯度
,首先需要計算它的歐幾里德梯度
,令
,
可得:
![]() |
通過黎曼梯度與歐幾里德梯度的關系式 ,即可導出黎曼梯度。由于黎曼梯度特殊的幾何特性,導致
(這里我們使用
來表示第k次迭代的字典張量)和
分別屬于兩個不同的切空間
和
,不能將它們直接進行運算。因此,需要采用矢量傳輸的方法,這里對于
和
,向量傳輸的定義由以下映射
給出:
。
黎曼共軛梯度方法在第k + 1 次迭代時使用如下遞歸來更新字典:
![]() |
其中 代表下降方向經矢量傳輸可得:(
)
![]() |
其中 。
1.3.2 黎曼稀疏編碼
對于給定 SPD 矩陣 X 和黎曼字典 B,黎曼稀疏編碼的目標是找到一組非負稀疏向量 ,
, 使得
在測地線距離 AIRM 下近似于矩陣
。因此,對于數據矩陣
,我們可以編寫稀疏編碼目標為:
![]() |
其中 是稀疏誘導函數。為簡單起見,我們使用的稀疏懲罰是
,其中
是正則化參數。由于我們規定
,所以直接用可微分的
來代替這個懲罰函數。在文獻[20]中已經證明在滿足集合條件
時,
是一個凸函數。
給定 SPD 矩陣 B、C 和 X,對函數
求導可得:(
)
![]() |
將 改寫為
,可以得到稀疏編碼優化問題的梯度
為:
![]() |
然后利用譜投影梯度(spectral projected gradient,SPG)方法[18]運行以下迭代來更新系數:
![]() |
其中步長 采用 Barzila-Borwein 步長[17],
表示簡化后的投影算子,定義為:
![]() |
1.4 分類的實現
在黎曼字典 B 中的原子被注釋的情況下,算法得到的稀疏向量可直接用于分類。假設字典 B 中有指定類別 i 的 n 個原子,則設 為指定類別的稀疏向量,其中
表示字典原子
的標簽和
表示離散狄拉克
函數。計算指定類別字典與對應指定類別稀疏向量的重構與輸入樣本的殘留誤差,殘差值最小的指定類別的標簽即為輸入樣本的預測標簽。輸入樣本 X 關于類別 i 的殘差計算公式為:
![]() |
2 結果與討論
2.1 評估指標
本文利用準確率(accuracy,Acc)、靈敏度(sensitivity,Se)和特異度(specificity,Sp)三個參數對算法的分類性能進行評估,并且采用混淆矩陣對分類結果進行分析。四個參數具體計算公式分別如式(16)(17)(18)所示:
![]() |
![]() |
![]() |
其中 TP 表示被正確檢測為房顫的信號數量;TN 表示被正確檢測為正常的信號數量;FN 表示房顫信號被錯誤檢測為正常的信號數量;FP 表示正常 ECG 信號被錯誤檢測為房顫發作的信號數量。
2.2 實驗數據
本文以麻省理工學院-貝斯以色列醫院(Massa-chusetts Institute of Technology-Beth Israel Hospital,MIT-BIH)房顫數據庫(MIT-BIH Atrial Fibrillation,MIT-BIH AF)和竇性心律數據庫(MIT-BIH Normal Sinus Rhythm,MIT-BIH NSR)中的數據作為實驗數據來源。MIT-BIH AF 數據庫共包含 23 例房顫患者(多數為陣發性房顫患者)的長期 ECG 數據記錄。單個記錄的持續時間約為 10 h,采樣頻率為 250 Hz,分辨率為 12 bit,采樣帶寬為 0.1~40 Hz。MIT-BIH NSR 數據庫包含 18 條長期 ECG 信號記錄數據,每條記錄持續時間約為 24 h,采樣頻率為 128 Hz。該數據庫中的受試者沒有明顯的心律失常,因此被認為是 NSR,在本文中用于驗證算法的特異度。與每條記錄一起的還有一個參考注釋文件,其中包含每個節拍的位置和類型。以上所有數據庫均可在網址 http:/www.physionet.org/phyobank/database/[21]上免費查閱下載。表1 歸納了數據庫的主要信息。

2.3 實驗結果與討論
對于字典學習,本文從 MIT-BIH AF 數據庫中選用前 6 個人的房顫數據樣本與正常心律數據樣本,分別進行房顫字典學習與正常心律字典學習,即字典原子庫是有注釋的。將學習到的兩類子類字典合并到一起構成一個大字典,求 MIT-BIH AF 數據庫和 MIT-BIH NSR 數據庫所有數據樣本在大字典上的稀疏表示。
2.3.1 實驗結果
本文算法檢測陣發性房顫片段取得了很好的效果。表 2 統計了滑動窗口 、字典原子數
和稀疏編碼正則化項
的變化對 MIT-BIH AF 數據庫檢測結果的影響。由表中可以看出,滑動窗口
的情況下,實驗結果準確率的變化范圍為 77.18%~96.13%,當
較小時結果隨著原子數
的增大而有所增加,當
增加到一定值時結果趨于平穩,而
較大時結果受原子數
的影響變小,最終得到的最佳實驗結果為靈敏度 97.70%、特異度 94.44%、準確率 96.13%。滑動窗口
的情況下,參數的變化對結果的影響沒有大范圍的波動,準確率保持在 93%~97% 的范圍內。在這種情況下得到了本文的最佳房顫檢測結果,其中參數為
、
,此時的靈敏度為 99.34%,特異度為 95.41%,總準確率為 97.45%。滑動窗口
的情況下,大部分結果也趨于平穩,準確率在 90%~95% 之間波動,最佳結果為靈敏度 98.45%、特異度 92.69%、總準確率 95.68%。由以上分析可得,在窗口短至 16 的情況下表現出的性能比窗口為 64 時更好,說明該算法對較短數據更加適用。由此可見,利用本文算法可以對早期陣發性房顫進行精確的檢測。

該模型還在 MIT-BIH NSR 數據庫上進行了驗證。根據表 1,此數據庫中沒有房顫發作,表明 TP=0 和 FN=0。因此,公式(16)等于公式(18),使得 Acc 等于 Sp。此數據庫用來驗證算法的特異度,由表 2 可得在不同滑動窗口下的最佳參數值,表 3 顯示了最佳參數值下的檢測結果。其中滑動窗口 時效果最好,特異度為 96.01%,其次是
的情況下特異度為 95.18%,
時效果最差為 89.96%。對于相同數據量,滑動窗口越小截取到的樣本總數越大,再加上更優的檢測性能導致滑動窗口越小檢測效果越好。

2.3.2 與其他檢測算法的比較
近年來,基于 RR 間期變異性的房顫檢測算法已經發展起來。這里,我們比較了過去 10 年中已發布的基于 RR 間期的算法,使用了相同的數據庫和相同的評估指標,對比結果如表 4[9-10,13,22-25]所示。其中大多數檢測算法使用數據窗口長度在 128 次心搏或 30 s(約 50 次心搏)以上,這種檢測算法會忽略比這更短暫的陣發性房顫發作。而本文采用的數據窗口長度為 33 次心搏,也就是 32 個 RR 間期,增加了短暫陣發性房顫發作的檢出率。這也是本文主要解決的問題,早期的精準檢測可以幫助患者在最佳治療期間及時醫治。

此外,也有與本文使用的數據窗口長度相似的算法被提出來。其中,文獻[25]繪制了 RR 間期與 RR 間期變化的網格地圖,然后通過對非空單元格計數來識別陣發性房顫,使用了與本文相同的 RR 間期數據長度,最終在 MIT-BIH AF 數據庫中達到了 94.30% 的靈敏度和 95.10% 的特異度,但在 MIT-BIH NSR 數據庫中的準確率僅為 84.10%,檢測效果遠不及本文算法。文獻[22]采用將卷積神經網絡和遞歸神經網絡相結合的端到端模型來檢測房顫,使用窗口大小為 31 次心搏,在 MIT-BIH AF 數據庫中靈敏度為 98.98%,特異度為 96.95%,在 MIT-BIH NSR 數據庫中的準確率為 95.01%。雖然文獻[22]獲得了較高的精度,但由于其采用深度學習的算法,需要大量的訓練數據才能訓練出性能較優的檢測模型,從而使得模型性能完全依賴于網絡參數的設置,并不能適用于臨床應用。而本文采用稀疏編碼的算法,將測試樣本轉化為稀疏表示形式,從而簡化了學習任務,參數設置少,因此算法性能僅依賴于實時輸入數據本身的結構。此外將稀疏編碼應用于流形空間使得本文算法能夠更好地提取到房顫患者 ECG 信號的特異性特征,因此對于 MIT-BIH AF 數據庫本文算法獲得的靈敏度最高(99.34%),提示對疾病的檢出較為敏感。
對于 MIT-BIH NSR 數據庫,本文算法在數據長度最小的情況下得到了較高的特異度(96.01%),其中結果略優于本文的算法所需數據長度較長,同理會忽略更短的早期陣發性房顫發作。綜合來看,本文所提房顫檢測算法在兩個數據庫中,不僅所需數據長度較短且取得了優秀的性能表現。因此,該算法能夠滿足陣發性房顫的檢測要求,可以實現早期房顫的精準檢測。
2.3.3 分析與討論
由結果分析我們發現在兩個數據庫上的特異度都有待提高,因此我們對正常心律數據做了進一步研究,發現在檢查原始 ECG 信號時顯示出多個有噪聲的部分,所以導致 R 峰自動檢測算法未能檢測到 R 峰。這將直接轉化為 RR 間期時間序列中的偽影,因為心室節律突然出現不規則,導致部分正常竇性心律被錯誤地歸類為房顫,如圖 2 所示。一種解決方案為在有噪聲的 ECG 信號中,采用性能更好的 R 峰檢測算法,比數據庫官方注釋的檢測算法更先進。另一種解決方案為對每段 ECG 信號中的噪聲水平進行分類,并拒絕超過估計閾值的部分。這樣,含有高噪音水平的片段就會被丟棄,而不是被錯誤地歸類為房顫片段。

a. RR 間期時間序列以及真實和預測的標簽;b. ECG 信號正常心律期間的錯誤預測
Figure2. Error of R peak auto-detectiona. RR interval time series and real and prediction labels; b. error prediction during the NSR period of ECG signal
據我們所知,本研究首次將黎曼流形學習引入陣發性房顫檢測中。流形空間作為歐氏空間中曲線、曲面等概念的推廣,通過黎曼測度來衡量點與點之間的距離,從而得到高維數據的低維流形結構,有助于我們尋找 ECG 信號中 RR 間期變化的本質,找到數據產生的內在規律。通過構建房顫信號字典與正常心律信號字典,字典也存在于黎曼流形空間中,用來稀疏地表示輸入信號具有更強的表達能力,使得算法檢測性能更高。此外,房顫和正常心律模式在不同受試者之間是不完全一樣的,即使同一受試者在不同的活動中,例如在清醒期和睡眠期,ECG 數據也會產生不同類型的波動。因此,如果能夠得到更多不同受試者和同一受試者不同活動下的 ECG 數據,以構建更完備的字典,就可以更加充分地表示房顫和正常心律模式。這種方案將得到更準確的稀疏表示,可以進一步改善本文算法的檢測精度。
3 結論
本文提出了一種基于黎曼流形稀疏編碼的陣發性房顫檢測算法,僅利用 RR 間期數據在分割窗口短至 16 的情況下就表現出良好的性能,是一種適用于早期陣發性房顫的快速、準確、有效的檢測算法。該算法通過構建協方差矩陣將 RR 間期子序列表示為黎曼流形上一點,采用黎曼稀疏編碼算法得到最優的稀疏表示,稀疏表示向量直接用來分類。相比于傳統的分類算法,降低了復雜度,且不需要很多的參數設置,可以實現長期的動態監測。我們未來的工作將專注于基于可穿戴設備陣發性房顫患者的長期記錄提出的技術進行前瞻性評估。
利益沖突聲明:本文全體作者均聲明不存在利益沖突。