基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法

1.本发明涉及单通道脑电信号伪迹去除技术领域，具体是一种基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法。

背景技术：

2.脑电信号(electroencephalogram，eeg)是大脑活动引起的电信号行为，可以反映脑神经细胞的电生理活动，可以获取到大量生理或心理活动以及一些疾病信息，是检测大脑功能和状况的有效手段之一。但在采集过程中，由于采集设备、外部环境及人体生理活动等影响，会不可避免的混入一些伪迹，如眼电伪迹、肌电伪迹等。这些伪迹会对后续脑电信号的处理造成很大影响，所以，我们需要开发一些有效的方法去除这些伪迹而且能够最大限度的保留有用信息。
3.目前来说，虽然没有一种特定的方法能够对所有类型的脑电信号都能有效去除伪迹，但是bss方法依旧是处理脑电信号最广泛、有效的方法。和其他方法一样，bss也有其局限性，此方法对源信号的通道数有要求，通道数较少会影响其效果，所以一般不会直接用于少通道或者单通道信号。而如今随着便携式脑电信号采集系统发展，很多情况下只能采集到少通道甚至单通道eeg，特别是生物医学方面，便携式脑电采集系统将更加适应社会需求。为了突破盲源分离无法直接应用于单通道eeg的限制，研究者们通常会引入信号分解方法，即先将原始eeg分解为多个信号分量，作为盲源分离的输入，再进行伪迹的识别与去除，最后重构出纯净的eeg。
4.在分解方法中，具有自适应分解的经验模态分解(emd)被广泛应用，但是用emd方法分解的信号会存在模态混叠现象，会造成伪迹去除不彻底或者误将有用信息去除的结果，而且，emd方法并没有坚实的理论基础，后续也有人提出了它的改进算法eemd、ceemd等方法，但还是会出现残存噪声等现象。
5.针对上述问题，本发明提出变分模态分解(variational mode decomposition，vmd)和sobi算法相结合的方法。既解决了emd分解方法的模态混叠现象，又克服了ica算法针对于处理不同伪迹的局限性。通过半模拟仿真实验表明，所提方法相较于其他现有去除方法具有更好的去除眼电、肌电伪迹效果，且能最大限度的保留有用信息。

技术实现要素：

6.本发明的目的是要解决单通道脑电信号常被多种伪迹污染而影响后续处理的问题，提供一种生物医学信号处理的途径。
7.为达到上述目的，本发明是按照以下技术方案实施：
8.一种基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法，包括以下步骤：
9.s1、采用变分模态分解(vmd)方法对被眼电伪迹(eog)和肌电伪迹(emg)同时污染的脑电混合信号y(t)进行分解，得到p个固有模态分量(imf)，表示为固有模态分量矩阵u(t)，u(t)＝[u1(t),u2(t),
…
,u
p
(t)]；
[0010]
s2、将p个imf作为盲源分离算法二阶盲辨识(sobi)的输入信号，进行进一步分离，得到p维信号分量矩阵s(t)，s(t)＝[s1(t),s2(t),
…
,s
p
(t)]。
[0011]
s3、计算经过s1和s2步骤后，得到的s(t)中的每个信号分量的模糊熵值，根据模糊熵值识别出伪迹分量，并将识别出伪迹分量置零，得到不含伪迹成分的信号矩阵
[0012]
s4、对信号分量矩阵s(t)中剩余的信号分量进行逆sobi变换，即对没有伪迹分量(没有置零)的信号分量进行sobi变换,得到不含伪迹成分的固有模态分量矩阵
[0013]
s5、对固有模态分量矩阵求和，重构出干净脑电信号x(t)。
[0014]
进一步的，所述s1中，具体分解步骤包括：
[0015]
s101、对被眼电伪迹(eog)和肌电伪迹(emg)同时污染的脑电混合信号y(t)建立变分约束模型：
[0016][0017]
其中，k为模态分解个数，k为总的模态分解个数，，uk(t)为第k个模态分量，ωk为uk(t)的中心频率,是每个模态分量信号相应的解析信号，是将乘子引入各个模态单边频谱中用以调整与之对应的中心频率ωk，将各模态分量的频谱调制到与其相对应的基频带中，使用l2范数对解调后信号进行高斯平滑，估计各个模态分量的带宽，并约束各模态分量的带宽之和最小，最后构造出带约束的变分模型，是模型约束，表示分解的各模态分量相加得到原始信号y(t)。
[0018]
s102、在约束条件下重构出{uk}和{ωk}，在式(1)引入二次惩罚函数项α和拉格朗日乘法算子λ将上述约束模型转化为非约束模型，即得到增广拉格朗日函数,表达式为：
[0019][0020]
s103、利用交替方向乘法器算法求解上述非约束模型得到拉格朗日函数的鞍点，从而得到窄带imf分量。具体过程如下：
[0021]
1)初始化参数
[0022]
2)n＝n+1，根据下式更新uk[0023][0024]
3)k＝k+1，重复步骤2，直到k＝k，终止上式，根据下式更新ωk[0025][0026]
4)k＝k+1，重复步骤3，直到k＝k，终止上式。
[0027]
5)更新拉格朗日乘数λ
[0028][0029]
式中，γ为噪声容忍度，满足信号分解的保真度要求。分别对应y(t)，λ(t)的傅里叶变换。
[0030]
6)重复步骤(2)-(5)，直到满足下式收敛条件
[0031][0032]
对于给定判断精度ε》0，结束循环，输出各imf的中心频率ωk。
[0033]
进一步的，所述s2中，具体分离步骤为：
[0034]
s201、首先对上述vmd分解得到的p维信号u(t)进行预白化处理，计算白化矩阵q，目的是去除通道间的相关性，提高分解效果。白化后的信号记为z(t)。
[0035]
z(t)＝qu(t)
ꢀꢀꢀ
(7)
[0036]
s202、计算z(t)的多个时延τ∈{τ1,τ2,
……
τk}的采样协方差矩阵：
[0037]
r(τ)＝e[z(t+τ)z(t)
t
]＝arz(τ)a
t
ꢀꢀꢀ
(8)
[0038]
s203、对上式计算的各个协方差矩阵r(τj)，进行联合近似对角化，计算正交矩阵v:
[0039]vt
r(τj)v＝djꢀꢀꢀ
(9)
[0040]
式中，{dj}是一组对角矩阵。
[0041]
s204、估计混合矩阵a和源信号s(t)：
[0042]
a＝q-1vꢀꢀꢀ
(10)
[0043]
s(t)＝a-1
u(t)＝v
t
qu(t)
ꢀꢀꢀ
(11)
[0044]
进一步的，所述s3中，具体计算步骤包括：
[0045]
s301、对上述sobi盲源分离算法分离得到的p维信号s(t)重构生成一组m维矢量为：
[0046][0047]
其中，i＝1,2,
…
,p-m，
[0048]
s302、定义为任意两个m维矢量之间的最大欧式距离，并用模糊函数μ定义相似度
[0049][0050]
式中，μ为关于的模糊隶属度函数，n为其边界的梯度，r为相似度容限。
[0051]
s303、定义：
[0052][0053]
s304、使维数增加1，对m+1为矢量重复步骤s301-s303，计算出m+1维度下的φ
m+1
(r)。
[0054]
s305、最后得到模糊熵的表达式：
[0055]
fe(m,r,p)＝lnφm(n,r)-lnφ
m+1
(n,r)
ꢀꢀꢀ
(15)
[0056]
式中，p为原始时间序列的长度。
[0057]
s306、通过上述步骤得到的每个信号分量的模糊熵值，利用模糊熵值识别出伪迹分量，并去除，得到不含伪迹成分的信号矩阵
[0058]
进一步的，所述s4中，具体包括：
[0059]
对上述得到的不含伪迹成分的信号矩阵按下面公式进行逆sobi变换：
[0060][0061]
式中,a为混合矩阵。
[0062]
进一步的，所述s5中，具体包括：
[0063]
根据以下求和公式重构出干净脑电信号：
[0064][0065]
式中，x(t)即为所求脑电信号。
[0066]
本发明的有益效果在于：脑电信号易被多种伪迹同时污染而影响后续处理，当前常见方法多是应对单种伪迹污染，且还存在一些问题，比如分解信号出现模态混叠现象、ica等盲源分离算法处理肌电伪迹效果不佳等。本发明提出一种基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法，能够解决经验模态分解出现的模态混叠现象，且对处理肌电伪迹的效果也较好，通过对仿真信号的实验结果表明，本发明的方法可以比较有效的去除脑电信号中的伪迹成分，并且可以最大限度的保留有用成分。
附图说明
[0067]
图1为vmd-sobi方法的流程图。
[0068]
图2为合成的原始脑电信号。
[0069]
图3为原始脑电信号经vmd-sobi方法后得到的信号分量图。
[0070]
图4为经vmd-sobi方法伪迹去除后的脑电信号。
[0071]
图5为经eemd-sobi方法伪迹去除后的脑电信号。
[0072]
图6为不同原始信噪比下的统计eemd-sobi和vmd-sobi两种方法的相对均方根误差。
[0073]
图7为不同原始信噪比下的统计eemd-sobi和vmd-sobi两种方法的相关系数。
具体实施方式
[0074]
为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
[0075]
本具体实施方式提出了一种基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法。
[0076]
总体流程如图1所示，具体步骤如下：
[0077]
步骤1：采用变分模态分解(vmd)方法对被眼电伪迹(eog)和肌电伪迹(emg)同时污染的脑电混合信号y(t)进行分解，得到多个固有模态分量(imf)，具体包括：
[0078]
1.1采用hilbert变换计算每个模态分量的解析信号，并获取其对应的单边频谱；
[0079]
1.2构建约束变分模型：
[0080][0081]
其中，k为模态分解个数，uk(t)为第k个模态分量，ωk为uk(t)的中心频率。
[0082]
1.3、引入二次惩罚函数项α和拉格朗日乘法算子λ将上述约束模型转化为非约束模型，表达式为：
[0083][0084][0085]
1.4、利用交替方向乘法器算法求解上述非约束模型得到拉格朗日函数的鞍点，从而得到窄带imf分量。具体过程如下：
[0086]
1)初始化参数
[0087]
2)n＝n+1，根据下式更新uk[0088][0089]
3)k＝k+1，重复步骤2，直到k＝k，终止上式，根据下式更新ωk[0090][0091]
4)k＝k+1，重复步骤3，直到k＝k，终止上式。
[0092]
5)更新拉格朗日乘数λ
[0093]
[0094]
式中，γ为噪声容忍度，满足信号分解的保真度要求。分别对应y(t)，λ(t)的傅里叶变换。
[0095]
6)重复步骤(2)-(5)，直到满足下式收敛条件
[0096][0097]
对于给定判断精度ε》0，结束循环，输出各imf的中心频率ωk。
[0098]
步骤2：分解得到的多个imf作为盲源分离算法二阶盲辨识(sobi)的输入信号，进行进一步分离，具体包括：
[0099]
2.1首先对上述vmd分解得到的p维信号u(t)进行预白化处理，计算白化矩阵q，目的是去除通道间的相关性，提高分解效果。白化后的信号记为z(t)。
[0100]
z(t)＝qu(t)
ꢀꢀꢀ
(7)
[0101]
2.2计算z(t)的多个时延τ∈{τ1,τ2,
……
τk}的采样协方差矩阵：
[0102]
r(τ)＝e[z(t+τ)z(t)
t
]＝arz(τ)a
t
ꢀꢀꢀ
(8)
[0103]
2.3对上式计算的各个协方差矩阵r(τj)，进行联合近似对角化，计算正交矩阵v:
[0104]vt
r(τj)v＝djꢀꢀꢀ
(9)
[0105]
式中，{dj}是一组对角矩阵。
[0106]
2.4估计混合矩阵a和源信号s(t)：
[0107]
a＝q-1vꢀꢀꢀ
(10)
[0108]
s(t)＝a-1
u(t)＝v
t
qu(t)
ꢀꢀꢀ
(11)
[0109]
如图2为合成的仿真信号，图3为经vmd-sobi方法分解分离后的信号分量。
[0110]
步骤3：将上述得到的多个信号分量分别计算模糊熵值，具体包括：
[0111]
3.1对上述sobi盲源分离算法分离得到的p维信号s(t)重构生成一组m维矢量为：
[0112][0113]
其中，i＝1,2,
…
,p-m，
[0114]
3.2定义为任意两个m维矢量之间的最大欧式距离，并用模糊函数μ定义相似度
[0115][0116]
式中，μ为关于的模糊隶属度函数，n为其边界的梯度，r为相似度容限。
[0117]
3.3定义：
[0118][0119]
3.4使维数增加1，对m+1为矢量重复步骤3.1-3.3，计算出m+1维度下的φ
m+1
(r)。
[0120]
3.5最后得到模糊熵的表达式：
[0121]
fe(m,r,p)＝lnφm(n,r)-lnφ
m+1
(n,r)
ꢀꢀꢀꢀ
(15)
[0122]
式中，p为原始时间序列的长度。
[0123]
如下表1为计算各个信号分量的模糊熵值。通过模糊熵值可以得出ic2和ic3为伪迹分量，也可以从图2可以看出ic2为眼电伪迹，ic3为肌电伪迹，所以使其置零。
[0124]
表1不同信号分量的模糊熵
[0125]
信号分量ic1ic2ic3模糊熵1.7
×
10-32×
10-4
8.84
×
10-2
[0126]
步骤4：对上述得到的不含伪迹成分的信号矩阵按公式(16)进行逆sobi变换：
[0127][0128]
式中,a为混合矩阵。
[0129]
步骤5：根据以下求和公式(17)重构出干净脑电信号：
[0130][0131]
式中，x(t)即为所求脑电信号。如图4为经本文所提方法伪迹去除后的脑电信号。
[0132]
对比实施例
[0133]
为说明上述实施例的方法的有效性和优越性，采用eemd-sobi方法对上述仿真数据进行伪迹去除，去除结果如图5所示。
[0134]
下表2是不同方法在原始信噪比-1db下的各自性能指标。
[0135]
表2
ꢀ‑
1db不同方法下的性能指标
[0136]
去噪方法ησγeemd-sobi5.25240.54610.8615本文方法7.51300.42110.9071
[0137]
图6和图7是四种方法在不同信噪比(-1.5db～1.5db)下的相对均方根误差和相对系数。
[0138]
图6和图7显示了eemd-sobi和vmd-sobi两种方法的性能对比，cc值越低，表明保留原始信号有效信号的能力越差，rrmse值越高表明去除伪迹的能力越差。可以看出，本文提出的vmd-sobi方法相较于已有方法eemd-sobi去除脑电信号中的eog和emg伪迹的效果更好，保留有用信息的能量也更强，尤其是在低信噪比的状态下。而且可以看到vmd的分解能力相较于eemd受信噪比的影响较小，证明了vmd-sobi处理eog和emg的优越性。

技术特征：

1.基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法，其特征在于，包括以下步骤：s1、采用变分模态分解vmd方法对被眼电伪迹eog和肌电伪迹emg同时污染的脑电混合信号y(t)进行分解，得到p个固有模态分量imf，将固有模态分量imf表示为固有模态分量矩阵u(t)，u(t)＝[u1(t)，u2(t)，...，u
p
(t)]；s2、将p个imf作为盲源分离算法二阶盲辨识sobi的输入信号，进行进一步分离，得到p维信号分量矩阵s(t)，s(t)＝[s1(t)，s2(t)，...，s
p
(t)]；s3、计算经过s1和s2步骤后，得到的s(t)中的每个信号分量的模糊熵值，根据模糊熵值识别出伪迹分量，并将识别出伪迹分量置零，得到不含伪迹成分的信号矩阵s4、对信号分量矩阵s(t)中剩余的信号分量，即没有置零的信号分量，进行逆sobi变换，得到不含伪迹成分的固有模态分量矩阵s5、对固有模态分量矩阵求和，重构出干净脑电信号x(t)。2.根据权利要求1所述的单通道脑电信号伪迹去除方法，其特征在于，所述步骤s1具体包括：s101、对被眼电伪迹(eog)和肌电伪迹(emg)同时污染的脑电混合信号y(t)建立变分约束模型：其中，k为模态分解个数，k为总的模态分解个数，u
k
(t)为第k个模态分量，ω
k
为u
k
(t)的中心频率，是每个模态分量信号相应的解析信号，是带约束的变分模型，是模型约束，表示分解的各模态分量相加得到原始信号y(t)；s102、在约束条件下重构出{u
k
}和{ω
k
}，在式(1)基础上引入二次惩罚因子α和拉格朗日乘子λ，将上述约束模型转化为非约束模型，表达式为：s103、利用交替方向乘法器算法求解上述非约束模型得到拉格朗日函数的鞍点，从而得到窄带imf分量，1)初始化参数2)n＝n+1，根据下式更新u
k
3)k＝k+1，重复步骤2，直到k＝k，终止上式，根据下式更新ω
k
4)k＝k+1，重复步骤3，直到k＝k，终止上式，5)更新拉格朗日乘数λ式中，γ为噪声容忍度，满足信号分解的保真度要求，n为迭代次数n＝{1，...，n}。分别对应y(t)，λ(t)的傅里叶变换，6)重复步骤(2)-(5)，直到满足下式收敛条件，循环迭代的收敛条件以各imf更新前后之差与更新前imf的比值作为系数：对于给定判断精度ε＞0，结束循环，输出各imf的中心频率ω
k
。3.根据权利要求2所述的单通道脑电信号伪迹去除方法，其特征在于，所述步骤s2具体包括：s201、首先对上述vmd分解得到的p维信号u(t)进行预白化处理，计算白化矩阵q，白化后的信号记为z(t)，z(t)＝qu(t)
ꢀꢀꢀꢀ
(7)s202、计算z(t)的多个时延τ∈{τ1，τ2，......τ
k
}的采样协方差矩阵r(τ)＝e[z(t+τ)z(t)
t
]＝ar
z
(τ)a
t
ꢀꢀꢀꢀꢀ
(8)s203、对上式计算的各个协方差矩阵r(τ
j
)，进行联合近似对角化，计算正交矩阵vv
t
r(τ
j
)v＝d
j
ꢀꢀꢀꢀ
(9)式中，{d
j
}是一组对角矩阵，s204、估计混合矩阵a和源信号s(t)：a＝q-1vꢀꢀꢀꢀꢀ
(10)s(t)＝a-1
u(t)＝v
t
qu(t)#(11) 。4.根据权利要求3所述的单通道脑电信号伪迹去除方法，其特征在于，所述步骤s3中计算信号模糊熵值的具体步骤包括：s301、对上述sobi盲源分离算法分离得到的p维信号s(t)重构生成一组m维矢量为其中，i＝1，2，...，p-m，{g(i)，g(i+1)，...，g(i+m-1)}代表从第i个点开始的连续m个g的值，g0(i)表示均值，s302、定义为任意两个m维矢量之间的最大欧式距离，即两者对应元素中差值最大的一个，并用模糊函数μ定义相似度值最大的一个，并用模糊函数μ定义相似度
式中，μ为关于的模糊隶属度函数，n为其边界的梯度，r为相似度容限，其中，(i，j＝1，...，p-m，i≠j)s303、定义函数s304、使维数增加1，对m+1为矢量重复步骤s301-s303，计算出m+1维度下的φ
m+1
(r)s305、得到模糊熵的表达式fe(m，r，p)＝lnφ
m
(n，r)-lnφ
m+1
(n，r)
ꢀꢀꢀꢀ
(15)式中，p为原始时间序列的长度；s306、通过上述步骤得到的每个信号分量的模糊熵值，利用模糊熵值识别出伪迹分量，并去除，得到不含伪迹成分的信号矩阵5.根据权利要求4所述的单通道脑电信号伪迹去除方法，其特征在于，所述步骤s4中包括：对上述得到的不含伪迹成分的信号矩阵按公式(16)进行逆sobi变换：式中，a为混合矩阵。6.根据权利要求5所述的单通道脑电信号伪迹去除方法，其特征在于，所述步骤s5中包括：根据以下求和公式(17)重构出干净脑电信号：。

技术总结

本发明涉及脑电信号伪迹去除技术领域，尤其是公开了一种基于变分模态分解和二阶盲辨识的脑电信号伪迹去除方法。本发明首先对被眼电伪迹（EOG）和肌电伪迹（EMG）两种伪迹同时污染的脑电信号进行VMD分解为多个固有模态函数（IMF），然后将这多个IMF作为SOBI的输入信号进一步盲源分离，使得脑电信号和伪迹信号分离开，接着计算每个分离得出的信号分量的模糊熵值，根据模糊熵值识别出伪迹分量，并将伪迹分量置零，最后做逆SOBI变换，重构出纯净脑电信号。通过对仿真信号的实验结果表明，本发明的方法可以比较有效的去除脑电信号中的伪迹成分，并且可以最大限度的保留有用成分。并且可以最大限度的保留有用成分。并且可以最大限度的保留有用成分。