CN111631710B - 一种状态相关的动态脑电信号中肌电伪迹的消除方法 - Google Patents

Abstract

本发明公开了一种状态相关的动态脑电信号中肌电伪迹的消除方法，其步骤包括：1、首先将收集到的脑电观测信号通过延时构造两个数据集；2、利用本发明提出的隐马尔科夫独立向量分析法进行动态的联合盲源分离，得到每个数据集在各个状态下的源信号矩阵和解混矩阵；3、选择脑电信号相对应的源信号矩阵和解混矩阵；4、依照自相关系数排序源信号矩阵中各个独立源成分，选择肌电噪声相关的独立源成分置零；5、盲源分离逆变换得到消除噪声后的干净脑电信号。本发明能在实际的动态环境中去除肌电噪声对脑电信号的影响，同时尽可能地保留脑电活动的信息不丢失，从而提高脑电信号的分析准确性，为脑电去噪提供一种新思路。

Description

一种状态相关的动态脑电信号中肌电伪迹的消除方法

技术领域

本发明涉及医学信号处理领域，尤其涉及复杂环境下神经信号中噪声去除的方法。

背景技术

通过测量神经元活动产生的电信号，脑电图(EEG)可以无创便捷的观测和记录大脑活动，已被广泛地应用于神经科学研究，疾病诊断和医疗健康监测中。然而由于脑电信号振幅较弱，其经常被源自眼，心脏和肌肉活动的各种非神经生理因素污染。在实际的长期医疗监测中，这些不可避免的噪声会明显干扰记录的脑电信号，从而对后续的分析产生不利的影响。因此，在不影响真实脑电数据的情形下，开发有效的噪声去除算法具有重要的意义。相比于眼电和心电噪声，肌电噪声由于其复杂的特性，例如高振幅，宽频谱等，使得对其有效的去除更为困难。

为了从脑电记录中去除肌电噪声，盲信号分离和联合盲信号分离的方法得到了极大的关注并取得了较好的效果。例如，独立成分分析(Independent component analysis，ICA)模型，作为一种盲信号分离方法被广泛的用于脑电去噪。它将脑电数据分解为独立的分量，并把自相关性较低的独立成分视为肌电噪声，将其置零，最后程重建源信号以达到去除噪声的目的。进一步为了满足多数据联合分析的需要，典型相关分析(canonicalcorrelation analysis，CCA)，独立分量分析(Independent vector analysis，IVA)等联合盲信号分离的方法被提出作为脑电去噪的可靠方法。然而，它们中的大多数方法是为理想条件下脑电降噪而设计的，不能满足实际的移动医疗监控中复杂变化的脑电降噪需求。实际情形下，各种肌电噪声伴随着脑电活动随时间发生变化，可能有多种混合模式出现在脑电记录中。而现存的脑电信号消噪方法多适用于静态情形，即模型假定在整个观测过程中，脑电记录的混合模式不发生变化。当在实际的移动医疗监护中应用传统静态脑电消噪方法时，肌电噪声无法被有效去除，从而对后续针对脑电记录的研究造成严重干扰。故发展可以在动态变化情形下消除噪声新方法是十分必要的需求。

发明内容

本发明为了克服现有技术的不足，提出了一种状态相关的动态脑电信号中肌电伪迹的消除方法，以期能在实际的动态环境中去除肌电噪声对脑电信号的影响，同时尽可能地保留脑电活动的信息不丢失，从而提高脑电信号的分析准确性，为脑电去噪提供一种新思路。

本发明为达到上述发明目的，采用如下技术方案：

本发明一种状态相关的动态脑电信号中肌电伪迹消除方法的特点是按如下步骤进行：

步骤一：由脑电测量设备采集并记录N通道的纯净脑电信号矩阵，记为：S_EEG＝[s₁,s₂,...,s_n,...,s_N]∈R^N×T，1≤n≤N，s_n为第n通道的纯净脑电信号，T表示总的采集时长；

采集收缩和舒张拳头时双手前臂的肌电源信号矩阵记为S_EMG∈R^N×T；

利用式(1)构造脑电观测信号矩阵X_EEG＝[x₁,x₂,...,x_n,...,x_N]，x_n为脑电观测信号矩阵X_EEG中第n通道观测信号，1≤n≤N：

X_EEG＝S_EEG+λ·X_EMG (1)

式(1)中，λ表示控制肌电噪声干扰脑电信号的程度；X_EMG∈R^N×T表示肌电观测信号，并有：

X_EMG＝[V₁S_EMG,V₂S_EMG,V₁S_EMG] (2)

式(2)中，V₁∈R^N×N，V₂∈R^N×N表示两个随机混合矩阵；

当时刻t∈[0,t₁]∪[t₂,T]时，肌电观测信号X_EMG由第一随机混合矩阵V₁作用于肌电源信号矩阵S_EMG产生，并表示为第一个混合模式；

当时刻t∈[t₁,t₂]时，肌电观测信号X_EMG由第二随机混合矩阵V₂作用于肌电源信号矩阵S_EMG产生，并表示为第二个混合模式；

利用式(3)得到信噪比SNR：

式(3)中，RMS(·)表示均方根函数；

步骤二：将所述脑电观测信号矩阵X_EEG通过延时K-1个点构造K个数据集X＝[X₁,X₂,...,X_k,...,X_K]，1≤k≤K，其中，X_k表示第k个数据集的脑电观测信号矩阵；

步骤三：利用隐马尔科夫独立向量分析法对K个数据集X进行动态联合盲源分离，得到不同混合模式下第k个数据集X_k的脑电源信号矩阵S_k以及解混矩阵W_k；

步骤四：对应于每一个混合模式，将第k个脑电源信号矩阵S_k中自相关系数较小的通道信号置零，得到肌电噪声分量置零后的脑电源信号矩阵

步骤五：利用式(4)得到第k个消噪后的脑电信号

式(4)中，(W_k)^-1为第k个逆混合矩阵W_k的逆矩阵。

本发明所述的状态相关的动态脑电信号中肌电伪迹消除方法的特点也在于，所述步骤三是按如下过程进行：

步骤3.1、令隐马尔科夫模型的转移状态记为Q＝{q₁,q₂,...,q_m,...,q_M}，q_m为第m个转移状态，1≤m≤M，M为转移状态总数；

K个数据集X实施维特比算法，得到转移状态序列

q_t′为时刻t所处的转移状态；

利用所述转移状态序列

标识出K个数据集X混合模式的变化；

步骤3.2、在每一个混合模式中，对K个数据集X进行联合盲源分离，得到K个解混矩阵W＝[W₁,W₂,...,W_k,...W_K]，W_k为在当前混合模式下第k个数据集X_k的解混矩阵；

利用式(5)得到在当前混合模式下第k个数据集X_k的脑电源信号矩阵S_k：

S_k＝W_k X_k，k＝1,2,...,K (5)。

本发明与传统静态盲分离方法如CCA和ICA，IVA相比的有益效果在于，不仅可以在动态变化环境中去除肌电噪声对脑电的影响，同时能很好地保留脑电信息，与动态算法HMM-ICA相比，也有更好的去噪能力。本发明方法性能的提高体现在：

1、本发明步骤二和步骤三中，通过延时得到了多个数据集，并利用维特比算法对观测数据集做状态划分，对划分的时间段内信号进行联合盲源分离提取噪声源信号。这一做法与传统的去噪方法如ICA，CCA和IVA相比，分段进行的联合盲信号分离满足了实际健康监护中长时间动态监测的需求，解决了复杂环境下变化的噪声去除的问题，提高了本发明方法的精度和噪声去除的准确性。

2、本发明步骤三中，隐马尔科夫的观测部分设置为独立向量分析模型，该方法对多数据集信号进行联合盲源分离时，分离出的源信号在同一数据集内相互独立，在不同的数据集中，相对应的源信号有着最大的相关性，即IVA同时使用了二阶统计量(SOS)与高阶统计量(HOS)。与ICA，CCA方法相比，IVA模型对肌电噪声的分离更加有效，减少了去噪过程中脑电信号信息的丢失。将IVA与状态估计相结合，提高了复杂环境下噪声去除的准确性，满足了实际应用中脑电信号预处理的需求。

附图说明

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

图2a为真实脑电信号19个通道10秒的示意图；

图2b为真实肌电信号10秒的示意图；

图2c为观测信号在SNR＝1.5，19通道10秒的示意图；

图3为动态盲分离得到的独立源信号和转移路径在SNR＝1.5，19通道10秒的示意图；

图4为本发明方法与其他方法消噪后的脑电信号与原始观测信号在SNR＝1.5，T7通道10秒的示意图；

图5为本发明方法与CCA，ICA，IVA和HMM-ICA方法在RRMSE与ACC去噪性能上的比较图。

具体实施方式

本实施例中，如图1所示，一种状态相关的动态脑电信号中肌电伪迹的消除方法：首先将收集到的脑电观测信号通过延时得到多个数据集，再用隐马尔科夫独立向量分析的方法对多个数据集进行联合盲源分离；得到观测信号的源信号矩阵和逆混合矩阵；依据自相关系数将部分噪声相关的独立源信号置零，保留脑电独立信号源；最后盲源分离逆变换重建后得到消除噪声后的脑电信号。

下面通过半仿真脑电信号为例，结合附图来说明具体的实施方式。

通过真实脑电和肌电信号数据，构造了半仿真数据集，详细介绍本发明的实施方式，并与ICA，CCA，IVA，HMM-ICA处理方法相比较。

步骤一：由脑电测量设备采集并记录N通道的纯净脑电信号矩阵，记为：S_EEG＝[s₁,s₂,...,s_n,...,s_N]∈R^N×T，1≤n≤N，s_n为第n通道的纯净脑电信号，T表示总的采集时长；

采集收缩和舒张拳头时双手前臂的肌电源信号矩阵记为S_EMG∈R^N×T；

利用式(1)构造脑电观测信号矩阵X_EEG＝[x₁,x₂,...,x_n,...,x_N]，x_n为脑电观测信号矩阵X_EEG中第n通道观测信号，1≤n≤N：

X_EEG＝S_EEG+λ·X_EMG (1)

式(1)中，λ表示控制肌电噪声干扰脑电信号的程度；X_EMG∈R^N×T表示肌电观测信号，并有：

X_EMG＝[V₁S_EMG,V₂S_EMG,V₁S_EMG] (2)

式(2)中，V₁∈R^N×N，V₂∈R^N×N表示两个随机混合矩阵；

当时刻t∈[0,t₁]∪[t₂,T]时，肌电观测信号X_EMG由第一随机混合矩阵V₁作用于肌电源信号矩阵S_EMG产生，并表示为第一个混合模式；

当时刻t∈[t₁,t₂]时，肌电观测信号X_EMG由第二随机混合矩阵V₂作用于肌电源信号矩阵S_EMG产生，并表示为第二个混合模式；

利用式(3)得到信噪比SNR：

式(3)中，RMS(·)表示均方根函数；

本实施例中，图2a表示由脑电测量设备采集并记录时间长度T＝10s，N＝19通道的干净脑电信号矩阵。为了获得多样的肌电信号，二十三名志愿者按照要求做三种类型的肌肉活动分别是(1)持续攥拳10s；(2)放松5s，再持续攥拳5s；(3)瞬间攥拳时间小于1s，再放松，数据记为S_EMG，如图2b所示。其数据大小同脑电信号相匹配。接下来，在不同时间段对脑电信号和肌电信号作用以不同的混合模式。具体的，设置两个混合矩阵V₁，V₂。其中，V₁作用在t∈[0,2.5]，以及t∈[7.5,10]，V₂作用在t∈[2.5,7.5]，根据式(1)构造观测信号X_EMG＝[V₁S_EMG,V₂S_EMG,V₁S_EMG]，X_EEG＝S_EEG+λ×X_EMG，如图2c所示。

步骤二：将脑电观测信号矩阵X_EEG通过延时K-1个点构造K个数据集X＝[X₁,X₂,...,X_k,...,X_K]，1≤k≤K，其中，X_k表示第k个数据集在时刻t的脑电观测信号；

步骤三：利用隐马尔科夫独立向量分析法对K个数据集X进行动态联合盲源分离，得到不同混合模式下第k个数据集X_k的脑电源信号矩阵S_k以及解混矩阵W_k：

步骤3.1、令隐马尔科夫模型的转移状态记为Q＝{q₁,q₂,...,q_m,...,q_M}，q_m为第m个转移状态，1≤m≤M，M为转移状态总数；

K个数据集X实施维特比算法，得到转移状态序列

q′_t为时刻t所处的转移状态；

利用所述转移状态序列

标识出K个数据集X混合模式的变化；

步骤3.2、在每一个混合模式中，对K个数据集X进行联合盲源分离，得到K个解混矩阵W＝[W₁,W₂,...,W_k,...W_K]，W_k为在当前混合模式下第k个数据集X_k的解混矩阵；

利用式(5)得到在当前混合模式下第k个数据集X_k的脑电源信号矩阵S_k：

S_k＝W_kX_k，k＝1,2,...,K (5)

本实施例中，独立向量分析模型(IVA)如下：

S_k＝W_kX_k,k＝1....K, (6)

式(6)中，W_k表示第k个数据集的解混合矩阵，S_k为分解出的源信号矩阵，记B＝{W_k},k＝1,2,...,K。模型保证S_k中各源相互独立，并且，数据集之间对应源保持相关性。根据多元分布的概率密度函数可以得到分解源信号与观测信号间的关系，即：

式(7)中，

表示第k个分解源信号矩阵的第i个源，

表示第k个观测数据集的第i个通道。sⁱ表示模型的第i个SCV。本发明对于动态环境的识别是由隐马尔科夫模型(HMM)来完成。具体的，它的状态转移部分设置如下，初始概率分布π，状态转移矩阵A，其中，转移矩阵尺寸设置为M*M,用以控制模型在M个状态之间转换。如前，将隐马尔科夫的观测部分设置为IVA模型，对应于模型的每一个隐藏状态，观测信号产生过程服从IVA模型。为了求解本发明方法所设置的参数，首先给出HMM模型的对数似然函数：

式(8)中，T表示观测数据长度，X(t)为在时刻t脑电观测信号的取值，记X(t)＝[X₁(t),X₂(t),...,X_K(t)]。

为整个模型的观测序列，取值是q′_i＝M,i＝1,2,...,T。为了求解模型设置的参数，需要优化辅助函数

其中记θ＝{π,A,B}。根据前向-后向算法可以求解模型的初始分布π，状态转移矩阵A。关于解混合矩阵求解过程如下，根据式(6)，隐藏状态q′_t下观测数据的多元概率密度函数可以表示为：

对上式解混矩阵求导，得到解混矩阵的更新公式：

式(10)中，

为在转移状态q_m下第k个数据集的解混矩阵，μ为学习率用于控制更新速率。γ_m[t]表示观测信号在时刻t处于隐藏状态q'_m的概率。Φ^k(s(t))＝[Φ^k(s¹(t)),Φ^k(s²(t)),...,Φ^k(s^N(t))]，且

利用隐马尔科夫独立向量分析法对2个数据集X进行动态联合盲源分离，经由维特比算法对观测信号划分状态，得到脑电观测信号矩阵X的脑电源信号矩阵S_k，k＝1,2。由于观测部分的多元高斯分布独立向量分析法考虑了数据集间相对应的源信号的二阶相关性，所以其分离出的源信号按照自相关系数从大到小自动排列。图3为经由动态盲源分离处理得到的独立源信号。图中第一行标示了隐藏状态的转移路径，不难发现，模型清晰的给出了对不同混合模式的鉴别。并在给出的分段里，完成动态联合盲源分离。

步骤四：对应于每一个混合模式，将第k个脑电源信号矩阵S_k中自相关系数较小的通道信号置零，得到肌电噪声分量置零后的脑电源信号矩阵

本实施例中，根据肌电噪声自相关系数较低的物理特性，可以将排在末尾的独立信号源置零，从而达到消除肌电噪声的目的。

步骤五：利用式(4)得到第k个消噪后的脑电信号

式(4)中，(W_k)^-1为第k个逆混合矩阵W_k的逆矩阵。

图4为经由盲源分离逆变换处理，重建得到的消噪脑电信号与其它方法的对比结果。灰色背景为原始脑电观测信号，黑色细线为重建脑电信号。可以直观的看出，本发明所提出的算法，在保留脑电信号的同时，对肌电噪声的去除更为彻底。

为了量化评估本发明的效果，为此将本发明方法与ICA，HMM-ICA，IVA和CCA四种算法相比较。选择相对均方根误差(RRMSE)和平均相关系数(ACC)两个性能指标作为评价指标。相对均方根误差的定义如下：

ACC用以计算干净脑电与去噪后脑电各通道自相关系数的平均值。RRMSE的值越小越好，越小表示去噪后的信号和原始的干净脑电信号差别越小；ACC的值越大越好，越大表示的去噪后的信号越接近干净的脑电信号，噪声去除的越干净。

按照上述步骤，重复仿真实验100次，画出了各种方法的RRMSE，ACC对比图，如图5所示，从图5中可以看出，在各个信噪比情形下，本发明方法显著强于其他比较方法，表现出了动态环境下优秀的消噪能力。

综上所述，本发明设计了一个新的算法，隐马尔科夫独立向量分析算法，能够解决肌电噪声的消除问题，而且可以尽量保留脑电信息。该方法适用于脑电信号的预处理，并且为后续的诊断和研究提供方便。与其他方法相比，在动态环境下能够取得更好的去噪效果，解决了实际健康监护中，复杂脑电信号的动态噪声去除的问题，对进一步发展长期健康监护具有重要意义。

Claims (2)

1.一种状态相关的动态脑电信号中肌电伪迹消除方法，其特征是按如下步骤进行：

步骤一：由脑电测量设备采集并记录N通道的纯净脑电信号矩阵，记为：S_EEG＝[s₁,s₂,...,s_n,...,s_N]∈R^N×T，1≤n≤N，s_n为第n通道的纯净脑电信号，T表示总的采集时长；

采集收缩和舒张拳头时双手前臂的肌电源信号矩阵记为S_EMG∈R^N×T；

利用式(1)构造脑电观测信号矩阵X_EEG＝[x₁,x₂,...,x_n,...,x_N]，x_n为脑电观测信号矩阵X_EEG中第n通道观测信号，1≤n≤N：

X_EEG＝S_EEG+λ·X_EMG (1)

式(1)中，λ表示控制肌电噪声干扰脑电信号的程度；X_EMG∈R^N×T表示肌电观测信号，并有：

X_EMG＝[V₁S_EMG,V₂S_EMG,V₁S_EMG] (2)

式(2)中，V₁∈R^N×N，V₂∈R^N×N表示两个随机混合矩阵；

当时刻t∈[0,t₁]∪[t₂,T]时，肌电观测信号X_EMG由第一随机混合矩阵V₁作用于肌电源信号矩阵S_EMG产生，并表示为第一个混合模式；

当时刻t∈[t₁,t₂]时，肌电观测信号X_EMG由第二随机混合矩阵V₂作用于肌电源信号矩阵S_EMG产生，并表示为第二个混合模式；

利用式(3)得到信噪比SNR：

式(3)中，RMS(·)表示均方根函数；

步骤二：将所述脑电观测信号矩阵X_EEG通过延时K-1个点构造K个数据集X＝[X₁,X₂,...,X_k,...,X_K]，1≤k≤K，其中，X_k表示第k个数据集的脑电观测信号矩阵；

步骤三：利用隐马尔科夫独立向量分析法对K个数据集X进行动态联合盲源分离，得到不同混合模式下第k个数据集X_k的脑电源信号矩阵S_k以及解混矩阵W_k；

步骤四：对应于每一个混合模式，将第k个脑电源信号矩阵S_k中自相关系数较小的通道信号置零，得到肌电噪声分量置零后的脑电源信号矩阵

步骤五：利用式(4)得到第k个消噪后的脑电信号

式(4)中，(W_k)^-1为第k个逆混合矩阵W_k的逆矩阵。

2.根据权利要求1所述的状态相关的动态脑电信号中肌电伪迹消除方法，其特征是，所述步骤三是按如下过程进行：

步骤3.1、令隐马尔科夫模型的转移状态记为Q＝{q₁,q₂,...,q_m,...,q_M}，q_m为第m个转移状态，1≤m≤M，M为转移状态总数；

K个数据集X实施维特比算法，得到转移状态序列

q′_t为时刻t所处的转移状态；

利用所述转移状态序列

标识出K个数据集X混合模式的变化；

步骤3.2、在每一个混合模式中，对K个数据集X进行联合盲源分离，得到K个解混矩阵W＝[W₁,W₂,...,W_k,...W_K]，W_k为在当前混合模式下第k个数据集X_k的解混矩阵；

利用式(5)得到在当前混合模式下第k个数据集X_k的脑电源信号矩阵S_k：

S_k＝W_kX_k，k＝1,2,...,K (5)。

Images