基于多通道的脑电信号中肌电噪声的消除方法与流程

本发明属于脑电信号处理技术领域，具体涉及一种基于总体平均经验模态分解和典型相关分析，从多通道脑电信号中自动识别肌电噪声并消除的方法，主要应用于人脑相关疾病和人脑功能的研究。

背景技术：

脑电信号是脑神经细胞的电生理活动在大脑皮层或头皮表面的总体反映。脑电图作为记录脑电信号的设备，已经和心电图、X射线检查一样，成为临床的重要检查手段。然而脑电信号作为微弱的电生理信号，经常受到如心电、眼电和肌电等多种噪声的干扰，影响后续对脑电分析的准确性。由于肌肉活动(如咬、嚼和皱眉)引起的肌电干扰，在脑电信号采集过程是难以避免的，并且其幅值大、频域分布广，造成肌电噪声是最难消除的干扰源。

过去通常用低通滤波器来去除肌电干扰。然而，若肌电干扰与感兴趣的脑电信号的频谱重叠，频率滤波器不仅会抑制肌电干扰，而且可能会滤掉有价值的脑电信息。

后来一些学者提出用独立成分分析来解决脑电中肌电噪声消除问题。独立成分分析(ICA)利用独立性将脑电信号分解成相互统计独立的分量。独立成分分析一般通过人为观察判定独立分量中是否包含肌电，置零这些含有肌电的分量后重建得到干净的脑电信号。实验证明独立成分分析在去除眼电和心电噪声具有良好的效果，然而在去除肌电噪声的效果并不是很理想，这是因为通过独立成分分析得到的大部分独立分量中既包含脑电也包含肌电。

为此，一些学者提出用典型相关分析(Clercq WD,VergultA,Vanrumste B,Van Paesschen W,Van Huffel S,Canonical correlation analysis applied to remove muscle artifacts from the electroencephalogram.IEEE Transactions on Biomedical Engineering,2006,53(12):2583-2587.)来解决脑电中肌电噪声消除问题。典型相关分析利用自相关性将脑电信号分解成互不相关的分量，通过求取分量的自相关系数，若低于设定阈值，则这些分量被判定是肌电噪声，置零这些分量后重建得到干净的脑电信号。由于肌电噪声的特性与白噪声相类似，因此肌电噪声相比脑电信号具有较低的自相关性，典型相关分析能把肌电噪声集中于少数典型变量中，通过设置自相关系数阈值，可以实现去除肌电干扰的目的。实验证明典型相关分析比独立成分分析具有更好的去肌电噪声效果。然而，独立成分分析和典型相关分析等方法将脑电信号分解成多个源分量后，被判定为肌电分量的源内仍然包含脑电成分，去除这些源分量虽然可以达到去噪的目的，然而不可避免地丢失了部分脑电信息。

近期，有学者提出了一种基于总体平均经验模态分解和独立成分分析的新型多通道去噪方法(KeZeng,DanChen,GaoxiangOuyang,LizheWang,XianzengLiu,XiaoliLi,AnEEMD-ICA approach to enhancing artifact rejection for noisy multivariate neural data.IEEE transactions on NeuralSystemsandRehabilitationEngineering,2016,24(6):630-638.)。实验证明该方法在多通道脑电信号去除肌电噪声方面，比独立成分分析具有更好的去噪效果。然而，该方法由于采用了独立成分分析进行盲源分离，其得到的独立分量内往往包含肌电和脑电二者信息，将这些分量置零将不可避免地丢失部分脑电信息。

技术实现要素：

本发明为了克服现有技术的不足之处，提出一种基于多通道的脑电信号中肌电噪声的消除方法，以期能去除肌电噪声对脑电信号的影响，并保留各本征模态分量中疑似脑电成分的信息不丢失，从而提高脑电信号分析的准确性。

本发明为解决技术问题，采用如下技术方案：

本发明一种基于多通道的脑电信号中肌电噪声的消除方法的特点是包括如下步骤：

步骤一：由脑电测量设备采集并记录t时刻N通道的脑电信号，记为：X(t)＝[x₁(t),x₂(t),…,x_n(t),…,x_N(t)]^T，x_n(t)为t时刻第n通道的脑电信号，T为矩阵的转置；1≤n≤N；

步骤二：应用总体平均经验模态分解将所述第n通道的脑电信号x_n(t)分解为P个本征模态分量，记为：I_n(t)＝[i_1(n)(t),i_2(n)(t),…,i_p(n)(t),…,i_P(n)(t)]^T；i_p(n)(t)为t时刻第n通道的脑电信号x_n(t)的第p个本征模态分量；1≤p≤P；从而获得t时刻N通道的脑电信号X(t)的本征模态分量矩阵，记为：I(t)＝[I₁(t),I₂(t),…,I_n(t),…,I_N(t)]^T；

步骤三：求取所述第n通道的脑电信号x_n(t)的第p个本征模态分量i_p(n)(t)的自相关系数值R_p(n)，当所述自相关系数R_p(n)低于阈值θ时，判定所述第p个本征模态分量i_p(n)(t)为含有肌电噪声的本征模态分量；从而从所述本征模态分量矩阵I(t)中挑选出所有含有肌电噪声的的本征模态分量，并组成含有肌电噪声的本征模态分量矩阵，记为M(t)＝[m₁(t),m₂(t),…,m_B(t)]^T；B表示含有肌电噪声的的本征模态分量的总数；

步骤四：用典型相关分析对所述含有肌电噪声的本征模态分量矩阵M(t)进行盲源信号的分离，得到混合矩阵A、解混矩阵W和源信号矩阵Y(t)＝[y₁(t),y₂(t),…,y_b(t),…,y_B(t)]^T；y_b(t)表示第b个典型变量，并有：M(t)＝AY(t)或Y(t)＝WM(t)；1≤b≤B；

步骤五：求取所述源信号矩阵Y(t)中的第b个典型分量y_b(t)的自相关系数值r_b，当所述自相关系r_b低于所设定的阈值e时，判定所述第b个典型分量y_b(t)为含有肌电噪声的典型分量；并将含有肌电噪声的典型分量置为零；从而将所述源信号矩阵Y(t)中所有含有肌电噪声的典型分量均置为零，得到不含有肌电噪声的源信号矩阵

步骤六：利用式(1)得到不含有肌电噪声的本征模态分量矩阵

$\tilde{M}\left(t\right)=A\tilde{Y}\left(t\right)---\left(1\right)$

步骤七：将所述不含有肌电噪声的本征模态分量矩阵中每个本征模态分量按照各自在所述本征模态分量矩阵I(t)中挑选前的位置，替换所述本征模态分量矩阵I(t)中对应的本征模态分量；从而得到去除噪后的本征模态分量矩阵I′(t)＝[I′₁(t),I′₂(t),…,I′_n(t),…,I′_N(t)]^T；

步骤八：利用式(2)得到去除噪后的第n通道的干净脑电信号从而获得去除噪后的N通道的的脑电信号

${\tilde{x}}_n\left(t\right)=\underset{p=1}{\overset{P}{\Σ}}{i}_{p\left(n\right)}^{\′}\left(t\right)---\left(2\right)$

式(2)中，i′_p(n)(t)表示第n通道的脑电信号x_n(t)的去除噪后的第p个本征模态分量。

本发明与传统多通道去噪方法相比，不仅能去除肌电对脑电的影响，同时能减少脑电信息在处理过程中的丢失；具体的有益效果体现在：

1、本发明步骤二和步骤三，对每一通道脑电信号进行总体经验平均模态分解，得到的多个通道的本征模态分量，通过求取各个本征模态分量的自相关系数，把自相关系数低于设定阈值的本征模态分量构成含噪本征模态分量矩阵，进行下一步处理。这种方式保留了只含脑电的本征模态分量，即保留了其中的脑电信息。与典型相关分析和独立成分分析对全部脑电处理的方式相比，能够更好地保留原始脑电，减少了因去噪过程带来的脑电信息丢失。

2、本发明步骤四、步骤五和步骤六，使用典型相关分析进行盲源信号分离，能够将肌电集中在少数的典型分量中。而现有基于总体平均经验模态分解和独立成分分析(EEMD-ICA)的多通道去噪方法则通过独立成分分析进行盲源信号分离，因为独立成分分析得到的独立分量中往往既包含肌电也包含脑电，置零含噪独立分量会损失脑电，保留的独立分量中会包含肌电。因此本发明相比于EEMD-ICA多通道去噪方法，不仅能更好地去除肌电，而且能更好地保留脑电信息。

附图说明

图1为本发明方法的主流程图；

图2a为19通道的干净模拟脑电信号示意图；

图2b为被肌电噪声干扰的19通道的模拟脑电信号(λ＝1)示意图；

图2c为通过本发明方法得到的部分含噪的本征模态分量示意图；

图2d为通过本发明方法得到的部分典型分量示意图；

图2e为通过本发明方法得到的消除肌电噪声后重建的脑电信号示意图；

图3a为19通道的模拟脑电信号示意图；

图3b为被肌电噪声干扰的19通道的模拟脑电信号(λ＝3)示意图；

图3c为通过本发明方法处理模拟脑电，所得到的去除肌电噪声后重建的脑电信号示意图；

图3d为通过EEMD-ICA多通道处理方法处理模拟脑电，所得到的去除肌电噪声后重建脑电信号示意图；

图3e为本发明方法与EEMD-ICA方法处理模拟脑电，其去噪性能指标相对均方根的比较图；

图3f为本发明方法与EEMD-ICA方法处理模拟脑电，其去噪性能指标相关系数的比较图；

图4a为19通道的干净真实脑电信号示意图；

图4b为被肌电噪声干扰的19通道混合的真实脑电信号示意图；

图4c为通过本发明方法处理真实脑电，所得到的去除肌电噪声后重建的脑电信号示意图；

图4d为通过EEMD-ICA的多通道处理方法处理真实脑电，所得到的去除肌电噪声后重建脑电信号示意图；

图4e为本发明方法与EEMD-ICA方法处理真实脑电，其去噪性能指标相对均方根误差的比较图；

图4f为本发明方法与EEMD-ICA方法处理真实脑电，其去噪性能指标相关系数的比较图；

图5a为21通道包含肌电和眼电噪声的癫痫脑电信号示意图；

图5b为通过本发明方法所得到的去除肌电噪声后重建的脑电信号示意图。

具体实施方式

如图1所示，基于多通道的脑电信号中肌电噪声的消除方法是：首先用总体平均经验模态对每一通道脑电信号进行分解，得到每一通道的本征模态分量；再通过自相关系数判定含噪声的本征模态分量，由含噪声的本征模态分量构成含噪本征模态分量矩阵；然后对含噪本征模态分量矩阵进行盲信号分离；最后用自相关系数判定含噪声的典型分量，置零噪声分量并重建信号。

下面分别通过纯模拟脑电信号、半模拟脑电信号与实测脑电信号为例，结合附图来说明具体的实施方式。

1.纯模拟脑电信号

在这一部分，将说明两个实例，第一个实例是将介绍本发明的具体实施方式，第二个实例是将本发明与EEMD-ICA多通道处理方法相比较。

(1)实例一

步骤一：由脑电测量设备采集并记录t时刻受肌电噪声干扰的N＝19个通道脑电信号X(t)＝[x₁(t),x₂(t),…,x_n(t),…,x₁₉(t)]^T，1≤n≤N；其中，x_n(t)为t时刻第n通道的脑电信号，T为矩阵的转置；如图2b所示的X(t)＝X_EEG(t)+λX_EMG(t)，λ＝1；且如图2a所示的X_EEG(t)＝[x_EEG1(t),x_EEG2(t),…,x_EEG19(t)]^T表示19通道的纯模拟脑电信号；X_EMG(t)＝[x_EMG1(t),x_EMG2(t),…,x_EMG19(t)]^T表示19通道的模拟肌电信号。模拟脑电X_EEG(t)每一个通道都是由5个2s长的EEG片段连接而成，而每一个EEG片段都是由4个频率在4-30Hz范围内的正弦信号叠加而成，并且这4个频率是随机产生的。而EMG信号则是由一段在20-60Hz频率范围内的19通道白噪声生成。信号采样频率为250Hz，采取共10秒的模拟信号，即总共有2500个采样点；

步骤二：应用总体平均经验模态分解将混合脑电信号X(t)内每一通道x_n(t)分解成为p＝11个本征模态分量I_n(t)＝[i_1(n)(t),i_2(n)(t),…,i_p(n)(t),…,i_11(n)(t)]^T；i_p(n)(t)为t时刻第n通道的脑电信号x_n(t)的第p个本征模态分量；1≤p≤P；对19个通道的脑电信号分别应用总体平均经验模态分解后得到19通道的本征模态分量矩阵I(t)＝[I₁(t),I₂(t),…,I_n(t),…,I₁₉(t)]^T；

步骤三：求取第n通道的脑电信号x_n(t)的第p个本征模态分量i_p(n)(t)的自相关系数值R_p(n)，设定阈值θ＝0.995，当自相关系数R_p(n)低于阈值θ时，判定第p个本征模态分量i_p(n)(t)为含有肌电噪声的本征模态分量；从而从本征模态分量矩阵I(t)中挑选出所有含有肌电噪声的的本征模态分量，并组成含有肌电噪声的本征模态分量矩阵，记为M(t)＝[m₁(t),m₂(t),…m_b(t),…,m₁₂₇(t)]^T；B＝127表示含有肌电噪声的的本征模态分量的总数；其中部分含噪本征模态分量如图2c；

步骤四：用典型相关分析对含有肌电噪声的本征模态分量矩阵M(t)进行盲源信号的分离，得到混合矩阵A、解混矩阵W和源信号矩阵Y(t)＝[y₁(t),y₂(t),…,y_b(t),…,y_B(t)]^T；y_b(t)表示第b个典型变量，并有：M(t)＝AY(t)或Y(t)＝WM(t)；1≤b≤B；

本实施例中，将含噪本征模态分量矩阵M(t)进行延时为1的处理，得到2个数据集：M(t)和J(t)＝M(t-1)，然后用典型相关分析对这两个数据集进行盲源信号的分离，由此得到M(t)的混合矩阵A、解混矩阵W和127个典型变量构成的源信号矩阵Y(t)＝[y₁(t),y₂(t),…,y_b(t),…,y₁₂₇(t)]^T，其中部分典型变量如图2d；

步骤五：求取源信号矩阵Y(t)中的第b个典型分量y_b(t)的自相关系数值r_b，当自相关系r_b低于所设定的阈值e时，判定第b个典型分量y_b(t)为含有肌电噪声的典型分量，此处e＝0.95；并将含有肌电噪声的典型分量置为零；从而将源信号矩阵Y(t)中所有含有肌电噪声的典型分量均置为零，得到不含有肌电噪声的源信号矩阵

步骤六：利用式(1)得到不含有肌电噪声的本征模态分量矩阵

$\tilde{M}\left(t\right)=A\tilde{Y}\left(t\right)---\left(1\right)$

步骤七：将不含有肌电噪声的本征模态分量矩阵中每个本征模态分量按照各自在本征模态分量矩阵I(t)中挑选前的位置，替换本征模态分量矩阵I(t)中对应的本征模态分量；从而得到去除噪后的本征模态分量矩阵I′(t)＝[I′₁(t),I′₂(t),…,I′_n(t),…,I′₁₉(t)]^T；

步骤八：利用式(2)得到去除噪后的第n通道的干净脑电信号从而获得去除噪后的N通道的的脑电信号

${\tilde{x}}_n\left(t\right)=\underset{p=1}{\overset{P}{\Σ}}{i}_{p\left(n\right)}^{\′}\left(t\right)---\left(2\right)$

式(2)中，i′_p(n)(t)表示第n通道的脑电信号x_n(t)的去除噪后的第p个本征模态分量。对19通道依次处理后，得到消除噪声的19通道脑电信号如图2e所示；

将去噪后的脑电信号与模拟的干净真实脑电信号X_EEG(t)进行对比，从图2a和图2e中可以清晰地观察到肌电噪声基本被完全消除，并且很好地保留了原干净真实脑电信号的细节信息，说明了本发明对多通道脑电信号中肌电噪声消除的有效性。

(2)实例二

为了量化评估本发明的效果，为此将发明方法与EEMD-ICA多通道处理方法相比较。

首先模拟N＝19的受肌电噪声干扰的多通道纯模拟脑电信号X₁(t)＝[x₁(t),x₂(t),…,x₁₉(t)]^T，1≤n≤N；其中图3b所示的X₁(t)为X_EEG(t)和X_EMG(t)的混合信号，即X₁(t)＝X_EEG(t)+λX_EMG(t)，图3a所示的X_EEG(t)＝[x_EEG1(t),x_EEG2(t),…,x_EEG19(t)]^T和X_EMG(t)＝[x_EMG1(t),x_EMG2(t),…,x_EMG19(t)]^T分别表示模拟的19通道干净脑电信号和肌电噪声。

信号采样频率为250Hz，采取共10秒的模拟信号，即总共有S₁＝2500个采样点，这里λ代表肌电噪声干扰脑电信号的强度，图3b中λ＝3。按照实例一中的步骤，本发明方法得到如图3c所示的去除肌电噪声重建的脑电信号EEMD-ICA多通道处理方法与EEMD-CCA多通道处理方式类似，两种方法设定相同的阈值。图3d所示为EEMD-ICA多通道处理方法去除肌电噪声重建的脑电信号

为了有效评估EEMD-CCA和EEMD-ICA两种算法的去噪效果，选择两个性能指标作为评价指标，分别是相对均方根误差(RRMSE)和相关系数(CC)。

信噪比定义：SNR＝RMS(X_EEG)/RMS(λ·X_EMG)，其中RMS表示均方根，和可以看出λ与信噪比SNR成反比，λ越大，信噪比越低。

此外，相对均方根误差定义如下：RRMSE值越小表明去噪效果越好。

第二个性能指标是相关系数CC。此处使用相关系数作为性能指标，可直观展示EEG信号去噪前后的结构相似度，相关系数的值越高说明数据恢复得越好。图3e、图3f展示了两种方法在不同信噪比下的消除肌电噪声和保留脑电的效果，从中均可以看出在本发明方法在不同信噪比稳定优于EEMD—ICA处理方法。

2.半模拟脑电信号

为了进一步验证本发明比EEMD-ICA具有更好的去噪效果，将使用由真实脑电信号和真实肌电信号混合得到的脑电信号，比较两种方法的去噪效果。

首先模拟N＝19的受肌电噪声干扰的多通道半模拟脑电信号X₂(t)＝[x₁(t),x₂(t),…,x₁₉(t)]^T，其中如图4b所示，λ＝1.5时，X₂(t)为X_EEG(t)和X_EMG(t)的混合信号：X₂(t)＝X_EEG(t)+λX_EMG(t)，如图4a所示的X_EEG(t)和X_EMG(t)分别表示模拟的19通道真实干净脑电信号和肌电噪声，且X_EEG(t)＝[x_EEG1(t),x_EEG2(t),…,x_EEG19(t)]^T，X_EMG(t)＝[x_EMG1(t),x_EMG2(t),…,x_EMG19(t)]^T。

信号采样频率为1000Hz，采取共10秒的模拟信号，即总共有S₂＝10000个采样点。按照实例一中的步骤，本发明方法得到如图4c所示的去除肌电噪声重建的脑电信号图4d所示为EEMD-ICA多通道处理方法去除肌电噪声重建的脑电信号用性能指标RRMSE和CC作为评价指标，图4e和图4f展示了两种方法在不同信噪比下的消除肌电噪声和保留脑电的效果，从中均可以看出在本发明方法在不同信噪比稳定优于EEMD-ICA处理方法。

3.实测脑电信号

在本实例的第三部分，使用真实脑电数据作为实验对象，使用EEMD-CCA多通道算法进行处理，评判发明方法的去噪效果。图5a是一段21通道的真实癫痫发作脑电信号，该信号的采样频率是250Hz，采样时间为10s，总共2500个采样点。从图5a中可以看出该段信号被眼电和肌电两种噪声所干扰。眼电可以在Fp1和Fp2通道的2.5、3.5、6和7.5s附近被明显观测到。同样肌电干扰可以在F7、T3、T5、C3、T1通道的0s-3.9s处以及F8、T4、F4、C4、P4的5s-10s处被观测到。癫痫发作可在通道T2、F8、T4、T6被观察到，一部分发作区域被肌电噪声严重干扰，这会影响到后续脑电信号的分析和对癫痫发作区域定位工作，因此消除肌电噪声干扰是十分有必要的。

按照实例一的步骤对上述真实脑电信号进行处理，使用本发明方法得到的去除肌电噪声后重建的脑电信号如图5b。通过对比图5a、图5b，可以发现本发明不仅将肌电噪声消除的非常干净，还能完整地保留脑电信号中的关键细节信息，例如，在通道F8、T4、T6中被肌电噪声干扰的癫痫发作脑电部分已经去除了肌电噪声，而且其中的眼电噪声被完整的保存下来。

综上，本发明能够解决多通道情况下肌电噪声消除的难题，另外本发明是一种完全自动的监测并消除肌电噪声的方法，避免了人为干预造成的影响。该方法适用于临床诊断和神经科学研究的多通道脑电设备中，与EEMD-ICA的多通道脑电处理方法相比，能够取得更好的去噪效果，对进一步研究大脑真实的电生理活动具有重要意义。