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

EEMD和TFPF聯合降噪法在齒輪故障診斷中的應用

2017-11-07 05:40:41寧少慧韓振南武學峰
振動、測試與診斷 2017年5期
關鍵詞:故障診斷模態振動

寧少慧, 韓振南, 武學峰, 趙 遠

(1.太原理工大學機械工程學院 太原,030024) (2.太原科技大學機械工程學院 太原,030024)

10.16450/j.cnki.issn.1004-6801.2017.05.024

EEMD和TFPF聯合降噪法在齒輪故障診斷中的應用

寧少慧1,2, 韓振南1, 武學峰3, 趙 遠1

(1.太原理工大學機械工程學院 太原,030024) (2.太原科技大學機械工程學院 太原,030024)

為了消除噪聲對齒輪傳動系統故障特征提取的影響,提出了一種基于集成經驗模態分解(ensemble empirical mode decomposition,簡稱EEMD)和時頻峰值濾波(time-frequency peak filtering,簡稱TFPF)相結合的降噪方法。針對TFPF算法在窗長的選擇方面受到限制的問題,采用了EEMD方法對其進行改進,使得信號在噪聲壓制和有效信號保真兩方面得到權衡;含噪聲的信號經過EEMD分解后,得到一系列頻率成分從高到低的本征模態函數(intrinsic mode functions,簡稱IMFs),計算出各IMFs間的相關系數,判斷需要濾波的IMFs。對不同的IMFs選擇不同的窗長進行TFPF濾波,把過濾后的IMFs和剩余的IMFs重構得到最終的降噪信號。用模擬仿真信號和齒輪齒根故障信號對該方法進行驗證,可見EEMD+TFPF能有效地去除噪聲,成功提取齒根裂紋故障特征。

時頻峰值濾波; 集成經驗模態分解; 齒根裂紋; 降噪

引 言

在齒輪傳動系統的故障診斷中,最常見的是通過分析齒輪箱體振動信號提取系統故障特征。當齒輪箱中的齒輪、軸或軸承等機械設備產生故障時,箱體振動信號中除微弱的故障信息外,還會夾雜著各種頻率的背景噪聲,而且信號頻帶與噪聲頻帶常會相互交錯,使得傳統的信號處理方法很難從包含強背景噪聲的箱體振動信號中提取到微弱的故障信息[1]。集成經驗模態分解(ensemble empirical mode decomposition,簡稱EEMD)不僅保持了經驗模態分解(empirical mode decomposition,簡稱EMD)的自適應分解特性,還在本質上解決了EMD的模態混疊問題,被廣泛用于時變的非線性、非平穩的齒輪傳動系統的故障診斷中[2]。但由于EEMD分解后會得到多個本征模態函數(intrinsic mode functions,簡稱IMFs),每個IMFs中都有可能包含著故障頻率成分,如何從EEMD分解的多個IMFs中提取故障頻率成為EEMD在故障診斷領域應用的關鍵。

文獻[3]把EEMD和相關系數結合,把相關系數大于0.5的IMFs保留,其他的IMFs去掉,在去噪聲的同時可能會把有用的信號也去掉了。文獻[4]把EEMD和排列熵結合,成功識別了高速列車轉向架的故障狀態。文獻[5]提出了基于EEMD、形態譜特征提取和模糊C均值聚類集成法對齒輪傳動系統的軸承做出了故障診斷。文獻[6]提出了EEMD時頻譜二值化方法,通過多尺度二進譜分析得到信號的權重譜,將其向時域累計得到權重向量,實現了微弱沖擊特征的增強,成功提取齒輪傳動系統中軸承故障特征。文獻[7-8]分別把EEMD與快速譜峭度圖和最小熵反褶積結合起來降噪,診斷齒輪系統的軸承故障。

時頻峰值濾波(time-frequency peak filtering,簡稱TFPF)是由Mesbah等[9]提出的一種信號消噪算法,特點是能夠在強噪聲環境中提取出有效信號,被廣泛應用在地震探測信號處理中[10-12]。時頻峰值濾波法基于時頻分析理論來消減隨機噪聲,通過頻率調制將含噪信號調制為解析信號,利用解析信號維納維爾分布沿瞬時頻率最為集中的特性,將其峰值作為信號的瞬時頻率來提取有效信號。文獻[13]把該方法與時頻分布結合應用到傳動系統的軸承故障診斷中。在TFPF算法中,降噪的同時還能保持有用的振動信號的關鍵問題是窗長的選擇。對于頻率成分復雜的箱體振動信號,選擇大的窗長在有效地去除噪聲同時會造成有效信號的幅值損失;選擇小窗長雖能很好地保真有效信號,但在降噪方面的力度卻不夠。而信號保真和噪聲壓制對能否成功提取出故障信號至關重要。針對此問題,筆者提出TFPF與EEMD聯合降噪,振動信號在經過EEMD分解后得到從高頻到低頻的一系列IMFs,對不同頻段的IMFs采用不同的窗長,解決了TFPF在窗長的選擇方面的問題,有效地提取強噪聲背景下的齒根裂紋的故障特征。

1 時頻峰值濾波(TFPF)算法

1 TFPF消除噪聲的基本原理

TFPF算法是基于偽Wigner-Ville分布的瞬時頻率估計,首先將含噪聲信號編碼成解析信號,再計算解析信號的偽Wigner-Ville分布,并將其峰值作為信號的瞬時頻率,從而消除隨機噪聲。

齒輪箱體振動信號x(t)中總會存在噪聲,把它表示為

x(t)=s(t)+n(t)

(1)

其中:s(t)為振動信號中的有效信號;n(t)為背景噪聲。

利用時頻峰值濾波去除齒輪箱體振動信號的噪聲的步驟如下。

1) 對包含噪聲的信號x(t)進行頻率調制,將其變為解析信號z(t)

(2)

其中:μ為頻率調制指數。

2) 計算解析信號z(t)的偽Wigner-Ville分布頻譜

(3)

其中:z*為z(t)的共軛。

3) 根據最大似然估計原理,求解析信號z(t)的PWVD分布頻譜的峰值,作為解析信號的瞬時頻率估計,即可得到原始有效信號的幅值估計

(4)

1.2 窗長的選擇

通過仿真信號來說明TFPF算法中窗長選擇的重要性。設定一個多調制源的仿真信號,采樣頻率為1 024 Hz,采樣點數是1 024,調制頻率分別為fn1=18 Hz,fn2=40 Hz;載波頻率為fz=200 Hz。在仿真信號中加入白噪聲n(t),仿真信號表達式為

x(t)= [1+cos(2πfn1t)+cos(2πfn2t)]·

cos(2πfzt)+n(t)

(5)

在信號中加入噪聲后,比較在不同的窗長下TFPF的過濾結果,如圖1所示。從圖1(a)看出,純凈信號已被所加的噪聲嚴重污染,合理地消除噪聲,盡可能恢復原始信號對故障診斷非常重要。圖1(b)和(c)是選擇不同窗長降噪的效果,含噪信號在經過TFPF濾波后,噪聲都有所減小,說明TFPF算法在噪聲壓制方面是非常有效的。圖1(b)為長窗長的降噪結果,可以看出,信號的去噪效果雖好,但在幅度方面會有所損失,尤其在波峰與波谷的位置。圖1(c)為短窗長的濾波效果,可以看到信號的波形和幅度損失小,但在噪聲壓制方面有所欠缺,濾波后依然殘留很多噪聲成分,信噪比將有所下降。

因此,TFPF算法中窗長的大小直接影響到信號保真和噪聲壓制的效果。振動信號經過EEMD分解后,對不同頻率的信號分量選用不同的窗長進行濾波,既能有效地去除隨機噪聲,也能保真有效信號。

圖1 TFPF法不同窗長的降噪結果Fig.1 Different window denoising signal using TFPF

2 集成經驗模態分解(EEMD)原理

集成經驗模態分解的實質是給原始信號加入極小幅度的白噪聲,利用了白噪聲頻譜均衡分布的特點和零均值特性,經過多次平均后將噪聲相互抵消,消除噪聲對原信號的影響,很好地解決了模態混疊問題。EEMD具體分解步驟如下。

1) 向原始信號x(t)中多次加入零均值、幅值標準差為常數的白噪聲nj(t)(j=1,2,…,M),j表示加入白噪聲的次數;

xj(t)=x(t)+nj(t)

(6)

其中:xj(t)表示第j次加入白噪聲后的信號。

2) 將添加了白噪聲的信號xj(t)通過EEMD算法分解為一系列本征模態函數ci,j(i=1,2,…,I),其中ci,j表示第j次加入白噪聲幅值后,分解得到的第i個IMF;如果j

3) 將每次得到的對應IMFs的集成平均值作為最后的IMF

(7)

其中:ci表示EEMD得到的第i個IMF;ci ,j表示第j次加入白噪聲幅值后,分解得到的第i個IMF;M表示總體平均的次數。

3 基于EEMD和TFPF的降噪方法

為了解決TFPF在窗長的選擇方面的局限性,把EEMD和TFPF相結合更適合齒輪傳動系統的噪聲消除,具體步驟如下。

1) 對含噪信號進行EEMD分解,得到頻率成分由高到低的不同的IMFs。但每個模態分量不是純噪聲模態或者是純信號模態,而是信號成分和噪聲成分相互交叉的模態,所以不能直接丟棄將噪聲主導的模態,或者直接將信號主導模態保留。要判斷哪些模態需要進行濾波處理。

2) 通過公式(8)計算各模態間的互相關系數并判斷需要濾波的模態。一般情況下,如果計算出的兩個相鄰模態間的互相關系數從一個較大的值往后都相對比較穩定,那么此值就可作為模態分界的閾值,這兩個模態中后一個模態及之后的模態就無需進行濾波處理,僅對其之前的模態分量選取合適窗長的TFPF進行降噪處理[14]。計算各IMFs間的互相關系數公式為

(8)

3) 選擇不同窗長的TFPF對需要降噪的IMFs分量進行處理,原則為:高頻分量主要包含的信號為噪聲信號,選用長窗長去噪;低頻分量主要成分是有用信號,因此選取短窗長。

4) 將去噪后的模態和剩余模態重構得到最終的濾波信號。為了提取故障頻率特征,對重構信號進行了循環解調分析,其過程如圖2所示。

圖2 含噪信號的降噪過程Fig.2 Denoising process of signal

4 仿真分析

用仿真信號驗證EEMD+TFPF的有效性。對式(5)的仿真信號采用EEMD+TFPF法進行降噪,信號的EEMD分解結果如圖3所示。

圖3 信號的EEMD分解Fig.3 Signal decomposition IMFs by EEMD

各相鄰IMF間的互相關系數計算結果如表1。

表1 仿真信號IMFs間的互相關系數

從表1可知,相鄰IMFs的互相關系數從IMF4后的值較穩定,故只需對IMF1,IMF2和IMF3選擇不同的窗長進行TFPF處理。為了說明EEMD+TFPF的降噪優勢,對含噪信號也進行EEMD和TFPF降噪。圖4(a)是EEMD降噪結果,可以看出在去噪的同時,也去掉了部分有效信號。圖4(b)是EEMD+TFPF的聯合降噪,既保真了信號幅值又使噪聲得到了有效的壓制。TFPF降噪在前面已經討論過。

信噪比和均方差是衡量降噪效果的重要指標,通過計算噪聲信號,TFPF,EEMD和EEMD+TF-PF的SNR和MSE,進一步說明TFPF+EEMD降噪的優越性。結果如表2所示:TFPF+EEMD降噪的信噪比最大,均方差最小,說明了TFPF+EEMD降噪的優越性。

圖4 兩種方法的降噪結果Fig.4 The denoising signal by EEMD and EEMD+TFPF

參數噪聲信號TFPF短窗降噪信號EEMD降噪信號TFPF+EEMD降噪信號信噪比/dB-1.17760.30802.02406.0157均方差 2.89181.95731.64490.5932

當齒輪傳動系統出現故障時,其箱體振動信號為調制信號[1],要從振動信號中提取故障頻率特征,需要對其進行解調分析。對EEMD+TFPF的降噪信號進行解調分析,結果如圖5所示,低頻處出現18,40及58 Hz是原信號的調制頻率。在高頻處二倍載波頻率400 Hz明顯突出,以調制頻率18和40 Hz為間隔的邊頻帶特征也被很清晰的解調出來。因此,EEMD+TFPF聯合降噪,既保持了有用信號,又最大程度地去除了噪聲。

圖5 仿真信號經過EEMD+TFPF降噪后的循環自相關分析Fig.5 Cyclic autocorrelation function of simulation signal by EEMD+TFPF denoising

5 齒輪傳動系統實測信號分析

5.1 齒輪傳動實驗系統

單級齒輪傳動實驗系統如圖6所示,主要包括電機、齒輪陪試箱、實驗齒輪箱、加速度傳感器及扭力桿等。在齒輪箱的軸承座上安裝了4個壓電加速度傳感器,從動輪為故障齒輪。實驗時,由電機驅動整個傳動系統的運轉,通過扭力桿加載,從扭矩測量儀觀察轉速和轉矩的大小。主動輪齒數為30,從動輪齒數為45。在從動輪的齒根處,人工加工長度為2 mm的裂紋(圖7),采用了動態數據采集分析系統進行信號采集,對采集到的信號采用TFPF,EEMD和EEMD+TFPF降噪法降噪。

圖6 齒輪傳動實驗臺及局部放大Fig.6 Single-stage spur gearbox test rig

圖7 齒根裂紋Fig.7 Gear tooth root cracks

5.2 EEMD+TFPF降噪分析

將集成經驗模態分解和時頻峰值濾波結合起來進行去噪,當齒輪傳動系統的負載為323 N·m,轉速為1 186 r/min,采樣頻率為8 kHz時,對系統采集到的箱體振動信號作EEMD分解,如圖8所示。高頻分量成分以噪聲為主,直接去掉雖然可以達到降噪目的,但有可能丟失存在高頻分量中的有效信號。也不能對所有模態進行TFPF濾波,因為在低頻模態含有純信號成分。所以要通過計算各IMFs間的相關系數來判斷哪些IMFs需要濾波(見表3)。

圖8 實測故障信號的EEMD分解Fig.8 Experimental signal decomposition IMFs by EEMD

圖9 EEMD和EEMD+TFPF降噪3種降噪結果Fig.9 The denoising signal by EEMD and TFPF and EEMD+TFPF

圖10 實驗信號經過EEMD+TFPF降噪后的循環自相關函數分析Fig.10 Cyclic autocorrelation function of experiment signal by EEMD+TFPF denoising

IMFs相關系數IMFs相關系數IMF1與IMF20.0455IMF6與IMF70.3519IMF2與IMF30.0830IMF7與IMF80.3953IMF3與IMF40.1032IMF8與IMF90.3793IMF4與IMF50.2736IMF9與IMF100.4003IMF5與IMF60.2966IMF10與IMF110.2682

從表3中可知,需要對前四個模態進行TFPF降噪處理,然后將去噪后的模態與不需去噪的模態重構得到最終的降噪信號,這使得TFPF降噪僅作用于以噪聲為主的高頻IMFs,改善了直接使用TFPF降噪的缺陷,降噪結果如圖9。為了說明EEMD+TFPF降噪的優越性,對3種降噪結果做了對比。

計算原始信號、EEMD和EEMD+TFPF的SNR和MSE,如表4所示。從對比結果可以看出,經過TFPF+EEMD降噪的信噪比最大,均方差最小,說明了TFPF+EEMD降噪的優越性。

表4實驗信號3種降噪結果的信噪比和均方差

Tab.4ThreedenosingSNRandMSEofexperimentalsignals

信號原始信號TFPF降噪信號EEMD降噪信號TFPF+EEMD降噪信號信噪比/dB19.246637.113238.213942.7782均方差0.34530.19220.16380.1150

裂紋是齒輪箱多種故障中比較難識別的故障,當輪齒齒根出現裂紋時,振動信號的頻率成分和幅值都會發生變化,因此調幅效應和調頻效應同時存在,頻譜上的邊頻成分由于具有不同的相位,使得信號的調制邊頻帶不再對稱。對采用EEMD+TFPF降噪后的信號進行循環自相關解調分析,結果如圖10所示,縱坐標表示幅值,用A表示,橫坐標為循環頻率,用α表示。根據解調原理,調幅調頻效應同時存在會導致信號的循環域低頻段出現調制源的1倍頻和2倍頻及以上成分,高頻段出現以嚙合頻率為中心頻率,以故障齒輪所在軸的轉頻為調頻的邊頻帶。圖10(b)中,出現了沖擊頻率13.2 Hz的1倍頻、2倍頻及3倍頻,說明此時發生了剛度變化而引得的調幅調頻同時存在。這種現象從圖10(c)也可以看出,由于調幅調頻同時存在,信號的邊頻帶不再對稱,嚙合頻率1 186 Hz的振幅也不再是最大。

6 結束語

通過將集成經驗模態分解(EEMD)與時頻峰值濾波法(TFPF)有效地結合,突出了兩種方法各自的優點,使TFPF降噪僅作用于含噪聲成分較多的IMFs,而不是在整個信號,解決了TFPF的窗長選擇在信號幅度和噪聲壓制上的矛盾,突破了TFPF方法窗長選擇的局限性,提高了分析的準確性。對降噪后的信號進行循環自相關解調分析,有效地提取了齒輪齒根裂紋的故障特征。EEMD+TFPF除了適用齒輪傳動系統故障診斷,還可以用于其他系統的故障診斷。

[1] 丁康,李巍華.齒輪及齒輪箱故障診斷實用技術[M].北京:機械工業出版社,2005:63-65.

[2] Wu Zhaohua, Huang Norden E. Ensemble empirical mode decomposition: A noise assisted data analysis method [J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41.

[3] 陳仁祥,湯寶平,呂中亮.基于相關系數的EEMD轉子振動信號降噪方法[J].振動、測試與診斷,2012,32(4):542-546.

Chen Renxiang, Tang Baoping, Lu Zhongliang. Ensemble empirical mode decomposition de-noising method based on correlation coefficients for vibration signal of rotor system[J]. Journal of Vibration, Measurement & Diagnosis, 2012,32(4):542-546. (in Chinese)

[4] 秦娜,蔣鵬,孫永奎,等. 基于EEMD排列熵的高速列車轉向架故障特征分析[J].振動、測試與診斷,2015,35(5):885-891.

Qin Na,Jiang Peng,Sun Yongkui, et al. Fault diagnosis of high speed train bogie based on EEMD and permutation entropy [J]. Journal of Vibration, Measurement & Diagnosis, 2015,35(5):885-891. (in Chinese)

[5] 鄭直,姜萬錄,胡浩松.基于EEMD形態譜和KFCM聚類集成的滾動軸承故障診斷方法研究[J].振動工程學報,2015,28(2):324-330.

Zheng Zhi,Jiang Wanlu,Hu Haosong. Research on rolling bearings fault diagnosis method based on EEMD morphological spectrum and kernel fuzzy C-means clustering [J]. Journal of Vibration Engineering, 2015,28(2):324-330. (in Chinese)

[6] 王 鵬,王太勇,董靖川.基于EEMD時頻譜二值化的振動信號微弱特征提取方法[J].天津大學學報:自然科學與工程技術版,2016,49(7):667-673.

Wang Peng,Wang Taiyong,Dong Jingchuan. Weak feature extraction of vibration signal based on binaryzation of EEMD time-frequency [J].Journal of Tianjin University:Science and Technology,2016,49(7):667-673. (in Chinese)

[7] 蔣超,劉樹林,姜銳紅,等.基于快速譜峭度圖的EEMD內稟模態分量選取方法[J].振動、測試與診斷,2015,35(6):1173-1178.

Jiang Chao,Liu Shulin,Jiang Ruihong, et.al. Feature extraction method of intrinsic mode function in EEMD based on fast kurtogram in machinery fault diagnosis [J]. Journal of Vibration, Measurement & Diagnosis, 2015,35(6):1173-1178. (in Chinese)

[8] 王志堅,韓振南,劉邱祖,等.基于MED-EEMD的滾動軸承微弱故障特征提取[J].農業工程學報.2014,30(23):70-78.

Wang Zhijian, Han Zhennan, Liu Qiuzu, et al. Weak fault diagnosis for rolling element bearing based on MED-EEMD[J]. Transactions of the Chinese Society of Agricultural Engineering, 2014, 30(23): 70-78. (in Chinese)

[9] Boashash B, Mesbah M. Signal enhancement by time frequency peak filtering[J] .IEEE Transactions on Signal Processing,2004, 52(4):929-937.

[10] Liu Yanping, Li Yue. Spatiotemporal time-frequency peak filtering method for seismic random noise reduction[J]. IEEE Geoscience and Remote Sensing Letters,2013,10(4):756-760.

[11] Zhang Jie, Li Yue,Wu Ning. Noise attenuation for seismic data by Hyperbolic-Trace time-frequency peak filtering[J]. IEEE Geoscience and Remote Sensing Letters, 2015,12(6):601-605.

[12] Zhang Chao , Li Yue, Lin Hongbo. Signal preserving and seismic random noise attenuation by Hurst exponent based time-frequency peak filtering[J].Geophysical Journal International, 2015,203(2):901-909.

[13] 楊平.基于時頻消噪 TFPF和時頻分布 MBD的軸承早期故障診斷[J]. 四川理工學院學報,2010,23(3):357-360.

Yang Ping. Rolling bearing incipient fault diagnosis based on TFPF and modified b-distribution[J]. Journal of Sichuan University of Science & Engineering,2010,23(3):357-360. (in Chinese)

[14] 劉彥萍.時空二維時頻峰值濾波方法壓制地震勘探隨機噪聲的研究[D].長春:吉林大學,2013.

國家自然科學基金資助項目(50775157);山西省基礎研究資助項目(2012011012-1)

2016-06-23;

2016-09-21

TH17; TH13

寧少慧,女,1978年10月生,博士、講師。主要研究方向為齒輪傳動系統故障診斷。曾發表《A novel fault diagnosis approach of gearbox using an embedded sensor fixed gear body》(《Journal of Vibroengineering》2016,Vol.18,No.7)等論文。

E-mail:nshzzl@126.com

韓振南,男,1958年2月生,教授、博士生導師。主要研究方向為機械傳動系統故障診斷。

E-mail:zhennan-han@hotmail.com

猜你喜歡
故障診斷模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
國內多模態教學研究回顧與展望
因果圖定性分析法及其在故障診斷中的應用
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 日韩毛片免费| 五月婷婷精品| 亚洲热线99精品视频| 亚洲中文字幕无码爆乳| 国产午夜无码片在线观看网站| 伊人久久大香线蕉影院| 国产人人射| 色综合久久久久8天国| 欧美日韩综合网| 亚洲精品无码在线播放网站| 国产网站黄| 性网站在线观看| 伊人久久综在合线亚洲2019| 国产jizz| 国产第一页亚洲| 内射人妻无套中出无码| 喷潮白浆直流在线播放| 成人日韩视频| 丁香婷婷久久| 玖玖精品视频在线观看| 日韩午夜片| 欧美日韩激情在线| 亚洲91在线精品| 天堂岛国av无码免费无禁网站| 中文无码精品A∨在线观看不卡| 免费看美女自慰的网站| 午夜国产小视频| 日韩在线中文| 国产成人啪视频一区二区三区| 91美女视频在线| 在线精品视频成人网| 青青草国产在线视频| 欧美成人手机在线观看网址| 九色综合视频网| 国产剧情无码视频在线观看| 人妻无码一区二区视频| 一级毛片网| h视频在线观看网站| 欧美亚洲中文精品三区| 国产免费精彩视频| 亚洲无码高清视频在线观看| 国产91麻豆免费观看| 欧美一级一级做性视频| 国产资源站| 欧美区一区二区三| 欧美日韩在线成人| 国产精品亚洲а∨天堂免下载| 精品国产网站| 国产在线97| 国产精选自拍| 免费一级成人毛片| 伦伦影院精品一区| 丁香六月激情综合| 国产激情无码一区二区APP | 国产精品自在在线午夜| 97视频在线观看免费视频| 国产亚洲精品在天天在线麻豆| 国产精品嫩草影院av| 免费三A级毛片视频| 91视频免费观看网站| 亚洲中文精品人人永久免费| 老司机午夜精品网站在线观看| 高清不卡一区二区三区香蕉| 成人精品在线观看| 亚洲日韩第九十九页| 亚洲精品成人片在线观看| 老司机久久99久久精品播放| 制服丝袜一区| 久草热视频在线| 国产乱肥老妇精品视频| 99视频全部免费| 欧美国产在线精品17p| 亚洲色图另类| 国产福利一区视频| 久久香蕉欧美精品| www成人国产在线观看网站| 欧美一区二区啪啪| 国产福利免费视频| 在线观看国产小视频| 综合久久五月天| 午夜三级在线| 国产人人射|