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

基于四元數(shù)理論與流形學(xué)習(xí)的多通道機械故障信號分類方法

2016-06-09 08:10:58易燦燦
武漢科技大學(xué)學(xué)報 2016年6期
關(guān)鍵詞:故障診斷分類故障

何 博, 呂 勇, 易燦燦,黨 章

(武漢科技大學(xué)機械自動化學(xué)院,湖北 武漢,430081)

?

基于四元數(shù)理論與流形學(xué)習(xí)的多通道機械故障信號分類方法

何 博, 呂 勇, 易燦燦,黨 章

(武漢科技大學(xué)機械自動化學(xué)院,湖北 武漢,430081)

提出一種基于增廣四元數(shù)矩陣奇異值分解與流形學(xué)習(xí)正交鄰域保持嵌入算法的多通道機械故障信號分類方法,通過引入四元數(shù)來耦合4個通道信號,并且利用四元數(shù)乘方的性質(zhì)對數(shù)據(jù)進行增廣處理,充分利用各通道信息并挖掘通道之間的相關(guān)性,從而減少因故障特征信息丟失對分類結(jié)果的影響。此外,針對傳統(tǒng)奇異譜分析提取特征參數(shù)的分類效果受噪聲影響較大的問題,引入正交鄰域保持嵌入算法對奇異值序列進行維數(shù)約簡,最后使用分類器完成故障分類。對仿真信號的分類結(jié)果表明,在強噪聲背景下,相較于單通道奇異譜分析方法和機械故障信號中常用的排列熵方法,本文提出的方法分類效果更好。將其應(yīng)用于更為復(fù)雜的實測軸承故障信號的分類與識別中,同樣有著較好的效果。

故障診斷;信號處理;四元數(shù);奇異值分解;流形學(xué)習(xí);故障分類

機械設(shè)備故障信號常常呈現(xiàn)出強非線性、非平穩(wěn)性,并且包含著與設(shè)備運行狀態(tài)相聯(lián)系的動力學(xué)特性[1],這使得許多傳統(tǒng)信號處理方法難以達到精確的故障診斷效果。奇異譜分析(singular spectrum analysis, SSA)是近年來興起的一種非線性信號處理方法,在機械故障診斷中被廣泛應(yīng)用[2-3]。然而,與目前絕大部分機械故障診斷方法一樣,SSA算法只針對單通道故障信號數(shù)據(jù),這必然會造成有效故障特征信息的遺漏,增大故障診斷結(jié)果的不確定性[4-6]。同時采集并分析多個通道的信號數(shù)據(jù)可以獲得更多包含故障特征的有用信息,從而提高故障診斷結(jié)果的準確性和置信度。Enshaeifar等[7]用四元數(shù)理論耦合四通道信號,并采用增廣四元數(shù)奇異值分解(augmented quaternion singular value decomposition, AQSVD)算法將信號分解為多個子空間信號,最后通過奇異譜分析選擇部分子空間進行重構(gòu)。該方法被成功地用于生物醫(yī)學(xué)信號處理,能準確實現(xiàn)信號的降噪及盲源分離。

流形學(xué)習(xí)近年來一直是模式識別領(lǐng)域的研究熱點,常常被用于高維數(shù)據(jù)的維數(shù)約簡,可挖掘出高維數(shù)據(jù)的低維結(jié)構(gòu)[8]。目前主要非線性流形學(xué)習(xí)算法包括局部線性嵌入算法、拉普拉斯特征映射算法、近鄰域保持嵌入算法等[9]。這些方法能夠較好地找到高維數(shù)據(jù)集所在流形的幾何特性以及非線性結(jié)構(gòu),并且挖掘原始數(shù)據(jù)的低維分布。正交鄰域保持嵌入(orthogonal neighborhood preserving embedding, ONPE)是一種根據(jù)鄰域保持嵌入改進的算法[10],該算法使得到的投影向量具有正交性,保留了數(shù)據(jù)的非線性特征,避免了局部子空間的結(jié)構(gòu)失真,對非線性高維數(shù)據(jù)具有較好的維數(shù)約簡以及分類效果。

為此,本文提出一種基于AQSVD和ONPE的多通道機械故障信號分類方法,并與傳統(tǒng)的排列熵算法及單通道奇異譜分析法相比較,以驗證其對強噪聲背景下故障信號提取及分類的有效性。

1 基于增廣四元數(shù)矩陣奇異值分解的故障信號特征提取

1.1 四元數(shù)的基本概念與性質(zhì)

四元數(shù)最早是由愛爾蘭數(shù)學(xué)家哈密頓提出的概念[11],它的數(shù)學(xué)意義是一種超復(fù)數(shù),一個四元數(shù)包括一個實部分量和3個虛部分量,其形式為[7]:

x=xa+ixb+jxc+kxd

(1)

其中,i、j、k表示虛部正交單位向量,具有如下性質(zhì):

ij=-ji=k,jk=-kj=i,ki=-ik=j

(2)

i2=j2=k2=ijk=-1

(3)

四元數(shù)域的另一個重要概念即為“自逆映射”,即四元數(shù)乘方,其計算公式如下[7]:

xi=-ixi=xa+ixb-jxc-kxd

xj=-jxj=xa-ixb+jxc-kxd

xk=-kxk=xa-ixb-jxc+kxd

(4)

1.2 基于增廣四元數(shù)矩陣奇異值分解的特征信號提取方法

對于一個四通道數(shù)據(jù),可通過四元數(shù)的概念將其耦合。本文將每個通道信號數(shù)據(jù)都視為四元數(shù)的其中一個分量,將4個通道數(shù)據(jù)耦合表示整體故障信息,即用一個超復(fù)數(shù)單通道信號來表示原始四通道信號,這樣就可用單通道計算方法來處理四通道信號,有效地利用了多通道信息。將4×N的輸入信號變?yōu)?×N的四元數(shù)序列,令該四元數(shù)序列為x=[x1,x2,…,xN],通過相空間重構(gòu),構(gòu)建一個Hankel矩陣W:

(5)

式中:N為信號源長度;L為窗口長度,1

不同于傳統(tǒng)奇異值分解算法,為了充分挖掘四通道之間的相關(guān)性,AQSVD算法中將對軌道矩陣W進行增廣處理。根據(jù)四元數(shù)的“自逆映射”(式(4)),可由W計算出3個乘方數(shù)軌道矩陣Wi、Wj和Wk。利用3個軌道矩陣對W進行增廣處理,得到增廣四元數(shù)矩陣Wa=[WT,WiT,WjT,WkT]T,Wa∈H4L×K,其中H為四元數(shù)域。然后計算增廣矩陣的協(xié)方差矩陣Ca∈Η4L×4L,計算公式如下:

Ca=E{WaWaH}=

(6)

式中:CWW為標準協(xié)方差矩陣,矩陣CWi、CWj和CWk是互補矩陣。

對協(xié)方差矩陣進行奇異值分解,由其特征向量和左右奇異值矩陣來表示,公式如下:

Ca=UΛ1/2VT

(7)

(8)

λj按逐個減小的方式排列為(λ1>λ2>…>λL),序列中的奇異值即為代表機械故障特征的參數(shù)。至此,特征提取已經(jīng)完成,奇異譜分解中的重構(gòu)部分將不被用到。

2 正交鄰域保持嵌入算法

2.1 正交化處理

設(shè)提取的原信號特征參數(shù)有n個,樣本維數(shù)為D,則原始特征集為X={x1,…,xn}∈D×n。流形學(xué)習(xí)進行維數(shù)約簡首先要計算出轉(zhuǎn)換矩陣。機械故障診斷中常常需要較多的特征參數(shù),即特征維度大于樣本數(shù),這樣得到的矩陣XXT為奇異矩陣,計算出來的轉(zhuǎn)換矩陣為非正交,不能保證數(shù)據(jù)集的幾何結(jié)構(gòu),因此需要對原特征集進行正交化處理。

將原始特征集X={x1,…,xn}∈D×n投影到PCA子空間,舍棄特征值為0的成分,得到正交處理后的數(shù)據(jù)X={x1,…,xn}∈d×n,式中d為降維后的維度,則矩陣XXT即為非奇異矩陣。將高維數(shù)據(jù)轉(zhuǎn)換到低維數(shù)據(jù)的轉(zhuǎn)換矩陣記為Apca。

2.2 鄰域權(quán)值矩陣構(gòu)造

利用K近鄰算法尋找數(shù)據(jù)集X中任意一點xi的k個近鄰點xj(j=1,…,k),計算其重構(gòu)權(quán)值矩陣Wij(即節(jié)點xj對節(jié)點xi的重構(gòu)貢獻度,兩節(jié)點距離越近則Wij越大),然后通過最小化重構(gòu)誤差函數(shù)構(gòu)造權(quán)值矩陣:

(9)

(10)

其中,xu、xv同樣為xi的近鄰點。權(quán)值矩陣Wij代表了原始數(shù)據(jù)集鄰域的結(jié)構(gòu)關(guān)系。

2.3 計算正交鄰域保持嵌入

正交領(lǐng)域保持嵌入的映射矩陣記為AONPE=(a1,a2,…,ad)。存在一種極限情況,即經(jīng)過PCA投影的數(shù)據(jù)集X映射到一維直線上,映射結(jié)果為y=(y1,y2,…,yn)∈1×n,y=aTX。由于數(shù)據(jù)從高維嵌入低維后其幾何分布不變,所以直線上每一數(shù)據(jù)都可視為其鄰域點的線性重構(gòu)。為了使重構(gòu)誤差最小,設(shè)計其最小懲罰函數(shù)為

(11)

定義過渡矩陣z和I:z=yT-WyT=(I-W)yT,那么有:

aTX(I-W)T(I-W)XTa

(12)

由于aTXXTa=1,令aTX(I-W)T(I-W)XTa=λ=λaTXXTa,可得:

(XXT)-1X(I-W)T(I-W)XTa=λa

(13)

其中λ為矩陣(XXT)-1X(I-W)T(I-W)XT的最小特征值,代表最小重構(gòu)誤差,其對應(yīng)的特征向量為a1。因此,ONPE的局部保持函數(shù)即為

(14)

ONPE的目標函數(shù)為

(15)

(16)

·

(XXT)-1X(I-W)T(I-W)XTak=λak

(17)

其中A(k-1)=[a1,a2,…,ak-1],S(k-1)=(A(k-1))T(XXT)-1A(k-1)(k=2,3,…,d)。

由式(17)迭代可求出映射矩陣AONPE,則低維特征集Y為

Y=(ApcaAONPE)TX

(18)

至此,實現(xiàn)了由高維數(shù)據(jù)集到低維數(shù)據(jù)集的正交鄰域保持投影。

3 模擬故障信號分析

3.1 模擬故障信號的設(shè)置

軸承故障信號常常表現(xiàn)為調(diào)制信號與噪聲信號的線性疊加,因此,本文將使用4種調(diào)制信號線性疊加并添加3 dB強背景噪聲來模擬多通道信號數(shù)據(jù)。調(diào)制信號的表達式如下:

ch1=0.2cos(2πf1t+5)+noise

ch2=0.25cos(2πf2t-15)+noise

ch3=0.25cos(2πf3t)+noise

ch4=0.25sin(2πf4t+15)+noise

(19)

式中:chi(i=1,2,3,4)表示通道編號;fi(i=1,2,3,4)為模擬故障信號的頻率;t為時間參數(shù);noise為3 dB的噪聲。

采樣頻率為1024 Hz,每個信號采樣長度均為20 000,則四通道仿真信號S為

(20)

為了驗證分類效果,設(shè)置3組模擬故障信號,其頻率如表1所示。

表1 3組模擬故障信號的頻率(單位:Hz)

由于多通道故障信號的各個通道之間存在一定關(guān)聯(lián)性,對于每組模擬故障信號,生成一個隨機矩陣C,使用隨機矩陣將原始多通道信號混合,則模擬四通道信號為

Xi=C×Si(i=1,2,3)

(21)

每組仿真故障信號的頻譜圖如1所示。由圖1中可見,在強噪聲背景下,各組模擬信號中設(shè)置的故障頻率均無法在頻譜圖上被找出。

將每類故障平均分為20組,則每類故障都生成了20個長度為1000的四通道信號樣本。為了驗證本文提出的基于AQSVD與ONPE方法的可靠性,將對其與排列熵算法及單通道奇異值分解算法的分類效果進行對比。

3.2 基于排列熵算法的仿真信號分類

排列熵算法是一種非線性數(shù)據(jù)處理方法,熵的大小可描述一維數(shù)據(jù)復(fù)雜度,并且可以檢測到系統(tǒng)的動態(tài)變化特征[12-13]。因此,排列熵算法近幾年被廣泛應(yīng)用于機械故障診斷中。

對于任意四通道信號樣本,計算每一通道信號的排列熵,令算法中嵌入維度為5,延時時間為1,結(jié)果如圖2所示。由圖2可見,無論在哪一通道,3組故障之間都沒有明顯界限,不同故障之間區(qū)分度不明顯,可見使用排列熵算法對仿真信號分類的效果較差,同時充分表明了排列熵算法在強噪聲背景下效果的局限性。

圖1 仿真信號頻譜圖

(a)通道1(b)通道2

(c)通道3(d)通道4

圖2 樣本排列熵值

Fig.2 Permutation entropy of samples

3.3 基于單通道奇異值分解算法的仿真信號分類

傳統(tǒng)單通道奇異值分解算法中,常使用峭度作為分類的特征參數(shù)[14]。對于任意四通道信號樣本,分別對其每一通道信號進行奇異值分解,得到奇異值序列,選取奇異值序列中貢獻度之和大于95%的前n個奇異值組成新的奇異值序列,計算該序列的峭度。

采用單通道奇異值分解算法對本文仿真信號進行分析,計算結(jié)果如圖3所示。由圖3可見,通道1和通道2中第三組仿真信號可以被準確區(qū)分開,而第一組與第二組數(shù)據(jù)混雜在一起,沒有明顯界限;通道3中也僅有第二組仿真信號可以被區(qū)分開;通道4是4個通道中效果最好的一個通道,其中第一組與其他兩組有著明顯界限,但第二、三組之間仍有少量樣本混疊。可見單通道算法明顯存在偶然性,其分類效果也并不理想。

(a)通道1

(b)通道2

(c)通道3

(d)通道4

圖3 樣本奇異值序列峭度

Fig.3 Kurtosis of sample singular value sequence

3.4 基于AQSVD與ONPE算法的仿真信號分類

對于任意四通道信號樣本,利用本文提出的基于AQSVD與ONPE的算法,將原始數(shù)據(jù)進行四元數(shù)耦合處理,生成一個四元數(shù)序列;繼而對序列進行相空間重構(gòu),形成Hankel矩陣,并依據(jù)四元數(shù)乘方的性質(zhì)對矩陣進行增廣處理;然后計算協(xié)方差矩陣,對協(xié)方差矩陣進行奇異值分解,得到奇異值序列。選取奇異值序列中貢獻度之和大于95%的前n個奇異值組成新的奇異值序列,對新的奇異值序列使用ONPE算法進行維數(shù)約簡,之后即可使用分類器進行故障分類。

圖4所示為使用ONPE算法將經(jīng)AQSVD步驟計算得出的奇異值序列降至二維后的樣本分布。由圖4中可看出,同組故障都有明顯類聚的特點,不同故障之間界限清晰,分類效果優(yōu)于前兩種方法。并且由于使用了基于四元數(shù)理論的AQSVD算法,避免了單通道信號遺漏故障信息且具有偶然性的缺點。

圖4 模擬信號的二維分布圖

Fig.4 2D distribution of the analog signals

4 實測故障信號分析

仿真信號的模擬環(huán)境比較理想,而機械設(shè)備的工作環(huán)境通常較為復(fù)雜,并包含了更多不確定因素。因此,為驗證本文所提出方法在實際機械故障信號分類中的有效性,采用美國辛辛那提大學(xué)智能維護系統(tǒng)中心提供的標準軸承故障測試數(shù)據(jù)[15]作進一步分析。所選3組故障信號的描述如表2所示。

由于第二組與第三組同為四通道信號,而第一組為八通道信號,故第一組取四個豎直方向的振動信號組成四通道數(shù)據(jù)。每組故障隨機取20組樣本,采樣點數(shù)為6000。

使用本文提出的基于AQSVD與ONPE的方法處理軸承故障數(shù)據(jù),并將奇異值序列維度降至二維,樣本的空間分布情況如圖5所示。由圖5可知,經(jīng)過基于AQSVD及ONPE算法后,3組軸承故障樣本在二維空間中的投影具有明顯的類聚特點,不同種類的故障樣本之間有著明顯界限。該方法無需對不同通道的信號進行篩選,有效簡化了分類步驟,并且降低了識別偶然性。由于對奇異值序列組成的特征集進行了維數(shù)約簡,在使用分類器進行分類時計算量將被大幅度削減。

表2 3組軸承故障信號的描述

圖5 實測信號的二維分布

Fig.5 2D distribution of the measured signals

5 結(jié)語

本文在傳統(tǒng)奇異值分解故障診斷方法的基礎(chǔ)上,提出了基于增廣四元數(shù)矩陣奇異值分解及正交鄰域保持嵌入的分類方法。針對單通道數(shù)據(jù)處理方法會造成有效信息遺漏這一缺點以及傳統(tǒng)奇異值分解算法在機械故障診斷應(yīng)用中的缺陷,提出了基于增廣四元數(shù)矩陣奇異值分解的算法,耦合四通道數(shù)據(jù)并充分挖掘各通道之間相關(guān)性,減少了有效信息的遺漏。針對傳統(tǒng)奇異譜分析提取的特征參數(shù)分類效果差的缺點,引入了正交鄰域保持嵌入算法對高維數(shù)據(jù)進行維數(shù)約簡并完成分類。對仿真信號以及軸承實測故障信號的分類結(jié)果表明,基于增廣四元數(shù)矩陣奇異值分解及正交鄰域保持嵌入的分類方法在機械故障診斷應(yīng)用中具有較好的效果,識別準確度要明顯高于傳統(tǒng)方法。

[1] ChenYS.Nonlineardynamicalprincipleofme-chanical fault diagnosis[J].Chinese Journal of Mechanical Engineering,2007,43(1):6.

[2] Rocco S C M. Singular spectrum analysis and forecasting of failure time series[J]. Reliability Engineering and System Safety,2013,114: 126-136.

[3] 盧德林,郭興明. 基于奇異譜分析的心音信號小波包去噪算法研究[J]. 振動與沖擊, 2013,32(18):63-69.

[4] Elsner J B, Tsonis A A. Singular spectrum analysis: a new tool in time series analysis[M].New York: Springer Science & Business Media,2013.

[5] VautardR,YiouP,GhilM.Singular-spectrumanal-ysis:a toolkit for short, noisy chaotic signals[J]. Physica D: Nonlinear Phenomena, 1992, 58(1-4): 95-126.

[6] Chouikhi S,Korbi I E, Ghamri-Doudane Y, et al. Fault tolerant multi-channel allocation scheme for wireless sensor networks[C]//2014 IEEE Wireless Communications and Networking Conference (WCNC). IEEE, 2014: 2438-2443.

[7] Enshaeifar S, Kouchaki S,Took C C, et al. Quaternion singular spectrum analysis of electroencephalogram with application in sleep analysis[J]. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2016, 24(1): 57-67.

[8] 蔣全勝,李華榮,黃鵬. 一種基于非線性流形學(xué)習(xí)的故障特征提取模型[J]. 振動與沖擊,2012,31(23):132-136.

[9] 王雷. 基于流形學(xué)習(xí)的滾動軸承故障診斷若干方法研究[D]. 大連:大連理工大學(xué), 2012.

[10]宋濤,湯寶平,李鋒. 基于流形學(xué)習(xí)和K-最近鄰分類器的旋轉(zhuǎn)機械故障診斷方法[J]. 振動與沖擊, 2013, 32(5):149-153.

[11]張榮輝,賈宏光,陳濤,等. 基于四元數(shù)法的捷聯(lián)式慣性導(dǎo)航系統(tǒng)的姿態(tài)解算[J]. 光學(xué)精密工程,2008,16(10):1963-1970.

[12]鄭近德,程軍圣,楊宇. 基于LCD和排列熵的滾動軸承故障診斷[J]. 振動、測試與診斷, 2014,34(5): 802-806.

[13]饒國強,馮輔周,司愛威,等. 排列熵算法參數(shù)的優(yōu)化確定方法研究[J]. 振動與沖擊, 2014,33(1):188-193.

[14]唐貴基,龐彬,劉尚坤. 基于奇異差分譜和平穩(wěn)子空間分析的滾動軸承故障診斷[J]. 振動與沖擊, 2015,34(11):83-87,115.

[15]Qiu H, Lee J, Lin J, et al. Wavelet filter-based weak signature detection method and its application on rolling element bearing prognostics[J]. Journal of Sound and Vibration, 2006, 289: 1066-1090.

[責(zé)任編輯 鄭淑芳]

A novel method for multi-channel mechanical fault signal classification based on quaternion and manifold learning

HeBo,LvYong,YiCancan,DangZhang

(College of Machinery and Automation, Wuhan University of Science and Technology, Wuhan 430081, China)

A novel method for multi-channel mechanical fault signal classification based on augmented quaternion matrix singular value decomposition and orthogonal neighborhood preserving embedding algorithm of manifold learning is proposed. Quaternion is used to couple four channel signals, and the nature of the quaternion power is employed for augmented processing of the data. The correlation between channels is made use of, and the information from each channel is employed to offset the negative influence of loss of characteristic information of faults on classification. Considering that the traditional classification method that uses singular spectrum analysis to extract characteristic parameters is seriously affected by noise, the orthogonal neighborhood preserving embedding algorithm is used to reduce the dimension of singular value sequence. Finally, the classifier is used to classify faults. The results show that, with the background of strong noise, the proposed method is superior to the traditional single-channel singular spectrum analysis method and the method of permutation entropy in fault classification. Applied to the complex identification and classification of real bearing fault signals, the proposed method shows good performance.

fault diagnosis; signal processing; quaternion; singular value decomposition; manifold learning; fault classification

2016-09-20

國家自然科學(xué)基金資助項目(51475339);武漢科技大學(xué)冶金裝備及其控制教育部重點實驗室開放基金資助項目(2015B11).

何 博(1990-),男,武漢科技大學(xué)碩士生.E-mail:sainthebo@163.com

呂 勇(1976-),男,武漢科技大學(xué)教授,博士生導(dǎo)師.E-mail:lvyong@wust.edu.cn

TH133.3;TH165.3

A

1674-3644(2016)06-0455-07

猜你喜歡
故障診斷分類故障
分類算一算
故障一點通
分類討論求坐標
數(shù)據(jù)分析中的分類討論
教你一招:數(shù)的分類
奔馳R320車ABS、ESP故障燈異常點亮
因果圖定性分析法及其在故障診斷中的應(yīng)用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 欧美亚洲国产精品久久蜜芽| 免费aa毛片| 欧洲成人在线观看| 亚洲国产亚洲综合在线尤物| 精品91在线| 97在线免费视频| 国产日本欧美亚洲精品视| 第一页亚洲| 亚洲美女久久| 97国产在线视频| 2022精品国偷自产免费观看| 自拍亚洲欧美精品| 亚洲成年人网| 国产9191精品免费观看| 在线视频亚洲欧美| 久久综合五月| 伊大人香蕉久久网欧美| 不卡视频国产| 成人精品免费视频| 久青草国产高清在线视频| 成人在线观看一区| 欧洲一区二区三区无码| 日韩免费中文字幕| 久久精品波多野结衣| 国产一二三区视频| 四虎精品国产AV二区| 欧美亚洲欧美| 日本精品视频一区二区 | 91九色国产porny| 亚洲国产高清精品线久久| 91午夜福利在线观看精品| 日韩精品资源| 极品尤物av美乳在线观看| 亚洲日韩Av中文字幕无码| 国产91在线免费视频| 欧美国产日韩在线观看| 国产欧美日韩另类精彩视频| 色悠久久久| 一级看片免费视频| 蜜臀av性久久久久蜜臀aⅴ麻豆 | 一本大道视频精品人妻 | 久久这里只有精品免费| 免费一级毛片完整版在线看| JIZZ亚洲国产| 四虎国产精品永久在线网址| 一级毛片在线直接观看| 亚洲欧美日韩另类| 国产精品亚欧美一区二区| 国产微拍精品| 六月婷婷激情综合| 欧美午夜视频在线| 国产老女人精品免费视频| 99久视频| 九九久久精品免费观看| 日韩不卡免费视频| 色精品视频| 57pao国产成视频免费播放| 国产一区二区视频在线| 国产精品手机视频| 97精品国产高清久久久久蜜芽| 日本一区二区不卡视频| 亚洲视频一区在线| 日本成人精品视频| 91精品人妻一区二区| 国产视频资源在线观看| 女人毛片a级大学毛片免费| 99精品一区二区免费视频| 国产91精品久久| 国产综合亚洲欧洲区精品无码| 色噜噜狠狠狠综合曰曰曰| 啪啪国产视频| 成人字幕网视频在线观看| 国产91高跟丝袜| 国产无遮挡裸体免费视频| 九九九国产| 久久a毛片| 国产91透明丝袜美腿在线| 精品国产福利在线| 久久精品国产国语对白| 亚洲天堂成人| 成人福利在线视频免费观看| 亚洲丝袜中文字幕|