融合生物动力学和NMF的下肢运动建模

戴海烽，张 超，张建海

(1.杭州电子科技大学计算机学院,浙江 杭州 310018；2.浙江中医药大学附属第一医院,浙江 杭州 310018)

0 引 言

越来越多可穿戴康复设备[1]引入了电刺激[2]反馈功能，电刺激反馈效果和肢体运动姿态与肌肉激活程度模型息息相关。肢体康复训练中，根据健康人正常肢体运动姿态的肌肉激活程度模型，对患者特定肌肉群实施有针对性的电刺激，为肢体运动进行正确、科学的电刺激反馈提供了指导性的参考。王江等[3]运用径向基(Radial Basis Function，RBF)神经网络对肌电信号进行建模仿真，该方法非常适合肌电信号非线性系统的建模，但RBF神经网络的理论基础是统计数学，缺乏生理意义。李强[4]运用肌纤维募集模型仿真第一背间肌运动神经元在不同激励水平下的表面肌电(Surface ElectroMyoGraphy，SEMG)信号和相应的收缩能力，研究发现仿真的动作电位信号与真实信号相似性极高。James等[5]通过研究肌梭反馈模型发现，反馈信号与肌肉长度之间具有非常高的相关性。虽然，这些生理学模型能够描述肌肉肌电与关节运动的物理关系，但是存在一定的片面性，拟合效果不如多模型融合方法。近年来，多模型融合方法开始运用于主动肌建模，白杰[6]使用一种由内驱力、高尔基腱、肌梭融合的模型拟合肌电，拟合误差与神经网络方法相当，但由于部分子模型为数学模型，拟合参数过多，导致3个子模型拟合的肌电成分不稳定，失去子模型该有的生理意义。在构建肢体运动姿态的肌肉激活模型时，多数以主动肌为主，辅助肌的作用同样不可忽视，对完善关节运动与肌肉群间的研究同样重要。协同分解[7]，例如非负矩阵分解(Nonnegative Matrix Factorization，NMF)[8]常常用于分析肌肉间隐含的协同关系。李飞[9]将NMF应用于小儿脑瘫的下肢步态分析中，为患者的康复训练提供评估和指导。本文使用运动配置策略以及膝关节运动相关损失函数改进融合生物动力学模型，采用NMF对辅助肌肌电进行建模，构建一个具有生理意义、全面性、肌电拟合效果好的下肢膝关节肌电模型。

1 融合生物动力学的下肢主动肌建模

在构建下肢运动姿态与肌肉激活关系模型时，为了使子模型生成较稳定的肌电，需选择模型系数较少或有限定范围的生理意义模型。同时，为了进一步加强各个子模型的稳定性，本文使用模型配置策略，根据下肢膝关节的运动特点使用特定的损失函数进行模型系数的搜索。

1.1 基本生理意义模型

1.1.1 意图-肌电子模型

意图-肌电子模型根据运动意图来控制肌肉收缩，肌肉收缩时产生肌肉肌电，即运动意图生成肌电，其过程如下：

ΦRTE=eα i
u(t)=wΦRTE

(1)

式中，ΦRTE为运动神经单元募集阈值估计(Raising Threshold Estimation，RTE)，w为放电增益，u(t)为最终的意图肌电。i为募集运动神经单元的指数系数，α为标量化的神经意图，一般无法直接测量，可通过其与角度偏差(Angle Error，AE)的关系求得，即α=kAAE，角度偏差可通过AAE=As+n+As-n-2As求得，即将n时刻前后的角度求和，再与2倍的当前时刻s下的角度As相减，k为神经意图与角度偏差的关系系数。

1.1.2 肌力-肌电子模型

肌力-肌电子模型需先计算主动肌产生的有用力矩信号，依据《中国成年人人体惯性参数国家标准》[10]中相关的回归方程，有用力矩的计算公式为：

(2)

式中，m，l，J分别为估算的小腿重量、重心和转动惯量，X1，X2分别为受试者的体重和身高。θ为关节角度，g为重力加速度。下肢有用力矩L分为2项[11]，分别与当前膝关节角度、角度加速度有关，其余力矩属于内力矩，起到调整运动姿态作用，可忽略。

肌力-肌电子模型选择逆hill模型。hill模型分为肌电-激活信号和激活信号-肌力，逆hill模型也可分为这2部分的相反过程。肌力-激活信号表示如下：

a(t)fe+fp=F/Fmax

(3)

式中，fe为肌肉主动收缩部分的系数，fp为肌肉被动收缩部分系数，可忽略不计，F为力矩信号，Fmax为力矩信号中的最大值，a(t)为所获得的激活信号。

获得激活信号后，再通过激活信号-肌电转化成肌电信号：

g(t)=vln[a(t)(eB-1)+1]/B

(4)

式中，g(t)为预测的肌电信号，v为拟合系数，B为非线性参数，取值范围为[-3，0]，该数值与动作的种类有关。

1.1.3 肌肉长度-肌电子模型

肌肉长度-肌电子模型选用肌梭模型，肌梭是肌肉中的肌肉长度感知器官，所以肌肉长度-肌电子模型是模仿肌梭反馈功能，感知肌肉长度变化并反馈与肌肉长度相关的肌电的模型。首先，针对股四头肌和股二头肌需要建立不同的肌肉长度-角度模型，其中股四头肌的肌肉长度-角度公式为：

LM=LM0+rθ

(5)

式中，LM为股四头肌肌肉长度，LM0为股四头肌初始肌肉长度，r为膝关节中心的旋转半径。

股二头肌的肌肉长度-角度公式为：

(6)

式中，LB为股二头肌肌肉长度，L11为股二头肌起始端到膝关节旋转中心的距离，L12为膝关节旋转中心至股二头肌止端的距离。

股四头肌肌梭模型处理过程如下：

uM-uM0=KMvp(LM-LM0)

(7)

式中，uM为股四头肌肌梭反馈的肌电信号，uM0为初始股四头肌肌梭反馈的肌电信号，v为肌肉长度变化速度，p为肌肉长度变化速度影响因子，通常取0.3，KM为待定系数。

类似地，股二头肌肌梭模型处理过程如下：

uB-uB0=KMvp(LB-LB0)

(8)

式中，uB为股二头肌肌梭反馈的肌电信号，uB0为初始股二头肌肌梭反馈的肌电信号，LB0为股二头肌初始肌肉长度。

1.2 运动配置策略

为了加强1.1节的3个子模型的稳定性，本文采用运动配置策略将各个子模型进行融合。首先，将下肢膝关节运动分成4个阶段，分别为膝关节伸展、膝关节伸展回归、膝关节屈曲和膝关节屈曲回归，如图1所示。其次，分析各个阶段的运动过程和肌电产生特点，以股四头肌为例，相比于其他过程，伸展阶段中，运动神经控制股四头肌收缩的意图最强，意图肌电成分所占比重最大；伸展和伸展回归阶段中，股四头肌主动收缩，控制小腿随膝关节中心旋转并对外做有用功，该过程生成的肌电大部分由主动肌力产生；屈曲和屈曲回归阶段中，股四头肌起协同肌作用，只随膝关节运动，使自身长度发生变化，起到调整运动姿势的作用，该过程生成的肌电基本由因肌肉长度变化产生的被动肌力。最后，依据膝关节运动特点，在伸展阶段配置3个子模型，在伸展回归阶段配置肌力-肌电子模型和肌肉长度-肌电子模型，在屈曲和屈曲回归配置肌肉长度-肌电子模型。

图1 膝关节运动中的4个阶段

与股四头肌类似，股二头肌也依据对应膝关节运动特点，在屈曲阶段配置3个子模型，在屈曲回归阶段配置肌力-肌电子模型和肌肉长度-肌电子模型，在伸展和伸展回归配置肌肉长度-肌电子模型。

1.3 膝关节运动相关损失函数

在融合3子模型中，为了不失去各子模型的生理意义，系数搜索时，需要依据膝关节运动的特点合理修改损失函数。以股四头肌为例，当只使用意图-肌电模型时，侧重运动意图对肌肉收缩的控制，即在伸展阶段，原始肌电与预测肌电间的误差最低；在只使用肌力-肌电模型时，着重考虑主动肌收缩发力对外做功时的情形；在只使用肌肉长度-肌电模型时，着重考虑非主动肌被动发力时的情况。当需要使用融合模型时，需要同时兼顾3个模型对4个阶段不同的参与度，在伸展阶段，意图、肌力、肌肉长度模型都在肌电预测方面具有一定的参与度；在伸展回归阶段，大脑处于放松阶段，意图模型不再参与预测肌电；在屈曲和屈曲回归阶段，股四头肌不再是主动肌，只有肌肉长度模型参与肌电预测。最后，通过最小化损失函数求解得到3个模型中的待定系数，使得模型拟合的肌电更趋近于各个子模型的生理意义。综合上述分析，融合模型的损失函数为：

(9)

式中，Sn(t)为各个子模型预测的肌电信号与真实信号间的损失函数，其中n为1表示意图-肌电子模型，2表示肌力-肌电子模型，3表示肌肉长度-肌电子模型，SF(t)为融合生物动力学模型与真实信号间的损失函数，F表示融合生物动力学模型。uSEMG为真实的表面肌电信号，un为各个子模型或者融合生物动力学模型预测的肌电。tn为膝关节运动阶段，t1表示意图控制运动阶段，在预测股四头肌肌电时，其表示膝伸展阶段，而预测股二头肌肌电时则表示膝屈曲阶段。t2表示肌肉主动收缩对外施展有用力矩的过程，在预测股四头肌肌电时，其表示膝伸展阶段和伸展回归阶段，而预测股二头肌肌电时则表示膝屈曲阶段和屈曲回归阶段。t3表示肌肉被动收缩而随关节运动产生被动肌电的过程，即协同过程，在预测股四头肌肌电时，其表示屈曲和屈曲回归阶段，而预测股二头肌肌电时则表示膝伸展阶段和伸展回归阶段。

2 基于NMF的下肢辅助肌建模

NMF常常用于分析肌肉间的协同情况[12]，与神经系统控制肌肉运动的分析过程相似，即在进行某个动作时，神经系统只控制一个更小的单位的肌肉协同模式而非一整块肌肉，分析式如下：

V=WH+E

(10)

式中，V为待分解的肌电矩阵，W为分解后的协同曲线，H为协同模式，E为噪声误差。使用迭代算法可使得噪声误差逐步趋近0，最基本的是文献[13]提出的基于欧氏距离的乘性迭代算法和基于K-L散度的加性迭代算法。本文采用基于欧氏距离乘性迭代的NMF。

运用NMF分解的主、副肌肉群间的协同效果如图2所示，8个肌肉通道分别对应阔筋膜张肌、股四头肌、股二头肌、内收肌、腓肠肌、胫骨前肌、比目鱼肌内侧、比目鱼肌外侧，协同曲线主模式分别为股四头肌、股二头肌以及阔筋膜张肌。

图2 NMF的分解效果

协同模式表示1组以某块主动肌为主要作用，其余肌肉为辅助作用的肌肉协同模式，对应的协同曲线表示该协同模式的能量随着时间变化的情况。分解肌电后，先将股四头肌与阔筋膜张肌对应协同模式上的特征值相除，再将该值乘以股四头肌肌电，将股四头肌肌电转化为阔筋膜张肌肌电，最后，叠加各自的协同模式，得到其他辅助肌肌电。

3 肌电拟合实验及分析

3.1 实验设备

实验选用8通道博瑞康无线肌电采集设备，采样频率为1 000 Hz，选取肌肉通道为阔筋膜张肌、内收肌、腓肠肌、胫骨前肌、比目鱼肌内侧、比目鱼肌外侧以及2块主动肌。角度传感器选用维特智能九轴角度传感器记录仪，采样频率为200 Hz，如图3所示。

图3 肌电设备和角度传感器

3.2 数据采集

选取4名健康男性受试者，均已签署知情同意书。采集受试者的身高、体质量、静息主动肌长度、大腿长度等生理参数。

下肢膝关节动作设计如图4所示。受试者先保持坐姿平衡，然后做出膝关节伸展运动，再返回至坐姿平衡状态，随后进行膝关节屈曲运动，再次返回平衡状态，共做5次。做完1个轮次后，休息5 min，避免肌肉疲劳。

图4 下肢膝关节运动实验过程

3.3 实验结果

在MATLAB实验环境下，采用本文方法对3.2节所述的膝关节运动实验中的主动/辅助肌肌电进行拟合。采用本文提出的融合生物动力学方法与文献[6]提出的主动肌的融合方法进行股四头肌肌电的拟合实验，文献[6]模型的拟合结果如图5所示，本文方法各子模型的拟合结果如图6所示。

图5 文献[6]模型拟合股四头肌肌电的结果

图6 本文方法子模型拟合股四头肌肌电的结果

从图5中可以看出，文献[6]融合模型的肌电拟合效果不错，但各子模型的肌电拟合不稳定，这是因为子模型的拟合参数过多，致使子模型的肌电成分不稳定，不符合生理意义。

从图6可以看出，意图-肌电模型在伸展阶段拟合情况较好，但其余运动阶段模型预测的肌电都高于原始肌电，这是因为除伸展阶段外的运动阶段过多考虑了意图的参与，拟合的肌电相对偏高。类似地，肌力-肌电模型在伸展阶段未考虑动作调整意图对肌电信号的影响，在屈曲和屈曲回归阶段，肌力-肌电模型又过多考虑主动肌力的参与，致使伸展阶段预测的肌电相对偏低，屈曲和屈曲回归阶段相对偏高。肌肉长度-肌电模型在伸展和伸展回归阶段未考虑主动肌力和动作调整意图对肌电的影响，致使肌肉长度-肌电模型预测的肌电在伸展和伸展回归阶段都相对偏低。

采用本文融合模型拟合股四头肌和股二头肌肌电的结果如图7所示。

图7 本文方法的各个子模型拟合主动肌肌电的结果

从图7可以看出，融合生物动力学模型中，各子模型的肌电拟合比较稳定，且3个子模型拟合的肌电叠加起来的外包络也能很好拟合原始肌电。相比文献[6]模型，本文采用模型配置策略，并运用运动相关损失函数来搜寻模型中的待定系数，使得各子模型拟合与各自模型意义相关的肌电，解决了子模型肌电拟合不稳定问题。

在主动肌肌电模型的基础上，采用本文提出的NMF求得的协同模式对各个辅助肌进行肌电拟合，得到辅助肌模型拟合结果如图8所示。

从图8可以看出，辅助肌肌电拟合情况良好，说明运用主动肌肌电和NMF算法分解的协同模式能有效预测辅助肌肌电，完善了下肢膝关节运动姿态与相应肌肉群激活程度的关系。

本文采用均方根误差(Root Mean Square Error，RMSE)作为定量评价拟合效果的指标。

(11)

式中，d为对比的两段肌电信号的长度。

本文各子模型及融合模型、文献[6]模型的主动肌肌电拟合效果如表1所示，各块辅助肌肌电拟合效果如表2所示。

表1 主动肌建模效果定量评价表 单位：%

表2 辅助肌建模效果定量评价表 单位：%

从表1可以看出，本文融合模型的RMSE值与文献[6]模型相差不大，说明两者的拟合效果相当；本文融合模型的RMSE值都低于单独使用子模型拟合的结果，说明融合模型拟合效果优于单个子模型。

将表2数据取平均值，得到辅助肌的平均RMSE值为4.62%。一般情况下，RMSE值低于8%时，预测的信号能很好地拟合原始信号[14]。因此，运用NMF对辅助肌进行肌肉激活程度建模可以得到良好的肌电拟合效果。

4 结束语

本文提出一种融合生物动力学和NMF的下肢运动建模方法，建立了下肢膝关节运动姿态与肌肉激活程度的关系模型。运用运动配置策略，采用运动相关损失函数，加强了子模型的稳定性，并建立基于NMF的辅助肌模型，完善了下肢关节运动与肌肉激活程度模型，为正确、有效进行功能电刺激提供借鉴。后续将尝试更多的分解方式以提升辅助肌的拟合效果，进一步优化下肢关节运动与肌肉激活程度模型。