董孝彤,曲新亮,魏守水△
(1. 山東大學(xué)控制科學(xué)與工程學(xué)院,濟南 250061;2.山東第一醫(yī)科大學(xué)第一附屬醫(yī)院,濟南 250014)
睡眠呼吸暫停會導(dǎo)致患者白天嗜睡及生活萎靡不振,是高血壓,心血管疾病的高風(fēng)險因素[1-2]。目前,對睡眠呼吸障礙的檢測主要是依據(jù)多導(dǎo)睡眠儀(PSG),但其價格昂貴且需醫(yī)院專業(yè)技師操作,無法推廣,因此,尋求價格低廉且易于推廣的方法一直是醫(yī)學(xué)工程界的共識[3]。睡眠呼吸暫停是指口鼻氣流量低于正常的90%,且持續(xù)時間超過10 s,低通氣是指氣流量低于正常的30%且血氧飽和度下降4%[4]。研究表明,睡眠呼吸暫停能引起心率的周期性變化,發(fā)生呼吸暫停時心動過緩,呼吸暫停結(jié)束時,心動過速且血壓升高。呼吸暫停發(fā)作時,由于缺氧及高碳酸血癥,交感神經(jīng)活動逐漸增加,在呼吸暫停結(jié)束時,交感神經(jīng)活動最為明顯,隨著呼吸暫停和呼吸恢復(fù),交感神經(jīng)活動強度迅速減弱。與健康人相比,即使在白天,睡眠呼吸暫?;颊叩纳窠?jīng)活動水平也很高,而且晝夜心率變異性(HRV)明顯減弱[5]。
因此,通過心電信號分析,實現(xiàn)睡眠呼吸暫停事件的檢測一直是熱門的研究課題。Penzel等[6]在PhysioNet發(fā)起的基于心電圖的呼吸暫停檢測和量化的挑戰(zhàn)賽總結(jié)中指出,表現(xiàn)最好的算法均使用了HRV的頻域特征或EDR(Electrocardiogram-derived Respiration)信號提取的特征。Correa等[7]對3種方法得到的EDR信號進(jìn)行頻域分析,使用閾值法進(jìn)行分類,實現(xiàn)Se=88.28%,Sp=88.66%。Bsoul等[8]對RR間期及EDR信號進(jìn)行時、頻域分析,使用徑向基核支持向量機(RBF-SVM)建立分類模型并根據(jù)ROC和F1的表現(xiàn)進(jìn)行特征選擇,實現(xiàn)Acc=86.24%。Manish等[9]使用雙正交反對稱小波濾波器組將心電信號片段分解為6個小波子帶,每個子帶信號計算模糊熵和對數(shù)能量特征,使用Mann -Whitney U檢驗進(jìn)行特征選擇,使用35折交叉驗證及最小二乘支持向量機(LS-SVM)實現(xiàn)Acc=90.11%,Se=90.87%,Sp=88.88%。梁九興等[10]使用HRV時頻特征及概率神經(jīng)網(wǎng)絡(luò)(PNN)構(gòu)建分類模型,實現(xiàn)Acc=75.97%,Se=82.51%,Sp=76.22%。Sadr等[11]對HRV與EDR信號進(jìn)行時頻分析并計算EDR信號的心肺耦合頻譜特征,實現(xiàn)Acc=86.5%。然而,很多研究致力于提高檢測的準(zhǔn)確率卻忽視了特征的穩(wěn)定性研究,本研究針對睡眠障礙檢測中的一些常用心電特征進(jìn)行穩(wěn)定性分析,并結(jié)合SVM建立睡眠呼吸暫停檢測模型。
本研究使用的數(shù)據(jù)庫為Physionet的Apnea-ECG數(shù)據(jù)庫。它包含70例睡眠記錄,其中,心電信號為V2導(dǎo)聯(lián),采樣頻率為100 Hz,時長約8 h。信號的每分鐘均由專家標(biāo)記是否包含睡眠呼吸暫停和低通氣事件,若存在,標(biāo)記為“A”,不存在,標(biāo)記為“N”。由于僅訓(xùn)練集的標(biāo)注公開,本研究使用35例訓(xùn)練集記錄,其中30例用作訓(xùn)練,5例用作測試。基于1 min時長實現(xiàn)睡眠呼吸暫停事件的檢測。
基于心電信號提取RR間期信號及EDR信號。使用HRV分析工具箱[12]檢測QRS波的位置,將相鄰QRS波位置差分,得到RR間期序列,并對RR間期進(jìn)行校正[13]。為了保留足夠的信息量,在分析中排除了校正后累計RR間期時長小于30 s的信號片段。
呼吸過程中,胸腔內(nèi)的電阻抗會發(fā)生變化,當(dāng)肺部的空氣量增加或減少時,電極相對于心臟的位置也會發(fā)生變化,因此呼吸運動的信息會反應(yīng)在ECG中。本研究使用QRS的面積估計EDR信號[7],圖1為真實呼吸信號與EDR信號。
圖1 真實呼吸信號與EDR信號
對RR間期序列及EDR信號進(jìn)行時域、頻域、非線性及相關(guān)性分析,共提取了45個特征,見表1。頻域分析中,使用三次樣條插值將RR間期和EDR信號重采樣到4 Hz,對重采樣后的RR間期與EDR信號進(jìn)行9層DB4小波變換,計算第4-9層近似系數(shù)(頻域范圍約為0.007~0.5 Hz)的方差(var_RR4- var_RR9、var_EDR4- var_EDR9)[8]。使用快速傅里葉變換(FFT)對重采樣后的RR間期與EDR信號估計功率譜密度(PSD),并分別計算極低頻功率VLF(0~0.04 Hz)、低頻功率LF(0.04~0.15 Hz)、高頻功率HF(0.15~0.4 Hz)及LF/HF。除上述特征外,EDR信號還計算了中心頻率(Fc)及主頻率(Fm)[7]兩個特征。計算公式如下,其中,fi表示頻率,N表示PSD的長度。
(1)
(2)
非線性分析中,RR間期的龐加萊圖可以反映RR間期的分布情況。其中,SD1為y=-x+2μ(μ為均值)方向上的方差,與高頻功率有關(guān)。SD2為y=x方向上的方差,與低頻及高頻功率均相關(guān)。兩者的比值SD1/SD2反映了長、短期HRV之間的平衡[14]。近似熵(sample entropy,ApEn)[15]是對不穩(wěn)定時間序列復(fù)雜度的一種度量,其思想是檢測一個時間序列中產(chǎn)生新子序列的概率。模糊測度熵(fuzzy measure entropy,F(xiàn)MEn)[16]使用模糊函數(shù)的隸屬度代替ApEn中使用的Heaviside函數(shù)作為向量相似判斷標(biāo)準(zhǔn)。同時,F(xiàn)MEn利用模糊局部測度熵和模糊全局測度熵反映生理信號中隱含的整體復(fù)雜性,彌補了模糊熵只關(guān)注局部復(fù)雜性的弱點。
相關(guān)性分析中,Bianchi等[17]在研究中發(fā)現(xiàn)RR間期序列及EDR信號的頻譜相關(guān)系數(shù)在睡眠呼吸暫停事件的識別中表現(xiàn)良好,因此,本研究計算了RR間期序列及EDR信號在極低頻段、低頻段及高頻段的相關(guān)系數(shù),分別用REv、REl及REf表示。除此之外,計算了RR間期序列及其延遲(1-5個樣本)序列的相關(guān)系數(shù)(RRcorr1- RRcorr5)[8]。
表1 基于RR間期及EDR信號的特征
特征選擇可以加快數(shù)據(jù)挖掘算法的速度,提高學(xué)習(xí)精度,增強模型的可理解性,但特征選擇結(jié)果往往因樣本差異或噪聲等因素而出現(xiàn)不一致。為了獲取穩(wěn)定的特征排序,本研究使用集成特征選擇方法進(jìn)行特征選擇[18]。將訓(xùn)練集隨機分成5組,分別使用最小冗余最大相關(guān)(mRMR)算法[19]進(jìn)行排序,使用穩(wěn)健排序聚合(RRA)[20]方法得到最終排序結(jié)果,我們將此方法稱為mRMR-RRA,整體結(jié)構(gòu)見圖2。其中,RRA通過計算每個特征在各次排名間的差異顯著性對特征排序,其復(fù)雜度與輸入維度大小成線性關(guān)系,對異常值、噪聲和誤差具有較強的魯棒性。mRMR是一種基于相關(guān)關(guān)系的過濾式特征選擇方法,它不僅考慮了特征與標(biāo)簽之間的相關(guān)性,而且考慮了特征之間的信息冗余。其中,相關(guān)性用互信息(MI)度量,兩個離散隨機變量X和Y的MI可以定義為:
圖2 mRMR-RRA方法的整體結(jié)構(gòu)
Fig.2The overall structure of the mRMR-RRA method
(3)
其中p(x,y)是x和y的聯(lián)合概率密度函數(shù),p(x)和p(y)分別是x和y的邊緣概率密度函數(shù)。MI越大,相關(guān)性越高。設(shè)S是所有特征的集合,Sm-1是已排序特征集合,根據(jù)最小冗余最大相關(guān)準(zhǔn)則,添加到Sm-1中的第m個特征滿足以下條件:
(4)
其中,x表示特征,c表示樣本標(biāo)簽。
使用斯皮爾曼排序相關(guān)系數(shù)(Ds)[21]對排序結(jié)果的穩(wěn)定性進(jìn)行評價。設(shè)r和r′分別為同一方法在不同數(shù)據(jù)上得到的排序向量,r(i)和r′(i)是特征i在排序向量r和r′中的位置,c為特征總數(shù),計算方法如下:
(5)
Ds的值在[-1,1]之間,Ds越接近1,排序結(jié)果越穩(wěn)定。對于n次排序結(jié)果,計算Ds的平均值:
(6)
使用RBF-SVM建立分類模型,通過10折交叉驗證選擇最佳特征數(shù)量,平臺為matlab2017b。為了加快訓(xùn)練速度,對特征進(jìn)行歸一化處理,為了提高分類性能,對分類結(jié)果進(jìn)行中值濾波后處理[22]。使用準(zhǔn)確率(Acc)、靈敏性(Se)、特異性(Sp)作為模型的評價指標(biāo)。
為了對比mRMR-RRA方法與mRMR方法的穩(wěn)定性,在訓(xùn)練集中隨機抽取數(shù)據(jù),每次抽取數(shù)據(jù)量為總數(shù)的80%,共抽取5次,計算每兩次間的Ds值和平均Ds值,見表2。從表中可以看出,使用mRMR-RRA方法,特征排序穩(wěn)定性明顯提高。
表2 mRMR- RRA方法和mRMR方法的排序穩(wěn)定性
圖3列出了各個特征在呼吸暫停信號與正常信號間的分布差異及最終排序結(jié)果。每個子圖中左側(cè)為正常睡眠信號,右側(cè)為呼吸暫停信號?!甊’表示特征來自RR間期,‘E’表示特征來自EDR信號。根據(jù)mRMR-RRA排序結(jié)果,特征從左到右,從上到下依次排列。排名前10的特征包括:RR間期的頻域分析(RRd6,RRd7,LF/HF),EDR信號的時域分析(var)、頻域分析(Fm)及非線性分析(FMEn),相關(guān)性分析中的RRcorr2,RRcorr3,RRcorr4,REh??梢钥闯觯判蚩壳暗奶卣髟趦深愰g的差異明顯大于排序靠后的特征,說明差異性大的特征表現(xiàn)較穩(wěn)定。雖然有些排序靠后的特征在兩類間的差異也很顯著,但由于其與排名靠前的特征相關(guān)性較高,因此排名靠后,如RRcorr2、RRcorr1等。
圖3正常睡眠信號與呼吸暫停信號間特征的分布差異(歸一化結(jié)果)
Fig.3Differences in the distribution of characteristics between normal sleep signals and apnea signals (Normalized results)
RR間期的頻域分析顯示,呼吸暫停信號的VLF,LF,LF/HF與正常睡眠信號相比偏高,這與睡眠呼吸暫停發(fā)生時,交感神經(jīng)活性增強的生理學(xué)結(jié)果吻合。觀察RRd9-RRd4的變化可以發(fā)現(xiàn),睡眠呼吸暫停與正常睡眠信號間的差異在低頻時更明顯。EDR的頻域分析與RR間期的相應(yīng)特征分布相似,說明EDR信號與RR間期信號有較好的相關(guān)性。睡眠呼吸暫停信號的Fc和Fm相對正常呼吸信號降低。非線性分析中,RR間期與EDR信號的熵值均顯示與ApEn相比,F(xiàn)MEn在兩類間的差異更顯著,并且由于交感神經(jīng)興奮會增加信號的確定度[23],因此呼吸暫停信號的熵值顯著降低。相關(guān)性分析中,RR間期及其延遲序列相關(guān)性在呼吸暫停信號中顯著大于正常睡眠信號,RR間期與EDR在高頻段的相關(guān)性更能反映兩類間的差異。
圖4為兩種方法10折交叉驗證準(zhǔn)確率隨特征數(shù)量的變化情況??梢钥闯觯琺RMR-RRA方法整體優(yōu)于mRMR方法,說明特征的穩(wěn)定性增加的同時,分類效果也得到了提升。當(dāng)特征數(shù)量為14時,mRMR-RRA方法達(dá)到局部最優(yōu),此時Acc=83.22%。隨著特征數(shù)量繼續(xù)增加,Acc呈小幅度波浪式變化。雖然特征數(shù)量大于39時,準(zhǔn)確率提高到了86%,但模型的復(fù)雜度過大。因此,在最終分類模型的訓(xùn)練中,使用了mRMR-RRA方法得到的前14個特征。
圖4 10折交叉驗證準(zhǔn)確率隨特征數(shù)量的變化
表3 分類結(jié)果
表3為分類模型在測試集上的表現(xiàn)。第一行顯示Se=71.30%,遠(yuǎn)遠(yuǎn)低于Sp=96.95%,說明假陰性(FN)較高,這在臨床中存在非常大的隱患。因此,在訓(xùn)練模型時重新設(shè)置代價矩陣,將FN的誤分類代價(Cost)設(shè)置為假陽性(FP)的兩倍,此時得到Se=86.71%,相比之前提升了15.41%,Acc=90.03%,相比之前提升了1.77%。
本研究針對心電信號檢測睡眠呼吸暫停事件中心電特征的穩(wěn)定性進(jìn)行研究。實驗結(jié)果表明,基于集成穩(wěn)定特征選擇策略,將mRMR特征選擇方法與RRA方法結(jié)合進(jìn)行特征穩(wěn)定性分析,與單獨使用mRMR方法相比,特征穩(wěn)定性及分類準(zhǔn)確率均有所提高。ECG指標(biāo)與EDR指標(biāo)有較好的相關(guān)性,結(jié)合ECG信號重構(gòu)呼吸信號值得進(jìn)一步研究。本研究給出的特征排序方法與結(jié)果,可以為心電信號識別睡眠障礙等相關(guān)研究提供一定的參考意義。