睡眠質量與人體健康密切相關,睡眠狀態的正確分期對于睡眠質量的評估、睡眠相關疾病的診斷分析都具有重要的意義。多導睡眠(PSG)信號常用于記錄和分析睡眠狀態,而有效的特征提取和表達方法是提高睡眠分期準確性的重要環節。本文提出采用協作表示(CR)方法重新表達腦電(EEG)信號的特征,并進一步采用核熵成分分析(KECA)對CR特征進行降維。為驗證本文提出算法的有效性,將CR-KECA特征與原始特征、CR特征和CR經過主成分分析(PCA)降維后的特征進行對比,實驗結果表明:基于CR-KECA方法的特征獲得了最好的睡眠狀態分期性能,其準確率、敏感度和特異性分別為68.74%±0.46%、68.76%±0.43%、92.19%±0.11%,明顯優于CR-PCA特征、CR特征和原始特征;同時,CR算法計算復雜度非常低,而且KECA降維后的特征維度更小,因而適于臨床海量睡眠數據的分析。
引用本文： 趙攀博, 施俊, 劉瀟, 江其坤, 顧宇. 基于協作表示與核熵成分分析的腦電睡眠狀態分期研究. 生物醫學工程學雜志, 2015, 32(4): 730-734. doi: 10.7507/1001-5515.20150133 復制
引言
睡眠是人體一種重要的生理現象，人的一生約有三分之一的時間是在睡眠中度過的。在睡眠過程中，睡眠不是一種固定不變的狀態，而是經歷著幾個相對穩定的階段，臨床上把這幾個階段定義為睡眠狀態的分期。根據美國睡眠醫學學會2007年公布的睡眠分期標準，把睡眠階段分成了5期，分別為：覺醒期、睡眠1期、睡眠2期、慢波睡眠期、非快速眼動期[1]。由于睡眠質量關系到人體健康，而睡眠分期的分析又與睡眠質量密切相關，因此，睡眠分期的研究分析對于預防和診斷睡眠相關疾病具有重要的臨床應用價值。
多導睡眠圖(polysomnography，PSG)在臨床上廣泛應用于記錄睡眠狀態，然后對所獲取的信號進行分析和診斷。PSG包括腦電圖(electroencephalogram，EEG)[2-4]、心電圖[4-6]、眼動圖[4, 7]和頰肌電圖[4]等多種信號。由醫師對長時間記錄的睡眠信號進行人工分析，不僅耗時耗力，而且易于引起漏診和誤診。因此，基于計算機的睡眠數據輔助分析是目前睡眠分析的研究熱點。
PSG信號的有效特征提取是實現睡眠分期的重要環節。目前，常用的特征提取方法包括時域、頻域和時頻域特征提取方法[3, 5, 8]。然而，基于現有特征提取方法的睡眠分期的效果還有待進一步的提高。
協作表示(collaborative representation，CR)是新提出的一種數據表示方法[9]，是一種與目前在信號處理領域取得非常成功應用的稀疏表示(sparse representation，SR)相對應的方法。SR旨在采用冗余字典中的少量字典原子的線性組合來表示信號，常采用l1-范數作為稀疏約束；而CR則是采用全部樣本的組合來表示信號，采用l2-范數作為約束條件，通過l2-范數求解系數并將不同類樣本協作共同用作訓練字典。研究還表明，是CR而非l1-范數稀疏約束使得基于SR的分類取得了巨大的成功[9]。此外，CR的計算復雜度遠遠小于基于l1-范數的SR。因此，基于CR的分類取得了與基于SR的分類相當的效果[9-11]。但是，目前CR主要是作為分類器來應用，鮮有用于特征提取，更未用于PSG信號處理。
另一方面，經過CR以后獲取的新的數據向量表示，其特征維度往往很高，而且由于CR是全部樣本的組合，其編碼系數在一定程度上還有冗余或者無關信息。因此，對CR獲得的編碼系數進行降維是降低其后續處理計算復雜度、進一步提高數據表達能力的有效方法。核熵成分分析(kernel entropy component analysis，KECA)是一種新的基于信息論的數據變換方法[12]。其特質在于尋找高維特征空間中能夠最大程度保持原始空間數據二次Renyi熵值的坐標軸作為投影方向。傳統的數據變換譜方法，如主成分分析(principal component analysis，PCA)等，所選擇的是數據變換后矩陣的前n個最大特征值所對應的特征作為降維后的特征，而KECA則與之不同，是根據熵值來確定降維后的特征，這就最大程度地保證了數據內在蘊含的結構屬性。因此，KECA已經在高維數據的降維、聚類等應用中表現出了甚至優于PCA的性能[12-14]。
本文針對基于PSG的睡眠狀態分期中的特征提取問題，采用CR對特征進行再表達，將CR編碼系數作為新的特征，然后通過KECA對CR特征進行降維，獲取更為緊致和更具表征能力的特征。
1 協作表示
對于SR和CR，都可以用以下一般模型來描述：
| $\hat \alpha = \arg {\min _\alpha }\left\{ {{{\left\| {{\rm{f}} - \boldsymbol{A}\alpha } \right\|}_{{l_q}}} + \lambda {{\left\| \alpha \right\|}_{{l_q}}}} \right\}$ |
其中，f是待分析信號，A是冗余字典，α是稀疏系數，λ是正則參數，p,q=1或者2，對于p，q不同的賦值可以得到不同的具體的協作表示模型。在SR中，p=1，q=2。
雖然l1-范數約束是SR常用的方法，但是l1-范數最小化問題求解復雜，計算復雜度高。Zhang等[9]用l2-范數即最小二乘正則化取代l1-范數，提出了CR的概念，其模型如下:
| $\hat \alpha = \arg {\min _\alpha }\left\{ {\left\| {{\rm{f}} - \boldsymbol{A}\alpha } \right\|_2^2 + \lambda \left\| \alpha \right\|_2^2} \right\}$ |
最小二乘正則項有兩個作用，一是使得最小二乘解穩定，二是求解時引入一定量的稀疏。 對式(2)求解問題，很容易求得唯一解：
| $\hat \alpha = \boldsymbol{P}{\rm{f}}$ |
其中，
| $\boldsymbol{P} = {\left( {{\boldsymbol{A}^{\rm{T}}}\boldsymbol{A} + \lambda \times \boldsymbol{I}} \right)^{ - 1}}{\boldsymbol{A}^{\rm{T}}}$ |
顯然，P獨立于信號f，可以作為投影矩陣。測試時，只需將信號f投影到P空間下，即可得到向量系數，這使得編碼效率大大提高。
在Zhang等[9]提出的協作表示模型中，冗余字典A往往直接采用訓練樣本替代。由于K-means聚類算法被廣泛應用于字典訓練，在本文中，采用K-means聚類算法來訓練字典A。
2 核熵成分分析
KECA的原理如下：
給定數據集S=x1,…,xN，p(x)是其概率密度函數，則其二次Renyi熵為：
| ${\rm{H}}\left( {\rm{p}} \right) = - \log \int {{{\rm{p}}^2}\left( x \right)} {\rm{d}}x$ |
由于對數函數的單調性，因此主要集中定性分析。采用如公式(6)所示的Parzen窗密度估計器來估計V(p)的核函數：
| ${\rm{\hat p}}\left( x \right) = \frac{1}{N}\sum\limits_{{x_t} \in S} {{{\rm{K}}_\sigma }\left( {x,{x_t}} \right)} $ |
其中，Kσ(x,xt)是Parzen窗或者稱之為以xt為中心、以σ為寬度的核函數。使用(x)的均值對V(p)進行估計，可以得到：
| ${\rm{\hat V}}\left( {\rm{p}} \right) = \frac{1}{N}\sum\limits_{{x_t} \in S} {{\rm{\hat p}}\left( {{x_t}} \right)} = \frac{1}{N}\sum\limits_{{x_t} \in S} {\frac{1}{N}\sum\limits_{{x_t} \in S} {{{\rm{K}}_\sigma }\left( {{x_t},{x_{t'}}} \right) = } } \frac{1}{{{N^2}}}{\boldsymbol{I}^{\rm{T}}}\boldsymbol{K}\boldsymbol{I}$ |
其中，I為(N×1)的單位向量，K為(N×N)的核矩陣，元素(t,t′)的核矩陣可表示為Kσ(xt,xt′)。那么，Renyi熵估計的核矩陣可以用特征值和特征向量來表示，即：K=EDET，其中，D是特征值為[λ1,λ2,…,λN]的對角矩陣，E是相應的特征向量，即：E=[e1,e2,…,eN]。那么，公式(7)可以表示為：
| ${\rm{\hat V}}\left( {\rm{p}} \right) = \frac{1}{{{N^2}}}\sum\limits_{i = 1}^N {{{\left( {\sqrt {{\lambda _i}} e_i^T\boldsymbol{I}} \right)}^2}} $ |
假設k(k<N)維數據通過Φ映射到子空間Uk，當且僅當子空間與Renyi熵建立聯系時，根據熵的大小將特征值和特征向量進行重新排序，產生KECA的映射Φeca
| ${\boldsymbol{\Phi} _{eca}} = {{\rm{\boldsymbol{U}}}_k}\boldsymbol{\Phi} = \boldsymbol{D}_k^{\frac{1}{2}}\boldsymbol{E}_k^T$ |
公式(9)的求解可以轉換成下面公式的最小化問題：
| ${\boldsymbol{\Phi} _{eca}} = \boldsymbol{D}_k^{\frac{1}{2}}\boldsymbol{E}_k^T:\min {\rm{\hat V}}\left( {\rm{p}} \right) - {{\rm{\hat V}}_k}\left( {\rm{p}} \right)$ |
結合公式(8)可知：
| ${{\rm{\hat V}}_k}\left( {\rm{p}} \right) = \frac{1}{{{N^2}}}\sum\limits_{i = 1}^k {{{\left( {\sqrt {{\lambda _i}} e_i^T\boldsymbol{I}} \right)}^2}} $ |
因此，可以發現對二次Renyi熵貢獻較高的是特定的特征值和相應的特征向量，而不是前面的特征值和相應的特征向量。顯然，KECA是通過選取對二次Renyi熵貢獻較高的特征值和相應的特征向量來盡可能地保持數據變換前后熵的不變性。
3 實驗與結果分析
3.1 數據庫
為了驗證本文提出的CR結合KECA的特征提取算法的有效性，本文采用St. Vincent’s University Hospital和University College Dublin提供的睡眠呼吸暫停數據庫[15]。這個數據庫包含25名疑似睡眠呼吸紊亂的被試者的整夜PSG數據，其中(C3-A2)通道的EEG信號被選擇用于本文的測試研究。將EEG信號每30秒分割成一個片段，每個片段對應于5種睡眠狀態的一個分期階段。
對每個片段進行如下計算來提取特征，獲取共計78個特征值：
(1)計算信號的均值、標準差、中值、幅值(信號最大值與最小值絕對值差)；
(2)計算信號的峰態系數和信息熵；
(3)計算信號譜密度線性擬合的負斜率；
(4)根據10個不同頻段(δ1,δ2,θ1,θ2,α1,α2,σ1,σ2,β,γ)運用快速傅里葉變換提取特征；
(5)對(4)中求得的信號最大值和最大值所在的頻帶作為特征；
(6)用傅里葉變換在五個頻段(α,β,δ,θ和高頻)提取特征；
(7)對(4)中得到的特征向量進行頻帶劃分，計算頻帶間比率作為特征；
(8)小波變換提取信號特征：計算六個頻段(α,β,δ,θ,高γ，低γ)的小波系數，再對每個頻段的小波系數計算最大值、最小值、標準差、均值、中值、四分位數、半四分位數、熵和能量。
3.2 實驗設計
為了驗證本文提出的CR結合KECA降維所提取特征(下文用CR-KECA特征表示)的性能，與如下3種特征提取算法進行對比實驗：①EEG信號的原始特征；②原始特征經過CR編碼以后獲得的特征(CR特征)；③對CR特征采用最常用的PCA降維以后獲取的特征(CR-PCA特征)。
由于所有樣本的片段數量巨大，我們隨機挑選了8 000個樣本數據用于本實驗，每一類(睡眠分期)的樣本數量平均分布。其中，隨機選擇5 000個樣本采用K-means聚類算法訓練字典，字典原子個數為312，即CR特征的維度為312。剩余的3 000個樣本中，隨機選擇1 000個樣本作為訓練集，2 000個樣本作為測試。上述實驗過程重復5次。分類器采用LIBSVM程序包中的線性支持向量機(support vector machine，SVM)[16]。
分類評價指標選擇分類精度(accuracy，ACC)、敏感度(sensitivity，SEN)、特異性(specificity,SPE)，計算公式如下：
| $\left\{ \begin{array}{l} {\rm{ACC = }}\frac{{{\rm{TP + TN}}}}{{{\rm{TP + FN + TN + FP}}}}\\ {\rm{SEN = }}\frac{{{\rm{TP}}}}{{{\rm{TP + FN}}}}\\ {\rm{SPE = }}\frac{{{\rm{TN}}}}{{{\rm{TN + FP}}}} \end{array} \right.$ |
其中，TP、FN、TN、FP分別表示真陽性、假陰性、真陰性和假陽性。五次實驗的平均結果作為最終的性能指標(均值±標準差)，并且進行t檢驗，當P值小于0.05時，差異有統計學意義。
3.3 結果及分析
表 1所示為CR-KECA特征在不同維度下的分類結果。其中，分類性能最高的結果出現在60維，其分類精度、敏感度和特異性分別為68.74%±0.46%、68.76%±0.43%、92.19%±0.11%。從表中還可以發現，從5維特征開始到30維，隨著維數的增加，準確率、敏感度和特異性均呈顯著提高，但是30維以后的特征維度下，結果相差不大。
表1
不同維度下的CR-KECA特征的分類精度、敏感度及特異性(%)
Table1.
Classification accuracy,sensitivity and specificity of CR-KECA feature with different dimensions (%)
為了直觀比較CR-KECA特征與其他三種特征的性能，圖 1依次顯示了原始特征、CR特征、CR-PCA特征和CR-KECA特征的平均分類精度、敏感度和特異性。由圖可以發現，CR-KECA的分類精度較原始特征、CR特征、CR-PCA特征分別提高了2.92%、1.47%、0.96%，敏感度分別提高了2.94%、1.49%、0.98%，特異性分別提高了0.73%、0.37%、0.24%。此外，CR-KECA特征與原始特征和CR特征之間在分類精度、敏感度和特異性這三個指標上的差異都具有統計學意義(P<0.05)，表明CR-KECA具有顯著的性能提高。另一方面，與原始的特征相比，CR特征的分類精度、敏感性和特異性較原始特征也有較大的提高，分別為1.45%、1.46%、0.36%。
圖1
四種特征分類的精度、敏感度和特異性結果
帶*的特征表示與CR-KECA特征相比，差異具有統計學意義(
the feature noted by * indicates that CR-KECA achieves significant improvement compared with this feature (
圖 1結果首先表明了CR可以作為特征提取方法對原始特征進行重表達得到新的特征。在基于EEG信號的睡眠分期實驗中，提取的原始特征包含了時域、頻域和時頻域的特征，但是其表征不同階段睡眠信號的性能還不夠，而基于CR特征的分類效果較原始特征的性能有顯著提高。這表明CR是一種有效的特征提取方法。其原因在于針對特定的目標信號，CR特征實際上是通過學習該目標信號內在的屬性而獲取的特征，而最初提取的78個原始特征是廣泛應用于各種生理信號的特征，并不具有CR特征的針對性，從而導致原始信號并不能很好地表征特定的睡眠狀態PSG信號，而通過學習得到的CR特征則具有很好的鑒別睡眠狀態的能力。
另一方面，由CR編碼所產生的系數向量，是基于冗余字典而實現的，對原始特征而言是一個升維的過程，所以CR特征是高維特征。而CR是將信號表示為冗余字典中所有字典原子的組合，因此，高維特征中必然還有一部分冗余或者無關的信息。因此，采用KECA對高維的CR特征進行降維，進一步提升了特征的表達能力。以應用最為廣泛的PCA為代表的傳統譜變換方法，其降維后的特征是選擇變換后數據的前n個最大特征值對應的特征；而KECA作為一種基于信息論的數據變換譜方法，其降維后的特征選擇是基于熵值，這就能最大程度地保留數據內在的結構屬性，因此，KECA降維后的特征，其分類效果也更優異。因而將CR與KECA相結合的特征提取方法應用于基于PSG的睡眠狀態分期，取得了最佳的分類效果。此外，CR-KECA特征提取算法具有較低的計算復雜度，時間效率高。一方面，SR中采用l1-范數作為約束條件，其計算復雜度非常高；而CR作為與SR相對應的數據編碼算法，采用l2-范數作為約束條件，從公式(4)可以發現，對于測試數據，只需將信號f投影到P空間下，即可得到編碼后的CR向量系數，計算非常直接，這樣就極大地降低了計算復雜度。另一方面，CR特征較原始特征而言是高維特征，不可避免地存在冗余信息，而通過KECA降維以后，不僅提高了特征表達的能力，同時使得新的特征具有更低的維度，這就降低了后續分類器的計算復雜度。上述兩方面的原因使得CR-KECA在特征表達能力上具有良好的效果和效率，適用于海量PSG數據的睡眠狀態分析。
由于KECA變換后的數據具有明顯的角度特性，在進一步的研究中，將采用基于角度相似度的分類器，以進一步提升分類效果。此外，在本文研究中，我們只是采用了單通道的EEG信號來驗證本文提出的特征提取算法，今后的工作還將研究如何將本文算法進一步擴展成針對多視圖數據的特征提取算法，以融合多個模態的PSG信號，進一步提高睡眠分期的效果。
4 結論
本文提出了一種新的特征提取方法，首先采用CR方法對原始特征進行重表達，然后再將高維的CR編碼特征利用KECA進行降維，最后將此方法應用于PSG信號特征提取。睡眠分期實驗結果表明原始特征經過CR編碼處理后，其特征表達能力顯著提高，而后續KECA降維則進一步提高了特征表達的性能，從而提升了睡眠分期的性能。本文提出的結合CR與KECA的特征提取算法還具有應用于其他生理信號分類的可行性。
引言
睡眠是人體一種重要的生理現象，人的一生約有三分之一的時間是在睡眠中度過的。在睡眠過程中，睡眠不是一種固定不變的狀態，而是經歷著幾個相對穩定的階段，臨床上把這幾個階段定義為睡眠狀態的分期。根據美國睡眠醫學學會2007年公布的睡眠分期標準，把睡眠階段分成了5期，分別為：覺醒期、睡眠1期、睡眠2期、慢波睡眠期、非快速眼動期[1]。由于睡眠質量關系到人體健康，而睡眠分期的分析又與睡眠質量密切相關，因此，睡眠分期的研究分析對于預防和診斷睡眠相關疾病具有重要的臨床應用價值。
多導睡眠圖(polysomnography，PSG)在臨床上廣泛應用于記錄睡眠狀態，然后對所獲取的信號進行分析和診斷。PSG包括腦電圖(electroencephalogram，EEG)[2-4]、心電圖[4-6]、眼動圖[4, 7]和頰肌電圖[4]等多種信號。由醫師對長時間記錄的睡眠信號進行人工分析，不僅耗時耗力，而且易于引起漏診和誤診。因此，基于計算機的睡眠數據輔助分析是目前睡眠分析的研究熱點。
PSG信號的有效特征提取是實現睡眠分期的重要環節。目前，常用的特征提取方法包括時域、頻域和時頻域特征提取方法[3, 5, 8]。然而，基于現有特征提取方法的睡眠分期的效果還有待進一步的提高。
協作表示(collaborative representation，CR)是新提出的一種數據表示方法[9]，是一種與目前在信號處理領域取得非常成功應用的稀疏表示(sparse representation，SR)相對應的方法。SR旨在采用冗余字典中的少量字典原子的線性組合來表示信號，常采用l1-范數作為稀疏約束；而CR則是采用全部樣本的組合來表示信號，采用l2-范數作為約束條件，通過l2-范數求解系數并將不同類樣本協作共同用作訓練字典。研究還表明，是CR而非l1-范數稀疏約束使得基于SR的分類取得了巨大的成功[9]。此外，CR的計算復雜度遠遠小于基于l1-范數的SR。因此，基于CR的分類取得了與基于SR的分類相當的效果[9-11]。但是，目前CR主要是作為分類器來應用，鮮有用于特征提取，更未用于PSG信號處理。
另一方面，經過CR以后獲取的新的數據向量表示，其特征維度往往很高，而且由于CR是全部樣本的組合，其編碼系數在一定程度上還有冗余或者無關信息。因此，對CR獲得的編碼系數進行降維是降低其后續處理計算復雜度、進一步提高數據表達能力的有效方法。核熵成分分析(kernel entropy component analysis，KECA)是一種新的基于信息論的數據變換方法[12]。其特質在于尋找高維特征空間中能夠最大程度保持原始空間數據二次Renyi熵值的坐標軸作為投影方向。傳統的數據變換譜方法，如主成分分析(principal component analysis，PCA)等，所選擇的是數據變換后矩陣的前n個最大特征值所對應的特征作為降維后的特征，而KECA則與之不同，是根據熵值來確定降維后的特征，這就最大程度地保證了數據內在蘊含的結構屬性。因此，KECA已經在高維數據的降維、聚類等應用中表現出了甚至優于PCA的性能[12-14]。
本文針對基于PSG的睡眠狀態分期中的特征提取問題，采用CR對特征進行再表達，將CR編碼系數作為新的特征，然后通過KECA對CR特征進行降維，獲取更為緊致和更具表征能力的特征。
1 協作表示
對于SR和CR，都可以用以下一般模型來描述：
| $\hat \alpha = \arg {\min _\alpha }\left\{ {{{\left\| {{\rm{f}} - \boldsymbol{A}\alpha } \right\|}_{{l_q}}} + \lambda {{\left\| \alpha \right\|}_{{l_q}}}} \right\}$ |
其中，f是待分析信號，A是冗余字典，α是稀疏系數，λ是正則參數，p,q=1或者2，對于p，q不同的賦值可以得到不同的具體的協作表示模型。在SR中，p=1，q=2。
雖然l1-范數約束是SR常用的方法，但是l1-范數最小化問題求解復雜，計算復雜度高。Zhang等[9]用l2-范數即最小二乘正則化取代l1-范數，提出了CR的概念，其模型如下:
| $\hat \alpha = \arg {\min _\alpha }\left\{ {\left\| {{\rm{f}} - \boldsymbol{A}\alpha } \right\|_2^2 + \lambda \left\| \alpha \right\|_2^2} \right\}$ |
最小二乘正則項有兩個作用，一是使得最小二乘解穩定，二是求解時引入一定量的稀疏。 對式(2)求解問題，很容易求得唯一解：
| $\hat \alpha = \boldsymbol{P}{\rm{f}}$ |
其中，
| $\boldsymbol{P} = {\left( {{\boldsymbol{A}^{\rm{T}}}\boldsymbol{A} + \lambda \times \boldsymbol{I}} \right)^{ - 1}}{\boldsymbol{A}^{\rm{T}}}$ |
顯然，P獨立于信號f，可以作為投影矩陣。測試時，只需將信號f投影到P空間下，即可得到向量系數，這使得編碼效率大大提高。
在Zhang等[9]提出的協作表示模型中，冗余字典A往往直接采用訓練樣本替代。由于K-means聚類算法被廣泛應用于字典訓練，在本文中，采用K-means聚類算法來訓練字典A。
2 核熵成分分析
KECA的原理如下：
給定數據集S=x1,…,xN，p(x)是其概率密度函數，則其二次Renyi熵為：
| ${\rm{H}}\left( {\rm{p}} \right) = - \log \int {{{\rm{p}}^2}\left( x \right)} {\rm{d}}x$ |
由于對數函數的單調性，因此主要集中定性分析。采用如公式(6)所示的Parzen窗密度估計器來估計V(p)的核函數：
| ${\rm{\hat p}}\left( x \right) = \frac{1}{N}\sum\limits_{{x_t} \in S} {{{\rm{K}}_\sigma }\left( {x,{x_t}} \right)} $ |
其中，Kσ(x,xt)是Parzen窗或者稱之為以xt為中心、以σ為寬度的核函數。使用(x)的均值對V(p)進行估計，可以得到：
| ${\rm{\hat V}}\left( {\rm{p}} \right) = \frac{1}{N}\sum\limits_{{x_t} \in S} {{\rm{\hat p}}\left( {{x_t}} \right)} = \frac{1}{N}\sum\limits_{{x_t} \in S} {\frac{1}{N}\sum\limits_{{x_t} \in S} {{{\rm{K}}_\sigma }\left( {{x_t},{x_{t'}}} \right) = } } \frac{1}{{{N^2}}}{\boldsymbol{I}^{\rm{T}}}\boldsymbol{K}\boldsymbol{I}$ |
其中，I為(N×1)的單位向量，K為(N×N)的核矩陣，元素(t,t′)的核矩陣可表示為Kσ(xt,xt′)。那么，Renyi熵估計的核矩陣可以用特征值和特征向量來表示，即：K=EDET，其中，D是特征值為[λ1,λ2,…,λN]的對角矩陣，E是相應的特征向量，即：E=[e1,e2,…,eN]。那么，公式(7)可以表示為：
| ${\rm{\hat V}}\left( {\rm{p}} \right) = \frac{1}{{{N^2}}}\sum\limits_{i = 1}^N {{{\left( {\sqrt {{\lambda _i}} e_i^T\boldsymbol{I}} \right)}^2}} $ |
假設k(k<N)維數據通過Φ映射到子空間Uk，當且僅當子空間與Renyi熵建立聯系時，根據熵的大小將特征值和特征向量進行重新排序，產生KECA的映射Φeca
| ${\boldsymbol{\Phi} _{eca}} = {{\rm{\boldsymbol{U}}}_k}\boldsymbol{\Phi} = \boldsymbol{D}_k^{\frac{1}{2}}\boldsymbol{E}_k^T$ |
公式(9)的求解可以轉換成下面公式的最小化問題：
| ${\boldsymbol{\Phi} _{eca}} = \boldsymbol{D}_k^{\frac{1}{2}}\boldsymbol{E}_k^T:\min {\rm{\hat V}}\left( {\rm{p}} \right) - {{\rm{\hat V}}_k}\left( {\rm{p}} \right)$ |
結合公式(8)可知：
| ${{\rm{\hat V}}_k}\left( {\rm{p}} \right) = \frac{1}{{{N^2}}}\sum\limits_{i = 1}^k {{{\left( {\sqrt {{\lambda _i}} e_i^T\boldsymbol{I}} \right)}^2}} $ |
因此，可以發現對二次Renyi熵貢獻較高的是特定的特征值和相應的特征向量，而不是前面的特征值和相應的特征向量。顯然，KECA是通過選取對二次Renyi熵貢獻較高的特征值和相應的特征向量來盡可能地保持數據變換前后熵的不變性。
3 實驗與結果分析
3.1 數據庫
為了驗證本文提出的CR結合KECA的特征提取算法的有效性，本文采用St. Vincent’s University Hospital和University College Dublin提供的睡眠呼吸暫停數據庫[15]。這個數據庫包含25名疑似睡眠呼吸紊亂的被試者的整夜PSG數據，其中(C3-A2)通道的EEG信號被選擇用于本文的測試研究。將EEG信號每30秒分割成一個片段，每個片段對應于5種睡眠狀態的一個分期階段。
對每個片段進行如下計算來提取特征，獲取共計78個特征值：
(1)計算信號的均值、標準差、中值、幅值(信號最大值與最小值絕對值差)；
(2)計算信號的峰態系數和信息熵；
(3)計算信號譜密度線性擬合的負斜率；
(4)根據10個不同頻段(δ1,δ2,θ1,θ2,α1,α2,σ1,σ2,β,γ)運用快速傅里葉變換提取特征；
(5)對(4)中求得的信號最大值和最大值所在的頻帶作為特征；
(6)用傅里葉變換在五個頻段(α,β,δ,θ和高頻)提取特征；
(7)對(4)中得到的特征向量進行頻帶劃分，計算頻帶間比率作為特征；
(8)小波變換提取信號特征：計算六個頻段(α,β,δ,θ,高γ，低γ)的小波系數，再對每個頻段的小波系數計算最大值、最小值、標準差、均值、中值、四分位數、半四分位數、熵和能量。
3.2 實驗設計
為了驗證本文提出的CR結合KECA降維所提取特征(下文用CR-KECA特征表示)的性能，與如下3種特征提取算法進行對比實驗：①EEG信號的原始特征；②原始特征經過CR編碼以后獲得的特征(CR特征)；③對CR特征采用最常用的PCA降維以后獲取的特征(CR-PCA特征)。
由于所有樣本的片段數量巨大，我們隨機挑選了8 000個樣本數據用于本實驗，每一類(睡眠分期)的樣本數量平均分布。其中，隨機選擇5 000個樣本采用K-means聚類算法訓練字典，字典原子個數為312，即CR特征的維度為312。剩余的3 000個樣本中，隨機選擇1 000個樣本作為訓練集，2 000個樣本作為測試。上述實驗過程重復5次。分類器采用LIBSVM程序包中的線性支持向量機(support vector machine，SVM)[16]。
分類評價指標選擇分類精度(accuracy，ACC)、敏感度(sensitivity，SEN)、特異性(specificity,SPE)，計算公式如下：
| $\left\{ \begin{array}{l} {\rm{ACC = }}\frac{{{\rm{TP + TN}}}}{{{\rm{TP + FN + TN + FP}}}}\\ {\rm{SEN = }}\frac{{{\rm{TP}}}}{{{\rm{TP + FN}}}}\\ {\rm{SPE = }}\frac{{{\rm{TN}}}}{{{\rm{TN + FP}}}} \end{array} \right.$ |
其中，TP、FN、TN、FP分別表示真陽性、假陰性、真陰性和假陽性。五次實驗的平均結果作為最終的性能指標(均值±標準差)，并且進行t檢驗，當P值小于0.05時，差異有統計學意義。
3.3 結果及分析
表 1所示為CR-KECA特征在不同維度下的分類結果。其中，分類性能最高的結果出現在60維，其分類精度、敏感度和特異性分別為68.74%±0.46%、68.76%±0.43%、92.19%±0.11%。從表中還可以發現，從5維特征開始到30維，隨著維數的增加，準確率、敏感度和特異性均呈顯著提高，但是30維以后的特征維度下，結果相差不大。
表1
不同維度下的CR-KECA特征的分類精度、敏感度及特異性(%)
Table1.
Classification accuracy,sensitivity and specificity of CR-KECA feature with different dimensions (%)
為了直觀比較CR-KECA特征與其他三種特征的性能，圖 1依次顯示了原始特征、CR特征、CR-PCA特征和CR-KECA特征的平均分類精度、敏感度和特異性。由圖可以發現，CR-KECA的分類精度較原始特征、CR特征、CR-PCA特征分別提高了2.92%、1.47%、0.96%，敏感度分別提高了2.94%、1.49%、0.98%，特異性分別提高了0.73%、0.37%、0.24%。此外，CR-KECA特征與原始特征和CR特征之間在分類精度、敏感度和特異性這三個指標上的差異都具有統計學意義(P<0.05)，表明CR-KECA具有顯著的性能提高。另一方面，與原始的特征相比，CR特征的分類精度、敏感性和特異性較原始特征也有較大的提高，分別為1.45%、1.46%、0.36%。
圖1
四種特征分類的精度、敏感度和特異性結果
帶*的特征表示與CR-KECA特征相比，差異具有統計學意義(
the feature noted by * indicates that CR-KECA achieves significant improvement compared with this feature (
圖 1結果首先表明了CR可以作為特征提取方法對原始特征進行重表達得到新的特征。在基于EEG信號的睡眠分期實驗中，提取的原始特征包含了時域、頻域和時頻域的特征，但是其表征不同階段睡眠信號的性能還不夠，而基于CR特征的分類效果較原始特征的性能有顯著提高。這表明CR是一種有效的特征提取方法。其原因在于針對特定的目標信號，CR特征實際上是通過學習該目標信號內在的屬性而獲取的特征，而最初提取的78個原始特征是廣泛應用于各種生理信號的特征，并不具有CR特征的針對性，從而導致原始信號并不能很好地表征特定的睡眠狀態PSG信號，而通過學習得到的CR特征則具有很好的鑒別睡眠狀態的能力。
另一方面，由CR編碼所產生的系數向量，是基于冗余字典而實現的，對原始特征而言是一個升維的過程，所以CR特征是高維特征。而CR是將信號表示為冗余字典中所有字典原子的組合，因此，高維特征中必然還有一部分冗余或者無關的信息。因此，采用KECA對高維的CR特征進行降維，進一步提升了特征的表達能力。以應用最為廣泛的PCA為代表的傳統譜變換方法，其降維后的特征是選擇變換后數據的前n個最大特征值對應的特征；而KECA作為一種基于信息論的數據變換譜方法，其降維后的特征選擇是基于熵值，這就能最大程度地保留數據內在的結構屬性，因此，KECA降維后的特征，其分類效果也更優異。因而將CR與KECA相結合的特征提取方法應用于基于PSG的睡眠狀態分期，取得了最佳的分類效果。此外，CR-KECA特征提取算法具有較低的計算復雜度，時間效率高。一方面，SR中采用l1-范數作為約束條件，其計算復雜度非常高；而CR作為與SR相對應的數據編碼算法，采用l2-范數作為約束條件，從公式(4)可以發現，對于測試數據，只需將信號f投影到P空間下，即可得到編碼后的CR向量系數，計算非常直接，這樣就極大地降低了計算復雜度。另一方面，CR特征較原始特征而言是高維特征，不可避免地存在冗余信息，而通過KECA降維以后，不僅提高了特征表達的能力，同時使得新的特征具有更低的維度，這就降低了后續分類器的計算復雜度。上述兩方面的原因使得CR-KECA在特征表達能力上具有良好的效果和效率，適用于海量PSG數據的睡眠狀態分析。
由于KECA變換后的數據具有明顯的角度特性，在進一步的研究中，將采用基于角度相似度的分類器，以進一步提升分類效果。此外，在本文研究中，我們只是采用了單通道的EEG信號來驗證本文提出的特征提取算法，今后的工作還將研究如何將本文算法進一步擴展成針對多視圖數據的特征提取算法，以融合多個模態的PSG信號，進一步提高睡眠分期的效果。
4 結論
本文提出了一種新的特征提取方法，首先采用CR方法對原始特征進行重表達，然后再將高維的CR編碼特征利用KECA進行降維，最后將此方法應用于PSG信號特征提取。睡眠分期實驗結果表明原始特征經過CR編碼處理后，其特征表達能力顯著提高，而后續KECA降維則進一步提高了特征表達的性能，從而提升了睡眠分期的性能。本文提出的結合CR與KECA的特征提取算法還具有應用于其他生理信號分類的可行性。
