999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

踝關(guān)節(jié)外骨骼人機(jī)耦合動(dòng)力學(xué)與助力性能分析1)

2023-01-15 12:32:50高鈺清靳葳徐鑒方虹斌
力學(xué)學(xué)報(bào) 2022年12期
關(guān)鍵詞:助力模型

高鈺清 靳葳 徐鑒 方虹斌

(復(fù)旦大學(xué)智能機(jī)器人研究院,上海 200433)

(復(fù)旦大學(xué)智能機(jī)器人教育部工程研究中心,上海 200433)

(復(fù)旦大學(xué)上海智能機(jī)器人工程技術(shù)研究中心,上海 200433)

引言

根據(jù)《第七次全國(guó)人口普查公報(bào)(第五號(hào))》的結(jié)果[1],我國(guó)60歲以上的人口數(shù)約為2.64 億人(2021 年5 月11 日),占全國(guó)人口的18.70%.由于人體衰老,老年人在日常行走和上下樓梯等場(chǎng)景會(huì)出現(xiàn)乏力癥狀,這使得老年人對(duì)助力設(shè)備的需求明顯增加.另一方面,士兵常需要攜帶武器裝備等負(fù)重進(jìn)行長(zhǎng)距離行軍,體力消耗巨大,迫切需要配備有效的行軍助力裝備.為滿足上述需求,近20年來,下肢助力外骨骼的設(shè)計(jì)、分析和性能優(yōu)化研究受到了廣泛關(guān)注.在下肢運(yùn)動(dòng)的過程中,踝關(guān)節(jié)為行走提供了最大的關(guān)節(jié)力矩[2],其在人體行走時(shí)輸出的正功率約為髖膝關(guān)節(jié)的總和[3].因此,踝關(guān)節(jié)外骨骼成為了相關(guān)研究的重點(diǎn)之一.

下肢踝關(guān)節(jié)外骨骼從驅(qū)動(dòng)方式上可以分為被動(dòng)外骨骼和主動(dòng)外骨骼兩種.被動(dòng)踝關(guān)節(jié)外骨骼具有結(jié)構(gòu)簡(jiǎn)單、輕量便攜等優(yōu)點(diǎn),典型代表包括由卡內(nèi)基梅隆大學(xué)Collins等[4]研發(fā)的無動(dòng)力踝關(guān)節(jié)外骨骼和范德堡大學(xué)Yandell等[5]研發(fā)的易于日常穿戴的被動(dòng)踝關(guān)節(jié)外骨骼.主動(dòng)踝關(guān)節(jié)外骨骼則具有更好的助力效果,典型代表包括由哈佛大學(xué)Walsh等[6-8]研制的繩驅(qū)動(dòng)外骨骼,加拿大女王大學(xué)Shepertycky等[9]研制的踝關(guān)節(jié)助力外骨骼,以及美國(guó)波士頓大學(xué)Aawd等[10]研制的踝關(guān)節(jié)助力外骨骼等.近幾年研究表明,在性能方面,被動(dòng)踝關(guān)節(jié)外骨骼的輔助效果不佳,而主動(dòng)踝關(guān)節(jié)外骨骼優(yōu)化了能量效率,具有更好的運(yùn)動(dòng)性能,值得進(jìn)一步研究[11].

為更好地了解踝關(guān)節(jié)外骨骼對(duì)人體行走的助力效果,需要從動(dòng)力學(xué)角度對(duì)人體-外骨骼耦合系統(tǒng)進(jìn)行建模和分析.目前的研究常基于外骨骼與人體下肢構(gòu)成的幾何關(guān)系,推導(dǎo)出外骨骼需要提供的力或力矩,通過對(duì)人體步態(tài)進(jìn)行劃分,形成外骨骼助力曲線,以達(dá)到根據(jù)人體步態(tài)周期助力的目的[12-14].但這種幾何關(guān)系推導(dǎo)并未從動(dòng)力學(xué)的角度解釋下肢踝關(guān)節(jié)外骨骼對(duì)人體行走的助力機(jī)制和效果.另一方面,針對(duì)機(jī)械臂或雙足、四足機(jī)器人,常基于拉格朗日方程進(jìn)行動(dòng)力學(xué)建模[15-16].這種建模方法能夠從能量角度推導(dǎo)出機(jī)器人各關(guān)節(jié)運(yùn)動(dòng)與關(guān)節(jié)力矩、外部力的關(guān)系,尤其適用于多自由度的復(fù)雜系統(tǒng).對(duì)于外骨骼與人體構(gòu)成的耦合系統(tǒng),由于自由度數(shù)高,主動(dòng)力和足-地交互作用力等因素復(fù)雜,相關(guān)的耦合動(dòng)力學(xué)建模以及基于動(dòng)力學(xué)分析的助力效果評(píng)估研究還不夠充分,這使得外骨骼的設(shè)計(jì)和研制缺乏動(dòng)力學(xué)依據(jù),無法事先對(duì)外骨骼設(shè)計(jì)進(jìn)行動(dòng)力學(xué)優(yōu)化.上述問題,要求建立包含外骨骼控制力和足-地交互作用力等因素的人體-外骨骼耦合動(dòng)力學(xué)模型,并以動(dòng)力學(xué)和生物力學(xué)量為指標(biāo),從多個(gè)角度分析踝關(guān)節(jié)外骨骼對(duì)人體行走的助力效果.

穿戴踝關(guān)節(jié)外骨骼的人體在行走過程中,人體自身動(dòng)力學(xué)、外骨骼施加給人體下肢的主動(dòng)助力以及足-地交互作用力這三方面因素對(duì)人體下肢-外骨骼耦合系統(tǒng)的動(dòng)力學(xué)行為具有重要影響,需要在建模過程中予以充分考慮.描述人體動(dòng)力學(xué)的模型主要包括[17]:骨骼模型、肌肉骨骼模型和神經(jīng)肌肉骨骼模型.對(duì)于骨骼模型,常根據(jù)捕獲的運(yùn)動(dòng)學(xué)數(shù)據(jù)和實(shí)測(cè)的足地交互力數(shù)據(jù),開展逆動(dòng)力學(xué)求解關(guān)節(jié)力矩[18-19],其可作為正向動(dòng)力學(xué)的輸入來估計(jì)人體運(yùn)動(dòng)[20-21],也有研究通過混合方法從正向動(dòng)力學(xué)仿真中得到相關(guān)數(shù)據(jù)來更新關(guān)節(jié)力矩[22].肌肉骨骼模型和神經(jīng)肌肉骨骼模型則考慮了神經(jīng)和肌肉對(duì)人體運(yùn)動(dòng)的影響[23-24].本文關(guān)注的是外骨骼對(duì)人體運(yùn)動(dòng)和關(guān)節(jié)力矩產(chǎn)生的影響,因此在理論建模上可暫不考慮肌肉和神經(jīng)的作用.下肢外骨骼對(duì)人體的主動(dòng)助力常通過繩驅(qū)[6-8]或氣動(dòng)[25]實(shí)現(xiàn),并以繩驅(qū)更為常見[26].足-地交互作用力包括地面反作用力和地面對(duì)足的摩擦力.常用的描述地面反作用力的模型有:赫茲接觸力模型[27]、Kelvin-Voigt 模型[28-29]、Hunt-Crossley模型[30]等;常用的描述摩擦力的模型有:庫倫摩擦模型[31]、Karnopp 摩擦模型[32]、Stribeck 摩擦模型[33]以及LuGre 摩擦模型[34]等.注意到,由Kelvin-Voigt模型和庫倫摩擦模型描述的足地交互作用力已廣泛應(yīng)用于多足機(jī)器人行走的動(dòng)力學(xué)分析[29,35-38],被證明可以較為準(zhǔn)確地表征地面反作用力和摩擦力.盡管上述三方面的研究各自已取得長(zhǎng)足發(fā)展,但綜合考慮這三方面因素的人體-外骨骼耦合動(dòng)力學(xué)研究才剛剛起步,相關(guān)研究水平仍然較低.

綜上所述,考慮人體動(dòng)力學(xué)、主動(dòng)助力和足-地交互作用力的人體-外骨骼耦合動(dòng)力學(xué)的研究尚處于起步階段,踝關(guān)節(jié)外骨骼對(duì)人體行走的助力效果尚未從動(dòng)力學(xué)角度得到深入分析,其難點(diǎn)包括人體-外骨骼耦合動(dòng)力學(xué)建模、主動(dòng)助力曲線生成以及助力效果的表征等.針對(duì)這些問題,本文將首先融合機(jī)器人正運(yùn)動(dòng)學(xué)描述方法和第二類拉格朗日方程,建立人體下肢-踝關(guān)節(jié)外骨骼耦合動(dòng)力學(xué)模型,其中足-地交互作用力由Kelvin-Voigt 模型和庫倫摩擦模型予以描述,人體關(guān)節(jié)力矩由PD 軌跡跟蹤控制生成并通過粒子群算法進(jìn)行控制參數(shù)調(diào)節(jié),主動(dòng)助力力矩通過上層控制器根據(jù)人體行走時(shí)步態(tài)周期的變化進(jìn)行實(shí)時(shí)調(diào)整.以此模型為基礎(chǔ)進(jìn)行動(dòng)力學(xué)數(shù)值計(jì)算,并以踝關(guān)節(jié)角度、踝關(guān)節(jié)力矩、踝關(guān)節(jié)功率和踝關(guān)節(jié)做功為評(píng)價(jià)指標(biāo),對(duì)踝關(guān)節(jié)外骨骼的助力效果進(jìn)行了分析和評(píng)價(jià).最后,本文在SCONE平臺(tái)中對(duì)人體-踝關(guān)節(jié)外骨骼耦合系統(tǒng)進(jìn)行肌肉骨骼動(dòng)力學(xué)預(yù)測(cè)仿真,并以肌肉激活度為指標(biāo),對(duì)外骨骼的助力效果進(jìn)行驗(yàn)證.

本文的創(chuàng)新貢獻(xiàn)在于:首先,進(jìn)一步完善了人體下肢-外骨骼耦合系統(tǒng)的動(dòng)力學(xué)建模方法,重點(diǎn)關(guān)注了足-地交互力、人體關(guān)節(jié)力矩和外骨骼助力力矩的生成方法,為人機(jī)耦合動(dòng)力學(xué)分析奠定了模型基礎(chǔ);其次,本文基于動(dòng)力學(xué)模型和肌肉骨骼模型分別對(duì)下肢踝關(guān)節(jié)外骨骼的助力效果進(jìn)行了分析和驗(yàn)證,通過關(guān)節(jié)角度、關(guān)節(jié)力矩、關(guān)節(jié)做功、肌肉激活度等多樣化指標(biāo),明確了踝關(guān)節(jié)外骨骼對(duì)行走的助力機(jī)制,提高了結(jié)論的可靠性.

1 人體下肢-踝關(guān)節(jié)外骨骼耦合動(dòng)力學(xué)模型

本研究所關(guān)注的繩驅(qū)踝關(guān)節(jié)外骨骼(圖1)的根本機(jī)理是模仿小腿肌肉為人體行走提供助力.由于足部生理結(jié)構(gòu)復(fù)雜,外骨骼提供的助力會(huì)對(duì)關(guān)節(jié)運(yùn)動(dòng)產(chǎn)生明顯影響,錯(cuò)誤的助力時(shí)機(jī)及大小均會(huì)降低外骨骼穿戴的舒適性和有效性,甚至影響人的正常行走.為此,有必要對(duì)人體下肢和踝關(guān)節(jié)外骨骼耦合系統(tǒng)進(jìn)行動(dòng)力學(xué)建模,并據(jù)此分析外骨骼對(duì)行走的助力效果.為此,本文首先介紹踝關(guān)節(jié)外骨骼系統(tǒng)及人體運(yùn)動(dòng)學(xué)模型,給出足-地交互力模型,并結(jié)合機(jī)器人正運(yùn)動(dòng)學(xué)描述手段建立人體-外骨骼耦合動(dòng)力學(xué)模型.

圖1 人體-踝關(guān)節(jié)外骨骼耦合系統(tǒng)及人體下肢運(yùn)動(dòng)學(xué)模型Fig.1 Human-ankle exoskeleton coupled system and the kinematic model of the human lower extremity

圖1 人體-踝關(guān)節(jié)外骨骼耦合系統(tǒng)及人體下肢運(yùn)動(dòng)學(xué)模型(續(xù))Fig.1 Human-ankle exoskeleton coupled system and the kinematic model of the human lower extremity(continued)

1.1 踝關(guān)節(jié)外骨骼系統(tǒng)及人體運(yùn)動(dòng)學(xué)模型

圖1(a)展示了人體-踝關(guān)節(jié)外骨骼耦合系統(tǒng),其中人體背負(fù)電機(jī)背包,背包內(nèi)包含外骨骼驅(qū)動(dòng)電機(jī)、電源和電子元件,從背包內(nèi)引出的兩條鮑登線連接于踝關(guān)節(jié)外骨骼的支撐模塊.為保證助力效果并減少外骨骼對(duì)人體的負(fù)面影響,工作過程中盡量保證外骨骼鮑登線的助力方向與小腿平行.

人體的運(yùn)動(dòng)在三維空間中可劃分為三個(gè)面:矢狀面、冠狀面和橫斷面.人體行走活動(dòng)主要發(fā)生在矢狀面內(nèi),冠狀面和水平面內(nèi)的運(yùn)動(dòng)較小,同時(shí)外骨骼的助力也僅發(fā)生在矢狀面內(nèi),因此冠狀面和橫斷面內(nèi)的運(yùn)動(dòng)在本研究中可以忽略.

在人體矢狀面內(nèi),下肢可以被視為一系列通過關(guān)節(jié)連接在一起的剛性桿,本文將這些剛性桿統(tǒng)稱為體節(jié).圖1(b)為本文所建立的人體運(yùn)動(dòng)學(xué)模型,坐標(biāo)系按Denavit-Hartenberg(DH)方法建立.為保證人體正常向前運(yùn)動(dòng),基坐標(biāo)系 {W} 設(shè)定在地面固定不動(dòng),同時(shí)在人體髖關(guān)節(jié)處建立坐標(biāo)系 {0},X0和Y0用來描述人體的向前運(yùn)動(dòng)和重心起伏.其他體節(jié)坐標(biāo)系以體節(jié)的起始點(diǎn)作為坐標(biāo)系原點(diǎn),Z軸垂直紙面向外,沿體節(jié)方向?yàn)閄軸正方向,并用右手法則確定Y軸.

具體地,體節(jié)1為軀干,體節(jié)2,3,4,以及5,6,7 分別為左側(cè)和右側(cè)的大腿、小腿和腳體節(jié),其對(duì)應(yīng)的DH 參數(shù)如表1 所示.其中,i=1,2,···,7為關(guān)節(jié)編號(hào),連桿長(zhǎng)度ai-1為沿Xi-1軸,從Zi-1移動(dòng)到Zi的距離;連桿扭角 αi-1為繞Xi-1軸,從Zi-1旋轉(zhuǎn)到Zi的角度;連桿偏距di為沿Zi軸,從Xi-1移動(dòng)到Xi的距離;關(guān)節(jié)角 θi為繞Zi軸,從Xi-1旋轉(zhuǎn)到Xi的角度.這里,約定 θ1=θbody,θ2=θleft_hip,θ3=θleft_knee,θ4=θleft_ankle,θ5=θright_hip,θ6=θright_knee,θ7=θright_ankle.

表1 人體體節(jié)DH 參數(shù)Table 1 DH parameters of the human body segments

根據(jù)DH 參數(shù)表可以求得每一個(gè)體節(jié)相對(duì)于基坐標(biāo)系 {W} 的齊次變換矩陣.以軀干為例

其中,X0和Y0表示髖關(guān)節(jié)(或坐標(biāo)系{0}的原點(diǎn))在基坐標(biāo)系W下的位置.

1.2 足-地交互力模型

足與地面的交互力主要包括地面反作用力以及摩擦力兩部分.為了簡(jiǎn)化計(jì)算,本文假設(shè)在人體行走的過程中,只有腳跟和腳尖兩點(diǎn)與地面存在接觸力和摩擦力.

地面反作用力模型由Kelvin-Voigt 模型來描述[29]

式中,Fdjn表示地面對(duì)腳的支持力,方向垂直于地面向上;下標(biāo)“d”取r 或l,分別表示右腳和左腳;“j”取h 或t,分別表示腳跟和腳尖,“n”表示反作用力;Kdjn為接觸剛度,Cdjn為接觸阻尼,δdj表示腳跟或者腳尖與地面間的接觸深度.

地面對(duì)足的摩擦力采用光滑化的庫倫摩擦模型描述[29]

式中,Fdjf表示摩擦力,方向沿水平方向;類似地,下標(biāo)“d”區(qū)分左右腳,下標(biāo)“j”區(qū)分腳跟和腳尖,“f ”表示摩擦力;C為常數(shù),表示摩擦力的光滑化程度;vj表示腳跟或腳尖相對(duì)地面的切向速度.

1.3 人體-外骨骼耦合動(dòng)力學(xué)模型

根據(jù)圖1(b)建立的運(yùn)動(dòng)學(xué)模型,人體體節(jié)的相對(duì)運(yùn)動(dòng)可以由7 個(gè)自由度 θi(i=1,2,···,7) 描述,人體的向前運(yùn)動(dòng)和重心起伏由X0,Y0兩個(gè)自由度描述.因此,本文描述的人體動(dòng)力學(xué)模型包含9 個(gè)自由度.在本研究中,考慮到繩驅(qū)外骨骼質(zhì)量遠(yuǎn)小于人體體節(jié)質(zhì)量,外骨骼本身的質(zhì)量予以忽略;外骨骼與人體的相互作用被簡(jiǎn)化為單向理想助力.圖1(b)和圖2 共同給出了人體-外骨骼耦合動(dòng)力學(xué)模型.結(jié)合機(jī)器人正運(yùn)動(dòng)學(xué)并基于拉格朗日方程,可以推導(dǎo)出模型的動(dòng)力學(xué)方程.

圖2 人體-外骨骼耦合動(dòng)力學(xué)模型Fig.2 Dynamics model of the human-exoskeleton system

具體地,人體任意體節(jié)上任意一點(diǎn)P相對(duì)于基坐標(biāo)系的位置矢量可以通過齊次變換矩陣表示

式中,wrp為點(diǎn)P在基坐標(biāo)系下的位置矢量,為坐標(biāo)系 {i}(i=1,2,···,7) 相對(duì)于基坐標(biāo)系的位姿變換矩陣,irP為點(diǎn)P在體節(jié)坐標(biāo)系 {i} 中的保持不變的位置矢量.由此,體節(jié)上任意一點(diǎn)P的速度可以表示為

速度的平方可寫為

其中,tr()表示矩陣的跡.設(shè)體節(jié)上任意質(zhì)點(diǎn)P的質(zhì)量為 dm,則該質(zhì)點(diǎn)的動(dòng)能為

積分得到任意體節(jié)的動(dòng)能可寫為

式中,Ii為偽慣量矩陣,其一般表達(dá)式為

在xoy平面內(nèi),假設(shè)慣量積為0,可以得到對(duì)應(yīng)每一個(gè)體節(jié)的偽慣量矩陣.以左足模型為例,其偽慣量矩陣由式(10)給出

對(duì)所有體節(jié)的動(dòng)能求和,可以得到總動(dòng)能為

各體節(jié)的勢(shì)能為

式中,mi是體節(jié)i的質(zhì)量,gT是重力行矢量

irci表示體節(jié)i的質(zhì)心在 {i} 坐標(biāo)系中的位置矢量.對(duì)所有體節(jié)的勢(shì)能求和,可以得到總勢(shì)能為

根據(jù)式(11)和式(14),系統(tǒng)的拉格朗日函數(shù)L可以表示為

將L對(duì)廣義坐標(biāo) θi以及X0,Y0求導(dǎo),可以得到對(duì)應(yīng)于廣義坐標(biāo) θi的動(dòng)力學(xué)方程

其中,θi(i=1,2,···,9) 由表1 給出.令X0=θ8,Y0=θ9,并將式(15)代入到式(16)計(jì)算,可得

將式(17)整理成矩陣形式,可得

式中,τhuman為人體關(guān)節(jié)力矩,由基于粒子群優(yōu)化算法的PD 控制給出(在2.2 節(jié)中予以具體闡述);τexo為外骨骼助力力矩,由上層控制器根據(jù)人體步態(tài)周期給出,其由外骨骼通過鮑登線施加于支撐模塊的力Fexo乘以作用力臂得到(在2.3 節(jié)中予以具體闡述); τf為足-地交互力所對(duì)應(yīng)的廣義力矢量,具體可以表示為

其中,Frhn,Flhn,Frtn,Fltn分別為右、左腳跟和右、左腳尖的地面支持力,由式(2) 具體給出;Frhf,Flhf,Frtf,Fltf分別為右、左腳跟與地面的摩擦力,和右、左腳尖與地面的摩擦力,由式(3)具體給出.Jf為將足-地交互力變換為廣義力矩的雅可比矩陣.和Jf的具體形式見附錄.

至此,本文得到了考慮足-地交互力、人體關(guān)節(jié)力矩和主動(dòng)助力的人體-外骨骼耦合動(dòng)力學(xué)模型.后續(xù)將基于該模型進(jìn)行動(dòng)力學(xué)分析.

2 人體-外骨骼耦合系統(tǒng)的動(dòng)力學(xué)分析

在上述動(dòng)力學(xué)模型的基礎(chǔ)上對(duì)人體-外骨骼耦合系統(tǒng)開展動(dòng)力學(xué)分析.數(shù)值計(jì)算基于Matlab R2020b 開展,其中用到了前期工作[29]中健康人體行走的運(yùn)動(dòng)學(xué)數(shù)據(jù)和人體參數(shù)(表2).

表2 人體節(jié)段長(zhǎng)度和重量值Table 2 Lengths and weights of human segments

2.1 人體步態(tài)階段劃分

人體運(yùn)動(dòng)研究中的關(guān)鍵一環(huán)是步態(tài)階段劃分.在外骨骼研究中,如果不能對(duì)步態(tài)階段進(jìn)行有效和準(zhǔn)確的劃分,外骨骼的助力效果將受到很大影響,甚至有可能對(duì)人體造成損傷[7].因此,在開展運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)分析之前,本研究首先對(duì)人體步態(tài)階段進(jìn)行劃分(圖3),其目的是有效設(shè)置外骨骼的助力時(shí)間與休眠時(shí)間,以期達(dá)到理想助力效果.

圖3 步態(tài)階段劃分Fig.3 Division of gait stages

圖3 展示了健康人行走時(shí)左側(cè)踝關(guān)節(jié)角度變化時(shí)間歷程[29].一個(gè)步態(tài)周期(one step)的劃分從單側(cè)腳跟著地開始,到下一次該側(cè)腳跟著地為止為一個(gè)步態(tài)周期(圖3 中淺色陰影區(qū)段).踝關(guān)節(jié)外骨骼在跖屈(plantar flexion)階段對(duì)人體進(jìn)行助力.踝關(guān)節(jié)跖屈的定義為腳后跟抬起然后腳尖向下的過程,在圖3 中用深色陰影表示.在步態(tài)階段劃分中踝關(guān)節(jié)跖屈的具體判斷方式為:腳跟離地瞬間(踝關(guān)節(jié)角度約為1.67 弧度)為跖屈開始時(shí)刻,小腿與足成90°時(shí)為踝關(guān)節(jié)跖屈的結(jié)束時(shí)刻.

2.2 基于粒子群優(yōu)化的人體關(guān)節(jié)控制力矩生成

在動(dòng)力學(xué)分析中,人體的關(guān)節(jié)力矩由PD 軌跡跟蹤控制實(shí)現(xiàn):即以人體期望軌跡與實(shí)際軌跡的差值為比例項(xiàng),使用差分近似導(dǎo)數(shù)構(gòu)成算法中的微分項(xiàng).具體地,基于PD 控制器的關(guān)節(jié)力矩可以表示為

其中,τj為關(guān)節(jié)j的力矩,Kj,P,Kj,D為比例和微分增益,ej(t)和分別表示了當(dāng)前時(shí)刻下關(guān)節(jié)j的位置誤差和速度誤差.

為體現(xiàn)人體對(duì)環(huán)境和外骨骼的自適應(yīng)性,PD 控制的參數(shù)Kj,P和Kj,D不采用固定值,而是通過具有出色搜索能力的啟發(fā)式算法自適應(yīng)地調(diào)節(jié).具體地,本研究選取粒子群優(yōu)化(PSO)算法來確定PD 參數(shù)以達(dá)到自適應(yīng)生成人體關(guān)節(jié)力矩的目的[39],這是因?yàn)槠淠軌蛟跊]有得知太多問題信息的情況下,有效地搜索龐大的解空間.

基于粒子群優(yōu)化算法和PD 軌跡跟蹤控制生成人體關(guān)節(jié)力矩的流程如圖4 所示,粒子群優(yōu)化算法的參數(shù)列于表3.其中,D為粒子維度,即待優(yōu)化PD控制參數(shù)的維度,前7 維為KP參數(shù),后7 維為KD參數(shù);m為種群規(guī)模,n為總迭代次數(shù),根據(jù)仿真所需時(shí)間給定; ω為慣性權(quán)重,體現(xiàn)了粒子保持前一運(yùn)動(dòng)狀態(tài)的能力;k為當(dāng)前迭代次數(shù),c1和c2為個(gè)體學(xué)習(xí)因子和群體學(xué)習(xí)因子,用于調(diào)節(jié)學(xué)習(xí)最大步長(zhǎng).

圖4 基于粒子群優(yōu)化的人體關(guān)節(jié)力矩PD 控制流程Fig.4 Block diagram of the PSO-based PD controller for human joint torques

表3 粒子群優(yōu)化算法參數(shù)值Table 3 Parameters of the PSO Algorithm

本研究中,動(dòng)力學(xué)數(shù)值計(jì)算通過歐拉法實(shí)現(xiàn),其計(jì)算步長(zhǎng)被設(shè)為0.1 ms,每個(gè)粒子的作用時(shí)間為2 ms.PSO 算法中速度更新和位置更新公式如下

其中,為第k次迭代中粒子i的速度矢量,為第k次迭代中粒子i的位置矢量,r1和r2為區(qū)間[0,1] 內(nèi)的隨機(jī)數(shù).Pi為粒子i的當(dāng)前個(gè)體最優(yōu)位置矢量,Pg為整個(gè)粒子群的當(dāng)前全局最優(yōu)位置矢量.r1,r2的引入增加了種群搜索的隨機(jī)性,提升了粒子群算法避免陷入局部最優(yōu)的能力.

粒子群優(yōu)化的損失函數(shù)定義為

其中,α和β為位置誤差和速度誤差的權(quán)重,損失函數(shù)表征了在ts時(shí)刻系統(tǒng)位置和速度的平均累計(jì)誤差.

基于圖4 的具體流程如下:根據(jù)每個(gè)粒子當(dāng)前的位置生成人體各個(gè)關(guān)節(jié)的力矩并作用于動(dòng)力學(xué)模型(18),每個(gè)粒子的作用時(shí)間為2 ms,計(jì)算當(dāng)前的損失函數(shù)f(ti) 并與之前的損失函數(shù)進(jìn)行比較,確定并更新該粒子的個(gè)體最優(yōu)位置Pi和粒子群的全局最優(yōu)位置Pg.當(dāng)m=100個(gè)粒子依次執(zhí)行完2 ms(共計(jì)0.2 s)后,基于當(dāng)前的個(gè)體最優(yōu)位置Pi和當(dāng)前全局最優(yōu)位置Pg,利用式(22)和式(23)更新每個(gè)粒子的位置和速度,并開始下一輪迭代.上述算法具有較好的收斂性和在線性,隨著迭代次數(shù)的增加,損失函數(shù)將很快收斂到最小值.

2.3 踝關(guān)節(jié)外骨骼助力曲線

為了取得理想的踝關(guān)節(jié)外骨骼控制效果,本文設(shè)計(jì)了外骨骼力矩控制器(即上層控制器),其通過人體步態(tài)周期確定踝關(guān)節(jié)外骨骼控制系統(tǒng)的期望助力力矩曲線,包括上升力矩曲線和下降力矩曲線.期望助力力矩曲線依據(jù)健康人行走時(shí)踝關(guān)節(jié)力矩曲線生成.

描述健康人踝關(guān)節(jié)力矩曲線的參數(shù)包括:start_time,fall_time,peak_time和peak_torque[39-40],其意義如圖5 所示.start_time和fall_time 分別表示踝關(guān)節(jié)力矩的上升時(shí)間和下降時(shí)間,peak_time 表示踝關(guān)節(jié)力矩到達(dá)峰值的時(shí)刻,peak_torque 表示踝關(guān)節(jié)的峰值力矩,它們常取值為37,45,15,90[40].首先,通過三次多項(xiàng)式擬合健康人踝關(guān)節(jié)力矩曲線.對(duì)于上升力矩曲線和下降力矩曲線,其兩端需分別滿足如下約束條件

圖5 踝關(guān)節(jié)外骨骼助力力矩曲線Fig.5 Assistive torque curve of the ankle exoskeleton

同時(shí),要求兩曲線在交點(diǎn)peak_time 處的斜率相等且為零,兩曲線在兩端的斜率也為0,即

求解以上約束可以給出健康人踝關(guān)節(jié)力矩的擬合曲線.圖5 表明,擬合曲線(點(diǎn)線)與真實(shí)踝關(guān)節(jié)力矩曲線(虛線)較為吻合.

為有效控制踝關(guān)節(jié)外骨骼,令實(shí)際施加的關(guān)節(jié)力矩為健康人踝關(guān)節(jié)力矩的70%,且假設(shè)左、右踝關(guān)節(jié)具有相同的助力力矩曲線.由此,可以得到外骨骼期望助力的上升力矩曲線和下降力矩曲線(圖5,實(shí)線),它們的表達(dá)式如下

擬合參數(shù)列于表4.為保證施加力矩后人體行走穩(wěn)定,本文規(guī)定,第一步和最后一步不施加外骨骼助力力矩.

表4 助力曲線的參數(shù)值Table 4 Parameters of the assistance curve

在本文基于動(dòng)力學(xué)模型的數(shù)值計(jì)算中,上述期望外骨骼力矩被直接施加于人體-外骨骼耦合動(dòng)力學(xué)模型.但在實(shí)際實(shí)驗(yàn)中,電機(jī)需要通過下層控制器跟蹤上層控制器的期望力矩曲線.這里,出于完整性考慮,對(duì)下層控制器算法進(jìn)行簡(jiǎn)述.在得到踝關(guān)節(jié)的期望助力力矩表達(dá)式后,可以依據(jù)所設(shè)計(jì)的外骨骼的幾何特征和電機(jī)的特性,計(jì)算電機(jī)的期望輸出轉(zhuǎn)矩,具體計(jì)算方法如下

其中,Fexo為電機(jī)轉(zhuǎn)軸處繩子的拉力,ra為關(guān)節(jié)支撐模塊的平面到踝關(guān)節(jié)軸的距離,η為繩子拉力的損耗系數(shù),rm為電機(jī)輪盤半徑,τm為電機(jī)的期望輸出轉(zhuǎn)矩.在仿真過程中,使用了由傳遞函數(shù)描述的電機(jī)模型,其描述了電機(jī)在零初始條件下的輸出量的拉普拉斯變換與輸入量的拉普拉斯變換之比[41].基于此模型可以計(jì)算出電機(jī)的實(shí)際輸出力矩,其與期望力矩的對(duì)比如圖6 所示.可以看出,期望力矩與實(shí)際輸出力矩差異并不明顯,表明PD 算法可以很好地跟蹤期望力矩.

圖6 踝關(guān)節(jié)外骨骼期望力矩與實(shí)際力矩對(duì)比Fig.6 Comparison of the actual torque with the desired torque of the ankle exoskeleton

2.4 踝關(guān)節(jié)外骨骼助力效果

首先根據(jù)所得到的運(yùn)動(dòng)學(xué)數(shù)據(jù)來評(píng)價(jià)外骨骼的助力效果[8].基于所建立的動(dòng)力學(xué)模型和健康人行走運(yùn)動(dòng)軌跡(行走速度2 km/h),得到了踝關(guān)節(jié)的實(shí)際運(yùn)動(dòng)軌跡.圖7 展示了左、右踝關(guān)節(jié)穿戴和未穿戴外骨骼兩種情況下的角度-時(shí)間歷程.圖7 表明,穿戴外骨骼顯著改變了左、右踝關(guān)節(jié)的角度規(guī)律.以左踝關(guān)節(jié)為例,從可以達(dá)到的最大角度來看,穿戴外骨骼要比未穿戴外骨骼的情況高約0.12 rad(6.88°);從最小角度來看,穿戴外骨骼也要比未穿戴外骨骼的情況高約0.17 rad(9.74°).由于得到的踝關(guān)節(jié)角度數(shù)據(jù)沒有呈現(xiàn)出明顯的分布規(guī)律,本文采用Wilcoxon 秩和檢驗(yàn)進(jìn)行差異顯著性分析.Wilcoxon秩和檢驗(yàn)屬于非參數(shù)檢驗(yàn)的一種,用于檢驗(yàn)兩個(gè)獨(dú)立的樣本是否符合同一分布,但是不要求被檢驗(yàn)的樣本具有相同的元素個(gè)數(shù),也不要求被檢驗(yàn)的樣本符合正態(tài)分布[42],因此,該檢驗(yàn)比較適合本文所提到的數(shù)據(jù).結(jié)果表明,穿戴外骨骼與未穿戴時(shí)況相比,左、右踝關(guān)節(jié)角度均具有顯著差異(p<0.01).造成這種顯著差異的原因是外骨骼在背屈過程中對(duì)踝關(guān)節(jié)施加了助力并牽引足部向上運(yùn)動(dòng),從而抬升了踝關(guān)節(jié)角度.

圖7 穿戴和未穿外骨骼情況下人體踝關(guān)節(jié)角度對(duì)比Fig.7 A comparison of human ankle angles for the cases with and without exoskeleton

除了運(yùn)動(dòng)學(xué)指標(biāo),也從人體踝關(guān)節(jié)力矩、踝關(guān)節(jié)功率和踝關(guān)節(jié)做功出發(fā),分析穿戴外骨骼對(duì)人體行走動(dòng)力學(xué)性能的影響.圖8 展示了左、右踝關(guān)節(jié)在穿戴和未穿戴外骨骼兩種情況下的力矩-時(shí)間歷程.從圖中結(jié)果可以看出,除了右腳邁出的第一步和左腳的最后一步由于沒有施加助力使得踝關(guān)節(jié)力矩?zé)o明顯變化,其他所有步態(tài)周期內(nèi),踝關(guān)節(jié)峰值力矩均有明顯的減少.Wilcoxon 秩和檢驗(yàn)(p<0.01)表明,穿戴外骨骼與未穿戴的情況相比,左、右踝關(guān)節(jié)力矩呈現(xiàn)出了顯著差異.從平均力矩的角度來看,左、右踝關(guān)節(jié)分別減少了40.76%和40.68%,充分展現(xiàn)了踝關(guān)節(jié)外骨骼的助力效果.注意到,左、右踝關(guān)節(jié)平均力矩減少程度有所不同,這一方面是由于健康人行走時(shí)左右踝關(guān)節(jié)角度的差異,另一方面粒子群優(yōu)化算法在生成人體關(guān)節(jié)力矩過程中具有一定的隨機(jī)性.

圖8 穿戴和未穿外骨骼情況下人體踝關(guān)節(jié)力矩對(duì)比Fig.8 A comparison of human ankle torques for the cases with and without exoskeleton

根據(jù)上述力矩結(jié)果,進(jìn)一步計(jì)算了左、右踝關(guān)節(jié)的功率,如圖9 所示.結(jié)果表明,穿戴外骨骼也可以顯著地降低兩側(cè)踝關(guān)節(jié)的峰值功率.此外,將踝關(guān)節(jié)功率對(duì)時(shí)間積分可以得到人體踝關(guān)節(jié)的做功.結(jié)果表明,穿戴外骨骼后,左、右踝關(guān)節(jié)在10s 內(nèi)的做功分別減少了37.47%和31.85%,實(shí)現(xiàn)了有效的行走助力.

圖9 穿戴和未穿外骨骼情況下人體踝關(guān)節(jié)功率對(duì)比Fig.9 A comparison of the human ankle power for the cases with and without exoskeleton

此外,為了展示踝關(guān)節(jié)外骨骼對(duì)不同行走狀態(tài)的適應(yīng)性,還針對(duì)5 種行走速度(2.5 km/h,3.5 km/h,4.5 km/h,5.5 km/h,6.5 km/h)分別進(jìn)行了動(dòng)力學(xué)仿真和助力效果分析.表5 給出了在6 種行走速度下,穿戴外骨骼后左、右踝關(guān)節(jié)平均力矩和做功的下降情況.結(jié)果表明,基于本文給出的助力力矩曲線設(shè)計(jì)方式,踝關(guān)節(jié)外骨骼對(duì)不同行走速度均具有較好的助力效果,可以實(shí)現(xiàn)至少24.84%的平均力矩降低和至少24.69%的做功降低.

表5 不同行走速度下外骨骼助力效果評(píng)估Table 5 Evaluation of exoskeleton assistance effect under different walking speed

上述分析表明,穿戴外骨骼后,人體踝關(guān)節(jié)角度變化規(guī)律仍能維持與健康人行走時(shí)一致,踝關(guān)節(jié)力矩、功率以及做功均實(shí)現(xiàn)了較大幅度的減小,說明本文提出的踝關(guān)節(jié)外骨骼及其助力方式是有效的,可以為人體行走提供助力.

3 基于SCONE 的肌肉-骨骼動(dòng)力學(xué)仿真

除了通過基于動(dòng)力學(xué)模型的數(shù)值計(jì)算展示踝關(guān)節(jié)外骨骼的助力效果,本研究還在SCONE 平臺(tái)中通過肌肉-骨骼動(dòng)力學(xué)模型預(yù)測(cè)仿真驗(yàn)證踝關(guān)節(jié)外骨骼的助力效果.以行走速度3.6 km/h為例進(jìn)行結(jié)果討論.

SCONE 是一款用于預(yù)測(cè)模擬人類、動(dòng)物和機(jī)器人運(yùn)動(dòng)的開源軟件.在SCONE 平臺(tái)中,執(zhí)行預(yù)測(cè)仿真所需要的一切都集成于一個(gè)場(chǎng)景中,其由以下幾部分組成:人體、動(dòng)物或機(jī)器人模型,一個(gè)為模型中的作動(dòng)器產(chǎn)生輸入信號(hào)的控制器,目標(biāo)任務(wù)以及優(yōu)化器.SCONE 使得肌肉骨骼模型有了強(qiáng)大的新應(yīng)用,如預(yù)測(cè)治療效果,優(yōu)化外骨骼設(shè)備的效率和功效等.更為重要的是,SCONE 使研究人員能夠研究模型和控制參數(shù)對(duì)整個(gè)運(yùn)動(dòng)的影響[43].本研究在SCONE 平臺(tái)中搭建了人體-外骨骼耦合系統(tǒng)的肌肉骨骼模型,并通過關(guān)鍵肌肉的激活程度來評(píng)價(jià)外骨骼對(duì)人體行走的助力效果.

本文搭建的人體-外骨骼耦合系統(tǒng)的肌骨模型如圖10所示,其由以下幾部分組成:電機(jī)背包、關(guān)節(jié)支撐模塊、人體骨骼模型、腿部肌肉,包括大腿肌肉、比目魚肌(soleus)和其他小腿肌肉.其中,比目魚肌受到重點(diǎn)關(guān)注,這是因?yàn)楸饶眶~肌的主要功能是作為踝關(guān)節(jié)的跖屈肌,它在走、跑、跳中對(duì)足的跖屈蹬地、固定踝關(guān)節(jié)和防止身體前傾等起到重要作用[44-46].模型在行走過程中足部受到地面反作用力,由Hunt-Crossley 模型予以描述并用黃色箭頭標(biāo)注,外骨骼同時(shí)還對(duì)人體踝關(guān)節(jié)施加由式(33)和式(34)給出的助力力矩.

圖10 SCONE 人體-外骨骼耦合系統(tǒng)的肌肉骨骼動(dòng)力學(xué)模型Fig.10 Musculoskeletal dynamic model of the human-exoskeleton coupled system in SCONE

在仿真過程中,通過評(píng)估SCONE 中目標(biāo)函數(shù)的收斂程度以及比目魚肌的激活程度來評(píng)價(jià)外骨骼的助力效果.具體地,肌肉骨骼動(dòng)力學(xué)模型的優(yōu)化目標(biāo)函數(shù)為

其中,fg是移動(dòng)能力指標(biāo),其限制步行速度和髖關(guān)節(jié)高度;fe為代謝耗能指標(biāo),評(píng)價(jià)行走過程中的代謝耗能;frf為地反力指標(biāo),避免行走過程中足部任一處的接觸力超過1.5 倍體重;flim_a和flim_k為約束指標(biāo),避免在行走過程中踝、膝關(guān)節(jié)的過伸或過屈.上述指標(biāo)的具體表示和解釋可見文獻(xiàn)[43],這里不做詳細(xì)解釋.wi為對(duì)應(yīng)指標(biāo)的權(quán)重,權(quán)重值設(shè)置如表6 所示.

表6 SCONE 中目標(biāo)函數(shù)的權(quán)重值Table 6 Weights of the objective function in SCONE

本文對(duì)無外骨骼助力(pure walking)和有外骨骼助力(exoskeleton on)兩種情況分別進(jìn)行了仿真預(yù)測(cè).圖11 給出了目標(biāo)函數(shù)收斂以及比目魚肌的激活程度的對(duì)比.結(jié)果表明,在SCONE 參數(shù)優(yōu)化過程中,式(37) 中的fg,frf,flim_a和flim_k四項(xiàng)在迭代至第20代左右時(shí)即收斂至零,這是因?yàn)檫@四項(xiàng)在目標(biāo)函數(shù)中起到加速收斂的作用以保證人體的正常、穩(wěn)定行走.在第20代以后,目標(biāo)函數(shù)的值由代謝耗能項(xiàng)fe主導(dǎo),其也逐漸收斂.圖11(a)表明,在第20代以后,穿戴外骨骼的情況要顯著優(yōu)于未穿戴外骨骼的情況,表現(xiàn)為更快的收斂速度和較低的目標(biāo)函數(shù)值.

圖11 SCONE 預(yù)測(cè)仿真結(jié)果Fig.11 Simulation results of SCONE

此外,圖11(b)表明,穿戴外骨骼可以降低比目魚肌激活度的峰值.Wilcoxon 秩和檢驗(yàn)也指出,比目魚肌的激活度在穿戴與未穿戴外骨骼情況下具有顯著差異(p<0.05).為了定量展示肌肉激活度,本文進(jìn)一步采用時(shí)域分析法來提取表面肌電信號(hào)特征,并采用RMS 指標(biāo)來量化肌肉激活度[32-34]

其中,N為采樣次數(shù),xi表示為一個(gè)歸一化的值,相當(dāng)于第i時(shí)刻的肌肉激活度(mV)除以肌肉可以產(chǎn)生的最大激活度.結(jié)果表明,穿戴和未穿戴外骨骼情況下比目魚肌的RMS 值分別為:0.097 8和0.104 3,即穿戴外骨骼實(shí)現(xiàn)了比目魚肌激活度約6.21%的下降.下降原因是外骨骼提供的力矩替代了部分比目魚肌的做功,減小了對(duì)比目魚肌的發(fā)力需求,從而降低了肌肉的激活度.

上述基于SCONE 平臺(tái)的肌肉骨骼模型和預(yù)測(cè)仿真分析從另一個(gè)角度驗(yàn)證了本文提出的繩驅(qū)踝關(guān)節(jié)外骨骼及其控制方式對(duì)人體行走的有效助力.通過評(píng)估比目魚肌的激活度以及肌電信號(hào)RMS 值,可以發(fā)現(xiàn)踝關(guān)節(jié)外骨骼對(duì)比目魚肌提供了一定的補(bǔ)償,從生理學(xué)角度解釋了踝關(guān)節(jié)外骨骼對(duì)人體行走助力的本質(zhì)機(jī)理.

4 結(jié)論

本文以自行設(shè)計(jì)的下肢繩驅(qū)踝關(guān)節(jié)外骨骼為研究對(duì)象,融合機(jī)器人正運(yùn)動(dòng)描述方法和拉格朗日方程,建立了人體-外骨骼耦合系統(tǒng)的人機(jī)動(dòng)力學(xué)模型.模型綜合考慮了以下因素:由Kelvin-Voigt 模型和庫倫摩擦模型描述的足-地交互力、由基于粒子群優(yōu)化和PD 控制生成的人體關(guān)節(jié)力和由上層控制器確定的外骨骼期望助力力矩.以此為基礎(chǔ)進(jìn)行動(dòng)力學(xué)數(shù)值計(jì)算,本文從人體踝關(guān)節(jié)角度、踝關(guān)節(jié)力矩、踝關(guān)節(jié)功率和踝關(guān)節(jié)做功多個(gè)角度出發(fā),系統(tǒng)分析了踝關(guān)節(jié)外骨骼的助力效果.本文也開展了基于SCONE 平臺(tái)的肌肉骨骼建模和預(yù)測(cè)仿真,從肌肉激活度的角度證實(shí)了踝關(guān)節(jié)外骨骼助力的有效性.本文得出的主要結(jié)論如下.

(1)基于動(dòng)力學(xué)模型的數(shù)值計(jì)算表明,穿戴外骨骼可以基本保持踝關(guān)節(jié)角度的變化規(guī)律,但會(huì)顯著影響踝關(guān)節(jié)的角度值(p<0.01).動(dòng)力學(xué)分析指出,穿戴外骨骼和未穿戴外骨骼情況下,人體行走動(dòng)力學(xué)性能具有顯著差異(p<0.01).在2.0km/h 到6.5 km/h 的人體行走步速下,穿戴外骨骼可以實(shí)現(xiàn)至少24.84%的人體踝關(guān)節(jié)的平均力矩下降和至少24.69%的踝關(guān)節(jié)做功降低.上述結(jié)果表明,該踝關(guān)節(jié)外骨骼可以對(duì)人體行走提供有效助力.

(2)基于SCONE 平臺(tái)的肌-骨動(dòng)力學(xué)預(yù)測(cè)仿真表明,在3.6 km/h 的步速下,比目魚肌的激活度在穿戴與未穿戴外骨骼情況下具有顯著差異(p<0.05),表現(xiàn)為激活度峰值的下降.定量地,穿戴外骨骼使得肌電信號(hào)的RMS 值下降了6.21%.上述結(jié)果從生理學(xué)的角度驗(yàn)證了本文提出的下肢踝關(guān)節(jié)和力矩控制器的助力效果.

綜上所述,本文的研究從動(dòng)力學(xué)角度入手,通過考慮足-地交互力、人體關(guān)節(jié)力矩以及外骨骼助力力矩,進(jìn)一步完善了人體下肢-外骨骼耦合系統(tǒng)的動(dòng)力學(xué)建模方法,為后續(xù)人機(jī)耦合動(dòng)力學(xué)分析奠定了基礎(chǔ).本文建立的動(dòng)力學(xué)模型也為后續(xù)踝關(guān)節(jié)外骨骼的實(shí)驗(yàn)研究提供了理論依據(jù),具有指導(dǎo)性意義.其次,本文通過動(dòng)力學(xué)數(shù)值計(jì)算和SCONE 肌肉骨骼動(dòng)力學(xué)預(yù)測(cè)仿真,從不同角度對(duì)踝關(guān)節(jié)外骨骼的助力效果進(jìn)行了分析和驗(yàn)證.在動(dòng)力學(xué)數(shù)值計(jì)算中,通過踝關(guān)節(jié)角度、踝關(guān)節(jié)力矩、踝關(guān)節(jié)功率和踝關(guān)節(jié)功率等指標(biāo),驗(yàn)證了踝關(guān)節(jié)外骨骼對(duì)人體行走助力的有效性;在SCONE 肌-骨動(dòng)力學(xué)預(yù)測(cè)仿真中,通過目標(biāo)函數(shù)的優(yōu)化結(jié)果以及肌肉激活度等指標(biāo),從生理學(xué)角度闡明了下肢踝關(guān)節(jié)外骨骼的助力機(jī)制.在后續(xù)工作中,將以本文的模型和方法為基礎(chǔ),進(jìn)一步開展下肢外骨骼的智能控制研究和實(shí)驗(yàn)研究.

附錄

對(duì)正文中的符號(hào)做下述簡(jiǎn)化

猜你喜歡
助力模型
一半模型
助力成功七件事
英語世界(2022年9期)2022-10-18 01:11:18
舉手之勞,助力我國(guó)實(shí)現(xiàn)碳中和
重要模型『一線三等角』
助力“一方水土養(yǎng)一方人”
金橋(2020年9期)2020-10-27 01:59:44
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
“自能”學(xué)習(xí)助力有機(jī)化學(xué)的學(xué)習(xí)
為更適宜的閱讀之地助力
商周刊(2017年17期)2017-09-08 13:08:58
健康助力 打贏脫貧攻堅(jiān)戰(zhàn)
3D打印中的模型分割與打包
主站蜘蛛池模板: 成人国产精品一级毛片天堂| 综合久久五月天| 亚洲人成网址| www.亚洲一区| 国产精品自在线拍国产电影| 99久久国产精品无码| 亚洲人精品亚洲人成在线| 久久伊人久久亚洲综合| 国产手机在线观看| 在线无码九区| 久久久黄色片| 国产av剧情无码精品色午夜| 88国产经典欧美一区二区三区| 欧美在线伊人| 亚洲色图欧美在线| 久久久久久久久18禁秘| 色视频国产| 国产精品美乳| 欧美亚洲国产一区| 狠狠亚洲婷婷综合色香| 亚洲人成网站日本片| 国产视频资源在线观看| 啪啪免费视频一区二区| 国产成人1024精品下载| 欧洲亚洲欧美国产日本高清| 亚洲欧美色中文字幕| 国产精品香蕉在线| 四虎成人在线视频| 亚洲成在线观看| 激情国产精品一区| 黄色网站不卡无码| 夜夜爽免费视频| 国产成人高清在线精品| 中文国产成人精品久久| 国产99视频精品免费视频7| 日本亚洲最大的色成网站www| 久久久久久尹人网香蕉| 国产精品成人免费综合| 亚洲永久色| 久久夜夜视频| 国内精品九九久久久精品| 狠狠久久综合伊人不卡| 爆操波多野结衣| 亚洲国产精品日韩av专区| 91成人精品视频| 亚洲av片在线免费观看| 激情六月丁香婷婷| 精品一区二区三区水蜜桃| 91精品视频在线播放| 人人看人人鲁狠狠高清| 欧美成人精品在线| 一区二区三区四区精品视频| 2021天堂在线亚洲精品专区| 992tv国产人成在线观看| 国产尤物在线播放| 欧美午夜小视频| 国产精品嫩草影院视频| 亚洲,国产,日韩,综合一区| 中文字幕在线欧美| 国产精品所毛片视频| 国产高潮流白浆视频| 99久久精品久久久久久婷婷| 色妞永久免费视频| 精品1区2区3区| 日韩国产无码一区| 91亚洲精品国产自在现线| 日韩高清在线观看不卡一区二区| 国产在线视频福利资源站| 色妞www精品视频一级下载| 亚洲色图综合在线| 国产色婷婷| 亚洲色欲色欲www网| 国产乱论视频| 嫩草影院在线观看精品视频| 宅男噜噜噜66国产在线观看| 国产在线日本| www亚洲精品| 在线免费a视频| 国产在线日本| 久久伊人操| 最新日韩AV网址在线观看| 香蕉久人久人青草青草|