全身麻醉是外科手术中保证患者安全的必不可少的部分，脑电图（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 指数有较高的相关性和较好的一致性，并且能较好地区分清醒、轻度麻醉、中度麻醉、深度麻醉四种麻醉水平，表明本文所采用的方法对麻醉深度监测是可行的。在以后的工作中，我们将增加患者的数量和开展药物的变异性影响分析，并增加多通道脑电信号对比采集，更全面地测试和优化上述麻醉深度计算方法，形成具有应用价值的独立算法引擎。
利益冲突声明：本文全体作者均声明不存在利益冲突。
