心音分割指對所獲取的心音信號按心動周期對收縮期、舒張期等進行分隔，是進行心音分類前的關鍵步驟。針對不依賴心電圖對心音信號直接分割準確度有限的難題，提出了一種基于持續時間隱馬爾可夫模型的心音分割算法。首先對心音樣本進行位置標注；然后采用自相關估計法對心音的心動周期持續時間進行估計，通過高斯混合分布對樣本的狀態持續時間進行建模；接著通過訓練集信號對隱馬爾可夫模型進行優化并建立基于持續時間的隱馬爾可夫模型（DHMM）；最后使用維特比算法對心音狀態進行回溯得出 S1、收縮期、S2、舒張期。使用 500 例心音樣本對本文算法性能進行測試，平均評估精度分數（F1）為 0.933，平均靈敏度為 0.930，平均精確率為 0.936。同其他算法相比，本文算法各項性能指標均有明顯提升，證實了該算法具有較高的魯棒性和抗噪聲性能，為臨床環境下所采集心音信號的特征提取與分析提供了一種新方法。
引用本文： 奎皓然, 潘家華, 宗容, 楊宏波, 粟煒, 王威廉. 基于持續時間隱馬爾可夫模型的心音分割算法. 生物醫學工程學雜志, 2020, 37(5): 765-774. doi: 10.7507/1001-5515.201911061 復制
引言
心音圖（phonocardiogram，PCG）在評估心血管異常的初始診斷篩選中起著至關重要的作用[1]。由于獲取 PCG 成本低、操作簡單，它在遠程健康診斷系統中有著重要的應用。它也被用來補充基于心電圖（electrocardiogram，ECG）的心臟診斷，以檢測心血管異常。隨著電子聽診器的問世，采集較高質量心音信號的問題得到解決[2]。由于與特定心臟病相關的雜音發生在心臟周期的某些時間點，因此準確地對心音圖的狀態（S1、收縮期、S2、舒張期）進行分割[3]是進一步分析心音信號的前提。
心臟雜音（murmur）是醫生對先天性心臟病進行初步診斷的重要信息[4]，根據雜音出現的位置可分為收縮期雜音與舒張期雜音。由于使用電子聽診器進行心音采集時容易引入其他非病理性雜音，例如聽診器與皮膚的摩擦音、患者的呼吸音、環境噪聲等[5]，這使得心音分割的過程變得十分困難。
現階段心音分割的主要算法包括：① 基于特征幅度閾值的分割方法，例如：Vepa 等[6]使用多層小波分解系數提取出基于能量的特征，通過閾值判定達到了 0.84 的準確率；Moukadem 等[7]對心音樣本進行 S 變換計算局部頻譜的香農能量，通過局部自適應閾值細化狀態邊界，使用 k 最近鄰（k-nearest neighbor，KNN）算法進行心音狀態分類，準確率可達 0.92。② 基于神經網絡的分割方法，例如：Sepehri 等[8]使用多層感知機（multilayer perceptron，MLP）對心音信號的短時能量進行訓練并驗證，在 60 個心音樣本測試集上獲得了 0.93 的準確率；Messner 等[9]使用深層遞歸神經網絡（deep recurrent neural network，DRNN）對心音頻譜進行序列檢測，獲得了 0.95 的準確率。③ 基于包絡的分割方法，例如：Yan 等[10]使用局部極值點對心音特征波形進行分段，取得了 0.93 的準確率；Sharma[11]提出基于多尺度的 Hilbert 包絡分割方法，在小樣本心音測試集下準確率可達 0.94。④ 基于隱馬爾可夫模型（hidden Markov model，HMM）的分割方法，例如：Gill 等[12]使用 HMM 對心音同態包絡進行建模，在 44 個測試集上取得了 0.97 的準確率；Oliveira[13]使用邏輯回歸函數作為初始概率建立 HMM 模型進行心音狀態分割，平均準確率達到 0.93。
上述研究中，基于特征幅度閾值的分割方法對心音進行維度變換，通過特征閾值進行心音分割。在理想狀況下，此類算法可達到較高準確度，其缺點為因非病理性雜音的干擾會影響特征的提取，不適用于臨床環境。基于神經網絡的分割方法主要使用已標記的心音信號對網絡進行訓練，從而得到狀態序列。由于神經網絡需要大量數據進行監督學習，針對心音小規模數據集進行訓練易造成過擬合現象，普適性較差。基于包絡的分割方法通過提取心音包絡進行心音分割，其特點為運算效率高，并有一定的抗噪聲性能。其缺點為由于形態學分割的局限性，分割精度遜色于其他幾類算法。基于 HMM 的方法是通過觀測序列與隱含狀態序列間的關系推斷出存在的合理狀態序列[14]。由于采集心音信號時心臟的狀態是未知的，但采集到的心音信號是可觀察的，所以基于 HMM 的方法比較符合心音分割任務的要求。現階段基于 HMM 的分割算法大多是對數據進行非監督學習，但這些算法使用的 HMM 存在一個問題，即沒有很好地對狀態的持續時間進行建模，所造成的現象是心音轉移至下一狀態的概率不受當前狀態持續時間的影響，這導致了基于 HMM 的算法用于心音狀態分割不夠準確。
針對上述問題，考慮到基于 HMM 算法對于心音狀態建模的優勢，本文使用高斯混合模型對心動周期持續時間進行建模，構建基于持續時間的隱馬爾可夫模型（duration hidden Markov model，DHMM），以探求一種魯棒性與準確性更高的心音分割算法。
1 模型與方法
1.1 數據來源
本文研究的心音數據來源于研究團隊采集構建的心音樣本庫。該心音樣本庫志愿者年齡分布于 6 個月至 18 歲之間，所有志愿者或其監護人均簽署知情同意書。本研究經云南大學醫學院倫理委員會審查通過。樣本庫中所有心音信號所屬分類疾病，均通過臨床專家使用超聲心動儀（Acuson Cypress，德國西門子股份有限公司）確診。采集所使用的心音傳感器為美國 ThinkLabs 公司的 THE ONE，其輸出為模擬音頻信號，采用實驗室自行研發的采集設備，可同步記錄心音心電信號，采樣頻率為 5 000 Hz。該采集設備能將模擬信號轉換為數字信號，便于儲存和處理。
1.2 心音信號的具體分割流程
本文心音分割算法的具體流程如圖1 所示，整個心音分割流程可概括為以下 7 個步驟。
圖1
心音分割算法流程圖
Figure1.
Flow chart of heart sound segmentation algorithm
（1）使用 R 峰檢測器與 T 波探測器對同步采集心音心電的信號進行位置標注；
（2）將標記的信號分為訓練集與測試集；
（3）對訓練集信號與測試集信號提取心動周期等特征；
（4）對訓練集信號進行高斯混合模型建模，為 DHMM 提供參數 
、
與 pj（d）；
（5）構建 DHMM，并使用訓練集信號對模型參數 λ 進行訓練；
（6）使用訓練后的 DHMM 對測試集信號進行建模，并通過改進的維特比算法尋找最優狀態序列；
（7）使用測試集信號心音心電位置的對應關系評估模型性能指標。
其中，步驟（1）和（2）中對測試集信號心音心電進行位置標注，其目的是為步驟（7）評估模型性能指標提供參考依據。在實際使用中，對于沒有任何參考標記的心音信號，經特征提取后能直接由訓練完成的 DHMM 進行心音分割。
1.3 信號標注
為了訓練基于 DHMM 的心音數據分割算法，必須使用外部參考信號標記 PCG 記錄中 S1 和 S2 的起始位置。使用 ECG 作為心音狀態 S1 和 S2 的參考標記位置，其中 ECG 的 R 峰對應 PCG 信號中 S1 的起始位置，而 S2 發生在 ECG 的 T 波末端處。因此首先需要在 ECG 中標記 R 峰與 T 波末端位置。本文采用了基于拋物線擬合的 R 峰檢測算法[15]與 Laguna 等[16]所提出的 T 波探測器，分別對 200 例訓練信號進行標記。訓練樣本中包括 80 例異常信號和 120 例正常信號，所有信號均來源于 1.1 節中的心音樣本庫。
由于 S1 的起始位置與 ECG 的 R 峰重合[17]，將一個 R 峰至下一個 R 峰之間 S1 平均持續時間段標記為 S1。S2 出現在 ECG 的 T 波末端附近[17]，且 S2 的振幅在 T 波末端附近將達到最大值。因此，設置搜索窗口在 T 波末端附近尋找 PCG 信號 Hilbert 包絡中的最大峰值，以此峰值作為 S2 中心，將其周圍 S2 平均持續時間段標記為 S2。收縮期與舒張期介于 S1 與 S2 的相對位置之間，將 S1 到 S2 之間的時間段標記為收縮期，S2 到下一個 S1 之間的時間段標記為舒張期。完成標記后的 PCG 如圖2 所示。
圖2
R 峰與 T 波末端分別作為定義 S1 與 S2 的參考位置
Figure2.
R peak and T wave terminal are the reference positions of S1 and S2 respectively
1.4 心動周期的提取
心音信號具有短時平穩性，研究中普遍認為在 10～30 ms 內心音信號可視作平穩[18]。采用自相關法[19]能較好地提取心動周期的長度，由于平滑的心音包絡更容易反映心動周期的特征，因此對心音信號的 Hilbert 包絡求自相關。將心音信號的 Hilbert 包絡記為
，其中 N 表示信號的總長度，則自相關計算公式可以表示為
 |
提取出的自相關包絡如圖3 所示，在圖中所示兩個紅點之間的間距為心動周期的持續時間（Dcycle），第一個紅點與綠點之間的間距稱為一個完整的收縮期（systole），其成分包括了 S1 與收縮期兩個狀態。其持續時間（Dsystole）對應公式如下
圖3
自相關法提取心動周期
Figure3.
Cardiac cycle was extracted by autocorrelation method
 |
完整的舒張期（diastole）包括 S2 與舒張期兩個狀態，其持續時間（Ddiastole）有
 |
由式（2）與式（3）可以得到心音各個狀態持續時間的約束關系，使得接下來對心動周期持續時間的高斯建模更加符合真實的心音持續時間分布。
1.5 心動周期持續時間的高斯建模
由于在平靜狀態下人類的心率大多介于 60～100 bpm 之間，則心動周期持續時長 heartcycle_duration 為
 |
其中 heartrate 代表心率。由式（4）可推出一個心動周期時長為 600～1 000 ms。圖4 是對 200 例訓練集心音信號單一心動周期的各狀態持續時間的統計直方圖，其中 n1、n2、n3、n4 分別對應 S1、收縮期、S2、舒張期的心音狀態持續時間出現次數。由于 n1、n2、n3 的持續時間較為集中，在相同的時間段內出現頻次有重疊部分，故繪制了對應狀態的包絡以便觀察。從圖中可看出，各個狀態的持續時間近似于正態分布。
圖4
對應狀態持續時間的出現頻次
Figure4.
The occurrence frequency of state duration
由于一個心動周期有四個狀態，故對心動周期的持續時間進行高斯混合模型建模[20]。定義 X = {x1，x2，
，xN}為觀測的隨機向量，又稱為不完全數據（incomplete-data），其中 x1，x2，
，xN 指觀測到的心音持續時間。（X，Z） = {（x1，z1），（x2，z2），
，（xN，zN）}稱為完全數據。Z 為隱隨機變量數據，表示對應的樣本 x 是屬于哪一個高斯分布。根據之前的分析，心音四個狀態的持續時間近似于正態分布。將心音各狀態的集合表示為 C（condition），有 C = （C1，C2，C3，C4），其中 C1 表示 S1，C2 表示收縮期，C3 表示 S2，C4 表示舒張期。可以確定 Z∈（C1，C2，C3，C4）。θ = {p1，p2，
，pk μ1，μ2，
，μk ∑1，∑2，
，∑k}為所需估計的參數，由于各狀態高斯模型的均值 μ 與方差∑可以通過訓練集樣本統計得出，故所需確定參數為 θ = {p1，p2，p3，p4}，其中 p1、p2、p3、p4 表示持續時間為 x 時分別屬于 C1、C2、C3、C4 高斯分布的概率。由此，可以得出離散概率分布如表1 所示
表1
高斯混合模型離散概率分布
Table1.
Discrete probability distribution of Gaussian mixture model
高斯混合模型的概率密度函數為
 |
式（5）中 N（x | μk，∑k）代表著 x 服從于第 Ck 類正態分布。式中 μk、∑k 通過訓練信號統計得出。使用最大期望算法（expectation-maximization algorithm，EM）計算模型最優參數 pk，EM 算法的遞推表達式為
 |
將式（6）中的 
稱為 Q 函數，有
 |
通過式（6）和（7），可以將其整理為
 |
在式（8）中，需要設定初始參數 θ（1），由于在采集過程中無法確定從哪個狀態開始采集，假設初始狀態等概率出現，則 
。根據遞推式可得到任意 t 時刻 θ（t）的值，因此只需求解未知參數 θ。使用拉格朗日乘數法確定 θ，得到遞推參數 pk（t + 1）：
 |
可由式（9）推出 p（t + 1） = {p1（t + 1），p2（t + 1），p3（t + 1），p4（t + 1）}，由此可以遞推各個時刻的最優參數。
1.6 心音信號的 DHMM 建模
DHMM 與傳統 HMM 最大的區別就是引入了持續時間的約束轉移條件。雖然傳統 HMM 很適合對心音信號狀態分割進行建模，但由于傳統 HMM 的限制性，所涉及的馬爾可夫過程是無記憶性的，處于下一狀態的概率僅取決于當前時刻的狀態，而與以前的狀態無關。這意味著心音不論在當前狀態持續多長時間，轉移到新狀態的概率都是相同的，導致每個狀態的持續時間呈幾何分布[21]，這顯然不符合基于時域的心音信號。雖然可以通過引入更多數量的子狀態來處理這一限制，但無形中引入了更多的參數，使模型變得復雜化，而且這只是緩解了癥狀，沒有從根本上解決問題。由此，引入持續時間密度 pj（d），其中 d（duration）表示狀態持續時間，概率密度函數 pj（d）表示在狀態 j 持續時間為 d 所發生的概率。j∈（C1，C2，C3，C4），故引入 pj（d）后，可以將模型簡化為四個狀態。
 |
其中 P（d）由本文第 1.5 節的高斯混合模型確定。
記 t 時刻模型所處的狀態為 ct，則有 ct∈（C1，C2，C3，C4）。設觀測值序列集合為 O，O = {O1，O2，O3，
，Ot}。定義 A 為狀態轉移概率矩陣，則
 |
aij 表示從狀態 Ci 到狀態 Cj 轉移時的轉移概率[22]。由于心音是規律性振動產生的準周期性生物信號，且總是按照 S1、收縮期、S2、舒張期循環出現。則狀態轉移概率集合
 |
DHMM 觀測序列與隱含狀態序列對應關系如圖5 所示，其中狀態間轉移概率由概率密度函數 pj（d）控制。
圖5
基于心動周期狀態的隱馬爾可夫模型
Figure5.
HMM based on cardiac cycle state
定義輸出觀測值概率集合 B = {bj（Ot）}，其中 bj（Ot）是在狀態 j 時觀測值符號 Ot 輸出的概率。
定義初始狀態概率集合為 π = {πi}，其中
 |
由于實際采集心音過程中可從心動周期的任何一個狀態開始采集，因此
 |
在使用 DHMM 之前，需要通過訓練集數據對其模型參數 λ 進行估計
 |
其中參數 
、
分別為觀測序列處于狀態 j 持續時間的平均值和方差，由 1.5 節中的高斯混合模型確定。
1.7 維特比算法尋找最優序列
心音分割的目的是找出所有已知觀測序列所對應的狀態序列，在已知模型參數 λ 的前提下，可以使用維特比算法（Viterbi algorithm）尋找心音信號的最佳狀態序列。
維特比算法是一種前向算法，用于計算前向概率 
。
表示時刻 t 時 ct 處于狀態 j 的瞬時概率。由于 DHMM 引入了持續時間的概念，故對傳統的維特比前向算法[21]進行改進使其適應 DHMM，改進后的前向概率 
表示為：
 |
給定持續時間 d，且前一狀態在時刻 
時為 i，則 ct 在狀態 j 結束的概率為
 |
其中 
是上一個狀態 i 在 
時刻結束的概率。
是在狀態 j 時持續 
到 t 期間觀測值輸出的概率。找到使式（16）最大化的參數 d 與前一狀態 i 來計算前向概率 
[23]：
 |
持續時間 d 限制在 
內且小于 t，其中 
和 
表示持續時間為 d 時所對應高斯混合模型的均值與方差。對于每個 t，存儲使式（17）最大化的持續時間 d 和前一狀態 i：
 |
在回溯之前，尋找最優狀態 j，使得結束狀態 
處于 j 的概率最大：
 |
實現回溯算法的具體程序流圖如圖6 所示。在圖6 中，首先初始化參數 
，將當前狀態設為觀測序列末尾所處的狀態。判定條件為 t > 1，目的為從序列末尾狀態回溯至序列首狀態。
讀取出 t 時刻最優狀態 
的持續時間 d*。
實現了將 
時刻至 
時刻之間的狀態標記為 
。
將 
時刻的狀態調整為 
的前一個最優狀態。
是下一次程序迭代的前提。當程序結束時，得到最優狀態序列集合 C* = {c1*，c2*，c3*，
，cT*}。
圖6
回溯算法程序流程圖
Figure6.
Program flow diagram of backtracking algorithm
2 實驗結果與分析
2.1 實驗環境與數據說明
實驗所用電腦 CPU 型號為 Intel（R）Core（TM）i7-8750H CPU @ 2.20 GHz，RAM 為 16.0 GB，仿真軟件為 MATLAB2019b。
本文算法與對比算法均使用本文測試集數據進行測試，測試集數據來源于 1.1 節所提及的心音數據庫。由于需要衡量各個算法在針對不同病理心音信號進行分割時的性能，故準備了 4 個獨立的測試集共 500 例心音信號對其進行評估，分別為 290 例正常（Normal）信號、70 例房間隔缺損（atrial septal defect，ASD）信號、70 例動脈導管未閉（patent ductus arteriosus，PDA）信號與 70 例室間隔缺損（ventricular septal defect，VSD）信號。
2.2 評價指標
為了評估模型對于分割心音信號的能力，使用 PCG-ECG 對照法評估模型在測試集中準確定位 S1 與 S2 的能力。S1 與 S2 的參考位置在 1.3 節中已經提及，所有參考位置均使用 R 峰檢測器與 T 波探測器進行標記。由于 R 峰檢測器與 T 波探測器標記位置時存在 100 ms 的檢測誤差。故若分割 S1 的起始位置在 ECG R 峰的 100 ms 內，則 S1 被標記為正確分段（TP）。同樣地，若 S2 中心位置在對應 T 波末端的 100 ms 內，則 S2 聲音被標記為正確分段（TP）。其他所有未被正確分段的位置標記為錯誤分段（FP）。圖7 為心音心電分割狀態對照圖，其中狀態序列（圖中紅色虛線所示）有 4 個狀態，狀態 1、2、3、4 分別對應 S1、收縮期、S2、舒張期。該圖可從心音、心電的位置對應關系客觀地評估模型分割性能指標。
圖7
心音-心電圖評估模型分割結果
Figure7.
Segmentation result of PCG-ECG evaluation model
對于心音分割這一任務，由于并不存在真正的負樣本，所以將靈敏度（sensitivity）定義為：
 |
其中 TP 表示分割正確的數量，All_signal 表示所有 S1、S2 位置的總數。recall 表示召回率，在數值上與 sensitivity 相等。
精確率（precision）定義為：
 |
其中 FP 表示分割錯誤的數量。
使用評估精度分數 F1[24]綜合評價 sensitivity 與 precision 指標：
 |
2.3 實驗結果
本文選取引言中提及的基于 Hilbert 的算法[11]與基于傳統 HMM 算法[12]作為對比算法。這兩類算法都有一定的抗噪聲性能，且運行效率較高，對含有部分非病理性噪音干擾的心音信號進行分割時均有較好表現。表2～5 反映了各算法在處理不同類型心音信號時的性能表現。
表2
對于正常心音信號各算法的評價指標（
）
Table2.
Evaluation index for each algorithm of normal heart sound signal（
）
表3
對于房間隔缺損心音信號各算法的評價指標（
）
Table3.
Evaluation index for each algorithm of ASD heart sound signal（
）
表4
對于動脈導管未閉心音信號各算法的評價指標（
）
Table4.
Evaluation index for each algorithm of PDA heart sound signal（
）
表5
對于室間隔缺損心音信號各算法的評價指標（
）
Table5.
Evaluation index for each algorithm of VSD heart sound signal（
）
圖8 為使用本文算法對各類心音進行分割的結果，該結果由程序自動生成。從圖中可觀察出對于病理性的心音信號，DHMM 仍能高精度地對其進行心音狀態分割。
圖8
對各類心音信號的分割結果
Figure8.
Segmentation result of various heart sound signals
2.4 討論
本文使用 4 個獨立的數據集對算法性能進行評估，評估指標結果均由程序自動運算得出，可以最大限度地降低主觀因素的影響。對比以往文獻所使用的測試數據與評估方法，本文實驗結果可以更為客觀地反映出各算法在處理不同類型心音信號時的性能表現。
從 2.3 節的表中可看出，對正常心音信號進行分割時，三種算法都能實現高精度的心音分割。其中文獻[12]基于傳統 HMM 心音分割算法使用同態包絡提取特征序列作為 HMM 的觀測序列，然后通過無監督學習對心音狀態進行分割。通過觀察表中數據可以發現，該算法總體標準差要大于其他兩種算法。出現這種現象的原因是此算法引入了大量子狀態用于緩解 HMM 無記憶性的特點，心音分割性能很大程度受限于轉移概率矩陣參數，這使得該算法分割性能在處理不同樣本時會產生較大變化。值得注意的是，該算法準確率同引言中描述相差較大，推測是由于文獻中僅使用 44 例數據進行測試，且分割精度評價指標不如本文嚴格。文獻[11]基于 Hilbert 包絡的心音分割算法運行效率最高，在對正常心音信號進行分割時性能表現穩定，但此算法在處理病理性心音信號時，分割精度明顯下滑。這是由形態學進行心音分割時的局限性所導致的，因為程序很難實現對復雜性病理心音信號進行包絡邊緣檢測。本文提出的基于 DHMM 的算法針對傳統 HMM 算法的缺點進行改進，對心音信號持續時間進行高斯混合模型建模，并將模型狀態數精簡為 4 個。從表中可看出，本文算法在對不同類型心音樣本進行處理時性能表現均優于其他兩種算法。與基于傳統 HMM 的算法相比，本文算法運行效率、泛化能力以及分割精度都有了大幅提升。基于心音狀態循環周期反復出現的特點，本文構建 DHMM，使得算法具有較高的魯棒性和抗噪聲性能，對復雜性病理心音信號與正常心音信號均能適用。
3 結論
基于 DHMM 的心音分割方法可以很好地對心音信號的準周期性、短時平穩性以及時域心音狀態分布進行建模。在使用訓練集信號完成 DHMM 訓練后，該算法能對無任何位置參考的心音信號進行心音狀態分割。由于算法具有較高的魯棒性與抗噪聲性能，更加符合臨床環境的使用需求。本算法在整體分割性能上相比傳統 HMM 分割算法有大幅提升，但受限于 DHMM 在序列開始和結束時對觀測序列建模能力的限制，對起始狀態與終止狀態的分割精度還存在少許不足。后續可對該問題進行研究，進一步提高算法精度，擴大算法適用范圍。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
心音圖（phonocardiogram，PCG）在評估心血管異常的初始診斷篩選中起著至關重要的作用[1]。由于獲取 PCG 成本低、操作簡單，它在遠程健康診斷系統中有著重要的應用。它也被用來補充基于心電圖（electrocardiogram，ECG）的心臟診斷，以檢測心血管異常。隨著電子聽診器的問世，采集較高質量心音信號的問題得到解決[2]。由于與特定心臟病相關的雜音發生在心臟周期的某些時間點，因此準確地對心音圖的狀態（S1、收縮期、S2、舒張期）進行分割[3]是進一步分析心音信號的前提。
心臟雜音（murmur）是醫生對先天性心臟病進行初步診斷的重要信息[4]，根據雜音出現的位置可分為收縮期雜音與舒張期雜音。由于使用電子聽診器進行心音采集時容易引入其他非病理性雜音，例如聽診器與皮膚的摩擦音、患者的呼吸音、環境噪聲等[5]，這使得心音分割的過程變得十分困難。
現階段心音分割的主要算法包括：① 基于特征幅度閾值的分割方法，例如：Vepa 等[6]使用多層小波分解系數提取出基于能量的特征，通過閾值判定達到了 0.84 的準確率；Moukadem 等[7]對心音樣本進行 S 變換計算局部頻譜的香農能量，通過局部自適應閾值細化狀態邊界，使用 k 最近鄰（k-nearest neighbor，KNN）算法進行心音狀態分類，準確率可達 0.92。② 基于神經網絡的分割方法，例如：Sepehri 等[8]使用多層感知機（multilayer perceptron，MLP）對心音信號的短時能量進行訓練并驗證，在 60 個心音樣本測試集上獲得了 0.93 的準確率；Messner 等[9]使用深層遞歸神經網絡（deep recurrent neural network，DRNN）對心音頻譜進行序列檢測，獲得了 0.95 的準確率。③ 基于包絡的分割方法，例如：Yan 等[10]使用局部極值點對心音特征波形進行分段，取得了 0.93 的準確率；Sharma[11]提出基于多尺度的 Hilbert 包絡分割方法，在小樣本心音測試集下準確率可達 0.94。④ 基于隱馬爾可夫模型（hidden Markov model，HMM）的分割方法，例如：Gill 等[12]使用 HMM 對心音同態包絡進行建模，在 44 個測試集上取得了 0.97 的準確率；Oliveira[13]使用邏輯回歸函數作為初始概率建立 HMM 模型進行心音狀態分割，平均準確率達到 0.93。
上述研究中，基于特征幅度閾值的分割方法對心音進行維度變換，通過特征閾值進行心音分割。在理想狀況下，此類算法可達到較高準確度，其缺點為因非病理性雜音的干擾會影響特征的提取，不適用于臨床環境。基于神經網絡的分割方法主要使用已標記的心音信號對網絡進行訓練，從而得到狀態序列。由于神經網絡需要大量數據進行監督學習，針對心音小規模數據集進行訓練易造成過擬合現象，普適性較差。基于包絡的分割方法通過提取心音包絡進行心音分割，其特點為運算效率高，并有一定的抗噪聲性能。其缺點為由于形態學分割的局限性，分割精度遜色于其他幾類算法。基于 HMM 的方法是通過觀測序列與隱含狀態序列間的關系推斷出存在的合理狀態序列[14]。由于采集心音信號時心臟的狀態是未知的，但采集到的心音信號是可觀察的，所以基于 HMM 的方法比較符合心音分割任務的要求。現階段基于 HMM 的分割算法大多是對數據進行非監督學習，但這些算法使用的 HMM 存在一個問題，即沒有很好地對狀態的持續時間進行建模，所造成的現象是心音轉移至下一狀態的概率不受當前狀態持續時間的影響，這導致了基于 HMM 的算法用于心音狀態分割不夠準確。
針對上述問題，考慮到基于 HMM 算法對于心音狀態建模的優勢，本文使用高斯混合模型對心動周期持續時間進行建模，構建基于持續時間的隱馬爾可夫模型（duration hidden Markov model，DHMM），以探求一種魯棒性與準確性更高的心音分割算法。
1 模型與方法
1.1 數據來源
本文研究的心音數據來源于研究團隊采集構建的心音樣本庫。該心音樣本庫志愿者年齡分布于 6 個月至 18 歲之間，所有志愿者或其監護人均簽署知情同意書。本研究經云南大學醫學院倫理委員會審查通過。樣本庫中所有心音信號所屬分類疾病，均通過臨床專家使用超聲心動儀（Acuson Cypress，德國西門子股份有限公司）確診。采集所使用的心音傳感器為美國 ThinkLabs 公司的 THE ONE，其輸出為模擬音頻信號，采用實驗室自行研發的采集設備，可同步記錄心音心電信號，采樣頻率為 5 000 Hz。該采集設備能將模擬信號轉換為數字信號，便于儲存和處理。
1.2 心音信號的具體分割流程
本文心音分割算法的具體流程如圖1 所示，整個心音分割流程可概括為以下 7 個步驟。
圖1
心音分割算法流程圖
Figure1.
Flow chart of heart sound segmentation algorithm
（1）使用 R 峰檢測器與 T 波探測器對同步采集心音心電的信號進行位置標注；
（2）將標記的信號分為訓練集與測試集；
（3）對訓練集信號與測試集信號提取心動周期等特征；
（4）對訓練集信號進行高斯混合模型建模，為 DHMM 提供參數 、 與 pj（d）；
（5）構建 DHMM，并使用訓練集信號對模型參數 λ 進行訓練；
（6）使用訓練后的 DHMM 對測試集信號進行建模，并通過改進的維特比算法尋找最優狀態序列；
（7）使用測試集信號心音心電位置的對應關系評估模型性能指標。
其中，步驟（1）和（2）中對測試集信號心音心電進行位置標注，其目的是為步驟（7）評估模型性能指標提供參考依據。在實際使用中，對于沒有任何參考標記的心音信號，經特征提取后能直接由訓練完成的 DHMM 進行心音分割。
1.3 信號標注
為了訓練基于 DHMM 的心音數據分割算法，必須使用外部參考信號標記 PCG 記錄中 S1 和 S2 的起始位置。使用 ECG 作為心音狀態 S1 和 S2 的參考標記位置，其中 ECG 的 R 峰對應 PCG 信號中 S1 的起始位置，而 S2 發生在 ECG 的 T 波末端處。因此首先需要在 ECG 中標記 R 峰與 T 波末端位置。本文采用了基于拋物線擬合的 R 峰檢測算法[15]與 Laguna 等[16]所提出的 T 波探測器，分別對 200 例訓練信號進行標記。訓練樣本中包括 80 例異常信號和 120 例正常信號，所有信號均來源于 1.1 節中的心音樣本庫。
由于 S1 的起始位置與 ECG 的 R 峰重合[17]，將一個 R 峰至下一個 R 峰之間 S1 平均持續時間段標記為 S1。S2 出現在 ECG 的 T 波末端附近[17]，且 S2 的振幅在 T 波末端附近將達到最大值。因此，設置搜索窗口在 T 波末端附近尋找 PCG 信號 Hilbert 包絡中的最大峰值，以此峰值作為 S2 中心，將其周圍 S2 平均持續時間段標記為 S2。收縮期與舒張期介于 S1 與 S2 的相對位置之間，將 S1 到 S2 之間的時間段標記為收縮期，S2 到下一個 S1 之間的時間段標記為舒張期。完成標記后的 PCG 如圖2 所示。
圖2
R 峰與 T 波末端分別作為定義 S1 與 S2 的參考位置
Figure2.
R peak and T wave terminal are the reference positions of S1 and S2 respectively
1.4 心動周期的提取
心音信號具有短時平穩性，研究中普遍認為在 10～30 ms 內心音信號可視作平穩[18]。采用自相關法[19]能較好地提取心動周期的長度，由于平滑的心音包絡更容易反映心動周期的特征，因此對心音信號的 Hilbert 包絡求自相關。將心音信號的 Hilbert 包絡記為，其中 N 表示信號的總長度，則自相關計算公式可以表示為
|
提取出的自相關包絡如圖3 所示，在圖中所示兩個紅點之間的間距為心動周期的持續時間（Dcycle），第一個紅點與綠點之間的間距稱為一個完整的收縮期（systole），其成分包括了 S1 與收縮期兩個狀態。其持續時間（Dsystole）對應公式如下
圖3
自相關法提取心動周期
Figure3.
Cardiac cycle was extracted by autocorrelation method
|
完整的舒張期（diastole）包括 S2 與舒張期兩個狀態，其持續時間（Ddiastole）有
|
由式（2）與式（3）可以得到心音各個狀態持續時間的約束關系，使得接下來對心動周期持續時間的高斯建模更加符合真實的心音持續時間分布。
1.5 心動周期持續時間的高斯建模
由于在平靜狀態下人類的心率大多介于 60～100 bpm 之間，則心動周期持續時長 heartcycle_duration 為
|
其中 heartrate 代表心率。由式（4）可推出一個心動周期時長為 600～1 000 ms。圖4 是對 200 例訓練集心音信號單一心動周期的各狀態持續時間的統計直方圖，其中 n1、n2、n3、n4 分別對應 S1、收縮期、S2、舒張期的心音狀態持續時間出現次數。由于 n1、n2、n3 的持續時間較為集中，在相同的時間段內出現頻次有重疊部分，故繪制了對應狀態的包絡以便觀察。從圖中可看出，各個狀態的持續時間近似于正態分布。
圖4
對應狀態持續時間的出現頻次
Figure4.
The occurrence frequency of state duration
由于一個心動周期有四個狀態，故對心動周期的持續時間進行高斯混合模型建模[20]。定義 X = {x1，x2，，xN}為觀測的隨機向量，又稱為不完全數據（incomplete-data），其中 x1，x2，，xN 指觀測到的心音持續時間。（X，Z） = {（x1，z1），（x2，z2），，（xN，zN）}稱為完全數據。Z 為隱隨機變量數據，表示對應的樣本 x 是屬于哪一個高斯分布。根據之前的分析，心音四個狀態的持續時間近似于正態分布。將心音各狀態的集合表示為 C（condition），有 C = （C1，C2，C3，C4），其中 C1 表示 S1，C2 表示收縮期，C3 表示 S2，C4 表示舒張期。可以確定 Z∈（C1，C2，C3，C4）。θ = {p1，p2，，pk μ1，μ2，，μk ∑1，∑2，，∑k}為所需估計的參數，由于各狀態高斯模型的均值 μ 與方差∑可以通過訓練集樣本統計得出，故所需確定參數為 θ = {p1，p2，p3，p4}，其中 p1、p2、p3、p4 表示持續時間為 x 時分別屬于 C1、C2、C3、C4 高斯分布的概率。由此，可以得出離散概率分布如表1 所示
表1
高斯混合模型離散概率分布
Table1.
Discrete probability distribution of Gaussian mixture model
高斯混合模型的概率密度函數為
|
式（5）中 N（x | μk，∑k）代表著 x 服從于第 Ck 類正態分布。式中 μk、∑k 通過訓練信號統計得出。使用最大期望算法（expectation-maximization algorithm，EM）計算模型最優參數 pk，EM 算法的遞推表達式為
|
將式（6）中的 稱為 Q 函數，有
|
通過式（6）和（7），可以將其整理為
|
在式（8）中，需要設定初始參數 θ（1），由于在采集過程中無法確定從哪個狀態開始采集，假設初始狀態等概率出現，則 。根據遞推式可得到任意 t 時刻 θ（t）的值，因此只需求解未知參數 θ。使用拉格朗日乘數法確定 θ，得到遞推參數 pk（t + 1）：
|
可由式（9）推出 p（t + 1） = {p1（t + 1），p2（t + 1），p3（t + 1），p4（t + 1）}，由此可以遞推各個時刻的最優參數。
1.6 心音信號的 DHMM 建模
DHMM 與傳統 HMM 最大的區別就是引入了持續時間的約束轉移條件。雖然傳統 HMM 很適合對心音信號狀態分割進行建模，但由于傳統 HMM 的限制性，所涉及的馬爾可夫過程是無記憶性的，處于下一狀態的概率僅取決于當前時刻的狀態，而與以前的狀態無關。這意味著心音不論在當前狀態持續多長時間，轉移到新狀態的概率都是相同的，導致每個狀態的持續時間呈幾何分布[21]，這顯然不符合基于時域的心音信號。雖然可以通過引入更多數量的子狀態來處理這一限制，但無形中引入了更多的參數，使模型變得復雜化，而且這只是緩解了癥狀，沒有從根本上解決問題。由此，引入持續時間密度 pj（d），其中 d（duration）表示狀態持續時間，概率密度函數 pj（d）表示在狀態 j 持續時間為 d 所發生的概率。j∈（C1，C2，C3，C4），故引入 pj（d）后，可以將模型簡化為四個狀態。
|
其中 P（d）由本文第 1.5 節的高斯混合模型確定。
記 t 時刻模型所處的狀態為 ct，則有 ct∈（C1，C2，C3，C4）。設觀測值序列集合為 O，O = {O1，O2，O3，，Ot}。定義 A 為狀態轉移概率矩陣，則
|
aij 表示從狀態 Ci 到狀態 Cj 轉移時的轉移概率[22]。由于心音是規律性振動產生的準周期性生物信號，且總是按照 S1、收縮期、S2、舒張期循環出現。則狀態轉移概率集合
|
DHMM 觀測序列與隱含狀態序列對應關系如圖5 所示，其中狀態間轉移概率由概率密度函數 pj（d）控制。
圖5
基于心動周期狀態的隱馬爾可夫模型
Figure5.
HMM based on cardiac cycle state
定義輸出觀測值概率集合 B = {bj（Ot）}，其中 bj（Ot）是在狀態 j 時觀測值符號 Ot 輸出的概率。
定義初始狀態概率集合為 π = {πi}，其中
|
由于實際采集心音過程中可從心動周期的任何一個狀態開始采集，因此
|
在使用 DHMM 之前，需要通過訓練集數據對其模型參數 λ 進行估計
|
其中參數 、 分別為觀測序列處于狀態 j 持續時間的平均值和方差，由 1.5 節中的高斯混合模型確定。
1.7 維特比算法尋找最優序列
心音分割的目的是找出所有已知觀測序列所對應的狀態序列，在已知模型參數 λ 的前提下，可以使用維特比算法（Viterbi algorithm）尋找心音信號的最佳狀態序列。
維特比算法是一種前向算法，用于計算前向概率 。表示時刻 t 時 ct 處于狀態 j 的瞬時概率。由于 DHMM 引入了持續時間的概念，故對傳統的維特比前向算法[21]進行改進使其適應 DHMM，改進后的前向概率 表示為：
|
給定持續時間 d，且前一狀態在時刻 時為 i，則 ct 在狀態 j 結束的概率為
|
其中 是上一個狀態 i 在 時刻結束的概率。是在狀態 j 時持續 到 t 期間觀測值輸出的概率。找到使式（16）最大化的參數 d 與前一狀態 i 來計算前向概率 [23]：
|
持續時間 d 限制在 內且小于 t，其中 和 表示持續時間為 d 時所對應高斯混合模型的均值與方差。對于每個 t，存儲使式（17）最大化的持續時間 d 和前一狀態 i：
|
在回溯之前，尋找最優狀態 j，使得結束狀態 處于 j 的概率最大：
|
實現回溯算法的具體程序流圖如圖6 所示。在圖6 中，首先初始化參數 ，將當前狀態設為觀測序列末尾所處的狀態。判定條件為 t > 1，目的為從序列末尾狀態回溯至序列首狀態。讀取出 t 時刻最優狀態 的持續時間 d*。 實現了將 時刻至 時刻之間的狀態標記為 。 將 時刻的狀態調整為 的前一個最優狀態。 是下一次程序迭代的前提。當程序結束時，得到最優狀態序列集合 C* = {c1*，c2*，c3*，，cT*}。
圖6
回溯算法程序流程圖
Figure6.
Program flow diagram of backtracking algorithm
2 實驗結果與分析
2.1 實驗環境與數據說明
實驗所用電腦 CPU 型號為 Intel（R）Core（TM）i7-8750H CPU @ 2.20 GHz，RAM 為 16.0 GB，仿真軟件為 MATLAB2019b。
本文算法與對比算法均使用本文測試集數據進行測試，測試集數據來源于 1.1 節所提及的心音數據庫。由于需要衡量各個算法在針對不同病理心音信號進行分割時的性能，故準備了 4 個獨立的測試集共 500 例心音信號對其進行評估，分別為 290 例正常（Normal）信號、70 例房間隔缺損（atrial septal defect，ASD）信號、70 例動脈導管未閉（patent ductus arteriosus，PDA）信號與 70 例室間隔缺損（ventricular septal defect，VSD）信號。
2.2 評價指標
為了評估模型對于分割心音信號的能力，使用 PCG-ECG 對照法評估模型在測試集中準確定位 S1 與 S2 的能力。S1 與 S2 的參考位置在 1.3 節中已經提及，所有參考位置均使用 R 峰檢測器與 T 波探測器進行標記。由于 R 峰檢測器與 T 波探測器標記位置時存在 100 ms 的檢測誤差。故若分割 S1 的起始位置在 ECG R 峰的 100 ms 內，則 S1 被標記為正確分段（TP）。同樣地，若 S2 中心位置在對應 T 波末端的 100 ms 內，則 S2 聲音被標記為正確分段（TP）。其他所有未被正確分段的位置標記為錯誤分段（FP）。圖7 為心音心電分割狀態對照圖，其中狀態序列（圖中紅色虛線所示）有 4 個狀態，狀態 1、2、3、4 分別對應 S1、收縮期、S2、舒張期。該圖可從心音、心電的位置對應關系客觀地評估模型分割性能指標。
圖7
心音-心電圖評估模型分割結果
Figure7.
Segmentation result of PCG-ECG evaluation model
對于心音分割這一任務，由于并不存在真正的負樣本，所以將靈敏度（sensitivity）定義為：
|
其中 TP 表示分割正確的數量，All_signal 表示所有 S1、S2 位置的總數。recall 表示召回率，在數值上與 sensitivity 相等。
精確率（precision）定義為：
|
其中 FP 表示分割錯誤的數量。
使用評估精度分數 F1[24]綜合評價 sensitivity 與 precision 指標：
|
2.3 實驗結果
本文選取引言中提及的基于 Hilbert 的算法[11]與基于傳統 HMM 算法[12]作為對比算法。這兩類算法都有一定的抗噪聲性能，且運行效率較高，對含有部分非病理性噪音干擾的心音信號進行分割時均有較好表現。表2～5 反映了各算法在處理不同類型心音信號時的性能表現。
表2
對于正常心音信號各算法的評價指標（）
Table2.
Evaluation index for each algorithm of normal heart sound signal（）
表3
對于房間隔缺損心音信號各算法的評價指標（）
Table3.
Evaluation index for each algorithm of ASD heart sound signal（）
表4
對于動脈導管未閉心音信號各算法的評價指標（）
Table4.
Evaluation index for each algorithm of PDA heart sound signal（）
表5
對于室間隔缺損心音信號各算法的評價指標（）
Table5.
Evaluation index for each algorithm of VSD heart sound signal（）
圖8 為使用本文算法對各類心音進行分割的結果，該結果由程序自動生成。從圖中可觀察出對于病理性的心音信號，DHMM 仍能高精度地對其進行心音狀態分割。
圖8
對各類心音信號的分割結果
Figure8.
Segmentation result of various heart sound signals
2.4 討論
本文使用 4 個獨立的數據集對算法性能進行評估，評估指標結果均由程序自動運算得出，可以最大限度地降低主觀因素的影響。對比以往文獻所使用的測試數據與評估方法，本文實驗結果可以更為客觀地反映出各算法在處理不同類型心音信號時的性能表現。
從 2.3 節的表中可看出，對正常心音信號進行分割時，三種算法都能實現高精度的心音分割。其中文獻[12]基于傳統 HMM 心音分割算法使用同態包絡提取特征序列作為 HMM 的觀測序列，然后通過無監督學習對心音狀態進行分割。通過觀察表中數據可以發現，該算法總體標準差要大于其他兩種算法。出現這種現象的原因是此算法引入了大量子狀態用于緩解 HMM 無記憶性的特點，心音分割性能很大程度受限于轉移概率矩陣參數，這使得該算法分割性能在處理不同樣本時會產生較大變化。值得注意的是，該算法準確率同引言中描述相差較大，推測是由于文獻中僅使用 44 例數據進行測試，且分割精度評價指標不如本文嚴格。文獻[11]基于 Hilbert 包絡的心音分割算法運行效率最高，在對正常心音信號進行分割時性能表現穩定，但此算法在處理病理性心音信號時，分割精度明顯下滑。這是由形態學進行心音分割時的局限性所導致的，因為程序很難實現對復雜性病理心音信號進行包絡邊緣檢測。本文提出的基于 DHMM 的算法針對傳統 HMM 算法的缺點進行改進，對心音信號持續時間進行高斯混合模型建模，并將模型狀態數精簡為 4 個。從表中可看出，本文算法在對不同類型心音樣本進行處理時性能表現均優于其他兩種算法。與基于傳統 HMM 的算法相比，本文算法運行效率、泛化能力以及分割精度都有了大幅提升。基于心音狀態循環周期反復出現的特點，本文構建 DHMM，使得算法具有較高的魯棒性和抗噪聲性能，對復雜性病理心音信號與正常心音信號均能適用。
3 結論
基于 DHMM 的心音分割方法可以很好地對心音信號的準周期性、短時平穩性以及時域心音狀態分布進行建模。在使用訓練集信號完成 DHMM 訓練后，該算法能對無任何位置參考的心音信號進行心音狀態分割。由于算法具有較高的魯棒性與抗噪聲性能，更加符合臨床環境的使用需求。本算法在整體分割性能上相比傳統 HMM 分割算法有大幅提升，但受限于 DHMM 在序列開始和結束時對觀測序列建模能力的限制，對起始狀態與終止狀態的分割精度還存在少許不足。后續可對該問題進行研究，進一步提高算法精度，擴大算法適用范圍。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
