一种手腕关节运动力矩的肌电估计方法

摘要

本发明公开了一种手腕关节运动力矩的肌电估计方法。该方法采集前臂6块肌肉在用户手腕关节弯曲伸展过程中的表面肌电信号，利用TKE算子计算肌肉收缩时间。利用等速肌肉力测试系统与肌电仪同步测量各肌肉等长收缩最大张力。建立腕关节正向肌骨模型，输入肌肉收缩时间与各肌肉等长收缩最大张力，输出腕关节弯曲伸展极限位置处估计力矩。利用三维运动捕捉系统获得腕关节弯曲伸展过程中运动学数据，并解算出极限位置处参考力矩。将模型估计力矩与参考力矩误差的平方作为目标函数，利用共轭梯度法完成腕关节正向肌骨模型的标定，实现表面肌电信号对手腕关节力矩的估计。本发明可应用于肌电假手、康复医疗、生机电一体化等领域。

说明书

【技术领域】

本发明属于智能假手与生机电一体化技术领域，涉及一种手腕关节运动力矩的肌电估计方法。

【背景技术】

目前，由于意外、疾病等原因造成的手部缺失患者人数已达千万。现有医疗水平还无法实现手部再生，故假肢是手部缺失患者恢复手部功能的唯一方式。肌电假手是基于人手臂表面肌电信号特征感知人体运动意图，控制假手实现相应操作，从而满足手部残疾患者日常生活所需。与传统假手相比，其功能更加完善，能够更好的提升患者使用的幻肢感。

大脑产生的兴奋信号经神经系统传递至肌肉纤维，从而产生动作电位。一方面动作电位相互叠加在人体表面形成表面肌电信号；另一方面动作电位沿着肌纤维向多个方向传播，引起肌纤维收缩，产生肌肉力并带动骨骼绕关节运动。又因为表面肌电信号的产生通常先于实际动作发生，所以表面肌电信号在某种程度上可以预测人体运动意图。肌电假手正是利用表面肌电信号与人体运动意图之间的关系来帮助手部缺失患者满足日常生活与工作所需。现有肌电假手在人体运动意图感知方面主要聚焦于假手抓握动作的分类，忽略抓握过程中力矩的研究。为了保证假手抓握的稳定性，除了提供合适的手部抓握力，对于腕关节的抓握力的估计也十分重要。基于此，研究前臂表面肌电信号对用户手腕关节运动力矩的估计是当前肌电假手研究领域的关键点，而其中极限位置处最大力矩的估计显得尤为重要。

【发明内容】

本发明为了增加表面肌电信号对于人体运动意图感知的方面，提供一种手腕关节运动力矩的肌电估计方法，该方法通过腕关节正向肌骨模型建立表面肌电信号与用户手腕关节弯曲伸展至极限位置时力矩之间映射关系，实现输入完整的弯曲伸展运动过程中前臂表面肌电信号，输出对应的用户手腕关节在弯曲伸展极限位置时的力矩值。

为达到上述目的，本发明采用以下技术方案予以实现：

一种手腕关节运动力矩的肌电估计方法，包括以下步骤：

步骤1：建立用户手腕关节正向肌骨模型；

步骤2：利用多关节等速肌力测试系统与表面肌电仪同步测量前臂各肌肉等长收缩最大张力；

步骤3：在利用三维运动捕捉系统采集用户手腕关节3种运动速度下腕关节运动学数据的同时，利用表面肌电仪同步采集人手前臂6块肌肉表面肌电信号，并进行预处理，利用TKE算子求出肌肉收缩时间；

步骤4：利用三维运动捕捉系统采集用户手腕关节运动学数据，建立用户手腕关节简化模型，进行逆动力学求解，得到用户手腕关节弯曲伸展极限位置处力矩；

步骤5：利用共轭梯度法完成腕关节正向肌骨模型标定。

与现有技术相比，本发明具有以下有益效果：

本发明易于操作并保证高精确性。多关节等速肌力测试系统可提供等长收缩张力测量模式，用于测量关节在多个固定角度下的峰力矩，由此得到的肌肉等长收缩最大张力精度高。三维运动捕捉系统具有高分辨率与高采样率，由此捕捉到的腕关节运动学数据可靠性高。表面肌电仪具有高采样率，可实时检测表面肌电信号反应肌肉收缩情况，映射肌肉力大小。三种仪器均具备同步采集接口，两两之间可实现同步采集。

建立的表面肌电信号估计用户手腕关节运动力矩的正向肌骨模型具有更好的实际应用性。本发明通过对Hill肌肉力模型进行改进，引入表面肌电信号，根据表面肌电信号起止时刻计算肌肉收缩时间，从而建立表面肌电信号与关节力矩之间的联系。并利用数值优化的思想寻求模型理想参数，避免传统肌肉力模型中生物力学参数不可直接测量带来的实际应用性差的现象。该模型的建立在医疗康复、人机工程等领域提供了新的思路。

【附图说明】

图1是人手前臂表面肌电信号对用户手腕关节运动力矩估计的整体方法框图；

图2是用户手腕关节动力学模型简图；

图3是6块肌肉在弯曲伸展整个过程中肌肉收缩时间的计算结果；

图4是腕关节正向肌骨模型标定方法图；

图5是共轭梯度法优化模型参数流程图。

【具体实施方式】

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，不是全部的实施例，而并非要限制本发明公开的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要的混淆本发明公开的概念。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的，其中为了清楚表达的目的，放大了某些细节，并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的，实际中可能由于制造公差或技术限制而有所偏差，并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。

本发明公开的上下文中，当将一层/元件称作位于另一层/元件“上”时，该层/元件可以直接位于该另一层/元件上，或者它们之间可以存在居中层/元件。另外，如果在一种朝向中一层/元件位于另一层/元件“上”，那么当调转朝向时，该层/元件可以位于该另一层/元件“下”。

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

下面结合附图对本发明做进一步详细描述：

参见图1，本发明建立腕关节正向肌骨模型，实现输入表面肌电信号输出估计的关节力矩。一方面采集用户手腕关节弯曲伸展运动过程中前臂肌肉的表面肌电信号，求取肌肉收缩时间，作为模型输入，基于模型给定的参数初始值，得到腕关节弯曲伸展至极限位置时的估计力矩。另一方面，利用三维运动捕捉系统采集到的腕关节运动学数据进行逆动力学解算，得到腕关节运动计算力矩。将计算力矩作为模型参考值，取计算力矩与估计力矩误差的平方作为目标函数，利用共轭梯度法完成模型离线参数学习，最终得到稳定的腕关节正向肌骨模型，实现再次输入表面肌电信号，输出对应的腕关节运动力矩。具体实施方案如下：

第一步：建立用户手腕关节正向肌骨模型，具体方法如下：

(1)基于表面肌电信号的单块肌肉力计算：

目前，应用比较广泛的肌肉力模型为Hill三元素模型，该模型由收缩单元、并联弹性单元和串联弹性单元混联构成。收缩单元代表肌节中的肌动蛋白微丝及肌球蛋白微丝。其在兴奋时可产生主动张力。串联弹性单元代表肌微丝、横桥闰盘及两端的腱结构。当收缩单元兴奋后，串联弹性单元使肌肉具有弹性。并联弹性单元代表肌束膜及肌纤维膜等结缔组织。当其被牵拉时会产生弹力，是一种被动张力。整块肌肉可认为是由许多这样的模型混联在一起构成。

Hill方程可描述为：

(a+F)(V+b)＝b(F₀+a)

简化为：

式中，a、b为实验参数，F₀是等长收缩最大张力，V是肌肉收缩速度，F₀、V是肌肉纤维初始长度、温度、周围环境的化学组成及浓度等的函数。

肌肉收缩速度V可表示为：

式中，t为肌肉收缩时间，Δl为肌肉收缩长度。一般肌肉达到最大力值所需时间为300～400ms，而许多运动中力作用时间比该时间短得多，直接测量困难。为了解决这一问题，可以采用表面肌电信号持续时间作为肌肉收缩时间，而表面肌电信号持续时间可以利用TKE算子求出。肌肉收缩长度作为实验参数进行优化。

则单块肌肉力最终简化表达式为：

(2)多块肌肉共同作用下腕关节力矩求解：

根据骨杠杆原理，单块肌肉收缩力在关节上产生的转动力矩为：

T_i＝r_i×F_i·cos(φ_i)i＝1,2,…,N

式中，r_i为关节转动中心到力作用点的位移矢量；F_i为肌肉力矢量；φ_i为肌肉羽状角；N为肌肉数量。

人体上肢前臂肌肉多为梭形肌紧附在前臂尺骨和桡骨，故羽状角φ_i与r_i较小，且运动过程中发生微小变化，故直接测量不便，为了使腕关节正向肌骨模型简单易用，将其转换为：

T_i＝k_i·F_i>

式中，k_i＝r_i·cos(φ_i)为比例系数。

确定每块肌肉对关节的力矩贡献并相加即可得到腕关节力矩，即：

至此得到腕关节正向肌骨模型，实现表面肌电信号估计腕关节运动力矩，即：

第二步：利用多关节等速肌力测试系统与表面肌电仪同步测量前臂各肌肉等长收缩最大张力，具体方法如下：

选取受试者右侧腕关节力矩为测试对象，受试者实验前进行适当热身，并用酒精擦拭右侧前臂。为了保证实验操作规范化，测试数据准确性，受试者按照多关节等速肌力测试系统操作手册进行姿态固定与关节对位，同时选取前臂桡侧腕长伸肌、桡侧腕屈肌、尺侧腕长伸肌、尺侧腕屈肌、指总伸肌和指浅屈肌为肌肉测点，利用表面肌电仪同步采集前臂6块肌肉表面肌电信号。

规定手掌与前臂保持水平为0°，用户手腕关节弯曲伸展范围为腕伸60°与腕屈75°之间。实验过程中，选取等长张力测量模式，分别测量受试者在腕伸30°、60°与腕屈30°、60°下屈肌和伸肌的峰力矩。测试时，每个角度重复3次，要求受试者1-2s内达到最大力量并保持5s，得到各角度下屈肌与伸肌的峰力矩。分别选取屈肌与伸肌在每个角度下所得最大峰力矩作为该角度下的峰力矩，再比较屈肌和伸肌在四种角度下峰力矩，选取最大值作为最大峰力矩，即屈肌和伸肌的等长收缩最大张力。同步采集人手前臂6块肌肉表面肌电信号，对表面肌电信号进行预处理之后提取表面肌电活跃强度特征值，并计算各通道表面肌电信号活跃强度最大幅值。计算桡侧腕长伸肌、尺侧腕长伸肌、指总伸肌的表面肌电信号活跃强度最大幅值之比，将伸肌最大峰力矩乘以各肌肉肌电活跃强度最大幅值占比计算出各伸肌的等长收缩最大张力。同理得到桡侧腕屈肌、尺侧腕屈肌、指浅屈肌的等长收缩最大张力。

第三步：利用表面肌电仪采集用户手腕关节3种运动速度下人手前臂6块肌肉表面肌电信号，并进行预处理，利用TKE算子求出肌肉收缩时间。具体方法如下：

(1)表面肌电信号的采集与预处理：

选取与用户手腕关节弯曲伸展运动相关的前臂6块肌肉，分别为桡侧腕屈肌、桡侧腕长伸肌、尺侧腕屈肌、尺侧腕长伸肌、指总伸肌、指浅屈肌。利用表面肌电仪采集其表面肌电信号，并对其进行带通滤波与陷波滤波，去除噪声干扰。

(2)肌肉收缩时间计算：

利用TKE算子检测表面肌电信号的起止时刻，由此计算表面肌电信号持续时间作为肌肉收缩时间。表面肌电信号可以表征为一串离散的数字信号，对于给定的离散信号，TKE算子ψ(n)可描述为：

其中，表示对表面肌电信号进行去均值处理，N为表面肌电信号序列总长度，M为背景噪声信号长度。

在表面肌电信号起始时刻检测过程中，需要根据TKE算子的均值与均方差设置合适的阈值。TKE算子均值与均方差计算公式分别如下：

得到表面肌电信号阈值为：

Th＝μ₀+jδ₀

其中，j为阈值乘子，通过调节j的取值来选取合适的阈值进行肌肉运动与否的判断。

j∈[5,7],j＝5

对表面肌电信号的TKE算子ψ(n)与阈值Th进行比较，得到二值化状态函数s(n)。

s(n)为一串0、1相间的表示肌肉收缩状态的序列。由于表面肌电信号受噪声干扰较严重，在动作未发生时将尖峰噪声信号误判为肌肉动作发生；而在肌肉持续收缩过程中，由于表面肌电信号不稳定，也可能将静息状态判为肌肉动作未发生。为了去除以上两种情况带来的误差，需要对s(n)做进一步处理：

首先，将s(n)中间隔小于T₁的两个1之间的0全部置为1，用于移除肌肉收缩频率过快引起表面肌电信号不稳定带来的误差。

然后，将s(n)中间隔小于T₂的两个0之间的1全部置为0，用于移除尖峰信号带来的干扰。

T₁表示有动作发生时，出现在动作持续区内部正常的非活动标识间隔；T₂表示无动作发生时，出现在非活动区内正常的伪活动标识间隔。利用TKE算子检测表面肌电信号起止时刻情况如参考图3所示。

则肌肉收缩时间为：

t＝max(n|s(n)＝1)-min(n|s(n)＝1),n＝1,2,…,M,…,N

第四步：在步骤三的过程中，利用三维运动捕捉系统同步采集用户手腕关节在3种运动速度下运动学数据，建立用户手腕关节简化模型，进行逆动力学求解，得到用户手腕关节弯曲伸展极限位置处力矩。具体方法如下：

(1)设计用户手腕关节运动学信息采集实验。受试者采取坐姿状态，大臂与小臂保持水平，手腕跟随指示视频做三种不同速度下弯曲伸展运动。可以将人手与手臂看成刚体，分别在其上粘贴三个不共线的mark点，并利用三维运动捕捉系统来获得用户手腕关节在空间中位姿变化的运动学数据。

(2)建立腕关节简化模型如参考图2所示，对腕关节位姿信息进行逆动力学求解，计算腕关节运动力矩，并选取弯曲伸展过程中最大值最为腕关节弯曲伸展极限位置处关节力矩值。将人体上肢简化成连杆形式，人手与腕关节看作是单一连杆通过铰链与竖直面相连接。自由度数相对较少，刚体数量只有一个，所以选用拉格朗日方法对用户手腕关节进行逆动力学求解。

拉格朗日动力学描述基于系统能量的概念。对于任何机械系统，拉格朗日函数L定义为系统总动能E_k与总势能E_p之差，即：

式中，q＝[q₁,q₂,…,q_n]是表示动能和势能的广义坐标，是相应的广义速度。

利用拉格朗日函数L，系统的动力学方程为：

式中，τ是n×1的关节驱动力矩矢量。由于势能Ep不含因而动力学方程变为：

对于人体上肢腕关节的简化模型如参考图2所示，设定手腕广义坐标为弯曲为正，伸展为负。

手部连杆的转动惯量为：

手部连杆动能和势能依次为：

式中，l₀为手掌质心与腕关节的距离。则：

因此腕关节弯曲/伸展力矩M为：

第五步：利用共轭梯度法完成腕关节正向肌骨模型标定，具体方法如下：

(1)确定目标函数：

腕关节正向肌骨模型最终表达式中肌肉等长收缩最大张力在第二步中已经确定，其余参数需要进一步优化。赋予a、b、k_i、Δl初始值，按照第四步的方法计算肌肉收缩时间并输入正向肌骨模型，输出对应的腕关节估计力矩。将第四步通过逆动力学解算出的计算力矩作为模型输出的参考值，求出估计力矩与参考力矩的误差的平方作为目标函数，则目标函数如下所示：

式中，n——整个样本中样本点个数；T_est(j)——由正向肌骨模型估计的第j极限位置处腕关节估计力矩，T_cal(j)——由运动学数据逆动力学解算出的第j极限位置处腕关节参考力矩。

将目标函数进一步简化成如下表达式：

Y＝f(k_i,Δl_m,a,b)

(2)基于共轭梯度法完成目标函数参数寻优：

目标函数与各参数之间存在非线性关系，所以采用共轭梯度法对目标函数进行优化，确定合适的参数使目标函数达到最小值。腕关节正向肌骨模型标定过程如参考图4所示。

设置迭代精度ε＝0.1和各参数初始值x⁰＝[k_i⁰,Δl⁰,a⁰,b⁰]^T，计算初始点处梯度，则第一次沿初始点负梯度方向进行搜索。

若||f(x⁰)||≤ε，则停止迭代，输出方程的解x^*＝x⁰。

否则确定新的搜索方向dⁿ⁺¹和步长β_n继续迭代：

从而得到新的搜索点：

xⁿ⁺¹＝xⁿ+β_ndⁿ

计算新的搜索点的梯度检查是否满足迭代精度要求，即：

共轭梯度法优化模型参数流程如参考图5所示。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。