基于串联弹性原理和增广上三角分解辨识算法的踝关节阻抗估计方法

徐 刚，向 馗，唐必伟，庞牧野

（武汉理工大学自动化学院，湖北武汉430070）

人类在与外界环境交互时可自主调节踝关节的“硬度”［1］，以从容应对外界环境的变化，这项能力被称为踝关节阻抗调节［2］。踝关节阻抗调节对维持身体平衡、吸收外界冲击和推进步态至关重要［3］。踝关节阻抗研究不仅能为人体下肢生理学研究提供线索，还能为仿生机器人控制策略研究提供重要依据。近年来，踝关节阻抗研究受到众多领域学者的关注［4］。

阻抗测量是研究踝关节阻抗调节的关键。在很多文献中，常将关节阻抗定义为关节位置到力矩的映射，以描述关节角度和关节力矩之间的动态关系［5］，一般用二阶模型描述关节阻抗［6］。通常情况下，学者通过对关节施加外部干扰（常见的干扰方式为随机激励法和脉冲法）来获取关节响应，并借助辨识方法获得关节阻抗模型，以研究关节阻抗调节机制。Tucker等［7］自主设计了一款可穿戴干扰膝关节装置。Shamaei等［8］和 Tucker［9］利用可穿戴干扰膝关节装置分别测量了膝关节准刚度和膝关节关节阻抗。Yagi等［10］基于弹性绳快速伸曲特性设计了一款便携式踝关节阻抗测量装置（impedance detection device for ankle joint，IDDAJ），该装置采用最小二乘法来估计人体在准静态条件下的踝关节阻抗，具有测量快速和方便携带等优势，但其驱动装置最大输出扭矩有限，无法提供大干扰力矩，且忽略了小腿肌肉激活程度对踝关节阻抗的影响。Hassan等［11］设计了一种穿戴式阻抗测量装置，用于测量人体在小腿肌肉放松状态和激活状态下的踝关节阻抗，但该装置无法输出大干扰力矩，测量结果与其他文献相差较大。Lee等［12］自主设计了AnkleBot装置，用于测量踝关节阻抗，该装置配备了2个Kollmorgen电机，以提供恒定大力矩和高频干扰力矩，并采用相关分析法来辨识踝关节（2个自由度）阻抗模型参数，但该装置的成本过高。

为了准确研究踝关节阻抗特性和降低设备耗费，笔者及同课题组成员［13］基于串联弹性原理，在IDDAJ中对称安装了1对扭簧，以实现分离恒定大力矩和高频干扰力矩的目的，提高装置的响应能力。但通过相关试验发现，前期设计的IDDAJ仍存在以下问题：1）扭簧易产生塑性形变，导致装置失效；2）机械结构复杂且机械间隙大；3）机械结构复杂使得装置阻抗模型的阶次升高，采用机理法得到的二阶模型可能会导致辨识精度降低。针对上述问题，笔者拟对前期设计的IDDAJ及辨识算法进行改进，改进的IDDAJ利用力矩传感器代替扭簧来检测力矩信号，以提高力矩测量精度；利用定制的橡胶弹性体代替扭簧，以降低设计、安装难度及避免塑性形变引发的精度下降问题；通过1根法兰式中心轴直接连接踏板和橡胶弹性体，以减小机械间隙；采用增广上三角分解辨识（augmented upper-diagonal decomposition identification，AUDI）算法［14］对采集的数据进行辨识处理以直接构建阻抗模型。最后利用所设计的IDDAJ对人体端坐状态下的踝关节阻抗进行研究。

1 踝关节阻抗测量装置设计

IDDAJ的机械结构如图1所示，包括踏板、调节螺柱、扰动电机、承载框架、橡胶弹性体和传感器等。踏板作为人体脚掌的承载平面，通过绷带与脚掌绑在一起；为了适应踝关节旋转中心到脚底的距离，通过调节螺柱微调踏板高度来使踏板旋转中心和踝关节旋转中心重合；扰动电机（型号为FAULHABER 3863A024C，额定功率为200 W）用于提供高频干扰力矩，可对踝关节施加随机扰动；橡胶弹性体（直径为46 mm、高度为20 mm的圆柱体）的弹性系数kl=0.14 Nm/°，其左端通过法兰式中心轴直接与踏板相连，右端固定在承载框架上；传感器包括力矩传感器（型号为SUNRISE M2210M）和光电编码器（型号为HEDS-9040t），分别用于采集扰动电机输出的力矩和踝关节响应角度。根据串联弹性原理可知，当IDDAJ的扰动电机未施加干扰时，橡胶弹性体为踏板提供静态支撑力矩；当扰动电机施加干扰时，橡胶弹性体提供恒定大力矩，扰动电机输出高频干扰力矩，实现了恒定大力矩和高频干扰力矩的分离。IDDAJ的工作原理如图2所示，橡胶弹性体为踏板和脚掌提供静态支撑力矩Ts，以抵消踏板和脚掌重力造成的影响，使初始状态下踝关节始终保持在其中心位置处；扰动电机为踏板和脚掌提供高频干扰力矩τ。

图1 踝关节阻抗测量装置机械结构示意图Fig.1 Schematic diagram of mechanical structure of IDDAJ

图2 踝关节阻抗测量装置的工作原理Fig.2 Working principle of IDDAJ

IDDAJ的控制原理如图3所示。以STM32F4单片机为控制器，利用PID（proportion integration differentiation，比例积分微分）算法控制扰动电机输出的高频干扰力矩；以200 Hz的频率同步采集力矩传感器和光电编码器的数据，借助STM32F4单片机串口，以921 600 b/s的传输速率将同步采集的数据发送至上位机并保存。同时，借助肌电仪（型号为ELONXI EMG 100-Ch-Y-RA）实时采集小腿肌肉的表面肌电信号（surface electromyography，sEMG），并发送给上位机。

2 踝关节和装置的阻抗模型及辨识算法

2.1 踝关节阻抗模型

踝关节的动力学特性可用二阶微分方程来描述［10］：

式中：θ和τ分别表示踝关节响应角度和高频干扰力矩；t为时间；I、B、K分别表示惯性、阻尼和刚度分量。

在零初始条件下，对式（1）进行拉氏变换，得到以下传递函数：

图3 踝关节阻抗测量装置的控制原理Fig.3 Control principle of IDDAJ

式中：ZA(s)表示踝关节阻抗；YA(s)表示踝关节导纳；s表示拉氏变量。

2.2 IDDAJ阻抗模型

在踝关节约束条件下，由IDDAJ和踝关节构成的测量系统框图如图4所示，其中反馈通道模型为橡胶弹性体的弹性系数k1；前向通道模型为踝关节和IDDAJ踏板的并联阻抗模型；θ0为IDDAJ踏板的初始角度；θ1是踏板的响应角度。

图4 踝关节阻抗测量装置测量系统框图Fig.4 Block diagram of detection system of IDDAJ

整个测量系统的闭环传递函数为：

式中：YCL(s)、ZCL(s)分别表示测量系统的导纳和阻抗；YA+I(s)、ZA+I(s)分别表示踝关节和IDDAJ踏板的联合导纳和阻抗。

由于IDDAJ的踏板和人体脚掌通过绷带绑定一起，两者在运动过程中的位移相同。此时，踝关节和IDDAJ踏板的联合阻抗可表示为：

式中：ZA(s)、ZI(s)分别表示踝关节和IDDAJ踏板的阻抗。

在未受到人体踝关节约束的条件下，对IDDAJ进行空载实验，并根据实验结果辨识得到IDDAJ的阻抗(ZI(s)+k1)，然后在同等条件利用IDDAJ对踝关节的阻抗进行测量，得到测量系统的阻抗(ZA+I(s)+k1)，再利用式（4）计算得到踝关节的阻抗。

2.3 辨识算法

AUDI算法具有计算简单、便于实现以及对未知阶数模型辨识快速和准确等优点。通过上三角分解，可直接得到n0（不大于给定的最大模型阶次）阶模型的待辨识参数［14］，避免了复杂的建模过程。本文采用该方法对IDDAJ和踝关节的阻抗模型进行参数辨识。在辨识过程中，将m阶时域模型离散化，可得［14］：

式中：ai、bj分别表示阻抗模型输出和输入的系数，i和j的取值由模型阶次决定，本文取1～m；z(k)、u(k)分别表示阻抗模型的输出和输入；v(k)表示均值为零的白噪声；k表示时刻。

将式（5）转换为矩阵形式，可得：

其中：

式中：θm和hm分别表示由输入和输出的系数构成的模型参数向量和由输入和输出构成的数据向量。

定义增广数据向量φm为：

利用输入和输出数据，构造k时刻的信息压缩矩阵cm(k)，表示为：

为辨识式（7）中的模型参数向量，利用上三角分解法分解信息压缩矩阵cm，以构建AUDI算法。信息压缩矩阵cm可分解为：

其中：

式（5）所示离散模型经双线性变换［15］后变为连续时域模型，可直接辨识相关参数。

3 踝关节阻抗测量实验与结果分析

3.1 实验数据获取

挑选8名健康男性（年龄为23～25岁，体重为60～75 kg）参与踝关节阻抗测量实验，所有受试者近5年内没有小腿关节受伤史，在实验前均签署了知情同意书。

踝关节阻抗测量实验分3步：预备实验、训练实验以及测量实验。预备实验的目的是获取肌肉最大主动激活程度（maximum voluntary contraction，MVC），以对采集的原始sEMG进行归一化处理；训练实验的目的是训练受试者控制肌肉激活程度；在进行测量实验时，将整个IDDAJ置于水平面上，通过弹性绳将人体大腿悬挂于一定高度，高度随受试者腿长进行调整，人体脚掌与IDDAJ踏板通过绷带连接；扰动电机对踝关节施加高频干扰力矩，采用AUDI算法离线辨识得到IDDAJ和踝关节的阻抗模型参数。踝关节阻抗测量实验现场如图5所示。鉴于胫骨前肌（tibialis anterior，TA）、腓肠肌（gastrocnemius，GAS）和比目鱼肌（soleus，SOL）的激活程度对踝关节阻抗有重要影响，且其sEMG易于获取，本文对TA、SOL和GAS的sEMG进行采集，并发送到上位机实时显示，便于受试者按要求控制肌肉的激活程度。参照国际标准，设置电极片贴附位置和进行事前皮肤处理［16］。

图5 踝关节阻抗测量实验现场Fig.5 Experimental site of ankle joint impedance detection

为保证实验数据的充分性和有效性，对每个受试者进行多次实验，舍弃受试者未按要求控制肌肉激活程度得到的实验数据，选取其中8组有效数据用于踝关节阻抗估计。采用拟合优度R2、均方根误差ERMS（root mean square error，RMSE）和百分方差占比PVAF（percentage variance accounted for，PVAF）来评价阻抗模型［14］：R2越接近1，ERMS越小，说明所构建模型的拟合效果越好，PVAF越大表明所构建模型的可靠性越高［17］。

3.2 实验数据处理

sEMG的处理流程如下：全波整流、线性包络提取以及归一化处理。全波整流指移除sEMG的均值和信号整流，线性包络提取指对sEMG进行滤波处理和提取其包络线，归一化处理指将sEMG的值归一化［18］。采用零相移巴特沃兹滤波器对力矩传感器和光电编码器采集的原始数据进行滤波处理，并对滤波得到的数据进行去均值处理，保留其动态分量。

3.3 实验结果分析

3.3.1 IDDAJ阻抗模型参数辨识

采用AUDI算法辨识得到不同阶次的IDDAJ阻抗模型的损失函数值，如图6所示。各阶模型的评价指标如表1所示。

图6 不同阶次的IDDAJ阻抗模型的损失函数值Fig.6 Loss function values of IDDAJ impedance models with different orders

表1 不同阶次的IDDAJ阻抗模型的评价指标Table 1 Evaluation indexs of IDDAJ impedance models with different orders

根据图6可得，三阶及三阶以上的IDDAJ阻抗模型的损失函数值基本相同；由表1可知，二阶和四阶模型的 R2＜0.9，PVAF＜90%，三阶模型的 R2＞0.9，ERMS＜0.6，PVAF＞90%，故采用三阶模型更为合适。辨识得到的三阶IDDAJ阻抗模型采用双线性变换法［15］进行连续化处理，可得：

通过辨识得到IDDAJ阻抗模型中的弹性系数为0.27 Nm/°，与文献［13］的测量值0.388 Nm/°相比约减小了30%，说明本文装置的柔顺性更好。本文三阶IDDAJ阻抗模型的PVAF的平均值为91.53%，比文献［13］的87.66%增大了3.87%，说明本文三阶模型的可靠性更高。

3.3.2 放松状态下踝关节阻抗模型参数辨识

放松状态下某位受试者小腿上TA、SOL和GAS的sEMG如图7所示。由图可知，TA、SOL和GAS的sEMG基本小于0.01，且波动较小（见表2）。结果表明该受试者小腿处于放松状态，且按照实验要求进行阻抗测量。

图7 放松状态下小腿肌肉的表面肌电信号Fig.7 sEMG of calf muscle in relaxed state

表2 放松状态下小腿肌肉表面肌电信号的平均值和标准差Table 2 Mean value and standard deviation of sEMG of calf muscle in relaxed state

根据8位受试者的踝关节阻抗测量数据，采用AUDI算法辨识踝关节阻抗模型的参数，结果如表3所示。从表3中可以看出，大部分R2＞0.9，ERMS在0.4左右波动，说明二阶线性时不变模型的拟合效果较好，绝大多数的PVAF＞90%，说明该二阶线性时不变模型可以较准确地估计放松状态下踝关节的阻抗；个别R2＜0.9，对应的PVAF＜90%，表示实际阻抗模型中含有较显著的非二阶模型成分，且该成分对阻抗模型的参数辨识具有一定程度的影响；基于8号受试者踝关节阻抗测量数据得到的惯性分量I与基于其他受试者得到的相比存在较大差异，这可能是因为小腿或大腿重力的共同影响。在实验过程中，通过调整大腿悬挂高度以减少小腿和大腿重力的影响，但高度调整主要依靠受试者的主观感受，无法准确检测实际效果，从而导致根据不同受试者的踝关节阻抗测量数据辨识得到的惯性分量I具有较大的差异。本文得到的踝关节阻抗模型的阻尼和刚度分量与文献［12］结果处于同一数量级，但整体数值小于文献［12］的结果。

表3 放松状态下踝关节阻抗模型的参数辨识结果及其评价指标Table 3 Parameter identification results and evaluation indexes of ankle joint impedance model in relaxed state

3.3.3 不同肌肉激活程度下踝关节阻抗模型参数辨识

不同肌肉激活程度下某位受试者小腿上TA、SOL和GAS的sEMG如图8所示。结果表明，在不同肌肉激活程度下，受试者小腿上各肌肉的sEMG在小范围内波动，这说明单独控制小腿某一部分肌肉的激活程度存在一定难度。基于2位受试者的踝关节阻抗测量数据，采用AUDI算法辨识得到4种肌肉激活程度下踝关节阻抗模型，结果如表4所示。由表4可知，该模型的ERMS＜0.4，大多数 PVAF＞80%，少数R2和PVAF较小，一方面是因为尽管对IDDAJ进行了一些优化改进，但IDDAJ仍存在机械间隙和摩擦力，这些非线性因素将会对辨识结果产生较大的影响；另一方面是因为肌肉激活程度变化会导致阻抗模型的参数发生变化。结果表明踝关节阻抗与小腿肌肉激活程度密切相关，阻尼分量B、刚度分量K均随着肌肉激活程度的增大而增大，其中刚度分量K的增大幅度大于阻尼分量B，惯性分量I基本保持不变。

不同肌肉激活程度下踝关节阻抗模型的伯德图如图9所示。当输入频率高于4 Hz时，踝关节阻抗模型的幅度变化主要由惯性分量I引起；当输入频率低于4 Hz时，踝关节阻抗模型的幅度随肌肉激活程度的增大而增大，但其斜率几乎为0，此时刚度分量K的影响占主导作用。本文所得结果在变化趋势上与文献［19］结果一致，但其数值比文献［19］结果稍大。

4 结论

为提高阻抗模型的辨识精度，本文主要从IDDAJ的机械结构设计和辨识算法两个方面进行改进。由于IDDAJ存在机械间隙和摩擦力，机理法建模对装置的测量精度有一定影响。本文采用AUDI算法辨识IDDAJ阻抗模型的最优阶次。本文辨识得到的三阶IDDAJ阻抗模型的百分方差占比PVAF的平均值为91.53%，相比于文献［13］增大了3.87%。借助IDDAJ，采用AUDI算法，获得不同肌肉激活程度下踝关节阻抗模型。结果表明，踝关节阻抗与小腿肌肉激活程度密切相关，阻尼分量B、刚度分量K均随肌肉激活程度的增大而增大，其中刚度分量K的增大幅度大于阻尼分量K，惯性分量I基本保持不变。本文所得结果与文献［12］和［19］保持一致。

本文实验方法可为IDDAJ的设计提供参考，也可为仿生智能机器人阻抗控制提供理论指导，实验结果可为踝关节神经系统疾病检测提供理论依据。

图8 不同肌肉激活程度下小腿肌肉的表面肌电信号Fig.8 sEMG of calf muscle under different activation levels of muscle

表4 不同肌肉激活程度下踝关节阻抗模型的参数辨识结果及其评价指标Table 4 Parameter identification results and evaluation indexes of ankle joint impedance model under different activation levels of muscle

图9 不同肌肉激活程度下踝关节阻抗模型伯德图Fig.9 Ankle impedance model Bode diagram under different activation levels of muscle