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

基于慣性傳感器的穿戴式步態分析系統設計與實現

2021-09-23 06:39:24李江慧連春快李玉榕
電氣技術 2021年9期
關鍵詞:信號

李江慧 連春快 李玉榕

基于慣性傳感器的穿戴式步態分析系統設計與實現

李江慧1,2連春快2李玉榕2

(1. 福州大學海洋學院,福州 350003; 2. 福州大學福建省醫療器械和醫藥技術重點實驗室,福州 350108)

步態分析是對下肢運動情況量化評估的重要手段。本文采用基于慣性傳感器的步態時空參數估計算法,根據穿戴在小腿上的慣性傳感器節點檢測步態事件點,再對穿戴在腳跟處的慣性傳感器加速度信號進行雙重積分計算步長,根據步態時相校正加速度信號以提高步長估計精度。本文設計并實現了基于慣性傳感器的步態數據節點采集,完成了硬件和軟件的設計。為了驗證系統的可靠性和準確性,對3個志愿者進行實驗。比較本文設計的步態分析系統和商用locometrix步態分析儀的結果,步長、步速和步頻的平均相對誤差分別為5.21%、4.51%和11.87%。結果表明利用本文研制的穿戴式步態分析系統能精確估計步態時空參數。

步態分析;慣性傳感器;步態時空參數; 智慧醫療

0 引言

步行是人類賴以生存的活動方式之一,人類通過深入研究步行的機制后發現,步行運動需要肌肉系統、神經系統及骨骼系統的協調配合[1],若其中任一環節出現問題,都可能導致步態異常甚至喪失步行能力。臨床上越來越多地使用步態分析對患者的下肢運動情況進行量化分析[2],并與正常人的步態數據對比來評估患者的病情,特別對于一些早期不容易被發現的下肢疾病的診斷更有意義。通過步態分析可以檢測出人眼難以察覺的微小步態異常,從而實現對疾病的盡早發現、盡早干預[3-4],避免病情的惡化。在患者進行康復訓練時,由步態分析量化出的結果可以作為醫生評估患者下肢康復情況的依據,有助于對患者制定針對性的康復訓練方案,從而實現高效康復[5]。

三維步態分析系統是目前精度最高且計算的參數最全面的系統[6],但它也存在一些問題:①系統操作復雜,非專業人員難以對其進行調試和設備維護,如實驗前對相機的校準、相機擺放位置的調試等,如果沒調試好會嚴重影響三維步態分析系統的測量精度;②三維步態分析系統對環境有很高的要求,需要有專門的場所來減少其他光源及噪聲的干擾,且場地較小,這限制了受試者的活動范圍;③設備成本高,難以大范圍推廣,患者使用三維步態分析系統做檢測的費用成本和時間成本都比較高[7-8]。因此,越來越多的研究者期望采用廉價傳感器代替三維步態分析系統,在盡量提高測量精度的同時降低設備成本,并朝可穿戴化方向發展。市場上已有多款商用化較成熟的采用微機電系統(micro-electro mechanical system, MEMS)傳感器的步態分析系統產品,比較著名的有荷蘭Xsens Technologies公司推出的Xsens MVN慣性動作捕捉系統[9]、北京諾亦騰科技有限公司生產的鑒步運動功能評估系統等。現有的基于慣性傳感器的穿戴式系統多為全身運動捕捉系統,功能復雜,采集節點較多,穿戴過程較為 繁瑣。

針對以上不足,本文采用價廉的傳感器代替三維步態分析系統[10-11],與現有其他類型系統相比,本文系統具有信息采集節點的體積較小、方便穿戴,所用的傳感器較為廉價、可降低設備成本,利于大范圍推廣,操作簡易方便,對算力要求較低等優勢,并實現可穿戴化[12]。設計步態參數計算算法時,首先根據穿戴在小腿上的慣性傳感器測得的角速度信號確定積分時間,然后對利用穿戴于腳跟處的慣性傳感器測得的角速度信號進行雙重積分得到步 長[13],并通過將站立相期間的足部加速度置零來減小積分誤差,提高步長估計精度。步頻通過一次步態實驗中的步數和時間的比值得到,步速通過一次步態實驗的位移和時間的比值得到。實驗表明,通過該步態分析系統,能得到準確的步態參數。

1 步態參數估計算法設計

1.1 步態時相

正常人的步行是一個規律性和重復性的運動過程,若以右足腳跟著地為起點,此時身體重心逐漸轉移到身體的右側,右腿進入支撐相,而左腿離開地面向前邁步,隨后左足腳跟著地,身體重心轉移到身體左側,然后右腿離開地面進入擺動相,直至右足腳跟再次著地,這些環節構成一個完整的步態周期。在一個步態周期內,通常可以根據足部與地面是否有接觸劃分為支撐相和擺動相,腳掌與地面有接觸為支撐相,反之則為擺動相,其中支撐相大約占步態周期的60%,擺動相大約占步態周期的40%[14]。

通常通過尋找特征明顯的臨界點區分步態時相,將這些臨界點稱為步態事件點。從支撐相轉變為擺動相的步態事件點為腳尖離地點,腳尖離地后腳掌將完全離開地面處于懸空狀態。從擺動相轉變為支撐相的步態事件點為腳跟著地點,從此時開始腳掌與地面有接觸。目前用于步態事件點檢測的方法有基于角速度信號、基于加速度信號、基于腳底壓力信號的步態事件點檢測方法等,其中基于腳底壓力信號的步態事件點檢測方法準確度較高。本文利用角速度信號進行步態事件點檢測,用基于腳底壓力信號檢測出的步態事件點作為對照。

1.2 步態事件點檢測

利用腳跟著地點可以劃分步態周期,定義同側下肢前后兩次腳跟著地點的間隔為一個步態周期。穿戴在小腿上的角速度傳感器測得的垂直于矢狀面的角速度分量,相比于另兩個軸的角速度分量有較大的幅值,且周期性特征明顯,常被用來判斷腳跟著地點。該角速度信號在一個步態周期中有明顯的兩個波谷和一個波峰。波峰對應的步態事件為擺動中相點,波峰前的波谷對應的步態事件為腳尖離地點,波峰后的波谷對應的步態事件為腳跟著地點,各步態事件點位置示意圖如圖1所示。基于這一特征可以對腳跟著地點和腳尖離地點進行檢測。

圖1 步態事件點位置示意圖

結合圖1可知,判斷腳跟著地點時對應的角速度信號需要滿足條件:①是擺動中相點出現后第一個小于零的極小值點;②角速度變化量為負。為了判斷腳跟著地點,需要先檢測擺動中相。擺動中相點至少要滿足以下兩個條件:①角速度應大于一個較大的角速度閾值;②角速度變化量為一個正值。

擺動中相點判斷:讀取并判斷小腿角速度信號是否大于預設閾值,若大于且當前角速度變化量為正值,則將此角速度存入變量中,當且僅當下一個采樣周期采集的角速度信號滿足擺動中相檢測條件且角速度大小大于時,用這個角速度覆蓋中的原值。不滿足上述任一條件時,退出本輪判斷,接著讀取下一個采樣周期的角速度數據,進行新一輪的判斷。重復上述步驟,直到取得滿足條件的最大的角速度信號,即為擺動中相點。本文中取值為50°/s。

腳跟著地點判斷:一個擺動中相點之后必有一個腳跟著地點,得到擺動中相點后,判斷角速度信號是否小于零,若小于零且角速度變化量為負值,將此角速度存入變量中,當且僅當下一個采樣周期采集的角速度信號滿足腳跟著地的檢測條件且角速度小于時,用這個角速度覆蓋中的原值。不滿足上述任一條件時,退出本輪判斷,接著讀取下一個采樣周期的角速度數據,進行新一輪的判斷。重復上述步驟,直到取得滿足條件的擺動中相點后的第一個小于零的極小值點,即為腳跟著地點。

由圖1可得,判斷腳尖離地時對應的角速度信號需要滿足條件:①是腳跟著地點出現后(大于一個時間閾值)的一個小于零的極小值點;②角速度變化量為負。為了檢測腳尖離地點,需要先檢測出擺動中相和腳跟著地點,因此無法檢測出第一個步態周期的腳尖離地點,故在計算時忽略第一個步態周期和最后一個步態周期。

腳尖離地點判斷:判斷得到當前周期的腳跟著地點后,記錄腳跟著地點出現的時刻,若判斷的角速度的時刻與腳跟著地點出現的時刻的差值大于一個時間閾值,則開始判斷角速度信號的幅值是否小于零,若小于零且角速度變化量為負值,將此角速度存入變量中,當且僅當下一個采樣周期采集的角速度信號滿足腳尖離地的檢測條件且角速度小于時,用這個角速度覆蓋中的原值。不滿足上述任一條件時,退出本輪判斷,接著讀取下一個采樣周期的角速度數據,進行新一輪的判斷。重復上述步驟,直到取得滿足條件的腳跟著地點后的一個小于零的極小值點,即為腳尖離地點。

1.3 姿態解算

1)坐標系轉換

2)基于四元數的姿態解算

假設空間矢量是不動的,通過旋轉坐標系得到空間矢量,可用單位四元數來描述這兩個空間坐標系間的轉換關系[16-17],即

式中:b為旋轉后的坐標系;′為四元數的共軛;r為原始的坐標系。

空間矢量和空間矢量用四元數表達為

代入式(2)得

坐標轉換矩陣中的各元素和坐標變換矩陣中的各元素聯系起來,可推導出

在實際中,傳感器坐標系是隨著人體運動而變化的,因而單位四元數中的各元素是隨時間變化的,為求解四元數建立如下的微分方程,即

式中,為傳感器坐標系上經過加速度信號和磁力信號補償的角速度矩陣。

將式(5)展開得

1.4 步態時空參數計算算法

1)步長計算

本文定義步長為在步行運動中,同側腿先后兩次腳跟著地點間的距離。本文對利用穿戴在腳踝處的慣性傳感器測得的加速度進行二重積分來計算步長[18],但腳踝處的慣性傳感器隨著腳踝運動,空間位置的變化使加速度傳感器測得的加速度的軸分量并非都是水平向前的,因此需要對慣性傳感器進行姿態解算,得到水平向前方向的加速度信號,在一個步態周期里對其進行二重積分便可得到步長。在腳跟著地后和腳尖離地前這段時間內,腳掌的位移可以忽略,在這期間腳踝處的加速度傳感器有較大誤差。為了減少步長計算誤差,本文在腳跟著地后和腳尖離地前這段時間內對加速度進行校正,將其置為0。

簡化在步行中的足部模型,腳踝在步行過程中只有矢狀面內的角度變化。慣性傳感器的穿戴及其坐標系軸向如圖2所示。利用四元數對慣性傳感器節點的姿態進行解算,求解傳感器的俯仰角,規定固定在腳踝處靜止時俯仰角為0°。則在擺動相水平加速度的分量計算式為

積分起始時刻為當前步態周期的腳尖離地點hs(-1),積分終止時刻為下一個步態周期的腳跟著地點to(),其中1, 2, 3,…,為步態周期數,則步長()的計算式為

2)步頻、步速計算

步頻為單位時間內行走的步態周期個數,為步行總的周期數與時間的比值,即

步速為單位時間內的步行長度,用平均速度來代替,為步行長度和時間的比值,即

2 數據采集節點設計

2.1 采集節點硬件設計

采集節點要求體積小,便于受試者穿戴,核心功能是能準確測量穿戴位置的角速度信號、加速度信號、磁場信號。為提高穿戴舒適性并適應不同實驗場地的需要,采集節點將數據存儲在安全數碼卡(secure digital memory card, SD Card)中,減少了數據傳輸線對受試者活動范圍的影響。采集節點的硬件主要包括主控芯片、MPU9250傳感信號采集電路、SD卡模塊、電源電路等。主控芯片選擇意法半導公司的STM32F411RET6芯片,該芯片的內核為32位的精簡指令集計算機(reduced instruction set computing, RISC),芯片的頻率可達100MHz,高主頻滿足采集節點的頻率需求,且很大地提高了對數據的計算、處理能力。核心部分包括MPU9250傳感信號采集電路和電源控制電路。

1)MPU9250外圍電路設計

本文所用的慣性傳感器為應美盛(InvenSense)的產品MPU9250,其能同步讀取慣性九軸數據,包括加速度、角速度及磁力信息,簡化了測量系統的搭建。MPU9250與主控芯片STM32F411RET6的接線及外圍電路如圖3所示。

2)電源模塊

采集節點由3.7V鋰電池供電,節點的電源開關為按鍵S1,采集節點電源控制電路如圖4所示。

圖3 MPU9250與主控芯片STM32F411RET6的接線及外圍電路

圖4 采集節點電源控制電路

電路主要器件是IRF71055雙路場效應晶體管、BAS70—05肖特基二極管和按鍵S1,V_BAT接入鋰電池,初始時場效應晶體管U7A的2端和場效應晶體管U7B的4端為低電平,即U7A和U7B處于截止狀態,VIN_BAT沒有接入V_BAT,節點處于關機狀態。長按S1時,VD4的2、3管腳導通,主控芯片通過KEY_PWR_MCU讀取到低電平,該低電平超過2s時和主控芯片連接的PWR_CTL_MCU發出高電平,則U7A導通,V_BAT、28和DGND間形成通路,U7B也導通,V_BAT向VIN_BAT供電,采集節點開機。

采集節點處于開機狀態時,再次長按S1,通過KEY_PWR_MCU讀取超過2s的低電平,松開S1后,PWR_CTL_MCU發出低電平,U7A截止。由于電路中沒有通路,U7B的4端同為低電平,U7B截止,切斷了向VIN_BAT的供電通路,采集節點關機。

鋰電池是內置于采集節點的,不便將鋰電池取出充電,所以設計基于TP4056芯片的恒定電流充電器,鋰電池充電器電路原理如圖5所示。

圖5 鋰電池充電器電路原理

外端電源通過VUSB端口向鋰電池充電,TP4056的PROG端口外接一個電阻后接到DGND,設置恒定電流的大小。BAT口和鋰電池相連進行充電,當充電完成時,/STDBY_TP4056_MCU為低電平。

2.2 采集節點軟件設計

程序首先對系統時鐘進行配置,對通用輸入輸出接口(general purpose input/output, GPIO)引腳、模-數轉換器(analog-to-digital converter, ADC)、通用同步/異步接收/發送器(general synchronous asynchronous receiver transmitter, USART)、SD卡等硬件模塊進行初始化。待系統啟動完畢后進入待機狀態,當收到來自按鍵的開始采集命令后進入一個while死循環中,每10ms讀取一次MPU9250的九軸運動數據到數據緩存區,每100ms將采集的數據保存進SD卡中。若收到停止采集命令,則退出循環,進入待機狀態。數據采集主程序流程如圖6所示。

3 基于Matlab的程序設計

穿戴在小腿和腳跟處的采集節點采集的數據以.xls的格式導出,然后在Matlab中進一步處理。根據第1節中所設計的步態參數計算算法,在Matlab環境中編程實現。對采集的數據進行滑動平均濾波后計算步態參數。依點讀取濾波處理后的小腿處的角速度數據及腳跟處的加速度數據,利用小腿處的角速度數據判斷步態事件點,對處于腳跟著地點前和腳尖離地點后的腳跟處加速度數據置零處理,利用其余數據對測量節點進行姿態解算,得到水平向前方向的加速度信號,然后對其進行二重積分便可得到步長。當后續每發生一次腳跟著地,則計數器加1,數據讀取完后,步頻由計數器當前的數值與步行時間的比值得到,步速由歷史步長的累計和與步行時間的比值得到,最后輸出計算結果,步態參數計算程序流程如圖7所示。

圖6 數據采集主程序流程

4 實驗驗證

4.1 實驗設置

為了驗證本文提出的步態時空參數的計算算法,選擇商用locometrix步態分析儀作為步態基本參數的參考標準值,locometrix步態分析儀是由法國國家自然科學研究中心(INRA)及Laval醫院康復醫學科的Bernard教授聯合研發的。locometrix步態分析儀采集腰部的三軸加速度信號,其數據采集頻率為100Hz,然后將所采集的數據導入自帶的上位機軟件進行分析,得到在一次步態實驗中的平均速度、平均步長、步頻等參數。

圖7 步態參數計算程序流程

圖8 傳感器穿戴位置

4.2 步態參數驗證

1)評估指標

利用固定在腳后跟處的傳感器節點計算出的各步態基本參數與locometrix步態分析儀導出的參考標準值,計算相對誤差(relative error, RE)為

式中:為計算值;r為參考標準值。

2)步態實驗結果

三名受試者的步態參數計算值、標準值及RE見表1。

健康受試者1的位移測試結果如圖9所示,沒有對加速度進行校正時,位移曲線呈線性增長,表明采取在站立相對加速度進行校正這一措施,對獲得準確的位移結果起到了很明顯的作用。

利用信號采集系統采集的數據經離線計算后,平均步長的計算值和參考標準值的平均相對誤差為5.21%,平均速度的計算值和參考標準值的平均相對誤差為4.51%,步頻的計算值和參考標準值的平均相對誤差為11.87%,以上數據表明,利用信號采集系統的數據計算出的步態基本時空參數是有效、可行的。

表1 三名受試者的步態參數計算值、標準值及RE

圖9 受試者1的位移測試結果

5 結論

本文實現了步態分析系統采集節點的硬件設計,以及步態參數估計算法的設計。其中,采集子節點的硬件是以STM32F411RET6為核心芯片,同時由傳感器電路、電源模塊、存儲模塊等構成;軟件包括下位機的數據采集、數據存儲等的程序設計。設計步態時空參數計算算法時,先根據穿戴在小腿上的慣性傳感器的垂直于矢狀面的角速度信號檢測步態事件點,對利用穿戴在腳跟處的慣性傳感器獲得的加速度信號進行雙重積分計算步長,根據步態時相進行校正提高步長的估計精度。經實驗驗證,本文設計的步態分析系統能夠準確估計步態參數。

[1] 錢競光, 宋雅偉, 葉強, 等. 步行動作的生物力學原理及其步態分析[J]. 南京體育學院學報(自然科學版), 2006(4): 1-7, 39.

[2] NONNEKES J, GOSELINK R J M, RUZICKA E, et al. Neurological disorders of gait, balance and posture: a sign-based approach (review)[J]. Nature Reviews Neurology, 2018, 14(3): 183-189.

[3] ZAGO M, SFORZA C, BONARDI D R, et al. Gait analysis in patients with chronic obstructive pulmonary disease: a systematic review[J]. Gait & Posture, 2018, 61: 408-415.

[4] KIRSCHBERG J, GORALSKI S, LAYHER F, et al. Normalized gait analysis parameters are closely related to patient-reported outcome measures after total knee arthroplasty[J]. Archives of Orthopaedic and Trauma Surgery, 2018, 138(5): 711-717.

[5] KIM H, MILLER L M, FEDULOW I, et al. Kinematic data analysis for post-stroke patients following bilateral versus unilateral rehabilitation with an upper limb wearable robotic system[J]. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2012, 21(2): 153-164.

[6] DONISI L, PAGANO G, CESARELLI G, et al. Bench marking between two wearable inertial systems for gait analysis based on a different sensor placement using several statistical approaches[J]. Measurement, 2021, 173: 108642.

[7] CHO Y S, JANG S H, CHO J S, et al. Evaluation of validity and reliability of inertial measurement unit- based gait analysis systems[J]. Annals of Rehabi- litation Medicine, 2018, 42(6): 872-883.

[8] ZHENG H, BLACK N D, HARRIS N D. Position- sensing technologies for movement analysis in stroke rehabilitation[J]. Medical & Biological Engineering & Computing, 2005, 43(4): 413-420.

[9] 荷蘭Xsens Technologies公司技術網站. XSENS Products技術資料[EB/OL]. [201912]. http://www.xsens.com.

[10] BOUDARHAM J, ROCHE N, PRADON D, et al. Variations in kinematics during clinical gait analysis in stroke patients[J]. PloS one, 2013, 8(6): e66421.

[11] KARASAWA Y, TERUYAMA Y, WATANABE T. A trial of making reference gait data for simple gait evaluation system with wireless inertial sensors[C]// 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 2013: 3427-3430.

[12] WU Yinfeng, SU Yiwen, FENG Renjian, et al. Wearable- sensor-based pre-impact fall detection system with a hierarchical classifier[J]. Measurement, 2019(140): 283-292.

[13] BAMBERG S J M, BENBASAT A Y, SCARBOROUGH D M, et al. Gait analysis using a shoe-integrated wire- less sensor system[J]. IEEE Transactions on Infor- mation Technology in Biomedicine, 2008, 12(4): 413-423.

[14] 盧軍. 助力外骨骼機器人隨動控制算法設計與實現[D]. 成都: 電子科技大學, 2016.

[15] 歐陽高詢. 雙目結構光深度獲取研究及光學平臺構建[D]. 西安: 西安電子科技大學, 2015.

[16] 王涵. 無人機航跡規劃及導航定位系統研究[D]. 杭州: 浙江大學, 2017.

[17] 李洪鳳, 林康, 李斌. 基于四元數的永磁動量球位置/電流雙閉環控制[J]. 電工技術學報, 2019, 34(2): 484-492.

[18] 安耕, 楊明靜. 基于壓力和慣性傳感器的步態分析驗證研究[J]. 電氣技術, 2020, 21(6): 45-49.

Design and implementation of wearable gait analysis system based on inertial sensor

LI Jianghui1,2LIAN Chunkuai2LI Yurong2

(1. Ocean School of Fuzhou University, Fuzhou 350003; 2. Fujian Key Lab of Medical Institute and Pharmaceutical Technology, Fuzhou University, Fuzhou 350108)

Gait analysis is an important method for the quantitative evaluation of lower limb movement. This paper uses the gait space-time parameter estimation algorithm based on inertial sensor. Firstly, the gait event point is detected according to the inertial sensor worn on the calf. Then the step length is calculated by the double integration of the acceleration signal of the inertial sensor worn on the heel, and the acceleration signal is calibrated according to the gait phase to improve the estimation accuracy of step length. The gait acquisition node based on inertial signals is designed and implemented in this paper. The hardware and software design are completed. To verify the reliability and accuracy of the system, three volunteers are recruited for the experiment. Compared with the results of the locometrix device, which is a kind of commercial gait analyzer, the average relative error of step length, speed and step frequency are 5.21%, 4.51% and 11.87% respectively. The results show that the developed gait analysis system can accurately estimate the gait space-time parameters.

gait analysis; inertial sensor; gait space-time parameters; smart medical care

2021-01-25

2021-03-10

李江慧(2000—),女,本科,主要研究方向為智能康復設備研發。

國家自然科學基金(61773124)

福建省自然科學基金(2019J01544)

猜你喜歡
信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
7個信號,警惕寶寶要感冒
媽媽寶寶(2019年10期)2019-10-26 02:45:34
孩子停止長個的信號
《鐵道通信信號》訂閱單
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
基于Arduino的聯鎖信號控制接口研究
《鐵道通信信號》訂閱單
基于LabVIEW的力加載信號采集與PID控制
Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: 亚洲精品黄| 国产探花在线视频| 在线无码av一区二区三区| 欧美日韩一区二区在线播放 | 国产美女精品人人做人人爽| 狠狠操夜夜爽| 亚洲日本中文字幕天堂网| 色婷婷狠狠干| 蜜桃视频一区二区三区| 日韩一区二区三免费高清| 国产一区二区三区夜色| 免费人成视频在线观看网站| 免费看一级毛片波多结衣| 国产精品13页| 欧美成人免费一区在线播放| 国产精品香蕉在线| 日韩天堂网| 亚洲国产精品人久久电影| 伊人久久青草青青综合| 最新国产你懂的在线网址| 国产91九色在线播放| 永久天堂网Av| 91亚瑟视频| 国产在线91在线电影| 伊人精品视频免费在线| 在线播放精品一区二区啪视频| 久久精品国产免费观看频道| 色偷偷一区二区三区| 91精品视频在线播放| 91午夜福利在线观看精品| 亚洲精品第一页不卡| 青青草原国产| 91毛片网| 伊人成人在线| 伊人91视频| 色婷婷狠狠干| 免费午夜无码18禁无码影院| 久久人妻xunleige无码| 欧美精品1区2区| 99视频在线免费| 亚洲综合婷婷激情| 99尹人香蕉国产免费天天拍| 国产精品第一区在线观看| 色香蕉影院| 99re这里只有国产中文精品国产精品 | 精品视频第一页| 日韩黄色在线| 在线国产欧美| 伊人久久综在合线亚洲2019| 久久中文字幕av不卡一区二区| 免费无码又爽又黄又刺激网站 | 亚洲一级色| 波多野结衣国产精品| 影音先锋丝袜制服| www.av男人.com| 国产熟女一级毛片| 国产高颜值露脸在线观看| 91蝌蚪视频在线观看| 午夜视频www| 国产成人喷潮在线观看| 国产欧美精品一区aⅴ影院| 爆乳熟妇一区二区三区| 亚洲精品自在线拍| 青青草91视频| 国产区人妖精品人妖精品视频| 久久国产精品影院| 色婷婷电影网| 国产日韩欧美视频| 国产欧美日韩另类| 色综合天天视频在线观看| 国产女主播一区| 精品无码专区亚洲| 日本国产在线| 欧美日韩国产一级| 这里只有精品在线| 中文国产成人久久精品小说| 国产成人精品免费视频大全五级| 国内精品视频区在线2021| 免费jjzz在在线播放国产| 亚洲一级毛片| 久久久受www免费人成| 婷婷亚洲综合五月天在线|