基于多路sEMG时序分析的人体运动模式识别方法

发布时间：2022-04-19

　　摘 要 针对主动康复训练中人体运动识别问题, 提出了一种基于多路表面肌电 (Surface electromyogram, sEMG) 时序特征的人体运动模式识别方法. 设计评估类周期 sEMG 信号波形相似度的方法来对多路 sEMG 信号进行特征选择; 以二维科荷伦自组织竞争网络 (Self-organization mapping net, SOM) 对多路信息进行编码; 最后, 建立描述各运动过程多路 sEMG 时序特征的隐马尔科夫模型 (Hidden Markov model, HMM), 基于最大似然估计法对多模型匹配进行综合判决获取识别结果. 并在对下肢踏车、椭圆、步行运动模式的识别实验中, 相对于经典线性及非线性算法, 识别率由 72.5 % 和 88.33 % 提高到 91.67 %, 验证了本文方法的有效性.

　　关键词 康复训练, 表面肌电, 人体运动意图识别, 自组织神经网络, 隐马尔科夫模型

　　中风是急性脑血管疾病的统称, 它不但可以造成永久性中枢神经损伤, 还可以引起多种并发症, 甚至死亡. 由于中风发病率高、致残率高、复发率高、死亡率高、并发症多, 是目前世界范围内成人致残的最主要原因以及致死的第二大原因, 世界医学界已经把它列为威胁人类健康的三大疾病之一[1] . 超过 80 % 的中风发病于 65 岁以上的老年人群体[2] , 并且有近 75 % 的幸存者致残[3] . 随着世界人口老龄化进程加剧, 中风造成的神经损伤性残疾已经给家庭和社会带来了沉重的经济和劳动力负担. 医学研究表明, 中枢神经系统的功能性康复依赖于其运动功能的频繁使用, 可以由长期的人体主动运动激励而来, 康复训练对中风后神经系统的恢复以及防止二次中风及其他后遗症具有非常积极的作用[4−5] , 尤其是结合人体运动意图的主动康复训练[6−7] . 因此, 准确、实时、可靠的人体运动意图识别方法的研究对实现基于患者主动意图的康复训练是十分必要的.

　　人体运动依赖于骨骼肌伸缩, 其过程中微伏级的生物电压在肌纤维细胞中产生并且可以在皮肤表面测得. 这种从皮肤表面、以非入侵方式测得的神经肌肉系统活动过程中的微弱电信号被称为表面肌电 (Surface electromyogram, sEMG)[8−9] 信号. 由于 sEMG 实时客观地反映着神经肌肉系统的运动功能指令信息, 因此, 可以被用来评估骨骼肌的运动状态从而实现对人体主动运动过程的解析[10] . 在主动康复训练研究中, 训练模式及功能是根据人体主动运动意图来规划的; 因而 sEMG 成为获取人体运动意图最主要的研究对象之一, 例如获取肢体姿势及运动模式[11−15]、实时肌肉力/力矩[16−18]、关节运动状态[19−21] 等. 根据识别对象的不同, 康复训练系统的功能规划也不同. 第一类, 使用 sEMG 用来识别肢体正在进行和即将进行的运动模式并规划训练系统完成一个特定的运动轨迹[13] , 因此, 还可以用于假肢控制[14] . 例如 Naik 等使用来源于 4 块前臂肌肉的 sEMG 和孪生支持向量机 (Twin support vector machine, Twin SVM) 方法成功对 7 种手指及腕部弯曲动作进行了分类, 得到了平均 84.83 % 的敏感性和 88.1 % 的特异性[11]. Momen 等使用模糊 C-means 聚类法对前臂 sEMG 特征空间进行分割[12] , 对 4 类及 5 类手势的平均识别准确率分别达到 92.7 ± 3.2 % 和 79.9 ± 16.8 %. Khushaba 等采用贝叶斯方法分析腕部的 2 路 sEMG 信号识别右手哪根手指处于弯曲状态, 达到近 90 % 的准确率[13] . 根据人体不同部位生理环境的差异, 对运动模式识别的效果的也不同, 通常来说, 上臂 sEMG 受脂肪、毛发等生理环境影响相对较小, 识别率通常高于其他部位. 人体躯干及下肢 sEMG 受较厚脂肪层、突出的皮肤表层及发毛等影响, 通常干扰较大, 目前少有研究对下肢相似运动模式进行识别. 例如 Miljkovi´c 等使用欧氏空间临近法对左右两块腰骶椎旁肌进行分析, 识别不同姿势下的坐、立两种过程, 仅获得了 60 % ∼ 70 % 的准确率[14] .

　　sEMG – 肌肉力/力矩的数学模型通常较为复杂[16−18] , 某些参数随生理条件和环境变化 (例如肌肉萎缩程度、疲劳程度、皮下脂肪含量、皮肤性质等) 并且不易测量, 因此, 较难投入实际应用. 通过 sEMG 实时估算肢体运动状态实时规划系统的运动轨迹是实现主动康复训练的另一个途径. 例如, Shrirao 等建立了手指屈伸过程中的伸趾浅肌 sEMG 与指关节角度的神经网络模型, 实验获得了较小的均方根误差及 200 ms 的延迟[20] . 我们的研究队伍分别建立了踏车和蹬踏运动过程下肢 7 块主要肌群 sEMG 和髋、膝、踝关节角度的 20 阶 BP (Back propagation) 神经网络模型, 在角度估算准确度上获得了平均 5 ◦ ∼ 6 ◦ 的均方根误差[21] . 然而, 由于传感器实际测得的 sEMG 受到生理、环境等多方干扰, 信噪比很低; 因此, 此类 sEMG – 肢体运动状态实时识别方法通常需要设定具体运动过程模式, 缺乏普适性; 所以目前尚需要配合运动过程模式识别方法共同使用[19−21] . 综上可知, 对运动过程模式的识别是实现主动康复训练系统的基础性必要条件, 而目前研究中的识别准确度尚需改善, 对相似运动过程的识别还存在较大的误判.

　　本文利用多路 sEMG 信息为主动康复训练系统提供一个可靠快速的人体肢体运动意图识别方法, 以下肢运动过程为例, 提出了一种基于多路 sEMG 时序特征分析法的人体下肢运动过程模式识别方法. 首先, 通过实验获取了来自 4 位健康个体的单侧下肢 8 块主要肌群在踏车、行走和椭圆运动 (可模拟步态) 三类相似运动过程中的多路 sEMG 数据样本, 使用了一种滑动巴特沃斯滤波法来减小 sEMG 数据处理的时间延迟; 并设计了一种针对 sEMG 与人体运动关系特点的类周期信号相似度评估方法对多路 sEMG 信号进行特征选择以减小计算的时间和空间复杂度. 其次, 使用自组织神经网络 (Self-organization mapping net, SOM) 对多路 sEMG 信号进行时序编码, 并建立描述不同运动过程中 sEMG 时序变化特点的隐马尔科夫模型 (Hidden Markov models, HMM). 最后, 通过实验验证了本方法可以实时、有效地区分人体相似运动模式.

　　本文章节安排如下: 第 1 节描述了实验数据样本获取过程及多路 sEMG 信号处理方法等; 第 2 节描述了多路 sEMG 特征选择及对其进行时序特征编码的方法; 第 3 节描述了使用隐马尔科夫模型建立描述不同运动过程中多路 sEMG 时序变化特征数学模型的方法; 第 4 节描述了实验过程并对结果进行讨论; 第 5 节总结全文. 1 多路 sEMG 信息获取与处理医学研究表明, 相对于其他运动训练方式, 踏车、行走和椭圆运动对人体下肢行走能力的恢复具有较为显著的积极作用[22−23] , 是目前康复训练中医师推荐的较为普遍的基本训练模式, 但是, 其运动过程中下肢主要肌群 sEMG 特征较为相似, 难以区分. 因此, 本文以踏车、行走和椭圆运动过程为例对人体下肢运动过程模式进行分类, 并且为改善时效性、减少计算复杂度对多路 sEMG 采用时域分析法.

　　1.1 信息获取为获取

　　人体下肢多路 sEMG 数据样本, 本文使用加拿大 Thought Technology 公司研制的 FlexComp Infiniti 表面肌电采集设备 (10 通道, 敏感度: < 0.1 µVRMS, 准确度: ± 0.3 µVRMS ± 4 %, 内置低通滤波器, 原始数据采样频率: 2 048 Hz) 来同步采集运动过程中单侧下肢多路 sEMG 信号.

　　sEMG 传感器贴片电极的配置位置对人体运动模式识别的影响是至关重要的[24] . 为充分分析人体下肢多路 sEMG 与踏车、行走及椭圆运动过程的影响, 本文将左侧下肢浅层主要肌群列入实验范围, 如表 1 所示. 为获取最佳传感数据, 本文实验使用Ag-AgCl 电极并放置于肌腹中心位置, 如图 1 所示 (原图出自文献 [25]), 并使正极与负极电极平行于肌纤维. 4 名健康志愿者参与了实验, 包括 3 名男性, 1 名女性, 年龄 27 ± 4 岁, 身高 173 ± 13 cm, 体重 57.5 ± 12.5 kg, 脂肪含量均为中等. 为增强本文方法的实用性, 实验中并没有去除皮肤表面的毛发、角质, 期望可以获得一个不受人体皮肤表面状况干扰的多路 sEMG— 运动模式识别方法. 实验于健身房利用踏车训练机、跑步机及椭圆机进行数据采集, 如图 2 所示; 并选取人体较为舒适的中等速度 (踏车 43 ± 2 圈/min, 行走 43 ± 2 步/min, 椭圆运动 43 ± 2 圈/min); 为不失普适性, 选取来自每个实验个体分别在踏车、行走、椭圆运动过程中完整的 10 个连续运动周期数据作为实验样本. 部分数据示例如图 3 所示, 已经经过放大.

　　1.2 信号欠采样

　　人体运动呈现较为明显的低频性, 因此, 在对运动模式进行识别研究前, 应对采集的多路 sEMG 信号进行欠采样. 较为常用的时域欠采样方法有直接欠采样、积分欠采样法和均方根 (Root mean square, RMS) 欠采样法等[27] . 对比前两者, RMS 欠采样法不但可以反映信号幅值在时域上的变化特点, 还取决于肌肉复合性因素和肌肉本身的生理、生化过程之间的内在联系, 可以实时反映肌肉活动状态[27] . 因此, 本文选用 RMS 欠采样法.

　　1.3 信号滤波

　　传感器采集到的 sEMG 信号是一种低信噪比的微弱信号, 它不只包含 sEMG, 还包含多种噪声: 不仅有毛发、皮肤、脂肪层以及其他生理因素产生的噪声, 还有工频干扰、直流基线噪声等. 尤其是人体下肢相对于上肢来说, 脂肪层、皮肤和毛发等方面的干扰更为强烈. 因此, 欠采样后的 sEMG 信号仍然很不平稳, 需要尽可能平稳地滤波. 然而经过滤波后的信号不可避免地会在时域范围内延迟于原信号. 为提高识别结果的实时性, 本文针对 sEMG 信号的特点使用了一种滑动巴特沃斯滤波法[28]

　　尽管这是一种可以最大限度地平滑原信号的滤波方法, 但是也会像其他滤波器一样在实际应用中对有限离散信号的滤波造成时间延迟. 滤波器的阶数 n 越高, 延迟越严重; 反之, 则信号平滑性变差. 为获得更平滑信号的同时减小时间延迟, 本研究组设计了一种滑动巴特沃斯滤波法. 其核心思想是通过当前时刻前的数据段信息获得当前时刻的滤波后信息元素, 进而减小信号滤波产生的时间延迟. 记当前时刻为 t, 则本方法过程如图 5 表示, 可以通过 4 个步骤实现:

　　步骤 1. 使用滑动时间窗口 (长度为 L, 且 L ≥ n) 获取一段从时刻 t − L + 1 到 t 的原信号;

　　步骤 2. 使用巴特沃斯滤波器对其进行滤波, 获得滤波后的信号段;

　　步骤 3. 取该滤波后的信号段尾部最后一个元素作为当前时刻 t 的滤波结果元素, 以获得无时间延迟的滤波结果;

　　步骤 4. 滑动时间窗口前移一步.

　　该方法的滤波效果如图 6 所示, 该方法不但可以很好地平滑信号, 还可以有效地减少信号滤波过程产生的延时.

　　2 多路 sEMG 特征提取与时序编码

　　2.1 多路 sEMG 特征选择方法

　　为提高计算效率及识别准确率, 对滤波后的多路 sEMG 信号进行分析以减少冗余特征是必要的. 在运动过程中采集到的 sEMG 信号除包含反映完成该运动过程所必需的运动功能指令信息之外, 还有诸多噪声、生理因素 (例如偶尔的肌肉抖动)、实验者主观因素影响下产生的和该运动无关的肌电信号等. 但是对于同一实验者而言, 在重复完成同种运动模式下的多次完整运动过程中, 反映完成该运动过程所必需的运动功能指令信息的 sEMG 分量之间应具有良好的重复性和周期性. 信号在多次类周期运动中的波形相似度越高, 则其与该种运动的关联性越强, 反映完成该种运动所需的运动功能指令信息分量越大, 干扰分量越小; 反之, 则关联性越差, 或者干扰信号的分量越大.

　　相关知识推荐：模式识别专业论文发sci容易吗

　　但是由于人体运动的随机性, 使得同种运动模式下各完整运动过程的完成时间之间不具备明确的关联关系, 因此, 类周期的 sEMG 信号波形的周期长度相似性[29] 并不能作为多路 sEMG 特征选择的依据. 由于 sEMG 波形幅值反映了肌肉力的大小, 并且其变化趋势反映了肌肉力的变化过程, 是决定运动模式类别的直接因素, 本文提出采用类周期 sEMG 信号波形的幅值相似度[29] 与变化趋势相似度相结合的 sEMG 波形相似度评估方法.

　　2.2 基于 SOM 的多路 sEMG 时序特征编码方法

　　为识别经信号处理后的 3 路 sEMG 信号序列 VR、VL、SM 的变化特点与运动过程的内在联系, 需要对其进行信息融合及编码. 由于人体运动具有过程性的特点, 单路 sEMG 时序变化特征反映了相应的肌肉力在运动过程中的时序变化, 多路 sEMG 时序则可以反映各肌群协调工作特点; 因此, 本文采用时序编码的方式将 3 路 sEMG 时序信息融合为一维时序码字, 用于表述人体运动过程中多路 sEMG 信息的时序特征.

　　通常实际采集到的 sEMG 幅值精度较高、值域范围较大, 因此, 对各时间点的 sEMG 幅值进行直接编码是不现实的, 需要对幅值的分布状况进行划分. 自组织神经网络 (Self-organization mapping net, SOM) 是基于无监督学习方法的神经网络的一种重要类型, 它可以通过其输入样本学会检测其规律性和输入样本相互之间的关系, 从而自适应地调整网络, 使其以后的响应与输入样本相适应[30−31] . 二维科荷伦 (Kohonen) 网络是 SOM 的一种形式, 它可以模拟细胞聚类的功能柱, 使输入样本通过竞争学习后, 功能相同的输入靠得较近, 不同的分得比较开, 由此将一些无规则的输入自动排开, 最后, 使权矢量的分布与输入样本的概率密度分布相似. 因此, 二维科荷伦网络能够表现出输入矢量在二维平面上的分布特性, 使 3 路 sEMG 信号投影到二维数据空间的某个区域上, 完成数据融合, 并通过二维输出层节点的编码信息完成对三维特征的编码. 可见, 这种方法不需要人为地对各路 sEMG 的幅值范围进行划分便可反映其分布情况, 在多路 sEMG 编码问题上简便而实用. 根据研究目标, 本文采用图 7 所示的二维科荷伦网络拓扑结构, 共分为两层:

　　1) 输入层: 有 r 个输入矢量, 即 r 个节点; 此处输入向量为经过信号处理后的 VR、VL、SM 信息, 因此 r = 3.

　　2) 输出层: 也称为竞争层, 为公平反映 r = 3 路 sEMG 输入的影响并降低识别算法的复杂度, 输出层采用 3 × 3 二维方阵型网络组成, 可输出 9 种编码码字, 即输出节点数量为 s = 9, 记编码值集合 V = {v1, v2, · · · , vs}; 所有的输入节点到输出节点之间均有权值链接; 按照竞争网络的特点其激活函数为二值型函数, 为更加接近实际神经细胞的特点, 兴奋区域采用获胜节点周围的六角形区域, 因此, 输出节点之间为六角形连接方式.

　　由此, 以 SOM 对 3 路滤波后的 sEMG 信号进行编码, 结合滑动巴特沃斯滤波法, 采用与滤波同频率的滑动时间窗口的方法即可实时生成 sEMG 时序特征序列.

　　3 基于隐马尔科夫模型的运动意图识别方法

　　人体运动是各种复杂随机过程的有机组合, 各运动过程具有过程性, 可由时序信息描述其特征, 进而达到对不同运动过程识别的目的. 隐马尔科夫模型 (Hidden Markov model, HMM) 具有很好的统计特性, 是一种根据随机过程的特点建立描述其时序信息特征的概率模型[32−33] , 目前已成功地应用于手势、语音、手写识别等领域[33] . 本文使用一阶 HMM 建立描述人体下肢运动过程 (踏车、行走及椭圆运动) 中多路 sEMG 信息的时序特征的数学模型, 并以此来对运动模式进行分类.

　　一阶 HMM 是一个双重随机过程, 此处, 它由一个描述下肢运动过程中各运动状态转移关系的、不能被 sEMG 传感器检测到的一阶马尔科夫链, 和一个用于描述运动状态与多路 sEMG 信息之间关系的一般随机过程组成, 即: 通过观测到的 sEMG 时序信息间接地分析隐含运动状态的转移. 结合本文研究, 定义人体下肢踏车、行走及椭圆运动过程中通过多路 sEMG 信息采集、处理后得到的编码码字时间序列为长 L 的观测值序列 O = {O1, O2, · · · , Ot , · · · , OL}, 其对应的隐含运动状态序列为 Q = {Q1, Q2, · · · , Qt , · · · , QL}, 其中 Ot 为 t 时刻的编码值, Qt 为 t 时刻的运动状态. 描述踏车、行走及椭圆运动的一阶 HMM 均可由一个 5 元组 λ = {M, N, π, A, B} 来表述, 如表 3 所示.

　　记描述踏车、行走及椭圆运动过程中多路 sEMG 信息时序特征的 HMM 分别为 λT 、λW 、λO. 以踏车运动为例, 以给定的该运动过程中多路 sEMG 时序特征信息作为 HMM 的观测值序列 O, 通过调整 λT 的参数使 P(O|λT ) 局部极大可训练得描述踏车运动过程多路 sEMG 时序特征的 HMM λT . 同理可得 λW 和 λO. 实际应用中, 对于给定的模型 λT 、λW 及 λO, 计算各自产生某个观测值序列 O (待识别运动过程中多路 sEMG 时序特征信息) 的输出概率 P(O|λT )、P(O|λW ) 及 P(O|λO), 便可估计 O 与模型 λT 、λW 、λO 的匹配程度, 基于最大似然估计法做出识别决策.

　　本文使用 Baum-Welch (B-W) 算法[32] 实现模型的训练和观测值输出概率的计算以提高时效性,降低计算复杂度; B-W 算法是一种广义期望最大化算法, 其核心思想是通过递归过程更新参数权重 (π, A, B), 进而得到可以更符合训练样本序列的模型参数. 为计算序列 O 在模型 λ 上的输出概率 P(O|λ) 需要定义两个变量:

　　1) 前向变量 αt(i): 给定参数 λ, t 时刻人体运动处于状态 ui , 此时已观测到部分经信号处理及编码后的多路 sEMG 信息特征序列 {O1, O2, · · · , Ot} 的概率, αt(i) = P(O1, O2, · · · , Ot , Qt = ui |λ);

　　2) 后 向 变 量 βt(i): 给 定 参 数 λ, t 时 刻人 体 运 动 处 于 状 态 ui , t 时 刻 之 后 观 测 到 的经 信 号 处 理 及 编 码 后 的 多 路 sEMG 信 息 特征 序 列 {Ot+1, Ot+2, · · · , OL} 的 概 率, βt(i) = P(Ot+1, Ot+2, · · · , OL, Qt = ui |λ).——论文作者：佟丽娜 1, 2 侯增广 1 彭 亮 1 王卫群 1 陈翼雄 1 谭 民 1