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

基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法

2015-10-29 02:35:55唐貴基王曉龍
中國機械工程 2015年11期
關鍵詞:故障診斷故障信號

唐貴基 王曉龍

華北電力大學,保定,071000

基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法

唐貴基王曉龍

華北電力大學,保定,071000

滾動軸承處于早期故障階段時,特征信號微弱,并且受環境噪聲影響嚴重,因此故障特征提取困難。針對這一問題,將最大相關峭度解卷積算法應用于軸承故障診斷,并通過包絡譜稀疏度來篩選最佳解卷積周期參數,提出了基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法。利用最佳參數相對應的最大相關峭度解卷積算法對原信號進行處理,得到解卷積信號后計算其包絡譜,通過分析包絡譜中幅值突出的頻率成分來判斷故障類型。早期故障仿真信號及實測全壽命數據分析結果表明,該方法可有效應用于軸承早期故障診斷。

滾動軸承;稀疏度;最大相關峭度解卷積;故障診斷

0 引言

滾動軸承是旋轉機械中的重要組成零部件,其狀態正常與否決定了整個系統的性能好壞。當軸承出現局部損傷或缺陷時,輕則使設備產生噪聲、振動異常,重則損壞設備,因此研究滾動軸承的狀態監測及故障診斷方法具有重要意義[1]。

當軸承處于早期故障階段時,由故障產生的沖擊成分通常比較微弱,并淹沒于強烈的背景噪聲中,因此想要準確提取出故障特征頻率信息相對困難,探尋軸承的早期故障診斷方法也一直是故障診斷領域的熱點和難點。文獻[2]利用小波相關濾波的降噪特性,將相關濾波降噪方法和包絡譜相結合,提出了基于小波相關濾波包絡分析的軸承早期故障特征提取方法;為確定由故障引起的共振調制邊頻帶,文獻[3]提出了共振解調結合小波包系數熵閾值降噪的綜合算法;為準確判斷故障類型,文獻[4]將雙重Q因子分析方法用于軸承早期故障診斷中,原信號被分解成高共振成分和低共振成分,并利用低Q因子來提取故障沖擊成分;文獻[5]采用基于互相關系數和峭度準則的經驗模態分解(EMD)降噪方法對信號進行預處理,再利用譜峭度選取最佳濾波參數,最后使用帶通濾波和包絡解調進行故障診斷。上述各種方法在軸承早期故障診斷應用中均取得了一定的效果。

1977年,Wiggins將最小熵概念引入到盲解卷積問題處理中,提出了最小熵解卷積方法,并被廣泛應用到地震波處理、超聲檢測及機械故障診斷等領域[6-7]。最近,McDonald等[8]提出了一種新的解卷積方法——最大相關峭度解卷積(MCKD),該方法以相關峭度最大化為目標,充分考慮了信號中沖擊成分的周期特性,旨在突出運算結果中的連續脈沖序列,本文嘗試將其引入到滾動軸承故障診斷領域。然而MCKD方法的運算結果受解卷積周期參數T影響嚴重,為避免實際故障信號處理過程中該參數盲目選擇對解卷積效果的影響,筆者以包絡譜稀疏度作為指導標準來選擇最佳的解卷積周期參數,并提出了基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法,模擬早期故障信號及全壽命周期加速試驗信號均驗證了該方法的正確性。

1 基本原理介紹[8]

1.1相關峭度和稀疏度

零均值信號y(n)(n=1,2,…,N)關于周期參數T的相關峭度表達式為

(1)

相關峭度是在峭度基礎上提出的概念,當T=0時,相關峭度即退化為峭度,該指標由于考慮了沖擊成分的周期特性,因此適合衡量信號中具有特定周期的脈沖序列所占的比重,與峭度相比,相關峭度更強調沖擊成分的連續性。

稀疏度是反應信號稀疏特性的統計參量,信號y(n)的稀疏度表達式為[9]

(2)

當信號的幅值分布相對均勻、差別不大時,信號的稀疏性較弱,稀疏度指標較小,而當信號中出現少量尖脈沖、個別位置幅值突出時,信號則呈現出較強的稀疏特性,稀疏度指標也隨之增大。

1.2最大相關峭度解卷積

假設yn為一個沖擊信號,hn為沖擊信號yn通過周圍環境及路徑時的傳輸衰減響應,xn為實際采集到的信號,以上過程可用公式表示為

xn=hnyn+en

(3)

為便于分析,在此先不考慮噪聲en的影響,則MCKD算法的本質是尋找一個有限脈沖響應(FIR)濾波器,通過輸出信號xn恢復輸入信號yn,即

(4)

其中,[f1f2…fL]T=f是長度為L的濾波器系數。

為使解卷積結果突出連續尖脈沖,該算法以信號的相關峭度KC為評定標準,并將其作為目標函數以求得最優結果,即

(5)

上述優化求解問題等價于求解方程:

(6)

由式(4)和式(5)可以得到如下表達式:

(7)

以矩陣的形式將式(7)重新表述為

4‖y‖-6‖β‖2X0y=0

(8)

β=[y1y1-Ty2y2-T…yNyN-T]T

整理式(8)可得

2‖β‖2X0y=‖y‖2(X0α0+XTα1)

(9)

由于已知:

(10)

則最終的濾波器系數可通過下式得到

(11)

1.3MCKD的數值實現過程

通過上述分析可歸納出最大相關峭度解卷積算法的具體數值實現步驟:

(1)首先確定一個解卷積周期參數T;

(3)設定好FIR濾波器的長度L,并假設初始濾波器系數f=[00…1-1…00]T;

(4)利用式(10)計算濾波后的信號y;

(5)通過y分別計算α0、α1和β;

(6)根據式(11)得到新的濾波器系數f;

(7)比較每次迭代過程的相關峭度變化量ΔKC(T)與閥值ε的大小,如果ΔKC(T)>ε,則返回步驟(4),否則結束整個循環,其中閥值ε是一個較小正數,用于控制迭代次數(本文設定的閥值ε=0.01)。

(8)最終的解卷積信號可由式(10)得到。

2 故障診斷流程

滾動軸承處于故障早期階段時,沖擊特征微弱,而最大相關峭度解卷積可實現信號的解卷積運算,突出信號中的連續脈沖,因此適合處理軸承早期故障信號。筆者研究后發現,在分析實際故障信號時,MCKD算法的處理結果受解卷積周期參數T影響嚴重,理論上,如果參數T與故障沖擊的間隔點相等,則可獲得最好的解卷積效果,然而實際故障信號的脈沖間隔是未知的,因此選擇最佳的解卷積周期參數、有效突出信號中隱含的周期性脈沖序列是利用MCKD算法進行軸承故障診斷的關鍵。

稀疏度能夠有效反應信號的稀疏特性,可用于評價故障信號解卷積處理的效果,然而時域信號的稀疏度容易受單個或少量大幅值脈沖的影響,因此直接以時域信號的稀疏度作為評價標準并不合適,而將信號轉換到頻域,計算包絡譜的稀疏度則可有效避免這一缺陷。在分析滾動軸承早期故障信號時,經過MCKD算法處理后的信號如果周期性沖擊特性不明顯,則解卷積信號包絡譜各頻率處的幅值相差不大,沒有峰值突出的成分,包絡譜稀疏度相對較小;如果信號中出現較明顯的周期性連續脈沖,則包絡譜的相應頻率處會出現較大譜峰,稀疏度也將隨之增大。鑒于上述分析,本文以包絡譜稀疏度為指導標準來搜尋最佳的解卷積周期參數,并提出基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法,具體的診斷流程如下:

(1)計算滾動軸承各部件理論故障特征頻率。根據軸承結構參數分別計算軸承內圈、外圈、滾動體和保持架的故障特征頻率fi、fo、fb和fc。

(2)設定最佳解卷積周期參數的搜尋中心。為加快參數的搜索速率,減少不必要的運算,設定Ti、To、Tb和Tc為4個搜尋中心點,其中Ti=|fs/fi|,To=|fs/fo|,Tb=|fs/fb|,Tc=|fs/fc|,fs為采樣頻率,|·|表示四舍五入取整。

(3)確定搜索范圍。由于實際故障特征頻率與理論計算值之間存在一定差別,如果僅在4個搜尋中心點處尋找最佳解卷積周期參數,很可能造成遺漏或產生偏差,而以中心點為中心確定一個搜索范圍,在搜索范圍內做進一步小范圍搜索則可有效避免此類現象發生。從計算成本角度考慮,為加快運算速率本文選擇在中心點左右5個點范圍內搜索,確定了[Ti-5,Ti+5]、[To-5,To+5]、[Tb-5,Tb+5]和[Tc-5,Tc+5]4個搜索范圍(每個范圍內有11個連續的數值點),在處理不同信號時,搜索范圍大小可根據具體情況作適當縮放調整。

(4)篩選最佳參數。以搜索范圍內的數值點為解卷積周期參數,按順序變化數值點,分別利用MCKD算法對故障信號進行處理,并計算所獲信號的包絡譜稀疏度,得到4條解卷積周期-包絡譜稀疏度曲線。MCKD算法處理原信號時解卷積周期參數越合理,所獲解卷積信號中包含的故障特征成分越多,解卷積信號包絡譜中特征頻率及其倍頻處譜峰越突出,包絡譜呈現出越強的稀疏特性,稀疏度值也相對較大,因此包絡譜稀疏度指標可有效評價解卷積周期參數設置的合理性,選擇曲線中稀疏度最大值點對應的參數即為最佳解卷積周期參數Tm。

(5)軸承故障診斷。設定MCKD算法的解卷積周期為Tm,對原故障信號進行處理,并計算解卷積信號的包絡譜,將包絡譜中幅值明顯的頻率成分與理論故障特征頻率進行對比,從而判斷滾動軸承的故障類型。

3 早期故障模擬信號分析

為驗證基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法的可靠性,本文通過向實測振動信號中添加高斯白噪聲來模擬早期故障信號。首先在千鵬故障(QPZZ)試驗平臺上利用內置電路(ICP)壓電加速度傳感器采集軸承內圈故障振動信號,采樣頻率為12 800 Hz,分析點數為8192,傳動軸轉速為1440 r/min,滾動軸承各部件的理論故障特征頻率如表1所示。

大部分文獻利用模擬信號驗證所提出的軸承早期故障診斷方法的有效性時,均采用故障模型來模擬早期故障信號,為更接近實際工況,本文直接向實測信號中添加較重噪聲來模擬軸承早期故障信號,早期故障信號信噪比rSN為-12 dB,計算公式為

rSN=20lg(ν/νn)

(12)

式中,ν、νn分別為實測振動信號和噪聲的有效值。

振動信號和模擬早期故障信號如圖1所示。

(a)振動信號

(b)模擬早期故障信號圖1 振動信號和模擬早期故障信號

通過對比可發現,添加白噪聲后,實測振動信號中明顯的規律性沖擊完全被噪聲所掩蓋,圖2所示是對模擬早期故障信號直接作包絡譜分析得到的結果,包絡譜中沒有出現峰值明顯的頻率成分。

圖2 模擬早期故障信號的包絡譜

利用本文提出的方法對信號進行分析,結果如圖3所示。首先需確定最佳解卷積周期參數的搜索范圍,根據采樣頻率及各部件故障特征頻率計算得到[70,80]、[105,115]、[208,218]和[1315,1325]4個搜尋范圍。以這4個范圍內的數值點為MCKD算法的解卷積周期參數,依照順序不斷變動參數,分別對模擬早期故障信號作最大相關峭度解卷積運算,并計算獲得信號的包絡譜稀疏度,得到圖3a所示的4條解卷積周期-包絡譜稀疏度曲線(T-S曲線)。分析這些曲線后發現最大稀疏度指標對應的解卷積周期為75,由此確定最佳解卷積周期參數Tm=75,設定MCKD算法的解卷積周期參數為75并對信號進行處理,結果如圖3b所示。與原信號對比后發現,解卷積信號中沖擊成分明顯增多,并呈現出一定的規律性。圖3c為解卷積信號的包絡譜,圖中內圈故障特征頻率fi及其倍頻、轉頻fr及其倍頻、特征頻率的轉頻調制邊帶等成分處譜線峰值明顯,由此可斷定,軸承內圈存在缺陷,診斷結果與實際情況相符。

(a)T-S曲線

(b)最大相關峭度解卷積信號

(c)解卷積信號的包絡譜圖3 本文方法的模擬信號分析結果

仿真信號分析表明,當信號中與故障相關的沖擊特征被較強背景噪聲淹沒時,如果直接對其做包絡譜分析,可能無法有效提取出特征頻率信息,而利用本文提出方法對原信號進行分析后,能夠獲得較理想的處理效果,實現軸承故障類型的準確判斷。

4 應用實例分析

實例分析數據來自于NSFI/UCR智能維護系統中心的滾動軸承全壽命周期加速試驗數據[10],試驗臺轉軸上同時安裝了4個軸承,轉速恒定在2000 r/min,約26 671 N的徑向載荷通過彈性系統加載到軸和軸承上,每個軸承的軸向和徑向各安裝了一個PCB公司生產的353B33型高靈敏度ICP加速度傳感器,圖4給出了軸承和傳感器的安裝位置,表2為試驗軸承各部件的理論故障特征頻率。

圖4 試驗臺示意圖

Hz

利用NI DAQCard-6062E采集卡采集振動信號,采樣頻率為20 kHz,共進行了3組全壽命加速試驗,其中第2組試驗持續了164 h,共采集984個文件,采樣間隔為10 min,每次采集20 480個點,試驗結束后發現1號軸承外圈出現局部損傷,說明第2組試驗中只有該軸承的數據為全壽命數據,本文便對該軸承的實測數據進行分析,圖5為軸承振動信號的均方根值趨勢圖。

圖5 1號軸承振動信號的均方根值趨勢圖

圖5中均方根值的變化趨勢基本反映了軸承從正常狀態到故障狀態的全過程。在7020 min之前,軸承運行相對穩定;當運行至7020 min時,振動信號均方根值發生微小突變,表明軸承狀態發生異常;在7020~9000 min時間段內,均方根值上下浮動,但變化幅度不大,說明軸承故障程度并不嚴重;當運行時間超過9000 min后,軸承振動信號均方根值呈連續遞增變化趨勢,并在9790 min時達到最大,此時軸承已達到壽命極限,出現嚴重故障。由于軸承早期故障的沖擊特征較微弱,容易被噪聲掩蓋而無法直接從原始波形中觀察到,所以振動信號的均方根值對早期故障并不敏感。

為驗證本文方法在軸承早期故障診斷應用上的可靠性,對5410 min時采集的信號進行分析,分析點數為8192,實測信號如圖6所示,通過波形觀察不到明顯的故障特征,對信號做進一步包絡譜分析,如圖7所示,也沒有發現幅值較高的頻率成分。

圖6 實測信號

圖7 實測信號的包絡譜

下面利用本文方法對信號進行處理,結果如圖8所示。首先根據采樣頻率fs=20 kHz、軸承內圈特征頻率fi=296.9 Hz、外圈特征頻率fo=236.4 Hz、滾動體特征頻率fb=140 Hz及保持架特征頻率fc=14.8 Hz得到最佳周期參數的搜索中心和范圍,其中獲得的4個搜索中心分別為Ti=67,To=85,Tb=143,Tc=1351,與搜索中心相對應的4個搜索范圍分別為[62,72]、[80,90]、[138,148]、[1346,1356],繼而繪制出圖8a中的T-S曲線。通過分析可知解卷積周期為85時稀疏度指標最大,則85即為最佳的解卷積周期參數Tm,設定MCKD算法的解卷積周期參數為85并對實測信號進行處理,所得結果如圖8b所示。進一步計算解卷積信號的包絡譜,結果如圖8c所示,包絡譜中外圈故障特征頻率fo及其諧波處幅值突出,由此可斷定軸承外圈存在局部損傷,理論分析結果與實際情況完全一致。

(a)T-S曲線

(b)最大相關峭度解卷積信號

(c)解卷積信號的包絡譜圖8 本文提出方法的實測信號分析結果

為進一步驗證基于包絡譜稀疏度和最大相關峭度解卷積的早期故障診斷方法的優勢,利用文獻[5]中基于EMD降噪和譜峭度的軸承早期故障診斷方法對實測信號進行分析,結果如圖9所示。圖9a為EMD降噪后得到的信號,通過計算降噪信號的快速峭度圖(結果如圖9b所示)獲得最佳帶通濾波中心及濾波帶寬,對濾波信號作包絡譜分析,結果如圖9c所示,圖中沒有出現任何頻率尖峰,無法提取出故障特征頻率信息。圖10是利用文獻[6]中基于最小熵解卷積的軸承微弱故障特征提取方法對原信號進行分析得到的結果,圖10b中僅外圈特征頻率2倍頻和3倍頻處譜線幅值偏高,但是峰值不夠突出,診斷效果與本文方法相比差距明顯。滾動軸承全壽命周期加速試驗數據對比分析結果表明,本文所述方法在實際軸承診斷應用中具有一定優勢,能夠有效提取出早期故障階段的微弱特征信息。

5 結語

(a)EMD降噪信號

(b)降噪信號的快速峭度圖

(c)濾波信號的包絡譜圖9 基于EMD降噪和譜峭度的實測信號分析結果

(a)最小熵解卷積信號

(b)解卷積信號的包絡譜圖10 基于最小熵解卷積的實測信號分析結果

本文提出了基于包絡譜稀疏度和最大相關峭度解卷積的滾動軸承早期故障診斷方法,仿真和試驗信號分析表明,該方法具有一定準確性和可靠性。筆者提出利用包絡譜稀疏度指標來篩選最佳解卷積周期參數,旨在起到一個拋磚引玉的效果。

[1]趙志宏,楊紹普,李韶華.基于Hilbert譜奇異值的軸承故障診斷[J].中國機械工程,2013,24(3):346-350.

Zhao Zhihong,Yang Shaopu,Li Shaohua.Bearing Fault Diagnosis Based on Hilbert Spectrum and Singular Value Decomposition[J].China Mechanical Engineering,2013,24(3):346-350.

[2]曾慶虎,邱靜,劉冠軍,等.基于小波相關濾波-包絡分析的早期故障特征提取方法[J].儀器儀表學報,2008,29(4):729-933.

Zeng Qinghu,Qiu Jing,Liu Guanjun,et al.Approach to Extraction of Incipient Fault Features Based on Wavelet Correlation Filter and Envelope Analysis[J].Chinese Journal of Scientific Instrument,2008,29(4):729-933.

[3]崔玲麗,康晨暉,胥永剛,等.滾動軸承早期沖擊性故障特征提取的綜合算法研究[J].儀器儀表學報,2010,31(11):2422-2427.

Cui Lingli,Kang Chenhui,Xu Yonggang,et al.Integrated Algorithm Research on Early Inpactive Fault Feature Extraction of Rolling Bearings[J].Chinese Journal of Scientific Instrument,2010,31(11):2422-2427.

[4]莫代一,崔玲麗,王婧.基于雙重Q因子的稀疏分解法在滾動軸承早期故障診斷中的應用[J].機械工程學報,2013,49(9):37-41.

Mo Daiyi,Cui Lingli,Wang Jing.Sparse Signal Decomposition Method Based on the Dual Q-factor and Its Application to Rolling Bearing Early Fault Diagnosis[J].Journal of Mechanical Engineering,2013,49(9):37-41.

[5]蘇文勝,王奉濤,張志新,等.EMD降噪和譜峭度法在滾動軸承早期故障診斷中的應用[J].振動與沖擊,2010,29(3):18-21.

Su Wensheng,Wang Fengtao,Zhang Zhixin,et al.Application of EMD Denoising and Spectral Kurtosis in Early Diagnosis of Rolling Element Bearings[J].Journal of Vibration and Shock,2010,29(3):18-21.

[6]Jiang Ruilong,Chen Jin,Dong Guangming,et al.The Weak Fault Diagnosis and Condition Monitoring of Rolling Element Bearing Using Minimum Entropy Deconvolution and Envelop Spectrum[J].Engineering Science Engineers,Part C:Journal of Mechanical Engineering Science,2013,227(5):1116-1129.

[7]江瑞龍.基于最小熵解卷積的滾動軸承故障診斷研究[D].上海:上海交通大學,2013.

[8]McDonald G L,Zhao Qing,Zuo M J.Maximum Correlated Kurtosis Deconvolution and Application on Gear Tooth Chip Fault Detection[J].Mechanical Systems and Signal Processing,2012,33:237-255.

[9]Tse P W,Wang Dong.The Design of a New Sparsogram for Fast Bearing Fault Diagnosis:Part 1 of the Two Related Manuscripts that Have a Joint Title as “Two Automatic Vibration-based Fault Diagnostic Methods Using the Novel Sparsity Measurement-Parts 1 and 2”[J].Mechanical Systems and Signal Processing,2013,40:499-519.

[10]Qiu Hai,Lee J,Lin Jing,et al.Wavelet Filter-based Weak Signature Detection Method and Its Application on Rolling Element Bearing Prognostics[J].Journal of Sound and Vibration,2006,289:1066-1090.

(編輯袁興玲)

Diagnosis Method for Rolling Bearing Incipient Faults Based on Sparsity of Envelope Spectrum and Maximum Correlated Kurtosis Deconvolution

Tang GuijiWang Xiaolong

North China Electric Power University,Baoding,Hebei,071000

Early fault features of rolling bearings are very weak and are affected by environment noise seriously,so it is difficult to draw fault features.Aiming at solving this problem,MCKD was tried to diagnose faults for bearings,and sparsity of envelope spectrum was used to select the optimal deconvolution period parameter,then incipient fault diagnosis method for rolling bearings was proposed based on sparsity of envelope spectrum and MCKD.MCKD method corresponding to the optimal parameter was used to process the original signals and the envelope spectrum of deconvolution signals was obtained,the bearing faults were judged by analyzing the envelope spectrum.Simulated incipient fault signals and full lifetime datasets of rolling bearings were used to examine the feasibility of this method and the results show the new method can be applied to diagnose the incipient fault effectively.

rolling bearing;sparsity;maximum correlated kurtosis deconvolution(MCKD);fault diagnosis

2014-02-17

中央高校基本科研業務費專項資金資助項目(13QN49);河北省自然科學基金資助項目(E2014502052)

TH133.3;TH17< class="emphasis_italic">DOI

:10.3969/j.issn.1004-132X.2015.11.006

唐貴基,男,1962年生。華北電力大學能源動力與機械工程學院教授、博士研究生導師。主要研究方向為機械故障診斷與狀態監測。王曉龍,男,1989年生。華北電力大學能源動力與機械工程學院博士研究生。

猜你喜歡
故障診斷故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 少妇高潮惨叫久久久久久| 中字无码精油按摩中出视频| 国产在线一区视频| 亚洲91精品视频| 国产浮力第一页永久地址| 综合久久五月天| 国产自视频| 国产www网站| 九九久久精品免费观看| 亚洲伊人久久精品影院| 亚洲Va中文字幕久久一区 | 54pao国产成人免费视频| www.91中文字幕| 91精品啪在线观看国产91| 日韩 欧美 小说 综合网 另类| 国产亚洲一区二区三区在线| 国产二级毛片| 毛片在线看网站| 3344在线观看无码| 亚洲精品色AV无码看| 免费av一区二区三区在线| 天天躁夜夜躁狠狠躁图片| 伊人久综合| 国产香蕉在线视频| 国产精品尤物在线| 日韩在线第三页| 精品人妻无码区在线视频| 免费毛片网站在线观看| 久久精品无码一区二区国产区| 精品一区二区三区自慰喷水| 亚洲国产成人精品一二区| 日韩精品成人在线| 欧美一级专区免费大片| 国外欧美一区另类中文字幕| 国产高清在线观看| 亚洲精品成人片在线观看| 亚洲无码在线午夜电影| 成人小视频在线观看免费| 91亚瑟视频| 一级片一区| 国产第一福利影院| 国内毛片视频| 国产精品自在自线免费观看| 天天做天天爱夜夜爽毛片毛片| 国产成人91精品| 亚洲国产日韩在线观看| 久久精品一品道久久精品| 久久精品波多野结衣| 国产福利在线免费| 国产精品无码在线看| 成人在线不卡视频| 国产欧美日韩综合一区在线播放| 欧美日韩理论| 欧美中日韩在线| 国产亚洲欧美在线专区| 精品国产香蕉在线播出| 国产香蕉国产精品偷在线观看| 在线观看91香蕉国产免费| 国产人在线成免费视频| 久久毛片网| 国产精品太粉嫩高中在线观看| 国产尤物在线播放| 五月天在线网站| 91在线丝袜| 91精品久久久久久无码人妻| 五月婷婷亚洲综合| 国产一区二区三区在线无码| 动漫精品啪啪一区二区三区| 久久久精品国产亚洲AV日韩| 国产精品网址你懂的| 国产地址二永久伊甸园| 亚洲人成影院在线观看| 啦啦啦网站在线观看a毛片| 欧美a级完整在线观看| 巨熟乳波霸若妻中文观看免费 | 成年人国产网站| 国产欧美视频一区二区三区| 高清乱码精品福利在线视频| 国产成人精品在线1区| 午夜综合网| 国产高清自拍视频| 亚洲国产精品日韩专区AV|