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

基于EEMD特征提取的滾動(dòng)軸承故障診斷

2022-06-22 07:24:10高淑芝李天池
關(guān)鍵詞:特征提取故障診斷模態(tài)

高淑芝, 李天池

(沈陽化工大學(xué) 信息工程學(xué)院, 遼寧 沈陽 110142)

滾動(dòng)軸承是關(guān)鍵機(jī)械基礎(chǔ)件之一,其失效是引起機(jī)械故障的主要原因,因此滾動(dòng)軸承的故障診斷技術(shù)極為重要[1].滾動(dòng)軸承失效的形式有疲勞失效、磨損失效、腐蝕失效、溫升失效、振動(dòng)失效、耦合失效,滾動(dòng)軸承故障診斷過程包括信息獲取、特征提取和模式識(shí)別[2].信息獲取的方法有很多,本篇采用振動(dòng)測(cè)量完成信號(hào)采集.

滾動(dòng)軸承的振動(dòng)信號(hào)通常具有非線性非平穩(wěn)性的特點(diǎn),對(duì)其處理方法有快速傅里葉變換(FFT)、小波變換等,但是這些方法都不能準(zhǔn)確地反應(yīng)信號(hào)本質(zhì)特征[3].經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)是針對(duì)該信號(hào)進(jìn)行處理的一種有效方法.該方法具有去局部時(shí)變特征的自適應(yīng)性,可將信號(hào)分解成若干個(gè)本征模態(tài)函數(shù)(IMF)之和,降低信號(hào)間特征信息的耦合.Pan 等[4]提出基于經(jīng)驗(yàn)?zāi)B(tài)分解的新概念,選擇合適的固有模態(tài)函數(shù)進(jìn)行后續(xù)包絡(luò)分析,利用EMD的帶通濾波特性,在IMF中捕獲待測(cè)結(jié)構(gòu)的共振頻帶,通過比較證實(shí)了該法的可行性.Cho等[5]針對(duì)機(jī)械故障特征提取等問題,采用振動(dòng)信號(hào),研究了基于EMD的信號(hào)降噪之法和基于SVM、最近鄰分類器的故障識(shí)別方法.可見利用EMD方法對(duì)信號(hào)進(jìn)行分解可以準(zhǔn)確有效地把握原始數(shù)據(jù)的特征信息,有利于信息特征挖掘[6-7].但EMD也有缺陷,即分解結(jié)果會(huì)出現(xiàn)模態(tài)混淆的現(xiàn)象.為此,Wu 等[8]提出新的方法,即集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)方法.不穩(wěn)定信號(hào)通過EEMD分解,可以得到平穩(wěn)的本征模態(tài)函數(shù),此時(shí)混淆度較低.Tabrizi等[9]提出了一種自動(dòng)檢測(cè)滾動(dòng)軸承微小缺陷的方法.首先,將小波包分解用于去除噪聲信號(hào);然后,利用EEMD技術(shù)提取信息特征向量;最后,利用支持向量機(jī)算法構(gòu)造超平面標(biāo)記樣本,對(duì)軸承狀態(tài)進(jìn)行檢測(cè).結(jié)果表明該方法能夠有效識(shí)別軸承的故障類別.針對(duì)齒輪振動(dòng)信號(hào)的非平穩(wěn)特征和現(xiàn)實(shí)中難以獲得大量典型故障樣本的實(shí)際情況,一種基于EEMD和支持向量機(jī)的齒輪故障診斷方法被提出,通過EEMD方法將非平穩(wěn)的原始加速度振動(dòng)信號(hào)分解成若干個(gè)平穩(wěn)的本征模態(tài)函數(shù),從包含有主要故障信息的IMF分量中提取出能量特征作為支持向量機(jī)的輸入,判斷齒輪的工作狀態(tài)和故障類型[10].

目前,熵理論開始應(yīng)用于機(jī)械故障診斷領(lǐng)域.部件發(fā)生不同損壞時(shí),不同頻帶內(nèi)的熵值也會(huì)變動(dòng).向丹等[11]研究了滾動(dòng)軸承故障診斷單一故障信號(hào)的局限性和故障特征的非線性,利用非線性動(dòng)力學(xué)參數(shù)熵作為特征,提出基于EMD熵特征融合的方法來解決滾動(dòng)軸承故障診斷問題.許凡等[12]針對(duì)滾動(dòng)軸承故障診斷中聚類模式識(shí)別方法需要預(yù)先設(shè)定聚類數(shù)目的問題,提出一種基于局部均值分解與基本尺度熵的相鄰傳播聚類故障診斷方法,并取得了良好效果.李兆飛等[13]提出基于過程功率譜熵與支持向量機(jī)理論的軸承故障分類方法,有望解決振動(dòng)故障隨機(jī)性和故障樣本提取困難的問題.可見,用EEMD分解后的熵判斷損壞程度是一種新思路.信號(hào)的起伏程度值用峭度描述:軸承正常狀況時(shí),其峭度值較小;隨著故障的發(fā)生,其峭度值會(huì)明顯增大.故IMF峭度值的大小也可作為故障特征的一種[14].特征提取是判別損傷情況的關(guān)鍵環(huán)節(jié),直接影響分類的準(zhǔn)確性.基于此,本文提出一種新穎的特征提取方法——將IMF的奇異熵、能量熵與峭度值融合,首先采用EEMD對(duì)原信號(hào)進(jìn)行處理,對(duì)分解后具有故障信息的信號(hào)求其奇異熵、能量熵與峭度值,再利用核主元分析(KPCA)進(jìn)行融合與降維[15],最后將融合數(shù)據(jù)輸入小樣本分類器SVM中進(jìn)行分類[16-17].實(shí)驗(yàn)結(jié)果表明此方法具有良好的故障診斷效果.

1 EEMD基本原理與熵特征

1.1 EEMD基本原理

EMD適用于不規(guī)律信號(hào)的處理,是Huang于1998年提出的方法[6,18-19],此方法將信號(hào)分解為若干本征模態(tài)函數(shù),以此來使瞬時(shí)頻率體現(xiàn)本身的物理意義.但分解時(shí)有模態(tài)混淆現(xiàn)象,即屬于前一個(gè)IMF的信號(hào)由于外部擾動(dòng)等影響,會(huì)將不同時(shí)間尺度的兩個(gè)IMF歸為一類,使信號(hào)頻域變形,影響判斷[20].為此提出EEMD來減小模態(tài)混疊的問題.EEMD方法中引入了頻率分布均勻的高斯白噪聲,使間斷信號(hào)變連續(xù),故能提高算法能力.算法步驟[20]為

(1) 原信號(hào)x(t),將標(biāo)準(zhǔn)差是常數(shù)、均值為0的白噪聲ni(t)引入,則有信號(hào)

xi(t)=x(t)+ni(t).

(1)

其中xi(t) 為第i次加入白噪聲后數(shù)據(jù).

(2) 對(duì)xi(t)進(jìn)行EMD分解,得到IMF,可知分量cij(t)和ri(t).cij(t)為在第i次加入白噪聲所得第j個(gè)IMF;ri(t)為xi(t)分解的余項(xiàng).

(3) 重復(fù)操作前兩步m次,每次添加隨機(jī)的高斯白噪聲序列,得到m個(gè)IMF,對(duì)m次分解得到的IMF進(jìn)行平均值求取,則EEMD的IMF為

(2)

其中cj(t)為分解后第j個(gè)IMF分量.則EEMD的分解結(jié)果為

(3)

其中r(t)為信號(hào)x(t)分解的余項(xiàng).

1.2 EEMD特征融合數(shù)據(jù)

1.2.1 峭度值

峭度是無量綱參數(shù),用以描述一組信號(hào)的波形[14].將原始信號(hào)經(jīng)過EEMD分解后求取前8個(gè)IMF的峭度值,則峭度的公式為

(4)

峭度值對(duì)沖擊信號(hào)十分敏感,適用于檢測(cè)軸承表面損傷情況.當(dāng)軸承正常工作時(shí),各IMF的峭度值在一定區(qū)間波動(dòng)且值較小;當(dāng)軸承出現(xiàn)損傷時(shí),沖擊變大,使含有故障信號(hào)的IMF分量峭度值不同程度地增大.因此,根據(jù)這一特征,可以將IMF分量的峭度值作為故障特征之一.

1.2.2 奇異熵

奇異熵可以用來描述數(shù)據(jù)成分的分布特征,因此可以用奇異熵來描述經(jīng)過EEMD分解信號(hào)的頻率組成成分特征[11],以此作為振動(dòng)信號(hào)特征數(shù)據(jù)之一.將EEMD分解后得到的IMF進(jìn)行奇異熵計(jì)算,那么EEMD奇異熵的表達(dá)式為

(5)

(6)

其中:σi為第i個(gè)奇異值;E為數(shù)據(jù)所有奇異值的和.

1.2.3 能量熵

當(dāng)滾動(dòng)軸承出現(xiàn)損傷時(shí),其振動(dòng)信號(hào)頻率分布也會(huì)發(fā)生改變,且其能量也會(huì)變化,為此引入能量熵作為損傷的特征數(shù)據(jù)之一[10].采用EEMD分解軸承信號(hào)得到IMF,假設(shè)IMF有M個(gè),并將本征模態(tài)函數(shù)對(duì)應(yīng)能量值求出為Ei(i=1,2,…,M),則有

(7)

(8)

2 KPCA-SVM算法原理

2.1 核主元分析(KPCA)算法原理

核主元分析的基本思想是將高維數(shù)據(jù)通過非線性映射降為低維空間的數(shù)據(jù),非線性映射函數(shù)Φ:Rm→F.其中F為映射的高維特征空間,可解釋為此高維空間是現(xiàn)有樣本在低維中的映射.利用KPCA提取數(shù)據(jù)特征[21],該方法巧妙地解決了樣本數(shù)據(jù)維數(shù)多、計(jì)算量大的特點(diǎn).

已知樣本數(shù)據(jù)集x1,x2,x3,…,xn∈Rm,利用Φ映射到高維F空間中,那么Φ(xi)的協(xié)方差矩陣C可求出為

(9)

對(duì)C進(jìn)行特征值分解,則特征值與特征向量為

λV=CV.

(10)

其中:λ為C的特征值;V為C的特征向量.

由此可定義核矩陣

K={kij}n×n.

(11)

其中kij=Φ(xi)·Φ(xj).

特征向量V在由Φ(xi)長成的空間內(nèi)有

(12)

將式(10)~式(12)代入式(9)中,求核矩陣K的特征值、特征向量.

nλα=Kα.

(13)

則樣本Φ(x)的投影方向?yàn)?/p>

(14)

80%

(15)

2.2 支持向量機(jī)(SVM)算法原理

圖1 最優(yōu)分類線Fig.1 Optimum separating line

2.3 特征提取與模式識(shí)別

將訓(xùn)練數(shù)據(jù)進(jìn)行EEMD分解,單一熵值在一定程度上有一定的局限性,很難完成多變信號(hào)的特征辨別.因此,采用奇異熵Ws、能量熵We與峭度值Kr相結(jié)合作為故障特征,融合的熵?cái)?shù)據(jù)能夠在一定程度上反映軸承的故障特征頻率.將得到EEMD分解后的IMF分別進(jìn)行奇異熵、能量熵與峭度值計(jì)算,將三組數(shù)據(jù)進(jìn)行組合得到特征數(shù)據(jù),將該數(shù)據(jù)輸入KPCA中進(jìn)行特征融合、特征提取與降維,再輸入到SVM中訓(xùn)練,最后將測(cè)試數(shù)據(jù)放入訓(xùn)練好的SVM中.SVM具有小樣本分類的特點(diǎn),可以進(jìn)行滾動(dòng)軸承的準(zhǔn)確分類,具體流程如圖2所示.

圖2 故障診斷流程Fig.2 Fault diagnosis flow chart

具體步驟如下:

(1) 在采樣頻率f下收集滾動(dòng)軸承正常、內(nèi)圈故障、滾動(dòng)體故障的振動(dòng)信號(hào)數(shù)據(jù),分別進(jìn)行N次采樣,故樣本一共獲得了3N個(gè)振動(dòng)信號(hào)數(shù)據(jù).

(3) 將F1、F2帶入KPCA中完成特征融合,將F1融合后數(shù)據(jù)存入SVM進(jìn)行訓(xùn)練,SVM1、SVM2組成多故障分類器,將F2融合后數(shù)據(jù)導(dǎo)入訓(xùn)練好的SVM中完成類別判斷.

3 試驗(yàn)與診斷步驟

為驗(yàn)證該方法的有效性,利用美國凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心的數(shù)據(jù)[23]進(jìn)行分析.試驗(yàn)裝置如圖3所示.采樣頻率為12 kHz,由1470 W三相感應(yīng)電動(dòng)機(jī)(左)、扭矩傳感器(中)和測(cè)功機(jī)(右)組成,轉(zhuǎn)速1772 r/min,利用傳感器收集數(shù)據(jù).用電火花加工的單點(diǎn)點(diǎn)蝕模擬軸承單點(diǎn)故障,選取3種軸承狀態(tài)模式:(1)正常狀態(tài),無點(diǎn)蝕,滾動(dòng)軸承正常運(yùn)行;(2)在滾動(dòng)軸承內(nèi)圈加單點(diǎn)點(diǎn)蝕故障,損壞點(diǎn)直徑為0.177 8 mm; (3)在滾動(dòng)體上加單點(diǎn)點(diǎn)蝕損傷,故障點(diǎn)直徑也為0.177 8 mm.

圖3 試驗(yàn)裝置Fig.3 Test bench

選取該實(shí)驗(yàn)平臺(tái)數(shù)據(jù)庫中的樣本:軸承正常狀態(tài)數(shù)據(jù)98.mat、內(nèi)圈故障數(shù)據(jù)106.mat、滾動(dòng)體故障數(shù)據(jù)119.mat(98.mat、106.mat和119.mat均為西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心的部分?jǐn)?shù)據(jù)名稱).從中各隨機(jī)抽取15組數(shù)據(jù),共45組數(shù)據(jù),每組有2048個(gè)數(shù)據(jù)點(diǎn),共92 160個(gè)數(shù)據(jù)點(diǎn).試驗(yàn)軸承數(shù)據(jù)集的描述如表1所示.

表1 軸承數(shù)據(jù)集描述Table 1 Description of bearing data set

滾動(dòng)軸承3種狀態(tài)的振動(dòng)信號(hào)時(shí)域圖如圖4所示.

圖4 滾動(dòng)軸承振動(dòng)信號(hào)時(shí)域圖Fig.4 Time domain diagram of rolling bearing vibration signal

將3種狀態(tài)的訓(xùn)練數(shù)據(jù)與測(cè)試數(shù)據(jù)分別進(jìn)行EEMD分解,得到各數(shù)據(jù)組的IMF分量,圖5為正常狀態(tài)下的一組IMF分量圖.圖5(a)為振動(dòng)信號(hào)及本征模態(tài)c1至c5的分解結(jié)果,圖5(b)為本征模態(tài)c6至c10以及余項(xiàng)r1的分解結(jié)果,可見主要信息在前幾個(gè)本征模態(tài)中,選取前8個(gè)IMF進(jìn)行數(shù)據(jù)特征求取.

圖5 正常狀態(tài)數(shù)據(jù)的EEMD分解結(jié)果Fig.5 EEMD decomposition results of normal state data

將得到的IMF進(jìn)行相應(yīng)的奇異熵、能量熵與峭度值計(jì)算,可以得到一個(gè)8×24的訓(xùn)練數(shù)據(jù)與7×24的測(cè)試數(shù)據(jù),由于數(shù)據(jù)篇幅較大,只列出每種運(yùn)行狀態(tài)的前兩組數(shù)據(jù)值,每組數(shù)據(jù)為經(jīng)EEMD分解后前3個(gè)IMF的奇異熵、能量熵、峭度值,如表2所示.

表2 部分IMF熵與峭度值數(shù)據(jù)Table 2 Partial IMFs entropy and kurtosis data

圖6為將正常數(shù)據(jù)、內(nèi)圈故障數(shù)據(jù)與滾動(dòng)體故障數(shù)據(jù)經(jīng)EEMD處理后的熵值及峭度值各拿出3組進(jìn)行歸一化后繪制的奇異熵、能量熵與峭度值分布圖,1~8個(gè)點(diǎn)為奇異熵,9~16個(gè)點(diǎn)為能量熵,16~24個(gè)點(diǎn)為峭度值.從圖6中可以看出3種軸承的數(shù)據(jù)分布具有一定的特征,但特征不夠明顯,因此利用KPCA進(jìn)行特征融合及提取特征值.

圖6 IMF奇異熵、能量熵與峭度值分布圖Fig.6 IMF singalar entropy,energy entropy and kurtosis vacue distribution

將軸承3種狀態(tài)分別提出5組IMF熵與峭度值,經(jīng)KPCA處理后的特征融合數(shù)據(jù)如圖7所示.從圖7中可以清楚地看出經(jīng)融合與特征提取后,每種狀態(tài)的特征更加明顯,更有利于故障分類.訓(xùn)練數(shù)據(jù)由原來的24×24維經(jīng)KPCA降維為24×4維,大大降低了維數(shù),并且有效地提取了特征值,訓(xùn)練數(shù)據(jù)經(jīng)KPCA處理后的部分特征值如表3所示.

圖7 經(jīng)KPCA處理后的特征融合數(shù)據(jù)Fig.7 Feature fasion data processed by KPCA

表3 KPCA特征值提取數(shù)據(jù)Table 3 KPCA feature extraction data

將KPCA降維后的訓(xùn)練數(shù)據(jù)輸入SVM進(jìn)行訓(xùn)練.本試驗(yàn)軸承共有3種故障類別:正常運(yùn)行(無故障)、內(nèi)圈故障與滾動(dòng)體故障.訓(xùn)練SVM1和SVM2兩個(gè)分類器:SVM1用于檢測(cè)軸承是否故障,正常運(yùn)行時(shí),SVM1輸出為+1,故障時(shí),SVM1輸出為-1;將SVM1輸出為-1的數(shù)據(jù)輸入訓(xùn)練好的SVM2,SVM2輸出為+1為軸承內(nèi)圈故障,輸出-1為軸承滾動(dòng)體故障.即先用訓(xùn)練好的SVM1檢測(cè)軸承是否故障,再將故障數(shù)據(jù)輸入訓(xùn)練好的SVM2判斷故障類型.數(shù)據(jù)篇幅較大,僅列出部分分類數(shù)據(jù),分類結(jié)果如表4、表5所示.

表4 SVM1檢測(cè)軸承是否故障Table 4 Detect bearing failure

表5 SVM2檢測(cè)軸承故障類型Table 5 Detect bearing failure type

判斷軸承是否故障,SVM1訓(xùn)練僅用了0.941 7 s.由表4可知:測(cè)試數(shù)據(jù)的軸承正常狀態(tài)數(shù)據(jù)經(jīng)SVM1后輸出都為+1,內(nèi)圈故障及滾動(dòng)體故障經(jīng)SVM1后輸出都為-1,可見分類正確率達(dá)100%.進(jìn)行故障類別分類,SVM2訓(xùn)練僅用0.779 2 s.由表5可知:將SVM1輸出為-1的數(shù)據(jù)輸入訓(xùn)練好的SVM2,若為內(nèi)圈故障,SVM2輸出為+1;若為滾動(dòng)體故障,SVM2輸出為-1.試驗(yàn)結(jié)果為:所有內(nèi)圈故障的測(cè)試數(shù)據(jù)經(jīng)SVM2后輸出均為+1,所有滾動(dòng)體故障的測(cè)試數(shù)據(jù)經(jīng)SVM2后輸出均為-1,可見分類正確率為100%.由表4和表5可知提出的方法故障診斷準(zhǔn)確率為100%,明顯高于僅有單一熵的原始特征故障診斷方法,可見此方法的可行性.

4 結(jié) 論

針對(duì)滾動(dòng)軸承振動(dòng)信號(hào)非線性非平穩(wěn)性的特點(diǎn),以及軸承故障特征提取對(duì)故障診斷起到?jīng)Q定性作用的問題,結(jié)合實(shí)際情況,采用EEMD對(duì)軸承振動(dòng)信號(hào)進(jìn)行處理,并提出將有故障特征的IMF進(jìn)行奇異熵、能量熵和峭度值相融合的特征提取方法,解決了以上問題,得出以下結(jié)論:

(1) 滾動(dòng)軸承采集的振動(dòng)信號(hào)為非線性非平穩(wěn)性信號(hào),因此采用EEMD的方法,將復(fù)雜信號(hào)分解為若干個(gè)平穩(wěn)的本征模態(tài)函數(shù),減少了模態(tài)混疊的問題.并求出前幾個(gè)本征模態(tài)函數(shù)的奇異熵、能量熵與峭度值,作為故障特征向量.

(2) 特征提取直接關(guān)系到分類的準(zhǔn)確性,因此,在故障分類前進(jìn)行特征提取極為重要.筆者提出用KPCA完成熵與峭度特征數(shù)據(jù)的融合,有效地提取故障數(shù)據(jù)特征,提高了分類的準(zhǔn)確性.并將KPCA處理后的數(shù)據(jù)輸入到SVM中,在本文應(yīng)用實(shí)例中得到了滿意的效果.

將滾動(dòng)軸承正常狀態(tài)、內(nèi)圈故障與滾動(dòng)體故障的振動(dòng)信號(hào)用于本研究進(jìn)行故障診斷,結(jié)果表明該方法能夠準(zhǔn)確無誤地對(duì)滾動(dòng)軸承狀態(tài)進(jìn)行分類,有望為實(shí)際應(yīng)用中滾動(dòng)軸承故障診斷提供一種新思路.

猜你喜歡
特征提取故障診斷模態(tài)
基于Gazebo仿真環(huán)境的ORB特征提取與比對(duì)的研究
電子制作(2019年15期)2019-08-27 01:12:00
一種基于LBP 特征提取和稀疏表示的肝病識(shí)別算法
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
因果圖定性分析法及其在故障診斷中的應(yīng)用
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
基于MED和循環(huán)域解調(diào)的多故障特征提取
由單個(gè)模態(tài)構(gòu)造對(duì)稱簡支梁的抗彎剛度
基于LCD和排列熵的滾動(dòng)軸承故障診斷
基于WPD-HHT的滾動(dòng)軸承故障診斷
高速泵的故障診斷
河南科技(2014年3期)2014-02-27 14:05:48
主站蜘蛛池模板: 看看一级毛片| 日本少妇又色又爽又高潮| 无码一区18禁| 亚洲欧美综合另类图片小说区| 亚洲国产亚洲综合在线尤物| 最新精品久久精品| 国产永久免费视频m3u8| 国产精品女主播| 日韩午夜片| 美女啪啪无遮挡| 国产精品久久久久久久久久久久| 国产va在线观看免费| 欧美69视频在线| 一级福利视频| 久久国产精品麻豆系列| 91精品国产福利| 欧美国产视频| 性视频一区| 欧美笫一页| 国产亚洲精| 91免费精品国偷自产在线在线| 久久久波多野结衣av一区二区| 国产免费a级片| 在线观看国产黄色| 国产精品一区二区不卡的视频| 欧美成人午夜视频| 国产成人乱码一区二区三区在线| 国产一在线| 伊人久久大香线蕉aⅴ色| 青青青国产视频手机| 国产浮力第一页永久地址| 亚洲品质国产精品无码| 精品一区二区无码av| 亚洲一本大道在线| 国产正在播放| 人人澡人人爽欧美一区| 国产一级做美女做受视频| 亚洲最新在线| 国产成人精品免费视频大全五级| 国产sm重味一区二区三区| 91麻豆国产在线| 91青青视频| 欧美天堂久久| 伊人久久久久久久| 在线观看无码a∨| av在线5g无码天天| 久久久受www免费人成| 一本一本大道香蕉久在线播放| 午夜日b视频| 亚洲第一成网站| 试看120秒男女啪啪免费| 91欧洲国产日韩在线人成| 久久人搡人人玩人妻精品一| 一区二区三区在线不卡免费| 国产精品白浆在线播放| 人妖无码第一页| 四虎国产永久在线观看| 国产波多野结衣中文在线播放| 亚洲综合国产一区二区三区| 久久成人18免费| 亚洲欧洲天堂色AV| 日韩无码黄色| 亚洲最大福利视频网| 国产女人在线| 国产美女91呻吟求| 免费一看一级毛片| 久久亚洲AⅤ无码精品午夜麻豆| 国产视频你懂得| 久久五月天国产自| 一区二区偷拍美女撒尿视频| 在线视频亚洲色图| www.国产福利| 日韩天堂视频| 永久免费av网站可以直接看的| 欧美日韩动态图| 亚洲国产精品不卡在线| 国产另类视频| 乱人伦99久久| 国产国语一级毛片| 亚洲一区网站| 久久综合九色综合97网| 久久99精品国产麻豆宅宅|