黃曉霞， 王盼盼

(上海海事大學(xué) 信息工程學(xué)院， 上海 200135)

基于近似熵快速算法的靜息態(tài)腦磁信號(hào)分析

黃曉霞*， 王盼盼

(上海海事大學(xué) 信息工程學(xué)院， 上海 200135)

為了探究靜息態(tài)精神分裂癥患者腦磁信號(hào)的非線性動(dòng)力學(xué)特性，提出了一種將小波變換和近似熵相結(jié)合的特征提取方法.該方法首先通過(guò)小波變換，將10個(gè)正常人和10個(gè)精神分裂癥患者的腦磁信號(hào)進(jìn)行6層小波分解，提取對(duì)應(yīng)于腦磁信號(hào)θ波段和α波段的小波系數(shù)，繼而計(jì)算和比較兩類人近似熵的分布情況.實(shí)驗(yàn)結(jié)果表明，相同情況下精神分裂癥患者M(jìn)EG信號(hào)的各腦區(qū)和各通道間的近似熵都普遍高于正常人，α波段的額葉和中央?yún)^(qū)域尤為突出.該結(jié)果為進(jìn)一步研究患者M(jìn)EG信號(hào)特征進(jìn)而建立相應(yīng)的分類診斷模型提供了思路.

腦磁信號(hào)； 小波變換； 近似熵； 精神分裂癥

精神分裂癥是一種多基因遺傳病，其臨床癥狀十分復(fù)雜，個(gè)體心理素質(zhì)和外部社會(huì)環(huán)境都會(huì)引發(fā)該類疾病.發(fā)病規(guī)律和病理原因至今也尚未明晰，尋找理想的腦疾病診斷工具對(duì)精神分裂癥進(jìn)行研究成為當(dāng)今社會(huì)關(guān)注的焦點(diǎn).腦磁圖(magnetoencephalography，MEG)能夠?qū)崿F(xiàn)對(duì)大腦皮層的精確定位，同時(shí)能以ms級(jí)的時(shí)間分辨率和毫秒級(jí)的空間分辨率，無(wú)損傷無(wú)副作用的探測(cè)腦磁生理信號(hào)[1].可以將腦磁圖作為理想的診斷工具之一，它同功能性磁共振成像(functionalmagneticresonanceimaging，fMR)、腦電圖(electroencephalography，EEG)相比有其特有優(yōu)越性，在臨床醫(yī)學(xué)上被廣泛應(yīng)用于精神分裂癥、癲癇灶定位和腦腫瘤治療等領(lǐng)域[2].

腦磁圖的現(xiàn)代分析方法大致可以分為時(shí)頻分析、非線性動(dòng)力學(xué)分析和人工神經(jīng)網(wǎng)絡(luò)等[3-4].刻畫非線性動(dòng)力學(xué)特性的常見(jiàn)參數(shù)有關(guān)聯(lián)維數(shù)(correlationdimension，CD)、Lorenz散點(diǎn)圖、李雅普諾夫指數(shù)(lyapunovexponent，LE)、Lempel-Ziv復(fù)雜度和熵等.MahsaTaghavi等[5]通過(guò)近似熵(approximateentropy，ApEn)定量分析精神分裂癥患者和正常人的EEG信號(hào)，發(fā)現(xiàn)精神分裂癥患者產(chǎn)生更高的ApEn.Kim等[6]應(yīng)用李雅普諾夫指數(shù)，針對(duì)25名精神分裂癥患者和10名正常人進(jìn)行研究，發(fā)現(xiàn)兩類人的李雅普諾夫指數(shù)在額葉及顳葉有明顯的差異.Takahashi等人[7]應(yīng)用多尺度熵(multi-scaleentropy，MSE)比較了精神分裂癥患者和正常對(duì)照組，也發(fā)現(xiàn)精神分裂癥患者在額葉和中央?yún)^(qū)域都表現(xiàn)出較高的復(fù)雜性.

上述所述的關(guān)聯(lián)維數(shù)、李雅普諾夫指數(shù)和多尺度熵等非線性動(dòng)力學(xué)參數(shù)在表征信號(hào)時(shí)需要較多的數(shù)據(jù)點(diǎn)，造成實(shí)際應(yīng)用困難.近似熵很好的解決了該類問(wèn)題，近似熵只需較短的數(shù)據(jù)點(diǎn)就可以估計(jì)隨機(jī)信號(hào)和確定性信號(hào)特性，具有很好的抗噪和抗干擾能力.MEG信號(hào)是一種隨機(jī)性很強(qiáng)的非平穩(wěn)信號(hào)，很適合利用近似熵進(jìn)行分析，目前采用近似熵對(duì)精神分裂癥MEG信號(hào)進(jìn)行研究的報(bào)道尚未見(jiàn)有.小波變換是信號(hào)處理中一種典型的時(shí)頻分析方法，具有多分辨率特性[8].它既克服了窗口不隨頻率變化的缺點(diǎn)又繼承了短時(shí)傅立葉變換的局部化思想，是處理非平穩(wěn)信號(hào)的理想工具.因此，本文重點(diǎn)討論將小波變換和近似熵結(jié)合，來(lái)研究精神分裂癥患者和正常人MEG信號(hào)在α波段和θ波段的近似熵規(guī)律，從而發(fā)現(xiàn)相應(yīng)特征.

1實(shí)驗(yàn)研究方法

1.1近似熵快速算法

近似熵這一概念是Pincus在20C90S初為了解決混沌現(xiàn)象中求解熵問(wèn)題而提出的[9]，但該算法在計(jì)算過(guò)程中存在冗余步驟，計(jì)算成本大，計(jì)算效率低.因此，Pincus根據(jù)實(shí)踐又給出了一種近似熵實(shí)用快速算法，使其計(jì)算速度較原算法提升了好幾倍.具體算法如下[10-11].

步驟1：定義一個(gè)采樣點(diǎn)個(gè)數(shù)為N的MEG信號(hào)時(shí)間序列為x(1)，x(2)…，x(n).

步驟2：給定閾值r，將N點(diǎn)時(shí)間序列組成N*N的二值距離矩陣D，D中的元素用dij表示，定義為第i行元素和第j列元素之間的距離.

(i=1，…，N;j=1，…，N;i≠j).

(1)

(2)

步驟5：當(dāng)N為有限值時(shí)，此時(shí)間序列的近似熵為：

ApEn=(N，m，r)=lim[φm(r)-φm+1(r)].

(3)

近似熵是通過(guò)一個(gè)非負(fù)數(shù)來(lái)度量維數(shù)變化時(shí)產(chǎn)生新模式的概率大小，近似熵越大，時(shí)間序列越復(fù)雜，產(chǎn)生新模式的概率也就越大.反之，產(chǎn)生新模式的概率也就越小.上述的近似熵快速算法通過(guò)構(gòu)造距離矩陣來(lái)代替原算法中矢量的構(gòu)造，同時(shí)用時(shí)間序列中數(shù)據(jù)點(diǎn)間的差值代替原算法中重復(fù)計(jì)算相同維數(shù)矢量之間的距離，避免了很多不必要的計(jì)算步驟.

1.2參數(shù)的選取

近似熵的變化與參數(shù)的選取有很大的聯(lián)系，由公式(3)可知，在計(jì)算近似熵時(shí)需要考慮3個(gè)參數(shù)的選取，即時(shí)間序列長(zhǎng)度N，模式維數(shù)m，相似容限r(nóng)，這幾個(gè)參數(shù)一旦確定后，在整個(gè)計(jì)算過(guò)程中都是固定不變的.

對(duì)于時(shí)間序列長(zhǎng)度N，一般取值在100～5 000就能很好的估計(jì)出比較穩(wěn)定的統(tǒng)計(jì)值和較小的ApEn偽差[12]，因此在實(shí)際應(yīng)用中N的取值都要控制在此范圍內(nèi).對(duì)于模式維數(shù)m的選取，文獻(xiàn)[12-13]認(rèn)為當(dāng)其通過(guò)聯(lián)合概率對(duì)信號(hào)進(jìn)行動(dòng)態(tài)重構(gòu)時(shí)，m=2比m=1能夠獲得更加詳細(xì)的信息.對(duì)于m>2的情況，則會(huì)使通過(guò)近似熵分析出來(lái)的時(shí)間序列丟失許多有用信息，因此在本實(shí)驗(yàn)中取m=2.

作為一個(gè)有效過(guò)濾，r的取值對(duì)于整個(gè)實(shí)驗(yàn)的結(jié)果至關(guān)重要，Pincus等人在對(duì)隨機(jī)過(guò)程和確定性過(guò)程的理論分析基礎(chǔ)上，總結(jié)出r在0.1～0.25 SD(SD為x(i)數(shù)據(jù)的標(biāo)準(zhǔn)差)之間能夠估計(jì)出比較有效的統(tǒng)計(jì)特性[10]，若r值太小，就不能得到理想的概率統(tǒng)計(jì)估計(jì)，若r值太大，則會(huì)使信號(hào)丟失很多有用信息[14].本文通過(guò)實(shí)驗(yàn)對(duì)多個(gè)r值進(jìn)行比較分析，如圖1所示，隨著相似容限r(nóng)值的增大，兩類人近似熵值逐漸減小，但兩者的區(qū)分度卻逐漸增大，當(dāng)r=0.2時(shí)精神分裂癥患者和正常人之間的區(qū)分度更加明顯.當(dāng)r值大于0.2時(shí)，兩類人的區(qū)分度又逐漸減小，所以本次實(shí)驗(yàn)取r=0.2.

圖1 不同相似容限下精神分裂癥患者和正常人之間的差異Fig.1 Differences between schizophrenia patients and controls under similar tolerance

2預(yù)處理過(guò)程

2.1數(shù)據(jù)來(lái)源

實(shí)驗(yàn)數(shù)據(jù)來(lái)源于美國(guó)國(guó)立精神衛(wèi)生研究所MEG核心實(shí)驗(yàn)室(National Institute of Mental Health MEG core Facility).該數(shù)據(jù)由加拿大VSM MedTech 公司研發(fā)的CTF-275超導(dǎo)量子干涉儀(SQUID)設(shè)備采集得到，包括精神分裂癥患者及正常人兩類人的靜息態(tài)腦磁數(shù)據(jù).數(shù)據(jù)采樣頻率為600 Hz，采樣時(shí)間4 min.本次實(shí)驗(yàn)選取了10名精神分裂癥作為患者組，10名健康正常人作為對(duì)照組，采用的腦磁信號(hào)數(shù)據(jù)包含有273個(gè)有效通道(F43和O13通道損傷，不能正常工作)，考慮近似熵算法在時(shí)間序列長(zhǎng)度為100～5 000時(shí)分析效果更佳，本次實(shí)剔除了受試者可能未完全進(jìn)入靜息態(tài)的前1 min數(shù)據(jù)和出現(xiàn)精神分散的后1 min數(shù)據(jù)，保留中間部分重采樣，重采樣頻率為150 Hz，采樣時(shí)間30s，濾波帶通為0.5～256 Hz.

2.2ICA去除偽跡

獨(dú)立成分分析(independent component analysis，ICA)是一種盲源分離技術(shù)[15]，在語(yǔ)音分離、特征提取、噪聲消除等方面有著廣泛應(yīng)用.腦磁信號(hào)的采集過(guò)程中，不可避免的會(huì)遭受其他干擾源的污染引入其他電生理信號(hào)成分，如心電(electrocardiogram，ECG)、眼電(electro-oculogram，EOG)、肌電(electromyogram，EMG)等.因此為了獲得相對(duì)純凈的腦磁信號(hào)，使用ICA將相互獨(dú)立的干擾源信號(hào)從觀察信號(hào)中分離出來(lái).

圖2 偽跡處理后的獨(dú)立分量時(shí)域波形圖Fig.2 Independent component time-domain waveforms after artifact processing

圖3 原始信號(hào)與去噪后信號(hào)的對(duì)比圖Fig.3 The comparison diagram of MEG signals before and after artifacts denoising

圖2是經(jīng)過(guò)ICA處理后挑選出來(lái)的偽跡獨(dú)立成分，根據(jù)有關(guān)偽跡去除的先驗(yàn)知識(shí)可知：圖2中分量1和分量4的頻率較高，幅值較大是肌電偽跡；分量2是眨眼偽跡，分量7是眼動(dòng)偽跡，分量8的波形呈現(xiàn)一定的規(guī)律性是心電偽跡，分量9的波形比較平滑未出現(xiàn)異常波動(dòng)是正常的腦磁信號(hào).為了不影響實(shí)驗(yàn)結(jié)果，需要將分量1、2、4、7、8過(guò)濾掉.圖3是第29個(gè)通道偽跡去除前后的對(duì)比圖，藍(lán)色為原始信號(hào)，紅色為偽跡去除后的信號(hào).從圖中可以看出偽跡去除后的波形更加符合腦磁信號(hào)的波形特征.

2.3小波變換提取節(jié)律波

利用緊支集正交小波(Daubechies小波)的db3(3是小波階數(shù))小波函數(shù)對(duì)經(jīng)過(guò)ICA處理后的腦磁信號(hào)進(jìn)行6層小波分解，即x(n)=A6+D6+D5+D4+D3+D2+D1.理論上[16-17]，在靜息態(tài)腦磁信號(hào)中α波和θ波所占比重較大，因此本文主要研究這兩個(gè)波段.表1給出了小波變換后，各層頻帶所對(duì)應(yīng)的節(jié)律波以及各節(jié)律波的范圍.α波對(duì)應(yīng)的頻帶范圍為8～13Hz， D5小波系數(shù)重構(gòu)后的頻帶范圍為8～16Hz，在α波的頻帶附近，因此提取D5頻帶對(duì)應(yīng)的能量均值作為特征量.同理，θ波對(duì)應(yīng)的頻帶范圍為4～8Hz，因此提取D6頻帶對(duì)應(yīng)的能量均值作為特征量.

表1 小波分解各層頻帶范圍

3實(shí)驗(yàn)結(jié)果與分析

3.1通道間差異

圖4分別比較了患者組和正常對(duì)照組θ波段和α波段273個(gè)通道近似熵均值的分布情況.觀察上圖可知，在θ波段患者組各通道近似熵均值分布在[0.2993，0.3702]，而正常對(duì)照組各通道近似熵均值分布在[0.2652，0.3190].在α波段患者組各通道近似熵均值分布在[0.4545，0.5392]，而正常對(duì)照組各通道近似熵均值分布在[0.4416，0.4657].從散點(diǎn)圖可看出，不管是在θ波段還是α波段患者組273個(gè)通道的近似熵均值均高于正常對(duì)照組.θ波段兩類人近似熵值在其平均線上下波動(dòng)幅度較大，且波動(dòng)幅度和變化趨勢(shì)較為相似，α波段正常組近似熵值集中分布在平均線附近，上下波動(dòng)幅度小，而患者組近似熵值分布明顯較正常組更加分散，上下波動(dòng)幅度大，可見(jiàn)在α波段患者組與正常對(duì)照組之間的差異性比θ波段更加明顯，個(gè)體差異性更大.

圖4 患者組與正常對(duì)照組各通道近似熵分布Fig.4 ApEn distribution of each channel between patients and controls

3.2腦區(qū)間差異

實(shí)驗(yàn)將全頭的273個(gè)有效通道劃分為13個(gè)腦區(qū)組，其中M表示腦磁信號(hào)，L表示左，R表示右，Z表示中線，F(xiàn)為額葉，C為中央?yún)^(qū)，P為頂葉，O為枕葉，T為顳葉，共組合成13個(gè)腦區(qū)，即LF、ZF、RF、LC、ZC、RC、LP、RP、LO、ZO、RO、LT和RT.如圖5所示，圖中藍(lán)色表示正常人，紅色表示精神分裂癥患者，橫軸為13個(gè)腦區(qū)，縱軸為每個(gè)腦區(qū)對(duì)應(yīng)的近似熵均值.從圖中顯示的結(jié)果來(lái)看，患者組θ波段各腦區(qū)近似熵主要分布在[0.3195， 0.3643]，而正常組主要分布在[0.2665， 0.3107].患者組α波段各腦區(qū)近似熵主要分布在[0.4604， 0.5185]，而正常組主要分布在[0.4524， 0.4628]分布相對(duì)均勻，兩者有明顯的差異性.由圖6可知θ波段的MLC、MLP、MRF、MZF和MZO腦區(qū)以及α波段的MLC、 MLF、MRC、MRF、MZC和MZF腦區(qū)差異更加明顯.其中，

圖5 患者組與正常組各腦區(qū)ApEn均值分布Fig.5 ApEn mean distribution of each brain area between patients and controls

圖6 θ波段和α波段通各腦區(qū)間的近似熵差異Fig.6 Differences of each brain area between θ rhythm and α rhythm

3.3近似熵的統(tǒng)計(jì)圖分析

圖7是根據(jù)兩類人θ波段和α波段273個(gè)通道近似熵的平均值繪制而成的全頭拓?fù)鋱D.在圖中熵值的變化通過(guò)顏色來(lái)反映，當(dāng)顏色由藍(lán)到紅逐漸變深，表明近似熵的值由小逐漸變大.觀察圖7可知，在θ波段和α波段患者組的紅色區(qū)域分布均比正常組分布更加廣泛，在α波的額葉和中央?yún)^(qū)域以及θ波的頂葉和枕葉區(qū)域差別更加明顯.因?yàn)榇竽X額葉與思考、記憶和個(gè)體情感需求密切相關(guān)，頂葉對(duì)于人的思維邏輯判斷起到了重要的作用，而中央?yún)^(qū)域和枕葉也是大腦結(jié)構(gòu)的重要組成成份，當(dāng)這些部分受到影響時(shí)，大腦活動(dòng)的邏輯性和連貫性就會(huì)受到限制.患者組和正常對(duì)照組相比在這些區(qū)域的近似熵普遍偏高，說(shuō)明精神分裂癥患者腦磁信號(hào)時(shí)間序列較無(wú)規(guī)則.

圖7 患者組與正常組近似熵均值拓?fù)鋱D(左為θ波，右為α波)Fig.7 Topology of ApEn mean between patients and controls (left is θ rhythm and right is α rhythm)

圖8 患者組與正常組近似熵的箱型圖Fig.8 Boxplot of ApEn between patients and controls

將θ波段和α波段兩類人所有通道上的近似熵值平均后，繪制如圖8所示的箱型圖，該圖顯示了兩組數(shù)據(jù)的分散情況.從圖中可以看出，兩種節(jié)律波中正常組的方盒長(zhǎng)度明顯比患者組要短，說(shuō)明正常組近似熵的正常值分布比較集中，而患者組比較分散.而且，患者組所有通道的近似熵均值也都高于正常對(duì)照組，這也從另一方面說(shuō)明了精神分裂癥患者腦磁信號(hào)的混亂程度要高于正常人.

4結(jié)論

本文將非線性動(dòng)力學(xué)方法和時(shí)頻分析方法相結(jié)合，比較了相同狀態(tài)下精神分裂癥患者和正常人的腦磁信號(hào)θ波段和α波段的近似熵特征，實(shí)驗(yàn)初步表明在兩種節(jié)律波中精神分裂癥患者的MEG信號(hào)近似熵值都顯著高于正常人，在α節(jié)律波的額葉和中央?yún)^(qū)域差別更加明顯，這對(duì)于診斷精神分裂癥患者具有一定的參考價(jià)值.由于精神分裂癥患者的臨床表現(xiàn)錯(cuò)綜復(fù)雜，個(gè)體差異，生理狀況都有可能影響實(shí)驗(yàn)結(jié)果.所以今后將從增加實(shí)驗(yàn)樣本容量，細(xì)分信號(hào)不同波段等方面入手，使得實(shí)驗(yàn)結(jié)果更具有普遍性.

致謝 感謝美國(guó)國(guó)立精神衛(wèi)生研究所 MEG 科研平臺(tái)提供的 MEG 研究數(shù)據(jù)，感謝Tzyy-Ping Jung 教授，Arnold J.Mandell 教授以及加州大學(xué)圣迭戈分校計(jì)算機(jī)系的 Chung-Kuan Cheng教授等對(duì)該課題研究的支持與幫助.

[1] 孫吉林. 腦磁圖(MEG)常見(jiàn)問(wèn)題[J]. 現(xiàn)代電生理學(xué)雜志， 2005， 12(4):245-249．

[2] DAI Y， ZHANG W， DICKENS D L， et al. Source connectivity analysis from MEG and its application to epilepsy source localization[J]. Brain Topography， 2012， 25(2):157-66．

[3] AFTANAS L I， LOTOVA N V， KOSHDAROV V I， et al. Nonlinear forecasting measurements of the human EEG during evoked emotions[J]. Brain Topography， 1997， 10(2):155-62．

[4] TAKEUCHI F， KAMADA K， KURIKI S. Time-frequency analysis of MEG signals during a visual language task[J]. Transactions of the Japanese Society for Medical & Biological Engineering Bme， 2006， 44:696-701．

[5] TAGHAVI M， BOOSTANI R， SABETI M， et al. Usefulness of approximate entropy in the diagnosis of schizophrenia[J]. Iranian Journal of Psychiatry & Behavioral Sciences， 2011， 5(5):62-70．

[6] KIM D J， JEONG J， CHAE J H， et al. An estimation of the first positive Lyapunov exponent of the EEG in patients with schizophrenia[J]. Psychiatry Research， 2000， 98(3):177-189．

[7] TAKAHASHI T， CHO R Y， MIZUNO T， et al. Antipsychotics reverse abnormal EEG complexity in drug-naive schizophrenia: a multiscale entropy analysis [J]. Neuroimage， 2010， 51(1):173-82．

[8] 張 毅， 羅明偉， 羅 元. 腦電信號(hào)的小波變換和樣本熵特征提取方法[J]. 智能系統(tǒng)學(xué)報(bào)， 2012， 7(4):339-344．

[9] PINCUS S M， MULLIGAN T， IRANMANESH A， et al. Older males secrete luteinizing hormone and testosterone more irregularly， and jointly more asynchronously， than younger males[J]. Proceedings of the National Academy of Sciences of the United States of America， 1996， 93(24):14100(1-5)．

[10] PINCUS， S M. Approximate entropy (ApEn) as a complexity measure[J]. Chaos， 1995， 5(1):110-117．

[11] AN X U， LIU J， PENG Q Y， et al. Fast algorithm of dynamic approximate entropy and its application in heart rate variance analysis[J]. Journal of Tongji University， 2005, 33(4):520-524.

[12] CHEN W T, WANG Z Z, REN X M. Characterization of surface EMG signals using improved approximate entropy[J].Journal of Zhejiang University Science(Life Science)2006, 10: 844-848.

[13] YENTES J M， HUNT N， SCHMID K K， et al. The appropriate use of approximate entropy and sample entropy with short data sets[J]. Annals of Biomedical Engineering， 2013， 41(2):349-365.

[14] PINCUS S M. Approximate entropy as a measure of system complexity[J]. Proceedings of the National Academy of Sciences of the United States of America， 1991， 88(6):2297-301．

[15] XIANFENG D U, ZHIJUN L I, FENGRONG B I, et al. Source separation of diesel engine vibration based on the empirical mode decomposition and independent component analysis[J].Chinese Journal of Mechanical Engineering, 2012(3):557-563.

[16] 湯 浩， 盧 青， 韓 莉， 等. 抑郁癥住院患者靜息態(tài)的腦磁頻譜分析[J]. 中華精神科雜志， 2012， 45(6):345-348．

[17] HINKLEY L B N， VINOGRADOV S， GUGGISBERG A G， et al. Clinical symptoms and alpha band resting-state functional connectivity imaging in patients with schizophrenia: implications for novel approaches to treatment[J]. Biological Psychiatry， 2011， 70(12):1134-1142．

Resting-state magnetic signals analysis based on the fast algorithm of approximate entropy

HUANG Xiaoxia, WANG Panpan

(Shanghai Maritime University， College of Information Engineering， Shanghai 200135)

In order to study the nonlinear dynamics of the schizophrenic patient’s MEG signals in resting-state, this paper presents a method of feature extraction which combined the wavelet variation with the approximate entropy. The brain magnetic signals of 10 controls and 10 patients are decomposed to six levels by wavelet decomposition and wavelet coefficient is extracted corresponding to the θ rhythm and α rhythm of MEG signals. Then the distribution of approximate entropy between two kinds of people are calculated and compared. The experiment results show that the entropy of each brain region and channel of the MEG signals in schizophrenic patients were generally higher than controls under the same situation, especially frontal and central regions in α rhythm. This result provides a guideline for the study of EEG signal characteristics of the patients and establishes the appropriate classification diagnostic model.

magnetoencephalography; wavelet transform; approximate entropy; schizophrenia

2016-10-21.

第48 批教育部留學(xué)回國(guó)人員科研啟動(dòng)基金項(xiàng)目．

1000-1190(2017)03-0309-08

R749.3

A

*E-mail: 1946883554@qq.com．