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

基于MF-SVD的滾動(dòng)軸承振動(dòng)信號(hào)故障特征提取方法研究

2017-07-05 10:37:10薩其日拉
石油化工自動(dòng)化 2017年3期
關(guān)鍵詞:特征提取振動(dòng)故障

薩其日拉

(東北石油大學(xué) 石油工程學(xué)院,黑龍江 大慶 163318)

基于MF-SVD的滾動(dòng)軸承振動(dòng)信號(hào)故障特征提取方法研究

薩其日拉

(東北石油大學(xué) 石油工程學(xué)院,黑龍江 大慶 163318)

極值域均值模式分解(EMMD)在抑制端點(diǎn)效應(yīng)、算法精度、計(jì)算耗時(shí)等方面均比經(jīng)驗(yàn)?zāi)J椒纸?EMD)算法和自適應(yīng)時(shí)變?yōu)V波分解(ATVFD)有著明顯的優(yōu)勢(shì),能夠有效地對(duì)旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)進(jìn)行故障特征提取,由于現(xiàn)場(chǎng)信號(hào)通常摻雜大量噪聲,嚴(yán)重影響了EMMD的分解精度。針對(duì)該問(wèn)題,提出了基于形態(tài)濾波-奇異值(MF-SVD)的去噪方法,并將其與EMMD相結(jié)合,建立了一種新的故障特征提取方法。實(shí)驗(yàn)結(jié)果表明: 該方法能夠有效、準(zhǔn)確地提取旋轉(zhuǎn)機(jī)械滾動(dòng)軸承內(nèi)圈損傷的故障特征,有著良好運(yùn)算速度和精確度。

極值域均值模式分解 形態(tài)濾波 奇異值分解 特征提取

1998年由N.E. Huang等人[1]提出一種新的時(shí)頻分析法HHT(hilbert-huang transform),能夠有效地對(duì)非線(xiàn)性、非平穩(wěn)信號(hào)進(jìn)行分析,同時(shí)具有良好的自適應(yīng)性,在機(jī)械故障診斷領(lǐng)域中HHT法已得到廣泛的應(yīng)用,并取得了良好的效果。但是N.E. Huang也指出該方法尚有不足之處,如三次樣條插值的“過(guò)包絡(luò)”和“欠包絡(luò)”、端點(diǎn)效應(yīng)以及模態(tài)混疊問(wèn)題[2-3]。近年來(lái),經(jīng)過(guò)專(zhuān)家學(xué)者的不斷研究,出現(xiàn)很多經(jīng)驗(yàn)?zāi)J椒纸?EMD)的改進(jìn)方法,其中較為先進(jìn)的為自適應(yīng)時(shí)變?yōu)V波分解法(ATVFD)和極值域均值模式分解法(EMMD)。經(jīng)研究表明,EMMD法與ATVFD法和EMD法相比,具有有效抑制端點(diǎn)效應(yīng)、高分解精度和更快的運(yùn)算速度的特點(diǎn)。

但是,現(xiàn)場(chǎng)采集的振動(dòng)信號(hào)通常摻雜大量噪聲,大部分為隨機(jī)噪聲和局部沖擊干擾,這些噪聲同樣進(jìn)行EMMD分解計(jì)算。因此,在對(duì)信號(hào)分析之前有必要對(duì)其進(jìn)行預(yù)處理,去除各種噪聲干擾影響,提高EMMD法的分解質(zhì)量和信號(hào)分析的可靠性和精度。

針對(duì)非線(xiàn)性、非平穩(wěn)的機(jī)械振動(dòng)信號(hào)有很多去噪方法,如沈路[4]等利用形態(tài)濾波對(duì)齒輪振動(dòng)信號(hào)濾波消噪提取了故障特征,脈沖噪聲干擾得到了有效抑制,隨機(jī)噪聲也得到了削弱,但是隨機(jī)噪聲還是比較明顯,信號(hào)依然受到噪聲干擾。曾作欽[5]等利用奇異值分解對(duì)信號(hào)進(jìn)行預(yù)處理后再進(jìn)行EMD法分解,取得比直接EMD法分解更好的效果,但是該方法對(duì)窄帶脈沖的去噪效果不甚理想。

因此,文中提出了一種基于形態(tài)濾波-奇異值(MF-SVD)和EMMD的振動(dòng)信號(hào)故障特征提取方法。該方法先將實(shí)測(cè)現(xiàn)場(chǎng)信號(hào)進(jìn)行基于MF-SVD法去噪,濾除隨機(jī)噪聲信號(hào)和局部強(qiáng)干擾噪聲,再進(jìn)行EMMD分解,得到固有模態(tài)函數(shù)(IMF)分量,從而提取故障特征。實(shí)驗(yàn)結(jié)果表明: 該方法能夠有效、準(zhǔn)確地提取滾動(dòng)軸承內(nèi)圈損傷的故障特征,還可抑制端點(diǎn)效應(yīng),有著良好運(yùn)算速度和精確度。

1 EMMD算法原理及實(shí)測(cè)信號(hào)中噪聲影響

1.1 EMMD算法原理

EMMD法是基于EMD法和ATVFD法的一種改進(jìn)算法,其計(jì)算局部均值既不同于EMD法的基于局部極值包絡(luò),也不同于ATVFD法的基于定積分中值定理原理的均值曲線(xiàn)[6]。

對(duì)于原信號(hào)數(shù)據(jù)求出所有N個(gè)局部極值點(diǎn),不需區(qū)分局部極大值和局部極小值,組成{e(ti)},其中ti(i=1, 2, …,N)為第i個(gè)局部極值點(diǎn)的時(shí)間位置。于是信號(hào)行x(t)的數(shù)據(jù)序列值,對(duì)ti與ti+1時(shí)刻極值點(diǎn)間的局部均值可計(jì)算如下:

(1)

考慮信號(hào)數(shù)據(jù)在極值點(diǎn)間一般是均勻變化的,其均值位置可視為在兩極值點(diǎn)的中點(diǎn)處,即tξ=(ti+ti+1)/2,這樣

(2)

同理,ti+1與ti+2時(shí)刻極值點(diǎn)間的局部均值mi+1為

(3)

這樣,mii+ti+1)/2為第i個(gè)極值點(diǎn)和第i+1個(gè)極值點(diǎn)信號(hào)的局部均值,而mi+1i+1+ti+2)/2為第i+1個(gè)極值點(diǎn)和第i+2個(gè)極值點(diǎn)信號(hào)的局部均值。現(xiàn)用2個(gè)相鄰局部均值mi和mi+1加權(quán)平均作為ti+1處極值點(diǎn)的局部均值mi(ti+1),即

(4)

(5)

利用式(4)可以得到信號(hào)各局部極值點(diǎn)處的局部均值m(0),m(t1),m(t2), …,m(tn),再使用3次樣條插值函數(shù)對(duì)其進(jìn)行曲線(xiàn)擬合得到局部均值m1(t),從而得到第1個(gè)IMF分量h1(t):

h1(t)=x(t)-m1(t)

(6)

這時(shí)h1(t)可能不滿(mǎn)足IMF定義,需要再次篩選:

h11(t)=h1(t)-m11(t)

(7)

重復(fù)式(7)直到滿(mǎn)足篩選停止準(zhǔn)則:

0.2

(8)

重復(fù)k次得到:

h1k(t)=h1(k-1)(t)-m1k(t)

(9)

此時(shí),h1k就是第1個(gè)IMF分量,表示為C1(t)=h1k(t),殘余信號(hào):

r1(t)=x(t)-C1(t)

(10)

以r1作為待分析信號(hào)重復(fù)以上步驟,直到所得剩余部分為單調(diào)信號(hào)或其值小于給定值時(shí),分解完畢,得到所有IMF分量C2(t),C3(t),C4(t), …,Cn(t)以及殘余分量Rn(t),這時(shí)信號(hào)可表示為所有IMF分量和殘余量之和:

(11)

實(shí)測(cè)信號(hào)中往往夾雜了大量的噪聲,而該噪聲是不需要分解的干擾量,不僅增加了EMD分解的層數(shù),使得分解時(shí)效性下降,而且過(guò)多的分解使得模態(tài)間累積泄漏加劇、邊界誤差不斷累加,嚴(yán)重時(shí)會(huì)使EMD分解失去實(shí)際的物理意義。

1.2 實(shí)測(cè)信號(hào)中噪聲的影響

某旋轉(zhuǎn)機(jī)械滾動(dòng)軸承振動(dòng)信號(hào)如圖1所示。分別對(duì)實(shí)測(cè)信號(hào)進(jìn)行EMD法和EMMD法分解,得到各IMF分量,由于篇幅限制,只取前4個(gè)分量,如圖2,圖3所示,相應(yīng)的希爾伯特邊際譜如圖4,圖5所示。

從圖3,圖5中可以看出EMMD法在分解精度、端點(diǎn)效應(yīng)抑制和計(jì)算耗時(shí)上均優(yōu)于EMD法,但是實(shí)測(cè)信號(hào)中大量噪聲同樣參加了計(jì)算過(guò)程,從分解出第1個(gè)IMF分量開(kāi)始,逐漸污染整個(gè)分解過(guò)程,尤其在50 Hz和120 Hz頻率附近,積累出能量很大的干擾量,影響了計(jì)算結(jié)果,需要對(duì)信號(hào)進(jìn)一步去噪。

圖1 實(shí)測(cè)滾動(dòng)軸承振動(dòng)信號(hào)示意

圖2 EMD法分解各IMF分量及殘余量曲線(xiàn)示意

2 MF和SVD相結(jié)合的去噪方法

2.1 形態(tài)濾波原理

與傳統(tǒng)的處理方法相比,數(shù)學(xué)形態(tài)學(xué)具有計(jì)算簡(jiǎn)單、實(shí)用性好、時(shí)延較小等特點(diǎn),為振動(dòng)信號(hào)去噪預(yù)處理提供了一種新的濾波方法[7]。

圖3 EMMD法分解各IMF分量及殘余量曲線(xiàn)示意

圖4 EMD算法邊際譜曲線(xiàn)示意

圖5 EMMD算法邊際譜曲線(xiàn)示意

數(shù)學(xué)形態(tài)學(xué)的基本運(yùn)算包括: 腐蝕運(yùn)算、膨脹運(yùn)算、開(kāi)運(yùn)算和閉運(yùn)算。設(shè)原始信號(hào)f(n)為定義在Df=(0, 1, 2, …,N-1)上的離散函數(shù),定義序列結(jié)構(gòu)元素g(n)為Dg=(0, 1, 2, …,M-1)上的離散函數(shù),且N≥M,則f(n)關(guān)于g(n)的的膨脹運(yùn)算和腐蝕運(yùn)算分別為

(f⊕g)(n)=max(f(n-m)+g(m))

(12)

(fΘg)(n)=min(f(n+m)-g(m))

(13)

f(n)關(guān)于g(n)的開(kāi)運(yùn)算和閉運(yùn)算分別為

(f°g)(n)=(fΘg⊕g)(n)

(14)

(f·g)(n)=(f⊕gΘg)(n)

(15)

其中m∈0, 1, 2, …,M-1;n∈0, 1, 2, …,N-1。

膨脹運(yùn)算除去負(fù)脈沖并平滑了正脈沖;腐蝕運(yùn)算除去正脈沖并平滑了負(fù)脈沖;開(kāi)運(yùn)算除去正脈沖并保留了負(fù)脈沖;閉運(yùn)算除去負(fù)脈沖并保留了正脈沖。為了同時(shí)去除正、負(fù)兩種噪聲,采用開(kāi)閉運(yùn)算的級(jí)聯(lián)形式,定義形態(tài)開(kāi)閉和閉開(kāi)濾波器:

Foc(f(n))=(f°g·g)(n)

(16)

Fco(f(n))=(f·g°g)(n)

(17)

為有效抑制脈沖干擾,本文采用開(kāi)閉和閉開(kāi)組合濾波器,定義:

y(n)=[Foc(f(n))+Fco(f(n))]/2

(18)

研究表明: 結(jié)構(gòu)元素越復(fù)雜,濾波能力越強(qiáng),計(jì)算量也越大[8]。為了確保EMMD法分解質(zhì)量,減少對(duì)信號(hào)中有用信息的損失,選擇扁平形結(jié)構(gòu)元素是最簡(jiǎn)單、也是最有效的方法,綜合考慮本文選擇扁平形結(jié)構(gòu)元素。

2.2 SVD原理

奇異值分解技術(shù)在故障診斷領(lǐng)域已有成功的應(yīng)用。已廣泛應(yīng)用于信號(hào)降噪或提取信號(hào)中的周期成分。文中利用SVD技術(shù),根據(jù)奇異值分布曲線(xiàn)確定降噪階次,進(jìn)行SVD降噪。

假設(shè)某故障機(jī)械系統(tǒng)的測(cè)試信號(hào)為xl(l=1, 2, 3, …),可由它重構(gòu)1個(gè)m×n維重構(gòu)相空間Dm[9]

若振動(dòng)信號(hào)中包含有一定的噪聲,則Dm可表示成Dm=D+W+V,其中D,W,V分別是光滑信號(hào)、故障信息、噪聲對(duì)應(yīng)的Dm中的軌跡矩陣,也可以將W和V視為對(duì)D矩陣的攝動(dòng)。現(xiàn)在僅知道Dm而不知道D和W,但是可以根據(jù)D和W的一些特點(diǎn)研究其奇異值情況對(duì)Dm做奇異值分解,Dm=USV′,U∈Rm×n,V′∈Rm×n,且UU′=I,VV′=I。S是m×n矩陣,對(duì)角線(xiàn)元素為λ1,λ2, …,λk,Dm的秩為k且k=min(m,n),一般取m?n,λ1≥λ2≥…≥λk稱(chēng)為Dm的奇異值,U和V分別表示左右奇異陣。

2.3 MF-SVD信號(hào)去噪方法

基于MF-SVD法和EMMD法的振動(dòng)信號(hào)故障特征提取方法,即先將實(shí)測(cè)現(xiàn)場(chǎng)信號(hào)進(jìn)行基于MF-SVD法去噪,結(jié)合形態(tài)濾波對(duì)振動(dòng)信號(hào)可以有效抑制脈沖噪聲干擾和奇異值分解針對(duì)隨機(jī)噪聲信號(hào)濾波效果明顯的優(yōu)點(diǎn),濾除隨機(jī)噪聲信號(hào)和局部強(qiáng)干擾噪聲,再進(jìn)行EMMD法分解,得到IMF分量,從而提取故障特征。

3 故障診斷實(shí)驗(yàn)

本文實(shí)驗(yàn)對(duì)象為某天然氣公司的制冷壓縮機(jī),是美國(guó)約克公司生產(chǎn)的RWB(II)-496型螺桿壓縮機(jī),單臺(tái)機(jī)組設(shè)計(jì)處理天然氣量為4 960 m3/h,壓縮介質(zhì)為丙烷,轉(zhuǎn)速為2 150 r/min,測(cè)點(diǎn)選擇在電機(jī)與螺旋桿之間連軸器靠近電機(jī)的支撐軸承上,軸承的型號(hào)為HRB—N205EM,參數(shù)見(jiàn)表1所列。

表1 軸承參數(shù)

在生產(chǎn)過(guò)程中發(fā)現(xiàn)測(cè)點(diǎn)在線(xiàn)振動(dòng)監(jiān)測(cè)裝置顯示振動(dòng)異常,檢修發(fā)現(xiàn)支撐軸承內(nèi)圈有中度損傷,其特征頻率應(yīng)為f=256.11 Hz。測(cè)得振動(dòng)信號(hào)如圖6所示,其頻譜如圖7所示,根據(jù)頻譜圖不能判斷出具體故障原因。

圖6 振動(dòng)信號(hào)時(shí)域波形曲線(xiàn)示意

圖7 振動(dòng)信號(hào)頻譜曲線(xiàn)示意

對(duì)信號(hào)直接進(jìn)行EMMD法分解,其邊際譜分析如圖8所示,根據(jù)邊際譜曲線(xiàn)所示,雖然頻率在260 Hz附近有較高的振幅,但在50 Hz和100 Hz附近同樣出現(xiàn)了較大的振幅,形成很大的干擾,而且邊界抖動(dòng)幅度大,這都說(shuō)明在直接EMMD法分解過(guò)程中,由于混有現(xiàn)場(chǎng)強(qiáng)噪聲,使整個(gè)分解過(guò)程受到污染,致使出現(xiàn)端點(diǎn)效應(yīng)和模態(tài)混疊問(wèn)題,影響了邊際譜的刻畫(huà),對(duì)于滾動(dòng)故障診斷來(lái)說(shuō)無(wú)法準(zhǔn)確判斷其故障類(lèi)型。

圖8 直接EMMD法分解邊際譜曲線(xiàn)示意

采用MF-SVD法去噪和EMMD法相結(jié)合的故障特征提取方法,對(duì)信號(hào)進(jìn)行處理。去噪后的信號(hào)如圖9所示,其IMF分量如圖10所示,其邊際譜如圖11所示。可以看出,邊際譜曲線(xiàn)頻率僅在255.6 Hz處幅度最大,其他頻率干擾消除,去噪效果明顯,有效地抑制了邊界效應(yīng)和模態(tài)混疊問(wèn)題,可以準(zhǔn)確判斷出故障類(lèi)型為滾動(dòng)軸承內(nèi)圈損傷。

圖9 MF-SVD去噪后信號(hào)時(shí)域曲線(xiàn)示意

圖10 EMMD法分解各IMF分量及 殘余信號(hào)曲線(xiàn)示意

4 結(jié)束語(yǔ)

筆者針對(duì)現(xiàn)場(chǎng)實(shí)測(cè)信號(hào)摻雜強(qiáng)噪聲影響EMMD法分解質(zhì)量和結(jié)果的問(wèn)題,提出了一種基于MF-SVD法與EMMD法分解相結(jié)合的故障特征提取方法。該方法先對(duì)實(shí)測(cè)現(xiàn)場(chǎng)信號(hào)采用MF-SVD法去噪,濾除隨機(jī)噪聲信號(hào)和局部強(qiáng)干擾噪聲,彌補(bǔ)了形態(tài)濾波去噪方法對(duì)于隨機(jī)噪聲信號(hào)處理的不足和奇異值分解對(duì)于脈沖干擾噪聲處理的不足,然后進(jìn)行EMMD法分解,得到IMF分量,從而提取故障特征。

圖11 去噪后經(jīng)EMMD法分解邊際譜曲線(xiàn)

實(shí)驗(yàn)結(jié)果表明: 該方法能夠去除隨機(jī)噪聲和局部強(qiáng)干擾噪聲,效果明顯,能夠有效、準(zhǔn)確地提取滾動(dòng)軸承內(nèi)圈損傷的故障特征,可抑制端點(diǎn)效應(yīng),有著良好的運(yùn)算速度和精確度。

[1] 王宏禹.信號(hào)處理相關(guān)理論綜合與統(tǒng)一法[M].北京: 國(guó)防工

業(yè)出版社,2005.

[2] Huang N E, Shen Z, Long S R. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J]. Proceedings of the Royal Society of London, 1998, 45(04): 903-995.

[3] Wu Zhaohua, Huang N E. Ensemble Empirical Mode Decomposition: a Noise-assisted Data Analysis Method[J]. Advances in Adaptive Data Analysis, 2009(01): 1-41.

[4] 沈路,楊富春,周曉軍,等.基于改進(jìn)EMD與形態(tài)濾波的齒輪故障特征提取[J].振動(dòng)與沖擊,2010,29(03): 154-157.

[5] 曾作欽.基于奇異值分解的信號(hào)處理方法及其在機(jī)械故障診斷中的應(yīng)用[D].廣州: 華南理工大學(xué),2011.

[6] 屈梁生,張西寧,沈玉娣.機(jī)械故障診斷理論與方法[M].西安: 西安交通大學(xué)出版社,2009.

[7] 李舜酩,李香蓮.振動(dòng)信號(hào)的現(xiàn)代分析技術(shù)與應(yīng)用[M].北京: 國(guó)防工業(yè)出版社,2008.

[8] 李兵,張培林,米雙山.機(jī)械故障信號(hào)的數(shù)學(xué)形態(tài)學(xué)分析與智能分類(lèi)[M].北京: 國(guó)防工業(yè)出版社,2011.

[9] 韓清凱,于曉光.基于振動(dòng)分析的現(xiàn)代機(jī)械故障診斷原理及應(yīng)用[M].北京: 科學(xué)出版社,2010.

[10] 劉佳寶,梁奕,徐漫江.一種過(guò)程數(shù)據(jù)趨勢(shì)特征提取方法.化工自動(dòng)化及儀表,2012,39(07): 850-853.

[11] 劉繼承,聶品磊,楊宏宇,等.基于形態(tài)濾波和HHT的滾動(dòng)軸承故障特征提取.化工自動(dòng)化及儀表,2014,41(05): 529-532,562.

Study on Fault Characteristics Extraction for MF-SVD Based Rolling Bearing Vibration Signal

Sachsla

(College of Petroleum Engineering, Northeast Petroleum University, Daqing, 163318, China)

In terms of calculation precision, restraining boundary effect and computation time, extremum field mean mode decomposition(EMMD) has obvious advantage comparing to empirical mode decomposition(EMD) and adaptive time varying filter decomposition(ATVFD). It can effectively extract machinery vibration signal fault characteristics. As the measured field signals are often mixed lots of noise, which influences EMMD decomposition quality seriously. Aiming at the problem, a new extracting fault feature method which is based on morphological filters-singular value decomposition (MF-SVD) denoising method, combined with EMMD, is put forward. The experimental results show fault features of rolling bearing inner ring damage can be effectively and accurately extracted. It has a good calculation speed and accuracy.

extremum field mean mode decomposition; morphological filters; singular value decomposition; feature extraction

黑龍江省長(zhǎng)江學(xué)者后備計(jì)劃資助(2012CJHB005)。

薩其日拉(1994—),男,蒙古族人,就讀于東北石油大學(xué)石油工程學(xué)院,研究方向?yàn)橛吞锕に嚰霸O(shè)備技術(shù)。

TP273

A

1007-7324(2017)03-0031-06

稿件收到日期: 2017-02-27,修改稿收到日期: 2017-04-27。

猜你喜歡
特征提取振動(dòng)故障
振動(dòng)的思考
振動(dòng)與頻率
故障一點(diǎn)通
基于Gazebo仿真環(huán)境的ORB特征提取與比對(duì)的研究
電子制作(2019年15期)2019-08-27 01:12:00
中立型Emden-Fowler微分方程的振動(dòng)性
一種基于LBP 特征提取和稀疏表示的肝病識(shí)別算法
奔馳R320車(chē)ABS、ESP故障燈異常點(diǎn)亮
故障一點(diǎn)通
江淮車(chē)故障3例
基于MED和循環(huán)域解調(diào)的多故障特征提取
主站蜘蛛池模板: 久操中文在线| 亚洲国产精品人久久电影| 国产成人精品高清不卡在线| 国产成人8x视频一区二区| 欧美一级色视频| 国产又黄又硬又粗| 午夜电影在线观看国产1区| 国模在线视频一区二区三区| 色偷偷男人的天堂亚洲av| 一级毛片在线免费视频| 天堂网国产| 亚洲中文字幕国产av| 狠狠做深爱婷婷综合一区| 99色亚洲国产精品11p| 欧美精品v日韩精品v国产精品| 久久久黄色片| 91毛片网| 色欲国产一区二区日韩欧美| 九九热免费在线视频| 一本大道视频精品人妻| 特级毛片免费视频| 58av国产精品| 性喷潮久久久久久久久| 久久先锋资源| 亚洲成a人在线观看| 欧美精品高清| 久久这里只有精品23| 99福利视频导航| 伊人无码视屏| 精品一區二區久久久久久久網站| 97一区二区在线播放| 日本亚洲欧美在线| 久久综合干| 人妻出轨无码中文一区二区| 看国产毛片| 日本免费福利视频| 国产成人一区| 国语少妇高潮| 国产黑丝视频在线观看| 伊人蕉久影院| 亚洲一区二区三区在线视频| 亚洲第七页| 亚洲天堂777| 亚洲无线国产观看| 中文字幕 日韩 欧美| 啪啪啪亚洲无码| 91香蕉国产亚洲一二三区 | 欧美日韩在线第一页| 欧美在线精品怡红院| 国产免费黄| 性色在线视频精品| 国产欧美视频在线| 国产成人一区二区| 1769国产精品视频免费观看| 91久久国产热精品免费| 五月激情婷婷综合| 黄色网在线| 国产流白浆视频| 人妻中文久热无码丝袜| 国产91全国探花系列在线播放| 国产日本欧美亚洲精品视| 色婷婷综合激情视频免费看 | 国产亚洲精品资源在线26u| 日韩欧美中文在线| 在线免费观看AV| 日本欧美午夜| 中国毛片网| 国产高清精品在线91| 亚洲精品无码抽插日韩| www.91在线播放| 噜噜噜久久| 伊大人香蕉久久网欧美| 91成人免费观看| 国产精品久久久精品三级| 亚洲嫩模喷白浆| 亚洲黄色片免费看| 天天婬欲婬香婬色婬视频播放| 伊人久久大香线蕉aⅴ色| 欧美成人综合在线| 久久久噜噜噜久久中文字幕色伊伊| 国产极品粉嫩小泬免费看| 国产高清国内精品福利|