CN101259015A - 一种脑电信号分析监测方法及其装置 - Google Patents

Abstract

本发明公开了一种脑电信号分析方法和装置。步骤如下：(1)采集原始的脑电信号；(2)去除原始脑电信号中的伪迹或噪声；(3)对预处理后的脑电信号进行相空间重构，得到一个递归图；(4)对递归图进行定量分析，计算其确定性数值DET。其中去除原始脑电信号中的伪迹，是采用小波变换和经验模式分解的混合方法。相空间重构是通过延迟时间的方法重建相空间向量，而得到递归图，再计算递归图对角结构上的递归点对所有递归点的比值，得到确定性数值DET。为医生提供良好的量化指标，在麻醉过程中提供准确的实时的麻醉深度的变化情况，或者用于脑功能变化、镇静与镇痛状况、认知功能的变化等研究分析。

Description

一种脑电信号分析监测方法及其装置

技术领域

本发明涉及脑电信号分析方法和装置，特别是涉及一种通过实时分析和监测脑电信号来描述大脑状态的方法和装置。

背景技术

脑电信号监测分析主要应用于脑功能变化、镇静与镇痛状况、认知功能的变化、麻醉深度监测和麻醉复苏变化的研究，它通常是将复杂的脑电信号转换成一个单一的指数，目前临床常用的脑电监测系统有脑电双频谱指数(BIS)、听觉诱发电位指数(AEI)、熵指数等等。

双频谱指数(BIS)分析方法是目前脑电监测麻醉镇静程度的一种常用方法(请参考Ira J Rampil，A primer for EEG signal processing in anesthesia，Anesthesiology，1998，89：980-1002)，迄今为止已有大量有关应用BIS进行麻醉深度评估的报道。这一方法在功率谱分析的基础上复合脑电相关函数谱分析技术，既测定脑电图(EEG)的线性(包括频率和功率)部分，又分析其非线性部分(包括位相和谐波)，使得双频谱指数分析在脑电图信号分析中具有重要价值。但是，BIS作为麻醉深度监测技术有其明显的局限性，许多因素都会对BIS产生影响，比如，BIS的脑电监护效果明显依赖于麻醉药的使用，如BIS与七氟醚的吸入浓度相关性良好，而与异氟醚吸入浓度则无相关性；N2O和氙气的麻醉对BIS值无影响，请参考Detsch O，Schneider G，Kochs E，et al.，Increasing isoflurane concentration may causeparadoxical increases in the EEG bispectral index in surgical patients，Br JAneasth，2000，84：33-37。

听觉诱发电位指数(AEI即AEPindex)是Mantzaridis根据听觉诱发电位形态学上的变化产生的数字化指标。请参考Mantzaridis H，Kenny GN，Auditory evoked potential index：a quantitative measure of changes in auditoryevoked potential during general anaesthesia.Anaesthesia，1997；52：1030-36。AEI的值从100-0的变化，为患者从清醒到脑电出现等电位的变化提供了良好的量化指标，能可靠地预测病人的意识消失，即AEI＜40。采用AEI的变化来判断手术期患者的记忆状态要比BIS更加有说服力，同时通过AEI的变化来反映患者从无意识状态到有意识状态的转变要比BIS更加灵敏。但是，听觉诱发电位指数的监测要求患者听力正常，因而对于伴发神经性耳聋或传导性耳聋的手术患者无法进行手术期的监测，还有儿童患者，其听力不同于成人的听力，在手术期监测以上的患者仍有一定限制。

近年来，非线性动力学分析开始运用于脑电信号的分析。常用的非线性分析手段包括：相图-相轨迹图；关联维数(Correlation Dimension，D₂)；相互维数(Mutual Dimension，D_m)；李亚普诺夫指数(LyapunovExponent，L₁)；柯尔莫哥诺夫熵(Kolmogorov Entropy，K₂)；复杂度；近似熵等等。中国专利申请03137747.5就是使用非线性动力学分析方法来监测和分析脑电信号。这些非线性动力学分析方法的主要缺点是对数据长度的依赖型和对噪声的敏感性。太长的数据会降低实时性，太短的数据则使结果不稳定；脑电信号中含有不同程度的噪声，这些噪声将直接影响结果。

发明简述

本发明针对现有技术的不足，提供了一种递归量化分析脑电信号的方法及其系统。

为了达到上述目的，本发明采用以下步骤：(1)采集原始的脑电信号；(2)去除原始脑电信号中的伪迹或噪声；(3)对去除伪迹预处理后的脑电信号进行相空间重构，得到一个递归图；(4)对递归图进行定量分析，计算其确定性数值DET。

所述步骤(2)去除原始脑电信号中的伪迹，本发明采用基于小波变换(wavelet transform)和经验模式分解(empirical mode decomposition，EMD)的混合方法。

所述步骤(3)是通过延迟时间的方法重建相空间，相空间向量x_i是基于过滤后的脑电信号u_i通过延迟时间的方法重建，如x_i＝(u_i，u_i+τ，...，u_i+(m-1)τ)，然后计算递归图(RP)R_i，j，如公式所示：

R_i，j＝Θ(ε-‖x_i-x_j‖)，i，j＝1，…，N.

其中N是状态空间向量的个数，ε是一个预定的截断距离，‖·‖是模(如欧几里得模，Euclidean norm)，Θ(x)是Heaviside函数。

所述步骤(4)中的DET是对角结构上的递归点对所有递归点的比值，如公式所示：

$\mathrm{DET}=\frac{{\∑}_{l=l_{\min }}^N\mathrm{lP}\left(l\right)}{{\∑}_{i,j}^N{R}_{i,j}},$

其中P(l)是递归图对角结构长度的频率分布，l_min是一个阈值。

附图说明

下面结合附图和具体实施例对本发明做进一步的详细描述。

图1是本发明的工作流程示意图；

图2A显示清醒状态下的一个原始的脑电信号；

图2B是对图2A中的原始脑电信号去除伪迹后得到的脑电信号；

图3显示一个脑电信号对麻醉反应的四个阶段，其中(a)是原始脑电信号，(b)是相应的去除了伪迹的滤波后的脑电信号；

图4显示图3中四个阶段的脑电信号的频率分布；

图5A显示一个原始的脑电信号序列，随着麻醉药物七氟醚吸入浓度增加而逐渐变化；

图5B显示七氟醚的浓度逐渐增加；

图5C显示确定性数值DET；

图6是一个利用本专利分析脑电信号来监测麻醉深度的系统。

具体实施例

图1所示为本发明的工作流程图。首先是步骤101，采集脑电信号。本实施例对9个在外科手术中的病人进行了脑电信号采集，其中有妇科、普通外科、整形外科手术。在本实施例中，采样频率是100Hz。采用Entropy^TM公司的复合电极记录前额和太阳穴处的脑电信号，在测量前，皮肤被清洗和干燥。

去除原始脑电信号中的伪迹(artefact)是非常重要的一个步骤，如步骤102。原始脑电信号中有一些典型的伪迹，包括肌电信号、眼运动和白噪声等。脑电活动是一个带宽很宽的信号(0～1000HZ)，同时外界的干扰信号也可能有各种频率成分，所以脑电信号就相当于是在各频率成分上都叠加有噪声的信号。如图2A所示，显示清醒状态下的一个原始脑电信号，它包括各种伪迹，如眼运动、眨眼、肌电信号等。这些脑电信号中的伪迹如果不去除，会导致许多错误的分析结果。

本发明采用一种新的去除伪迹的滤波方法，是基于小波变换(wavelettransform)和经验模式分解(empirical mode decomposition，EMD)方法。小波变换可以在多尺度空间分解脑电信号，然后针对分解的脑电信号进行阈值去噪，重构出一个新的信号。这种方法对脑电信号去噪，阈值的选择是主要困难。融合经验模式分解方法(EMD)可以提高小波滤波器的功能。EMD方法是基于数据本身的特性，采用自适应的方法逐步地分解信号，此过程非常相似于脑电信号叠加的反过程，因此，结合小波变换和EMD方法很适合于脑电信号的分解，最终可以去掉噪声信号。

在一个4秒时间的脑电信号上进行6阶小波变换，使用Sym8小波核，将得到在频带0.5～1.5，1.5～3，3～6，6～12，12～25和25～50Hz上的脑电信号。去除低频带0～3Hz和高频带25～50Hz的信号，因为这些频带的信号含有太多伪迹，如在频带25～50Hz上的脑电信号通常含有大量的肌电信号。为了去除在各个频带成分上的白噪声信号，本实施例使用贝叶斯估计(Bayesianestimation)为小波分解后不同尺度下的信号设定不同的阈值，超过阈值的信息将保留下来。阈值的标准是基于贝叶斯形式和假设的先验分布。最后，使用离散小波反变换来重构脑电信号序列。为了进一步去除肌电信号伪迹和其它噪声，对保留的脑电信号采用经验模式分解(EMD)方法。对于一个给定的脑电信号s(t)，EMD方法可以将其分解成多个固有模式函数(intrinsic mode functions，IMFs)的线性组合C_n(n＝1，2，...，N)，如公式(1)所示：

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{C}_i\left(t\right)+r_{N+1}\left(t\right),---\left(1\right)$

其中N是固有模式函数(IMFs)的个数，r_N+1(t)是脑电信号的余量，表示一个趋势。分解过程依赖于信号本身的振动特性，也就是说，固有模式函数(IMFs)的个数是由振动特性来决定的。在EMD分解之后，会发现肌电信号伪迹主要集中在第一个固有模式函数(最高频成分)和r(t)趋势信息(最低频成分)上。于是在本实施例中，删除第一个固有模式函数和余量信息。然后将其它的固有模式函数迭加起来得到一个新的脑电信号序列u_i。

图2B所示为在图2A的原始脑电信号中去除伪迹后得到的新的脑电信号。再参见图3，图3显示脑电信号对麻醉反应的四个阶段。其中(a)是原始脑电信号，(b)是相应的去除了伪迹的滤波后的脑电信号。阶段I(a)和(b)显示在清醒状态下眼运动伪迹的去除，阶段II(a)和(b)显示去除原始脑电信号中的慢波，在阶段III和IV，脑电信号中主要包括麻醉过程中的低频波(5～12Hz)。图4显示图3中四个阶段过滤后的频率分布，在清醒阶段I，脑电信号主要包括两部分，一是5～12Hz的波，另一是20Hz以上的快波；在麻醉的初始阶段II，在β波段功率上有一个急升(12～20Hz)，而且在低频上也有一个幅度的升高；随着麻醉药物七氟醚浓度的增加，病人从轻度麻醉进入深度麻醉阶段III，12～25Hz的β波的幅度显著下降，进一步的麻醉深度使得低频波的幅度增加，如阶段IV。

如图1，在去除脑电信号中的伪迹(步骤102)之后，对新的脑电信号进行相空间重构，如步骤103。为了得到递归图，首先计算一个如公式(2)所示的N×N矩阵：

R_i，j＝Θ(ε-‖x_i-x_j‖)，i，j＝1，…，N.(2)

其中N是状态空间向量的个数，ε是一个预定的截断距离，‖·‖是模(如欧几里得模，Euclidean norm)，Θ(x)是Heaviside函数。相空间向量x_i可以基于过滤后的脑电信号u_i通过延迟时间的方法重建，如x_i＝(u_i，u_i+τ，...，u_i+(m-1)τ)。在本实施例中，4秒的脑电信号(400个数据电)被用来建立相空间，相空间中的维数m和延迟时间τ分别由最近邻域法和互相关函数来确定。在本实验中m＝12，τ＝5。截断距离ε定义一个以x_i为中心的球体，如果x_i落在球体内部的话，即状态靠近x_j的话，R_i，j＝1，否则R_i，j＝0。从而第i时刻在第j时刻的递归状态可以通过一个2-维方阵用黑点(R_i，j＝1)或者白点(R_i，j＝0)表示，R_i，j在此被称为递归图(RP)。

在重构相空间之后，对递归图上单个的点及其对角线结构分布进行统计分析，得到递归定量分析，如图1中的步骤104。递归图的量化分析就称为递归定量分析(RQA)。RQA可以量化递归图的结构，包括一个动态系统递归的数量和周期。在本实施例中，引入一个确定性(determinism，DET)概念，参数DET是一个动态系统的确定性量度，是递归图中对角线大于l_min的点对所有递归点的比值，如公式(3)所示：

$\mathrm{DET}=\frac{{\∑}_{l=l_{\min }}^N\mathrm{lP}\left(l\right)}{{\∑}_{i,j}^N{R}_{i,j}},---\left(3\right)$

其中P(l)是递归图对角结构长度的频率分布，则 ${\∑}_{l=l_{\min }}^N\mathrm{lP}\left(l\right)$ 是递归图中对角线大于l_min的点的个数，l_min是一个阈值，在此实施例中，l_min＝2。DET是基于递归图中的对角线来度量系统的可预测性，对于一随机或浑电系统，其递归图中的对角线很短或没有，而对于一确定性的系统其对角线显得很长，没有单个或孤立的对角线。

图5是一个麻醉深度监测的例子。图5A显示一个原始的脑电信号序列，随着麻醉药物七氟醚吸入浓度增加而逐渐变化；图5B显示七氟醚浓度的逐渐增加，将图5A的原始脑电信号输入滤波器，去除其中的噪声及伪迹后，可以计算出脑电信号的确定性数值DET，如图5C所示。随着七氟醚浓度的增加，过滤后的脑电信号的确定性数值DET显著增加，直到达到一个高度。

DET量化值是一个动态系统的预测性量度，它意味着神经活动随着麻醉药物浓度的增加变得更可预测(低复杂性)，这是DET量化值的物理含义。

尽管许多研究表明脑电信号多数呈现非线性特征，而且最近有很多方法也使用非线性动力学分析方法来实现麻醉监测，但是实际上脑电信号不仅仅是非线性的，有时也会呈现出线性的特征。本发明充分考虑这个可能性，对脑电信号的特征不做任何假设，既可以适用于含噪声的线性或非线性脑电信号，又可以适用于短的数据。这一特点使得本发明方法优于现有方法，比如专利申请03137747.5的方法，就没有分析数据长度的影响。实际上不同长度的数据会产生不同的结果，通常数据越长，结果越稳定，但是实时性下降。本发明对数据长度并不是很敏感，可以适用于短数据，这对临床使用很有帮助。

图6是一个利用本专利分析脑电信号来监测麻醉深度的系统。首先利用电极601采集脑电信号，然后通过模拟/数字转换器602转换成数字信号，以便进一步处理。数字滤波器603将采集到的脑电信号中的伪迹去除，这里数字滤波器可以单独是经验模式分解(EMD)模块605，也可以包括小波变换模块604和EMD模块605，小波变换模块604使用上述方法在多尺度空间分解脑电信号，然后针对分解的脑电信号进行阈值去噪，EMD模块605使用上述自适应的方法逐步地分解信号，最终可以去掉噪声信号。去除了伪迹的脑电信号再经过相空间重建模块606，建立相空间向量，递归定量分析RQA模块607对相空间向量矩阵进行计算，得到确定型数值DET，并在显示模块608上直观地显示出来，为医生提供良好的量化指标，在麻醉过程中提供准确的实时的麻醉深度的变化情况，或者用于脑功能变化、镇静与镇痛状况、认知功能的变化等研究分析。

Claims (10)

1. 一种脑电信号分析方法，包括以下步骤：

(1)采集原始的脑电信号；

(2)去除原始脑电信号中的伪迹或噪声；

(3)对去除了伪迹的脑电信号进行相空间重构，得到一个递归图；

(4)对递归图进行定量分析，计算其确定性数值DET。

2. 根据权利要求1所述的方法，其特征在于：其中所述步骤(2)是使用小波变换方法在多尺度空间分解脑电信号，得到多个频带上的信号，去掉3Hz以下的和25Hz以上的信号；再使用贝叶斯估计对小波变换后不同尺度下的信号设定不同的阈值，超过阈值的信息将保留下来，使用所述离散小波反变换来重构一个新的脑电信号序列。

3. 根据权利要求2所述的方法，其特征在于：所述阈值的标准是基于贝叶斯形式和假设的先验分布。

4. 根据权利要求3所述的方法，其特征在于还包括：对所述新的脑电信号s(t)，将其分解成多个固有模式函数的线性组合C_n(n＝1，2，...，N)，如以下公式所示：

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{C}_i\left(t\right)+r_{N+1}\left(t\right),$

其中N是固有模式函数的个数，r_N+1(t)是脑电信号的余量；删除第一个固有模式函数和余量信息；然后，将其他的固有模式函数迭加起来得到一个新的过滤后的脑电信号序列u_i。

5. 根据权利要求1所述的方法，其特征在于，在其中所述步骤(2)，对所述原始脑电信号s(t)，将其分解成多个固有模式函数的线性组合C_n(n＝1，2，...，N)，如以下公式所示：

$s\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{C}_i\left(t\right)+r_{N+1}\left(t\right),$

其中N是固有模式函数的个数，r_N+1(t)是脑电信号的余量；删除第一个固有模式函数和余量信息；然后，将其他的固有模式函数迭加起来得到一个新的过滤后的脑电信号序列u_i。

6. 根据权利要求4或5所述的方法，其特征在于，在所述步骤(3)，对所述过滤后的脑电信号u_i通过延迟时间的方法重建相空间，相空间向量x_i＝(u_i，u_i+τ，...，u_i+(m-1)τ)，然后如下公式计算N×N矩阵的递归图：

R_i，j＝Θ(ε-||x_i-x_j||)，i，j＝1，...，N.

其中N是状态空间向量的个数，ε是一个预定的截断距离，||·||是模，Θ(x)是Heaviside函数。

7. 根据权利要求6所述的方法，其特征在于，在所述步骤(4)中，计算所述递归图R_i，j对角结构上的递归点对所有递归点的比值，即确定性数值DET，如下公式：

$\mathrm{DET}=\frac{{\mathrm{\Σ}}_{l=l_{\min }}^N\mathrm{lP}\left(l\right)}{{\mathrm{\Σ}}_{i,j}^N{R}_{i,j}},$

其中P(l)是递归图对角结构长度的频率分布，l_min是一个阈值。

8. 一种脑电信号分析装置，包括：

脑电信号采集电极；

模拟/数字转换器，用于将模拟脑电信号转换成数字信号；

数字滤波器，用于去除脑电信号中的伪迹；

其特征在于，还包括：

相空间重建模块，用于建立脑电信号的相空间向量；和

递归定量分析模块，用于计算确定性数值DET。

9. 根据权利要求8所述的脑电信号分析装置，其特征在于，所述数字滤波器包括小波变换模块和经验模式分解模块。

10. 根据权利要求9所述的脑电信号分析装置，其特征在于，所述数字滤波器是经验模式分解。

Images