非平整路面下被动行走机器人的步态动力学演化

Evolution of gait dynamics of passive walking robot on rough terrain

摘要

为了探索非平整路面环境对被动行走机器人的步态演化影响，采用正弦函数模型刻画路面的凹凸特性，并通过引入随机变化的路面幅度和频率来模拟真实行走环境下道路的起伏状态。通过运用分岔图和李雅普诺夫指数等数值分析工具，对比不同路面不平度对机器人步态稳定性的影响，并深入探讨全局分岔中的边界激变事件引发的机器人跌倒行为。此外，搭建机器人样机并进行ADAMS动力学仿真，验证了模型的行走能力。研究表明，在正弦波形的非平整路面上，机器人的步态为拟周期运动，并随系统参数的改变经环面倍化分岔逐渐演变为混沌；而在随机非平整路面上，路面对机器人步态的激励效应导致其偏离极限环轨迹，系统的稳定性下降；在上述行走环境中，产生于鞍结分岔的不稳定轨道引发的双边界激变事件是步态吸引子最终消失的重要原因。

Abstract

In order to study the excitation effect of gait displacement caused by rough terrain environment on passive walking robot， the sinusoidal function model is used to describe the concave and convex characteristics of the road surface， and the random variation of amplitude and frequency is introduced to simulate the unevenness in the real road environment. By using bifurcation diagram and Lyapunov exponent， the influence of road unevenness on the gait stability of robot is compared and analyzed， and the boundary crisis event in global bifurcation is deeply studied. In addition， a robot prototype is built and ADAMS walking simulation is performed to verify the walking capability of the model. It is shown that， on an uneven sinusoidal road surface， the robot's gait exhibits quasi-periodic motion and transitions into chaos with the change of system parameter， occurring via the torus multiplication bifurcation. The excitation effect of the robot gait caused by the rough terrain makes the robot deviate from the limit cycle trajectory， and this results in system degradation in terms of stability. In the walking environments above， the double boundary crisis event triggered by the unstable orbit generated by the saddle-node bifurcation is the main reason for the disappearance of the gait attractor.

Graphical abstract

关键词

机器人技术 / 被动行走 / 非平整路面 / 分岔 / 李雅普诺夫指数 / 双边界激变

Key words

robotics / passive walking / rough terrain / bifurcation / Lyapunov exponents / double boundary crisis

0 引 言

被动行走机器人无须外部控制和驱动系统，仅凭借自身重力就能够在斜面上稳定行走，其步态具有能量利用率高和行走自然的优点，为揭示类人行走机理和设计高效节能的双足机器人提供了崭新的视角。由于被动机器人的稳定步态在相空间中呈现为一个封闭的环，因此该行走方式也被称为“极限环行走”^［1］。在机器人的行走过程中，一个完整的步态由腿的连续摆动阶段和足地碰撞的离散事件两部分组成^［2］，所以被动行走机器人系统属于典型的脉冲混杂系统，其动力学特性异常复杂。

迄今为止，关于被动行走机器人的研究主要集中在平坦的斜坡环境下^［3-7］。Goswami等^［8］和Garcia等^［9］构建了罗盘模型和最简模型，揭示了步态失稳的分岔路径，即从倍周期分岔到混沌。此后，研究者们基于这两类基础模型，广泛研究了被动行走机器人的步态动力学特性。Wisse等^［10］在最简模型中添加上肢体，研究上肢对步态稳定性的影响规律；Narukawa等^［11］为罗盘模型增加了膝关节和上体躯干等结构，研究了此模型的步态动力学特性，并利用躯干驱动实现了平地上的稳定行走；柳宁^［12］和Kim等^［13］分别建立了圆弧足和平板足模型，分析了足部形状参数对机器人步态的影响；Wu等^［14］和Safartoobi等^［15］分别探究了髋关节和柔性腿对步态稳定性的影响。

然而，在实际步行环境中，路面的不平整性为行走步态带来了不可避免的外界干扰^［16］，但目前关于非平整路面导致机器人步态产生位移激励的研究还相对较少。Afshar等^［17］基于正弦函数构造了粗糙路面，分析了该路面对机器人步态稳定性的影响，揭示了路面粗糙度对步态的影响规律。但是，使用正弦函数构造的地面轮廓仍是结构化的路面，无法刻画实际的地表轮廓；Su和Dingwell^［18］研究了具有随机扰动的不规则路面对最简模型步态的影响，结果表明，较小扰动将导致步态偏离极限环，较大扰动则会增加步态的变异性。注意到，这些研究多集中于对步态的特征分析，而路面不平整性对步态的演化分析以及演化中的局部及全局分岔行为则鲜有深入研究。

综上所述，对被动行走机器人在非平整路面环境下的步态动力学演化行为研究是极其必要的。本文旨在利用正弦函数构建凹凸路面特征，并通过随机数描述路面不平度^［19］和起伏频率用于模拟非平整路面环境。首先，综合利用分岔图、李雅普诺夫指数等数值工具，对比分析路面不平度大小对机器人步态稳定性的影响；其次，重点研究了导致步态吸引子消失的双边界激变现象，以揭示路面不平整性对机器人摔倒事件的影响规律；最后，使用ADAMS软件对机器人的行走过程进行了动力学仿真。本文的研究结果不仅有助于深入地理解被动行走机器人在自然环境下的步态动力学特性，而且为仿人双足机器人的优化和设计提供了理论参考。

1 罗盘模型及非平整路面建模

1.1　罗盘模型

如图1所示，本文以罗盘模型为研究对象，该模型由两条相同的刚性直腿和一个髋关节组成。其中，腿质量为m，质心到髋关节的距离为a，质心到足部的距离为b，腿的长度为l（l=a+b），髋关节的质量为m_{\mathrm{H}}。{\theta }_{\mathrm{s}}和{\theta }_{\mathrm{n}\mathrm{s}}分别表示支撑腿和摆动腿与竖直方向的夹角，两角度均以逆时针方向为正。将机器人跨距角记为\mathrm{2}\alpha，斜坡倾角记为\phi。图1中各参数的取值如表1所示。对该模型作如下假设：①行走过程中髋关节处无阻尼、无摩擦；②忽略运动中位时的摆动腿擦地现象；③机器人的足地碰撞视为完全非弹性，满足角动量守恒定理；④足地碰撞瞬时完成，且支撑腿与地面间无相对滑动和弹起。

当为机器人设置合适的运动初始条件时，该机构即可依赖自身重力沿斜面向下保持稳定行走。一个完整的周期步态可分为单腿支撑的摆动阶段和足地碰撞的步态切换阶段，下面将分别建立这两个阶段的动力学方程。

1.1.1　摆动阶段

在此阶段，支撑腿绕支撑足做倒立摆运动，摆动腿则绕髋关节做单摆运动。由于仅有重力做功，故系统机械能守恒。因此，利用第二类拉格朗日方法可建立模型在此阶段的动力学方程为：

\mathit{M}(\mathit{\theta })\ddot{\mathit{\theta }}+\mathit{C}(\mathit{\theta },\dot{\mathit{\theta }})\dot{\mathit{\theta }}+\mathit{G}(\mathit{\theta })=\mathrm{0}（1）

式中：向量\mathit{\theta }=[{\theta }_{\mathrm{s}},{\theta }_{\mathrm{n}\mathrm{s}}{]}^{\mathrm{T}}，其余各矩阵的具体表达式如下所示：

\mathit{M}(\mathit{\theta })=

\left[\begin{array}{cc}m{b}^{\mathrm{2}}& -mbl\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{s}}-{\theta }_{\mathrm{n}\mathrm{s}})\\ -mbl\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{s}}-{\theta }_{\mathrm{n}\mathrm{s}})& m{a}^{\mathrm{2}}+(m_{\mathrm{H}}+m)l^{\mathrm{2}}\end{array}\right]，

\mathit{C}(\mathit{\theta },\dot{\mathit{\theta }})=

\left[\begin{array}{cc}\mathrm{0}& mlb\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{s}}-{\theta }_{\mathrm{n}\mathrm{s}}){\dot{\theta }}_{\mathrm{s}}\\ -mlb\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{s}}-{\theta }_{\mathrm{n}\mathrm{s}}){\dot{\theta }}_{\mathrm{n}\mathrm{s}}& \mathrm{0}\end{array}\right]，

\mathit{G}(\mathit{\theta })=\left[\begin{array}{c}mgb\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{n}\mathrm{s}}\\ -(ma+m_{\mathrm{H}}l+ml)g\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{s}}\end{array}\right]。

1.1.2　碰撞阶段

在碰撞前后，机器人的腿部位姿不变，角速度发生突变，整个系统关于摆动腿与地面接触点的角动量守恒，而支撑腿关于髋关节的角动量守恒。因此，根据角动量守恒定理建立此阶段的代数方程为：

{\mathit{Q}}^{-}(\alpha ){\dot{\mathit{\theta }}}^{-}={\mathit{Q}}^{+}(\alpha ){\dot{\mathit{\theta }}}^{+}（2）

式中：上标“-”和“+”分别为标识足地碰撞前、后系统的状态。

各矩阵的具体表达式如下所示：

{\mathit{Q}}^{-}(\alpha )=\left[\begin{array}{cc}-mab& -mab+(m_{\mathrm{H}}{l}^{\mathrm{2}}+\mathrm{2}mal)\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\alpha )\\ \mathrm{0}& -mab\end{array}\right]，

{\mathit{Q}}^{+}(\alpha )=

\left[\begin{array}{cc}m{b}^{\mathrm{2}}-mbl\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\alpha )& K\\ m{b}^{\mathrm{2}}& mbl\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\alpha )\end{array}\right]。

式中：K=(m{l}^{\mathrm{2}}+m{a}^{\mathrm{2}}+m_{\mathrm{H}}{l}^{\mathrm{2}})-mbl\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\alpha )。

1.2　非平整路面建模

1.2.1　正弦波路面建模

图2为正弦波路面的结构图，将图1中的平坦斜面视为构造正弦波的基准面（以虚线标识），它与水平面的夹角仍记为\phi。在此基准上建立坐标系，x轴正方向沿斜面向下，y轴正方向垂直斜面向上。此时地表轮廓函数可以表示为：

y=A\mathrm{s}\mathrm{i}\mathrm{n}(\omega x)（3）

式中：A为正弦波的幅值，控制了路面不平度的大小，用于量化路面激励的强度，mm；ω为正弦波的频率，控制了路面凹凸变化的周期性，rad/mm。通过调整这两个参数，可以模拟不同类型的路面特征。

当摆动腿与地面接触时，根据模型碰撞时刻的位姿关系，得：

h_{\mathrm{s}}+A\mathrm{s}\mathrm{i}\mathrm{n}(\omega x_{\mathrm{s}})=h_{\mathrm{n}\mathrm{s}}+A\mathrm{s}\mathrm{i}\mathrm{n}(\omega x_{\mathrm{n}\mathrm{s}})（4）

式中：x_{\mathrm{s}}、x_{\mathrm{n}\mathrm{s}}分别为支撑腿和摆动腿的横坐标；h_{\mathrm{s}}、h_{\mathrm{n}\mathrm{s}}分别为y轴方向上支撑腿和摆动腿足端到髋关节的距离，其满足如下关系式：

\left\{\begin{array}{l}{h}_{\mathrm{s}}=l\cdot \mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{s}}+\phi )\\ h_{\mathrm{n}\mathrm{s}}=l\cdot \mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{n}\mathrm{s}}+\phi )\end{array}\right.（5）

显然，通过明确双腿足端的横坐标值，即可确定机器人在碰撞时刻的位姿。鉴于此，设置机器人初始行走时的支撑腿位于坐标原点处，即x_{\mathrm{s}}=\mathrm{0}，根据几何关系，此时摆动腿的足端横坐标x_{\mathrm{n}\mathrm{s}}可表示为：

x_{\mathrm{n}\mathrm{s}}=x_{\mathrm{s}}+l\left(\left|\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{s}}+\phi )\right|+\left|\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{n}\mathrm{s}}+\phi )\right|\right)（6）

当完成一次足地碰撞后，将摆动腿足端的横坐标赋给支撑腿：{x_{\mathrm{s}}}^{(i+\mathrm{1})}={x_{\mathrm{n}\mathrm{s}}}^{(i)}，其中i为机器人行走的步数。将式（5）（6）代入式（4）中，即可得出机器人在正弦波路面上行走的触地检测条件：

\begin{array}{c}l[\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{s}}+\phi )-\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{n}\mathrm{s}}+\phi )]=A\mathrm{s}\mathrm{i}\mathrm{n}(\omega x_{\mathrm{s}})-\\ A\mathrm{s}\mathrm{i}\mathrm{n}[\omega (x_{\mathrm{s}}+l(\left|\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{s}}+\phi )\right|+\left|\mathrm{s}\mathrm{i}\mathrm{n}({\theta }_{\mathrm{n}\mathrm{s}}+\phi )\right|))]\end{array}（7）

1.2.2　随机路面建模

为了模拟更加真实的行走环境，探究路面不平整性对机器人产生的步态位移激励效应，构建了随机性的非平整路面模型，如图3所示。此路面基于图2中的正弦波路面，不同的是，此时不平度A和起伏频率ω将在一定范围内随机取值。

定义\sigma为从[\mathrm{0,1}]区间中选择的随机数，将机器人在第i步足地碰撞后的\sigma取值为{\sigma }_i，此时的路面随机幅值A_i和频率{\omega }_i为：

\left\{\begin{array}{l}{A}_i={\sigma }_i\cdot A^{\text{'}}\\ {\omega }_i={\sigma }_i\cdot {\omega }^{\text{'}}\end{array}\right.（8）

式中：A^{\text{'}}和{\omega }^{\text{'}}分别为路面函数的幅值和频率变化范围的上限。

由于被动行走机器人的吸引域非常狭窄^［20］，对外界干扰的抵抗能力相对较弱，当路面频率变化过快时，机器人可能行走失败。因此，为确保机器人能够稳定行走，需限制路面频率变化的上限{\omega }^{\text{'}}\le \mathrm{3}^［17］。此时非平整路面函数表达式为：

y_i=A_i\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{i}x)（9）

机器人在随机路面上行走的触地检测条件为：

l[\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{s}}+\phi )-\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{n}\mathrm{s}}+\phi )]\le

A_i\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_i{x}_{\mathrm{s}})-A_i\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_i{x}_{\mathrm{n}\mathrm{s}})（10）

由于道路结构变化的不规律性，在机器人的行走过程中，非平整路面将对机器人的步态产生位移激励。步态位移激励指路面的不平整性对机器人步态产生的直接影响，这包括由于路面的凹凸、坡度的变化或其他不规则特征引起的机器人行走周期演化、腿部和足部的位移变化等。本文将从系统的局部分岔和全局分岔两方面出发，对非平整路面下机器人的步态演化做出详细研究。

2 步态演化中的局部分岔

非线性系统常见的局部分岔类型包括3种：倍周期分岔（Period-doubling bifurcation， PDB）、鞍结分岔（Saddle-node bifurcation， SNB）和霍普分岔。当前的研究表明，多次的PDB是导致被动行走机器人由周期步态演化至混沌进而失稳的主要原因。为绘制机器人步态的分岔图，首先采用Newton-Raphson法求得罗盘模型在斜面角度\phi =\mathrm{0.052}\mathrm{4}（单位为rad，在下文中省略）时的周期步态吸引子：({\theta }_{\mathrm{s}},{\theta }_{\mathrm{n}\mathrm{s}},{\dot{\theta }}_{\mathrm{s}},{\dot{\theta }}_{\mathrm{n}\mathrm{s}})=（0.218 6，-0.323 4，-1.091 8，-0.377 2）。使用该值作为步态的行走初始值，利用四阶Runge-Kutta法计算步态随斜面角度\phi变化的响应特性。此外，考虑到机器人在足地碰撞时刻的雅克比矩阵不存在，本研究采用正交扰动法^［21］对步态演化中的最大和第二大李雅普诺夫指数（Lyapunov exponents， LEs）进行计算，将它们依次记为{\lambda }_{\mathrm{1}}和{\lambda }_{\mathrm{2}}。

2.1　正弦波路面上的步态演化

本节将探究正弦波路面不平度A的大小对机器人步态演化的影响，此时固定频率ω=3，路面结构由式（3）决定。图4（a）中的分岔曲线展示了A=\mathrm{0.0}时罗盘模型的步态通向混沌的经典倍周期级联路径^［8］，此时机器人最高可行走斜面角度\phi =\mathrm{0.090}\mathrm{82}。分岔曲线对应的LEs如图4（b）所示，观察到，{\lambda }_{\mathrm{1}}在\phi =\mathrm{0.076}\mathrm{64}时到达0值，此位置发生第一次PDB；在时，{\lambda }_{\mathrm{1}}\le \mathrm{0}，这表明机器人一直处于周期行走状态；而当\phi \ge \mathrm{0.087}\mathrm{96}后，{\lambda }_{\mathrm{1}}>0，此时机器人开始保持混沌步态行走。这些指数的变化情况与图4（a）中的步态演化具有良好的一致性。

图5（a）为A=\mathrm{0.1}时机器人步态的分岔图，与图4（a）对比可知，相较于平坦路面上的行走，当略微出现路面起伏时，机器人的行走步态很难再保持原有的周期性；当放大图5（a）中的矩形框时，可以观察到明显的锁相行为。将此分岔图对应的LEs图绘制在图5（b）中，观察到，当时，{\lambda }_{\mathrm{1}}始终在0值附近振荡，而{\lambda }_{\mathrm{2}}为负值，这表明当机器人受到来自行走路面的固定激励后，其步态表现为拟周期。由此图可知，第一次环面倍化分岔（Torus-doubling bifurcation， TDB）发生于\phi =\mathrm{0.076}\mathrm{55}，此时{\lambda }_{\mathrm{2}}到达0值。为更加直观地了解机器人的行走状态，在分岔图上依次选择\phi=0.052 4、0.077 52、0.086 45及0.088 15这4组参数，绘制这些参数下机器人的步态在Poincaré截面的投影面({\theta }_{\mathrm{n}\mathrm{s}},{\dot{\theta }}_{\mathrm{n}\mathrm{s}})上的相图，结果如图6（a）~（d）所示。从这4个子图中可以看到，随着\phi的增大，步态在相空间中的轨迹开始为一个二维环面，吸引子在Poincaré截面上呈现为一个T¹不变圈^［22］（见图6（a））；随着TDB的发生，吸引子演化为2T¹不变圈（见图6（b））；而当经历锁相后，2T¹不变圈又演化为稳定的长周期运动（见图6（c））；当在图6（d）所处的参数下，2T¹不变圈破裂导致步态通向混沌，此时的吸引子表现出显著的分形结构。由上述分析可知，在正弦波路面激励下，系统原本的倍周期级联现象消失，步态的失稳路径变为环面倍化分岔伴随锁相通向混沌。

当路面不平度A增大至0.3时，机器人步态的分岔图如图7（a）所示。相比于图5（a）中A=\mathrm{0.1}的情况，此时机器人的步态混乱程度明显增强。图7（b）为图7（a）对应的LEs图，由此图可知，机器人的失稳路径与A=\mathrm{0.1}时基本一致，路面不平度的增大并未导致步态演化方式发生改变。同时，TDB的发生位置也与图5（b）相同，均为\phi=0.076 55。

2.2　随机路面上的步态演化

为揭示随机路面上的步态位移激励效应，对机器人在路面不平度为A^{\text{'}}=\mathrm{0.1}和A^{\text{'}}=\mathrm{0.3}时的行走进行数值模拟，此时固定频率{\omega }^{\text{'}}=3，路面结构由式（9）决定。

图8（a）为A^{\text{'}}=0.1时机器人步态的分岔图。观察到，相比于正弦波路面上A=\mathrm{0.1}的步态演化（图5（a）），随机路面使机器人的步态偏离了原本的拟周期环面，光滑的分岔曲线变得粗糙，机器人规律性的行走状态被打破，行走步态的混乱性增加。此外，在TDB附近步态吸引子的散布区间明显增大，路面不平整性在该位置造成的步态位移激励效果增强。图8（b）为分岔曲线对应的LEs，当\phi \le \mathrm{0.086}\mathrm{82}时，{\lambda }_{\mathrm{1}}为0而{\lambda }_{\mathrm{2}}为负值，表明此时机器人的步态仍为拟周期；而当\phi \ge \mathrm{0.086}\mathrm{82}时，{\lambda }_{\mathrm{1}}为正而{\lambda }_{\mathrm{2}}为0，机器人处于混沌步态。注意到，由于路面的随机性，此时指数曲线的上下波动情况非常明显。

图9（a）为A^{\text{'}}=\mathrm{0.3}时机器人步态的分岔图，此时随机路面的激励导致步态在相空间中的轨迹严重偏离拟周期环面，分岔曲线变得极其粗糙。图9（b）为图9（a）对应的LEs图，此时指数曲线的上下波动更加剧烈。图10（a）~（c）依次为图9（a）中对应\phi=0.052 4，0.08以及0.086时吸引子的Poincaré截面图。观察到，随着参数\phi的增大，吸引子在截面上由一团离散点（见图10（a））倍化至两团（见图10（b）），实际上它们均为对拟周期步态的偏离，行走步态在此期间发生了TDB；在图10（c）中，原本应该具有清晰分形结构的混沌吸引子变得无法分辨，表明系统受到较强的外界环境干扰，机器人更加容易跌倒。

综上可知，基于正弦波路面构建的随机路面环境使机器人步态失去周期以及拟周期特性，演化为杂乱无章的双足行走状态。

3 步态演化中的全局分岔

全局分岔包括边界激变（Boundary crisis， BC）、内部激变（Interior crisis， IC）以及吸引子合并激变^［23］等动力学行为。其中，BC被认为是导致混沌吸引子突然消失或产生的直接原因，它通常伴随着系统拓扑结构的重大变化。本文将对被动行走机器人步态的双边界激变现象^［24］进行研究，以探索机器人在非平整路面上的跌倒以及运动失稳机制。

3.1　正弦波路面上的双边界激变现象

对正弦波路面不平度A=0.0及0.1这两种路面环境进行研究，\omega仍固定为3。图11为A=0.0时的步态分岔图，与2.1节中的图4（a）不同的是，除占据斜面范围最广的A₁分岔曲线外，在图11中还搜索到一支与A₁共存的、由鞍结分岔（SNB）产生的周期3步态演化而来的分岔曲线A₂。这两条分岔曲线上的吸引子都经由倍周期级联路径通向混沌，并最终消失。此外，SNB同时产生了一条不稳定周期3轨道（Period-three unstable periodic orbit， p-3 UPO），在图11中以虚线标识。随着斜面角度\phi的增加，p-3 UPO分别与A₂和A₁上的混沌吸引子发生碰撞，依次引发了两次BC事件，该事件相继造成两个步态吸引子的湮灭。

在A=\mathrm{0.1}时，步态分岔曲线如图12所示。此图中，分岔曲线A₂的吸引子出现在\phi \in [\mathrm{0.067}\mathrm{88,0.068}\mathrm{81}]的参数区间内，与A₁分岔线上的拟周期吸引子共存，同样经历了由环面倍化分岔到混沌的步态失稳路径。与图11类似，图12中的A₂同样是由SNB产生的吸引子演化而来的，其吸引子类型为3个吸引的不变圈，并且SNB也同时产生了不稳定的拟周期轨道^［25］（Unstable quasi-periodic orbit， UQO）。随着分岔参数的变化，UQO随之延伸并依次与步态A₂、A₁相撞，同样引发了两次BC。显然，与平坦路面（A=0.0）相比，一方面，路面不平度的增加并未使机器人的共存步态消失，反而由于路面激励的作用，共存的周期3步态变为拟周期步态；另一方面，机器人在受到路面激励时仍通过双边界激变引发机器人的跌倒，其跌倒机制与A=0.0时的情况相同。这种相似的全局分岔机制对不同行走环境下的机器人步态控制是有利的。

3.2　随机路面上的双边界激变现象

图13为随机路面不平度A^{\text{'}}=\mathrm{0.1}下的机器人步态分岔图。当斜面角度\phi增大至约为0.067 99时，被动行走系统中出现了SNB现象，进而产生了与图13中A₁分岔曲线共存的拟周期步态分岔A₂和对应的UQO。由于随机路面的激励，A₂与UQO的边界都不太光滑。随着参数\phi的增大，UQO先与步态A₂的三带混沌碰撞，发生了第一次BC，该事件导致了步态A₂的消失。随后，UQO仍未消失而是继续延伸，直到在\phi \approx \mathrm{0.090}\mathrm{18}位置处再次与分岔曲线A₁上的混沌吸引子相撞，发生第二次BC，进而造成步态A₁最终消失，机器人跌倒。显然，UQO在非平整路面下的机器人跌倒事件中发挥着重要作用。

综合图12与图13可知，在正弦波路面不平度A=\mathrm{0.1}与随机路面不平度A^{\text{'}}=\mathrm{0.1}下，机器人的步态周期、分岔类型以及跌倒机理都基本一致。因此，仅需通过比较不同A值下的步态演化，即可对随机路面下的行走状态做出相应判断。图14中的A₁、A₂和A₃分别为正弦波路面不平度 A等于0.0、0.1和0.3时机器人步态分岔曲线的后半部分，图中标注了它们对应的不稳定轨道以及发生第二次BC的参数值\phi。由图14可知，随着A的增大，UQO的响应范围逐渐加大，但其演化趋势总体上保持不变。同时，A的增大也导致混沌吸引子发生BC的位置提前，如A=0.0时，分岔曲线A₁最高可延伸至\phi=0.090 82处，而当A=0.1与0.3时，A₂和A₃最高可到达的斜面角度仅为0.090 14与0.087 45。由此可见，路面不平度A的增大将导致机器人在更小的斜面角度下摔倒。

3.3　ADAMS仿真与行走验证

本节展示了一款简易的点足罗盘模型样机及其行走路面，并将它们导入ADAMS 软件中开展动力学仿真。如图15所示，机器人主体由成对的内外腿组成，为了防止侧向扭矩，采用连接件保持成对腿的同步运动。使用半圆球模拟机器人的点足特征，此结构被放大到图15的右上侧。为了防止运动中的摆动腿擦地现象，行走斜面被布置为棋格状。设置各脚块表面平坦无起伏，用于模拟不平度A=0时的路面环境。当斜面倾角\phi=0.06 rad，此时机器人的步态吸引子({\theta }_{\mathrm{s}},{\theta }_{\mathrm{n}\mathrm{s}},{\dot{\theta }}_{\mathrm{s}},{\dot{\theta }}_{\mathrm{n}\mathrm{s}})=（0.223 8，-0.343 6，-1.115 8，-0.309 3），式（11）为机器人的步长公式，将{\theta }_{\mathrm{s}}和{\theta }_{\mathrm{n}\mathrm{s}}的值代入此式，可算得两个相邻脚块间的距离约为0.559 8 m。

d=\mathrm{2}l\mathrm{s}\mathrm{i}\mathrm{n}[({\theta }_{\mathrm{s}}-{\theta }_{\mathrm{n}\mathrm{s}})/\mathrm{2}]（11）

仿真开始前，为机器人各关节添加转动副，并参照表1中的参数为机器人设置质量、质心位置等结构信息；为模拟完全非弹性碰撞，设置足地碰撞时的恢复系数为0；令机器人的双腿垂直于斜面，仅对摆动腿添加1.610 775 rad/s的初始角速度，经验证，此时机器人能够在斜面上保持稳定行走。

图16（a）为机器人在ADAMS仿真环境下的行走过程，此时足部的落脚点大致位于脚块的中间位置。图16（b）与16（c）分别为双腿角度以及角速度的时间响应图，虚线为数值仿真曲线，实线为ADAMS仿真曲线。注意到，经由ADAMS仿真获得的行走步态，其角度及角速度都低于数值仿真结果，这是由于在ADAMS环境中为机器人的腿部设置了相应的转动惯量引起的。由步态的时间响应可知，使用两种不同的仿真方式得到的机器人稳定行走步态均为周期1，且角度及角速度的变化趋势以及每步间的相位差是接近一致的，这证明了机器人行走能力以及理论分析的正确性。

4 结　论

（1）路面不平整性产生的步态位移激励导致机器人步态呈现非常规的失稳路径，即由T¹不变圈（拟周期运动）经环面倍化分岔到2T¹不变圈，再经历环面振荡、锁相等动力学行为演化为混沌步态。

（2）非平整路面激励导致机器人的步态在相空间中的轨迹偏离了极限环，步态的混乱程度增强，系统的稳定性下降。随着路面不平度的增大，机器人的步态演化机制并未改变，但是机器人能够行走的最大坡度显著降低。

（3）机器人在非平整路面行走时，系统通过鞍结分岔产生共存的拟周期步态和不稳定的拟周期轨道，后者将与混沌吸引子相撞引发步态吸引子的湮灭，这被称为边界激变。值得强调的是，不稳定轨道引发了两次激变，这种行为将导致步态吸引子消失，造成机器人摔倒。

参考文献

[1] Moon J S, Lee M, Bae J, et al. Analysis of period-1 passive limit cycles for flexible walking of a biped with knees and point feet[J]. Robotica, 2016, 34(11): 2486-2498.

[2] Znegui W, Gritli H, Belghith S. A new Poincaré map for investigating the complex walking behavior of the compass-gait biped robot[J]. Applied Mathematical Modelling, 2021, 94: 534-557.

[3] Wu Y, Qiao S, Yao D. A hybrid chaotic controller integrating hip stiffness modulation and reinforcement learning-based torque control to stabilize passive dynamic walking[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2023, 237(3): 673-691.

[4] Mcgeer T. Passive dynamic walking[J]. The International Journal of Robotics Research, 1990, 9(2): 62-82.

[5] Iqbal S, Zang X, Zhu Y, et al. Bifurcations and chaos in passive dynamic walking: a review[J]. Robotics and Autonomous Systems, 2014, 62(6): 889-909.

[6] 钱志辉, 周亮, 任雷, 等. 具有仿生距下关节和跖趾关节的完全被动步行机[J]. 吉林大学学报: 工学版, 2018, 48(1): 205-211.

[7] Qian Zhi-hui, Zhou Liang, Ren Lei, et al. Completely passive walking machine with bionic subtalar jointand matatarsal phalangeal joint[J]. Journal of Jilin University(Engineering and Technology Edition), 2018, 48(1): 205-211.

[8] 高家昌, 高建设, 陶征. 被动步行平足机器人动力学参数研究[J]. 机械传动, 2023, 46(12): 22-30.

[9] Gao Jia-chang, Gao Jian-she, Tao Zheng. Research on dynamic parameters of a passive walking robot with flat feet[J]. Journal of Mechanical Transmission, 2023, 46(12): 22-30.

[10] Goswami A, Thuilot B, Espiau B. A study of the passive gait of a compass-like biped robot: symmetry and chaos[J]. International Journal of Robotics Research, 1998, 17(12): 1282-1301.

[11] Garcia M, Chatterjee A, Ruina A, et al. The simplest walking model: stability, complexity, and scaling[J]. Journal of Biomechanical Engineering, 1998, 120(2): 281-288.

[12] Wisse M, Schwab A L, van der Helm F C T. Passive dynamic walking model with upper body[J]. Robotica, 2004, 22(6): 681-688.

[13] Narukawa T, Takahashi M, Yoshida K. Level-ground walk based on passive dynamic walking for a biped robot with torso[C]∥Proceedings 2007 IEEE International Conference on Robotics and Automation, Rome, Italy, 2007: 3224-3229.

[14] 柳宁. 双足被动行走器动力学仿真与实验研究[D]. 北京:清华大学航天航空学院, 2009.

[15] Liu Ning. Simulation and experimental study on the dynamics of passive bipedal walker[D]. Beijing: School of Aerospace Engineering, Tsinghua University, 2009.

[16] Kim J, Choi C H, Spong M W. Passive dynamic walking with symmetric fixed flat feet[C]∥2007 IEEE International Conference on Control and Automation, Guangzhou, China, 2007: 24-30.

[17] Wu B, Zhao M. Bifurcation and chaos of a biped robot driven by coupled elastic actuation[C]∥Proceeding of the 11th World Congress on Intelligent Control and Automation, Shenyang, China, 2014: 1905-1910.

[18] Safartoobi M, Dardel M, Daniali H M. Gait cycles of passive walking biped robot model with flexible legs[J]. Mechanism and Machine Theory, 2021, 159: No.104292.

[19] Iida F, Tedrake R. Minimalistic control of biped walking in rough terrain[J]. Autonomous Robots, 2010, 28: 355-368.

[20] Afshar P N, Ren L. Dynamic stability of passive bipedal walking on rough terrain: a preliminary simulation study[J]. Journal of Bionic Engineering, 2012, 9(4): 423-433.

[21] Su J L S, Dingwell J B. Dynamic stability of passive dynamic walking on an irregular surface[J]. Journal of Biomechanical Engineering, 2007, 129(6): 802-810.

[22] Bogsjö K, Podgórski K, Rychlik I. Models for road surface roughness[J]. Vehicle System Dynamics, 2012, 50(5): 725-747.

[23] Rao X B, Gao J S, Ding S L, et al. Multistability of gaits, the basin of attraction and its external topology in the simplest passive walking model on stairs[J]. Chaos, Solitons & Fractals, 2023, 172:No. 113592.

[24] Balcerzak M, Dąbrowski A, Stefański A, et al. Spectrum of Lyapunov exponents in non-smooth systems evaluated using orthogonal perturbation vectors[C]∥MATEC Web of Conferences, Lodz, Poland, 2018, 148: No.10003.

[25] 丁旺才. 多自由度碰撞振动系统的环面分岔与混沌研究[D]. 成都: 西南交通大学力学与工程学院, 2004.

[26] Ding Wang-cai. Research on torus bifurcations and chaos of multi-degree-of-freedom vibro-impact systems[D]. Chengdu: School of Mechanics and Engineering, Southwest Jiaotong University, 2004.

[27] Singh J P, Roy B K. Crisis and inverse crisis route to chaos in a new 3-D chaotic system with saddle, saddle foci and stable node foci nature of equilibria [J]. Optik, 2016, 127(24): 11982-12002.

[28] Gritli H, Belghith S, Khraief N. Cyclic-fold bifurcation and boundary crisis in dynamic walking of biped robots[J]. International Journal of Bifurcation and Chaos, 2012, 22(10): No.1250257.

[29] Zhang Y, Luo G. Detecting unstable periodic orbits and unstable quasiperiodic orbits in vibro-impact systems[J]. International Journal of Non-Linear Mechanics, 2017, 96: 12-21.

基金资助

国家自然科学基金项目(59106225)