CN105615890A - 人体下肢行走关节角度肌电连续解码方法 - Google Patents

Abstract

本发明公开一种人体下肢行走关节角度肌电连续解码方法，利用光学运动捕捉系统记录人体行走过程中下肢体表光学标记点的运动轨迹，通过人体下肢运动学建模精确计算下肢关节运动角度；同步采集人体行走过程中与下肢运动相关的八块主力肌肉的表面肌电信号，通过滤波、整流预处理提取其活动强度信息，通过主成分分析方法提取描述表面肌电信号强度的最优独立特征向量集；采用基因表达式编程符号回归分析方法，构建从表面肌电信号特征(自变量)到矢状面关节运动角度(因变量)的非线性回归模型，预测下肢运动轨迹。本发明主要应用于医疗康复机械的设计制造。

Description

人体下肢行走关节角度肌电连续解码方法

技术领域

本发明涉及医疗康复器械，具体讲，涉及人体下肢行走关节角度肌电连续解码方法。

背景技术

当今社会，由于车祸、工程意外、或脑中风等疾病引起的骨折、肢体运动功能缺失及其并发症的发病率有逐年上升的趋势。其中，下肢运动功能的缺失对患者日常生活的影响尤为严重。实验表明，特定的功能训练能促进中枢神经的重组和代偿，恢复患者的肢体运动功能，下肢康复机器人由此发明和发展。

人体的表面肌电信号(sEMG)是一种生物电信号，能够直接反映人体的运动意图，非常适合作为康复机器人机器人的控制信号来源，即人机交互系统的有效的接口。前期通过sEMG识别人体行为意图的研究大多集中在对人体动作分类，此类方法只能预测少数离散肢体动作，应用预测结果控制机器人，无法实现机器人关节像人体关节一样连续自如运动，而保证机器人运动与人体运动高度匹配是实现多类机器人系统安全控制的先决条件，因此通过sEMG估计连续的人体关节运动量是当前理论与应用研究的关键点。

肌电解码预测角度方面有关专利：

1.发明专利。基于相关系数的肌电假肢控制源导联优化方法，201210192585.2。

专利1利用人工神经网络建立关联sEMG和关节运动量的回归模型，神经网络的结构形式对预测结果影响很大，却没有规范方法定义合适的网络结构，并且网络参数确定过度依赖于训练数据，当测试数据和训练数据有差异时，易造成网络输出偏离真实值。

肌电解码预测角度方面有关论文：

1.BuchananTS,LloydDG,ManalK,etal.Neuromusculoskeletalmodeling：Estimationofmuscleforcesandjointmomentsandmovementsfrommeasurementsofneuralcommand[J].JournalofAppliedBiomech,2004,20(4):367-395

2.FleischerC,HommelG.Ahuman-exoskeletoninterfaceutilizingelectromyography[J].IEEETransactionsonRobotics,2008,24(4):872-882

文献1和文献2结合肌肉生理学建立以sEMG为输入的关节动力学模型，进而计算关节力矩/角加速度等量，但是构建的模型结构复杂，含有许多无法测量的生理学参数，难以实际应用。

发明内容

本发明针对现有技术存在的上述不足，提供一种人体下肢行走关节角度肌电连续解码方法，该方法可通过采集患者下肢表面肌电信号来实时准确稳定预测下肢关节角度，患者可通过穿戴下肢康复机器人实现下肢运动功能康复的主动训练，从而增加康复的积极性，逐渐拥有自主下肢运动的能力。

本发明是通过以下技术方案实现的：

一种人体下肢行走关节角度肌电连续解码方法，包括与细以下步骤：

步骤一、利用光学运动捕捉系统记录人体下肢行走运动三维坐标，通过人体下肢运动学建模方法计算下肢髋关节、膝关节和踝关节的运动角度；

步骤二、通过16路肌电仪同步采集运动期间人体下肢运动相关8块肌肉的表面肌电信号，8块肌肉具体为：股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧和胫骨前肌；

步骤三、对采集的8块肌肉的表面肌电信号进行滤波、整流预处理提取其活动强度信息，并使用PCA方法对肌电强度信息降维，提取肌电信号活跃度主分量；

步骤四、以步骤三选择所得的肌电信号活跃度主分量作为输入量和以步骤一计算所得关节运动角度作为输出量构建训练集，使用GEP符号回归分析方法，建立非线性下肢关节角度预测模型；

步骤五、采集人体下肢行走运动相关8块肌肉的表面肌电信号，8块肌肉具体为：股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧和胫骨前肌；对采集的8块肌肉的表面肌电信号进行滤波、整流预处理提取其活动强度信息，并使用PCA方法对肌电强度信息降维，提取肌电信号活跃度主分量；将提取的肌电信号活跃度主分量输入所述非线性下肢关节角度预测模型中，获得人体下肢髋关节、膝关节和踝关节的运动角度预测值。

进一步的，步骤一中计算关节运动角度具体包括：

所求角度均为各个关节在人体矢状面上的运动角度，其计算公式如下：

θ₁＝acos(T₁k_Thigh·k_Pelvis)

θ₂＝acos(T₂k_Calf·k_Thigh)

θ₃＝acos(T₃k_Foot·k_Calf)

式中：θ₁、θ₂和θ₃分别为髋关节、膝关节和踝关节的屈/伸角；T₁、T₂和T₃分别为从大腿坐标系到骨盆坐标系、小腿坐标系到大腿坐标系和足部坐标系到小腿坐标系的转化矩阵；k_Pelvis、k_Thigh、k_Calf和k_foot分别为骨盆坐标系、大腿坐标系、小腿坐标系和足部坐标系竖直轴上的单位矢量。

进一步的，步骤三具体包括：

利用20-500Hz带通滤波器和50Hz工频陷波器对肌电信号进行降噪和伪迹消除处理，之后对信号进行全波整流和归一化处理，得到肌肉活跃强度X＝{x_ij}_m×p作为PCA降维的输入；选择贡献率之和大于给定提取效率的前h个肌肉活跃度主分量α＝[α₁,α₂,…,α_h]作为GEP算法条件属性，其h取值为2～6；p＝8表示表面肌电信号通道数，m表示采样点数，给定提取效率为90％。

进一步的，步骤四中GEP符号回归方法具体包括：

用GEP算法的多基因染色体表示问题的可行解，包括可行解的形式和参数，然后利用遗传进化在解空间中进行搜索，最后根据适应度值进行优胜劣汰，求得最优解，即拟合样本肌肉活跃度主分量与关节角度最好的非线性预测模型。

进一步的，所述遗传进化具体为：

设定函数集{+，-，*，/，Q}和终结符集{α₁α₂…α_h？}，其中Q表示平方运算，α表示肌电信号活跃度主分量，？代表常数量；定义3段式基因结构，包括头部、尾部和常数域，头部元素来自函数集或终极符集，尾部元素来自终结符集，头部长度设为9，尾部长度由固定关系计算得为10，常数域长度与尾部长度相等，染色体由6个基因随机产生；初始化生成500组未知染色体形成初始子种群，即种群规模为500，进化代数计数器的初始值设为1，最大遗传代数设为200；

设定动态适应度函数

F_i'＝F_avg(F_i-F_min)/(F_avg-F_min)

$E_i=\underset{j=1}{\overset{n}{\Σ}}{(P_{ij}-T_j)}^2$

$F_i=1000\×\frac{1}{1+E_i}$

式中：E_i为预测值与目标值均方误差平方；P_ij表示输入第i个染色体第j组样本特征时染色体表达式值即关节角度预测值；T_j为对应的关节角度目标值，即利用光学运动捕捉系统记录人体下肢行走运动三维坐标所计算的关节角度；F_i为以均误差平方作为染色体适应度评价方法的评价函数；F_avg是个体平均适应度；F_min是个体最小适应度；

个体选择采用二次选择策略，一次选择排除由于除法运算中出现除数为零的无效基因，二次选择策略为轮盘赌法，同时选择部分最优个体3组直接进入下一代。

进一步的，GEP算法中突变概率为0.051，单点重组率为0.3，双点重组概率为0.3，基因重组率为0.1，插串迁移概率为0.1，根插串迁移概率为0.1，基因迁移概率为0.1。

进一步的，利用遗传进化在解空间中进行搜索的终止条件设定为达到最大迭代次数200次或者连续20代最好适应值未变化。

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

本发明提取下肢8块肌肉肌电信号的肌肉活跃度主分量，结合运动学建模所得下肢关节运动角度，使用GEP符号回归分析方法，建立非线性回归模型，预测下肢运动轨迹。本发明实现人体下肢步态动作的连续精细识别，建立表面肌电信号与下肢关节运动角度之间的回归模型，达到控制下肢康复机器人主动训练的目的。本发明所建模型稳定可靠，预测角度实施准确，在康复工程、人机接口、虚拟世界等方面都有潜在的应用价值。

附图说明

图1是人体下肢行走关节角度肌电连续解码方法框图；

图2是下肢标记点位置和坐标系布置示意图；其中，1为左髂前上棘；2为右髂前上棘；3为左髂后上棘；4为右髂后上棘；5为左膝外侧；6为左膝内侧；7为左踝外侧；8为左踝内侧；9为右膝外侧；10为右膝内侧；11为右踝外侧；12为右踝内侧；13为左足跖骨；14为左脚后跟；15为右脚趾骨；16为右脚后跟；

图3是关节运动角度计算结果示意图；

图4是表面肌电信号活跃度主分量示意图；

图5是GEP算法基本流程图；

图6是关节运动角度预测结果示意图。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于了解，下面结合具体图示，进一步阐述本发明技术方案的具体实现方式。

请参阅图1至图6所示，本发明一种人体下肢行走关节角度肌电连续解码方法，包括下列步骤：

采用红外三维运动捕捉系统记录人体运动三维坐标，进而计算下肢髋关节、膝关节和踝关节的运动角度；利用16通道肌电采集仪采集人体在行走过程中8块肌肉的肌电信号，他们是：股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧、胫骨前肌，经过预处理和降维后提取肌电信号活跃度主分量作为特征参数；随后，利用基因表达式编程符号回归分析方法(GEP)，构建从表面肌电信号特征到矢状面关节运动角度的非线性回归模型。人体下肢行走关节角度肌电连续解码方法的示意图如图1。结果发现，本发明的方法模型可靠，角度预测正确率高，有利于下肢康复机器人的训练和使用。

1、关节角度计算与肌电特征提取的方法：

光学运动捕捉系统记录到人体下肢行走过程中16个体表关键标记点在世界坐标系下的三维运动轨迹，其下肢标记点的位置和坐标系布置如图2所示。其世界坐标系X轴方向指向人体正前方，Z轴竖直向上，两者构成平面为矢状面，而X轴与Y轴构成水平面，Y轴与Z轴构成额状面。局部参考坐标系以小腿坐标系为例，首先，基于标记点的位置坐标确定膝和踝关节中心点的位置坐标，对于膝关节和踝关节，以关节内外侧标记点的中点作为关节中心，即

$p_{kjc}=\frac{p_5+p_6}{2}$

$p_{ajc}=\frac{p_7+p_8}{2}$

式中：p₅、p₆、p₇和p₈分别为图2中下肢标记点5、6、7和8点的位置坐标；p_kjc为膝关节中心位置坐标；p_ajc为踝关节中心位置坐标。其次，利用关节中心位置坐标和标记点位置坐标确定小腿的局部参考坐标系，规定坐标系位于小腿质心位置，k轴由p_ajc指向p_kjc，i轴垂直于p_ajc、p_kjc和p₆三点组成的平面向前，j轴垂直于k轴和i轴，小腿坐标系单位矢量为

$k_{Calf}=\frac{p_{kjc}-p_{ajc}}{\left|p_{kjc}-p_{ajc}\right|}$

$i_{Calf}=\frac{(p_{ajc}-p_{kjc})\×(p_6-p_{hjc})}{|(p_{ajc}-p_{kjc})\×(p_6-p_{hjc})|}$

j_Calf＝k_Calf×i_Calf

式中：k_Calf、i_Calf和j_Calf分别小腿坐标系k轴、i轴和j轴的单位矢量。

其骨盆坐标系、大腿坐标系和足坐标系建立方法与小腿坐标系一致。针对建立的局部坐标系，定义3个关节的屈伸角度：髋关节屈/伸角θ₁为k_Thigh在骨盆坐标系矢状面的投影与竖直轴的夹角；膝关节屈/伸角θ₂为k_Calf在大腿坐标系矢状面的投影与竖直轴的夹角；踝关节屈/伸角θ₃为k_foot在大腿坐标系矢状面的投影与竖直轴的夹角。计算公式如下

θ₁＝acos(T₁k_Thigh·k_Pelvis)

θ₂＝acos(T₂k_Calf·k_Thigh)

θ₃＝acos(T₃k_Foot·k_Calf)

式中：T₁、T₂和T₃分别为从大腿坐标系到骨盆坐标系、小腿坐标系到大腿坐标系和足部坐标系到小腿坐标系的转化矩阵。k_Pelvis、k_Thigh和k_foot分别为骨盆坐标系、大腿坐标系和足部坐标系k轴上的单位矢量。图3给出某次行走过程中计算的关节运动角度。

采集的肌电信号利用带通滤波器(20-500Hz)和50Hz工频陷波器进行降噪和伪迹消除处理，并对信号进行全波整流和归一化处理，得到肌肉活跃强度X＝{x_ij}_m×p(p＝8表示表面肌电信号通道数，m表示采样点数)作为主成分分析方法(PCA)降维的输入。

其利用PCA实现肌电活跃主分量提取的过程为

(1)计算肌电信号强度指标数据的相关系数矩阵，R＝(r_ij)_p×p

$r_{ij}=\frac{\mathrm{cov}(x_i,x_j)}{\sqrt{D\left(x_i\right)}\sqrt{D\left(x_j\right)}}$

式中：cov(x_i,x_j)为强度指标x_i和x_j的协方差；D(x_i)和D(x_j)分别为x_i和x_j的方差。

(2)计算相关系数矩阵R的特征值λ以及对应的特征向量，同时对特征向量进行正交归一化得到单位正交化特征向量A；

(3)计算每一个特征值的贡献率，根据给定的提取效率(>90％)，选定h个肌电信号活跃度主分量α＝[α₁,α₂,…,α_h]，其h取值2～6。图4给出某次行走过程中处理得到的肌肉活跃度主分量。

2、GEP符号回归方法建模原理

基因表达式程序设计(Geneticexpressionprogramming，GEP)是由葡萄牙学者Ferreira承袭遗传算法与遗传规划算法而新发展出来的一种进化算法，GEP利用适应度函数对于种群个体的选择加上演化的机制，将程式自行演化以达成特定的演化目标，适用于解决各种不同的领域的问题。在功能上，GEP能发现揭示问题本质的规则、公式、以及描述问题解答过程的程序，通过简单紧凑的编码和简单的遗传操作，能够解决复杂应用问题。从形式上，GEP是一种进化计算方法，采用了个体表现型和遗传编码不同的设计思路，采用等长线性符号编码，遗传操作非常简单，基本类似于标准遗传算法。GEP符号回归方法的流程图如图5所示，其具体步骤如下：

步骤一、产生初始群体，解码并计算适应度；

步骤二、执行精英策略，判断结果是否符合停机条件(进化是否达到一定代数或者连续n代最好适应值未变化，n优选取20)，若满足，则输出相应的表达式并结束运算；

步骤三、按某一策略执行选择操作，形成交配池；

步骤四、种群进化(分别以某一概率执行变异、插入、迁移、重组等遗传算子)；

步骤五、产生新种群，解码并计算适应度，跳转至步骤二重复执行上面步骤。

所述步骤一具体包括：

步骤11、初始种群创建中，设定函数集{+，-，*，/，Q}和终结符集{α₁α₂…α_h？}，其中Q表示平方运算，α表示肌电信号活跃度主分量，？代表常数量；定义3段式基因结构，包括头部、尾部和常数域，头部元素来自函数集或终极符集，尾部元素来自终结符集，常数域长度与尾部长度相等，头部与尾部长度满足以下关系

Tl＝Hl×(N-1)+1

式中：Hl为基因头部长度(事先确定的值)；Tl为基因尾部长度；N表示函数符号集中参数个数；

本发明预测关节角度的数据对象是多元时间序列数据，染色体构成采用多基因染色体的方式。初始化生成M组未知染色体形成初始子种群，即种群规模为M，进化代数计数器Gen的初始值设为1，最大遗传代数为MaxGen；

步骤12、设定动态适应度函数

F_i'＝F_avg(F_i-F_min)/(F_avg-F_min)

$E_i=\underset{j=1}{\overset{n}{\Σ}}{(P_{ij}-T_j)}^2$

$F_i=1000\×\frac{1}{1+E_i}$

式中：E_i为预测值与目标值均方误差平方；P_ij表示输入第i个染色体第j组样本特征时染色体表达式值即关节角度预测值；T_j为对应的关节角度目标值，即利用光学运动捕捉系统记录人体下肢行走运动三维坐标所计算的关节角度；F_i为以均误差平方作为染色体适应度评价方法的评价函数；F_avg是个体平均适应度；F_min是个体最小适应度。

所述步骤二停机条件中适应度的收敛精度即|P_ij-T_j|采用多次残差回归模型修正的控制方法，即利用残差数据建立残差回归模型，以此修正回归模型，提高精度。

所述步骤三个体选择采用二次选择策略，一次选择排除由于除法运算中出现除数为零的无效基因二次选择策略为轮盘赌法，同时选择部分最优个体m组直接进入下一代。

所述步骤四种群进化所使用的遗传算子有选择算子、变异算子、变换和插入序列元素以及重组算子。

3、预测结果论证

截取3种不同行走速度下完整步态周期数据，组成样本数据，并进行关节角度计算和表面肌电信号活跃强度主分量提取，以计算关节角度为决策属性(因变量)，以主分量构成独立特征向量为条件属性(自变量)，编写GEP算法程序，建立函数关系模型。其中GEP算法参数设置如表1所示

表1GEP参数设置

通过肌电连续解码，使用GEP符号回归方法预测关节运动角度的预测结果如图6所示，从图6可以看出，该方法能够有效稳定、实时精确的预测下肢关节运动角度，可用于下肢康复机器人的主动训练控制。

本发明是一种人体下肢行走关节角度肌电连续解码方法，以实现下肢康复机器人在实际使用过程中的精确控制。该方法提取股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧、胫骨前肌八块肌肉行走过程中的肌电信号，结合下肢关节运动角度，利用GEP符号回归方法建立非线性回归模型，达到精确控制康复机器人的目的。该发明所建模型可靠，角度预测正确率高，有利于下肢康复机器人的训练和使用，具有可观的社会效益和经济效益。

Claims (7)

1.一种人体下肢行走关节角度肌电连续解码方法，其特征在于，包括以下步骤：

步骤一、利用光学运动捕捉系统记录人体下肢行走运动三维坐标，通过人体下肢运动学建模方法计算下肢髋关节、膝关节和踝关节的运动角度；

步骤二、通过16路肌电仪同步采集运动期间人体下肢运动相关8块肌肉的表面肌电信号，8块肌肉具体为：股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧和胫骨前肌；

步骤三、对采集的8块肌肉的表面肌电信号进行滤波、整流预处理提取其活动强度信息，并使用PCA方法对肌电强度信息降维，提取肌电信号活跃度主分量；

步骤四、以步骤三选择所得的肌电信号活跃度主分量作为输入量和以步骤一计算所得关节运动角度作为输出量构建训练集，使用GEP符号回归分析方法，建立非线性下肢关节角度预测模型；

步骤五、采集人体下肢行走运动相关8块肌肉的表面肌电信号，8块肌肉具体为：股直肌、股外侧肌、股二头肌、半腱肌、比目鱼肌、腓肠肌内侧、腓肠肌外侧和胫骨前肌；对采集的8块肌肉的表面肌电信号进行滤波、整流预处理提取其活动强度信息，并使用PCA方法对肌电强度信息降维，提取肌电信号活跃度主分量；将提取的肌电信号活跃度主分量输入所述非线性下肢关节角度预测模型中，获得人体下肢髋关节、膝关节和踝关节的运动角度预测值。

2.据权利要求1所述的人体下肢行走关节角度肌电连续解码方法，其特征在于，步骤一中计算关节运动角度具体包括：

所求角度均为各个关节在人体矢状面上的运动角度，其计算公式如下：

θ₁＝acos(T₁k_Thigh·k_Pelvis)

θ₂＝acos(T₂k_Calf·k_Thigh)

θ₃＝acos(T₃k_Foot·k_Calf)

式中：θ₁、θ₂和θ₃分别为髋关节、膝关节和踝关节的屈/伸角；T₁、T₂和T₃分别为从大腿坐标系到骨盆坐标系、小腿坐标系到大腿坐标系和足部坐标系到小腿坐标系的转化矩阵；k_Pelvis、k_Thigh、k_Calf和k_foot分别为骨盆坐标系、大腿坐标系、小腿坐标系和足部坐标系竖直轴上的单位矢量。

3.据权利要求1所述的人体下肢行走关节角度肌电连续解码方法，其特征在于，步骤三具体包括：

利用20-500Hz带通滤波器和50Hz工频陷波器对肌电信号进行降噪和伪迹消除处理，之后对信号进行全波整流和归一化处理，得到肌肉活跃强度X＝{x_ij}_m×p作为PCA降维的输入；选择贡献率之和大于给定提取效率的前h个肌肉活跃度主分量α＝[α₁,α₂,…,α_h]作为GEP算法条件属性，其h取值为2～6；p＝8表示表面肌电信号通道数，m表示采样点数，给定提取效率为90％。

4.据权利要求1所述的人体下肢行走关节角度肌电连续解码方法，其特征在于，步骤四中GEP符号回归方法具体包括：

用GEP算法的多基因染色体表示问题的可行解，包括可行解的形式和参数，然后利用遗传进化在解空间中进行搜索，最后根据适应度值进行优胜劣汰，求得最优解，即拟合样本肌肉活跃度主分量与关节角度最好的非线性预测模型。

5.据权利要求4所述的人体下肢行走关节角度肌电连续解码方法，其特征在于，所述遗传进化具体为：

设定函数集{+，-，*，/，Q}和终结符集{α₁α₂…α_h？}，其中Q表示平方运算，α表示肌电信号活跃度主分量，？代表常数量；定义3段式基因结构，包括头部、尾部和常数域，头部元素来自函数集或终极符集，尾部元素来自终结符集，头部长度设为9，尾部长度由固定关系计算得为10，常数域长度与尾部长度相等，染色体由6个基因随机产生；初始化生成500组未知染色体形成初始子种群，即种群规模为500，进化代数计数器的初始值设为1，最大遗传代数设为200；

设定动态适应度函数

F′_i＝F_avg(F_i-F_min)/(F_avg-F_min)

$E_i=\underset{j=1}{\overset{n}{\Σ}}{(P_{ij}-T_j)}^2$

$F_i=1000\×\frac{1}{1+E_i}$

式中：E_i为预测值与目标值均方误差平方；P_ij表示输入第i个染色体第j组样本特征时染色体表达式值即关节角度预测值；T_j为对应的关节角度目标值，即利用光学运动捕捉系统记录人体下肢行走运动三维坐标所计算的关节角度；F_i为以均误差平方作为染色体适应度评价方法的评价函数；F_avg是个体平均适应度；F_min是个体最小适应度；

个体选择采用二次选择策略，一次选择排除由于除法运算中出现除数为零的无效基因，二次选择策略为轮盘赌法，同时选择部分最优个体3组直接进入下一代。

6.据权利要求4所述的人体下肢行走关节角度肌电连续解码方法，其特征在于：

GEP算法中突变概率为0.051，单点重组率为0.3，双点重组概率为0.3，基因重组率为0.1，插串迁移概率为0.1，根插串迁移概率为0.1，基因迁移概率为0.1。

7.根据权利要求4所述的人体下肢行走关节角度肌电连续解码方法，其特征在于：

利用遗传进化在解空间中进行搜索的终止条件设定为达到最大迭代次数200次或者连续20代最好适应值未变化。