CN113576463B - 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统 - Google Patents

肌电信号驱动的膝关节肌骨模型接触力估计方法及系统 Download PDF

Info

Publication number
CN113576463B
CN113576463B CN202110877393.4A CN202110877393A CN113576463B CN 113576463 B CN113576463 B CN 113576463B CN 202110877393 A CN202110877393 A CN 202110877393A CN 113576463 B CN113576463 B CN 113576463B
Authority
CN
China
Prior art keywords
muscle
model
mtu
force
ith
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110877393.4A
Other languages
English (en)
Other versions
CN113576463A (zh
Inventor
李玉榕
张文萱
苏志超
连章汇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fuzhou University
Original Assignee
Fuzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fuzhou University filed Critical Fuzhou University
Priority to CN202110877393.4A priority Critical patent/CN113576463B/zh
Publication of CN113576463A publication Critical patent/CN113576463A/zh
Application granted granted Critical
Publication of CN113576463B publication Critical patent/CN113576463B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/112—Gait analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/1036—Measuring load distribution, e.g. podologic studies
    • A61B5/1038—Measuring plantar pressure during gait
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/1121—Determining geometric values, e.g. centre of rotation or angular range of movement
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/389—Electromyography [EMG]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/389—Electromyography [EMG]
    • A61B5/397—Analysis of electromyograms
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/45—For evaluating or diagnosing the musculoskeletal system or teeth
    • A61B5/4538—Evaluating a particular part of the muscoloskeletal system or a particular medical condition
    • A61B5/4585—Evaluating the knee
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters

Abstract

本发明涉及一种肌电信号驱动的膝关节肌骨模型接触力估计方法及系统，该方法包括：建立接触力估计模型，包括信号处理模型、OpenSim个体模型、肌肉激活度模型、肌肉收缩力学模型以及外部内收力矩模型，信号处理模型对输入数据进行预处理，OpenSim个体模型将通用模型缩放至个体生理参数，肌肉激活度模型用于获取肌肉激活度，肌肉收缩力学模型用于获取肌肉力矩，外部内收力矩模型用于获取外部内收力矩；对接触力估计模型进行参数优化；获取肌电信号数据、关节角度信号数据和足底压力信号数据，将其作为输入信号输入到参数优化后的接触力估计模型中，计算得到膝关节接触力。该方法及系统使用方便，且有利于提高接触力估计的准确性。

Description

肌电信号驱动的膝关节肌骨模型接触力估计方法及系统
技术领域
本发明属于步态分析领域，具体涉及一种肌电信号驱动的膝关节肌骨模型接触力估计方法及系统。
背景技术
膝关节是人类下肢活动不可避免的一部分，膝关节骨性关节炎(Osteoarthritis,OA)作为一种常见的慢性关节疾病，尤其在中老年群体中尤为常见，在我国患病人数已经超过3600万，50岁以上患病率以10年为单位成倍增长，具有很高的致残率，进行全膝关节置换手术(Total knee arthroplasty,TKA)的患者数量以10万人/年的速度增长。膝关节骨性关节炎的预防与康复与人们的生活质量息息相关，其发展前景决定了能有多少患者重新回归原有的生活状态。
膝骨关节炎致病因素包括生物力学因素、生理因素等，异常的膝关节接触力是膝骨关节炎的发生和进展的主要机械危险因素，且伴随着整个膝关节骨性关节炎的发展过程。根据现有研究表明，膝关节内部高负载和低负载都有可能会造成膝关节损伤。就算进行了TKA手术，若仍旧长期维持不正确的步态，膝关节生物力学的改变对于关节假体会造成不同程度的磨损，进而造成病症的反复或加重，假体的更换也带来了经济负担。因此，估计体内膝关节接触力对膝关节的假体设计、TKA术后运动功能评价及康复治疗都具有重要意义。
正常步态需要中枢神经系统、周围神经系统及骨骼肌肉的动态整合，保证适宜的关节运动和适宜的力量，否则就会出现异常步态。动态肌电图是步态分析非常重要的组成部分，用于检测步行时肌肉活动与步态的关系。由于病理步态分析的发展，临床对于明确步态障碍关键肌肉的需求日益提高，动态肌电图的价值越来越突出。肌肉力量的强弱、关节的灵活性和稳定性以及步态的调整和反应速度等是人为可控因素中最有效降低关节运动损伤的因素。绝大部分运动者有能力复制出正确步态模式应有的力学机制，转言之，膝关节接触力也可以作为步态合理性判断的重要指标。
现有的接触力间接估计方法主要包括有限元建模、机理建模、黑箱建模以及肌骨模型四种，其中，有限元建模仅考虑膝关节局部力学环境，忽略整个下肢解剖结构及运动学影响，例如肌肉对关节力学和运动学的贡献；使用机理建模进行接触力，计算复杂，涉及众多未知生理参数，就算传感器得到的关节角度和关节力矩相同，每个个体也拥有独特的神经解决方案；应用黑箱模型时，用户无法清晰认识到自身的运动方式，即肌肉激活程度及肌力负荷分配方式。肌骨模型可以同时预测和评估关节力学负荷和肌肉作用，因此，基于此种方式进行膝关节接触力的估计是合理的。
步态分析可以提供关于人体运动的动力学和运动学的重要信息，已经成为人体运动科学分析工程师和医疗专业人员评估被试的力量和运动特性的重要研究方法，诸如运动捕捉系统，测力平台或压力传感器和肌电等工具对于进行步态分析是必要的。
现有技术成果多集中于使用黑箱模型研究步频、步长、足姿态参数以及落地模式等外部参数对膝关节接触力的影响，或利用AnyBody等商用生物力学仿真软件或OpenSim等开源肌骨建模软件，建立个体化模型以仿真人体运动。人体的运动学分析及动力学分析主要基于专用的步态分析实验室。通过定点摄像头及运动捕捉系统记录并计算被试对象的运动轨迹，利用三维测力板测量步行时垂直、水平和侧向作用力，将运动学数据和人体骨骼肌肉模型结合，通过多体动力学算法计算到骨骼力、关节力与力矩、肌肉力等生物力学数值。
现有技术中，为了评估每个被试施加的地面反作用力(Ground reaction force,GRF)，最常见的方法有两种。
一是在步态分析中使用至少一个测力平台(板)，以收集其动力学信息——六轴的力和力矩。这些测力板通常与地板齐平地安放，不会影响到被试者正常的走路或跑步，使得被试踏入板上时不会尝试改变步态。但是，不能保证踩踏在测力板上脚步的合理性是步态分析中公认的弱点。为了保证收集的信息可用，被试在测试过程中同时只能有一只脚踩在同一块测力板上，如果受试者错过了测力板，只能将脚的一部分踩在测力板上，或者将两只脚均触到同一个测力板上，则该次试验应被忽略，直到收集到一系列成功的试验，这是极其繁琐的过程。在测试很难维持正常步态的被试时，这个问题特别突出，例如对于交叉步态的患者，几乎不可能有单脚踩踏测力板的情况出现。
针对这个问题现有的几个解决方案均存在一定的缺点，一种常用方法是通过指示受试者从具体位置特定脚开始，令受试者瞄准力板。使用这种方法，受试者需要了解测力平台可能在哪里，并调整步态以获得成功的试验，这将会导致不准确的步态数据。为了防止这种情况，另一种方法是从个人选择的位置开始。然而，Oggero等人审查了一系列步态分析案例，受试者自行选择其起始位置时，即使是60cm的测力平板，只有25％的被试需要3次以下的机会才能成功，多达42％的患者并不能保证试验成功。整个实验过程对被试来说容易产生疲劳，更难获得良好的数据。
对于两只脚同时踩在同一块测力板上时如何分析数据这一问题，通过测力跑台下装有并排的测力板，可以避免这个问题，即使在双脚同时触地期间，每个测力板也分别记录每只脚的地面反作用力。但使用测力跑台的一个问题是如何在跑台上适应传统地面步行的步态，这在老年人群和残疾人士中尤为突出，他们可能没有任何跑步机使用经验。当个体刚开始在跑步机上行走时，会缩短步长和改变正常节奏，大多数受试者通常需要1-3min时间才能适应机器并且做到目测类似于地面行走步态。
综上，已公开发表或公开使用的技术存在的问题和缺陷如下：
1、数据方面：现有的膝关节内部接触力间接估计主要基于肌骨模型进行仿真，输入数据为：通过运动捕捉系统获得的运动轨迹、通过测力平台得到的三位力/力矩以及通过肌电传感器测得的肌电信号，前两项输入数据依赖于高端步态实验室，受场地大小限制，在一般性的研究场所或医疗场景中难以满足。
2、校准方面：公开通用肌骨模型由多个人体数据的平均值进行初始化，对于每块肌肉的最大等长力、肌腱松弛长度、最佳肌纤维长度以及羽状角4个肌肉参数均使用经验值，不能很好地匹配到用户个体，进而会对膝关节接触力的估计造成误差，影响对假体磨损程度的判断。
3、操作方面：现有利用OpenSim此类开源软件进行膝关节接触力仿真估计的案例中，多使用官方提供的GUI界面进行模型配置，进而实现离线操作。使用此种开发者界面进行膝关节接触力估计仅适用于具有一定专业知识的科研人群，对于其他一般性受众使用较为困难，无法实现便捷操作。
发明内容
本发明的目的在于提供一种肌电信号驱动的膝关节肌骨模型接触力估计方法及系统，该方法及系统使用方便，且有利于提高接触力估计的准确性。
为实现上述目的，本发明采用的技术方案是：一种肌电信号驱动的膝关节肌骨模型接触力估计方法，包括以下步骤：
建立接触力估计模型，包括信号处理模型、OpenSim个体模型、肌肉激活度模型、肌肉收缩力学模型以及外部内收力矩模型，所述信号处理模型包括肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块、步态周期划分模块以及时间归一化处理模块，所述肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块分别用于对肌电信号数据、关节角度信号数据和足底压力信号数据进行预处理，所述步态周期划分模块用于进行步态周期划分，所述时间归一化处理模块用于进行时间归一化处理，以统一步态周期长度；所述OpenSim个体模型用于将模型库中的通用模型缩放至个体特定的生理参数；所述肌肉激活度模型用于通过肌电信号获得肌肉激活度；所述肌肉收缩力学模型用于获取肌肉力矩；所述外部内收力矩模型用于获取外部内收力矩；
对接触力估计模型进行参数优化；
获取肌电信号数据、关节角度信号数据和足底压力信号数据，将其作为输入信号输入到参数优化后的接触力估计模型中，计算得到膝关节接触力。
进一步地，通过无线表面肌电传感器获取原始肌电信号数据，通过角度传感器获取原始关节角度信号数据，通过放置于足底的压力传感器获取原始足底压力信号数据。
进一步地，所述肌电信号数据预处理模块按如下方法对肌电信号数据进行预处理：首先，采用4阶巴特沃斯零相移带通滤波器，去除直流分量；并且，进行全波整流将信号取正，以便于后期数据处理；而后，为模拟肌肉特性，进行4阶巴特沃斯零相移低通滤波；由此产生的线性包络相对于从整个试验记录得到的滤波处理过的峰值肌电信号值进行归一化；最后，重采样至与另外两种信号相同的指定频率；
所述关节角度数据预处理模块按如下方法对关节角度信号数据进行预处理：对原始数据进行平均滤波，以消除在数据采集时产生的抖动；最后，重采样至与另外两种信号相同的指定频率；
所述足底压力数据预处理模块按如下方法对足底压力信号数据进行预处理：提取压力传感器测得的力在额状面上垂直于地面的分量作为个体足底压力数据；对原始数据进行平均滤波，以消除在数据采集时产生的抖动；最后，重采样至与另外两种信号相同的指定频率。
进一步地，按如下方法建立肌肉激活度模型：
首先建立表面电极募集到肌肉的动作电位总和，即表面肌电信号sEMG与反映肌肉状态的神经激活度之间的关系，利用归一化处理后的肌电信号求解神经激活度，然后再建立神经激活度与肌肉激活度之间的非线性模型，求解得到肌肉激活度；
使用临界阻尼的线性离散二阶差分系统来描述sEMG和神经激活度之间的关系：
ui(t)＝αei(t-d)-β1ui(t-1)-β2ui(t-2)
其中，ui(t)、ui(t-1)、ui(t-2)分别是第i块肌肉在第t、t-1、t-2个时刻经过递归滤波器计算后的神经激活度，ei(t-d)表示第i块肌肉在第t-d个时刻经过滤波整流后的肌电信号，i＝1,…,n，d表示机电延迟，α是肌肉的增益系数，β1和β2是肌肉的递归系数；为了使方程的稳定解和肌肉激活度位于[0,1]之间，α、β1、β2的约束为：
β1＝C1+C2、β2＝C1×C2且α-β1-β2＝1，|C1|＜1、|C2|＜1
通过非线性关系建立肌肉激活和神经激活度之间的关系：
Figure BDA0003190823250000051
其中，ai(t)为第i块肌肉在第t个时刻肌肉激活度，Ai为第i块肌肉的非线性因子，取值范围为(-3,0)。
进一步地，所述肌肉收缩力学模型包括肌纤维长度拟合模块、肌纤维速度求解模块、归一化处理模块、肌肉力臂计算模块、建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块、肌肉力计算模块以及肌肉力矩计算模块。
进一步地，所述肌纤维长度拟合模块采用膝关节角度进行三次多项式拟合得到肌肉肌腱单元MTU的肌纤维长度；在t时刻，对于跨越一个关节的第i个MTU，将与膝关节在额状面的内收-外展运动、在矢状面的屈曲-伸展运动、在横截面的内旋-外旋运动相关的MTU长度分别拟合为三次多项式：
Figure BDA0003190823250000052
Figure BDA0003190823250000053
Figure BDA0003190823250000054
其中，θi/frontal(t)、θi/sagittal(t)、θi/cross(t)为第i个MTU在膝关节额状面、矢状面及横截面随时间变化的关节角度，b10-13、b20-23、b30-33均为常系数；
在t时刻，第i个MTU的肌纤维全长为：
Figure BDA0003190823250000055
所述肌纤维速度求解模块将第i个MTU的肌纤维长度对时间求导，获得在t时刻，第i个MTU随时间变化的肌纤维速度为：
Figure BDA0003190823250000061
所述归一化处理模块对肌纤维长度和肌纤维速度进行归一化处理，其计算公式为：
Figure BDA0003190823250000062
其中，lst为第i个MTU的肌腱松弛长度，lom为第i个MTU的最佳肌纤维长度，
Figure BDA0003190823250000063
为第i个MTU的羽状角；
所述肌肉力臂计算模块将肌纤维长度对相应活动面的关节角度求偏导，获得在膝关节额状面、矢状面及横截面上第i个MTU的肌肉力臂分别为：
Figure BDA0003190823250000064
Figure BDA0003190823250000065
Figure BDA0003190823250000066
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的主动力和肌纤维长度的数学模型如下：
Figure BDA0003190823250000067
其中，δ1为常系数，li m(t)为第i个MTU归一化之后的肌纤维长度；
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的主动力和肌纤维速度的数学模型如下：
Figure BDA0003190823250000068
其中，δ2为常系数；
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的被动力与肌纤维长度变化的数学模型如下：
Figure BDA0003190823250000071
其中，δ3为常系数；
所述肌肉力计算模块将每个MTU都建模为Hill肌肉模型，肌纤维与非线性弹性肌腱单元串联；第i个MTU产生的肌肉力为：
Figure BDA0003190823250000072
其中，各变量均对于第i个MTU而言，Fi MTU(t)为产生的肌肉力，是肌肉激活ai(t)、肌纤维长度li m(t)和肌纤维收缩速度vi m(t)的函数，Fi max为最大自主收缩力，
Figure BDA0003190823250000073
为羽状角，f(li m)为肌纤维长度和肌肉主动力之间的关系、f(vi m)为肌纤维速度和肌肉主动力之间的关系，fP(li m)为肌纤维长度和肌肉被动力之间的关系；
所述肌肉力矩计算模块建立以全部MTU在指定活动面引起的膝关节力矩如下：
Figure BDA0003190823250000074
其中，n为MTU数量，i表示第i条MTU，
Figure BDA0003190823250000075
表示t时刻在指定活动面上由全部关注的MTU引起的肌肉力矩之和，
Figure BDA0003190823250000076
为第i条MTU作用在指定活动面的肌肉力，
Figure BDA0003190823250000077
为对应的肌肉力臂；肌肉力臂是支点到肌拉力线的垂直距离，随着肌拉力角的增大而增大，肌拉力角为肌肉在运动骨上的附着点到关节中心的连线与肌拉力线之间的夹角。
进一步地，按如下方法建立外部内收力矩模型：建立待训练的长短时记忆神经网络LSTM，然后将预处理后的膝关节角度训练数据和足底压力训练数据作为输入，将利用OpenSim模型仿真软件中的逆向动力学分析工具ID计算出来的内收力矩作为金标准，对LSTM网络进行训练，得到外部内收力矩模型，用于估计由外部引起的内收力矩。
进一步地，采用遗传算法对接触力估计模型进行参数优化，第i条MTU需要进行优化的参数包括C1i、C2i、Ai这3个肌肉激活系数以及最大等长力Fmax、肌腱松弛长度lst、最佳肌纤维长度lom、羽状角
Figure BDA0003190823250000081
这4个肌肉参数。
进一步地，采用参数优化后的接触力估计模型计算膝关节接触力，包括肌肉力矩计算、外部内收力矩计算以及膝关节内、外侧接触力计算；
首先进行肌肉力计算，在指定关节活动面，将产生的肌肉力与相应肌肉力臂相乘分别求得i条MTU在单个肌肉作用下的产生的膝关节内、外侧肌肉力矩；i条MTU产生的肌肉力矩求和以确定全部MTU在关节处产生的肌肉总力矩；
对于额状面的膝关节内收-外展力矩：
对于膝关节外侧接触点(LC)：
Figure BDA0003190823250000082
对于膝关节内侧接触点(MC)：
Figure BDA0003190823250000083
其中，n为MTU总数量，Fi MTU(t)为第i个MTU在t时刻产生的力量；
Figure BDA0003190823250000084
Figure BDA0003190823250000085
分别为第i个MTU相对于内/外髁突接触点的内收-外展力臂，
Figure BDA0003190823250000086
Figure BDA0003190823250000087
分别为t时刻内、外侧的额状面肌肉力矩；
然后将未经过训练的膝关节角度与足底压力数据输入到外部内收力矩模型中进行计算，得出由外部接触力施加的膝关节内收力矩；
最后通过额状面力平衡关系计算内、外侧接触力：
Figure BDA0003190823250000088
Figure BDA0003190823250000089
其中，dIC为髁突长度，
Figure BDA00031908232500000810
分别为t时刻内、外侧的肌肉力矩，
Figure BDA00031908232500000811
分别为内、外侧的外部内收力矩，FMC/frontal(t)、FLC/frontal(t)分别为膝关节在额状面上的内、外侧接触力。
本发明还提供了一种肌电信号驱动的膝关节肌骨模型接触力估计系统，包括存储器、处理器以及存储于存储器上并能够被处理器运行的计算机程序指令，当处理器运行该计算机程序指令时，能够实现上述的方法步骤。
与现有技术相比，本发明具有以下有益效果：提供了一种肌电信号驱动的膝关节肌骨模型接触力估计方法及系统，该方法及系统使用方便，所需输入数据通过穿戴设备即可获得，使用足底压力数据代替地面反作用力(GRF)数据，以关节角度数据代替标记轨迹数据作为内收力矩模型的输入数据，两种数据可以直接利用加速度计、压力传感器进行测量，摆脱了步态实验室的约束，解决了原始信号计算复杂的问题，拓宽了接触力评估的适用面。此外，通过对模型参数进行优化，使基于模型的力矩、接触力估算更加准确，降低了接触力估计误差。本发明可用于人机交互控制系统中，通过足底压力、膝关节角度及肌电信号进行膝关节屈伸力矩，内收力矩，内、外侧接触力的估计。在此基础上，用户可以明确对应步态模式下的肌肉激活模式、所用肌力大小以及膝关节接触力大小，进而根据建议自主进行步态调整，缓解对植入假体的磨损。
附图说明
图1是本发明实施例的方法实现流程图；
图2是本发明实施例中三类信号的预处理流程图；
图3是本发明实施例中“拧紧旋转”时PZ最大值示意图；
图4是本发明实施例中肌肉激活度模型的工作流程图；
图5是本发明实施例中肌肉收缩力学模型的工作流程图；
图6是本发明实施例中外部内收力矩模型的建模流程图；
图7是本发明实施例中模型参数优化的流程图；
图8是本发明实施例中右腿膝关节内、外侧受力示意图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
应该指出，以下详细说明都是示例性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
如图1所示，本实施例提供了一种肌电信号驱动的膝关节肌骨模型接触力估计方法，包括以下步骤：
1)建立接触力估计模型，包括信号处理模型、OpenSim个体模型、肌肉激活度模型、肌肉收缩力学模型以及外部内收力矩模型，所述信号处理模型包括肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块、步态周期划分模块以及时间归一化处理模块，所述肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块分别用于对肌电信号数据、关节角度信号数据和足底压力信号数据进行预处理，所述步态周期划分模块用于进行步态周期划分，所述时间归一化处理模块用于进行时间归一化处理，以统一步态周期长度；所述OpenSim个体模型用于将模型库中的通用模型缩放至个体特定的生理参数；所述肌肉激活度模型用于通过肌电信号获得肌肉激活度；所述肌肉收缩力学模型用于获取肌肉力矩；所述外部内收力矩模型用于获取外部因素引起的内收力矩。
2)对接触力估计模型进行参数优化。
3)获取肌电信号数据、关节角度信号数据和足底压力信号数据，将其作为输入信号输入到参数优化后的接触力估计模型中，计算得到膝关节接触力。
在本实施例中，通过无线表面肌电传感器获取原始肌电信号数据，通过角度传感器获取原始关节角度信号数据，通过放置于足底的压力传感器获取原始足底压力信号数据。
第一步：建立接触力估计模型
1、信号处理模型
此模型包括五个模块，分别是：肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块、步态周期划分模块以及时间归一化处理模块。
①肌电信号数据预处理模块
原始肌电信号数据来自无线表面肌电传感器，首先，采用4阶巴特沃斯零相移带通滤波器(截止频率为0.08-0.9Hz)，去除直流分量；并且，进行全波整流将信号取正，便于后期数据处理；随后，为模拟肌肉特性，进行4阶巴特沃斯零相移低通滤波(截止频率为随步态周期变化而变化)；由此产生的线性包络相对于从整个试验记录得到的滤波处理过的峰值肌电信号值进行归一化。最后，重采样至与另外两种信号相同的指定频率。
②关节角度数据预处理模块
原始关节角度数据来自角度传感器，需要对原始数据进行平均滤波，以便消除在数据采集时产生的抖动，并重采样至与另外两种信号相同的指定频率。
③足底压力数据预处理模块
原始足底压力数据来自放置于足底的压力传感器，提取压力传感器测得的力在额状面上垂直于地面的分量作为个体足底压力数据。需要对原始数据进行平均滤波，以便消除在数据采集时产生的抖动，并重采样至与另外两种信号相同的指定频率。
三类信号的处理流程如图2所示。
④步态周期划分模块
三种信号数据预处理完成后，可根据应用情况决定是否调用步态周期划分模块。步态周期摆动相末期向支撑相转化期间，腿部近似伸直，胫骨围绕自身轴线向外轻微转动，发生膝关节“拧紧旋转”现象。此时，足底压力中心与膝关节中心的连线在人体额状面的投影距离P达到最大值，矢状面法向Z上的分量也达到最大，如图3所示。滤波等处理完成后，可以将PZ的最大值出现的位置或膝关节外旋角度最大值处统一作为该组步态数据的起点，完成步态周期的划分。
⑤时间归一化处理模块
步态周期划分后，需要对数据进行时间归一化处理，统一步态周期长度，为后续步态周期内膝关节接触力的计算分析奠定基础。
2、OpenSim个体模型
模型仿真需要利用开源OpenSim软件进行个体建模，将模型库中的通用模型缩放至个体特定的生理参数，如身高、体重等，以便接下来通过逆向运动学和逆向动力学分析，得到膝关节内收力矩的金标准数据，进而进行模型参数优化。
3、肌肉激活度模型
基于肌电信号的神经-肌肉-骨骼系统模型，需要通过肌电信号生成肌肉激活。表面肌电信号(Surface electromyography，sEMG)是表面电极募集到肌肉的动作电位总和，肌肉激活度是反映肌肉状态的模型，用来反映肌肉的激活水平。首先要建立sEMG与神经激活度之间的关系，归一化后的肌电信号取值范围为[0,1]，1为完全激活状态，利用归一化处理后的肌电信号求解神经激活度，然后再建立神经激活度与肌肉激活度之间的非线性模型求解得到肌肉激活度。肌肉激活度数据用于后续计算，处理流程如图4所示。
使用临界阻尼的线性离散二阶差分系统来描述sEMG和神经激活度之间的关系：
ui(t)＝αei(t-d)-β1ui(t-1)-β2ui(t-2)
其中，ui(t)、ui(t-1)、ui(t-2)分别是第i块肌肉在第t、t-1、t-2个时刻经过递归滤波器计算后的神经激活度，ei(t-d)表示第i块肌肉在第t-d个时刻经过滤波整流后的肌电信号，i＝1,…,n，在本实施例中，n＝6；d表示机电延迟，α是肌肉的增益系数，β1和β2是肌肉的递归系数。为了实现方程的稳定解和肌肉激活度位于[0,1]之间，α、β1、β2的约束为：
β1＝C1+C2、β2＝C1×C2且α-β1-β2＝1，|C1|＜1、|C2|＜1
通过非线性关系建立肌肉激活和神经激活度之间的关系：
Figure BDA0003190823250000121
其中，ai(t)为第i块肌肉在第t个时刻肌肉激活度，Ai为第i块肌肉的非线性因子，取值范围为(-3,0)。在本实施例中，初始值取-1.5。
4、肌肉收缩力学模型
此模型包括7个模块，分别是：肌纤维长度拟合、肌纤维速度求解、归一化处理、肌肉力臂计算、建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型、肌肉力计算以及肌肉力矩计算，应用流程如图5所示。
①肌纤维长度拟合
肌肉带动关节运动，在这个过程中，肌纤维的长度会随着关节角度的变化而变化。然而实际中很难直接测得肌纤维长度，只能通过诸如磁共振成像(MRI)系统或超声波扫描仪之类的医疗设备来测量，实验室环境中一般通过肌肉的运动轨迹来确定li mt。
在本实施例中，肌纤维长度拟合模块采用三个活动面的膝关节角度进行三次多项式拟合得到肌肉肌腱单元(Muscletendonunit，MTU)的肌纤维长度li mt。
在t时刻，对于跨越一个关节的第i个MTU，将与膝关节在额状面的内收-外展运动、在矢状面的屈曲-伸展运动、在横截面的内旋-外旋运动相关的MTU长度分别拟合为三次多项式：
Figure BDA0003190823250000122
Figure BDA0003190823250000123
Figure BDA0003190823250000124
其中，θi/frontal(t)、θi/sagittal(t)、θi/cross(t)为第i个MTU在膝关节额状面、矢状面及横截面随时间变化的关节角度，b10-13、b20-23、b30-33均为常系数；
在t时刻，第i个MTU的肌纤维全长为：
Figure BDA0003190823250000131
对于跨越两个关节的MTU，为了简化模型，将跨越踝关节或髋关节的双关节肌肉也用单关节肌肉的近似表达式来表示。采用这样的处理方式在实际操作过程中可以只采集膝关节的相关角度，减少传感器的使用，在建模过程中降低了模型的计算复杂度。
运用OpenSim对个体模型进行人体运动模拟仿真后，导出各MTU随时间变化的的肌纤维长度值li mt(t)。对导出的数据使用三次多项式对li mt(t)进行拟合，以得到三次多项式的系数b10-13、b20-23、b30-33。在拟合过程中关节角度为弧度制。
②肌纤维速度求解
动作过程中，关节角度为时间的函数，肌纤维速度可以看作肌纤维长度对时间求导。在计算过程中关节角速度为弧度制。将第i个MTU的肌纤维长度对时间求导，获得第i个MTU随时间变化的肌纤维速度：
Figure BDA0003190823250000132
③归一化处理
肌肉的力学特性主要为肌纤维长度-肌张力关系以及肌纤维速度-肌张力关系，需要得到归一化的肌纤维长度和肌纤维速度。
对肌纤维长度和肌纤维速度进行归一化处理具体使用公式为：
Figure BDA0003190823250000133
其中，lst为第i个MTU的肌腱松弛长度，lom为第i个MTU的最佳肌纤维长度，
Figure BDA0003190823250000134
为第i个MTU的羽状角；归一化后的肌纤维长度范围应满足[0.5,1.3]，归一化后的肌纤维速度取值范围为[-1,1]，才能符合生理运动规律。对于模型的初值，生理参数选用临床标准值。
④肌肉力臂计算
对于肌肉力臂而言，各力臂可以看作肌纤维长度对相应活动面关节角度的一阶偏导数，获得在膝关节额状面、矢状面及横截面上第i个MTU的肌肉力臂分别为：
Figure BDA0003190823250000141
Figure BDA0003190823250000142
Figure BDA0003190823250000143
⑤建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型
在得到归一化后的肌纤维长度和肌纤维速度之后，建立肌纤维长度-肌肉张力的关系、肌纤维速度-肌肉张力的数学关系模型。
肌肉通过肌腹的收缩产生了主动力，归一化后主动力会跟随肌纤维长度变化进行变化，肌纤维归一化长度为1对应的是肌纤维在机体内静止时的长度，其余对应的是对应与该静止状态下肌纤维的各种长度，收缩或者伸展时对应的肌力。肌纤维归一化长度在[0,1]时，随着肌纤维长度增加，肌肉主动力逐渐增加。肌纤维归一化长度[1,2]时，肌纤维伸展，肌肉主动力逐渐减小。建立第i个MTU的主动力和肌纤维长度的数学模型：
Figure BDA0003190823250000144
其中，δ1为常系数，取值为1.35；li m(t)为第i个MTU归一化之后的肌纤维长度，取值范围为[0,2]，正常活动情况下的范围为[0.5,1.3]。
肌肉产生的总肌力为主动力和被动力的合力。根据主动力和肌纤维长度的关系可知，当肌纤维归一化长度为1时，即肌纤维在体内静止时对应的肌纤维速度为0，随之肌纤维收缩到伸展，肌肉力也随之增加。肌纤维速度达到最大时，归一化肌肉力达到1.8。因此建立的第i个MTU的主动力和肌纤维速度的数学模型为：
Figure BDA0003190823250000145
其中，δ2为常系数，取值为1.8。
肌肉的被动张力是肌腱拉伸产生的，肌腱本身并不能收缩产生肌力，在肌纤维收缩产生力时，肌腱会被动受到肌肉纤维的拉伸，以此来增加肌肉的长度，将力传至骨骼带动骨骼运动。当归一化的肌纤维长度在[1,2]时，肌腱被拉伸产生被动肌力，被动力随着肌纤维的伸长而逐渐增大。第i个MTU的被动力与肌纤维长度变化的数学模型建立如下：
Figure BDA0003190823250000151
其中，δ3为常系数，取值为1.001。
⑥肌肉力计算
每个MTU都建模为Hill肌肉模型，肌纤维与非线性弹性肌腱单元串联。第i个MTU产生的肌肉力为：
Figure BDA0003190823250000152
其中，各变量均对于第i个MTU而言，Fi MTU(t)为产生的肌肉力，是肌肉激活ai(t)、肌纤维长度li m(t)和肌纤维收缩速度vi m(t)的函数，Fi max为最大自主收缩力，
Figure BDA0003190823250000153
为羽状角，f(li m)为肌纤维长度和肌肉主动力之间的关系、f(vi m)为肌纤维速度和肌肉主动力之间的关系，fP(li m)为肌纤维长度和肌肉被动力之间的关系。
⑦肌肉力矩计算
关节力矩为选定肌肉共同作用的结果。对于膝关节而言，胫骨和股骨在肌肉拉力的作用下，沿着膝关节进行运动。因此，以全部MTU在指定活动面引起的膝关节力矩为：
Figure BDA0003190823250000154
其中，n为MTU数量，i表示第i条MTU，
Figure BDA0003190823250000155
表示t时刻在指定活动面上由全部关注的MTU引起的肌肉力矩之和，
Figure BDA0003190823250000156
为第i条MTU作用在指定活动面的肌肉力，
Figure BDA0003190823250000157
为对应的肌肉力臂；肌肉力臂是支点到肌拉力线的垂直距离，随着肌拉力角的增大而增大，肌拉力角为肌肉在运动骨上的附着点到关节中心的连线与肌拉力线之间的夹角。以肌肉引起的在额状面的膝关节内收-外展力矩为例：
Figure BDA0003190823250000161
5、外部内收力矩模型
长短时记忆神经网络(Long Short-Term Memory Network,LSTM)是一种记忆性很强的深度神经网络，适合于处理和预测时间序列中间隔和延迟较长的重要事件，对步态信息此类具有时序性的非线性模型拟合效果很好。
建立待训练的LSTM网络，将预处理后的膝关节角度数据和足底压力数据作为输入，将利用OpenSim模型仿真软件中的逆向动力学分析工具(Inverse dynamics，ID)计算出来的额状面力矩分量(内收力矩)作为金标准，对LSTM网络进行训练，得到外部内收力矩模型，用于对外部因素(如地面)导致的内收力矩的估计。建模流程如图6所示。
第二步：模型参数优化
本实施例采用遗传算法进行，目的是对第i条MTU识别1组肌电激活系数C1i、C2i、Ai和4项肌肉参数最大等长力Fi max、肌腱松弛长度li st、最佳肌纤维长度li om、羽状角
Figure BDA0003190823250000162
以便最大程度降低力矩估计的误差，如图7所示。
公开的通用肌骨模型是由多个人体数据的平均值进行初始化的，对于每条MTU的最大等长力Fi max、肌腱松弛长度li st、最佳肌纤维长度li om以及羽状角
Figure BDA0003190823250000163
4项肌肉参数均使用经验值。因此，最初时刻的模型肌肉参数不能很好地匹配用户个体，需要进行模型校准。
Fi max一般允许50％-200％的变化范围。li st基于肌纤维长度和关节位置之间的测量关系，一般约束在初始值的±15％内。
Figure BDA0003190823250000164
会随着瞬时肌纤维长度而变化，优化方程为：
Figure BDA0003190823250000165
羽状角的范围一般在[5°,30°]。最佳肌纤维长度变化与肌肉激活相关，当肌肉激活最大的时候，就会对应最佳肌纤维长度。
li om_after(t)＝li om(γ(1-ai(t))+1)
其中，
Figure BDA0003190823250000166
在建立基于Hill模型的膝关节力矩估计模型时，选取股直肌、股中间肌、股外侧肌、外侧腓肠肌、半膜肌以及股二头肌长头6块肌肉作为分析对象，首先分别对各MTU的4项肌肉参数使用经验值，然后对每条MTU使用遗传算法优化C1i、C2i、Ai，i＝1,…,6这3个参数，直至目标函数达到最小，其中C1i、C2i的约束范围为(-1,1)，Ai的约束范围均为(-3,0)。
通过遗传算法进行模型参数优化的目的是对每条MTU识别1组肌电激活系数和4项肌肉参数，利用这些肌肉参数和系数，单个MTU产生的肌肉力乘以相对于关节各自的力臂并进行加和时，得到的预测力矩更接近从逆动力学ID计算得到的肌肉力矩，目标函数J为：
Figure BDA0003190823250000171
其中，k为采样点编号，共m个采样点，
Figure BDA0003190823250000172
为通过OpenSim中的肌肉分析工具获得模型内收力矩，
Figure BDA0003190823250000173
为通过OpenSim中的肌肉分析工具获得模型屈伸力矩，两者作为目标函数中的金标准，通过模型计算出的
Figure BDA0003190823250000174
随着优化过程的不断进行逐渐逼近
Figure BDA0003190823250000175
和
Figure BDA0003190823250000176
创建初始种群z。种群规模一般情况下在20-160之间，种群太大虽然更容易得到全局最优解，但是运行时间也会相应变长，本方法最终选择100。对于建立的初始种群，种群里面的每个个体都包括C1i、C2i、Ai共3×i个元素，针对每个元素，使用距离方式(<0.5赋值0，否则为1)进行编码。
对初始种群的二进制进行解码，然后带入目标函数中，得到100个目标函数值(即适应度)J1、J2、…、J100，并计算他们的平均值：
Figure BDA0003190823250000177
其中，k为当前个体编号，j为遗传代数编号。通过轮盘赌法，根据适应度大小，选择复制到下一代的个体。具体做法为：计算种群中每个个体的适应度，然后计算各个个体的选择概率：
Figure BDA0003190823250000181
根据适应度大小排序得到的概率值，通过轮盘赌法产生下一代。该过程中还涉及交叉和变异的操作，交叉是产生新个体的主要方式，通过父代某部分互换得来。变异提高了种群的多样性，能够防止过度早熟。遗传算法的收敛有两种方式，其一是满足最大进化代数；其二是多次连续代之间的适应度的差值小于某个设置的常数。
使用平均绝对误差MAE和相关系数CC来评估系统在关节力矩估计方面的性能。
Figure BDA0003190823250000182
Figure BDA0003190823250000183
其中，m代表采样点个数。
第三部：膝关节接触力估计
当前两步处理完成后，即可进行第三步的膝关节接触力估计。此步骤包括三个方面，分别是：肌肉力矩计算、外部内收力矩计算以及膝关节内、外侧接触力计算。
①肌肉力矩计算
使用参数优化后的模型进行肌肉力计算，在指定关节活动面，将产生的肌肉力与相应肌肉力臂相乘分别求得i条MTU在单个肌肉作用下的产生的膝关节内、外侧肌肉力矩；i条MTU产生的肌肉力矩求和以确定全部MTU在关节处产生的肌肉总力矩。以额状面的膝关节内收-外展力矩为例：
对于膝关节外侧接触点(LC)：
Figure BDA0003190823250000184
对于膝关节内侧接触点(MC)：
Figure BDA0003190823250000185
其中，n为MTU总数量，Fi MTU(t)为第i个MTU在t时刻产生的力量；
Figure BDA0003190823250000186
Figure BDA0003190823250000187
分别为第i个MTU相对于内/外髁突接触点的内收-外展力臂，
Figure BDA0003190823250000188
Figure BDA0003190823250000189
分别为t时刻内、外侧的额状面肌肉力矩；
②外部内收力矩计算
将未经过训练的膝关节角度与足底压力数据输入到训练好的LSTM外部内收力矩模型中进行计算，得出由外部接触力施加的膝关节内收力矩。
③通过额状面力平衡关系计算内、外侧接触力
MTU产生的肌肉力矩和外部内收力矩之差除以髁突长度得到膝关节内、外侧接触力。膝关节内、外侧接触力的受力如下图8所示(以右腿为例)。
Figure BDA0003190823250000191
Figure BDA0003190823250000192
其中，dIC为髁突长度，
Figure BDA0003190823250000193
分别为t时刻内、外侧的肌肉力矩，
Figure BDA0003190823250000194
分别为内、外侧的外部内收力矩，FMC/frontal(t)、FLC/frontal(t)分别为膝关节在额状面上的内、外侧接触力。
本发明基于肌电信号、关节角度和足底压力数据，利用OpenSim搭建个体适应性肌骨模型，利用MATLAB进行数据处理、参数优化和计算分析，通过肌电驱动模型的方式实现对膝关节接触力的预测。
本发明通过API接口将OpenSim建模能力与MATLAB数据处理方面的优势相结合，将肌电信号数据、采集自穿戴式加速度计的膝关节角度数据及采集自放置于鞋垫式压力传感器的足底压力数据作为输入信号，利用遗传算法迭代调整肌肉参数和激活系数，对最大等长力、肌腱松弛长度、最佳肌纤维长度以及羽状角4项肌肉参数进行优化，实现对个体肌骨模型的校准，进而求得肌肉力矩，结合通过LSTM模型求得的外部内收力矩，对膝关节处的力平衡关系进行分析实现对人体内部膝关节内、外侧接触力的估计。
不同的肌肉激活模式驱动肌骨模型产生不同的运动变化，也必定会引起关节接触力的改变，因此在使用本方法进行膝关节接触力估计之前，需要获得对应时段的肌电信号、关节角度和足底压力数据。本方法拟利用MATLAB将建模与运算封装在一个.m文件中，应用时只需将肌电信号、关节角度和足底压力，这三种可以直接通过传感器测量得到的信号导入至系统，即可实现对应时刻膝关节内、外侧接触力的估计，并将结果可视化。
本实施例还提供了一种肌电信号驱动的膝关节肌骨模型接触力估计系统，包括存储器、处理器以及存储于存储器上并能够被处理器运行的计算机程序指令，当处理器运行该计算机程序指令时，能够实现上述的方法步骤。
进行膝关节内、外侧接触力在线估计时，通过加速度计、压力传感器及肌电传感器获得用户的关节角度、足底压力以及选定肌肉的肌电信号数据。随后，只需要对目标时段的三类输入数据输入至系统，进行适当的滤波预处理(是否进行步态周期划分可根据需要决定)，在完成参数优化的模型中即可实现膝关节屈伸力矩、膝关节内收力矩以及膝关节内、外侧接触力的估计，如图1所示。根据可视化界面的分析提示，用户可以明确对应步态模式下的肌肉激活模式、所用肌力大小以及膝关节接触力大小，进而判断对关节假体的磨损影响，根据建议自主进行步态调整。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (7)

1.一种肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，包括以下步骤：
建立接触力估计模型，包括信号处理模型、OpenSim个体模型、肌肉激活度模型、肌肉收缩力学模型以及外部内收力矩模型，所述信号处理模型包括肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块、步态周期划分模块以及时间归一化处理模块，所述肌电信号数据预处理模块、关节角度数据预处理模块、足底压力数据预处理模块分别用于对肌电信号数据、关节角度信号数据和足底压力信号数据进行预处理，所述步态周期划分模块用于进行步态周期划分，所述时间归一化处理模块用于进行时间归一化处理，以统一步态周期长度；所述OpenSim个体模型用于将模型库中的通用模型缩放至个体特定的生理参数；所述肌肉激活度模型用于通过肌电信号获得肌肉激活度；所述肌肉收缩力学模型用于获取肌肉力矩；所述外部内收力矩模型用于获取外部内收力矩；
对接触力估计模型进行参数优化；
获取肌电信号数据、关节角度信号数据和足底压力信号数据，将其作为输入信号输入到参数优化后的接触力估计模型中，计算得到膝关节接触力；
所述肌肉收缩力学模型包括肌纤维长度拟合模块、肌纤维速度求解模块、归一化处理模块、肌肉力臂计算模块、建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块、肌肉力计算模块以及肌肉力矩计算模块；所述肌纤维长度拟合模块采用膝关节角度进行三次多项式拟合得到肌肉肌腱单元MTU的肌纤维长度；在t时刻，对于跨越一个关节的第i个MTU，将与膝关节在额状面的内收-外展运动、在矢状面的屈曲-伸展运动、在横截面的内旋-外旋运动相关的MTU长度分别拟合为三次多项式：
Figure FDA0003526048210000011
Figure FDA0003526048210000012
Figure FDA0003526048210000013
其中，θi/frontal(t)、θi/sagittal(t)、θi/cross(t)为第i个MTU在膝关节额状面、矢状面及横截面随时间变化的关节角度，b10-13、b20-23、b30-33均为常系数；
在t时刻，第i个MTU的肌纤维全长为：
Figure FDA0003526048210000021
所述肌纤维速度求解模块将第i个MTU的肌纤维长度对时间求导，获得在t时刻，第i个MTU随时间变化的肌纤维速度为：
Figure FDA0003526048210000022
所述归一化处理模块对肌纤维长度和肌纤维速度进行归一化处理，其计算公式为：
Figure FDA0003526048210000023
其中，lst为第i个MTU的肌腱松弛长度，lom为第i个MTU的最佳肌纤维长度，
Figure FDA0003526048210000024
为第i个MTU的羽状角；
所述肌肉力臂计算模块将肌纤维长度对相应活动面的关节角度求偏导，获得在膝关节额状面、矢状面及横截面上第i个MTU的肌肉力臂分别为：
Figure FDA0003526048210000025
Figure FDA0003526048210000026
Figure FDA0003526048210000027
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的主动力和肌纤维长度的数学模型如下：
Figure FDA0003526048210000028
其中，δ1为常系数，li m(t)为第i个MTU归一化之后的肌纤维长度；
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的主动力和肌纤维速度的数学模型如下：
Figure FDA0003526048210000029
其中，δ2为常系数；
所述建立肌纤维长度-肌肉张力与肌纤维速度-肌肉张力关系模型模块建立第i个MTU的被动力与肌纤维长度变化的数学模型如下：
Figure FDA0003526048210000031
其中，δ3为常系数；
所述肌肉力计算模块将每个MTU都建模为Hill肌肉模型，肌纤维与非线性弹性肌腱单元串联；第i个MTU产生的肌肉力为：
Figure FDA0003526048210000032
其中，各变量均对于第i个MTU而言，Fi MTU(t)为产生的肌肉力，是肌肉激活ai(t)、肌纤维长度li m(t)和肌纤维收缩速度vi m(t)的函数，Fi max为最大等长收缩力，
Figure FDA0003526048210000033
为羽状角，f(li m)为肌纤维长度和肌肉主动力之间的关系、f(vi m)为肌纤维速度和肌肉主动力之间的关系，fP(li m)为肌纤维长度和肌肉被动力之间的关系；
所述肌肉力矩计算模块建立以全部MTU在指定活动面引起的膝关节力矩如下：
Figure FDA0003526048210000034
其中，n为MTU数量，i表示第i条MTU，
Figure FDA0003526048210000035
表示t时刻在指定活动面上由全部关注的MTU引起的肌肉力矩之和，
Figure FDA0003526048210000036
为第i条MTU作用在指定活动面的肌肉力，
Figure FDA0003526048210000037
为对应的肌肉力臂；肌肉力臂是支点到肌拉力线的垂直距离，随着肌拉力角的增大而增大，肌拉力角为肌肉在运动骨上的附着点到关节中心的连线与肌拉力线之间的夹角；
采用参数优化后的接触力估计模型计算膝关节接触力，包括肌肉力矩计算、外部内收力矩计算以及膝关节内、外侧接触力计算；
首先进行肌肉力计算，在指定关节活动面，将产生的肌肉力与相应肌肉力臂相乘分别求得i条MTU在单个肌肉作用下的产生的膝关节内、外侧肌肉力矩；i条MTU产生的肌肉力矩求和以确定全部MTU在关节处产生的肌肉总力矩；
对于额状面的膝关节内收-外展力矩：
对于膝关节外侧接触点(LC)：
Figure FDA0003526048210000041
对于膝关节内侧接触点(MC)：
Figure FDA0003526048210000042
其中，n为MTU总数量，Fi MTU(t)为第i个MTU在t时刻产生的力量；
Figure FDA0003526048210000043
Figure FDA0003526048210000044
分别为第i个MTU相对于内/外髁突接触点的内收-外展力臂，
Figure FDA0003526048210000045
Figure FDA0003526048210000046
分别为t时刻内、外侧的额状面肌肉力矩；
然后将未经过训练的膝关节角度与足底压力数据输入到外部内收力矩模型中进行计算，得出由外部接触力施加的膝关节内收力矩；
最后通过额状面力平衡关系计算内、外侧接触力：
Figure FDA0003526048210000047
Figure FDA0003526048210000048
其中，dIC为髁突长度，
Figure FDA0003526048210000049
分别为t时刻内、外侧的肌肉力矩，
Figure FDA00035260482100000410
分别为内、外侧的外部内收力矩，FMC/frontal(t)、FLC/frontal(t)分别为膝关节在额状面上的内、外侧接触力。
2.根据权利要求1所述的肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，通过无线表面肌电传感器获取原始肌电信号数据，通过角度传感器获取原始关节角度信号数据，通过放置于足底的压力传感器获取原始足底压力信号数据。
3.根据权利要求1所述的肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，所述肌电信号数据预处理模块按如下方法对肌电信号数据进行预处理：首先，采用4阶巴特沃斯零相移带通滤波器，去除直流分量；并且，进行全波整流将信号取正，以便于后期数据处理；而后，为模拟肌肉特性，进行4阶巴特沃斯零相移低通滤波；由此产生的线性包络相对于从整个试验记录得到的滤波处理过的峰值肌电信号值进行归一化；最后，重采样至与另外两种信号相同的指定频率；
所述关节角度数据预处理模块按如下方法对关节角度信号数据进行预处理：对原始数据进行平均滤波，以消除在数据采集时产生的抖动；最后，重采样至与另外两种信号相同的指定频率；
所述足底压力数据预处理模块按如下方法对足底压力信号数据进行预处理：提取压力传感器测得的力在额状面上垂直于地面的分量作为个体足底压力数据；对原始数据进行平均滤波，以消除在数据采集时产生的抖动；最后，重采样至与另外两种信号相同的指定频率。
4.根据权利要求1所述的肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，按如下方法建立肌肉激活度模型：
首先建立表面电极募集到肌肉的动作电位总和，即表面肌电信号sEMG与反映肌肉状态的神经激活度之间的关系，利用归一化处理后的肌电信号求解神经激活度，然后再建立神经激活度与肌肉激活度之间的非线性模型，求解得到肌肉激活度；
使用临界阻尼的线性离散二阶差分系统来描述sEMG和神经激活度之间的关系：
ui(t)＝αei(t-d)-β1ui(t-1)-β2ui(t-2)
其中，ui(t)、ui(t-1)、ui(t-2)分别是第i块肌肉在第t、t-1、t-2个时刻经过递归滤波器计算后的神经激活度，ei(t-d)表示第i块肌肉在第t-d个时刻经过滤波整流后的肌电信号，i＝1,...,n，d表示机电延迟，α是肌肉的增益系数，β1和β2是肌肉的递归系数；为了使方程的稳定解和肌肉激活度位于[0,1]之间，α、β1、β2的约束为：
β1＝C1+C2、β2＝C1×C2且α-β1-β2＝1，|C1|＜1、|C2|＜1
通过非线性关系建立肌肉激活和神经激活度之间的关系：
Figure FDA0003526048210000051
其中，ai(t)为第i块肌肉在第t个时刻肌肉激活度，Ai为第i块肌肉的非线性因子，取值范围为(-3,0)。
5.根据权利要求1所述的肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，按如下方法建立外部内收力矩模型：建立待训练的长短时记忆神经网络LSTM，然后将预处理后的膝关节角度训练数据和足底压力训练数据作为输入，将利用OpenSim模型仿真软件中的逆向动力学分析工具ID计算出来的内收力矩作为金标准，对LSTM网络进行训练，得到外部内收力矩模型，用于估计由外部引起的内收力矩。
6.根据权利要求1所述的肌电信号驱动的膝关节肌骨模型接触力估计方法，其特征在于，采用遗传算法对接触力估计模型进行参数优化，第i条MTU需要进行优化的参数包括C1i、C2i、Ai这3个肌肉激活系数以及最大等长收缩力Fmax、肌腱松弛长度lst、最佳肌纤维长度lom、羽状角
Figure FDA0003526048210000061
这4个肌肉参数。
7.一种肌电信号驱动的膝关节肌骨模型接触力估计系统，其特征在于，包括存储器、处理器以及存储于存储器上并能够被处理器运行的计算机程序指令，当处理器运行该计算机程序指令时，能够实现如权利要求1-6任一项所述的方法步骤。
CN202110877393.4A 2021-07-31 2021-07-31 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统 Active CN113576463B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110877393.4A CN113576463B (zh) 2021-07-31 2021-07-31 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110877393.4A CN113576463B (zh) 2021-07-31 2021-07-31 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统

Publications (2)

Publication Number Publication Date
CN113576463A CN113576463A (zh) 2021-11-02
CN113576463B true CN113576463B (zh) 2022-05-10

Family

ID=78253316

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110877393.4A Active CN113576463B (zh) 2021-07-31 2021-07-31 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统

Country Status (1)

Country Link
CN (1) CN113576463B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114431832B (zh) * 2021-12-27 2023-10-20 同济大学 一种肌肉能量消耗的量化分析方法
CN114469142A (zh) * 2022-01-06 2022-05-13 中南大学 一种基于人体肌肉动力学模型和肌电信号的肌力解码方法
CN115687898B (zh) * 2022-12-30 2023-07-11 苏州大学 基于多模态信号的步态参数自适应拟合方法
CN116978564B (zh) * 2023-07-27 2024-04-23 广东省第二中医院（广东省中医药工程技术研究院） 一种股四头肌肌肉质量评估方法、系统、设备及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105288933A (zh) * 2015-11-20 2016-02-03 武汉理工大学 并联下肢康复机器人自适应训练控制方法及康复机器人
CN108324503A (zh) * 2018-03-16 2018-07-27 燕山大学 基于肌骨模型和阻抗控制的康复机器人自适应控制方法
CN109044352A (zh) * 2018-06-22 2018-12-21 福州大学 一种确定预测人体关节力矩的人工智能输入变量的方法
CN109346176A (zh) * 2018-08-27 2019-02-15 浙江大学 一种基于人体动力学建模和表面肌电信号修正的肌肉协同分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8864846B2 (en) * 2005-03-31 2014-10-21 Massachusetts Institute Of Technology Model-based neuromechanical controller for a robotic leg

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105288933A (zh) * 2015-11-20 2016-02-03 武汉理工大学 并联下肢康复机器人自适应训练控制方法及康复机器人
CN108324503A (zh) * 2018-03-16 2018-07-27 燕山大学 基于肌骨模型和阻抗控制的康复机器人自适应控制方法
CN109044352A (zh) * 2018-06-22 2018-12-21 福州大学 一种确定预测人体关节力矩的人工智能输入变量的方法
CN109346176A (zh) * 2018-08-27 2019-02-15 浙江大学 一种基于人体动力学建模和表面肌电信号修正的肌肉协同分析方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo;David G. Lloyd等;《Journal of Biomechanics》;20031231(第36期);第765–776页 *
基于Hill 肌肉模型的人体关节力矩智能预测;熊保平等;《北京生物医学工程》;20210228;第40卷(第1期);第11-23页 *
基于sEMG变阻抗控制的肘关节康复外骨骼系统研究;张腾;《哈尔滨工业大学硕士学位论文》;20210115;第1-65页 *
基于肌电信号的膝关节肌肉力分布模型;李翰君;《北京体育大学博士（毕业）学位论文》;20110915;第1-84页 *

Also Published As

Publication number Publication date
CN113576463A (zh) 2021-11-02

Similar Documents

Publication Publication Date Title
CN113576463B (zh) 肌电信号驱动的膝关节肌骨模型接触力估计方法及系统
Chen et al. Surface EMG based continuous estimation of human lower limb joint angles by using deep belief networks
Fleischer et al. Predicting the intended motion with EMG signals for an exoskeleton orthosis controller
US9975249B2 (en) Neuromuscular model-based sensing and control paradigm for a robotic leg
Erdemir et al. Model-based estimation of muscle forces exerted during movements
Zheng et al. Gait phase estimation based on noncontact capacitive sensing and adaptive oscillators
Buchanan et al. Estimation of muscle forces and joint moments using a forward-inverse dynamics model
US8821417B2 (en) Method of monitoring human body movement
Alamdari et al. A review of computational musculoskeletal analysis of human lower extremities
Liu et al. Lower extremity joint torque predicted by using artificial neural network during vertical jump
Guo et al. Mechanical energy and power flow of the upper extremity in manual wheelchair propulsion
Hossain et al. Deepbbwae-net: A cnn-rnn based deep superlearner for estimating lower extremity sagittal plane joint kinematics using shoe-mounted imu sensors in daily living
Ling et al. Lower limb exercise rehabilitation assessment based on artificial intelligence and medical big data
Nuckols et al. Automated detection of soleus concentric contraction in variable gait conditions for improved exosuit control
Dinovitzer et al. Accurate real-time joint torque estimation for dynamic prediction of human locomotion
Hahn et al. A neural network model for estimation of net joint moments during normal gait
HARAGUCHI et al. Prediction of ground reaction forces and moments and joint kinematics and kinetics by inertial measurement units using 3D forward dynamics model
Wang et al. Qualitative evaluations of gait rehabilitation via EMG muscle activation pattern: Repetition, symmetry, and smoothness
Raison On the quantification of joint and muscle efforts in the human body during motion
Weng et al. A gait stability evaluation method based on wearable acceleration sensors
Khanna Modelling and Simulation of Ideal and Torque-limited Devices to Assist Elderly Walking
Smith Modeling human dynamics for powered exoskeleton control
Slaboda et al. The use of splines to calculate jerk for a lifting task involving chronic lower back pain patients
Liang Human Locomotion Synergy for Robotic Assistive Devices
McCabe Utilizing Neural Networks and Wearables to Quantify Hip Joint Angles and Moments During Walking and Stair Ascent

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant