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

基于TVD和MSB的滾動(dòng)軸承故障特征提取

2019-06-13 09:57:10朱丹宸張永祥朱群偉
振動(dòng)與沖擊 2019年8期
關(guān)鍵詞:特征提取故障信號(hào)

朱丹宸 ,張永祥,趙 磊,朱群偉

(海軍工程大學(xué) 動(dòng)力工程學(xué)院,武漢 430033)

滾動(dòng)軸承在旋轉(zhuǎn)機(jī)械中得到廣泛使用,其工作狀態(tài)直接影響著整個(gè)機(jī)器設(shè)備的運(yùn)行穩(wěn)定性,因此,對滾動(dòng)軸承進(jìn)行狀態(tài)監(jiān)測和故障診斷就具有重要意義。就目前而言,振動(dòng)信號(hào)分析是軸承故障診斷的常用手段。滾動(dòng)軸承故障一般包括內(nèi)圈故障、外圈故障以及滾動(dòng)體故障,這些故障產(chǎn)生的振動(dòng)信號(hào)往往是非平穩(wěn)、非線性的,同時(shí)周期性故障沖擊一般都淹沒在強(qiáng)背景噪聲之中[1],如何降低原始振動(dòng)信號(hào)中的噪聲成分,提取出軸承故障特征就成了主要研究內(nèi)容。

TVD(Total Variation Denoising)是一種基于最優(yōu)問題的降噪方法,最早被用于去除圖像中的噪聲成分[2],其本質(zhì)上是一種信號(hào)的稀疏表示方法,隨后,TVD方法得到了廣泛的研究并在一維信號(hào)處理中得到了運(yùn)用。由于TVD是通過將代價(jià)函數(shù)最小化來獲得最終輸出,文獻(xiàn)[3]提出利用Majorization-Minimization(MM)算法對TVD進(jìn)行求解從而構(gòu)成TV-MM算法,然而,該算法中,參數(shù)λ的選取對最終濾波結(jié)果有著較大影響。Zhang等[4]將變分模態(tài)分解以及TV-MM算法進(jìn)行結(jié)合,利用加權(quán)峭度求解最優(yōu)參數(shù)λ,并證明該方法優(yōu)于一些經(jīng)典的去噪算法。Yi等[5]提出了一種二階TVD算法,將懲罰函數(shù)引入算法之中,通過排列熵和相關(guān)系數(shù)選取參數(shù)λ,提高了降噪效果。

近來,Alwodai等[6]利用MSB(Modulation Signal Bispectrum)分析電機(jī)中的電流信號(hào),展現(xiàn)出比傳統(tǒng)功率譜更好的頻譜特性,實(shí)現(xiàn)了電機(jī)軸承的故障診斷。文獻(xiàn)[7]同樣利用電流信號(hào),通過MSB分析,實(shí)現(xiàn)了齒輪磨損狀況的監(jiān)測。Rehab等[8]利用MSB從包絡(luò)信號(hào)中提取出故障特征,體現(xiàn)出比傳統(tǒng)包絡(luò)分析更強(qiáng)的抗噪與故障監(jiān)測能力。文獻(xiàn)[9]提出了一種更為穩(wěn)定的MSB解調(diào)方法,實(shí)現(xiàn)了滾動(dòng)軸承的故障特征提取,且分析效果優(yōu)于傳統(tǒng)的Kurtogram算法。

受到強(qiáng)背景噪聲的影響和信號(hào)自身特點(diǎn)的制約,TVD能夠消除部分噪聲,但同時(shí)也有一些殘留噪聲影響故障特征提取效果。本文提出將TVD與MSB相結(jié)合,利用MSB分析TVD處理后的濾波信號(hào),進(jìn)一步抑制噪聲的干擾,提取出軸承故障特征。首先,二階TVD被用于處理采樣得到的振動(dòng)信號(hào),利用包絡(luò)譜相關(guān)峭度選取參數(shù)λ,獲得較優(yōu)濾波結(jié)果;然后利用MSB分析濾波后的信號(hào),通過故障特征頻率成分占比p選取最佳的5個(gè)載波頻率切片進(jìn)行平均得到復(fù)合切片譜,提取出軸承故障特征;最后,通過分析復(fù)合切片譜,判斷故障類型。仿真和實(shí)驗(yàn)分析表明,該方法能夠有效抑制隨機(jī)噪聲,提高故障特征提取效果。

1 二階全變分去噪

1.1 基本原理

全變分去噪可以看成是一個(gè)數(shù)值優(yōu)化過程,包含二次數(shù)據(jù)保真項(xiàng)和凸正則化項(xiàng),常用的全變分過程的基礎(chǔ)是通過一階或者二階差分實(shí)現(xiàn)對原始信號(hào)的稀疏表示。本文選用二階差分定義信號(hào)的全變分。假設(shè)一維信號(hào)x(n),(0≤n≤N-1),定義信號(hào)x(n)的二階全變分為

(1)

式中:D2x為二階差分,其中D2可以用矩陣形式表示為

(2)

‖·‖p為p階范數(shù),且‖x‖p可以表示為

(3)

當(dāng)p=1以及p=2時(shí),信號(hào)的一階與二階范數(shù)可分別表示為

‖x‖1=|x1|+|x2|+…+|xN|

(4)

(5)

假設(shè)信號(hào)y(n)由有效成分x(n)與白噪聲w(n)組成,即

y(n)=x(n)+w(n)

(6)

由此,二階全變分去噪模型可以通過一個(gè)選優(yōu)過程表示

(7)

式中:F(x)為目標(biāo)函數(shù);λ>0為正則化參數(shù)。

在獲得目標(biāo)函數(shù)之后,本文利用MM算法求解式(7),結(jié)果可以表示為

(8)

利用式(8)就可以實(shí)現(xiàn)信號(hào)的全變分去噪,但通過觀察式(7)和式(8)可知,正則化參數(shù)λ對降噪過程會(huì)產(chǎn)生一定的影響,當(dāng)λ→0時(shí),降噪后的結(jié)果x將無限趨近于原始信號(hào)y,當(dāng)λ→∞時(shí),信號(hào)的保真度將會(huì)下降。因此,如何選取最為合適的參數(shù)λ是接下來將要研究的內(nèi)容。

1.2 最優(yōu)參數(shù)的選取方法

軸承故障信號(hào)在時(shí)域范圍內(nèi)包含周期性的沖擊特征,所以相關(guān)峭度(Correlated Kurtosis,CK)作為一種衡量信號(hào)中沖擊特性的指標(biāo),由于它既保留了峭度的特性,也包含了相關(guān)函數(shù)的特性,常被用于軸承故障診斷之中,其計(jì)算公式為[10]

(9)

式中:xn為信號(hào)序列;N為信號(hào)的長度;T為感興趣信號(hào)的長度;M為偏移的周期個(gè)數(shù)。同時(shí),軸承故障沖擊特性在頻域范圍內(nèi)會(huì)表現(xiàn)為故障特征頻率及其倍頻。由此,本文利用頻域相關(guān)峭度自適應(yīng)選取全變分去噪中的參數(shù)λ。頻域相關(guān)峭度可以通過原信號(hào)的包絡(luò)譜定義為

(10)

式中:E(x)n為信號(hào)序列xn的包絡(luò)譜;此時(shí)T選取為故障特征頻率。當(dāng)信號(hào)的包絡(luò)譜中包含明顯的故障特征頻率成分時(shí),頻域相關(guān)峭度較大,相反,當(dāng)故障特征頻率成分沒有在信號(hào)的包絡(luò)譜中得到明顯體現(xiàn)時(shí),頻域相關(guān)峭度較小。

2 調(diào)制信號(hào)雙譜

2.1 基本原理

MSB作為一種考慮邊頻帶的雙譜分析方法能對調(diào)制信號(hào)進(jìn)行有效分析,該方法可以較好地抑制隨機(jī)噪聲和非周期成分的干擾,清晰反映出信號(hào)中的調(diào)制成分。MSB在頻域范圍內(nèi),利用原始信號(hào)x(t)的離散傅里葉變換形式X(f)可以表示為

BMS(fc,fx)=
E〈X(fc+fx)X(fc-fx)X*(fc)X*(fc)〉

(11)

式中:E〈·〉為數(shù)學(xué)期望;fc為載波頻率;fx為調(diào)制頻率。對MSB進(jìn)行歸一化可得

(12)

2.2 最優(yōu)切片譜選取

為了進(jìn)行有效的故障特征提取,選擇合適的分析頻段是十分重要的,在利用MSB進(jìn)行分析時(shí),也即是選取合適的fc切片位置,而一般來說,故障特征較為明顯的切片位置有多個(gè),選取多個(gè)切片可以綜合各切片處的故障特征信息,同時(shí)也可以減小存在于單個(gè)切片中的干擾成分。由此,本文提出利用前三階故障特征頻率占比p來選取最佳的5個(gè)fc切片,并對其進(jìn)行平均得到復(fù)合切片譜,從而提取出軸承故障特征。

設(shè)f1,f2,f3分別為軸承故障的特征頻率及其二倍頻和三倍頻,則在fc切片位置,特征頻率占比p可以定義為

(13)

(14)

3 故障特征提取流程

基于TVD與MSB的滾動(dòng)軸承故障特征提取步驟如下,流程如圖1所示:

步驟1獲取振動(dòng)信號(hào),在0.1~5內(nèi),間隔為0.1選取λ并計(jì)算相應(yīng)的TVD處理結(jié)果;

步驟2以最大包絡(luò)譜相關(guān)峭度為準(zhǔn)則確定參數(shù)λ并計(jì)算最優(yōu)去噪結(jié)果;

步驟3對步驟2得到的去噪信號(hào)進(jìn)行MSB分析;

步驟4根據(jù)特征頻率成分占比p選取MSB中的5個(gè)最優(yōu)載波頻率切片;

步驟5對步驟4中得到的5個(gè)切片譜進(jìn)行平均,得到復(fù)合切片譜,提取出軸承故障特征,判斷故障形式。

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

為了驗(yàn)證本文前一部分所述故障特征提取方法是正確而又有效的,先利用仿真信號(hào)進(jìn)行研究分析,以滾動(dòng)軸承內(nèi)圈故障為例構(gòu)造仿真信號(hào)

(15)

圖1 滾動(dòng)軸承故障特征提取流程圖Fig.1 Flow chart of fault feature extraction for rolling element bearings

式中:軸的轉(zhuǎn)頻fr為42 Hz;Ai為以1/fr為周期的幅值調(diào)制;n(t)為隨機(jī)白噪聲;r為系統(tǒng)的阻尼系數(shù)且設(shè)定為0.05;兩個(gè)相鄰沖擊的間隔T為1/185 s;ti為第i各周期內(nèi),由滾動(dòng)體滑移引起的延遲,ti=0.01T;B(t)為諧波分量;fm=100;φA=π/2;φw=0;A0,B0,C0分別為2,0.5,0.5;系統(tǒng)自然頻率fn為3 200 Hz;采樣頻率fs為32 768 Hz;仿真時(shí)間為1 s。

為了分析TVD的去噪效果,將隨機(jī)白噪聲分別設(shè)為0,-5 db與-10 db,分別得到仿真信號(hào)1、仿真信號(hào)2和仿真信號(hào)3,其時(shí)域波形如圖2所示。由于受到噪聲的干擾,周期性沖擊成分在三組仿真信號(hào)中都不能得到體現(xiàn),且所添加白噪聲越強(qiáng),干擾越明顯。

利用TVD處理仿真信號(hào),根據(jù)圖3(a)、圖3(c)和圖3(e)給出的三組λ值變化曲線,選取的最優(yōu)參數(shù)λ分別為0.7,1.3和2.1,得到的降噪結(jié)果如圖3(b)、圖3(d)和圖3(f)所示。對比圖3和圖2,雖然經(jīng)過TVD處理后的時(shí)域波形中能夠看出噪聲的減小,但周期性沖擊成分仍然難以得到體現(xiàn)。

圖4對比了三組仿真信號(hào)在降噪前后的包絡(luò)譜,分析仿真信號(hào)1,降噪前后的包絡(luò)譜中都能明顯識(shí)別沖擊頻率185 Hz以及其倍頻370 Hz,555 Hz,軸承故障特征得到了有效提取,且圖4(b)中,圍繞故障特征頻率及其倍頻,間隔為轉(zhuǎn)頻的調(diào)制邊頻帶更為明顯,說明利用TVD能夠進(jìn)一步降低噪聲提高故障特征提取效果;分析仿真信號(hào)2,原始信號(hào)的包絡(luò)譜無法識(shí)別故障特征頻率,通過TVD處理后,降噪信號(hào)的包絡(luò)譜能夠識(shí)別故障特征頻率及其倍頻成分,但白噪聲對故障特征提取的干擾是較為明顯的;分析仿真信號(hào)3,降噪前后信號(hào)的包絡(luò)譜都不能識(shí)別故障特征頻率及其倍頻。綜上所述,當(dāng)信號(hào)中的背景噪聲較強(qiáng)時(shí),TVD的降噪效果并不理想,需要利用其他分析手段進(jìn)一步加以處理。

圖2 仿真信號(hào)的時(shí)域波形Fig.2 Time domain waveform of simulation signals

圖3 TVD處理后的結(jié)果Fig.3 Time domain waveform after TVD processing

圖4 原始仿真信號(hào)與去噪信號(hào)的包絡(luò)譜Fig.4 Envelope spectrum of original simulation signal and denoised simulation signal

由于在強(qiáng)背景噪聲下,僅利用TVD處理仿真信號(hào)不能實(shí)現(xiàn)較好的故障特征提取,接下來將利用MSB進(jìn)一步分析圖3(d)和圖3(f)中的信號(hào)。對于圖3(d)中的信號(hào),通過比較載波頻率各切片的p值,選取288 Hz,658 Hz,3 004 Hz,2 920 Hz和472 Hz處(以p值從大到小進(jìn)行排序)的切片構(gòu)成復(fù)合切片譜,結(jié)果如圖5(a)所示;對于圖3(f)中的信號(hào),選取472 Hz,288 Hz,658 Hz,2 920 Hz和2 776 Hz處的切片構(gòu)成復(fù)合切片譜,結(jié)果如圖5(b)所示。

圖5中能夠明顯識(shí)別故障特征頻率及其倍頻成分,白噪聲的干擾得到了明顯的抑制,且故障特征提取效果要明顯優(yōu)于圖4(d)和圖4(f)的結(jié)果。

圖5 本文算法處理后的結(jié)果Fig.5 Results of the proposed method

為了進(jìn)行對比,利用文獻(xiàn)[11]提出的Fast Kurtogram算法分析仿真信號(hào)2,結(jié)果如圖6所示。快速峭度圖確定的最佳分析頻帶中心頻率為16 213 Hz,帶寬為341 Hz,其平方包絡(luò)譜如圖6(b)所示,從圖6(b)可知,不能明顯識(shí)別故障特征頻率及其倍頻成分,說明對于該仿真信號(hào),快速峭度圖方法失效。由此可見,本文的算法在抑制白噪聲,提取弱故障沖擊方面具有一定的優(yōu)勢。

圖6 Fast Kurtogram算法分析仿真信號(hào)2的結(jié)果Fig.6 Analysis results of simulation signal 2 by Fast Kurtogram

5 實(shí)測信號(hào)分析

為了驗(yàn)證所提出的方法在滾動(dòng)軸承內(nèi)圈和外圈故障特征提取中的有效性,使用實(shí)驗(yàn)室滾動(dòng)軸承故障模擬平臺(tái)(見圖7)進(jìn)行軸承故障模擬。測點(diǎn)選取在機(jī)匣外側(cè)以增加振動(dòng)信號(hào)的復(fù)雜程度,測量方向?yàn)閺较颉?/p>

圖7 實(shí)驗(yàn)平臺(tái)Fig.7 The test rig

5.1 軸承內(nèi)圈故障分析

軸承內(nèi)圈故障信號(hào)測量時(shí)所用的軸承為6 010,該軸承的內(nèi)徑為50 mm,外徑為80 mm,滾動(dòng)體直徑為9 mm,滾動(dòng)體個(gè)數(shù)為13,接觸角α=0°。測試時(shí)的轉(zhuǎn)速為3 000 r/min,采樣頻率為32 768 Hz,通過計(jì)算可得故障特征頻率的理論值為370.01 Hz。

為了提高分析效果,先對測量得到的振動(dòng)信號(hào)進(jìn)行標(biāo)準(zhǔn)化處理。圖8(a)為標(biāo)準(zhǔn)化處理后,軸承內(nèi)圈故障信號(hào)的時(shí)域波形,且從圖8(b)的頻譜中無法識(shí)別軸承故障特征頻率及其倍頻成分。

圖8 軸承內(nèi)圈故障信號(hào)Fig.8 The inner ring fault signal

利用本文的算法處理該內(nèi)圈故障信號(hào),由圖9(a)中的λ值變化曲線,可以確定TVD處理時(shí)的參數(shù)λ為1.6,得到的降噪后時(shí)域波形如圖9(b)所示。從圖9(b)可知,軸承故障產(chǎn)生的周期性沖擊得到了一定的體現(xiàn),之后利用MSB分析降噪后的故障信號(hào),比較各載波頻率切片的p值,選取2 154 Hz,5 292 Hz,3 990 Hz,4 040 Hz和2 522 Hz(以p值從大到小進(jìn)行排序)處的切片構(gòu)成復(fù)合切片譜,得到的結(jié)果如圖9(c)所示。

圖9(c)中軸承的轉(zhuǎn)頻50 Hz及其倍頻,軸承內(nèi)圈故障特征頻率為368 Hz(與理論值相接近)以及其倍頻736 Hz,1 104 Hz能夠得到明顯識(shí)別,圍繞內(nèi)圈故障特征頻率及其倍頻,間隔為轉(zhuǎn)頻的調(diào)制邊頻帶也能得到體現(xiàn),噪聲對分析的干擾得到明顯抑制。通過與故障特征頻率的理論值相對比,可以判斷該軸承存在內(nèi)圈故障。

圖9 本文算法處理軸承內(nèi)圈故障信號(hào)的結(jié)果Fig.9 Results of inner ring fault signal by proposed method

為了進(jìn)行對比,與仿真分析類似,利用文獻(xiàn)[11]提出的Fast Kurtogram算法分析該故障信號(hào),結(jié)果如圖10所示。快速峭度圖確定的最佳分析頻帶中心頻率為341 Hz,帶寬為682 Hz,其平方包絡(luò)譜如圖10(b)所示。從圖10(b)可知,不能明顯識(shí)別故障特征頻率及其倍頻成分,快速峭度圖方法失效。由此可見,本文的算法在提取弱故障沖擊方面具有一定的優(yōu)勢。

圖10 Fast Kurtogram算法處理軸承內(nèi)圈故障信號(hào)的結(jié)果Fig.10 Results of inner ring fault signal by Fast Kurtogram

5.2 軸承外圈故障分析

軸承外圈故障信號(hào)測量時(shí)所用的軸承為7 010 C,該軸承的內(nèi)徑為50 mm,外徑為80 mm,滾動(dòng)體直徑為8.7 mm,滾動(dòng)體個(gè)數(shù)為19,接觸角a=15°。測試時(shí)的轉(zhuǎn)速為3 000 r/min,采樣頻率為32 768 Hz,通過計(jì)算可得故障特征頻率的理論值為413.6 Hz。圖11(a)為標(biāo)準(zhǔn)化處理后,軸承外圈故障信號(hào)的時(shí)域波形,受背景噪聲的干擾,周期性故障沖擊不能得到識(shí)別,同時(shí),圖11(b)的頻譜中,軸承故障特征頻率及其倍頻成分也沒有得到體現(xiàn)。

圖11 軸承外圈故障信號(hào)Fig.11 The outer ring fault signal

利用本文的算法處理該外圈故障信號(hào),由圖12(a)中的變化曲線可得TVD處理時(shí)選取的參數(shù)λ為2.5,得到的降噪后時(shí)域波形如圖12(b)所示。相比原始故障信號(hào),從圖12(b)可知,背景噪聲得到了明顯減少,之后利用MSB分析降噪后的故障信號(hào),選取載波頻率2 502 Hz,4 168 Hz,3 334 Hz,4 584 Hz和2 918 Hz處的切片構(gòu)成復(fù)合切片譜,得到的結(jié)果如圖12(c)所示。

圖12(c)中軸承外圈故障特征頻率416 Hz(與理論值相接近)以及其倍頻832 Hz,1 248 Hz能夠得到明顯識(shí)別,噪聲的干擾得到明顯抑制。通過與故障特征頻率的理論值相對比,可以判斷該軸承存在外圈故障。

圖12 本文算法處理軸承外圈故障信號(hào)的結(jié)果Fig.12 Results of outer ring fault signal by proposed method

進(jìn)行對比,圖13為Fast Kurtogram算法分析該故障信號(hào)的結(jié)果,快速峭度圖確定的最佳分析頻帶中心頻率為14 336 Hz,帶寬為1 365 Hz,其平方包絡(luò)譜如圖13(b)所示,由圖13(b)可知,外圈故障特征頻率416 Hz及其三倍頻1 248 Hz的幅值十分突出,但受到噪聲干擾,故障特征頻率的二倍頻832 Hz能得到識(shí)別但不明顯,分析效果欠佳,通過比較,體現(xiàn)了本文算法的有效性。

圖13 Fast Kurtogram算法處理軸承外圈故障信號(hào)的結(jié)果Fig.13 Results of outer ring fault signal by Fast Kurtogram

6 結(jié) 論

通過仿真和實(shí)驗(yàn)室實(shí)測信號(hào)驗(yàn)證表明,將TVD和MSB結(jié)合進(jìn)行軸承故障特征提取是可行的。本文得到的結(jié)論主要有:

(1)TVD能夠有效突出信號(hào)中的沖擊成分,減少噪聲干擾,但受到信號(hào)本身特點(diǎn)的制約與強(qiáng)背景噪聲的影響,TVD的效果并不是十分理想。通過MSB進(jìn)行解調(diào)可以有效抑制信號(hào)中的隨機(jī)噪聲干擾。將MSB與TVD相結(jié)合,可以利用兩種方法的優(yōu)勢,提高故障特征提取的有效性。相比于經(jīng)典的快速峭度圖方法,本文提出的算法在抑制背景噪聲,提取弱故障沖擊特征方面更為有效。

(2)針對TVD和MSB解調(diào)時(shí)的不確定性,利用頻域相關(guān)峭度選取TVD中的參數(shù)λ,通過特征頻率成分占比p選取最佳切片位置,確保了結(jié)果的準(zhǔn)確性。

(3)本文所提出的方法也存在一定的缺陷,TVD作為一種基于信號(hào)稀疏表示的方法,其分析結(jié)果會(huì)受到信號(hào)自身特點(diǎn)的影響,因此,如何根據(jù)信號(hào)自身特點(diǎn),選取更為合適的稀疏表示方法,獲得更理想的降噪效果,值得進(jìn)一步研究。

猜你喜歡
特征提取故障信號(hào)
信號(hào)
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點(diǎn)通
基于Gazebo仿真環(huán)境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
電子制作(2018年11期)2018-08-04 03:25:42
一種基于LBP 特征提取和稀疏表示的肝病識(shí)別算法
奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
基于LabVIEW的力加載信號(hào)采集與PID控制
故障一點(diǎn)通
江淮車故障3例
主站蜘蛛池模板: AV色爱天堂网| 91破解版在线亚洲| 国产精品久线在线观看| 国产精品刺激对白在线| 日韩小视频在线观看| 欧美精品成人| 欧美福利在线| 午夜毛片免费看| 白浆视频在线观看| 久久综合丝袜日本网| 波多野结衣视频一区二区| 欧美劲爆第一页| 国产一级毛片高清完整视频版| 色综合中文综合网| 亚洲欧美另类色图| 国产女人水多毛片18| 四虎成人精品在永久免费| 日韩欧美国产综合| 在线视频亚洲色图| 一本一道波多野结衣一区二区 | 污视频日本| www.狠狠| 99久久国产综合精品2020| 欧类av怡春院| 亚洲av无码人妻| 高清无码手机在线观看| 日韩无码真实干出血视频| 国产精品久久久久久久伊一| 色135综合网| 手机在线国产精品| 青青草a国产免费观看| 超清无码一区二区三区| 久久久久青草线综合超碰| 国产拍揄自揄精品视频网站| 香蕉国产精品视频| 在线观看精品国产入口| 国产h视频免费观看| 亚洲综合中文字幕国产精品欧美 | 成人福利在线观看| 久久婷婷六月| 91视频首页| 亚洲第一区在线| 欧美啪啪精品| 999精品色在线观看| 精品91自产拍在线| 91小视频在线观看| 成人一区专区在线观看| 白浆视频在线观看| 国产网站免费观看| 伊伊人成亚洲综合人网7777| 亚洲 欧美 中文 AⅤ在线视频| 手机成人午夜在线视频| 国产靠逼视频| 国产精品九九视频| 乱码国产乱码精品精在线播放| 欧美一区二区啪啪| 五月婷婷综合网| 国产精品部在线观看| AV片亚洲国产男人的天堂| 亚洲中文字幕无码爆乳| 亚洲bt欧美bt精品| 先锋资源久久| 国产大全韩国亚洲一区二区三区| 国产欧美精品一区二区 | 香蕉久久国产超碰青草| 91人人妻人人做人人爽男同| 青青青草国产| 国产主播福利在线观看| 欧美日韩国产成人高清视频| 97精品国产高清久久久久蜜芽| 欧美国产精品不卡在线观看| 人妻一本久道久久综合久久鬼色| 首页亚洲国产丝袜长腿综合| 无码日韩人妻精品久久蜜桃| 亚洲中文在线看视频一区| 狼友视频一区二区三区| 亚洲精品视频免费观看| 手机在线免费不卡一区二| 久久人人爽人人爽人人片aV东京热 | 国产十八禁在线观看免费| 亚洲国产成人超福利久久精品| 欧美成一级|