疑似腦死患者腦電信號的動態(tài)近似熵分析

楊 琨 曹建庭，2，3* 王如彬 朱惠莉

1(華東理工大學(xué)認(rèn)知神經(jīng)動力學(xué)研究所，上海 200237)

2(日本埼玉工業(yè)大學(xué)大學(xué)院工學(xué)研究科，日本埼玉縣深谷市 396-0293)

3(日本理化研究所腦科學(xué)研究中心.日本埼玉縣和光市 351-0198)

4(華東醫(yī)院呼吸科，上海 200040)

引言

腦死亡，是指人的全腦和腦干功能完全的、不可逆轉(zhuǎn)的喪失。腦死亡是人生命的終結(jié)，所以對腦死亡判定的研究具有重要的意義。在1968年的第22屆世界醫(yī)學(xué)大會上，美國哈佛醫(yī)學(xué)院腦死亡定義審查特別委員會提出，將“腦功能完全不可逆地喪失”作為新的死亡標(biāo)準(zhǔn)，并制定了世界上第一個腦死亡判定標(biāo)準(zhǔn):“1.深度昏迷;2.自主呼吸停止;3.腦干反射消失;4.腦電消失或平坦(2 μV以下)”［1－2］。盡管不同國家對腦死亡的判定標(biāo)準(zhǔn)略有不同，如日本的判定標(biāo)準(zhǔn)除了上述4項(xiàng)外還增加了瞳孔散大的項(xiàng)目，但是一些基本的判定項(xiàng)目是一致的［3］。

現(xiàn)有的腦死判定程序比較耗時，而且存在一定的風(fēng)險。例如，在實(shí)施“有無自主呼吸”的檢查時，需要拔去患者的呼吸器，由此而導(dǎo)致患者死亡的例子也曾經(jīng)發(fā)生過［4］。還有，“確認(rèn)平坦腦電”的檢查一般至少需要12 h，個別也有超過200 h才能確認(rèn)的。在此，引入一種更為簡單有效的方法，即利用腦電預(yù)測系統(tǒng)對患者進(jìn)行腦死亡預(yù)判定。在預(yù)判定過程中，對還具有腦電節(jié)律的患者可中止判定程序并返回對其進(jìn)行治療，對未提取出腦電節(jié)律的患者仍可以按現(xiàn)有的腦死判定程序進(jìn)行判定。引入腦電預(yù)測系統(tǒng)的最大優(yōu)點(diǎn)在于可避免腦死亡的誤判，并且能減輕醫(yī)生和患者的負(fù)擔(dān)［5－6］。

復(fù)雜度是用來衡量時間序列隨機(jī)性程度的一個指標(biāo)。如果一個序列的復(fù)雜度高，那么這個序列的隨機(jī)程度也較高。常見的復(fù)雜度指標(biāo)包括倫伯爾－齊夫定義的 CLZ復(fù)雜度、近似熵(approximate entropy，ApEn)、C0復(fù)雜度等。作為一個隨機(jī)性復(fù)雜度，近似熵描述的是一個時間序列在模式上的自相似程度，被應(yīng)用于分析睡眠腦電等研究中。例如，對正常人睜眼、閉眼、淺睡、深睡等4種狀態(tài)腦電的近似熵研究顯示，人在深睡時沒有意識，大腦的活動相對不是很“復(fù)雜”，其復(fù)雜度很低;而在清醒時，大腦接收大量外界的信息，其近似熵比閉眼時的要大［7－9］。

在對疑似腦死患者進(jìn)行腦死亡判定的研究中，首先對所有患者的腦電信號進(jìn)行獨(dú)立成分分析(ICA)［5］，把其結(jié)果與醫(yī)生通過臨床診斷患者的結(jié)果進(jìn)行對照，發(fā)現(xiàn)無論從數(shù)據(jù)分析的角度還是從臨床判定的角度，對腦死者的判定結(jié)果是一致的。雖然獨(dú)立成分分析具有較強(qiáng)的去噪和成分提取的功能，但對于未提取到腦電節(jié)律的部分患者的判定缺乏嚴(yán)謹(jǐn)，所以又引入能量指標(biāo)［10］和復(fù)雜度指標(biāo)，對疑似腦死患者的腦電信號進(jìn)行分析。在復(fù)雜度方面，首先運(yùn)用了近似熵(ApEn 法)、C0復(fù)雜度［11－12］、NSSE 法［13］、DFA 法［14－15］4 種不同的復(fù)雜度，對采集到的部分患者的腦電數(shù)據(jù)進(jìn)行處理，經(jīng)過概率驗(yàn)證(P值驗(yàn)證等)后，從宏觀上得出了昏迷患者和腦死亡者的兩類腦電確實(shí)具有特征差異的結(jié)論［6］。隨著采集的腦電數(shù)據(jù)量的增加，對所有35例患者、共計61次腦電數(shù)據(jù)進(jìn)行近似熵分析，進(jìn)一步確認(rèn)和衡量腦死者和昏迷患者腦電隨機(jī)性程度的差異［16］。

通過對所有被測腦電信號的靜態(tài)近似熵的分析，雖然從宏觀上得到了昏迷患者與腦死者兩類腦電信號的差異，但是在微觀方面，靜態(tài)近似熵值僅是一個值，不能實(shí)時反映疑似腦死患者狀態(tài)的變化。本研究引入動態(tài)近似熵，從而在微觀方面能反映腦電信號的近似熵變化，也能觀察到疑似腦死患者生命體征的變化。引入動態(tài)近似熵分析法后，發(fā)現(xiàn)個別疑似腦死患者腦電信號的分析結(jié)果和期望值有一定偏差，其原因在于搶救室里各種儀器設(shè)備對腦電記錄存在干擾。為減弱噪聲的影響，引入小波分析法先對腦電信號進(jìn)行去噪前處理，然后再進(jìn)行動態(tài)近似熵分析。

1 動態(tài)近似熵

靜態(tài)近似熵分析法被用于反映一個時間序列在模式上的自相似程度，即當(dāng)維數(shù)變化時序列中產(chǎn)生新模式的可能性大小。ApEn越大，產(chǎn)生新模式的機(jī)會越大，因而序列越復(fù)雜［8］。靜態(tài)近似熵分析只是對其中的某段數(shù)據(jù)求值，不能反映整個測量時間段近似熵的變化過程。在此，引入一個時間變量t，對現(xiàn)有的靜態(tài)近似熵分析法在時間軸上進(jìn)行擴(kuò)展，從而得到動態(tài)近似熵。

對于一個時間長度為T的時間序列s(t)，定義其動態(tài)近似熵為 ApEn(m，r，t)，其中 m為矢量維數(shù)，r為閾值，即選定的相似容限，t為時間。引入一個寬度為t′的時間窗口，窗口沿著時間軸 t不斷移動，進(jìn)入窗口的時間序列可表示為

式中:t=0，…，T－t′;N為時間窗口對應(yīng)的序列長度，見圖1。

在序列 Xt中，按照順序構(gòu)造兩個 m維矢量Vt(i)和 vt(j)，則有

式中，i，j≤N － m+1。

D［vt(i)，vt(j)］表示矢量 vt(i)和 vt(j)之間的距離，則有

圖1 動態(tài)近似熵的窗口移動Fig.1 The diagram of dynamic ApEn

對于給定閾值 r及矢量 vt(i)，當(dāng) j≤N－m+1時，計算出與vt(i)的距離小于r的矢量vt(j)個數(shù)，則Cm

i(r，t)表示對于所有 j≤N－m+1時的矢量，vt(j)與vt(i)距離小于r的概率，有

對于每個 m+1，重復(fù)上述步驟，得到 φm+1(r，t)，則動態(tài)近似熵為

在一般情況下，時間序列取100～5 000點(diǎn)即可，m=2，r取時間序列標(biāo)準(zhǔn)差的 0.1 ～0.2 倍［7－8］。由于近似熵是把隨機(jī)性作為對復(fù)雜型的度量，即衡量一個時間序列的隨機(jī)性程度，因此對于隨機(jī)序列和正弦波來說，正弦波的隨機(jī)性要小于隨機(jī)序列。具體到疑似腦死患者腦電信號的分析中，由于昏迷患者還有腦活動，測得的腦電信號是具有一定周期性的腦電節(jié)律與噪聲的疊加，其近似熵較小;而腦死者的腦電信號是各種干擾源產(chǎn)生的噪聲總和，其近似熵應(yīng)該較大。根據(jù)這個原理，可以利用近似熵分析法來識別是否腦死。

2 小波分析

采用小波分析法對信號去噪，其本質(zhì)在于小波變換對信號和噪聲的瞬時特性表現(xiàn)不一樣。本研究所使用的腦電數(shù)據(jù)是從重癥監(jiān)護(hù)室中獲得的，采集到的腦電受到外部各種干擾源的影響。對于一個含噪聲e(t)的原始信號f(t)，其基本模型可以表示為

式中，e(t)為噪聲，σ稱為噪聲強(qiáng)度。

在最簡單的情況下，可以假設(shè)e(t)為高斯白噪聲，且 σ=1。小波法去噪的目的，就是要抑制e(t)以恢復(fù)s(t)。

2.1 小波變換

小波變換的基本思想是用一簇函數(shù)表示或逼近待分析信號或函數(shù)。這一簇函數(shù)被稱為小波函數(shù)系，是由基小波或母小波的函數(shù)通過平移和伸縮而構(gòu)成的一組正交基［17－18］。將待分析信號在小波函數(shù)系上做分解，即可得到一維連續(xù)小波變換，有

式中:a，b∈R，R為實(shí)數(shù)域，a被稱為伸縮因子(又稱尺度因子)且 a≠0，b為平移因子;ψ(a，b)(t)表示由基小波生成的小波函數(shù)系表示小波函數(shù)系的共軛。Wf(a，b)為待分析信號f(t)在固定小波函數(shù)ψ(a，b)(t)上的分量，表達(dá)了 f(t)與該小波函數(shù)的相似程度。

在實(shí)際應(yīng)用中，為了方便計算機(jī)處理，尺度因子和平移因子需要進(jìn)行離散化處理，則離散的小波變換可表示為

式中:j，k是將 a，b離散化后的參數(shù)，R為實(shí)數(shù)域，ψj，k(t)為離散小波函數(shù)，時間t仍為連續(xù)變量。

2.2 小波去噪方法

利用小波變化進(jìn)行去噪的方法有很多，本研究利用閾值法去噪，即對含噪信號f(t)作離散小波變換，所得到的小波系數(shù)一部分是與信號對應(yīng)的小波系數(shù)，另一部分是與噪聲對應(yīng)的小波系數(shù)，具有不同的統(tǒng)計特性，通常信號的系數(shù)值要大于噪聲的系數(shù)值。需要找到一個合適的閾值λ，將絕對值小于閾值的小波系數(shù)舍棄，將絕對值大于閾值的小波系數(shù)保留或收縮，得到小波系數(shù)的估計值;根據(jù)這個估計值進(jìn)行信號重構(gòu)，從而得到去噪后的信號［19－20］。

3 患者腦電數(shù)據(jù)采集和分析計算流程

研究所用到的腦電數(shù)據(jù)全部是從醫(yī)院的重癥監(jiān)護(hù)室里測得的，所使用的腦電儀是便攜式NEUROSCAN ESI-32系統(tǒng)。在這個系統(tǒng)中，一共有7個電極被放在昏迷患者的前額，分別是 Fp1，F(xiàn)p2，F(xiàn)3，F(xiàn)4，F(xiàn)7，F(xiàn)8 和地極，A1，A2 作為參考電極放在耳朵上，見圖2。腦電信號的采樣頻率是1 000 Hz，電極阻抗低于8 000 Ω［5］。

圖2 電極安放系統(tǒng)［14］Fig.2 The electrode layout［14］

在2004年6月 ～2006年3月期間，共采集了35名患者的腦電信號。這些患者包括19名昏迷患者和16名腦死者，年齡從17～85歲不等。根據(jù)這些患者的健康狀況，記錄了不同的腦電數(shù)據(jù)。圖3是兩例疑似腦死患者的腦電信號。

加入小波前處理的動態(tài)近似熵算法流程如下，其中N為時間窗口 t＇對應(yīng)的序列長度，r為閾值，維數(shù)m初值為2，nt為時間窗口的起始點(diǎn)，wname為選擇的小波函數(shù)，lev為小波分解層數(shù)。

步驟1:讀取腦電信號 data，設(shè)置初值 N，r，m，nt，wname，lev。

步驟2:由 wname，lev，對腦電信號 data采用小波分析法去噪，得到去噪后的信號datalb。

步驟3:從datalb中提取起始值為nt，長度為N的時間序列datnew。

步驟4:給定初值 i=1，sum=0。

步驟 5:由 i，得到 xi=datnew(i:i－ m+1)。

步驟6:由每個不大于N－m+1的j值得到xj=datnew(j:j－m+1)，計算出 D(xi，xj)≤r的個數(shù) n，c=n/(N － m+1)，sum=sum+logc。

步驟 7:i=i+1，若 i≤N － m+1，重復(fù)步驟 5、6，直到 i>N －m+1，則 φm=sum/(N －m+1)。

步驟8:m=m+1時，重復(fù)步驟4～步驟7，得到φm+1，則 apen=φm－φm+1。

步驟9:移動窗口，使 nt=nt+N，重復(fù)步驟3～步驟8，得到一系列apen，并作圖。

4 數(shù)據(jù)分析結(jié)果

4.1 昏迷患者和腦死者的動態(tài)近似熵分析

圖3 腦電預(yù)測系統(tǒng)所測的疑似腦死患者的腦電信號。(a)患者A的腦電信號;(b)患者B的腦電信號Fig.3 The EEG recordings of two patients.(a)the EEG recording of patient A;(b)the EEG recording of patient B

對所測得的35例患者、共計61次腦電信號進(jìn)行動態(tài)近似熵分析，發(fā)現(xiàn)一部分腦電信號的近似熵較小，在0附近較小的范圍內(nèi)變化，如圖4(a)所示的患者A動態(tài)近似熵的分析結(jié)果;另一部分腦電信號的近似熵值較大，在1附近較小的范圍內(nèi)變化，如圖4(b)所示的患者B動態(tài)近似熵的分析結(jié)果。根據(jù)近似熵的定義，較小的近似熵值說明腦電信號中存在周期性的腦電節(jié)律，即患者還有腦活動，處于昏迷狀態(tài);而較大的近似熵值說明腦電信號中存在隨機(jī)性的噪聲，患者可能為腦死者。將這些數(shù)據(jù)分析的結(jié)果同醫(yī)生對患者病狀的診斷進(jìn)行對照，發(fā)現(xiàn)用動態(tài)近似熵分析法的結(jié)果和醫(yī)生的判定是一致的。

圖4 利用動態(tài)近似熵對患者的腦電信號進(jìn)行分析的結(jié)果。(a)對患者A的腦電信號的分析結(jié)果;(b)對患者B的腦電信號的分析結(jié)果Fig.4 Dynamic ApEn results for two patients.(a)Dynamic ApEn results for patient A;(b)Dynamic ApEn results for patient B

采集的腦電信號還包括了一例體征狀態(tài)改變(即由昏迷變?yōu)槟X死狀態(tài))患者的腦電信號。2005年10月18日，患者意識不明、瞳孔散大，但有微弱的對光反射，處于深昏迷狀態(tài);10月19日，患者狀況惡化，各種腦干反射均消失，被醫(yī)生判定為腦死亡。在患者昏迷狀態(tài)以及被判定為腦死狀態(tài)時，對其的腦電信號進(jìn)行了多次測量，利用文中第1節(jié)提出的動態(tài)近似熵進(jìn)行分析，得到如圖5所示的結(jié)果。圖5左側(cè)的近似熵值較小，在整個過程中，總在0附近較小的范圍內(nèi)變化;而圖5右側(cè)的近似熵值較大，大部分時間都在1附近較小的范圍內(nèi)變化。根據(jù)近似熵的定義，圖5左側(cè)的腦電信號存在周期性腦電節(jié)律，患者處在昏迷狀態(tài)。而圖5右側(cè)的腦電信號存在隨機(jī)噪聲，患者可能處于腦死狀態(tài)?？梢姡P者對數(shù)據(jù)分析的結(jié)果同醫(yī)生的臨床判定結(jié)論是一致的。包括上述結(jié)果，35個病例中的31個病例只需利用開發(fā)的動態(tài)復(fù)雜度便可得到類似的識別結(jié)果。

4.2 與小波分析結(jié)合的動態(tài)近似熵分析

雖然從圖4和圖5中可以看出，不管是處于不同狀態(tài)的不同患者，還是體征狀態(tài)改變的同一患者，處于昏迷狀態(tài)時的近似熵總是在0附近的較小范圍內(nèi)變化，而處于腦死狀態(tài)時的近似熵則在1附近的較小范圍內(nèi)變化。但是，在圖5中，腦死狀態(tài)的腦電信號在300 s處的動態(tài)近似熵忽然降低，和左側(cè)昏迷狀態(tài)的動態(tài)近似熵相差不大。如果僅對300 s前后的這段腦電信號進(jìn)行傳統(tǒng)的動態(tài)近似熵分析，對患者狀態(tài)的判定會出現(xiàn)失誤。為了避免誤判，在對所測的腦電信號實(shí)施動態(tài)近似熵分析之前，使用小波分析法去除由測量環(huán)境中儀器設(shè)備產(chǎn)生的高頻周期噪聲。圖6是采用小波分析法濾波后，昏迷患者A和腦死亡者B的動態(tài)近似熵比較。對比圖4和圖6可以看出，噪聲干擾較小時，對腦電信號的動態(tài)近似熵的分析結(jié)果影響不大。

圖5 利用動態(tài)近似熵對患者C的腦電信號進(jìn)行分析的結(jié)果Fig.5 Dynamic ApEn results for patient C

圖6 利用小波分析法和動態(tài)近似熵對患者的腦電信號進(jìn)行分析的結(jié)果。(a)利用小波分析法和動態(tài)近似熵對患者A的腦電信號的分析結(jié)果;(b)利用小波分析法和動態(tài)近似熵對患者B的腦電信號的分析結(jié)果Fig.6 Dynamic ApEn with a wavelet denosing technique for two patients.(a)Dynamic ApEn with a wavelet denosing technique for patient A;(b)Dynamic ApEn with a wavelet denosing technique for patient B

用小波法去噪后，對昏迷變?yōu)槟X死狀態(tài)的患者腦電信號進(jìn)行動態(tài)近似熵比較，見圖7。對比圖5和圖7，在300 s處，采用小波分析法去噪后，高頻周期噪聲的影響降低，腦死狀態(tài)腦電信號的近似熵總在1附近變化。包括上述結(jié)果在內(nèi)，35個病例中的4個病例在利用小波分析與動態(tài)近似熵結(jié)合的方法后，其識別效果更好。因此，在進(jìn)行動態(tài)近似熵分析之前，對測得的數(shù)據(jù)進(jìn)行小波去噪前處理是非常必要的。

圖7 利用小波分析法和動態(tài)近似熵對患者C的腦電信號進(jìn)行分析的結(jié)果Fig.7 Dynamic ApEn with a wavelet denosing technique for patient C

為了進(jìn)一步分析圖5中患者的腦死狀態(tài)腦電信號在第300 s前后的噪聲成分，對導(dǎo)聯(lián)Fp1上這段長約103 s的信號采用小波分析法去噪，濾除的噪聲如圖8(a)所示;對噪聲進(jìn)行傅里葉變換，得出其頻譜圖，如圖8(b)所示。從圖8(b)中可以看出，被濾除的噪聲帶有很強(qiáng)的150 Hz的成分，而這一成分的干擾主要是來自外部設(shè)備的電源干擾。

圖8 患者C在腦死狀態(tài)300s前后濾除的噪聲。(a)濾除的噪聲;(b)對噪聲進(jìn)行傅立葉變換的結(jié)果Fig.8 The interfering noise in patient C’s EEG data around the 300th second.(a)the interfering noise;(b)Fourier transformation of the interfering noise

5 結(jié)論

在本研究中，首先在靜態(tài)近似熵分析法中引入時間參量，并把它擴(kuò)展成動態(tài)近似熵分析法，用來識別昏迷患者與腦死者的腦電信號，以及觀察患者病狀的變化過程。由于腦電信號中存在由儀器設(shè)備產(chǎn)生的高頻周期噪聲，所以首先采用小波分析法對腦電信號進(jìn)行去噪前處理，然后再對去噪后的腦電信號進(jìn)行動態(tài)近似熵分析。在整個測量和分析過程中，昏迷患者的近似熵小于腦死者的近似熵，昏迷患者的近似熵在接近于0的很小范圍內(nèi)變化，而腦死者的近似熵則在1附近的很小范圍內(nèi)變化。共對35例患者進(jìn)行了數(shù)據(jù)分析和識別，得出的結(jié)論與醫(yī)生的臨床判定是完全一致的。今后的研究重點(diǎn)是高噪聲下腦死亡數(shù)據(jù)的分析與識別。

［1］Ad hoc committee of the Harvard medical school to examine the definition of brain death.A definition of irreversible coma［J］.JAMA，1968，205(6):337－340.

［2］陳忠華.腦死亡:現(xiàn)代死亡學(xué)［M］.北京:科學(xué)出版社，2004.

［3］Wijdicks EFM.Brain death worldwide:accepted fact but no global consensus in diagnostic criteria ［J］.Neurology，2002，58:20－25.

［4］Marks SJ，Zisfein J.Apneic oxygenation in apnea tests for brain death:a controlled trial［J］.Archives of Neurology，1990，47(10):1066－1068.

［5］Cao Jianting.Analysis of the quasi-brain-death EEG data based on a robust ICA approach ［J］.Lecture Notes in Computer Science，2006，10:1240－1247.

［6］Cao Jianting，Chen Zhe.Advanced EEG signal processing in brain death diagnosis［M］//Signal Processing Techniques for Knowledge Extraction and Information Fusion.Berlin:Springer，2008:275－298.

［7］Pincus SM.Approximate entropy(ApEn)as a measure of system complexity［J］.Proc Natl Acad Sci USA，1991，88:110－117.

［8］顧凡及，梁培基.神經(jīng)信息處理［M］.北京:北京工業(yè)大學(xué)出版社，2007.

［9］張涇周，馬穎穎，李婷，等.基于復(fù)雜性測度的睡眠腦電分期處理方法研究［J］.中國生物醫(yī)學(xué)工程學(xué)報，2009，28(3):367－371.

［10］Shi Qiwei，Yang Juhong，Cao Jianting，et al.EMD based power spectral pattern analysis forquasi-brain-death EEG ［J］.Emerging Intelligent Computing Technology and Applications.With Aspects of Artificial Intelligence，2009，5755:814－823.

［11］Gu Fanji，Meng Xin，Shen Enhua，et al.Can we measure consciousness with EEG complexities?［J］International Journal of Bifurcation and Chaos，2003，13:733－742.

［12］Chen Fang，Xu Jinghua，Gu Fanji，et al.Dynamic process of information transmission complexity in human brains［J］.Biol Cybern，2000，83:355－366.

［13］RobertsSJ， Penny W， Rezek L. Temporaland spatial complexity measures for EEG-based brain-computer interfacing［J］.Medical and Biological Engineering and Computing，1998，37(1):93－99.

［14］Peng CK，Buldyrev SV，Havlin S，et al.Mosaic organization of DNA nucleotides［J］.Phys Rev E，1994，49:1685 －1689.

［15］Little M，McSharry P，Moroz I，et al.Nonlinear，biophysicallyinformed speech pathology detection［C］//Proceedings of ICASSP’06.New York:IEEE Publishers，2006:1080 －1083.

［16］Yang Kun，Shi Qiwei，Cao Jianting，et al.Analyzing EEG of quasi-brain-death based on approximate entropy measures［C］// Advances in Cognitive Neurodynomics. Berlin:Springer，2009(in press)

［17］徐晨，趙瑞珍，干小兵.小波分析·應(yīng)用算法［M］.北京:科學(xué)出版社，2004.

［18］高志，余嘯海.Matlab小波分析工具箱原理與應(yīng)用［M］.北京:國防工業(yè)出版社，2004.

［19］吳鎮(zhèn)揚(yáng).數(shù)字信號處理的原理與實(shí)現(xiàn)［M］.南京:東南大學(xué)出版社，2001.

［20］張翀，余紅英.基于小波分析的心電信號降噪研究［J］.山西大學(xué)學(xué)報(自然科學(xué)版)，2008，31(S2):57 －62.