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

基于時序的水田平地機俯仰角預測建模與試驗

2018-06-21 09:29:26趙潤茂羅錫文唐靈茂
農業工程學報 2018年11期
關鍵詞:模型

趙潤茂,胡 煉,2※,羅錫文,2,唐靈茂,周 浩,杜 攀,賀 靜

(1. 華南農業大學南方農業機械與裝備關鍵技術教育部重點實驗室,廣州 510642;2. 南方糧油作物協同創新中心,長沙 410128)

0 引 言

農田高低不平引起農業機械行駛姿態變化,導致機具作業質量與效率下降。為實現機具位姿精準控制,水田激光平地機采用IMU(inertial measurement unit)實時檢測拖拉機橫滾角偏差并自動調節平地鏟水平傾角,以使平地鏟不受車身姿態影響而始終保持水平[1-4];馬敏等[5]采用紅外光電傳感器對煙草高度精確仿形,實現煙草打頂高度的精確控制以及打頂后抑芽劑的精量噴施;李君等[6]采用超聲波冠形探測方法,設計了一種懸掛式電動柔性疏花機,其伺服控制系統能滿足疏花機平面位置的精度要求;王松林[7]設計了檢測噴桿高度的接觸式傳感裝置,通過液壓控制回路實現噴霧高度的自動調節。上述所有自動仿形系統均采用傳感器測量被控對象實時狀態,計算偏差、施加控制、測量反饋,形成閉環,反饋輸入信號來源于傳感器當前時刻測量值,系統偏差計算始終滯后于被控對象的一個運動周期[8],影響即時控制響應。但通過對農機運動狀態進行預測,根據預測信息設計控制律,可提高農具控制精度和農機作業效率[9]。

預測控制(model predictive control,MPC)的關鍵是預測模型,可根據對象的歷史信息和未來輸入,預測其未來狀態或輸出[10]。國內外諸多學者對此開展了大量研究。Mozaffari等[11]通過進化最小二乘算法預測混合動力轎車車速,以實現動力優化控制;Wu等[12]采用卡爾曼濾波器預測未來短時內車輛位置;Gallieri等[13]以AR(9)為預測模型設計了海藻收割機高程前饋控制系統;Koo等[14]以船體和飛機動力學方程為預測模型,實現了艦載機降落的預測控制;趙曉莉等[15]采用新陳代謝GM(1,1)模型預測玉米葉片的生長長度;張軍等[16]采用灰色理論對溫室系統溫度的不精確性、時變性和多擾動等進行預測補償及節能控制;此外,還有大量用于多旋翼無人機[17-19]和民用飛機[20]高程與姿態預測控制的預測模型研究。可見,只要是具有預測功能的信息集合,不論其具體表現形式如何均可作為預測模型,包括經典的狀態方程和傳遞函數等[21],但未見有為實現農機精準作業控制而進行的相關預測模型研究。考慮到農業機械田間運動姿態的隨機、不精確、時變和多擾動等特性,精確的農機姿態模型難以獲得,本文以水田平地機為研究對象,提出一種基于時間序列分析的、面向水田平地機平地鏟預測控制律設計的預測模型在線辨識及其參數估計方法,并通過試驗分析驗證所建模型的準確性。

1 基于時間序列的水田平地機俯仰角預測模型

將平地機田間行駛時的俯仰角變化看作一個隨機過程 Yt(t > 0),其在每一采樣時刻的樣本實現y1,y2,…,yt, …,yn構成一個俯仰角時間序列。

1.1 俯仰角時序模型結構識別

為識別描述平地機田間行駛俯仰角的時序模型是自回歸(auto-regressive,AR)、滑動平均(moving-average,MA)或自回歸滑動平均(auto-regressive and movingaverage,ARMA)模型,分別計算俯仰角數據樣本的自相關函數(autocorrelation function,ACF)和偏自相關函數(partial autocorrelation function,PACF),根據其呈現的“拖尾”或“截尾”特征進行判別[22],如表1所示。

表1 時間序列ACF和PACF的一般特征Table 1 ACF and PACF characteristics for time series

ACF和PACF的計算式[23]分別為式(1)和式(2)。式中n為樣本數據量,y為樣本均值,k為遲滯階數,ρ?k為俯仰角序列ACF的估計,φ?kk為俯仰角序列PACF的估計,

由表 1和田間預試驗可知,平地機田間行駛俯仰角時序的 ACF和 PACF均呈現“拖尾”現象,因此采用ARMA表達平地機俯仰角動態變化,模型結構為式(3)。

式中 A ( q ) = 1 + a1q-1+…+,

C( q ) = 1 + c1q-1+…+q-nc,y( t)為t時刻模型輸出,na為自回歸項階數,nc為滑動平均項系數,q為后移算子, e( t)為系統噪聲, e( t)~N ID(0, σ2)。

1.2 俯仰角ARMA模型定階

采用基于信息量準則的 AIC(akaike information criterion)定階方法[24]估計平地機俯仰角預測模型ARMA的階次(na, nc),式(4)。

式中N為樣本個數;為殘差方差。

AIC定階需要進行大范圍內的階數搜索,并對每次搜索得到的模型計算AIC信息量[25];隨著搜索階數的遞增,式(4)中第一項數值下降,第二項數值上升,直至某階AIC(na, nc)最小即為最佳階數。

1.3 模型參數在線辨識

模型結構確定后,需根據平地機俯仰角傳感信息估計式(3)中的參數 a1, a2, …,c1, c2,… cnc值。水田地況復雜多變,采用批處理最小二乘法(least squares,LS)離線辨識得到的定常參數模型無法長期描述平地機俯仰角的動態變化;而始終利用一定長度的最新數據進行在線參數估計雖可有效跟蹤系統時變,但運算速度無法滿足預測控制實時性要求,尤其是ARMA中MA部分的參數估計需要大量非線性計算[26]。為保證平地機俯仰角的時變模型參數實時有效估計又避免每次傳感數據更新均需重復估計參數的耗時運算[27],采用平地機最新俯仰角數據對上一步估計出的參數進行校正,即遞推最小二乘算法(recursive least square,RLS)[28]。RLS參數遞推方程為式(5)~(7)。

式中( t) = [ a1, … ,,c1, … cn]為t時刻模型擬辨識參數;c

y?( t)為t時刻平地機俯仰角估計值; y( t)為t時刻平地機俯仰角測量值;φ?( t ) = [- y ( t- 1 ),… ,- y ( t- na),( t - 1),… ,e?( t - nc)]T, e? ( t) = y( t) -( t?(t ),為t時刻模型殘差;K( t)為校正系數。

2 水田平地機俯仰角預測模型構建

2.1 數據采集與預處理

試驗在華南農業大學岑村實驗基地進行,將 AHRS(XSENS, Mti-300)安裝于以久保田2ZGQ-6D1為配套動力的1PJ-3.0型水田激光平地機座椅后方,采集平地機行駛俯仰角數據,采樣頻率為20 Hz。傳感數據通過串口通信RS232發至上位機實時處理,平地機俯仰角在線建模在MATLAB 2016b平臺中實現。

初始采集得到的數據由于可能包含趨勢性或周期性等因素,無法保證時間序列的統計特征不隨時間推移而變化[29],故不能直接用于時序建模。因此,采用差分法去除平地機俯仰角時序中非平穩因素,建模后通過逆差分還原序列[30]。圖1為任意采集的一段時長為200 s的平地機水田行駛俯仰角度值,由圖可見,在水田中平地機行駛姿態起伏變化頻繁,為非平穩隨機過程,通過一階差分使序列平穩化,差分處理后俯仰角數據圍繞零均值上下波動,序列基本平穩。序列平穩性采用MATLAB提供的檢驗函數adftest進行ADF單位根檢驗,函數返回值為1表明序列平穩,可以進行時序建模。

圖1 平地機水田行駛俯仰角序列Fig.1 Time series for pitch angle of leveler in paddy field

2.2 預測模型結構識別

差分后序列ACF與PACF值的計算結果如圖2所示。圖2表明,俯仰角序列的ACF與PACF值絕大部分落在2倍標準差外,其衰減過程震蕩且收斂緩慢,可確定ACF與 PACF“拖尾”,依據表 1,選擇 ARMA(na,nc)(na>0,nc>0)模型對平地機俯仰角序列建模[23]。從na= 1,nc= 1開始逐次升階擬合前20階平地機俯仰角的ARMA(na,nc)模型并計算AIC值。對于高階模型定階,選用(2n,2n–1)的建模方案[26]進行最佳階數搜索以減小運算次數,AIC計算結果如表2所示,最小AIC值為–5.272 9,此時模型階數na=18,nc=17。再在(18,17)附近各進行一步搜索(19,18)與(17,16)模型,計算結果AIC(19,18)=–5.269 5,AIC(17,16)= –5.271 9。綜上可確定最佳模型為:ARMA(18,17)。

圖2 差分序列ACF與PACFFig.2 ACF and PACF of differential time series

2.3 預測模型參數估計

采用RLS算法在線辨識模型參數,令式(5)~(7)中初值?(0)=0,P(0)=μI,其中I為35階單位矩陣,μ=1× 1 06;可以證明當μ取足夠大的正數時,經過不多次的遞推,初值影響就會消失[27]。輸出俯仰角 ARMA(18,17)模型的向前1步遞推值并與AHRS實測值比較,如圖3所示,誤差分析表明最大絕對誤差MAE(maximum absolute error)為 3.527 1°,均方根誤差 RMSE(root mean square error)為 1.274 9°。

表2 AIC計算結果Table 2 Results of AIC calculation

圖3 RLS辨識模型輸出結果Fig.3 Output results of model identified by RLS

對圖3的分析可看出,約50 s后模型絕對誤差開始增大,在近150 s處達到最大為–3.527 1°,這是因為隨著數據量的增加,P ( t)、K ( t)變小,產生“數據飽和”現象,從而對θ?(t)的修正能力變弱,導致時變參數估計失敗。為改善RLS在時變系統中的跟蹤性能,提出采用遺忘因子[28,31](FF,forgetting factor)λ(0<λ≤1)對在線傳感數據加權以強調新息對參數估計的貢獻,保證RLS可跟蹤系統參數變化。FFRLS辨識算法見式(8)~式(10)。

FFRLS算法中遺忘因子λ的取值對俯仰角模型精度有重要影響。若λ選取過小,模型過多重視新息,模型輸出較實測值跳動震蕩大;若λ選取過大,即建模更注重序列歷史長期趨勢,模型輸出則較為平滑,對平地機俯仰角值的高頻變化表達不夠。一般地,若時間序列的連續T0個采樣趨于穩定常值,則有式(11)[32]。

對于時變與非線性較強的平地機俯仰角序列,T0=3,此時0.67λ=。為更加精確確定λ值,分別在0.61λ≤≤時進行仿真試驗并以模型誤差RMSE與MAE評價λ取值優劣,結果如表3所示。依據表3可見,對于當前序列,當遺忘因子0.7λ=時,模型輸出RMSE為0.026 1°,MAE為0.027 0°,誤差最小。

2.4 預測模型適用性檢驗

在將模型用于平地機俯仰角預測前,應對其適用性進行檢驗,主要包括殘差分析與模型仿真。殘差分析檢驗殘差是否為白噪聲,若為白噪聲,說明姿態模型擬合良好,否則說明殘差序列仍有相關信息未被提取,需要重新建立模型[33];通過模型仿真可對模型輸出誤差進行直觀分析與比較。根據式(1),ARMA模型的1步遞推輸出為式(12)所示。

表3 λ取值對模型擬合殘差的影響Table 3 Influence of different λ values on residual error

寫成最小二乘形式為式(13)所示。

對ARMA(18,17)模型進行殘差自相關檢驗,如圖4所示。圖4表明高于1階滯后的殘差自相關函數值均落在非顯著相關的99%置信區間內,說明殘差序列值之間不存在任何相關性,為白噪聲,因此當前 ARMA(18,17)模型適用于平地機俯仰角序列建模。

圖4 殘差相關性檢驗Fig.4 Correlation test of residual

3 試驗與結果分析

為進一步驗證 ARMA(18,17)模型對描述平地機俯仰角的普適性與魯棒性,于2017年5月31日在華南農業大學岑村實驗基地(北緯 23.163429°,東經113.370198°)任選取包含坑洼、坡度、路障的自然水泥路面以模擬階躍、斜坡、正弦激勵信號源。駕駛 1PJ-3.0型水田激光平地以作業速度(1.2 m/s)分別完成過溝,爬坡,前進-后退越障、前進-掉頭-前進越障,每組試驗重復3次,如圖5所示。利用AHRS實時采集上述地況下平地機行駛俯仰角度值并同步在線建模,采樣率為 20 Hz,依次比較模型輸出與AHRS實測值,如圖6a~6d所示,圖6e為平地機水田(旋耕后靜置48 h)作業時連續200 s俯仰角測量值與模型輸出對比。

從圖6看出,在4種模擬典型地況及水田環境下,算法的一步遞推輸出結果與AHRS實測值趨勢變化相同,所建模型可有效跟蹤平地機俯仰角動態變化,具有較高準確性。誤差分析結果如表 4所示,各模擬地況下(圖6a~圖6d)模型輸出MAE與RMSE值均小于0.1°;水田作業時(圖6e)模型輸出MAE與RMSE值小于0.2°,誤差較大是由于復雜環境中模型時變性強造成。綜上所述,用ARMA(18,17)對平地機俯仰角建模并通過FFRLS實時估計和遞推更新模型參數的方法是可行的,所建模型輸出精度滿足平地機姿態預測控制律設計的要求。

圖5 不同地況模擬試驗Fig.5 Algorithm testing under different situations

圖6 不同地況下ARMA(18,17)預測與AHRS實測對比Fig.6 Comparison between output value of ARMA(18, 17)and AHRS measurement

表4 不同模擬地況下ARMA(18,17)輸出誤差Table 4 Output errors of ARMA (18,17) under different situations

4 結 論

1)基于時間序列分析方法,建立了水田平地機俯仰角序列的 ARMA(18,17)模型,通過模型擬合殘差單位根ADF檢驗,驗證了所建模型適用性。

2)為滿足預測控制系統對預測模型辨識的實時性要求,采用遺忘因子遞推最小二乘法對模型參數在線估計并遞推更新,通過模型仿真試驗得到其遺忘因子取值0.7。分別采集幾種模擬典型地況及實際水田作業時平地機行駛俯仰角數據在線建模,將模型輸出值與 AHRS實測值比較,結果表明:模型輸出與AHRS實測值變化趨勢一致,MAE與RMSE計算值均不超過0.2°,說明采用ARMA(18,17)對平地機俯仰角建模并通過FFRLS在線估計模型參數的方法是可行的,滿足控制系統對平地機俯仰角預測的誤差要求。

面向水田平地機姿態預測控制律設計,本文提出了平地機俯仰角預測模型在線辨識方法。由于水田平地過程中,平地鏟運動并不影響車身俯仰角變化,因此本文研究平地機俯仰角在線建模時,未分析平地鏟與動力車身之間的相對運動關系。后續在設計MPC控制器時,將針對性開展鏟姿態與車身俯仰角之間的運動關聯分析,并結合硬件平臺測試系統動態響應,進一步對建模算法精度和執行效率進行檢驗、分析與優化。

[1] 胡煉,林潮興,羅錫文,等. 農機具自動調平控制系統設計與試驗[J]. 農業工程學報,2015,31(8):15-20.Hu Lian, Lin Chaoxing, Luo Xiwen, et al. Design and experiment on auto leveling control system of agricultural implements[J].Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2015, 31(8): 15-20. (in Chinese with English abstract)

[2] 胡煉,羅錫文,林潮興,等. 1PJ-4.0型水田激光平地機設計與試驗[J]. 農業機械學報,2014,45(4):146-151.Hu Lian, Luo Xiwen, Lin Chaoxing, et al. Development of 1PJ-4.0 Laser Leveler Installed on a Wheeled Tractor for Paddy Field[J]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(4): 146-151. (in Chinese with English abstract)

[3] 趙祚喜,羅錫文,李慶,等. 基于MEMS慣性傳感器融合的水田激光平地機水平控制系統[J]. 農業工程學報,2008,24(6):119-124.Zhao Zuoxi, Luo Xiwen, Li Qing, et al. Leveling control system of laser-controlled land leveler for paddy field based on MEMS inertial sensor fusion[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE),2008, 24(6): 119-124. (in Chinese with English abstract)

[4] 胡煉,羅錫文,林潮興,等. 一種作業機具自動調平系統及其控制方法:CN104206063A[P]. 2014-12-17.

[5] 馬敏,張曉輝,宋濤,等. 智能煙草打頂抑芽機控制系統[J]. 農業工程學報,2011,27(2):129-135.Ma Min, Zhang Xiaohui, Song Tao, et al. Control system of intelligentized tobacco chopping and sucker controlling machine[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2011, 27(2): 129-135. (in Chinese with English abstract)

[6] 李君,徐巖,許績彤,等. 懸掛式電動柔性疏花機控制系統設計與試驗[J]. 農業工程學報,2016,32(18):61-66.Li Jun, Xu Yan, Xu Jitong, et al. Design and experiment of control system for suspended electric flexible thinner[J].Transactions of the Chinese Society of Agricultural Engineering(Transactions of the CSAE), 2016, 32(18): 61-66. (in Chinese with English abstract)

[7] 王松林. 噴桿高度自動調節系統設計與試驗研究[D]. 楊凌:西北農林科技大學,2014.Wang Songlin. Design and Experiments on Boom Height Automatic Adjusting System[D]. Yangling: Northwest A&F University, 2014. (in Chinese with English abstract)

[8] 林守利. 隨機海浪作用下船上物體保持水平的預測控制研究[D]. 沈陽:沈陽工業大學,2003.Lin Shouli. The Study on Predictive Control of an Object's Level on a Ship in Random Ocean Waves[D]. Shenyang:Shenyang University of Technology, 2003. (in Chinese with English abstract)

[9] Suh J, Yi K, Jung J, et al. Design and evaluation of a model predictive vehicle control algorithm for automated driving using a vehicle traffic simulator[J]. Control Engineering Practice, 2016, 51: 92-107.

[10] 席裕庚. 預測控制(第2版)[M]. 北京:國防工業出版社,2013.

[11] Mozaffari L, Mozaffari A, Azad N L. Vehicle speed prediction via a sliding-window time series analysis and an evolutionary least learning machine: A case study on San Francisco urban roads[J]. Engineering Science and Technology, an International Journal, 2015, 18(2): 150-162.

[12] Wu Chaozhong, Peng Liqun, Huang Zhen, et al. A method of vehicle motion prediction and collision risk assessment with a simulated vehicular cyber physical system[J]. Transportation Research Part C: Emerging Technologies, 2014, 47: 179-191.

[13] Gallieri M, Ringwood J, Giantomassi A, et al. Predictor design for altitude control of a seaweed harvester[J]. IFAC Conference on Control Applications in Marine Systems,Germany, 2010, 43(20): 301-306.

[14] Koo S, Kim S, Suk J. Model predictive control for UAV automatic landing on moving carrier deck with heave motion[J]. IFAC Proceedings Volumes (IFAC-Papers Online), 48(5): 59-64.

[15] 趙曉莉,馬新明,王舉才,等. 基于新陳代謝 GM(1,1)模型的玉米葉長動態預測[J]. 農業工程學報,2013,29(10):183-189.Zhao Xiaoli, Ma Xinming, Wang Jucai, et al. Dynamic prediction on leaf length of maize based on metabolic model GM(1, 1)[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(10):183-189. (in Chinese with English abstract)

[16] 張軍,張侃諭. 溫室溫度控制系統不確定性與干擾的灰色預測補償算法[J]. 農業工程學報,2013,29(10):225-233.Zhang Jun, Zhang Kanyu. Grey prediction compensation algorithm for the uncertainty and interference of greenhouse temperature control system[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE),2013, 29(10): 225-233. (in Chinese with English abstract)

[17] Bangura M, Mahony R. Real-time model predictive control for quadrotors[J]. IFAC Proceedings Volumes, 2014: 11773-11780.

[18] Alexis K, Nikolakopoulos G, Koveos Y, et al. Switching model predictive control for a quadrotor helicopter under severe environmental flight conditions[J]. IFAC Proceedings Volumes, 2011, 44(1): 11913-11918.

[19] Alexis K, Nikolakopoulos G, Tzes A. Switching model predictive attitude control for a quadrotor helicopter subject to atmospheric disturbances[J]. Control Engineering Practice,2011, 19: 1195-1207.

[20] Dong Yiqun, Zhang Youmin, Ai Jianliang. Full-altitude attitude angles envelope and model predictive control-based attitude angles protection for civil aircraft[J]. Aerospace Science and Technology, 2016, 55: 292-306.

[21] 鄒濤,丁寶蒼,張端. 模型預測控制工程應用導論[M]. 北京:化學工業出版社,2010.

[22] Cryer, Jonathan D, Chan, et al. Time Series Analysis: with Applications in R[M]. Berlin: Springer Science & Business Media, 2008.

[23] 王燕. 應用時間序列分析[M]. 北京:中國人民大學出版社,2015.

[24] 彭秀艷,趙希人,高奇峰. 船舶姿態運動實時預報算法研究[J]. 系統仿真學報,2007,19(2):267-271.Peng Xiuyan, Zhao Xiren, Gao Qifeng. Research on real-time prediction algorithm of ship attitude motion[J]. Journal of System Simulation, 2007, 19(2): 267-271. (in Chinese with English abstract)

[25] 龔嘯. 基于 ARMA遞推算法的電力系統低頻振蕩模式在線辨識研究[D]. 重慶:重慶大學,2011.Gong Xiao. Research on Online Identification for Low Frequency Oscillation Modes in Power System Based on ARMA Recursive Method[D]. Chongqing: Chongqing University, 2011. (in Chinese with English abstract)

[26] 潘迪特,吳憲民. 時間序列及系統分析與應用[M]. 北京:機械工業出版社,1988.

[27] 高奇峰. 艦載機起降指導系統仿真研究[D]. 哈爾濱:哈爾濱工程大學,2006.Gao Qifeng. The Simulated Research for Guidance System of Carrier Plane Taking Off & Landing[D]. Harbin: Harbin Engineering University, 2006. (in Chinese with English abstract)

[28] Ljung L. System identification:theory for the user[M].Beijing: Tisinghua University Press, 2002.

[29] 鐘美霞. 不同采樣頻率時間序列的多步預測比較研究[D].廣州:廣東工業大學,2013.Zhong Meixia. A Comparative Study of the Multi-step Prediction of Time Series with Different Sampling Frequencies[D].Guangzhou: Guangdong University of Technology, 2013. (in Chinese with English abstract)

[30] 楊叔子,吳雅,軒建平. 時間序列分析的工程應用[M]. 武漢:華中科技大學出版社,2007.

[31] 龐中華,崔紅. 系統辨識與自適應控制[M]. 北京:北京航空航天大學出版社,2013.

[32] Mathworks. Matlab Product Documentation[Z]. Natick: 2016.

[33] 蘇孟輝. 基于時間序列的網絡流量監測系統的設計與實現[D]. 廣州:華南理工大學, 2015.Su Menghui. Design and Implementation of Network Traffic Monitoring System Based on Time Series[D]. Guangzhou:South China University of Technology, 2015. (in Chinese with English abstract)

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 人人爽人人爽人人片| 中文字幕亚洲另类天堂| 波多野结衣在线se| 国产超碰在线观看| 婷婷综合缴情亚洲五月伊| 国产成人亚洲综合A∨在线播放 | 欧美视频在线播放观看免费福利资源| 亚洲经典在线中文字幕| 国产视频入口| 国产成人综合久久精品下载| 亚洲成人网在线播放| 自拍欧美亚洲| 欧美中文一区| 亚洲精品无码高潮喷水A| 尤物国产在线| 国产精品亚洲精品爽爽| 毛片一区二区在线看| 国产高清在线丝袜精品一区| 国产乱人伦偷精品视频AAA| 午夜福利在线观看成人| 98精品全国免费观看视频| 乱人伦中文视频在线观看免费| 一区二区三区国产精品视频| 美女无遮挡免费网站| 久久久国产精品免费视频| 欧美成人午夜视频| 日韩人妻精品一区| 综合人妻久久一区二区精品 | a亚洲视频| 欧美亚洲国产精品第一页| 久久久亚洲色| 最新国产麻豆aⅴ精品无| 特级毛片8级毛片免费观看| 亚洲精品色AV无码看| 2018日日摸夜夜添狠狠躁| 亚洲天堂视频网站| 欧美国产中文| 亚洲天堂高清| 免费又黄又爽又猛大片午夜| 久久黄色影院| 国产美女丝袜高潮| 美女被操黄色视频网站| 国产在线高清一级毛片| 亚洲中文字幕久久无码精品A| 亚洲第一成年网| 色欲色欲久久综合网| 日本不卡免费高清视频| 国模视频一区二区| V一区无码内射国产| 国产丝袜精品| 午夜福利亚洲精品| 亚洲天堂视频在线播放| 亚洲精品动漫在线观看| 午夜福利在线观看成人| 国产精品无码制服丝袜| 午夜国产精品视频| 92午夜福利影院一区二区三区| 国产精品毛片一区视频播| 国产人人干| 最新亚洲人成网站在线观看| 日韩精品无码一级毛片免费| 2020国产在线视精品在| 99国产精品免费观看视频| 激情六月丁香婷婷| 精品一区二区三区水蜜桃| 国产在线欧美| 国产女人18毛片水真多1| 亚洲男人的天堂久久香蕉 | 九九九九热精品视频| …亚洲 欧洲 另类 春色| 亚洲精品欧美重口| 91区国产福利在线观看午夜| 91精品国产一区| 国产亚洲视频中文字幕视频| 久久天天躁狠狠躁夜夜2020一| 少妇精品在线| 91系列在线观看| 大乳丰满人妻中文字幕日本| 伊人成人在线| 亚洲黄色成人| 午夜小视频在线| 特级欧美视频aaaaaa|