CN109523486B - 噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法 - Google Patents

Abstract

本发明公开了一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法，其步骤包括：1、建立噪声环境下多通道脑电信号模型；2、基于多通道脑电信号的协稀疏性和低秩性，建立噪声环境下多通道脑电信号对应的优化模型；3、采用交替方向乘子法求解所述优化模型。本发明考虑了按照噪声分布的特点分为稠密噪声和稀疏噪声两种类型的噪声，能提高多通道脑电信号在压缩感知的过程中对噪声干扰的鲁棒性，从而能提高重构后的多通道脑电信号的精度和准确性。

Description

噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法

技术领域

本发明涉及信号处理领域，具体涉及一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法。

背景技术

脑电图(Electroencephalogram,EEG)信号是最常用的生物医学信号之一，它在医疗保健、脑-计算机接口等方面有着重要的应用。连续的脑电图监测通常需要采样和传输大量的数据，对设备的硬件要求很高。

实践中，脑电记录每天很容易产生1GB的数据，传输所需的能量非常高。传统的数据压缩方法是在传输前对数据进行压缩，由于大量的样本数据在压缩过程中被丢弃，严重浪费了资源。

为了应对这一挑战，提出了压缩感知技术(Compressed sensing，CS)，此时模拟信号在奈奎斯特(Nyquist)采样速率下不再是第一次采样，而是在压缩过程中被丢弃，直接以较低的采样率获得压缩信号，并通过非线性算法从压缩数据中恢复信号。

目前从压缩信号中重构多通道脑电信号的方法大多数基于l₀范数的优化问题，如基追踪法(basis pursuit)提出用l₁范数代替l₀范数来解决最优化问题，用线性规划方法来求解、正交匹配寻踪法(orthogonal matching pursuit，OMP)提出贪婪迭代的方式选择压缩测量矩阵的列，使得在每次迭代中所选择的列与当前的冗余向量最大程度地相关，从测量向量中减去相关部分并反复迭代，直到迭代次数达到稀疏度K，强制迭代停止。迭代硬阈值法(iterative hard thresholding，IHT)通过迭代执行硬阈值函数来解决最优化问题；2013年美国加州大学Z.Zhang课题组发现某些脑电信号在任何稀疏域都不是稀疏的，提出块稀疏贝叶斯学习法(block sparse Bayesian learning，BSBL)用于脑电信号的重构；2015年中国电子科技大学Yipeng Liu课题组提出，考虑到多通道脑电信号相互关联，多通道信号形成的矩阵具有低秩特性，建立了多通道脑电信号协稀疏低秩模型(simultaneouscosparsity and low-rank，SCLR)。

然而，现有的方法很少考虑噪声或仅考虑传输过程中产生的高斯噪声的影响。在实际情况下，噪声是不可避免的因素，根据噪声分布的特点，可以分成稠密噪声和稀疏噪声，当处理从复杂噪声环境中采集得到的压缩信号，以往的信号重构的方法性能都会下降。

发明内容

本发明针对现有技术的不足之处，提供一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法，以期能提高多通道脑电信号在压缩感知的过程中对噪声干扰的鲁棒性，从而能提高重构后的多通道脑电信号的精度和准确性。

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

本发明一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法的特点在于，包括以下步骤：

步骤1，利用式(1)建立噪声环境下多通道脑电模型：

Y＝ΦX+N+S (1)

式(1)中，X表示重构后的多通道脑电信号矩阵，且X∈R^m×n，n表示脑电信号的通道数，m表示每个通道脑电信号的数据长度；Φ表示压缩测量矩阵，且Φ∈R^k×m，k表示压缩后的数据长度，

表示信号的压缩率；N表示所述多通道脑电信号中的稠密噪声矩阵，并以高斯噪声信号矩阵为代表，且N∈R^k×n；S表示所述多通道脑电信号中的稀疏噪声矩阵，并以脉冲噪声信号矩阵为代表，且S∈R^k×n；Y表示压缩后的被噪声干扰的多通道脑电信号矩阵，且Y∈R^k×n；

步骤2，利用式(2)建立所述多通道脑电信号重构方法所对应的优化模型：

式(2)中，min表示最小化算子；Ω表示二阶差分矩阵生成的协稀疏分析字典，且Ω∈R^m×m；

表示重构后的多通道脑电信号矩阵X的核范数，其中σ_i表示重构后的多通道脑电信号矩阵X的第i个奇异值，且i＝1,...,r，r表示奇异值的数目；||N||_F表示稠密噪声矩阵N的弗罗贝尼乌斯范数，且

||S||₁表示稀疏噪声矩阵S中所有行和列元素绝对值的总和；λ作为正则化参数，用于平衡张量秩函数和l₀范数之间的贡献；α，β作为正则化参数，用于表示两种类型噪声在模型里的贡献作用；s.t.表示约束条件；

步骤3，采用交替方向乘子法求解所述优化模型，获得重构后的多通道脑电信号矩阵X：

步骤3.1，引入两个辅助变量V₁和V₂，对所述优化模型进行重写，得到如式(3)所示的重写后的优化模型：

利用式(4)得到所述重写后的优化模型所对应的增广拉格朗日函数：

式(4)中，

表示堆叠的辅助变量矩阵，

表示堆叠的约束条件矩阵，μ表示惩罚系数，

表示堆叠的拉格朗日乘子的缩放矩阵，分别表示式(3)中的三个约束条件对应的拉格朗日乘子的缩放矩阵；

步骤3.2，定义当前迭代次数为k，并初始化k＝0；初始化第k次迭代的变量，包括：X^k、N^k、S^k、

步骤3.3，利用式(5)更新第k+1次迭代的重构后的多通道脑电信号矩阵X^k+1：

式(5)中，I表示单位矩阵；

步骤3.4，利用式(6)更新第k+1次迭代的高斯噪声信号矩阵N^k+1：

步骤3.5，利用式(7)更新第k+1次迭代的脉冲噪声信号矩阵S^k+1：

式(7)中，令变量

则

是脉冲噪声信号软收缩算子，

表示脉冲噪声信号门限，sgn(x)表示x的符号函数，max(·)表示取较大值算子；

步骤3.6，利用式(8)更新第k+1次迭代的第一个辅助变量

式(8)中，令函数

表示当函数f(V₁)取最小值时所对应的V₁的取值，

是第一个辅助变量V₁的软收缩算子；

步骤3.7，利用式(9)更新第k+1次迭代的第二个辅助变量

式(9)中，

是奇异值收缩算子，令

令

Σ是重构后的多通道脑电信号矩阵

的奇异值构成的对角阵；U和W分别是正交列矩阵；

步骤3.8，利用式(10)更新第k+1次迭代的拉格朗日乘子的第一缩放矩阵

第二缩放矩阵

和第三缩放矩阵

步骤3.9，利用式(11)更新第k+1次迭代的原始误差u^k+1和第k+1次迭代的对偶误差d^k+1：

式(11)中，

表示对角单位矩阵，

表示负对角单位矩阵；

步骤3.10，判别收敛条件：

如果

且

则表示得到在噪声环境下重构后的多通道脑电信号矩阵X，其中ε表示收敛门限，否则，令k+1赋值给k，并回转执行步骤3.3。

与现有技术相比，本发明的有益效果在于：

1、本发明考虑了脑电信号实际采样传输过程中据噪声的分布特点划分了两类噪声，一类是不同信噪比的高斯噪声代表的稠密噪声，一类是不同噪声强度的脉冲噪声代表的稀疏噪声，即包含高斯噪声和脉冲噪声的复杂噪声环境，在此基础上建立多通道脑电模型，从而模拟出受到复杂噪声环境中多通道脑电信号的压缩感知，为解决复杂噪声下多通道脑电信号重构建立了数学模型。

2、本发明依据多通道脑电的协稀疏和低秩的结构特点，结合高斯噪声的概率分布和脉冲噪声的稀疏特性，建立了重构方法对应的优化模型，并采用交替方向乘子法求解该优化模型，从而适用于实际中的复杂噪声下多通道脑电信号重构问题，提高了重构信号的鲁棒性和精度。

附图说明

图1是本发明的流程图；

图2是本发明的重构信号均方根误差图；

图3是本发明的重构信号结构相似指数图。

具体实施方式

本实施例中，如图1所示，一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法主要由3个步骤组成：1、建立建立噪声环境下多通道脑电模型；2、基于多通道脑电信号的低秩性和稀疏性建立重构方法对应的优化模型；采用交替方向乘子法求解优化模型；具体的说，是按以下步骤进行：

步骤1，为了克服多通道脑电信号在复杂噪声环境下压缩感知受到噪声干扰的问题，本方法将多通道脑电压缩感知时的噪声环境划分为两大类，即稠密噪声和稀疏噪声。稠密噪声是指大部分的多通道脑电信号被噪声污染，主要包括高斯噪声。稀疏噪声是指一小部分的多通道脑电信号被噪声污染，主要包括脉冲噪声。

利用式(1)建立复杂噪声环境下多通道脑电信号模型

Y＝ΦX+N+S (1)

式(1)中，X表示干净多通道脑电信号，且X∈R^m×n，n表示脑电信号的通道数，m表示每个通道信号的数据长度；Φ表示压缩测量矩阵，且Φ∈R^k×m，k表示压缩后的数据长度，

表示信号的压缩率；N表示多通道脑电信号中的稠密噪声，且N∈R^k×n；S表示多通道脑电信号中的稀疏噪声，且S∈R^k×n；Y表示压缩后的被噪声干扰的多通道脑电信号，且Y∈R^k×n。

步骤2，研究表明，脑电信号很难被稀疏表示，并且传统稀疏矩阵Ψ要求列是非相关的；与传统的稀疏模型或块稀疏信号模型不同，协稀疏信号模型利用协稀疏分析字典Ω(Ω＝Ψ^-1)的相关性可以提高多通道脑电信号重构时的分辨率；根据脑电信号分段近似线性的特点，优化模型中将二阶差分矩阵作为协稀疏分析字典Ω，可以由l₀范数来表征，但基于l₀范数的优化问题通常是NP难问题，为此采用l₁范数来替代l₀范数；

同时，由于多通道脑电信号之间存在的相关性，促使在多通道脑电信号的压缩感知中可以利用数据的低秩结构，由于秩函数的优化问题也是NP难问题,因此优化模型中用核范数作为秩函数的凸代理函数；根据噪声的特点可以将噪声分成稠密噪声和稀疏噪声；稠密噪声是指污染大部分的多通道脑电信号的噪声，主要包括高斯噪声，根据高斯噪声的概率分布，故模型中可用F-范数来表征稠密噪声。稀疏噪声指污染小部分的多通道脑电信号的噪声，主要包括脉冲噪声，故稀疏噪声具有稀疏特性，类似地也可以用l₁范数来表征。于是在噪声环境下基于多通道脑电信号的协稀疏性和低秩性，利用式(2)建立多通道脑电信号重构方法对应的优化模型：

式(2)中，由于稀疏噪声S具有稀疏特性，可以由l₀范数来表征，但基于l₀范数的优化问题通常是NP难问题，为此采用l₁范数来替代l₀范数；min表示最小化算子；F-范数表征稠密噪声；Ω表示二阶差分矩阵生成的协稀疏分析字典，且Ω∈R^m×m；

表示X的核范数，其中σ_i(i＝1,...,r)表示矩阵的奇异值，r表示奇异值的数目；||N||_F表示稠密噪声矩阵N，

其中m，n表示稠密噪声矩阵N的行与列；||S||₁＝∑_i,j||S_i,j||表示稀疏噪声矩阵S中第i行第j列元素绝对值的总和；λ、α、β表示正则化参数，实施例中已调优；s.t.表示约束条件。

步骤3，由于交替方法乘子法广泛应用于求解约束条件下的优化问题，并取得了令人满意的效果，故本方法采用交替方法乘子法求解优化模型，获得重构多通道脑电信号X：

步骤3.1，引入两个辅助变量V₁和V₂,对优化模型进行重写，得到如式(3)所示的重写后的优化模型：

式(3)中，交替方向乘子法利用辅助变量将全局问题转换成多个容易求解的子问题，通过协调子问题的解来得到全局问题的解。

利用式(4)得到重写后的优化模型所对应的增广拉格朗日函数：

式(4)中，

表示堆叠的辅助变量矩阵，

表示堆叠的约束条件矩阵,μ表示惩罚系数，

表示堆叠的拉格朗日乘子的缩放矩阵，分别表示式(3)中的三个约束条件对应的拉格朗日乘子的缩放矩阵。

步骤3.2，定义当前迭代次数为k，并初始化k＝0；初始化第k次迭代的变量，包括：X^k、N^k、S^k、

步骤3.3，利用式(5)更新第k+1次迭代的重构的多通道脑电信号X^k+1：

步骤3.4，利用式(6)更新第k+1次迭代的高斯噪声信号N^k+1：

步骤3.5，利用式(7)更新第k+1次迭代的脉冲噪声信号S^k+1：

式(7)中，令

S_β/μ[x]＝sgn(x)max(|x|-β/μ,0)是软收缩算子，且β/μ表示门限，sgn(x)表示x的符号函数，max(·)表示取较大值算子；

步骤3.6，利用式(8)更新第k+1次迭代的第一个辅助变量

式(8)中，令函数

表示当函数f(V₁)取最小值时所对应的V₁的取值，

是第一个辅助变量V₁的软收缩算子；

步骤3.7，利用式(9)更新第k+1次迭代的第二个辅助变量

式(9)中，

是奇异值收缩算子，令

令

是重构后的多通道脑电信号矩阵

的奇异值构成的对角阵；U和W分别是正交列矩阵；

步骤3.8，利用式(10)更新第k+1次迭代的拉格朗日乘子的第一缩放矩阵

第二缩放矩阵

和第三缩放矩阵

步骤3.9，利用式(11)更新第k+1次迭代的原始误差u^k+1和第k+1次迭代的对偶误差d^k+1：

式(11)中，

表示对角单位矩阵，

负对角单位矩阵，I表示单位矩阵；

步骤3.10，判别收敛条件

如果

且

则表示得到在噪声环境下多通道脑电信号压缩感知重构的X，其中ε表示收敛门限，否则，令k+1赋值给k，并回转执行步骤3.3。此外，μ的选择对收敛速度有较大的影响，更新μ使得原始误差和对偶误差的范数的比值保持在一定的范围内，并且最终都收敛于0。

具体实施中，将进行半模拟实验来验证本发明提出算法(NRCS)的有效性，对比算法有同时贪婪分析追踪法(SGAP)，块稀疏贝叶斯学习法(BSBL)，同时稀疏低秩法(SCLR)并采用均方根误差(RMSE)和结构相似指数(MCC)来度量重构信号的精度。一般来说，RMSE越小，MCC越大，重构的精度越高。

将采用美国麻省理工学院发布的Physiobank数据库中的CHB-MIT scalp EEG数据集，在该数据集中随机抽选一组数据作为干净脑电数据，这组多通道脑电信号有23个通道，采样频率256Hz，选择100个片段进行实验，每段数据包含256×23个样本点，然后加入不同类型的噪声(1)高斯噪声：所有通道的信号都加入零均值的高斯白噪声，噪声的信噪比是30dB；(2)脉冲噪声，所有通道的信号加入脉冲噪声，噪声强度从1％到10％。

仿真结果分析：

图2和图3中横坐标表示脉冲噪声强度，纵坐标分别表示MSE和MCC。叉号曲线表示本发明随着噪声强度增大的性能指标变化曲线，加号曲线表示BSBL随着噪声强度增大的性能指标变化曲线，圆圈曲线表示SGAP随着噪声强度增大的性能指标变化曲线，星号曲线表示SCLR随着噪声强度增大的性能指标变化曲线。由图2可见，随着脉冲噪声强度增大，对比方法的均方根误差随之增大，而由图3可见，对比方法的信号结构相似指数随着脉冲噪声强度增强而下降，说明其他方法对于脉冲噪声变化的环境非常敏感，而本发明在脉冲噪声增强的噪声环境中重构的信号均方误差基本保持稳定，结构相似指数也保持较高的结果。

综上，由半模拟实验得到的结果表明，相较于对比方法，本发明提出的方法对多种复杂噪声条件下的多通道脑电信号重构更加鲁棒，并能提升重构的精度。

Claims (1)

1.一种噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法，其特征在于，包括以下步骤：

步骤1，利用式(1)建立噪声环境下多通道脑电模型：

Y＝ΦX+N+S (1)

式(1)中，X表示重构后的多通道脑电信号矩阵，且X∈R^m×n，n表示脑电信号的通道数，m表示每个通道脑电信号的数据长度；Φ表示压缩测量矩阵，且Φ∈R^k×m，k表示压缩后的数据长度，

表示信号的压缩率；N表示所述多通道脑电信号中的稠密噪声矩阵，且N∈R^k×n；S表示所述多通道脑电信号中的稀疏噪声矩阵，且S∈R^k×n；Y表示压缩后的被噪声干扰的多通道脑电信号矩阵，且Y∈R^k×n；

步骤2，利用式(2)建立所述多通道脑电信号重构方法所对应的优化模型：

式(2)中，min表示最小化算子；Ω表示二阶差分矩阵生成的协稀疏分析字典，且Ω∈R^{m ×m}；

表示重构后的多通道脑电信号矩阵X的核范数，其中σ_i表示重构后的多通道脑电信号矩阵X的第i个奇异值，且i＝1,...,r，r表示奇异值的数目；||N||_F表示稠密噪声矩阵N的弗罗贝尼乌斯范数，且

||S||₁表示稀疏噪声矩阵S中所有行和列元素绝对值的总和；λ作为正则化参数，用于平衡张量秩函数和l₀范数之间的贡献；α，β作为正则化参数，用于表示两种类型噪声在模型里的贡献作用；s.t.表示约束条件；

步骤3，采用交替方向乘子法求解所述优化模型，获得重构后的多通道脑电信号矩阵X：

步骤3.1，引入两个辅助变量V₁和V₂，对所述优化模型进行重写，得到如式(3)所示的重写后的优化模型：

利用式(4)得到所述重写后的优化模型所对应的增广拉格朗日函数：

式(4)中，

表示堆叠的辅助变量矩阵，

表示堆叠的约束条件矩阵，μ表示惩罚系数，

表示堆叠的拉格朗日乘子的缩放矩阵，分别表示式(3)中的三个约束条件对应的拉格朗日乘子的缩放矩阵；

步骤3.2，定义当前迭代次数为k，并初始化k＝0；初始化第k次迭代的变量，包括：X^k、N^k、S^k、V₁ ^k、

步骤3.3，利用式(5)更新第k+1次迭代的重构后的多通道脑电信号矩阵X^k+1：

式(5)中，I表示单位矩阵；

步骤3.4，利用式(6)更新第k+1次迭代的高斯噪声信号矩阵N^k+1：

步骤3.5，利用式(7)更新第k+1次迭代的脉冲噪声信号矩阵S^k+1：

式(7)中，令变量

则

是脉冲噪声信号软收缩算子，

表示脉冲噪声信号门限，sgn(x)表示x的符号函数，max(·)表示取较大值算子；

步骤3.6，利用式(8)更新第k+1次迭代的第一个辅助变量V₁ ^k+1：

式(8)中，令函数

表示当函数f(V₁)取最小值时所对应的V₁的取值，

是第一个辅助变量V₁的软收缩算子；

步骤3.7，利用式(9)更新第k+1次迭代的第二个辅助变量

式(9)中，

是奇异值收缩算子，令

令

Σ是重构后的多通道脑电信号矩阵

的奇异值构成的对角阵；U和W分别是正交列矩阵；

步骤3.8，利用式(10)更新第k+1次迭代的拉格朗日乘子的第一缩放矩阵

第二缩放矩阵

和第三缩放矩阵

步骤3.9，利用式(11)更新第k+1次迭代的原始误差u^k+1和第k+1次迭代的对偶误差d^k+1：

式(11)中，

表示对角单位矩阵，

表示负对角单位矩阵；

步骤3.10，判别收敛条件：

如果

且

则表示得到在噪声环境下重构后的多通道脑电信号矩阵X，其中ε表示收敛门限，否则，令k+1赋值给k，并回转执行步骤3.3。

Images