楊 琨 曹建庭，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項外還增加了瞳孔散大的項目，但是一些基本的判定項目是一致的［3］。

現(xiàn)有的腦死判定程序比較耗時，而且存在一定的風(fēng)險。例如，在實施“有無自主呼吸”的檢查時，需要拔去患者的呼吸器，由此而導(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)點在于可避免腦死亡的誤判，并且能減輕醫(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)行獨立成分分析(ICA)［5］，把其結(jié)果與醫(yī)生通過臨床診斷患者的結(jié)果進(jìn)行對照，發(fā)現(xiàn)無論從數(shù)據(jù)分析的角度還是從臨床判定的角度，對腦死者的判定結(jié)果是一致的。雖然獨立成分分析具有較強(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)過概率驗證(P值驗證等)后，從宏觀上得出了昏迷患者和腦死亡者的兩類腦電確實具有特征差異的結(jié)論［6］。隨著采集的腦電數(shù)據(jù)量的增加，對所有35例患者、共計61次腦電數(shù)據(jù)進(jìn)行近似熵分析，進(jìn)一步確認(rèn)和衡量腦死者和昏迷患者腦電隨機(jī)性程度的差異［16］。

通過對所有被測腦電信號的靜態(tài)近似熵的分析，雖然從宏觀上得到了昏迷患者與腦死者兩類腦電信號的差異，但是在微觀方面，靜態(tài)近似熵值僅是一個值，不能實時反映疑似腦死患者狀態(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點即可，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ù)域，a被稱為伸縮因子(又稱尺度因子)且 a≠0，b為平移因子;ψ(a，b)(t)表示由基小波生成的小波函數(shù)系表示小波函數(shù)系的共軛。Wf(a，b)為待分析信號f(t)在固定小波函數(shù)ψ(a，b)(t)上的分量，表達(dá)了 f(t)與該小波函數(shù)的相似程度。

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

式中:j，k是將 a，b離散化后的參數(shù)，R為實數(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為時間窗口的起始點，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)?？梢?，筆者對數(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)失誤。為了避免誤判，在對所測的腦電信號實施動態(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ī)生的臨床判定是完全一致的。今后的研究重點是高噪聲下腦死亡數(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ù)字信號處理的原理與實現(xiàn)［M］.南京:東南大學(xué)出版社，2001.

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