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

一種多層經(jīng)驗(yàn)小波變換與多指標(biāo)交叉融合的列車輪對(duì)軸承故障診斷研究

2022-07-22 14:29:52王大鵬鄧飛躍
關(guān)鍵詞:模態(tài)故障信號(hào)

王大鵬,李 忞,鄧飛躍

(1.中車大連機(jī)車車輛有限公司, 遼寧 大連 116000;2.石家莊鐵道大學(xué) 機(jī)械工程學(xué)院,石家莊 050043)

0 引言

鐵路運(yùn)輸是我國(guó)交通運(yùn)輸?shù)囊粋€(gè)重要組成部分,對(duì)保障國(guó)家經(jīng)濟(jì)發(fā)展、人民生活需求具有十分重要的意義。輪對(duì)軸承是列車轉(zhuǎn)向架的核心部件之一,其健康狀態(tài)直接影響車輛的安全運(yùn)營(yíng)。由于惡劣的運(yùn)行工況,列車輪對(duì)軸承極易產(chǎn)生各類故障。相對(duì)于列車車體、控制系統(tǒng)、輪軌等大型部件運(yùn)行安全的高度重視,列車軸承的關(guān)注度較少,目前仍缺乏較為有效的故障診斷方法[1]。為了降低車輛故障性停機(jī)率,避免事故發(fā)生,亟待開展有效的列車輪對(duì)軸承故障診斷研究。

滾動(dòng)軸承故障振動(dòng)信號(hào)具有典型的非平穩(wěn)特征,傳統(tǒng)的時(shí)域、頻域方法無(wú)法同時(shí)兼顧信號(hào)在時(shí)頻域內(nèi)的局部化特征。為此,基于時(shí)頻分解的信號(hào)處理方法取得了很大發(fā)展[2],但目前仍存在一些問(wèn)題有待解決:如短時(shí)傅里葉變換(short fourier transform,STFT)只適合分析緩變信號(hào)[3];Wigner-Ville(WVD)分布易產(chǎn)生交叉項(xiàng)干擾[4];小波變換(wavelet transform,WT)使用固定的基函數(shù),自適應(yīng)性欠缺[5];經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法可以自適應(yīng)分解信號(hào)獲取模態(tài)分量,但缺乏完備的數(shù)學(xué)理論基礎(chǔ),存在模態(tài)混疊、端點(diǎn)效應(yīng)等缺點(diǎn)[6]。

經(jīng)驗(yàn)小波變換(empirical wavelet transform,EWT)是近年來(lái)出現(xiàn)的一種新的時(shí)頻分解方法[7],它基于EMD自適應(yīng)性和小波理論框架,克服了EMD存在的不足,得到了國(guó)內(nèi)外學(xué)者廣泛應(yīng)用,李志農(nóng)等[8]通過(guò)EWT分析轉(zhuǎn)子故障信號(hào),識(shí)別出了不同程度的碰磨故障;Cao等[9]利用EWT分別提取出滾動(dòng)軸承單一故障與復(fù)合故障特征頻率;Kedadouche等[10]則證實(shí)EWT方法要優(yōu)于EMD、EEMD方法,分解模態(tài)分量更為準(zhǔn)確,效率也更高。

雖然EWT方法具有較好性能,但在實(shí)際應(yīng)用過(guò)程中依然存在一些問(wèn)題:① 傳統(tǒng)EWT方法根據(jù)頻譜局部極值點(diǎn)來(lái)劃分頻譜區(qū)間,易受雜頻成分干擾,導(dǎo)致信號(hào)頻譜邊界分割并不準(zhǔn)確;② 劃分信號(hào)頻譜區(qū)間,需要確定故障稀疏度最優(yōu)的頻帶范圍,雖然已有的峭度、平滑因子、峰度系數(shù)等多個(gè)統(tǒng)計(jì)量指標(biāo)可以定量評(píng)價(jià)故障稀疏度的大小,但每個(gè)單一指標(biāo)均存在一定程度的應(yīng)用局限性[11]。

為解決上述問(wèn)題,提出了一種基于多層經(jīng)驗(yàn)小波變換(multi-layer EWT,MLEWT)與多指標(biāo)交叉融合的列車輪對(duì)軸承故障診斷方法。該方法在劃分頻譜邊界過(guò)程中,不再以局部極值點(diǎn)作為頻譜區(qū)間分割的依據(jù),而是通過(guò)設(shè)定頻譜區(qū)間個(gè)數(shù)來(lái)對(duì)信號(hào)頻譜進(jìn)行更為細(xì)致的多層分解,同時(shí)提出了一種多個(gè)統(tǒng)計(jì)量指標(biāo)交叉融合方法來(lái)評(píng)價(jià)、確定信號(hào)中故障稀疏度最優(yōu)的模態(tài)分量信號(hào)。所提方法可以對(duì)信號(hào)進(jìn)行多層次分解,并自適應(yīng)尋找最優(yōu)的模態(tài)分量信號(hào)提取軸承故障特征,通過(guò)仿真信號(hào)與實(shí)驗(yàn)分析,證實(shí)了所提方法的有效性。

1 多層經(jīng)驗(yàn)小波變換(MLEWT)

1.1 傳統(tǒng)經(jīng)驗(yàn)小波變換

給定信號(hào)x(t),其Fourier譜X(f)表示為:

(1)

(2)

(3)

原信號(hào)可重構(gòu)為:

(4)

(5)

有關(guān)EWT方法的詳細(xì)過(guò)程可參考文獻(xiàn)[7],不再贅述。根據(jù)上述過(guò)程,可以發(fā)現(xiàn)EWT方法中,信號(hào)頻譜劃分區(qū)間的個(gè)數(shù)等于分解后模態(tài)分量信號(hào)的個(gè)數(shù),并且頻譜分割后不同區(qū)間范圍對(duì)應(yīng)以不同緊支撐頻率為中心的模態(tài)函數(shù)。因此,如何分割信號(hào)頻譜邊界直接影響EWT分解的結(jié)果。

1.2 MLEWT

受Antoni等[12]所提譜峭度方法的啟示,所提的MLEWT方法不再以信號(hào)頻譜中局部極值點(diǎn)作為頻譜區(qū)間分割的邊界,而是通過(guò)設(shè)定分解層數(shù),來(lái)實(shí)現(xiàn)對(duì)信號(hào)頻譜更為細(xì)致的多層分解。具體分解過(guò)程如下:

1) 設(shè)定劃分信號(hào)頻譜的區(qū)間個(gè)數(shù)為k,則整個(gè)信號(hào)頻譜區(qū)間[0,π]平均分割為k個(gè)連續(xù)區(qū)間(如圖1所示),表達(dá)式如下:

圖1 信號(hào)頻譜劃分

(6)

2) 頻譜區(qū)間劃分后,每個(gè)區(qū)間的頻譜帶寬表示如下:

B=fs/(2k), (k=2,4,…,N)

(7)

式中:fs為所分析信號(hào)的采樣頻率。通過(guò)設(shè)置k值大小,能夠?qū)崿F(xiàn)對(duì)頻譜區(qū)間范圍的調(diào)節(jié),k值越大,信號(hào)頻譜分割后區(qū)間個(gè)數(shù)越多,對(duì)信號(hào)頻帶劃分就越細(xì)致,最終實(shí)現(xiàn)對(duì)信號(hào)頻譜區(qū)間的多層劃分,如圖2所示。

圖2 信號(hào)頻譜多層劃分

MLEWT方法克服了傳統(tǒng)EWT方法根據(jù)頻譜局部極值點(diǎn)劃分頻譜區(qū)間存在的缺陷,不再受頻譜中強(qiáng)背景噪聲所引起的大幅值雜頻成分的干擾,通過(guò)設(shè)置k值大小,可以對(duì)信號(hào)整個(gè)頻譜區(qū)間進(jìn)行更為靈活、細(xì)致的劃分。該方法在繼承傳統(tǒng)EWT方法優(yōu)勢(shì)基礎(chǔ)上,通過(guò)多層分解可以得到多個(gè)模態(tài)分量信號(hào),進(jìn)而獲取包含故障稀疏度特征最為明顯的分量信號(hào)。

2 多指標(biāo)交叉融合

MLEWT方法分解信號(hào)后,將會(huì)得到多個(gè)子分量信號(hào),每個(gè)信號(hào)對(duì)應(yīng)的頻譜區(qū)間范圍各不相同,因此如何確定區(qū)間個(gè)數(shù)k以及如何選取最優(yōu)模態(tài)分量信號(hào),將是本節(jié)所要解決的問(wèn)題。

2.1 不同統(tǒng)計(jì)量評(píng)價(jià)指標(biāo)

當(dāng)滾動(dòng)軸承出現(xiàn)故障,損傷點(diǎn)與其他元件表面接觸時(shí),振動(dòng)信號(hào)會(huì)出現(xiàn)瞬態(tài)沖擊成分,造成信號(hào)某些統(tǒng)計(jì)量指標(biāo)發(fā)生較為明顯的變化。當(dāng)前應(yīng)用較為廣泛的統(tǒng)計(jì)量指標(biāo)有峭度、平滑因子、稀疏值和峰值系數(shù)等。

峭度是描述信號(hào)波形尖峰度的一個(gè)無(wú)量綱參數(shù)[13],數(shù)學(xué)表達(dá)式為:

(8)

式中:u為信號(hào)x(t)的均值;σ為信號(hào)x(t)的指標(biāo)值。峭度指標(biāo)對(duì)沖擊信號(hào)較為敏感,當(dāng)信號(hào)中沖擊成分所占比重較多時(shí),峭度值會(huì)明顯增大,但同時(shí)峭度對(duì)信號(hào)中局部極值點(diǎn)也非常敏感,因此僅采用峭度指標(biāo),難免會(huì)出現(xiàn)判斷失誤。

平滑因子表達(dá)式為:

(9)

式中:GM(·)為信號(hào)x(t)的幾何均值;AM(·)為信號(hào)x(t)的算術(shù)均值。與峭度指標(biāo)相反,信號(hào)中沖擊成分越多時(shí),平滑因子越小,并且平滑因子對(duì)信號(hào)中局部極值不敏感。但是,信號(hào)中噪聲成分對(duì)平滑因子影響較大,導(dǎo)致其可能會(huì)無(wú)法區(qū)分真實(shí)的沖擊分量和噪聲分量[11]。

稀疏值和峰度系數(shù)也是常用于評(píng)價(jià)信號(hào)故障稀疏度的統(tǒng)計(jì)量指標(biāo),它們分別定義如下:

(10)

(11)

式中:RMS(x(t))為根均方函數(shù)。

當(dāng)信號(hào)分量中出現(xiàn)故障沖擊成分時(shí),稀疏值會(huì)明顯增大,但文獻(xiàn)[14]指出采用稀疏值定量描述信號(hào)的沖擊性效果欠佳。峰度系數(shù)與峭度相似,其數(shù)值越大,表明信號(hào)中沖擊分量越大,但它對(duì)沖擊成分的敏感度較低,只有出現(xiàn)較為明顯的沖擊信號(hào)時(shí)才較為有效。雖然上述每個(gè)統(tǒng)計(jì)量指標(biāo)都能在一定程度上評(píng)價(jià)滾動(dòng)軸承信號(hào)中故障稀疏度的大小,但均存在一定的局限性。

2.2 多指標(biāo)交叉融合

多指標(biāo)交叉融合的基本思想是通過(guò)選取多個(gè)表征信號(hào)中故障稀疏度大小的統(tǒng)計(jì)量指標(biāo),將其融合在一起,綜合考慮它們的整體變化情況,克服單一指標(biāo)存在的缺陷,進(jìn)而形成一個(gè)最終指標(biāo),以此確定蘊(yùn)含故障稀疏度最優(yōu)的經(jīng)驗(yàn)?zāi)B(tài)分量。本文中所提出的融合上述4種指標(biāo)的多指標(biāo)交叉融合方法具體步驟如下:

1) 信號(hào)經(jīng)過(guò)MLEWT處理后,針對(duì)每個(gè)固定的k值得到k個(gè)模態(tài)分量信號(hào),分別計(jì)算每個(gè)模態(tài)分量信號(hào)的峭度、平滑因子、稀疏值、峰值系數(shù)指標(biāo)。其中,因?yàn)槠交蜃訑?shù)值大小與信號(hào)故障稀疏程度成反比,所以取其倒數(shù)。將每一類指標(biāo)作為行向量,構(gòu)建指標(biāo)矩陣如下:

(12)

式中:i=k(k=2,3,…,N)。

2) 分別對(duì)矩陣中每一類指標(biāo)進(jìn)行歸一化處理:

(13)

3) 計(jì)算矩陣A中每個(gè)行向量的Shannon熵值:

(14)

Shannon熵可以較好地評(píng)價(jià)某一序列的復(fù)雜程度,其值大小反映了各個(gè)模態(tài)分量在同一統(tǒng)計(jì)量指標(biāo)下的概率分布情況。Shannon熵值越小,各分量信號(hào)所對(duì)應(yīng)的統(tǒng)計(jì)量指標(biāo)數(shù)值彼此越接近,分布越均勻,彼此差別不大;反之,則說(shuō)明彼此差別較大,出現(xiàn)了個(gè)別較為突出的指標(biāo)值,這代表了某個(gè)模態(tài)分量信號(hào)具有較好的故障稀疏度特征。

4) 計(jì)算不同類型指標(biāo)下總熵值和:

(15)

計(jì)算每一類指標(biāo)對(duì)應(yīng)的權(quán)重系數(shù)為:

sj=Ej/H(j=1,2,3,4)

(16)

如果通過(guò)某一類指標(biāo)識(shí)別出表征故障稀疏度大小最優(yōu)的模態(tài)分量信號(hào),則該類指標(biāo)的Shannon熵值會(huì)較大,相應(yīng)的權(quán)重系數(shù)也會(huì)較大;如果在該類指標(biāo)下未找到最優(yōu)分量,則Shannon熵值較小,權(quán)重系數(shù)也相對(duì)較小。雖然某個(gè)統(tǒng)計(jì)量指標(biāo)可能會(huì)因?yàn)檎`判導(dǎo)致數(shù)值較大,但所選的4個(gè)評(píng)價(jià)指標(biāo)彼此補(bǔ)充平衡,只有故障稀疏度最優(yōu)的分量信號(hào)所對(duì)應(yīng)的權(quán)重系數(shù)才能達(dá)到整體數(shù)值最大。

5) 計(jì)算多指標(biāo)交叉融合后各個(gè)模態(tài)分量的最終評(píng)價(jià)指標(biāo)為:

(17)

取maxC(i),(i=1,2,…,k)對(duì)應(yīng)的分量信號(hào)作為信號(hào)MLEWT處理后故障稀疏度大小最優(yōu)的模態(tài)分量。上述建立的多指標(biāo)交叉融合評(píng)價(jià)方法,綜合考慮了4種統(tǒng)計(jì)量指標(biāo)所表征的故障稀疏度大小,避免了單一指標(biāo)可能造成的誤判,有效提高了識(shí)別最優(yōu)模態(tài)分量信號(hào)的可靠性。所提方法的算法流程如圖3所示。

圖3 算法流程框圖

3 仿真分析

滾動(dòng)軸承故障仿真信號(hào)x(t)包含故障脈沖成分、機(jī)械系統(tǒng)其他部件產(chǎn)生的頻率成分和高斯白噪聲組成的背景噪聲成分,數(shù)學(xué)表達(dá)式如下:

x(t)=exp(-1 300t0)(cos(2πf1t))+

0.1sin(2πf2t)+0.15cos(2πf3t)+

xn(t)

(18)

式中:t0=mod(n/fs,fm),n=0,1,2,…,8 191,采樣點(diǎn)長(zhǎng)度為8 192,采樣頻率fs=12 000 Hz,故障特征頻率fm=100 Hz。系統(tǒng)載波頻率f1=3 500 Hz,其他部件產(chǎn)生的振動(dòng)頻率為f2=130 Hz、f3=430 Hz,xn(t)是添加信噪比為-3 dB的高斯白噪聲。

仿真信號(hào)的時(shí)域波形、頻譜及包絡(luò)譜如圖4所示,從時(shí)域波形中很難識(shí)別周期性的故障脈沖成分;從頻譜圖中可以看出低頻周期性頻率成分,但無(wú)法識(shí)別軸承故障特征頻率;包絡(luò)譜中也無(wú)法提取出軸承故障特征頻率。根據(jù)所提方法,對(duì)仿真信號(hào)進(jìn)行MLEWT處理后,對(duì)每個(gè)分量信號(hào)進(jìn)行多指標(biāo)交叉融合計(jì)算,結(jié)果如圖5所示。

圖4 軸承故障仿真信號(hào)

從圖5中可知,利用所提方法對(duì)仿真信號(hào)進(jìn)行了累加的多層分解,通過(guò)計(jì)算得到的多指標(biāo)交叉融合最大值對(duì)應(yīng)的區(qū)間個(gè)數(shù)為k=8,選取第4個(gè)子信號(hào)(紅色虛線框內(nèi))為最優(yōu)的模態(tài)分量信號(hào)。該信號(hào)的時(shí)域波形及包絡(luò)解調(diào)分析結(jié)果如圖6所示。從圖6中可以清晰地提取出故障特征頻率100 Hz及其多倍頻成分,這說(shuō)明通過(guò)所提方法,軸承故障特征信息被準(zhǔn)確提取出來(lái)。

圖6 仿真信號(hào)故障診斷結(jié)果

4 實(shí)驗(yàn)分析

實(shí)驗(yàn)是在鐵道車輛輪對(duì)軸承跑合實(shí)驗(yàn)臺(tái)上完成的,測(cè)試軸承為雙列圓錐滾子軸承,軸承的具體幾何參數(shù)如表1所示。需要說(shuō)明的是,測(cè)試軸承是從某一型號(hào)列車拆卸下來(lái)的,軸承在長(zhǎng)期使用過(guò)程中外圈出現(xiàn)了剝落故障(如圖7所示),此故障為真實(shí)的軸承故障形式,并非人為加工。實(shí)驗(yàn)過(guò)程中,軸承轉(zhuǎn)速為465 r/min,采樣頻率為12 800 Hz。通過(guò)計(jì)算,軸承外圈故障特征頻率為66.75 Hz。

表1 測(cè)試軸承幾何參數(shù)

圖7 測(cè)試軸承故障

輪對(duì)軸承故障振動(dòng)信號(hào)的時(shí)域波形、頻譜及其包絡(luò)譜如圖8所示。時(shí)域波形中沒有發(fā)現(xiàn)明顯的故障沖擊成分,包絡(luò)譜中也未提取出較為明顯的外圈故障特征頻率,因此難以診斷該測(cè)試軸承出現(xiàn)故障。采用本文方法對(duì)該信號(hào)進(jìn)行分析,結(jié)果如圖9所示,從圖9中可知,當(dāng)劃分區(qū)間k=12時(shí),最大多指標(biāo)交叉融合值對(duì)應(yīng)的模態(tài)分量信號(hào)為第11個(gè)(紅色虛線框內(nèi))。提取該信號(hào)進(jìn)行包絡(luò)解調(diào)分析,結(jié)果如圖10所示,包絡(luò)譜中測(cè)試軸承外圈缺陷頻率及其2、3倍頻被清晰地提取了出來(lái),通過(guò)該方法可以準(zhǔn)確診斷出測(cè)試軸承的外圈故障。

圖8 測(cè)試軸承信號(hào)

圖9 內(nèi)圈故障信號(hào)頻譜的分割

圖10 測(cè)試軸承故障診斷結(jié)果

為了進(jìn)一步證明所提方法的有效性,采用快速譜峭度方法和傳統(tǒng)EWT方法分別分析輪對(duì)軸承故障信號(hào)與所提方法。快速譜峭度結(jié)果如圖11所示,最大峭度值對(duì)應(yīng)的中心頻率為fc=5 600 Hz,帶寬為Bw=1 600 Hz。根據(jù)這些參數(shù)構(gòu)建帶通濾波器,濾波后信號(hào)及其包絡(luò)譜如圖12所示。

圖11 測(cè)試軸承快速譜峭度圖

圖12 測(cè)試軸承快速譜峭度診斷結(jié)果

從圖12中可以發(fā)現(xiàn)外圈缺陷頻率,但無(wú)法提取出故障頻率的倍頻成分。傳統(tǒng)EWT方法中,信號(hào)頻譜區(qū)間分割采用locmaxmin方法,分割的個(gè)數(shù)為6,結(jié)果如圖13所示。從圖13中可知,由于頻譜中雜頻的干擾,信號(hào)頻譜區(qū)間的劃分主要集中在低頻段內(nèi),信號(hào)的中高頻段并未進(jìn)行有效的分割。選取其中最優(yōu)模態(tài)分量信號(hào)進(jìn)行包絡(luò)解調(diào)分析,結(jié)果如圖14所示,從圖14中可以提取出軸承外圈故障頻率,沒有明顯的倍頻成分。通過(guò)與圖10結(jié)果進(jìn)行對(duì)比,可以清楚發(fā)現(xiàn)本文中所提方法可以準(zhǔn)確提取出軸承故障特征頻率及多倍頻成分,診斷結(jié)果更為可靠,提取軸承故障特征效果更好。

圖13 測(cè)試軸承EWT結(jié)果

圖14 測(cè)試軸承EWT診斷結(jié)果

5 結(jié)論

1) MLEWT方法不再以傳統(tǒng)EWT方法中根據(jù)頻譜局部極值點(diǎn)來(lái)劃分頻譜區(qū)間,而是通過(guò)設(shè)定區(qū)間個(gè)數(shù)k,對(duì)信號(hào)整個(gè)頻譜區(qū)間進(jìn)行更為靈活的分割。該方法在繼承傳統(tǒng)EWT方法優(yōu)勢(shì)基礎(chǔ)上,通過(guò)多層分解可以得到多個(gè)模態(tài)分量信號(hào)。

2) 多指標(biāo)交叉融合方法綜合考慮峭度、平滑因子、稀疏值和峰值系數(shù)4個(gè)統(tǒng)計(jì)量指標(biāo),通過(guò)交叉融合后的最終指標(biāo)來(lái)評(píng)價(jià)信號(hào)中故障稀疏度的大小,克服了采用單一指標(biāo)引起的誤判,并以此自適應(yīng)搜尋MLEWT后包含故障稀疏特征最明顯的模態(tài)分量信號(hào)。

3) 通過(guò)對(duì)仿真信號(hào)和實(shí)測(cè)列車軸承故障信號(hào)進(jìn)行分析,證實(shí)所提方法能夠有效提取軸承故障特征信號(hào),診斷效果優(yōu)于譜峭度和傳統(tǒng)EWT方法,具有一定的工程應(yīng)用價(jià)值。

猜你喜歡
模態(tài)故障信號(hào)
信號(hào)
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點(diǎn)通
基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
基于LabVIEW的力加載信號(hào)采集與PID控制
國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
故障一點(diǎn)通
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
江淮車故障3例
主站蜘蛛池模板: 亚洲狠狠婷婷综合久久久久| 激情综合激情| 91九色国产在线| 亚洲无码一区在线观看| 久久精品嫩草研究院| 免费国产不卡午夜福在线观看| 美女内射视频WWW网站午夜| 91国内在线观看| 日韩精品中文字幕一区三区| 国产乱视频网站| 亚洲日本中文综合在线| 久久久精品国产SM调教网站| 亚洲人成人伊人成综合网无码| 丝袜高跟美脚国产1区| 黄色网站不卡无码| 999精品色在线观看| 国内精品视频在线| 三区在线视频| 少妇高潮惨叫久久久久久| 亚洲精品午夜天堂网页| 成人在线亚洲| 国产精品性| 日韩人妻精品一区| 男人天堂伊人网| 欧美一级在线看| hezyo加勒比一区二区三区| 一本大道香蕉高清久久| 国产一级无码不卡视频| 成人夜夜嗨| 国产成人精品综合| 99re免费视频| 久久毛片网| 国产亚洲美日韩AV中文字幕无码成人| 精品五夜婷香蕉国产线看观看| 9丨情侣偷在线精品国产| 日韩欧美中文| 综合五月天网| 91精品久久久久久无码人妻| 精品欧美日韩国产日漫一区不卡| 欧美日本在线播放| 中文字幕永久在线看| 国内精品手机在线观看视频| 99精品在线视频观看| 色综合成人| 幺女国产一级毛片| 成年av福利永久免费观看| 日韩av无码精品专区| 国产女人在线观看| 欧美综合激情| 亚洲a级在线观看| www.亚洲色图.com| 久久国产精品影院| 欧美成人午夜视频| 精品国产免费观看| 小蝌蚪亚洲精品国产| 香蕉久人久人青草青草| 激情成人综合网| 天天婬欲婬香婬色婬视频播放| 久久精品国产精品青草app| 无码中文字幕乱码免费2| 欧美日韩专区| 久久精品国产在热久久2019| 中文一级毛片| 在线高清亚洲精品二区| 国产手机在线观看| 国产精品免费露脸视频| 国产精品欧美亚洲韩国日本不卡| 97人妻精品专区久久久久| 最新痴汉在线无码AV| 亚洲精品777| 国产在线98福利播放视频免费 | 欧美亚洲一区二区三区导航| 国产美女精品一区二区| 伊人久久久久久久| 国产一区二区网站| 91在线高清视频| 黑人巨大精品欧美一区二区区| 国产成人喷潮在线观看| 欧美黄网在线| 无码精品福利一区二区三区| 毛片免费在线视频| 亚洲成人高清在线观看|