基于盲辨识模型和极限学习机的下肢假肢路况识别

2019-12-23刘作军宋寅卯吴青娥

振动与冲击 2019年23期

刘磊，杨鹏，刘作军，宋寅卯，吴青娥

(1.郑州轻工业大学建筑环境工程学院,郑州 450002； 2.河北工业大学人工智能与数据科学学院,天津 300130)

截止到 2014 年底, 我国 60 岁以上老龄人口已达2.12亿。随着交通事故、自然灾害等意外伤害的频繁发生以及脑血管疾病、糖尿病和骨关节病等慢性疾病的不断蔓延，根据世界卫生组织报告的最新结果显示，截肢患者达到 1725万, 下肢截肢约占肢体截肢总人数的51%。研发高性能智能下肢假肢是代偿截肢者基本运动功能的重要途径。20世纪90年代，智能假肢开始兴起并取得了一定的研究成果。德国奥托博克公司的 Genium 智能仿生假肢、冰岛奥索公司的锐欧仿生磁控膝关节(REHO KNEE)、英国英中耐公司的艾伦 (Elan) 仿生电子踝脚等智能假肢产品和样机逐步升级优化[1]。我国在假肢研究方面起步较晚，金德闻等[2]研发了一种六连杆假肢膝关节。台湾德林假肢有限公司推出三款多连杆大腿智能假肢产品。Chen等[3]设计出基于动态行走机理的“机器人假肢”。假肢控制首先是能够准确识别截肢者的路况，目前常用的路况识别信号源有运动学信息和肌电信号(Electromyography，EMG)。由于运动学信息需要某些特定的传感器，如足底压力、膝关节角度会在不同的路况中发生周期性变化，通过这些传感器获取的信号可以判断路况信息，但传感器在实际使用期间会磨损或者产生异位等现象造成路况识别准确率较低。肌电信号是产生肌肉力的电信号根源，是更接近原始运动意图的信号，通过EMG可以在未做出动作前获取主动运动意图，相对于采集足底压力、速度等运动力学信息的传统方法具有明显优势。因此本文选择EMG作为路况识别信息源。

肌电信号假肢路况识别研究中，特征提取和分类识别是关键的2个环节，特征提取方法主要有时域法、时域-频域法、频域法、非线性特征法等[4]。上述特征提取未考虑肌电信号源信号未知的特点，为了克服这一缺点，本文建立辅助变量盲辨识模型，并将双通道肌电信号的线性组合作为辅助变量克服了传统辅助变量必须为已知源信号的难题，采用模型系数作为EMG特征。高性能的分类器可以确保其控制的准确性。肌电信号分类方法主要包括BP神经网络、简约支持向量机、极限学习机、隐马尔科夫模型分类器、K 最近邻分类器等方法。吴剑锋[13]提出的基于表面肌电信号的简约支持向量机可识别六种路况，平均识别率达到 90%。佟丽娜等[5]使用隐马尔科夫模型分类器识别方法识别起立、坐下等8种路况，分类精度为91.67%。Huang等[6]同时使用了判别分析和支持向量机的分类识别方法识别平地行走等5种路况，识别率提高到 91.25%。李影[7]使用小波包能量作为特征，采用ELM分类识别路况，实验结果表明ELM方法优于BP神经网络、简约支持向量机、隐马尔科夫模型等方法，但是ELM方法中输入层到隐含层节点的权值和偏置都是根据样本分布随机赋值，如果选取不合适会影响识别结果。烟花算法(Fireworks Algorithm，FA)是一种模拟烟花爆炸过程的全局优化算法，为了克服这一缺点，本文采用烟花算法优化极限学习机进行路况分类。

本文首先将双通道肌电信号的线性组合作为辅助变量，将模型系数作为信号的特征，并与常用EMG特征比较，最后采用烟花算法优化极限学习机进行路况分类，并与ELM、BP等方法进行对比，实验结果表明，该方法对平地行走、上楼、下楼、上坡、下坡、跑步进行有效识别，识别准确率、鲁棒性优于ELM等方法。

1 实验方法

1.1 辅助变量盲辨识模型

本文根据肌电信号的产生机理[8]，采用2块肌肉(股内侧肌、股直肌)肌电信号，对肌电信号建立单输入多输出的有限冲击响应(Finite Impulse Response，FIR)模型。本文建立的FIR模型是将肌电信号产生的整个过程等效为一个FIR系统，系统的输入为大脑皮层中控制运动部分的运动神经元发放的电脉冲序列，电脉冲经脊髓、运动神经传给肌肉的整个传导过程等效为系统的一个传递函数H，输出信号为实验采集的EMG。图1为双通道FIR模型结构图。

图1 双通道FIR模型

图1中：u为运动神经元引发的激励信号，H(1)和H(2)为两路肌纤维传导通道的传递函数；x为系统输出；N为观测噪声；y为实际采集的肌电信号。

不同路况会引发不同的肌肉收缩模式，这个差异是通过传递函数H(1)和H(2)来体现的，本文将辨识传递函数的系数作为信号的特征。

图1中，输出向量y(t)=[y1(t),y2(t),…,ym(t)]T，输入向量u为未知向量，N(t)=[N1(t),N2(t),…,Nm(t)]T是互不关联的随机噪声向量，传递函数H形式

H(i)(z)=bi(0)+bi(1)z-1+bi(2)z-2+…+

bi(ni)z-ni

(1)

bi(l)是需要辨识的参数。

由图1可得式(2)：

y(i)=H(i)(z)u(t)+N(i)(t),i=1,2，…

(2)

因为u是未知的，因此，对于任意非零的常数α，都有以下关系成立：

(3)

即式(2)和(3)就会有相同的输出，为了克服传递函数H都不唯一的缺点。本文参考文[9]，先对传递函数H进行归一化处理，假设传递函数均为首1多项式，即bi(0)=1,i=1,2。

由模型可得

y(1)(t)-N(1)(t)=H(1)(z)u(t)，

y(2)(t)-N(2)(t)=H(2)(z)u(t)

(4)

式(4)两端相除并且十字相乘后可将式(4)化成差分方程的形式

(5)

定义向量

λT=[-y(2)(t-1),-y(2)(t-2),…，

-y(2)(t-n1),-y(1)(t-1),-y(1)(t-2),…，

-y(1)(t-n2)]T

(6)

σ=[b1(1),b1(2),…,b1(n1),b2(1),b2(2),…,

b2(n2)]T

(7)

(8)

则式(4)可以写成

y(2)(t)-y(1)(t)=λTσ+e(t)

(9)

式中:噪声e(t)无法确定其模型，本文采用辅助变量盲辨识的方法，该方法不用考虑噪声的具体形式，但是，该系统的输入变量是未知的，不符合传统辅助变量为输入变量的要求。辅助变量盲辨识方法的关键就在于构造辅助变量，且辅助变量要满足一定的限制条件：所选择的辅助变量要与噪声信号N(1)(t)和N(2)(t)无关，但是与y(1)(t)和y(2)(t)相关。

本文采用双通道信号的线性组合构造辅助变量。设

则式(9)可写成

Y=Aσ+E

(10)

这里n为数据的长度，若存在一个与A同维数的矩阵F使得下式成立

(11)

式中:Q是非奇异矩阵，在式(10)两端同时乘以，便可得到的无偏估计为

(12)

本文所构造的辅助变量形式为：

(13)

1.2 烟花极限学习机分类模型

文献[10-11]介绍了ELM分类算法，FA主要包括爆炸算子、变异操作。式(14)给出了初始化每个火花个数计算公式。

(14)

式中：ωi指第i个烟花产生的火花数量，i∈(1,N),m是常数，用来限制火花数量，f(xi)表示个体xi的适应度值，ymax是当前种群中适应度值最差个体的适应度值。ε是一个极小的常数。

FA采用爆炸烟花的方式多次迭代选择出理想的输入权值和隐含层偏差矩阵构成ELM分类器。

烟花爆炸范围的计算公式如下：

(15)

式中：Si表示第i个烟花的爆炸范围；Z表示最大的爆炸幅度。ymax是当前种群中适应度值最好个体的适应度值。f(xi)与ε与式(14)相同。

位移操作是对烟花的每一维进行位移，表达式为

(16)

式中：random(0,Ai)表示在幅度Ai内生成的均匀随机数。

高斯变异提高了算法的鲁棒性，表示第i个个体在第k维上的位置，计算方法如下

(17)

式中：g是服从均值为1，方差为1的高斯分布的随机数，g=N(1,1)。

2 基于盲辨识模型和极限学习机假肢路况识别

采用烟花算法进化极限学习机的下肢假肢路况识别方法首先选取下肢表面肌电信号作为信息源，提出肌电信号与运动关系的辅助变量盲辨识模型，将模型系数作为信号特征，分类器ELM算法中随机产生的输入权值和隐含层偏差只有少部分是比较优越的，为了克服该问题，本文提出了烟花算法进化极限学习机算法。图2给出了算法流程图，算法具体步骤如下：

(1) 按照1.1节方法提取路况特征值，特征值的维度为d，p为隐含层节点数，训练样本的均方根误差(Root Mean Square error，RMSE)作为烟花算法的适应度函数。

(2) 通过实验比较后选择合适的激活函数。

(3) 初始化N个烟花，确定烟花位置Hi，根据适应度函数计算出最优适应度值B。

(4) 根据式(14)计算出每一个烟花产生子代火花的个数wi(i=1,2,…,N)，i=1。

(5) 根据式(15)计算出每一个烟花的爆炸范围Gi。

(6) 根据式(16)对子代火花进行位移变异，对从k维中随机选择一定维度进行维度变异，在变异后的火花通过适应度函数选择最优的火花Qi。

(7) 随机选择m个烟花按照式(17)高斯变异，依据评价函数选出最优的火花Mi。

(8) 从Hi、Qi、Mi中选择最优的N个火花位置作为下一代火花爆炸位置。

(9)i=i+1，判断i是否达到最大迭代次数，不成立跳转到(3)。

(10) 得到最优火花位置，即最优的输入权值矩阵和隐含层偏差构建ELM分类器，分类6种路况。

图2 算法流程图

3 实验结果及分析

本文主要是对上楼、下楼、上坡、下坡、平地行走、跑步6种路况识别。受试者为10位假肢穿戴者，表1记录了受试者的年龄、残肢长度比例、假肢型号以及穿戴时间等信息。本文参考文献[4-5]方法选定股内侧肌和股直肌作为肌肉群，楼梯有10级台阶，台阶高度为广泛应用的150 mm；斜坡是无障碍通道的10°斜坡；平地实验在宽阔的室外走廊中进行，跑步在跑步机上完成。采用Delsys公司的Trigno肌电采集设备。该系统具有16个无线传感器，Trigno肌电采集设备采样频率为1 000 Hz。图3为平地行走实验，图4为上坡实验，为了避免速度对路况识别结果的影响，受试者所有路况均匀速行走。在实时处理表面肌电信号的过程中，本文使用300 ms小段数据来进行表面肌电信号的分析，后文做了实验证明。10位受试者分别采集了6种路况的EMG数据各30组。将每种路况的1 800组数据900份用来训练，900份用来测试。后文研究了训练样本所占比例不同情况下算法的识别准确率。表2给出了路况分类期望输出。

表1 受试者基本资料

表2 路况分类期望输出

表3给出了不同传递函数阶数ELM分类器路况识别准确率，当传递函数的阶数n1=2，n2=5时，路况识别准确率最高。表4为其中1位受试者6种路况下肌电信号通过辅助变量盲辨识方法辨识出的FIR模型系数。

表3 不同传递函数阶数下ELM分类器路况识别准确率

Tab.3 Road condition recognition accuracy of ELM classifier under different transfer function orders

n1=2n1=3n1=4n1=5n1=6n2=282.7587.2881.5579.2677.90n2=379.1581.3482.6488.5180.08n2=488.9590.1684.8285.2783.85n2=591.7682.4989.0990.8582.16n2=687.8388.7290.1684.8281.92

表4 不同路况下的FIR模型系数

在肌电信号路况识别问题的研究中，许多的研究方法是提取时域、频域和时频域特征。常用的时域特征有偏度(A1)、绝对值积分(A2)、方差(A3)、直方图(A4)、峭度(A5)、过零点数(A6)、均方根值(A7)、4阶自回归系数(A8)等。表面肌电信号的频域分析主要借助于功率谱分析方法,常用的频域特征有平均功率频率(A9)、中值频率(A10)和功率谱比值(A11)等，常用的时频域特征有小波变换系数(A12)、小波系数绝对值最大值(A13)。为了验证本文所提出的FIR模型系数特征用于路况识别的可行性，将这些单一特征值作为ELM的输入，再进行分类。从图5中可以看出：时频域特征识别准确率高于频域特征，频域特征高于时域特征。从方差角度来看，频域特征值显示了更好的稳定性，所以采用频域特征进行路况识别稳定性优于时频特征和时域特征。FIR模型系数(A14)特征表现出最好的分类准确率，FIR模型系数特征更能表征不同路况。

图5 采用不同的单一特征时的平均识别准确率和方差对比

Fig.5 Average recognition accuracy and variance comparison when different single features are used

表5 基于FIR模型系数特征的路况识别的混淆矩阵

表5给出了基于FIR模型系数特征的路况识别的混淆矩阵,L1、L2、L3、L4、L5、L6分别表示平地行走、跑步、上楼、下楼、上坡、下坡。由表5的混淆矩阵不难看出,平地行走和跑步路况比较容易识别,两者的识别准确率都超过了93.28%。其次为上坡和下坡2种路况，上坡和下坡的识别准确率均达到了91.57%。FIR模型系数类区分性更好，容易区分。较大差异，类间分离性较好，更具有代表性。上楼和下楼2种路况涉及多块相邻肌群的收缩，其识别准确率相对较低。上坡易误判为平地行走或下坡，产生误判主要是因为上坡和平地行走、下坡路况对应的某些肌肉所起的作用具有相似性。

本文建立的烟花极限学习机(FA-ELM)的识别模型中，该算法的参数设置设置为m=72,Z=30。常用的激活函数为Sine、Sigmoidal、Hardlim函数，图6为3种激活函数下6种路况的平均识别准确率的实验结果，Sigmoidal函数相对于Sine与Hardlim函数识别准确率高，波动小，因此本文选用sigmond函数。本文以 RMSE 作为误差的衡量标准，图7给出了ELM和FA-ELM训练误差和测试误差随隐含层节点数目变化的曲线图。

图6 3种激活函数下6种路况的平均识别准确率

Fig.6 Average recognition accuracy of 6 road conditions under 3 activation function

图7 误差随隐含层节点数目变化曲线

从图7可知，FA-ELM在隐含层节点数为10的时候RMSE达到最小值，而ELM为17时RMSE达到最小值。当隐含层节点数逐渐增加到50时，ELM的 RMSE 先减小后逐渐稳定。而FWAELM在隐含层节点数为10时达到最小，到50的过程中基本是稳定的。

设FA-ELM的隐含层节点数为10，迭代次数i=1，3，…，29，图8给出了FA-ELM算法中测试误差与训练误差随迭代次数的变化曲线。从图8可知，随着迭代次数增加，FA-ELM的测试误差和训练误差依次变小，迭代次数为19时，训练误差基本保持不变；测试误差达到最小值。

图8 迭代次数对烟花极限学习机误差的影响

迭代次数越多直接导致训练时间延长，为了权衡迭代次数和训练时间，本文设置FA-ELM的最优迭代次数为19，FA-ELM的隐含层节点p=10。

在训练分类器时，训练样本大小影响到识别准确率，为了分析训练样本数对识别准确率的影响。在900组训练集中，选择60组作为验证集，然后从余下的840组中按照抽取比例在10%～90%之间变化作为训练集样本，最后采用验证集检验算法的性能。

图9给出了训练样本数对BP、ELM、FA-ELM3种算法识别准确率的影响。有标签训练样本占训练集比例为10%时，本文方法的分类准确率为84.94%，比ELM算法高了 2.90%，比BP算法高了4.93%。随着有标签训练样本数量比例的上升，ELM等算法的分类精度逐渐提高，但是有标签训练样本占训练集比例为40%时，ELM分类精度有所下降之后又逐渐提高。当训练样本占比为90%时，FA-ELM的分类准确率为93.18%，仍然高于ELM算法(88.48%)，也高于BP算法(85.61%)。

图9 有标签训练所占比例不同情况下算法的识别准确率

图10给出了数据长度与识别准确率的关系，从图9的实验结果可以看出，数据长度的选择会对分类结果产生影响，对于数据长度小于200 ms，路况识别准确率都不高(小于88%)。当每次识别利用的数据长度100 ms时，随着数据长度的增大，识别准确率有所提高。当数据长度为300 ms达到最高的识别准确率，随着数据长度的继续增加，准确度并没有继续递增，而是略微的起伏变化。综合考虑系统的延迟时间和识别准确率，每次特征提取利用的数据长度为300 ms。

图10 识别准确率随数据长度变化曲线图

为了比较FA-ELM的性能，本文从识别准确率、训练时间和测试时间与BP网络进行对比。设置BP网络学习速率为0.1，最大迭代次数为2 000，BP神经网络收敛条件为输出均方差小于0.000 01。BP神经网络学习函数为trainlm。表6给出了基于FIR模型系数特征FA-ELM、ELM、BP 神经网络对6种路况的平均识别准确率和测试时间。从表 6可以看出，基于烟花算法进化ELM的分类方法比 ELM和BP神经网络分类平均准确率分别高1.82%和4%。采用ELM分类器上楼、下楼的识别准确率较低86.39%、82.18%。烟花算法进化的ELM对上楼、下楼的分类准确率分别为88.93%、88.26%。通过观察表6可以看出经过烟花算法优化ELM分类器的下肢假肢路况识别中，烟花算法进化ELM分类器的识别率高于没有经过参数优化的ELM分类器，除了上坡、下坡以外，每种路况达到了90%以上。实验结果表明，通过特征值的合理选择以及分类器的优化，路况识别准确率明显提升。

表6 FIR模型系数特征时路况识别情况

由表6可知，基于FIR模型系数单一特征，采用BP神经网络训练时间最长，采用烟花算法进化ELM与ELM时间相比，前者约是后者的2倍，基于FIR模型系数特征的烟花算法进化ELM测试时间约是基于FIR模型系数特征的ELM测试时间的三分之一，BP神经网络分类时间3种算法中耗时最长。

本文取最优迭代次数为19，FA-ELM的隐含层节点为10，做10次训练测试实验，测试结果如表7所示。利用烟花算法训练出来的极限学习机分类器的训练精度与测试精度差别较小。为了直观反映烟花算法进化极限学习机的鲁棒性，表7中的ELM和FA-ELM测试精度的标准差分别为0.182 9、0.024 5，由此说明FA-ELM鲁棒性优于ELM。