CN110355761B - 一种基于关节刚度和肌肉疲劳的康复机器人控制方法 - Google Patents

Abstract

本发明公开了一种基于关节刚度和肌肉疲劳的康复机器人控制方法。本发明采集受试者的关节角度信号和表面肌电信号，结合正、逆动力学原理，通过遗传算法辨识肌肉的个性化生理参数，建立个性化关节骨骼肌肉模型，计算运动过程中的关节刚度信息；计算运动过程中表面肌电信号的中值频率，并利用其相对变化值获取受试者的疲劳信息；采用关节刚度信息和运动疲劳信息对阻抗模型的参数进行自适应调整，同时通过饱和函数对刚度和阻尼参数进行约束，实现康复机器人自适应阻抗控制方法。本发明同时考虑与关节相关的刚度信息和肌肉相关的疲劳信息，将两者引入康复机器人控制，能够兼顾关节和肌肉在康复训练过程中的安全性，实现安全有效的康复训练。

Description

一种基于关节刚度和肌肉疲劳的康复机器人控制方法

技术领域

本发明属于基于表面肌电信号的康复机器人控制领域，尤其一种基于关节刚度和肌肉疲劳的康复机器人控制方法。

背景技术

康复机器人主动康复控制的实现需要考虑机器人和患者之间的人机交互作用。通过分析患者在康复训练过程中产生的生物传感信息，可以获取患者的运动意图，并由此控制机器人相应的运行状态，可为患者提供主动安全的康复训练。交互力信号和生理信号是典型的人机交互信息。交互力信号的获取通常需要依托康复机器人的机械结构，不够方便灵活，且不能及时反馈患者信息，可能导致闭环控制的时延。生理信号不仅可以判断患者的运动意图，还可用于评估康复训练过程中肢体的健康状态，使得基于生理信号的康复机器人主动康复控制得到广泛的研究。目前基于生理信号的康复机器人主动控制大多关注人体的运动意图的识别，而对于康复信息的提取如关节刚度、肌肉疲劳状态等缺少进一步研究。

人体关节刚度是关节的一个重要物理特征，其在运动过程中变化很大，在康复训练过程中康复机器人需对这种变化做出反应，以提供更加安全舒适的控制方式；使用康复机器人对患者进行康复训练能够在一定程度上恢复其肢体运动功能，但随着运动进行，患者易出现肌肉疲劳或疲劳加深的状况，过度训练造成的肌肉疲劳会对患者造成二次损伤，因此需在康复过程中实时获取患者的肌肉疲劳信息，并根据疲劳信息调节康复机器人的运行状态。基于此，在得到患者的关节刚度信息和运动疲劳信息之后需要设计可靠的康复机器人控制方法，使之能安全地帮助患者完成康复训练。

中国专利CN109645962.A公开了一种基于疲劳感知的机器人辅助康复人机协作训练方法，通过表面肌电信号进行疲劳特征提取，并检测患者的疲劳状态：正常/轻度疲劳/重度疲劳，根据疲劳状态和决策控制器切换康复训练模式，从而避免过度训练对患者造成损伤。该控制方法仅是将患者的疲劳状态当作控制开关来切换康复训练的状态，缺乏人机交互过程中疲劳信息的持续反馈，无法在两种状态之间充分利用疲劳信息进行人机交互控制，且忽略了关节信息在康复训练中的变化，无法保障康复训练中的关节安全。中国专利CN106109174.A公开了一种肌电反馈式阻抗自适应的康复机器人控制方法，使用表面肌电信号辨识关节的屈伸状态，通过关节角度和肌肉活动水平自适应调整阻抗参数，同时依据疲劳程度调节初始期望静态平衡力，实现康复机器人的自适应控制。该控制方法使用患侧镜像健侧进行肌肉活动水平的计算来调节阻抗参数，仅适用于单侧损伤的患者，而且疲劳程度的分级难以定量分析，使得该方法缺乏普适性。中国专利 CN109718059.A公开了一种手部康复机器人自适应控制方法及装置，其中的控制方法通过运动想象脑电信号获取手部康复轨迹，同时利用表面肌电信号获取上肢活动度从而建立变阻抗控制模型。该控制方法中的参数调整仅考虑康复机器人的自适应问题，未考虑安全性的影响，使得该方法存在安全隐患；而且需要同时获取脑电信号和表面肌电信号进行实时处理，此过程的实施较为困难。

通过上述分析可知，目前主动康复控制器的设计大多依据肌肉的信息，忽略关节在康复过程中的影响，且没有兼顾两者的安全性。为了克服传统康复机器人主动控制的局限性，提高康复训练的有效性和安全性，本发明利用表面肌电信号获取人体的关节刚度信息和运动疲劳信息，在此基础上设计基于位置的阻抗控制模型，实现基于刚度估计和疲劳反馈的康复机器人自适应阻抗控制，为患者的康复训练提供自然、有效和安全的控制手段。

发明内容

为了解决上述技术问题，本发明提出了一种基于关节刚度和肌肉疲劳的康复机器人控制方法。

本发明的技术方案为一种基于关节刚度和肌肉疲劳的康复机器人控制方法，具体包括以下步骤：

步骤1：采用运动采集与分析系统对表面肌电信号和关节角度信号进行同步采集，并分别对表面肌电信号和关节角度信号进行预处理，得到预处理之后的表面肌电信号和实际关节角度信号；

步骤2：结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型，将预处理后表面肌电信号EMG_ap和实际关节角度信号θ分别输入关节刚度估计模型和关节逆向动力学模型，利用遗传算法辨识出与受试者相匹配的生理参数，得到个性化关节刚度估计模型；

步骤3：采集关节动作最大自主收缩时的表面肌电信号，通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度，结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息；

步骤4：采集康复训练过程中的人机交互力矩，建立阻抗模型，设定康复机器人的刚度参数和阻尼参数的基础值，实现基于位置的阻抗控制模型；

步骤5：设定康复机器人运行的期望轨迹，将关节刚度信息和运动疲劳信息引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数，得到实际的运行轨迹。

作为优选，步骤1中所述采用运动采集与分析系统对表面肌电信号进行同步采集为EMG；

步骤1中所述采用运动采集与分析系统对关节角度信号进行同步采集并进行计算为：

在关节转动时固定肢体的末端、关节转动中心和关节转动时转动肢体的末端分别贴上标记点，对应记为点A(x_A,y_A,z_A)、点B(x_B,y_B,z_B)和点C(x_C,y_C,z_C)，记

为 B点和A点构成的向量，

为B点和C点构成的向量；

通过运动采集与分析系统中标记点获取关节角度信号，通过标记点构成的向量计算关节角度信号θ_a：

步骤1中所述对表面肌电信号进行预处理为：

采用阶数为5阶，截止频率为20Hz的巴特沃斯高通滤波器对表面肌电信号进行预处理，以减少基线漂移和伪影噪声的影响，并得到预处理后表面肌电信号EMG_ap；

步骤1中所述对关节角度信号进行预处理为：

θ＝θ_a-θ₀

其中，θ为实际关节角度信号，θ₀为标记点的初始角度误差。

作为优选，步骤2中所述结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型：

引入步骤1中预处理后表面肌电信号EMG_ap取绝对值并进行归一化处理，得到归一化表面肌电信号e(t)如下：

其中，EMG_rea＝|EMG_ap|为取绝对值之后的表面肌电信号，EMG_mvc为最大自主收缩期间记录的处理过的表面肌电信号峰值，EMG_res为静息状态下的表面肌电信号；

然后使用二阶离散线性模型得到神经激活度u(t)如下：

u(t)＝0.9486·e(t-de)+0.052·u(t-1)-0.000627·u(t-2)

其中，de＝80ms。最后获取肌肉激活度a(t)如下：

引入步骤1中实际关节角度信号θ，通过肌肉路径方程和个性化骨骼肌肉几何学模型计算肌肉长度l_mt和肌肉力力臂长度d如下：

β＝σ+θ

d＝l_PRl_QR sinβ/l_PQ

其中，P为肌肉起点，Q为肌肉止点，R为关节转动中心，σ为初始状态肌肉起点和止点与转动中心之间的角度；

由Hill肌肉肌腱模型结合给定的生理参数计算肌肉力F如下：

其中，

是最佳肌肉纤维长度下的最大收缩力，a为肌肉激活度，l为标准化的肌肉纤维长度，v为标准化的肌肉纤维收缩速度，F_A(l)为标准化的主动力-肌纤维长度函数，F_V(v)为标准化的主动力-肌纤维收缩速度函数，F_p(l)为标准化的被动力-肌纤维长度函数，α称为羽状角，人体肌肉在运动的过程中，羽状角的变化引起的肌肉的输出力的变化很小，所以在计算肌肉力时，可以忽略因羽状角的变化而引起的肌肉力的变化。

使用肌肉力和肌肉力力臂长度计算出关节力矩τ_a如下：

引入步骤1中实际关节角度信号θ，使用关节刚度估计模型计算关节转动的刚度信息K_a如下：

其中，i为参与计算的第i块肌肉，τ_a为关节力矩。

引入步骤1中实际关节角度信号θ，使用逆向动力学模型计算参考力矩τ_r如下：

其中，g为重力加速度，I、M、l分别为关节转动时转动肢体的转动惯量、质量以及重力力臂长度。

引入步骤1中实际关节角度信号θ，通过参考力矩计算参考刚度K_r如下：

以刚度K_a和参考刚度K_r之间的误差最小化为优化目标，为每块肌肉设置对应的最佳肌肉纤维长度下的最大收缩力

最佳肌纤维长度l_mopt和肌腱长度调节因子s_t的初始值及寻优范围，通过遗传算法的选择、交叉和变异操作对最佳肌肉纤维长度下的最大收缩力、最佳肌纤维长度和肌腱长度调节因子进行寻优，从而得到优化后参数即受试者个性化生理参数，以构建个性化关节刚度估计模型如下：

其中，

l_mopt和s_t对每块肌肉有与之对应的寻优范围，h为数据点，H为总的数据点数。

作为优选，步骤3中所述采集关节动作最大自主收缩时的表面肌电信号为EMG_m：

步骤3中所述通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度为：

假设有Q种不同关节动作，对这Q种动作采集的最大自主收缩时的表面肌电信号分别提取均方根值(Root Mean Square，RMS)作为肌肉活动水平如下：

其中，x_i,q代表了第q个动作的表面肌电信号序列，n_q表示第q个序列长度。

采用肌肉协同理论计算第i肌肉对第q个动作的贡献度WD_iq如下：

其中，N为肌肉数量，W_i为第i块肌肉在第q个动作的肌肉协同矩阵，并取所有关节动作的贡献度均值作为第i块肌肉的平均贡献度c_i如下：

步骤3中所述结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息为：

计算第i块肌肉的中值频率反映肌肉的疲劳信息，并通过其相对变化值结合平均贡献度获取运动疲劳信息p如下：

其中，p_i为第i块肌肉的疲劳信息，p_0i为第i块肌肉的初始疲劳信息。

作为优选，步骤4中所述建立阻抗模型为:

采集人机交互力矩τ，通过经验设定刚度参数的基础值即K₀和阻尼参数得基础值即B₀，由于康复机器人的运行速度通常比较匀速缓慢，加速度项通常可以忽略不计，可以得到阻抗模型如下：

其中，

和q分别表示速度和轨迹，

和q_d分别表示期望速度和轨迹；

步骤4中所述基于位置的阻抗控制模型为：

作为优选，步骤5中所述设定康复机器人运行的期望轨迹为：

引入步骤1中的关节刚度估计模型计算关节刚度信息k，引入步骤3中的运动疲劳信息p，结合饱和函数，将关节刚度信息k和运动疲劳信息p对刚度参数和阻尼参数的基础值进行调节，得到自适应变化的刚度参数K和阻尼参数B如下：

其中，B₀为阻尼的基础值，K₀为刚度参数的基础值，B_l0为阻尼参数的下界， B_h0为阻尼参数的上界，K_l0为刚度参数的下界，K_h0为刚度参数的上界，p为疲劳信息，k为刚度信息，η为疲劳信息的权重因子，υ为刚度信息的权重因子；

步骤5中所述引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数，得到实际的运行轨迹为：

引入步骤4中基于位置的阻抗控制模型，将自适应变化的刚度参数和阻尼参数替换刚度参数和阻尼参数的基础值，得到自适应阻抗控制模型如下：

将上式

中的机器人运行速度

进行积分操作，可得到实际的运行轨迹q。

与现有技术相比，本发明采用上述方案能够产生这样的有益效果：

本发明采集人体表面肌电信号和关节角度信号，建立个性化关节刚度估计模型，同时获取康复训练过程中的运动疲劳信息。用关节刚度信息调节康复机器人的刚度参数，用运动疲劳信息调节阻尼参数，由此提供一种基于刚度估计和疲劳反馈的康复机器人自适应阻抗控制方法。随着关节刚度的增大，康复机器人的刚度自适应地变小，确保在关节刚度较大的角度也能够让患者通过自己的努力调整运动姿态，防止康复机器人刚度过大对关节有害；同时随着疲劳的加深，阻尼参数也自适应地减小，使得患者能够更加自由地通过自己的主动力影响康复机器人的运动，在保证康复训练积极性的同时能够提高康复训练的安全性。

附图说明

图1为本发明基于关节刚度和肌肉疲劳的康复机器人自适应控制方法的原理图；

图2为脚踝的表面肌电信号电极分布示意图；

图3为脚踝的运动采集与分析系统标记点的空间位置；

图4为脚踝屈伸刚度估计模型生理参数辨识原理图。

2-1为腓肠内肌电极位置、2-2为腓肠外肌电极位置、2-3为胫骨前肌电极位置、2-4为比目鱼肌电极位置、3-1为膝关节转动中心标记点、3-2为脚踝转动中心标记点、3-3为足尖标记点、3-4为脚踝屈伸角度、3-5为运动采集与分析系统标记原点。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合图1至图4，以脚踝屈伸动作为例，介绍本发明的具体实施方式为：

结合图1，本发明所述的基于刚度估计和疲劳反馈的脚踝康复机器人自适应阻抗控制方法采用基于位置的阻抗模型，首先设定参考轨迹，在康复训练过程中采集表面肌电信号和关节角度信号，获取关节刚度信息和运动疲劳信息，分别对刚度参数和阻尼参数进行自适应调节；使用力矩传感器采集人机交互力矩信息，结合自适应刚度参数和自适应阻尼参数重塑参考轨迹生成实际轨迹；将实际轨迹输入康复机器人位置控制器中，使康复机器人在患者主动运动意图之下带动患者肢体运动，实现安全有效的康复训练。

步骤1：采用运动采集与分析系统对表面肌电信号和关节角度信号进行同步采集，并分别对表面肌电信号和关节角度信号进行预处理，得到预处理之后的表面肌电信号和实际关节角度信号；

结合图2，与脚踝屈伸运动相关的肌肉群主要包括腓肠内肌(MedialGastrocnemius，MG)，腓肠外肌(Lateral Gastrocnemius，LG)，胫骨前肌 (TibialisAnterior，TA)，比目鱼肌(Soleus，SO)，选取这四块肌肉作为表面肌电信号的采集点，具体位置如图2中2-1、2-2、2-3和2-4所示。

步骤1中所述采用运动采集与分析系统对表面肌电信号进行同步采集为 EMG；

步骤1中所述采用运动采集与分析系统对关节角度信号进行同步采集并进行计算为：

结合图3，为了计算脚踝屈伸角度，以3-5(记为O点)为坐标系原点建立一个三维空间，从而得到三个标记点的空间坐标。脚踝转动的固定肢体为小腿，转动肢体为足，记3-1为小腿末端标记点A(x_A,y_A,z_A)，记3-2为脚踝转动中心标记点B(x_B,y_B,z_B)，记3-3为足部末端标记点C(x_C,y_C,z_C)，记

为B点和A点构成的向量，

为B点和C点构成的向量，通过运动采集与分析系统中标记点A、 B和C获取脚踝屈伸角度3-4(记为θ_a)可通过这两个向量计算而得：

步骤1中所述对表面肌电信号进行预处理为：

采用阶数为5阶，截止频率为20Hz的巴特沃斯高通滤波器对表面肌电信号进行预处理，以减少基线漂移和伪影噪声的影响，并得到预处理后表面肌电信号EMG_ap；

步骤1中所述对关节角度信号进行预处理为：

θ＝θ_a-θ₀

其中，θ为实际关节角度信号，θ₀为标记点的初始角度误差，也即足部与小腿垂直时标记点计算得到的初始角度。

步骤2：结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型，将预处理后表面肌电信号EMG_ap和实际关节角度信号θ分别输入关节刚度估计模型和关节逆向动力学模型，利用遗传算法辨识出与受试者相匹配的生理参数，得到个性化关节刚度估计模型；

作为优选，步骤2中所述结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型：

引入步骤1中预处理后表面肌电信号EMG_ap取绝对值并进行归一化处理，得到归一化表面肌电信号e(t)如下：

其中，EMG_rea＝|EMG_ap|为取绝对值之后的表面肌电信号，EMG_mvc为最大自主收缩期间记录的处理过的表面肌电信号峰值，EMG_res为静息状态下的表面肌电信号；

然后使用二阶离散线性模型得到神经激活度u(t)如下：

u(t)＝0.9486·e(t-de)+0.052·u(t-1)-0.000627·u(t-2)

其中，de＝80ms。最后获取肌肉激活度a(t)如下：

将肌肉和骨骼简化为几何线条，使用数学方法建立骨骼肌肉几何模型，根据人体的形态学参数获取肢体的相关参数，结合肌肉路径评估方程，可得到运动过程中的肌肉长度和肌肉力力臂长度。引入步骤1中实际关节角度信号θ，通过肌肉路径方程和个性化骨骼肌肉几何学模型计算肌肉长度l_mt和肌肉力力臂长度d 如下：

β＝σ+θ

d＝l_PRl_QR sinβ/l_PQ

其中，P为肌肉起点，Q为肌肉止点，R为关节转动中心，σ为初始状态肌肉起点和止点与转动中心之间的角度；

由Hill肌肉肌腱模型结合给定的生理参数计算肌肉力F如下：

其中，

是最佳肌肉纤维长度下的最大收缩力，a为肌肉激活度，l为标准化的肌肉纤维长度，v为标准化的肌肉纤维收缩速度，F_A(l)为标准化的主动力-肌纤维长度函数，F_V(v)为标准化的主动力-肌纤维收缩速度函数，F_p(l)为标准化的被动力-肌纤维长度函数，α称为羽状角，人体肌肉在运动的过程中，羽状角的变化引起的肌肉的输出力的变化很小，所以在计算肌肉力时，可以忽略因羽状角的变化而引起的肌肉力的变化。

在得到肌肉力和相应的力臂之后便可计算肌肉力矩，多块肌肉的力矩之和就可得到关节力矩，取脚踝背伸动作的方向为正方向，使用肌肉力和肌肉力力臂长度计算出脚踝屈伸的关节力矩τ_a如下：

其中F_i表示肌肉i所产生的的肌肉力，d_i表示肌肉i的力臂长度。引入步骤 1中实际关节角度信号θ，使用关节刚度估计模型计算关节转动的刚度信息K_a如下：

引入步骤1中实际关节角度信号θ，使用逆向动力学模型计算参考力矩τ_r如下：

其中，g为重力加速度，I、M、l分别为关节转动时转动肢体的转动惯量、质量以及重力力臂长度。

引入步骤1中实际关节角度信号θ，通过参考力矩计算参考刚度K_r如下：

结合图4，以刚度K_a和参考刚度K_r之间的误差最小化为优化目标，为腓肠内肌，腓肠外肌，胫骨前肌和比目鱼肌分别设置对应的最佳肌肉纤维长度下的最大收缩力

最佳肌纤维长度l_mopt和肌腱长度调节因子s_t的初始值及寻优范围，通过遗传算法的选择、交叉和变异操作对每块肌肉的最佳肌肉纤维长度下的最大收缩力、最佳肌纤维长度和肌腱长度调节因子进行寻优，从而得到优化后参数即受试者个性化生理参数，以构建个性化关节刚度估计模型如下：

其中，

l_mopt和s_t对每块肌肉有与之对应的寻优范围，h为数据点，H为总的数据点数。

步骤3：采集关节动作最大自主收缩时的表面肌电信号，通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度，结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息；

步骤3中所述采集关节动作最大自主收缩时的表面肌电信号为EMG_m：

步骤3中所述通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度为：

肌肉协同理论将肌肉活动状态表示为肌肉协同元和激活系数的线性组合，得到的肌肉协同矩阵可表示肌肉对该动作的贡献度，如下：

V_N×T≈W_N×K×H_K×T

其中，V_N×T是表面肌电信号数据集矩阵，N是选取的肌肉数量，T是时间样本数，K是协同元数量，W_N×K是具有K个协同元的肌肉协同矩阵，H_K×T为肌肉激活系数矩阵。

脚踝屈伸有背伸和跖屈共Q＝2种动作，对这2种动作采集的最大自主收缩时的表面肌电信号分别提取均方根值(Root Mean Square，RMS)作为肌肉活动水平如下：

其中，x_i,q代表了第q个动作的表面肌电信号序列，n_q表示第q个序列长度。

对脚踝背伸或者跖屈动作协同元数目取1，此时肌肉协同矩阵变成列向量，采用肌肉协同理论计算第i肌肉对第q个动作的贡献度WD_iq如下：

其中，N＝4为肌肉数量，W_i为第i块肌肉在第q个动作的肌肉协同矩阵，并取所有关节动作的贡献度均值作为第i块肌肉的平均贡献度c_i如下：

步骤3中所述结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息为：

运动疲劳的计算选择的特征是中值频率(Median Frequency，MF)，计算如下：

其中，P(f)表示表面肌电信号的功率谱密度，f₁和f₂分别表示信号带宽的最低频率和最高频率。进一步可得运动疲劳信息如下：

其中，MF_i为第i块肌肉的疲劳信息，MF_oi为第i块肌肉的初始疲劳信息，c_i为第i块肌肉在屈伸过程中的平均贡献度。

步骤4：采集康复训练过程中的人机交互力矩，建立阻抗模型，设定康复机器人的刚度参数和阻尼参数的基础值，实现基于位置的阻抗控制模型；

作为优选，步骤4中所述建立阻抗模型为:

式(15)中

和q分别表示实际加速度、速度和轨迹，

和q_d分别表示期望加速度、速度和轨迹，τ表示人机交互力矩，有力矩传感器采集获得， M₀、B₀和K₀分别表示惯性、阻尼和刚度参数的基础值。

采集人机交互力矩τ，通过经验设定刚度参数的基础值即K₀和阻尼参数得基础值即B₀，由于康复机器人的运行速度通常比较匀速缓慢，加速度项通常可以忽略不计，可以得到阻抗模型如下：

进一步转换，步骤4中所述基于位置的阻抗控制模型为：

步骤5：设定康复机器人运行的期望轨迹，将关节刚度信息和运动疲劳信息引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数，得到实际的运行轨迹。

作为优选，步骤5中所述设定康复机器人运行的期望轨迹为：

结合图1，采用饱和函数对刚度和阻尼参数进行约束，随着关节刚度的增大，康复机器人的刚度自适应地变小，；同时随着疲劳的加深，阻尼参数也自适应地减小，引入步骤1中的关节刚度估计模型计算关节刚度信息k，引入步骤3中的运动疲劳信息p，结合饱和函数，将关节刚度信息k和运动疲劳信息p对刚度参数和阻尼参数的基础值进行调节，得到自适应变化的刚度参数K和阻尼参数B 如下：

其中，B₀为阻尼的基础值，K₀为刚度参数的基础值，B_l0为阻尼参数的下界， B_h0为阻尼参数的上界，K_l0为刚度参数的下界，K_h0为刚度参数的上界，p为疲劳信息，k为刚度信息，η为疲劳信息的权重因子，υ为刚度信息的权重因子；

步骤5中所述引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数，得到实际的运行轨迹为：

引入步骤4中基于位置的阻抗控制模型，将自适应变化的刚度参数和阻尼参数替换刚度参数和阻尼参数的基础值，得到自适应阻抗控制模型如下：

将上式

中的机器人运行速度

进行积分操作，可得到实际的运行轨迹q。

本发明中的阻尼参数和刚度参数的基础值、上界和下界，以及疲劳信息和刚度信息的权重因子可根据具体的实验平台和实验人员的经验设定。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (2)

1.一种基于关节刚度和肌肉疲劳的康复机器人控制方法，其特征在于，包括以下步骤：

步骤1：采用运动采集与分析系统对表面肌电信号和关节角度信号进行同步采集，并分别对表面肌电信号和关节角度信号进行预处理，得到预处理之后的表面肌电信号和实际关节角度信号；

步骤2：结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型，将预处理后表面肌电信号EMG_ap和实际关节角度信号θ分别输入关节刚度估计模型和关节逆向动力学模型，利用遗传算法辨识出与受试者相匹配的生理参数，得到个性化关节刚度估计模型；

步骤3：采集关节动作最大自主收缩时的表面肌电信号，通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度，结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息；

步骤4：采集康复训练过程中的人机交互力矩，建立阻抗模型，设定康复机器人的刚度参数和阻尼参数的基础值，实现基于位置的阻抗控制模型；

步骤5：设定康复机器人运行的期望轨迹，将关节刚度信息和运动疲劳信息引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数，得到实际的运行轨迹；

步骤2中所述结合肌肉激活模型、Hill肌肉肌腱模型和骨骼肌肉几何模型建立关节刚度估计模型：

引入步骤1中预处理后表面肌电信号EMG_ap取绝对值并进行归一化处理，得到归一化表面肌电信号e(t)如下：

其中，EMG_rea＝|EMG_ap|为取绝对值之后的表面肌电信号，EMG_mvc为最大自主收缩期间记录的处理过的表面肌电信号峰值，EMG_res为静息状态下的表面肌电信号；

然后使用二阶离散线性模型得到神经激活度u(t)如下：

u(t)＝0.9486·e(t-de)+0.052·u(t-1)-0.000627·u(t-2)

其中，de＝80ms；最后获取肌肉激活度a(t)如下：

引入步骤1中实际关节角度信号θ，通过肌肉路径方程和骨骼肌肉几何模型计算肌肉长度l_mt和肌肉力力臂长度d如下：

β＝σ+θ

d＝l_PRl_QRsinβ/l_PQ

其中，P为肌肉起点，Q为肌肉止点，R为关节转动中心，σ为初始状态肌肉起点和止点与转动中心之间的角度；

由Hill肌肉肌腱模型结合给定的生理参数计算肌肉力F如下：

其中，

是最佳肌肉纤维长度下的最大收缩力，a为肌肉激活度，l为标准化的肌肉纤维长度，v为标准化的肌肉纤维收缩速度，F_A(l)为标准化的主动力-肌纤维长度函数，F_V(v)为标准化的主动力-肌纤维收缩速度函数，F_p(l)为标准化的被动力-肌纤维长度函数，α称为羽状角；

使用肌肉力和肌肉力力臂长度计算出关节力矩τ_a如下：

引入步骤1中实际关节角度信号θ，使用关节力矩对实际关节角度信号求导数计算关节转动的刚度K_a如下：

其中，i为参与计算的第i块肌肉，τ_a为关节力矩；

引入步骤1中实际关节角度信号θ，使用逆向动力学模型计算参考力矩τ_r如下：

其中，g为重力加速度，I、M、l分别为关节转动时转动肢体的转动惯量、质量以及重力力臂长度；

引入步骤1中实际关节角度信号θ，通过参考力矩计算参考刚度K_r如下：

以刚度K_a和参考刚度K_r之间的误差最小化为优化目标，为每块肌肉设置对应的最佳肌肉纤维长度下的最大收缩力

最佳肌肉纤维长度l_mopt和肌腱长度调节因子s_t的初始值及寻优范围，通过遗传算法的选择、交叉和变异操作对最佳肌肉纤维长度下的最大收缩力、最佳肌肉纤维长度和肌腱长度调节因子进行寻优，从而得到优化后参数即受试者个性化生理参数，以构建个性化关节刚度估计模型如下：

其中，

l_mopt和s_t对每块肌肉有与之对应的寻优范围，h为数据点，H为总的数据点数；

步骤3中所述采集关节动作最大自主收缩时的表面肌电信号为EMG_m：

步骤3中所述通过肌肉协同理论分析受试者肌肉在运动过程中的平均贡献度为：

假设有Q种不同关节动作，对这Q种动作采集的最大自主收缩时的表面肌电信号分别提取均方根值RMS作为肌肉活动水平如下：

其中，x_i,q代表了第q个动作的表面肌电信号序列，n_q表示第q个序列长度；

采用肌肉协同理论计算第i块肌肉对第q个动作的贡献度WD_i，q如下：

其中，N为肌肉数量，W_i,q为第i块肌肉在第q个动作的肌肉协同矩阵，并取所有关节动作的贡献度均值作为第i块肌肉的平均贡献度c_i如下：

步骤3中所述结合运动过程中计算的每块肌肉的中值频率，得到运动疲劳信息为：

计算第i块肌肉的中值频率反映肌肉的疲劳信息，并通过第i块肌肉的中值频率反映肌肉的疲劳信息相对变化值结合平均贡献度获取运动疲劳信息p如下：

其中，p_i为第i块肌肉的疲劳信息，p_0i为第i块肌肉的初始疲劳信息；

步骤4中所述建立阻抗模型为:

采集人机交互力矩τ，通过经验设定刚度参数的基础值即K₀和阻尼参数的基础值即B₀，得到阻抗模型如下：

其中，

和q分别表示速度和轨迹，

和q_d分别表示期望速度和轨迹；

步骤4中所述基于位置的阻抗控制模型为：

步骤5中所述将关节刚度信息和运动疲劳信息引入基于位置的阻抗控制模型中调节刚度参数和阻尼参数为：

引入步骤2中的个性化关节刚度估计模型计算关节刚度信息k，引入步骤3中的运动疲劳信息p，结合饱和函数，将关节刚度信息k和运动疲劳信息p对刚度参数和阻尼参数的基础值进行调节，得到自适应变化的刚度参数K和阻尼参数B如下：

其中，B₀为阻尼参数的基础值，K₀为刚度参数的基础值，B_l0为阻尼参数的下界，B_h0为阻尼参数的上界，K_l0为刚度参数的下界，K_h0为刚度参数的上界，p为疲劳信息，k为刚度信息，η为疲劳信息的权重因子，υ为刚度信息的权重因子；

步骤5中所述实际的运行轨迹，具体计算方法为：

引入步骤4中基于位置的阻抗控制模型，将自适应变化的刚度参数和阻尼参数替换刚度参数和阻尼参数的基础值，得到自适应阻抗控制模型如下：

将上式

中的机器人速度

进行积分操作，可得到实际的运行轨迹q。

2.根据权利要求1所述的基于关节刚度和肌肉疲劳的康复机器人控制方法，其特征在于：步骤1中所述采用运动采集与分析系统对表面肌电信号进行同步采集为EMG；

步骤1中所述采用运动采集与分析系统对关节角度信号进行同步采集并进行计算为：

在关节转动时固定肢体的末端、关节转动中心和关节转动时转动肢体的末端分别贴上标记点，对应记为点A(x_A,y_A,z_A)、点B(x_B,y_B,z_B)和点C(x_C,y_C,z_C)，记

为点B和点A构成的向量，

为点B和点C构成的向量；

通过运动采集与分析系统中标记点获取关节角度信号，通过标记点构成的向量计算关节角度信号θ_a：

步骤1中所述对表面肌电信号进行预处理为：

采用阶数为5阶，截止频率为20Hz的巴特沃斯高通滤波器对表面肌电信号进行预处理，以减少基线漂移和伪影噪声的影响，并得到预处理后表面肌电信号EMG_ap；

步骤1中所述对关节角度信号进行预处理为：

θ＝θ_a-θ₀

其中，θ为实际关节角度信号，θ₀为标记点的初始角度误差。

Images