睡眠呼吸暫停綜合征患者的睡眠分期算法研究

魯柯柯 祁 霞 張建保 王 剛＃* 閆相國＃*

1(生物醫學信息工程教育部重點實驗室，西安交通大學生命科學與技術學院， 健康與康復科學研究所，西安 710049)2(國家醫療保健器具工程技術研究中心， 廣州 510500)3(神經功能信息學與康復工程民政部重點實驗室，西安 710049)

引言

睡眠呼吸暫停綜合征(sleep apnea syndrome，SAS)是發生在睡眠過程中呼吸出現短暫停止現象的一種睡眠障礙性疾病[1]，與高血壓、冠心病以及心律失常等心血管疾病密切相關[2]。 臨床上通過睡眠分期來分析SAS 患者的睡眠結構和睡眠質量，為診斷睡眠疾病提供有效信息。

臨床按照美國睡眠醫學學會(The American Academy of Sleep Medicine，AASM)標準和多導睡眠圖(polysomnography，PSG)，以30 s 為間隔將夜間睡眠活動劃分為5 個階段:清醒(wake，W)期、非快速眼動(non-rapid eye movement，NREM)Ⅰ(N1)期、Ⅱ(N2) 期、 Ⅲ(N3) 期和快速眼動(rapid eye movement，REM)期，其中N1 期、N2 期統稱為淺睡(light sleep，LS)期，N1 期、N2 期和N3 期統稱為NREM 期[3]。 SAS 患者相比于健康人N3 期和REM期在整晚睡眠中占比較低，且睡眠呼吸暫停主要發生在N3 期[4]。 因此針對SAS 患者通過睡眠分期分析其睡眠結構具有重要意義。

PSG 作為臨床診斷睡眠疾病和睡眠分期的金標準，記錄了夜間睡眠腦電圖(electroencephalogram， EEG)、 眼電圖( electro-Oculogram，EOG)、心電圖(electrocardiography，ECG)等多種生理信號[5]，然而因其價格昂貴及受試者需在專業的睡眠實驗室佩戴大量傳感器，很難適用于普通的家庭以及病房監護。 因此，提出一種能夠利用便于檢測的生理信號進行自動睡眠分期的方法具有實際應用價值。

心率和呼吸信號與睡眠階段高度相關[6-7]，并且測量技術較為成熟。 Wang 等[8]從7 名SAS 患者的ECG 信號中以每5 min 為間隔提取時域、頻域和非線性特征共24 個，使用序列前向選擇(sequential forward seelction，SFS)進行特征選擇，最后使用基于決策樹的支持向量機(decision-tree-based support vector machines，DTB-SVM) 分類器實現W-REMNREM 三分類的睡眠分期任務，采用留一交叉驗證法測試得到73.51%的平均準確率。 Willemen 等[9]提取了25 例SAS 患者心電、呼吸信號的時域及頻域特征，使用線性判別分析(linear discriminant analysis，LDA)在W-REM-NREM 三分類睡眠分期任務上獲得了70%的準確率，Cohen's Kappa 系數為0.41。 黃文漢等[10]通過提取23 例SAS 患者心電和呼吸信號的時域特征、頻域特征和耦合特征，使用隨機森林(random forest，RF)在W-REM-NREM 三分類睡眠分期任務上獲得了71.9%的準確率(Cohen's Kappa ＝0.36)。

上述基于心電和呼吸信號的自動睡眠分期研究取得了較好的分期效果，但都是以三分類為主，應用范圍和場景具有較大局限性。 而且基于小規模數據集訓練測試的方法或模型往往缺乏一定的泛化能力，當用于其它數據集時，準確率可能只會下降10%左右，而Cohen's Kappa 系數卻會下降20%左右[11]。 同時，睡眠是一個結構化過程，不同的睡眠階段之間彼此緊密聯系，但是已有的針對SAS 患者進行睡眠分期的研究忽略了時間維度信息。 另外值得注意的是，SAS 嚴重影響心電信號的穩定性，睡眠分期數據較普通人而言更加不連續和不平衡[12-13]，現有的基于健康人的睡眠分期方法[14]用于SAS 患者時分類性能嚴重下降。 針對上述問題，本研究以R-R 間期序列為基礎，利用較大規模數據集，提出一種基于門控循環單元(gate recurrent unit，GRU)網絡進行多分類自動睡眠分期的方法。

1 材料和方法

1.1 SHRSV 數據集

所采用方法的總體描述如圖1 所示。 首先從R-R 間期序列中提取心率變異性(heart rate variability，HRV)、 呼吸幅度變異性(respiratory amplitude variability， RAV) 和呼吸率變異性(respiratory rate variability，RRV)信號；然后從這些信號中提取時域、頻域和非線性特征共55 個，并且根據分類器類型，對原六分類的睡眠狀態注釋進行合并處理，分別得到4 種對應不同分類器的睡眠標簽；最后基于GRU 網絡分別構建區分W-S 的二分類器、W-REM-NREM 的三分類器、W-REM-LS-SWS的四分類器、W-REM-N1-N2-N3 的五分類器。

圖1 自動睡眠分期整體框圖Fig.1 An overall description of the method proposed

所采用的訓練集和測試集來自SHRSV(sleep heart rate and stroke volume，SHRSV)數據庫[15]，該數據庫包含274 例SAS 患者[男229 例，女45 例，平均年齡(61±10)歲]，其中142 例為輕度SAS 患者(呼吸暫停低通氣指數(apnea hyponea index，AHI)＜15)，132 例為中重度患者(AHI＞15)。 數據庫僅由R-R 間期(heart interbeat interval)序列、每搏輸出量(stroke volume，SV)數據和睡眠狀態注釋組成。 R-R間期序列是從采樣率為500 Hz 的心電信號中采用自動心率分析軟件提取，并經過手動檢視和校正；睡眠狀態注釋是由專家按照早期的Rechtschaffen 和Kales(R＆K)睡眠分期標準將睡眠活動分為W、REM、N1、N2、N3、N4 等共6 個階段[16]。

1.2 特征提取

提取的特征主要包括R-R 間期序列、HRV、RAV 和RRV 的時域特征，HRV 和RRV 的頻域及非線性特征。 首先從R-R 間期序列得到心率信號HR，然后利用三次樣條插值對HR 進行插值處理，得到采樣率10 Hz 的等時間間隔HRV 信號。

根據心率變化受呼吸過程調制這一特點，使用一個截止頻率為0.1～0.4 Hz 的三階巴特沃茲帶通濾波器處理HRV 信號，即可得到頻率調制RAV 信號。 RRV 信號提取方法與HRV 相似，首先對RAV信號使用最大斜率法[17]檢測呼吸波位置，并計算相鄰呼吸波峰位置的時間間隔，即可得到P-P 呼吸間隔序列；然后從P-P 呼吸間期序列計算呼吸率信號，對其利用三次樣條插值進行插值處理，得到采樣率10 Hz 的等時間間隔RRV 信號。

提取特征之前，首先對HRV、RAV 和RRV 進行均值為0，方差為1 的標準化處理，并對整晚數據按30 s 為單位進行劃分。 考慮到計算頻域特征需要一定的長度，尤其是獲取信號超低頻帶(very low frequence，VLF)的可靠信息至少需要5 min[18]，同時為了保證特征數據的連續性和平穩性，時域、頻域和非線性的特征計算均以當前30 s 為中心，前后各延伸2.25 min，共5 min 為一段，最后將5 min 數據計算的特征作為中心30 s 信號的特征。

從R-R 間期序列、HRV、RAV 和RRV 信號中提取的17 個時域特征如表1 所示，對于R-R 間期序列，提取了特征RRn 和PNN50，RRn 表示該段R-R間期序列中最常出現的R-R 間期數值，PNN50 表示相鄰R-R 間期之差超過50 ms 的個數與總竇性心搏個數的比值，其余的15 個時域特征在文獻[14]中有詳細說明。 此外，HRV 和RRV 的頻域特征能夠深層次的反映心率和呼吸率的變化規律。 表2 是這一部分的特征說明，其中特征SA、SB、SC、SA/C、SB/C，分別表示0.05 ～0.40 Hz、0.005 ～0.010 Hz、0.01 ～0.05 Hz 頻段的能量、0.05～0.40 Hz 與0.01～0.05 Hz能量的比值、0.005 ～0.010 Hz 與0.01 ～0.05 Hz 能量的比值，其余HRV 和RRV 分別對應的5 個頻域特征來自文獻[14]使用的特征。

表1 時域特征提取Tab.1 The time domain features from R-R intervals sequence， HRV， RAV and RRV

表2 頻域特征提取Tab.2 The frequency domain features from HRV and RRV

考慮到心臟是一種復雜的非線性混沌系統，時域、頻域特征很難完整揭示反映其本質特征的關鍵信息[19]，進一步提取了原信號S，即HRV 和RRV 信號的多尺度熵(multiscale entropy，MSE)[20]、過零檢測(zero crossing analysis，ZCL)、赫斯特指數和多重分形譜寬度[21]4 種非線性特征。 具體特征如表3所示，其中Sl和Sh為小波分解和重構后的低頻和高頻信號，分別視為HRV 和RRV 信號LF 和HF 分量的估計。

表3 非線性特征提取Tab.3 The nonlinear features from HRV and RRV

1.3 分類網絡

考慮到不同應用場景的需求，研究基于Tensorflow 2.1.0 軟件庫的Keras 2.3.1 框架構建了4 個分類器，其中Python 版本為3.7。 每個分類器共有7 層，分別是輸入層、Masking 層、全連接層、批標準化層、雙向GRU 層、Dropout 層和輸出層，最大區別僅在于輸出層的神經元數量不同，代表不同的類別數。

此外，每個人每晚睡眠時間長度的不同導致提取的特征長度也不一致，而keras 構建網絡模型時，輸入層需要指定特征序列的長度和個數。 研究選擇所有數據中最長的數據長度，即睡眠時間長度為9 h 2 min 30 s，以30 s 為一個時間單位，共劃分為1 085個時間單位，對于數據長度不足1 085 的數據，采用keras 的Masking 策略，即用0 填充。 圖2 為分類器的結構。

圖2 分類網絡結構Fig.2 Classification network structure

1)輸入層:輸入特征數據，指明輸入維度，即時間維度為1085，特征維度為55。

2)Masking 層:屏蔽輸入特征數據中補0 的時間單位，保證該時間單位在后續的所有層被跳過，不參與計算。

3)全連接層:將每個時間單位的特征線性映射到另一個特征空間，該層神經元數量為輸入特征維度的2 倍。 其中，每個神經元節點均與上一層每個神經元節點相連接。

4)批標準化層:具有正則化的作用，使得每一批輸入特征數據的分布重新置為均值接近為0，方差接近為1，這一層的作用能夠減少中間協變量的轉換，提高網絡訓練的速度[14]。

5)雙向GRU 層:設置100 個GRU 網絡單元，當前時間單位的輸出不僅依賴于當前時間單位的特征，還依賴于之前以及之后時間單位的特征。 設置序列標志位為true，意味著該層的輸出與輸入的時間維度一致。

6)Dropout 層:采用Dropout 策略[22]，隨機以0.5 的失活概率丟棄雙向GRU 層的輸出節點，這樣可以增強模型的泛化能力，防止網絡陷入過擬合。

7)輸出層:本質為一個全連接層，將前一層的輸出映射到輸出特征空間中，該層神經元個數為分類器的輸出類別數。 使用softmax 激活函數使得神經元輸出為概率值，各神經元對應的概率值之和為1，概率值最大的即為預測的類別。

1.4 模型訓練

使用交叉熵代價函數作為損失函數，用來衡量訓練過程中預測結果與實際結果之間的誤差。通過賦予損失函數類別權重系數來降低類別樣本數不均衡對模型訓練的影響，類別權重系數定義為

式中，wi為第i個類別的權重系數，i取值范圍為1～n，n為類別總數；N為樣本總數；Ni為第i個類別樣本數；α為權重指數，取值范圍為0～1。

從SHRSV 數據集的274 例數據中隨機選取20%作為測試集，剩余數據按照4 ∶1 隨機劃分訓練集和驗證集，即訓練集176 例，驗證集44例，測試集54 例，同時保證各數據集中AHI＜15以及AHI＞15 的數據在總數據集中所占比例基本相同。 訓練過程中的batch_size 為5，訓練次數epoch 為50，初始學習率為0.001，選擇收斂速度快且效果較優的Adam 優化方法[23]用來更新權重。

實驗環境為一臺處理器為Intel(R) Core(TM)i5-6400 CPU ＠ 2.71 GHz、RAM 為16 GB、Win10 64 位操作系統、裝有Matlab R2019b，Tensorflow-CPU 2.1.0 的DELL 臺式機。

1.5 評價指標

1.5.1 統計學評價

使用混淆矩陣可視化分類結果，并計算總體準確率和Cohen's Kappa 系數k，其中k的取值范圍在-1～1 之間，具體統計學意義為:k＜0 表示分類器預測結果與真實結果沒有一致性；0.01≤k≤0.2 表示有極低的一致性；0.21≤k≤0.4 表示有一般一致性；0.41≤k≤0.6 表示有中等一致性；0.61≤k≤0.8 表示有較高一致性；0.81≤k≤1 表示幾乎完全一致[24]。

1.5.2 睡眠結構評價指標

睡眠分期的主要目的是分析評估睡眠結構，進而判斷人體的睡眠質量。 共使用了11 個評價指標用以分析睡眠結構，對睡眠質量進行評價[25]，如表4 所示。

表4 睡眠結構評價指標Tab.4 Sleep structure evaluation index

2 結果

一晚睡眠數據的特征提取時間與30 s 時間單位特征的提取時間呈線性關系，其中30 s 時間單位的特征提取時間為0.81 s；模型的訓練時間隨類別數的增加而增加，二分類、三分類、四分類、五分類網絡的訓練時間分別約為32 min 10 s、42 min 31 s、47 min 32 s、53 min 23 s，測試時間隨類別數的變化微乎其微，每個類別測試所需時間均約1 s。

2.1 權重指數α 的影響

圖3 顯示了當損失函數被賦予類別權重系數wi時，三分類驗證集的準確率和Cohen's Kappa 系數隨權重系數wi中的權重指數α的變化過程。 可以看出當α取0.3 時Cohen's Kappa 系數最大且準確率無明顯降低，此時整體效果最優。

圖3 準確率和Cohen's Kappa 系數隨α 的變化Fig.3 The change of accuracy and Cohens Kappa statistic k with different α

2.2 模型在測試集上的表現

表5～8 分別給出了測試集二分類、三分類、四分類和五分類的混淆矩陣、準確率和Cohen's Kappa系數k，其中4 個分類器的準確率分別為85.06%、75.44%、63.80%、62.13%，Cohen's Kappa 系數達到了0.54、0.49、0.41、0.41，均達到了與臨床PSG 結果中等的一致性。

表5 二分類睡眠分期的混淆矩陣Tab.5 Confusion matrix of two-class sleep staging

表6 三分類睡眠分期的混淆矩陣Tab.6 Confusion matrix of three-class sleep staging

表7 四分類睡眠分期的混淆矩陣Tab.7 Confusion matrix of four-class sleep staging

表8 五分類睡眠分期的混淆矩陣Tab.8 Confusion matrix of five-class sleep staging

圖4 是測試集中某一受試者臨床分期結果與文中方法預測結果的對比，可以看出只有少部分時間單位睡眠階段被分錯。 該受試者預測結果與專業醫師的分析結果相比，準確率達到了86.39%，Cohen's Kappa 系數為0.80，與臨床PSG 判讀結果有較高的一致性，而臨床中不同醫師判讀結果的一致性僅有76.8%[26]。

圖4 模型分期結果與臨床醫師分期結果對比.模型分期結果(上)；專業醫師分期結果(下)Fig.4 The comparison of predicted sleep stage by GRU network with clinical technician of one subject from the test data.Sleep stages by GRU network (the top).Sleep stages by clinical analysis (the bottom).

用表4 列出的11 個評價指標，對54 個測試用例進行統計分析。 使用SPSS 軟件分析發現相關評價指標不滿足正態分布和方差齊性，因此采用曼-惠特尼U 檢驗規則檢驗[27]預測結果和臨床PSG 判讀結果之間的統計學差異。 表9 分別列出了預測結果和臨床PSG 判讀結果的均值、標準差和顯著性水平P值，由于P值均大于0.05，表明預測結果和臨床PSG 判讀結果之間的差異無統計學意義。

表9 睡眠結構評估指標的統計學分析Tab.9 Statistical analysis of sleep structure evaluation indicators

2.3 特征選擇的可行性

為了從提取的特征集中篩選對當前分類任務有益的和剔除冗余的特征，提出一種更加簡化的基于深度學習模型的封裝式方法，具體過程為:對于原始特征集合F＝{f1，f2，…，fm}，首先將每個特征fi獨立的依次輸入分類器，并依據驗證集上的Cohen's Kappa 系數由大到小對F進行排序得到新的特征集合F' ＝{f1'，f2'，f3'，…，f n'} 。 然后將最優的特征f1'放入空特征子集Fs，接著將F'中剩下的n -1 個特征依次加入Fs， 若加入特征f i'后得到的Cohen's Kappa 系數比之前好，則保留特征f i'， 否則去除。最后特征子集Fs作為最終輸入分類器訓練和測試的特征集合。

表10 列出了不同粒度分類器篩選的特征集合，54 例測試數據的二分類、三分類、四分類和五分類的準確率分別是 85.10%、 76.16%、 69.07%、64.44%，Cohen's Kappa 系數是0.54、0.51、0.45、0.42。 可以看出三、四和五分類的分類效果相比于特征選擇之前得到了提升，同時減少了需要的特征規模，節約了計算時間，充分說明了特征選擇的有效性。

表10 選擇的特征集合Tab.10 Selected set of features

2.4 國內外同類研究對比

表11 列出了課題組前期Wei 等[14]提出的基于健康人的睡眠自動分期算法應用于SAS 患者數據后的表現，可以看出不同粒度睡眠分期任務的效果大大降低，尤其是五分類準確率降低了14.66%，Cohen's Kappa 系數降低了0.22。 而表12 是本文模型或算法與國內外同類研究的對比結果，其中由于公開數據庫UCD[28]中數據較少，直接采用文中訓練好的模型對25 個SAS 患者進行預測，而對于健康人數據的處理過程與上文描述一致，由此發現準確率和Cohen's Kappa 系數均有所提高，這也表明所提出的自動睡眠分期算法的可行性與有效性。

表11 健康人的自動睡眠分期算法在不同數據集中的表現Tab.11 Performance of automatic sleep staging algorithm in healthy subjects in different data sets

表12 國內外同類研究方法的對比Tab.12 Comparison of similar research methods at home and abroad

3 討論

本研究從SAS 患者的R-R 間期序列中提取了HRV、RAV 和RRV 等3 種信號，進一步地從這些信號中提取了時域、頻域及非線性特征共55 個，基于GRU 網絡實現了不同分類粒度的睡眠自動分期。其中所有睡眠分期任務的Cohen's Kappa 系數均在0.41～0.61 范圍內，意味著與臨床PSG 結果有中等一致性，可是與課題組前期關于健康人的睡眠分期研究[14]相比，所有睡眠分期任務的表現較差，一方面是由于SAS 患者與健康人相比，睡眠分期類別數呈現更加不平衡的特點，本研究雖通過賦予損失函數類別權重系數來緩解睡眠分期類別樣本數嚴重不均衡對分類器學習過程的影響，但是從表5 ～8 中依然能看出二分類中有較多的W 期被錯分為S，三分類中有較多的W 期和REM 期被分為NREM 期，四分類中錯分為LS 期的數據較多，五分類中錯分為N2 期數據較多，這也表明采用的樣本權重雖然在一定程度上能夠緩解，但并不能完全解決樣本不均衡的問題。 另外，睡眠過程中的呼吸暫停事件也會模糊W 期和S 期、REM 和NREM 之間的決定界限[9]，這也是上述分期任務中被錯分的主要原因。

睡眠呼吸暫停是睡眠障礙中的一類普遍問題，相應的患者對在家庭或社區可進行睡眠監測的需求也日益強烈[10]，而心電信號由于其測量技術成熟，已經成為可穿戴睡眠監測設備常用的生理信號，主要用來對受試者進行睡眠分期和睡眠結構評價。 可是表11 中列舉的一些針對SAS 患者的研究，準確率和臨床PSG 金標準結果相比還是有很大的差距，這也一直是制約可穿戴睡眠監測設備的推廣使用的因素。 此外，這些研究只關注了準確率和Cohen's Kappa 系數兩個評價指標，而可穿戴睡眠監測設備除了對睡眠進行分期外，還包括睡眠質量的評估，它們并沒有從表4 列舉的睡眠結構指標角度去綜合評估算法的性能。 而本研究計算了11 個睡眠結構評估指標并利用統計學方法進行分析結果如表9 所示，所有指標均與臨床結果之間的差異無統計學意義，由此發現本研究使用GRU 網絡對SAS患者進行睡眠自動分期的有效性，雖然準確率與臨床金標準PSG 的結果還有很大差距，但總的來說對于睡眠結構分析和睡眠質量評估還是比較準確的，一定程度上可移植于可穿戴睡眠監測設備，用于用戶睡眠結構和睡眠質量的分析評估。

國內外同類研究往往集中于健康人，針對SAS患者的研究甚少，可是從表11 看出當健康人的算法用于SAS 患者時效果明顯變差，這表明算法的泛化或推廣能力較低，主要原因是SAS 嚴重影響心電信號的穩定性且Wei 等[14]的算法沒有考慮睡眠分期類別樣本數不均衡的影響。 然而從表12 看出本研究提出的睡眠自動分期算法和已訓練好的睡眠分期模型在不同人群中有良好的普適性，值得注意的是，Willemen 和黃文漢等使用的訓練集和測試集均來自于同一個數據集，一定程度會導致結果相較于訓練集和測試集彼此獨立的結果大幅升高，而SHRSV 和UCD 數據集來自于不同的采集設備與醫師標注，會引起數據分布存在很大差異。

此外，雖然本研究在SAS 患者的數據基礎上提出了一種自動睡眠分期算法，但是當用于健康人數據集時，表11 和表12 的結果對比表明，所有睡眠分期任務的表現均得到提升，尤其是二、三和四分類的Cohen's Kappa 系數大于0.61，與臨床PSG 結果達到較高的一致性。 這可能是因為本研究相比于課題組前期Wei 等[14]提取了更多的特征，尤其是HRV 和RRV 信號頻域的SB/C特征量即0.005～0.01 Hz 的能量與0.01 ～0.05 Hz 的能量的比值，它可以較好區分睡眠呼吸暫停和正常片段。 另一方面，考慮到人體心臟系統的非線性特性，本研究在時域、頻域特征的基礎上提取了一些非線性特征，這些特征的重要性在表10 中也有所體現。 綜上，本研究基于SAS 患者的數據提出的自動睡眠分期算法具有較好的穩健性和普適性，用于可穿戴睡眠監測設備時可以適用于不同人群，反之部分設備會首先進行睡眠呼吸暫停事件的判斷，然后加載與之對應的睡眠分期算法或模型完成睡眠分期和睡眠質量的評價，這無疑是增加了設備的分析時間和成本。

當然，本研究也存在一定的局限性。 一方面，提取特征過程復雜和特征數量較多，一定程度會限制它的推廣能力，除可以通過本研究提出的一種基于深度學習模型的特征選擇方法來降低特征數量外，未來將研究更有效的特征和復雜度更低的特征提取方法；另一方面，睡眠結構指標與臨床結果之間的差異無統計學意義，基本上滿足可穿戴睡眠監測設備中睡眠質量評估的需求，但是當考慮用于臨床某些場合，準確率仍亟待提高，未來也將圍繞改進模型來提高分類效果。

4 結論

睡眠分期是睡眠監測和睡眠質量評價中的一項重要內容。 本研究針對睡眠呼吸暫停綜合征患者提出了一種基于GRU 網絡的自動睡眠分期方法。在睡眠結構評估指標方面，本方法與臨床分析結果無統計學差異，基本可以滿足睡眠質量評估的需求，適用于可穿戴睡眠監測和家庭護理監護。