外骨骼设备是用于增强、帮助以及恢复人体运动能力的机电设备^[1-4]. 外骨骼与人体运动的协调是外骨骼控制要解决的关键问题之一，通过正确预识别路况信息来判断人体运动意图，可以辅助外骨骼的控制决策. 目前的研究成果可以将路况识别的方法分为两大类：1）通过传感器采集人体运动过程中的速度、加速度^[5]等物理信息进行识别；2）通过外置传感器采集的距离信号^[6]识别. 由于物理信号伴随运动过程产生，对人体意图的识别具有一定的滞后性. 本文选择采用轻巧的、适于可穿戴设备装配的二维激光测距仪识别路况. 路永乐等^[7]采用支持向量机的方法判断人体运动模式，利用核映射的方法解决非线性分类问题，但处理多维输入数据时耗费的时间较长；神经网络分类^[8]方法需要对数据提取特征后再进行训练；K均值（K-means）聚类^[9]作为一种无监督学习方法，没有有效应用样本的标记信息.

膝关节外骨骼是外骨骼设备的一类分支，随着老龄化的加剧以及膝关节的脆弱特性，人们对膝关节助力外骨骼的需求不断增加. 美国Yobotics公司开发的膝关节外骨骼Roboknee^[10]，通过安装在大小腿之间的串联弹性驱动器为膝关节提供力矩，采用比例放大控制方法控制力拒. Rifai等^[11]基于拉格朗日法建立小腿-足部-外骨骼系统的动力学方程，采用有界控制策略对膝关节外骨骼进行控制. Chen等^[12]设计多模式弹性驱动器的膝关节外骨骼机械腿，对驱动器进行动力学建模后，采用运动状态机的控制策略进行控制. 韩亚丽等^[13]构建人机动力学模型，通过肌电信号对人体意图进行预判，采用基于导纳原理的等效惯量补偿控制方法，研究膝关节外骨骼机械腿. 以上基于动力学模型的控制方法已经较成熟和完善，然而由于人机系统的复杂性，其动力学模型和参数难以精确地确定，模型中不可避免地存在参数和结构误差. 无模型自适应控制（model free adaptive control，MFAC）作为数据驱动控制（data driven control，DDC）的一种方法，基于系统的输入输出数据对被控对象进行控制，已成功应用于化工精馏塔、通信拓扑跟踪、泊车系统、小型无人直升机控制^[14-16]等问题中.

本文针对膝关节外骨骼行走控制问题，判别人体运动意图，采用二维激光测距仪采集地面数据，选用学习向量量化^[17-18]（learning vector quantization，LVQ）方法，利用样本的类别标记信息，有监督地迭代学习，得到不同分类簇的原型向量，进而实现有效分类. 设计无模型自适应控制器，融合识别出的不同地形的先验力矩信息对膝关节外骨骼控制，搭建联合仿真平台，对控制方法进行验证.

膝关节外骨骼的控制系统结构如图1所示. 采用二维激光测距仪，扫描人体运动方向的地面距离进行识别，针对不同地形选择相应的期望角度轨迹，以协调人机系统在不同路况下的控制. 基于膝关节的输入输出角度建立动态线性化模型，依据输入准则函数设计无模型自适应控制器，融合先验力矩控制外骨骼. 先验力矩为不同路况下仿真得到的关节力矩，在控制器中添加一定比例的先验力矩，可以利用人体正常行走时力矩间的数据变化趋势，提高控制的准确性.

生活中常见的地形可以分为平地、上下楼梯、上下斜坡5种类型. 为了便于人体移动携带，选用图2所示的体积小、角度分辨率高的URG-04LX-UG01二维激光测距仪对地形识别，传感器的特征参数如表1所示.

传感器的扫描角度为0~240°，如图3（a）所示. 当把扫描仪垂直水平面放置在人体腰间时，扫描仪180°位置所扫描的对应点为人体站立点. 图3（b）中的 $A$点为安置于人体腰间的传感器位置点， $O$点为地面的人体站立点， ${h_1}$为人体腰部至足底的距离， ${\alpha _1}{\text{、}}{\alpha _2}$分别为传感器向前、向后扫描的角度， ${h_2}$、 ${h_3}$分别为传感器识别的人体运动方向及反向的地面距离. 地面识别距离的计算公式为

(1) $ {h_2} = {h_1}\tan {\alpha _1}{\text{，}} \\ {h_3} = {h_1} \tan {\alpha _2}{\text{，}} \\ BC = {h_2} + {h_3}. \\ $

测量受试者的 ${h_1}$为960 mm，选取 ${\alpha _1}{\text{、}}{\alpha _2}$分别为61°、25°，计算可得传感器识别的地面距离 $BC$为2 179 mm，满足地形预识别距离要求. 研究的数据采集范围为二维激光测距仪扫描的119°~205°内246个点的数据.

对采集后的数据进行预处理，剔除相邻点间距离相差较大的点，以提高数据的质量.

下肢运动数据采集部分使用VICON步态分析系统，在人体髋关节、膝关节、大腿、小腿、踝关节和足跟6个部位左、右侧各粘贴一个红外反光球作为标记点，如图4所示. 人体运动时步态分析系统的红外拍摄头会捕捉红外反光球的运动轨迹，通过计算可得下肢各关节的角度数据. 将计算得到的各关节角度数据作为期望轨迹.

通过对二维激光传感器测量范围的选取，将对5种地形的识别转化为对每帧246维数据的识别问题.

如图5所示为采集的5种地形数据. 可以看出，不同地形之间的距离特征相差明显，依据距离分类的方法可以有效分类. 在收集的数据 $({{x}},y)$中， ${{x}}$为n维的距离数据， $y \in (1,2,\cdots,5)$分别表示5种地形类别. 选用LVQ分类方法，对于n维的输入样本，LVQ通过有监督的学习得到一组n维的原型向量，每个原型向量表示一种类别. 在输入样本后，计算输入样本与每一个原型向量的欧式距离，选取与输入样本距离最近的原型向量所属类别为该输入样本的类别. LVQ的算法步骤如下.

1）预设原型向量个数为q，设置各原型向量的地形类别标记 ${t_1},{t_2},\cdots, {t_q}$，其中 $t \in (1,2,\cdots,5)$. 选取学习率 $\alpha \in (0,1)$，初始化一组原型向量 ${{{p}}_{{1}}},{{{p}}_{{2}}},\cdots,{{{p}}_{{q}}}$.

2）从采集数据中随机选取第 $j$个样本 $({{{x}}_{{j}}},{y_j})$，计算 ${{{x}}_{{j}}}$与每一个原型向量 ${{{p}}_{{i}}}$的距离，可得

(2) ${d_{ji}} = {\left\| {{{{x}}_{{j}}} - {{{p}}_{{i}}}} \right\|_2}; \;\;i = 1,2,\cdots,q.$

找到与 ${{{x}}_{{j}}}$距离最近的原型向量 ${{{p}}_{{k}}}$，其中 $k = \arg {\min _{i \in \{ 1,2,\cdots,q\} }}{d_{ji}}$.

3）判断选取的样本的类别标记是否与原型向量类别标记一致来更新原型向量 ${{{p}}_{{k}}}$. 若两者的类别标记相同，则令 ${{{p}}_{{k}}}$向 ${{{x}}_{{j}}}$的方向靠拢，否则更新 ${{{p}}_{{k}}}$远离 ${{{x}}_{{j}}}$的方向.

(3) $ \begin{array}{l} \left. \begin{array}{l} {{p}}' = {{{p}}_{{k}}} + \alpha ({{{x}}_{{j}}} - {{{p}}_{{k}}}){\text{，}}{y_j} = {t_k}{\text{;}}\; \\ {{p}}' = {{{p}}_{{k}}} - \alpha ({{{x}}_{{j}}} - {{{p}}_{{k}}}){\text{，}}{y_j} \ne {t_k}{\text{;}}\; \\ \end{array} \right\} \\ \;\; {{ p}_k} = {{ p}'}. \end{array}$

其中， ${y_j}$表示 ${{{x}}_{{j}}}$的地形类别标记， ${t_k}$表示 ${{{p}}_{{k}}}$的类别标记.

4）判断是否超过迭代次数. 若是，则输出原型向量；否则，返回2）.

经过学习训练得到一组原型向量后，即可实现对路况分类. 对于任意的输入样本 ${{x}}$，它将被划分至与该输入样本距离最近的原型向量所代表的类别中.

人机外骨骼模型如图6所示，其中二维激光测距仪放置在人体腰间位置. 采用拉格朗日动力学分析法，可得下肢动力学系统通用模型：

(4) ${{M}}({{\theta }}){{\ddot \theta }} + C({{\theta ,\dot \theta }}){{\dot \theta }} + G({{\theta }}) = {{\tau }} + {{{\tau }}_{\rm{h}}}.$

式中： ${{\theta }}$为关节角度向量， ${{M}}({{\theta }})$为正定惯量矩阵， ${{C}}({{\theta ,\dot \theta }})$为哥式力和离心力相关项矩阵， ${{G}}({{\theta }})$为重力相关项矩阵， ${{\tau }}$为外骨骼关节力矩， ${{{\tau }}_{\rm{h}}}$为人体作用于外骨骼的力矩. 其中，各项参数是考虑穿戴下肢外骨骼对人体模型影响后的参数.

式（4）可以整理为

(5) ${{\ddot \theta }} = {{M}}{({{\theta }})^{ - 1}}( - {{C}}({{\theta ,\dot \theta }}){{\dot \theta }} - {{G}}({{\theta }}) + {{\tau }} + {{{\tau }}_{\rm{h}}}).$

研究 $\theta $为膝关节角度时的控制情况，因而系统转换为单输入单输出系统. 定义外骨骼的膝关节力矩 $\tau $为输入 $u$，膝关节角度 $\theta $为输出 $y$. 其中，膝关节角度定义为股骨纵轴延长线与胫骨纵轴平行线的夹角，屈曲为正，伸展为负.

对于一个离散系统而言，在k时刻式（5）中的 $C\left( k \right)$、 $G\left( k \right)$、 $M\left( k \right)$是特定的值，人体作用于外骨骼的力矩 ${\tau _{\rm h}}(k)$是有界的，式（5）可以离散化为

(6) $\ddot y(k) = \frac{{ - C(k)\dot y(k) - G(k) + u(k) + {\tau _{\rm h}}(k)}}{{M(k)}}.$

当离散系统的采样时间足够小时，设采样时间为 $T$， $\dot y\left( k \right)$、 $\ddot y\left( k \right)$可以表示为

(7) $\dot y(k) = \frac{{y(k + 1) - y(k)}}{T},$

(8) $\ddot y(k) = \frac{{y(k + 2) - 2y(k + 1) + y(k)}}{{{T^2}}}.$

将 $\dot y\left( k \right)$、 $\ddot y\left( k \right)$代入式（6），可得

(9) $\begin{split} y(k + 1) =& 2y(k) - y(k - 1) + \frac{{ - C(k - 1)T(y(k) - y(k - 1))}}{{M(k - 1)}} + \\ &\frac{{{T^2}( - G(k - 1) + u(k - 1) + {\tau _{\rm h}}(k - 1))}}{{M(k - 1)}}.\end{split} $

式（9）满足如下条件：

1）对变量 $y(k)$和 $u(k)$存在连续偏导数；

2）满足广义Lipschitz条件^[19]，即对于任意时刻 ${k_1} \ne {k_2}$， ${k_1},{k_2} \geqslant 0$和 $u({k_1}) \ne u({k_2})$， $y({k_1}) \ne y({k_2})$时有

$\frac{{\left\| {y({k_1} + 1)} \right. - \left. {y({k_2} + 1)} \right\|}}{{\left\| {{{[y({k_1}),u({k_1})]}^{\rm T}} - } \right.\left. {{{[y({k_2}),u({k_2})]}^{\rm T}}} \right\|}} \leqslant b {\text{.}}$

式中： b为一个常数， $b > 0$. 为了方便下文叙述，定义2个相邻时刻的输出变化、输入变化分别为

(10) $\Delta y(k) = y(k) - y(k - 1),$

(11) $\Delta u(k) = u(k) - u(k - 1).$

对于满足条件1）、2）的系统，当 $\left\| {[y(k),\left. {u(k){]^{\rm T}}} \right\|} \right. $ $ \ne 0$时，必定存在伪梯度的时变参数向量 ${{\varPhi }}(k) = $ $[{\phi _1}(k),{\phi _2}(k)]$，使得系统（9）可以转化为全格式动态线性化模型^[20]：

(12) $\Delta y(k + 1) = {\phi _1}(k)\Delta y(k) + {\phi _2}(k)\Delta u(k).$

证明：系统（9）可以写成

(13) $y(k + 1) = f(y(k),y(k - 1),u(k),u(k - 1)).$

由式（10）、（13），可得

(14) $ \begin{split} \Delta y(k + 1) =& f(y(k),y(k - 1),u(k),u(k - 1)) -\\ & f(y(k - 1),y(k - 1),u(k - 1),u(k - 1)) +\\ & f(y(k - 1),y(k - 1),u(k - 1),u(k - 1))- \\ & f(y(k - 1),y(k - 2),u(k - 1),u(k - 2)). \end{split} $

令

(15) $ \begin{split} \varphi (k) = &f(y(k - 1),y(k - 1),u(k - 1),u(k - 1)) - \\ &f(y(k - 1),y(k - 2),u(k - 1),u(k - 2)){\text{，}} \end{split} $

由条件1）可知， $y(k)$和 $u(k)$存在连续偏导数，同时根据Cauchy微分中值定理，可以将式（14）转化为

(16) $\Delta y(k + 1) = \frac{{\partial f}}{{\partial y(k)}}\Delta y(k) + \frac{{\partial f}}{{\partial u(k)}}\Delta u(k) + \varphi (k).$

式中： $\displaystyle{{\partial f}}/{{\partial y\left( k \right)}}$、 $\displaystyle{{\partial f}}/{{\partial u\left( k \right)}}$为 $y(k + 1)$关于 $y(k)$和 $u(k)$的偏导数在 $[y(k),y(k - 1),u(k),u(k - 1)]$ 和 $ [y(k - 1),y(k - 1),$ $u(k - 1),u(k - 1)]$之间某一点处的值.

k时刻时，考虑含有变量 ${{\eta }}(k) = [{\eta _1}(k),{\eta _2}(k)]$的方程：

(17) $\varphi (k) = {{\eta }}(k){[\Delta y(k),\Delta u(k)]^{\rm T}}.$

由于 $[\Delta y(k),\Delta u(k)] \ne 0$，式（17）至少存在一个解. 令

(18) ${{\varPhi }}(k) = {{\eta }}(k) + \left[\frac{{\partial f}}{{\partial y(k)}},\frac{{\partial f}}{{\partial u(k)}}\right].$

式（14）可以写为

$ \begin{split} \Delta y(k + 1) = {{\varPhi }}(k){[\Delta y(k),\Delta u(k)]^{\rm T}}= {\phi _1}(k)\Delta y(k) + {\phi _2}(k)\Delta u(k). \end{split} $

式（12）证毕.

为了防止控制算法产生稳态的跟踪误差和过大的控制输入，采用如下控制输入准则函数^[21]设计控制器：

(19) $ \begin{split} J({u_{\rm MFA}}(k)) =& {({y^*}(k + 1) - y(k + 1))^2} + \\ & \lambda {({u_{\rm MFA}}(k) - {u_{\rm MFA}}(k - 1))^2}. \end{split} $

式中： $\lambda $为用来限制控制输入量的权重因子， ${u_{\rm MFA}}(k)$为无模型自适应控制输入， ${y^*}(k + 1)$为期望输出信号. 将式（12）代入式（19），对 ${u_{\rm MFA}}(k)$求导并令求导结果为零，可得控制输入为

(20) $ \begin{split} {u_{\rm MFA}}(k) =& {u_{\rm MFA}}(k - 1) + \frac{{\rho {\phi _2}(k)}}{{\lambda + {{\left| {{\phi _2}(k)} \right|}^2}}}({y^*}(k + 1) - \\ & y(k)) - \frac{{\rho {\phi _1}(k){\phi _2}(k)}}{{\lambda + {{\left| {{\phi _2}(k)} \right|}^2}}}\Delta y(k). \end{split} $

式中： $\rho $为步长因子，选取 $\rho $为 $(0,1]$是为了使控制算法具有更大的灵活性. 式（20）中，为了得到 ${u_{\rm MFA}}(k)$，须求解出伪梯度时变向量 ${{\varPhi }}(k) = [{\phi _1}(k),{\phi _2}(k)]$. 系统的数学模型是未知的， ${{\varPhi }}(k)$是一个时变向量，所以 ${{\varPhi }}(k)$的精确值很难获取. 采用修改后的投影估计算法^[22]估计伪梯度向量 ${{\varPhi }}(k)$：

(21) $ \begin{split} J({{\varPhi }}(k)) =& {\left| {\Delta y(k) - {{\varPhi }}(k){{[\Delta y(k - 1),\Delta u(k - 1)]}^{\rm T}}} \right|^2} + \\ & \mu {\left\| {{{\varPhi }}(k) - {{\hat \varPhi }}(k - 1)} \right\|^2}. \end{split} $

式中：μ为权重因子， $\mu > 0$； ${{\hat \varPhi }}(k)$为对 ${{\varPhi }}(k)$的估计值，对 ${{\varPhi }}(k)$进行微分并令求导的公式结果为零，可得

(22) $ \begin{split} {{\hat \phi }_1}(k) =& {{\hat \phi }_1}(k - 1) + \\ &\frac{{\eta \Delta y(k - 1)(\Delta y(k) - {{\hat \varPhi }}(k - 1){{[\Delta y(k - 1),\Delta u(k - 1)]}^{\rm T}})}}{{\mu + {{\left\| {[\Delta y(k - 1),\Delta u(k - 1)]} \right\|}^2}\;\;\;\;}} \end{split} {\text{，}}$

(23) $ \begin{split} {{\hat \phi }_2}(k) =& {{\hat \phi }_2}(k - 1) + \\ &\frac{{\eta \Delta u(k - 1)(\Delta y(k) - {{\hat \varPhi }}(k - 1){{[\Delta y(k - 1),\Delta u(k - 1)]}^{\rm T}})}}{{\mu + {{\left\| {[\Delta y(k - 1),\Delta u(k - 1)]} \right\|}^2}}}. \end{split} $

选取步长因子 $\eta \in (0,2]$是为了使控制算法具有更好的灵活性.在估算出 ${\hat \phi _1}(k)$和 ${\hat \phi _2}(k)$后，将 ${\hat \phi _1}(k)$、 ${\hat \phi _2}(k)$代入式（20），可以求解出 ${u_{\rm MFA}}(k)$.

将采集的各关节角度数据导入ADAMS仿真软件，利用驱动函数CUBSPL函数对人体模型进行驱动，通过仿真得到膝关节先验力矩. 如图7所示为仿真结束后的后处理模块中得到的膝关节力矩.

为了避免先验力矩对控制力矩的过多影响，为无模型自适应控制输入 ${u_{\rm MFA}}$和膝关节先验力矩 ${u_{\rm prior}}$分配权重，分别为 $1 - F$和 $F$，其中 $F \in [0,1]$. $F$越大，先验力矩对控制力矩的影响越大.

最终可得控制器形式如下：

(24) $u(k) = (1 - F){u_{\rm MFA}}(k) + F{u_{\rm prior}}(k).$

将平地、上下楼梯、上下斜坡地形分别标记为1~5类. 其中上楼梯地形包含平地至上楼梯、上楼梯以及上楼梯至平地3种形式，下楼梯、上、下斜坡为3种形式，即5种地形共13种形式. 实验采集5种地形情况下各2 400组数据，并将每种地形的前1 500组数据用于训练，后900组数据用于检验，即用7 500组数据进行训练，4 500组数据进行检验.

设置学习率 $\alpha $=0.01，原型向量个数为15，设置迭代次数为500次. 因为原型向量初始化具有随机性，重复训练300次，选取效果最好的一组原型向量作为最终的识别向量. 在训练结束后，通过测试数据检验LVQ的识别率，与支持向量机、K均值算法、奇异值识别3种方法进行比较. 比较结果如表2、图8所示.

外骨骼助行对地形识别的准确率和识别时间有较高的要求，通过上述比较发现，在LVQ、SVM、K-means和SVD 4种方法中，LVQ方法的识别率和识别时间综合结果更优. LVQ和K-means这2种方法是基于数据间的距离特征进行分类，但LVQ设置了15个原型向量，每种原型向量表示一种地形的聚类簇，5种地形一共有13种形式，LVQ的多原型向量可以学习得到多种特征，K-means算法只能选择分为5类，每一类的均值向量无法有效地代表同一地形的不同距离特征.

人体是一个复杂的系统，在人体运动特性和主要参数不变的前提下，对人体模型进行一定程度的简化. 将人体简化为头部、上躯干、左右手臂、左右大腿、左右小腿和左右脚10个体段，计算各体段的参数，如表3所示，在SolidWorks软件中建立相应的模型并导入至ADAMS软件中.

按照人体实际行走过程中的功能关系对模型添加关节副，其中头部和上躯干之间添加固定副固定；两手臂与躯干间添加旋转副；左右腿的髋关节、膝关节和踝关节添加旋转副. 人体的正常行走主要是在矢状面内的运动，故在人体上躯干处添加平面副，使得人体模型在矢状面内行走.

将采集的各关节角度数据依次导入ADAMS软件中为Spline函数曲线，并在下肢各关节处添加Motion驱动，通过样条函数CUBSPL将Spline曲线与关节关联起来. 采用关节角度驱动ADAMS中的人体模型行走，如图9所示为仿真的行走步态图.

在ADAMS的后处理模块中，将仿真模型的各关节角度导出，与实际角度进行对比，如图10~12所示. 经过计算可知，仿真的关节角度与采集的期望角度误差分别为：髋2.03%，膝1.63%，踝3.6%，求得的误差在允许范围内^[23]，表明了建立的人体模型的正确性.

ADAMS仿真软件可以将实际中的器件抽象成逻辑关系表达，使用2个“Link”杆来建立简化的下肢膝关节外骨骼，如图6所示. 分别在杆与上躯干和两杆间添加旋转副，模拟实际的转动情况. 现实中，膝关节处的动力是由电动机提供的，故添加膝关节力矩来模拟实际电机的输出力矩.

在MATLAB中建立系统的控制平台，在ADAMS中定义膝关节力矩为模型的输入变量，定义膝关节角度为输出变量. 将人体模型接口导入MATLAB中，在仿真时MATLAB计算的输入力矩和ADAMS的输出角度通过接口可以实时交互信息. 将控制算法写入MATLAB中的S-Function函数中作为控制器，ADAMS的模型接口作为被控对象进行仿真实验.

选取模拟路况为平地地形，使用MATLAB计算得到的控制力矩驱动膝关节，由期望关节角度驱动其他各关节. 如图13所示为通过仿真得到的人体行走步态图，其中控制器的各参数设置为 $\rho = 1,\eta = 1$， $\lambda = 0.8$， $\mu = 2$， $F = 0.1$，即在控制器中为先验力矩分配10%的权重.

选择仿真时间为人体正常行走的3个步态周期，如图14、15所示为得到的仿真效果和跟踪误差对比图. 由于模型初始化的原因，在仿真开始时刻膝关节角度从42°开始收敛逼近期望角度. 根据足部是否与地面接触，可以将步态周期划分为支撑期和摆动期两阶段. 仿真结果显示，在支撑期阶段输出角度与期望角度误差小于3°，基于MFAC的控制方法难以在摆动期阶段使得膝关节角度快速跟随期望角度的变化，融合先验力矩信息的MFAC控制方法能够在摆动期阶段实现良好的跟踪效果.

如图16所示为仿真先验力矩与控制力矩的对比. 可以看出，在摆动期内控制力矩与仿真先验力矩的差距较大，在支撑期内两者力矩较趋近.

对膝关节力矩的控制目的是使外骨骼可以最大化地对人体正常行走提供助力，减小行走时人体对外骨骼膝关节的力矩. 仿真结果表明，该控制算法能够仅在摆动期中期需要人体助力，验证了该控制算法的有效性.

针对外骨骼与人体运动的协调问题，本文提出利用二维激光测距仪预识别路况信息来判断人体的运动意图，采用LVQ识别方法，利用不同地形数据间的距离特征进行分类. 基于数据驱动的无模型自适应控制方法建立动态线性化模型设计控制器，融合先验力矩信息对膝关节外骨骼控制. 实验结果表明，平地、上下楼梯和上斜坡4种地形的识别率可以达到100%，对下斜坡地形的识别率达到95%. MATLAB和ADAMS的联合仿真结果验证了控制方法的有效性，融合先验力矩的MFAC控制方法在摆动期内可以更快速地跟随关节角度的变化. 下一步的工作重点是将路况识别和控制系统应用到实际系统中进行优化.