[发明专利]基于多元变分模态分解的脑电多域特征提取方法

权利要求书

1.基于多元变分模态分解的脑电多域特征提取方法，其特征在于：包括以下步骤：

步骤(1)：采集受试者的多通道运动想象脑电信号，接着利用带通滤波器对脑电信号进行8-30Hz的低通滤波，用于脑电信号特征提取；

步骤(2)：对多通道数据进行MVMD分解；具体为：

对于包含C个数据通道的输入数据X(t)，记为[x₁(t),x₂(t),…x_C(t)]，MVMD算法的具体步骤如下：

①、预先定义k个窄带IMF分量u_k(t)，使

其中，u_k(t)＝[u₁(t),u₂(t),…u_C(t)]；

利用希尔伯特变换的算法分别得到向量u_k(t)中每一个元素的解析表示，记为以计算单边频谱，然后与指数项相乘来调整其对应的中心频率ω_k(t)，使每个模态的频谱被调制到与它相对应的基频带上，通过谐波转换后的梯度函数的L₂范数来估计出各模态u_k(t)的带宽；

由于单个的频率分量ω_k被用作为整个向量的谐波混频，因此需要在多元振荡u_k(t)中找到在多个通道的共同频率分量ω_k，考虑的所有通道，使每个通道所分解的IMF分量之和能够再现输入数据，且最小化模态函数带宽之和，受约束的相关优化问题变为：

其中，u_k,c(t)表示相应的通道数c和向量u_k(t)中的元素，表示对应的u_k,c(t)中每一个元素的解析表示，表示与时间相关的偏导数；

③、对变分问题求解，构造增广的拉格朗日表示为：

其中，α为二次惩罚因子，其保障存在高斯噪声时信号的重构精度，λ_c(t)为拉格朗日乘子，使得约束条件保持严格性；

④、为了解决这个转变后的非约束性的变分问题，应用乘子交替方向法来更新和以寻找增广拉格朗日表达式的鞍点；具体做法为：

1)模态更新

模态更新的最小化问题相当于：

利用Parseval傅里叶变换，将其变换到频域，与VMD算法的模式更新方法类似，得到模式更新为：

2)中心频率更新

由于增广拉格朗日函数的后两项不依赖于ω_k，则中心频率更新所对应的优化问题可以简化为：

同样利用Parseval傅里叶变换，将该优化转换到频域，得到在傅里叶域中的等效问题为：

得到中心频率的更新方法：

通过更新关系自适应地分解信号的频带，得到k个窄带IMF分量；此外，由于MVMD方法能同时计算多个通道数据，保证了通道间频率的一致性，使信号的分析更趋于稳定；

步骤(3)：将式(1)的各分量u_k(t)进行希尔伯特谱分析，定义为u_k(t)与1/πt的卷积，即：

其中，μ为积分变量；

由u_k(t)作为实部，它的希尔伯特变换作为虚部，得到其解析信号和瞬时幅值为：

U_k(t)＝u_k(t)+j₀H[u_k(t)]                  (10)

其中，j₀为虚数单位；

然后计算每个通道的平均瞬时能量值为：

其中n为采样点数目，为第i个采样点对应的瞬时能量值；取瞬时能量均值作为时域特征，记为F₁；

步骤(4)：

①根据采样点数对IMF分量作一维离散时间序列{u₁,u₂,…,u_i,…,u_n}，然后进行多尺度化处理来构造粗粒化形式的时间序列{y^(τ)}：

其中，τ为尺度因子，序列长度为M＝int(n/τ)；

②构造一组m维的向量：

X_m(i)＝{y_i+b:0≤b≤m-1}               (14)

③计算X(i)和X(j)间距离:

d[X(i),X(j)]＝max|y_(i+b)-y_(i+b)|             (15)

其中，0≤b≤m-1；i,j＝1,2...M-m+1；i≠j，M为序列长度；

④计算d[X(i),X(j)]中小于阈值r的个数，也是模板匹配的个数，记为B^m(i)，在计算B^m(i)与距离总数的比值，记作：

⑤计算的平均值C^m(r)：

⑥将维数增加至m+1，重复②～⑤步骤，计算和C^m+1(r)；

得到时间序列长度为M时的样本熵值：

SampEn(m,r,M)＝-ln[C^m+1(r)/C^m(r)]            (18)

重复以上步骤得到多尺度下的样本熵特征，记为F₂；

步骤(5)：将每个分量看作为单个通道的脑电信号，并以原始信号的采样频率对其采样，然后合并分量构造新的信号矩阵，矩阵维数为k×n，k表示IMF分量个数，n表示采样点数目；以C3,C4通道为例，矩阵表示为：

记X₁和X₂分别为按式(16)构造的两类运动想象任务下的信号矩阵，构造空间滤波器，首先求信号矩阵协方差，并归一化，得到：

其中，trace(·)为矩阵的迹，然后求出混合空间的协方差矩阵R,并进行特征值分解：

其中，为两类任务实验的平均协方差矩阵，U为矩阵R的特征向量，Λ为对应的特征值矩阵；进一步可求出白化矩阵P：

接着对R₁与R₂做白化处理以及主成分分解：

则矩阵S₁与矩阵S₂的特征向量相等，即B₁＝B₂，且Λ₁+Λ₂＝I，其中，I为单位矩阵；两类矩阵的特征值之和为1，因此当一类的特征值最大时另一类最小，从而可以实现分类的作用；则构造空间滤波器：

W＝B^TP                            (24)

将式(19)所示的信号矩阵经W滤波，得到特征矩阵，接着取其前q个和后q个特征值构成信号Z_P(P＝1,…,2q)，求取特征向量：

得到表示空域的特征，记为F₃；

步骤(6)：分别求出时域瞬时能量均值特征F₁，多尺度熵特征F₂，以及重构信号的CSP特征F₃，最后得到多域的结合的特征向量F＝{F₁,F₂,F₃}，并进行归一化；将训练集特征输入SVM分类器训练分类模型后，对测试集特征进行分类得到测试结果。