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

基于奇異值分解和變分模態分解的軸承故障特征提取

2016-12-12 11:34:25趙洪山郭雙偉
振動與沖擊 2016年22期
關鍵詞:特征提取模態振動

趙洪山, 郭雙偉, 高 奪

(華北電力大學(保定) 電氣與電子工程學院,河北 保定 071003)

?

基于奇異值分解和變分模態分解的軸承故障特征提取

趙洪山, 郭雙偉, 高 奪

(華北電力大學(保定) 電氣與電子工程學院,河北 保定 071003)

為了有效提取軸承故障,提出了基于變分模態分解和奇異值分解降噪的故障特征提取方法。通過對故障信號進行變分模態分解,獲得其本征模態函數。基于峭度指標,選擇包含故障信息的本征模態函數進行信號重構。利用奇異值分解降噪技術對重構信號進行處理,提高信噪比。最后對降噪信號進行包絡解調提取故障特征頻率。與常見的故障特征提取方法相比,該方法能有效辨別滾動軸承的典型故障,突出故障特征,提高滾動軸承的故障診斷效果。

變分模態分解;奇異值分解;滾動軸承;故障特征提取

軸承是機械傳動系統的核心部件[1],在功率傳遞的過程中發揮著至關重要的作用,如新能源風力發電機主軸承、齒輪箱的各級軸承等。一旦軸承發生故障,傳動系統的正常運行會受到極大影響。在軸承的故障中,元件表面損傷最為常見。軸承元件發生表面損傷后,會和與之配合的元件表面發生撞擊,產生頻帶很寬的脈沖力,覆蓋軸承的固有頻率,從而激起系統的高頻固有振動,導致軸承故障信號的頻譜產生多個共振峰,共振峰中包含了故障信息。

共振解調法[2]通過對振動信號進行窄帶濾波得到包含故障信息的故障頻帶,利用包絡解調進行故障特征提取。該方法能有效提取滾動軸承的局部故障特征,但若共振頻段選擇不合理,會造成故障特征無法提取,進而影響故障診斷結果。小波變換[3]能夠對振動信號進行分解,得到每一頻帶內振動信號的變化規律。該方法具有多分辨率的優點,可以由粗到細地逐步觀察信號,從而實現故障特征的有效提取。但基函數與閾值的選擇對檢測結果有較大影響,對設計者的要求也較高。以Wingner-Ville分布為基礎的雙線性時頻分析方法也被應用到故障特征的提取中[4-5],該方法提供了時間域與頻率域的聯合分布信息,清楚地描述了信號頻率隨時間變化的關系,具有很好的頻譜分辨率。但由于交叉項的存在,容易造成信息丟失與頻率混疊,不利于對瞬態信號的探測。經驗模態分解[6](Empirical Mode Decomposition, EMD)可以處理非線性、非平穩信號,能自適合地將非平穩信號分解為若干平穩的本征模態函數(IMF),克服了小波變換和自適應時頻分析方法的不足,具有自適應、正交性和完備性的特點,在故障特征提取方面受到廣泛關注[7-8]。然而,EMD方法本身存在著一些不足,如模態混疊[9]、存在端點效應、受采樣頻率影響較大等[10]。對此,DRAGOMIRETSKIY等[11]提出一種自適應信號處理新方法—變分模態分解(Variational Mode Decomposition, VMD),既保留了EMD對非平穩信號自適應分解的優點,又彌補了EMD方法的不足。

奇異值分解(Singular Value Decomposition,SVD)本身具有極好的穩定性和不變性。通過構造信號的Hankel矩陣,并對矩陣進行奇異值分解,選取恰當的奇異值進行信號重構能有效地消除信號中的隨機成分,最大限度地保留有用信息,剔除無用信息,提高信號的信噪比。

本文結合變分模態分解(VMD)和奇異值分解提出一種新的軸承故障特征提取方法。故障信號經VMD分解為若干IMF,基于峭度指標,選取包含故障信息的IMF進行信號重構。利用SVD對重構信號進行降噪處理,最后對降噪信號進行包絡解調,提取軸承的故障特征。

1 變分模態分解(VMD)

變分模態分解 (VMD)是一種新的信號分解估計方法。在VMD算法中,每一個本征模態函數(IMF)均被視為調幅-調頻信號uk(t):

uk(t)=Ak(t)cos(φk(t))

(1)

VMD的分解過程是一個變分問題的求解過程。假設每個IMF具有有限帶寬,變分問題可表示為尋求k個模態函數uk(t),使得所有模態函數的估計帶寬之和最小,并且滿足各模態之和等于原始輸入信號f的約束條件。具體的分解步驟如下:

(1)對每一個模態分量信號uk(t),利用Hilbert變換計算與之對應的解析信號,得到其單邊頻:

(δ(t)+j/πt)·uk(t)

(2)

(2)加入指數項e-jωkt調整每一種模態函數對應解析信號的預估中心頻率,將每個模態的頻譜轉移至基帶:

[(δ(t)+j/πt)·uk(t)]e-jωkt

(3)

(3)利用H1高斯平滑估計移頻后解析信號的帶寬,得到受約束的變分問題如下:

s.t. ∑uk=f

(4)

式中:{uk}:{u1,u2,…uk},{ωk}:={ω1,ω2,…ωk}

(4)引入二次懲罰因子α和拉格朗日乘子算子λ(t)構造擴展拉格朗日表達式Γ(uk,ωk,λ):

Γ(uk,ωk,λ)=α∑k‖?t[(δ(t)+

(5)

通過反復迭代,尋找擴展拉格朗日表達式的“鞍點”求解最小值,獲得最優解。最優解為本征模態函數{uk}及各自的中心頻率{ωk}。

2 奇異值分解(SVD)

對于離散數字信號序列{hi,i=1,2,…,N},根據相空間重構理論[12]可得到L×K階的Hankel矩陣:

(6)

式中:N為信號長度,K=N-L+1,延時值為1,H為軌道矩陣。

軌道矩陣H經奇異值分解得到:

H=USVT

(7)

式中:U、V為正交矩陣:VT為V的轉置;S為對角陣,S=diag(σ1,σ2,…σn)σ1,σ2,…σn為奇異值,V∈RK×K。

由奇異值理論知[13]:信號中的有用分量對應前k個奇異值,噪聲分量對應后面較小的奇異值。利用前k個奇異值進行矩陣重構可以在Forbeious范數意義下實現對H的最佳逼近,從而降低噪聲,提高信噪比。

奇異值的差分譜為相鄰兩奇異值的差值:

di=σi-σi+1,(i=1,2,…,t)

(8)

式中:t=min(L,K)-1。根據差分譜的定義知,兩個相鄰的奇異值相差愈大,在差分譜中對應的峰值也愈大,所表現出的特征也愈明顯。選取合適的差分譜譜峰[14]對應的奇異值對信號進行重構可實現信號的降噪處理。由于差分序列的元素數值較大,本文采用變量Sk表示差分譜,其計算公式如下:

(9)

3 基于VMD與SVD降噪的故障特征提取

3.1 VMD分解的優點

EMD采用循環包絡篩分方法處理信號,實現信號從高頻到低頻的自適應劃分,本質上是一組頻率由高到低的帶通濾波器。由于EMD算法本身的缺陷,在劃分信號時會產生模態混疊,強噪聲時更加明顯。此外,模態混疊還會導致包絡解調的頻帶中產生大量無關頻帶,極不利于軸承故障特征的提取。VMD算法利用遞歸迭代計算變分模型最優解來確定每個IMF的頻率中心與帶寬,IMF的頻率中心及帶寬在變分模型的迭代求解中不斷變化,自適應地實現信號頻域剖分與各IMF的分離。每個模態分量的頻帶緊緊圍繞在中心頻率附近,不會出現模態混疊現象,包絡解調時亦沒有無關頻帶,便于故障特征的提取。

3.2 峭度指標對故障信號的篩選

軸承信號經過VMD分解,故障信息被分解到IMF中。若要進行故障特征提取,需選擇合適的指標將包含故障信息的IMF篩選出來。

峭度被用來表示樣本的密度函數圖形頂峰的凸平度,其計算公式如下:

K=E(x-μ)4/σ4

(10)

式中:μ為信號x的均值;σ為信號x的標準差。

當軸承處于正常運行狀態時,振動信號的幅值分布近似于正態分布,其峭度值約等于3。信號中存在較多沖擊成分時,即包含較多故障信息時,信號的峭度會明顯變大[15]。峭度值愈大,信號中沖擊成分所占的比重愈多,故障特征信息也越易提取[16]。

3.3 算法步驟

基于變分模態分解和奇異值分解降噪進行故障特征提取的具體步驟如下:

(1)對原始振動信號進行VMD分解,得到若干IMF。

(2)利用峭度指標對沖擊信號的指示作用,計算不同IMF的峭度值,選擇峭度較大的IMF進行信號重構。

(3)對重構信號的Hankel矩陣進行奇異值分解。

(4)求取奇異值的差分譜,選擇合適的差分譜譜峰,利用其對應的奇異值進行信號重構。

(5)對步驟4得到的信號進行包絡解調,確定故障是否發生,以及故障發生部位。

4 實驗數據分析

4.1 實驗裝置

以美國凱斯西儲大學軸承數據中心的故障數據作為研究對象,對軸承進行故障診斷,并判斷故障類型及其發生的位置。

測試平臺包括驅動電機、測力計、轉矩傳感器和電子控制裝置。電機轉軸由被測軸承支撐,兩端分別為驅動端和風扇端,軸承型號為:SKF6025-2RS。軸承的振動信號由加速度傳感器采集,電機驅動端和風扇端的12點鐘方向(即徑向載荷方向)分別安裝一個加速度傳感器。故障則是通過電火花技術在外圈和內圈上加工凹坑來模擬。

滾動軸承試驗裝置如圖1所示,圖中標注了試驗裝置的各個部件,三個安裝在不同位置的加速度傳感器依次編號為①,②,③,分別表示電機驅動端的傳感器、電機風扇端的傳感器、試驗臺基座上的傳感器。

圖1 試驗測試平臺Fig.1Experiment platform

4.2 實驗故障數據分析

試驗采集的軸承外圈故障振動信號波形與頻譜如圖2(a)和圖2(b)所示,此時轉速傳感器測得的驅動端的轉速為 1 796 r/min,信號采樣頻率為12 000 Hz,采樣點個數為8 000,外圈故障特征頻率為106.5 Hz。

圖2 軸承外圈故障振動信號的波形與頻譜Fig.2 Vibration signal waveform and spectrum of bearing outer ring fault

軸承外圈振動信號VMD分解的各模態頻譜如圖3所示。從圖3可知,原始振動信號經VMD分解得到若干IMF。

圖3 外圈故障信號VMD分解各模態頻譜Fig.3 Vibration signal spectrum of each mode by VMD decomposition of the bearing outer ring fault

不同IMF頻率成分不同,每個IMF緊緊圍繞在某一頻率中心。圖3中各IMF分量的峭度值如表1所示。

表1 各IMF分量的峭度值Tab.1 The kurtosis of each IMF

由表1可知,IMF2和IMF3的峭度值較大,包含的故障信息較多。因而利用IMF2和IMF3得到重構信號u(t),u(t)的波形和頻譜如圖4所示。

圖4 重構信號u(t)的時域波形與頻譜Fig.4 Reconstructed signal u(t) waveform and spectrum

利用奇異譜分解對信號u(t)進行降噪以提高信噪比。對信號u(t)信號的Hankel矩陣進行奇異值分解,奇異值的分布如圖5所示。

圖5 信號u(t)軌道矩陣的奇異值Fig.5 Singular value of u(t) Hankel matrix

計算信號u(t)奇異值的差分譜,結合信號u(t)軌道矩陣的奇異值,利用前14個奇異值對信號進行重構,獲得信號v(t)。信號v(t)的波形與頻譜如圖6(a)、圖6(b)所示。

與圖4(b)中的頻譜相比,可以觀測到經過SVD降噪處理,無關噪聲頻帶被去除,有用信號的頻帶則被最大限度保留。經SVD降噪處理信號和初始信號的信噪比如表2所示。

表2 初始信號與處理信號的信噪比Tab.2 The original and reconstructed signal-to-interference ratio

通過對比可知,信噪比獲得較大提高。

利用Hilbert變換對信號v(t)進行包絡解調,去除高頻振動頻率成分,得到低頻包絡信號k(t)。低頻包絡信號k(t)的頻譜如圖6(c)所示。

圖6 信號v(t)的波形、頻譜、包絡譜Fig.6 Vibration signal waveform, spectrum and envelop spectrum of signal v(t)

轉軸驅動端的轉速為1 796 r/min,轉頻為29.63 Hz。在圖6(c)中,頻率28.5 Hz的位置處存在峰值,十分接近理論轉頻29.63 Hz。軸承處于工作狀態時,滾子和滾道會產生輕微的相對滑動,因而理輪轉頻與實際提取的轉頻存在細微差異。圖6(c)中的105 Hz處可觀測到明顯峰值,對應軸承外圈的故障特征頻率。同時在210 Hz、315 Hz的位置也存在較為明顯的峰值,與故障特征頻率的二倍頻、三倍頻相對應,和軸承外圈的故障特征相吻合。從圖6(c)可知,經過VMD分解和SVD降噪處理,故障特征頻率及其倍頻變得非常突出,軸承的故障特征明顯。

圖7 軸承內圈故障信號的波形、頻譜、包絡譜Fig.7 Vibration signal waveform, spectrum and envelop spectrum of bearing inner ring fault

利用同樣方法處理內圈故障信號w(t),w(t)的波形與頻譜如圖7(a)和圖7(b)所示。此時,驅動端的轉頻為1 794 r/min,理論轉頻為29.9 Hz,經計算理論故障特征頻率為160 Hz。利用上述方法對信號w(t)處理得到的信號包絡譜如圖7(c)所示。

在圖7(c)中,頻率30 Hz、58.5 Hz、88.5 Hz處可觀測到較為明顯的峰值,分別對應轉頻及其二倍頻與三倍頻。此外,頻率160.5 Hz處亦可觀測到明顯波峰,對應著故障特征頻率。此外,在160.5 Hz附近的130.5 Hz和189 Hz處存在較為明顯的峰值,對應著1倍轉頻的調制邊頻帶,與軸承內圈的故障特征相符。

圖8 EMD分解后前3個IMF分量的頻譜Fig.8 Vibration signal spectrum of the first three IMF by EMD decomposition

4.3 與EMD分解的對比分析

采用EMD算法對上述內圈故障信號w(t)進行分析。取前3個對故障特征提取影響較大的本征模態函數進行分析,這 3個IMF的時域波形與頻譜如圖8、圖9所示。

圖9 EMD分解的前3個IMF分量的包絡譜Fig.9 The envelop spectrum of the first three IMF by EMD decomposition

由圖8可知,與VMD分解得到的IMF相比,EMD分解得到的IMF分量會出現頻譜混疊,不同分量之間頻帶相差不明顯,而且會相互影響。

在圖9 EMD分解的IMF包絡譜中,可在故障特征頻率160.5 Hz、轉頻30 Hz及其三倍頻88.5 Hz處觀測到較為明顯的幅值,驗證了本文提出方法在故障特征提取時的準確性。但與VMD分解得到的包絡譜相比,故障特征頻率被大量無關頻帶包圍,不利于故障特征的識別,突出了本文提出方法的優越性。

5 結 論

針對軸承故障信號所在頻帶不易選擇和易受噪聲干擾的問題,提出了基于變分模態分解(VMD)和奇異值分解降噪的軸承故障特征提取方法。該方法通過對振動信號進行VMD分解獲得本征模態函數,利用峭度指標對包含故障信息的本征模態函數進行選擇,有效解決了共振解調提取故障頻率時最優頻帶的選擇問題。經過SVD降噪,去除了背景噪聲的干擾。最后對降噪信號進行包絡解調提取故障特征,取得了良好的效果。

[1] 陳雪峰,李繼猛,程航,等. 風力發電機狀態監測和故障診斷技術的研究與進展[J].機械工程學報,2011,47(9):45-52. CHEN Xuefeng,LI Jimeng,CHENG Hang, et al.Research and application of condition monitoring and fault diagnosis technology in wind turbines [J].Power System Technology,2011,47(9):45-52.

[2] 周智,朱永生,張優云,等. 基于EEMD和共振解調的滾動軸承自適應故障診斷[J]. 振動與沖擊,2013,32(2):76-80. ZHOU Zhi,ZHU Yongsheng,ZHANG Youyun,et al.Adaptive fault diagnosis of rolling bearings based on EEMD and demodulated resonance [J].Journal of Vibration and Shock,2013,32(2):76-80.

[3] 張家凡,易啟偉,李季. 復解析小波變換與振動信號包絡解調分析[J]. 振動與沖擊,2010,29(9):93-96. ZHANG Jiafan,YI Qiwei,LI Ji.Complex analytic wavelet transform and vibration signals envelope-demodulation analysis[J].Journal of Vibration and Shock,2010,29(9):93-96.

[4] 王新晴,馬瑞恒,王耀華,等. 基于一種新的時頻分布的機械故障診斷[J]. 機械工程學報,2003,39(7):150-153. WANG Xinqing,MA Ruihuan,WANG Yaohua,et al.Mechanical fault diagnosis based on a new time-frequency distribution [J].Power System Technology,2003,39(7):150-153.

[5] 馬瑞恒,王新晴. 基于一種新的時頻分布的機械故障診斷[J]. 振動與沖擊,2003,22(3):59-62. MA Ruihuan,WANG Xinqing.Mechanical fault diagnosis based on a new time-frequency distribution [J].Journal of Vibration and Shock,2003,22(3):59-62.

[6] HUANG N E,SHEN Z,LONG S R,et al.The empirical mode decomposition and the Hilbert spectrum for non-liner and non-stationary time series analysis[J]. Proceedings of the Royal Society ,1998,(454):903-993.

[7] 朱文龍,周建中,肖劍,等. 獨立分量分析-經驗模態分解特征提取在水電機組振動信號中的應用[J]. 中國電機工程學報,2013,33(29):95-101. ZHU Wenlong,ZHOU Jianzhong,XIAO Jian, et al.An ICA-EMD feature extraction method and its application to vibration signals of hydroelectric generating units [J].Proceedings of the CSEE,2013,33(29):95-101.

[8] 楊江天,趙明元. 改進雙譜和經驗模態分解在牽引電機軸承故障診斷中的應用[J]. 中國電機工程學報,2012,32(18):116-122. YANG Jiangtian,ZHAO Mingyuan, et al.Fault diagnosis of traction motor bearings using modified bispectrum and empirical mode decomposition [J].Proceedings of the CSEE,2012,32(18):116-122.

[9] 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.

[10] RILLING G,FLANDRIN P.On the influence of sampling on the empirical mode decomposition[C]//Proceedings of IEEE Conference on Acoustics,Speech and Signal Processing. Toulouse,France,2006.

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

[12] 王建國,李健,萬旭東. 基于奇異值分解和局域均值分解的滾動軸承故障特征提取方法[J].機械工程學報,2015,51(3):104-110. WANG Jianguo,LI Jian,WAN Xudong.Fault feature extraction method of rolling bearings based on singular value decomposition and local mean decomposition [J].Power System Technology,2015,51(3):104-110.

[13] 湯寶平,蔣永華,張詳春. 基于形態奇異值分解和經驗模態分解的滾動軸承故障特征提取方法[J]. 機械工程學報,2010,46(5):37-42. TANG Baoping,JIANG Yonghua,ZHANG Xiangchun.Feature extraction method of rolling bearing fault based on singular value decomposition-morphology filter and empirical mode decomposition[J]. Journal of Mechanical Engineering, 2010,46(5):37-42.

[14] 趙學智,葉邦彥,陳統堅. 奇異值差分譜理論及其在車床主軸箱故障診斷中的應用[J]. 機械工程學報,2010,46(1):100-108. ZHAO Xuezhi,YE Bangyan,CHEN Tongjian.Difference spectrum theory of singular value and its application to the fault diagnosis of headstock of lathe[J]. Journal of Mechanical Engineering, 2010,46(1):100-108.

[15] 蘇文勝,王奉濤,張志新,等.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.

[16] 胡愛軍,馬萬里,唐貴基.基于集成經驗模態分解和峭度準則的滾動軸承故障特征提取方[J].中國電機工程學報,2012,32(11):106-111. HU Aijun,MA Wanli,TANG Guiji.An ICA-EMD feature extraction method and its application to vibration signals of hydroelectric generating units [J].Proceedings of the CSEE,2012,32(11):106-111.

Fault feature extraction of bearing faults based on singular value decomposition and variational modal decomposition

ZHAO Hongshan, GUO Shuangwei, GAO Duo

(School of Electrical and Electronic Engineering, North China Electric Power University, Baoding 071003, China)

In order to extract fault features of rolling bearings effectively, a method based on variational mode decomposition and singular value decomposition was proposed. The Intrinsic Mode Function (IMF) was obtained by variational mode decomposition. The IMF containing fault information was selected to reconstruct the signal according to the index of kurtosis. The singular value decomposition was used to reduce noise and increase the ratio of signal-to-noise. Then the fault features were extracted by using envelope spectrum analysis. Compared with common fault features extraction methods, the proposed method can distinguish typical faults, highlight fault features and improve diagnostic effect.

variational mode decomposition; singular value decomposition; rolling bearing; fault features extraction

國家自然科學基金項目(51277074)

2015-08-13 修改稿收到日期:2015-11-18

趙洪山 男,博士,教授,1965年生

郭雙偉 男,碩士,1989年生

TH17

A

10.13465/j.cnki.jvs.2016.22.027

猜你喜歡
特征提取模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
中立型Emden-Fowler微分方程的振動性
一種基于LBP 特征提取和稀疏表示的肝病識別算法
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
基于MED和循環域解調的多故障特征提取
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: www.狠狠| 国产成人综合日韩精品无码首页| 国产精品视屏| 亚洲an第二区国产精品| 青青草原国产| 亚洲人成亚洲精品| 久久国产亚洲欧美日韩精品| 国产在线观看第二页| av在线手机播放| 国产女同自拍视频| 国产精品区网红主播在线观看| 大香伊人久久| 亚洲精品手机在线| 国产乱子伦精品视频| 欧美在线视频a| 亚洲国产成人超福利久久精品| 免费看久久精品99| 在线看片中文字幕| 日韩高清在线观看不卡一区二区 | 青青草原国产av福利网站| 日韩精品毛片| 曰韩免费无码AV一区二区| 免费福利视频网站| 四虎国产精品永久在线网址| 国产浮力第一页永久地址| 日韩黄色精品| 热这里只有精品国产热门精品| 欧美精品在线看| 一级毛片在线免费视频| 欧美啪啪网| 亚洲欧美极品| 欧美成人午夜影院| 久久中文字幕2021精品| 亚洲国产中文精品va在线播放 | 国精品91人妻无码一区二区三区| 日韩欧美国产成人| 日韩午夜伦| 99青青青精品视频在线| 免费一级毛片不卡在线播放 | 久草网视频在线| 秋霞午夜国产精品成人片| 人妻中文久热无码丝袜| 国产日韩欧美成人| 自拍欧美亚洲| 亚洲欧美激情小说另类| 欧美国产菊爆免费观看| 国产午夜在线观看视频| 99热亚洲精品6码| 国产成人免费| 亚洲动漫h| 久久国产黑丝袜视频| 亚洲综合九九| 欧美精品1区2区| 亚洲精品无码不卡在线播放| 国产男人天堂| 国产亚洲精品精品精品| 激情综合网激情综合| AV片亚洲国产男人的天堂| 成年人国产视频| 国产精品毛片一区视频播| 99在线视频免费| 亚洲人网站| 日本国产在线| 精品视频在线观看你懂的一区| 成人蜜桃网| 丰满的少妇人妻无码区| 国产理论最新国产精品视频| 亚洲成人网在线观看| 中文字幕免费播放| 无码人中文字幕| 国产va在线观看| 国产成人精品一区二区| 亚洲一区色| 国产免费久久精品99re丫丫一| 亚洲侵犯无码网址在线观看| 亚洲精品亚洲人成在线| 四虎永久在线| 欧美国产日韩在线观看| 色综合中文| 永久在线播放| 欧美日韩精品在线播放| 丁香亚洲综合五月天婷婷|