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

基于線調頻小波路徑追蹤的滾動軸承故障特征提取

2017-07-18 11:49:23程衛東劉東東趙德尊
振動與沖擊 2017年13期
關鍵詞:故障信號

程衛東, 劉東東, 趙德尊

(北京交通大學 機械與電子控制工程學院,北京 100044)

基于線調頻小波路徑追蹤的滾動軸承故障特征提取

程衛東, 劉東東, 趙德尊

(北京交通大學 機械與電子控制工程學院,北京 100044)

針對齒輪噪源以及變轉速的工作條件雙重干擾下的滾動軸承的故障診斷,提出了一種基于線調頻小波路徑追蹤的滾動軸承故障特征提取方法。對混合信號進行Hilbert變換得到包絡信號,并在滿足采樣定理的條件下對包絡信號進行降采樣;對降采樣包絡信號應用線調頻小波路徑追蹤算法提取軸承的瞬時故障特征頻率趨勢線,再對軸承的瞬時故障特征頻率趨勢線和轉速曲線進行等時間重采樣,并求各時間點的計算瞬時故障特征系數;根據計算瞬時故障特征系數與軸承的外圈、內圈和滾動體的故障特征系數進行比較,完成故障診斷。通過處理仿真信號和實測信號證明了該方法在不對混合信號進行濾波和使用階比分析的情況下,不僅能檢測滾動軸承是否出現故障,而且能確定故障的發生位置。

齒輪噪源;變轉速;軸承故障診斷;線調頻小波;計算瞬時故障特征系數

齒輪噪源干擾和變轉速的工作模式一直是困擾滾動軸承故障診斷的兩大難題[1]。齒輪較大幅值的振動往往會掩蓋軸承的故障沖擊特征,齒輪嚙合振動的基頻、諧頻譜線成分也影響軸承故障引起的共振頻帶的獲取[2]。自適應噪聲消除技術、周期與隨機量分離技術、時域同步平均技術以及適應性噪聲消除技術均曾被用來去除齒輪噪源的干擾[3-4]。文獻[5]在不同的位置安裝兩個傳感器,分別采集噪聲背景下的軸承信號和背景噪聲,將這兩個信號分別作為自適應去噪系統的主輸入和參考輸入使用自適應消噪技術去噪,但是,這種去噪方法需要安裝兩個傳感器,受到安裝成本以及安裝條件的限制。文獻[6]指出周期與隨機量分離算法要求保證軸承的轉速不存在波動。張西寧等[7]運用同步平均的方法對信號進行處理,當軸承轉速發生變化時,同步平均不能完全去除齒輪嚙合、不對中產生的周期性噪聲成分[8]。通過以上分析發現,現有的對于齒輪噪源影響的軸承故障的診斷基本都是先使用濾波技術去掉齒輪噪源的影響,再對剩余信號進行分析。濾波算法不僅復雜,效率較低,而且還可能需要安裝額外的輔助設備,使設備檢測成本提高。

目前,對于變轉速工況下的滾動軸承的故障診斷常采用階比跟蹤技術[9-11],將時域信號轉化為角域信號,再進行處理。常用的階比分析技術有硬件階比跟蹤和計算階比跟蹤兩種。與硬件階比跟蹤技術相比計算階比跟蹤不需要采樣率合成器、抗混疊跟蹤濾波器、轉速計等硬件設備[12],因而引起了學者的廣泛研究。然而,計算階比分析也存在諸多弊端。例如,Saavedra等[13]證明計算階比跟蹤算法中插值算法存在不可避免的誤差。等角度重采樣時從時域到角域轉變會涉及到大量高次方程,計算效率低,甚至出現無解情況,同樣會存在誤差。另外,階比跟蹤技術還存在包絡畸變、幅值誤差等問題。

近年,Candès等[14]提出了線調頻小波路徑追蹤算法,該算法與常用提取故障特征頻率的峰值搜索算法相比,不僅能更精確的提取信號各頻率成分變化趨勢線,而且還能根據需要選擇提取不同頻率范圍內的頻率成分變化趨勢線。本文利用這一特性并結合分析信號在時頻域的特點提出了不使用濾波技術和階比跟蹤技術的基于線調頻小波路徑追蹤的滾動軸承故障特征提取方法。

具體方法:首先,對原始振動信號進行Hilbert變換得到包絡信號,并在滿足采樣定理的條件下對包絡信號進行降采樣;然后,對降采樣包絡信號應用線調頻小波路徑追蹤算法提取故障軸承的瞬時故障特征頻率(Instantaneous Fault Characteristic Frequency,IFCF)趨勢線,再對軸承的IFCF趨勢線和轉速曲線進行等時間重采樣,并求各時間點的計算瞬時故障特征系數;最后,根據計算瞬時故障特征系數與軸承的外圈、內圈和滾動體的故障特征系數進行比較,完成故障診斷。通過處理仿真信號和實驗信號來證明該方法在不對原始信號進行濾波和使用階比分析的情況下,不僅能檢測滾動軸承是否出現故障,而且能確定故障發生位置。

1 滾動軸承的故障特征頻率及特征系數

當滾動軸承的外圈、內圈或者滾動體出現凹坑或裂紋時,故障點與其接觸點會產生沖擊,由于系統的衰減特性,形成高幅值而快速衰減的沖擊響應。根據軸承結構,尺寸參數以及轉頻可以計算軸承的各部分的故障特征頻率FCF(Fault Characteristic Frequency)[15-16]

外圈故障特征頻率fo

(1)

內圈故障特征頻率fi

(2)

滾動體故障特征頻率fb

(3)

式中:N為滾動體個數;fr為軸承內圈轉頻;d為滾動體直徑;D為節圓直徑;φ為接觸角。

根據各部分的故障特征頻率便可以得到故障特征系數(Fault Characteristic Coefficient, FCC)

(4)

(5)

(6)

由式(4)、式(5)和式(6)可以看出FCC的大小只與滾動軸承的參數有關,不會隨著轉速以及其他因素的變化而變化,根據這個特性,在滾動軸承的故障診斷的過程中我們只需要求出FCC并與使用軸承各部分的FCC比較,不僅能判斷出軸承是否發生故障還能確定故障位置。由于FCF以及轉速曲線的獲取都涉及插值估計等函數,因此,求出的每個時間點的FCC不可能為定值,應該與實際FCC有一定誤差。本文將使用線調頻小波路徑追蹤獲取的瞬時故障特征頻率與對應點轉速的比值稱為計算瞬時故障特征系數(Calculating Instantaneous Fault Characteristic Coefficient, CIFCC),進而通過分析求得的CIFCC來判斷軸承故障的情況。

2 線調頻小波路徑追蹤算法

線調頻小波路徑追蹤原理是將信號時間長度以二進制劃分為動態的時間支撐區,從定義的chirplet原子庫中尋找各動態支撐區相關性最大的原子并采用最佳路徑連接的算法逐一連接,自適應的尋找信號中相關性最大的信號分量的瞬時頻率變化趨勢線。具體算法如下:

建立chirplet原子庫

(7)

式中:I為動態支撐區間,I=[kN2-j,(k+1)N2-j]×Δt;k為動態區間的序號,k=0,1,…,2j-1;N為待分析信號的采樣點數;j動態支撐區間的尺度系數,j=0,1,…,log2N-5;1I(t)為矩形窗函數,即t∈I時,1I(t)=1,t?I時,1I(t)=0;aμ為調頻率系數;bμ為頻率偏置系數;根據采樣定理,aμt+bμ為動態支撐區間內瞬時頻率,應小于fs/2。

定義分析信號為信號與原子庫中所有原子的內積,即

(8)

式中:Δt為采樣時間間隔;n為t時刻采樣點數。

在建立的原子庫中尋找一組原子連接,使其既能覆蓋整個分析信號的長度,又能使分析信號在時間支撐區域上有最大的投影系數,支撐區間內最大投影系數為

(9)

定義最大投影系數kI在時間支撐區I內代表的信號分量為dI(t),得

(10)

連接最大投影系數kI使得信號分量dI(t)信號具有最大的總能量,即

(11)

式中,Π應覆蓋整個時間長度且不重疊。

由此獲得最佳路徑連接算法:

(1) 初始化 取e(i)=0,pred(i)=0,其中,e(i)為選取i-1個動態時間支撐區的最大投影系數對應的總能量和,pred(i)為連接到第i個動態時間支撐區的前置支撐區的序號。

e(i)+d(i)>e(k)

(12)

則令

(13)

該連接方法能確保基原函數的組合與待分析信號的相似度最高,而動態支撐區Π內的瞬時頻率為aμt+bμ,線性連接各個動態支撐區的瞬時頻率就獲得了待分析信號分解信號的瞬時頻率。

3 基于線調頻小波路徑追蹤的滾動軸承故障特征提取方法

從第2節的線調頻小波路徑追蹤的算法可以看出,線調頻小波路徑追蹤算法與峰值搜索算法不同,應用線調頻小波路徑追蹤算法可以提取選定頻率范圍內的復雜信號中的頻率成分的變化趨勢線。齒輪的嚙合頻率為齒輪所在軸的整齒數倍,為了防止齒輪加工的過程中出現“根切現象”,齒輪的齒數一般大于等于17。而如式(4)~式(6)所示,軸承的故障特征頻率與轉頻的比值為故障特征系數,故障特征系數一般很小,所以,齒輪嚙合頻率及其倍頻不會影響故障特征頻率的提取。正常齒輪的振動無調幅調頻現象,不會出現嚙合頻率被轉頻及其高倍頻調制的現象,對齒輪振動信號進行Hilbert變換產生的包絡信號不會在低頻部分產生轉頻及其倍頻[17]。因此,齒輪的嚙合頻率與倍頻及其包絡信號均不會對軸承低頻的包絡信號中的故障特征頻率的提取產生干擾。鑒于此,提出了基于線調頻小波路徑追蹤的滾動軸承的故障特征提取方法。

圖1給出了具體算法如下:

(1) 對原始振動信號進行Hilbert變換得到包絡信號,并在滿足采樣定理的條件下對包絡信號進行降采樣;

(2) 根據實驗軸承的各部分的FCC和轉速曲線計算FCF的變化范圍,并應用線調頻小波路徑追蹤算法對降采樣包絡信號提取故障軸承的IFCF趨勢線;

(3) 為了使IFCF趨勢線的點與轉速曲線點在時間上一一對應,對IFCF趨勢線和轉速曲線進行等時間重采樣(這里的重采樣區別于階比分析需要的等角度重采樣,不涉及信號的時域到角域的轉換,只是用到簡單的一次朗格朗日插值函數,而且一般重采樣100個點就能滿足需要,避免了引言所說的等角度重采樣需要解大量高次方程等缺陷);

(4) 根據重采樣的IFCF趨勢線和轉速曲線求各時間點的計算瞬時故障特征系數;

(5) 根據計算瞬時故障特征系數的分布范圍與軸承的外圈、內圈和滾動體的FCC進行比較,判斷軸承是否發生故障。

圖1 基于線調頻小波路徑追蹤的滾動軸承故障診斷方法流程圖

Fig.1 Flowchart of rolling bearing fault diagnosis based on chirplet path pursuit

4 仿真信號分析

為了驗證提出診斷方法的有效性,對滾動軸承的仿真信號進行分析。模擬滾動軸承的外圈出現故障,故障特征系數C0=3,滾動軸承內圈的轉速如圖5(實線)所示隨時間成線性變化為v(t)=10t+5(rps),軸承轉過角度與時間之間的關系

(14)

式中,θi=i/C0(i=1,2,…,k),由θi=i/C0和式(14)可以求得每個沖擊對應的時標Ti。

該混合信號主要由故障軸承的沖擊成分xbearing、齒輪嚙合振動產生信號成分xgear和白噪聲成分n(t)三部分構成,如式(15)所示。

x(t)=xbearing+xgear+n(t)

(15)

本文改進了勻轉速工況下滾動軸承的振動信號的仿真模型[18],構造變轉速工況下故障滾動軸承振動信號的仿真模型

φω)×u(t-iT-τi)

(16)

齒輪振動信號的仿真模型

φm)

(17)

式中:Ai與Am分別為軸承和齒輪信號第i和m個沖擊時的最大幅值,簡化該幅值與時間成線性關系;fr為共振頻率;fz(t)為齒輪轉頻隨時間變化規律;參考軸承轉速的時間浮動,定義滑移常數λ;軸承和齒輪信號的初相角φω、φm均設定為0;齒輪箱主動輪和外圈故障軸承同軸,齒數為Z,單級傳動,其參數設定如下表所示。

根據以上仿真模型和參數得到如圖2所示的時域波形圖。對仿真信號做短時傅里葉變換(Short-time Fourier Transform,STFT)得到如圖3所示的時頻圖,從時頻圖可以看到齒輪嚙合頻率及其倍頻與軸承共振頻率產生交叉,不能通過帶通濾波直接獲得共振頻率。圖4給出了譜峭度圖,從譜峭度圖中看出最大峭度值出現在尺度為7,中心頻率為119 Hz處,而不是在中心頻率為2 000 Hz(共振頻率)處,這也直觀反應了利用譜峭度濾波的限制。分析得知,通過濾波處理以獲得共振頻率有一定難度。因此,這里試圖避開使用濾波處理仿真信號。

表1 振動信號仿真模型參數

根據軸承轉速變化范圍5~25 Hz(如圖5所示)和C0=3,求得故障特征頻率的變化范圍為15~75 Hz,為了提高線調頻小波路徑追蹤算法處理的效率,對經過Hilbert變換得到的包絡信號進行降采樣,降采樣后的采樣頻率為1 024 Hz,通過計算,降采樣后的采樣頻率仍滿足采樣定理。設置頻率提取范圍為10~80 Hz,應用線調頻小波追蹤算法得到如圖5所示的IFCF趨勢線。對IFCF趨勢線與轉速曲線進行重采樣處理(取100個點)并求各時間點的CIFCC,做如圖6所示的CIFCC分布范圍的柱狀圖,81%的點數落在了3±0.015范圍內并且100%的點落在了3±0.045范圍內,恰好對應設置的外圈的FCC(C0=3)。因為同一軸承的外圈、內圈以及滾動體的FCC不同,所以可以判斷滾動軸承是否發生故障以及故障發生的位置。對仿真信號的分析驗證了本文提出方法的有效性。

圖2 仿真信號時域波形圖

圖3 仿真信號時頻譜

Fig.3 The time-frequency representation of simulated signal

圖4 仿真信號譜峭度圖

圖5 仿真信號轉頻與IFCF趨勢線

圖6 仿真信號的計算瞬時故障特征系數分布圖

5 實測信號驗證

利用滾動軸承實驗臺測取的外圈故障滾動軸承和內圈故障滾動軸承的振動信號對本算法進一步驗證。其中,與軸承同軸的齒輪齒數為15,實驗軸承的參數如表2所示。

表2 滾動軸承結構參數

采集系統的采樣率為24 000 Hz,截取時間段為2 s的數據,轉速如圖10所示,變化范圍41~94 r/min。圖7和圖8給出了實驗信號的時域波形圖和時頻圖,從時頻圖可以明顯的看到齒輪嚙合頻率及其倍頻與軸承共振頻率交叉,不能直接通過帶通濾波直接提取共振頻率。圖9給出了仿真信號的譜峭度圖,從譜峭度圖中看出最大峭度值出現在尺度7,中心頻率為2 390 Hz處,而不是在中心頻率為6 000 Hz(共振頻率)處,這也直觀反應了利用譜峭度濾波的限制。

對經Hilbert變換得到包絡信號進行降采樣,降采樣后的采樣頻率為2 048 Hz。設置頻率提取范圍為220~510 Hz,應用線調頻小波追蹤算法得到如圖10所示的IFCF趨勢線。對IFCF趨勢線與轉速曲線進行重采樣處理(取100個點)并求各時間點的CIFCC,做如圖11所示的CIFCC分布范圍的柱狀圖,81%的點數落在了3.58±0.02范圍內并且100%的點落在了3.58±0.06范圍內,恰好對應軸承的外圈的故障特征系數C0=3.58,由此可以判斷軸承的外圈出現故障。

通過處理內圈故障滾動軸承振動信號進一步驗證本文算法的有效性,圖12為內圈故障軸承振動信號的時域波形圖,利用本文算法對該振動信號進行處理,得到如圖13所示的CIFCC分布柱狀圖,觀察發現,100%的點落在了5.42±0.04范圍內,恰好對應軸承的內圈的故障特征系數C0=5.42,由此可以判斷軸承的內圈出現故障。

圖7 外圈故障滾動軸承信號時域波形圖

圖8 外圈故障滾動軸承信號時頻譜

Fig.8 The time-frequency representation of outer faulty bearing signal

圖9 外圈故障滾動軸承信號譜峭度圖

圖10 外圈故障軸承信號轉頻與IFCF趨勢線

Fig.10 The curve of outer faulty bearing signal rotating frequency and FCF

圖11 外圈故障軸承信號的計算瞬時故障特征系數分布圖

圖12 內圈故障滾動軸承信號時域波形圖

圖13 內圈故障軸承信號的計算瞬時故障特征系數分布圖

6 結 論

對混合信號經過Hilbert變換和降采樣處理后得到降采樣包絡信號,然后利用線調頻小波路徑追蹤算法可以從降采樣包絡信號中提取軸承的IFCF趨勢線。根據IFCF趨勢線與轉速曲線可以得到CIFCC,通過分析CIFCC的大小分布能夠判斷軸承發生故障的位置。整個處理過程沒有使用濾波算法和階比跟蹤技術,避免了濾波算法使用的局限性和因安裝輔助設備產生的成本的提高以及階比跟蹤技術的不足,為齒輪噪源影響下的變轉速工況下滾動軸承的故障診斷提供了一種新思路。

[1] ZHAO Dezun, LI Jianyong, CHENG Weidong. Feature extraction of faulty rolling element bearing under variable rotational speed and gear interferences conditions[J]. Shock and Vibration, 2015(3): 1-9.

[2] 王天楊,李建勇,程衛東. 基于改進的自適應噪聲消除和故障特征階比譜的齒輪噪源干擾下變轉速滾動軸承故障診斷[J]. 振動與沖擊, 2014, 33(18): 7-13.

WANG Tianyang,LI Jianyong,CHENG Weidong. Fault diagnosis of rolling bearing under a variable rotational speed and gear vibration noise based on revised ANC algorithm and FCO spectrum[J]. Journal of Vibration and Shock, 2014, 33(18): 7-13.

[3] ANTONI J, RANDALL R B. Unsupervised noise cancellation for vibration signals:part II-a novel frequency-domain algorithm[J]. Mechanical Systems and Signal Processing, 2004, 18(1): 103-117.

[4] ANTONI J, RANDALL R B. Unsupervised noise cancellation for vibration signals:Part I-evaluation of adaptive algorithms[J]. Mechanical Systems and Signal Processing, 2004, 18(1): 89-101.

[5] CHATURVEDI G K, THOMAS D W. Bearing fault detection using adaptive noise canceling[J]. Journal of Mechanical Design, 1982, 104(2): 280-289.

[6] RANDALL R B, SAWALHI N. A comparison of methods for separation of deterministic and random signals[J].Condition Monitoring, 2011,1(1):11-43.

[7] 張西寧,吳婷婷,徐進杰. 變轉速齒輪箱振動信號監測的無鍵相時域同步平均方法[J].西安交通大學學報, 2012, 46(6): 111-116.

ZHANG Xining, WU Tingting, XU Jinjie,et al. Time-domain synchronous average method without key-phase signal for vibration monitoring of variable speed gearbox[J]. Journal of Xi’an Jiaotong University, 2012, 46(6): 111-116.

[8] BORGHESANI P, RANDALL R B. Application of cepstrum pre-whitening for the diagnosis of bearing faults under variable speed conditions[J]. Mechanical Systems and Signal Processing, 2013, 36(2): 370-384.

[9] BONNARDOT F, BADAOUI M E, RANDALL R B, et al. Use of the acceleration signal of a gearbox in order to perform angular resampling(with limited speed fluctuation)[J]. Mechanical Systems and Signal Processing, 2005, 19(4): 766-785.

[10] 郭瑜,秦樹人. 無轉速計旋轉機械升降速振動信號零相位階比跟蹤濾波[J]. 機械工程學報,2004, 40(3): 50-54.

GUO Yu, QIN Shuren. Tacholess order tracking filtering for run-up or coast down vibration signal of rotating machinery based on zero-phase distortion digital filtering[J]. Chinese Journal of Mechanical Engineering, 2004, 40(3): 50-54.

[11] 郭瑜,秦樹人. 旋轉機械非穩定信號的偽轉速跟蹤階比分析[J].振動與沖擊,2004, 23(1): 61-69.

GUO Yu, QIN Shuren. Pseudo-speed tracking order analysis for non-stationary[J]. Journal of Vibration and Shock, 2004, 23(1): 61-69.

[12] CHENG Weidong, GAO R X, WANG Jinjiang. Envelope deformation in computed order tracking and error in order analysis[J]. Mechanical Systems and Signal Processing, 2014, 48(1/2): 92-102.

[13] SAAVEDRA P N,RODRIGUEZ C G. Accurate assessment of computed order tracking[J]. Shock and Vibration, 2006, 13(1): 13-21.

[15] RANDALL R B, ANTONI J. Rolling element bearing diagnostics-a tutorial[J]. Mechanical Systems and Signal Processing, 2011, 25(2): 485-520.

[16] ANTONI J, RANDALL R B. The spectral kurtosis: application to the vibratory surveillance and diagnostics of rolling machines[J]. Mechanical Systems and Signal Processing, 2006, 20(2): 308-331.

[17] 彭富強,于德介,武春燕. 基于多尺度線調頻基稀疏信號分解的包絡解調方法及其在齒輪故障診斷中的應用[J]. 機械工程學報, 2010, 46(12): 1-12.

PENG Fuqiang, YU Dejie, WU Chunyan. AM-FM signal extraction method by the sparse signal decomposition based on multi-scale chirplet and its application to gear fault diagnosis[J]. Journal of Mechanical Engneering, 2010, 46(12): 1-12.

[18] LIANG M, BOZCHALOOI I S. An energy operator approach to joint application of amplitude and frequencydemodulations for bearing fault detection[J]. Mechanical Systems and Signal Processing, 2010, 24(5): 1473-1494.

Fault feature extraction method for rolling bearings based on chirplet path tracing

CHENG Weidong, LIU Dongdong, ZHAO Dezun

(School of Mechanical Electronic and Control Engineering, Beijing Jiaotong University, Beijing 100044, China)

Aiming at fault diagnosis of rolling element bearings under disturbances of variable rotating speed and gear vibration noise, a fault feature extraction method for rolling bearings was proposed based on the chirplet path tracing method here. Firstly, the original vibration signals were transformed with Hilbert transformation to get envelope signals, the envelope ones were down sampled under the condition of satisfying the sampling theorem. Then, the down sampled envelope signals were used to extract the instantaneous fault characteristic frequency (IFCF) trend curve of the faulty bearing with the chirplet path tracing method. The instantaneous fault characteristic frequency (IFCF) trend curve and the rotating speed curve of the bearing were resampled with equal time intervals and the calculating instantaneous fault characteristic coefficient (CIFCC) at each time point was gained. Finally, the fault diagnosis was completed by comparing CIFCCs with the fault characteristic coefficients (FCCs) of the outer ring, inner ring and rolling element. Through processing simulated signals and test ones without filtering and order analysis, it was shown that the method can not only detect if there is any fault in rolling bearings, but also determine the fault’s location.

gear vibration noise; variable rotating speed; fault diagnosis; chirplet; calculating instantaneous fault characteristic coefficient(CIFCC)

國家自然科學基金資助項目(51275030)

2016-03-21 修改稿收到日期:2016-05-23

程衛東 男,博士,副教授,1967年生

TH212;TH213.3

A

10.13465/j.cnki.jvs.2017.13.024

猜你喜歡
故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
孩子停止長個的信號
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
故障一點通
故障一點通
故障一點通
主站蜘蛛池模板: 亚洲精品桃花岛av在线| 欧美第二区| 亚洲精品国产首次亮相| 新SSS无码手机在线观看| 伊人国产无码高清视频| 99热国产在线精品99| 美女亚洲一区| 国产精品久久久久无码网站| 国产日本视频91| 国产精品成人啪精品视频| 欧美国产精品不卡在线观看| 亚洲系列无码专区偷窥无码| 亚洲香蕉伊综合在人在线| 亚洲色图欧美视频| 国产国产人成免费视频77777 | 国产精品成| 一级成人a毛片免费播放| 亚洲综合18p| 99激情网| 日韩色图在线观看| 国产精品精品视频| 国产乱码精品一区二区三区中文| 久久精品中文字幕少妇| 亚洲VA中文字幕| 国产精品男人的天堂| 免费人成黄页在线观看国产| 亚洲无码37.| 伊人久久大香线蕉成人综合网| 成年看免费观看视频拍拍| 欧美中文一区| 亚洲国产综合精品一区| 国产新AV天堂| 日本三区视频| 国产在线精品美女观看| 国产精品永久久久久| 国产精欧美一区二区三区| 青青草久久伊人| 中文字幕无码av专区久久 | 国产精品9| 亚洲欧美日韩视频一区| 丁香婷婷在线视频| 国产资源免费观看| 国产极品粉嫩小泬免费看| 国产成人免费手机在线观看视频| 国产综合在线观看视频| av免费在线观看美女叉开腿| 91成人在线免费观看| 97免费在线观看视频| 亚洲成人在线网| 国产呦视频免费视频在线观看| 波多野结衣AV无码久久一区| 最新精品久久精品| 美女啪啪无遮挡| 毛片视频网址| 88av在线播放| 午夜国产不卡在线观看视频| 日韩小视频在线观看| 国产色婷婷| 国产AV无码专区亚洲A∨毛片| 日韩av在线直播| 国产精品福利尤物youwu| 久久国产黑丝袜视频| 青青操国产| 亚洲一区二区三区中文字幕5566| 午夜不卡福利| 亚洲不卡影院| 亚洲无码视频一区二区三区| 久久夜色撩人精品国产| 亚洲欧美成人在线视频| 中国国语毛片免费观看视频| 亚洲国产看片基地久久1024| 亚洲 欧美 偷自乱 图片 | 丁香六月激情综合| 2018日日摸夜夜添狠狠躁| 中文字幕1区2区| 99久久无色码中文字幕| 亚洲成A人V欧美综合| 欧美曰批视频免费播放免费| 欧美黄网在线| 制服丝袜在线视频香蕉| 大香伊人久久| 99re在线视频观看|