一种基于多元经验模式分解的肌电信号同步处理方法与流程

本发明属于多通道(或多变量)非平稳非线性信号同步分析的技术领域，涉及一种基于多元经验模式分解的多通道肌电信号同步处理方法。

技术背景

肌电信号(Electromyography，EMG)是肌肉收缩过程中人体神经肌肉系统的生理活动在机体皮肤表面产生的表征，是一种复杂的非线性非平稳生理电信号。由于SEMG是机体神经肌肉系统的直接表征形式，因此挖掘SEMG信号的生理机制可在一定程度上揭示机体神经肌肉系统的功能特性，在生理学上可为研究人体神经肌肉系统的生理反馈机制提供理论依据；在临床医用上也可为诊断与治疗神经肌肉系统功能性疾病(如运动神经元疾病和肌萎缩等肌肉疾病等)提供有效帮助；在康复工程上，包括人体运动康复、肌能评价等方面也能提供指导性意见。

现有的SEMG信号分析与处理方法主要是针对单一通道下信号进行分析，主要包括以下五种方法。(1)小波分析，该方法借助小波基函数，通过连续微分，使得被分解信号分量的时频特征有较好的独立性，在实时处理SEMG信号上具有优势，但由于小波分解具有交叉项，并不能完整的描述多成分、非平稳SEMG信号的自然属性，更不能合理解释运动神经元特征规律。(2)高阶统计量，该方法根据线性时不变模型对SEMG信号进行建模，能较好评估肌肉恒力收缩(静态)SEMG信号的非线性、非平稳性和非高斯性程度，但缺点是该方法不能对动态负荷或肌肉力作用下的SEMG特性进行分析。(3)人工神经网络，该方法虽然能较好地描述SEMG信号，但也存在一定的局限性，例如人工神经网络的训练过程往往需要大量的先验数据。(4)独立成分分析，该方法假设噪声信息和运动神经元信息属于不同的源信号，通过构造信号混合模型进行反解运算，从复杂信号中提取信号模式，能处理非高斯性的信号，但是该方法处理后得到的成分模式存在丢失，往往无法逆向重构原始信号。(5)标准经验模式分解和集合经验模式分解，上述两种方法都是基于希尔伯特-黄变换，不同于以传统的傅立叶变换为基础的线性和稳态频谱分析方法的是，这些方法是基于信号序列时间尺度局部特性的分解方法，能够直接的、后验的和自适应的处理SEMG信号。目前，标准经验模式分解或集合经验模式分解方法在非线性信号分析中已得到广泛认可。

然而，上述方法由于仅利用了SEMG信号在时间维度上的特性进行分析，并未考虑多通道表面肌电(Multi-channel Electromyography,MEMG)信号的空间信息，在处理MEMG信号时忽略了机体不同位置的肌肉群之间的信息。目前，多种生理参数监护与分析和复杂生物反馈系统的建模更能较为合理的揭示机体运动生理反应的机制问题，已逐步成为前沿研究热点。因此，传统的基于单一通道的SEMG信号分析由于缺乏在空间尺度上生理信息的分析，严重制约了上述技术的发展。针对这一问题，本专利公开了一种基于多元经验模式分解的多通道肌电信号同步处理方法，该方法将采集的MEMG信号通过低差异的哈默斯利序列将数据映射到多维向量空间中，进而在多变量条件下基于经验模式分解自适应的分解成若干个固有模态函数(Intrinsic Mode Functions,IMF)，有效的利用了所包含的MEMG时空间信息。本发明能确保MEMG分解后通道间IMFs阶数的一致性，同时也能实现分解后通道间对应各IMF信息具有较强的相关性。

技术实现要素：

本发明主要目的在于克服现有技术的不足，提供一种基于多元经验模式分解的多通道肌电信号同步处理方法，利用MEMG信号所包含的时空间信息，将采集的MEMG信号通过低差异的哈默斯利序列将数据映射到多维向量空间中，进而在多变量条件下基于经验模式分解自适应的分解成若干个固有模态分量，有效确保了每一通道IMFs阶数一致性，其次与标准经验模式分解或集合经验模式分解方法相比，多元经验模式分解(MEMD)得到的通道间IMFs信息具有较好的相关性。

为了实现以上目的，本发明的技术方案是一种基于多元经验模式分解的肌电信号同步处理方法，主要包括以下步骤：

步骤1：通过p,p≥3个通道肌电信号采集仪器同步拾取人体在特定运动动作条件下典型肌肉群处的肌电信号，获取肌电信号样本数据，构成一个p维原始信号s(t)，s(t)＝｛s₁(t),s₂(t),…,s_p(t)｝，其中t是时刻值；

步骤2：通过低差异的哈默斯利序列，将步骤1所获取的s(t)映射到一个p-1维球体上形成新的数据点集u_p＝{u₁,u₂,…,u_p}；

步骤3：将步骤2得到的点集，沿着方向向量集其中，p-1维球体上角度点集θ_v＝{θ_v1,θ_v2,…,θ_p-1|θ_v∈R^p}，方向向量v＝1,2,…,V分别进行投影并得到一组投影集

步骤4：计算步骤3得到的投影集的极值及对应的时刻

步骤5：将步骤4得到的极值矩阵插值，得到包络曲线集

步骤6：将步骤5得到的包络曲线集求平均，得到平均包络曲线，

步骤7：通过公式(5)提取成分d(t)：

$d\left(t\right)=s\left(t\right)-\stackrel{\‾}{a}\left(t\right)---\left(5\right)$

此时，若d(t)满足筛选条件，则将d(t)被定义为固有模态分量c(t)，同时将s(t)-d(t)视为s(t)，重复以上步骤1至7；若d(t)不满足筛选条件，则将d(t)视为s(t)，重复以上步骤1至7得到残余分量r(t)，停止循环计算；

通过以上步骤重复筛选后，原始信号s(t)将被分解为N个固有模态分量，{c_N(t)},N＝1,2,…,N和残余分量r(t)之和，如公式(6)所示：

$s\left(t\right)={\mathrm{\Σ}}_1^N{c}_N\left(t\right)+r\left(t\right)---\left(6\right)$

其中，筛选条件是其中a(t)是包络平均曲线的幅值，而σ是调节参数，其值根据实际情况决定且σ∈(0,1)。

进一步的，所述步骤2的获得低差异的哈默斯利序列的具体方法为：

首先将获得的多通道肌电信号数据构建Halton序列，即假设整个质数集由小到大排列后，前n个质数为{x₁,x₂,…,x_n}，则某一维度上的Halton序列的第i个样本为：

$r_i^x=\frac{a_0}{x}+\frac{a_1}{x^2}+\frac{a_2}{x^3}+...+\frac{a_s}{x^{s+1}}---\left(1\right)$

其中，基数x为整数，a_j(j＝1,2,…,s)是样本索引i按x进制展开后的j位上的数，s表示样本索引i按x进制展开后的总个数，因此样本索引i通过基数x的表示方法如下：

i＝a₀+a₁×x+a₂×x²+…+a_s×x^s (2)

将质数分别带入公式(2)，得到n个样本的Halton序列如下：

$(r_i^{x_1},r_i^{x_2},...,r_i^{x_n})---\left(3\right)$

由Halton序列得到n个样本的哈默斯利序列如下：

$(\frac{i}{n},r_i^{x_1},r_i^{x_2},...,r_i^{x_{n-1}})---\left(4\right).$

进一步的，所述步骤3的具体方法为：

将哈默斯利序列采样得到的(p-1)维球体上形成新的数据点集u_p＝{u₁,u₂,…,u_p}沿着(p-1)维球体的所有方向向量分别投影到球体上，得到一组投影集此时，投影集长度为单位长度，不具有物理意义，投影集方向向量值代表了u_p中各个元素。

进一步的，所述步骤4的方法为：

通过设定检索数据长度，并计算过零点个数，得出投影集对应的极值及时刻

进一步的，所述本发明肌电信号包括通过介入和非介入两种方式拾取的肌电电信号。

本发明一种基于多元经验模式分解的肌电信号同步处理方法，根据上述步骤各个通道分解出的固有模态分量包含了MEMG信号的时空间信息；基于多元经验模式分解的多通道肌电信号同步处理方法，能确保MEMG分解后通道间IMFs阶数的一致性，其次与标准经验模式分解或集合经验模式分解方法相比，该方法分解得到的通道间对应各IMF信息具有较强的相关性。

附图说明

图1为本发明的流程图；

图2下肢腿部实验协议；图2中(a)为坐立伸腿示意图，(b)为站立屈曲示意图，(c)为行走示意图；

图3表面肌电电极放置位置；图3中(a)为人体正面示意图，(b)为人体背面示意图；

图4第二类数据正常组被试1的坐立伸腿实验的MEMD结果中的股内侧肌部分，图4中(a)为第1-4阶IMF分量时域图，(b)为第5-8阶IMF分量时域图，(c)为第9-12阶IMF分量时域图，(d)为第13-16阶IMF分量时域图。

具体实施方式

下面结合附图详细说明，本发明在针对多通道肌电信号(MEMG)处理时(相比目前在肌电信号分析与处理领域中广泛采用的集合经验模式分解(EEMD)方法)具体实施方法以及性能分析结果。

典型肌肉群包括且不限于斜方肌、背阔肌、肩胛提肌、菱形肌、竖脊肌、胸大肌、腹外斜肌、腹直肌、胸小肌、肋间外肌、腹内斜肌、肩胛提肌、三角肌、肱二头肌、喙肱肌、三角肌、拇短伸肌、桡侧腕短伸肌、桡侧腕长伸肌、小指伸肌、指伸肌、肱桡肌、掌长肌、旋前圆肌、尺侧腕屈肌、桡侧腕屈肌、臀大肌、股二头肌、腓肠肌、趾长伸肌和比目鱼肌；

为验证本专利所提出的基于多元经验模式分解的多通道肌电信号同步处理方法的性能，首先实施了实验并采集了11个正常男性被试的下肢运动SEMG数据。在实验过程中，每位被试者被要求完成三种腿部运动动作，包括坐立状态下的腿部伸展(简称“坐立伸腿”)、站立状态下的腿部屈曲(简称“站立屈曲”)以及行走，如图2所示。为了消除个体差异的影响，每位被试至少执行4次实验动作。

多通道肌电信号采集设备的采样率为1000Hz，采样精度14位，单位mV。设备共有8个数字通道和4个模拟通道。实验过程中采用了其中5个通道，包括4个数字通道用于采集表面SEMG信号和1个模拟通道用于采集大腿与小腿间的夹角角度。原始数据被直接存储在设备的内置存储卡中，并通过蓝牙实时返回到数据采集系统的软件设备。第1至4通道的SEMG信号依次对应以下肌肉：股直肌(RF)、股二头肌(BF)、股内侧肌(VM)和半腱肌(ST)。肌电电极分布如图3所示。

将采集得到的四通道肌电信号数据，依次按照步骤2-7进行多元经验模式分解，分解后得到IMFs。以被试1的坐立伸腿实验为例，其股内侧肌处的肌电信号分解后得到16阶IMFs，具体结果如图4所示。此外，股直肌、股二头肌和半腱肌三个通道的EMG信号，也同时分解得到16阶IMFs。

为了验证本专利所公开方法在处理多通道肌电信号的性能优劣，我们采用了两个性能评估指标，即固有模态函数(IMF)阶数一致性以及模式队列程度，同时与目前在肌电信号分析与处理领域中广泛采用的集合经验模式分解(EEMD)方法相比较，对上述分解结果进行评估。固有模态函数(IMF)阶数一致性是比较各个通道分解出的IMFs阶数数量的差异性。而模式队列程度能够评估不同通道的同阶IMF分量之间的相关性程度，对于多通道EMG信号此相关性越高则每一通道对应IMF一致性越好，模式队列的程度越高。对于上述IMFs，其平均相关性系数的计算方法如下。首先，分别对各个IMF分量进行标准化并计算标准化IMF的功率谱；其次，进一步计算第n1(n1＝1,2,…,N)通道与第n2(n2＝1,2,…,N)通道的第m(m＝1,2,…,M)阶IMF分量功率谱之间的相关性系数c(n1,n2)，最终得到一个N阶相关性系数矩阵C_mC_m：

当n₁≠n₂n₁≠n₂时计算矩阵C_mC_m元素的平均值ρ_mρ_m，该值即为该被试对应实验的第m阶IMF的平均相关性系数。依次计算各阶IMF的平均相关性系数计算该组系数的平均值则可得到该被试在相应实验动作下IMFs的平均相关性系数ρ。

性能指标1：固有模态函数(IMF)阶数一致性

表1被试四通道SMEG信号经EEMD和MEMD处理所得IMFs的总阶数

由表1可知，对于不同通道的IMFs总阶数，EEMD不一定相等，而MEMD方法所得IMFs其不同通道的总阶数相等。

性能指标2：模式队列程度

根据分解的IMFs功率谱分析，计算坐立伸腿、站立屈曲、行走三种实验动作的平均相关性系数，用以衡量IMFs的模式队列性。为了进一步对比，这里对正常组所有被试的平均相关性系数进行统计。对于单个被试某个动作分解所得的一组IMFs，首先计算不同通道间同阶IMFs的平均相关性系数，再剔除80％能量对应频率小于20Hz的IMFs的平均相关性系数，计算剩余系数之间的平均值。最后求取正常组被试间相同实验动作在相同分解方法下对应系数的均值与标准差，结果如表2所示。

表2正常组EEMD和MEMD所得IMFs的平均相关性系数表

由表2可知，MEMD的平均相关性系数相近，对于三种实验动作中被试的EMG信号，MEMD的平均相关性系数外都在0.78以上，而三种实验动作中EEMD的平均相关性系数在0.69以下。从被试的角度分析，对于坐立伸腿实验中81.82％的被试、站立屈曲实验中90.91％的被试和行走实验中90.91％的被试，MEMD的平均相关性系数高于EEMD；对于坐立伸腿实验中90.91％的被试、站立屈曲中90.91％的被试和行走实验中63.64％的被试，MEMD的平均相关性系数高于EEMD。综上所述，除行走实验中MEMD的平均相关性系数相近，总体上MEMD的平均相关性系数高于EEMD。