CN108309290A - 单通道脑电信号中肌电伪迹的自动去除方法 - Google Patents

Abstract

本发明公开了一种用于单通道脑电信号中肌电伪迹的去除方法，将脑电信号通过SSA分解得到P个信号分量；将P个信号分量按行拼接成一个P维数据矩阵；将P维数据矩阵进行时间延迟处理得到若干个数据矩阵；利用MCCA对若干个数据矩阵进行盲源分离，得到源估计矩阵S和混合矩阵A；识别源估计矩阵中与肌电伪迹相关的源；去除源估计矩阵中的肌电伪迹，将识别为肌电伪迹的源置零，得到消除肌电伪迹后的源估计矩阵S’，并通过重构得到去除肌电伪迹后的多通道脑电信号X′＝A*S′；将多通道脑电信号X’的各行求和，即可最终得到去除了肌电伪迹后的单通道脑电信号x’。本发明去除肌电伪迹的同时尽可能保留了脑电信息，提高了脑电信号分析的准确性。

Description

单通道脑电信号中肌电伪迹的自动去除方法

技术领域

本发明涉及脑电信号处理技术领域，特别涉及一种单通道脑电信号中肌电伪迹的自动去除方法。

背景技术

脑电图(Electroencephalograph,EEG)是大脑神经细胞的电生理活动在头皮上的反映。由于具有高时间分辨率、非侵入式、低成本和适合长期监测等特点，脑电图已被广泛用于研究脑功能和病理性脑机制。但是EEG在记录大脑活动的同时，也记录了大脑以外的其他活动产生的电信号，而这些被记录的不是由大脑活动产生的电信号称为伪迹。常见的伪迹包括眼电(Electrooculogram，EOG)、肌电(Electromyography，EMG)、心电(Electrocardiography，ECG)以及由电路或外部设备带来的噪声等。相比于其他类型的伪迹信号，比如心电(ECG)和眼电(EOG)，头部肌肉收缩(即肌电(EMG)信号)对EEG产生的干扰更难以消除，其主要原因可能在于肌电伪迹具有高振幅、宽频谱以及生理分布更广等特点，导致一些简单的伪迹去除方法失效。伪迹的存在对后续脑电信号分析造成了极大的影响，可能会导致重要信息的丢失。因此，伪迹去除是神经信息处理应用中最重要的预处理步骤之一，对神经科学研究以及实际临床诊断都具有重要意义。

目前的肌电伪迹去除方法主要分为两类：伪迹拒绝和伪迹校正。伪迹拒绝是将包含伪迹的脑电时段整个去除，这种方法的主要缺点是它也会消除重要的脑电信息，导致数据丢失；而伪迹校正是指采用各种方法去除伪迹成分的同时尽可能的保留脑电信息。伪迹校正技术主要包括：1)滤波，通过使用特定频段的滤波器来去除脑电信号中的肌电伪迹；然而由于肌电信号频谱分布范围广，和脑电信号有较大的重合，所以在去除肌电伪迹的同时也会损失一部分脑电信号。2)小波变换(wavelet transform，WT)，通过小波变换将观测信号分解为一系列基函数，然后对系数进行阈值处理以去除肌电伪迹，最后通过小波重构得到纯净的EEG信号；但是小波变换方法需要大量的实验来选择合适的小波基函数和分解层次。3)盲源分离(Blind Source Separation，BSS)，盲源分离类算法用于从多个观测到的混合信号中分析出未知的信号源以及混合矩阵，然后将伪迹相关的信号源去除，最后通过混合矩阵就可以得到纯净信号；常见的BBS算法有：主成分分析(Principal componentanalysis，PCA)、独立成分分析(Independent component analysis，ICA)和典型相关分析(canonical correlation analysis，CCA)。4)经验模态分解(empirical modedecomposition，EMD)，通过EMD将观测到的EEG信号分解得到若干本征模式分量(IMF)，然后将与肌电伪迹相关的IMF去除再重建信号即可得到纯净EEG。EMD是一种数据驱动的方法，用于对非平稳、非线性、随机过程进行处理，因此非常适用于脑电信号分析和处理，但是EMD计算复杂度高，一般需要耗时数分钟，不太适合实时处理。

随着近年来可穿戴设备的快速发展以及未来个人医疗的需求，如何有效去除单通道脑电数据中的伪迹信号变得尤为重要。

发明内容

本发明是为了克服现有技术的不足，提供一种单通道脑电信号中肌电伪迹的自动去除方法，从观测脑电信号中去除肌电伪迹的同时尽可能保留脑电信息，避免脑电信息的丢失，从而提高脑电信号分析的准确性。

本发明采用如下技术方案：单通道脑电信号中肌电伪迹的自动去除方法，包括以下步骤：

步骤一：将通过单通道脑电电极传感器采集得到的脑电信号通过奇异谱分析算法SSA分解得到P个信号分量；

步骤二：将步骤一中得到的P个信号分量按行拼接成一个P维数据矩阵；

步骤三：将拼接成的P维数据矩阵进行时间延迟处理，得到若干个数据矩阵；

步骤四：利用多重集典型相关分析MCCA对步骤三中得到的若干个数据矩阵进行盲源分离，得到源估计矩阵S和混合矩阵A；

步骤五：识别源估计矩阵中与肌电伪迹相关的源；

步骤六：去除源估计矩阵中的肌电伪迹，将识别为肌电伪迹的源置零，得到消除肌电伪迹后的源估计矩阵S’，并通过重构得到去除肌电伪迹后的多通道脑电信号X′＝A*S′；

步骤七：将多通道脑电信号X’的各行求和，即可最终得到去除了肌电伪迹后的单通道脑电信号x’。

从以上技术方案可知，本发明首先将混合了肌电伪迹的单通道脑电信号经过奇异谱分析(singular spectrum analysis，SSA)算法分解得到若干信号分量；接着将得到的所有信号分量拼接成一个数据矩阵并进行延迟处理得到若干个数据矩阵；然后利用多重集典型相关分析(MCCA)进行盲源分离，得到若干典型变量；最后，将自相关系数小于选定阈值的分量当成是肌电伪迹，并重构得到去除肌电伪迹后的干净脑电信号。与已有技术相比，本发明所取得的有益效果包括：

1、使用基于奇异谱分析(SSA)算法和多重集典型相关分析(MCCA)算法的肌电伪迹去除方法，无需参考信号，只需单通道脑电信号即可有效去除肌电伪迹，同时保留大量有用的脑电信号，适用于近年来快速发展的穿戴式单通道脑电设备，并且去除过程中无需人工参与即可自动去除肌电伪迹。

2、相比基于EMD的方法，其运行时间更短，具有较佳的实时性，对日常个人脑电分析技术具有重要意义。

附图说明

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

图2为半模拟脑电信号图，其中(a)为真实的干净脑电信号，(b)为模拟的肌电伪迹信号，(c)为两者的合成信号；

图3为利用MCCA进行分解后得到的估计源信号；

图4为原真实的干净脑电信号与采用本发明方法去除肌电伪迹后脑电信号的对比图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明，但本发明的实施方式不限于此。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些详细描述在这里将被忽略。

本实施例以半模拟脑电信号为例，其流程图参见图1，具体包括以下步骤：

步骤一：将通过单通道脑电电极传感器采集得到的脑电信号通过奇异谱分析(SSA)算法分解得到P个信号分量；

具体为：设所采集到的单通道含肌电伪迹的脑电信号为x＝[x(1),x(2),...,x(T)]，其中序列x(t)(t＝1,2,…,T)表示脑电信号x在t时刻的值，脑电信号x经SSA分解后产生P个信号分量，即其中每个分量x_i的长度也都为T。

本步骤一中利用SSA算法将单通道含肌电伪迹的脑电信号分解为P个信号分量的过程为：

(1)计算轨迹矩阵(trajectory matrix)X

对于给定的一维时间序列x＝(x₁,x₂,…,x_N)，根据窗口长度L计算其轨迹矩阵X，轨迹矩阵X为L x(N‐L+1)阶：

(2)奇异值分解(SVD)

定义矩阵S＝XX^T，计算矩阵S的特征值和特征变量，并将特征值按降序排列：λ₁,λ₂,…,λ_N(λ₁≥λ₂≥,…,≥λ_N≥0)，则其对应的特征向量为u₁,u₂,…,u_L。我们因此可以将轨迹矩阵X表示为d个初等矩阵的和：X＝Ｘ₁+X₂+…+X_d，其中d＝argmax_i{λ_i＞0}，

(3)分组

将初等矩阵X_i的下标{1,2,…,d}分成p个不相交的子集I₁,I₂,…,I_p，设I＝{i₁,i₂,…,i_m}，那么合成矩阵计算集合I₁,I₂,…,I_p的每个合成矩阵，那么轨迹矩阵X可表示为：

(4)对角平均

将步骤(3)中分组得到的每个矩阵利用对角平均公式转换为长度为N的新序列则初始序列x被分解为p个序列之和：初始序列x指单通道含肌电伪迹的脑电信号。

本实施例中，设x(t)＝x_EEG(t)+x_EMG(t)，这里x_EEG(t)和x_EMG(t)分别代表真实采集的不含伪迹的“纯净”EEG信号和模拟合成的EMG伪迹信号，信号采样频率为200Hz，信号长度为3s，信噪比为‐5dB。图2(a)、图2(b)、图2(c)分别为x_EEG(t)、x_EMG(t)和混合信号x(t)的波形。然后设定窗口长度L＝8，利用SSA算法分解x(t)得到8个信号分量。

步骤二：将步骤一中得到的P个信号分量按行拼接成一个P维数据矩阵；

具体为：将P个信号分量x_i(i＝1,2,…,P)拼接成一个P x T的数据矩阵X，即X＝[x₁,x₂,…,x_P]^T；本实施例中，步骤一分解得到8个信号分量，本步骤将8个信号分量按行拼接成一个8x 600的数据矩阵X(t)。

步骤三：将拼接成的P维数据矩阵进行时间延迟处理，得到若干个数据矩阵；

具体为：将数据矩阵X进行时间延迟处理得到M个数据矩阵：

X_m(t)＝X(t-(m-1)τ),m＝1,2,…,M.

本实施例将数据矩阵X(t)进行时间延迟处理，单位延迟时间τ＝70，得到M＝3个数据矩阵：

X₁(t)＝X(t),X₂(t)＝X(t-70),X₃(t)＝X(t-140)

步骤四：利用多重集典型相关分析(MCCA)对步骤三中得到的若干个数据矩阵进行盲源分离，得到源估计矩阵和混合矩阵；

具体为：利用MCCA对步骤三得到的M个数据矩阵X_m(m＝1,2,…,M)进行盲源分离，得到源估计矩阵S和混合矩阵A。本实施例中，M＝3，所获得的源估计矩阵中各个源的波形如图3所示。

本步骤利用MCCA进行盲源分离的过程为：

对于给定的M个P维数据矩阵X_m(m＝1,2,…,M)，每个数据矩阵可表示为X_m(t)＝A_my_m(t)或y_m(t)＝W_mX_m(t)，其中A_m＝W_m ^-1，这里A_m、W_m和y_m(t)分别称为X_m(t)的混合矩阵、解混矩阵和典型变量矩阵，解混矩阵可表示为通过最大化M个P维数据矩阵中相应盲信号源的总体相关性来实现各个数据矩阵的盲源分离，其目标函数分为P步：

第一步：

第二步到第P步：受限于

其中， 为的标准差。

通过以上P步，就可得到数据矩阵X_m(t)的混合矩阵A_m、解混矩阵W_m和典型变量矩阵y_m(t)，其中m＝1,2,…,M。

步骤五：识别源估计矩阵中与肌电伪迹相关的源；

具体为：计算源估计矩阵S中每个源s_i(t)的自相关系数，将自相关系数小于选定阈值的源判定为肌电伪迹，因为肌电信号的自相关系数一般较小。其中自相关系数的计算公式为：

本实施例中，发现s₂(t)、s₃(t)、s₄(t)、s₆(t)、s₇(t)、s₈(t)的自相关系数小于阈值0.8，故将这些源判定为肌电伪迹。

步骤六：去除源估计矩阵中的肌电伪迹，并通过重构得到去除肌电伪迹后的多通道脑电信号；具体为：将步骤五中被判定为肌电伪迹的源置零，得到消除肌电伪迹后的源估计矩阵S’，然后通过步骤四中得到的混合矩阵A获得去除肌电伪迹的多通道脑电信号X’，即X′＝A*S′。

本实施例中，将源s₂(t)、s₃(t)、s₄(t)、s₆(t)、s₇(t)、s₈(t)置零，以得到消除肌电伪迹后的源估计矩阵S’。

步骤七：将多通道信号X’的各行求和，即可最终得到去除了肌电伪迹后的单通道脑电信号x’，如图4中第二个波形所示，第一个波形为x_EEG(t)。从图中可以明显地观察到原肌电伪迹x_EMG(t)基本被完全去除，并且很好地保留了原干净真实脑电信号x_EEG(t)的细节信息，说明了本发明对单通道脑电信号中肌电伪迹去除的有效性。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.单通道脑电信号中肌电伪迹的自动去除方法，其特征在于，包括以下步骤：

步骤一：将通过单通道脑电电极传感器采集得到的脑电信号通过奇异谱分析算法SSA分解得到P个信号分量；

步骤二：将步骤一中得到的P个信号分量按行拼接成一个P维数据矩阵；

步骤三：将拼接成的P维数据矩阵进行时间延迟处理，得到若干个数据矩阵；

步骤四：利用多重集典型相关分析MCCA对步骤三中得到的若干个数据矩阵进行盲源分离，得到源估计矩阵S和混合矩阵A；

步骤五：识别源估计矩阵中与肌电伪迹相关的源；

步骤六：去除源估计矩阵中的肌电伪迹，将识别为肌电伪迹的源置零，得到消除肌电伪迹后的源估计矩阵S’，并通过重构得到去除肌电伪迹后的多通道脑电信号X′＝A*S′；

步骤七：将多通道脑电信号X’的各行求和，即可最终得到去除了肌电伪迹后的单通道脑电信号x’。

2.根据权利要求1所述的单通道脑电信号中肌电伪迹的自动去除方法，其特征在于，步骤一中利用SSA将脑电信号分解为P个信号分量的过程为：

(1)计算轨迹矩阵X

对于给定的一维时间序列x＝(x₁,x₂,…,x_N)，根据窗口长度L计算其轨迹矩阵X，轨迹矩阵X为L x(N‐L+1)阶：

(2)奇异值分解

定义矩阵S＝XX^T，计算矩阵S的特征值和特征变量，并将特征值按降序排列：λ₁,λ₂,…,λ_N，其中λ₁≥λ₂≥,…,≥λ_N≥0，则其对应的特征向量为u₁,u₂,…,u_L；将轨迹矩阵X表示为d个初等矩阵的和：X＝X₁+X₂+…+X_d，其中d＝argmax_i{λ_i＞0}，

(3)分组

将初等矩阵X_i的下标{1,2,…,d}分成p个不相交的子集I₁,I₂,…,I_p，设I＝{i₁,i₂,…,i_m}，合成矩阵计算集合I₁,I₂,…,I_p的每个合成矩阵，轨迹矩阵X表示为：

(4)对角平均

将步骤(3)中分组得到的每个矩阵利用对角平均公式转换为长度为N的新序列则初始序列x被分解为p个序列之和：初始序列x指单通道含肌电伪迹的脑电信号。

3.根据权利要求1所述的单通道脑电信号中肌电伪迹的自动去除方法，其特征在于，步骤四中利用MCCA进行盲源分离的过程为：

对于给定的M个P维数据矩阵X_m，每个数据矩阵可表示为X_m(t)＝A_my_m(t)或y_m(t)＝W_mX_m(t)，其中A_m＝W_m ^-1，这里A_m、W_m和y_m(t)分别称为X_m(t)的混合矩阵、解混矩阵和典型变量矩阵，解混矩阵表示为通过最大化M个P维数据矩阵中相应盲信号源的总体相关性来实现各个数据矩阵的盲源分离，其目标函数分为P步：

第一步：

第二步到第P步：受限于

其中， 为X_m(t)的标准差；

通过以上P步，得到数据矩阵X_m(t)的混合矩阵A_m、解混矩阵W_m和典型变量矩阵y_m(t)，其中m＝1,2,…,M。

4.根据权利要求1所述的单通道脑电信号中肌电伪迹的自动去除方法，其特征在于，步骤三中将P维数据矩阵X进行时间延迟处理，得到M个数据矩阵：

X_m(t)＝X(t-(m-1)τ),m＝1,2,…,M

其中，τ为单位延迟时间。

5.根据权利要求1所述的单通道脑电信号中肌电伪迹的自动去除方法，其特征在于，步骤五的过程为：计算源估计矩阵S中每个源s_i(t)的自相关系数，其中自相关系数的计算公式为：

将自相关系数小于选定阈值的源判定为肌电伪迹。