CN103690163B - 基于ica和hht融合的自动眼电干扰去除方法 - Google Patents

Abstract

本发明公开了一种基于ICA和HHT融合的自动眼电干扰去除方法。首先，用独立成分分析将采集到的包含眼电信号的脑电信号分解为若干独立成分。希尔伯特-黄变换由两部分构成：经验模态分解和希尔伯特谱分析，然后，对每一个独立成分通过经验模态分解提取趋势项，计算其统计特征，从而确定包含眼电信号的独立成分，利用希尔伯特谱分析，保留这些独立成分中不属于眼电信号的高频成分，同时去除属于眼电信号的低频成分；本发明使脑电信号信号中不包含眼电信号的段不受影响，去除眼电信号后，在没有眼电信号发生的区域，更接近于原始的脑电信号。

Description

基于ICA和HHT融合的自动眼电干扰去除方法

技术领域

本发明具体涉及一种基于ICA和HHT融合的自动眼电干扰去除方法。

背景技术

在过去的几十年中，对脑电波(Electroencephalogram，EEG)信号的分析与研究在神经科学、心理学以及生物医学等领域起到了重要的作用。然而，由于EEG信号的电压非常小，通常只有50微伏左右，在记录过程中经常会受到干扰源的干扰，从而产生一些伪迹，这些伪迹经常妨碍脑电波的分析与研究。在伪迹中，眼电(Electrooculogram，EOG)信号是最主要的一种干扰，它会随机地出现在EEG信号中，并且幅度较大，EOG分为水平EOG和垂直EOG。水平EOG是由于眼部随着关注目标的移动而产生，可以通过被试在实验中减少眼部的水平运动而减少甚至消除，对于垂直EOG，是由于眨眼产生，如果要求被试在实验中，尽量减少眨眼的次数，被试会产生不适感，从而对实验造成影响，此外，还可以通过舍弃EEG信号中被EOG干扰的部分，但是这种方法可能会引起重要信息的丢失。因此，如何去除EOG的干扰，寻找行之有效的检测并提取出真实可靠的EEG信号的方法，具有重要意义。

目前检测并去除EOG的方法主要有如下三种：基于回归分析的方法，这种方法出现的最早，并被广泛应用于流行的EEG诊断软件中；基于主成分分析(Principal Component Analysis，PCA)的方法，由于该方法要求的条件比较严格，因此应用不多；基于独立成分分析(Independent Component Analysis，ICA)的方法，EEG信号是自发EEG信号与各种伪迹的线性混合，满足信号源独立的条件，因此EEG信号适合用该方法去除EOG干扰。

基于回归分析的方法会错误地将包含的EEG信号去除。基于PCA的方法去除EOG干扰后，会使没有EOG干扰的EEG信号段发生变形，从而影响对EEG信号的分析。基于ICA的方法无法实现自动去除含有EOG干扰的独立成分(Independent Component，IC)，并且由于其分解IC过程是一个估计的过程，因此分解得到的IC并不准确，会将部分高频的EEG成分去除，从而影响对EEG信号的分析。

发明内容

针对上述基于独立成分分析(ICA)去除眼电信号(EOG)方法的不足，本发明提供一种基于独立成分分析(ICA)和希尔伯特-黄变换(HHT)融合的自动眼电干扰去除方法，该方法是一种不需要人为干预的自动眼电信号(EOG)去除方法，同时使去除眼电信号(EOG)后的脑电信号(EEG)，在没有眼电信号(EOG)的部分更接近于原始的脑电信号(EEG)。

本发明的目的是通过以下技术方案实现的：一种基于独立成分分析(ICA)和希尔伯特-黄变换(HHT)融合的自动眼电干扰去除方法，具体步骤如下：

步骤一：通过独立成分分析(ICA)将原始脑电信号(EEG)分解为若干独立成分(IC)；

步骤二：对所有独立成分(IC)分段并加窗，段长为10秒，在每一段内的所有独立成分(IC)上加窗，窗长为一次眨眼的时长，为400毫秒，为了减小眼电发生在窗的边缘部分对处理过程产生的影响，采用窗间交叠的办法，交叠部分为半个窗长；

步骤三：通过经验模态分解(EMD)提取窗内各独立成分(IC)的趋势项，经验模态分解(EMD)是基于时域局部特征的，能够把复杂数据分解成少量的固有模态函数(IMF)，每个固有模态函数(IMF)要满足如下两个条件：(1)在整个数据序列中，极值点的数量与过零点的数量相等，或相差一个；(2)在任一时间点上，信号的局部最大值和局部最小值定义的包络均值为零；

步骤四：计算段内每个独立成分(IC)的统计特征，包括段内各窗趋势项最大值与最小值之差的和与方差，计算公式如下所示：

range(j)_i，w＝max(j)_i，w-min(j)_i，w

${\mathrm{diffsum}\left(j\right)}_i=\underset{i=1}{\overset{W}{\mathrm{\Σ}}}\mathrm{range}{\left(j\right)}_{i,w}$

$\mathrm{diffvar}{\left(j\right)}_i=\frac{1}{W}\underset{i=1}{\overset{W}{\mathrm{\Σ}}}{(\mathrm{range}{\left(j\right)}_{i,w}-\mathrm{meanrange}{\left(j\right)}_i)}^2$

其中，j∈1....N表示独立成分分析(ICA)分解得到的独立成分(IC)，N为独立成分分析(ICA)分解得到的独立成分(IC)数量，max(j)_i，w、max(j)_i，w和range(j)_i，w分别表示第j个独立成分(IC)的第i段的第w个窗的趋势项的最大值、最小值以及它们的差，meanrange(j)_i表示第j个独立成分(IC)的第i段中所有窗的趋势项的最大值与最小值之差的均值，diffsum(j)_i和diffvar(j)_i分别表示第j个独立成分(IC)的第i段中的所有窗的趋势项的最大值与最小值之差的和以及方差，由此形成一个特征向量(diffsum(j)_i，diffvar(j)_i)；

步骤五：对每一段中提取的统计特征通过K-均值聚类算法聚成两类，通过计算这两类类间可分性测度，判断其中一类的独立成分(IC)是否包含眼电信号；

步骤六：由于信号并不是由单一的频率成分构成的，而眼电信号(EOG)通常被认为是低频成分，因此，通过计算瞬时频率，保留高于眼电信号(EOG)频率的成分；通过对包含眼电信号(EOG)的独立成分(IC)在被判断为含有眼电信号(EOG)的段中，经过经验模态分解(EMD)分解得到的固有模态函数(IMF)做希尔伯特谱分析(HSA)，计算其瞬时频率，保留高于4Hz的频率成分，合成一个新的独立成分(IC)；

步骤七：将所有独立成分(IC)通过独立成分分析(ICA)得到的混合矩阵，重新合成一个新的脑电信号(EEG)，即为去除眼电信号(EOG)后的脑电信号(EEG)。

本发明还具有如下特征：

1、以上方法的步骤三所述的趋势项具体提取过程如下所示：

(1)根据信号x(t)的极大点和极小点，获取其上包络r_up(t)和下包络r_low(t)；

(2)考察x(t)与上下包络均值之差d(t)：

$d\left(t\right)=r\left(t\right)-\frac{(r_{\mathrm{up}}\left(t\right)+r_{\mathrm{low}}\left(t\right))}{2}$

(3)将d(t)视为新的x(t)重复(1)和(2)，直到d(t)满足固有模态函数(IMF)条件，记为c₁(t)＝d(t)，c₁(t)即为一个固有模态函数(IMF)，

r(t)＝x(t)-c₁(t)

r(t)即为新的x(t)

(4)重复(1)、(2)和(3)，直到残差为零或单调函数时停止，最后得到关于x(t)的表达式：

$x\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_j\left(t\right)+r_N\left(t\right)$

式中，N代表分解得到的固有模态函数(IMF)的个数，残差r_N(t)即为该段信号的趋势项。

2、以上方法的步骤五的判断其中一类的独立成分(IC)是否包含眼电信号的计算过程如下所示：

$J=\frac{J_b}{J_w}$

其中J_b和J_w分别表示类内和类间的均方距离，计算公式如下所示，

$J_w=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{P}_j{J}_j$

其中，T表示矩阵转置，P_j表示样本属于第j类的先验概率，通过第j类中样本数量占样本总数量的比值来估计，为第j类的样本均值，为所有样本的均值，是第j类样本的均方距离，计算公式如下所示，

其中，N_j是第j类中样本的数量，是第j类的第i个样本向量，设定一个阈值J_th，如果J＞J_th，则表示聚类均值大的一类中的独立成分(IC)包含眼电信号(EOG)信号，否则，该段不包含眼电信号(EOG)信号；对所有判断为包含眼电信号(EOG)的段以及所属独立成分(IC)进行统计，对J_th的设置要确保含有眼电的独立成分(EOG-IC)的检测有较高的准确率，这里设为6.0；本发明具有如下有益的效果：

1、将独立成分(IC)分段考虑，包含眼电信号(EOG)的段和不包含眼电信号(EOG)的段分别处理，使脑电信号(EEG)中不包含眼电信号(EOG)的段不受影响；

2、通过经验模态分解(EMD)获取各段独立成分(IC)的趋势项，然后获取这些趋势项的统计特征，整个过程不需要人为干预，从而实现了自动地判断哪个独立成分(IC)包含眼电信号(EOG)，以及出现在哪一段；

3、通过对包含眼电信号(EOG)的独立成分(IC)由经验模态分解(EMD得到的固有模态函数(IMF)，利用希尔伯特谱分析(HSA)，计算这些独立成分(IC)各频率成分的瞬时频率，保留高于眼电信号(EOG)频率的频率成分，相比其他方法，在脑电波信号中包含眼电信号(EOG)的段，去除眼电信号(EOG)后，在没有眼电信号(EOG)发生的区域，更接近于原始的脑电信号(EEG)。

附图说明

图1为本发明的执行过程示意图；

图2为本发明的包含EOG的窗和不包含EOG的窗的EEG信号，以及它们的趋势项示意图；

图3为本发明的对各IC的所有段的统计特征聚类的示意图；

图4为本发明的原始EEG信号和去除EOG后EEG信号的示意图；

图5为本发明的Fp1导联去除EOG前后的EEG信号示意图。

具体实施方式

下面结合附图对本发明作进一步的说明：

实施例1

按照本发明的技术方案，结合HHT中包含的两部分：EMD和HSA，分别在趋势项特征提取和瞬时频率计算的应用，实现对原始EEG信号经过ICA后，自动判断哪些分量包含EOG以及这些分量中的哪一段包含EOG，并且尽可能的减小去除EOG操作对不包含EOG部分的影响，为从EEG信号中去除EOG提供了一种新的方法。发明的具体实施方法如下：

步骤一：通过ICA将原始EEG信号分解为若干IC。一般认为，采集EEG数据时使用多少导联，ICA便将其分解为多少个IC。

步骤二：对所有IC分段并加窗。段长选择要适中，太短不利于后面统计特征的提取，太长必然会每一段都包含EOG信号，不利于去除EOG操作对减小不包含EOG部分的影响，同时还要考虑被试在该段时间内可以不眨眼并且没有不适感，通常认为这个时间段为10秒比较合适。在每一段内的所有IC上加窗，窗长为一次眨眼的时长，通常为400毫秒。为了减小眼电发生在窗的边缘部分对处理过程产生的影响，我们采用窗间交叠的办法，交叠部分为半个窗长。

步骤三：通过EMD提取窗内各IC的趋势项(如图2所示)。

EMD是基于时域局部特征的，可以把复杂数据分解成少量的IMF，每个IMF要满足如下两个条件：(1)在整个数据序列中，极值点的数量与过零点的数量相等，或最多相差不能多于一个；(2)在任一时间点上，信号的局部最大值和局部最小值定义的包络均值为零。分解是自适应别且高效的，特别适合用来分析非平稳非线性的时变过程，能清晰地分辨出交叠复杂数据的内蕴模式。

趋势项可以反映一定长度内信号的变化趋势，具体提取过程如下所示：

(1)根据信号x(t)的极大点和极小点，获取其上包络r_up(t)和下包络r_low(t)；

(2)考察x(t)与上下包络均值之差d(t)：

$d\left(t\right)=r\left(t\right)-\frac{(r_{\mathrm{up}}\left(t\right)+r_{\mathrm{low}}\left(t\right))}{2}---(1-1)$

(3)将d(t)视为新的x(t)重复(1)和(2)，直到d(t)满足IMF条件，记为c₁(t)＝d(t)，c₁(t)即为一个IMF，

r(t)＝x(t)-c₁(t)            (1-2)

r(t)即为新的x(t)

(4)重复(1)、(2)和(3)，直到残差为零或单调函数时停止。

最后得到关于x(t)的表达式：

$x\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_j\left(t\right)+r_N\left(t\right)---(1-3)$

式中，N代表分解得到的IMF的个数，残差r_N(t)即为该段信号的趋势项。

步骤四：计算段内每个IC的统计特征，包括段内各窗趋势项最大值与最小值之差的和与方差，计算公式如下所示：

range(j)_i，w＝max(j)_i，w-min(j)_i，w        (1-4)

${\mathrm{diffsum}\left(j\right)}_i=\underset{i=1}{\overset{W}{\mathrm{\Σ}}}\mathrm{range}{\left(j\right)}_{i,w}---(1-5)$

$\mathrm{diffvar}{\left(j\right)}_i=\frac{1}{W}\underset{i=1}{\overset{W}{\mathrm{\Σ}}}{(\mathrm{range}{\left(j\right)}_{i,w}-\mathrm{meanrange}{\left(j\right)}_i)}^2---(1-6)$

其中，j∈1....N表示ICA分解得到的IC，N为ICA分解得到的IC数量，max(j)_i，w、max(j)_i，w和range(j)_i，w分别表示第j个IC的第i段的第w个窗的趋势项的最大值、最小值以及它们的差，meanrange(j)_i表示第j个IC的第i段中所有窗的趋势项的最大值与最小值之差的均值，diffsum(j)_i和diffvar(j)_i分别表示第j个IC的第i段中的所有窗的趋势项的最大值与最小值之差的和以及方差，由此形成一个特征向量(diffsum(j)_i，diffvar(j)_i)。

步骤五：对每一段中提取的统计特征通过K-均值聚类算法聚成两类(如图3所示)，我们通过计算这两类类间可分性测度，判断其中一类的IC是否包含眼电信号，计算过程如下所示。

$J=\frac{J_b}{J_w}---(1-7)$

其中J_b和J_w分别表示类内和类间的均方距离，计算公式如下所示。

$J_w=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{P}_j{J}_j---(1-9)$

其中，T表示矩阵转置，P_j表示样本属于第j类的先验概率，可以通过第j类中样本数量占样本总数量的比值来估计，为第j类的样本均值，为所有样本的均值，是第j类样本的均方距离，计算公式如(1-10)所示。

其中，N_j是第j类中样本的数量，是第j类的第i个样本向量。

设定一个阈值J_th，如果J＞J_th，则表示聚类均值大的一类中的IC包含EOG信号，否则，该段不包含EOG信号。对所有判断为包含EOG的段以及所属IC进行统计。对J_th的设置要确保EOG-IC的检测有较高的准确率，这里设为6.0。

步骤六：由于信号并不是由单一的频率成分构成的，而EOG通常被认为是低频成分，因此，通过计算瞬时频率，保留高于EOG频率的成分。我们通过对包含EOG的IC在被判断为含有EOG的段中，经过EMD分解得到的IMF做HSA，计算其瞬时频率。一般认为，EOG频率低于4Hz，因此我们保留高于4Hz的频率成分，合成一个新的IC。

步骤七：将所有IC通过ICA得到的混合矩阵，重新合成一个新的EEG信号，即为去除EOG后的EEG信号(如图4和图5所示)。

Claims (3)

1.基于独立成分分析和希尔伯特-黄变换融合的自动眼电干扰去除方法，其特征在于：具体步骤如下：

步骤一：通过独立成分分析将原始脑电信号分解为若干独立成分；

步骤二：对所有独立成分分段并加窗，段长为10秒，在每一段内的所有独立成分上加窗，窗长为一次眨眼的时长，为400毫秒，为了减小眼电发生在窗的边缘部分对处理过程产生的影响，采用窗间交叠的办法，交叠部分为半个窗长；

步骤三：通过经验模态分解提取窗内各独立成分的趋势项，经验模态分解是基于时域局部特征的，能够把复杂数据分解成少量的固有模态函数，每个固有模态函数要满足如下两个条件：(1)在整个数据序列中，极值点的数量与过零点的数量相等，或相差一个；(2)在任一时间点上，信号的局部最大值和局部最小值定义的包络均值为零；

步骤四：计算段内每个独立成分的统计特征，包括段内各窗趋势项最大值与最小值之差的和与方差，计算公式如下所示：

range(j)_i，w＝max(j)_i，w-min(j)_i，w

$\mathrm{diffsum}{\left(j\right)}_i=\underset{i=1}{\overset{W}{\mathrm{\Σ}}}\mathrm{range}{\left(j\right)}_{i,w}$

$\mathrm{diffvar}{\left(j\right)}_i=\frac{1}{W}\underset{i=1}{\overset{W}{\mathrm{\Σ}}}{(\mathrm{range}{\left(j\right)}_{i,w}-\mathrm{meanrange}{\left(j\right)}_i)}^2$

其中，j∈1....N表示独立成分分析分解得到的独立成分，N为独立成分分析分解得到的独立成分数量，max(j)_i，w、min(j)_i，w和range(j)_i，w分别表示第j个独立成分的第i段的第w个窗的趋势项的最大值、最小值以及它们的差，meanrange(j)_i表示第j个独立成分的第i段中所有窗的趋势项的最大值与最小值之差的均值，diffsum(j)_i和diffvar(j)_i分别表示第j个独立成分的第i段中的所有窗的趋势项的最大值与最小值之差的和以及方差，由此形成一个特征向量(diffsum(j)_i，diffvar(j)_i)；

步骤五：对每一段中提取的统计特征通过K-均值聚类算法聚成两类，通过计算这两类类间可分性测度，判断其中一类的独立成分是否包含眼电信号；

步骤六：由于信号并不是由单一的频率成分构成的，而眼电信号通常被认为是低频成分，因此，通过计算瞬时频率，保留高于眼电信号频率的成分；通过对包含眼电信号的独立成分在被判断为含有眼电信号的段中，经过经验模态分解得到的固有模态函数做希尔伯特谱分析，计算其瞬时频率，保留高于4Hz的频率成分，合成一个新的独立成分；

步骤七：将所有独立成分通过独立成分分析得到的混合矩阵，重新合成一个新的脑电信号，即为去除眼电信号后的脑电信号。

2.根据权利要求1所述的基于独立成分分析和希尔伯特-黄变换融合的自动眼电干扰去除方法，其特征在于，步骤三所述的趋势项具体提取过程如下所示：

(1)根据信号x(t)的极大点和极小点，获取其上包络r_up(t)和下包络r_low(t)；

(2)考察x(t)与上下包络均值之差d(t)：

$d\left(t\right)=r\left(t\right)-\frac{(r_{\mathrm{up}}\left(t\right)+r_{\mathrm{low}}\left(t\right))}{2}$

(3)将d(t)视为新的x(t)重复(1)和(2)，直到d(t)满足固有模态函数条件，记为c₁(t)＝d(t)，c₁(t)即为一个固有模态函数，

r(t)＝x(t)-c₁(t)

r(t)即为新的x(t)；

(4)重复(1)、(2)和(3)，直到残差为零或单调函数时停止，最后得到关于x(t)的表达式：

$x\left(t\right)\text{=}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_j\left(t\right)+r_N\left(t\right)$

式中，N代表分解得到的固有模态函数的个数，残差r_N(t)即为窗内各独立成分的趋势项。

3.根据权利要求1所述的基于独立成分分析和希尔伯特-黄变换融合的自动眼电干扰去除方法，其特征在于，步骤五的判断其中一类的独立成分是否包含眼电信号的计算过程如下所示：

$J=\frac{J_b}{J_w}$

其中J_b和J_w分别表示类内和类间的均方距离，计算公式如下所示，

$J_w=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{P}_j{J}_j$

其中，T表示矩阵转置，P_j表示样本属于第j类的先验概率，通过第j类中样本数量占样本总数量的比值来估计，为第j类的样本均值，为所有样本的均值，J_j尽第j类样本的均方距离，计算公式如下所示，

其中，N_j是第j类中样本的数量，是第j类的第i个样本向量，设定一个阈值J_th，为了确保眼电信号检测的准确率，J_th＝6.0，如果J＞J_th，则表示聚类均值大的一类中的独立成分包含眼电信号，否则，该段不包含眼电信号；对所有判断为包含眼电信号的段以及所属独立成分进行统计。