一种基于表面肌电分解的运动单位识别方法

摘要

本发明公开了一种基于表面肌电分解的运动单位识别方法，包括：1、预采集指定肌肉运动下的表面肌电信号，构建用以初始化的离线数据集；2、采用离线状态的逐步变量剥离算法分解离线数据集，解算每个运动单位对应的解混向量并保存；3、在线处理阶段，实时获取表面肌电信号并对其分窗，利用已获取的解混向量直接得到不同的运动单位源信号；4、利用自适应修正阈值的方法提取分析窗内源信号对应的运动单位发放序列；5、随着实时数据的输入，执行步骤3和4，待没有新的数据输入时，依次连接不同分析窗的分解结果并输出。本发明能以在线的方式识别表面肌电信号中的运动单位活动，并实时获取到运动单位的发放序列。

说明书

技术领域

本发明属于肌电信号处理技术领域，具体涉及一种基于肌电分解技术的运动单位识别方法，主要应用于高效智能的神经指令解码和运动意图捕获。

背景技术

表面肌电信号(Surface Electromyography,SEMG)作为一种伴随运动产生的电生理信号，反映驱动肌肉收缩运动的神经指令。其形成过程由多个运动单位(Motor Unit,MU)受中枢调控激活发放产生的运动单位动作电位(Motor UnitAction Potential,MUAP)沿肌纤维传播在检测电极处的时间、空间叠加而成的，直接反应了用户的运动意图和神经驱动信息。特别的，随着肌电采集技术的进步，包含有丰富的时空信息的高密度表面肌电为肌电控制应用注入了新的活力。更细致准确的解析肌电运动信息，特别是从刻画肌肉整体活动的宏观层面深入至探究支配肌肉活动的神经指令的微观层面，是肌电控制技术革新所面临的根本性问题。因此，现有的肌电控制技术不再满足于各项宏观时频域特征的简单应用，转而聚焦于蕴含在神经驱动下由大量的信源产生的丰富活动信息，寻找更符合运动底层生理规律的控制策略。在解析微观神经驱动信息的过程中，需要获取神经肌肉系统最基本的组成单位——运动单位的活动信息，因此肌电分解是不可或缺的一步。肌电分解可以看作肌电形成的逆过程，将混叠的表面肌电还原为其最基本的组成部分——MUAP序列的形式，帮助我们从运动单位的微观层面获取神经驱动信息。MUAP序列中包含了运动单位的发放和波形信息，其中发放信息是最为重要的，直接展示了神经指令下发和传达的时刻，表征为运动单位的0-1发放序列。大量的学者利用肌电分解算法识别出肌肉活动中运动单位的信息，做出了很多深入的研究，也驱动着肌电控制技术朝着解码微观神经指令的方向逐步发展。

在过去的十年内，最具代表性的肌电分解算法包括卷积核抵消算法(ConvolutionKernel Compensation,CKC)和逐步变量剥离算法(Progressive FastICAPeel-off,PFP)。CKC算法粗略地估计运动单位发放序列和表面肌电信号之间的互相关矢量，然后通过估计的互相关矢量和EMG信号的相关矩阵来计算该运动单位的发放模式。PFP算法通过盲源分离经典算法快速独立成分分析(Fast Independent Component Analysis,FastICA)算法求解出一系列运动单位源信号，提取出其中包括的发放序列，通过带时域约束的FastICA对结果进行检测和修正，并通过剥离策略将已得到的运动单位从原始信号中剥离出来，以寻找到更多的运动单位。两种肌电分解策略在离线的过程中取得了较好的成果，其分解得到的运动单位活动信息是解码表面肌电的神经指令的重要依据，直接反映了用户的运动意图和运动状态，提供了重要的人机接口，应用潜力巨大。

然而，目前的以肌电分解技术为基础的运动单位识别方法大都是在离线的过程中处理，鲜有在在线情况的研究，但是在科技飞速发展的现代社会，对信号的实时处理分析有着巨大的需求，离线状态的运动单位识别显然存在很严重的应用局限性。表面肌电信号呈现为动作电位波形的高度混叠，从中分解出代表微观神经指令的单个运动单位信息已然是领域内公认的难题。对于离线的表面肌电分解技术，需要对较长时间的信号段进行复杂的迭代计算求出解混向量，依赖专家检测的过程来提升分解结果的准确性，计算复杂性较高，确保性能的前提下难以满足肌电控制任务所需的实时获取分解结果的要求，这对于在线运动单位识别的是一个巨大的困难与挑战。

发明内容

本发明为了克服现有运动单位识别方法的不足之处，提出一种基于表面肌电分解的运动单位在线识别方法，以期能够在实时输入肌电信号的情形下，以在线的方式识别表面肌电信号中的运动单位活动，并实时获取到运动单位的发放序列，同时保持较高的精度，从而达到与离线分解相当的效果。

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

本发明一种基于表面肌电分解的运动单位识别方法的特点是按如下步骤进行：

步骤1、利用肌电测量设备和M通道高密度电极在待测肌肉处预采集t时刻肌肉收缩任务的表面肌电信号数据，记为x(t)＝[x

步骤2、采用逐步变量剥离算法分解所述离线数据集的肌电信号：

步骤2.1、利用式(1)对所述离线数据集进行扩展，得到扩展信号

式(1)中，K是延迟因子；

定义并初始化运动单位发放集合γ及其对应的解混向量集合ε为空集；

定义并初始化残差信号

步骤2.2、执行FastICA算法，设置算法执行的总次数为n：

步骤2.2.1、初始化i＝1；

步骤2.2.2、初始化j＝0；

步骤2.2.3、定义并初始化残差信号

步骤2.2.4、利用式(2)计算第i个解混向量在第j+1次的计算结果w

式(2)中，G'表示非多项式函数G的一阶导数；G”表示G的二阶导数；E是数学期望；

步骤2.2.5、如果

W'＝(WW

步骤2.2.6、计算运动单位源信号

步骤2.3：采用图像处理的大津法获取第i个运动单位源信号y

步骤2.3.1、将区间(0,max(y

步骤2.3.2、初始化h＝1；

步骤2.3.3、令th

步骤2.3.4、根据式(4)计算得到类间方差F

F

步骤2.3.5、若h＜M，则令h+1赋值给h，返回至步骤2.3.3；否则进入步骤2.3.6；

步骤2.3.6、找到最大的F

步骤2.4、将y

步骤2.5、通过带时域约束的FastICA算法检测分解结果的可靠性：

步骤2.5.1、定义迭代计算的次数为d，并初始化d＝1；定义并初始化第d次迭代的拉格朗日乘子μ

步骤2.5.2、通过式(5)得到第d+1次迭代的第i个解混向量w

式(5)中，g'

步骤2.5.3、对第d+1次迭代的第i个解混向量w

步骤2.5.4、根据式(6)、式(7)和式(8)分别更新其他参数：

μ

δ

ξ

式(6)、式(7)和式(8)中，μ

步骤2.5.5、更新第d+1次迭代的第i个运动单位源信号

步骤2.5.6、若

步骤2.5.7、令第i个运动单位源信号y

若c

步骤2.6、遍历运动单位发放集合γ并分别计算任意两个不同的发放序列之间的相关系数，若大于C，则表示两个发放序列重复，并随机选取一个从运动单位发放集合γ和解混向量集合ε中删除；否则，两个发放序列均保留；经此操作后，可得到共有q个发放序列的运动单位发放集合γ和共有q个解混向量的解混向量集合ε；

步骤2.7、将运动单位发放集合γ中q个发放序列表示为R＝[r

步骤2.8、将更新后的残差信号

步骤3、实时采集肌电信号，并将获取的肌电数据流按照时间轴分为一系列窗长为L，步长为D的分析窗；在每个分析窗内，根据式(1)得到扩展后的肌电数据X，并利用卷积模型

步骤4、采用自适应的阈值选取策略提取分析窗内第i个运动单位源信号s

步骤5、每当新的肌电数据输入时，按照步骤3和步骤4的过程进行处理和保存；

当没有新的数据输入时，利用式(9)将所保存的P个分析窗中的发放序列按照时间顺序进行拼接，得到总的分解结果Spike并作为最终的结果：

Spike＝[Win

式(9)中，Win

本发明所述的基于表面肌电分解的运动单位识别方法的特点也在于，所述步骤4中的自适应阈值选择策略是按如下过程进行：

步骤4.1、按照步骤2.3中的大津法计算第i个运动单位源信号s

步骤4.2、将s

步骤4.3、初始化k＝1；

步骤4.4、令中间变量th＝peak

步骤4.5、设置向量a，用来存储运动单位发放序列b

cos

式(10)中，std(a)表示向量a的标准差，mean(a)表示向量a的平均值；

步骤4.6、若k＜l，则令k+1赋值给k，返回步骤4.4；否则找到最大的cos

本发明能够实时识别表面肌电包含的运动单位，且能够以自适应的方式提取运动单位的发放序列，达到离线肌电分解相当的准确率，具体的有益效果体现在：

1、本发明步骤二对离线数据集进行了分解，获取到了不同运动单位对应的解混向量。基于运动单位结构与位置不变性，步骤三直接将解混向量集合作用于卷积模型求解新获取的肌电信号中的运动单位源信号，从而避免了庞大的迭代计算的步骤，将在线处理的部分计算工作转移到离线的过程中，为在线处理加速。

2、本发明步骤四中通过优化运动单位辨别与提取方法，采用自适应的阈值选取策略，克服了离线肌电分解需要专家检测的弊端，保证了所提的方法能够兼顾高性能和实时性。

3、本发明步骤三中通过设置随时间滑动的分析窗，保证了肌电数据实时输入和分解结果的实时输出，满足了在线处理数据的需求；步骤六通过相关系数的大小链接不同分析窗内的分解结果，保证了分解结果在时间上的连续性。

附图说明

图1为本发明实施例提供的方法流程示意图；

图2为本发明实施例提供的实验示意图；

图3为本发明实施例提供的在线处理和离线分解结果对比示意图；

图4为本发明实施例提供的在线处理结果示意图。

具体实施方式

本实施例中，一种基于表面肌电分解的运动单位识别方法，如图1所示，包括如下步骤：

步骤1、利用肌电测量设备和M通道高密度电极在待测肌肉处预采集t时刻肌肉收缩任务的表面肌电信号数据，记为x(t)＝[x

具体实施中包括：

(1)募集d位受试者，引导每位受试者将待测手放于3D打印模块上固定，拇指套入连接着力度传感器的戒指环，设备采集手掌大鱼际肌的M通道高密度表面肌电信号，如图2中的(a)部分所示。阵列电极阵列中单个电极触点的直径为p，电极中心间距为q。示例性的，可以设置：d＝8，M＝64，p＝2mm，q＝4mm。受试者记为S1-S8；

(2)逐一采集受试者执行拇指收缩任务的连续肌电信号。肌肉收缩任务设置为拇指肌力在2s的时间内从0逐渐增加至30％最大收缩力(在每次实验前测量不同受试者的拇指最大收缩力)，然后维持这个力度约3s，最后直接下降至0，如图2中的(b)部分所示。每个手势任务采集5-10次，保证有足够的数据用来测试；

(3)将所采集的数据作为离线数据保存，通常选取5次收缩任务的数据，以保证能够获取到足够的解混向量。故离线数据集时长T

步骤2、采用逐步变量剥离算法分解所述离线数据集的肌电信号：

步骤2.1、利用式(1)对所述离线数据集进行扩展，得到扩展信号

式(1)中，K是延迟因子；

定义并初始化运动单位发放集合γ及其对应的解混向量集合ε为空集；

定义并初始化残差信号

步骤2.2、执行FastICA算法，设置算法执行的总次数为n：

步骤2.2.1、初始化i＝1；

步骤2.2.2、初始化j＝0；

步骤2.2.3、定义并初始化残差信号

步骤2.2.4、利用式(2)计算第i个解混向量在第j+1次的计算结果w

式(2)中，G'表示非多项式函数G的一阶导数；G”表示G的二阶导数；E是数学期望；

步骤2.2.5、如果

W'＝(WW

步骤2.2.6、计算运动单位源信号

步骤2.3：采用图像处理的大津法获取第i个运动单位源信号y

步骤2.3.1、将区间(0,max(y

步骤2.3.2、初始化h＝1；

步骤2.3.3、令th

步骤2.3.4、根据式(4)计算得到类间方差F

F

步骤2.3.5、若h＜M，则令h+1赋值给h，返回至步骤2.3.3；否则进入步骤2.3.6；

步骤2.3.6、找到最大的F

步骤2.4、将y

步骤2.5、通过带时域约束的FastICA算法检测分解结果的可靠性：

步骤2.5.1、定义迭代计算的次数为d，并初始化d＝1；定义并初始化第d次迭代的拉格朗日乘子μ

步骤2.5.2、通过式(5)得到第d+1次迭代的第i个解混向量w

式(5)中，g'

步骤2.5.3、对第d+1次迭代的第i个解混向量w

步骤2.5.4、根据式(6)、式(7)和式(8)分别更新其他参数：

μ

δ

ξ

式(6)、式(7)和式(8)中，μ

步骤2.5.5、更新第d+1次迭代的第i个运动单位源信号

步骤2.5.6、若

步骤2.5.7、令第i个运动单位源信号y

若c

步骤2.6、遍历运动单位发放集合γ并分别计算任意两个不同的发放序列之间的相关系数，若大于C，则表示两个发放序列重复，并随机选取一个从运动单位发放集合γ和解混向量集合ε中删除；否则，两个发放序列均保留；经此操作后，可得到共有q个发放序列的运动单位发放集合γ和共有q个解混向量的解混向量集合ε；

步骤2.7、将运动单位发放集合γ中q个发放序列表示为R＝[r

步骤2.8、将更新后的残差信号

本实施例中，利用了较为代表性的传统的离线分解算法——PFP算法，来识别离线数据集中包含的运动单位。其中，步骤2.1的矩阵扩展部分，延迟因子K选择为5，这与实验的采样率有关；步骤2.2为了尽可能迭代得到较多的解混向量，一般n会取得较大，示例性的，可以设置为n＝200；步骤2.2.3中的非多项式函数通常设置为G(x)＝log(cosh(x))；步骤2.2.5中的

步骤3、实时采集肌电信号，并将获取的肌电数据流按照时间轴分为一系列窗长为L，步长为D的分析窗；在每个分析窗内，根据式(1)得到扩展后的肌电数据X，并利用卷积模型

本实施例中，实时采集肌电信号的实验方案和设备与离线数据集的采集过程保持一致，这里就不再赘述；为了保证窗长设置的足够合理，既能满足在线数据处理的实时性需求，又可以保证分解结果的可靠，示例性的，可以设置为L＝5s，D＝1s；

步骤4、采用自适应的阈值选取策略提取分析窗内第i个运动单位源信号s

步骤4.1、按照步骤2.3中的大津法计算第i个运动单位源信号s

步骤4.2、将s

步骤4.3、初始化k＝1；

步骤4.4、令中间变量th＝peak

步骤4.5、设置向量a，用来存储运动单位发放序列b

cos

式(9)中，std(a)表示向量a的标准差，mean(a)表示向量a的平均值；

步骤4.6、若k＜l，则令k+1赋值给k，返回步骤4.4；否则找到最大的cos

本实验中，采用了间隔变异系数作为运动单位序列质量的衡量指标，通常认为正常的运动单位是以较为稳定的时间间隔募集和发放，故其对应的发放变异系数应当尽可能的小，因而这里寻找到发放变异系数最小的结果作为最终的结果；

步骤5、每当新的肌电数据输入时，按照步骤3和步骤4的过程进行处理和保存；

当没有新的数据输入时，利用式(10)将所保存的P个分析窗中的发放序列按照时间顺序进行拼接，得到总的分解结果Spike并作为最终的结果：

Spike＝[Win

式(10)中，Win

如前所述，实时采集的肌电数据以窗长L、步长D划分分析窗，故两个相邻的分析窗中有L-D的重叠段，因此在连接结果的时候需要去除这部分重复段，仅取每个分析窗结尾为D长度的部分进行连接；此外，在肌电分解领域，通常认为一个解混向量对应一个运动单位，故在本发明直接采用矩阵拼接的形式连接最终的结果。

为了量化评估本发明的效果，通过8位不同受试者的实验数据将本发明方法与传统的离线分解方法(即步骤2内容)相比较，所有的数据处理过程在软件MATLAB 2020a上运行。

对比实验中，在线处理了20s的表面肌电信号，采集过程中共执行了4次的拇指收缩任务。在线过程的20s的信号除了完成了在线识别运动单位发放序列，还应用步骤二所详述的离线分解算法进行处理(将数据分为4个5s的数据段分别分解，每个数据段对应于1次收缩任务)。图3是其中一位受试者一个5s数据段的分解结果示意图。从图上可以看出，在线处理的结果和离线的结果很好地匹配上了，存在较少的错误或者遗漏的发放。图4展示了8位受试者的在线处理结果。可以发现，离线分解算法共获取到了10.31±1.46个运动单位，在线处理过程共获取到了11.38±2.50个运动单位，共有6.56±1.83个运动单位被成功匹配上了，这些运动单位称之为离线和在线过程共有的运动单位。在这些共有的运动单位中，匹配度(Matching Rate,MR)平均为(92.65±3.77)％，错误发现率(False DiscoveryRate,FDR)为0.066±0.030，假阴性率(False Negative Rate,FNR)为0.074±0.044。可以发现，本发明在不同的受试者上均能与离线的分解算法达到90％以上的匹配度，错误率维持在一个较小的状态，且运行速度相较于离线分解得到大幅的提升，平均每个分析窗的处理时间从62.44秒缩减到0.53秒，证明了所提出的在线识别运动单位方法的可行性和可靠性。

综上所述，本发明能够解决肌电分解难以在线运行的问题，将离线分解的结果作为先验知识指导在线的处理过程，并且通过自适应的阈值选取策略，对结果进行修正和完善，以提升整体的精度。与传统的离线分解算法相比，本发明能够更加快速地获取到运动单位的活动信息，在肌电控制方面有着更广泛的应用面，对运动控制和医疗康复等领域具有重要意义。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例并可以通过相应的软件和代码实现。上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例的方法。