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

基于圖像形狀特征和LLTSA的故障診斷方法

2016-06-23 09:38:30張前圖房立清
振動與沖擊 2016年9期
關鍵詞:故障診斷

張前圖, 房立清

(1.駐356廠軍代室,昆明 650114; 2.軍械工程學院 火炮工程系,石家莊 050003)

基于圖像形狀特征和LLTSA的故障診斷方法

張前圖1,2, 房立清2

(1.駐356廠軍代室,昆明650114; 2.軍械工程學院 火炮工程系,石家莊050003)

摘要:針對滾動軸承故障診斷問題,提出了一種基于圖像形狀特征和線性局部切空間排列(LLTSA)的故障診斷方法。首先采用SDP(Symmetrized Dot Pattern)方法對時域信號進行變換,得到極坐標空間下的雪花圖像,在分析圖像特點的基礎上,從圖像處理的角度初步提取出圖像的形狀特征;然后利用LLTSA對初步提取的特征進行維數約簡以提取低維特征;最后采用支持向量機(SVM)對低維特征進行分類評估。滾動軸承的故障診斷實驗表明圖像形狀特征能夠表征軸承的狀態,經LLTSA約簡后特征數據的復雜度得到降低,且具有更好的聚類效果,而SVM對軸承4種狀態的識別率也達到了100%,驗證了該方法的有效性。

關鍵詞:SDP;形狀特征;線性局部切空間排列;支持向量機;故障診斷

滾動軸承是機械設備中最常見、最容易損壞的部件之一,其正常與否決定著機械設備狀態的好壞。振動信號作為滾動軸承故障信息的載體,包含著軸承運行狀態的重要信息,對其進行分析處理是實現軸承故障診斷的常用方法。目前對振動信號進行分析主要從時域、頻域以及時頻域等方面入手,采用不同的信號處理方法,提取能夠表征軸承運行狀態的特征,取得了不錯的效果[1-3]。近年來,也有不少學者[4-5]采用相應的方法將振動信號可視化,進而以圖像處理的方式,提取出能夠表征其運行狀態的特征并加以分類識別。目前從圖像處理的角度進行故障診斷時,圖像的來源大多都是對信號進行時頻分析后得到的時頻圖,但時頻變換計算復雜,且從圖像提取的特征往往是一些統計特征量,選擇沒有規律性且計算復雜,易丟失重要的時頻信息。

SDP[6-8](Symmetrized Dot Pattern)作為一種圖像生成方法,相比于時頻分析具有理論簡單、運算速度快的特點,能夠將一維的時間序列通過相應的計算公式變換到極坐標空間下的對稱雪花圖,不同信號之間的差異可以通過雪花花瓣形狀的不同得以體現。根據這一特點,可以將滾動軸承振動信號進行SDP變換,然后從圖像處理的角度提取出反映花瓣形狀的特征,采用相應的模式識別方法,以達到故障診斷目的。然而提出的特征中往往包含有冗余信息,影響故障診斷的精度,因此有必要采用維數約簡算法對特征進行二次提取。線性局部切空間排列[9](Liner Local Tangent Space Alignment,LLTSA)作為一種新型流行學習算法,相比于局部保持投影[10](Locality Preserving Projection,LPP)等算法具有更好的降維效果[11],能夠充分挖掘高維非線性數據中的低維特征,具有很好的非線性復雜信息處理能力。

基于以上分析,本文提出了一種基于圖像形狀特征和LLTSA的滾動軸承故障診斷方法。即對滾動軸承振動信號進行SDP變換,然后從SDP圖像中提取形狀特征并采用LLTSA進行維數約簡,最后將提取的低維特征輸入支持向量機(Support Vector Machine,SVM)進行分類識別。滾動軸承的故障診斷實例驗證了該方法的有效性。

1SDP方法和圖像形狀特征

1.1SDP方法的基本原理

在信號的離散數據采樣序列中,xn是在時刻n時的幅值,xn+l是時刻n+l時的幅值,通過SDP變換公式,可以將其變換到極坐標空間P(r(n),Θ(n),φ(n))中的點。圖1[7]為SDP方法的基本原理圖。圖中r(n)為極坐標的半徑;Θ(n)為極坐標逆時針沿初始線旋轉的角度;φ(n)為極坐標順時針沿初始線旋轉的角度。這三個變量的具體計算公式如下:

(1)

(2)

(3)

式中,xmin為采樣序列的最小值;xmax為采樣序列最大值;θ為鏡像對稱平面旋轉角;g為角度放大因子;l為時間間隔參數。

圖1 SDP方法原理圖Fig.1 Principle of SDP method

1.2SDP的性質及參數選擇

從SDP的計算公式可以看出,SDP方法的重點就在于極坐標中點位置的確定。在極坐標空間中,時域信號兩個時間間隔為l的幅值決定了點的位置,假定l=1,如果信號中主要包含高頻成分,則時域波形中i處幅值xi和i+1處幅值xi+1的差異就較大,而由SDP得到的極坐標空間中的點就會有較小的偏轉角度和較大的半徑,反之亦然[12]。因此,信號頻率間的差異主要體現在極坐標中點分布位置和曲率的不同。圖2為3個頻率分別為100 Hz、200 Hz和400 Hz的正弦仿真信號以及高斯白噪聲信號的SDP圖像(θ=60°,g=40°,l=5)。從中可以看出,隨著頻率的增大,SDP圖像的點的位置發生了較大變化,使得由點連接成的花瓣逐漸變得飽滿,重疊部分也開始增加,SDP圖像點的分布位置直觀反映了正弦信號頻率的變化。對比正弦信號和高斯白噪聲信號的SDP圖像可知,從周期信號到非線性信號的變化同樣可以通過SDP圖像的差異得到體現。

圖2 仿真信號SDP圖像Fig.2 SDP image of simulation signals

在SDP方法中,θ,g,l這三個參數的選擇尤為重要。通常選取θ=60°,其鏡像對稱平面為0°、60°、120°、180°、240°、300°。在極坐標中,這六個平面組成雪花狀的六邊形。如果鏡像對稱平面旋轉角θ過大,在極坐標中顯示的圖像不發生重疊,圖像對稱性較差;θ過小則圖像的重疊較多,導致特征不明顯。當θ為60°時,能夠對信號特征清晰描述,便于特征的提取。通過大量實驗發現[8],不同信號之間的細微區別主要依靠g和l的選取,一般g要小于θ,取20°~50°為宜,l的值在1~10的范圍內最佳。

1.3圖像形狀特征參數

圖像的特征識別大致可歸納為三類[13]:顏色或灰度的統計特征;紋理、邊緣特征;形狀特征。由于不同信號SDP圖像的差異主要體現在幾何形狀的不同上,因此采用形狀特征對SDP圖像進行描述。圖像形狀特征主要包括區域面積、區域邊界、區域質心、與區域具有相同二階中心矩的橢圓、方向角等參數,圖3為以上各參數的示意圖。上述幾種形狀特征參數定義如下[14]:

(1) 區域面積:描述了區域的大小,計算區域面積的一種簡單方法就是對所屬區域的像素進行計算。圖3(a)中圖像區域由5個白色像素點組成,按定義可知其區域面積為5像素。

(2) 區域質心:根據區域內所有點計算出來的坐標就是質心的坐標。圖3(a)中灰色圓點即為白色像素圖像區域質心。質心坐標的計算公式如下:

(4)

式中A為區域面積。

(3) 區域邊界:是指包含圖像區域最小的外接矩形,圖3(a)中白色矩形為白色像素圖像的區域邊界。矩形的大小由矩形左上角頂點坐標以及矩形的長和寬決定。

(4) 與區域具有相同二階中心矩的橢圓:圖3(b)中橢圓即為與白色像素區域具有相同二階中心矩的橢圓;圖3(c)中長直線為橢圓長軸,短直線為橢圓短軸,長直線上兩個紅點為橢圓的焦點。

(5) 方向角:是指橢圓長軸與水平線的夾角。圖3(c)中長直線與水平虛線的夾角即為方向角。

圖3 形狀特征參數示意圖Fig.3 Sketch map of shape feature

2線性局部切空間排列

線性局部切空間排列的主要思想就是尋找一個轉換矩陣A將RD空間中具有N個點的含噪數據集XORG(即本文中的故障樣本集)映射為Rd空間數據集Y=[y1,…,yN],即:

Y=ATXORGHN(d

(5)

式中HN=I-eeT/N為中心矩陣,I為單位矩陣,e為所有元素為1的N維列向量。Y為XORG潛在的d維非線性流行。

LLTSA主要有以下幾個步驟[15]:

(1) 主成分分析投影。使用PCA將數據集映射到主體子空間,用APCA表示PCA的轉換矩陣,同時為了表示方便,仍采用X表示PCA子空間數據集。

(2) 確定領域。采用K-近鄰法(KNN)搜索數據點xi的領域,即首先采用歐式距離構造所有數據點的距離矩陣,然后通過分析距離矩陣找到數據點xi的k個同類近似點xj。

(3) 獲取局部信息。計算由XiHk(Xi=[xi1,…,xik])的d個最大特征值對應d個特征矢量構成的矩陣Vi。其中Hk=I-eeT/k。

(4) 構造排列矩陣。通過局部累加構造矩陣B為:

B(Ii,Ii)←B(Ii,Ii)+WiWiT

(i=1,2,…,N)

(6)

初始化B=0,式中Ii={ii,…,ik}表示xi的k個近鄰點的索引集,Wi=Hk(i-ViViT)(i=1,2,…,N)。

(5) 計算映射。計算廣義特征問題的特征值和特征矢量:

XHNBHNXTα=λXHNXTα

(7)

與特征值λ1,λ2,…,λn(λ1<λ2…<λn)對應的特征矢量解為α1,α2,…,αd,則ALLTSA=(α1,α2,…,αd),因此轉換矩陣A=APCA×ALLTSA,X→Y=ATXORGHN。

3形狀特征-LLTSA-SVM軸承故障診斷方法

滾動軸承發生不同故障后,故障信號中會出現不同頻率成分的周期性沖擊,使得信號的幅值和頻率成分發生改變,進而影響極坐標空間中點的分布情況。而前面定義的形狀特征又主要由點的分布情況決定,點分散的越開,重疊的點越少,則區域面積和區域邊界就會變大,使得區域質心位置和橢圓位置發生改變,而方向角也會隨著橢圓位置的改變而變化,反之亦然。因此,滾動軸承振動信號中頻率成分的改變引起的點分布位置的變化,就建立了軸承振動信號與形狀特征之間的一種內在聯系,使得從振動信號SDP圖像中提取的形狀特征能夠對軸承振動信號進行表征。同時,LLTSA具有較強的非線性信息處理能力,通過它對形狀特征進行維數約減,能夠剔除形狀特征中的冗余信息,可以獲得一個維數低、敏感度高且聚類效果好的低維特征。另一方面,支持向量機(SVM)是建立在統計學基礎上的機器學習方法,具有小樣本學習、泛化能力強等優點,是目前故障診斷領域比較常用的方法[2, 16]。因此,把SVM作為診斷方法的一部分,將其用于低維特征的分類識別,以達到故障診斷的目的。

基于圖像形狀特征-LLTSA-SVM的滾動軸承故障診斷方法的流程如圖4所示,具體步驟如下:

圖4 形狀特征和LLTSA軸承故障診斷方法流程圖Fig.4 Flow chart ofbearing fault diagnosis smethod based image shape feature and LLTSA

(1) 對滾動軸承故障信號首先進行降噪處理,以突出故障信號的特征和提高信噪比,然后對降噪后信號進行SDP變換,得到SDP雪花圖像。進行SDP變換時,采用圖像相關分析確定g和l的取值,即g在20°~50°的范圍內步長為10°,l在1~10的范圍內步長為1,分別固定g和l的值,由于需要將不同故障狀態進行區分,因此將各故障狀態的雪花圖進行兩兩相關分析,取各相關系數之和最小時對應的g和l作為最優值。圖像相關性的計算公式為:

(2) 根據前面定義的圖像形狀特征,對各SDP圖像進行特征提取,得到一個維數為11的圖像形狀特征集。

(3) 將訓練樣本和測試樣本特征集同時輸入LLTSA進行訓練,求得轉換矩陣。

(4) 對訓練樣本和測試樣本進行維數約簡,分別得到d維的非線性主流行(即經過約簡的特征集),其中1≤d≤11。

(5) 將約簡后的訓練樣本輸入SVM中進行訓練,而后將約簡后的測試樣本輸入訓練好的SVM中進行測試,得到診斷結果。

4實驗驗證與結果分析

本文采用美國CWRU(Case Western Reserve University)軸承數據中心網站的滾動軸承故障數據來驗證本文所提方法的有效性。測試軸承為深溝球軸承,軸承型號為SKF6205,軸承損傷是通過電火花加工技術人為設置,并通過加速度傳感器采集振動信號。實驗中,實驗轉速為1 797 r/min,載荷為0,故障損傷直徑為0.177 8 mm,采樣頻率為12 kHz,對于軸承正常(N)、內圈故障(IF)、滾動體故障(BF)和外圈故障(OF)4種狀態,每種狀態測取20組數據樣本。通過多次實驗發現,數據樣本點數為4 096時,得到的圖像離散點少且花瓣飽滿,故樣本長度確定為4 096。

采用SDP方法對時域信號進行變換,利用圖像相關分析確定g和l的取值,通過計算,發現g=40°以及l=2時得到相關系數之和最小。圖5為軸承4種狀態在最優g和l值下的SDP圖像,從中可以看出,軸承4種狀態的SDP圖像具有明顯的差異,通過肉眼就可以直觀地分辨出軸承的狀態。通過對大量實驗數據的分析,發現軸承4種狀態的SDP圖像具有很好的重復性。

圖5 軸承4種狀態的SDP圖像Fig.5 SDP image of bearing with four kinds of fault

由于SDP圖像是通過對稱鏡像得到,取其中一部分花瓣進行分析就能夠表達整個圖像的特點。因此,選擇第一象限兩個花瓣進行分析。首先,只繪制第一象限兩個花瓣的極坐標圖并將其直接轉換為二值圖像,如圖6所示。然后,通過Matlab中regionprops函數實現對軸承4種狀態第一象限花瓣形狀特征的提取,能夠得到4個20×11的特征矩陣。因為二值圖像是由第一象限花瓣極坐標圖直接轉換得到,極坐標中心與二值圖像左下角的相對位置不會變化,而regionprops函數默認二值圖像左下角為0像素點,使得形狀特征的計算基準是一致的,這就保證了軸承不同故障狀態形狀特征的可比性。表1給出了軸承4種狀態下2個樣本信號的11維形狀特征。

圖6 第一象限花瓣二值圖像Fig.6 Binary image of petals in first quadrant

在得到形狀特征后,通過LLTSA將11維的特征向量進行維數約簡。設置鄰域參數K=12,目標維數d=3。作為比較,采用同樣參數設置的LTSA(局部切空間排列)和LPP對特征集進行降維處理。圖7為3種降維方法低維特征提取效果圖。從圖7中可以看出,LLTSA提取的低維特征能夠將軸承的4種狀態完全區分開,并且低維特征具有較好的聚類效果;LTSA雖然基本上能將4種狀態區分開,但仍然有少部分滾動體故障樣本混入到外圈故障和內圈故障樣本中,并且不同故障類型的距離也較近;LPP能夠將正常和外圈故障進行分離,并且各狀態的聚類效果較好,但也有部分內圈故障樣本和部分滾動體故障樣本出現了的混疊。以上分析表明LLTSA的性能要優于LTSA和LPP。

表1 軸承4種狀態下的形狀特征

圖7 3種流行學習算法特征約簡結果對比Fig.7 Comparison of dimension reduction effects of three manifold algorithms

為了進一步說明本文方法的有效性,采用SVM對提取低維特征進行分類識別。作為比較,同樣采用SVM對未經約簡的原始特征、經LTSA約簡的低維特征以及經LPP約簡的低維特征進行分類識別。實驗將各狀態20組樣本隨機分成兩組,10組作為訓練樣本,剩余10組作為測試樣本。試驗中SVM的核函數選用效果較優的徑向基核函數(RBF),懲罰參數C=10,核參數g=1,軸承正常、內圈故障、滾動體故障以及外圈故障的訓練標簽依次為1、2、3和4。表2為各方法故障特征的識別結果。

表2 支持向量機識別結果/%

從表2可以看出,根據SDP圖像提取的形狀特征雖然能夠較為準確的診斷出軸承的工作狀態,但識別率不是很高,這說明提取的形狀特征中含有一些重疊的故障信息,降低了故障分類的準確率,從而證明需要進行維數約簡以提取低維特征。而LTSA-SVM的識別率相比于不進行維數約簡的識別率有了提高,但和圖7(b)的分析結果一樣,除軸承正常狀態被完全正確識別外,其余三種狀態都有錯誤識別。在LPP-SVM中,有2個滾動體故障被錯分為內圈故障,識別率得到了一定程度的提高。在LLTSA-SVM中,軸承各狀態的識別率均達到了100%,效果優于LPP-SVM和LTSA-SVM方法,這說明本文所提方法是有效可行的。

5結論

本文在SDP方法、圖像處理和LLTSA的基礎上提出了基于圖像形狀特征和LLTSA的滾動軸承故障診斷方法。結合SDP對不同信號的可視化能力、形狀特征對圖像的表征能力以及LLTSA對非線性特征的提取能力等特點,對滾動軸承進行故障診斷。實驗結果表明,提取的SDP圖像形狀特征能夠反映軸承狀態的變化,而經LLTSA維數約簡后的低維特征聚類效果更好,提高了故障特征的分辨率,簡化了特征數據,增強了故障模式識別的準確率。將提取后的低維特征輸入SVM進行識別,對軸承4種狀態的識別率均達到了100%,驗證了本文方法的有效性。

參 考 文 獻

[1] Li Xu,Zheng A’nan,Zhang Xu-nan,et al. Rolling element bearing fault detection using support vector machine with improved ant colony optimization[J]. Measurement,2013,46:2726-2734.

[2] 唐貴基,鄧飛躍,何玉靈,等. 基于時間-小波能譜熵的滾動軸承故障診斷研究[J]. 振動與沖擊,2014,33(7):68-72,91.

TANG Gui-ji,DENG Fei-yue,HE Yu-ling,et al. Rolling element bearing fault diagnosis based on time-wavelat energy spectrum entropy[J]. Journal of Vibration and Shock,2014,33(7):68-72,91.

[3] 程利軍,張英堂,李志寧,等. 基于時頻分析及階比跟蹤的曲軸軸承故障診斷研究[J]. 振動與沖擊,2012,31(19):73-78.

CHENG Li-jun,ZHANG Ying-tang,LI Zhi-ning,et al.Fault diagnosis of an engine’s main bearing based on time-frequency analysis and order tracking[J]. Journal of Vibration and Shock,2012,31(19):73-78.

[4] 竇唯,劉占生. 旋轉機械故障診斷的圖形識別方法研究[J]. 振動與沖擊,2012,31(17):171-175.

DOU Wei,LIU Zhan-sheng. A fault diagnosis method based on graphic recognition for rotating machine[J]. Journal of Vibration and Shock,2012,31(17):171-175.

[5] 張立國,楊瑾,李晶,等. 基于小波包和數學形態學結合的圖像特征提取方法[J]. 儀器儀表學報,2010,31(10):2285-2290.

ZHANG Li-guo,YANG Jin,LI Jing,et al. Image characteristic extraction method based on wavelat packet andmathematicalm orphology[J]. Chinese Journal of Scientific Instrument,2010,31(10):2285-2290.

[6] Shibata K,Takahashi A,Shirai T. Fault diagnosis of rotating machinery through visualization of sound signals[J]. Mechanical Systems and Signal Processing,2000,14(2):229-241.

[7] Wu Jian-da,Chuang Chao-qin. Fault diagnosis of internal combustion engines using visual dot patterns of acoustic and vibration signals[J]. NDT Int.,2005,38(8):605-614.

[8] Delvecchio S,D’Elia G,Mucchi E,et al. Advanced signal processing tools for the vibratory surveillance of assembly faults in diesel engine cold tests[J]. ASME Journal of Vibration and Acoustics,2010,132:1-10.

[9] Kouropteva O,Okun O,Pietikainen M. Supervised locally linear embedding algorithm for pattern recognition[J]. Pattern Recognition and Image Analysis,2003,2652(9):386-394.

[10] He X F,Yan S C,Hu Y X,et al. Face recognition using laplacianfaces[J]. IEEE Trans.Pattern Analysis and Machine Intelligence,2005,27(3):328-340.

[11] Zhang T H,Yang J,Zhao D L,et al. Linear local tangent space alignment and application to face recognition[J]. Neurocomputing,2007,70:1547-1553.

[12] 楊成,馮燾,王中方,等. 基于SDP法診斷發動機的異響[J]. 聲學技術,2010,29(5):523-527.

YANG Cheng,FENG Tao,WANG Zhong-fang,et al. Abnormal noise diagnosis of motorcycle engines based on Symmetrized Dot Pattern method[J]. Technical Acoustics,2010,29(5):523-527.

[13] 任玲輝,劉凱,張海燕. 基于圖像處理技術的機械故障診斷研究進展[J]. 機械設計與研究,2011,27(5):21-24.

REN Ling-hui,LIU Kai,ZHANG Hai-yan. Progress on mechanical fault diagnosis based on image processing[J]. Machine Design and Research,2011,27(5):21-24.

[14] 秦襄培,鄭賢中. Matlab數字圖像處理寶典[M]. 北京:電子工業出版社,2011.

[15] 李鋒,湯寶平,陳法法. 基于線性局部切空間排列維數化簡的故障診斷[J]. 振動與沖擊,2012,31(13):36-40.

LI Feng,TANG Bao-ping,CHEN Fa-fa. Fault diagnosis model based on dimension reduction using linear local tangent space alignment[J]. Journal of Vibration and Shock,2012,31(13):36-40.

[16] 向丹,葛爽. 一種基于小波包樣本熵和流行學習的故障特征提取模型[J]. 振動與沖擊,2014,33(11):1-5.

XIANG Dan,GE Shuang. A model of fault feature extraction based on wavelet packet sample entropy and manifold learning[J]. Journal of Vibration and Shock,2014,33(11):1-5.

Fault diagnosis method based on image shape features and LLTSA

ZHANG Qian-tu1,2, FANG Li-qing2

(1. Military Representative Office in NO.356 Factory, Kunming 650114, China;2. Department of Artillery Engineering, Ordnance Engineering College, Shijiazhuang 050003, China)

Abstract:Aiming at fault diagnosis problems of rolling bearing, a fault diagnosis method based on image shape features and linear local tangent space alignment (LLTSA) was proposed. Firstly, a time waveform was transformed into a snowflake image in polar coordinate space with the symmetrized dot pattern (SDP) method. The image shape features were initially extracted on the basis of analyzing the characteristics of the image. Secondly, LLTSA was introduced to compress the initial high-dimensional features into low-dimensional ones with a better discrimination. Finally, a support vector machine (SVM) was employed to classify and evaluate low-dimensional features. The test results of rolling bearing fault diagnosis showed that the image shape features can characterize bearing states, and LLTSA can reduce the complexity of feature data and enhance their clustering effect; furthermore, a relatively higher identification rate of four bearing states reaches 100% with SVM, the effectiveness of the proposed method was verified.

Key words:symmetrized dot pattern (SDP); shape feature; linear local tangent space alignment (LLTSA); support vector machine (SVM); fault diagnosis

基金項目:軍內科研項目

收稿日期:2015-01-07修改稿收到日期:2015-05-03

通信作者房立清 男,博士后,博士生導師,1969年生

中圖分類號:TH165.3

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.09.028

第一作者 張前圖 男,碩士生,1991年生

猜你喜歡
故障診斷
基于包絡解調原理的低轉速滾動軸承故障診斷
一重技術(2021年5期)2022-01-18 05:42:10
ILWT-EEMD數據處理的ELM滾動軸承故障診斷
水泵技術(2021年3期)2021-08-14 02:09:20
凍干機常見故障診斷與維修
基于EWT-SVDP的旋轉機械故障診斷
數控機床電氣系統的故障診斷與維修
電子制作(2018年10期)2018-08-04 03:24:46
基于改進的G-SVS LMS 與冗余提升小波的滾動軸承故障診斷
因果圖定性分析法及其在故障診斷中的應用
改進的奇異值分解在軸承故障診斷中的應用
基于LCD和排列熵的滾動軸承故障診斷
基于KPCA和PSOSVM的異步電機故障診斷
主站蜘蛛池模板: 99爱视频精品免视看| 一本色道久久88| 亚洲欧洲AV一区二区三区| 99草精品视频| 2020最新国产精品视频| 毛片最新网址| 久久精品无码一区二区国产区| 欧美三级日韩三级| 国产一级毛片高清完整视频版| 国产剧情无码视频在线观看| 亚洲无码免费黄色网址| 国产网站一区二区三区| 国产女人18水真多毛片18精品| 熟女日韩精品2区| 在线观看av永久| 国产精品视频猛进猛出| 欧美五月婷婷| 国产va在线| 国产成人精品一区二区不卡 | 中国一级毛片免费观看| 亚洲综合色婷婷| 免费一级毛片| 91视频青青草| 毛片免费在线视频| 尤物成AV人片在线观看| 成年人视频一区二区| 欧亚日韩Av| 91福利免费视频| 91小视频版在线观看www| 99精品影院| 伊人狠狠丁香婷婷综合色| 都市激情亚洲综合久久| 欧美日韩国产在线人成app| 久久永久精品免费视频| 中文无码毛片又爽又刺激| 久久女人网| 国产无人区一区二区三区| 成人蜜桃网| 日韩国产精品无码一区二区三区 | 国内精品久久人妻无码大片高| 无码一区18禁| 99热这里只有精品2| 四虎永久免费在线| 成人福利在线看| 欧美在线中文字幕| 青青极品在线| Jizz国产色系免费| 免费在线一区| 美女毛片在线| 国产一区二区网站| 97狠狠操| 亚洲精品第五页| 欧美性天天| av在线手机播放| 91色爱欧美精品www| 久久久精品国产SM调教网站| 免费国产高清精品一区在线| 99尹人香蕉国产免费天天拍| 亚洲黄色网站视频| 色综合日本| 国产成人91精品| 亚洲黄色片免费看| 午夜一级做a爰片久久毛片| av一区二区三区在线观看| a在线亚洲男人的天堂试看| 国产在线视频导航| 免费AV在线播放观看18禁强制| 国产国语一级毛片| 黄色在线网| 区国产精品搜索视频| 国产一区二区三区在线观看视频| 日韩精品一区二区三区免费在线观看| 男女男免费视频网站国产| 日韩欧美亚洲国产成人综合| 久久国产精品波多野结衣| 黄色片中文字幕| 日韩一区二区三免费高清| 国外欧美一区另类中文字幕| 伊人久久大香线蕉综合影视| 久久精品无码中文字幕| 亚洲欧美一区二区三区图片| 国产成人无码综合亚洲日韩不卡|