一种用于人体下肢表面肌电信号的辨识方法
【专利摘要】本发明涉及一种用于人体下肢表面肌电信号的辨识方法。本发明主要步骤为：实时采集对应肌肉块在活动动作刺激下的表面EMG信号；对采集的EMG信号进行预处理，获得消除伪迹信号后的EMG信号；使用离散小波变换方法对获得的EMG信号进行分解，分解后获得低频系数向量和高频系数列；采用时、频域结合的滤波方法，将获得的小波分量进行奇异值分解，并将分解得到的奇异值构成特征矩阵；采用支持向量机对特征样本进行训练，并生成支持向量机分类器用于对盲样本进行分类识别。本发明的有益效果为，相对于传统技术，本发明提出了一种新的下肢EMG信号预处理方法。
【专利说明】
一种用于人体下肢表面肌电信号的辨识方法
技术领域
[0001 ] 本发明设计用于人体下肢的生物电信号中的表面肌电信号(E1 e c tr omy ograph i c Signal，EMG)的伪迹消除、特征提取与辨识技术领域。
【背景技术】
[0002] 据统计数据显示，我国在21世纪初已正式迈入老龄化社会，且老龄化进程超过其 他国家，预计2020年我国老年人口会达到2.48亿，2050年将达到4亿。老年人及肢体残疾人 在人口结构中正在迅速扩大，上述人群的显著特点是其日常行为活动均需提供辅助帮助。 瘫痪是导致上述人群丧失运动能力的最常见原因之一，尤其是下身瘫痪，它涉及到肢体、躯 干的部分的或完全的功能丧失。目前临床专家普遍认为，肢体运动康复被认为是一个有效 的解决方案，其要求受影响的下肢肢体主动地参与积极锻炼。然而，对于肢体功能存在障碍 的人群，通常难以顺利完成例如站立、下蹲、行走等下肢肢体运动。因此，通过EMG信号的研 究能够帮助探索运动神经与肌肉组织的生物反馈机制，预测与感知肢体运动态势，评估老 年人、残疾人、以及亚健康人群的肌肉活动能力，为发展适用于老龄人和肢体残疾人士的康 复下肢运动辅助设备提供理论依据与应用基础。
[0003] 表面EMG信号是一种分析人体与日常生活活动相关的下肢肢体运动动作的生理信 息来源，该信号是通过电极引导、放大而从肌肉表面记录下来的神经肌肉系统活动时的非 平稳一维时间序列生物电信号，它能够反映肌肉强度和人的运动特征。人的运动意图通常 通过神经兴奋刺激肌肉细胞的收缩与舒张实现，由于不同的肢体运动中肌肉收缩的模式不 同，导致相应的表面肌电信号的特征也具有差异，在人的意图控制下产生的表面EMG信号能 较好地反映出肢体运动或动作特征，评估人体运动意图。人体运动意图的实现主要通过肢 体完成，由于下肢EMG信号相对于上肢更为复杂，同时受到更大的噪声干扰，因此对人体运 动的研究主要集中在上肢EMG信号，而对下肢EMG信号及其识别的研究需要进一步完善。目 前在EMG信号的识别中，常用的方法都是基于传统的分类或聚类算法，如支持向量机、神经 网络（Neural Network Algorithm，简称 "NNA"）、线性鉴别分析（Linear Discriminant Analysis，简称"LDA"）等。其中，LDA算法即可以识别单一的某种动作，也可以将多个动作加 上标签作为特殊的一类进行识别。在下肢EMG信号的研究上，2016年，John A. Spanias等人 使用LDA算法，研究了只用EMG信号进行分类、以及将EMG信号和器械传感器返回的其他类型 的数据一起进行分类的方法；2014年，AJ Young等人使用Sensor Time History的方法对 EMG信号进行分类，但该方法仅考虑了整个过程中信号的时间跨度。
【发明内容】
[0004] 本发明的目的，就是针对上述问题，提供一种用于人体下肢表面肌电信号的辨识 方法。
[0005] 通过实验分析，发现下肢EMG信号尺度通常非常微弱，同时由于硬件局限性和肢体 移动等原因，也十分容易受到工频、基线漂移和高斯白噪声干扰，因此实验所采集到的下肢 EMG信号伪迹严重，基于原始数据直接进行特征提取和分类识别是不现实的。传统伪迹消除 方法是对原始信号进行带通滤波处理(下肢EMG信号的主要成分集中分布在20~500Hz频段 上）。本发明中分别对工频干扰、基线漂移进行陷波滤波与低通滤波，根据本发明提出的工 频干扰噪声因子和基线漂移噪声因子分别对下肢EMG信号进行评估，当噪声的能量超出阈 值时，采用陷波和低通滤波，否则采用滤波前信号。而对于高斯白噪声，首先计算下肢EMG信 号的过零点数，通过过零点数对下肢EMG时域信号进行区间划分，再通过高斯白噪声因子逐 次对每一区间的噪声进行评估，当噪声的能量超出阈值时，采用滤波结果，否则采用滤波前 的信号。在表面EMG信号具有非平稳性、非线性特点的前提下，本发明将离散小波变换 (Discrete Wavelet Transform，DWT)应用到下肢EMG信号的特征提取中，它是将时域、频域 分析相结合对表面EMG信号的时间和频率所包含的信息进行分析。而传统的伪迹消除(如带 通滤波器）、特征提取(如傅立叶变换、时域和频域分析)方法仅仅独立地在时域或者频域中 分析数据，并将EMG信号视为平稳或者短时平稳信号进行处理，因此传统的方法并不能准确 地获取人体下肢肢体运动动作的EMG生理反应特征。奇异值分解（Singular Value DeC〇mp〇Siti〇n，SVD)是一种有效的代数特征提取方法，由于奇异值特征在描述信号数值上 比较稳定，且具有转置不变性、旋转不变形、位移不变性等重要性质，因此奇异值特征可以 作为信号的一种有效的代数特征描述。最后，本发明将离散小波变换(DWT)和奇异值分解 (SVD)相结合，基于新的下肢EMG信号预处理方法，进一步提出时、频域结合的滤波方法，根 据上述得到的特征矩阵，采用支持向量机(Support Vector Machine，SVM)的方法对信号进 行分类识别。
[0006] 本发明的技术方案为:一种用于人体下肢表面肌电信号的辨识方法，其特征在于， 包括以下步骤：
[0007] a.将一次性肌电电极粘放至下肢肌肉组织表皮上，实时采集对应肌肉块在活动动 作刺激下的表面EMG信号；
[0008] b.对步骤a中采集的EMG信号进行预处理，获得消除伪迹信号后的EMG信号;所述预 处理方法包括工频干扰滤波、基线漂移滤波和高斯白噪声滤波;本步骤中，工频干扰滤波、 基线漂移滤波和高斯白噪声滤波的顺序可以进行任意排列；
[0009] C.使用离散小波变换方法对步骤b获得的EMG信号进行分解，分解后获得低频系数 向量cAl和高频系数向量cDl;使用离散小波变换的方法对低频系数向量cAl进行分解，获得 低频系数向量cA2和高频系数cD2;继续重复使用离散小波变换的方法对低频系数向量进行 分解，直至获得低频系数向量cA5和5个高频系数列cDl，cD2，cD3，cD4，cD5为止；
[0010] d.采用时、频域结合的滤波方法，将步骤c中获得的小波分量进行奇异值分解，并 将分解得到的奇异值构成特征矩阵；
[0011] e.以步骤d中获得的特征矩阵为样本，采用支持向量机对特征样本进行训练，并生 成支持向量机分类器用于对盲样本进行分类识别。
[0012] 2、根据权利要求1所述的一种用于人体下肢表面肌电信号的辨识方法，其特征在 于，步骤b中所述预处理的具体方法包括：
[0013] bl.工频干扰滤波;具体方法为:将步骤a中采集得到的EMG信号c(t)作为原始输入 信号CPLI(t)进行陷波滤波，得到滤波结果为a(t)，定义工频干扰因子e PLI，则ePLI可通过如下 公式1计算：
[0014]
[0015] 其中，var是信号方差算子，用于计算时间序列的方差，通过工频干扰因子ePLI对滤 波结果进行修正如下公式2所示：
[0016]
[0017] 其中，sPLI(t)是滤除工频干扰噪声的最终结果，公式2表明如果工频干扰噪声的能 量占比超过原始信号能量的10%，则采用陷波滤波器的滤波结果；
[0018] b2.基线漂移滤波;具体方法为:将步骤bl获得的信号SPLI(t)作为原始输入信号CBW (t)进行低通滤波，滤波结果为d(t)，定义基线漂移因子为εΒν，则εΒνΡΤ通过如下公式3计算：
[0019]
[0020]根据基线漂移因子对滤波得到的基线漂移噪声d(t)进行修正得到b(t)，其表达 式如公式4所示：
[0021]
[0022]最后，得到去除基线漂移噪声后的信号SBw(t)如下公式5所示：
[0023] SBw(t) = CBw(t)-b(t) (公式 5);
[0024] b3.高斯白噪声滤波;具体方法为:将步骤b2得到的信号sBW(t)作为原始输入信号 CWCN(t)进行滤波;对时域信号CTCN(t)基于过零点进行区间划分，再通过阈值判断方法对信 号CTCN(t)进行高斯白噪声滤波;具体为:假设信号内包含j个过零点 21(1 = 1，2，...，」），则 对于区间Zi<tj<Zi+1内的信号有：
[0025]
[0026] 其中，T为信号c (t)的阈值，t m为区间内的极值点；阈值T由如下公式7和公式8确 定：
[0027]
[0028] o=median( | CwGN(t) | :t = l，2, · · ·，L)/0.6745 (公式8)
[0029 ]此处，σ为输入信号ctcn (t)的噪声水平;L为信号ctcn (t)的长度，特别地对于离散信 号L即为离散点点数;median为信号中位数算子，用于获取时间序列的中位数。
[0030] 本发明的有益效果为，相对于传统技术，本发明提出了一种新的下肢EMG信号预处 理方法，同时基于伪迹消除后的下肢EMG信号特点，首先提出结合DWT和SVD，对下肢典型运 动动作的EMG信号如行走支撑、行走行走、下蹲和站立等进行特征提取，最后通过SVM对所提 出的特征矩阵进行分类辨识。
【附图说明】
[0031] 图1为本发明的流程图；
[0032]图2为利用小波函数对信号进行5层小波分解过程的逻辑框图。
【具体实施方式】
[0033]下面结合附图，详细描述本发明的技术方案：
[0034]如图1所示，本发明用于人体下肢表面肌电信号的辨识方法，主要包括以下步骤： [0035]步骤1)下肢动作原始EMG信号获取：
[0036] 将一次性肌电电极粘放至下肢肌肉组织表皮上，实时采集对应肌肉块在活动动作 刺激下的表面EMG信号c(t);
[0037] 步骤2)原始信号预处理(伪迹消除方法）：
[0038] 对采集得到的原始EMG信号进行伪迹消除，其包括对工频干扰噪声、高斯白噪声、 基线漂移噪声的滤除。本步骤中的三种噪声的滤波没有固定的先后顺序，这里按工频干扰 噪声的滤波、高斯白噪声的滤波、基线漂移的滤波这一顺序进行描述;其次，在叙述上本文 中的"去噪"与"对噪声的滤波"含义相同。
[0039] a)工频干扰的去噪:将采集得到的表面EMG信号，c(t)，作为本次去噪过程的原始 输入信号，cPLI(t)，进行陷波滤波，得到滤波结果，a(t)，同时定义工频干扰因子，e PLI，如公 式(1)所示：
[0040] (1)
[0041] 其中，var是信号方差算子，用于计算时间序列的方差。通过工频干扰因子ePLI对滤 波结果进行修正如下：
[0042] (2)
[0043] 其中，SPLI(t)是滤除工频干扰噪声的最终结果。公式(2)表明如果工频干扰噪声的 能量占比超过原始信号能量的10%，则采用陷波滤波器的滤波结果。
[0044] b)基线漂移的去噪:将进行工频干扰滤波后得到的信号，sPLI(t)，作为本次去噪过 程的原始输入信号，CBW(t)，进行低通滤波，滤波结果为d(t)。类似地，定义基线漂移因子， ￡bw，如公式（3)所不：
[0045] … (3)
[0046] 根据基线漂移因子，eBw，对滤波得到的基线漂移噪声，d(t)，进行修正得到，b(t)， 其表达式如公式(4)所示：
[0047] (4)
[0048] 最后，得到去除基线漂移噪声后的信号SBw(t)如下：
[0049] SBw(t) = CBw(t)-b(t) (5)
[0050] C)高斯白噪声的去噪:将进行基线漂移噪声滤波后得到的信号SBW(t)作为本次去 噪过程的原始输入信号CTCN(t)进行滤波。对时域信号CTCN(t)基于过零点进行区间划分，再 通过阈值判断方法对信号C TCN(t)进行高斯白噪声滤波。假设，信号内包含j个过零点21(1 =
[0051] 1，2,. . .，j)，则对于区间Zi<tj<Zi+1内的信号有：
(6)[0052]其中，T为信号c (t)的阈值，t m为区间内的极值点。阈值T由公式(7)和(8)确定：[0053]
(7)
[0054] o=median( | cwGN(t) | : t = 1，2，· · ·，L)/0·6745 (8)
[0055]此处，σ为输入信号CWGN( t)的噪声水平;L为信号cwgn( t)的长度，特别地对于离散信 号L即为离散点点数;median为信号中位数算子，用于获取时间序列的中位数。
[0056] 步骤3)离散小波变换(DWT):
[0057] 小波变换是将时域分析和频域分析相结合形成的一种新的分析方法，反映的是表 面肌电信号在时间和频率两个维度上所呈现的变化，因此该方法相对于以上两种方法而言 理论上应具有一定的优势，能够充分利用表面肌电信号所包含的信息。小波分析是一种窗 口的大小固定、形状可变的时频局部化信号分析方法，即在低频部分具有较高的频率分辨 率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率。
[0058]给定一个长度为N的信号S，离散小波分解(DWT)最多可以把信号分解成log2N个频 率级。第一步分解开始于信号S，分解后分解系数由两部分组成:低频系数向量cAl和高频系 数向量cDl，前者也被称为近似(Approximation)成分，后者也被称为细节(Detail)成分。向 量cAl是由信号S与低通分解滤波器经过卷积运算得到的，向量cDl是由信号S与高通分解滤 波器经过卷积运算得到的。在下一步分解中，用同样的方法把低频系数cAl分成两部分，即 把上面的信号S用cAl代替，分解后返回尺度2的低频系数cA2和高频系数cD2,同时尺度1的 高频系数cDl保持不变，依此类推继续分解。最后得到一个低频系数列cA5和五个高频系数 列cDl，〇02,〇03,〇04,〇05，如图2所不。
[0059]基于对下肢运动的EMG信号的分解得到的分量，我们在对其进行带通滤波处理，使 得。01，。02，。03，。04^5的频带范围分别是256-512抱、128-256抱、64-128抱、32-64抱以及 16-32HZ。
[0060] 步骤4)奇异值分解(SVD):
[0061]根据步骤3中的方法，下肢运动的EMG信号可以通过DWT方法进行分解，从步骤3中 得到的一个低频系数列cA5和五个高频系数列cDl，cD2，cD3，cD4，cD5中，选取5个数据来构 成矩阵表达式，如选择:cD2，cD3，cD4，cD5和cA5构成矩阵，得到5 X I的矩阵表达如下，Mij = [d2i;d3i;d4i ;d5i;a5i]，i = l，2, · · ·，1，其中I为下肢运动动作样本数，d2i、d3i、d4i、d5i和 a5i对应DWT分解获得的系数，为了获取能在SVM分类中能够使用的特征值，这里将使用SVD 对Mij进行分解获取奇异值。
[0062] 对于M，NeCmXn，如果存在m阶酉矩阵U和η阶酉矩阵V使得等式UTMV=N成立，则认为 矩阵Μ和N酉等价。如果ME >矽，则矩阵MTM的特征值应满足以下关系：
[0063] λχ λ2 ^： · · · ^： λΓ λΓ+ι = · · · = λη=〇 (9)
[0064] 于是可得，矩阵Μ的奇异值为
[0065] 老
则一定存在m阶酉矩阵U和η阶酉矩阵V，使得矩阵Μ，矩阵U和
[0066] 矩阵V满足公式(10):
(10)
[0067] Σ =^38(0^02, · · .or)，其中〇i(i = l，2, · · ·，r是矩阵 Μ 的非零奇异值；公式（10) 等价于公式(11):
[0068] UU
[0069] 这里，
[0070] U=MMT (12)
[0071] V=MTM (13)
[0072] 其中，公式（11)也被称为矩阵Μ的奇异值分解(SVD)。
[0073] 步骤5)SVM分类：
[0074]为了对下肢运动进行分类，通过DWF分解得到的这些矩阵的奇异值，被用作下肢 EMG信号的时频域特征矩阵。对于被选定的数据段，矩阵由以下公式确定：
[0075]
[0076] 其中，i是每次实验(Trial)的数据的奇异值个数，j是EMG样本个数。
[0077] Vapnik等人创建了支持向量机(SVM)方法。SVM的基本思想是将数据映射到高维空 间中，并找出具有最大决策边缘的超平面。SVM公式所表征的风险最小化原则已被证明优于 其他传统的经验性的风险最小化原则。为了获得一个鲁棒性好的实时分类器，在本发明中， 用于SVM分类所需的特征值将从时域或频域两方面获取，或使用时频域相结合的分析方法。
[0078] 假设存在1个样本的训练向量^和相对应的分类标签7^(~71)，...，（~71)￡1^ X {-1，1}，则SVM本质上是求解二次最优化问题：
[0079]
[0080] yi(w · Xi+b)^Ki，Gi>0(i = l，2, · · ·，1) (15)
[0081] 其中，C是常数是松弛变量，w是权重向量，b是偏差值。对于实例X，其判别函数 如下所5·
[_] (]6)
[0083]其中，Ns是支持向量的个数，cti是正的拉格朗日乘子，G(x，Xi)是核函数。SVM中最重 要的规则参数也由训练集的十倍交叉验证过程所决定。
[0084]本专利采用线性内核函数K处理特征矩阵M^，支持向量机的分类函数f(x)可描述 为：
[_5] (17)
[0086] η是样本个数，幻是第j个样本，η是第j个样本的SVM输出，K是用于数据变换的线性 内核函数，A是对偶优化问题的拉格朗日乘子。
[0087]样本数据首先通过归一化处理（零均值+单位标准偏差处理），然后通过十倍交叉 验证技术评估分类辨识效果。本发明能够使人体下肢肢体运动动作的识别准确率得到有效 提1?。
【主权项】
1. 一种用于人体下肢表面肌电信号的辨识方法，其特征在于，包括以下步骤： a. 将一次性肌电电极粘放至下肢肌肉组织表皮上，实时采集对应肌肉块在活动动作刺 激下的表面EMG信号； b. 对步骤a中采集的EMG信号进行预处理，获得消除伪迹信号后的EMG信号;所述预处理 方法包括工频干扰滤波、基线漂移滤波和高斯白噪声滤波； c. 使用离散小波变换方法对步骤b获得的EMG信号进行分解，分解后获得低频系数向量 cAl和高频系数向量cDl;使用离散小波变换的方法对低频系数向量cAl进行分解，获得低频 系数向量cA2和高频系数cD2;继续重复使用离散小波变换的方法对低频系数向量进行分 解，直至获得低频系数向量cA5和5个高频系数列cDl，cD2，cD3，cD4，cD5为止； d. 采用时、频域结合的滤波方法，将步骤c中获得的小波分量进行奇异值分解，并将分 解得到的奇异值构成特征矩阵； e. 以步骤d中获得的特征矩阵为样本，采用支持向量机对特征样本进行训练，并生成支 持向量机分类器用于对盲样本进行分类识别。2. 根据权利要求1所述的一种用于人体下肢表面肌电信号的辨识方法，其特征在于，步 骤b中所述预处理的具体方法包括： bl.工频干扰滤波;具体方法为:将步骤a中采集得到的EMG信号c(t)作为原始输入信号 CPLI(t)进行陷波滤波，得到滤波结果为a(t)，定义工频干扰因子ePLI，则机河通过如下公式 1计算：其中，var是信号方差算子，用于计算时间序列的方差，通过工频干扰因子ePLI对滤波结 果进行修正如下公式2所示：其中，sPLI(t)是滤除工频干扰噪声的最终结果，公式2表明如果工频干扰噪声的能量占 比超过原始信号能量的10%，则采用陷波滤波器的滤波结果； b2.基线漂移滤波;具体方法为:将步骤bl获得的信号sPu(t)作为原始输入信号cBW(t) 进行低通滤波，滤波结果为d(t)，定义基线漂移因子为εΒν，则εΒνΡΤ通过如下公式3计算：根据基线漂移因子ε?对滤波得到的基线漂移噪声d(t)进行修正得到b(t)，其表达式如 公式4所示：最后，得到去除基线漂移噪声后的信号SBW(t)如下公式5所示： SBW(t) = CBW(t)-b(t)(公式5); b3.高斯白噪声滤波；具体方法为:将步骤b2得到的信号SBW(t)作为原始输入信号CTCN (t)进行滤波;对时域信号CTCN(t)基于过零点进行区间划分，再通过阈值判断方法对信号 CWCN⑴进行高斯白噪声滤波;具体为:假设信号内包含j个过零点21(1 = 1，2，...，」），则对 于区间Zi<tj<Zi+1内的信号有：其中，T为信号C(t)的阈值，U为区间内的极值点；阈值T由如下公式7和公式8确定：〇=mediarU I cwGN〈t川：t = U，· · ·，L)/0.6745(公式8) 此处，σ为输入信号ctcn( t)的噪声水平;L为信号ctcn( t)的长度，特别地对于离散信号L 即为离散点点数;median为信号中位数算子，用于获取时间序列的中位数。
【文档编号】A61B5/04GK105997064SQ201610326948
【公开日】2016年10月12日
【申请日】2016年5月17日
【发明人】张羿, 温悦欣, 张向刚, 秦开宇
【申请人】成都奥特为科技有限公司