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

基于小波包的EITD 風力發電機組齒輪箱故障診斷

2015-06-06 07:28:24鄢小安
動力工程學報 2015年3期
關鍵詞:故障診斷故障信號

向 玲, 鄢小安

(華北電力大學 能源動力與機械工程學院,河北保定071003)

風能作為一種清潔能源,近年來引起了人們的廣泛關注.隨著風力發電機組的增加,其故障出現得越來越頻繁,對風力發電機組的故障診斷也就顯得越來越重要[1-3].齒輪箱是風力發電機組的重要組成構件之一,其結構復雜,在運行過程中容易出現故障,且一旦其發生故障,將直接影響整個系統的運行狀態和性能[4].因此,開展對風力發電機組齒輪箱的故障診斷研究,及時診斷出系統的早期故障并進行維修,對確保風力發電機組安全穩定運行具有重要意義.

小波變換最早由Morlet提出,其劃分尺度是按二進制變化的,在大尺度(低頻)時,頻率分辨率高,時間分辨率低,在小尺度(高頻)時,頻率分辨率低,時間分辨率高,但對于一些復雜的信號,小波變換不能很好地適應信號的特點[5-6].于是Wickerhauser[7]提出了小波包變換(WPT),小波包變換可以將頻帶進行多層次劃分,不僅可以分解低頻部分,還可以分解高頻部分,提高信號的時頻分辨率,是一種自適應的非線性分析方法.因此,利用小波包變換分析方法進行故障特征提取有著廣泛的應用價值.

Frei等[8]提出了一種自適應時頻分析方法——固有時間尺度分解(ITD)方法,該方法可將一個復雜的非線性非平穩信號分解為若干個固有旋轉分量(PRC)和一個殘余項之和.雖然文獻[8]中對ITD方法與經驗模態分解(EMD)方法進行了比較,ITD方法在端點處理和分解速度上有著明顯優勢,但文獻[8]中沒有闡述ITD 方法及PR 分量的物理意義;同時,ITD 方法使用線性變換方法進行信號分解,有可能使得到的PR 分量出現毛刺而失真.因此,筆者基于EMD 中的三次樣條插值方法和ITD 方法中的線性變換方法,提出了一種集成固有時間尺度分解(EITD)方法.目前國內學者[9-11]已將固有時間尺度分解方法用于故障診斷中并取得了一定進展.

為了提高風力發電機組齒輪箱故障診斷的正確性,筆者提出了一種基于小波包的EITD 風力發電機組齒輪箱故障診斷方法:首先利用EITD 對齒輪加速度振動信號進行分解;然后選擇相關性最大的PR 分量進行小波包分解,得到一系列小波包系數,分別計算各小波包系數的能量,選擇能量比重較大的小波包系數重構原PR 分量;最后計算重構PR分量的關聯維數,對齒輪的工作狀態和故障類型進行識別.

1 EITD方法

1.1 ITD 方 法

ITD 方法是一種新的時頻分析方法,該方法與EMD 和局部均值分解(LMD)方法一樣具有自適應性.對于任意給定的一個非線性非平穩信號Xt,應用ITD 方法能夠將其分解為若干個瞬時頻率具有物理意義的PRC和一個單調趨勢項之和.該方法具體分解過程如下:

(1)定義一個算子L 用于提取低頻基線信號,使得從原信號中去除一個基線后得到的剩余信號成為一個固有旋轉分量,信號Xt的一次分解為

式中:Lt和Ht分別為基線信號和固有旋轉分量.

(2)確定信號Xt(t≥0)的局部極值Xk及其對應的時刻τk{k=1,2,…,M},其中M 為極值總數.定義τ0=0,為便于分析,設LXt=Lt,HXt=Ht.

其中

式中:Lk和Lk+1分別為第k 個和第k+1個基線控制點;0<α<1,通常α取0.5.

1.2 EITD方法

根據ITD 的原理可知,ITD 方法對基線的定義是基于信號的線性變換而來的,該方法可能會導致分解信號波形出現毛刺而失真.EMD 方法使用三次樣條插值擬合上下包絡,存在過包絡和欠包絡等缺陷.考慮到2 種方法存在的不足,提出了EITD 方法,其具體分解過程如下:

(1)確定原信號Xt所有局部極值點,與ITD方法一樣,根據式(2)、式(3)提取各基線控制點Lk.

(2)采用鏡像對稱延拓法對時間序列信號進行端點處理,獲得左右兩端極值點(τ0,X0)和(τM+1,XM+1).令k分別為0和M-1,按照式(2)和式(3)求出L1與LM的值.然后使用三次樣條插值來擬合所有的Lk,得到基線信號L1(t).

(3)將基線從原信號中分離出來,得到h1(t):

理想情況下,h1(t)為一個固有旋轉分量,即h1(t)=βPRC1,若h1(t)不滿足固有旋轉分量條件,即基線,則將h1(t)作為原信號重復以上過程,循 環k 次,直 到h1k(t)為 一 個PR 分 量,即h1k(t)=βPRC1.實 際 中,可 以 設 置 一 個 變 量Δ,當時迭代結束.

(4)將βPRC1從原信號中分離出來,可得到一個新信號r1(t),即

(5)將r1(t)作為原信號重復以上過程,得到原信號Xt的第2個滿足PRC 條件的分量βPRC2,重復循環n-1次,得到Xt的第n 個滿足PRC條件的分量βPRCn,直到rn(t)為一單調函數或常數為止.至此原信號Xt被分解成n 個固有旋轉分量βPRCn 和一個單調函數rn(t)之和,即

應用EITD 方法的一個關鍵問題就是如何選擇PR 分量的迭代終止判據Δ,這關系到分解的效果和迭代的次數,目前常用的迭代終止判據有標準差法和三參數法[7]等.筆者在PR 分量定義的基礎上,采用三參數法確定Δ 的取值.經多次仿真實驗分析發現,當信號基線控制點Lk+1與其對應極值點Xk+1的比值小于某取值范圍,即0.005 Xk+1<Δ<0.05 Xk+1時,分解得到的分量滿足PR 分量條件,此時迭代次數通常為4~10次.

1.3 端點效應的評價指標

由ITD 方法的理論可知,ITD 方法與EMD 和LMD 方法一樣,也存在端點效應,這種端點效應將影響原信號分解的精度,給原信號增加一些虛假成分,使得各PRC的總能量發生相應改變.原信號經ITD 方法分解后得到的每個PRC 的能量之和應該等于原信號的能量,因此可以通過比較原信號經ITD 方法分解前后產生的能量變化來評價端點效應對原信號的影響程度.

為了更好地比較,一般先求取原信號Xt和經ITD 方法分解得到的各分量的有效值,即

式 中:E 為 信 號 有 效 值;Xi為 信 號 序 列;n 為 信 號 的采樣點數.當Xi=Xt時,E=Ex;當時,E=Ep.

根據式(9)比較各PRC有效值的總和與原信號有效值,獲得一個評價指標θ:

式中:Ex為原信號Xt的有效值;Ep為第p 個PRC分量的有效值;k+1為PRC 分量的總數,包括分解的殘余項.

根據定義,θ≥0,且θ的值越大,表示分解精度越低,端點效應的影響越大;θ=0,表示端點效應對ITD 方法沒有影響.

為評價端點效應對EITD 方法的影響程度,考察如下仿真信號

該信號由2個調幅-調頻分量組成,其時域波形如圖1所示.對信號x(t)分別通過未經過鏡像延拓處理的ITD(端點處理前)、經過鏡像延拓處理的ITD(端點處理后)和EITD 方法分解,所得分解結果如圖2~圖4所示.其中前2個分量代表分解后的真實分量,最后一個分量r為殘余項.為了更好地比較3種方法的分解效果,可考察它們的端點效應評價指標θ,分解前后2個PR 分量與原信號的互相關系數ρ1、ρ2 及分解速度.具體分解效果評價指標如表1所示.

圖1 信號x(t)的時域波形Fig.1 Time domain waveform of the signal x(t)

由圖2~圖4和表1可明顯發現端點處理后的ITD 和EITD 分解結果精確地體現了原信號的頻率組成成分,而端點處理前ITD 分解結果的第2 個PR 分量出現了明顯的變形.與圖3 相比,圖4 中EITD 分解結果的第2個PR 分量表現得更光滑,而ITD 方法分解得到的第2個PR 分量因有毛刺而失真.由表1可知,ITD 方法的端點效應評價指標θ較大,分解速度較慢,前2個分量與原信號的相關性較小.由此可知,EITD方法在基線的計算中采用一次三次樣條插值,這樣既可以保證分解得到的分量更光滑,且分解速度更快,又避免了ITD 方法采用線性變換求取基線使得波形出現失真的現象.可見,EITD 方法能有效地解決端點效應問題,并且能快速準確地分解信號.

圖2 信號x(t)的ITD分解結果(端點處理前)Fig.2 ITD result generated from the signal x(t)before endpoint processing

圖3 信號x(t)的ITD分解結果(端點處理后)Fig.3 ITD result generated from the signal x(t)after endpoint processing

圖4 信號x(t)的EITD分解結果Fig.4 EITD result generated from the signal x(t)

2 小波包變換

小波包變換是將信號在小波包函數系上展開,也就是求信號與小波包函數的內積,可以更加精確地對信號進行分析與重構,屬于一種自適應的非線性分析方法.小波包分解不僅對低頻部分進行分解,而且對高頻部分進行二次分解,提高了信號的時頻分辨率.小波包分解的算法為

式中:xj,m(a)為信號x(a)在尺度j上的小波包分解系數;j為尺度;m 為頻帶;a 為信號x(a)中各點的時域位置.

為便于理解,以一個3層分解的小波包分解樹來說明其分解過程,如圖5所示.

圖5 小波包分解過程Fig.5 Process of the wavelet packet decomposition

圖5 中,x(a)為原信號,x1,0(a)為小波包分解的第1層低頻系數,x1,1(a)為小波包分解的第1層高頻系數,依次類推,第3層按低頻到高頻共分解為8個小波包系數,即總信號可表示為

3 基于小波包的EITD 風力發電機組齒輪箱故障診斷

針對風力發電機組齒輪箱故障診斷的復雜性,提出了一種基于小波包的EITD 風力發電機組齒輪箱故障診斷方法,其診斷流程如圖6所示.該方法可具體描述為:

(1)分別選取正常、第2級平行軸齒輪點蝕故障及第3級平行軸齒輪磨損和斷齒混合故障的振動加速度信號作為研究對象.

(2)利用EITD 對每個振動信號進行分解,得到若干個PR 分量,計算各PR 分量與原振動信號的互相關系數,并據此選出相關性最大的PR 分量作為主PR 分量,進行第i層小波包分解并得到一系列小波包系數(WPC).

(3)計算各小波包系數的能量分布并選取能量較大的系數重構原PR 分量.計算表達式如下

式中:E(βWPCi,γ)表示各小波包系數的能量;K 為小波包空間位置標識;γ表示頻帶.

(4)分別計算重構PR 分量的關聯維數,對齒輪的工作狀態和故障類型進行識別.值得注意的是,延遲時間和嵌入維數的選擇對關聯維數的計算結果影響很大,延遲時間選取過大或過小都無法反映系統 的運 動 特 征[12-13].筆 者 采 用 關 聯 積 分 法[14](C-C法)來確定重構相空間中合適的延遲時間.

圖6 基于小波包的EITD故障診斷流程圖Fig.6 Flowchart of the fault diagnosis based on EITD-WPT

4 應用實例

以張家口蔚縣風電廠1.5 MW 風力發電機組齒輪箱為研究對象.利用SCADA 系統采集齒輪箱加速度振動信號,采樣頻率為32 768 Hz,采樣點數為16 384.圖7所示為該風力發電機組齒輪箱傳動系統的結構圖.該傳動系統具有3級傳動裝置,第1級為行星輪系,采用內嚙合傳動;第2級和第3級為平行軸系,采用外嚙合傳動.太陽輪齒數z0為20;齒輪箱內齒圈齒數z1為100;第2級傳動大齒輪齒數z3為100,小齒輪齒數z5為23;第3級傳動大齒輪齒數z4為93,小齒輪齒數z6為25.齒輪箱在轉速1 400r/min下發生故障,主要包括第2級平行軸(中間軸)小齒輪點蝕故障,第3級平行軸(高速軸)小齒輪磨損和斷齒混合故障.故障齒輪如圖8和圖9所示.

圖7 風力發電機組齒輪箱傳動系統的結構圖Fig.7 Structural diagram of the wind turbine gearbox transmission system

圖8 點蝕故障齒輪Fig.8 Gear with pitting fault

圖9 磨損和斷齒混合故障齒輪Fig.9 Gear with wear and broken teeth mixed fault

圖10(a)~圖10(c)分別為齒輪在正常、第2級平行軸齒輪點蝕故障、第3級平行軸齒輪磨損和斷齒混合故障下的加速度振動信號.由圖10 可知,3種狀態下的時域波形不同,由于含有大量噪聲,使得信號的沖擊調制特征被覆蓋,因此不能很好地提取特征信息.風電機組齒輪振動信號經EITD 分解所得的主要PR 分量與常見的自適應時頻分解方法(如EMD 方法)分解得到的IMF 分量相比,保留了更多的頻率信息并攜帶豐富的故障信息,因此可以提取主要PR 分量的特征信息對風電機組齒輪進行故障診斷.考慮到風電機組齒輪振動信號往往受到強噪聲信號的干擾和背景噪聲的影響,可采用小波包變換優良的降噪性能對主要PR 分量進行分解并重構,最后計算重構PR 分量的關聯維數來有效提取齒輪的故障信息.限于篇幅,以第2級平行軸小齒輪點蝕故障振動信號為例,圖11所示為其EITD 分解結果,表2列出了各PR 分量與振動信號的互相關系數.根據最大相關系數準則選取PR2分量作為主PR 分量,并據此進一步采用小波包變換法進行分析與重構.圖12所示為主PR 分量經3層小波包分解后各小波包系數的能量分布情況,其中8個頻帶的能量值分別為55.15、23.31、5.11、13.61、0.03、0.28、1.38和1.14,據此選擇前4個能量比重大的小波包系數重構PR2分量作為降噪處理.

圖10 3種狀態下齒輪振動信號的時域波形Fig.10 Time domain waveform of gear vibration signal under three conditions

圖11 點蝕故障下齒輪振動信號的EITD分解結果Fig.11 EITD result generated from gear vibration signal with existence of pitting fault

表2 各PR 分量與齒輪點蝕振動信號的互相關系數Tab.2 Correlation coefficients between PRCs and the gear pitting fault signal

圖12 點蝕故障下齒輪振動信號的能量分布Fig.12 Energy distribution of gear vibration signal with existence of pitting fault

在計算關聯維數提取特征信息前,首先采用文獻[14]中的C-C法和文獻[15]中的方法(Cao方法)來確定延遲時間和嵌入維數.根據C-C 法的原理,由圖13可知,C-C法的統計量均值S(t)的第1個零點對應的時間為5s,即最佳延遲時間為5s.根據Cao提出的方法[15]來確定最小嵌入維數,由圖14可明顯看出,在嵌入維數為10 時參數E1和E2的值趨近于1,即最小嵌入維數為10.隨后分別計算嵌入維數在10~24 內重構PR2分量的關聯維數.同理按照此過程,依次計算齒輪在正常、第3級平行軸齒輪磨損和斷齒混合故障下重構主PR 分量的關聯維數,計算結果如圖15所示.由圖15可明顯看出,不同工作狀態下的關聯維數不同,在一定嵌入維數(10<m<24)范圍內可以有效分辨出齒輪的各個工作狀態和故障類型,表明EITD 方法與小波包分解降噪相結合計算關聯維數實現了風力發電機組齒輪故障特征與噪聲的分離,且具有較好的識別效果.

圖13 C-C法求延遲時間Fig.13 Time delay determined by correlation integral method

圖14 根據Cao方法求嵌入維數Fig.14 Embedding dimension determined by Cao method

圖15 齒輪振動信號經EITD-WPT處理后在不同嵌入維數下的關聯維數Fig.15 Correlation dimension of gear vibration signal obtained by EITD-WPT under different embedding dimensions

為了驗證本文方法的有效性,首先對3種工作狀態下的齒輪振動信號進行EMD 分解,并直接計算其主IMF分量的關聯維數,如圖16所示.由于強背景噪聲等因素對故障特征信息的干擾與耦合,圖16中3種故障類型出現了相互交叉,很難辨識出各自的工作狀態.圖17所示為3種工作狀態下的齒輪振動信號經EMD 方法和小波包分解聯合處理后,對包含故障特征信息最豐富的IMF 分量進行計算得到的不同嵌入維數下的關聯維數.由圖17可知,雖然嵌入維數大于15時3種工作狀態下的關聯維數有著較大區分間隔,但在嵌入維數為10~14時出現了重疊現象,造成齒輪的工作狀態和故障類型分辨模糊.可見,經小波包的EMD 方法處理后能夠減少一些噪聲的干擾,但在強噪聲干擾下獲得的主IMF分量的關聯維數并不能精確提取齒輪的故障特征,且類型識別不明顯.

圖16 齒輪振動信號經EMD處理后在不同嵌入維數下的關聯維數Fig.16 Correlation dimension of gear vibration signal obtained by EMD under different embedding dimensions

圖17 齒輪振動信號經EMD-WPT處理后在不同嵌入維數下的關聯維數Fig.17 Correlation dimension of gear vibration signal obtained by EMD-WPT under different embedding dimensions

5 結 論

(1)在ITD 方法和EMD 方法的三次樣條插值基礎上提出了EITD 方法.EITD 方法不僅具有端點效應小、分解速度快等優勢,而且避免了PR 分量失真現象.將該方法和小波包變換相結合,實現了風力發電機組齒輪箱故障的精確診斷.

(2)風力發電機組齒輪振動信號受噪聲的影響較大,經EMD 方法分解后直接計算主IMF分量的關聯維數和由EMD 與小波包聯合降噪處理后計算IMF分量的關聯維數都不能很好地提取齒輪故障特征信息并識別工作狀態,而振動信號經EITD 自適應時頻分解及小波包變換對主PR 分量降噪后,再計算PR 分量的關聯維數來提取故障特征,可以實現齒輪狀態和故障類型的有效識別.

[1] 許小剛,王松嶺,劉錦廉.基于小波包能量分析及改進支持向量機的風機機械故障診斷[J].動力工程學報,2013,33(8):606-612.XU Xiaogang,WANG Songling,LIU Jinlian.Mechanical fault diagnosis of fan based on wavelet packet energy analysis and improved support vector machine[J].Journal of Chinese Society of Power Engineering,2013,33(8):606-612.

[2] YANG W,TAVNER P J,WILKINSON M R.Condition monitoring and fault diagnosis of a wind turbine synchronous generator drive train[J].IET Renewable Power Generation,2009,3(1):1-11.

[3] 陳雪峰,李繼猛,程航,等.風力發電機狀態監測和故障診斷技術的研究與進展[J].機械工程學報,2011,47(9):45-52.CHEN Xuefeng,LI Jimeng,CHENG Hang,et al.Research and application of condition monitoring and fault diagnosis technology in wind turbines[J].Journal of Mechanical Engineering,2011,47(9):45-52.

[4] 龍泉,劉永前,楊勇平.基于粒子群優化BP神經網絡的風電機組齒輪箱故障診斷方法[J].太陽能學報,2012,33(1):120-125.LONG Quan,LIU Yongqian,YANG Yongping.Fault diagnosis method of wind turbine gearbox based on BP neural network trained by particle swarm optimization algorithm[J].Acta Energiae Solaris Sinica,2012,33(1):120-125.

[5] 唐貴基,蔡偉.應用小波包和包絡分析的滾動軸承故障診斷[J].振動、測試與診斷,2009,29(2):201-204.TANG Guiji,CAI Wei.Rolling bearings fault diagnosis by using wavelet packet and envelope analysis[J].Journal of Vibration,Measurement &Diagnosis,2009,29(2):201-204.

[6] XIANG Ling,HU Aijun.Comparison of methods for different time-frequency analysis of vibration signal[J].Journal of Software,2012,7(1):68-74.

[7] WICKERHAUSER M V.Acoustic signal compression with wavelet packets[M]//CHUI C K.Wavelets analysis and its applications.Boston,US:Academic Press,1992:670-700.

[8] FREI M G,OSORIO I.Intrinsic time-scale decomposition:time-frequency-energy analysis and realtime filtering of non-stationary signals[J].Proceedings of the Royal Society A,2007,463:321-342.

[9] 安學利,蔣東翔,陳杰,等.基于ITD 和LS-SVM 的風力發電機組軸承故障診斷[J].電力自動化設備,2011,31(9):10-13.AN Xueli,JIANG Dongxiang,CHEN Jie,et al.Bearing fault diagnosis based on ITD and LS-SVM for wind turbine[J].Electric Power Automation Equipment,2011,31(9):10-13.

[10] 楊宇,李杰,潘海洋,等.VPMCD 和改進ITD 的聯合智能診斷方法研究[J].振動工程學報,2013,26(4):608-616.YANG Yu,LI Jie,PAN Haiyang,et al.Research on combined intelligent diagnostic method based on VPMCD and improved ITD[J].Journal of Vibration and Engineering,2013,26(4):608-616.

[11] 段禮祥,張來斌,岳晶晶.基于ITD 和模糊聚類的齒輪箱故障診斷方法[J].中國石油大學學報:自然科學版,2013,37(4):133-139.DUAN Lixiang,ZHANG Laibin,YUE Jingjing.Fault diagnosis method of gearbox based on intrinsic time-scale decomposition and fuzzy clustering[J].Journal of China University of Petroleum:Edition of Natural Science,2013,37(4):133-139.

[12] 王炳成,任朝暉,聞邦椿.基于非線性多參數的旋轉機械故障診斷方法[J].機械工程學報,2012,48(5):63-69.WANG Bingcheng,REN Zhaohui,WEN Bangchun.Fault diagnoses method of rotating machines based on nonlinear multi-parameters[J].Journal of Mechanical Engineering,2012,48(5):63-69.

[13] 高建強,唐樹芳,姜華偉,等.基于關聯維數的鼓泡流化床風帽壓力波動特性研究[J].動力工程學報,2012,32(4):268-272.GAO Jianqiang,TANG Shufang,JIANG Huawei,et al.Characteristic study on pressure fluctuation of wind caps in a bubbling fluidized-bed based on correlation dimension[J].Journal of Chinese Society of Power Engineering,2012,32(4):268-272.

[14] KIM H S,EYKHOLT R,SALAS J D.Nonlinear dynamics,delay times,and embedding windows[J].Physica D:Nonlinear Phenomena,1999,127(1/2):48-60.

[15] CAO Liangyue.Practical method for determining the minimum embedding dimension of a scalar time series[J].Physica D:Nonlinear Phenomena,1997,110(1/2):43-50.

猜你喜歡
故障診斷故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 亚洲乱亚洲乱妇24p| 亚洲成人网在线播放| 日本三区视频| 亚洲综合二区| 国产福利一区视频| 中文字幕资源站| 97国产在线视频| 婷婷色中文| 国产人妖视频一区在线观看| 日韩精品一区二区三区swag| 久久女人网| 国产小视频a在线观看| 在线视频一区二区三区不卡| 一级不卡毛片| 黄色三级毛片网站| 亚洲欧美日韩另类| 亚洲国产精品国自产拍A| 国产微拍一区二区三区四区| 中文字幕无线码一区| 女同久久精品国产99国| 欧美区一区| 午夜视频日本| hezyo加勒比一区二区三区| 国产成人免费高清AⅤ| 少妇精品网站| 无码国内精品人妻少妇蜜桃视频| 国产手机在线小视频免费观看| 国产高清无码第一十页在线观看| 99视频精品全国免费品| 国产精品99久久久久久董美香| 久久精品视频一| 欧美一区二区三区欧美日韩亚洲 | 久久久久久高潮白浆| 国产精品尤物在线| 日本黄色不卡视频| 国产精品护士| 中文字幕在线看视频一区二区三区| 手机在线免费不卡一区二| 老司机精品久久| 2020最新国产精品视频| 中文无码精品A∨在线观看不卡| 野花国产精品入口| 亚洲精品手机在线| 亚洲一区二区三区国产精品 | aa级毛片毛片免费观看久| 91视频日本| 日韩欧美高清视频| 久久精品电影| 91青青草视频在线观看的| 亚洲人成网18禁| 全免费a级毛片免费看不卡| 亚洲中文制服丝袜欧美精品| 国产精品污视频| 色婷婷亚洲十月十月色天| 国产精品内射视频| 色婷婷成人网| 中文字幕无线码一区| 亚洲h视频在线| 日韩无码视频网站| 九九九久久国产精品| 青草视频网站在线观看| 亚洲日本中文字幕乱码中文| 内射人妻无套中出无码| 国产成人精品一区二区三区| 亚洲精品人成网线在线| 日本久久网站| 中文字幕在线永久在线视频2020| 久久人人爽人人爽人人片aV东京热| 国产精品太粉嫩高中在线观看 | 全午夜免费一级毛片| 色综合久久综合网| 四虎国产永久在线观看| 免费国产高清精品一区在线| 性视频久久| 97综合久久| 国产亚洲精| 九九精品在线观看| 亚洲国产成熟视频在线多多 | 日韩高清欧美| 亚洲—日韩aV在线| 国产区人妖精品人妖精品视频| 国产第一福利影院|