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

五軸重型特種車輛的狀態(tài)估計及驗證

2023-03-20 02:21:28于傳強舒洪斌劉志浩陳漸偉
振動與沖擊 2023年5期
關鍵詞:模型

于傳強, 舒洪斌, 劉志浩, 陳漸偉, 高 楊

(1. 火箭軍工程大學 導彈工程學院, 西安 710025; 2. 火箭軍裝備部駐長治地區(qū)代表室, 山西 長治 046000)

隨著軍事化需求的推進,特種車輛行駛環(huán)境復雜多變,加上特種車輛載重大、質(zhì)心高等特點,在機動過程中易造成車輛操縱失穩(wěn)等事故。為了更好的判斷車輛行駛過程中的狀態(tài)變化,減小事故發(fā)生的可能性,借助有限的傳感器,基于動力學模型和參數(shù)估計方法,來獲取盡可能符合實際的狀態(tài)參數(shù),是一種經(jīng)濟有效的方法[1-3]。

近年國內(nèi)外的研究學者通過無跡卡爾曼濾波(unscented Kalman filter, UKF)及其改進算法對車輛的運動狀態(tài)參數(shù)進行研究成為一主要的估計方法,并取得了一些成果。Antonov等[4]基于建立的動力學模型采用UKF估計車輛的輪胎滑移率和質(zhì)心側(cè)偏角,在進行試驗后驗證了該估計算法具有較高的精度和較好的收斂性。趙萬忠等[5]利用UKF算法對車輛狀態(tài)進行了實時估計,并將Carsim和Simulink 聯(lián)合仿真進行驗證,結(jié)果表明,采用UKF算法響應快,估計精度要高于擴展卡爾曼濾波(extended kdman filter, EKF)。Doumiati等[6]對基于EKF和UKF的車輛質(zhì)心側(cè)偏角估計器進行了性能對比,結(jié)果證明無跡卡爾曼濾波在估計精度上具有一定優(yōu)勢。Gu等[7]以兩軸分布式驅(qū)動車輛為研究對象,建立非線性 7 自由度車輛動力學模型,采用UKF算法對縱向速度、側(cè)向速度和橫擺角速度進行了聯(lián)合估計,試驗結(jié)果表明該算法估計精度更高。上文所述皆為四輪小型車的狀態(tài)估計,且試驗結(jié)果表明無跡卡爾曼濾波的估計精度更高,針對多軸重型車輛運動狀態(tài)估計的研究則較少。

因此,本文以五軸重型車輛為研究對象,設計了基于無跡卡爾曼濾波(UKF)的狀態(tài)估計器,并以經(jīng)過實車驗證的13自由度仿真平臺開展狀態(tài)估計研究,來驗證估計算法的有效性和估計精度。

1 車輛動力學模型搭建及驗證

1.1 整車模型搭建

分析五軸車輛的動力學特性,基于車輛動力學理論,建立較為完善的整車動力學模型,為算法驗證提供高保真的仿真平臺,以進行不同工況下的仿真驗證。

建立包括縱向、側(cè)向、橫擺和每個車輪的轉(zhuǎn)動,總共十三自由度的車輛模型。

1.1.1 車體動力學模型

如圖1所示為車輛在xoy平面的運動,忽略了車輛的俯仰運動和側(cè)傾運動,且不考慮車輛的滾動阻力、迎風阻力和空氣阻力等,得出五軸特種車輛在縱向、側(cè)向以及橫擺運動的運動平衡方程如下[8]

圖1 車輛整體模型Fig.1 Overall vehicle model

縱向運動微分方程

(1)

側(cè)向運動微分方程

(2)

橫擺運動微分方程

l2(Fy2l+Fy2r)-l3(Fy3l+Fy3r)-

l4(Fy4l+Fy4r)-l5(Fy5l+Fy5r)

(3)

式中:m為整車質(zhì)量;g為重力加速度;li為各軸到質(zhì)心的距離;B為車輛輪距;vx為縱向速度;vy為側(cè)向速度;Iz為車輛繞z軸的轉(zhuǎn)動慣量;wz為車身橫擺角速度;Fxil和Fxir為第i軸左、右側(cè)輪胎的地面縱向力,F(xiàn)yil和Fyir為第i軸左、右側(cè)輪胎的地面?zhèn)认蛄?。在輪胎坐標系上,第i軸左、右側(cè)縱向力表示為Fxwil和Fywir,側(cè)向力表示為Fywil和Fywir。車輛坐標系與輪胎坐標系之間的輪胎力的關系式為

(4)

式中,α為車輪側(cè)偏角。

1.1.2 車輪運動模型

考慮車輪在平面內(nèi)的旋轉(zhuǎn)自由度,建立車輪的運動模型如圖2所示。不考慮地面?zhèn)认蛄Φ挠绊?,在車輪旋轉(zhuǎn)運動的方向上,車輪受到驅(qū)動力矩Ti、制動力矩Tbi、反力矩Tdi[9];由于輪胎的彈性和地面的摩擦作用,垂向力的作用點偏前于輪胎接地印跡幾何中心一個距離還對輪心形成一個滾動阻力矩

圖2 車輪受力模型Fig.2 Wheel force model

Tfi=FziΔi=Fzifrrw

(5)

式中:Fzi為地面對輪胎的作用力;fr為滾動阻力系數(shù);rw為輪胎滾動半徑。

將作用于車輪上的力和力矩對輪心取力矩,得到車輪滾動方程

(6)

式中:Iw為車輪轉(zhuǎn)動慣量;ωi為車輪旋轉(zhuǎn)角速度。

1.1.3 輪胎模型

輪胎作為車輛與地面的一個交互部件,它的建模精度直接影響著整車模型仿真結(jié)果,魔術(shù)公式(magic formula)輪胎模型能夠較為精確的表征輪胎的非線性,采用魔術(shù)公式輪胎模型來描述輪胎的動力學特性。它的表達公式為

(7)

式中:y為輪胎縱向力或側(cè)向力;x為輪胎側(cè)偏角或滑移率;Sv為垂直偏移量;Sh為水平偏移量;B、C、D、E分別為剛度因子、形狀因子、峰值因子、曲率因子,其數(shù)值可通過參數(shù)辨識得到。由于該類車型所使用的重載子午輪胎具有垂向載荷范圍大、扁平比大的特點,因此本文采用文獻[10]中針對重載工況下的魔術(shù)公式修正模型。由式(7)可知,輪胎力的計算需要側(cè)偏角和滑移率的輸入,接下來計算這兩個值。

因為本文中的建模忽略側(cè)傾和俯仰運動對速度產(chǎn)生的影響, 由于車輛橫擺運動的緣故,各輪輪心的速度可由車輛質(zhì)心速度和橫擺角速度表示,則各個輪心的縱向速度和側(cè)向速度為

(8)

vyi=vy±liwz

(9)

輪心沿縱軸方向的速度為

vi=vxicosδi+vyisinδi

(10)

得各個輪胎的側(cè)偏角αi和滑移率λi為

(11)

(12)

式中:δi為各軸車輪轉(zhuǎn)角;ωi為各車輪旋轉(zhuǎn)角速度。

1.1.4 轉(zhuǎn)向系統(tǒng)運動學模型

建立如圖3所示的阿克曼轉(zhuǎn)向模型,在阿克曼轉(zhuǎn)向模型中,為了使所有車輪處于純滾動或只有極小滑動的狀態(tài),要求所有車輪具有同一個瞬時轉(zhuǎn)向中心,該五軸車的1、2、4、5軸為轉(zhuǎn)向軸,瞬時轉(zhuǎn)向中心位于三軸的延長線上。為簡化轉(zhuǎn)向模型,假定車輛轉(zhuǎn)向系統(tǒng)的傳動比不變?yōu)閕sw,傳動比為方向盤轉(zhuǎn)角和一軸左側(cè)輪胎轉(zhuǎn)角的比值,方向盤轉(zhuǎn)角為δsw,根據(jù)圖3得出各個車輪轉(zhuǎn)角的幾何關系:

圖3 阿克曼轉(zhuǎn)向模型Fig.3 Ackerman steering model

(13)

1.1.5 整車模型集成

通過Matlab/Simulink實現(xiàn)的整車集成模型,輸入方向盤信號給轉(zhuǎn)向模塊,轉(zhuǎn)向模塊通過Matlab Function函數(shù)表示出前兩軸和后兩軸的8個轉(zhuǎn)向輪的轉(zhuǎn)角并輸出給車輛運動模塊。車輛運動模塊的輸入包括 8 輪轉(zhuǎn)角和 10輪轉(zhuǎn)矩,輸出為10輪轉(zhuǎn)速和縱向速度、縱向加速度、側(cè)向速度、側(cè)向加速度和橫擺角速度等,若在后續(xù)研究中需要其他參數(shù),可在Simulink中增加相應的輸出。

1.2 實車試驗驗證

為了驗證所搭建模型的準確性,利用所具有的特種車輛行駛狀態(tài)監(jiān)測系統(tǒng)進行實車驗證。實車試驗方案為:方向盤轉(zhuǎn)角傳感器采集轉(zhuǎn)向信息、SPEEDBOX-INS車輛姿態(tài)測量系統(tǒng)、衛(wèi)星導航傳感器等實時獲得車輛運行過程中車速、軌跡、高程、車輛姿態(tài)、橫擺角速度等信息,通過CAN總線輸出至擁有多通道、多參數(shù)的DEWE43數(shù)據(jù)采集系統(tǒng),用Race Technology和DEWEsoftX3軟件對試驗數(shù)據(jù)進行記錄分析處理,具體試驗過程如圖4所示。

圖4 實車試驗Fig.4 Real vehicle test

車輛模型驗證的過程如圖5所示。將試驗測得的方向盤轉(zhuǎn)角δexp、縱向速度Vx,exp作為整車模型的輸入,其中速度控制模型根據(jù)作為期望速度的縱向速度Vx,exp和當前模型實際車速Vx,mod,將對應的驅(qū)動力或制動力傳給整車模型的車輪運動模塊。將仿真模型得到的縱向速度Vx,mod、橫擺角速度wz,mod、側(cè)向加速度ay,mod與試驗數(shù)據(jù)進行比較。因為搭建的模型目的是作為狀態(tài)觀測的仿真平臺,所需觀測的縱向速度是各類控制算法中重要的輸入?yún)?shù),且橫擺角速度和質(zhì)心側(cè)偏角狀態(tài)量是以橫向動力學模型為基礎,橫擺角速度與側(cè)向加速度二者有研究用來驗證模型的準確度[11]。根據(jù)試驗數(shù)據(jù)與仿真數(shù)據(jù)的差值,通過試錯調(diào)整模型參數(shù),最終確定最小均方根誤差下的模型。

圖5 驗證模型Fig.5 Validate the model

基于現(xiàn)有的試驗條件,安裝好試驗器材后在某段高原公路與砂石路上進行實車試驗,選取的驗證路況分別為近乎直線行駛路段和帶有兩處較大拐彎路段,將選取路段的數(shù)據(jù)信息用于車輛模型參數(shù)調(diào)整與驗證。

路況1該路段為近乎直道行駛路段,主要為砂石路,較為顛簸。實車行駛軌跡、方向盤轉(zhuǎn)角分別為圖6、圖7;模型與實際值的縱向速度、橫擺角速度和側(cè)向加速度對比圖分別為圖8(a)~(c)。

圖6 實車行駛軌跡圖Fig.6 Real vehicle driving trajectory map

圖7 方向盤轉(zhuǎn)角Fig.7 Steering wheel angle

(a) 縱向速度對比

(b) 橫擺角速度對比

(c) 側(cè)向加速度對比圖8 模型驗證對比Fig.8 Model validation comparison

路況2該路段為帶有兩處較大彎道行駛路段,主要為柏油路與砂石路。實車行駛軌跡、方向盤轉(zhuǎn)角分別為圖9、圖10;模型與實際值的縱向速度、橫擺角速度和側(cè)向加速度對比圖分別為圖11(a)~(c)。

圖9 實車行駛軌跡圖Fig.9 Real vehicle driving trajectory map

圖10 方向盤轉(zhuǎn)角Fig.10 Steering wheel angle

(a) 縱向速度對比

(b) 橫擺角速度對比

(c) 側(cè)向加速度對比圖11 模型驗證對比Fig.11 Model validation comparison

將實車所得的橫擺角速度和側(cè)向加速度去噪處理,實車數(shù)據(jù)與模型輸出數(shù)據(jù)的均方根誤差如表1所示,由圖8、圖11以及表1可知,實際值與模型輸出的縱向速度、橫擺角速度以及側(cè)向加速度對比,二者會有一些誤差,但總體上變化趨勢與實際相符,較為吻合。其原因是數(shù)學模型中忽略了懸架的動態(tài)載荷且模型中輪胎的轉(zhuǎn)角和實際轉(zhuǎn)角會有一些誤差,進而導致簧上質(zhì)量的載荷轉(zhuǎn)移有些誤差,且實車試驗環(huán)境為開闊的高原公路與砂石路,起關鍵因素的輪胎會與平原的胎壓有所不同,魔術(shù)公式輪胎模型的參數(shù)擬合會有所差異,間接影響輪胎與地面相應的作用力,同樣會受風速尤其是側(cè)向風速、路面不平整度等一些因素的影響。該實車試驗驗證結(jié)果能夠說明該模型可用于狀態(tài)估計器的仿真驗證。最終調(diào)整驗證模型過程中,在原有車輛數(shù)據(jù)的基礎上,主要對轉(zhuǎn)向系傳動比、車輛質(zhì)量繞z軸轉(zhuǎn)動慣量這幾個參數(shù)做微調(diào)。在調(diào)整驗證過程中,轉(zhuǎn)向系傳動比的變化對橫擺角速度與側(cè)向加速度的變化影響較大,實車轉(zhuǎn)向系傳動比隨著檔位的變化而變化,模型中將其作為一個定值,所以傳動比是調(diào)整變化較大的一個參數(shù),按照均方根誤差數(shù)值最小的調(diào)整基準,模型最終確定的主要參數(shù)數(shù)值如表2所示。

表1 均方根誤差Tab.1 Root mean square error

表2 車輛參數(shù)Tab.2 Vehicle parameters

2 UKF觀測器設計

2.1 五軸車觀測方程

所建立的車輛模型自由度多,且會因模型的復雜度等影響估計的實時性,因此在建立狀態(tài)方程時將13自由度的模型降為三自由度車輛模型,作為多軸車輛狀態(tài)估計模型。

為研究車輛行駛的操縱穩(wěn)定性,3自由度車輛模型考慮了縱向、側(cè)向、橫擺方向的運動,車輛模型的動力學平衡方程為[12]

縱向運動微分方程

(14)

側(cè)向運動微分方程

(15)

橫擺運動微分方程

(16)

非線性系統(tǒng)的狀態(tài)空間方程可以表示為

(17)

式中:xk+1,zk,uk分別為系統(tǒng)狀態(tài)量、觀測量和輸入量;wk和vk分別為過程噪聲和觀測噪聲。

根據(jù)所建立的車輛動力學模型,方向盤的轉(zhuǎn)角及縱向加速度可以通過傳感器精確測量,二者作為輸入變量為

u=[δax]T

(18)

式中:δ為方向盤轉(zhuǎn)角;ax為縱向加速度。

三維狀態(tài)變量為

x(t)=[wzβvx]T

(19)

式中:wz為橫擺角速度;β為質(zhì)心側(cè)偏角;vx為縱向速度。

橫向加速度ay為觀測量[13]。

2.2 系統(tǒng)估計方程

利用UKF濾波算法進行狀態(tài)估計,建立系統(tǒng)估計方程。

UKF利用相似分布原理,計算的Sigma點集跟原分布均值和協(xié)方差一樣,帶進非線性系統(tǒng)后進行無跡變換,進而得出估計量:

步驟1計算Sigma采樣點χ(i)并求出采樣點的權(quán)值ω(i)

(20)

(21)

步驟2用式(20)重新求得Sigma點集和其相應權(quán)值并計算Sigma點的進一步預測。

χ(i)k=

(22)

(23)

進而求出狀態(tài)量預測和協(xié)方差矩陣的預測

(24)

(25)

步驟3將步驟2計算出的狀態(tài)預測量通過無跡變換計算出預測量的采集樣點

(26)

計算Sigma點的預測觀測量

(27)

求觀測的預測均值,得出預測的均值及協(xié)方差

(28)

(29)

(30)

步驟4計算卡爾曼增益矩陣

(31)

步驟5更新狀態(tài)量和協(xié)方差矩陣[14]

(32)

(33)

根據(jù)以上公式可得出三自由度整車狀態(tài)方程及測量方程

(34)

(35)

式中:ki為側(cè)偏剛度;k1i為各軸輪胎轉(zhuǎn)角與一軸的比例值。

3 觀測器虛擬驗證

在Simulink中將建立的十三自由度整車動力學模型與UKF整合,如圖12所示,利用經(jīng)過試驗驗證過的高保真模型進行仿真驗證。車輛先加速到20 m/s,在到達20 m/s后分別對方向盤進行角階躍輸入、角脈沖輸入和正弦輸入驗證。

圖12 估計算法仿真驗證Fig.12 Simulation verification of estimation algorithm

工況1在角階躍輸入工況中,車輛勻加速到20 m/s后,在35 s時對方向盤轉(zhuǎn)角進行角階躍輸入如圖13(a)所示,縱向速度、橫擺角速度和質(zhì)心側(cè)偏角的估計值和實際值如圖13(b)~(d)所示。

(a) 方向盤轉(zhuǎn)角

(b) 縱向速度

(c) 橫擺角速度

(d) 質(zhì)心側(cè)偏角圖13 角階躍輸入的UKF估計結(jié)果Fig.13 UKF estimation result of angular step input

工況2在角脈沖輸入工況中,車輛勻加速到20 m/s后,在35 s時對方向盤轉(zhuǎn)角進行角脈沖輸入如圖14(a)所示,縱向速度、橫擺角速度和質(zhì)心側(cè)偏角的估計值和實際值如圖14(b)~(d)所示。

(a) 方向盤轉(zhuǎn)角

(b) 縱向速度

(c) 橫擺角速度

(d) 質(zhì)心側(cè)偏角圖14 角脈沖輸入的UKF估計結(jié)果Fig.14 UKF estimation result of angular pulse input

工況3在正弦輸入工況中,車輛勻加速到20 m/s后,在35 s時對方向盤轉(zhuǎn)角進行正弦輸入如圖15(a)所示,縱向速度、橫擺角速度和質(zhì)心側(cè)偏角的估計值和實際值如圖15(b)~(d)所示。

(a) 方向盤轉(zhuǎn)角

(c) 橫擺角速度

(d) 質(zhì)心側(cè)偏角圖15 正弦輸入的UKF估計結(jié)果Fig.15 UKF estimation result of sinusoidal input

三種工況下的均方根誤差如表3所示。由圖13、圖14、圖15和表3可知,在車輛直線勻加速行駛階段,三種工況下的縱向速度幾乎是吻合的,經(jīng)計算此階段的縱向速度均方根誤差僅為0.008 4 m/s,橫擺角速度和質(zhì)心側(cè)偏角的估計會隨著速度的增大而有些偏差,均方根誤差為0.001 6 rad/s和0.000 17 rad。

表3 均方根誤差Tab.3 Root mean square error

在速度穩(wěn)定后的勻速行駛階段,方向盤轉(zhuǎn)角發(fā)生突變時,縱向速度和質(zhì)心側(cè)偏角能夠保持較好的跟隨性,且三種工況下的質(zhì)心側(cè)偏角均方根誤差均值為0.003 8 rad。在這個階段,三種工況下的橫擺角速度誤差在角脈沖輸入中出現(xiàn)最大差值為0.048 rad/s,但很快會趨于穩(wěn)定并誤差縮小。在正弦輸入工況下,汽車行駛狀態(tài)時刻發(fā)生變化,橫擺角速度的估計結(jié)果能夠?qū)嶋H值保持良好的跟隨性。橫擺角速度的誤差主要在峰值和峰谷處,其原因在于UKF估計器中初始協(xié)方差數(shù)值的設置,協(xié)方差大小的取值影響著橫擺角速度和質(zhì)心側(cè)偏角的估計精度,考慮到質(zhì)心側(cè)偏角在實際測量中困難更大,設備價格更高的原因,在本文中的估計器中所設置的大小更符合質(zhì)心側(cè)偏角的估計,對于橫擺角速度相對而言差一些,這需要在后期做進一步的參數(shù)辨識以及測量噪聲和過程噪聲自適應的改進。

本文還對30 s勻加速過程中進行正弦輸入,即加速階段的轉(zhuǎn)向驗證,縱向速度、橫擺角速度和質(zhì)心側(cè)偏角的均方根誤差分別為0.027 9 m/s、0.008 3 rad/s和0.004 9 rad。

綜上結(jié)果表明該狀態(tài)觀測器精度較高,可實現(xiàn)對縱向速度、橫擺角速度和質(zhì)心側(cè)偏角參數(shù)的動態(tài)估計。

4 結(jié) 論

(1) 針對五軸重型車輛搭建整車13自由度模型,依據(jù)現(xiàn)有的試驗條件將試驗測得的參量作為模型的輸入,將可以反應車輛狀態(tài)的關鍵參數(shù)-縱向速度、橫擺角速度與側(cè)向加速度作為驗證對比的對象,在綜合分析后到試驗結(jié)果表明所搭建模型能夠作為算法驗證的仿真平臺。

(2) 面向五軸重型車輛設計了基于無跡卡爾曼濾波的狀態(tài)估計器,通過13自由度模型分別在方向盤角階躍輸入、角脈沖輸入和正弦輸入工況下進行驗證。驗證結(jié)果表明,五軸重型車輛的行駛狀態(tài)發(fā)生改變時,所設計的估計器能夠較精確的估計出狀態(tài)參數(shù),具有良好的跟隨性,算法精度高,可為后續(xù)的車輛穩(wěn)定性控制提供可靠的參數(shù)。

(3) 本文所搭建的未考慮車輛應用的實際場景,會導致仿真結(jié)果與試驗結(jié)果有偏差,后續(xù)要結(jié)合五軸重型車輛的作業(yè)環(huán)境搭建仿真精度更高的整車模型。且在實車試驗驗證中,因試驗條件以及試驗環(huán)境的限制,現(xiàn)只能對車輛的橫擺角速度與側(cè)向加速度進行驗證,在今后的試驗中,需著重考慮質(zhì)心側(cè)偏角、輪胎力等參數(shù)的對比驗證。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 高清不卡一区二区三区香蕉| 91在线精品免费免费播放| 自拍亚洲欧美精品| 亚洲天堂网在线观看视频| 国产亚洲精品自在线| 国产精品毛片一区视频播| 日韩成人在线网站| 亚洲国产欧美自拍| 国产三区二区| 香蕉99国内自产自拍视频| 天堂久久久久久中文字幕| 成年女人a毛片免费视频| 欧美亚洲国产视频| 亚洲成aⅴ人片在线影院八| 性网站在线观看| 韩日免费小视频| 国产后式a一视频| 中文无码日韩精品| 色综合五月婷婷| 欧洲成人在线观看| A级毛片无码久久精品免费| 伊人婷婷色香五月综合缴缴情| A级毛片高清免费视频就| 国产精品午夜福利麻豆| 综合色区亚洲熟妇在线| 被公侵犯人妻少妇一区二区三区| 天堂av综合网| 最新国产午夜精品视频成人| 中文字幕波多野不卡一区| 一级在线毛片| 四虎精品黑人视频| 亚洲色精品国产一区二区三区| 亚洲黄色网站视频| jizz国产在线| 亚洲欧美成人在线视频| 色爽网免费视频| 日本AⅤ精品一区二区三区日| a天堂视频| 午夜国产理论| 国产欧美视频在线观看| 亚洲国产欧美国产综合久久| 国产在线视频欧美亚综合| 在线观看无码a∨| 园内精品自拍视频在线播放| 日韩专区第一页| 国产成人一区在线播放| 亚洲国产成人久久77| 日本午夜视频在线观看| 久草性视频| 久久国产精品电影| 91无码国产视频| 国产xx在线观看| 91久久夜色精品| 国产剧情一区二区| 久草热视频在线| 亚洲精品国产综合99久久夜夜嗨| 有专无码视频| 久久这里只有精品23| 亚洲欧美色中文字幕| 五月天综合网亚洲综合天堂网| 午夜国产精品视频| 在线国产91| 亚洲三级影院| 亚洲AV无码不卡无码| 亚洲精品男人天堂| 999国内精品久久免费视频| 九九热在线视频| 日本一区二区三区精品国产| 免费高清a毛片| 人妻丰满熟妇αv无码| 欧美天堂在线| 国产福利免费观看| 国产精品极品美女自在线| 99偷拍视频精品一区二区| 国产美女免费| 欧美成人精品在线| 91日本在线观看亚洲精品| 91精品啪在线观看国产60岁 | 日韩免费毛片| 9丨情侣偷在线精品国产| 欧美中文字幕在线二区| 九九这里只有精品视频|