亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于階頻譜相關(guān)組合切片能量和SVM的滾動軸承故障診斷

        2017-12-23 02:58:46汪治安夏均忠白云川劉鯤鵬呂麒鵬
        軍事交通學(xué)院學(xué)報 2017年12期
        關(guān)鍵詞:階次外圈特征向量

        汪治安,夏均忠,白云川,劉鯤鵬,呂麒鵬

        (1.陸軍軍事交通學(xué)院 研究生管理大隊,天津 300161;2.陸軍軍事交通學(xué)院 軍用車輛系,天津 300161)

        基于階頻譜相關(guān)組合切片能量和SVM的滾動軸承故障診斷

        汪治安1,夏均忠2,白云川2,劉鯤鵬1,呂麒鵬1

        (1.陸軍軍事交通學(xué)院 研究生管理大隊,天津 300161;2.陸軍軍事交通學(xué)院 軍用車輛系,天津 300161)

        針對變轉(zhuǎn)速工況下滾動軸承故障診斷問題,提出基于階頻譜相關(guān)組合切片能量(SEOFSC)和支持向量機(SVM)的滾動軸承故障診斷方法。首先研究階頻譜相關(guān)(OFSC)對滾動軸承故障特征提取的原理;其次針對OFSC的不足,計算振動信號階頻譜相關(guān)組合切片以提高故障特征提取效率,對其譜頻率軸積分得到SEOFSC實現(xiàn)故障特征降維;最后選取切片簇中能量極大值構(gòu)成特征向量輸入到“一對一”多分類SVM中進行故障模式識別。通過對比試驗表明,該方法可以實現(xiàn)變轉(zhuǎn)速工況下滾動軸承快速、準確的故障診斷,具有一定的工程實用價值。

        滾動軸承;故障診斷;角度/時間循環(huán)平穩(wěn);階頻譜相關(guān)組合切片能量;支持向量機(SVM)

        v

        滾動軸承廣泛應(yīng)用于各類機械設(shè)備中,其工作狀態(tài)直接影響整臺機械設(shè)備的性能。在高速重載的環(huán)境下,滾動軸承容易失效[1]。對滾動軸承的故障診斷方法的研究一直是國內(nèi)外研究的熱點[2-3]。

        實際運行過程中,大多數(shù)機械通常工作在速度波動較大的變轉(zhuǎn)速工況下[4],滾動軸承故障振動信號出現(xiàn)復(fù)雜的頻率、幅值和相位調(diào)制等現(xiàn)象[5-6]。為此,Abboud等[7-8]在循環(huán)平穩(wěn)的基礎(chǔ)上,提出了角度/時間循環(huán)平穩(wěn)(angle/time cyclostationary, AT-CS),利用階頻譜相關(guān)(order-frequency spectral correlation, OFSC)表征變轉(zhuǎn)速下滾動軸承故障特征,為循環(huán)平穩(wěn)在變轉(zhuǎn)速工況下的滾動軸承故障診斷奠定了基礎(chǔ)。

        OFSC存在兩點不足[9-10]:①計算量大。為較好地從信號中提取滾動軸承故障特征,在估計信號的階頻譜相關(guān)時,需要保證較高的階次分辨率和足夠的階次范圍,從而使得遍歷的循環(huán)階次數(shù)目增多,計算量較大。②特征表達不直觀,維數(shù)較高。OFSC是通過三維圖的方式表征變轉(zhuǎn)速工況下滾動軸承故障特征,在噪聲干擾下,難以從圖中準確識別故障特征階次。

        Vapnik等[11]提出的支持向量機(support vector machine, SVM)算法為基于結(jié)構(gòu)風(fēng)險最小化方法的統(tǒng)計學(xué)習(xí)理論,具有很強的學(xué)習(xí)和泛化能力,能夠較好地解決小樣本、非線性、高維數(shù)等模式識別問題[12]。在滾動軸承的故障診斷中,由于采集的故障樣本數(shù)量有限,因此SVM被廣泛應(yīng)用于滾動軸承故障診斷。明陽等[13]將譜相關(guān)密度切片分析與SVM相結(jié)合實現(xiàn)滾動軸承故障診斷。張小龍等[14]提取ITD分解后PR分量的Lempelk-Ziv復(fù)雜度構(gòu)建特征向量,將其輸入到SVM中實現(xiàn)滾動軸承不同類型故障模式識別。

        本文在研究階頻譜相關(guān)的基礎(chǔ)上,提出基于階頻譜相關(guān)組合切片能量(sliced energy of order-frequency spectral correlation,SEOFSC)和支持向量機的滾動軸承故障診斷方法。通過計算信號的階頻譜相關(guān)組合切片,減少計算量,提高特征提取效率;對階頻譜相關(guān)組合切片的譜頻率軸積分得到SEOFSC,降低故障特征維數(shù);從SEOFSC中選擇局部能量極大值組成特征向量輸入到SVM中進行模式識別,實現(xiàn)滾動軸承故障診斷。

        1 特征提取

        1.1 階頻譜相關(guān)(OFSC)

        變轉(zhuǎn)速工況下,滾動軸承點蝕故障具有如下特點:故障產(chǎn)生的沖擊不再具有時域周期性,由于軸承是旋轉(zhuǎn)對稱結(jié)構(gòu),相鄰兩次沖擊間隔的角度不變,因此其故障沖擊仍具有角域周期性;同時,故障沖擊依然隨時間衰減,故其脈沖響應(yīng)在時域描述最佳。稱此類信號為角度/時間循環(huán)平穩(wěn)信號(AT-CS)。

        為實現(xiàn)變轉(zhuǎn)速工況下的滾動軸承故障特征提取,需要建立一個既能夠表征信號的角域周期性,又能夠描述時域脈沖響應(yīng)的統(tǒng)計量。將傳統(tǒng)循環(huán)統(tǒng)計量從時域轉(zhuǎn)化到角度/時間域,此時其角度/時間自相關(guān)函數(shù)(angle/time autocorrelation function,ATCF)具有周期性:

        R2x(τ,θ)=E{x(t(θ))·x(t(θ))-τ)*}=

        R2x(τ,t(θ))

        (1)

        式中:τ為時間延遲;E為集總平均運算。

        時間t和角位移θ存在下列關(guān)系:

        (2)

        式中ω(t)為瞬時角速度。

        定義ATCF的二次傅里葉變換為階頻譜相關(guān)(OFSC):

        (3)

        式中α為循環(huán)階,對于滾動軸承故障診斷,可視為參考軸每旋轉(zhuǎn)一周故障沖擊產(chǎn)生的次數(shù)。

        對于AT-CS信號,其OFSC具有如下性質(zhì):

        (4)

        階頻譜相關(guān)是譜頻率f和循環(huán)階α的函數(shù),在α—f平面內(nèi),存在一系列垂直于α軸的譜線,譜線對應(yīng)信號的循環(huán)特征階次及其諧波成分,譜線的幅值主要集中在信號的響應(yīng)頻帶上。根據(jù)階頻譜相關(guān)的這個特點,可以有效提取滾動軸承振動信號中的AT-CS成分,進行故障診斷。

        1.2 階頻譜相關(guān)組合切片能量

        為克服階頻譜相關(guān)在計算效率和特征表達方面的不足,分別從兩方面對其進行改進。

        為降低階頻譜相關(guān)的計算量,將組合切片思想與階頻譜相關(guān)估計相結(jié)合。根據(jù)滾動軸承的幾何參數(shù),選擇外圈、內(nèi)圈、滾動體的理論故障特征階次及其前3個倍頻,計算相應(yīng)的階頻譜相關(guān)組合切片,通過切片之間的能量對比能夠初步判斷軸承的故障狀態(tài),稱該方法為階頻譜相關(guān)組合切片分析。為降低故障特征的維數(shù),在滾動軸承的故障診斷中,只提取能夠反映調(diào)制特征的階次信息,而忽略載頻信息。對階頻譜相關(guān)組合切片方法進一步改進,將信號的階頻譜相關(guān)組合切片的譜頻率軸積分得到階頻譜相關(guān)組合切片能量,選取切片簇中能量極大值作為該特征循環(huán)頻率對應(yīng)的切片能量,將故障特征從三維降低到二維。

        基于階頻譜相關(guān)組合切片能量(SEOFSC)的故障特征提取流程如圖1所示,其主要步驟包括:

        (1)由滾動軸承的瞬時轉(zhuǎn)速曲線積分得到角位移曲線;

        (2)根據(jù)滾動軸承的結(jié)構(gòu)參數(shù)計算理論故障特征階次,并以理論故障特征階次及其前3個倍頻作為組合切片中心,選擇一定的切片寬度,并由此確定切片區(qū)間;

        (3)利用平均循環(huán)周期圖法計算振動信號在切片區(qū)間的OFSC,得到階頻譜相關(guān)組合切片;

        (4)選取合適的積分頻帶對階頻譜相關(guān)組合切片的譜頻率軸積分得到SEOFSC,降低故障特征的維數(shù);

        (5)對于每一個切片簇,選取其中能量極大值作為該切片中心對應(yīng)的能量值,將所有切片中心對應(yīng)的能量值構(gòu)成故障特征向量。

        圖1 軸承故障特征提取流程

        2 支持向量機(SVM)

        2.1 SVM分類原理

        對于線性可分的訓(xùn)練樣本集,Q={(xi,yi)|xi∈Rn,yi∈{-1,1},i=1,2,…,l},存在超平面H使得訓(xùn)練樣本中的兩類輸入分別位于H的兩側(cè)。

        H:wTx+b=0

        (5)

        式中:wT為分類面權(quán)向量;b為分類閾值。

        滿足式(5)條件的超平面并不唯一,需要尋找一個最優(yōu)超平面,使得所有樣本到超平面的距離最大,盡量減小分類誤差,實現(xiàn)結(jié)構(gòu)風(fēng)險最小化。

        (6)

        式(6)是一個帶線性不等式約束的二次規(guī)劃問題,為此,定義Lagrange函數(shù):

        (7)

        式中βi為Lagrange系數(shù)。

        f(x)=sgn {wTx+b}

        (8)

        以上分析是對于線性可分樣本,對于線性不可分樣本,利用核函數(shù)K(xi,yi)將樣本非線性映射到高維空間,在高維空間尋找最優(yōu)超平面。為允許一定程度的錯分,引入非負松弛變量ξi和懲罰因子C,將最優(yōu)分類面的約束條件變?yōu)?/p>

        (9)

        式中,松弛變量ξi指出了離群點的離群程度,值越大,離群越遠。懲罰因子C則決定了在多大程度上重視離群點帶來的危害,其值越大,重視程度越高。

        此時,最優(yōu)分類函數(shù)為

        (10)

        典型核函數(shù)有線性核函數(shù)、多項式核函數(shù)和高斯徑向基核函數(shù),3種核函數(shù)的表達形式如下:

        (11)

        2.2 多分類SVM

        SVM本質(zhì)上是兩分類的分類器,實際應(yīng)用中軸承故障類型多于兩類,因此需要構(gòu)造合適的多分類器。通過組合多個二分類SVM,可以實現(xiàn)多分類器構(gòu)造,常用方法有“一對多”“一對一”“二叉樹”等方法。其中“一對一”方法,是指訓(xùn)練時為任意兩類樣本構(gòu)造一個分類器,N個類別所需SVM分類器數(shù)量為N(N-1)/2。將待分類樣本輸入到每個兩類分類器中,統(tǒng)計每個分類器的分類結(jié)果,用“投票法”決定待分類樣本的所屬類別。這種方法具有訓(xùn)練和識別速度快的特點,但是隨著類別數(shù)的增加,需要構(gòu)造的二分類數(shù)目會急劇上升(類別數(shù)的平方級別),使得訓(xùn)練過程速度變慢。由于滾動軸承的故障類型有限,需要構(gòu)造的二分類數(shù)目較少,因此本文采用“一對一”的方法實現(xiàn)SVM多分類模式識別。

        3 軸承故障診斷

        基于OFSC組合切片能量與SVM的軸承故障診斷流程如圖2所示。

        圖2 軸承故障診斷流程

        采集不同故障類型的滾動軸承振動信號及其轉(zhuǎn)速信號,將其分為訓(xùn)練樣本和測試樣本。計算各個樣本的階頻譜相關(guān)組合切片能量,選取每個切片簇中能量極大值作為切片對應(yīng)的能量值,構(gòu)成特征向量。將訓(xùn)練樣本的特征向量及其對應(yīng)的故障類型輸入到多分類SVM中進行訓(xùn)練,得到SVM模型。將測試樣本的特征向量輸入到訓(xùn)練好的SVM模型中,根據(jù)SVM模型識別的結(jié)果與測試樣本對應(yīng)的故障類型進行比對,判斷SVM故障診斷的準確性。

        4 試驗驗證

        試驗裝置主要由三相驅(qū)動電機、變頻控制器、驅(qū)動軸、從動軸等組成[15]。通過變頻器控制電機轉(zhuǎn)頻在10~20 Hz之間,采樣頻率為50 kHz,采樣時長為21 s。

        試驗軸承安裝在從動軸上,其主要技術(shù)參數(shù)見表1。分別在3個軸承上加工外圈、內(nèi)圈、滾動體故障,故障截面形狀為正方形,深度為1 mm,邊長為3 mm。根據(jù)表1中的技術(shù)參數(shù)計算得到滾動軸承的理論故障特征階次,外圈故障特征階次BPO=3.592,內(nèi)圈故障特征階次BPI=5.409,滾動體故障特征階次BPB=0.399。

        表1 試驗軸承技術(shù)參數(shù)

        不同故障軸承試驗的轉(zhuǎn)速曲線和振動信號時域波形,分別如圖3、圖4所示。3種故障狀態(tài)的轉(zhuǎn)頻變化均在10~15 Hz,其對應(yīng)的轉(zhuǎn)速差最大可達300 r/min,說明轉(zhuǎn)速變化較快,且都經(jīng)歷了加速和減速過程,速度變換也較為復(fù)雜。3種故障振動信號的幅值受到轉(zhuǎn)速的調(diào)制,從中無法識別出故障類型。

        圖3 試驗軸承運行轉(zhuǎn)速曲線

        以外圈故障特征為例,介紹基于階頻譜相關(guān)組合切片能量的故障特征提取過程。對瞬時轉(zhuǎn)速曲線積分得到角位移曲線如圖5所示。

        圖5 外圈故障軸承角位移曲線

        根據(jù)滾動軸承理論故障特征階次及其二、三倍頻確定切片中心αT=[BPB, 2BPB, 3BPB,BPO,BPI, 2BPO, 3BPO,2BPI,3BPI](按升序排列),以Δα=0.05作為切片寬度,確定組合切片區(qū)間。計算組合切片區(qū)間對應(yīng)的階頻譜相關(guān),得到外圈階頻譜相關(guān)組合切片(如圖6所示)。

        通過觀察圖6中階頻譜相關(guān),發(fā)現(xiàn)其能量主要集中在10~500 Hz,故以此為積分頻帶對階頻譜相關(guān)組合切片的譜頻率軸積分得到階頻譜相關(guān)組合切片能量(如圖7所示)。

        圖6 外圈故障階頻譜相關(guān)組合切片

        圖7 外圈故障階頻譜相關(guān)組合切片能量

        從階頻譜相關(guān)組合切片能量中提取特征向量。分別提取各切片簇中能量極大值作為切片中心對應(yīng)的能量值,得到外圈故障信號的故障特征向量,同理,可得到內(nèi)圈和滾動體故障對應(yīng)的特征向量,3種故障對應(yīng)的特征向量見表2。對于某一具體故障類型的特征向量,其理論故障特征階次及其諧波處對應(yīng)的幅值較大,不同故障類型之間的特征向量差異性較大。

        表2 故障特征向量

        單個樣本的規(guī)律難以說明特征向量的穩(wěn)定性,為研究故障特征向量的一致性,提取不同故障多個樣本的特征向量進行了比較分析。分別提取外圈、內(nèi)圈和滾動體故障各30組樣本的特征向量,得到不同樣本特征向量的分布圖(如圖8所示)。同時計算同一故障的不同特征向量之間的平均相關(guān)系數(shù),外圈、內(nèi)圈、滾動體故障特征向量之間的平均相關(guān)系數(shù)分別為0.98、0.83、0.97,說明同種故障向量之間的故障特征向量高度相關(guān),結(jié)合圖8可以看出,對于同一種故障類型,其故障特征向量具有較好的魯棒性。

        圖8 不同故障特征向量分布

        試驗所用訓(xùn)練樣本從數(shù)據(jù)集中隨機抽取,外圈、內(nèi)圈、滾動體故障數(shù)據(jù)各選30組,每組100 000個數(shù)據(jù)點,共90組訓(xùn)練樣本,用于訓(xùn)練SVM。用同樣的方法得到90組測試樣本,用于測試SVM模型識別故障類型的準確性。

        用類別Ⅰ表示外圈故障,類別Ⅱ表示內(nèi)圈故障,類別Ⅲ表示滾動體故障。為測試SVM分類的準確性,分別用30、60、90組測試樣本對其進行驗證,其預(yù)測準確率見表3。說明通過多分類SVM可以準確識別滾動軸承的故障類型。

        表3 不同樣本數(shù)量識別結(jié)果

        5 結(jié) 論

        (1)針對階頻譜相關(guān)在計算效率和特征表達上的不足,分別從兩方面對其改進,提出基于階頻譜相關(guān)組合切片能量的特征提取方法,將組合切片思想融入到階頻譜相關(guān)估計,得到階頻譜相關(guān)組合切片,對其譜頻率軸積分得到相應(yīng)的階頻譜相關(guān)組合切片能量。較之于階頻譜相關(guān)方法,具有計算效率高,特征維數(shù)低的特點。

        (2)選取階頻譜相關(guān)組合切片簇中能量極大值構(gòu)成特征向量,該特征向量不受轉(zhuǎn)速變化的影響,具有較好的魯棒性。

        (3)利用支持向量機結(jié)構(gòu)風(fēng)險最小化的優(yōu)點,將階頻譜相關(guān)組合切片能量與“一對一”多分類SVM相結(jié)合實現(xiàn)變轉(zhuǎn)速工況下滾動軸承故障診斷,實測數(shù)據(jù)表明該方法能夠準確識別滾動軸承的外圈、內(nèi)圈、滾動體故障。

        [1] 夏均忠,趙磊,白云川,等.基于Teager能量算子和ZFFT的滾動軸承故障特征提取[J].振動與沖擊,2017,36(11):106-110.

        [2] 汪治安,夏均忠,但佳壁,等.循環(huán)平穩(wěn)在滾動軸承故障診斷中的應(yīng)用[J].軍事交通學(xué)院學(xué)報,2017,19(6):25-30.

        [3] RANDALL R B,ANTONI J. Rolling element bearing diagnostics: a tutorial[J].Mechanical Systems and Signal Processing, 2011,25(2): 485-520.

        [4] ABBOUD D, ANTONI J, SIEG-ZIEBA S, et al. Envelope analysis of rotating machine vibrations in variable speed conditions: a comprehensive treatment[J]. Mechanical Systems and Signal Processing, 2017, 84:200-226.

        [5] 林京,趙明.變轉(zhuǎn)速下的機械設(shè)備動態(tài)信號分析方法的回顧與展望[J].中國科學(xué)(技術(shù)科學(xué)),2015,45(7):669-686.

        [6] ANTONI J, GRIFFATON J, ANDRE H, et al. Feedback on the surveillance 8 challenge: vibration-based diagnosis of a safran aircraft engine[J]. Mechanical Systems and Signal Processing, 2017,97:1112-144.

        [7] ABBOUD D, ANTONI J, ELTABACH M, et al. Angle/time cyclostationarity for the analysis of rolling element bearing vibrations[J]. Measurement, 2015,75:29-39.

        [8] SOPHIE B, DIDER R, ANTONI J, et al. Non-intrusive rattle noise detection in non-stationary conditions by an angle/time cyclostationary approach[J]. Journal of Sound and Vibration, 2016, 366:501-513.

        [9] ANTONI J,XIN G,HAMZAOUI N. Fast computation of the spectral correlation[J]. Mechanical Systems and Signal Processing, 2017, 92: 248-277.

        [10] 畢果,陳進.組合切片分析在滾動軸承故障診斷中的應(yīng)用研究[J].機械科學(xué)與技術(shù),2009,28(2):182-185.

        [11] VAPNIL V N.The nature of statistical learning theory[M].New York:Springer-Verlag,1995:92.

        [12] LIU Y, BI J W, FAN Z P. A method for multi-class sentiment classification based on an improved one-vs-one (OVO) strategy and the support vector machine (SVM) algorithm[J]. Information Sciences, 2017, 394-395:38-52.

        [13] 明陽,陳進.基于譜相關(guān)密度切片分析和SVM的滾動軸承故障診斷[J].振動與沖擊, 2010,29(1):196-199.

        [14] 張小龍,張氫,秦仙蓉,等.基于ITD復(fù)雜度和PSO-SVM的滾動軸承故障診斷[J].振動與沖擊, 2016,35(24):102-107.

        [15] MISHRA C,SAMANTARAY A K,CHAKRABORTY G,et al. Rolling element bearing defect diagnosis under variable speed operation through angle synchronous averaging of wavelet de-noised estimate[J]. Mechanical Systems and Signal Processing, 2016,72-73:206-222.

        RollingBearingFaultDiagnosisBasedonSEOFSCandSVM

        WANG Zhian1, XIA Junzhong2, BAI Yunchuan2, LIU Kunpeng1, LYU Qipeng1

        (1.Postgraduate Training Brigade, Army Military Transportation University, Tianjin 300161, China;2.Military Vehicle Department, Army Military Transportation University, Tianjin 300161, China)

        Considering the fault diagnosis problem of rolling bearing under variable speed condition, the paper proposes a fault diagnosis method of rolling bearing based on sliced energy of order-frequency spectral correlation (SEOFSC) and support vector machine (SVM). Firstly, it studies the principle of order-frequency spectral correlation (OFSC) on bearing fault feature extraction. Then, it calculates combined slices of OFSC to improve the efficiency of fault feature extraction, and obtains SEOFSC by integrating the spectral frequency axis to reduce the fault feature dimension. Finally, it selects maximal values from slice cluster to compose feature vector and input it into “one-to-one” multi-classification SVM for fault pattern recognition. The experimental indicate that the proposed method can realize rapid and accurate fault diagnosis of the rolling bearing under the variable speed condition, and it has some practical value in engineering.

        rolling bearing; fault diagnosis; angle/time cyclostationary; sliced energy of order-frequency spectral correlation(SEOFSC); support vector machine(SVM)

        2017-08-30;

        2017-09-12.

        汪治安(1992—),男,碩士研究生;

        夏均忠(1967—),男,博士,教授,碩士研究生導(dǎo)師.

        10.16807/j.cnki.12-1372/e.2017.12.007

        TH133.33;TP206+.3

        A

        1674-2192(2017)12- 0029- 06

        (編輯:張峰)

        猜你喜歡
        階次外圈特征向量
        二年制職教本科線性代數(shù)課程的幾何化教學(xué)設(shè)計——以特征值和特征向量為例
        深溝球軸承外圈表面凹坑缺陷分析
        哈爾濱軸承(2022年1期)2022-05-23 13:13:16
        克羅內(nèi)克積的特征向量
        角接觸球軸承外圈鎖口高度自動檢測規(guī)改進
        哈爾濱軸承(2020年2期)2020-11-06 09:22:34
        階次分析在驅(qū)動橋異響中的應(yīng)用
        基于Vold-Kalman濾波的階次分析系統(tǒng)設(shè)計與實現(xiàn)*
        一類特殊矩陣特征向量的求法
        基于齒輪階次密度優(yōu)化的變速器降噪研究
        價值工程(2017年28期)2018-01-23 20:48:29
        EXCEL表格計算判斷矩陣近似特征向量在AHP法檢驗上的應(yīng)用
        雙溝球軸承外圈冷輾擴數(shù)值模擬與試驗研究
        軸承(2010年2期)2010-07-28 02:25:56
        亚洲视频一区二区久久久| 黑人巨大跨种族video| 国产精品污www一区二区三区| 亚洲 日韩 在线精品| 亚洲天堂一二三四区在线| 老女老肥熟女一区二区| 2021久久精品国产99国产精品| 日日摸夜夜欧美一区二区| av有码在线一区二区| 亚洲天堂成人av在线观看| 伊人激情av一区二区三区| 欧美韩日亚洲影视在线视频| 日韩极品免费在线观看| 日韩女优av一区二区| 白天躁晚上躁麻豆视频| 久久综合视频网站| 中文字幕东京热一区二区人妻少妇| 秋霞在线视频| 丰满女人又爽又紧又丰满| 无码一区二区三区在| av免费网站免费久久网| 无码人妻丰满熟妇区五十路| 亚洲男同志gay 片可播放| 看黄色亚洲看黄色亚洲 | 91丝袜美腿亚洲一区二区| 亚洲日韩av无码中文字幕美国 | 97人妻无码免费专区| 24小时在线免费av| 欧美俄罗斯40老熟妇| 加勒比无码专区中文字幕| 少妇太爽高潮在线播放| 国产成人精品亚洲日本在线观看 | 大陆啪啪福利视频| 亚洲乱码中文字幕视频| 高清不卡一区二区三区| 91人妻无码成人精品一区91| 久久精品国产亚洲av热东京热| 欧美 日韩 人妻 高清 中文| 亚洲a∨无码一区二区| 国产区高清在线一区二区三区| 变态另类人妖一区二区三区 |