CN111543990B - 肌肉协同的提取方法及基于肌肉协同的运动功能分析评估方法 - Google Patents

Abstract

本发明提供肌肉协同的提取方法及基于肌肉协同的运动功能分析评估方法。人体运动过程中采集多块肌肉的表面肌电信号，预处理后得到多通道肌电数据矩阵D，利用多元曲线分辨‑交替最小二乘对不同协同数下的方差占比进行分析，方差占比大于等于80％的最小协同数k为最终协同数，再利用多元曲线分辨‑交替最小二乘从肌电数据矩阵D中提取出k个肌肉协同。与现有基于非负矩阵分解的肌肉协同提取方法相比，该方法实现了肌肉协同的稳定和准确提取，尤其对于非稀疏肌肉协同提取具有更好的重复性和稳定性，为肌肉协同用于人体运动功能的分析评估奠定了坚实的基础。

Description

肌肉协同的提取方法及基于肌肉协同的运动功能分析评估 方法

技术领域

本发明属于生物信号处理技术与传感器技术领域，尤其涉及肌肉协同的提取方法及基于肌肉协同的运动功能分析评估方法。

背景技术

中枢神经系统如何控制身体的肌肉骨骼系统完成日常各种复杂的肢体动作一直是众多学者研究的焦点和热点。而肌肉协同理论能够较好地解释中枢神经系统对于肌肉骨骼系统的控制机制。

一个肌肉协同是由多块不同激活强度的肌肉组成，而一块肌肉可在多个肌肉协同中被激活。人体的各种复杂肢体动作能够通过少数几个肌肉协同在时域上的叠加控制完成。因此，肌肉协同可用于分析、评估人体的运动功能。

目前，肌肉协同主要通过矩阵分解算法分析多通道肌电数据计算得到。非负矩阵分解是应用最为广泛的肌肉协同提取方法，该方法首先随机产生初始协同和激活矩阵，并由重构矩阵与原始数据矩阵的残差构建目标函数，最后通过迭代优化的方式不断缩小重构矩阵与原始数据矩阵之间的残差，直到达到收敛条件停止迭代，实现肌肉协同的提取。非负矩阵分解具有原理简单、容易实现、分解结果非负性的优点。该方法的计算过程存在稀疏性约束，因此对于稀疏协同具有较好的分析效果。但是对于肌肉群普遍肌张力保持较高水平的个体而言，例如脑卒中、脊髓损伤等患者，其肌肉协同难以保证稀疏性。此外，人体很多多自由度运动的肌肉协同也是非稀疏的。在这些情况下非负矩阵分解法容易陷入局部最优，重复性较差，难以获得准确的肌肉协同。因此，稳定、可靠的肌肉协同提取方法具有重要的意义。

发明内容

针对上述技术现状，本发明人经过大量的试验探索后发现，采用多元曲线分辨-交替最小二乘法提取的肌肉协同具有唯一性，大大提高了肌肉协同提取的稳定性和准确性，对运动功能的分析评估具有重要意义。

即，本发明的技术方案为：

一种肌肉协同提取方法，包括如下步骤：

(1)人体运动过程中采集多块肌肉的表面肌电信号，预处理后得到多通道肌电数据矩阵D(m×n)，m为肌电数据的通道数，n为肌电数据的变量数；

(2)通过分析多通道肌电数据矩阵D(m×n)，提取肌肉协同；

其特征是：所述的步骤(2)中，采用多元曲线分辨-交替最小二乘法提取肌肉协同，即，针对肌电数据矩阵D(m×n)，利用多元曲线分辨-交替最小二乘法在不同协同数下的方差占比进行分析，方差占比大于80％的最小协同数k为最终协同数，2≤k≤m，再利用多元曲线分辨-最小二乘法从肌电数据矩阵D中提取k个肌肉协同。

2、如权利要求1所述的基于肌肉协同的运动功能分析评估方法，其特征是：所述步骤(2)包括如下：

(2-1)确定第1个纯变量

对于变量j，1≤j≤n，其平均值μ_j和均方差σ_j分别为：

其中，d_ij是肌电数据矩阵D中第i行、第j列的元素；

根据(式1)(式2)计算第1个纯度谱，第1个纯度谱中的第j(1≤j≤n)个元素为p_1,j：

其中，α是矫正参数，对信号噪声有一定的抑制作用，其值为1到5的常数；

在第1个纯度谱中，假设第t₁个元素纯度最大，则变量t₁(1≤t₁≤n)为第1个纯变量；

(2-2)确定第i个纯变量，2≤i≤k

首先利用变量长度λ对肌电数据矩阵D进行归一化，肌电变量j的长度λ_j为：

归一化后的肌电数据矩阵D(λ)中第i行、第j列元素为：

d(λ)_i,j＝d_i,j/λ_j (式5)

计算COO(Correlation around the origin)矩阵E：

E＝(1/m)D(λ)^TD(λ) (式6)

根据COO矩阵E，纯变量j的相关行列式为：

其中，t₁为第一个纯变量，ω_2,j作为权重系数用于计算第2个纯度谱，根据(式1)(式2)(式7)计算第2个纯度谱，第2个纯度谱中的第j个元素为p_2,j，1≤j≤n：

在第2个纯度谱中，假设第t₂个元素值最大，则变量t₂为第2个纯变量，1≤t₂≤n；

以此类推，对于i≥2时，相关行列式表示为：

对于i≥2时，第i个纯度谱中的第j个元素为p_i,j，2≤i≤k，1≤j≤n：

在第i个纯度谱中，假设第t_i个元素值最大，则变量t_i为第i个纯变量，1≤t_i≤n；

(2-3)假设协同数k＝2；

(2-4)从肌电数据矩阵D中提取出所有纯变量(t₁,t₂,…,t_k)对应位置的元素，组成初始协同矩阵S(m×k)，然后利用交替最小二乘(式11和式12)迭代优化，使目标函数(式13)值最小，

C＝(S^TS)^-1(S^TD) (式11)

S＝(DC^T)(CC^T)^-1 (式12)

再利用式14计算方差占比VAF(Variance account for)：

VAF＝1-(||D-SC||²/||D-mean(D)||²) (式14)

(2-5)当步骤(2-4)计算得到的方差占比小于80％时，协同数k增加1，重复步骤(2-1)、(2-2)、(2-4)，直到计算得到的方差占比大于或等于80％，对应的协同数k即为最终协同数，在该协同数下计算得到的矩阵S为最终的肌肉协同。

所述肌肉不限，包括上肢肌肉、下肢肌肉等。

本发明中的肌肉协同计算方法基于多元曲线分辨-交替最小二乘，包括纯变量计算、协同初始化矩阵、迭代优化等，首先计算各肌电信号变量的纯度谱，将成分贡献率较高的变量提取出来构建初始化协同矩阵，避免了传统非负矩阵分解法的随机初始化，然后基于初始化协同矩阵，利用交替最小二乘法优化分析结果，使多元曲线分辨的重构数据与原始数据更接近。与现有的非负矩阵分解法提取的肌肉协同相比，具有如下有益效果：

(1)对于非稀疏肌肉协同的提取具有更好的准确性和稳定性；

(2)无需重复运行程序寻找最小残差的分解结果，因此运行速度快于非负矩阵分解法；

(3)利用本发明提取得到的肌肉协同进行人体运动功能分析、评估、预测、控制时可提高相应的准确性与稳定性。

附图说明

图1是本发明实施例1中受试者A的数据矩阵D1经非负矩阵分解法重复计算20多次得到的肌肉协同统计盒须图。

图2是本发明实施例1中受试者A的数据矩阵D1经非负矩阵分解法以及多元曲线分辨-交替最小二乘法重复提取的肌肉协同统计图。

具体实施方式

下面结合实施例与附图对本发明作进一步详细描述，需要指出的是，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。

实施例1：

(1)六名受试者A、B、C、D、E、F(平均年龄28.5岁)进行规定的测试动作：六名受试者(右利手)右手握住测试杆进行折线运动，每名受试者重复该运动六次；运动过程中，利用Trigno Wireless System设备对7块肌肉：三角肌前束、三角肌后束、肱二头肌、肱三头肌、桡侧腕伸肌、肱肌和肱桡肌进行肌电数据采集，采样频率为2kHz。将采集到的7个通道肌电原始数据分别进行20Hz～450Hz的带通滤波、整流、1Hz的低通滤波(取包络线)、归一化等数据预处理，各受试者得到六次运动的数据矩阵D(7×1000)分别为D_n1、D_n2、D_n3、D_n4、D_n5、D_n6，其中n指代A、B、C、D、E或者F。

(2)利用非负矩阵分解法分析受试者A第一次运动得到的数据矩阵D_A1，提取肌肉协同。

(3)利用多元曲线分辨-交替最小二乘法分析该数据矩阵D_A1，提取肌中肉协同。具体如下：

(3-1)确定第1个纯变量

对于变量j(1≤j≤1000)，其平均值μ_j和均方差σ_j分别为：

其中，d_ij是肌电数据矩阵D_A1中第i行、第j列的元素；

根据(式1)(式2)计算第1个纯度谱，第1个纯度谱中的第j(1≤j≤1000)个元素为p_1,j：

其中，α是矫正参数，其值为5；

在第1个纯度谱中，假设第t₁个元素纯度最大，则变量t₁(1≤t₁≤1000)为第1个纯变量；

(3-2)确定第i(2≤i≤k)个纯变量，

首先利用变量长度λ对肌电数据矩阵D_A1进行归一化，肌电变量j的变量长度λ_j为：

归一化后的肌电数据矩阵D_A1(λ)中第i行、第j列元素为：

d(λ)_i,j＝d_i,j/λ_j(式5)

计算COO(Correlation around the origin)矩阵E：

E＝(1/m)D(λ)^TD(λ) (式6)

根据COO矩阵E，纯变量j的相关行列式为：

其中，t₁为第一个纯变量，ω_2,j作为权重系数用于计算第2个纯度谱，根据(式1)(式2)(式7)计算第2个纯度谱，第2个纯度谱中的第j(1≤j≤1000)个元素为p_2,j：

在第2个纯度谱中，假设第t₂个元素值最大，则变量t₂(1≤t₂≤1000)为第2个纯变量；

以此类推，对于i≥2时，相关行列式表示为：

对于i≥2时，第i(2≤i≤k)个纯度谱中的第j(1≤j≤1000)个元素为p_i,j：

在第i个纯度谱中，假设第t_i个元素值最大，则变量t_i(1≤t_i≤1000)为第i个纯变量；

(3-3)假设协同数k＝2；

(3-4)从肌电数据矩阵D_A1中提取出所有纯变量(t₁,t₂,…,t_k)对应位置的元素，组成初始协同矩阵S(m×k)，然后利用交替最小二乘(式11和式12)迭代优化，使目标函数(式13)值最小。收敛条件为:1)迭代次数小于1000次；2)目标函数(残差)迭代变化率小于0.01％。

C＝(S^TS)^-1(S^TD) (式11)

S＝(DC^T)(CC^T)^-1 (式12)

再利用式14计算方差占比VAF(Variance account for)：

VAF＝1-(||D-SC||²/||D-mean(D)||²) (式14)

(3-5)当步骤(3-4)计算得到的方差占比小于80％时，协同数k增加1，重复步骤(3-1)、(3-2)、(3-4)，直到计算得到的方差占比大于或等于80％，对应的协同数k即为最终协同数，在该协同数下计算得到的矩阵S为最终的肌肉协同；

经计算，本实施例中，最终协同数k＝3。

图1是采集得到的受试者A的同一数据矩阵D_A1经步骤(2)中的非负矩阵分解法20次重复计算得到的肌肉协同统计盒须图，可见，3个协同波动范围较大，平均相关系数分别为0.834、0.813、0.739，盒长分别为：0.183、0.180、0.288。

图2中左图是采集得到的受试者A的同一数据矩阵D_A1经步骤(2)中的非负矩阵分解法重复计算得到的肌肉协同结构统计图；图2中右图是采集得到的受试者A的同一数据矩阵D_A1经步骤(3)多元曲线分辨-交替最小二乘法重复计算得到的肌肉协同结构统计图。由图2可见，利用非负矩阵分解法得到的结果波动性较大，而利用多元曲线分辨交替最小二乘法得到的解具有唯一性，因此稳定性更好。

由于同一个人的肌肉协同是稳定的，因此利用同一受试者不同运动次数间的肌肉协同一致性来评价两种算法。即，进行如下步骤(4)与步骤(5)：

(4)利用非负矩阵分解法分别计算六名受试者的六次运动肌电数据矩阵D_n1、D_n2、D_n3、D_n4、D_n5、D_n6，其中n指代A、B、C、D、E或者F，提取肌肉协同，然后取平均值得到一致性相关系数。

(5)利用多元曲线分辨交替最小二乘法分别计算六名受试者的六次运动肌电数据矩阵D_n1、D_n2、D_n3、D_n4、D_n5、D_n6，其中n指代A、B、C、D、E或者F，提取肌肉协同，然后取平均值得到一致性相关系数。

在步骤(4)中，利用非负矩阵分解法计算得到的六名受试者六次运动的一致性相关系数为0.916、0.843、0.784；而在步骤(5)中，利用多元曲线分辨-交替最小二乘法计算得到的六名受试者六次运动的一致性相关系数为0.957，0.903，0.852，可见多元曲线分辨-交替最小二乘法在肌肉协同提取上稳定性更好。

因此，本实施例中，采用多元曲线分辨-交替最小二乘法分析数据矩阵D提取肌肉协同，然后根据肌肉协同分析、评估六名受试者的运动功能。

以上所述的实施例对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的具体实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改和改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.肌肉协同的提取方法，包括如下步骤：

(1)人体运动过程中采集多块肌肉的表面肌电信号，预处理后得到多通道肌电数据矩阵D(m×n)，m为肌电数据的通道数，n为肌电数据的变量数；

(2)通过分析多通道肌电数据矩阵D(m×n)，提取肌肉协同；

其特征是：所述的步骤(2)中，采用多元曲线分辨-交替最小二乘法提取肌肉协同，即，针对肌电数据矩阵D(m×n)，利用多元曲线分辨-交替最小二乘法对不同协同数下的方差占比进行分析，方差占比大于80％的最小协同数k为最终协同数，2≤k≤m，再利用多元曲线分辨-最小二乘法从肌电数据矩阵D中提取出k个肌肉协同；

所述步骤(2)包括如下：

(2-1)确定第1个纯变量

对于变量j，1≤j≤n，其平均值μ_j和均方差σ_j分别为：

其中，d_ij是肌电数据矩阵D中第i行、第j列的元素；

根据(式1)(式2)计算第1个纯度谱，第1个纯度谱中的第j(1≤j≤n)个元素为p_1,j：

其中，α是矫正参数，对信号噪声有一定的抑制作用，其值为1到5的常数；

在第1个纯度谱中，假设第t₁个元素纯度最大，则变量t₁(1≤t₁≤n)为第1个纯变量；

(2-2)确定第i个纯变量，2≤i≤k

首先利用变量长度λ对肌电数据矩阵D进行归一化，肌电变量j的长度λ_j为：

归一化后的肌电数据矩阵D(λ)中第i行、第j列元素为：

d(λ)_i,j＝d_i,j/λ_j (式5)

计算COO(Correlation around the origin)矩阵E：

E＝(1/m)D(λ)^TD(λ) (式6)

根据COO矩阵E，纯变量j的相关行列式为：

其中，t₁为第一个纯变量，ω_2,j作为权重系数用于计算第2个纯度谱，根据(式1)(式2)(式7)计算第2个纯度谱，第2个纯度谱中的第j个元素为p_2,j，1≤j≤n：

在第2个纯度谱中，假设第t₂个元素值最大，则变量t₂为第2个纯变量，1≤t₂≤n；

以此类推，对于i≥2时，相关行列式表示为：

对于i≥2时，第i个纯度谱中的第j个元素为p_i,j，2≤i≤k，1≤j≤n：

在第i个纯度谱中，假设第t_i个元素值最大，则变量t_i为第i个纯变量，1≤t_i≤n；

(2-3)假设协同数k＝2；

(2-4)从肌电数据矩阵D中提取出所有纯变量(t₁,t₂,…,t_k)对应位置的元素，组成初始协同矩阵S(m×k)，然后利用式11的激活矩阵C和式12交替最小二乘迭代优化，使式13的目标函数f(S,C)值最小，

C＝(S^TS)^-1(S^TD) (式11)

S＝(DC^T)(CC^T)^-1 (式12)

再利用式14计算方差占比VAF(Variance account for)：

VAF＝1-(||D-SC||²/||D-mean(D)||²) (式14)

(2-5)当步骤(2-4)计算得到的方差占比小于80％时，协同数k增加1，重复步骤(2-1)、(2-2)、(2-4)，直到计算得到的方差占比大于或等于80％，对应的协同数k即为最终协同数，在该协同数下计算得到的矩阵S为最终的肌肉协同。

2.如权利要求1所述的肌肉协同的提取方法，其特征是：所涉及的肌肉包括上肢肌肉与下肢肌肉。

3.基于肌肉协同的运动功能分析评估方法，其特征是：利用权利要求1或2所述的提取方法提取得到的肌肉协同分析评估个体的运动功能。

Images