1. 引言

人与结构之间的相互作用研究具有广泛的应用背景，如人与多层结构，人与网架结构，以及人与人行桥等。在这些研究中一个比较基本的问题是如何确定人产生的激励力。以人与人行桥相互作用的研究为例，在人–桥相互作用的研究中通常采用半经验的方法来确定行人产生的激励力的大小 [1] [2] [3] ：根据对千禧桥的封闭测试，Dallard等 [4] 假设行人产生的侧向力和桥面运动的速度成正比。而Nakamura等 [5] 则认为仅当桥面运动的速度比较小时行人产生的侧向力与桥面的侧向速度成正比，当桥面的速度较大时行人产生的侧向激励力会趋向一个有限值而不会无限加大。Roberts [6] 认为人和桥之间的相互作用力总是简谐变化且频率总在1 Hz附近。尽管通过半经验方法得到的这些行人激励力的假设模型在实际应用中发挥了很大的作用 [7] [8] [9] ，但对机理性研究往往还不够深入。例如，为何人行桥的过大侧向振动总是集中在1 Hz附近，千禧桥中著名的跳跃现象的机理等，并不能根据这些经验侧向力模型得到很好的解释。因此有必要从理论上对行人产生的侧向力进行建模和分析。

当行人行走在静止的刚性平面上时由于身体的重心的振动会使得人在水平和竖直方向上都产生一个随着时间变化的激励力。正常行走时激励力主要集中在竖直方向，其大小在体重附近波动，而侧向力的大小约为体重的4%，频率为竖直方向步频的一半，约在0.7~1.2 Hz之间。在生物力学领域倒摆模型常用于研究人体的平衡性 [10] 。一些研究倾向于寻找能直接计算行人行走动力行为的公式 [11] ，而另外一些研究考虑用系统方法研究人体的姿态控制 [12] 。对人体进行真实的生物力学建模非常复杂，很难直接应用到人与结构的相互作用研究中。更为实际的做法是抓住人与结构相互作用过程中的主要特性，对行人的动力学建模进行简化。根据实验观察，Erlicher等 [13] 认为行人行走时身体具有自我保持性。在此观察基础上较为合理并可行的做法是将人体简化为倒摆。Macdonald [10] 采用倒摆模型模拟人的侧向力讨论了人–桥之间的相互作用，得到了和实测较为吻合的结果。在他们的研究中并未考虑倒摆支点的持续运动，因此无法讨论人的行走和行走过程中人体的稳定性。但这两点在人与结构相互作用的研究中往往非常重要，因为结构的振动会影响人的稳定性，而人体稳定性的改变又会影响到行走产生的激励力，从而反过来影响结构的振动。为更合理地用倒摆模型来模拟人的行走，本文考虑倒摆的支点在水平和竖直方向同时进行简谐运动。根据非线性动力学稳定性理论分析支点运动下倒摆的稳定条件，并在参数平面上利用摄动方法计算分界线(transition curves)的位置，由此讨论参数变化对行人行走稳定性的影响。我们的分析表明行人行走时水平和竖向方向步频之比，振幅之比都对人体的稳定性具有很大的影响，这和对行人正常行走的观测结果是吻合的。

本文剩下部分结构安排如下：在第二节中建立倒摆模型；在第三和四节中分别讨论倒摆系统的通解和稳定性分界线；在第五节中讨论参数变化对倒摆稳定性的影响；结论在第六节中给出。

2. 倒摆模型

本文用倒摆模拟行人在静止平面上的行走，如图1所示。摆球代表行人的质心，摆杆代表行人的腿。当人行走时质心上下左右波动，因此我们考虑摆杆的支点上下左右进行简谐运动来模拟人的行走。摆球的质量，摆杆的长度和支点竖向简谐运动的频率分别为m，L和ω。支点摆动的水平和竖直方向的幅值假设分别为A，B。n为水平和竖直振动频率的比值，ϕ为水平和竖直方向频率的相位差，g为重力加速度。图1中所示倒摆的控制方程可由Lagrange方法得到。假设 $\theta$是摆杆运动时任意时刻与竖直方向所成的角度， $\left(x\left(t\right),y\left(t\right)\right)$是为任意时刻倒摆支点的位置坐标，则总有和 $y\left(t\right)=A\sin \left(\omega t\right)$。摆球的位置坐标可以写为 $\left(x-L\sin \theta ,y+L\cos \theta \right)$。对位置坐标求导，可以得到摆球的速度

$\upsilon =\sqrt{{\stackrel{\dot{}}{x}}^2-2L\cos \theta \stackrel{\dot{}}{x}\stackrel{\dot{}}{\theta }+y^2-2L\stackrel{\dot{}}{y}\stackrel{\dot{}}{\theta }\sin \theta +L^2{\stackrel{\dot{}}{\theta }}^2}$。(1)

倒摆系统的Lagrangian可以写为

$M=T-V=\frac{1}{2}m{\upsilon }^2-mg\left(y+L\cos \theta \right)$，(2)

将(1)和(2)带入下面的Lagrange方程

$\frac{\text{d}}{\text{d}t}\frac{\partial M}{\partial \stackrel{\dot{}}{\theta }}-\frac{\partial M}{\partial \theta }=0$，(3)

有

$L\ddot{\theta }+\left(-g+A{\omega }^2\sin \left(\omega t\right)\right)\sin \theta +B{n}^2{\omega }^2\cos \left(n\omega t+\varphi \right)\cos \theta =0$。(4)

将时间重新标度 $\omega t\to t$，(4)可以写为

$\ddot{\theta }+\left(\delta +\epsilon \sin t\right)\sin \theta +\nu {\epsilon }^2\sin \left(nt+\varphi \right)\cos \theta =0$，(5)

其中，

$\delta =-\frac{g}{L{\omega }^2}$， $\epsilon =\frac{A}{L}$， $\nu =\frac{n^{2}BL}{A^2}$。

对于正常行走的行人 $0<B<A<L$和 $0<n<1$成立，所以 $\epsilon$是一个小量。考虑到 $\left|\theta \right|<1$，(5)可以近似写为

$\ddot{\theta }+\left(\delta +\epsilon \sin t\right)\left(\theta -\frac{1}{6}{\theta }^3\right)+\nu {\epsilon }^2\sin \left(nt+\varphi \right)\left(1-\frac{1}{2}{\theta }^2\right)=0$。(6)

3. 倒摆控制方程的通解

如果 $\epsilon \to 0\left(A\to 0\right)$，(6)退化为

$\ddot{\theta }+\delta \left(\theta -\frac{{\theta }^3}{6}\right)=0$，(7)

因为 $\delta <0$，(7)没有周期解，如果 $\left|\theta \right|<1$， $\theta =0$是唯一不稳定的平衡点。 $\epsilon \to 0$意味着行人停止行走并保持站立。本文重点关注 $\delta =0$， $\epsilon =0$附近(6)的解的稳定性。因此可以通过在(7)的解附近摄动来获得(6)的解。在 $\epsilon -\delta$参数平面上，(7)的解稳定性区域由两条分界线包围而成，一条从 $\delta =0$发出，而另一条从 $\delta >0$发出。下面考虑(7)在 $\delta >0$情况下如下形式的通解

$\theta \left(t\right)=\epsilon {\theta }_1\left(t\right)+{\epsilon }^2{\theta }_2\left(t\right)+\cdots$。(8)

将(8)带入(6)并且分别令 $\epsilon$的各次幂系数为0，有

$\begin{array}{l}{\ddot{\theta }}_1+\delta {\theta }_1=0,\\ {\ddot{\theta }}_2+\delta {\theta }_2=-\sin t{\theta }_1-\nu \sin \left(nt+\varphi \right).\end{array}$(9)

(9)中第一式的通解可以写为

${\theta }_1=a\cos \left(\omega t\right)+b\sin \left(\omega t\right)$，(10)

其中 ${\omega }^2=\delta >0$，a和b是依赖初始条件的常数。将(10)带入(9)的第二式，可以得到其特解

$\begin{array}{c}{\theta }_2=-\frac{1}{2\left(2\omega -1\right)}\left[b\cos \left(\left(\omega -1\right)t\right)-a\sin \left(\left(\omega -1\right)t\right)\right]-\frac{\nu \sin \left(nt+\varphi \right)}{{\omega }^2-n^2}\\ -\frac{1}{2\left(2\omega +1\right)}\left[b\cos \left(\left(\omega +1\right)t\right)-a\sin \left(\left(\omega +1\right)t\right)\right].\end{array}$(11)

$\delta >0$时(6)的二阶近似解可以写为

$\begin{array}{c}\theta =\epsilon \left(a\cos t+b\sin t\right)+{\epsilon }^2\{-\frac{1}{2\left(2\omega -1\right)}\left[b\cos \left(\left(\omega -1\right)t\right)-a\sin \left(\left(\omega -1\right)t\right)\right]\\ -\frac{1}{22\left(2\omega +1\right)}\left[b\cos \left(\left(\omega +1\right)t\right)-a\sin \left(\left(\omega +1\right)t\right)\right]-\frac{\nu \sin \left(nt+\varphi \right)}{{\omega }^2-n^2}\}\end{array}$(12)

显然，系统(6)的共振条件为 $\omega =\pm 1/2$和 $\omega =\pm n$(或者等价条件 $\delta =1/4$和 $\delta =n^2$)。当 $n=1/2$时，系统(6)中出现强共振。四种不同参数情况下(12)的解析解和直接数值结果的比较结果如图2所示，表明了(12)的有效性。需要注意的是，当 $\epsilon$增大时(12)的精度会下降，如果要提高精度，需要计算 $\epsilon$的高阶项。

4. 分界线和周期解

在这一节中讨论 $\epsilon -\delta$参数平面上(6)解的稳定性区域。根据上一节的讨论，稳定性区域由从δ-轴上 $\delta =0$， $\delta =n^2\left(n\ne 1/2\right)$和 $\delta =1/4$分别出发的分界线所包围而成。下面分别进行讨论。

首先假设参数δ能写成如下展开形式

$\delta ={\delta }_0+{\delta }_1\epsilon +{\delta }_2{\epsilon }^2+\cdots$，(13)

其中 ${\delta }_i$是常数，可通过避免久期项来计算。将(8)和(13)带入(6)并且令 $\epsilon$的各次幂系统为0，有

${\ddot{\theta }}_1+{\delta }_0{\theta }_1=0$，(14)

${\ddot{\theta }}_2+{\delta }_0{\theta }_2=-\left({\delta }_1+\sin t\right){\theta }_1-\nu \sin \left(nt+\varphi \right)$，(15)

${\ddot{\theta }}_3+{\delta }_0{\theta }_3=-\left({\delta }_1+\sin t\right){\theta }_2-\left({\delta }_2{\theta }_1-\frac{1}{6}{\delta }_0{\theta }_1^3\right)$。(16)

4.1. δ₀= 0

将 ${\delta }_0=0$带入(14)，有 ${\theta }_1=a$，其中 $a$为常数依赖初始条件。此时(15)变为

${\ddot{\theta }}_2+a\left({\delta }_1+\sin t\right)+\nu \sin \left(nt+\varphi \right)=0$。(17)

如果令 ${\delta }_1$为0， ${\theta }_2$是周期的。(17)的特解可以写为

${\theta }_2=a\sin t+\frac{\nu }{n^2}\sin \left(nt+\varphi \right)$。(18)

将(18)带入(16)并消去久期项，有

${\theta }_3=-\frac{a}{8}\cos \left(2t\right)+\frac{\nu \cos \left[\left(n-1\right)t+\varphi \right]}{2n^2{\left(n-1\right)}^2}-\frac{\nu \cos \left[\left(n+1\right)t+\varphi \right]}{2n^2{\left(n+1\right)}^2}.$(19)

从δ-轴上 ${\delta }_0=0$点出发的分界线可以近似表示为

$\delta =-\frac{1}{2}{\epsilon }^2$，(20)

分界线上的周期解可以解析地表示为

$\theta =a\epsilon +\left[a\sin t+\frac{\nu }{n^2}\sin \left(nt+\varphi \right)\right]{\epsilon }^2+\left[-\frac{1}{8}a\cos \left(2t\right)+\frac{\nu \cos \left(\left(n-1\right)t+\varphi \right)}{2n^2{\left(n-1\right)}^2}-\frac{\nu \cos \left(\left(n+1\right)t+\varphi \right)}{2n^2{\left(n+1\right)}^2}\right]{\epsilon }^3.$(21)

4.2. δ₀= n²(n ≠ 1/2)

将 ${\delta }_0=n^2$带入(14)，得到(14)的通解

${\theta }_1=a\cos \left(nt\right)+b\sin \left(nt\right)$。(22)

将(22)带入(15)，(15)变为

${\ddot{\theta }}_2+n^2{\theta }_2=-\left({\delta }_1+\sin t\right)\left[a\cos \left(nt\right)+b\sin \left(nt\right)\right]-\nu \sin \left(nt+\varphi \right)$。(23)

要 ${\theta }_2$为周期解，必须要消去久期项，因此有

$\begin{array}{l}\nu \cos \left(\varphi \right)+b{\delta }_1=0,\\ \nu \sin \left(\varphi \right)+a{\delta }_1=0.\end{array}$(24)

a，b是任意依赖初始条件的常数，(24)的非平凡解可以写为

$a=0,\varphi =0,{\delta }_1=-\frac{\nu }{b},$

或者

$b=0,\varphi =\frac{\text{π}}{2},{\delta }_1=-\frac{\nu }{a}.$

4.2.1. $a=0,\varphi =0,{\delta }_1=-\frac{\nu }{b},$

在这种情况下(23)的特解可以表示为

${\theta }_2=-\frac{b\cos \left(\left(n-1\right)t\right)}{2\left(2n-1\right)}-\frac{b\cos \left(\left(n+1\right)t\right)}{2\left(2n+1\right)}.$(25)

通过和前面相同的步骤，可以得到

$\begin{array}{l}{\delta }_2=\frac{4b^2{n}^4-b^2{n}^2+4}{8\left(4n^2-1\right)},\\ {\theta }_3=-\frac{\nu \cos \left(\left(n-1\right)t\right)}{2{\left(2n-1\right)}^2}+\frac{\nu \cos \left(\left(n+1\right)t\right)}{{\left(2n+1\right)}^2}+\frac{b^3}{192}\sin \left(3nt\right)-\frac{b\sin \left(\left(n-2\right)t\right)}{16\left(2n-1\right)\left(n-1\right)}-\frac{b\sin \left(\left(n+2\right)t\right)}{16\left(2n+1\right)\left(n+1\right)}.\end{array}$(26)

此时分界线可以表示为

$\delta =n^2-\frac{\nu }{b}\epsilon +\frac{4b^2{n}^4-b^2{n}^2+4}{8\left(4n^2-1\right)}{\epsilon }^2.$(27)

4.2.2. $b=0,\varphi =\pi /2,{\delta }_1=-\nu /a$

在这种情况下分界线和分界线上的周期解可以分别表示为

$\delta =n^2-\frac{\nu }{a}\epsilon +\frac{4a^2{n}^4-a^2{n}^2+4}{8\left(4n^2-1\right)}{\epsilon }^2,$(28)

$\begin{array}{c}\theta =a\cos \left(nt\right)\epsilon +[\frac{a\sin \left(\left(n-1\right)t\right)}{2\left(2n-1\right)}+\frac{a\sin \left(\left(n+1\right)t\right)}{2\left(2n+1\right)}]{\epsilon }^2+[\frac{\nu \sin \left(\left(n-1\right)t\right)}{2{\left(2n-1\right)}^2}\\ -\frac{\nu \sin \left(\left(n+1\right)t\right)}{2{\left(2n+1\right)}^2}-\frac{a\cos \left(\left(n-2\right)t\right)}{\left(2n-1\right)\left(n-1\right)}-\frac{a\cos \left(\left(n+2\right)t\right)}{\left(2n+1\right)\left(n+1\right)}-\frac{a^3}{192}\cos \left(3nt\right)]{\epsilon }^3.\end{array}$(29)

4.3. δ₀= 1/4

和前面过程类似，可以分为两种情况，下面分别讨论。

4.3.1. a = b

在这种情况下分界线和分界线上的周期解可以分别表示为

$\delta =\frac{1}{4}-\frac{1}{2}\epsilon +\left(\frac{b^2}{16}-\frac{1}{8}\right){\epsilon }^2,$(30)

$\begin{array}{c}\theta =b\left[\cos \left(\frac{1}{2}t\right)+\sin \left(\frac{1}{2}t\right)\right]\epsilon +\frac{1}{4}\left[-b\cos \left(-\frac{3}{2}t\right)+b\sin \left(-\frac{3}{2}t\right)+\frac{16\nu \sin \left(nt+\varphi \right)}{4n^2-1}\right]{\epsilon }^2\\ +[\frac{b\left(b^2+6\right)}{96}\left(\cos \left(-\frac{3}{2}t\right)+\sin \left(-\frac{3}{2}t\right)\right)-\frac{b}{48}\left(\cos \left(\frac{5}{2}t\right)+\sin \left(\frac{5}{2}t\right)\right)\\ +\frac{8\nu \cos \left(\left(n-1\right)t+\varphi \right)}{\left(2n+1\right)\left(2n-3\right){\left(2n-1\right)}^2}-\frac{8\nu \cos \left(\left(n+1\right)t+\varphi \right)}{\left(2n-1\right)\left(2n+3\right){\left(2n+1\right)}^2}-\frac{8\nu \sin \left(nt+\varphi \right)}{{\left(4n^2-1\right)}^2}]{\epsilon }^3.\end{array}$(31)

4.3.2. a = −b

在这种情况下分界线和分界线上的周期解可以分别表示为

$\delta =\frac{1}{4}\text{+}\frac{1}{2}\epsilon +\left(\frac{b^2}{16}-\frac{1}{8}\right){\epsilon }^2,$(32)

$\begin{array}{c}\theta =b\left[-\cos \left(\frac{1}{2}t\right)+\sin \left(\frac{1}{2}t\right)\right]\epsilon -\frac{1}{4}\left[b\cos \left(\frac{3}{2}t\right)+b\sin \left(\frac{3}{2}t\right)-\frac{16\nu \sin \left(nt+\varphi \right)}{4n^2-1}\right]{\epsilon }^2\\ -[\frac{b\left(b^2+6\right)}{96}\left(\cos \left(\frac{3}{2}t\right)+\sin \left(\frac{3}{2}t\right)\right)+\frac{b}{48}\left(\cos \left(\frac{5}{2}t\right)+\sin \left(\frac{5}{2}t\right)\right)\\ +\frac{8\nu \cos \left(\left(n-1\right)t+\varphi \right)}{\left(2n+1\right)\left(2n-3\right){\left(2n-1\right)}^2}-\frac{8\nu \cos \left(\left(n+1\right)t+\varphi \right)}{\left(2n-1\right)\left(2n+3\right){\left(2n+1\right)}^2}+\frac{8\nu \sin \left(nt+\varphi \right)}{{\left(4n^2-1\right)}^2}]{\epsilon }^3.\end{array}$(33)

根据上面的解析表达式，在分界线上有无穷多的周期解。以从 $\delta =0$出发的分界线为例，在(20)中令 $\delta =-0.02$，则有 $\epsilon =0.2$。另外，令 $\nu =0.01$，n = 0.4， $\varphi =0.15$，并且取初始条件为 $\theta \left(0\right)=0$和 $\stackrel{\dot{}}{\theta }\left(0\right)\text{=0}\text{.001}$。直接采用数值方法对(6)进行积分得到倒摆摆角的时程图，相图和时程曲线的傅里叶谱，如图3所示。从图3(c)可以分辨出傅里叶谱中的频率分量为 $\omega =0.4,0.6,1,1.4\cdots$，这和解析解(21)的结果是吻合的。接下来选取 $\delta =-0.02$， $\nu =0.01$， $n=0.8$， $\varphi =0.15$， $b=0.016$， $\epsilon =0.482$，此时(31)为近似解析解。直接对(6)进行积分得到时程图，相图和时程曲线的傅里叶谱，如图4所示。从图4(c)可以看到，周期解有频率分量 $\omega =0.2,0.5,0.8,1.5$，这和(31)的结果是吻合的。

5. 稳定性区域分析和讨论

本节中用(20)，(27)，(28)，(30)和(32)来研究参数平面中的稳定性区域随参数改变的变化。令 $\nu$， $\varphi$，a和 $b$的值分别为0.01，0.15，0.1，−0.2。n是倒摆支点水平与竖向振动频率之比。令n取0.25，0.45，0.55和0.75，在 $\epsilon -\delta$参数平面上分别画相应的分界线并标注稳定性区域，如图5所示。从图中可以看出n值的变化对分界线以及稳定性区域都有很大的影响。在条件 $\epsilon >0$和 $\delta <0$下，如果 $n<1/2$，随着n值的增大参数平面上稳定性区域不断增加；但当 $n>1/2$后n值再增加不会改变稳定性区域的面积。此外，当 $n<1/2$时，稳定性区域是不连续的。图5表明 $n>1/2$对于倒摆的稳定性而言是有利的。

需要指出的是，当 $n>1/2$时依赖初始条件的常数a和b都对稳定性区域的面积有影响。在对(6)直接数值积分时如果选择a = 0.5，b = −0.01，且其它参数与图6中保持一致，稳定性区域如图6所示。

对比图5和6，当 $n<1/2$时初始条件对(6)的解的稳定性区域面积和连续性都有重要影响。而当 $n>1/2$时，初始条件几乎不再有影响。

$\nu$是倒摆支点水平和竖向振幅之比。为分析 $\nu$对稳定性区域的影响，分别令 $\nu$取0.05，0.1，0.15。其它参数和图5中参数保持一致。 $\nu$对稳定区域的影响如图7所示。从图5和图7可以看出，当 $n<1/2$时 $\nu$对稳定性区域有很大的影响；当 $n>1/2$时 $\nu$的变化几乎没有影响。当 $\nu$很小并接近0时， $n<1/2$时的稳定性区域的面积可能超过 $n>1/2$时，但 $n<1/2$时稳定性区域是不连续的(图5(b))。当 $n<1/2$时，随着 $\nu$的增加，稳定性区域的不连续性开始消失，稳定性区域的面积显著下降。对比图5和图7， $\nu$的增加总是使得不稳定性增加，此时不管n值的大小，通过减小 $\epsilon$和δ来增加稳定性总是有效的。根据图5~7的分析对倒摆系统如果要获得好的稳定性需要 $n>1/2$。

考虑到真实行人的重心高度和竖向步频，重点关注δ的值在−0.03到−0.13范围内变化， $\epsilon$的值在0到0.4范围内变化。行人行走时重心上下振动，与此同时身体不可避免地会左右摇摆，因此侧向摆动幅度不会为0。根据上面对倒摆的稳定性分析，如果侧向步频与竖向步频之比小于1/2，身体保持稳定性的区域在参数平面上是不连续的。这意味着若要保持身体平衡可能需要不间断进行跳跃性调整，这对于行走而言是不合适的。因此从稳定性角度需要保持侧向和竖向步频之比大于1/2。而从能量消耗角度来看，步频之比大于1/2并不经济。因此，人行走时保持侧向和竖向步频之比在1/2附近是最为合理的。这和对行人行走时的观测结果是吻合的。另外，当行人的侧向和竖向步频之比小于1/2，则身体保持稳定所对应的参数区域的面积和连续性都极易被侧向和竖向步幅之比以及初始条件所影响。这可以用于解释为何一个喝醉的人或者刚开始学走路的婴幼儿非常容易摔倒。因为他们不能保证侧向步频与竖向步频之比足够大，并且难以及时调整步频和身体重心的高度。根据图5~7可以推知当行人的侧向摆幅比较大的，降低身体重心和加快竖向步频总可以有效地增强身体的平衡。这和对行人行走的实际观测是一致的。弯腰、快速调整落脚速度都是人体对脚下不稳情况下的自然反应。

6. 结论

本文采用倒摆模型研究行人行走的动力行为。考虑倒摆支点在水平和竖直方向做简谐运动来模拟行人行走时身体质心在侧向和竖向的起伏运动。利用摄动方法通过计算参数平面上的分界线得到了倒摆系统的稳定性条件，并用数值计算结果验证了理论分析的正确性。本文利用摆动稳定性分析论证了人行走时保持侧向和竖向步频之比在1/2附近是最为合理的。根据本文讨论，用支点作简谐运动的倒摆来模拟人的行走可以用于分析行人的动力学行为。通过将倒摆和结构耦合可以用于人–结构相互作用的分析，这可以避免直接确定行人产生的激励力带来的分析困难。这为理论分析人–结构耦合作用提供了一条新的思路。

基金项目

本文工作得到了国家自然科学基金项目(11472160, 51708349)的资助。

参考文献