四剛體人體行走模型及其人行荷載分析

2022-09-03 09:11:00李永甫楊浩文

振動工程學報 2022年4期

李永甫，吳斌，包宇，楊浩文

（1.武漢理工大學土木工程與建筑學院，湖北武漢 430070；2.哈爾濱工業大學土木工程學院，黑龍江哈爾濱 150090）

引言

在建筑中，人體完成步行、跳躍、奔跑、屈伸、上下樓梯等動作時，施加于結構上的動力作用可引起工程結構的振動，嚴重時會影響結構的正常使用［1?3］。例如英國倫敦千禧橋，由于結構設計時未重視人致荷載對結構的影響［4］，低估了人與結構間的相互作用［5］，在開放3 天后橋梁就因為嚴重的人致振動問題而不得不關停。

以往對結構進行動力響應分析時，一般將人體看作外部動荷載或附加質量［1］。但大量的試驗數據表明，即便是靜止人體，對結構的作用也不能單純等效為靜荷載。土木工程結構中人不僅作為外部荷載激勵結構，同時能夠改變所在結構的頻譜特性［6］。Nimmen 等［7］發現人體對結構動力特性的影響與姿態有關；劉澤龍等［8］采用IDA 方法對教學樓進行地震動響應分析，指出相較于空載和靜止人群工況，結構在運動人群工況下動力響應最大。在設計人群密集的建筑結構時，需要對人與結構相互作用進行充分考慮，確定一個合理的人體動力模型愈發重要。

基于人體行走時的物理特性所建立的運動模型能反映人與結構間的相互作用。學者們對單質點、倒立擺、多剛體等模型進行過研究。Bocian等［9］使用倒立擺模型進行人行橋系統動力計算，與測力跑步機上實測足底力結果吻合較好；Qin 等［10］為支撐桿件添加時變剛度阻尼，得到可用于結構分析的倒立擺模型。前述模型通過直接定義人體動力屬性代替行走過程中各體段姿態改變對運動的影響，進而固定了模型的運動形式。其中經典的倒立擺模型，剛性桿件無法產生變形，使得質心的運動軌跡與實際值相差較大，進而造成了受力性能存在較大差異。

人體正常運動時，骨骼以及附屬肌肉顫動對運動影響有限。將人體骨骼及外部組織視為剛體，人體可劃分為多個剛體組成的運動系統。多剛體模型可以很好地再現人體運動軌跡。張夢詩等［11］通過實測步行軌跡，利用十五剛體模型對人步行進行軌跡重構并對比實測步行荷載，驗證了多剛體模型模擬步行荷載的可行性；Mazzoleni 等［12］使用此類模型完成人體跳躍荷載的重構。上述剛體模型通過代入實測軌跡，得到擬合結果較好的步行荷載曲線，研究思路與時域力荷載模型類似，但未建立模型動力方程。

基于此，本文擬建立簡化人體多剛體模型，類比人體行走驅動模式，在各關節施加扭矩，并利用反饋線性化控制方法調節模型扭矩輸入，使模型能完成擬人行走。模型忽略人體運動時變化幅度較小的軀干，采用四剛體建模，對象分別為下肢的大小腿。并建立擺動階段的動力方程，求解模型行走時程。從模型行走軌跡和力學性能出發，判斷模型行走擬人程度，并驗證模型對不同身高、體重人體的適用性。

1 模型建立

1.1 動力學方程

正常人體行走單個步行周期為一側下肢完成從足跟著地到再次落地的過程，包括該下肢的支撐相與擺動相。在忽略人體行走左右腿區別后，模型單個行走周期定義為單腿完成一次擺動與碰撞切換過程。將人體沿矢狀面簡化為平面內多剛體模型，以描述人正常行走。

本文對人體進行簡化，軀干及以上簡化成單個質點，大腿和小腿簡化為四個均勻分布的剛性桿，編號為1～4，四個剛性桿之間以及支撐點與地面之間連接為鉸接，驅動力簡化為集中在連接節點①～③處的力矩。各物理量定義如圖1所示，其中大腿與小腿等長，長度為l，質量分別為mt，ms，集中于剛體質心，軀干與髖關節的集中質量為m1。根據歐拉?拉格朗日方程，建立模型擺動階段運動方程如下：

圖1 人體簡化模型及扭矩作用形式Fig.1 Four rigid body model and torque action form for human walking

式中q=[θ，q2，q3，q4]T為角位移，以逆時針轉動為正；θ表示髖關節節點③與模型觸地點連線相對于y軸的角度，為欠驅動自由度；qa=[q2，q3，q4]T，分別為剛性桿j（j=2，3，4）相對于剛性桿j-1 的轉動角度，為驅動自由度；M(q)∈R4×4為慣性矩陣，M11為慣性矩陣中第一個元素；h1(q，)∈R1×1，h2(q，)∈R3×1包含方程中離心力、科氏力和重力，其表達式見附錄，以及文獻［13］；u=[u1，u2，u3]T為模型各關節處獨立驅動力矩；B∈R4×3為外力矩位置矩陣，其數值根據廣義力定義求得：

1.2 碰撞切換

多剛體系統發生碰撞過程時，常以模型碰撞后的角速度作為未知量，假定碰撞前后模型位形不變，根據沖量定理建立方程。但是這種方法計算量大，且迭代不易收斂。實際上人作為一種智慧生物，其在支撐腿蹬地與擺動腿足部完全觸地過程中存在著復雜的機制以使人體能夠穩定地切換到下一步態。模擬人體這種生物力學機制大大超出了本文的范圍，作者更關心切換的目標，即觸地碰撞后任意一肢的速度等于另外一肢的上一步初始速度。本文假定這個目標能實現，不討論如何實現這個目標。

2 模型行走

區別于傳統結構動力學已知外界激勵求解結構響應的問題，行走模型運動方程（1）中驅動力矩u并不能預先知道，并且隨意假定u也不能使模型完成擬人行走，甚至可能導致數值模擬失敗。為了克服這個困難，本文提出如下方法：對驅動自由度qa在整個擺動階段的時程進行規劃（在文中稱為規劃軌跡qa，d），代入式（1）可得對應此規劃軌跡驅動力u的時程（此后稱為規劃力矩ud）；考慮人與結構相互作用以及人運動過程中可能受到的外部擾動（例如結構振動、地震作用等），模型實際軌跡qa與規劃軌跡qa，d之間可能產生誤差，因此需要反饋線性化控制方法對驅動力矩u進行動態調整，以此完成模型行走仿真。

2.1 步態規劃

為使多剛體模型進行擬人行走，首先需對系統進行步態規劃。本文使用勞斯規約法進行步態規劃。勞斯規約法是一種幾何降維法，針對對稱性循環動力系統，采用單個運動周期內具有單調特性的循環變量作為系統基本變量，降低系統維度以達到簡化求解目的［14］。雙足行走動力系統進行向前行走行為時，支撐腿與地面間夾角θ在單行走周期內具有單調遞減性，如圖1所示。利用此特性，將平面內行走的多變量運動系統用θ進行描述。

對于模型平面步行的擺動階段，將驅動自由度qa時程假定為變量θ的函數，通常采用如下形式的貝塞爾函數：

與初始定義的qa相區分，式中d 表示規劃量；αa，k為規劃系數，控制軌跡形狀。以q2，d為例，q3，d，q4，d形式可類比：

式中s=0 和s=1 對應模型周期的初始和終止時刻，對qa的邊界條件有：

由式（4）～（7）可得：

將驅動自由度規劃軌跡qa，d代入式（1）后，可得到關于變量θ的單自由度運動方程：

注意，以上方程不僅需要滿足初始條件，即：

還須滿足模型單個行走周期的終止條件：

式中tf為常微分方程（9）求解的終止時刻，其判斷依據為擺動腿末端與地面發生接觸。

2.2 單步行走狀態

擺動腿末端判定與地面接觸后，多剛體模型當前行走周期結束，模型碰撞前的狀態即為此行走周期的終止狀態擺動腿觸地判定條件為ys為擺動腿末端的y坐標（剛體4 端點）。如圖2所示，擺動腿觸地后，支撐腿與擺動腿互換，相應的位形坐標轉換關系可表示為：

圖2 前后步位形切換Fig.2 Shape between front and rear steps

其中：

假定人體以穩定的周期步態行走，即下一周期內運動軌跡與上一周期的軌跡相同，各個行走周期的起始狀態間也存在對應關系，可表述為：

由式（12）和（14）可得單個行走周期內終止位形和初始位形的關系：

單個行走周期內，模型的起始位形與終止位形需滿足式（15）。由于變量qa與θ間引入式（3）規劃的軌跡關系，原四剛體行走模型微分方程組簡化為僅關于θ的常微分方程（9）。求解式（9）可得欠驅動角度θ的軌跡，歸一化后代入式（3）得到此解對應驅動自由度規劃軌跡。根據貝塞爾曲線性質，驅動自由度qa初始與終止條件能自動滿足，無需進行驗證。但是，qa的軌跡與欠驅動初始和終止狀態相關。

變量θ的初始與終止位形須滿足式（15），本文假定已知欠驅動初始角度，按式（15）確定欠驅動終止角度。將和列為未知量，以按方程（9）計算得到的終止值與預設終止值相等為目標進行迭代，從而得到滿足初始和終止邊界條件的欠驅動角位移時程，進而根據貝塞爾函數式（3）得到驅動自由度角位移的規劃時程qa，d。

至此欠驅動行走步態規劃可轉化為如下問題。在給定的[θi，θf]T下，尋找特定的關于θ的起始與終止速度量，使方程（9）存在滿足邊界條件的可行解，并能同時滿足一系列關于行走的不等式約束：包括欠驅動角速度始終小于零，即<0；以及支撐足與地面不發生滑動，即|tanθ| <μ，上述約束能保證多剛體模型行走時符合日常觀測。

2.3 反饋線性化控制

人體在靜止地面上行走時，若外界環境不存在擾動，則其能穩定地進行周期性行走。但實際上，人體行走時存在受到外界擾動的可能性，尤其是在人?結構相互作用體系中，結構的動力響應會影響到人體行走行為。在這種情形下，動力方程組（1）右側存在擾動力項，即關于欠驅動自由度θ的動力方程（9）右側存在擾動力。對人體行走這一動力系統引入控制器，動態調整系統的輸入扭矩u，以減小各驅動自由度與預設的規劃軌跡間因外部擾動導致的誤差。

如圖3所示為模型控制器的設計思路，動力系統的輸出為誤差函數y=qa-qa，d。由于輸出y與輸入u間的強非線性，本文使用反饋線性化方法，將u設計為由規劃軌跡相關的前饋項，以及與y相關的反饋項構成，使系統的實際輸出與規劃軌跡的誤差y，關于y的一階、二階導數均為零，詳見文獻［15］。

圖3 控制器結構Fig.3 Controller structure

2.4 人-結構相互作用

本文的目的是建立人體行走的多剛體模型，人與結構相互作用并不是本文的研究重點，但本文的基本方法可以直接考慮人與結構的相互作用。下面僅從概念上對人與結構相互作用進行討論。以地面為絕對坐標系建立動力學方程，可得到人與結構耦合的運動方程：

式中Ms，Cs，Ks，分別表示結構質量、阻尼和剛度矩陣列；x表示人致荷載作用下結構動力響應；Fh為模型行走足底力，作用點為人體與結構接觸處；足底力的反作用力轉化為廣義力后變為零，即模型運動方程右側第二項為零。結構發生振動時，人體模型相對坐標系原點運動，模型作為運動于結構上物體，各質心處存在由結構產生的慣性力，Fi=[Fxi，Fyi]T表示模型所受慣性力，i=1，2，…，5，更具體的可見文獻［13］。結構的振動將會影響到模型的實際行走軌跡，同時會改變四剛體模型對樓面的作用力，體現了人?結構間的相互作用。

3 模型仿真

步行是人通過下肢雙足的交替動作移動機體的活動。步態分析利用生理學知識和力學概念等，對不同人步行時的姿態、行為特征進行分析對比。常用的步態分析參數包括行走時間參數如行走頻率、行走時間、雙足支撐期占比；時空參數如下肢關節運動角度時程，行走步長；動力學參數如人行激勵力、關節輸入扭矩等。通常針對人與結構行為進行力學分析時，更關注人致荷載在建筑結構中可能引起的共振現象，對人行走頻率更為敏感。相應地，在建立荷載模型時，其模型參數也表達為以人行頻率作為自變量的函數。因此，本文進行模型行走軌跡和力學性能的結果驗證時，同樣以人行頻率作為自變量。

對四剛體雙足行走模型進行仿真驗證，行走模型參數如表1所示。設定模型質量參數時，參考中國男性青年人體體段質心位置［16］，以及中國男性青年質量分布［17］。為保持模型質心位置與正常人體一致，參考二自由度行走模型［10］，將大小腿長度定為0.6 m，模型行走時質心位于1.02 m 附近，詳見圖1。表2展示步速為1.12 m/s 時，一種可行的初始和終止狀態值。

表1 模型物理參數Tab.1 Parameters of four rigid body model

表2 步速1.12 m/s 的初始和終止狀態值Tab.2 Initial and terminal state values of walking speed at 1.12 m/s

3.1 模型仿真步態對比

圖4為多剛體模型在其行走周期內的運動軌跡。為與人體步行軌跡進行對比分析，圖中展示了多剛體模型兩個行走周期?？梢钥闯?，與實際人行走時運動軌跡相比，多剛體模型能夠較好地模擬人體的運動特征：右膝屈曲角度增加至最大值，左腿完成抬腿動作；右膝屈曲角度逐漸減小，此時模型質心位于右腿上方；在右腿站立相后期，大腿與小腿間相對角度接近于零，擺動腿左腿接近觸地；擺動腿切換后，右腿邁步相開始，其膝關節屈曲角度變大，并加速向前擺動；右肢達到膝關節最大屈曲位置，此時擺動腿右腿剛與支撐腿左腿發生交錯；擺動相后期，髖關節達到最大屈曲角度，右腿觸地。

圖4 逐幀人體模型行走圖Fig.4 Schematic drawing of human walking during one cycle frame by frame

圖5為各自由度在單個行走周期內的時程，與事先約定的相同，欠驅動角度θ在仿真時呈現出單調性；q2形狀類似于正弦曲線；q3為兩大腿間相對角度，在行走時，擺動腿向前，則q3角逐步增大，在觸地前擺動腿回落，此時q3角減??；q4始終小于零。

圖5 各個自由度q 在一個行走周期內的時程曲線Fig.5 Time history for rotation angle of each degreeof-freedom q given a walking cycle

建模時對人體步態進行若干簡化，實際上多剛體模型的行走在忽略了左右腿區別后，周期為實際人體行走周期的一半。圖6展示了通過捕捉髖關節、膝關節處標記點運動軌跡，利用圖像識別技術得到支撐腿的完整周期角度時程。對比多剛體模型行走時的時程曲線，證實了模型行走軌跡的擬人性。

圖6 周期行走角度時程對比Fig.6 Periodic walking angle comparison

上述表達為針對人行走步態的主觀描述，表3為模型步態量化指標。以一名青年男性為受試對象，在平整且大小合適的鞋底放置測力鞋墊，對人體周期行走時的力學性能進行測試，更詳細的方案見文獻［18］。在步態參數的描述中，步長、步頻、步速滿足兩種數據即說明當前步態與對應實測人行走步態擬合較好。荷載峰值因子由荷載時程峰值與人體自重計算得到，表征模型力學性能，觀察到剛體模型與實測和荷載模型在相同頻率情況下吻合較好。

表3 步態參數對比Tab.3 Gait parameter comparison

3.2 模型足底力驗證

許多學者嘗試直接對人行足底力進行測量，本文采用文獻［19］中由中國人體步行數據擬合得到的傅里葉級數荷載模型，來表示單步步行豎向荷載曲線：

式中n為模型階數；αn，φn稱為第n階動載因子和相位角，詳見文獻［19］。

如圖7所示為在給定四種不同步行頻率的情況下，單步行走的豎向荷載規則化時程曲線的對比?？梢杂^察到，剛體模型的步頻與步長能較好地與實測模型吻合，且豎向荷載曲線呈現出“M”型，與實測差距較小，總體在可接受范圍內。豎向荷載峰值與兩類模型對比吻合較好，但后續的峰谷及第二峰值有一定差距。以圖7（a）為例，對比實測值數據，豎向荷載峰值相差在10%以內，圖中峰谷相差25%，第二峰值相差5%。