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

一種改進的峭度圖方法及其在復雜干擾下軸承故障診斷中的應用

2017-12-27 10:48:29顧曉輝楊紹普劉永強廖英英
振動與沖擊 2017年23期
關鍵詞:故障信號方法

顧曉輝, 楊紹普, 劉永強, 廖英英

(1.石家莊鐵道大學 交通運輸學院,石家莊 050043;2.河北省交通工程結構力學行為演變與與控制重點實驗室,石家莊 050043)

一種改進的峭度圖方法及其在復雜干擾下軸承故障診斷中的應用

顧曉輝1,2, 楊紹普1,2, 劉永強2, 廖英英2

(1.石家莊鐵道大學 交通運輸學院,石家莊 050043;2.河北省交通工程結構力學行為演變與與控制重點實驗室,石家莊 050043)

快速峭度圖是一種常用的滾動軸承故障診斷方法,但由于峭度指標對沖擊過于敏感,在干擾較復雜的工況中,該方法往往無法正確識別出最優的共振頻帶進行包絡解調。然而,解調信號的包絡譜對噪聲具有一定的免疫能力,而且包絡譜中通常會清晰的出現故障特征頻率及其倍頻成分,呈現出典型的周期性脈沖特點。因此,提出應用相關峭度定量地刻畫窄帶信號的包絡譜幅值,即以頻域相關峭度值生成峭度圖,用于最優頻帶的自適應地識別。同時,基于相關峭度的指向性,可以將該方法應用于軸承的復合故障診斷。最后通過實驗分析,驗證了該方法對軸承微弱故障和復合故障診斷的有效性。

峭度圖;頻域相關峭度;包絡分析;滾動軸承;故障診斷

滾動軸承作為旋轉機械中的關鍵部件,在機床、電機、機車車輛等機械設備中發揮著不可替代的作用。但通常滾動軸承所處的工作環境也是最復雜、最惡劣的,在服役過程中極易在內、外圈滾道或滾動體上產生點蝕、剝落等局部缺陷。故障的產生和發展輕則影響設備的工作精度,重則導致機毀人亡的重大災難事故。因此,研究有效的滾動軸承故障診斷方法十分必要。

振動信號常用于軸承的故障診斷研究,通常認為故障沖擊會激起軸承元件或軸承座等機械結構的共振響應,通過識別共振頻帶進行包絡解調,即可有效分離出信噪比較高的故障沖擊進而通過沖擊間隔診斷出故障發生位置[1]。為了自適應地識別共振頻帶,Antoni創造性地提出了譜峭度(Spectral Kurtosis)理論[2-3]和快速峭度圖(Fast Kurtogram)方法[4],在提出至今的10年間得到了充分的肯定和廣泛的關注[5]。與此同時,眾多學者在應用的過程中也相繼提出了多種優秀的改進方案以提升原方法的診斷能力。根據不同的出發點,大致可以歸為以下三類:① 更精確的帶通濾波器:基于內積變換原理,文獻[6-8]分別應用Db小波、提升多小波、諧波小波替換Antoni的短時傅里葉變換(STFT)和半解析有限沖擊響應(FIR)濾波器組,以改善分解效果,最大程度地從染噪信號中提取與小波基函數相似的故障特征。② 更完善的頻帶分割:Antoni采用金字塔式的1/3-二叉樹來劃分頻率-頻率分辨率平面,將信號分解到不同的頻帶上。文獻[9]認為這種自上而下的劃分方法易將最優的共振頻帶分割,由此提出一種自適應的相鄰頻帶融合方法。文獻[10]在其基礎上進一步提出了一種多尺度的頻帶融合方法。文獻[11]提出了一種基于典型故障頻率的頻帶劃分方法。③ 更魯棒的指標:Antoni應用的峭度是一種刻畫非高斯性的四階統計量,峭度對故障沖擊敏感但同時也無法克服對沖擊性噪聲敏感的缺陷,特別是在故障形成的早期階段。為彌補以上不足,文獻[12-13]通過引入新的統計模型對軸承的故障響應進行更準確的描述。文獻[14]提出了一種譜峭度平均的改進方法。文獻[15]從故障沖擊的周期性考慮,在診斷中采用相關峭度[16]代替峭度來識別共振頻帶。此外,文獻[17-18]提出采用包絡譜幅值的峭度來代替時域信號的峭度,即認為共振頻帶對應的解調信號的包絡譜是稀疏的。文獻[19]在其基礎上提出了Sparsogram方法。更進一步,文獻[20]結合時域的沖擊性和頻域的稀疏性提出了一種時-頻集總峭度。近期,Antoni結合故障信號的沖擊性和循環平穩性,提出了一種新的Infogram方法[21]。

由以上分析可知,如何準確地刻畫解調信號是識別共振頻帶的關鍵所在。由于滾動體滑動等隨機因素的影響,軸承中產生的重復性故障沖擊在時域并不具有嚴格的周期性[22],但在頻域即對應的包絡譜中,故障特征頻率的周期性通常可以得到保持?;谏鲜鲈?,本文從前文分析的第三個角度出發,提出了一種改進的峭度圖方法,采用相關峭度而非峭度刻畫包絡譜幅值來克服復雜干擾對診斷結果的影響。此外,相關峭度具有明確的指向性,通過輸入不同的故障特征頻率,可以將峭度圖方法擴展至軸承的復合故障診斷,使其具有更廣泛的應用價值。

1 Fast Kurtogram和Infogram方法簡介

1.1 Fast Kurtogram方法

一個條件非平穩過程的Wold-Cramér分解可以表示為

(1)

式中,H(t,f)為過程y(t)的時頻復包絡。在此基礎上,Antoni給出了譜峭度(SK)的定義:

(2)

式中,〈·〉表示平均,|·|表示取模?;诮y計理論,過程y(t)的非高斯性越強,四階累積量越大,也就是峭度越大,即SK可以用來度量y(t)在頻率f處的峰值特性。因此,利用SK,我們能夠識別出信號中的非平穩分量并可以確定瞬態分量所在的頻帶。

由前文分析可知,H(t,f)可以通過多種方式得到,為了便于對比新方法與Infogram方法的診斷效果,本文同樣選用STFT計算被分析信號的時頻復包絡:

(3)

式中,w(τ)為時窗。對于采樣長度為N、采樣頻率為Fs的離散非平穩信號y(n)的STFT可以表示為

(4)

式中,P為時窗的移動步長。

當STFT的時窗長度Nw取不同的值時,可以構造出頻率分辨率Δf不同的帶通濾波器組,即Δf~Fs/Nw。由此,基于式(2)和式(4),我們可以得到SK在(f,Δf)平面內的分布,即如圖1所示的峭度圖。

圖1 基于STFT的峭度圖Fig.1 The STFT-based kurtogram

1.2 Infogram方法

在信噪比較高的情況下,應用Fast Kurtogram方法通??梢詼蚀_地找到軸承的共振頻帶進行包絡解調,但在故障沖擊不明顯或信號中存在高峰值的脈沖干擾時,該方法易于失效。這是因為峭度對離群野值過于敏感,僅能刻畫故障信號的沖擊性而無法表征故障沖擊的循環平穩性。為此,Antoni基于時域的譜負熵和頻域的譜負熵提出了Infogram方法。其中,時域的譜負熵定義為

(5)

式中,SE為平方包絡,即SE=H(kP,f)2。

頻域的負熵定義為

(6)

式中,SES為平方包絡譜即SES=DFT[SE(n;f,Δf)]。

ΔIe和ΔIE分別可以表征時域的平方包絡信號和頻域的平方包絡譜幅值的能量流動特性,即分別用來表征故障信號的沖擊特性和循環平穩特性。此外,根據熵的不確定原理,Antoni進一步提出了平均譜負熵用來統一表征沖擊特性和循環平穩特性,即:

ΔI1/2=ΔIe/2+ΔIE/2

(7)

然而,Infogram方法采取的時、頻域平均并沒有從根本上克服沖擊性干擾的影響,但給如何設計更魯棒的指標指明了方向,即如何用統一表征沖擊特性和循環平穩特性。

2 基于頻域相關峭度的改進峭度圖方法

考慮如下仿真信號:

(8)

式中:Ak為軸承故障沖擊的幅值;T0為故障沖擊的間隔;τk為滾動體滑動造成的隨機偏差;θ(t)為高幅值的沖擊性干擾;n(t)為均值為0、標準差為1的高斯白噪聲;h(t)為沖擊響應函數:

(9)

式中:bw為帶寬參數;f0為中心頻率。圖2為噪聲強度δ分別取0、0.1、0.2、0.3、0.4和0.5時,各仿真信號的時域波形及包絡譜。

通過對比可知,即使在噪聲很強的情況下,在頻域的包絡譜中依然可以找到故障特征頻率及其高階倍頻。因此,作者認為利用McDonald提出的相關峭度指標定量的描述包絡譜中的周期性故障特征頻率,即利用頻域相關峭度(FDCK)可以統一故障信號的沖擊特性和循環平穩特性。FDCK的定義為

FDCK(T)=

(10)

式中,T為待檢測的故障特征頻率,根據不同的輸入,由FDCK生成的峭度圖可以識別不同故障對應的共振頻帶。

3 軸承微弱故障的診斷分析

為了驗證新方法對軸承微弱故障診斷的有效性,應用QPZZ-Ⅱ旋轉機械振動及故障模擬實驗平臺進行驗證。測試軸承型號為6 205,采用電火花加工技術在外圈滾道上加工了直徑為0.2 mm的點蝕故障。實驗臺及加工的故障形貌如圖3所示。

(a) δ=0

(b) δ=0.1

(c) δ=0.2

(d) δ=0.3

(e) δ=0.4

(f) δ=0.5圖2 不同噪聲強度下的仿真信號及其包絡譜Fig.2 The simulated signals and the corresponding envelope spectrums under different noise

實驗時,轉速設定為1 478 r/min,采樣頻率為25.6 kHz,應用加速度傳感器采集1 s的數據,故障信號的時域波形及頻譜如圖4所示。根據軸承尺寸,計算得到軸承的內圈故障特征頻率fi、外圈故障特征頻率fo、滾動體故障特征頻率fb分別為133.40 Hz、88.31 Hz和58.05 Hz。

首先,應用本文方法對圖4信號進行分析,為了確保結果的可靠性,以最底層的帶寬不小于2倍的最大故障特征頻率為原則選取最大的分解層數,得到的峭度圖如圖5所示。從中可知,最大的FDCK對應的中心頻率和帶寬分別為2 200 Hz和400 Hz,該窄帶信號的時域波形及平方包絡譜如圖6所示。時域波形中出現了不太明顯的故障沖擊,但從平方包絡譜中可以很容易的找到外圈故障特征頻率fo及其倍頻,表明該方法可以成功識別出軸承的微弱外圈故障。

圖3 實驗臺及故障形貌Fig.3 The test bench and the artificial fault

(a) 時域波形

(b) 頻譜圖4 故障軸承振動信號及其頻譜Fig.4 Original vibration signal and its spectrum of the faulty bearing

圖7為應用Infogram方法對相同的信號進行分析的結果。受噪聲的影響,表征沖擊性特征的時域譜負熵ΔIe確定的中心頻率和帶寬分別為200 Hz和400 Hz,顯然這不是軸承的共振頻帶。表征循環平穩特征的頻域譜負熵ΔIE確定的中心頻率和帶寬分別為2 200 Hz和400 Hz,與本文方法識別的結果相同。然而,最終的平均譜負熵ΔI1/2仍然未能克服噪聲的干擾,確定的確定的中心頻率和帶寬同樣分別為200 Hz和400 Hz,對應的窄帶信號及包絡譜如如8所示。

圖5 基于FDCK的峭度圖Fig.5 FDCK-based kurtogram

(a) 窄帶信號的時域波形

(b) 平方包絡譜圖6 本文方法得到的窄帶信號及其平方包絡譜

Fig.6 Filtered signal by the proposed method and its squred envelope spectrum

(a) 基于ΔIe的infogram

(b) 基于ΔIE的infogram

(c) 基于ΔI1/2的infogram圖7 三種InfogramFig.7 Three kinds infograms

(a) 窄帶信號的時域波形

(b) 平方包絡譜圖8 由Infogram方法得到的窄帶信號及其平方包絡譜Fig.8 Filtered signal by Infogram method and its squred envelope spectrum

4 軸承復合故障的診斷分析

軸承的復合故障相互耦合、相互干擾,具有一定的診斷難度。基于FDCK的指向性,本文進一步將新方法應用于復合故障的分析,并采用如圖9所示的貨車輪對軸承跑合實驗臺驗證了診斷效果。測試軸承型號為197726,軸承故障為服役過程中自然形成,除了內圈滾道出現明顯的剝落外,在外圈滾道出現了幾條輕微的壓痕。

圖9 輪對軸承實驗臺及故障形貌Fig.9 The test bench of wheel set bearing and the compound faults

實驗時,轉速設定為465 r/min,采樣頻率為25.6 kHz,應用加速度傳感器采集1 s的數據,故障信號的時域波形及頻譜如圖10所示。根據軸承尺寸,計算得到軸承的內圈故障特征頻率fi、外圈故障特征頻率fo、滾動體故障特征頻率fb分別為88.25 Hz、66.75 Hz和27.08 Hz。

圖11為設定T=66.75時,應用本文方法得到的峭度圖。從中可知,最大的FDCK對應的中心頻率和帶寬分別為9 700 Hz和200 Hz,該窄帶信號的時域波形及平方包絡譜如圖12所示。時域波形中出現了多個清晰的故障沖擊,從包絡譜中可以很容易的找到外圈故障對應的特征頻率fo。因此,可以判斷軸承外圈存在故障。

(a) 時域波形

(b) 頻譜圖10 故障軸承振動信號及其頻譜Fig.10 Original vibration signal and its spectrum of the faulty bearing

圖11 T=fo時基于FDCK的峭度圖Fig.11 FDCK-based kurtogram when T=fo

(a) 窄帶信號的時域波形

(b) 平方包絡譜圖12 T=fo時本文方法得到的窄帶信號及其平方包絡譜

Fig.12 Filtered signal by the proposed method and its squred envelope spectrum whenT=fo

同理,圖13為設定T=88.25時,應用本文方法得到的峭度圖。從中可知,最大的FDCK對應的中心頻率和帶寬分別為3 500 Hz和200 Hz,該窄帶信號的時域波形及平方包絡譜如圖14所示。時域波形中出現了更密集的故障沖擊,并且從包絡譜中可以很容易的找到內圈故障對應的特征頻率fi及其倍頻和以轉頻fr為間隔的調制現象,即可以判斷軸承內圈同樣存在故障。

圖15~16為應用Infogram方法對該復合故障信號的分析結果,由于負熵并不具有指向性,該方法僅能識別出軸承中的外圈故障。

圖13 T=fi時基于FDCK的峭度圖Fig.13 FDCK-based kurtogram when T=fi

(a) 窄帶信號的時域波形

(b) 平方包絡譜圖14 T=fi時本文方法得到的窄帶信號及其平方包絡譜Fig.14 Filtered signal by the proposed method and its squred envelope spectrum when T=fi

(a) 基于ΔIe的infogram

(b) 基于ΔIE的infogram

(c) 基于ΔI1/2的infogram圖15 三種InfogramFig.15 Three kinds infograms

(a) 窄帶信號的時域波形

圖16 由Infogram方法得到的窄帶信號及其平方包絡譜Fig.16 Filtered signal by Infogram method and its squred envelope spectrum

5 結 論

本文從包絡譜中的故障特征頻率分布出發,提出了一種新的指標——頻域相關峭度,改進了峭度圖方法。新方法可以有效消除復雜干擾噪聲的影響,自適應地識別早期微弱故障對應的最優共振頻帶,提高了包絡分析中中心頻率和帶寬參數選擇的準確性。此外,基于頻域相關峭度的指向性,新方法同時可以應用于軸承的復合故障診斷。根據不同的故障特征頻率識別不同的故障對應的共振頻帶。通過對兩組實驗信號進行分析,驗證了新方法的有效性和相對于Infogram方法的優越性。

[1] MCFADDEN P D, SMITH J D. Vibration monitoring of rolling element bearings by the high-frequency resonance technique—a review[J]. Tribology International, 1984, 17(1): 3-10.

[2] ANTONI J. The spectral kurtosis: a useful tool for characterising non-stationary signals[J]. Mechanical Systems and Signal Processing, 2006, 20(2): 282-307.

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

[4] ANTONI J. Fast computation of the kurtogram for the detection of transient faults[J]. Mechanical Systems and Signal Processing, 2007, 21(1): 108-124.

[5] WANG Y, XIANG J, MARKERT R, et al. Spectral kurtosis for fault detection, diagnosis and prognostics of rotating machines: A review with applications[J]. Mechanical Systems and Signal Processing, 2016, 66: 679-698.

[6] LEI Y, LIN J, HE Z, et al. Application of an improved kurtogram method for fault diagnosis of rolling element bearings[J]. Mechanical Systems and Signal Processing, 2011, 25(5): 1738-1749.

[7] 王曉冬, 何正嘉, 訾艷陽. 滾動軸承故障診斷的多小波譜峭度方法[J]. 西安交通大學學報, 2010, 44(3): 77-81.

WANG Xiaodong, HE Zhengjia, ZI Yanyang. Spectral kurtosis of multiwavelet for fault diagnosis of rolling bearing[J]. Journal of Xi’an Jiaotong University, 2010, 44(3): 77-81.

[8] 田福慶, 羅榮, 李萬, 等. 改進的諧波小波包峭度圖及其應用[J]. 上海交通大學學報, 2014, 48(1): 39-44.

TIAN Fuqing, LUO Rong, LI Wan, et al. Improved harmonic wavelet packet kurtogram and its application[J]. Journal of Shanghai Jiaotong University, 2014, 48(1): 39-44.

[9] WANG Y, LIANG M. An adaptive SK technique and its application for fault detection of rolling element bearings[J]. Mechanical Systems and Signal Processing, 2011, 25(5): 1750-1764.

[10] LI C, CABRERA D, DE OLIVEIRA J V, et al. Extracting repetitive transients for rotating machinery diagnosis using multiscale clustered grey infogram[J]. Mechanical Systems and Signal Processing, 2016, 76: 157-173.

[11] 馬新娜, 楊紹普. 典型譜峭圖在共振解調方法中的應用[J]. 振動. 測試與診斷, 2015, 35(6): 1140-1144.

MA Xinna, YANG Shaopu. Study and application of demodulated resonance based on typi-kurtogram[J]. Journal of Vibration, Measurement & Diagnosis, 2015, 35(6): 1140-1144.

[12] YU G, LI C, ZHANG J. A new statistical modeling and detection method for rolling element bearing faults based on alpha-stable distribution[J]. Mechanical Systems and Signal Processing, 2013, 41(1): 155-175.

[13] OBUCHOWSKI J, WYA, ZIMROZ R. Selection of informative frequency band in local damage detection in rotating machinery[J]. Mechanical Systems and Signal Processing, 2014, 48(1): 138-152.

[14] 代士超, 郭瑜, 伍星, 等. 基于子頻帶譜峭度平均的快速譜峭度圖算法改進[J]. 振動與沖擊, 2015, 34(7): 98-102.

DAI Shichao, GUO Yu, WU Xing, et al. Improvement on fast kurtogram algorithm based on sub-frequency-ban spectral kurtosis average[J]. Journal of Vibration and Shock, 2015, 34(7): 98-102.

[15] ZHANG X, KANG J, XIAO L, et al. A new improved kurtogram and its application to bearing fault diagnosis[J]. Shock and Vibration, 2015(3): 1-22.

[16] MCDONALD G L, ZHAO Q, 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.

[17] BBRSZCZ T, JABOSKI A. A novel method for the optimal band selection for vibration signal demodulation and comparison with the Kurtogram[J]. Mechanical Systems and Signal Processing, 2011, 25(1): 431-451.

[18] WANG D, PETER W T, TSUI K L. An enhanced Kurtogram method for fault diagnosis of rolling element bearings[J]. Mechanical Systems and Signal Processing, 2013, 35(1): 176-199.

[19] PETER W T, WANG D. 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(2): 499-519.

[20] CHEN B Q, ZHANG Z S, ZI Y Y, et al. Detecting of transient vibration signatures using an improved fast spatial-spectral ensemble kurtosis kurtogram and its applications to mechanical signature analysis of short duration data from rotating machinery[J]. Mechanical Systems and Signal Processing, 2013, 40(1): 1-37.

[21] ANTONI J. The infogram: Entropic evidence of the signature of repetitive transients[J]. Mechanical Systems and Signal Processing, 2016, 74: 73-94.

[22] HO D, RANDALL R B. Optimisation of bearing diagnostic techniques using simulated and actual bearing fault signals[J]. Mechanical Systems and Signal Processing, 2000, 14(5): 763-788.

Animprovedkurtogrammethodanditsapplicationinfaultdiagnosisofrollingelementbearingsundercomplexinterferences

GU Xiaohui1,2, YANG Shaopu1,2, LIU Yongqiang2, LIAO Yingying2

(1. School of Traffic and Transportation, Shijiazhuang Tiedao University, Shijiazhuang 050043, China;2. Key Laboratory of Mechanical Evolution and Control of Traffic Structure in Hebei, Shijiazhuang 050043, China)

Fast Kurtogram is one of the most useful methods in fault diagnosis of rolling element bearings. However, in some cases of complex interferences, it cannot exactly recognize the optimal resonance frequency band for envelope demodulation due to that the kurtosis index is too sensitive to impulsive noise. In fact, the envelope spectrum of demodulated signals in frequency domain has a certain immunity ability to noise, the bearing fault characteristic frequency and its harmonics often appear clearly with typical periodic impulse features in the envelope spectrum. Here, the frequency domain correlated kurtosis was proposed to quantitatively describe envelope spectrum amplitudes of narrow-band signals and generate the kurtogram. Moreover, the proposed method can be applied in the compound fault detection based on the directivity of correlated kurtosis. Two cases of real bearing fault signals were employed to verify the effectiveness and robustness of the proposed method in bearing weak fault diagnosis and compound fault diagnosis.

kurtogram; frequency domain correlated kurtosis; envelope analysis; rolling element bearing; fault diagnosis

國家自然科學基金(11227201;U1534204;11472179;11572206;11302137;11372197);河北省自然科學基金(A2016210099)

2016-11-08 修改稿收到日期:2017-02-09

顧曉輝 男,博士生,1988年11月生

楊紹普 男,教授,博士生導師,1962年10月生

TH113.1;TH133

A

10.13465/j.cnki.jvs.2017.23.028

猜你喜歡
故障信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
故障一點通
主站蜘蛛池模板: 久久久四虎成人永久免费网站| 国产精品久久久久婷婷五月| 亚洲精品色AV无码看| 国产成人成人一区二区| 亚洲男人在线天堂| 欧洲亚洲一区| 亚洲人成色77777在线观看| 秘书高跟黑色丝袜国产91在线| 免费观看国产小粉嫩喷水| 亚洲精品你懂的| 亚洲香蕉伊综合在人在线| 国产95在线 | 亚洲人成日本在线观看| 国产h视频在线观看视频| 日韩AV无码一区| 最新亚洲人成网站在线观看| 日韩中文字幕亚洲无线码| 欧美人与动牲交a欧美精品| 久久夜色精品国产嚕嚕亚洲av| 直接黄91麻豆网站| 一级毛片不卡片免费观看| 亚洲aaa视频| 国产日产欧美精品| 午夜日b视频| 中文字幕亚洲无线码一区女同| 久久香蕉国产线看观看式| 中文字幕1区2区| 国产情精品嫩草影院88av| 亚洲天堂网在线播放| 国产成人资源| 久久婷婷六月| 国产xx在线观看| 亚洲第一区在线| 婷婷激情亚洲| 久久亚洲美女精品国产精品| 亚洲人成人伊人成综合网无码| 亚洲综合激情另类专区| 国产真实自在自线免费精品| 超碰免费91| 亚洲AV无码乱码在线观看裸奔| 日本精品视频| 中文字幕无码中文字幕有码在线| 欧美 亚洲 日韩 国产| 99r在线精品视频在线播放| 一级毛片不卡片免费观看| 99精品在线看| 亚洲人成人无码www| 一级片免费网站| 日韩一区二区在线电影| 第一区免费在线观看| 偷拍久久网| 欧美成人免费午夜全| 色呦呦手机在线精品| 久久综合伊人 六十路| 久久天天躁狠狠躁夜夜躁| 国产亚洲欧美另类一区二区| 国产精品流白浆在线观看| 欧美日在线观看| 红杏AV在线无码| 国产香蕉一区二区在线网站| 秘书高跟黑色丝袜国产91在线 | 精品国产自| 无码一区中文字幕| 久久综合五月| 欧美、日韩、国产综合一区| 国产资源免费观看| 婷婷色中文网| 一本一道波多野结衣av黑人在线| 久久综合结合久久狠狠狠97色| 色哟哟精品无码网站在线播放视频| 中文一级毛片| 欧美精品综合视频一区二区| 国产精品美女免费视频大全| 国产精品无码AⅤ在线观看播放| 又爽又大又黄a级毛片在线视频| 精品欧美一区二区三区久久久| A级毛片无码久久精品免费| 波多野结衣国产精品| 91年精品国产福利线观看久久 | 久久人搡人人玩人妻精品一| 99精品免费在线| 国产原创第一页在线观看|