全身麻醉是外科手術中保證患者安全的必不可少的部分，腦電圖（EEG）能反映大腦活動狀況，包含豐富的信息，因此已廣泛應用于監測麻醉深度。本文提出了一種將小波變換與人工神經網絡（ANN）相結合的方法來估計麻醉深度。利用離散小波變換（DWT）將腦電信號進行分解，根據分解得到的近似系數與細節系數計算 9 種特征參數，并對這 9 種特征參數進行克魯斯卡爾-沃利斯統計檢驗，結果表明這 9 種特征參數在清醒、輕度麻醉、中度麻醉和深度麻醉這四種不同麻醉水平間的差異均有統計學意義（P < 0.001）。將這 9 種特征參數作為 ANN 的輸入，以雙譜指數（BIS）作為參考輸出，使用 8 例全麻手術的患者數據對該方法進行了評估。該方法在 7∶3 留出法中對測試集四種不同麻醉水平的分類準確度為 85.98%，與 BIS 的相關系數為 0.977 0。結果表明，該方法能較好地區分四種不同麻醉水平，對于麻醉深度監測具有廣闊的應用前景。
引用本文： 袁思念, 葉繼倫, 張旭, 周晶晶, 檀雪, 李若薇, 鄧鑄強, 丁耀茂. 基于小波變換與人工神經網絡的麻醉深度計算方法研究. 生物醫學工程學雜志, 2021, 38(5): 838-847. doi: 10.7507/1001-5515.202007003 復制
引言
在手術過程中，全身麻醉是保證患者安全的必要和重要的環節，麻醉劑量過大可能會延長患者的恢復時間，而麻醉不足可能會導致患者產生術中記憶甚至在術中蘇醒，給患者帶來生理或者心理上的創傷[1]，所以在手術過程中進行準確的無創麻醉深度監測是非常必要的。由于麻醉藥物主要作用于中樞神經系統，在麻醉過程中，大腦活動的變化表現為腦電圖（electroencephalogram，EEG）波形或信號的振幅和頻率變化[2]，處理麻醉過程中的 EEG 信號有助于監測患者的麻醉深度。
美國 Aspect 公司于 1992 年開發的雙譜指數（bispectral index，BIS）監護儀由于具備實用性和有效性，獲得了市場的廣泛認可，是麻醉深度監測的可靠比較標準[3]，目前，臨床上主要采用基于 BIS 監護儀的麻醉深度評價方法。雖然 BIS 是一種已經應用于臨床數十年的指標，在降低術中意識發生率和指導麻醉藥物使用方面都取得了一定的成功，但是其算法并未公開。國內麻醉深度監護儀的需求較大，但是目前市場上可供選擇的產品卻非常有限[3]。由于從國外引進的儀器價格昂貴，國內制造的麻醉深度監護儀沒有技術上的改進，而麻醉深度監護儀最關鍵的技術就是其算法，所以，研究麻醉深度的計算方法是具有重要意義的。
早期基于傅立葉理論的頻譜分析方法假設數據在短時間內是平穩的，由于 EEG 信號在人體生物系統中呈現出非平穩、非線性的特征，這些基于頻率的方法可能會忽略有價值的信息[4]。近年來，許多研究者也采用非線性分析方法處理 EEG 信號，如希爾伯特-黃變換[5]、復雜性分析[6]、熵[4-5]和去趨勢波動分析[7]等。這些分析雖然提供了傳統 EEG 分析無法獲得的信息，但是也存在不足之處：希爾伯特-黃變換由于其計算量大不能應用于計算機上的實時分析[3]；復雜性分析、熵等非線性動力學方法對于初始條件敏感[8]；去趨勢波動分析方法在深度麻醉狀態下的分析效果不太好[9]。
小波變換作為時頻域分析工具，能夠同時反映信號的時域和頻域特征，很好地適用于非平穩的信號的分析，因而被廣泛用于生物醫學信號的分析[10-13]。EEG 信號通常被認為是非平穩的，而小波變換適合分析非平穩非線性的信號，基于 Mallat 快速算法的離散小波變換（discrete wavelet transform，DWT）是 EEG 信號分析的常用算法。模式識別與機器學習在 EEG 信號分析中也得到了廣泛的應用，一些研究人員使用支持向量機（support vector machine，SVM）和人工神經網絡（artificial neural network，ANN）做麻醉深度的回歸與分類[1,4,13-14]，而且 ANN 在模型開發中具有強大的自適應能力、學習能力、非線性映射能力、魯棒性和容錯性的優勢[13]。
本文使用 DWT 分析了麻醉中的腦電信號，并通過分解得到的小波系數提取了 δ 波、θ 波、α 波、β 波、γ 波的相對能量比，腦電總能量，肌電能量，β 率以及小波熵（wavelet entropy，WE）9 種特征參數，使用克魯斯卡爾-沃利斯統計檢驗來確定所使用的參數對四種麻醉水平是否具有顯著差異，并結合 ANN，將 BIS 作為參考輸出來評估麻醉深度，從靈敏度和分類準確度兩個方面對該方法的性能進行了評價。此外，我們還將 ANN 與 SVM 兩種模型進行了性能比較。
1 材料和方法
1.1 數據獲取
本研究的實驗數據均來自高州人民醫院麻醉科，使用美國 Aspect medical 公司生產的 BISX 外置模塊及數據記錄系統記錄了 8 名成人手術全過程全身麻醉狀態下的 EEG 信號。麻醉藥物基本選用芬太尼、異苯酚、卡肌寧、丙泊酚等靜脈麻醉藥物，麻醉方式均為靜吸麻醉方式。采集過程中 BISX 模塊濾波器模式設置為高通 0.5 Hz、低通 70 Hz、帶阻 50&60 Hz，以 256 Hz 的采樣率記錄，總共記錄約 30.5 h。數據已獲醫院授權。
1.2 EEG 信號處理流程
考慮到采集過程中存在電刀等干擾，所以使用自適應方差閾值法去除 EEG 信號中的尖銳噪聲，然后再將數據集分割成樣本；將分割好的樣本通過 DWT 提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，得到特征集；最后將特征集輸入模型進行訓練，得到可使用的能夠較為準確地計算麻醉深度指數的模型。
1.3 腦電數據預處理
1.3.1 自適應方差閾值法去噪[15 ]
在采集手術中的 EEG 信號時，由于外界環境電刀等帶來的干擾，會出現異常腦電數據段，表現為尖銳的噪聲。為了較好地去除大的波動干擾，又不對正常 EEG 信號造成影響，本文采用自適應調整的方法確定閾值，具體操作如下：選擇當前 60 s 作為一個時間段，將當前 60 s 的腦電數據分割成 60 個 1 s 的時間段，求得 60 個方差值，再利用這 60 個方差值的均值與中值計算得到自適應閾值，閾值設計公式如下：
 |
式（1）中，NewT 是計算得到的方差閾值，Mean 與 Med 分別是最近 60 s 數據段 60 個方差的均值與中值，OldT 是上一個 60 s 數據段的閾值。每次檢測平移 1 s 時間段，并且更新閾值，實現方差閾值的自動調整。如果計算出該段數據的方差大于計算的閾值，則認定該段數據是較大的波動干擾，將其去除。
1.3.2 數據集分割
使用步長 5 s、大小為 60 s 的滑動窗，將 8 例腦電數據分割成每段 60 s 的 EEG 數據長度的樣本，每個樣本包含 15 360 個采樣點，共獲得 21 821 個樣本。整個數據集根據 BIS 分為四個級別。BIS 是根據麻醉深度顯示 0～100 范圍內的無量綱值，0 表示大腦基本沒有活動，100 表示完全清醒，范圍 0～100 可以細分為 4 個狀態，表 1 對不同麻醉狀態的 BIS 值范圍作了說明，中度麻醉狀態被認為是在全身麻醉下進行手術的合適水平[10]。
表1
不同麻醉狀態的 BIS 指數范圍與對應級別
Table1.
BIS index range and levels for different anesthesia states
圖 1 為不同麻醉階段對應的 60 s 長的 EEG 信號。由圖 1 可以看出，隨著麻醉深度的增加，EEG 信號幅度逐漸增大，頻率逐漸減小。在深度麻醉的狀態下，容易出現爆發抑制信號。
圖1
四種不同麻醉水平的 EEG 信號
Figure1.
EEG of four levels of anesthesia
1.4 使用 DWT 進行分解
給定基本小波函數：
 |
其中，變量a、b分別是尺度因子和平移因子，且
。通過小波基把任意
空間中的函數
展開，則信號的連續小波變換為：
 |
取
，
，j, 
，Z 是自然數，DWT 定義為：
 |
通過對 a、b 進行采樣，小波變換可基于 Mallat 算法進行離散化，如圖 2 所示。
圖2
Mallat 算法分解架構
Figure2.
The Mallat algorithm decomposition architecture
其中，
為高通濾波器，將濾掉信號的低頻分量而輸出高頻分量，
為低通濾波器，濾掉信號的高頻分量而輸出低頻分量，根據下式將腦電信號分解成細節系數和近似系數：
 |
 |
式中 
與 
分別表示 EEG 信號分解得到的近似系數與細節系數。
DWT 分析的兩個重要考慮因素是選擇合適的母小波和分解層數。由于 Daubechies 小波被廣泛用于與 EEG 信號相關的應用[15]，本研究采用 db4 小波為母小波，使用 DWT 將 EEG 信號分解為不同的頻帶。小波分解的水平取決于五個腦電波（δ 波、θ 波、α 波、β 波、γ 波）以及肌電圖（electromyogram，EMG）信號的頻帶范圍。根據 EEG 信號頻帶分類，選取 5 階分解階數。因此，60 s 一段的 EEG 信號被分解成五個細節系數 D1-D5 和最后一個近似系數的 A5，其分別對應的頻帶如表 2 所示。
表2
EEG 信號頻帶分解
Table2.
EEG frequency band decomposition
1.5 特征提取
在不同的麻醉深度下，EEG 信號具有不同的波形幅度與頻率分布，可以通過 DWT 分解得到的小波系數，從中提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數。其中，δ%、θ%、α%、β%、γ% 用于量化 EEG 中的頻率成分的變化，EEGTP、EMGTP 代表 EEG 與 EMG 功率大小的變化，β 率有利于放大 EEG 信號中高頻成分與低頻成分的轉移程度，WE 則從熵的角度反映信號譜能量在各個子空間分布的有序或無序程度，因此，使用這些特征參數，可以更全面地量化 EEG 信號中的變化。
小波系數的能量給出了信號強度的信息，細節系數與近似系數的能量由以下方程給出：
 |
 |
這里 
是分解級別，
是對應的小波系數下標。
（1）由于 EEG 的能量基本集中在 0～64 Hz 頻帶范圍內，所以定義 EEG 總能量為：
 |
（2）為了避免計算出的數值過大，對計算出的 EEG 總能量作對數處理，定義參數對數 EEG 能量 EEGTP 為：
 |
（3）定義參數 δ%、θ%、α%、β%、γ% 分別為：
 |
 |
 |
 |
 |
（4）由于 EMG 信號主要分布在 64～128 Hz 頻帶內，定義參數肌電能量 EMGTP 為：
 |
（5）β 率定義為高頻的 β 波能量與低頻的 α 波能量之比，β 率計算公式為：
 |
（6）WE 的計算公式如下：
① 首先計算分解得到的小波總能量：
 |
② 計算小波子帶的相對能量
 |
 |
③ 計算 WE
 |
1.6 克魯斯卡爾-沃利斯統計檢驗
由于利用 DWT 提取的腦電特征數據樣本量較少，數據不服從正態分布，所以本文采用克魯斯卡爾-沃利斯統計檢驗方法進行統計分析，以檢驗所提取的特征在清醒、輕度麻醉、中度麻醉和深度麻醉不同麻醉水平下的差異是否具有統計學意義。克魯斯卡爾-沃利斯統計檢驗方法是檢驗樣本是否來自相同分布的一種非參數方法[16]，它被用來確定兩組或多組自變量和連續或有序因變量之間的差異是否具有統計學意義，測定這些組在統計上是否不同。本文還使用了箱線圖描述每個特征在四種不同麻醉水平下的數據分布情況。
1.7 回歸與分類
ANN 是一種靈活的非參數的并行計算模型，基于模擬人腦的神經結構而發展起來[14]。ANN 通常由多個相互連接的多層節點組成，分為輸入層、隱含層和輸出層，輸入層中的每個節點接收外部信息，輸出層中的每個節點生成模型解決方案并輸出最終結果。本文使用的是誤差反向傳播神經網絡即 BP 神經網絡來輸出最終的麻醉深度指數，BP 神經網絡是一種多層前饋 ANN，也是目前應用最廣泛的 ANN 之一[1]。
SVM 是基于統計學習理論的一種機器學習方法[17]，它通過適當的非線性映射將輸入向量映射到一個高維的特征空間，使得數據總能被一個超平面分割成兩類，并尋找一個最優超平面。SVM 不僅能將兩類數據正確分開，而且使分類間隔最大，使其在有限樣本下具有了較強的泛化能力。在本研究中，我們使用 LIBSVM 工具箱完成 SVM 的相關分析。
1.8 分類性能分析
由于腦電特征數據樣本有限，我們使用固定的抽樣間隔抽取數據集的 70% 為訓練集，剩下的 30% 作為測試集來評估分類方法的性能，并且使用了 10 折交叉驗證評估了分類方法的泛化性能。根據測試集輸出的結果，利用靈敏度和分類精度可以對分類器的性能進行定量分析。靈敏度是一種麻醉狀態被正確地識別的數目與相應麻醉狀態的總數的比率。分類精度是指所有被正確識別的麻醉狀態的數目與實際麻醉狀態的總數的比率。
此外，本文還計算了 BIS 與所提出方法的指標之間的皮爾遜相關系數來評價所提出的方法，采用 Bland-Altman 分析來直觀反映兩種測量結果的一致性界限。最后，使用配對樣本t檢驗來檢驗最后 ANN 輸出的指數與 BIS 的差異有無統計學意義，判斷它們是否來自于分布相同的總體。
2 結果
2.1 自適應方差閾值法去除噪聲
方差閾值的計算方法如式（1）所述，如果發現該段數據的方差大于計算的閾值，則認定該段數據是較大的波動干擾，應將其去除。圖 3a 是一段 60 s 帶噪聲的原始 EEG 信號。圖 3b 中，藍線是計算得到的 60 s 的 EEG 數據的方差，紅線是計算得到的自適應閾值，從圖 3b 中可以看出，EEG 信號波動較大的地方計算得到的方差也較大，而通過自適應方差閾值法計算出的閾值則較為平穩，所以，只要判定其方差是否大于閾值，就可知道此數據段是否為噪聲。圖 3c 是使用自適應方差閾值法去除噪聲后的 EEG 信號，從圖 3c 中可以看出，波動較大的尖銳噪聲可以被很好地去除。
圖3
自適應方差閾值法去噪
a.帶噪聲的 EEG 信號；b.方差與閾值；c.去噪后的 EEG 信號
Figure3. Adaptive variance threshold denoising in EEGa. EEG with noise; b. variance and threshold; c. EEG after denoising
2.2 DWT 與特征提取
圖 4 顯示了 60 s 清醒狀態下的 EEG 信號 DWT 五級分解得到小波細節系數 D1～D5 與近似系數 A5 重構后的信號，從圖 4 中可以看出，DWT 能很好地將 EEG 各個頻段的信號區分開來，并且分解后的信號包含了原始 EEG 信號的全部信息。
圖4
清醒狀態下 60 s 的 EEG 信號 DWT 分解重構后的信號
Figure4.
The reconstructed signal of DWT decomposition of the 60 s EEG signal in the awake state
利用前面所述方法，對獲得的 21 821 個樣本提取特征，得到 21 821 個特征樣本，每個特征樣本包含 9 個特征。本文采用克魯斯卡爾-沃利斯統計檢驗方法檢驗提取的特征對四種麻醉狀態的鑒別能力，結果如表 3 所示，克魯斯卡爾-沃利斯統計檢驗得到的P值均小于 0.001，表明九個特征在四種不同麻醉水平下的差異具有統計學意義。
表3
克魯斯卡爾-沃利斯統計檢驗結果
Table3.
Kruskal-Wallis statistical test results
每個特征在不同麻醉深度下兩兩比較的結果中，只有兩組特征在Ⅰ vs. Ⅳ的差異沒有統計學意義：特征 θ% 的P值為 0.126 > 0.05，特征 α% 的P值為 0.170 > 0.05，剩余組的兩兩比較結果P值均小于 0.001，表明剩余組的差異均具有統計學意義。
不同麻醉深度下九種特征的箱線圖如圖 5 所示。從圖 5a～5e 中可以看出，δ% 反映的是 EEG 信號中低頻分量所占 EEG 總能量的比例，其值隨著麻醉深度的增加而增加。β% 與 γ% 在清醒階段的值較高，反映的是 EEG 信號中的高頻分量所占 EEG 總能量的比例，其值均與麻醉深度成反比。當受試者處于清醒狀態時，由于焦慮和壓力，代表高頻成分的 β 波與 α 波會更加活躍。然而，當麻醉時，代表低頻成分的 δ 波與 θ 波則變得活躍，因此，在中度麻醉和深度麻醉狀態時，δ% 的值較高。
圖5
九種特征參數在四種麻醉水平下的箱線圖
a. δ%；b. θ%；c. α%；d. β%；e. γ%；f. EEGTP；g. EMGTP；h. β 率；i. WE
Figure5. The boxplots of nine characteristic parameters under four levels of anesthesiaa. δ%; b. θ%; c.α%; d. β%; e. γ%; f. EEGTP; g. EMGTP; h. β ratio; i. WE
EEGTP 反映的是 EEG 信號總能量，從圖 5f 可以看出，其值隨著麻醉深度的增加而增加，表明隨著麻醉深度的增加，EEG 的總能量也在逐漸增加。
EMGTP 反映的是肌電能量，有研究表明 EMG 對高 BIS 的診斷有相當大的幫助，尤其是對于確定 BIS > 80 是至關重要的[18]。如圖 5g 所示，隨著麻醉藥物與肌松劑的使用，麻醉深度逐漸增加，肌電能量則逐漸減少，表現為其值與麻醉深度成反比。
β 率反映的是高低頻能量比，從圖 5h 中可以看出，隨著麻醉深度的增加，EEG 中的頻率能量分布從高頻逐漸轉向低頻，所以 β 率的值也與麻醉深度成反比。
小波熵算法是從信號所包含的能量信息角度對信號復雜性的一種度量，對時間序列的非線性動力學特征可以有很好的表達[9]。從圖 5i 可以看出，WE 在清醒和輕度麻醉狀態有較高的值，在中度麻醉與深度麻醉中則有較低的值。
2.3 分類結果
在本研究中，從每 60 s 腦電信號中提取的九種特征被設置為 ANN 的輸入，以區分清醒、輕度、中度和深度麻醉狀態。本研究采用的 ANN 結構分為四層：輸入層有 9 個節點，第一個隱含層 9 個節點，第二個隱含層 16 個節點，輸出層 1 個結點。網絡中所有神經元單元的激活函數均為“tansig”函數。ANN 結構如圖 6 所示。
圖6
人工神經網絡結構圖
Figure6.
The structure diagram of artificial neural network
從 8 例病例中共提取出 21 821 個樣本，每個樣本包含 9 個特征參數（δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE）作為 ANN 的輸入，采用 7∶3 留出法測試分類方法，將整個數據集按照 7∶3 的比例分為訓練集與測試集，得到訓練集共 15 274 個樣本，測試集共 6 546 個樣本，以 BIS 指數為參考進行回歸訓練。我們還比較了 ANN 和 SVM 的性能，結果如表 4 所示。ANN 模型的分類準確度較高，達到 85.98%（配對樣本t 檢驗，t = 1.526，P = 0.127 > 0.05，表明 ANN 輸出的指數與 BIS 指數差異無統計學意義，可以認為來自于分布相同的總體），SVM 的分類準確度為 79.13%。同時，ANN 對四種麻醉狀態的靈敏度均高于 SVM。
表4
ANN 與 SVM 四種麻醉水平分類靈敏度與準確度比較
Table4.
Comparison of the sensitivity and accuracy calculated by the two methods of ANN and SVM
圖 7 為 7∶3 留出法的測試集通過訓練得到的 ANN 模型輸出的結果，紅色線段為 ANN 輸出指數，藍色線段為 BIS 指數，兩者之間的皮爾遜相關系數為 0.977 0，對四種麻醉水平的分類準確率為 85.98%。Bland-Altman 分析得到的結果如圖 8 所示，Bland-Altman 圖的橫坐標為 ANN 輸出指數與 BIS 指數的均值，縱坐標為 ANN 輸出指數與 BIS 指數的差值，最終計算得到的偏差為 ? 0.082 2，一致性的區間為 8.639 0 和 ? 8.803 5，表明偏差很小，一致性非常好。
圖7
測試集的 ANN 輸出指數與對應的 BIS 指數
Figure7.
The ANN output index and corresponding BIS index of the test set
圖8
測試集的 ANN 輸出指數與 BIS 指數的 Bland–Altman 圖
Figure8.
Bland–Altman diagram of ANN output index and BIS index of the test set
此外我們還采用 10 折交叉驗證法來驗證該模型的泛化性能，將整個數據集使用分層抽樣劃分成 10 個子集，每個子集包含 2 182 個樣本，使用 9 個子集的并集作為訓練集，余下的一個子集作為測試集，以 BIS 指數為參考進行回歸訓練，最后得到了與 7∶3 留出法相似的結果，表明 ANN 模型泛化性能較好。結果如表 4 所示。
3 討論
麻醉深度是外科手術中十分重要的指標。Benzy 等[9]使用了 WE 作為麻醉深度的度量，得到了較好的結果；Gu 等[1]使用排列熵與 95% 功率譜邊緣頻率、BetaRatio、SynchFastSlow 三個功率譜參數，結合 ANN 得到的模型對麻醉深度的監測也有較好的效果；劉軍等[19]引入樣本熵結合決策樹，嘗試了麻醉意識深度監測中腦電參數的選取，得到可完備地描述麻醉狀態變化過程的一套方案；王峰等[20]研究了小波變換和獨立分量分析相結合的去噪算法，并評價了 WE、樣本熵和復雜度這三種特征參數對 EEG 信號的狀態區分效果，得到結果為 WE 的狀態變化率最高。
在過去的這些研究中，有些研究人員嘗試使用基于腦電圖的特征結合機器學習的方法來評估麻醉深度，然而，他們只使用了一種基于 EEG 的特征，或者只是基于傅里葉變換的方法提取特征。因為腦電信號是非平穩非線性的，我們認為，利用 DWT 對腦電信號進行分解提取特征，引入 WE，與 ANN 相結合，能更全面地提取出 EEG 的有效信息，對麻醉深度的度量結果會更好。
在本研究中，我們提出了一種通過離散小波變換提取出 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，使用克魯斯卡爾-沃利斯統計檢驗對 9 種參數進行統計檢驗，并結合 ANN 來評估麻醉深度的方法，8 例全身麻醉患者的數據集被用來驗證本研究所提出的方法。最后的分類結果表明，使用 ANN 模型輸出的指標與 BIS 指數有較高的相關性和較好的一致性，并且分類準確度比 SVM 要更好。
本研究的局限性是，雖然我們計算得到的麻醉深度指數可以很好地區分清醒狀態和其他麻醉狀態，但由于樣本量小，沒有考慮到人體 EEG 的高變異性，而且，本研究也沒有測試藥物的變異性。這些局限性將在我們今后的工作中予以克服。
4 結論
本文使用小波變換提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，并結合 ANN 來評估麻醉深度。結果表明，該方法計算出的麻醉深度指數與 BIS 指數有較高的相關性和較好的一致性，并且能較好地區分清醒、輕度麻醉、中度麻醉、深度麻醉四種麻醉水平，表明本文所采用的方法對麻醉深度監測是可行的。在以后的工作中，我們將增加患者的數量和開展藥物的變異性影響分析，并增加多通道腦電信號對比采集，更全面地測試和優化上述麻醉深度計算方法，形成具有應用價值的獨立算法引擎。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
在手術過程中，全身麻醉是保證患者安全的必要和重要的環節，麻醉劑量過大可能會延長患者的恢復時間，而麻醉不足可能會導致患者產生術中記憶甚至在術中蘇醒，給患者帶來生理或者心理上的創傷[1]，所以在手術過程中進行準確的無創麻醉深度監測是非常必要的。由于麻醉藥物主要作用于中樞神經系統，在麻醉過程中，大腦活動的變化表現為腦電圖（electroencephalogram，EEG）波形或信號的振幅和頻率變化[2]，處理麻醉過程中的 EEG 信號有助于監測患者的麻醉深度。
美國 Aspect 公司于 1992 年開發的雙譜指數（bispectral index，BIS）監護儀由于具備實用性和有效性，獲得了市場的廣泛認可，是麻醉深度監測的可靠比較標準[3]，目前，臨床上主要采用基于 BIS 監護儀的麻醉深度評價方法。雖然 BIS 是一種已經應用于臨床數十年的指標，在降低術中意識發生率和指導麻醉藥物使用方面都取得了一定的成功，但是其算法并未公開。國內麻醉深度監護儀的需求較大，但是目前市場上可供選擇的產品卻非常有限[3]。由于從國外引進的儀器價格昂貴，國內制造的麻醉深度監護儀沒有技術上的改進，而麻醉深度監護儀最關鍵的技術就是其算法，所以，研究麻醉深度的計算方法是具有重要意義的。
早期基于傅立葉理論的頻譜分析方法假設數據在短時間內是平穩的，由于 EEG 信號在人體生物系統中呈現出非平穩、非線性的特征，這些基于頻率的方法可能會忽略有價值的信息[4]。近年來，許多研究者也采用非線性分析方法處理 EEG 信號，如希爾伯特-黃變換[5]、復雜性分析[6]、熵[4-5]和去趨勢波動分析[7]等。這些分析雖然提供了傳統 EEG 分析無法獲得的信息，但是也存在不足之處：希爾伯特-黃變換由于其計算量大不能應用于計算機上的實時分析[3]；復雜性分析、熵等非線性動力學方法對于初始條件敏感[8]；去趨勢波動分析方法在深度麻醉狀態下的分析效果不太好[9]。
小波變換作為時頻域分析工具，能夠同時反映信號的時域和頻域特征，很好地適用于非平穩的信號的分析，因而被廣泛用于生物醫學信號的分析[10-13]。EEG 信號通常被認為是非平穩的，而小波變換適合分析非平穩非線性的信號，基于 Mallat 快速算法的離散小波變換（discrete wavelet transform，DWT）是 EEG 信號分析的常用算法。模式識別與機器學習在 EEG 信號分析中也得到了廣泛的應用，一些研究人員使用支持向量機（support vector machine，SVM）和人工神經網絡（artificial neural network，ANN）做麻醉深度的回歸與分類[1,4,13-14]，而且 ANN 在模型開發中具有強大的自適應能力、學習能力、非線性映射能力、魯棒性和容錯性的優勢[13]。
本文使用 DWT 分析了麻醉中的腦電信號，并通過分解得到的小波系數提取了 δ 波、θ 波、α 波、β 波、γ 波的相對能量比，腦電總能量，肌電能量，β 率以及小波熵（wavelet entropy，WE）9 種特征參數，使用克魯斯卡爾-沃利斯統計檢驗來確定所使用的參數對四種麻醉水平是否具有顯著差異，并結合 ANN，將 BIS 作為參考輸出來評估麻醉深度，從靈敏度和分類準確度兩個方面對該方法的性能進行了評價。此外，我們還將 ANN 與 SVM 兩種模型進行了性能比較。
1 材料和方法
1.1 數據獲取
本研究的實驗數據均來自高州人民醫院麻醉科，使用美國 Aspect medical 公司生產的 BISX 外置模塊及數據記錄系統記錄了 8 名成人手術全過程全身麻醉狀態下的 EEG 信號。麻醉藥物基本選用芬太尼、異苯酚、卡肌寧、丙泊酚等靜脈麻醉藥物，麻醉方式均為靜吸麻醉方式。采集過程中 BISX 模塊濾波器模式設置為高通 0.5 Hz、低通 70 Hz、帶阻 50&60 Hz，以 256 Hz 的采樣率記錄，總共記錄約 30.5 h。數據已獲醫院授權。
1.2 EEG 信號處理流程
考慮到采集過程中存在電刀等干擾，所以使用自適應方差閾值法去除 EEG 信號中的尖銳噪聲，然后再將數據集分割成樣本；將分割好的樣本通過 DWT 提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，得到特征集；最后將特征集輸入模型進行訓練，得到可使用的能夠較為準確地計算麻醉深度指數的模型。
1.3 腦電數據預處理
1.3.1 自適應方差閾值法去噪[15 ]
在采集手術中的 EEG 信號時，由于外界環境電刀等帶來的干擾，會出現異常腦電數據段，表現為尖銳的噪聲。為了較好地去除大的波動干擾，又不對正常 EEG 信號造成影響，本文采用自適應調整的方法確定閾值，具體操作如下：選擇當前 60 s 作為一個時間段，將當前 60 s 的腦電數據分割成 60 個 1 s 的時間段，求得 60 個方差值，再利用這 60 個方差值的均值與中值計算得到自適應閾值，閾值設計公式如下：
|
式（1）中，NewT 是計算得到的方差閾值，Mean 與 Med 分別是最近 60 s 數據段 60 個方差的均值與中值，OldT 是上一個 60 s 數據段的閾值。每次檢測平移 1 s 時間段，并且更新閾值，實現方差閾值的自動調整。如果計算出該段數據的方差大于計算的閾值，則認定該段數據是較大的波動干擾，將其去除。
1.3.2 數據集分割
使用步長 5 s、大小為 60 s 的滑動窗，將 8 例腦電數據分割成每段 60 s 的 EEG 數據長度的樣本，每個樣本包含 15 360 個采樣點，共獲得 21 821 個樣本。整個數據集根據 BIS 分為四個級別。BIS 是根據麻醉深度顯示 0～100 范圍內的無量綱值，0 表示大腦基本沒有活動，100 表示完全清醒，范圍 0～100 可以細分為 4 個狀態，表 1 對不同麻醉狀態的 BIS 值范圍作了說明，中度麻醉狀態被認為是在全身麻醉下進行手術的合適水平[10]。
表1
不同麻醉狀態的 BIS 指數范圍與對應級別
Table1.
BIS index range and levels for different anesthesia states
圖 1 為不同麻醉階段對應的 60 s 長的 EEG 信號。由圖 1 可以看出，隨著麻醉深度的增加，EEG 信號幅度逐漸增大，頻率逐漸減小。在深度麻醉的狀態下，容易出現爆發抑制信號。
圖1
四種不同麻醉水平的 EEG 信號
Figure1.
EEG of four levels of anesthesia
1.4 使用 DWT 進行分解
給定基本小波函數：
|
其中，變量a、b分別是尺度因子和平移因子，且。通過小波基把任意空間中的函數展開，則信號的連續小波變換為：
|
取，，j, ，Z 是自然數，DWT 定義為：
|
通過對 a、b 進行采樣，小波變換可基于 Mallat 算法進行離散化，如圖 2 所示。
圖2
Mallat 算法分解架構
Figure2.
The Mallat algorithm decomposition architecture
其中， 為高通濾波器，將濾掉信號的低頻分量而輸出高頻分量， 為低通濾波器，濾掉信號的高頻分量而輸出低頻分量，根據下式將腦電信號分解成細節系數和近似系數：
|
|
式中 與 分別表示 EEG 信號分解得到的近似系數與細節系數。
DWT 分析的兩個重要考慮因素是選擇合適的母小波和分解層數。由于 Daubechies 小波被廣泛用于與 EEG 信號相關的應用[15]，本研究采用 db4 小波為母小波，使用 DWT 將 EEG 信號分解為不同的頻帶。小波分解的水平取決于五個腦電波（δ 波、θ 波、α 波、β 波、γ 波）以及肌電圖（electromyogram，EMG）信號的頻帶范圍。根據 EEG 信號頻帶分類，選取 5 階分解階數。因此，60 s 一段的 EEG 信號被分解成五個細節系數 D1-D5 和最后一個近似系數的 A5，其分別對應的頻帶如表 2 所示。
表2
EEG 信號頻帶分解
Table2.
EEG frequency band decomposition
1.5 特征提取
在不同的麻醉深度下，EEG 信號具有不同的波形幅度與頻率分布，可以通過 DWT 分解得到的小波系數，從中提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數。其中，δ%、θ%、α%、β%、γ% 用于量化 EEG 中的頻率成分的變化，EEGTP、EMGTP 代表 EEG 與 EMG 功率大小的變化，β 率有利于放大 EEG 信號中高頻成分與低頻成分的轉移程度，WE 則從熵的角度反映信號譜能量在各個子空間分布的有序或無序程度，因此，使用這些特征參數，可以更全面地量化 EEG 信號中的變化。
小波系數的能量給出了信號強度的信息，細節系數與近似系數的能量由以下方程給出：
|
|
這里 是分解級別，是對應的小波系數下標。
（1）由于 EEG 的能量基本集中在 0～64 Hz 頻帶范圍內，所以定義 EEG 總能量為：
|
（2）為了避免計算出的數值過大，對計算出的 EEG 總能量作對數處理，定義參數對數 EEG 能量 EEGTP 為：
|
（3）定義參數 δ%、θ%、α%、β%、γ% 分別為：
|
|
|
|
|
（4）由于 EMG 信號主要分布在 64～128 Hz 頻帶內，定義參數肌電能量 EMGTP 為：
|
（5）β 率定義為高頻的 β 波能量與低頻的 α 波能量之比，β 率計算公式為：
|
（6）WE 的計算公式如下：
① 首先計算分解得到的小波總能量：
|
② 計算小波子帶的相對能量
|
|
③ 計算 WE
|
1.6 克魯斯卡爾-沃利斯統計檢驗
由于利用 DWT 提取的腦電特征數據樣本量較少，數據不服從正態分布，所以本文采用克魯斯卡爾-沃利斯統計檢驗方法進行統計分析，以檢驗所提取的特征在清醒、輕度麻醉、中度麻醉和深度麻醉不同麻醉水平下的差異是否具有統計學意義。克魯斯卡爾-沃利斯統計檢驗方法是檢驗樣本是否來自相同分布的一種非參數方法[16]，它被用來確定兩組或多組自變量和連續或有序因變量之間的差異是否具有統計學意義，測定這些組在統計上是否不同。本文還使用了箱線圖描述每個特征在四種不同麻醉水平下的數據分布情況。
1.7 回歸與分類
ANN 是一種靈活的非參數的并行計算模型，基于模擬人腦的神經結構而發展起來[14]。ANN 通常由多個相互連接的多層節點組成，分為輸入層、隱含層和輸出層，輸入層中的每個節點接收外部信息，輸出層中的每個節點生成模型解決方案并輸出最終結果。本文使用的是誤差反向傳播神經網絡即 BP 神經網絡來輸出最終的麻醉深度指數，BP 神經網絡是一種多層前饋 ANN，也是目前應用最廣泛的 ANN 之一[1]。
SVM 是基于統計學習理論的一種機器學習方法[17]，它通過適當的非線性映射將輸入向量映射到一個高維的特征空間，使得數據總能被一個超平面分割成兩類，并尋找一個最優超平面。SVM 不僅能將兩類數據正確分開，而且使分類間隔最大，使其在有限樣本下具有了較強的泛化能力。在本研究中，我們使用 LIBSVM 工具箱完成 SVM 的相關分析。
1.8 分類性能分析
由于腦電特征數據樣本有限，我們使用固定的抽樣間隔抽取數據集的 70% 為訓練集，剩下的 30% 作為測試集來評估分類方法的性能，并且使用了 10 折交叉驗證評估了分類方法的泛化性能。根據測試集輸出的結果，利用靈敏度和分類精度可以對分類器的性能進行定量分析。靈敏度是一種麻醉狀態被正確地識別的數目與相應麻醉狀態的總數的比率。分類精度是指所有被正確識別的麻醉狀態的數目與實際麻醉狀態的總數的比率。
此外，本文還計算了 BIS 與所提出方法的指標之間的皮爾遜相關系數來評價所提出的方法，采用 Bland-Altman 分析來直觀反映兩種測量結果的一致性界限。最后，使用配對樣本t檢驗來檢驗最后 ANN 輸出的指數與 BIS 的差異有無統計學意義，判斷它們是否來自于分布相同的總體。
2 結果
2.1 自適應方差閾值法去除噪聲
方差閾值的計算方法如式（1）所述，如果發現該段數據的方差大于計算的閾值，則認定該段數據是較大的波動干擾，應將其去除。圖 3a 是一段 60 s 帶噪聲的原始 EEG 信號。圖 3b 中，藍線是計算得到的 60 s 的 EEG 數據的方差，紅線是計算得到的自適應閾值，從圖 3b 中可以看出，EEG 信號波動較大的地方計算得到的方差也較大，而通過自適應方差閾值法計算出的閾值則較為平穩，所以，只要判定其方差是否大于閾值，就可知道此數據段是否為噪聲。圖 3c 是使用自適應方差閾值法去除噪聲后的 EEG 信號，從圖 3c 中可以看出，波動較大的尖銳噪聲可以被很好地去除。
圖3
自適應方差閾值法去噪
a.帶噪聲的 EEG 信號；b.方差與閾值；c.去噪后的 EEG 信號
Figure3. Adaptive variance threshold denoising in EEGa. EEG with noise; b. variance and threshold; c. EEG after denoising
2.2 DWT 與特征提取
圖 4 顯示了 60 s 清醒狀態下的 EEG 信號 DWT 五級分解得到小波細節系數 D1～D5 與近似系數 A5 重構后的信號，從圖 4 中可以看出，DWT 能很好地將 EEG 各個頻段的信號區分開來，并且分解后的信號包含了原始 EEG 信號的全部信息。
圖4
清醒狀態下 60 s 的 EEG 信號 DWT 分解重構后的信號
Figure4.
The reconstructed signal of DWT decomposition of the 60 s EEG signal in the awake state
利用前面所述方法，對獲得的 21 821 個樣本提取特征，得到 21 821 個特征樣本，每個特征樣本包含 9 個特征。本文采用克魯斯卡爾-沃利斯統計檢驗方法檢驗提取的特征對四種麻醉狀態的鑒別能力，結果如表 3 所示，克魯斯卡爾-沃利斯統計檢驗得到的P值均小于 0.001，表明九個特征在四種不同麻醉水平下的差異具有統計學意義。
表3
克魯斯卡爾-沃利斯統計檢驗結果
Table3.
Kruskal-Wallis statistical test results
每個特征在不同麻醉深度下兩兩比較的結果中，只有兩組特征在Ⅰ vs. Ⅳ的差異沒有統計學意義：特征 θ% 的P值為 0.126 > 0.05，特征 α% 的P值為 0.170 > 0.05，剩余組的兩兩比較結果P值均小于 0.001，表明剩余組的差異均具有統計學意義。
不同麻醉深度下九種特征的箱線圖如圖 5 所示。從圖 5a～5e 中可以看出，δ% 反映的是 EEG 信號中低頻分量所占 EEG 總能量的比例，其值隨著麻醉深度的增加而增加。β% 與 γ% 在清醒階段的值較高，反映的是 EEG 信號中的高頻分量所占 EEG 總能量的比例，其值均與麻醉深度成反比。當受試者處于清醒狀態時，由于焦慮和壓力，代表高頻成分的 β 波與 α 波會更加活躍。然而，當麻醉時，代表低頻成分的 δ 波與 θ 波則變得活躍，因此，在中度麻醉和深度麻醉狀態時，δ% 的值較高。
圖5
九種特征參數在四種麻醉水平下的箱線圖
a. δ%；b. θ%；c. α%；d. β%；e. γ%；f. EEGTP；g. EMGTP；h. β 率；i. WE
Figure5. The boxplots of nine characteristic parameters under four levels of anesthesiaa. δ%; b. θ%; c.α%; d. β%; e. γ%; f. EEGTP; g. EMGTP; h. β ratio; i. WE
EEGTP 反映的是 EEG 信號總能量，從圖 5f 可以看出，其值隨著麻醉深度的增加而增加，表明隨著麻醉深度的增加，EEG 的總能量也在逐漸增加。
EMGTP 反映的是肌電能量，有研究表明 EMG 對高 BIS 的診斷有相當大的幫助，尤其是對于確定 BIS > 80 是至關重要的[18]。如圖 5g 所示，隨著麻醉藥物與肌松劑的使用，麻醉深度逐漸增加，肌電能量則逐漸減少，表現為其值與麻醉深度成反比。
β 率反映的是高低頻能量比，從圖 5h 中可以看出，隨著麻醉深度的增加，EEG 中的頻率能量分布從高頻逐漸轉向低頻，所以 β 率的值也與麻醉深度成反比。
小波熵算法是從信號所包含的能量信息角度對信號復雜性的一種度量，對時間序列的非線性動力學特征可以有很好的表達[9]。從圖 5i 可以看出，WE 在清醒和輕度麻醉狀態有較高的值，在中度麻醉與深度麻醉中則有較低的值。
2.3 分類結果
在本研究中，從每 60 s 腦電信號中提取的九種特征被設置為 ANN 的輸入，以區分清醒、輕度、中度和深度麻醉狀態。本研究采用的 ANN 結構分為四層：輸入層有 9 個節點，第一個隱含層 9 個節點，第二個隱含層 16 個節點，輸出層 1 個結點。網絡中所有神經元單元的激活函數均為“tansig”函數。ANN 結構如圖 6 所示。
圖6
人工神經網絡結構圖
Figure6.
The structure diagram of artificial neural network
從 8 例病例中共提取出 21 821 個樣本，每個樣本包含 9 個特征參數（δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE）作為 ANN 的輸入，采用 7∶3 留出法測試分類方法，將整個數據集按照 7∶3 的比例分為訓練集與測試集，得到訓練集共 15 274 個樣本，測試集共 6 546 個樣本，以 BIS 指數為參考進行回歸訓練。我們還比較了 ANN 和 SVM 的性能，結果如表 4 所示。ANN 模型的分類準確度較高，達到 85.98%（配對樣本t 檢驗，t = 1.526，P = 0.127 > 0.05，表明 ANN 輸出的指數與 BIS 指數差異無統計學意義，可以認為來自于分布相同的總體），SVM 的分類準確度為 79.13%。同時，ANN 對四種麻醉狀態的靈敏度均高于 SVM。
表4
ANN 與 SVM 四種麻醉水平分類靈敏度與準確度比較
Table4.
Comparison of the sensitivity and accuracy calculated by the two methods of ANN and SVM
圖 7 為 7∶3 留出法的測試集通過訓練得到的 ANN 模型輸出的結果，紅色線段為 ANN 輸出指數，藍色線段為 BIS 指數，兩者之間的皮爾遜相關系數為 0.977 0，對四種麻醉水平的分類準確率為 85.98%。Bland-Altman 分析得到的結果如圖 8 所示，Bland-Altman 圖的橫坐標為 ANN 輸出指數與 BIS 指數的均值，縱坐標為 ANN 輸出指數與 BIS 指數的差值，最終計算得到的偏差為 ? 0.082 2，一致性的區間為 8.639 0 和 ? 8.803 5，表明偏差很小，一致性非常好。
圖7
測試集的 ANN 輸出指數與對應的 BIS 指數
Figure7.
The ANN output index and corresponding BIS index of the test set
圖8
測試集的 ANN 輸出指數與 BIS 指數的 Bland–Altman 圖
Figure8.
Bland–Altman diagram of ANN output index and BIS index of the test set
此外我們還采用 10 折交叉驗證法來驗證該模型的泛化性能，將整個數據集使用分層抽樣劃分成 10 個子集，每個子集包含 2 182 個樣本，使用 9 個子集的并集作為訓練集，余下的一個子集作為測試集，以 BIS 指數為參考進行回歸訓練，最后得到了與 7∶3 留出法相似的結果，表明 ANN 模型泛化性能較好。結果如表 4 所示。
3 討論
麻醉深度是外科手術中十分重要的指標。Benzy 等[9]使用了 WE 作為麻醉深度的度量，得到了較好的結果；Gu 等[1]使用排列熵與 95% 功率譜邊緣頻率、BetaRatio、SynchFastSlow 三個功率譜參數，結合 ANN 得到的模型對麻醉深度的監測也有較好的效果；劉軍等[19]引入樣本熵結合決策樹，嘗試了麻醉意識深度監測中腦電參數的選取，得到可完備地描述麻醉狀態變化過程的一套方案；王峰等[20]研究了小波變換和獨立分量分析相結合的去噪算法，并評價了 WE、樣本熵和復雜度這三種特征參數對 EEG 信號的狀態區分效果，得到結果為 WE 的狀態變化率最高。
在過去的這些研究中，有些研究人員嘗試使用基于腦電圖的特征結合機器學習的方法來評估麻醉深度，然而，他們只使用了一種基于 EEG 的特征，或者只是基于傅里葉變換的方法提取特征。因為腦電信號是非平穩非線性的，我們認為，利用 DWT 對腦電信號進行分解提取特征，引入 WE，與 ANN 相結合，能更全面地提取出 EEG 的有效信息，對麻醉深度的度量結果會更好。
在本研究中，我們提出了一種通過離散小波變換提取出 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，使用克魯斯卡爾-沃利斯統計檢驗對 9 種參數進行統計檢驗，并結合 ANN 來評估麻醉深度的方法，8 例全身麻醉患者的數據集被用來驗證本研究所提出的方法。最后的分類結果表明，使用 ANN 模型輸出的指標與 BIS 指數有較高的相關性和較好的一致性，并且分類準確度比 SVM 要更好。
本研究的局限性是，雖然我們計算得到的麻醉深度指數可以很好地區分清醒狀態和其他麻醉狀態，但由于樣本量小，沒有考慮到人體 EEG 的高變異性，而且，本研究也沒有測試藥物的變異性。這些局限性將在我們今后的工作中予以克服。
4 結論
本文使用小波變換提取 δ%、θ%、α%、β%、γ%、EEGTP、EMGTP、β 率、WE 這 9 種特征參數，并結合 ANN 來評估麻醉深度。結果表明，該方法計算出的麻醉深度指數與 BIS 指數有較高的相關性和較好的一致性，并且能較好地區分清醒、輕度麻醉、中度麻醉、深度麻醉四種麻醉水平，表明本文所采用的方法對麻醉深度監測是可行的。在以后的工作中，我們將增加患者的數量和開展藥物的變異性影響分析，并增加多通道腦電信號對比采集，更全面地測試和優化上述麻醉深度計算方法，形成具有應用價值的獨立算法引擎。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
