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

VMD和ICA聯合降噪方法在軸承故障診斷中的應用

2017-07-18 11:49:23馬增強柳曉云張俊甲王建東
振動與沖擊 2017年13期
關鍵詞:故障診斷模態故障

馬增強, 柳曉云, 張俊甲, 王建東

(石家莊鐵道大學 電氣與電子工程學院,石家莊 050043)

VMD和ICA聯合降噪方法在軸承故障診斷中的應用

馬增強, 柳曉云, 張俊甲, 王建東

(石家莊鐵道大學 電氣與電子工程學院,石家莊 050043)

針對振動信號易受噪聲干擾的影響、故障特征提取困難的問題,提出一種基于變分模態分解(Variational Mode Decomposition, VMD)和獨立分量分析(Independent Component Analysis, ICA)相結合的去噪方法。該方法首先利用VMD算法將振動信號分解成若干不同頻率的本征模態分量(Intrinsic Mode Function, IMF),有效的抑制了LMD分解中存在的模態混疊現象和端點效應等問題,然后依據峭度準則選取相應分量進行重構,引入虛擬噪聲通道;最后利用FastICA將重構后信號再次進行去噪處理,分離出有效的故障特征分量,從而識別故障類型。將該方法應用到滾動軸承故障數據中,并與LMD-ICA方法作對比,結果表明,提出方法不僅能夠有效的解決去噪過程中丟失故障信息以及由于模態混疊導致噪聲不能完全去除的問題,還能更清晰、準確地提取出故障特征頻率。

變分模態分解;獨立分量分析;降噪;滾動軸承;故障診斷

滾動軸承是旋轉機械中不可或缺的部件之一,據統計,在旋轉機械中,約有30% 的機械故障由滾動軸承引起,滾動軸承質量的好壞對機械設備工作狀況有很大影響。因此,對滾動軸承的故障診斷具有重要的意義。

機械系統中,由于受復雜背景噪聲以及其它干擾源影響,導致故障特征難以精確地提取出來。一些學者在振動信號降噪方面做了大量研究。Wu等[1]提出了經驗模態分解(Empirical Mode Decomposition,EMD)算法,這是目前廣泛使用的一種自適應信號處理方法,把非平穩信號分解成不同頻段的模態分量進而轉化成平穩信號進行分析。蘇文勝等[2]將EMD應用到滾動軸承早期故障診斷中,實現了軸承故障診斷。吳小濤等[3]針對 EMD 的模態混疊,用集成經驗模態分解(Ensemble Empirical Mode Decomposition, EEMD)清晰地提取到故障特征信息,該方法可自適應地實現非平穩振動信號降噪。

Lin[4]提出了局部均值分解(Local Mean Decomposition,LMD)。陳亞農等[5]將LMD算法應用到滾動軸承故障診斷中,并提取出故障特征頻率。鞠萍華等[6]通過LMD得到若干個PR(Production Functions)分量,并對每一個分量進行能量算子解調,從而進行故障診斷,并指出了該方法優于EEMD診斷方法。

EEMD或者LMD被廣泛地應用在故障特征提取方面,然而兩者都屬于遞歸“篩選”模態,存在端點效應和模態混疊現象,且受采樣頻率影響,分解誤差較大。文獻[7]提出一種自適應信號處理新方法-變分模態分解(Variational Mode Decomposition,VMD),該方法通過迭代搜尋變分模型最優解從而確定每個分解分量的頻率中心和帶寬,可以自適應地實現信號的頻域剖分和各分量的有效分離。相比EEMD和LMD的遞歸“篩選”模態,VMD將信號分解轉化為非遞歸、變分模態分解模態,表現出更好的噪聲魯棒性。

獨立分量分析(Independent Component Analysis, ICA)是一種基于樣本高階統計信息的新的信號處理方法[8-10],目前已廣泛應用在通信信號處理、語音信號降噪等各個領域。為了進一步提高降噪效果,在以上研究的基礎上,本文提出一種基于VMD與ICA聯合降噪的方法,并將該方法應用在滾動軸承的故障診斷中。通過對實測數據分析結果表明,本文提出方法能夠有效地實現特征信號提取,相比其他方法降噪效果更好。

1 VMD-ICA聯合降噪方法基本原理

1.1 變分模態分解原理

VMD是一種新的信號分解估計方法,是基于經典維納濾波、希爾伯特變換和混頻的變分問題求解過程,該方法通過迭代搜尋變分模型最優解從而確定每個分解分量的頻率中心和帶寬,可以自適應地將信號分解成具有稀疏特性的分量[11-13]。

假設每個模態是具有中心頻率的有限帶寬,中心頻率和帶寬在分解過程中不斷更新,VMD分解是尋求K個估計帶寬之和最小的模態函數uk(t),k∈{1,2,…,K},各模態之和為輸入信號f。通過以下步驟確定每個模態函數的帶寬:

(1) 為了獲得模態函數的解析信號,對每個模態函數uk(t)進行希爾伯特變換,即

(1)

式中:t為時間,是大于0的正數;δ(t)為沖擊函數;{uk(t)}={u1(t),u2(t),…,uK(t)}為分解得到的K個IMF分量。

(2) 對各模態解析信號預估中心頻率e-jωkt進行混合,將每個模態的頻譜調制到相應的基頻帶,即

(2)

式中,{ωk}={ω1,…,ωK}為各IMF分量uk(t)的中心頻率。

(3) 計算以上解調信號的梯度的平方L2范數,估計出各模態分量的帶寬。對應的約束變分模型表達式為

(3)

為了求取上述約束變分問題,引入二次懲罰因子α與Lagrange乘法算子λ(t),其中α為足夠大的正數,可在高斯噪聲存在的情況下保證信號的重構精度,Lagrange算子使得約束條件保持嚴格性,擴展的Lagrange表達式為

L({uk(t)},{ωk},λ(t))=

(4)

(5)

式中,i∈{1,2,…,K}且i≠k。利用傅里葉等距變換,將式(5)轉變到頻域后用ω-ωk代替ω,并將得到的結果轉換為負頻率區間積分的形式,則優化問題的解為

k∈{1,2,…,K}

(6)

根據同樣的過程,求得中心頻率更新方式

(7)

(2) 執行循環n=n+1;

(8)

式中,τ為時間常數,常取為0。

(5) 給定判別精度ε>0,重復上述步驟,直到滿足迭代停止條件

(9)

1.2 峭度準則

峭度是反映隨機變量分布特性的數值統計量,是歸一化的4階中心矩。

(10)

式中:μ、σ分別為信號的均值和標準差;E(x)為變量x的期望值。峭度系數對軸承早期故障比較敏感,當軸承沒有任何故障時振動信號近似服從正態分布,而正態分布的峭度系數約等于3;當軸承逐步發生故障,工作表面出現損傷引起機械沖擊時,其振動信號概率密度偏離正態分布,峭度系數值增加較快。所以可以推斷出,峭度值越大說明這些IMF分量中包含有更多的沖擊成分,對這些分量進行重構,則合成的信號峭度值會明顯提高,包含的故障特征信息越多、故障特征越明顯。

1.3 ICA算法基本原理

獨立分量分析是一種盲信源分離技術。該方法是指在先驗信息很少的情況下,僅通過觀測信號來估計有用源信號的一種新型信號處理方法。模型為

假設x(t)=(x1(t),x2(t),…,xm(t))T為m個觀測信號,s(t)=(s1(t),s2(t),…,sn(t))T為n個未知的相互獨立的源信號,x(t)中各分量可由s(t)各獨立信源線性表示,用矩陣可表示為

x(t)=As(t)

(11)

式中:x(t)為觀測矩陣;s(t)為信號源矩陣;A為一個未知的滿秩的m×n的混合系數矩陣(m≥n)。ICA的目的是要在信號源矩陣s(t)和系數矩陣A均未知的情況下,利用統計學途徑獲得分離矩陣W,使得分離矩陣滿足

y(t)=Wx(t)=WAs(t)

(12)

式中,y(t)為源信號s(t)的估計。獨立分量分析是最新提出的算法,其目的是使非高斯數據可以被線性表示。從線性變換角度來看,源信號即為相互獨立的非高斯信號,也可稱之為線性空間的基信號,而觀測信號則是源信號的線性表示。ICA就是在源信號和系數矩陣A均未知的情況下,從混合信號中估計出有用的輸入信號。該算法在統計獨立的意義下對混合信號進行分離,從而提取出混合信號中的有用信號,相比傳統的自適應濾波方法效果更好。

2 VMD-ICA聯合降噪方法研究

由于軸承故障信號具有能量小、頻帶分布寬、易受其它震源干擾等問題,單獨使用VMD算法對信號進行處理可能導致故障特征信息不能被完全提取出來。而獨立分量分析目前的研究大都是基于多觀測通道(即要求振動傳感器的個數大于等于振動源信號。本文首先將故障信號進行VMD分解得到若干IMFs分量,然后依據峭度準則將IMF分量重構,引入虛擬噪聲通道,由于引入信號包含原始信號先驗信息,所以避免了盲目引入噪聲信號導致處理結果不佳的影響。接下來對重構信號進行FastICA解混,最后對降噪后的信號進行能量算子解調,并將理論值與解調后的頻率進行對比,識別故障類型。算法流程如圖1所示,主要步驟為:

圖1 VMD-ICA聯合降噪算法流程圖

步驟1 獲取振動信號,初始化模態數K=2,懲罰因子α和帶寬τ使用默認值,α=2 000,τ=0。

步驟3 判斷中心頻率是否相近,如果相近,則確定模態數為K=K-1,否則以模態數K=K+1繼續進行步驟2。

步驟4 計算各個IMF分量的峭度值,依據峭度準則將IMF分量進行分組重構,得到觀測信號s1(t),并構造虛擬通道信號s2(t)。

步驟5 利用FastICA算法對步驟4得到的觀測信號s1(t)和虛擬通道信號s2(t)進行解混,得到聯合降噪后信號y1(t)。

步驟6 對步驟5得到的信號進行能量算子解調,提取故障特征頻率,并與軸承特征頻率的理論值進行對比,從而識別故障類型。

3 實驗分析

為了驗證本文提出方法在滾動軸承故障特征提取中的有效性,采用實際軸承內圈故障信號進行驗證,實驗平臺如圖2所示的QPZZ-Ⅱ旋轉機械故障試驗臺。本實驗測試參數及軸承技術參數如表1所示,由軸承特征頻率理論計算公式可得到內圈故障頻率f=37.5 Hz。滾動軸承內圈故障振動信號時域及頻域波形如圖3所示。

表1 軸承測試參數與技術參數

圖2 QPZZ-Ⅱ旋轉機械故障試驗臺

(a)內圈故障振動信號波形(b)內圈故障振動信號頻譜

圖3 內圈故障信號時域波形及頻譜

Fig.3 Waveform and spectrum of the bearing inner fault

本文將從故障信號分解效果方面、信號去噪方面以及故障信號頻率檢測精度方面進行分析,并通過和LMD-ICA方法進行對比說明本文方法的有效性和優越性。

3.1 故障信號分解效果對比

為了防止出現欠分解和過分解現象,根據不同的K值對應的各個模態中心頻率來確定分解次數。由表2可知當K=5時中心頻率相差較大,出現了欠分解現象;而K=7時出現了中心頻率相近的模態分量,發生了過分解現象。因此最終選取模態數K=6,分解結果如圖4(a)所示。為了進行比較,將同一故障信號進行了LMD處理,結果如圖4(b)所示。

表2 不同K值各個模態分量的中心頻率

(a) VMD模態分量波形及頻譜

(b) LMD模態分量波形及頻譜

圖4 軸承故障信號分解結果

Fig.4 Decomposition results of the original signal

從分解時域圖中可以看出,VMD分解在一定程度上克服了LMD中存在的模態混疊和端點效應問題。而且從頻譜圖中可以看出,LMD的頻段分離不是很分明,尤其是第一頻譜分量,在整個頻率范圍內幾乎都有分量,而VMD的各個頻段的分離效果則較好。

3.2 重構信號去噪效果對比

3.2.1 實測信號分析

計算VMD分解后各個IMF分量的峭度值,如表3所示。由表3可知,IMF1、IMF2峭度值較大,說明其保留了原始信號中最多的沖擊特征,故提取這兩個IMF分量進行原始信號重構,剩余分量重構作為虛擬噪聲通道,由于噪聲通道的構成包含了原始信號本身部分的先驗信息,避免了盲目選擇噪聲信號引起處理效果不佳的影響。重構信號如圖5(a)、圖5(b)所示,將重構信號利用FastICA方法進行分離,分離結果如圖5(c)、圖5(d)所示。將同一故障信號進行LMD分解后,計算各個分量的峭度值如表4所示,選取前三個分量進行重構,并利用FastICA方法進行分離,圖6(a)~圖6(d)為用LMD-ICA方法對信號進行處理得到的結果。

表3 IMFs分量峭度值

(a)重構故障信號(b)重構噪聲信號(c)獨立分量IC1(d)獨立分量IC2

圖5 本文所提方法分析結果

Fig.5 Analysis results with the proposed method

表4 IMFs分量峭度值

圖5、圖6定性地反映了以上兩種方法的去噪效果,從圖中可以看出,FastICA算法可有效地從混合信號中提取出有用信號,將噪聲信號和有用信號分離開,由于ICA沒有任何關于原始信號的先驗知識,具有不確定性,所以解混后信號幅值可能不同,但是并不影響對信號的分析。

(a)重構故障信號(b)重構噪聲信號(c)獨立分量IC1(d)獨立分量IC2

圖6 LMD-ICA算法分析結果

Fig.6 Analysis results with LMD-ICA

為了定量地分析重構故障信號在以上兩種方法降噪后的去噪效果,選取峭度值(K)、峰值信噪比(PSNR)作為降噪后效果的評價指標。峭度值說明了重構信號中包含的故障特征信息,峭度值越大則說明信號中故障信息越多;峰值信噪比則反映了信號的去噪能力,峰值信噪比越大則說明去噪效果越好。

表5 降噪結果對比

3.2.2 仿真信號分析

在以上研究的基礎上,構造如下仿真信號。其中,幅值A=1;衰減系數K=800;系統共振頻率ωr=2×pi×1 000;u(t)為單位階躍函數;設定信號的故障特征頻率fr=128;滾珠和滾道之間微小滑動對故障特征頻率的影響因子τi為0.01/fr-0.02/fr之間的隨機數,采樣頻率為25 600 Hz,n(t)為白噪聲。故障仿真信號的波形及頻譜如圖 7(a)、圖7(b)所示。

(13)

相關系數反映了去噪后信號與原始信號之間整體波形的相似程度,在仿真信號中加入不同的信噪比,用相似度指標進一步比較以上兩種方法的去噪性能,曲線圖如圖8所示。從圖中可以看出,在不同信噪比情況下本文提出方法得到的去噪信號與原始信號相似度更高。尤其是在信噪比較低的情況下,本文提出方法仍能有效地降低噪聲的干擾能力,在去噪方面效果更好。

(a)故障仿真信號波形(b)故障仿真信號頻譜

圖7 仿真信號波形及頻譜

Fig.7 Waveform and spectrum of simulated signals

圖8 信噪比-相似度曲線圖

3.3 故障信號頻率檢測精度對比

將降噪后的信號分別進行能量算子解調得到相應能量譜,如圖9(a)、圖9(b)所示。圖10(a)、圖10(b)為用LMD-ICA方法得到的能量譜。依據文獻[14],本文定義故障頻率檢測精度為

(a)獨立分量IC1能量譜(b)獨立分量IC2能量譜

圖9 本文所提方法分析結果

圖10 LMD-ICA分析結果

Fig.10 The analysis results of LMD-ICA

4 結 論

振動信號中常受到噪聲的影響,這些干擾導致故障特征提取困難,基于此本文提出了基于VMD和FastICA聯合降噪的方法,并應用于軸承故障特征提取中,通過實測數據驗證了提出方法的實用性和有效性,結論如下:

(1) 通過VMD將信號進行分解,引入了虛擬噪聲通道,既解決了盲目引入噪聲導致處理結果不佳的影響,又解決了ICA算法不能對單通道信號進行處理的問題。

(2) 通過峭度值選擇相應分量重構,避免了盲目選擇分量,造成故障特征信息丟失的弊端,同時對比了將重構信號進行ICA解混后的降噪效果,驗證了本文提出方法在降噪方面的優越性。

(3) 相較于LMD-ICA方法,本文提出方法解決了去噪過程中丟失故障信息以及由于模態混疊導致噪聲不能完全去除的問題。

[1] WU Z H, HUANG N E. Ensemble empirical mode decomposition: a noise-assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41.

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

SU Wensheng, WANG Fengtao, ZHANG Zhixin, et al. Application of EMD denoising and spectral kurtosis in early fault diagnosis of rolling element bearings[J]. Journal of Vibration and Shock, 2010, 29(3): 18-21+201.

[3] 吳小濤,楊錳,袁曉輝,等. 基于峭度準則EEMD及改進形態濾波方法的軸承故障診斷[J]. 振動與沖擊,2015,34(2):38-44.

WU Xiaotao, YANG Meng, YUAN Xiaohui, et al. Bearing fault diagnosis using EEMD and improved morphological filtering method based on kurtosis criterion[J]. Journal of Vibration and Shock, 2015,34(2):38-44.

[4] LIN Jinshan. Fault diagnosis of gearboxes based on the local mean decomposition method[J]. Intelligent Information Technology Application Association, 2011,141:559-563.

[5] 陳亞農,郜普剛,何田,等.局部均值分解在滾動軸承故障綜合診斷中的應用[J].振動與沖擊,2012,31(3):73-78.

CHEN Yanong, GAO Pugang, HE Tian, et al. Roller bearing comprehensive fault diagnosis based on LMD[J]. Journal of Vibration and Shock, 2012, 31(3): 73-78.

[6] 鞠萍華,秦樹人,趙玲.基于LMD的能量算子解調方法及其在故障特征信號提取中的應用[J].振動與沖擊,2011,30(2):1-4.

JU Pinghua, QIN Shuren, ZHAO Ling. Energy operator demodulating approach based on LMD and its application in extracting characteristics of a fault signal[J]. Journal of Vibration and Shock, 2011, 30(2): 1-4.

[7] DRAGOMIRETSKIY K, ZOSSO D. Variational mode decomposition[J]. IEEE Transactions on Signal Processing, 2013, 62(3): 531-544.

[8] HYVARINEN A, KARHUNEN J, OJA E. Independent component analysis[M]. New York: John Wiley, 2001.

[9] 張俊紅,李林潔,馬文朋,等. EMD-ICA聯合降噪在滾動軸承故障診斷中的應用[J]. 中國機械工程,2013,24(11):1468-1472.

ZHANG Junhong, LI Linjie, MA Wenpeng, et al. Application of EMD-ICA to fault diagnosis of rolling bearings[J]. China Mechanical Engineering,2013,24(11):1468-1472.

[10] 姚家馳,向陽,李勝楊,等. 基于VMD-ICA-CWT的內燃機噪聲源識別方法[J]. 華中科技大學學報(自然科學版),2016,44(7):20-24.

YAO Jiachi, XIANG Yang, LI Shengyang, et al. Noise sources identification method of internal combustion engine based on VMD-ICA-CWT[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition),2016,44(7):20-24.

[11] 武英杰,甄成剛,劉長良.變分模態分解在風電機組故障診斷中的應用[J].機械傳動,2015,39(10):129-132.

WU Yingjie, ZHEN Chenggang, LIU Changliang. Application of variational mode decomposition in wind power fault diagnosis[J].Journal of Mechanical, 2015, 39(10): 129-132.

[12] 趙洪山,郭雙偉,高奪.基于奇異值分解和變分模態分解的軸承故障特征提取[J].振動與沖擊,2016,35(22):183-188.

ZHAO Hongshan, GUO Shuangwei, GAO Duo. Fault feature extraction of bearing faults based on singlular value decomposition and variational modal decomposition[J].Journal of Vibration and Shock,2016,35(22):183-188.

[13] 唐貴基,王曉龍.參數優化變分模態分解方法在滾動軸承早期故障診斷中的應用[J].西安交通大學學報,2015,49(5):73-81.

TANG Guiji, WANG Xiaolong. Parameter optimized variational mode decomposition method with application to incipient fault diagnosis of rolling bearing[J]. Journal of Xi’an Jiaotong University, 2015, 49(5): 73-81.

[14] 熊偉,呂英華. 視頻紅信號頻譜識別中波峰搜索算法研究[J]. 軟件,2012,11:190-193.

XIONG Wei, Lü Yinghua. Peak finding in the identify of video red signal spectrum[J]. Software,2012, 11: 190-193.

Application of VMD-ICA combined method in fault diagnosis of rolling bearings

MA Zengqiang, LIU Xiaoyun, ZHANG Junjia, WANG Jiandong

(Electrical and Electronics Engineering, Shijiazhuang Railway University, Shijiazhuang 050043, China)

In order to solve the problem of fault features of rolling bearings being difficult to extract, the signal analysis method of variational mode decomposition(VMD) combined with independent component analysis(ICA) was proposed here. At first, VMD was used to decompose a multi-component vibration signal into a number of quasi-orthogonal intrinsic mode functions (IMFs) to effectively suppress problems, such as, mode mixing, end effect and so on existing in LMD algorithm. Then the kurtosis criterion was used to choose the corresponding IMFs to be reconstructed, and to induce a virtual noise channel. Finally, the reconstructed signal was denoised again with the fastICA, the effective fault feature components were extracted and the fault types were identified. In order to verify the effectiveness of the proposed method, this method was applied in fault diagnosis of rolling bearings, its effect was compared with that of LMD-ICA algorithm. The results demonstrated that the proposed method can not only solve problems of losing fault information and not enough eliminating noise due to mode mixing in denoising process, but also extract fault feature frequencies more clearly and correctly.

rolling bearing; variational mode decomposition (VMD); noise reduction; independent component analysis (ICA); fault diagnosis

國家自然科學基金(11227201;11372199;11572206);河北省自然科學基金(A2014210142)

2016-07-18 修改稿收到日期:2016-10-17

馬增強 男,博士,教授,1975年生

柳曉云 女,碩士生,1991年生

TH165+.3

A

10.13465/j.cnki.jvs.2017.13.032

猜你喜歡
故障診斷模態故障
故障一點通
奔馳R320車ABS、ESP故障燈異常點亮
國內多模態教學研究回顧與展望
因果圖定性分析法及其在故障診斷中的應用
故障一點通
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 亚洲国产一区在线观看| 亚洲精品中文字幕无乱码| 好紧好深好大乳无码中文字幕| 亚洲天堂首页| 中国一级毛片免费观看| 99re视频在线| 日韩精品视频久久| 精品久久久久成人码免费动漫| 成人午夜久久| 国产一级在线播放| 无码视频国产精品一区二区| AV无码无在线观看免费| 亚洲人妖在线| 亚洲IV视频免费在线光看| 欧美啪啪视频免码| 老司机久久99久久精品播放| 91无码人妻精品一区二区蜜桃| 午夜国产小视频| 国产香蕉一区二区在线网站| 亚洲日韩每日更新| 欧美日韩北条麻妃一区二区| 国产又粗又爽视频| 日韩av在线直播| 欧美成人一级| 一区二区三区高清视频国产女人| 亚洲一级毛片在线观播放| 亚洲天堂福利视频| 精品国产福利在线| 亚洲天堂在线视频| 国产麻豆aⅴ精品无码| 农村乱人伦一区二区| 亚洲天堂色色人体| 日韩在线视频网站| 欧美性精品| 国产一区二区人大臿蕉香蕉| 538国产视频| 国产精品女人呻吟在线观看| 久久99久久无码毛片一区二区| 视频二区中文无码| 国禁国产you女视频网站| 免费人成视网站在线不卡| 91视频99| 亚洲精品视频免费看| 久草中文网| 国产成人精品一区二区三区| 国产精品污污在线观看网站| 亚洲日本一本dvd高清| 婷婷亚洲天堂| 欧美在线中文字幕| 精品亚洲国产成人AV| 久久国产乱子| 国产精品尤物在线| 国产精品所毛片视频| 国产精品视频久| 99久久免费精品特色大片| 五月激激激综合网色播免费| 国产精品无码影视久久久久久久| 黄色国产在线| 国产成人AV男人的天堂| 日韩中文无码av超清| 狠狠色综合久久狠狠色综合| 原味小视频在线www国产| 欧美一级爱操视频| 色噜噜狠狠狠综合曰曰曰| 久久91精品牛牛| 国产成人精品高清不卡在线| 人妻丰满熟妇AV无码区| 久久精品亚洲中文字幕乱码| 中文字幕欧美日韩| 亚洲va在线观看| 国产浮力第一页永久地址| 国产一线在线| 亚洲成a人在线播放www| 综合人妻久久一区二区精品 | 99在线视频网站| 亚洲人在线| 国产精品久久久久久久久| 日韩免费毛片| 狼友av永久网站免费观看| 男人天堂亚洲天堂| 国产日韩精品欧美一区喷| 毛片基地美国正在播放亚洲 |