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

基于短時傅里葉變換和卷積神經網絡的軸承故障診斷方法

2018-10-20 02:12:52秦仙蓉孫遠韜
振動與沖擊 2018年19期
關鍵詞:故障診斷故障信號

李 恒,張 氫,秦仙蓉,孫遠韜

(同濟大學 機械與能源工程學院,上海 201804)

滾動軸承是大多數機械傳動裝置的關鍵性零件,同時又極易損壞,若發生故障,甚至對整個機械設備造成嚴重破壞。由于從機械設備現場采集到的軸承振動信號,往往具有非平穩的特點,且噪聲干擾嚴重,給故障模式識別帶來了巨大的困難。近年來,隨著監測的設備群規模變大,每個裝備需要的測點增多,每個測點的采樣頻率高,從開始服役到壽命終止的數據收集歷時長,致使機械設備故障診斷領域也進入了“大數據”時代[1]。如何利用“大數據”來提高強噪聲、非平穩信號的故障診斷精度,將成為今后的研究熱點。

目前,針對強噪聲和非平穩信號的軸承故障診斷,主要集中于特征提取方法的研究,希望找到能反映軸承故障模式,并對工況和噪聲不敏感的特征。張小龍等[2]提出通過ITD-形態濾波濾除高頻噪聲信號后,用Teager能量譜提取有效故障特征的方法;錢林等[3]提出了變分模態分解(VMD)和自適應形態學濾波結合的方法,對故障信號進行降噪處理后提取故障特征頻率;Unal等[4]利用Hibert-Huang變換分析滾動軸承振動信號,并通過傅里葉變換(FFT)提取故障特征。上述方法存在兩個主要問題:①人工提取對噪聲和工況不敏感的故障特征較為困難,且需要一定的先驗知識;②將故障診斷分為信號預處理、特征提取和故障分類三個孤立的階段,破壞了各個階段的耦合關系,并造成了一部分故障信息的損失。因此,將預處理、特征提取和特征分類統一到一個框架下的端到端的故障診斷方式逐漸成為研究熱點。目前,部分學者提出基于深度學習的端到端的故障診斷方法,不僅能夠自適應的提取故障特征,而且能直接處理原始的故障信號。其中,Tamilselvan等[5]提出了一種基于深度置信網絡(DBN)的故障診斷方法,驗證了深度學習方法在故障模式識別方面的優勢及可行性。

卷積神經網絡(Convolution Neural Network,CNN)是一種有監督的深度學習算法,被廣泛的應用于模式識別領域。魏東等[6]提出了基于CNN的輸電線路內外故障模式的識別方法。CNN不同于其他深度學習算法的三大突出特點:局部感受野、權值共享以及池化。不僅降低了網絡的復雜度,也減少了過擬合的風險,使得構建處理海量數據的深度學習框架成為可能。

Lu等[7]將CNN應用到滾動軸承的故障診斷中去,并驗證了該方法在噪聲環境下了魯棒性。雖然在診斷結果上有一定的提高,但是存在著兩個主要問題:①以時間序列重構的二維矩陣作為卷積神經網絡的輸入,丟失了故障信號的頻率信息;②樣本量少且CNN模型簡單,沒有體現CNN在處理大數據方面的優勢。因此,本文提出了一種短時傅里葉變換(Short-time Fourier Transform,STFT)和CNN相結合的故障診斷方法,同時考慮了時域和頻域信息,并利用CNN對故障特征進行的自適應提取,實現端到端的故障診斷。最后實驗驗證了CNN在處理“大數據”方面的優勢,以及針對強噪聲非平穩信號的魯棒性。

1 短時傅里葉變換

短時傅里葉變換是針對時變、非平穩信號的一種聯合時頻分析方法。能將一維的故障振動信號變換成適應于CNN處理的二維矩陣:一種包含時域和頻域信息的特征譜。

1.1 短時傅里葉變換的定義

STFT的基本思想是采用固定長度的窗函數對時域信號進行截取,并對截取得到的信號進行傅里葉變換,得到時刻t附近很小時間段上的局部頻譜。通過窗函數在整個時間軸上的平移,最終變換得到每一時間段上局部頻譜的集合,因此,STFT是關于時間和頻率的二維函數。基本運算公式如下

(1)

式中:f(t)為時域信號;g(t-τ)為中心位于τ時刻的時間窗口。由此可以看出,STFT就是將信號f(t)乘以一個以τ為中心的窗函數g(t-τ)所做的傅里葉變換。

1.2 參數選擇

對于STFT,影響其變換效果的因素有兩個:窗函數類型和窗的寬度。選取好的窗函數能有效減少頻譜泄露和譜間干擾[8],而窗的寬度則影響著時域和頻域的相對分辨率:窗寬則頻率分辨率高;窗窄則時間分辨率高。因而需要根據特定信號的特征來合理選擇窗函數寬度。時域和頻域的分辨率可由下列公式計算得到

(2)

(3)

其中T和F分別表示時頻譜圖的時間和頻率分辨率,[·]表示向下取整運算。Nx是參與STFT的樣本長度,Nw是窗口寬度,No是窗口重疊寬度,Nf是參與傅里葉變換的點數,一般取Nf等于Nx。根據圖像分辨率的定義,將時頻譜的分辨率P定義為

P=T×F

(4)

P值大小決定了時頻譜所包含的信息量,這些信息既包括有用的故障信息,也包含噪聲干擾。因此合理的T和F,能夠使故障信號顯著,噪聲干擾減少。

在軸承故障診斷實踐中,故障特征頻率是一種重要的參考指標。針對滾動軸承的振動信號,可以以故障特征頻率的一倍頻為依據選擇樣本長度,以二倍頻為依據選擇STFT的窗口寬度。滾動軸承的故障頻率可以由下列公式計算得到

(5)

(6)

(7)

式中:fBPI、fBPO和fBS分別為滾動體過內圈頻率、滾動體過外圈頻率和滾動體自轉頻率;fH=n/60為內圈旋轉頻率;n為轉速;d和D分別為滾動體和節圓直徑;Z為滾子數目;α為接觸角度。

2 卷積神經網絡

卷積神經神經網絡CNN是一種前饋神經網絡,近年來在圖像識別和目標檢測方面有著出色的表現,能很好的處理過擬合問題,使更大規模的深度學習得以實現[9]。自從Lecun等[10]提出LeNet-5之后,卷積神經網絡的基本結構就被確定下來,主要由卷積層、池化層和全連接層組成,如圖1所示。

圖1 卷積神經網絡結構Fig.1 Structure of convolution neural network

在本文中CNN的輸入為時頻譜,經過卷積層、池化層和全連接層后,輸出為故障類型的類標。CNN的學習過程主要分為前向傳播過程和反向參數更新過程。

2.1 前向傳播階段

卷積神經網絡的前向傳播過程主要為三個部分:卷積層、池化層和全連接層。

2.1.1 卷積層

在卷積層上,用多個卷積核與輸入圖像進行卷積,加上偏置后通過一個激活函數就可以得到一系列特征圖[11],卷積過程數學表達式為

(8)

f(x)=max(0,lg(1+ex))

(9)

2.1.2 池化層

池化層常接在卷積層后面,對特征圖進行降維,同時在一定程度上保持特征尺度的不變性。常用的池化方法有:最大值池化(Max Pooling)、平均值池化(Mean Pooling)、隨機池化(Stochastic Pooling)等。池化層一般只進行降維操作,沒有參數,不需要進行權值更新。

在池化層上,對卷積層輸出的特征圖在每個不重疊的大小為n×n區域進行池化操作,選取每個區域上的最大值或者平均值,最終使輸出圖像在兩個維度上都縮小了的n倍。

2.1.3 全連接層

輸入圖像經過多個卷積層和池化層的交替傳播之后,依靠全連接層網絡針對提取的特征進行分類。在全連接層上,輸入是所有特征圖展開的一維特征向量經加權求和并通過激活函數后得到的

yk=f(wkxk-1+bk)

(10)

式中:k為網絡層的序號;yk為全連接層的輸出;xk-1是展開的一維特征向量;wk為權重系數;bk為偏置項。激活函數f(·)常使用Softmax,是一種適用于分類任務的激活函數。

2.2 反向參數更新

對于一個具體的分類任務,CNN的訓練目標是最小化網絡的損失函數,因此選擇一個合適的損失函數十分重要。常見的損失函數有均方誤差函數、交叉熵函數、負對數似然函數等[13],本文選用效果較好的交叉熵作為損失函數,表達式如下

(11)

式中:n為該類故障的樣本數;t為預測值;y為真實值。在訓練過程中,最小化損失函數的方法是梯度下降法。通過對式(11)求一階偏導數,可以逐層更新卷積神經網絡的可學習參數(w和b)

(12)

(13)

式中:w′和b′為更新后的權重和偏置;w和b為現有的權重和偏置;η為學習速率參數,用來控制權值更新的步長。η太大,會使網絡陷入局部最優;η過小,則會增加網絡的訓練時間。

2.3 基于CNN的故障診斷流程

基于CNN的故障診斷方法能夠將信號預處理、故障特征提取和故障模式分類融合到一起,實現故障特征自適應提取及智能診斷的具體流程,如圖2所示。采集的振動信號經STFT后,被分成訓練集和測試集。首先將訓練集輸入到CNN中進行參數學習,并用梯度下降法不斷更新權重(w)和偏置(b)。然后將訓練好的參數應用于測試集,得到故障診斷結果。

3 試驗驗證

本節將針對滾動軸承實測的振動信號,驗證所提出方法的可行性和有效性,并討論該方法在變轉速工況和強噪聲環境下的魯棒性。

3.1 數據集描述

數據集1 美國凱斯西儲大學所公開的軸承試驗數據。試驗用電動機轉速為1 792 r/min,采樣頻率為12 kHz。本文將驅動端軸承座的振動信號分為10種故障模式,依次為正常狀態以及內圈、外圈和滾動體分別在損傷斑點為0.18 mm、0.36 mm、0.54 mm三種直徑下的故障模式。每類故障的振動信號被分為400個樣本,共建立了一個具有4 000個樣本的故障數據集。對每種故障類標采用“獨熱(one-hot)”的編碼方式,即設置一個向量,其維數與故障類別數相同,除了某一位數字是1以外其余各維數字都是0。例如將正常軸承信號的類標設置為0,則可以編碼成([1,0,0,0,0,0,0,0,0,0])。

圖2 基于CNN的故障診斷流程圖Fig.2 Flowchart of the CNN-based method for fault diagnosis

滾動軸承型號為SKF6205,參數見表1。由式(5)~式(7)計算出故障特征頻率,如表2所示。根據一倍頻和二倍頻所對應的特征周期確定樣本長度為0.085 s(1 024個采樣點),STFT的窗口寬度為0.01 s(128個采樣點)。

表1 滾動軸承SKF6205參數Tab.1 Parameters of the SKF6205

表2 SKF6205故障特征頻率Tab.2 Fault frequencies of the SKF6205

STFT采用64點和114點兩種不同的重疊點數,得到兩種不同分辨率的時頻譜。由式(2)~式(4)可以計算得到分辨率分別為65×15和65×65,將這兩種時頻譜都輸入到CNN進行故障識別,以正確率為指標選擇其一作為輸入。圖3為內圈和外圈故障情況下,不同分辨率的時頻譜。

數據集2 由QPZZ-II旋轉機械故障試驗平臺(見圖4)采集到的軸承實驗數據。實驗使用東華測試的便攜式數據采集分析系統(DH5901)采集軸承座處的振動信號,采樣頻率為12.8 kHz。分別測試軸承在轉速600 r/min、900 r/min、1 200 r/min的運行狀態,并設置一個600~1 200 r/min的升速過程模擬變轉速非平穩運行工況。實驗用軸承型號為N205E和NU205E,分別為外圈可拆卸式和內圈可拆卸式,滾子數目為13個,其他參數與SKF6025相同,如表1所示。計算轉速900 r/min情況下的特征頻率,如表3所示。試驗前用線切割加工的方式,分別在內圈、外圈和滾動體上加工出不同損傷程度的凹槽,如圖5所示。

圖3 不同分辨率下的故障時頻譜Fig.3 Time-frequency spectrums of faults at different resolutions

圖4 QPZZ-II旋轉機械故障試驗平臺Fig.4 QPZZ-II rotating machinery test rig

表3 N205E/NU205E故障特征頻率Tab.3 Fault frequencies of the N205E/NU205E

圖5 滾動軸承內圈和滾動體的故障形式Fig.5 The faults of inner ring and rolling element

試驗安排了9種模式的故障:正常、內圈和外圈輕-中-重三種程度的故障、滾動體輕度和重度故障。在每類故障中選取400個樣本,建立一個具有3 600個樣本的數據集。同樣可以根據故障特征頻率確定樣本長度為1 024個采樣點,STFT的窗口寬度為128的采樣點。

3.2 卷積神經網絡模型建立

目前,對如何確定CNN的結構尚沒有明確的方法,很大程度上取決于經驗。卷積層中卷積核的大小和數目、池化層的大小及層數、全連接層的神經元數目等參數的選擇都影響著CNN的識別效果。本文在LetNet-5基礎上,改變一些網絡參數,從中尋找適應于滾動軸承故障診斷的CNN模型。其參數如表4所示,其中C表示卷積層,P表示池化層,F表示全連接層。

表4 LetNet-5模型的參數Tab.4 The parameters of LetNet-5

將卷積層C1和C2的卷積層數目、卷積核大小(C1和C2相同)以及全連接層F的神經元數目作為變量,并定義取值范圍為

{′C1′:2,4,8,16,32,64;

′C2′:2,4,8,16,32,64;

′F′:128,256,512,1024,2 048;

′卷積核大小′:2×2,4×4,5×5,6×6;}

這些參數可以組合成576種CNN模型,以遍歷的方式,同時考慮識別正確率和訓練時間兩個指標,在這些組合中確定最優的網絡結構。在實驗過程中,發現去掉LeNet-5模型中的C3層,并在輸入層和卷積層C1之間增加一個池化層P0之后,能在取得高故障識別正確率的同時大大減少訓練時間。結果如表5所示,實驗平臺為windows10+Python,CPU采用英特爾Core i5-6500。

表5 最優網絡結構的識別結果Tab.5 Fault identification results of optimal network structure

取正確率最高的情況,確定輸入為65×65的時頻譜,最優網絡結構,如表6所示。為一個8層的CNN模型“In-P0-C1-P1-C2-P2-F-Out”,卷積核大小為4×4。可將時頻譜分辨率修改為64×64(消去一行一列),都圓整成2的倍數形式,能減少一定訓練時間,且對最后的正確率影響不大。

表6 卷積神經網絡的結構Tab.6 The structure of convolution neural network

確定輸入和網絡結構之后,本文分別使用不同的輸入數據格式和網絡結構,如表7所示。設計了兩組對比試驗:①驗證輸入為二維圖譜的方法比輸入為一維時間序列具有更好的故障識別效果。對比實驗采用“In-H1-H2-F-Out”的5層BP神經網絡結構直接對一維時間序列進行訓練,這里的H1和H2代表隱層。時間序列樣本長度取1 024個采樣點;②驗證以STFT后的時頻譜作為CNN的輸入比Lu等直接對時間序列重構得到的時間圖作為輸入具有更高的故障識別正確率。時間圖是將時間序列按行重新排列得到的二維矩陣,本文將1 024個采樣點重構成32×32的時間圖作為輸入。其使用的7層CNN結構為:“In-C1-P1-C2-P2-F-Out”。

表7 輸入數據格式設置Tab.7 The formats setting of the input data

對于“3.1”節所建立的數據集,同樣采用遍歷的方法,得到表7中所設置的三種網絡結構,能很好的兼顧計算時間以及訓練集的收斂速度。

在數據集1中,有10類故障,每類故障有400個樣本,總共4 000個樣本,分成3 200個樣本的訓練集和800個樣本的測試集。同時在數據集2中選取900 r/min轉速下的9類故障信號組成具有3 600個樣本的數據集,并分成3 000個樣本的訓練集和600個樣本的測試集。訓練過程采用“抓包”的形式,即每步從訓練集中抓取一個子集輸入網絡做進行訓練,并在每次訓練之后將得到的參數應用于測試集,將識別出的故障類標與正確類標比較,并按下式計算每個子集的正確率

(14)

式中:ui為子集中第i個樣本的識別結果,若識別正確,ui等于1,識別錯誤,ui等于0;N為子集樣本數,本文中為50。

圖6表示采用兩個數據集進行軸承故障診斷模擬實驗的結果。每幅圖的橫坐標代表訓練步數,縱坐標代表故障識別正確率。由圖6(a)和圖6(d)可以輸入為時間序列雖然在訓練集上都達到了100%的故障識別正確率,但測試集的正確率僅為60%和30%。圖6(b)和圖6(e)所示輸入為時間圖,故障識別正確率有了一定的提高,達到了90%和60%,且訓練集和測試集的差距減少。如圖6(c)和圖6(f)所示,將時間圖替換成STFT后的時頻譜,則不僅在訓練集和測試集上都取得了非常高的識別精度,而且收斂速度也得到了明顯的提升。

輸入為時間序列

(a)

(b)

(c)

(d)

(e)

(f)

表8給出了三種輸入格式下,CNN所需的訓練時間以及對單個樣本進行識別的時間。在實際的故障診斷中,由于CNN的網絡參數是事先訓練好的,所以故障診斷的時間效率應該用識別時間作為評價標準。從表8可知,雖然三種輸入形式在訓練時間上存在較大差距,但是在對單個樣本的識別時間上,差距不大,都在1 s左右。

表8 三種形式輸入的時間效率Tab.8 The time efficiency of the three forms of input

從實驗結果中可以看出,輸入為時間序列和時間圖時,在數據集1上取得的效果要優于數據集2,這是由于在QPZZ-II旋轉機械故障試驗平臺進行實驗時,軸承裝配間隙比較大,噪聲干擾嚴重,振動信號呈非平穩性。而使用STFT與CNN結合的識別時頻圖的方法在兩個數據集上都取得了很高的故障識別正確率,說明此方法具有一定的魯棒性,能夠處理非平穩和強噪聲干擾的信號。

3.3 方法魯棒性驗證

將時頻譜輸入CNN經訓練后得到的網絡參數,直接用于魯棒性驗證試驗。在數據集2中,選取變轉速運行時的振動信號樣本一共1 000個樣本,分成10組,進行故障識別,并與“3.2”節中恒定轉速900 r/min下的故障識別正確率進行對比,結果如圖7所示。

由圖7可知,對于恒定轉速900 r/min的軸承振動信號數據,“3.2”節中訓練得到的卷積神經網絡參數能很好的識別出故障類別,但在變轉速情況下的故障識別正確率卻差了很多,這是因為所使用的網絡參數是在900 r/min的恒定轉速數據集上訓練得到的,CNN提取的特征不能很好的排除轉速變化的影響。因此,希望能通過構造一個更加龐大的數據集,包含各種轉速下的故障信息,并輸入到CNN中進行學習,得到對變轉速工況不敏感的故障特征。試驗采用包括600 r/min,900 r/min和1 200 r/min轉速下的所有的故障信號,建立了具有10 800個樣本的數據集。經訓練得到網絡參數,應用在測試集上。得出的識別結果,如圖8所示。

圖7 變轉速與恒定轉速識別結果對比Fig.7 Results of fault recognition between variable speed and constant speed

圖8 大數據集下變轉速與恒定轉速識別結果對比Fig.8 Comparison of variable speed and constant speed in bigger dataset

使用一個包含各種轉速的故障數據集進行訓練得到的網絡參數,能很好的識別出變轉速和恒定轉速情況下的故障模式。說明CNN是一種基于數據驅動的故障診斷方法,數據量越大,故障模式越完整,診斷效果越好。同樣,在解決噪聲干擾時,也可以通過增加數據量的方式提高模型的魯棒性。圖9為變轉速工況下正常和外圈重度故障的振動數據分別加上10 dB和20 dB的高斯噪聲后的時域圖像。使用與圖8一樣的網絡參數進行故障識別。得到的結果,如圖10所示。

從圖10可知,在識別受噪聲干擾的故障信號時,正確率明顯下降。同樣使用增加數據的方法,構造一個包含噪聲干擾在內的數據集,輸入到CNN中進行訓練,提取出對噪聲不敏感的故障特征,進行故障識別,結果如圖11所示。

表9總結了不同數據格式在各種環境下的故障識別正確率。可以看出,如果采集到的軸承振動信號能包含各種工況、各種環境下的各類故障模式,那么基于STFT和CNN的故障診斷方法,能夠取得非常好的故障識別率。這就建立了魯棒性和數據量之間的關系,說明此種方法是一種適應于處理“大數據”的故障診斷方法。

正常 外圈重度故障

圖9 不同噪聲下的振動信號時域圖像Fig.9 Time domain images of vibration signal with different noises

圖10 不同噪聲下的故障識別結果Fig.10 Results of fault recognition under noises

圖11 帶噪聲數據集故障識別結果Fig.11 Results of fault recognition under noises with bigger dataset

4 結 論

本文提出了一種基于STFT和CNN的滾動軸承故障診斷方法。將故障信號的預處理,特征提取和分類結合在一起,一定程度上實現了端到端的軸承故障診斷。首先根據滾動軸承故障特征頻率選擇合適的樣本長度和STFT的窗口寬度,然后用遍歷的方法,在LetNet-5模型的基礎上選擇合適的網絡參數,構建適應于軸承故障診斷的CNN模型。通過滾動軸承故障模擬試驗驗證了:

表9 不同數據格式在各種環境下的故障識別正確率Tab.9 Fault recognition accuracy of under various situations %

(1)將故障信號通過STFT后得到時頻譜能夠很好的區分出不同故障模式,作為CNN的輸入在訓練集上得到網絡參數,能夠在測試集上取得非常好的效果。

(2)CNN的魯棒性受到數據量的影響,針對強噪聲和非平穩信號,如果增加對應環境下的數據進行訓練,能夠顯著提高算法的故障識別正確率。本文的研究為如何利用機械設備監測“大數據”來提高故障診斷效果提供了新的思路,即通過構建能夠處理大量數據的深度學習框架,提高強噪聲干擾,非平穩信號的故障識別精度。

猜你喜歡
故障診斷故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 国产成人亚洲精品无码电影| 福利在线一区| 无码国产伊人| 18禁高潮出水呻吟娇喘蜜芽| 国产激情无码一区二区APP| 国产在线观看精品| 狠狠综合久久| 香蕉久久国产超碰青草| 亚洲欧洲日产国产无码AV| 午夜不卡视频| 欧美在线一二区| 色综合婷婷| 粉嫩国产白浆在线观看| 国产精品青青| 久久国语对白| 国产欧美日本在线观看| 欧美午夜理伦三级在线观看| 欧美黄网站免费观看| 国产真实乱人视频| 国产96在线 | 久久人搡人人玩人妻精品一| 久久一本日韩精品中文字幕屁孩| 亚洲精品另类| 理论片一区| 日韩欧美视频第一区在线观看| 国产精品xxx| 精品国产女同疯狂摩擦2| 91久久国产成人免费观看| 九九精品在线观看| 91在线视频福利| 亚洲最大看欧美片网站地址| 欧美一级高清免费a| 青青操国产| 欧美国产精品拍自| 高清欧美性猛交XXXX黑人猛交 | 成·人免费午夜无码视频在线观看| 国产精品99久久久| 丁香五月婷婷激情基地| 操操操综合网| 国产粉嫩粉嫩的18在线播放91 | 天天综合网亚洲网站| 亚洲国产欧美国产综合久久| 日本道综合一本久久久88| 好紧好深好大乳无码中文字幕| 国产在线观看91精品亚瑟| 91九色国产porny| 岛国精品一区免费视频在线观看| 99re精彩视频| 黄色网址免费在线| 3344在线观看无码| 久久精品只有这里有| 夜夜高潮夜夜爽国产伦精品| 国产成人综合网| 在线综合亚洲欧美网站| 伊人无码视屏| 99热这里只有精品5| 91精品国产91欠久久久久| 国产精品99r8在线观看| 国产黄色爱视频| 中文字幕在线观看日本| 国产主播喷水| 99视频在线免费| 色一情一乱一伦一区二区三区小说| 日韩国产高清无码| 国产xxxxx免费视频| 六月婷婷综合| 国产丝袜丝视频在线观看| 精品三级在线| 四虎永久在线精品国产免费| 欧美成人A视频| 91精品国产麻豆国产自产在线| 国产女人18水真多毛片18精品| 国产精品女人呻吟在线观看| 亚洲无码高清一区二区| 一级毛片视频免费| 国产欧美日韩免费| 看av免费毛片手机播放| 亚洲欧洲日本在线| 午夜国产理论| 成年女人a毛片免费视频| 好吊色妇女免费视频免费| 国产人人射|