基于粒子群優(yōu)化-支持向量機(jī)的睡眠呼吸暫停檢測(cè)

張大可， 馬雋， 王立英， 王鋼*， 蔡靖， 孫玉冰

(1.北華大學(xué)電氣與信息工程學(xué)院， 吉林 132021； 2. 北華大學(xué)基礎(chǔ)醫(yī)學(xué)院， 吉林 132021； 3. 清華大學(xué)附屬垂楊柳醫(yī)院， 北京 100020； 4.吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院， 長春 130061)

睡眠呼吸暫停(sleep apnea，SA)是指在連續(xù)7 h睡眠中發(fā)生30次以上的呼吸暫停，每次氣流中止10 s以上(含10 s)，或平均每小時(shí)低通氣次數(shù)(呼吸紊亂指數(shù))超過5次，而引起慢性低氧血癥及高碳酸血癥的臨床綜合征，是一種睡眠時(shí)呼吸停止的睡眠障礙[1-2]，主要表現(xiàn)為睡眠過程中反復(fù)發(fā)生上氣道塌陷、阻塞，睡眠打鼾和白天嗜睡，常見于中年以上肥胖男性，對(duì)患者的心血管系統(tǒng)和內(nèi)分泌系統(tǒng)等有著嚴(yán)重的損害，增大腫瘤患病的可能性[3]。

目前關(guān)于睡眠呼吸暫停的診斷主要是利用多導(dǎo)睡眠儀(polysomnography，PSG)[4]對(duì)患者的腦電、心電(electrocardiogram, ECG)、血氧飽和度或眼電等生理信號(hào)進(jìn)行監(jiān)測(cè)，然而PSG價(jià)格昂貴，且佩戴電極繁雜，佩戴過程煩瑣，影響患者正常睡眠[5]，監(jiān)測(cè)過程費(fèi)時(shí)費(fèi)力，往往需要數(shù)位醫(yī)生共同分析監(jiān)測(cè)結(jié)果，無法給予患者及時(shí)有效的診斷反饋。

Guilleminault等[6]最早利用單導(dǎo)ECG信號(hào)對(duì)睡眠呼吸暫停實(shí)現(xiàn)了診斷。近年來，研究者們發(fā)現(xiàn)，睡眠呼吸暫停患者在呼吸暫停發(fā)作時(shí)，會(huì)產(chǎn)生明顯的心率周期性變化現(xiàn)象。梁九興等[7]利用心率變異性(heart rate variability, HRV)的時(shí)、頻域特征構(gòu)建概率神經(jīng)網(wǎng)絡(luò)(probabilistic neural network，PNN)模型，用于檢測(cè)SA，使用10折交叉驗(yàn)證實(shí)現(xiàn)準(zhǔn)確率達(dá)75.97%。Cafer等[8]通過心電信號(hào)估算心電伴生呼吸(derived respiration，EDR)信號(hào)，進(jìn)而利用EDR信號(hào)估算受試者瞬時(shí)心率序列，對(duì)3 min信號(hào)對(duì)應(yīng)的瞬時(shí)心率序列進(jìn)行小波分解，再利用非線性自回歸型人工神經(jīng)網(wǎng)絡(luò)分類器構(gòu)建分類模型，準(zhǔn)確率可達(dá)78.3%。范文兵等[9]通過多尺度熵證明ECG信號(hào)多尺度熵指數(shù)能夠區(qū)分SA。葛靖等[10]采用長短時(shí)記憶模型，通過心電的間接信號(hào)提取特征，識(shí)別睡眠呼吸暫停，準(zhǔn)確率可達(dá)87.4%。然而，利用上述方法進(jìn)行睡眠呼吸暫停檢驗(yàn)或分類的準(zhǔn)確率不高，導(dǎo)致誤檢的可能性增加。

現(xiàn)采用粒子群優(yōu)化-支持向量機(jī)(particle swarm optimization-support vector machine，PSO-SVM)方法[11]，分析處理受試者的心電信號(hào)，實(shí)現(xiàn)對(duì)睡眠呼吸暫停的準(zhǔn)確檢測(cè)，以期不僅減輕醫(yī)務(wù)工作者的負(fù)擔(dān)，也及時(shí)分析患者的檢測(cè)結(jié)果并得出結(jié)論，提高工作效率。

1 方法

1.1 總體框架

提出的睡眠呼吸暫停檢測(cè)方法的主要計(jì)算流程如圖1所示，首先讀取數(shù)據(jù)文件，根據(jù)注釋文件對(duì)睡眠呼吸暫停狀態(tài)進(jìn)行分類，再根據(jù)期望數(shù)據(jù)要求對(duì)心電信號(hào)進(jìn)行預(yù)處理，去除噪聲和因其他原因產(chǎn)生的雜亂信號(hào)，進(jìn)而利用峰值檢波提取心電信號(hào)的R波，并計(jì)算心率和心率變異性(heart rate variablity，HRV)，對(duì)上述信號(hào)進(jìn)行特征提取，再利用方差閾值法對(duì)特征進(jìn)行篩選，最后利用粒子群算法優(yōu)化支持向量機(jī)，從而實(shí)現(xiàn)對(duì)睡眠呼吸暫停的檢測(cè)。

圖1 總體框圖Fig.1 Overall block diagram

1.2 數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理主要包含兩方面：對(duì)心電信號(hào)進(jìn)行切割分段；對(duì)心電信號(hào)進(jìn)行心率變異性的計(jì)算。

首先對(duì)ECG按照60 s時(shí)間窗分段，對(duì)應(yīng)睡眠呼吸暫停或正常睡眠標(biāo)簽；在分段后，利用峰值檢測(cè)尋找R波位置，并記錄R波產(chǎn)生時(shí)間。將R波產(chǎn)生時(shí)間相減，得到受試者RR間期信號(hào)，進(jìn)而得到受試者的心率變化，再將相減后的R波時(shí)間做差分，即可得到心電信號(hào)的心率變異性HRV序列[12]。

1.3 特征提取

根據(jù)心電信號(hào)和心率變異性的定義和特點(diǎn)，以及睡眠呼吸暫停產(chǎn)生過程中患者心率波動(dòng)，對(duì)心率變異性的時(shí)域、頻域、時(shí)頻域和非線性特征[13]進(jìn)行了提取，時(shí)域特征主要包含均值(MEAN)、總體標(biāo)準(zhǔn)差(SDNN)、均值標(biāo)準(zhǔn)差(SDANN)和差值均方的平方(r-MSSSD)等[14]，頻域特征主要包含HRV頻譜特征、功率譜密度(PSD)、信號(hào)總功率(TP，f≤0.4 Hz)、甚低頻段功率(VLF，0.003 3 Hz≤f≤0.04 Hz)、低頻段功率(LF，0.04 Hz≤f≤0.15 Hz)、高頻段功率(HF，0.15 Hz≤f≤0.4 Hz)和低頻段與高頻段功率比(LF/HF)[15]，時(shí)頻域特征主要包含瞬時(shí)能量和瞬時(shí)中位頻率[16]，非線性特征主要包含LZ復(fù)雜度、近似熵(ApEn)、樣本熵(SampEn)、信息熵(InEn)、邊際譜熵(MSEn)和能量譜熵(ESEn)[17]。部分特征參數(shù)的計(jì)算公式如下所示。

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式(7)中：pm為第m個(gè)邊際譜在整個(gè)邊際譜中所占的百分比，即

(8)

(9)

(10)

式中：hx為邊際譜；H(m,t)為Hilbert譜；Pk為第k個(gè)能量譜在整個(gè)能量譜中所占的百分比，即

(11)

(12)

式中：ex(k)為能量譜；H2(k,t)為能量譜密度。

1.4 特征選擇

在進(jìn)行特征提取后，由于上述特征間存在冗余，且存在與睡眠呼吸暫停相關(guān)性小的特征，因此采用方差閾值法，計(jì)算各個(gè)特征的方差，刪除方差低于預(yù)設(shè)閾值的特征。

方差閾值法(variance thresholding, VT)是一種利用方差實(shí)現(xiàn)特征提取的方法，通過程序設(shè)置允許通過的方差最小值，即可對(duì)方差小于該值的特征進(jìn)行去除，從而篩選與檢測(cè)對(duì)象相關(guān)性大的特征，提高檢測(cè)的準(zhǔn)確率。方差表示數(shù)據(jù)分布的可變性，公式為

(13)

若一組數(shù)據(jù)的方差為0，則表示該組數(shù)據(jù)在預(yù)測(cè)過程中無意義，只能增加模型的復(fù)雜性，而不會(huì)增加模型的預(yù)測(cè)能力，方差閾值法就是將無法提高模型預(yù)測(cè)能力的特征篩除的過程。

然而，由于不同特征的分布不同，導(dǎo)致計(jì)算所得方差的尺度也不盡相同，無法直接利用方差判別其有效性，因此首先對(duì)特征進(jìn)行歸一化處理，即每一個(gè)特征feati除以特征的均值，得到歸一化特征z_feat，再求取其方差，實(shí)現(xiàn)特征選擇。

(14)

(15)

1.5 SA檢測(cè)模型建立

由于ECG信號(hào)的非線性特征，以及SA檢測(cè)標(biāo)簽僅分為“呼吸正常”和“呼吸紊亂”兩類，在特征選擇后，利用支持向量機(jī)(SVM)對(duì)睡眠呼吸暫停進(jìn)行訓(xùn)練和檢測(cè)。SVM是一種基于間隔最大化原則實(shí)現(xiàn)二分類的機(jī)器學(xué)習(xí)方法[18-19]，其目標(biāo)函數(shù)和約束條件可表示為

(16)

式(16)中：L(a)為拉格朗日函數(shù)；ai為拉格朗日系數(shù)；xi為數(shù)據(jù)向量；yi為xi對(duì)應(yīng)的數(shù)據(jù)類別；K(xixj)為核函數(shù)。基于徑向基核函數(shù)的非線性SVM具有較好的性能，因此擬選擇徑向基核函數(shù)；C為懲罰因子，C越大則懲罰越重，且泛化能力會(huì)隨之降低；高斯核超函數(shù)g表示對(duì)徑向范圍的控制。

由于懲罰因子C和高斯核超函數(shù)g會(huì)嚴(yán)重影響SVM的分類能力，采用粒子群算法(PSO)實(shí)現(xiàn)參數(shù)尋優(yōu)，以提高SVM分類器的泛化性能。粒子群算法是一種基于迭代的優(yōu)化算法，結(jié)構(gòu)簡單，具有全局尋優(yōu)能力，能夠以高概率和高收斂速度找到最優(yōu)解。PSO算法中根據(jù)式(17)和式(18)更新粒子的速度和位置[19-20]。

(17)

(18)

利用粒子群算法優(yōu)化支持向量機(jī)的程序流程圖如圖2所示。

圖2 PSO-SVM程序流程圖Fig.2 PSO-SVM program flow chart

2 結(jié)果與分析

2.1 數(shù)據(jù)源

所采用的數(shù)據(jù)來源于2000年心臟病學(xué)計(jì)算會(huì)議和PhysioNet聯(lián)合進(jìn)行的一項(xiàng)競(jìng)賽中首次公開的Physionet Apnea-ECG數(shù)據(jù)庫[21]。數(shù)據(jù)庫中共包含70組采樣頻率為100 Hz的心電圖，總持續(xù)時(shí)間為401～578 min，正常呼吸時(shí)間在11～535 min變化，呼吸紊亂時(shí)間在0～534 min變化。A類為呼吸暫停組，即患者組，定義為呼吸紊亂超過100 min的受試者，共包含16名29～63歲受試者的40條記錄，其中15名男性，1名女性；C類為對(duì)照組，即健康組，定義為呼吸紊亂少于5 min的受試者，共包含11名年齡為27～42歲的受試者的20條記錄，其中6名男性，5名女性。

在數(shù)據(jù)庫提供的測(cè)試集文件中，列舉了20條患者組和10條健康組的ECG信號(hào)和標(biāo)簽，ECG信號(hào)存儲(chǔ)于.dat文件，ECG信號(hào)的標(biāo)簽存儲(chǔ)于.apn文件，即記錄每分鐘是否存在睡眠呼吸暫停事件，若存在，則標(biāo)記為A；否則標(biāo)記為N。

2.2 數(shù)據(jù)預(yù)處理

部分原始心電信號(hào)如圖3所示，由于數(shù)據(jù)庫定義每分鐘心電信號(hào)存在一個(gè)標(biāo)簽，對(duì)心電信號(hào)按照60 s的時(shí)間窗分段，每段數(shù)據(jù)包含60×100個(gè)數(shù)據(jù)點(diǎn)，對(duì)應(yīng)睡眠呼吸暫停(A)或正常睡眠(N)標(biāo)簽。取每位受試者測(cè)試時(shí)間前423 min即423×60×100=253 800個(gè)采樣點(diǎn)的測(cè)試數(shù)據(jù)，按60 s分為423段，分別對(duì)應(yīng)每分鐘的SA標(biāo)簽，取14位受試者共計(jì)5 922組樣本數(shù)據(jù)。

圖3 原始心電信號(hào)Fig.3 Original ECG signal

圖4為利用峰值檢波檢測(cè)R波位置結(jié)果；圖5為利用RR間期計(jì)算受試者1 min心率變化；圖6為利用RR間期計(jì)算受試者1 min心率變異性波動(dòng)；圖7為心率變異性頻譜圖。

圖4 心電信號(hào)R波定位Fig.4 R wave location of ECG

圖5 受試者1 min內(nèi)心率波動(dòng)Fig.5 1 min heart rate fluctuation of subjects

圖6 受試者1 min內(nèi)心率變異性波動(dòng)Fig.6 1 min HRV fluctuation of subjects

圖7 HRV頻譜圖Fig.7 HRV spectrum

2.3 特征提取

圖8為心電信號(hào)RR間期的均值、標(biāo)準(zhǔn)差、均值標(biāo)準(zhǔn)差和標(biāo)準(zhǔn)差均值的變化情況；圖9為某受試者HRV功率譜密度；圖10為HRV信號(hào)總功率、低頻段功率、高頻段功率和功率比的變化情況；圖11為HRV邊際譜熵和能量譜熵的變化情況。可以看出，受試者處于不同睡眠呼吸階段，其心電信號(hào)各項(xiàng)特征存在明顯的波動(dòng)變化，可用于進(jìn)行睡眠呼吸暫停的標(biāo)記。

圖8 時(shí)域特征Fig.8 Time domain features

圖9 HRV功率譜密度對(duì)數(shù)Fig.9 Logarithm of HRV power spectral density

圖10 頻域特征Fig.10 Frequency domain characteristics

圖11 非線性特征Fig.11 Nonlinear characteristics

2.4 特征選擇

通過方差閾值法，按照各特征的歸一化方差排序，共篩選了10個(gè)特征放入檢測(cè)算法中，其中包含4個(gè)時(shí)域特征(均值、標(biāo)準(zhǔn)差、均值標(biāo)準(zhǔn)差、差值均方的平方)、3個(gè)頻域特征(信號(hào)總功率、低頻段功率、高頻段功率)、1個(gè)時(shí)頻域特征(瞬時(shí)中位頻率)和2個(gè)非線性特征(邊際譜熵、能量譜熵)。

2.5 模型評(píng)價(jià)

利用粒子群算法對(duì)支持向量機(jī)的參數(shù)進(jìn)行選取，進(jìn)而實(shí)現(xiàn)睡眠呼吸暫停的檢測(cè)，采用準(zhǔn)確率(Acc)衡量SVM檢測(cè)表現(xiàn)[22-23]，公式為

(18)

式(18)中：TP表示真正結(jié)果數(shù)；TN表示真負(fù)結(jié)果數(shù)；FP表示假正結(jié)果數(shù)；FN表示假負(fù)結(jié)果數(shù)。

圖12為粒子群算法適應(yīng)度曲線，經(jīng)過300次迭代，針對(duì)SVM的參數(shù)優(yōu)化的最優(yōu)解的準(zhǔn)確率可達(dá)97.89%，優(yōu)化后的懲罰因子C=100，高斯核超參數(shù)g=0.01。

圖12 PSO適應(yīng)度曲線Fig.12 PSO fitness curve

圖13為利用PSO-SVM實(shí)現(xiàn)睡眠呼吸暫停的部分結(jié)果與其真實(shí)值的對(duì)比情況，測(cè)試結(jié)果良好。

圖13 部分訓(xùn)練結(jié)果對(duì)比Fig.13 Comparison of some training results

PSO-SVM與SVM相比，實(shí)現(xiàn)了SVM最優(yōu)參數(shù)的選取，減小了算法的時(shí)間復(fù)雜度[24]，假設(shè)支持向量個(gè)數(shù)為N，輸入向量維度為dl，訓(xùn)練樣本點(diǎn)個(gè)數(shù)為1，則基于徑向基核函數(shù)SVM的計(jì)算復(fù)雜度為O(N3+N2l+Ndll+O(dl)N)；假設(shè)存在M個(gè)N維粒子，進(jìn)行T次迭代，則PSO-SVM的時(shí)間復(fù)雜度為O(TMN+dll2+O(dl)N)。圖14繪制了利用PSO-SVM和SVM檢測(cè)SA的接受者操作特征(receiver operating characteristics，ROC)曲線，展示了PSO-SVM的良好準(zhǔn)確率，曲線下的面積(area under curve，AUC)為96.94%。表1中列舉了不同方式SVM算法的執(zhí)行時(shí)間，以本文訓(xùn)練數(shù)據(jù)規(guī)模為基礎(chǔ)，PSO-SVM可以滿足對(duì)準(zhǔn)確率和執(zhí)行效率的要求。

圖14 PSO-SVM與SVM的ROC曲線對(duì)比圖Fig.14 ROC curve comparison between PSO-SVM and SVM

表1為單獨(dú)利用SVM不同核函數(shù)檢測(cè)SA和利用PSO-SVM檢測(cè)SA的結(jié)果參數(shù)，可以看出，利用粒子群算法對(duì)支持向量機(jī)進(jìn)行優(yōu)化，提高了支持向量機(jī)檢測(cè)的準(zhǔn)確率，可達(dá)94.0%。

表1 訓(xùn)練結(jié)果對(duì)比Table 1 Comparison of training results

3 結(jié)論

利用粒子群算法優(yōu)化支持向量機(jī)對(duì)睡眠呼吸暫停進(jìn)行檢測(cè)，通過單一的心電信號(hào)，降低了多種信號(hào)采集處理的復(fù)雜性和困難程度；利用方差閾值法實(shí)現(xiàn)ECG信號(hào)的特征選擇，減少檢測(cè)模型的輸入，提高檢測(cè)效率；通過PSO-SVM對(duì)睡眠呼吸暫停建立預(yù)測(cè)模型，使檢測(cè)準(zhǔn)確率提升至94.0%，優(yōu)于單一的SVM檢測(cè)模型。實(shí)驗(yàn)結(jié)果表明，單一的心電信號(hào)通過求取心率變異性及其他特征，可以用于睡眠呼吸暫停的檢測(cè)，且分類性能良好。

不足之處在于數(shù)據(jù)來源單一且測(cè)量時(shí)間較早，與當(dāng)前國內(nèi)患者可能存在一定的信號(hào)差異；此外，由于數(shù)據(jù)標(biāo)簽限制，僅實(shí)現(xiàn)了睡眠呼吸暫停的二分類檢測(cè)，無法根據(jù)模型結(jié)果判斷其嚴(yán)重程度。在未來的工作中，可以進(jìn)一步擴(kuò)展數(shù)據(jù)來源，深度優(yōu)化計(jì)算模型，提高泛化能力和檢測(cè)范圍。