睡眠狀況是評價人體健康狀態的重要指標。本文提出一種基于枕下式的無擾睡眠監測系統，通過無擾獲取的心率信號測算心率變異性（HRV），并結合隱馬爾可夫模型（HMM），在對用戶無擾無接觸的環境下求解睡眠分期。針對現有 HMM 睡眠分期存在的問題，提出采用集合經驗模態分解（EEMD）消除 HRV 個體差異導致的分期誤差，再求解相應的睡眠分期。試驗選取廣州醫學院呼吸疾病研究所 10 例不同年齡及性別的無睡眠障礙的院內正常受試者，并與多導睡眠圖（PSG）睡眠分期結果相比較。研究結果證明本文所提無擾式睡眠監測方案可實現 S1～S4 睡眠分期，正確率超過 60%，且性能優于現有 HMM 睡眠分期方案。
引用本文： 李翔, 劉勇, 陳澎彬, 吳潔偉, 張涵. 基于隱馬爾可夫模型的枕下無擾式新型睡眠監測方案. 生物醫學工程學雜志, 2018, 35(2): 280-289. doi: 10.7507/1001-5515.201703059 復制
引言
睡眠分期結果可以對睡眠質量進行評估，是相關疾病診斷的重要指標[1-4]。20 世紀 60 年代美國加州大學腦研究所通過對腦電（electroencephalogram，EEG）、2 路頜肌電（electromyogram，EMG）以及 2 路眼電（electrooculogram，EOG）信號的分析研究，將睡眠狀態分為 6 期，即清醒期（wake，W）、快速眼動期（rapid eye movement，REM）、非快速眼動期（non-rapid eye movement，NREM），其中 NREM 期包括睡眠 Ⅰ 期（sleep stage 1，S1）、睡眠 Ⅱ 期（sleep stage 2，S2）、睡眠 Ⅲ 期（sleep stage 3，S3）、睡眠 Ⅳ 期（sleep stage 4，S4）。目前，EEG 已成為判斷睡眠分期的金標準，但是在此類監測方式下，受試者需長時間接觸多路電極，其自然睡眠會受到影響，不適于家居環境下的睡眠監護；而且研究所需的龐大數據量會耗費巨大的人力和時間成本。因此，探索其他生理信號在睡眠過程中的變化規律，并據此規律研究更加準確、高效的睡眠分期方法，是當前的研究熱點[5-8]。
心率的變化與自主神經系統的調控息息相關，自主神經系統的活動在不同睡眠期的表現也各不相同，現有研究[9-11]表明心率變異性（heart rate variability，HRV）中隱含著潛在的睡眠信息，具體表現為在上述 6 個分期階段心率變化呈現出不同的特點，通過合適的信號處理方法并結合先驗知識就可以反過來推算睡眠結構。Kortelainen 等[12]利用頻譜分析的方法對 RR 間期（RR interval，RRI）信號進行小波變換，將 RRI 通過頻譜變換分解成超低頻（very low frequency，VLF）、低頻（low frequency，LF）和高頻（high frequency，HF）三個頻譜區，經過對頻譜區特征的對比發現，LF/HF、VLF、VLF/HF 的均值與睡眠期存在線性關系，而且如果降低低頻信號的功率譜，能大大提高 S3 和 S4 期睡眠的識別率。莊志等[13]利用隱馬爾可夫模型（hidden Markov model，HMM）對 HRV 進行分析，在信號處理的過程中采用標準化的方式消除個體差異，對于睡眠分期的識別準確率可以接受。王金海等[14]利用支持向量機（support vector machin，SVM）實現了基于 HRV 信息的睡眠自動分期方法，采用主成分分析方法（principal component analysis，PCA）來降低數據冗余度，用少數幾個不相關的主成分替代原來相關互聯的眾多指標，獲得了較高的預測準確度。針對 HRV 的獲取方式，Wisse 等[15]利用壓力傳感的鉑電極傳感器獲得 HRV 信號，再結合體動信息通過 HMM 來預估睡眠分期。Ebrahimi 等[16]以睡眠心臟健康研究數據庫 HRV 信號為依據，提取特征值信息，再計算不同頻帶的能量和熵等指標對睡眠分期進行預估，能夠區分覺醒期和深睡期。
然而，目前睡眠分期方案存在以下 3 方面的問題：① 針對 HRV 信號的獲取，均需要用戶與多路電極直接接觸，用戶體驗度低[17-20]；② 國際公認的睡眠分期將睡眠狀態分為 6 期，現有睡眠分期方案僅分析對比睡眠/覺醒期，沒有針對 S1～S4 不同間期的睡眠分期對比驗證；③ 現有分期方案中，未考慮因用戶個體化差異帶來的分期誤差，分期輸出結果可信度相對較低。
針對上述問題，本文率先設計了一種枕下無擾式睡眠監控方案，主要貢獻可歸納為兩方面：
（1）針對 HRV 的獲取，提出采用壓電傳感器置于枕下的方式，無擾地采集用戶心率及 HRV 信息。所設計的采集平臺相較傳統院內、院外智能心率監測設備具有如下優勢：① 采樣率高達 1 000 Hz，達到醫療精準級別（目前院內心電采樣頻率 300～500 Hz），而市面可穿戴設備受功率制約采樣率極低，無法用于 HRV 的測量；② 新型非接觸無擾式檢測，不需與檢測對象直接接觸，相較現有院內專業采集設備需用戶與多路電極直接接觸的約束，具有明顯的用戶體驗度優勢。
（2）利用無擾心率檢測平臺進行睡眠分期的檢測，提出基于 HMM 的睡眠分期方案。與現有方案不同，為了消除因個體化差異導致的訓練輸出誤差，提出利用集合經驗模態分解（ensemble empirical mode decomposition，EEMD）去除基線漂移等歸一化預處理方案，消除 HRV 的個體化差異，同時，對歸一化后的信號進行限幅量化編碼，減少數據冗余度。所提優化方案將會在用戶體驗度和睡眠分期結果的準確性方面有顯著提升。
1 枕下無擾式實驗平臺
HRV 信號的采集通過枕下無擾式實驗平臺來實現，系統整體設計由枕下壓電傳感、信號采集卡以及終端處理三部分組成，設計簡圖如圖 1 所示。傳感器置于一塊平板中，放于枕頭底下，集成了印制電路板（printed circuit board，PCB）的采集卡放置于床邊，采集信號經過處理后在終端顯示。各部分的具體實現及功能如下：
（1）枕下壓電傳感內置采集模塊（系統的“感受裝置”），主要通過壓電傳感器來實現，通過比較后我們選擇的傳感器是壓電陶瓷蜂鳴片，它具有厚度薄、體積小、穩定性好、靈敏度高、耗能低、耐高低溫、壽命長、耐高壓、無需供電等優點，可以使信號采集模塊在枕下發揮功效。當壓電陶瓷蜂鳴片被施加直流電壓時，產生交替重復的機械彎曲變形，從而在空氣中產生聲波；然而當有外力施加時，壓電陶瓷片會受外力影響產生電波，將機械能轉化為電信號。因此我們可以使用壓電陶瓷蜂鳴片作為一個無源傳感器。
（2）采集卡集成了信號調理模塊和 A/D 轉換模塊，各模塊的具體實現及功能為：① 信號調理模塊，先利用高通濾波電路（截止頻率 0.2 Hz）、低通濾波電路（截止頻率 159 Hz）和兩級放大電路（放大倍數分別為 11 倍和 1.21 倍）將采集信號濾波、放大和抬升，并經過鉗位電路使輸出電壓被鉗位在 — 0.7～4 V 之間，起到保護后級模數轉換器（analog to digital converter，ADC）芯片的作用；②A/D 轉換模塊，采用 TI 公司 12 位 ADC 芯片 ADS1013，標準 I2C 接口，低功耗，轉換速率為 1.6 ksps。將通過信號調理模塊后的模擬信號進行 A/D 轉換，將所得數字信號通過串口實時發送至上位機進行處理。
（3）終端處理利用 LabVIEW 應用軟件通過串口接收 A/D 轉換后的數字信號，經過處理分析，最終以波形圖的形式實時顯示所采集的原始生命體征信號。
上述枕下式睡眠監測系統采用單通道方式采集數據，采樣頻率為 1 000 Hz，波特率為 115 200，采集到的原始生命體征信號如圖 2 所示。由圖 2 可知，所得生命體征信號混雜著不同頻率、不同振幅的生物電信號以及不可避免的工頻干擾、設備儀器噪聲、人為運動等噪聲信號，我們首先通過時頻聯合分析法設計巴特沃斯帶阻濾波器、陷波濾波器以及低通濾波器等濾除噪聲并初步分離信號，再通過形態濾波去除基線漂移，最終從原始生命體征信號中分離出脈搏信號、呼吸信號與體動信號，如圖 3 所示。
圖1
系統總體設計簡圖
Figure1.
The overall design of the system
圖2
由枕下式睡眠監測系統采集的原始生命體征信號
Figure2.
The original vital signals collected by the sleep monitoring system under pillow
2 信號預處理
利用 HMM 對 HRV 進行分析，識別 HRV 在不同睡眠分期的模式變化，從而推算出相應的睡眠分期，需要注意的是，現有基于 HRV 的睡眠分期方案[21-22]，忽略了個體化差異會對睡眠分期分析產生影響。因此，所提分期模型僅適用于共性個體，難以實現針對年齡、性別存在明顯差異化個體的精準睡眠分期。本文首先對利用枕下式睡眠監測系統采集的生命體征信號以及獲取的 RRI 進行預處理，在最大程度上消除由于用戶性別、年齡等個體化差異帶來的分期影響。
2.1 RRI 提取
由枕下式睡眠監測系統采集的原始生命體征信號中分離出如圖 3 所示的脈搏信號、呼吸信號以及體動信號后，通過時頻聯合分析設計濾波器、峰值檢測、協方差匹配、自回歸模型算法等計算出準確的呼吸率、脈率，得到呼吸、脈搏的波形，并得到準確的的 RRI。整個 RRI 獲取的算法步驟如下所示：
圖3
經枕下式睡眠監測系統采集的原始生命體征信號中分離出的信號
a. 脈搏信號；b. 呼吸信號；c. 體動信號
Figure3. Signals separated from the original vital signals collected by the sleep monitoring system under pillowa. pulse signals; b. breath signals; c. movement signals
Algorithm
要求：從圖 2 所示的原始生命體征信號中分離出的脈搏信號、呼吸信號，按照每 10 s 劃分為一組
初始化：j = 1，k = 1，pulse ≠ φ
1：for i = 1 to Nn
2： corr（i） = cov（p（i），pulse）
3：end for
4：while（corr（j） ≠ φ）
5： Pmax = max{window（corr（j））}
6： ?P∈Pmax
7： if Tjmin < P < T jmax
8： 保留 P 并通過差分方程獲得初步 RRI
9： else
10： 去除 P
11： end if
12：end while
13：將得到的初步 RRI 信號按照每 6 s 劃分一組
14：while（RRI（n） ≠ φ）
15： ?R∈RRI（n）
16： Y（n） = a * R（n – 1） + a2 * R（n – 2） + a3 * R（n – 3）
17：其中 a 由以下公式獲得 a + a2 + a3 = 1
18： if R（n）*（1 – 30%）<|Y（n） – R（n）|<R（n）*（1 + 30%）
19： Y（n） 保留
20： else
21： Y（n） 去除
22： end if
23：end while
pulse 是從原始脈搏信號中選取的一段脈搏波形脈沖，將脈搏信號每 10 s 劃分成一組，p（i） 表示第 i 組的脈搏信號，Nn 是分成的組數。Corr（i）是 p（i）和 pulse 進行協方差匹配的結果。然后建立一個移動窗口，將窗口在相關函數中移動，選取每個窗口內的最大值為 Pmax。
6～12 行是去除不在閾值
和
之間的無關信號點，計算出初步 RRI，它表示將 RRI 以 6 s 為間期進行分組。14～23 行建模一個間期內有效峰值信息的 AR 模型，根據第 16 行公式由前三個點的值估計第四個點的值 Y（n），最后將 Y（n）與 6 s 間期內的第四個點的實際值 R（n）進行比較，將約束范圍設定為 30%，即當|Y（n）– X（n）|的值在 X（n）的 ± 30% 以內時，記錄此 RRI；如不滿足，視為失敗檢測。
2.2 差異化均衡
為消除因個體化差異導致的睡眠分期誤差，本文提出基于 EEMD[23-24]的 RRI 差異化預處理方案，可有效消除現有 EMD 模態混淆問題，具體實現步驟如下：
（1）提取并刪減異常 RRI，此類信息視為干擾，不做統計，包括過長（如 > 1 500 ms）、過短（如<300 ms），以及早搏前后的不規則 RRI。考慮到醫院多導睡眠圖（polysomnography，PSG 每 30 s 產生一個睡眠分期值，本文將排除異常干擾后的 RRI 設定為 30 s 一組，計算每一組的均值和方差。在排除異常 RRI 后，提出基于 EEDM 的個體差異均衡，具體步驟如下：
步驟 1：在原始信號 S（t）中，分別添加均值為零的白噪聲信號
和
，即：
 |
 |
其中
表示添加噪聲信號的幅值，i = 1，2，
，Ne，Ne 表示添加白噪聲的對數。
步驟 2：分別對
和
進行 EMD 分解，得到 IMF 分量序列和
，對上述得到的分量
進行集成：
 |
步驟 3：重復以上步驟直至 IMF 分量不是異常信號。
步驟 4：將已分解的前 p — 1 個分量從原始信號中分離出來，即：
 |
步驟 5：對剩余信號
進行 EEMD，將得到的所有 IMF 分量按高頻到低頻排列。
上述步驟中，EEMD 通過對原始信號多次加入不同的白噪聲進行 EMD 分解，將多次分解的結果進行平均即得到最終的 IMF。我們把去除異常值的 RRI 以 30 s 為周期分組計算出均值和方差以后，再將所得的 RRI 均值進行 EEMD，從而得到從高頻到低頻依次順序排列的 IMF 分量，最后 5 層低頻的 IMF 分量最能代表 RRI 均值的整體趨勢，于是將最后 5 層疊加來重構 RRI 均值的低頻從而得到其趨勢圖，結果如圖 4 所示（由于 PSG 每 30 s 產生一個睡眠分期值，因此圖中橫軸物理量采用與 PSG 相同的方式，標記為 u，即每 30 s 產生一個 RRI 均值）。
（2）在上述 EEMD 重構低頻得到趨勢信號后，將采用以下步驟進一步消除個體差異：
步驟 1：將 RRI 均值和 RRI 方差全部除以上述重構的趨勢信號，對 RRI 均值和方差歸一化。
步驟 2：為了減少觀察量的取值數目，對歸一化 RRI 均值和方差進行非均勻量化處理。預先設定好量化級數 S，將上述歸一化的值進行限幅編碼，從而得到與量化級數 S 相等的觀察量的取值數目。
步驟 3：在保證 HMM 準確度的前提下，盡量減少 HMM 輸入參數的個數，將由量化后的 RRI 均值和方差組成的二維訓練集通過線性關系映射到一維數組上，從而得到一組標識著 RRI 特征值的經過量化和歸一化的一維訓練集 RIS（i）。
3 基于 HMM 的睡眠分期
通過對枕下式睡眠監測系統 HRV 信息的差異化預處理，采用 HMM 對 HRV 進行分析，識別 HRV 在不同睡眠分期的模式變化，從而推算出相應的睡眠分期。HMM 是一種用參數表示的用于描述隨機過程統計特性的概率模型，是一個雙重隨機過程，由兩個部分組成：馬爾可夫鏈和一般隨機過程。其中，馬爾可夫鏈用來描述狀態的轉移，用轉移概率描述；一般隨機過程用來描述狀態與觀察序列間的關系，用觀察值概率描述。HMM 由初始狀態概率向量 pi、狀態轉移概率矩陣 A 和觀測概率矩陣 B 決定，pi 和 A 決定狀態序列，B 決定觀測序列。因此，HMM 的 lambda 可以用三元符號表示，即：lambda = （pi，A，B）。
總結起來，HMM 相關的算法主要分為三類，稱為 HMM 的三個基本問題：評估問題、解碼問題和學習問題，分別描述如下。
（1） 評估問題：已知觀察序列 O 和模型 lambda = （pi，A，B），如何有效計算在給定模型 lambda 條件下產生觀察序列 O 的條件概率即 P（O|lambda）；
圖4
重構低頻所得 RRI 趨勢信號
Figure4.
The trend RRI signals of reconstructing the low frequency
（2） 解碼問題：已知觀察序列 O 和模型 lambda = （pi，A，B），如何選擇相應的在某種意義上最佳（最好能解釋觀察序列）的狀態序列；
（3） 學習問題：如何調整模型參數 lambda = （pi，A，B），以使條件概率 P（O|lambda）最大。
本文采用 Baum-Welch[25]算法實現基于 HRV 分析的 HMM 睡眠分期問題。HMM 的輸入如下：
（1） PSG 采集的訓練集的睡眠分期結果：W、REM、NREM（S1，S2，S3，S4）；
（2） 前文預處理部分將 RRI 均值和方差經過歸一化以及量化后的一組標識著 RRI 特征值的一維訓練集 RIS（i）；
（3） 睡眠分期狀態總數 M；
（4） 觀察序列狀態總數 N。
在建立由以上輸入對應的 HMM 的 lambda = （pi，A，B） 時，睡眠分期狀態總數 M 定為 4，觀察序列狀態總數 N 為 9（由上述量化級數 S 而定），建立離散 HMM 的過程，為三個基本問題中的學習問題，即對于一個訓練樣本 RIS（i），確定 HMM 的參數 lambda = （pi，A，B），使得概率 P（RIS（i）|lambda）最大。具體步驟如下：
步驟 1：首先定義兩個輔助變量，第一個變量定義為 t 時狀態 i 和 t + 1 時狀態 j 的概率，即：
 |
第二個變量定義為后驗概率，即在給定觀察狀態序列和 HMM 的情況下，t 時狀態 i 的概率，即：
 |
步驟 2：根據公式（7）、（8）、（9）計算重估的 HMM 參數，即：
 |
 |
 |
步驟 3：重復步驟 1 和步驟 2，不斷地重新估計 HMM 的參數，在多次迭代以后得到 HMM 的一個最大似然估計，即為學習問題的最優解。
得到 HMM 模型以后，就得到狀態轉移矩陣 A 以及觀測概率矩陣 B，將測試集作為觀察序列 O，選擇在某種意義上最佳的狀態序列來解釋觀察序列，這個最佳的狀態序列即為最終測試集睡眠分期的估計結果，即對應 HMM 模型的第二個經典問題解碼問題。本文采用維特比算法來實現，具體步驟如下：
步驟 1：定義一個部分概率
和部分最優路徑。
為到達某個中間狀態的概率，部分最優路徑為到達每一個中間狀態和終止狀態的最可能的路徑。用
來表示在 t 時刻到達狀態 i 的所有可能的路徑中概率最大的序列的概率，部分最優路徑即為到達這個最大概率的路徑。
步驟 2：計算 t=1 時刻的部分概率。當 t=1 時，到達某個狀態最大可能的路徑還不存在，但是我們可以直接使用在 t=1 時刻某個狀態的概率和這個狀態可觀察序列 k1 的轉移概率
，即：
 |
步驟 3：計算 t > 1 時刻的部分概率。根據公式（11）由 t — 1 時刻的部分概率求 t 時刻的部分概率，其中
表示從狀態 j 到狀態 i 的概率，
表示狀態 i 到狀態
的概率。
 |
步驟 4：根據公式（12）利用一個后向指針
來記錄導致某個狀態最大局部概率的前一個狀態，即：
 |
步驟 5：重復步驟 3、4，即可找到最可能生成這個可觀察狀態序列的隱藏狀態序列，這個隱藏狀態序列即為某種意義上最佳的解釋觀察序列的狀態序列，即為最終測試集睡眠分期的估計結果。
通過以上 HMM 建模和識別的過程，最終得到測試集的睡眠分期的估計值。
4 結果
4.1 枕下采集準確性對比
本文將所提枕下式睡眠監測系統應用于廣州醫學院呼吸疾病研究所，同時與 PSG 測試數據進行同步比對，以驗證所提方案的可行性。
首先，驗證枕下非接觸壓電傳感方案獲取 RRI 的準確性，與 PSG 采集到的 HRV 信號的對比如圖 5 所示。由圖 5 所示，所設計的枕下心率檢測方案與 PSG 的心電信息在時間上保持嚴格同步（枕下式監測系統與 PSG 設備開機時間存在 10 s 左右的時延）。特別地，可以看到在圖 5a 中間部分心電信號發生了畸變，而與之相對應的圖 5b 中也同步出現了嚴重的波形變化。此外，心率所定位的 R 峰與 PSG 心電保持了良好的擬合度。
為了進一步驗證獲取 RRI 的準確性，按照前述 RRI 獲取的方法，得到一段 10 min 的 RRI，如圖 6 所示。對獲取的 RRI 進行傅里葉變換后的頻譜如圖 7 所示。由圖 6、7 可知，所設計的枕下式睡眠監測系統所獲取的 RRI 與 PSG 心電信息高度一致，誤差在 ± 50 ms 以內，經過傅里葉變換后得到的幅度譜和相位譜吻合率也非常之高，進一步表明本文所設計的枕下式睡眠監測系統可實現高準確度的心率 RRI 捕捉，可用于進一步基于 HRV 的深入分析。
圖5
同步采集 HRV 信號
a. PSG 采集；b. 傳感器采集
Figure5. Synchronous acquisition of HRV signalsa. by PSG; b. by sensor
圖6
PSG 與壓電傳感方案采集的 10 分鐘 RRI 信號
Figure6.
RRI signals of 10 minutes by PSG and sensor
4.2 睡眠分期
為了驗證所提分期方案的準確性，將本試驗結果與 PSG 分期金標準進行對比。為方便分析，試驗預設中將 NREM-1、2、3 合并為一種狀態，即共輸出四種睡眠分期的狀態，并采用 leave-one-out cross-validation 的交叉驗證方式。試驗數據由廣州中科新知科技有限公司授權，選取廣州醫學院呼吸疾病研究所 10 例不同年齡及性別的無睡眠障礙的院內正常受試者，10 名測試對象中，男、女各 5 例，最大年齡差異為 53 歲，其中年齡在 18 歲以下的有 3 人（1 男 2 女），在 18～50 歲之間的有 4 人（2 男 2 女），50 歲以上的有 3 人（2 男 1 女）。試驗過程為，采用前述枕下式睡眠監測系統分別采集 10 位測試對象夜間 HRV 信息，受試者同時采用 PSG 記錄當晚的睡眠分期結果。選擇第一位受試者作為測試集，剩余作為訓練集，通過 HMM 后得到估計的睡眠分期結果，然后將其與 PSG 采集的睡眠分期結果進行對比，結果如圖 8 所示，各分期時長及準確率如表 1 所示。經統計測算，分期準確率達 64.48%，證明本文非接觸的枕下睡眠監測分期與 PSG 睡眠分期相比，具備較高的擬合度。
依次改變測試集，分別將第 2，3，
，10 組作為測試集，剩余作為訓練集，經大量的重復性測試，得到的睡眠分期的準確率如表 2 所示。可以看出，本文所提出的 RRI 預處理消除個體化差異后的 HMM 睡眠分期方案準確率遠優于未歸一化組的 HMM 睡眠分期，以 PSG 為金標準，睡眠分期準確性超過 60％。上述數據對比試驗證明了本文所提方案的有效性。
圖7
對 RRI 進行傅里葉變換后的頻譜
Figure7.
Fourier transform spectrum of RRI
圖8
睡眠分期結果對比
Figure8.
Comparison of sleep staging results
5 總結與討論
本文設計了一種置于枕下的睡眠監測方案，采用枕下無擾采集裝置精確識別 HRV，在此基礎上，提出 EEMD 消除 HRV 個體化差異后，結合 HMM 進行睡眠分期監測。經過驗證，所提方案在對用戶無擾非接觸環境下獲取的 RRI 與 PSG 心電信息具有較高的擬和度，時域誤差范圍在 ± 50 ms 以內。以差異化預處理后的 HRV 為依據，采用 HMM 估計所得睡眠分期結果與 PSG 結果相比，在四期睡眠的劃分標準下，睡眠分期 S1～S4 正確匹配率達 61.51% ± 3%。為了驗證本文算法的有效性，同時驗證了文獻[12]中未經差異化預處理的 HMM 睡眠分期算法。經對比分析，本文所提方案睡眠分期的平均準確率具有較大的優勢，反映出不同個體差異化對睡眠分期建模具有較大的影響。
針對本文方案進一步的改進，包括：① 借助枕下傳感挖掘不同分期內呼吸、體動等參數特征，挖掘與 RRI 的相關性；借助此類相關特性，預期可以提高分析正確率，同時消除數據冗余，提高分析速度。② 從實驗和理論上尋找新的歸一化方法，在不同尺度下消除個體差異的影響，預期將進一步提升 HMM 建模的普適度與準確率。
表1
PSG 與枕下式睡眠監測系統睡眠分期結果對比表
Table1.
Comparison of the result of each sleep staging between PSG and sleep monitoring system under pillow
表2
枕下式睡眠監測系統數據歸一化前后睡眠分期結果對 比表
Table2.
Comparison of sleep staging results before and after data normalization based on sleep monitoring system under pillow
引言
睡眠分期結果可以對睡眠質量進行評估，是相關疾病診斷的重要指標[1-4]。20 世紀 60 年代美國加州大學腦研究所通過對腦電（electroencephalogram，EEG）、2 路頜肌電（electromyogram，EMG）以及 2 路眼電（electrooculogram，EOG）信號的分析研究，將睡眠狀態分為 6 期，即清醒期（wake，W）、快速眼動期（rapid eye movement，REM）、非快速眼動期（non-rapid eye movement，NREM），其中 NREM 期包括睡眠 Ⅰ 期（sleep stage 1，S1）、睡眠 Ⅱ 期（sleep stage 2，S2）、睡眠 Ⅲ 期（sleep stage 3，S3）、睡眠 Ⅳ 期（sleep stage 4，S4）。目前，EEG 已成為判斷睡眠分期的金標準，但是在此類監測方式下，受試者需長時間接觸多路電極，其自然睡眠會受到影響，不適于家居環境下的睡眠監護；而且研究所需的龐大數據量會耗費巨大的人力和時間成本。因此，探索其他生理信號在睡眠過程中的變化規律，并據此規律研究更加準確、高效的睡眠分期方法，是當前的研究熱點[5-8]。
心率的變化與自主神經系統的調控息息相關，自主神經系統的活動在不同睡眠期的表現也各不相同，現有研究[9-11]表明心率變異性（heart rate variability，HRV）中隱含著潛在的睡眠信息，具體表現為在上述 6 個分期階段心率變化呈現出不同的特點，通過合適的信號處理方法并結合先驗知識就可以反過來推算睡眠結構。Kortelainen 等[12]利用頻譜分析的方法對 RR 間期（RR interval，RRI）信號進行小波變換，將 RRI 通過頻譜變換分解成超低頻（very low frequency，VLF）、低頻（low frequency，LF）和高頻（high frequency，HF）三個頻譜區，經過對頻譜區特征的對比發現，LF/HF、VLF、VLF/HF 的均值與睡眠期存在線性關系，而且如果降低低頻信號的功率譜，能大大提高 S3 和 S4 期睡眠的識別率。莊志等[13]利用隱馬爾可夫模型（hidden Markov model，HMM）對 HRV 進行分析，在信號處理的過程中采用標準化的方式消除個體差異，對于睡眠分期的識別準確率可以接受。王金海等[14]利用支持向量機（support vector machin，SVM）實現了基于 HRV 信息的睡眠自動分期方法，采用主成分分析方法（principal component analysis，PCA）來降低數據冗余度，用少數幾個不相關的主成分替代原來相關互聯的眾多指標，獲得了較高的預測準確度。針對 HRV 的獲取方式，Wisse 等[15]利用壓力傳感的鉑電極傳感器獲得 HRV 信號，再結合體動信息通過 HMM 來預估睡眠分期。Ebrahimi 等[16]以睡眠心臟健康研究數據庫 HRV 信號為依據，提取特征值信息，再計算不同頻帶的能量和熵等指標對睡眠分期進行預估，能夠區分覺醒期和深睡期。
然而，目前睡眠分期方案存在以下 3 方面的問題：① 針對 HRV 信號的獲取，均需要用戶與多路電極直接接觸，用戶體驗度低[17-20]；② 國際公認的睡眠分期將睡眠狀態分為 6 期，現有睡眠分期方案僅分析對比睡眠/覺醒期，沒有針對 S1～S4 不同間期的睡眠分期對比驗證；③ 現有分期方案中，未考慮因用戶個體化差異帶來的分期誤差，分期輸出結果可信度相對較低。
針對上述問題，本文率先設計了一種枕下無擾式睡眠監控方案，主要貢獻可歸納為兩方面：
（1）針對 HRV 的獲取，提出采用壓電傳感器置于枕下的方式，無擾地采集用戶心率及 HRV 信息。所設計的采集平臺相較傳統院內、院外智能心率監測設備具有如下優勢：① 采樣率高達 1 000 Hz，達到醫療精準級別（目前院內心電采樣頻率 300～500 Hz），而市面可穿戴設備受功率制約采樣率極低，無法用于 HRV 的測量；② 新型非接觸無擾式檢測，不需與檢測對象直接接觸，相較現有院內專業采集設備需用戶與多路電極直接接觸的約束，具有明顯的用戶體驗度優勢。
（2）利用無擾心率檢測平臺進行睡眠分期的檢測，提出基于 HMM 的睡眠分期方案。與現有方案不同，為了消除因個體化差異導致的訓練輸出誤差，提出利用集合經驗模態分解（ensemble empirical mode decomposition，EEMD）去除基線漂移等歸一化預處理方案，消除 HRV 的個體化差異，同時，對歸一化后的信號進行限幅量化編碼，減少數據冗余度。所提優化方案將會在用戶體驗度和睡眠分期結果的準確性方面有顯著提升。
1 枕下無擾式實驗平臺
HRV 信號的采集通過枕下無擾式實驗平臺來實現，系統整體設計由枕下壓電傳感、信號采集卡以及終端處理三部分組成，設計簡圖如圖 1 所示。傳感器置于一塊平板中，放于枕頭底下，集成了印制電路板（printed circuit board，PCB）的采集卡放置于床邊，采集信號經過處理后在終端顯示。各部分的具體實現及功能如下：
（1）枕下壓電傳感內置采集模塊（系統的“感受裝置”），主要通過壓電傳感器來實現，通過比較后我們選擇的傳感器是壓電陶瓷蜂鳴片，它具有厚度薄、體積小、穩定性好、靈敏度高、耗能低、耐高低溫、壽命長、耐高壓、無需供電等優點，可以使信號采集模塊在枕下發揮功效。當壓電陶瓷蜂鳴片被施加直流電壓時，產生交替重復的機械彎曲變形，從而在空氣中產生聲波；然而當有外力施加時，壓電陶瓷片會受外力影響產生電波，將機械能轉化為電信號。因此我們可以使用壓電陶瓷蜂鳴片作為一個無源傳感器。
（2）采集卡集成了信號調理模塊和 A/D 轉換模塊，各模塊的具體實現及功能為：① 信號調理模塊，先利用高通濾波電路（截止頻率 0.2 Hz）、低通濾波電路（截止頻率 159 Hz）和兩級放大電路（放大倍數分別為 11 倍和 1.21 倍）將采集信號濾波、放大和抬升，并經過鉗位電路使輸出電壓被鉗位在 — 0.7～4 V 之間，起到保護后級模數轉換器（analog to digital converter，ADC）芯片的作用；②A/D 轉換模塊，采用 TI 公司 12 位 ADC 芯片 ADS1013，標準 I2C 接口，低功耗，轉換速率為 1.6 ksps。將通過信號調理模塊后的模擬信號進行 A/D 轉換，將所得數字信號通過串口實時發送至上位機進行處理。
（3）終端處理利用 LabVIEW 應用軟件通過串口接收 A/D 轉換后的數字信號，經過處理分析，最終以波形圖的形式實時顯示所采集的原始生命體征信號。
上述枕下式睡眠監測系統采用單通道方式采集數據，采樣頻率為 1 000 Hz，波特率為 115 200，采集到的原始生命體征信號如圖 2 所示。由圖 2 可知，所得生命體征信號混雜著不同頻率、不同振幅的生物電信號以及不可避免的工頻干擾、設備儀器噪聲、人為運動等噪聲信號，我們首先通過時頻聯合分析法設計巴特沃斯帶阻濾波器、陷波濾波器以及低通濾波器等濾除噪聲并初步分離信號，再通過形態濾波去除基線漂移，最終從原始生命體征信號中分離出脈搏信號、呼吸信號與體動信號，如圖 3 所示。
圖1
系統總體設計簡圖
Figure1.
The overall design of the system
圖2
由枕下式睡眠監測系統采集的原始生命體征信號
Figure2.
The original vital signals collected by the sleep monitoring system under pillow
2 信號預處理
利用 HMM 對 HRV 進行分析，識別 HRV 在不同睡眠分期的模式變化，從而推算出相應的睡眠分期，需要注意的是，現有基于 HRV 的睡眠分期方案[21-22]，忽略了個體化差異會對睡眠分期分析產生影響。因此，所提分期模型僅適用于共性個體，難以實現針對年齡、性別存在明顯差異化個體的精準睡眠分期。本文首先對利用枕下式睡眠監測系統采集的生命體征信號以及獲取的 RRI 進行預處理，在最大程度上消除由于用戶性別、年齡等個體化差異帶來的分期影響。
2.1 RRI 提取
由枕下式睡眠監測系統采集的原始生命體征信號中分離出如圖 3 所示的脈搏信號、呼吸信號以及體動信號后，通過時頻聯合分析設計濾波器、峰值檢測、協方差匹配、自回歸模型算法等計算出準確的呼吸率、脈率，得到呼吸、脈搏的波形，并得到準確的的 RRI。整個 RRI 獲取的算法步驟如下所示：
圖3
經枕下式睡眠監測系統采集的原始生命體征信號中分離出的信號
a. 脈搏信號；b. 呼吸信號；c. 體動信號
Figure3. Signals separated from the original vital signals collected by the sleep monitoring system under pillowa. pulse signals; b. breath signals; c. movement signals
Algorithm
要求：從圖 2 所示的原始生命體征信號中分離出的脈搏信號、呼吸信號，按照每 10 s 劃分為一組
初始化：j = 1，k = 1，pulse ≠ φ
1：for i = 1 to Nn
2： corr（i） = cov（p（i），pulse）
3：end for
4：while（corr（j） ≠ φ）
5： Pmax = max{window（corr（j））}
6： ?P∈Pmax
7： if Tjmin < P < T jmax
8： 保留 P 并通過差分方程獲得初步 RRI
9： else
10： 去除 P
11： end if
12：end while
13：將得到的初步 RRI 信號按照每 6 s 劃分一組
14：while（RRI（n） ≠ φ）
15： ?R∈RRI（n）
16： Y（n） = a * R（n – 1） + a2 * R（n – 2） + a3 * R（n – 3）
17：其中 a 由以下公式獲得 a + a2 + a3 = 1
18： if R（n）*（1 – 30%）<|Y（n） – R（n）|<R（n）*（1 + 30%）
19： Y（n） 保留
20： else
21： Y（n） 去除
22： end if
23：end while
pulse 是從原始脈搏信號中選取的一段脈搏波形脈沖，將脈搏信號每 10 s 劃分成一組，p（i） 表示第 i 組的脈搏信號，Nn 是分成的組數。Corr（i）是 p（i）和 pulse 進行協方差匹配的結果。然后建立一個移動窗口，將窗口在相關函數中移動，選取每個窗口內的最大值為 Pmax。
6～12 行是去除不在閾值
和
之間的無關信號點，計算出初步 RRI，它表示將 RRI 以 6 s 為間期進行分組。14～23 行建模一個間期內有效峰值信息的 AR 模型，根據第 16 行公式由前三個點的值估計第四個點的值 Y（n），最后將 Y（n）與 6 s 間期內的第四個點的實際值 R（n）進行比較，將約束范圍設定為 30%，即當|Y（n）– X（n）|的值在 X（n）的 ± 30% 以內時，記錄此 RRI；如不滿足，視為失敗檢測。
2.2 差異化均衡
為消除因個體化差異導致的睡眠分期誤差，本文提出基于 EEMD[23-24]的 RRI 差異化預處理方案，可有效消除現有 EMD 模態混淆問題，具體實現步驟如下：
（1）提取并刪減異常 RRI，此類信息視為干擾，不做統計，包括過長（如 > 1 500 ms）、過短（如<300 ms），以及早搏前后的不規則 RRI。考慮到醫院多導睡眠圖（polysomnography，PSG 每 30 s 產生一個睡眠分期值，本文將排除異常干擾后的 RRI 設定為 30 s 一組，計算每一組的均值和方差。在排除異常 RRI 后，提出基于 EEDM 的個體差異均衡，具體步驟如下：
步驟 1：在原始信號 S（t）中，分別添加均值為零的白噪聲信號
和
，即：
|
|
其中
表示添加噪聲信號的幅值，i = 1，2，
，Ne，Ne 表示添加白噪聲的對數。
步驟 2：分別對
和
進行 EMD 分解，得到 IMF 分量序列和
，對上述得到的分量
進行集成：
|
步驟 3：重復以上步驟直至 IMF 分量不是異常信號。
步驟 4：將已分解的前 p — 1 個分量從原始信號中分離出來，即：
|
步驟 5：對剩余信號
進行 EEMD，將得到的所有 IMF 分量按高頻到低頻排列。
上述步驟中，EEMD 通過對原始信號多次加入不同的白噪聲進行 EMD 分解，將多次分解的結果進行平均即得到最終的 IMF。我們把去除異常值的 RRI 以 30 s 為周期分組計算出均值和方差以后，再將所得的 RRI 均值進行 EEMD，從而得到從高頻到低頻依次順序排列的 IMF 分量，最后 5 層低頻的 IMF 分量最能代表 RRI 均值的整體趨勢，于是將最后 5 層疊加來重構 RRI 均值的低頻從而得到其趨勢圖，結果如圖 4 所示（由于 PSG 每 30 s 產生一個睡眠分期值，因此圖中橫軸物理量采用與 PSG 相同的方式，標記為 u，即每 30 s 產生一個 RRI 均值）。
（2）在上述 EEMD 重構低頻得到趨勢信號后，將采用以下步驟進一步消除個體差異：
步驟 1：將 RRI 均值和 RRI 方差全部除以上述重構的趨勢信號，對 RRI 均值和方差歸一化。
步驟 2：為了減少觀察量的取值數目，對歸一化 RRI 均值和方差進行非均勻量化處理。預先設定好量化級數 S，將上述歸一化的值進行限幅編碼，從而得到與量化級數 S 相等的觀察量的取值數目。
步驟 3：在保證 HMM 準確度的前提下，盡量減少 HMM 輸入參數的個數，將由量化后的 RRI 均值和方差組成的二維訓練集通過線性關系映射到一維數組上，從而得到一組標識著 RRI 特征值的經過量化和歸一化的一維訓練集 RIS（i）。
3 基于 HMM 的睡眠分期
通過對枕下式睡眠監測系統 HRV 信息的差異化預處理，采用 HMM 對 HRV 進行分析，識別 HRV 在不同睡眠分期的模式變化，從而推算出相應的睡眠分期。HMM 是一種用參數表示的用于描述隨機過程統計特性的概率模型，是一個雙重隨機過程，由兩個部分組成：馬爾可夫鏈和一般隨機過程。其中，馬爾可夫鏈用來描述狀態的轉移，用轉移概率描述；一般隨機過程用來描述狀態與觀察序列間的關系，用觀察值概率描述。HMM 由初始狀態概率向量 pi、狀態轉移概率矩陣 A 和觀測概率矩陣 B 決定，pi 和 A 決定狀態序列，B 決定觀測序列。因此，HMM 的 lambda 可以用三元符號表示，即：lambda = （pi，A，B）。
總結起來，HMM 相關的算法主要分為三類，稱為 HMM 的三個基本問題：評估問題、解碼問題和學習問題，分別描述如下。
（1） 評估問題：已知觀察序列 O 和模型 lambda = （pi，A，B），如何有效計算在給定模型 lambda 條件下產生觀察序列 O 的條件概率即 P（O|lambda）；
圖4
重構低頻所得 RRI 趨勢信號
Figure4.
The trend RRI signals of reconstructing the low frequency
（2） 解碼問題：已知觀察序列 O 和模型 lambda = （pi，A，B），如何選擇相應的在某種意義上最佳（最好能解釋觀察序列）的狀態序列；
（3） 學習問題：如何調整模型參數 lambda = （pi，A，B），以使條件概率 P（O|lambda）最大。
本文采用 Baum-Welch[25]算法實現基于 HRV 分析的 HMM 睡眠分期問題。HMM 的輸入如下：
（1） PSG 采集的訓練集的睡眠分期結果：W、REM、NREM（S1，S2，S3，S4）；
（2） 前文預處理部分將 RRI 均值和方差經過歸一化以及量化后的一組標識著 RRI 特征值的一維訓練集 RIS（i）；
（3） 睡眠分期狀態總數 M；
（4） 觀察序列狀態總數 N。
在建立由以上輸入對應的 HMM 的 lambda = （pi，A，B） 時，睡眠分期狀態總數 M 定為 4，觀察序列狀態總數 N 為 9（由上述量化級數 S 而定），建立離散 HMM 的過程，為三個基本問題中的學習問題，即對于一個訓練樣本 RIS（i），確定 HMM 的參數 lambda = （pi，A，B），使得概率 P（RIS（i）|lambda）最大。具體步驟如下：
步驟 1：首先定義兩個輔助變量，第一個變量定義為 t 時狀態 i 和 t + 1 時狀態 j 的概率，即：
|
第二個變量定義為后驗概率，即在給定觀察狀態序列和 HMM 的情況下，t 時狀態 i 的概率，即：
|
步驟 2：根據公式（7）、（8）、（9）計算重估的 HMM 參數，即：
|
|
|
步驟 3：重復步驟 1 和步驟 2，不斷地重新估計 HMM 的參數，在多次迭代以后得到 HMM 的一個最大似然估計，即為學習問題的最優解。
得到 HMM 模型以后，就得到狀態轉移矩陣 A 以及觀測概率矩陣 B，將測試集作為觀察序列 O，選擇在某種意義上最佳的狀態序列來解釋觀察序列，這個最佳的狀態序列即為最終測試集睡眠分期的估計結果，即對應 HMM 模型的第二個經典問題解碼問題。本文采用維特比算法來實現，具體步驟如下：
步驟 1：定義一個部分概率
和部分最優路徑。
為到達某個中間狀態的概率，部分最優路徑為到達每一個中間狀態和終止狀態的最可能的路徑。用
來表示在 t 時刻到達狀態 i 的所有可能的路徑中概率最大的序列的概率，部分最優路徑即為到達這個最大概率的路徑。
步驟 2：計算 t=1 時刻的部分概率。當 t=1 時，到達某個狀態最大可能的路徑還不存在，但是我們可以直接使用在 t=1 時刻某個狀態的概率和這個狀態可觀察序列 k1 的轉移概率
，即：
|
步驟 3：計算 t > 1 時刻的部分概率。根據公式（11）由 t — 1 時刻的部分概率求 t 時刻的部分概率，其中
表示從狀態 j 到狀態 i 的概率，
表示狀態 i 到狀態
的概率。
|
步驟 4：根據公式（12）利用一個后向指針
來記錄導致某個狀態最大局部概率的前一個狀態，即：
|
步驟 5：重復步驟 3、4，即可找到最可能生成這個可觀察狀態序列的隱藏狀態序列，這個隱藏狀態序列即為某種意義上最佳的解釋觀察序列的狀態序列，即為最終測試集睡眠分期的估計結果。
通過以上 HMM 建模和識別的過程，最終得到測試集的睡眠分期的估計值。
4 結果
4.1 枕下采集準確性對比
本文將所提枕下式睡眠監測系統應用于廣州醫學院呼吸疾病研究所，同時與 PSG 測試數據進行同步比對，以驗證所提方案的可行性。
首先，驗證枕下非接觸壓電傳感方案獲取 RRI 的準確性，與 PSG 采集到的 HRV 信號的對比如圖 5 所示。由圖 5 所示，所設計的枕下心率檢測方案與 PSG 的心電信息在時間上保持嚴格同步（枕下式監測系統與 PSG 設備開機時間存在 10 s 左右的時延）。特別地，可以看到在圖 5a 中間部分心電信號發生了畸變，而與之相對應的圖 5b 中也同步出現了嚴重的波形變化。此外，心率所定位的 R 峰與 PSG 心電保持了良好的擬合度。
為了進一步驗證獲取 RRI 的準確性，按照前述 RRI 獲取的方法，得到一段 10 min 的 RRI，如圖 6 所示。對獲取的 RRI 進行傅里葉變換后的頻譜如圖 7 所示。由圖 6、7 可知，所設計的枕下式睡眠監測系統所獲取的 RRI 與 PSG 心電信息高度一致，誤差在 ± 50 ms 以內，經過傅里葉變換后得到的幅度譜和相位譜吻合率也非常之高，進一步表明本文所設計的枕下式睡眠監測系統可實現高準確度的心率 RRI 捕捉，可用于進一步基于 HRV 的深入分析。
圖5
同步采集 HRV 信號
a. PSG 采集；b. 傳感器采集
Figure5. Synchronous acquisition of HRV signalsa. by PSG; b. by sensor
圖6
PSG 與壓電傳感方案采集的 10 分鐘 RRI 信號
Figure6.
RRI signals of 10 minutes by PSG and sensor
4.2 睡眠分期
為了驗證所提分期方案的準確性，將本試驗結果與 PSG 分期金標準進行對比。為方便分析，試驗預設中將 NREM-1、2、3 合并為一種狀態，即共輸出四種睡眠分期的狀態，并采用 leave-one-out cross-validation 的交叉驗證方式。試驗數據由廣州中科新知科技有限公司授權，選取廣州醫學院呼吸疾病研究所 10 例不同年齡及性別的無睡眠障礙的院內正常受試者，10 名測試對象中，男、女各 5 例，最大年齡差異為 53 歲，其中年齡在 18 歲以下的有 3 人（1 男 2 女），在 18～50 歲之間的有 4 人（2 男 2 女），50 歲以上的有 3 人（2 男 1 女）。試驗過程為，采用前述枕下式睡眠監測系統分別采集 10 位測試對象夜間 HRV 信息，受試者同時采用 PSG 記錄當晚的睡眠分期結果。選擇第一位受試者作為測試集，剩余作為訓練集，通過 HMM 后得到估計的睡眠分期結果，然后將其與 PSG 采集的睡眠分期結果進行對比，結果如圖 8 所示，各分期時長及準確率如表 1 所示。經統計測算，分期準確率達 64.48%，證明本文非接觸的枕下睡眠監測分期與 PSG 睡眠分期相比，具備較高的擬合度。
依次改變測試集，分別將第 2，3，
，10 組作為測試集，剩余作為訓練集，經大量的重復性測試，得到的睡眠分期的準確率如表 2 所示。可以看出，本文所提出的 RRI 預處理消除個體化差異后的 HMM 睡眠分期方案準確率遠優于未歸一化組的 HMM 睡眠分期，以 PSG 為金標準，睡眠分期準確性超過 60％。上述數據對比試驗證明了本文所提方案的有效性。
圖7
對 RRI 進行傅里葉變換后的頻譜
Figure7.
Fourier transform spectrum of RRI
圖8
睡眠分期結果對比
Figure8.
Comparison of sleep staging results
5 總結與討論
本文設計了一種置于枕下的睡眠監測方案，采用枕下無擾采集裝置精確識別 HRV，在此基礎上，提出 EEMD 消除 HRV 個體化差異后，結合 HMM 進行睡眠分期監測。經過驗證，所提方案在對用戶無擾非接觸環境下獲取的 RRI 與 PSG 心電信息具有較高的擬和度，時域誤差范圍在 ± 50 ms 以內。以差異化預處理后的 HRV 為依據，采用 HMM 估計所得睡眠分期結果與 PSG 結果相比，在四期睡眠的劃分標準下，睡眠分期 S1～S4 正確匹配率達 61.51% ± 3%。為了驗證本文算法的有效性，同時驗證了文獻[12]中未經差異化預處理的 HMM 睡眠分期算法。經對比分析，本文所提方案睡眠分期的平均準確率具有較大的優勢，反映出不同個體差異化對睡眠分期建模具有較大的影響。
針對本文方案進一步的改進，包括：① 借助枕下傳感挖掘不同分期內呼吸、體動等參數特征，挖掘與 RRI 的相關性；借助此類相關特性，預期可以提高分析正確率，同時消除數據冗余，提高分析速度。② 從實驗和理論上尋找新的歸一化方法，在不同尺度下消除個體差異的影響，預期將進一步提升 HMM 建模的普適度與準確率。
表1
PSG 與枕下式睡眠監測系統睡眠分期結果對比表
Table1.
Comparison of the result of each sleep staging between PSG and sleep monitoring system under pillow
表2
枕下式睡眠監測系統數據歸一化前后睡眠分期結果對 比表
Table2.
Comparison of sleep staging results before and after data normalization based on sleep monitoring system under pillow
