針對所有原始腦電信號都受到低頻和尖峰噪聲干擾的問題, 提出了小波變換和獨立分量分析相結合的去噪算法; 對預處理后的腦電數據, 進行小波熵、近似熵和復雜度這三種特征參數的數值表征, 并進一步通過特征參數的狀態變化率來判斷腦電信號的狀態區分效果。麻醉與非麻醉的腦電數據處理結果表明, 三種特征參數的狀態變化率分別達到50.5%、21.6%和19.5%, 其中小波熵的狀態變化率最高, 這些特征參數可作為基于腦電信號分析的麻醉深度量化研究的基礎。
引用本文： 王鋒, 李曉歐. 基于腦電信號的麻醉特征參數分析. 生物醫學工程學雜志, 2015, 32(1): 13-18, 31. doi: 10.7507/1001-5515.20150003 復制
引言
麻醉是通過抑制中樞神經系統而達到意識消失、止痛的目的，在外科手術中是不可或缺的。精確、實時的麻醉深度監護對于提高麻醉質量、保障手術安全，以及減少麻醉并發癥和控制麻醉藥品用量都具有重要的意義。腦電信號是腦內部億萬個神經元興奮或抑制而記錄在頭皮表面的電生理信號，它能夠直接反映大腦意識的興奮或抑制程度，對其進行處理和分析，無論是在臨床上對腦疾病的診斷和治療，還是對腦認知功能的研究都是十分重要的[1-2]。因此，利用EEG信號分析大腦受抑程度是最為直接的，并可作為一種非插入式的監測手段，從而來反饋、改變手術過程中麻醉藥物的使用。
在腦電信號采集過程中，不可避免地會受到各種干擾，例如工頻干擾、基線漂移、眼動和肌電等。而且腦電信號本身是非線性、非平穩且十分微弱的，狀態不同，噪聲的成分和類型也不同，導致了信號被淹沒其中。均值法、中值法等傳統的去噪算法只是對腦電信號在空間域上進行處理，自適應差，局部性強，結果不能令人滿意。小波變換具有多尺度分析的特點，能夠聚焦到信號的任意細節進行多分辨率的時頻分析，對處理像腦電這種非平穩信號的干擾是非常適合的[3]。獨立分量分析(independent component analysis, ICA)方法需要滿足一定分離條件，即干擾信號與腦電信號之間是統計獨立的，從而將腦電信號按照統計獨立的原則通過優化算法分解為若干獨立來源成分，這樣就將去噪問題轉化為獨立分量分離問題，特別適合處理那些突發的，尖銳的噪聲。所以，采用兩者的結合，可以有效去除低頻和尖峰等高頻噪聲[4-5]。
早期通過腦電信號監測麻醉深度主要是依靠腦電波形的時域特征, 隨著快速傅里葉變換(fast fourier transform, FFT)理論的不斷發展，越來越多的腦電頻域特征被用來反映麻醉深度，腦電雙頻指數(bispectral index, BIS)就是其中之一。但迄今為止尚未尋找到對不同藥物、不同麻醉方法都具有普遍性的，且適用于各類麻醉手術的深度實時監測的特征指標[6]。人體的動力學特性是混沌的，而熵和復雜度都是非線性動力學分析方法, 比其他生理參數能挖掘出更多的信息，可以反映大腦在麻醉過程中清醒和麻醉狀態的不同。所以，熵和復雜度分析方法都非常適合處理腦電信號[7-8]。
本文在小波變換和ICA的基礎上，根據噪聲的獨立性和在頻域中分布范圍不同的特點，選擇多種算法相結合的方式進行有效地濾除處理。并利用小波熵(wavelet entropy, WE)、近似熵(approximate entropy, ApEn)和復雜度這三種特征參數對清醒與麻醉腦電狀態進行分析。目的是尋找多指標、多方法的綜合性特征參數，建立麻醉深度指標的數學模型。
1 腦電信號的預處理
1.1 基于小波變換去噪
通過小波基把任意L2(R)空間中的函數f(t)展開，其連續小波變換定義為
| $ \begin{align} & W{{T}_{f}}(a, b)=\left\langle f\left(t \right){{\psi }_{a, b}}\left(t \right)\right\rangle=\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\left| a \right|}^{-1/2}}\int\limits_{R}{f\left(t \right)\text{ }\!\!\times\!\!\text{ }\psi \left(\frac{\overline{t-b}}{a} \right.}\text{d}t \\ \end{align} $ |
其中，變量a、b分別是尺度因子和平移因子，且a≠0。將其離散化，取，, j, k∈Z, Z是自然數，離散小波變換定義為
| $ W{{T}_{f}}(j, k)=\int\limits_{R}{f\left(t \right)\text{ }\!\!\times\!\!\text{ }\overline{{{\psi }_{j, k}}\left(t \right)}}\text{d}t $ |
(1) 小波分解與重構法
對腦電信號在不同尺度j下進行分解，不同的a、b值對應不同的帶寬和中心頻率, 從而將整個頻帶劃分為一個個子頻帶，進而在不影響有用信號的前提下將某些干擾信號的頻帶置零。例如用來處理基線漂移噪聲，其頻率在0~0.5 Hz范圍內，根據采樣定理，認為信號的最高頻率是采樣頻率fs的一半，根據Mallat算法，將信號進行L級分解后整個頻帶就被分解成L+1個的子頻帶，從而可將其中所在頻帶置零，達到消除基線漂移噪聲的目的[3]。
(2) 小波變換閾值法
針對的是疊加了高斯白噪聲的腦電信號，根據白噪聲在任何正交小波變換中仍然是白噪聲且有著相同的幅度，而信號的小波系數值必然大于那些能量分散且幅值較小的噪聲的小波系數值。因此，選擇合適的閾值，對小波系數進行軟閾值或者硬閾值處理，就可以從被白噪聲污染的腦電信號中恢復出原始信號[3-4]。
(3) 小波變換模極大法
針對的是腦電信號中的奇異點，即突變點。Lip指數是表征信號局部奇異點特征的一種量度，即在尺度s下，?t∈δt0, 若有｜Wf(s, t)｜≤｜Wf(s, t0)｜，則稱t0為小波變換在尺度s下的局部模極大值點, 信號f(t)的Lip指數與小波變換模極大值滿足:
| $ \text{log}_{2}^{{{W}_{2}}jf(t)}\le \text{log}_{2}^{k}+ja $ |
即信號和噪聲的小波系數a在不同尺度j上具有不同的傳播特性。隨著尺度的增大，噪聲所對應的模極大值迅速衰減。信號的模極大值出現三種情況：對緩變信號，模極大值逐漸增大；對階躍信號，模極大值保持不變；對脈沖信號，所對應的正、負極值組成的脈沖對的幅值將同時變小。因此，觀察不同尺度間小波變換模極大值變化的規律，剔除那些由噪聲引起的模極大值并重構信號，實現去噪[9]。
(4) 小波熵法
把腦電信號在尺度j下分解的每一個高頻細節作為一個單獨的信號，將其中的高頻小波系數先分成幾個相等的小區間并計算對應的小波熵，然后根據小波熵進行最優閾值的確定，最后利用最高一層小波分解的低頻系數和經過閾值量化的不同尺度的高頻小波系數，組成重構信號所需的系數[10]。
本文就是在小波變換的基礎上，選擇了這四種算法，對采集的原始腦電信號中存在的低頻和高頻噪聲進行有效地去噪處理，并從去噪效果的評價來進行算法的結合。
1.2 基于獨立分量分析法去噪
ICA需要源信號滿足獨立性和非高斯性，而腦電信號與瞬時、尖銳的噪聲之間滿足了觀察數目m不小于源信號數目n，且源信號S 的各分量之間是相互統計獨立且最多只允許有一個是高斯分布的條件[2, 5]。
假設源信號S=[s1, s2, …, sn]T和觀測信號X=[x1, x2, …, xn]T之間的線性組合關系為X=AS，其中，A為非奇異混合，列滿秩的矩陣。過程就是不知道混合矩陣A的情況下, 尋找線性映射分離矩陣W，從觀測信號X中提取源信號S的估計值Y=WX。
本文采用的基于負熵最大的FastICA算法，是一種在每一步迭代中有大量樣本數據參與的快速尋優迭代算法，可以實現順序地提取獨立源，使得收斂更加快速、穩健。這樣，在經過小波變換去噪后，可以繼續進行ICA去噪，去除可能存在的尖銳噪聲。
1.3 基于諧波小波變換的波形預處理
諧波小波變換(harmonic wavelet transform, HWT)的定義為
| $ \begin{align} & h\text{ }(m, n, k)=(n-m)\int\limits_{-\infty }^{+\infty }{\text{f}}(t)W_{m.n.k}^{*}\left(t-\frac{k}{n-m} \right)\text{d}t, \\ & k\in Z, m, n\in \boldsymbol{\text{R} } \\ \end{align} $ |
式中，m, n決定諧波小波變換的尺度，信號的分解層數為s時，則第s層對應2s個帶寬為的子頻帶，其中，fn是乃奎斯特頻率，且。
在式(4)中:
| $ \frac{\exp \left\{ i2n\pi \left[t-\frac{k}{n-m} \right] \right\}-\exp \left\{ i2m\pi \left[t-\frac{k}{n-m} \right] \right\}}{i2\left(n-m \right)\pi \left[t-\frac{k}{n-m} \right]} $ |
其在頻域的表達式為
| $ {{W}_{m.n.k}}\left(w \right)=\left\{ \begin{align} & \frac{1}{\left(n-m \right)}\exp \left(-iw\frac{k}{n-m} \right), 2m\pi \le w\pi 2n\pi \\ & 0, \text{elsewhere} \\ \end{align} \right. $ |
HWT的特點是通過FFT和IFFT實現運算，從式(6)中看出，諧波小波只在相應的頻域段內具有恒定的幅值，而在頻段外都為0，具有嚴格的“盒形”帶通濾波特性。所以，通過變化m、n的值來調節帶寬大小和中心頻率，以匹配不同頻帶的信號，可以實現將信號互不交疊而又無遺漏地分解到相互獨立的頻段中。這個特點很適合刻畫腦電信號的特征，并從劃分出的α、β、δ、θ波來進一步反映去噪效果[11]。
2 腦電信號特征提取
2.1 小波熵
腦電信號經二進正交小波變換，在第j分解尺度k時刻的高頻分量系數為cDj, k，低頻分量系數為cAj, k，進行單支分量重構后
| $ f\left(t \right)={{\text{D}}_{1}}+{{\text{D}}_{2}}+{{\text{A}}_{2}}=\sum\limits_{j=1}^{M}{{{\text{D}}_{j}}+{{A}_{M}}} $ |
根據正交小波變換的特性，即在某時間窗內信號總功率等于各個分量功率之和。所以設Pj=，∑jpj =1, 則WE定義為
| $ \text{WE}=-\sum\limits_{j}{{{p}_{j}}\cdot \text{ln}{{\text{p}}_{j}}} $ |
其中，poj由各分量的系數cDj, k和cAj, k，或者由重建小波分量Dj和AM求得[12]。
2.2 復雜度分析
采用Lempel—Ziv算法，即不斷比較某一字符串是否是另一字符串的子串，若是，則復雜度維持不變，否則加1。通過這種比較和添加的方式，不僅可以表示一個腦電時間序列隨其長度的增長而出現新模式的速率，而且算法簡單而且實時性好, 不依賴于實驗對象的差異[13-14]。
首先采用二值化的方法對腦電信號進行粗粒化處理，即重構，形成“0-1”序列；然后對重構序列按Lempel—Ziv算法進行比較和添加，得到各個子串的界定，并計算復雜度d(n)；最后進行歸一化。當n→∞時，幾乎所有的d(n)都會趨向于一定值，即
| $ \underset{n\to \infty }{\mathop{\text{lim}}}\, d\left(n \right)=b\left(n \right)=n/lo{{g}_{a}}\left(n \right) $ |
其中，a為字符串中不同種字符的個數。
用C(n)來測度時間序列的復雜性變化，用b(n)對d(n)進行歸一化，即:
| $ C(n)=d(n)/b(n) $ |
2.3 近似熵
ApEn分析無需對時間序列粗粒化且只需很短的數據即可達到穩定的值，具有一些其它復雜性參數所不具備的特點，表征了腦電信號中產生新模式的概率，概率越大，序列的復雜性越大，ApEn值越大[15-16]。
將序列X(i)按照順序組成m維矢量(m為模式維數)，并計算其與其它矢量X(j)之間的距離
| $ \begin{align} & d[X(i)-X(j)]=\underset{k=0, \cdots, m-1}{\mathop{max}}\, \left[\left| x\left(i=k \right)-x\left(j+k \right)\right| \right], \\ & \left(i, j=1, 2, \cdots, N-m+1 \right)\\ \end{align} $ |
根據d≤r(r為相似容限)，統計每個i的個數Nm(i)與總的個數(N-m+1)的比值Crm(i)，然后對每個Cim(i)取對數并對所有i求平均，即fm(r)=。將維數增加到m+1，重復上述步驟，計算出Crm+1(i)和fm+1(r)。由于數據點數N是有限的，所以ApEn只是一個近似值，記為
| $ \begin{align} & \text{ApEn}\left(m, r, N \right)=\underset{n\to \infty }{\mathop{\text{lim}}}\, \left[{{f}^{m}}\left(r \right)-{{f}^{m+1}}\left(r \right)\right]\approx \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{f}^{m}}\left(r \right)-{{f}^{m+1}}\left(r \right)\\ \end{align} $ |
其中，m一般情況下取2，r取0.1~0.25 SD，SD為原始數據的標準差，即
| $ \text{SD}=\sqrt{\frac{1}{N-1}\sum\limits_{n-1}^{N}{{{\left[x\left(n \right)-\frac{1}{N}\sum\limits_{n-1}^{N}{x\left(n \right)} \right]}^{2}}}} $ |
ApEn分析對偶爾產生的瞬態強干擾具有良好的屏蔽作用，因為其勢必造成距離大于給定的相似容限距離r，在距離檢測時被忽略，從而對數據段的計算具有很強的穩定性。
3 結果與分析
3.1 腦電數據獲取
腦電信號的采集、記錄采用上海諾誠公司的8道腦電記錄系統，采樣頻率為128 Hz，采樣精度12位。實驗中采用雙極導聯方式，取作用電極為Fp1和Fp2，并以前額正中間為GND，頭頂為REF，左右耳垂為A1、A2。在采集中，各導聯的阻抗應小于5 kΩ。采集對象是上海市第五人民醫院麻醉科提供的全麻手術，年齡在30~50歲之間的男女患者共7例，以及同時采集相同數量的清醒腦電數據做對比。在提取麻醉腦電數據的同時，還同步記錄了心電呼吸、體溫等生理參數，以及過程中出現的藥物、體位和手術操作等變化，作為判定麻醉深度參考依據。
3.2 去噪及分類實現
截取含4 096個采樣點的數據段，通過信噪比(signal noise ratio, SNR)和均方根誤差(root mean square error, RMSE)來判斷腦電信號的去噪效果，即信號信噪比越高, 均方根誤差越小, 則估計信號就越接近于原始信號, 去噪效果越好。
| $ \begin{align} & \text{SNR}=10\log \left\{ \frac{\sum\limits_{n}{{{f}^{2}}\left(t \right)}}{\sum\limits_{n}{{{\left[f\left(t \right)-\hat{f}\left(t \right)\right]}^{2}}}} \right\} \\ & \text{RMSE=}\sqrt{\frac{1}{n}\sum\limits_{n}{{{\left[f\left(t \right)-\hat{f}\left(t \right)\right]}^{2}}}} \\ \end{align} $ |
觀察表 1，其中在數據源中，A1~A7是麻醉數據，NA1~NA7是非麻醉數據。可以判別出麻醉與非麻醉腦電信號所含噪聲成分的比例是不一樣的，導致對應去噪方式的效果區別很大。從中說明小波閾值法可以對腦電信號進行整體去噪，特別是白噪聲之類的存在于整個頻率段內的干擾，而其余的去噪算法則對特定噪聲有比較明顯的作用。在麻醉狀態下，噪聲成分區分比較明顯，其中，小波模極大值法與小波熵法的效果比較一致，由于小波模極大值法計算量大，信號重構復雜，所以本文選擇了小波熵法。與此同時，對于腦電信號中存在的眼動和肌電噪聲，本文采用了ICA法濾除(見圖 1、2)。
圖1
去噪時域圖
(a)原始腦電；(b)小波熵去噪；(c)小波分解與重構法去噪；(d)小波閾值法去噪；(e)ICA法去噪
Figure1. Denoising results of time domain(a) raw EEG; (b) wavelet entropy denoising; (c) wavelet denoising; (d) wavelet threshold de-noising; (e) ICA denoising
圖2
去噪頻域圖
(a)原始腦電；(b)小波熵去噪；(c)小波分解與重構法去噪；(d)小波閾值法去噪；(e)ICA法去噪
Figure2. Denoising results of frequency domain(a) raw EEG; (b) wavelet entropy denoising; (c) wavelet denoising; (d) wavelet threshold de-noising; (e) ICA denoising
表1
去噪性能比較
Table1.
De-noising performance
另外，根據采樣頻率fs和腦電一般劃分為4個基本不重疊的特征波段α、β、δ、θ，利用諧波小波對各個頻帶進行劃分后，從圖 3可以看出劃分出來的特征波十分明顯，進一步反映了之前去噪處理的有效性。
圖3
諧波小波分解結果圖
(a)提取的α、β、δ、θ特征波時域圖; (b)提取的α、β、δ、θ特征波頻譜圖
Figure3. Harmonic wavelet decomposition results(a) extraction ofα, β, δ, θcharacteristic wave time domain; (b) extraction ofα, β, δ, θcharacteristic wave spectrum
將預處理后的腦電信號進行WE、ApEn和KC三種特征向量的計算，從表 2可以看出，清醒狀態下的WE、ApEn和KC值與麻醉狀態下的比，區分明顯。這三個特征參數可以將兩種狀態給區分開來，從而可以用來判斷麻醉中的狀態改變。
表2
兩組狀態下的WE、ApEn和KC值
Table2.
WE, ApEn, and KC values under two sets of state
其中，通過觀察變化率
| $ \begin{align} & \Delta \text{WE}=\frac{\text{0}\text{.166 8-0}\text{.082 6}}{\text{0}\text{.166 8}}=0.505 \\ & \Delta \text{ApEn}=\frac{\text{0}\text{.425 9-0}\text{.333 8}}{\text{0}\text{.425 9}}=0.216 \\ & \Delta \text{KC}=\frac{\text{0}\text{.223 1-0}\text{.179}6}{\text{0}\text{.223 1}}=0.195 \\ \end{align} $ |
WE的狀態變化率可以達到50%，而ApEn和KC兩者在20%上下，說明WE對腦電信號狀態變化的敏感性優于ApEn和KC。三者中，WE最適合應用于麻醉深度指數量化的研究中，而且算法簡便，數據量短，實時性強。比較ApEn與KC，雖然判斷性能比較接近，但ApEn算法復雜，計算量大，實時性較差。
4 結語
本文通過小波變換與ICA法的結合，實現了原始腦電信號的有效去噪，并利用諧波小波包分類、提取了腦電信號中的4個基本特征波，從結果上進一步反映了去噪效果的可靠性。預處理后，利用復雜度、近似熵、小波熵這三種特征算法來計算腦電數據段，區分出了兩種狀態下的腦電信號，達到了預期效果。因此，本文為麻醉深度量化的研究提供了參考依據。
引言
麻醉是通過抑制中樞神經系統而達到意識消失、止痛的目的，在外科手術中是不可或缺的。精確、實時的麻醉深度監護對于提高麻醉質量、保障手術安全，以及減少麻醉并發癥和控制麻醉藥品用量都具有重要的意義。腦電信號是腦內部億萬個神經元興奮或抑制而記錄在頭皮表面的電生理信號，它能夠直接反映大腦意識的興奮或抑制程度，對其進行處理和分析，無論是在臨床上對腦疾病的診斷和治療，還是對腦認知功能的研究都是十分重要的[1-2]。因此，利用EEG信號分析大腦受抑程度是最為直接的，并可作為一種非插入式的監測手段，從而來反饋、改變手術過程中麻醉藥物的使用。
在腦電信號采集過程中，不可避免地會受到各種干擾，例如工頻干擾、基線漂移、眼動和肌電等。而且腦電信號本身是非線性、非平穩且十分微弱的，狀態不同，噪聲的成分和類型也不同，導致了信號被淹沒其中。均值法、中值法等傳統的去噪算法只是對腦電信號在空間域上進行處理，自適應差，局部性強，結果不能令人滿意。小波變換具有多尺度分析的特點，能夠聚焦到信號的任意細節進行多分辨率的時頻分析，對處理像腦電這種非平穩信號的干擾是非常適合的[3]。獨立分量分析(independent component analysis, ICA)方法需要滿足一定分離條件，即干擾信號與腦電信號之間是統計獨立的，從而將腦電信號按照統計獨立的原則通過優化算法分解為若干獨立來源成分，這樣就將去噪問題轉化為獨立分量分離問題，特別適合處理那些突發的，尖銳的噪聲。所以，采用兩者的結合，可以有效去除低頻和尖峰等高頻噪聲[4-5]。
早期通過腦電信號監測麻醉深度主要是依靠腦電波形的時域特征, 隨著快速傅里葉變換(fast fourier transform, FFT)理論的不斷發展，越來越多的腦電頻域特征被用來反映麻醉深度，腦電雙頻指數(bispectral index, BIS)就是其中之一。但迄今為止尚未尋找到對不同藥物、不同麻醉方法都具有普遍性的，且適用于各類麻醉手術的深度實時監測的特征指標[6]。人體的動力學特性是混沌的，而熵和復雜度都是非線性動力學分析方法, 比其他生理參數能挖掘出更多的信息，可以反映大腦在麻醉過程中清醒和麻醉狀態的不同。所以，熵和復雜度分析方法都非常適合處理腦電信號[7-8]。
本文在小波變換和ICA的基礎上，根據噪聲的獨立性和在頻域中分布范圍不同的特點，選擇多種算法相結合的方式進行有效地濾除處理。并利用小波熵(wavelet entropy, WE)、近似熵(approximate entropy, ApEn)和復雜度這三種特征參數對清醒與麻醉腦電狀態進行分析。目的是尋找多指標、多方法的綜合性特征參數，建立麻醉深度指標的數學模型。
1 腦電信號的預處理
1.1 基于小波變換去噪
通過小波基把任意L2(R)空間中的函數f(t)展開，其連續小波變換定義為
| $ \begin{align} & W{{T}_{f}}(a, b)=\left\langle f\left(t \right){{\psi }_{a, b}}\left(t \right)\right\rangle=\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\left| a \right|}^{-1/2}}\int\limits_{R}{f\left(t \right)\text{ }\!\!\times\!\!\text{ }\psi \left(\frac{\overline{t-b}}{a} \right.}\text{d}t \\ \end{align} $ |
其中，變量a、b分別是尺度因子和平移因子，且a≠0。將其離散化，取，, j, k∈Z, Z是自然數，離散小波變換定義為
| $ W{{T}_{f}}(j, k)=\int\limits_{R}{f\left(t \right)\text{ }\!\!\times\!\!\text{ }\overline{{{\psi }_{j, k}}\left(t \right)}}\text{d}t $ |
(1) 小波分解與重構法
對腦電信號在不同尺度j下進行分解，不同的a、b值對應不同的帶寬和中心頻率, 從而將整個頻帶劃分為一個個子頻帶，進而在不影響有用信號的前提下將某些干擾信號的頻帶置零。例如用來處理基線漂移噪聲，其頻率在0~0.5 Hz范圍內，根據采樣定理，認為信號的最高頻率是采樣頻率fs的一半，根據Mallat算法，將信號進行L級分解后整個頻帶就被分解成L+1個的子頻帶，從而可將其中所在頻帶置零，達到消除基線漂移噪聲的目的[3]。
(2) 小波變換閾值法
針對的是疊加了高斯白噪聲的腦電信號，根據白噪聲在任何正交小波變換中仍然是白噪聲且有著相同的幅度，而信號的小波系數值必然大于那些能量分散且幅值較小的噪聲的小波系數值。因此，選擇合適的閾值，對小波系數進行軟閾值或者硬閾值處理，就可以從被白噪聲污染的腦電信號中恢復出原始信號[3-4]。
(3) 小波變換模極大法
針對的是腦電信號中的奇異點，即突變點。Lip指數是表征信號局部奇異點特征的一種量度，即在尺度s下，?t∈δt0, 若有｜Wf(s, t)｜≤｜Wf(s, t0)｜，則稱t0為小波變換在尺度s下的局部模極大值點, 信號f(t)的Lip指數與小波變換模極大值滿足:
| $ \text{log}_{2}^{{{W}_{2}}jf(t)}\le \text{log}_{2}^{k}+ja $ |
即信號和噪聲的小波系數a在不同尺度j上具有不同的傳播特性。隨著尺度的增大，噪聲所對應的模極大值迅速衰減。信號的模極大值出現三種情況：對緩變信號，模極大值逐漸增大；對階躍信號，模極大值保持不變；對脈沖信號，所對應的正、負極值組成的脈沖對的幅值將同時變小。因此，觀察不同尺度間小波變換模極大值變化的規律，剔除那些由噪聲引起的模極大值并重構信號，實現去噪[9]。
(4) 小波熵法
把腦電信號在尺度j下分解的每一個高頻細節作為一個單獨的信號，將其中的高頻小波系數先分成幾個相等的小區間并計算對應的小波熵，然后根據小波熵進行最優閾值的確定，最后利用最高一層小波分解的低頻系數和經過閾值量化的不同尺度的高頻小波系數，組成重構信號所需的系數[10]。
本文就是在小波變換的基礎上，選擇了這四種算法，對采集的原始腦電信號中存在的低頻和高頻噪聲進行有效地去噪處理，并從去噪效果的評價來進行算法的結合。
1.2 基于獨立分量分析法去噪
ICA需要源信號滿足獨立性和非高斯性，而腦電信號與瞬時、尖銳的噪聲之間滿足了觀察數目m不小于源信號數目n，且源信號S 的各分量之間是相互統計獨立且最多只允許有一個是高斯分布的條件[2, 5]。
假設源信號S=[s1, s2, …, sn]T和觀測信號X=[x1, x2, …, xn]T之間的線性組合關系為X=AS，其中，A為非奇異混合，列滿秩的矩陣。過程就是不知道混合矩陣A的情況下, 尋找線性映射分離矩陣W，從觀測信號X中提取源信號S的估計值Y=WX。
本文采用的基于負熵最大的FastICA算法，是一種在每一步迭代中有大量樣本數據參與的快速尋優迭代算法，可以實現順序地提取獨立源，使得收斂更加快速、穩健。這樣，在經過小波變換去噪后，可以繼續進行ICA去噪，去除可能存在的尖銳噪聲。
1.3 基于諧波小波變換的波形預處理
諧波小波變換(harmonic wavelet transform, HWT)的定義為
| $ \begin{align} & h\text{ }(m, n, k)=(n-m)\int\limits_{-\infty }^{+\infty }{\text{f}}(t)W_{m.n.k}^{*}\left(t-\frac{k}{n-m} \right)\text{d}t, \\ & k\in Z, m, n\in \boldsymbol{\text{R} } \\ \end{align} $ |
式中，m, n決定諧波小波變換的尺度，信號的分解層數為s時，則第s層對應2s個帶寬為的子頻帶，其中，fn是乃奎斯特頻率，且。
在式(4)中:
| $ \frac{\exp \left\{ i2n\pi \left[t-\frac{k}{n-m} \right] \right\}-\exp \left\{ i2m\pi \left[t-\frac{k}{n-m} \right] \right\}}{i2\left(n-m \right)\pi \left[t-\frac{k}{n-m} \right]} $ |
其在頻域的表達式為
| $ {{W}_{m.n.k}}\left(w \right)=\left\{ \begin{align} & \frac{1}{\left(n-m \right)}\exp \left(-iw\frac{k}{n-m} \right), 2m\pi \le w\pi 2n\pi \\ & 0, \text{elsewhere} \\ \end{align} \right. $ |
HWT的特點是通過FFT和IFFT實現運算，從式(6)中看出，諧波小波只在相應的頻域段內具有恒定的幅值，而在頻段外都為0，具有嚴格的“盒形”帶通濾波特性。所以，通過變化m、n的值來調節帶寬大小和中心頻率，以匹配不同頻帶的信號，可以實現將信號互不交疊而又無遺漏地分解到相互獨立的頻段中。這個特點很適合刻畫腦電信號的特征，并從劃分出的α、β、δ、θ波來進一步反映去噪效果[11]。
2 腦電信號特征提取
2.1 小波熵
腦電信號經二進正交小波變換，在第j分解尺度k時刻的高頻分量系數為cDj, k，低頻分量系數為cAj, k，進行單支分量重構后
| $ f\left(t \right)={{\text{D}}_{1}}+{{\text{D}}_{2}}+{{\text{A}}_{2}}=\sum\limits_{j=1}^{M}{{{\text{D}}_{j}}+{{A}_{M}}} $ |
根據正交小波變換的特性，即在某時間窗內信號總功率等于各個分量功率之和。所以設Pj=，∑jpj =1, 則WE定義為
| $ \text{WE}=-\sum\limits_{j}{{{p}_{j}}\cdot \text{ln}{{\text{p}}_{j}}} $ |
其中，poj由各分量的系數cDj, k和cAj, k，或者由重建小波分量Dj和AM求得[12]。
2.2 復雜度分析
采用Lempel—Ziv算法，即不斷比較某一字符串是否是另一字符串的子串，若是，則復雜度維持不變，否則加1。通過這種比較和添加的方式，不僅可以表示一個腦電時間序列隨其長度的增長而出現新模式的速率，而且算法簡單而且實時性好, 不依賴于實驗對象的差異[13-14]。
首先采用二值化的方法對腦電信號進行粗粒化處理，即重構，形成“0-1”序列；然后對重構序列按Lempel—Ziv算法進行比較和添加，得到各個子串的界定，并計算復雜度d(n)；最后進行歸一化。當n→∞時，幾乎所有的d(n)都會趨向于一定值，即
| $ \underset{n\to \infty }{\mathop{\text{lim}}}\, d\left(n \right)=b\left(n \right)=n/lo{{g}_{a}}\left(n \right) $ |
其中，a為字符串中不同種字符的個數。
用C(n)來測度時間序列的復雜性變化，用b(n)對d(n)進行歸一化，即:
| $ C(n)=d(n)/b(n) $ |
2.3 近似熵
ApEn分析無需對時間序列粗粒化且只需很短的數據即可達到穩定的值，具有一些其它復雜性參數所不具備的特點，表征了腦電信號中產生新模式的概率，概率越大，序列的復雜性越大，ApEn值越大[15-16]。
將序列X(i)按照順序組成m維矢量(m為模式維數)，并計算其與其它矢量X(j)之間的距離
| $ \begin{align} & d[X(i)-X(j)]=\underset{k=0, \cdots, m-1}{\mathop{max}}\, \left[\left| x\left(i=k \right)-x\left(j+k \right)\right| \right], \\ & \left(i, j=1, 2, \cdots, N-m+1 \right)\\ \end{align} $ |
根據d≤r(r為相似容限)，統計每個i的個數Nm(i)與總的個數(N-m+1)的比值Crm(i)，然后對每個Cim(i)取對數并對所有i求平均，即fm(r)=。將維數增加到m+1，重復上述步驟，計算出Crm+1(i)和fm+1(r)。由于數據點數N是有限的，所以ApEn只是一個近似值，記為
| $ \begin{align} & \text{ApEn}\left(m, r, N \right)=\underset{n\to \infty }{\mathop{\text{lim}}}\, \left[{{f}^{m}}\left(r \right)-{{f}^{m+1}}\left(r \right)\right]\approx \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{f}^{m}}\left(r \right)-{{f}^{m+1}}\left(r \right)\\ \end{align} $ |
其中，m一般情況下取2，r取0.1~0.25 SD，SD為原始數據的標準差，即
| $ \text{SD}=\sqrt{\frac{1}{N-1}\sum\limits_{n-1}^{N}{{{\left[x\left(n \right)-\frac{1}{N}\sum\limits_{n-1}^{N}{x\left(n \right)} \right]}^{2}}}} $ |
ApEn分析對偶爾產生的瞬態強干擾具有良好的屏蔽作用，因為其勢必造成距離大于給定的相似容限距離r，在距離檢測時被忽略，從而對數據段的計算具有很強的穩定性。
3 結果與分析
3.1 腦電數據獲取
腦電信號的采集、記錄采用上海諾誠公司的8道腦電記錄系統，采樣頻率為128 Hz，采樣精度12位。實驗中采用雙極導聯方式，取作用電極為Fp1和Fp2，并以前額正中間為GND，頭頂為REF，左右耳垂為A1、A2。在采集中，各導聯的阻抗應小于5 kΩ。采集對象是上海市第五人民醫院麻醉科提供的全麻手術，年齡在30~50歲之間的男女患者共7例，以及同時采集相同數量的清醒腦電數據做對比。在提取麻醉腦電數據的同時，還同步記錄了心電呼吸、體溫等生理參數，以及過程中出現的藥物、體位和手術操作等變化，作為判定麻醉深度參考依據。
3.2 去噪及分類實現
截取含4 096個采樣點的數據段，通過信噪比(signal noise ratio, SNR)和均方根誤差(root mean square error, RMSE)來判斷腦電信號的去噪效果，即信號信噪比越高, 均方根誤差越小, 則估計信號就越接近于原始信號, 去噪效果越好。
| $ \begin{align} & \text{SNR}=10\log \left\{ \frac{\sum\limits_{n}{{{f}^{2}}\left(t \right)}}{\sum\limits_{n}{{{\left[f\left(t \right)-\hat{f}\left(t \right)\right]}^{2}}}} \right\} \\ & \text{RMSE=}\sqrt{\frac{1}{n}\sum\limits_{n}{{{\left[f\left(t \right)-\hat{f}\left(t \right)\right]}^{2}}}} \\ \end{align} $ |
觀察表 1，其中在數據源中，A1~A7是麻醉數據，NA1~NA7是非麻醉數據。可以判別出麻醉與非麻醉腦電信號所含噪聲成分的比例是不一樣的，導致對應去噪方式的效果區別很大。從中說明小波閾值法可以對腦電信號進行整體去噪，特別是白噪聲之類的存在于整個頻率段內的干擾，而其余的去噪算法則對特定噪聲有比較明顯的作用。在麻醉狀態下，噪聲成分區分比較明顯，其中，小波模極大值法與小波熵法的效果比較一致，由于小波模極大值法計算量大，信號重構復雜，所以本文選擇了小波熵法。與此同時，對于腦電信號中存在的眼動和肌電噪聲，本文采用了ICA法濾除(見圖 1、2)。
圖1
去噪時域圖
(a)原始腦電；(b)小波熵去噪；(c)小波分解與重構法去噪；(d)小波閾值法去噪；(e)ICA法去噪
Figure1. Denoising results of time domain(a) raw EEG; (b) wavelet entropy denoising; (c) wavelet denoising; (d) wavelet threshold de-noising; (e) ICA denoising
圖2
去噪頻域圖
(a)原始腦電；(b)小波熵去噪；(c)小波分解與重構法去噪；(d)小波閾值法去噪；(e)ICA法去噪
Figure2. Denoising results of frequency domain(a) raw EEG; (b) wavelet entropy denoising; (c) wavelet denoising; (d) wavelet threshold de-noising; (e) ICA denoising
表1
去噪性能比較
Table1.
De-noising performance
另外，根據采樣頻率fs和腦電一般劃分為4個基本不重疊的特征波段α、β、δ、θ，利用諧波小波對各個頻帶進行劃分后，從圖 3可以看出劃分出來的特征波十分明顯，進一步反映了之前去噪處理的有效性。
圖3
諧波小波分解結果圖
(a)提取的α、β、δ、θ特征波時域圖; (b)提取的α、β、δ、θ特征波頻譜圖
Figure3. Harmonic wavelet decomposition results(a) extraction ofα, β, δ, θcharacteristic wave time domain; (b) extraction ofα, β, δ, θcharacteristic wave spectrum
將預處理后的腦電信號進行WE、ApEn和KC三種特征向量的計算，從表 2可以看出，清醒狀態下的WE、ApEn和KC值與麻醉狀態下的比，區分明顯。這三個特征參數可以將兩種狀態給區分開來，從而可以用來判斷麻醉中的狀態改變。
表2
兩組狀態下的WE、ApEn和KC值
Table2.
WE, ApEn, and KC values under two sets of state
其中，通過觀察變化率
| $ \begin{align} & \Delta \text{WE}=\frac{\text{0}\text{.166 8-0}\text{.082 6}}{\text{0}\text{.166 8}}=0.505 \\ & \Delta \text{ApEn}=\frac{\text{0}\text{.425 9-0}\text{.333 8}}{\text{0}\text{.425 9}}=0.216 \\ & \Delta \text{KC}=\frac{\text{0}\text{.223 1-0}\text{.179}6}{\text{0}\text{.223 1}}=0.195 \\ \end{align} $ |
WE的狀態變化率可以達到50%，而ApEn和KC兩者在20%上下，說明WE對腦電信號狀態變化的敏感性優于ApEn和KC。三者中，WE最適合應用于麻醉深度指數量化的研究中，而且算法簡便，數據量短，實時性強。比較ApEn與KC，雖然判斷性能比較接近，但ApEn算法復雜，計算量大，實時性較差。
4 結語
本文通過小波變換與ICA法的結合，實現了原始腦電信號的有效去噪，并利用諧波小波包分類、提取了腦電信號中的4個基本特征波，從結果上進一步反映了去噪效果的可靠性。預處理后，利用復雜度、近似熵、小波熵這三種特征算法來計算腦電數據段，區分出了兩種狀態下的腦電信號，達到了預期效果。因此，本文為麻醉深度量化的研究提供了參考依據。
