亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        抗差自適應濾波的導向鉆具動態(tài)姿態(tài)測量方法

        2016-04-19 09:08:40汪躍龍程為彬
        中國慣性技術學報 2016年4期
        關鍵詞:振動測量

        高 怡,汪躍龍,程為彬

        (1. 西安石油大學 陜西省油氣井測控技術重點實驗室,西安 710065;2. 西安石油大學 電子工程學院,西安 710065)

        抗差自適應濾波的導向鉆具動態(tài)姿態(tài)測量方法

        高 怡1,2,汪躍龍1,2,程為彬1,2

        (1. 西安石油大學 陜西省油氣井測控技術重點實驗室,西安 710065;2. 西安石油大學 電子工程學院,西安 710065)

        針對鉆井工程中,導向鉆井工具在近鉆頭振動條件下,由于底部鉆具振動對導向工具姿態(tài)測量造成嚴重干擾,使得姿態(tài)參數(shù)失真,導致姿態(tài)測量不準確的問題。為了消除或削弱振動對空間姿態(tài)測量的不利影響,快速解算出準確的鉆具姿態(tài)參數(shù),提出一種抗差自適應濾波的動態(tài)姿態(tài)測量方法,利用等價權函數(shù)和自適應因子合理的分配信息,有效地濾除鉆具振動對動態(tài)姿態(tài)測量的影響。仿真實驗和實鉆井數(shù)據(jù)試驗證明,采用該算法能濾除近鉆頭振動干擾信號,井斜角控制在5.5°左右,工具面角測量誤差小于6°,有效地提高了導向鉆井工具姿態(tài)參數(shù)的動態(tài)測量精度,解決動態(tài)測量時姿態(tài)參數(shù)測不準的問題。

        導向鉆井;動態(tài)測量;鉆具振動;抗差估計;自適應濾波

        隨著石油工業(yè)的不斷發(fā)展和油氣勘探開發(fā)難度的增大,國內外競爭也愈來愈激烈。載體空間姿態(tài)測量在石油工業(yè)領域中的地位也越來越顯赫。這就意味著,對姿態(tài)參數(shù)(井斜角、工具面角和方位角)測量的實時性、精確性以及連續(xù)、動態(tài)測量的要求也越來越高[1-2]。井下動態(tài)空間姿態(tài)參數(shù)的實時測量是導向工具實現(xiàn)實時導向控制的前提條件。然而,在實際鉆井工程中,由于底部鉆具振動對導向工具姿態(tài)測量造成嚴重干擾,使得姿態(tài)參數(shù)嚴重失真,導致姿態(tài)測量不準確。目前普遍采用隨鉆測量(Measurement While Drilling,MWD)技術對井下姿態(tài)進行測量。該方法雖然能得到準確的姿態(tài)參數(shù),但要求姿態(tài)測量時必須停止鉆進,屬于靜態(tài)測量,存在時效低、成本高和風險大等問題。如何在不停鉆條件下獲得實時性強、精度高的連續(xù)動態(tài)姿態(tài)測量,是目前該方向的研究熱點及難點之一。

        文獻[3]提出了一種采用卡爾曼濾波的新狀態(tài)空間法連續(xù)實時測量鉆井軌跡,針對旋轉導向鉆井和自動垂直鉆井建立一套捷聯(lián)測量系統(tǒng)[4],采用卡爾曼濾波估計系統(tǒng)狀態(tài),但是該方法沒有考慮測量傳感器信號中包含的大量振動加速度。楊全進等提出一種新的旋轉導向鉆井動態(tài)位置測量方法[5],介紹了常用的石油鉆井定向傳感器的結構和測量參數(shù),提供了一種在鉆柱旋轉情況下的動態(tài)位置測量技術,詳細分析了如何計算實時精確的重力工具面角。采用無跡卡爾曼濾波方法消除傳感器數(shù)據(jù)中的干擾噪聲,提高旋轉導向鉆具姿態(tài)測量的準確性[6]。徐寶昌等在此基礎上又提出了一種旋轉導向系統(tǒng)有色噪聲的改進無跡卡爾曼濾波方法。該方法利用四元數(shù)理論構造觀測方程和時變狀態(tài)方程,實時解算出鉆具姿態(tài),有效的濾除了加速度傳感器數(shù)據(jù)中的有色噪聲[7]。但是該算法限制條件較多,只是將鉆具軸向振動信號作為主要干擾進行分析。Elgizawy和Jurkov等提出一種基于慣性導航系統(tǒng)的測量井斜與方位角的傳感器封裝,通過仿真驗證基于慣性導航的傳感器封裝可以實現(xiàn)高精度連續(xù)MWD測量[8-10]。Yanshun等人基于簡易慣性單元開發(fā)了一種相似的MWD裝置[11]。由于井下工況復雜,鉆柱旋轉向下運動常伴隨著各種振動現(xiàn)象,鉆具振動嚴重影響了傳感器的測量精度,導致測量誤差較大。

        在上述研究的基礎上,針對導向鉆井工具動態(tài)測量受鉆具振動的影響而導致測量不準確的問題,提出一種抗差自適應濾波的動態(tài)空間姿態(tài)測量方法。通過分析鉆具振動對姿態(tài)測量的影響,并吸收抗差估計和自適應濾波的優(yōu)點,利用抗差等價權矩陣自適應的確定量測信息,通過自適應因子調整狀態(tài)模型信息對狀態(tài)參數(shù)的整體貢獻,從而消除鉆具振動對動態(tài)姿態(tài)測量的影響,獲得實時性強、精度高的姿態(tài)參數(shù),提高鉆井效率,降低鉆井風險。

        1 鉆具振動對姿態(tài)測量的影響

        鉆井振動與鉆柱及其組成部分的動力特性相關,通過監(jiān)測鉆柱的動力特性,削弱振動對鉆井測量的影響,可以有效的提高鉆井效率,降低鉆井事故,因此,鉆井的振動問題成為目前的研究熱點之一。

        在實際鉆井過程中,鉆頭切削巖層、鉆柱與井壁的碰撞等會使鉆具產生橫向振動、縱向振動和扭轉振動等,這些振動嚴重的影響了測量傳感器輸出信號的正確性。而橫向振動又是石油鉆井下部鉆具(鉆頭及近鉆頭鉆鋌和鉆桿)的主要振動形式之一[12]。章楊烈指出橫向振動對鉆具的損害大于縱向振動和扭轉振動[13]。采用常用的牙輪鉆頭進行鉆進,鉆具橫向振動具有較強的隨機性,其強度屬于劇烈分區(qū),近鉆頭橫向振動強度隨鉆壓、鉆速增大而增大。鉆頭振動信號是一種非穩(wěn)態(tài)連續(xù)隨機振動,具有較寬的頻帶,影響因素復雜。

        鉆柱的橫向振動的傳播方程如下:

        式中:ρ表示鉆柱質量密度,A表示鉆柱橫截面積,ω表示鉆柱微元偏離井眼軸線的橫向位移,z為軸向坐標,S為源于X方向的剪切力,fF為外部力。

        式中:力矩M為

        設彈性模量為E,橫截面積對中性軸的慣性矩為I。

        將式(2)和(3)代入到式(1)中,得到關于橫向振動的傳播方程:

        根據(jù)上述動力學方程即可對鉆具橫向振動的傳播進行數(shù)值模擬。

        在導向鉆井系統(tǒng)的動態(tài)空間姿態(tài)測量時,由于三軸磁通門數(shù)據(jù)中包含磁干擾,三軸加速度計測量數(shù)據(jù)中包含大量的振動加速度,嚴重影響組合鉆具井斜及方位姿態(tài)測量精度。為了研究井下振動對姿態(tài)測量信號的影響,通過分析橫向振動信號特性,根據(jù)橫向振動的傳播方程,井下鉆具橫向振動加速度響應曲線如圖1所示。

        圖1 鉆具橫向振動加速度響應曲線Fig.1 Transverse vibration acceleration response curve of drilling tool

        2 抗差自適應濾波動態(tài)測量姿態(tài)解算

        抗差自適應濾波的基本思想是當觀測值存在異常時,對觀測值采用抗差估計原則,能夠控制觀測異常的影響;當動力學模型存在異常誤差時,將動力學模型信息作為一個整體,采用統(tǒng)一的自適應因子調整動力學模型信息對狀態(tài)參數(shù)的整體貢獻[14-16]。

        2.1 導向鉆井工具姿態(tài)測量系統(tǒng)動態(tài)建模

        選取地理坐標系以“東北天(ENU坐標系)”為順序構成的右手直角坐標系,即E軸沿當?shù)厮矫嬷赶驏|,N軸沿當?shù)厮矫嬷赶虮?,U軸沿當?shù)卮咕€指向上。再以鉆具三個基本軸建立鉆具坐標系(XYZ坐標系),取Z軸和鉆具的軸線方向一致,X軸和Y軸在同一水平面內且互相垂直,Z軸跟二者垂直并構成右手直角坐標系。

        在以上選取的坐標系下,各姿態(tài)角如圖2所示,其中,H為水平面,V為鉆孔彎曲平面,P代表鉆具橫截面。方位角ψ為磁北方向沿逆時針方向到Z軸在水平面的投影間的夾角,其范圍在 0°~360°之間,井斜角θ為鉆進軸Z軸與水平面所成的夾角,規(guī)定向下為正,反之為負,其范圍為-90°~90°,工具面向角γ則為鉆孔橫截面內由鉆孔高邊到Y軸所成的角度,范圍在 0°~360°之間。這樣,我們就準確的定義了井下鉆具的方位角ψ、井斜角θ和工具面向角γ,且角度的正向都符合右手系原則。

        根據(jù)導航理論中的歐拉定理,載體在空間中的姿態(tài)可用相對于地理坐標系有限次的轉動來表示,每次轉動的角度即為航向角、俯仰角和橫滾角。同理,井下鉆進過程中,鉆具在空間的任一姿態(tài)都可以用相對于地理坐標系的一系列旋轉來表示,旋轉的角度為方位角、傾角和工具面向角。在XYZ坐標系中安裝三軸加速度計和三軸磁通門,如圖3所示。

        圖2 坐標系中表示的姿態(tài)角Fig.2 Attitude angles in a coordinate system

        圖3 加速度計和磁通門組合測量安裝示意圖Fig.3 Combined measurements of accelerometer and fluxgate

        設Gx、Gy和Gz分別為地球重力加速度g向鉆具坐標系上投影:

        由式(5)可以解算得到井斜角和工具面角為

        根據(jù)上述理論,建立導向鉆井工具姿態(tài)測量的動態(tài)數(shù)學模型,給出狀態(tài)方程和量測方程。

        假定tk時刻的動力學模型方程為

        式中:xk和xk-1分別為tk和tk-1時刻的n維狀態(tài)參數(shù)向量,為n×n維狀態(tài)轉移矩陣;wk為p維動力學模型誤差向量,其數(shù)學期望為0,協(xié)方差矩陣為

        顯然,wk為高斯白噪聲序列。

        設tk時刻的量測方程為

        式中:yk為tk時刻的m維觀測向量;Hk為m×n維測量矩陣,也稱為觀測矩陣;vk為m維觀測誤差向量,其數(shù)學期望為0,協(xié)方差矩陣為

        顯然,vk為高斯白噪聲序列。在i=k時,的協(xié)方差矩陣分別為互不相關。

        狀態(tài)向量為

        式(9)表明直接將鉆具姿態(tài)參數(shù)作為狀態(tài)向量,而非姿態(tài)誤差作為狀態(tài)。

        當鉆具坐標系相對地理坐標線加速度為零的情況下,三軸加速度計的輸出記為即此時鉆具坐標系的加速度輸出也就是地球的重力加速度g在鉆具坐標系上的投影。三軸磁通門的輸出記為

        將加速度計、磁通門測量值與當?shù)刂亓?、磁場強度在鉆具坐標系下投影的差值作為量測,共6維,則量測方程如下

        式中:ya三軸加速度計測量時輸出數(shù)據(jù)矩陣,ym為磁強計測量值;G為重力場,G=[00g]T;M為磁力場,M=[M10M2]T,M1和M2分別為地理坐標系中地磁場北向分量和垂直分量的測量值,M1=ξcosθ,M2=ξsinθ,ξ為當?shù)氐卮艔姸?;C為姿態(tài)旋轉矩陣,。

        2.2 導向鉆井工具抗差自適應濾波的動態(tài)姿態(tài)測量

        采用抗差自適應濾波器實現(xiàn)各傳感器數(shù)據(jù)融合,計算導向鉆具工具姿態(tài)測量的框圖如圖4所示,其中,ψ、θ和γ分別為方位角、井斜角和工具面角。

        圖4 姿態(tài)參數(shù)估計框圖Fig.4 Attitude parameter estimation

        考慮模型式(7)和(9),系統(tǒng)狀態(tài)預測向量為

        式中:Vk為觀測殘差向量,為預測殘差向量。Vk和的協(xié)方差矩陣分別為

        合理地選擇自適應因子不但能夠自適應地平衡動力學模型預測信息與量測信息的權比,而且能夠控制動力學模型擾動異常對濾波解的影響?;陬A測殘差誤差判別統(tǒng)計量的抗差自適應因子函數(shù)為

        等價權矩陣為

        式中:k0∈(1,1.5),k1∈(3,8),Vk為觀測值的殘差向量。

        由式(14)~(16)和式(19)~(21)可構造如下極值:

        式中:Kk為增益矩陣,根據(jù)矩陣恒等式,可表示為

        對量測信息采用抗差估計,自適應的確定觀測噪聲協(xié)方差矩陣,并利用自適應因子調節(jié)狀態(tài)噪聲的協(xié)方差矩陣,因此,可以有效的控制量測異常和動態(tài)模型噪聲異常對空間狀態(tài)參數(shù)估值的影響。

        3 實驗及分析

        3.1 數(shù)值計算仿真及分析

        實驗室地理條件為北緯34.24°,東經108.99°,地球自轉角速度為 15 (°)/h,磁傾角為 55.4°,磁場強度為52.5 μT,地球重力加速度為9.8 m/s2。在實驗室條件下,根據(jù)測斜校驗裝置測量得到一組理想的實驗數(shù)據(jù)。

        將模擬振動試驗中得到的三個加速度輸出信號經過解算后送到抗差自適應濾波器中進行濾波處理,得到濾波前后的井斜角和工具面角對比曲線分別如圖 5和圖6所示。

        圖5 數(shù)值計算濾波前后井斜角對比曲線Fig.5 Contrast curves of deviation angles before and after filtering by numerical calculation

        圖6 數(shù)值計算濾波前后工具面角對比曲線Fig.6 Contrast curves of tool face angles before and after filtering by numerical calculation

        從圖5和圖6中可以非常直觀地看到,經過抗差自適應濾波處理后的井斜角和工具面角波動范圍明顯變小,濾波后的數(shù)據(jù)曲線更加平滑,橫向振動干擾信號基本濾除,濾波效果明顯。

        3.2 實鉆井數(shù)據(jù)試驗及分析

        為了進一步驗證算法的有效性,采用實鉆井數(shù)據(jù)進行試驗和分析,分別對提出的抗差自適應濾波算法與卡爾曼濾波(KF)算法進行對比。

        實驗數(shù)據(jù)來源于 2012年陜北某井的實鉆井過程中的某一段數(shù)據(jù),鉆井深度為1876m,鉆速為45m/h,鉆壓1~10 t,溫度40℃。

        從圖7和圖8可以看出,濾波前井斜角和工具面角波動較大,采用卡爾曼濾波器后井斜角誤差仍較大,而經過抗差自適應濾波器后,井斜角控制在5.5°左右,與靜態(tài)測量下測得的井斜角基本相同。工具面角經過濾波后曲線更加平滑,干擾信號基本濾除,工具面角誤差小于 6°。實驗結果表明,提出的抗差自適應濾波能夠抑制橫向振動對導向鉆井工具動態(tài)姿態(tài)測量的影響。

        圖7 實鉆井數(shù)據(jù)濾波前后井斜角對比曲線Fig.7 Contrast curves of deviation angles before and after filtering by real drilling data calculation

        圖8 實鉆井數(shù)據(jù)濾波前后工具面角對比曲線Fig.8 Contrast curves of tool face angles before and after filtering by real drilling data calculation

        4 結 論

        針對近鉆頭強振動對姿態(tài)參數(shù)所造成的測不準問題,提出采用抗差自適應濾波方法來減小近鉆頭強振動對導向鉆井工具姿態(tài)參數(shù)測量的影響。從橫向振動傳播方程著手,根據(jù)振動干擾信號特性,設計了針對性的模擬振動實驗。分析振動加速度為姿態(tài)測量的影響,通過抗差自適應濾波器有效地濾除近鉆頭振動信號,解決導向鉆井工具姿態(tài)參數(shù)的動態(tài)測量問題,提高了導向鉆井工具姿態(tài)測量的準確性。

        這里只是考慮了鉆具的橫向振動對導向鉆井工具姿態(tài)測量的影響,沒有考慮縱向振動和扭轉振動的影響,下一步工作考慮如何對這兩種振動進行模擬,采用更合適的濾波算法消除干擾信號,以期獲得更精確的動態(tài)測量姿態(tài)參數(shù)。

        (References):

        [1] 薛啟龍, 丁青山, 黃蕾蕾. 旋轉導向鉆井技術最新進展及發(fā)展趨勢[J]. 石油機械, 2013, 41(7): 1-6. Xue Qi-long, Ding Qing-shan, Huang Lei-lei. The latest progress and development trend of rotary steering drilling technology[J]. China Petroleum Machinery, 2013, 41(7): 1-6.

        [2] 汪躍龍, 王海皎, 康思民, 等. 導向鉆井穩(wěn)定控制平臺的反饋線性化控制[J]. 石油學報, 2014, 35(5): 952-957. Wang Yue-long, Wang Hai-jiao, Kang Si-min, et al. Output feedback linearization of servo platform for rotary steering drilling system[J]. Acta Petrolei Sinica, 2014, 35(5): 952-957.

        [3] Xue Qi-long, Leung H, Wang Ruihe, et al. Continuous real-time measurement of drilling trajectory with new state space models of kalman filter[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(1): 144-154.

        [4] 薛啟龍, 王瑞和, 孫峰, 等. 捷聯(lián)式旋轉導向井斜方位動態(tài)解算方法[J]. 中國石油大學學報(自然科學版), 2012, 36(2): 93-97, 107. Xue Qi-long, Wang Rui-he, Sun Feng, et al. Dynamic solution approach to inclination and azimuth of strapdown rotary steerable system[J]. Journal of China University of Petroleum, 2012, 36(2): 93-97, 107.

        [5] Yang Quan-jin. A new method for dynamic position measurement while drilling string rotating[J]. Applied Mechanics and Materials, 2012(152-154): 1102-1105.

        [6] 楊全進, 徐寶昌, 左信, 等. 旋轉導向鉆具姿態(tài)的無跡卡爾曼濾波方法[J]. 石油學報, 2013, 34(4): 1168 -1175. Yang Quan-jin, Xu Bao-chang, Zuo Xin, et al. An unscented Kalman filter method for attitude measurement of rotary steerable drilling assembly[J]. Acta Petrolei Sinica, 2013, 34(4): 1168-1175.

        [7] 徐寶昌, 楊全進, 蔣海旭. 旋轉導向系統(tǒng)有色噪聲的改進無跡卡爾曼濾波方法[J]. 中國石油大學學報(自然科學版), 2015, 39(2): 157-163. Xu Bao-chang, Yang Quan-jin, Jiang Hai-xu. Improved unscented Kalman filtering method for colored noises of rotary steerable system[J]. Journal of China University of Petroleum, 2015, 39(2): 157-163.

        [8] Elgizawy M, Noureldin A, Georgy J, et al. Wellbore serveying while drilling based on kalman filtering[J]. American Journal of Engineering and Applied Sciences, 2010, 3(2): 240-259.

        [9] Elgizawy M, Noureldin A, El-Sheimy N. Continuous wellbore surveying while drilling utilizing MEMS gyroscopes based on Kalman filtering[C]//Proceedings of 2010 SPE Annual Technical Conference & Exhibit. Florence, Italy, 2010: 5416-5428.

        [10] Jurkov A S, Cloutier J, Pecht E, et al. Experimental feasibility of the in drilling alignment method for inertial navigation in measurement while drilling[J]. IEEE Transactions on Instrumentation and Measurement, 2011, 60(3): 1080-1090.

        [11] Zhang Yanshun, Wanf Shuwei, Fang Jiancheng. Measurement while drilling instrument based on predigested inertial measurement unit[J]. IEEE Transactions on Instrumentation and Measurement, 2012, 61(12): 3295-3302.

        [12] 祝效華, 胡志強. 基于鉆頭破巖鉆井的下部鉆具橫向振動特性研究[J]. 振動與沖擊, 2014, 33(17): 90-93. Zhu Xiao-hua, Hu Zhi-qiang. Lateral vibration characteristics analysis of a bottom hole assembly based on interaction between bit and rock[J]. Journal of Vibration and Shock, 2014, 33(17): 90-93.

        [13] 章楊烈. 鉆柱運動學與動力學[M]. 北京: 石油工業(yè)出版社, 2001. Zhang Yang-lie. Kinematics and dynamics of drill string [M]. Beijing: Petroleum Industry Press, 2001.

        [14] Yang Y X, Zhang S. Adaptive fitting of systematic errors in navigation[J]. Journal of Geodesy, 2005, 79(1-3): 43-49.

        [15] Yang Y, Cui X. Adaptively robust filter with multi adaptive factors[J]. Survey Review, 2008, 40(309): 260-270.

        [16] Gao S, Zhong Y, Li W. Robust adaptive filtering method for SINS/SAR integrated navigation system[J]. Aerospace Science and Technology, 2011, 15(6): 425-430.

        [17] 孫 霄,董景新,顧啟泰. 近鉆頭測斜器最優(yōu)八位置標定法[J]. 中國慣性技術學報, 2014, 22(1): 5-8. Sun Xiao, Dong Jing-xin, Gu Qi-tai. Optimal 8-position calibration for at-bit inclinometer[J]. Journal of Chinese Inertial Technology, 2014, 22(1): 5-8.

        Robust adaptive filtering method for dynamic attitude measurement of steering drilling

        GAO Yi1,2, WANG Yue-long1,2, CHENG Wei-bin1,2
        (1. Shaanxi Key Laboratory of Measurement and Control Technology for Oil and Gas Wells, Xi’an Shiyou University, Xi’an 710065, China; 2. School of Electronic Engineering, Xi’an Shiyou University, Xi’an 710065, China)

        In drilling process, there are serious interferences in the attitude measurement for the steering drilling tool under conditions of near-bit vibration of steering drilling tool, which is due to the bottom drill’s string vibration and would lead to inaccurate attitude measurement. In order to eliminate or weaken these adverse vibration influences on the space attitude measurement and calculate the accurate attitude parameters of drilling tool, a dynamic attitude measurement algorithm with robust adaptive filter is proposed which can effectively eliminate the influence of drilling tool vibration on dynamic attitude measurement and distribute information reasonably by equivalent weight function and adaptive factor. Experimental results and comparison analysis demonstrate that the proposed robust adaptive filtering can eliminate the vibration disturbance signal of the near-bit, the deviation angle can be controlled to about 5.5°, and the tool face angle error is less than 6°. It improves the dynamic measurement accuracy of attitude parameters and solves the problem of uncertainty in dynamic attitude measurement of steering drilling tool.

        steering drilling; dynamic measurement; drill string vibration; robust estimation; adaptive filter

        TP301.6

        :A

        2016-04-15;

        :2016-07-28

        國家自然科學基金(51604226,61174191);陜西省教育廳重點實驗室科研計劃項目(16JS090);中國石油科技創(chuàng)新基金(2015D-5006-0307);西安石油大學青年科技創(chuàng)新基金項目(2015BS44)

        高怡(1978—),女,博士,講師,從事油氣井測控技術、控制理論與控制工程。E-mail: gy@xsyu.edu.cn

        1005-6734(2016)04-0437-06

        10.13695/j.cnki.12-1222/o3.2016.04.004

        猜你喜歡
        振動測量
        振動的思考
        科學大眾(2023年17期)2023-10-26 07:39:14
        噴水推進高速艇尾部振動響應分析
        This “Singing Highway”plays music
        把握四個“三” 測量變簡單
        滑動摩擦力的測量和計算
        滑動摩擦力的測量與計算
        振動攪拌 震動創(chuàng)新
        中國公路(2017年18期)2018-01-23 03:00:38
        中立型Emden-Fowler微分方程的振動性
        測量的樂趣
        測量
        女同久久精品国产99国产精| 久久精品无码中文字幕| 日韩精品久久久一区| 久久精品日本美女视频| 91九色播放在线观看| 五月av综合av国产av| 欧美日韩亚洲国内综合网| 免青青草免费观看视频在线| 小黄片免费在线播放观看| 午夜精品久久久久久久| 免费无码午夜福利片69| 男性一插就想射是因为啥| 青青草视频在线播放观看| 亚洲色偷偷综合亚洲avyp| 国产精品欧美成人| 国产未成女年一区二区| 成人国产av精品麻豆网址| 国产精品国产三级国产av品爱网| 亚洲欧美日韩一区二区三区在线| 国产精品亚洲综合色区韩国| 乳乱中文字幕熟女熟妇| 亚洲一区av在线观看| 丝袜足控一区二区三区| 日本第一区二区三区视频| 精品国产精品三级在线专区| 又黄又硬又湿又刺激视频免费| 日韩精品区欧美在线一区| 亚洲中文字幕在线精品2021| 无码人妻久久久一区二区三区| 午夜精品久久久久成人| 国产精品久久久久9999小说| 国产成人综合久久精品推| 一区二区三区精品偷拍av| 99久久精品费精品国产一区二| 久久人人爽人人爽人人av| 2020久久精品亚洲热综合一本| 国产黄色一区二区在线看| 亚洲av永久无码精品放毛片| 久久精品国产91久久性色tv| 国产高清大片一级黄色| 伊人大杳焦在线|