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

自適應無參經驗小波變換及其在轉子故障診斷中的應用

2016-09-08 06:57:45鄭近德潘海洋潘紫微羅潔思
中國機械工程 2016年16期
關鍵詞:故障診斷故障信號

鄭近德 潘海洋 潘紫微 羅潔思

1.安徽工業大學,馬鞍山,243032  2.廈門理工學院,廈門,361024

?

自適應無參經驗小波變換及其在轉子故障診斷中的應用

鄭近德1潘海洋1潘紫微1羅潔思2

1.安徽工業大學,馬鞍山,2430322.廈門理工學院,廈門,361024

為了實現經驗小波變換中Fourier譜的自適應分割,提出了自適應無參經驗小波變換(APEWT)方法。同時,為了克服希爾伯特變換解調的不足,更精確地估計信號的時頻分布,提出了改進歸一化希爾伯特變換(INHT)。通過分析仿真信號將APEWT和INHT方法與經驗模態分解(EMD)、總體平均經驗模態分解(EEMD)和局部特征尺度分解等方法進行對比,結果表明了APEWT和INHT方法的優越性。最后,將基于APEWT和INHT的時頻分析方法應用于轉子局部碰磨故障診斷,試驗數據分析結果表明,所提出的方法不僅能夠有效地診斷轉子局部碰磨故障,而且診斷效果優于EMD和EEMD方法。

轉子故障;碰磨;經驗模態分解;希爾伯特變換

0 引言

時頻分析由于能夠同時提供信號的時域和頻域局部信息而在機械故障診斷領域得到了廣泛的應用[1]。機械故障診斷中常用的時頻分析方法主要有Wigner-Ville分布、小波變換、經驗模態分解(empiricalmodedecomposition,EMD)[1-3],以及在EMD基礎上發展的局部均值分解(localmeandecomposition,LMD)和局部特征尺度分解(localcharacteristic-scaledecomposition,LCD)等[4-7]。來五星等[8]將Wigner-Ville時頻分布應用于齒輪故障診斷;Yan等[9]對小波分析在旋轉機械故障診斷中的應用進行了總結,并對小波變換在故障診斷的應用前景進行了預測;于德介等[10]將EMD引入到機械故障診斷領域,并對其存在的問題進行了改進和完善;雷亞國[11]研究了基于總體平均經驗模態分解(ensembleempiricalmodedecomposition,EEMD)改進的希爾伯特-黃變換(Hilbert-Huangtransform,HHT)在機械故障診斷中的應用。然而,這些時頻分析方法都有各自的固有缺陷[12]。

最近,在小波變換的基礎上,Gilles等[13-14]提出了一種新的非平穩信號分析方法——經驗小波變換(empiricalwavelettransform,EWT)。EWT通過對Fourier譜進行分割,在每個分割的區間建立小波正交基,從而能夠將一個多分量信號分解為若干個具有緊支集頻譜的調幅調頻信號之和。但是,EWT方法的一個關鍵問題是Fourier譜的分割。為了實現Fourier譜的自適應分割,筆者提出了一種新的非平穩信號分解方法——自適應無參經驗小波變換(adaptiveparameterlessEWT,APEWT)方法。

采用APEWT對多分量信號進行分解,得到若干個單分量信號,再對每個單分量信號進行解調即可得到原始信號的完整時頻分布。為了更精確地估計信號的時頻分布,本文提出一種改進的歸一化希爾伯特變換(improvednormalizedHilberttransform,INHT)方法。INHT方法避開了希爾伯特變換(Hilberttransform,HT)方法構造解析信號,瞬時特征的估計精度更高。

1 自適應無參經驗小波變換

本文提出的APEWT方法主要包含兩部分:Fourier譜的自適應分割和濾波器組的建立。APEWT的關鍵是對原始信號的Fourier譜進行自適應分割。文獻[15]提出將Fourier譜轉化成尺度空間表示(函數的尺度空間表示參見文獻[16]),將Fourier譜的自適應分割問題轉化為如何在尺度空間表示中找到“有意義”模態,并將此問題轉化為一個二類聚類問題。

(1)

(2)

為方便,假設τn與ωn成正比,即τn=γωn,0<γ<1。β(x)是一任意區間[0,1]連續的任意函數,依據Meyer小波,選擇β(x)=x4(35-84x+70x2-20x3)。

然后,采用類似經典的小波變換方法,通過對信號和經驗小波作內積運算的方式得到經驗小波變換系數,通過對信號和尺度函數作內積運算得到逼近系數,進而得到原始信號的各階模態。

本質上,對任意一個復雜的多分量信號,APEWT建立一組自適應的帶通濾波器,通過對信號進行帶通濾波的方式將其分解為若干個具有緊支集頻譜的單分量信號之和,再對得到的單分量信號進行解調,即可得到原始信號完整的時頻分布。

2 改進的歸一化希爾伯特變換

HT是一種常用且有效的信號解調方法。但由于Bedrosian定理的限制,HT會出現無法解釋的負頻率。為了克服此問題,文獻[18]提出了歸一化HT(normalized Hilbert transform,NHT),NHT首先對EMD分解得到的內稟模態函數(intrinsic mode function,IMF)進行歸一化,再對歸一化后的IMF作HT。NHT的對象是歸一化IMF,不再受Bedrosian定理的限制,但由于受Nuttall定理的局限,跟真實頻率之間仍有一定的誤差,而且HT和NHT都有嚴重的端點效應。

為了克服HT和NHT方法的缺陷,本文提出了改進的NHT方法——INHT方法。NHT和INHT都是基于經驗調幅調頻分解(empirical AM-FM demodulation,EAD)[12,18]方法而提出的,EAD方法能夠將一個單分量信號(如IMF)分解為調幅部分和調頻部分的乘積,調幅部分即為單分量信號的瞬時幅值。

INHT的具體步驟如下:

(1)對于單分量信號c(t),首先可采用EAD將其寫作調幅部分和調頻部分的乘積,即:c(t)=a(t)c1(t),其中,a(t)是c(t)的瞬時幅值,c1(t)是歸一化信號,幅值為1。

3 仿真信號分析

令仿真信號

s(t)=s1(t)+s2(t)+s3(t)+s4(t)

(3)

t∈[0,1]

其中,s1(t)是頻率為200 Hz的載波被頻率為指數函數調制構成的調幅調頻信號,s2(t)=(1+0.3cos(2π(5t))sin(2π(80t)+2π(5t2)),s3(t)=2(t2+1)cos(2π(40t)),s4(t)=2t2。各個分量的時域波形如圖1中虛線所示。

圖1 仿真信號s(t)的各個成分及APEWT分解結果

采用APEWT方法對s(t)進行分解,結果如圖1中實線所示,其中Fourier譜的分割中采用基于K均值聚類的尺度空間法。通過與真實分量(虛線)對比可以發現,APEWT方法得到的四個分量C1、C2、C3、C4分別對應真實分量s1、s2、s3、s4,對應的實線和虛線幾乎重合,它們對應的相關系數分別為0.9929、0.9935、0.9987和0.9973,因此,APEWT方法得到的分量與真實分量的相關性和吻合度很高,分解效果比較理想。

圖2 仿真信號s(t)的EMD分解結果

為了對比,再分別采用EMD、EEMD[17]和最近提出的LCD[6,12]對s(t)進行分解。其中,EMD方法的分解結果如圖2所示。為節約篇幅,EEMD和LCD方法的分解未給出,下文將給出它們的時頻譜。由圖2可以看出,EMD的分解結果發生了嚴重的模態混疊,而且C4、C5、C6為擬合產生的虛假分量。對應的前三個分量和剩余趨勢項與對應真實分量s1、s2、s3、s4的相關系數分別為0.9275、0.5408、0.8358和0.9990,分解結果不如APEWT的結果理想。EEMD由于添加的白噪聲并不能完全抵消會產生高頻的噪聲分量,而且仍出現了一個模態被分解為相鄰的幾個IMF分量的情況;LCD也產生了較多的虛假分量,而且分解的分量與真實分量的吻合度較低。因此,APEWT方法得到的分量與真實分量的吻合度和相似性最高,分解效果最為理想。

采用APEWT方法對原始信號進行分解,得到若干個瞬時頻率具有物理意義的單分量信號,再對它們進行解調即可得到原始信號的時頻分布。

分別采用HT、NHT、INHT1和INHT2四種方法估計單分量信號s2(t)和s3(t)的瞬時頻率,結果如圖3所示,其絕對誤差分別如圖4所示。由圖3可知,四種方法估計的s2(t)的瞬時頻率都比較精確,差別很小。但由圖4的絕對誤差中可以看出,HT和NHT方法估計結果有明顯的端點效應;而INHT2估計的結果與真實值非常接近,絕對誤差水平非常小,結果比較理想。對于s3(t),分別采用四種方法估計的瞬時頻率中,HT端點效應嚴重,且由兩端向中間傳播,導致估計值波動較大;而NHT方法的端點效應較小,但仍存在;INHT1克服了二者的端點效應,估計結果比較精確,但估計值比真實值偏小;INHT2方法不僅克服了HT和NHT方法的端點效應,而且絕對誤差也非常小,和真實值非常吻合。綜上分析,本文提出的估計瞬時頻率的INHT方法克服了HT方法的端點效應,估計值更加精確,尤其是INHT2方法,估計效果要優于其他方法。最后,為了對比,分別采用EMD-HT,EEMD-NHT,LCD-INHT1和APEWT-INHT2四種方法估計信號s(t)的時頻譜,結果分別如圖5a~圖5d所示。由圖5可以看出,圖5a中EMD方法分解得到的偽分量較多,HT方法有明顯的端點效應,而且出現了負頻率;圖5b中EEMD方法得到的分量較多,仍存在模態混疊現象;LCD-INHT1方法中,LCD分解得到的分量也較多;而本文提出的APEWT-INHT2方法得到的時頻譜最接近原信號的理想時頻譜。因此,上述分析表明了基于APEWT和INHT2的時頻分析方法的優越性。

(a) 不同方法估計的s2(t)的瞬時頻率

(b)不同方法估計的s3(t)的瞬時頻率圖3 四種不同方法估計的s2(t)和s3(t)的瞬時頻率

(a) 不同方法估計的s2(t)的瞬時頻率的絕對誤差

(b)不同方法估計的s3(t)的瞬時頻率的絕對誤差圖4 四種不同方法估計的s2(t)和s3(t)的瞬時頻率絕對誤差

(a)EMD-HT時頻譜  (b)EEMD-NHT時頻譜

(c)LCD-INHT 1時頻譜 (d)APEWT-INHT 2時頻譜圖5 四種不同方法得到的信號s(t)的時頻譜

4 實測信號分析

局部碰磨是轉子系統常見的非線性振動故障,當轉子發生局部碰磨時,轉子在旋轉的過程中動靜件會周期性地發生摩擦,其碰磨信號往往表現為調幅特征,由于具有調幅特征的碰磨信號非常微弱,因此,如何從信號中提取包含故障特征信息的調幅信號是轉子局部碰磨故障診斷的關鍵[19]。

為了驗證APEWT和INHT方法的有效性,將其應用于轉子動靜碰磨故障實驗數據的分析。采用Z-T3型轉子振動模擬實驗臺進行轉子系統故障實驗,實驗裝置參見文獻[19]。在采樣頻率為2048Hz、轉頻為50Hz條件下,采集到存在局部碰磨故障的轉子徑向位移信號,其時域波形及頻譜如圖6所示,由時域波形和頻譜僅能看出工頻,而與故障有關的碰磨特征被背景信號和噪聲淹沒。

(a)轉子位移信號

(b)轉子位移信號的頻譜圖6 存在局部碰磨故障的轉子位移信號與頻譜

采用APEWT對存在故障的轉子徑向位移信號進行分解,結果如圖7所示。第一個分量的包絡譜和第二、三、四個分量的頻譜分別如圖8所示。綜合兩圖可以看出,第一個分量出現了高頻成分的頻率調制現象,具有明顯的調幅特征,由其包絡譜可以看出,調幅的頻率為50Hz,剛好等于工頻,這是由于轉子每轉動一周,動靜件就碰磨一次造成的。第二、三、四個分量分別為工頻的三倍頻、二倍頻和一倍頻,而且,高頻諧波分量(三倍頻、二倍頻)的幅值較大,說明碰磨的程度較重[20]。因此,通過對轉子局部碰磨信號進行APEWT分解,將碰磨信號、背景信號和噪聲分離,可將高頻碰磨信號從強大的背景信號中提取出來。上述分析表明,APEWT方法能夠有效地提取轉子碰磨故障特征信息。

圖7 存在局部碰磨故障轉子位移信號的APEWT分解結果

(a)第一個分量的包絡譜

(b)第二個分量的頻譜

(c)第三個分量的頻譜

(d)第四個分量的頻譜圖8 APEWT分解第一個分量的包絡譜和其他分量的頻譜

為了對比,采用EMD方法對上述存在碰磨故障的轉子徑向位移信號進行分解,結果如圖9所示。由圖9可以看出,EMD第一個分量雖然具有調幅特征,但是局部發生了模態混疊,第一個分量的調制頻率不僅有轉頻及其二倍頻,而且還包含有低頻調制,對故障診斷造成干擾,而且EMD還無法分解出二倍轉頻和三倍轉頻信息。再采用EEMD方法對上述轉子徑向位移信號進行分解,前四個IMF分量和其對應的殘余分量結果如圖10所示,其中添加噪聲的標準差大小為0.15,總體平均次數為100。由圖10可以看出,雖然EEMD能夠得到具有調幅特征的碰磨信號,但無法分解出二倍轉頻和三倍轉頻,第二個分量仍是一個多分量信號,物理意義不明確。因此,APEWT的診斷效果要優于EMD和EEMD方法。

圖9 存在局部碰磨故障轉子位移信號的EMD分解結果

圖10 存在局部碰磨故障轉子位移信號的EEMD分解結果

圖11為正常轉子徑向位移信號(圖中第一條曲線)及其APEWT分解,結果(圖中第2~6條曲線)所示。第一個分量的包絡譜和頻譜如圖12所示。綜合兩圖可以看出,正常轉子位移信號的第一個分量的調幅特征不明顯,通過其頻譜發現,第一個分量為工頻的四倍頻,而第二、三、四個分量分別為工頻的三倍頻、二倍頻和一倍頻,沒有明顯碰磨故障特征信息,與存在碰磨故障的信號分析結果區別明顯。

圖11 正常轉子位移信號及其APEWT分解結果

(a) 第一個分量的包絡譜

(b)第一個分量的頻譜圖12 正常轉子信號APEWT分解的第一個分量的包絡譜和頻譜

5 結論

(1)提出了一種自適應無參經驗小波變換(APEWT)方法,通過仿真試驗信號將其與EMD、LCD和EEMD等方法進行了對比,結果表明,APEWT方法得到的分量更精確。

(2)針對希爾伯特變換等現有解調方法的不足,提出一種新的瞬時頻率估計方法——改進的歸一化希爾伯特變換。將其與希爾伯特變換和歸一化希爾伯特變換進行了對比,結果表明了本文提出的方法端點效應更小,精確性更高。

(3)將提出的基于APEWT和INHT的時頻分析方法應用于正常轉子和轉子局部碰磨故障的信號分析,通過與EMD和EEMD等方法進行對比,結果表明APEWT方法不僅能有效地診斷轉子局部碰磨故障,而且診斷效果優于EMD和EEMD等方法。

APEWT方法也有不足之處,如分解結果依賴于頻譜的劃分,也存在輕微的端點效應,筆者正針對這些問題展開進一步研究。

[1]何正嘉, 訾艷陽, 孟慶豐, 等. 機械設備非平穩信號的故障診斷原理及應用[M]. 北京: 高等教育出版社, 2001.

[2]吳強, 孔凡讓, 何清波,等. 基于小波變換和ICA的滾動軸承早期故障診斷[J]. 中國機械工程, 2012, 23(7): 835-840.

WuQiang,KongFanrang,HeQingbo,etal.EarlyFaultDiagnosisofRollingElementBearingBasedonWaveletTransformandIndependentComponentAnalysis[J].ChinaMechanicalEngineering, 2012, 23(7): 835-840.

[3]鄭近德, 程軍圣, 曾鳴, 等. 基于偽極值點假設的經驗模態分解及其在轉子故障診斷中的應用[J]. 中國機械工程, 2014,25(18): 2467-2472.

ZhengJinde,ChengJunsheng,ZengMing,etal.Pseudo-extrema-basedEMDandItsApplicationtoRotorFaultDiagnosis[J].ChinaMechanicalEngineering, 2014,25(18): 2467-2472.

[4]SmithJS.TheLocalMeanDecompositionandItsApplicationtoEEGPerceptionData[J].JournaloftheRoyalSocietyInterface, 2005, 2(5):443-454.

[5]張亢, 程軍圣, 楊宇. 基于局部均值分解的階次跟蹤分析及其在齒輪故障診斷中的應用[J]. 中國機械工程, 2011, 22(14):1732-1736.ZhangKang,ChengJunsheng,YangYu.OrderTrackingAnalysisbasedonLocalMeanDecompositionandItsApplicationstoGearFaultDiagnosis[J].ChinaMechanicalEngineering, 2011, 22(14):1732-1736.

[6]程軍圣, 張亢, 楊宇. 局部均值分解方法及其在滾動軸承故障診斷中的應用[J]. 中國機械工程, 2009,20(22):2711-2717.

ChengJunsheng,ZhangKang,YangYu.LocalMeanDecompositionandItsApplicationstoRollingBearingFaultDiagnosis[J].ChinaMechanicalEngineering, 2009,20(22):2711-2717.

[7]鄭近德, 程軍圣, 楊宇. 部分集成局部特征尺度分解:一種新的基于噪聲輔助數據分析方法[J]. 電子學報, 2013, 41(5):1030-1035.

ZhengJinde,ChengJunsheng,YangYu.PartlyEnsembleLocalCharacteristic-ScaleDecomposition:aNewNoiseAssistedDataAnalysisMethod[J].ActaElectronicaSinica. 2013, 41(5):1030-1035.

[8]來五星, 軒建平, 史鐵林, 等.Wigner-Ville時頻分布研究及其在齒輪故障診斷中的應用[J]. 振動工程學報,2003, 16(2):147-250.

LaiWuxing,XuanJianping,ShiTielin,etal.ResearchofWigner-VilleTimeFrequencyandApplicationinDetectingGearPinionFault[J].JournalofVibrationEngineering, 2003,16(2):147-250.

[9]YanRQ,GaoRX,ChenX.WaveletsforFaultDiagnosisofRotaryMachines:aReviewwithApplications[J].SignalProcessing, 2014, 96:1-15.

[10]于德介, 程軍圣, 楊宇. 機械故障診斷的Hilbert-Huang變換方法[M]. 北京:科學出版社, 2006.

[11]雷亞國. 基于改進Hilbert-Huang變換的機械故障診斷[J]. 機械工程學報, 2011, 47(5):71-77.

LeiYaguo.MachineryFaultDiagnosisBasedonImprovedHilbert-HuangTransform[J].JournalofMechanicalEngineering, 2011, 47(5):71-77.

[12]程軍圣, 鄭近德, 楊宇. 基于局部特征尺度分解的經驗包絡解調方法及其在機械故障診斷中的應用[J]. 機械工程學報, 2012, 48(19):87-94.

ChengJunsheng,ZhengJinde,YangYu.EmpiricalEnvelopeDemodulationApproachBasedonLocalCharacteristic-scaleDecompositionandItsApplicationstoMechanicalFaultDiagnosis[J].JournalofMechanicalEngineering, 2012, 48(19):87-94.

[13]GillesJ.EmpiricalWaveletTransform.IEEETransactionsonSignalProcessing[J]. 2013, 61(16):3999-4010.

[14]GillesJ,TranG,OsherS. 2DEmpiricalTransforms:Wavelets,RidgeletsandCurveletsRevisited[J].SIAMJournalonImagingSciences, 2014, 7(1):157-186.

[15]GillesJ,HEALK.AParameterlessScale-spaceApproachtoFindMeaningfulModesinHistograms-applicationtoImageandSpectrumSegmentation[J].InternationalJournalofWaveletsMultiresolutionandInformationProcessing, 2014, 12(6):1450044(1-17).

[16]LindebergT.Scale-spaceTheoryinComputerVision[M].Berlin:Springer, 1994.

[17]WuZH,HuangNE.EnsembleEmpiricalModeDecomposition:aNoise-assistedDataAnalysisMethod[J].AdvancesinAdaptiveDataAnalysis, 2011, 1(1):1-41.

[18]HuangNE,WuZ,LongSR,etal.OnInstantaneousFrequency[J].AdvancesinAdaptiveDataAnalysis, 2009, 1(2):177-229.

[19]程軍圣. 基于Hilbert-Huang變換的旋轉機械故障診斷方法研究[D]. 長沙:湖南大學, 2005.

[20]胡愛軍, 安連鎖, 唐貴基. 轉子碰摩故障振動時頻特征的實驗研究[J]. 動力工程學報, 2007, 27(4):482-486.

HuAijun,AnLiansuo,TangGuiji.ExperimentalStudyonTime-frequencyCharacteristicsofRubbingInitiatedVibrationofRotor[J].JournalofPowerEngineering, 2007, 27(4):482-486.

(編輯盧湘帆)

AdaptiveParameterlessEmpiricalWaveletTransform(EWT)andItsApplicationstoFaultDiagnosisofRotorSystem

ZhengJinde1PanHaiyang1PanZiwei1LuoJiesi2

1.AnhuiUniversityofTechnology,Maanshan,Anhui,243032 2.XiamenUniversityofTechnology,Xiamen,Fujian,361024

TofulfillanadaptiveseparationofFourierspectruminEWT,anadaptiveparameterlessEWT(APEWT)methodwasproposedherein.ToovercometheshortcomingsofHilberttransformandestimatemoreaccuratetime-frequencydistributionofagivensignal,animprovednormalizedHilberttransform(INHT)wasputforward.TheproposedAPEWTandINHTwerecomparedwithempiricalmodedecomposition(EMD),ensembleEMD(EEMD)andlocalcharacteristic-scaledecomposition(LCD)methodsandtheanalysisresultsdemonstratetheeffectivenessoftheproposedmethod.Finally,APEWTandINHTbasedtime-frequencyanalysismethodwereappliedtolocalrubbingfaultdiagnosisofarotorsystem,andtheanalysisresultsofexperimentaldataindicatethattheproposedmethodmayfulfillrotorrubbingfaultdiagnosiseffectivelyandhavebettereffectivenessthanthatofEMDandEEMDmethods.

rotorfault;rubbing;empiricalmodedecomposition;Hilberttransform

2016-01-08

國家自然科學基金資助項目(51505002);安徽省高校自然科學研究重點項目(KJ2015A080);福建省自然科學基金資助項目(2014J05065)

TH165.3; TN911.7

10.3969/j.issn.1004-132X.2016.16.016

鄭近德,男,1986年生。安徽工業大學機械工程學院講師、博士。主要研究方向為非線性動力學、動態信號處理和機械故障診斷等。發表論文近40篇。潘海洋,男,1989年生。安徽工業大學機械工程學院助教、碩士。潘紫微,男,1956年生。安徽工業大學機械工程學院教授。羅潔思,女,1985年生。廈門理工學院機械與汽車工程學院副教授。

猜你喜歡
故障診斷故障信號
信號
鴨綠江(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在线| 国产呦视频免费视频在线观看| 91精品国产一区| 综合色区亚洲熟妇在线| 伊人久久久久久久| 午夜在线不卡| 国产va视频| 人与鲁专区| 国产一级α片| 国产一区二区三区夜色| 老司机aⅴ在线精品导航| 天堂成人在线| 爆乳熟妇一区二区三区| 国产精品任我爽爆在线播放6080| 免费在线a视频| 国产在线精品美女观看| 中文无码伦av中文字幕| 国产精品主播| 伊人久久婷婷| 少妇高潮惨叫久久久久久| 激情网址在线观看| 国产精品视频导航| 免费看黄片一区二区三区| 国产无吗一区二区三区在线欢| 成人福利在线免费观看| 成人精品午夜福利在线播放| 天天摸天天操免费播放小视频| 午夜激情婷婷| 夜色爽爽影院18禁妓女影院| 亚洲最新在线| 久久综合色视频| 看你懂的巨臀中文字幕一区二区| 欧美19综合中文字幕| 国产亚洲精久久久久久久91| 国产精品自在在线午夜| www亚洲天堂| 青青网在线国产| 91区国产福利在线观看午夜| 国产亚洲视频在线观看| 天天摸夜夜操| 国产一二视频| 凹凸国产熟女精品视频| 777国产精品永久免费观看| 国产日韩欧美成人| 无码免费试看| 精品1区2区3区| аⅴ资源中文在线天堂| 日韩区欧美国产区在线观看| 97综合久久| 国产精品欧美在线观看| 欧美在线一级片| 久久久久久国产精品mv| 国产丝袜无码一区二区视频| 国产精品熟女亚洲AV麻豆| 在线观看亚洲天堂| 992Tv视频国产精品| 日本亚洲国产一区二区三区| 日本国产精品一区久久久| 亚洲日韩AV无码一区二区三区人| 亚洲国产日韩在线观看| 国产精品久久自在自线观看| 毛片手机在线看| 国产黄色片在线看| 2048国产精品原创综合在线| 日韩免费毛片| 亚洲成人在线网| 国产一级裸网站| 成人在线第一页| 国产精品久久久久婷婷五月| 欧美翘臀一区二区三区| 欧美午夜理伦三级在线观看| 一本色道久久88| 成人福利在线视频| 操国产美女| 无码日韩精品91超碰| 国产人成网线在线播放va|