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

基于多尺度子帶樣本熵和LPP的軸承故障診斷方法

2016-11-24 06:36:26王廣斌杜謀軍韓清凱李學軍
振動與沖擊 2016年20期
關鍵詞:故障診斷故障信號

王廣斌, 杜謀軍, 韓清凱, 李學軍

(湖南科技大學 機械設備健康維護湖南省重點實驗室,湖南 湘潭 411201)

?

基于多尺度子帶樣本熵和LPP的軸承故障診斷方法

王廣斌, 杜謀軍, 韓清凱, 李學軍

(湖南科技大學 機械設備健康維護湖南省重點實驗室,湖南 湘潭 411201)

軸承損傷是機械設備損傷的主要原因之一,其產生的振動信號具有微弱、非平穩和非線性的特點。針對不能準確從微弱信號中提取故障特征的問題,提出使用多尺度子帶樣本熵,首先對信號進行小波包分解得到多尺度信號,再將每一個多尺度信號進行子帶分解得到多尺度子帶信號,再求其樣本熵得到多尺度子帶樣本熵,該方法能深入挖掘微弱信號的本質特征;針對非平穩信號能量密度分布不均的問題,提出使用平滑偽Wigner-Ville分布,其可對非平穩信號的瞬時對稱相關函數進行時頻聚集處理,使信號的能量均勻分布;針對不能準確的挖掘非線性數據的主流形的問題,提出使用局部保持投影(LPP,Locality Preserving Projection),LPP在投影過程中保持了最優的數據局部鄰域關系,可以準確的挖掘非線性數據的主流形。文中分別采用四組正常、內圈故障、滾珠故障和外圈故障信號作為原始數據來驗證該方法的有效性,實驗結果證明該方法能有效地對信號故障進行分離和識別。

軸承損傷;特征提??;多尺度子帶樣本熵;平滑偽Wigner-Ville分布;LPP分布

隨著制造業的發展,機械水平的不斷提高,軸承的使用也越來越廣泛,但軸承發生故障的頻率也在不斷增加,由于軸承故障不易察覺,從而容易引起不必要的損失。軸承故障診斷的發展,一直以來都受到了國內外的廣泛關注,能否準確的從軸承故障信號中提取出準確的特征參數,是軸承故障診斷的關鍵。隨著人們對故障診斷方法的研究,越來越多的方法被提了出來,以2000年在Science雜志上發表的兩篇等距映射(Isometric Mapping, ISOMAP)和局部線性嵌入(Locally Linear Embedding, LLE)算法論文[1-2]為起點興起的流形學習算法研究,使得故障診斷方法獲得了新的動力。經典的流形學習方法有等距特征映射(Isometric Mapping, ISOMAP)[1]、拉普拉斯特征映射(Laplacian Eigenmaps, LE)[2]、局部線性嵌入(Locally Linear Embedding, LLE)[3]、局部切空間排列(Local Tangent Space Alignment, LTSA)[4]等。丁曉喜等[5]利用小波包分解和LPP相結合的方法,有效的提取了軸承故障信號潛在的特征信息,并對軸承故障及故障不同的損傷程度進行了診斷;張曉濤[6]利用多尺度正交PCA與LPP相結合的方法,消除投影分量間的冗余信息,是處理之后的齒輪箱故障信號內含的故障特征得到增強,提高了故障的識別率;李國芳[7]利用2DPCA與LPP相結合的方法對人臉特征進行了提?。秽嵔碌萚8]利用多尺度熵和支持向量機相結合的方法,對軸承故障特征進行提取,并有效的對軸承進行故障診斷;臧懷剛等[9]利用EMD和平滑偽Wigner-Ville相結合的方式,充分提取信號的平滑偽Wigner-Ville譜熵,對軸承進行故障診斷。

本文提出一種基于多尺度子帶樣本熵和LPP相結合的軸承故障診斷方法,該方法利用了多尺度子帶樣本熵對微弱信號的時頻域的挖掘能力,利用了平滑偽Wigner-Ville分布不僅可以使非平穩信號的能量密度均勻化,還可以消除在子帶分解時產生的模態混疊效應,再結合LPP對非線性數據主流形的挖掘能力,最后經過實驗對比,驗證該方法的有效性。

1 多尺度子帶樣本熵

樣本熵是時間序列復雜度的一種度量,其是為了彌補近似熵匹配自身的缺陷,且具有比近似熵更高的精度[10],但是如果數據過于復雜,則普通的樣本熵提取則不能深入數據內部對其進行準確的提取。為了解決該問題在這里提出了多尺度子帶樣本熵這個概念,首先將信號用小波包分解成多尺度信號,初步對信號特征值進行提取,然后選取合適的子帶因子對信號進行子帶分解得到多尺度子帶,最后求出每一個多尺度子帶信號的樣本熵,因為該方法能深入數據內部,并對數據的特征值和特征向量進行層層分解和挖掘,所以能有效的將復雜微弱數據的故障特征準確的提取出來。多尺度子帶樣本熵具體步驟為:

(1) 原始信號為X,X=[x1,x2,x3,…,xn]T;

(2) 對信號X=[x1,x2,x3,…,xn]T進行小波包分解,即:

(1)

(2)

式中:hk為正交高通濾波器系數;gk為正交低通濾波器系數;Y0(t)=φ(t)為尺度函數;Y1(t)=φ(t)為小波函數,由此可知函數{Yn(t)|n∈z+}稱為正交尺度函數φ(t)的正交小波包;

(3) 信號X(t)經過小波包j層分解與重構后,可以得到2j個小波包分解與重構序列S(j,k)(k=0,1,2,…,2j-1),S(j,k)為小波包對信號X(t)進行j層分解的第k個節點序列,這種小波包分解與重構可以看作是小波包對信號的一種劃分,則定義這種劃分的測度為:

(3)

式中:φF(j,k)(i)是φF(j,k)的第i個值(i=1,2,3,…,N),N是原始信號的長度,在本文中j=3,則k=0,1,2,…,7,且φF(j,k)是φ(j,k)的傅里葉變換;

(4) 依據數據經小波包分解后節點序列S(j,k)的數據長度M,選取樣本熵最佳子帶因子a(a=1,2,3,…);

(5) 將數據S(j,k)按照最佳子帶因子分為a個數據子帶Si(j,k)(i=1,2,…,a);

(6) 以S1(j,k)為例,給定維數m,將子帶信號S1(j,k)按照序號組成一組m維向量,即:

S1(j,k)(i)=[S1(j,k)(i),…,

S1(j,k)(i+m+1)]

(4)

(7) 計算S1(j,k)(i)與S1(j,k)(n)之間的距離din,則din為:

din=d[S1(j,k)(n)-S1(j,k)(i)]=

max|S1(j,k)(n+l)-S1(j,k)(i+l)|

(5)

式中:l=0,1,2,…,m-1。

(6)

(9) 增加維數至m+1,重復步驟(6)~(8),得到Bm+1(r);

(10) 若D為有限值,計S1(j,k)的樣本熵為c1,則其樣本熵c1為:

(7)

(11) 重復步驟(6)~(10),分別計算出S1(j,k),S2(j,k),S3(j,k),…,Sa(j,k)的樣本熵c1,c2,c3,…,c(a);

(12) 將樣本熵c1,c2,c3,…,c(a)進行特征矩陣構造,并進行歸一化處理,記作C(j,k),則該節點序列的子帶樣本熵為C(j,k);

(13) 在本文中,分別計算各個節點序列S(j,k)(j=3),則k=0,1,2,…,7,記其樣本熵為C(3,0),C(3,1),C(3,2),C(3,3),C(3,4),C(3,5),C(3,6)和C(3,7),并將其進行特征矩陣構造,進行歸一化處理,則可得該故障信號的多尺度子帶樣本熵為T,則:

T=[C(3,0),C(3,1),C(3,2),C(3,3),

C(3,4),C(3,5),C(3,6),C(3,7)]

(8)

從以上步驟可以看出,樣本熵的值與m和r密切相關,因此,對m和r的值的確定非常重要。其中m=1或m=2,0.1Std≤r≤0.25Std(Std是原始數據的標準差),本文取m=2和r=0.2Std。

2 平滑偽Wigner-Ville分布

在故障信號的分析過程中,對信號進行時頻分析是一個很重要的環節,時頻分析能直觀的了解到信號在時頻域上的分布,Wigner-Ville分布就是一種比較好的時頻分析方法,它能有效的分解出信號的時頻特征?,F有信號為s(t),則其Wigner-Ville分布如下所示[11-12]:

(9)

從上式可以看出,Wigner-Ville可以看成是信號的瞬時相關函數的傅里葉變換,其還具有信號在時間-頻率平面上的能量密度分布的含義,所以其結果能反應信號的時頻特性。Wigner-Ville擁有許多優良特性,比如其具有時移不變性、頻移不變性、時域有界性和頻域有界性等,使其在許多領域得到成功應用,雖然其有很多有用特性,但它也存在交叉干擾項等缺陷,在式(1)的定義中,函數x(t)出現了兩次,說明這是一種雙線性變換,不滿足疊加原理。

為了解決Wigner-Ville存在的交叉項的問題,提出了平滑偽Wigner-Ville分布,定義兩個平滑窗函數為h(τ)和g(u),利用平滑窗函數可以抑制Wigner-Ville帶來的交叉項效應,利用這兩個平滑窗函數改進Wigner-Ville分布可定義為:

SPWVDx(t,f)=

(10)

3 局部保持投影(LPP)

局部保持投影(LPP)是經典流形學習方法的一種,它是一種線性降維算法[14-15],是以拉普拉斯特征投影的線性近似為理論基礎的,可以有效的保留數據內在的非線性結構和數據子空間的局部特性,因而得到廣泛應用。

假設有一個d維的高維矩陣X=[x1,x2,x3,…,xn],其通過一個d維的映射矩陣W后得到一組向量矩陣Y=[y1,y2,y3,…,yn],其中,Y=WTX。

LPP的優化函數為:

(11)

式中:S為相似性矩陣,用k近鄰點法定義為:

近鄰點k的選取直接關系到Sij的大小,若xi和xj相互都在對方的k近鄰域內,則可確定Sij的大小,其中t>0。

對優化函數進行代數變換得:

WTX(D-S)XTW=WTXLXTW

(13)

式中:X=[x1,x2,x3,…,xn],D是n×n的對角矩陣,對角元素Dii=∑jSij,拉普拉斯矩陣L=D-S,存在約束條件:WTXLXTW=1。則W可通過求解下式所示的廣義特征值而得到:

XLXTW=μXDXTW

(14)

則可知,映射函數W是式(14)的前m個特征值μ1,μ2,μ3,…,μm(μm-1<μm)所對應的特征向量w1,w2,w3,…,wm構成,則W=[w1,w2,w3,…,wm]。

4 基于多尺度子帶樣本熵和LPP的故障特征提取

針對軸承故障信號具有非平穩、非線性和微弱的特點,本文提出了多尺度子帶樣本熵、平滑偽Wigner-Ville分布和LPP為一體的故障診斷方法,因為多尺度子帶樣本熵可以很好的挖掘出微弱故障信號的本質特征,平滑偽Wigner-Ville分布可使信號的能量密度均勻分布,局部保持投影(LPP)可以準確的挖掘非線性數據的主流形,所以該方法能準確的對軸承故障進行診斷與識別。該方法的具體步驟如下:

(1) 分別選取四種故障在相同故障尺寸不同負載下的四組故障信號并構成信號矩陣;

(2) 選取合適的小波基函數和分解層數,求出故障信號的多尺度子帶樣本熵;

(3) 對故障信號的多尺度子帶樣本熵進行希爾伯特變換,求出其解析信號;

(4) 選擇合適的窗函數及其長度,求出解析信號的平滑偽Wigner-Ville分布,挖掘故障信號的時頻分布特征;

(5) 對經過平滑偽平滑偽Wigner-Ville分布之后的信號進行LPP降維,挖掘故障信號的主流形結構,并得到故障信號的特征值和特征向量;

(6) 根據故障信號的特征值和特征向量,對軸承故障進行識別。

該故障診斷方法的流程圖如圖1所示。

圖1 故障診斷方法流程圖Fig.1 Fault diagnosis flowchart

5 實驗驗證分析

5.1 數據來源

本文針對滾動軸承各種常見的故障的數據用該方法進行故障診斷,并用其他方法加以對比,從而驗證該方法的有效性。數據來自于美國Case Western Reserve University軸承數據中心提供的公開數據,軸承型號6205-2RS JEM SKF深溝球軸承,分別選取正常、內圈故障、外圈故障和滾動體故障四種類型的數據,每一種故障分別選取四組負載和轉速不一樣的數據,軸承故障尺寸為0.021″,電機負載分別為0 HP 、1 HP、2 HP和3 HP,轉速分別為1 797 r/min、1 772 r/min、1 750 r/min和1 730 r/min, 信號的采樣頻率為12 000 Hz。

5.2 故障特征提取

軸承故障具有三種形式,分別為內圈故障、外圈故障和滾珠故障,每一種故障都具有不同形式的故障尺寸和故障頻率等,以滾珠故障為例,滾珠故障信號的時頻域圖(見圖2)。

圖2 滾珠故障信號的時域頻譜圖Fig.2 Time domain diagram and spectrum of ball fault signal

求取每一個故障信號的多尺度子帶樣本熵,初步提取故障信號的故障特征。在本文中,小波基函數為db3小波基函數,分解層數為3層,可求得8個尺度信號;再取每個子帶信號的最佳子帶因子,在本文中,將每一個尺度信號分成4個子帶信號,再求取每一個子帶的樣本熵,將求得的多尺度子帶樣本熵進行矩陣構造,將其作為特征矩陣向量。表1為不同故障信號的多尺度子帶樣本熵,從表1可以看出,四種狀態的多尺度子帶樣本熵各有不同,正常情況下的多尺度子帶樣本熵的平均值最大,外圈故障和內圈故障次之,滾珠故障的多尺度子帶樣本熵平均值最小,在每個故障的多尺度子帶樣本熵集合中,其大小都不一樣,表明在數據內部各子帶的結構、復雜度和所包含的信息等不一樣,這些異同都可以由多尺度子帶樣本熵表現出來。

表1 四種故障的多尺度子帶樣本熵

圖3為四種故障的多尺度子帶樣本熵的走勢圖,從圖中可以更直觀的看出同種故障信號在不同尺度不同頻段的多尺度子帶樣本熵不一樣,同一狀態下的多尺度子帶樣本熵變化平穩,且隨著多尺度子帶序列的推移慢慢減小,且各種狀態的多尺度子帶樣本熵不相同,可以表明同種故障信號內部的復雜性。因此可以看出多尺度子帶樣本熵不僅可以表現出微弱故障信號的本質特征,又可以將信號的變化特征完整的提取出來,所以采用多尺度子帶樣本熵可以對不同故障進行區分。因為軸承在實際運轉過程中,其工作環境復雜,背景噪聲較大,所以就導致采集的軸承故障數據復雜,且具有非線性、微弱和不平穩的特性,所以,為了更好的的表征其最準確的特征值和特征向量,挖掘其準確的主流形,我們就可以采用多尺度子帶樣本熵。

圖3 四種不同故障的多尺度子帶樣本熵Fig.3 The multi-scale sample entropy of four different fault

圖4 故障信號經過LPP分解的識別效果圖Fig.4 Identification effect of fault signal after LPP decomposition

圖5 故障信號經過多尺度子帶樣本熵和平滑偽Wigner-Ville分布的識別效果圖Fig.5 Recognition renderings of multi-scale sub-band sample entropy and smoothed pseudo Wigner-Ville distribution

圖6 故障信號經平滑偽Wigner-Ville分布和LPP分解的識別效果圖Fig.6 Recognition renderings of smoothed pseudo Wigner-Ville distribution and LPP decomposition

圖7 故障信號經過多尺度子帶樣本熵和LPP分解的識別效果圖Fig.7 Recognition renderings of multi-scale sub-bandsample entropy and LPP decomposition

對比圖4~圖7可知,從圖4可以看出,內圈故障可以和其他三種故障分開,但是正常、滾珠故障和外圈故障不能完全分開,故障相互融合, 無法進行故障識別;從圖5可以看出,正常、內圈故障、滾珠故障和外圈故障都無法分開,相互重疊,無法進行故障識別;從圖6中可以看出,正常故障可以識別,但內圈故障、滾珠故障和外圈故障相互重疊,無法進行故障識別;從圖7可以看出,四種故障內聚性較好,且故障沒有重疊的地方,可以很好的進行故障識別,能達到故障識別的要求,說明基于多尺度子帶樣本熵和LPP分布的故障診斷方法可以有效的對軸承故障進行故障診斷。圖8~圖10為故障識別效果的不同平面的識別圖。

圖8 為XY平面的故障效果識別圖Fig.8 Recognition renderings of XY plane

圖9 為XZ平面的故障效果識別圖Fig.9 Recognition renderings of XZ plane

圖10 為YZ平面的故障效果識別圖Fig.10 Recognition renderings of YZ plane

從圖8~圖10中可以更好地看出,正常、內圈故障、滾珠故障與外圈故障在不同的平面上都沒有相互重疊的部分,各種故障相互分開,且內聚性好,識別效果好。

6 結 論

本文基于多尺度子帶樣本熵、平滑偽Wigner-Ville分布和LPP分布的基礎上,結合多尺度子帶樣本熵對微弱信號時頻域的挖掘能力、平滑偽Wigner-Ville分布對非平穩信號的能量密度均分的能力和LPP對非線性信號特征提取的優越性,對滾動軸承故障進行診斷。結果表明基于多尺度子帶樣本熵和LPP分布的軸承故障診斷方法對軸承的故障診斷有很好的效果,該方法能夠深入挖掘故障數據的主流行結構,從而可以更好地獲得其故障特征。

[1] SEUNG H S, DANIEL D L.The manifold ways of perception[J].Science(S0036-8075),2000,290(5500):2268-2269.

[2] BELKIN M, NIYOGI P. Laplacian eigenmaps for dimensionality reduction and data representation[J]. Neural Computation,2003,15(6):1373-1396.[3] ROWEIS S, SAUL L. Nonlinear dimensionality reduction by locally linear embedding[J]. Science(S0036-8075),2000,290(5500):2323-2326.

[4] ZHANG Zhenyue, ZHA Hongyuan. Principal manioflds and nonlinear dimension reduction via tangent space alignmnet[J]. SIAM Journal on Scientific Computing,2005,26(1):313-338.

[5] 丁曉喜,何清波.基于WPD和LPP的設備故障診斷方法研究[J].振動與沖擊,2014,33(3):89-93.

DING Xiaoxi, HE Qingbo. Machine fault diagnosis based on WPD and LPP[J]. Journal of Vibration and Shock, 2014,33(3):89-93.

[6] 張曉濤.基于多尺度正交PCA-LPP流形學習算法的故障特征增強方法[J].振動與沖擊,2015,34(13):66-70.

ZHANG Xiaotao. Fault feature enhancement method based on multi-scale orthogonal PCA-LPP manifold learning algorithm[J].Journal of Vibration and Shock, 2015,34(13):66-70.

[7] 李國芳.基于2DPCA和流形學習LPP算法的人臉特征提取應用[J].電腦知識與技術,2014,10(31):7438-7441.

LI Guofang. Face feature extraction based on 2DPCA and LPP manifold learning algorithm[J].Computer Knowledge and Technology,2014,10(31):7438-7441.

[8] 鄭近德,程軍圣,胡思宇.多尺度熵在轉子故障診斷中的應用[J].振動、測試與診斷,2013,33(2):294-297.

ZHENG Jinde, CHENG Junsheng, HU Siyu. Rotor fault diagnosis based on multi-scale entropy[J]. Journal of Vibration, Measurement & Diagnosis,2013,33(2):294-297.

[9] 臧懷剛,王石云,李玉奎. EMD和平滑偽Wigner-Ville譜熵的軸承故障診斷[J]. 噪聲與振動控制,2014,34(5):145-149.

ZANG Huaigang, WANG Shiyun, LI Yukui. Bearing fault diagnosis based on EMD and smoothed pseudo Wigner-Ville spectrum entropy[J].Noise and Vibration Control,2014,34(5):145-149.

[10] RICHMAN J S, MOORMAN J R. Physiological time series analysis using approximate entropy and sample entropy[J]. Am J Physiol Circ Heart Physio,2000,278:2039-2049.

[11] 齊曉軒,郭婷婷,賈志勇. 基于Fast-ICA的Wigner-Ville分布交叉項消除方法[J].計算機工程,2015,41(8):71-75.

QI Xiaoxuan, GUO Tingting, JIA Zhiyong. Approach of eliminating cross-term of Wigner-Ville distribution based on Fast-ICA[J]. Computer Engineering,2015,41(8):71-75.

[12] GUPTA S, DEGRANDE G, LOMBAERT G. Experimental validation of a numerical model for subway induced vibrations[J]. Journal of Sound and Vibration, 2009,321(3/4/5):786-812.

[13] 樂葉青. 基于Wigner-Ville分布的電能質量擾動的分析[D].浙江:浙江大學,2007:11-21.

[14] LI Wei, PRASAD S, FOWLER J E, et al. Locality-preserving dimensionality reduction and classification for hyperspectral image analysis[J]. IEEE Transactions on Geoscience and Remote Sensing,2012,50(4):1185-1198.

[15] LU Guifu, LIN Zhong, JIN Zhong. Orthogonal complete discriminant locality preserving projections for face recognition[J].Neural Processing Letters,2011,33(3):235-250.

A bearing fault diagnosis method based on multi-scale sub-band sample entropy and LPP

WANG Guangbin, DU Moujun, HAN Qingkai, LI Xuejun

(Hunan Provincial Key Laboratory of Health Maintenance for Mechanical Equipment, Hunan University of Science and Technology, Xiangtan 411201, China)

One of the primary causes of machinery malfunction is bearing damage. Its vibration signal is weak, non-stationary, and nonlinear. This work proposed a multi-scale sub-band sample entropy concept aiming at the problem of eigenvalues and eigenvectors that cannot be accurately extracted from the weak signals. Firstly, the multi-scale signal was obtained by wavelet packet decomposition. Then, sub-band decomposing of the scales was realized. Finally, the sample entropy of each sub-band was solved. This method can dig deep into the essential characteristics of the weak signal. It presents that adopting smooth pseudo Wigner-Ville to solve the problem of uneven distribution of energy density of non-stationary signal. It can be used to deal with the time and frequency aggregation of the instantaneous symmetric correlation function of non-stationary signals, and make the signal energy distributed evenly. It proposed that using LPP (Locality Preserving Projection) decomposition to settle the issue of getting accurate image of the mainstream of the nonlinear data. LPP in the process of projection keeps the relationship between the optimal data of the local neighborhood. The main manifold can be accurately excavated from nonlinear data. Signals from a group of normal, inner ring fault, balls fault and outer ring fault bearings were used to verify the method. Experimental results prove that the method can effectively separate and identify bearing failure.

bearing damage; feature extraction; multi-scale sub-band sample entropy; smoothed pseudo Wigner-Ville distribution; LPP decomposition

國家自然科學基金項目(51575178;11572125)

2015-11-25 修改稿收到日期:2016-01-28

王廣斌 男,博士,副教授,1974年12月生

杜謀軍 男,碩士生,1992年12月生

E-mail:admjt163@163.com

TP277

A

10.13465/j.cnki.jvs.2016.20.012

猜你喜歡
故障診斷故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 全部免费特黄特色大片视频| 九九九久久国产精品| 白丝美女办公室高潮喷水视频| 国产在线98福利播放视频免费| 亚洲天堂精品在线| 成人在线第一页| 欧美不卡在线视频| 国产污视频在线观看| 国产精品极品美女自在线看免费一区二区| 国产一区二区三区在线无码| 欧美综合一区二区三区| 精品精品国产高清A毛片| 国产亚洲高清视频| 免费人成在线观看成人片| 日韩 欧美 国产 精品 综合| 99精品福利视频| 久久女人网| 国产精品综合久久久| 五月天丁香婷婷综合久久| 亚洲清纯自偷自拍另类专区| 麻豆国产原创视频在线播放| 欧美精品一二三区| 激情综合网激情综合| 欧美激情视频在线观看一区| 男人天堂亚洲天堂| 国产高清无码麻豆精品| 国产成人午夜福利免费无码r| 国产亚洲一区二区三区在线| 久草中文网| 日韩精品免费一线在线观看| 美女无遮挡被啪啪到高潮免费| 精品無碼一區在線觀看 | 91麻豆国产精品91久久久| 黄片一区二区三区| 青草精品视频| 亚洲大尺度在线| 色婷婷色丁香| 免费A级毛片无码免费视频| 国产黄色爱视频| 国产网站在线看| 中美日韩在线网免费毛片视频| 伊人无码视屏| 国产成人做受免费视频| 国产欧美视频一区二区三区| 伊人久久婷婷| 99在线免费播放| 欧美精品三级在线| 欧美精品在线看| 日韩在线网址| 欧美激情视频二区| 高清大学生毛片一级| 全部无卡免费的毛片在线看| 亚洲AV无码久久精品色欲| 国产精品视频白浆免费视频| 日韩AV无码免费一二三区| 国产成人福利在线视老湿机| 国产精品妖精视频| 71pao成人国产永久免费视频| 亚洲午夜国产片在线观看| 国产无码高清视频不卡| 无码人妻热线精品视频| 无码区日韩专区免费系列| 成人免费视频一区二区三区 | 精品国产电影久久九九| 国产一区二区三区精品久久呦| 国产免费福利网站| 国产亚洲精品无码专| 久久免费视频6| 91色在线观看| 免费高清a毛片| 国产福利一区二区在线观看| 久久特级毛片| 男人天堂亚洲天堂| 亚洲免费福利视频| 国产成人亚洲无码淙合青草| 久久精品一品道久久精品| 亚洲swag精品自拍一区| 波多野一区| 婷婷在线网站| 日韩午夜伦| 亚洲国产理论片在线播放| 久久精品人人做人人爽电影蜜月 |