心律失常的自動檢測對于早期預防與診斷心血管疾病具有重要意義。傳統的心律失常診斷受限于專家知識，缺乏多維的特征表示能力，算法復雜，不適用于穿戴式心電監測設備。本研究提出基于自回歸移動平均（ARMA）模型擬合的特征提取方法，以不同類型心拍作為模型輸入，利用 ARMA 模型對心律失常信號以合適的階數進行系數擬合，完成心電特征提取，并將特征向量分別輸入支持向量機（SVM）與 K 近鄰分類器（KNN）進行心電自動分類。所提算法采用 MIT-BIH 心律失常數據庫與房顫數據庫為數據集進行驗證，結果表明：ARMA 模型擬合系數組成的特征工程結合支持向量機分類器得到查全率為 98.2%，查準率為 98.4%，F1 指數為 98.3%。該算法具有較高的性能，滿足臨床診斷需求，算法復雜度低，可采用低功耗嵌入式處理器進行實時運算，適用于穿戴式心電監測設備的實時預警。
引用本文： 閆惠君, 莫思特, 黃華, 劉彥. 基于自回歸移動平均模型擬合的高效心律失常自動分類研究. 生物醫學工程學雜志, 2021, 38(5): 848-857, 868. doi: 10.7507/1001-5515.202101054 復制
引言
心律失常是心血管疾病（cardiovascular disease，CVD）的重要表現，由心律失常所引發的心血管疾病每年導致全球約 1 790 萬人死亡，占總死亡人數的 31%[1]，其早期的檢測與治療非常關鍵。穿戴式心電監測設備的飛速發展也正是為了實時監測個體心律失常信號特征，這一應用需求依賴高效便捷的心電檢測算法。心電圖（electrocardiogram，ECG）記錄了心臟至體表的電生理活動，通過其形態和節奏能有效評估心血管系統的狀態，是臨床檢測心律失常的必要手段之一。
心律失常的自動檢測（automatic detection of arrhythmia，ADA）是指使用計算機設備代替繁瑣耗時的人工 ECG 數據分析進行心律失常的診斷篩查[2]。隨著近年來機器學習的飛速發展，ADA 受到了更多學者的關注與研究。ECG 信號微弱不穩定的特性使得其特征提取成為研究的關鍵[3]。Raj 等[4]提出基于離散余弦變換的正交 stockwell 變換方法，利用支持向量機（support vector machine，SVM）分類器對來自 MIT-BIH 心律失常數據庫 16 類心律失常信號進行分類，但缺乏合適的 R 峰檢測算法，需要對特征集的大小進行最佳選擇以減少分類器的訓練時間。Mazaheri 等[5]提出在心電預處理后對 ECG 信號進行分段，提取其形態特征、頻域特征和非線性指標，利用元啟發式優化算法刪除多余或不相關的特征實現降維，但研究所應用的 NSGA II 特征選擇方法處理時間較長，難以實現實時穿戴式心電監測。
由于心跳的能量集中于 ECG 波形的 QRS 復合波中，R 波特有的陡峭特性有利于 ECG 的特征提取。Merah 等[6]設計了一種基于平穩小波變換的方法，省略了變換算法中每個級別的系數抽取，可更好地執行離群值檢測。Jung 等[7]提出利用稀疏判別提取重要的小波系數來自動檢測室性早搏的方法，通過定位適當的基準點以確定 ECG 信號中準確的 QRS 波群。基于經驗模式分解（empirical mode decomposition，EMD）的方法可以有效去除 ECG 信號中的基線漂移[8]，是心電預處理中的重要步驟。EMD 和離散小波變換（discrete wavelet transform，DWT）的組合方法為 ECG 降噪提供了更好的時間分辨率[9]。Lu 等[10]提出了一種基于 EMD 和 SVM 分類器集成的新方法，使用具有可變范圍編碼和動態搜索策略的改進遺傳算法優化統計特征子集，以獲得更好的分類精度，但缺少尋找分類精度的評價過程。
可見，針對 ECG 信號的特征提取方面使用時頻域及形態學特征的方法較為普遍。目前研究 ECG 信號自動分類的主流，一是對心電信號作較為復雜的變換提取相關特征參數，二是利用深度學習的端到端算法直接對原始信號進行分類[11-12]。但是，這些方法準確率的保證對原始 ECG 信號質量要求較高且計算過程復雜，消耗功耗大，且易受到 P 波、T 波等雜波的干擾，高效低功耗的算法研究少，因此現有研究成果難以應用于便攜式心電設備。
而直觀的自回歸移動平均模型（auto regressive moving average，ARMA）作為表征依賴數據最常用的模型之一，在語音處理、生物醫學、電力消耗預測和計量經濟學等領域有著悠久的傳統，但在心電分析方面的研究較少[13-15]。例如，Aguila 等[16]提出假最近鄰算法（false nearest neighbor，FNN）與 ARMA 模型相結合的時間序列建模對心電信號進行預測。Heurtefeux 等[17]采用 ARMA 模型，通過計算赤池信息準則對心電嵌入式系統進行了評價，建立了 p 階和 q 階以達到模型復雜度和準確性之間的權衡，相對較低的階數被認為在受約束的嵌入式系統中實現是可行的。
為了使用 ARMA 模型進行心電預測，需要選擇合適的模型參數。針對上述問題，本研究在 QRS 復合波檢測上利用雙微分法與自適應閾值準確判斷 R 峰位置，以此作為心拍分割基準點。在特征提取階段，應用 ARMA 模型將相關時序轉化為獨立時序，進行統計處理分析。上述步驟基于 ECG 信號特性，探求一種魯棒性與準確性更高的 ECG 分類算法，在極低的算法復雜度下可以對 ECG 實現高效分類，算法簡單，適用于穿戴式心電檢測，有利于對心血管疾病患者進行日常監護和預警，促進實時在線心電監護的發展。
1 實驗數據與方法
1.1 實驗數據
本文的數據來自 MIT-BIH 的心律失常數據庫（MIT-BIH Arrhythmia Database，MIT-BIH-AR），該數據庫包含從 1975 至 1979 年在心律不齊實驗室中的案例，選擇 48 條信號的雙導聯數據，每條 ECG 信號記錄時間長度均超過 30 min，采樣頻率：360 Hz，采樣分辨率：11 bit，采樣帶寬：0.1～40 Hz。考慮到房顫作為最常見的心律失常疾病，具有重要的研究意義，選擇 MIT-BIH 房顫數據庫（MIT-BIH Atrial Fibrillation Database，MIT-BIH-AF）中的房顫信號作為補充數據。將正常 ECG 信號與不同種類的心律失常信號交叉組成訓練集與測試集，以驗證房顫檢測算法的特異度。MIT-BIH-AR 數據庫與 MIT-BIH-AF 數據庫中正常與非正常心搏的具體數量、記錄數及持續時間如表1 所示。
表1
MIT-BIH 心電數據庫統計
Table1.
Statistics sorted out from MIT-BIH database
1.2 心電自動分類算法流程
本研究提出的算法流程如圖1 所示，首先對 ECG 信號進行預處理，利用小波閾值方法濾除噪聲、基線等干擾，采用雙微分法檢測 QRS 波群并提取 R 峰。根據 R 峰位置，在波峰前后各取一定數量的樣本點分割心電心拍，根據心拍注釋分類不同的心律失常信號，完成輸入數據的準備。利用 ARMA 模型與單變量特征選擇擬合對輸入的心拍數據進行系數擬合并降維，與心拍標注類型組成特征向量輸入到 SVM 和 K 近鄰（K-nearest neighbor，KNN）兩種分類器中，用于預測輸入樣本標簽。最后，評估模型識別的查準率、查全率等參數，為心律失常的自動檢測方法提供參考。
圖1
心電信號自動檢測流程圖
Figure1.
Flow diagram of automatic detection of ECG signal
1.3 心電預處理與 R 波提取
原始心電信號多數存在基線漂移（≤ 5 Hz）、肌電噪聲（10～300 Hz）和工頻噪聲（50 Hz）等干擾信號，心電預處理需要濾除噪聲干擾以提高后續檢測的準確性。本研究利用小波閾值去噪方法，對信號進行正交小波變換，得到一組小波分解系數
，變換后的小波系數中含有大量有用 ECG 信號的數據，ECG 噪聲的圖像分布情況與小波分解系數分布相同，因此小波閾值去噪方法可有效濾除噪聲。由于小波系數大小不一，將絕對值小的設為 0，對絕對值大的采用軟閾值方法進行保存或壓縮，處理后的小波系數要更加平滑，得到其估計值
，利用小波系數重構新的 ECG 信號，達到去噪效果。針對 ECG 信號的低頻特點，閾值λ的選擇根據 sqtwolog 準則確定[18]，能消除更多噪聲，其計算表達式與軟閾值函數表達式分別如式（1）和式（2）所示，其中n代表小波系數向量長度。
 |
 |
對濾波后的波形進行雙微分處理，以獲得心電波形的陡度與高度特性：通過在樣本點的的左右兩側設置閾值區間，分別尋找最大平均斜率（
，
）與最小平均斜率（
，
），交叉求差后求取最大值
，如式（3）所示。設置微分閾值
，其值的確定根據前 8 個峰值的平均值計算。雙微分處理后利用低通濾波消除雙峰現象，滑動窗口積分使得 ECG 絕對振幅增大，寬度根據經驗設為 17 個采樣點，此時預處理后的信號更便于檢測 R 峰。
 |
檢測 R 波時選擇自適應閾值的設計，其選擇準則為：① 能在各種類型的記錄下良好地區分目標與非目標，達到較高準確率與可用度；② 選擇經驗參數需固定，針對不同心電記錄能跟隨波形實時穩健地變化；③ 自適應閾值可防錯檢、漏檢。例如利用“不應期”的含義，在檢測到一個 QRS 波后短時間內不再檢出 QRS 波，防止錯檢；同時可以進行波形檢測回溯，以防漏檢。基于以上三點，設計自適應機制如式（4）~（5）所示（THR0為低閾值，THR1為高閾值）。
 |
 |
其中，peak指檢測到的峰值，peak_buffer指存儲了當前峰值之前 8 個峰值，mean（）指均值。THR1_lim與THR0_lim為兩個經驗常數，代表兩個閾值變化的下限，在這里分別取值為 0.3 和 0.23。R 波提取的流程如圖2 所示，對 ECG 進行預處理檢測結果分別如圖3、圖4 所示。
圖2
R 波提取流程
Figure2.
Process of R wave extraction
圖3
ECG 預處理結果
Figure3.
Results of ECG preprocessing
圖4
“雙微分”處理后標記 R 波波峰提取示意圖
Figure4.
Schematic diagram of extraction of marked R-wave peaks after "double differential" processing
1.4 心拍分割
心律失常信號的異常主要體現在心電圖所蘊含波形形態中，完整的心電記錄包含多個心拍，直接全局處理不利于精細化識別，常通過截取單個心拍來分析判別，分割長度對于最終信號的分類準確度至關重要。根據先驗知識可知，波形異常主要表征在 PR 間期以及 QRS 波群的異常，因此在預處理的基礎上，通常以陡峭的 R 波波峰為關鍵點，前后分別截取 100 和 150 個樣本點共 250 個樣本點作為樣本心拍，心拍注釋采用 MIT-BIH 數據庫中的專家標注，根據標注中的心律失常類型進行后續的特征提取與分類工作。
本文共選取正常心拍（N）、左束支傳導阻滯心拍（L）、右束支傳導阻滯心拍（R）、室性早搏心拍（V）以及房顫信號（AF）五類心電信號。從 MIT-BIH-AR 數據庫的 48 條記錄中的 MLII 通道數據中進行心拍劃分，總計獲取 96 871 個樣本心拍，其中正常心拍（N）有 74 466 個，左束支傳導阻滯心拍（L）有 8 051 個，右束支傳導阻滯心拍（R）有 7 239 個，室性早搏心拍（V）有 7 115 個。從 MIT-BIH-AF 數據庫的 24 條記錄中的 MLII 通道數據中獲取 79 654 個心拍作為房顫（AF）類別。根據 MIT-BIH 數據庫的心拍注釋，對每類心電信號選取 2 500 個心拍，共 12 500 個心拍，利用 1.3 節 R 波提取后的位置進行心拍的分割。圖5 為所截取各類型心拍示意圖。
圖5
不同 ECG 心拍提取示意圖
Figure5.
Schematic diagram of different ECG heartbeat extraction
1.5 特征工程
ARMA 模型是目前應用最廣泛的線性平穩時間序列預測模型，其建模技術具有容易實現、處理速度快和能同時對多類 ECG 進行分類的優點，此外，它還能對觀測數據進行自主平滑和濾波，進一步達到去噪效果，從而估計出觀測數據的真值[19]。ARMA（p，q）的模型結構如式（6）所示。
 |
蘊含在心電數據的取值大小可由 ARMA（p，q）模型的模型參數
刻畫，心電數據的頻率特征可由 
和 
刻畫。基于模型適用于一維信號的平穩特性，可以通過對心電數據建立 ARMA（p，q）模型，提取擬合的模型參數作為每個樣本心拍的特征值。ARMA 模型可細分為 AR 模型和 MA 模型，通過對一維時序信號的觀測值反映出其基本特征。算法具體步驟如下：
步驟一：基于心電信號的 ARMA 線性輸出方程式可寫為：
 |
方程（7）可分解為 AR 過程和 MA 過程：
AR 過程：
 |
MA 過程：
 |
其中，G(x)為 ECG 信號的時間序列；ai和bj分別為 AR 系數和 MA 系數；p和q分別為 AR 階次和 MA 階次；R(x)為 AR 誤差的序列；n(x)為 ARMA 預測誤差，即 ARMA 建模誤差。
步驟二：在建模中，階次p和q的選擇至關重要，本研究選擇最小描述長度（minimum description length，MDL）準則與敏感性函數來確定模型階次[20]。MDL 準則如式（10）所示：
 |
式中，
為建模誤差的方差，N為數據窗口的寬度，Np為模型的階次，在這里Np=p（AR 過程）或Np=q（MA 過程）。另外，MDL 的敏感性函數可以反映 ARMA 模型階次對擬合效果的影響，其表示見式（11）：
 |
其中，
是 MDL 的變化量，
反映了模型階次Np對 MDL 的影響。
步驟三：通過 MDL 準則確定 ARMA 模型階數范圍在 2~4 之間，利用殘差檢驗驗證所選擇階數，如果殘差是隨機正態分布且不自相關，則說明所有有用信號都已被提取到 ARMA 模型中。同時，利用德賓-沃森統計量（Durbin-Watson Statistic，DWS）檢驗模型誤差項。
步驟四：根據 MDL 準則確定 ARMA（4, 2）模型系數作為心電信號系數。利用 R 語言中的 arima 函數擬合 ARMA 模型系數，并對得到系數矩陣進行降維。系數矩陣由擬合得到的 ARMA 系數xi和對應心拍標簽組成，使用主成分分析（principal component analysis，PCA）降維方法可有效降低由原始心拍數據維度過大造成的 ARMA 系數復雜問題。
對于線性 PCA 計算過程如下：首先，計算擬合系數均值
并根據式（12）計算協方差矩陣 C。
 |
式（13）給出了特征向量V和特征值E的計算。其中，特征值由高到低排序組成特征向量，各特征向量中，特征值最高的特征向量為數據的主成分。
 |
線性 PCA 技術的重要優點是其簡單的幾何解釋和彌補后續分類算法過擬合的計算特性作為 ECG 的特征矩陣[21]。另外，為了對比計算，使用 kernel PCA 用非線性映射完成 PCA 的線性運算[22]。核主成分分析計算核矩陣的主特征向量，而非線性方法中計算協方差矩陣，核矩陣系數 
的計算如式（14）所示。
 |
其中，n是數據樣本的數量，k為核函數，例如高斯函數或多項式函數等。
整體流程如圖6 所示。
圖6
ARMA 模型擬合系數流程圖
Figure6.
Flow diagram of ARMA model fitting coefficient
殘差計算結果如圖7 所示，圖中分別顯示了殘差計算結果圖、QQ 分位圖以及自相關和偏自相關圖。可以發現殘差檢驗分布滿足正態分布，且信號圖形趨近于白噪聲波形，通過自相關和偏自相關圖可以發現模型階數選擇符合實際提取心拍。DWS 值越接近 2，則越說明殘差不存在一階相關性，本研究計算結果為 1.999 998，因此殘差結果符合要求。
圖7
ARMA 模型階數殘差檢驗結果
Figure7.
The results of ARMA model order residual test
利用選定階數預測心拍序列，可得圖8 所示波形，此時模型擬合的仿真信號與原始信號最為接近。故本研究選用 ARMA（4，2）模型，只要確定各參數的取值，對應的模型也就完全確定，將其用來擬合數據、考察數據內在的統計特征以及作最佳預測時，在數學上的分析處理是便捷高效的。模型序列的譜密度為有理譜密度，而連續譜密度可以用有理譜密度來逼近，并能達到理想的近似程度。
圖8
ARMA (4, 2)模型預測心拍序列結果
Figure8.
Prediction results of heart beat sequence by ARMA (4, 2) model
1.6 機器學習方法分類與比較
本研究使用了 SVM 與 KNN 兩種分類器。SVM 是一種監督學習的二分類模型，其主要目的是找到一種超平面對樣本數據進行分割，從而轉化為求解凸二次規劃的問題。然而在大部分的分類任務中，單純的線性關系完成分類實例的準確分割是較為困難的，SVM 模型利用非線性核函數將分類任務從線性不可分的低維空間轉換到線性可分的高維空間，然后通過算法分析找出最佳超平面，最后返回到原始空間以獲得初始空間的解[23]。
KNN 通常基于某種距離度量找出訓練集中與其最靠近的 k 個訓練樣本，利用“投票法”選擇 k 個樣本中出現最多的類別標記作為預測結果。它是“懶惰學習”的代表，訓練時間開銷為零，雖然其結構與運算簡單，但它的泛化錯誤率不超過貝葉斯最優分類器錯誤率的兩倍，因此在許多分類任務中的應用也是非常廣泛的[23]。
本研究選擇使用 K 次交叉驗證方法（K=10）無重復地抽樣劃分輸入數據作為訓練集和驗證集，訓練集用于訓練模型，驗證集評估模型性能，以此來降低選擇模型過擬合的現象。另外，通過使用網格化搜索函數來優化分類模型，提高模型的穩定性以及泛化性能。
2 實驗評估與分析
2.1 評估指標
本文利用查全率（以符號
表示）、查準率（以符號
表示）、評估精度分數
來評估心電信號的分類效果，由真陽性（true positive，TP）、假陰性（false negative，FN）、假陽性（false positive，FP）的數量以及數據實例決定，指標計算公式如式（15）~（17）所示：
 |
 |
 |
在多分類問題中，對于第i類心律失常樣本，TP 表示該類被正確分類的樣本數量，FN 表示第i類的樣本被錯誤地分類為第j類的數量，FP 表示第j類的樣本被錯誤地分類為第i類的數量。為了便于表達，以符號
表示五種類型 ECG 的平均查全率，以符號
表示五種類型 ECG 的平均查準率，以符號
表示五種類型 ECG 的平均
指數。
2.2 實驗結果分析
部分正常心電信號樣本心拍的 ARMA 模型擬合系數降維后如表2 所示，表中的
、
系數分別表征了輸入心拍的取值大小以及頻率特征，且結構簡單。將其與心電類型標簽組合作為分類器模型的輸入，表3 為各類別心電信號的分類結果統計，KNN 分類器對 MIT-BIH AR 數據集分類得到查全率、查準率以及
指數分別為 98.1%、97.7% 以及 97.9%，而 SVM 分類器的結果分別為 98.2%、98.4% 以及 98.3%。在 M 類分類問題中，混淆矩陣可以顯示分類器獲得分類結果的詳細分布，而且是一個 M 維方陣。圖9 為 KNN 與 SVM 分類算法結果混淆矩陣。從結果可以看出，文中提出的預處理以及特征工程算法在 MIT-BIH 數據庫上有較好的分類性能。從分類器特性上看，KNN 對噪聲有較好的魯棒性，但模型有欠擬合或過擬合傾向，SVM 在建立了一個成功的模型基礎上，對新樣本的分類較為準確，目前應用更加廣泛。對比 1.5 小節中兩種 PCA 降維技術，它們應用在 KNN 分類器上的效果差異并不大，這是由于降維操作實際上是對數據進行了一定程度的近似，丟失了部分信息所導致的。從綜合衡量查全率與查準率的
指數看出，SVM 分類器要優于 KNN 分類器。
表2
部分正常心電信號樣本心拍的 ARMA 模型擬合系數
Table2.
ARMA model fitting coefficient of cardiac beats in partial normal ECG signal samples
表3
機器學習方法分類效果統計
Table3.
Classification accuracy statistics of machine learning methods
圖9
KNN 分類算法與 SVM 分類算法結果混淆矩陣
Figure9.
Confusion matrix of SVM and KNN classification algorithm results
在處理 ECG 數據方面，本研究避免了特異性處理特殊心律失常的復雜步驟，采用簡單的過濾和標記技術實現算法的低復雜特性。另外 ARMA 模型應用 R 語言的 arima 函數來實現，其數據收集與回歸過程通過測試模型擬合時間驗證得到計算復雜度為O(n)，可以近似看作線性模型。從結果來看，提出方法的查準率與查全率均可以達到很好的分類結果，且特有的低復雜度適用于便攜式的心電實時監測。
表4 顯示了近年來關于心律失常自動檢測工作[10,24-28]與本文提出算法在計算復雜度與評估指標方面的對比，所得分類效果均基于 MIT-BIH 數據庫。相較于計算復雜度相同的算法，本文所提出算法在各個評估指標上均具有顯著優勢；相較于 Oh 等[24]、Yildirim 等[25]提出的基于長短期記憶人工神經網絡（long short-term memory，LSTM）的復雜模型，本文提出算法的復雜度更低。另外，考慮算法的實際應用方面，選擇高精度低功耗 ECG 傳感芯片 KS1081 作為參考芯片，計算不同算法在相同運算情況下的功耗。KS1081 芯片的單通道功耗約為 120 μW，平均執行速度為 532 MIPS（每秒處理的百萬級語言指令數）32 位定點運算，可滿足 16 位采樣精度數據運算需求，且可實現在長時間連續信號監測情況下提升設備續航時間的功能。故應用該芯片的心電方案可以實現在手指、手腕處使用金屬電極獲取高精度的 ECG 信號，適合手環、手表等便攜設備的應用。
表4
不同算法分類評估指標對比
Table4.
Comparison of classification effects of different algorithms
本文算法對于心率很快的單個心拍（心拍時長 0.5 s，心率每秒 120 次）單循環步驟為 16 次，循環次數為 10 次，總運算量約為 320 MIPS 定點運算，約為嵌入式低功耗心電芯片的 60%，基本功耗約為 72 μW。而遺傳算法和卷積編碼器，如表4 所示，其算法復雜度分別為O(n2)與O(n3)，故循環運算次數分別為 100 次與 1 000 次，對同樣的單個心拍進行運算時，總運算量可達 3 200 MIPS 定點運算和 32 000 MIPS 定點運算，故功耗也是低功耗心電芯片的 6 倍及 60 倍，達到了 720 μW 與 7 200 μW。所以復雜度高的算法雖然評估指標較高，但其運算占用資源相對來說也是非常大的，無法應用在便攜式心電設備中。
基于以上分析，本文算法是一種快速、復雜度低和輕量級的算法，可以實現對心電信號的準確檢測分類，算法指標滿足臨床檢測需求，可以適用于低功耗嵌入式處理芯片，可為便攜式心電監測設備心率失常實時監測提供科學有效的算法。
3 結論
本文提出了一種基于 ARMA 模型擬合的高效心律失常自動分類算法，可以很好地對心電序列進行類模型擬合與分析，在小波閾值去噪基礎上，通過雙微分方法與自適應閾值檢測 R 波并劃分心拍，對常見的心律失常信號進行了有效的特征表征。在使用數據集完成 10 層交叉訓練后，該算法能對有標記的心電信號進行自動檢測與分類。相較于目前心律失常分類算法大多存在復雜度較高且功耗較大的問題，該算法具有快速、復雜度低且準確率高的性能，完全滿足嵌入式設備的實時運算需求，并可實現基于便攜式設備的高效、高可靠、低功耗的心律失常分類運算，具有較高的現實意義與價值。
盡管本文算法在整體分類性能上相比于傳統或使用復雜網絡結構的心電分類算法有較大提升，但受限于 ARMA 模型的建模容量限制，對于更加多樣性、復雜化的心律失常分類精度還存在些許不足。未來可針對這些問題深入探究，以提出分類魯棒性更好、適用范圍更廣的算法。此外，本算法為一種通用算法，對于具體嵌入式硬件電路，可以根據芯片結構進一步優化。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
心律失常是心血管疾病（cardiovascular disease，CVD）的重要表現，由心律失常所引發的心血管疾病每年導致全球約 1 790 萬人死亡，占總死亡人數的 31%[1]，其早期的檢測與治療非常關鍵。穿戴式心電監測設備的飛速發展也正是為了實時監測個體心律失常信號特征，這一應用需求依賴高效便捷的心電檢測算法。心電圖（electrocardiogram，ECG）記錄了心臟至體表的電生理活動，通過其形態和節奏能有效評估心血管系統的狀態，是臨床檢測心律失常的必要手段之一。
心律失常的自動檢測（automatic detection of arrhythmia，ADA）是指使用計算機設備代替繁瑣耗時的人工 ECG 數據分析進行心律失常的診斷篩查[2]。隨著近年來機器學習的飛速發展，ADA 受到了更多學者的關注與研究。ECG 信號微弱不穩定的特性使得其特征提取成為研究的關鍵[3]。Raj 等[4]提出基于離散余弦變換的正交 stockwell 變換方法，利用支持向量機（support vector machine，SVM）分類器對來自 MIT-BIH 心律失常數據庫 16 類心律失常信號進行分類，但缺乏合適的 R 峰檢測算法，需要對特征集的大小進行最佳選擇以減少分類器的訓練時間。Mazaheri 等[5]提出在心電預處理后對 ECG 信號進行分段，提取其形態特征、頻域特征和非線性指標，利用元啟發式優化算法刪除多余或不相關的特征實現降維，但研究所應用的 NSGA II 特征選擇方法處理時間較長，難以實現實時穿戴式心電監測。
由于心跳的能量集中于 ECG 波形的 QRS 復合波中，R 波特有的陡峭特性有利于 ECG 的特征提取。Merah 等[6]設計了一種基于平穩小波變換的方法，省略了變換算法中每個級別的系數抽取，可更好地執行離群值檢測。Jung 等[7]提出利用稀疏判別提取重要的小波系數來自動檢測室性早搏的方法，通過定位適當的基準點以確定 ECG 信號中準確的 QRS 波群。基于經驗模式分解（empirical mode decomposition，EMD）的方法可以有效去除 ECG 信號中的基線漂移[8]，是心電預處理中的重要步驟。EMD 和離散小波變換（discrete wavelet transform，DWT）的組合方法為 ECG 降噪提供了更好的時間分辨率[9]。Lu 等[10]提出了一種基于 EMD 和 SVM 分類器集成的新方法，使用具有可變范圍編碼和動態搜索策略的改進遺傳算法優化統計特征子集，以獲得更好的分類精度，但缺少尋找分類精度的評價過程。
可見，針對 ECG 信號的特征提取方面使用時頻域及形態學特征的方法較為普遍。目前研究 ECG 信號自動分類的主流，一是對心電信號作較為復雜的變換提取相關特征參數，二是利用深度學習的端到端算法直接對原始信號進行分類[11-12]。但是，這些方法準確率的保證對原始 ECG 信號質量要求較高且計算過程復雜，消耗功耗大，且易受到 P 波、T 波等雜波的干擾，高效低功耗的算法研究少，因此現有研究成果難以應用于便攜式心電設備。
而直觀的自回歸移動平均模型（auto regressive moving average，ARMA）作為表征依賴數據最常用的模型之一，在語音處理、生物醫學、電力消耗預測和計量經濟學等領域有著悠久的傳統，但在心電分析方面的研究較少[13-15]。例如，Aguila 等[16]提出假最近鄰算法（false nearest neighbor，FNN）與 ARMA 模型相結合的時間序列建模對心電信號進行預測。Heurtefeux 等[17]采用 ARMA 模型，通過計算赤池信息準則對心電嵌入式系統進行了評價，建立了 p 階和 q 階以達到模型復雜度和準確性之間的權衡，相對較低的階數被認為在受約束的嵌入式系統中實現是可行的。
為了使用 ARMA 模型進行心電預測，需要選擇合適的模型參數。針對上述問題，本研究在 QRS 復合波檢測上利用雙微分法與自適應閾值準確判斷 R 峰位置，以此作為心拍分割基準點。在特征提取階段，應用 ARMA 模型將相關時序轉化為獨立時序，進行統計處理分析。上述步驟基于 ECG 信號特性，探求一種魯棒性與準確性更高的 ECG 分類算法，在極低的算法復雜度下可以對 ECG 實現高效分類，算法簡單，適用于穿戴式心電檢測，有利于對心血管疾病患者進行日常監護和預警，促進實時在線心電監護的發展。
1 實驗數據與方法
1.1 實驗數據
本文的數據來自 MIT-BIH 的心律失常數據庫（MIT-BIH Arrhythmia Database，MIT-BIH-AR），該數據庫包含從 1975 至 1979 年在心律不齊實驗室中的案例，選擇 48 條信號的雙導聯數據，每條 ECG 信號記錄時間長度均超過 30 min，采樣頻率：360 Hz，采樣分辨率：11 bit，采樣帶寬：0.1～40 Hz。考慮到房顫作為最常見的心律失常疾病，具有重要的研究意義，選擇 MIT-BIH 房顫數據庫（MIT-BIH Atrial Fibrillation Database，MIT-BIH-AF）中的房顫信號作為補充數據。將正常 ECG 信號與不同種類的心律失常信號交叉組成訓練集與測試集，以驗證房顫檢測算法的特異度。MIT-BIH-AR 數據庫與 MIT-BIH-AF 數據庫中正常與非正常心搏的具體數量、記錄數及持續時間如表1 所示。
表1
MIT-BIH 心電數據庫統計
Table1.
Statistics sorted out from MIT-BIH database
1.2 心電自動分類算法流程
本研究提出的算法流程如圖1 所示，首先對 ECG 信號進行預處理，利用小波閾值方法濾除噪聲、基線等干擾，采用雙微分法檢測 QRS 波群并提取 R 峰。根據 R 峰位置，在波峰前后各取一定數量的樣本點分割心電心拍，根據心拍注釋分類不同的心律失常信號，完成輸入數據的準備。利用 ARMA 模型與單變量特征選擇擬合對輸入的心拍數據進行系數擬合并降維，與心拍標注類型組成特征向量輸入到 SVM 和 K 近鄰（K-nearest neighbor，KNN）兩種分類器中，用于預測輸入樣本標簽。最后，評估模型識別的查準率、查全率等參數，為心律失常的自動檢測方法提供參考。
圖1
心電信號自動檢測流程圖
Figure1.
Flow diagram of automatic detection of ECG signal
1.3 心電預處理與 R 波提取
原始心電信號多數存在基線漂移（≤ 5 Hz）、肌電噪聲（10～300 Hz）和工頻噪聲（50 Hz）等干擾信號，心電預處理需要濾除噪聲干擾以提高后續檢測的準確性。本研究利用小波閾值去噪方法，對信號進行正交小波變換，得到一組小波分解系數，變換后的小波系數中含有大量有用 ECG 信號的數據，ECG 噪聲的圖像分布情況與小波分解系數分布相同，因此小波閾值去噪方法可有效濾除噪聲。由于小波系數大小不一，將絕對值小的設為 0，對絕對值大的采用軟閾值方法進行保存或壓縮，處理后的小波系數要更加平滑，得到其估計值，利用小波系數重構新的 ECG 信號，達到去噪效果。針對 ECG 信號的低頻特點，閾值λ的選擇根據 sqtwolog 準則確定[18]，能消除更多噪聲，其計算表達式與軟閾值函數表達式分別如式（1）和式（2）所示，其中n代表小波系數向量長度。
|
|
對濾波后的波形進行雙微分處理，以獲得心電波形的陡度與高度特性：通過在樣本點的的左右兩側設置閾值區間，分別尋找最大平均斜率（，）與最小平均斜率（，），交叉求差后求取最大值，如式（3）所示。設置微分閾值，其值的確定根據前 8 個峰值的平均值計算。雙微分處理后利用低通濾波消除雙峰現象，滑動窗口積分使得 ECG 絕對振幅增大，寬度根據經驗設為 17 個采樣點，此時預處理后的信號更便于檢測 R 峰。
|
檢測 R 波時選擇自適應閾值的設計，其選擇準則為：① 能在各種類型的記錄下良好地區分目標與非目標，達到較高準確率與可用度；② 選擇經驗參數需固定，針對不同心電記錄能跟隨波形實時穩健地變化；③ 自適應閾值可防錯檢、漏檢。例如利用“不應期”的含義，在檢測到一個 QRS 波后短時間內不再檢出 QRS 波，防止錯檢；同時可以進行波形檢測回溯，以防漏檢。基于以上三點，設計自適應機制如式（4）~（5）所示（THR0為低閾值，THR1為高閾值）。
|
|
其中，peak指檢測到的峰值，peak_buffer指存儲了當前峰值之前 8 個峰值，mean（）指均值。THR1_lim與THR0_lim為兩個經驗常數，代表兩個閾值變化的下限，在這里分別取值為 0.3 和 0.23。R 波提取的流程如圖2 所示，對 ECG 進行預處理檢測結果分別如圖3、圖4 所示。
圖2
R 波提取流程
Figure2.
Process of R wave extraction
圖3
ECG 預處理結果
Figure3.
Results of ECG preprocessing
圖4
“雙微分”處理后標記 R 波波峰提取示意圖
Figure4.
Schematic diagram of extraction of marked R-wave peaks after "double differential" processing
1.4 心拍分割
心律失常信號的異常主要體現在心電圖所蘊含波形形態中，完整的心電記錄包含多個心拍，直接全局處理不利于精細化識別，常通過截取單個心拍來分析判別，分割長度對于最終信號的分類準確度至關重要。根據先驗知識可知，波形異常主要表征在 PR 間期以及 QRS 波群的異常，因此在預處理的基礎上，通常以陡峭的 R 波波峰為關鍵點，前后分別截取 100 和 150 個樣本點共 250 個樣本點作為樣本心拍，心拍注釋采用 MIT-BIH 數據庫中的專家標注，根據標注中的心律失常類型進行后續的特征提取與分類工作。
本文共選取正常心拍（N）、左束支傳導阻滯心拍（L）、右束支傳導阻滯心拍（R）、室性早搏心拍（V）以及房顫信號（AF）五類心電信號。從 MIT-BIH-AR 數據庫的 48 條記錄中的 MLII 通道數據中進行心拍劃分，總計獲取 96 871 個樣本心拍，其中正常心拍（N）有 74 466 個，左束支傳導阻滯心拍（L）有 8 051 個，右束支傳導阻滯心拍（R）有 7 239 個，室性早搏心拍（V）有 7 115 個。從 MIT-BIH-AF 數據庫的 24 條記錄中的 MLII 通道數據中獲取 79 654 個心拍作為房顫（AF）類別。根據 MIT-BIH 數據庫的心拍注釋，對每類心電信號選取 2 500 個心拍，共 12 500 個心拍，利用 1.3 節 R 波提取后的位置進行心拍的分割。圖5 為所截取各類型心拍示意圖。
圖5
不同 ECG 心拍提取示意圖
Figure5.
Schematic diagram of different ECG heartbeat extraction
1.5 特征工程
ARMA 模型是目前應用最廣泛的線性平穩時間序列預測模型，其建模技術具有容易實現、處理速度快和能同時對多類 ECG 進行分類的優點，此外，它還能對觀測數據進行自主平滑和濾波，進一步達到去噪效果，從而估計出觀測數據的真值[19]。ARMA（p，q）的模型結構如式（6）所示。
|
蘊含在心電數據的取值大小可由 ARMA（p，q）模型的模型參數刻畫，心電數據的頻率特征可由 和 刻畫。基于模型適用于一維信號的平穩特性，可以通過對心電數據建立 ARMA（p，q）模型，提取擬合的模型參數作為每個樣本心拍的特征值。ARMA 模型可細分為 AR 模型和 MA 模型，通過對一維時序信號的觀測值反映出其基本特征。算法具體步驟如下：
步驟一：基于心電信號的 ARMA 線性輸出方程式可寫為：
|
方程（7）可分解為 AR 過程和 MA 過程：
AR 過程：
|
MA 過程：
|
其中，G(x)為 ECG 信號的時間序列；ai和bj分別為 AR 系數和 MA 系數；p和q分別為 AR 階次和 MA 階次；R(x)為 AR 誤差的序列；n(x)為 ARMA 預測誤差，即 ARMA 建模誤差。
步驟二：在建模中，階次p和q的選擇至關重要，本研究選擇最小描述長度（minimum description length，MDL）準則與敏感性函數來確定模型階次[20]。MDL 準則如式（10）所示：
|
式中，為建模誤差的方差，N為數據窗口的寬度，Np為模型的階次，在這里Np=p（AR 過程）或Np=q（MA 過程）。另外，MDL 的敏感性函數可以反映 ARMA 模型階次對擬合效果的影響，其表示見式（11）：
|
其中，是 MDL 的變化量，反映了模型階次Np對 MDL 的影響。
步驟三：通過 MDL 準則確定 ARMA 模型階數范圍在 2~4 之間，利用殘差檢驗驗證所選擇階數，如果殘差是隨機正態分布且不自相關，則說明所有有用信號都已被提取到 ARMA 模型中。同時，利用德賓-沃森統計量（Durbin-Watson Statistic，DWS）檢驗模型誤差項。
步驟四：根據 MDL 準則確定 ARMA（4, 2）模型系數作為心電信號系數。利用 R 語言中的 arima 函數擬合 ARMA 模型系數，并對得到系數矩陣進行降維。系數矩陣由擬合得到的 ARMA 系數xi和對應心拍標簽組成，使用主成分分析（principal component analysis，PCA）降維方法可有效降低由原始心拍數據維度過大造成的 ARMA 系數復雜問題。
對于線性 PCA 計算過程如下：首先，計算擬合系數均值并根據式（12）計算協方差矩陣 C。
|
式（13）給出了特征向量V和特征值E的計算。其中，特征值由高到低排序組成特征向量，各特征向量中，特征值最高的特征向量為數據的主成分。
|
線性 PCA 技術的重要優點是其簡單的幾何解釋和彌補后續分類算法過擬合的計算特性作為 ECG 的特征矩陣[21]。另外，為了對比計算，使用 kernel PCA 用非線性映射完成 PCA 的線性運算[22]。核主成分分析計算核矩陣的主特征向量，而非線性方法中計算協方差矩陣，核矩陣系數 的計算如式（14）所示。
|
其中，n是數據樣本的數量，k為核函數，例如高斯函數或多項式函數等。
整體流程如圖6 所示。
圖6
ARMA 模型擬合系數流程圖
Figure6.
Flow diagram of ARMA model fitting coefficient
殘差計算結果如圖7 所示，圖中分別顯示了殘差計算結果圖、QQ 分位圖以及自相關和偏自相關圖。可以發現殘差檢驗分布滿足正態分布，且信號圖形趨近于白噪聲波形，通過自相關和偏自相關圖可以發現模型階數選擇符合實際提取心拍。DWS 值越接近 2，則越說明殘差不存在一階相關性，本研究計算結果為 1.999 998，因此殘差結果符合要求。
圖7
ARMA 模型階數殘差檢驗結果
Figure7.
The results of ARMA model order residual test
利用選定階數預測心拍序列，可得圖8 所示波形，此時模型擬合的仿真信號與原始信號最為接近。故本研究選用 ARMA（4，2）模型，只要確定各參數的取值，對應的模型也就完全確定，將其用來擬合數據、考察數據內在的統計特征以及作最佳預測時，在數學上的分析處理是便捷高效的。模型序列的譜密度為有理譜密度，而連續譜密度可以用有理譜密度來逼近，并能達到理想的近似程度。
圖8
ARMA (4, 2)模型預測心拍序列結果
Figure8.
Prediction results of heart beat sequence by ARMA (4, 2) model
1.6 機器學習方法分類與比較
本研究使用了 SVM 與 KNN 兩種分類器。SVM 是一種監督學習的二分類模型，其主要目的是找到一種超平面對樣本數據進行分割，從而轉化為求解凸二次規劃的問題。然而在大部分的分類任務中，單純的線性關系完成分類實例的準確分割是較為困難的，SVM 模型利用非線性核函數將分類任務從線性不可分的低維空間轉換到線性可分的高維空間，然后通過算法分析找出最佳超平面，最后返回到原始空間以獲得初始空間的解[23]。
KNN 通常基于某種距離度量找出訓練集中與其最靠近的 k 個訓練樣本，利用“投票法”選擇 k 個樣本中出現最多的類別標記作為預測結果。它是“懶惰學習”的代表，訓練時間開銷為零，雖然其結構與運算簡單，但它的泛化錯誤率不超過貝葉斯最優分類器錯誤率的兩倍，因此在許多分類任務中的應用也是非常廣泛的[23]。
本研究選擇使用 K 次交叉驗證方法（K=10）無重復地抽樣劃分輸入數據作為訓練集和驗證集，訓練集用于訓練模型，驗證集評估模型性能，以此來降低選擇模型過擬合的現象。另外，通過使用網格化搜索函數來優化分類模型，提高模型的穩定性以及泛化性能。
2 實驗評估與分析
2.1 評估指標
本文利用查全率（以符號表示）、查準率（以符號表示）、評估精度分數來評估心電信號的分類效果，由真陽性（true positive，TP）、假陰性（false negative，FN）、假陽性（false positive，FP）的數量以及數據實例決定，指標計算公式如式（15）~（17）所示：
|
|
|
在多分類問題中，對于第i類心律失常樣本，TP 表示該類被正確分類的樣本數量，FN 表示第i類的樣本被錯誤地分類為第j類的數量，FP 表示第j類的樣本被錯誤地分類為第i類的數量。為了便于表達，以符號表示五種類型 ECG 的平均查全率，以符號表示五種類型 ECG 的平均查準率，以符號表示五種類型 ECG 的平均指數。
2.2 實驗結果分析
部分正常心電信號樣本心拍的 ARMA 模型擬合系數降維后如表2 所示，表中的、系數分別表征了輸入心拍的取值大小以及頻率特征，且結構簡單。將其與心電類型標簽組合作為分類器模型的輸入，表3 為各類別心電信號的分類結果統計，KNN 分類器對 MIT-BIH AR 數據集分類得到查全率、查準率以及指數分別為 98.1%、97.7% 以及 97.9%，而 SVM 分類器的結果分別為 98.2%、98.4% 以及 98.3%。在 M 類分類問題中，混淆矩陣可以顯示分類器獲得分類結果的詳細分布，而且是一個 M 維方陣。圖9 為 KNN 與 SVM 分類算法結果混淆矩陣。從結果可以看出，文中提出的預處理以及特征工程算法在 MIT-BIH 數據庫上有較好的分類性能。從分類器特性上看，KNN 對噪聲有較好的魯棒性，但模型有欠擬合或過擬合傾向，SVM 在建立了一個成功的模型基礎上，對新樣本的分類較為準確，目前應用更加廣泛。對比 1.5 小節中兩種 PCA 降維技術，它們應用在 KNN 分類器上的效果差異并不大，這是由于降維操作實際上是對數據進行了一定程度的近似，丟失了部分信息所導致的。從綜合衡量查全率與查準率的指數看出，SVM 分類器要優于 KNN 分類器。
表2
部分正常心電信號樣本心拍的 ARMA 模型擬合系數
Table2.
ARMA model fitting coefficient of cardiac beats in partial normal ECG signal samples
表3
機器學習方法分類效果統計
Table3.
Classification accuracy statistics of machine learning methods
圖9
KNN 分類算法與 SVM 分類算法結果混淆矩陣
Figure9.
Confusion matrix of SVM and KNN classification algorithm results
在處理 ECG 數據方面，本研究避免了特異性處理特殊心律失常的復雜步驟，采用簡單的過濾和標記技術實現算法的低復雜特性。另外 ARMA 模型應用 R 語言的 arima 函數來實現，其數據收集與回歸過程通過測試模型擬合時間驗證得到計算復雜度為O(n)，可以近似看作線性模型。從結果來看，提出方法的查準率與查全率均可以達到很好的分類結果，且特有的低復雜度適用于便攜式的心電實時監測。
表4 顯示了近年來關于心律失常自動檢測工作[10,24-28]與本文提出算法在計算復雜度與評估指標方面的對比，所得分類效果均基于 MIT-BIH 數據庫。相較于計算復雜度相同的算法，本文所提出算法在各個評估指標上均具有顯著優勢；相較于 Oh 等[24]、Yildirim 等[25]提出的基于長短期記憶人工神經網絡（long short-term memory，LSTM）的復雜模型，本文提出算法的復雜度更低。另外，考慮算法的實際應用方面，選擇高精度低功耗 ECG 傳感芯片 KS1081 作為參考芯片，計算不同算法在相同運算情況下的功耗。KS1081 芯片的單通道功耗約為 120 μW，平均執行速度為 532 MIPS（每秒處理的百萬級語言指令數）32 位定點運算，可滿足 16 位采樣精度數據運算需求，且可實現在長時間連續信號監測情況下提升設備續航時間的功能。故應用該芯片的心電方案可以實現在手指、手腕處使用金屬電極獲取高精度的 ECG 信號，適合手環、手表等便攜設備的應用。
表4
不同算法分類評估指標對比
Table4.
Comparison of classification effects of different algorithms
本文算法對于心率很快的單個心拍（心拍時長 0.5 s，心率每秒 120 次）單循環步驟為 16 次，循環次數為 10 次，總運算量約為 320 MIPS 定點運算，約為嵌入式低功耗心電芯片的 60%，基本功耗約為 72 μW。而遺傳算法和卷積編碼器，如表4 所示，其算法復雜度分別為O(n2)與O(n3)，故循環運算次數分別為 100 次與 1 000 次，對同樣的單個心拍進行運算時，總運算量可達 3 200 MIPS 定點運算和 32 000 MIPS 定點運算，故功耗也是低功耗心電芯片的 6 倍及 60 倍，達到了 720 μW 與 7 200 μW。所以復雜度高的算法雖然評估指標較高，但其運算占用資源相對來說也是非常大的，無法應用在便攜式心電設備中。
基于以上分析，本文算法是一種快速、復雜度低和輕量級的算法，可以實現對心電信號的準確檢測分類，算法指標滿足臨床檢測需求，可以適用于低功耗嵌入式處理芯片，可為便攜式心電監測設備心率失常實時監測提供科學有效的算法。
3 結論
本文提出了一種基于 ARMA 模型擬合的高效心律失常自動分類算法，可以很好地對心電序列進行類模型擬合與分析，在小波閾值去噪基礎上，通過雙微分方法與自適應閾值檢測 R 波并劃分心拍，對常見的心律失常信號進行了有效的特征表征。在使用數據集完成 10 層交叉訓練后，該算法能對有標記的心電信號進行自動檢測與分類。相較于目前心律失常分類算法大多存在復雜度較高且功耗較大的問題，該算法具有快速、復雜度低且準確率高的性能，完全滿足嵌入式設備的實時運算需求，并可實現基于便攜式設備的高效、高可靠、低功耗的心律失常分類運算，具有較高的現實意義與價值。
盡管本文算法在整體分類性能上相比于傳統或使用復雜網絡結構的心電分類算法有較大提升，但受限于 ARMA 模型的建模容量限制，對于更加多樣性、復雜化的心律失常分類精度還存在些許不足。未來可針對這些問題深入探究，以提出分類魯棒性更好、適用范圍更廣的算法。此外，本算法為一種通用算法，對于具體嵌入式硬件電路，可以根據芯片結構進一步優化。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
