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

基于Laplace小波卷積和BiGRU的少量樣本故障診斷方法

2023-01-03 04:36:48路顏萍
振動與沖擊 2022年24期
關鍵詞:故障診斷特征信號

羅 浩, 何 超, 陳 彪, 路顏萍, 張 欣, 張 利

(1. 遼寧大學 信息學院,沈陽 110036;2. 東北大學 材料科學與工程學院,沈陽 110819)

機械旋轉設備故障診斷已步入智能化、自動化時代。其中,滾動軸承扮演著重要角色,其健康狀態直接影響設備的穩定性、可靠性[1]。然而,滾動軸承往往在復雜條件下工作,受到材料退化、溫濕度等多種要素影響,故障后輕則影響工廠效益,重則導致人員傷亡。因此,對其健康狀態監測是極具研究意義的課題。

過去基于信號分析、群智能進化、機器學習的故障診斷方法不斷涌現[2-4]。然而相關算法依賴專家先驗知識,難以處理高維數據;群智能優化的尋優結果難以穩定且時間復雜度高。隨后,深度學習智能故障診斷如火如荼?;蚨嘟嵌确治鲂D機械振動信號特征,或改進深度優化算法,或設計新穎高效的網絡結構[5]。

近些年,少量樣本故障診斷成為研究熱點。有的利用模型特征提取優勢和正則化策略,例如Han等[6]提出基于雙向長短時記憶(bidirectional long short-term memory, BiLSTM)和膠囊網絡少量樣本故障診斷方法。振動信號經卷積神經網絡(convolutional neural network,CNN)和BiLSTM去噪融合后,膠囊網絡對少量樣本故障診斷有著良好性能。Saufi等[7]提出基于譜峰度濾波和粒子群優化堆疊稀疏自編碼器少量樣本診斷方法,當每個故障類型訓練數為100時,取得較高精度。有的依據數據分布生成多而高質量樣本,例如Li等[8]提出條件Wasserstein生成對抗網絡(conditional Wasserstein generative adversarial networks, CWGAN),利用源域大量標簽數據訓練CWGAN生成大量樣本。針對有限標簽目標域數據,加載預訓練參數并微調CWGAN實現少量樣本遷移,取得良好效果。有的應用元學習、遷移學習等新興機器學習技術實現診斷,例如Zhang等[9]提出基于孿生神經網絡的小樣本故障診斷方法,輸入相同或不同類的樣本對,計算兩者特征向量的L1距離,判斷是否同屬一類來訓練模型,最后將支持集與查詢集作為特征對并計算相似度來實現故障診斷。Wang等[10]在此基礎上,利用全連接來模擬特征對相似性度量,并添加正則化手段提升效果。Wu等[11]比較了基于特征遷移、微調、元關系網絡三者間小樣本遷移性能,得出當樣本較少或源域、目標域相似性較大時,元遷移占據主導地位;反之,特征遷移優勢逐漸明顯。

針對模型特征提取優勢和正則化策略,除上述膠囊網絡外,沈濤等[12]探究CNN-LSTM(convolutional neural network-long short-term memor)對少量樣本表達能力,選擇訓練集比例0.3來探討模型應對各種復雜工況能力。然而單尺度單向的CNN-LSTM,沒有充分挖掘模型性能,并且批歸一化和Dropout等都略常規。Yang等[13]提出一種基于1DCNN和BiGRU框架損害檢測方法,將多傳感器信號組成信號矩陣輸入模型中,融合兩者提取的特征實現故障診斷,但未探究少量訓練樣本下BiGRU效果。另外,雙路Laplace小波卷積和雙向門控循環單元(dual Laplace wavelet convolution bidirectional gated recurrent unit, DLWCB)每類訓練數是Saufi等研究的60%為60個,泛化難度更大。

在上述分析基礎上,研究內容總結如下:

(1) 從拓撲結構和正則化出發,提出一種少量樣本故障診斷思路,設計一種端對端融合雙路Laplace小波卷積核和BiGRU的故障診斷方法。首先提出Laplace小波卷積核,并和平均能量池化共同作用于所設計模型,提取多尺度特征,大卷積核具有強魯棒性[14]。然后,BiGRU學習隱藏信息進一步提取高級特征。設計全局平均池化(global average pooling, GAP)增強通道間聯系并提高BiGRU的特征利用率。

(2) 引入故障診斷不常見的手段標簽平滑正則化(label smoothing regularization, LSR)和流形正則化并結合批歸一化(batch normalization, BN)、PReLU動態激活、AdamP等來改善DLWCB泛化性,緩解少量樣本下模型的過擬合現象。

(3) 針對噪聲樣本,提出具有參數傳遞的故障診斷框架,經少量樣本微調模型便具備更高的效率。探究了BiGRU和GAP在少量樣本故障診斷應用中的性能。實驗表明,該方法具有較高診斷效率。

1 相關理論基礎

1.1 卷積神經網絡

卷積神經網絡(convolutional neural network,CNN)一般由濾波塊和分類診斷塊兩個模塊組成。一般故障診斷CNN,如圖1所示。

圖1 故障診斷CNN結構Fig.1 CNN for fault diagnosis

信號處理領域,一般應用具有相同內核的一維卷積計算信號的延遲累積。輸出y如式(1)所示

(1)

式中:kw和bw為卷積核和偏置;xt-w+1為輸入振動信號信號; *為卷積; PReLU(·)為神經網絡激活函數。

池化對特征進行選擇并降低模型參數量防止過擬合,故障診斷中常用最大池化。激活函數可以增強神經網絡的表示和學習能力,提高計算效率。BN不僅可以提高優化效率,而且由于其隨機選擇批次,可以增強模型的泛化能力。全連接將振動信號分布式特征表示映射到樣本標記空間,最后應用SoftMax激活函數進行故障診斷。

1.2 雙向門控循環神經網絡

門控循環神經網絡(gate recurrent unit,GRU)由更新門zt和重置門rt組成。BiGRU由雙向多個GRU組成,如圖2、圖3所示,信號通過不同且不共享參數的隱藏層將前后向輸出連接到相同層以提取過去和未來的特征。如式(2)~式(4)所示

(2)

(3)

(4)

式中:wt,vt為前向和后向的狀態權重矩陣;xt為輸入信號;bt為偏置。

然而,深度故障診斷中,一般取最后一個隱藏神經元細胞輸出作為BiGRU學習到的振動信號特征并作為后續輸入,忽略了其他GRU細胞,在DLWCB中,將在BiGRU后連接GAP解決。

圖2 GRU細胞Fig.2 GRU cell

圖3 DLWCB智能故障診斷模型Fig.3 DLWCB intelligent fault diagnosis model

2 Laplace小波卷積

2.1 Laplace小波卷積定義

受Morlet小波[15-16]啟發,將Laplace小波思想融入卷積核中,提出并定義Laplace小波卷積核。

在時域中,小波基本字典ψu,s(t)定義如式(5)

(5)

式中:ψ(·)為小波基函數;t為時間;s為尺度因子,u為平移因子;s,u為自適應可調節參數。

考慮到滾動軸承的機械振動信號屬于實信號,故采用實Laplace小波基函數分析信號,如式(6)所示

(6)

式中:f為信號頻率;ξ為黏性阻尼比;τ為時間參數;A為小波歸一化函數。

由式(5)、式(6)可得實Laplace小波卷積字典ψL,u,s(t)如式(7)所示

(7)

將式(7)代入式(1),得到Laplace小波卷積核輸出yL,如(8)式所示

(8)

2.2 Laplace小波卷積參數量

暫不考慮偏置,普通卷積參數量為卷積核尺寸與卷積核個數的乘積,而Laplace小波卷積核僅s,u需調整,參數量是卷積核個數的二倍。DLWCB中,第一層參數量為50×2=100,而普通卷積核參數量為50×18=900。顯然卷積核尺寸越大,參數減少越明顯。

2.3 參數s和u更新

s和u更新依據反向傳播算法,在第一層Laplace小波卷積中,參數更新可以描述為

(9)

(10)

同時,根據鏈式求導法則可以得到兩個參數的偏導數如式(11)和式(12)所示,并將式(11)、式(12)代入式(9)、式(10)中對兩個參數進行更新。

(11)

(12)

3 DLWCB故障診斷模型

CNN-RNN(convolutional neural network-recurrent neural network)已取得一定應用[17-18]。但少量樣本下其性能表現卻鮮有研究,而且優化算法和訓練方式較為常規,潛在性能未能得到進一步挖掘。為此,提出一種少量樣本智能故障診斷方法——DLWCB,具體結構見圖3。DLWCB由數據增強層、雙路Laplace小波卷積層、特征融合層、BiGRU、GAP和診斷層組成。

圖3結構中,將Laplace小波思想融入卷積過程中,并應用于首個卷積層。

GAP則解決了1.2節所提問題,通過計算各個GRU提取的特征矩陣作全局平均池化,綜合考慮各個GRU的特征輸出,降低了診斷層的訓練壓力并提高特征利用率。DLWCB相關結構細節參數,如表1所示。

表1 DLWCB結構細節Tab.1 Structure details of DLWCB圖層類型

4 實驗結果與分析

4.1 數據描述與實驗配置

本節將結合兩個軸承故障診斷案例,以證明所提方法性能、效率和應用價值。實驗應用pytorch 1.8.0,運行在GTX970M GPU上。采用平均功率池化,PReLU動態激活函數,AdamP反向傳播算法[19]和帶有流形正則項的標簽平滑損失函數[20]來訓練DLWCB。

數據集S1-凱斯西儲大學滾動軸承數據[21]。采樣頻率12 kHz,負載0~2.237 1 kW。電動煙花引起單點故障。加速度傳感器位于電機外殼驅動端,收集加速度數據。據負載不同,劃分ABCD 4組,如表2所示。

表2 數據集S1描述Tab.2 Data set description of S1

數據集S2-為驗證算法的有效性和可靠性,搭建的軸承故障實驗平臺,如圖4所示。該平臺由電機、齒輪箱、聯軸器、軸承座、操作柜和操作臺等部分組成。采集頻率為50 kHz,使用加速度計采集1 000 r/min下振動信號。采集單元為德國Bruel & Kjaer Vibro公司的一款用于振動分析的VDAU-6000,具有16個可以實現同步數據采集的輸入通道。根據不同位置的不同狀態將振動信號分為4類。如表3所示。

圖4 機械故障模擬實驗臺Fig.4 Mechanical fault simulation test bench

表3 數據集S2描述Tab.3 Data set description of S2

實驗中,S1,S2均采用滑動窗口采集樣本,窗口大小為400,樣本長度為1 024。兩個數據集學習率分別為0.001,0.000 4,最大迭代次數epoch=150, 丟棄率Dropout=0.2, 權重衰減系數為0.000 1。采取早停策略(patience=10)。訓練樣本比例是0.1(20)~0.5(100)。

4.2 少量樣本問題描述

少量樣本的概念很難統一準確定義,這與數據和診斷模型復雜度有關。通常,機器學習中“少量樣本”會導致模型過擬合。當訓練數據不足時,網絡學習到特征表示是有限的,只能很好地擬合訓練數據,導致訓練集精度高,測試集精度低。當訓練數據充足時,算法可以有效提取整個數據集特征,在訓練和測試集上都取得良好性能。因此,其在兩個案例中描述如下:訓練數據較少,以訓練一個足夠泛化的模型,可以在測試集上實現有效故障分類。

4.3 batch_size參數選擇

批處理(batch_size, b)會影響模型訓練效率和泛化性。針對數據集S1-B,訓練集比例為0.3(60),僅改變b,實驗結果如圖5所示。

當b≤100時,b越大,DLWCB收斂時需要的epoch越多,單個epoch較短的訓練時間并不能使得整體收斂速度提高。至于收斂損失,lb=32,80,100<0.56,相對其他情況收斂損失更低。三者驗證集準確率基本相同。另外,從時間來看,tb=32=75.71 s,然而tb=80,100>100 s,當b=32時,達到與b=80,100相似的準確率所需時間更少。綜上,b=32具有最佳診斷效率。

圖5 不同批次訓練結果比較Fig.5 Comparison of training results of different batch_size

4.4 Laplace小波卷積核性能

CNN的第一個卷積層影響著整個模型的性能[22]。因此將第一層普通卷積核替換為Laplace小波卷積核,與信號卷積達到小波分析的效果。這樣使得時域信號轉換到頻域,更有效地識別軸承狀態[23]。實驗中,f=100,ξ=0.03,τ=0.1,A=0.08,s∈[1,100],u∈[0,100],s,u向量維度為輸出通道數大小。

采用數據集S2,每類訓練樣本20~100,對比DLWCB和DCB(dual convolution bidirectional gated recurrent unit)(第一層未使用Laplace小波卷積)的性能表現,相關結果如表4所示。

表4 數據集S2性能表現Tab.4 Performance of S2 data set

隨著訓練樣本數增加,模型性能逐漸提升。由于采用正則化策略,當訓練樣本數大于60時,取得較高精度。即使這樣,DLWCB相較于DCB準確率更高、方差更小、更穩定。當樣本數少于40時,由于訓練樣本數太少導致正則化方法更難泛化網絡。此時,DLWCB在收斂損失上更低;當樣本數為20或40時,DLWCB達到87.56%和94.15%,相較于DCB提升17.33%和7.79%。但DLWCB的時間復雜度較高。

4.5 少量樣本參數選擇

在S1-B,b=32下,進行5次實驗取平均值,得到訓練時間和測試集準確率,如圖6所示。

圖6 S1不同訓練集比例下性能Fig.6 Performance of different training proportion in S1

結合表4和圖6,軸承數據集S2由于信號特征不夠明顯,達到與S1相似性能需更多訓練樣本。

針對數據集S2,據4.2節所述,同時考慮正則化的影響,采用訓練集比例0.3(60)建立少量樣本,因為此時訓練樣本數對模型的影響還沒達到最大,可以體現正則化對模型的作用;至于數據集S1,當訓練集比例大于0.3(60)時,隨著訓練樣本數增加,準確率提升不足1%且時間基本呈線性增加;當訓練集比例為0.3(60)時,DLWCB達到99%以上。權衡兩者,選擇0.3(60)為訓練集比例,建立少量樣本。

同樣,針對軸承數據集S1,在0.3(60)下,得到是否應用Laplace小波卷積的準確率,如表5所示。

表5 Laplace小波卷積Tab.5 LW convolution

少量樣本下,訓練樣本數對模型的影響較低,此時觀察Laplace小波卷積的作用具有較大參考價值。從表5可以看到相較于原始卷積,Laplace卷積最終收斂損失減少了約0.02,準確率提升2%。

綜上,針對少量樣本故障診斷,將信號處理領域中Laplace小波思想融入卷積核中也是一種應對思路。

4.6 正則化分析

采用數據集S2,如4.5節分析,將訓練集比例設置為0.3(60),建立少量樣本。驗證正則化方法的收斂性。

4.6.1 AdamP收斂性分析

幾種梯度下降算法在驗證集上的表現,如圖7所示。SGDM收斂損失最大且收斂速度最慢。Adam,AdamP有著較快的收斂速度。RMSprop損失收斂曲線波動相對較大。反觀AdamP則兼具平緩、較快收斂,在epoch=57時收斂于0.64,這使得DLWCB更具穩定性。

圖7 不同算法損失值對比圖Fig.7 Comparison of loss values of different algorithms

4.6.2 PReLU和ReLU泛化性分析

ReLU和PReLU兩種激活函數在不同epoch下的收斂曲線,如圖8所示。結果顯示了在每個訓練階段后訓練過程的訓練準確率和誤差??梢钥闯?,與ReLU相比,PReLU-DLWCB準確率更高并且收斂損失更低。由圖8(b)可知,PReLU相較于ReLU訓練過程更加平滑,在epoch=72時,收斂到最低損失。表明PReLU-DLWCB更能學習到信號特征且更具穩定性。

圖8 不同激活函數下DLWCB 性能對比Fig.8 Comparison of different activation functions

4.6.3 目標函數泛化性分析

針對少量樣本,僅改變損失函數,LSR與CrossEntropy, Focal Loss[24]和GHMC(gradient harmonizing mechanism classification)[25]損失收斂曲線,如圖9所示。當訓練集比例為0.3(60)時,CrossEntropy達到97.32%,GHMC達到97.32%,LSR達到98.04%。在數據集S2下,LSR在eopch=72時最早完成收斂,收斂速度更快,準確率曲線波動較小。

圖9 不同損失函數下DLWCB準確率對比圖Fig.9 Accuracy comparison of under different loss functions

顯然,少量訓練樣本下,LSR促使模型更快、更好地學習數據特征,具備良好的穩定性和訓練效率。

4.7 DLWCB消融實驗、通用性分析

為進一步分析DLWCB各部分作用,使用軸承數據集S1-B,在訓練樣本比例0.3(60)下與DCNN,BiGRU,BiLSTM進行性能對比。對比結果如表6所示。

表6 DLWCB消融實驗Tab.6 DLWCB ablation experiment

在對比實驗過程中,單層LSTM(long short-term memor)和GRU處理少量的高維原始故障樣本振動信號的效果不理想,故對原始振動信號采用主成分分析(principal component analysis, PCA)降維。圖6中,DLWCB精度達到99.86%。表6表明,GRU由于參數較少,相比于LSTM提升約0.17%,BiGRU相比于GRU提升約1.3%。

另外,BiGRU使得DLWCB學習到振動信號不同隱藏位置特征,準確率提升約8%。DCNN學習到少量樣本高低頻信號,提升約10%。兩者對于少量樣本都有不錯的性能表現。

同樣,為驗證DLWCB在不同負載下通用性,以數據集S1-B為訓練集。將其遷移至數據集S1下A, C,D進行故障診斷,得到混淆矩陣如圖10所示。

原始振動信號中,負載越低,信號中包含脈沖信息越少;負載越高,則信號包含噪聲越多。由圖10可知,DLWCB在不同負載下仍具有較高識別率,在A,C,D下準確率分別是99.62%,99.50%,98.35%,說明DLWCB具有較好域適應能力。同樣,觀察到在不同負載下,對于4,5,8部分樣本識別不敏感,還需進一步研究改善,但整體上對正常和故障樣本都得到明確區分,各類診斷準確率在98%以上。

圖10 不同負載下域適應實驗Fig.10 Domain adaptation under different loads

上述相關實驗和分析表明,DLWCB具有較高的故障識別率和定位表現,可以為軸承實際工作中的故障快速定位診斷和維護提供一定指導。

4.8 可視化分析

為進一步揭示所提方法的特征提取和識別過程,針對軸承數據集S2,使用T-SNE將DLWCB提取的特征降至二維。相關結構特征表示,如圖11所示,不同灰度描述不同狀態。圖11(a)~圖11(e)依次為輸入振動信號,大卷積核提取的低頻特征,小卷積核提取的高頻特征,BiGRU只保存最后一個GRU提取的特征,以及添加GAP后BiGRU保存各個GRU輸出特征。

可以發現,圖11(b)、圖11(c)對部分混亂樣本進行初步學習,雖然大卷積核特征分離不明顯,但可以有效應對噪聲;小卷積核卷積越深,特征分離越明顯。經特征融合后,相較于只關注最后一個隱藏GRU單元(圖11(d)),關注所有GRU單元使得DLWCB利用了更豐富的信息,提高了特征利用率,更有利于其從少量樣本中分離出各類故障(圖11(e))。這正是設計GAP的結果,反映出利用GAP針對少量樣本故障診斷具備先進性能。

圖11 不同結構輸出T-SNE可視化圖Fig.11 Visualization of T-SNE with different structures

綜上所述,DLWCB可以更好地將特征從不同類中分離出來,也暗示DLWCB具有較強泛化性能。

4.9 抗噪魯棒性分析

實際工廠信號多含有噪聲,為此,分析不同信噪比條件下,DLWCB抗噪能力。將高斯白噪聲與原始樣本混合形成新的含噪聲復合樣本,如式(13)所示。

(13)

式中:Psignal為原始信號功率;Pnosie為噪聲功率。

與以往研究直接將含噪信號輸入模型不同,DLWCB訓練和測試均使用圖12所示的故障診斷框架。此框架由離線預訓練和在線測試組成。離線框架用于訓練DLWCB所需預訓練參數,在線框架主要用于泛化性測試和快速訓練以應對噪聲樣本的故障診斷。

圖12 基于參數傳遞的DLWCB故障診斷框架Fig.12 DLWCB fault diagnosis framework based on parameter transfer

實驗將對數據集S1-B和S2下的預訓練參數進行微調。由于權重參數接近,這將減少DLWCB訓練時間,提高診斷效率,以便快速應對各種噪聲,實現故障定位。

將信噪比(signal-noise ratio,SNR)=-4~10 dB的高斯白噪聲添加到原始信號中,其他設置一致,加載訓練比例為0.3(60)下預訓練參數,使用相同比例含噪樣本微調DLWCB[26]。

針對軸承數據集S1-B,表7中,當SNR≥2 dB時,DLWCB取得較高的診斷精度,訓練時間在10~40 s。與圖7比較,診斷時間縮短約1/2。隨著SNR降低,信號中噪聲占比越來越高,對DLWCB的魯棒性提出更高的要求,診斷精度有所降低,但依然達到90%左右。一方面是由于Laplace小波卷積和大卷積核對噪聲強抗性,過濾部分噪聲并有效地提取信號特征;另一方面則是BiGRU提取少量噪聲樣本下不同位置隱藏信息,并通過GAP使得DLWCB關注到更多有用信息。

表7 S1下DLWCB不同SNR抗噪聲能力Tab.7 Anti-noise capability of different SNR of DLWCB in S1

與F方案相比,加載預訓練參數的E方案對含噪樣本診斷效率取得提升。信噪比越高,F所需診斷時間越長,但是E僅僅35 s就可以完成診斷;并且在各種SNR下,E都取得約1%~5%領先,可以看出在噪聲環境下,此方法可以取得較高的效率。

此外,針對數據集S2也做相同實驗,數據處理方法與S1相同,也是利用含噪聲樣本微調。以進一步證明所提方法的泛化性和可靠性。除學習率變為0.000 4外,其他設置相同,得到相關結果如表8所示。

表8 S2下DLWCB不同SNR抗噪聲能力Tab.8 Anti-noise capability of different SNR of DLWCB in S2

由表8可知,隨著信噪比增加,DLWCB的準確率逐漸提升,其在干擾場景下也具有相當高的性能表現,最高達到98.21%,表明了該方法具備良好的魯棒性。同樣可以看到,加載預訓練參數在時間和準確率上都取得提升。以上兩種數據集的實驗充分說明E方案具備較高診斷效率。

5 結 論

針對少量樣本的故障診斷問題,提出一種端到端融合Laplace小波卷積的DLWCB故障診斷模型,無需額外降噪算法,便可實現高效故障診斷。首先,提出Laplace小波卷積核并評估其性能。其次,通過LSR多重目標函數、引入AdamP、PReLU等在故障診斷領域不常用正則化進一步提高泛化性。最后,使用可視化手段初步理解DLWCB。實驗主要從準確率和時間復雜度等方面驗證該方法的可靠性,表明在少量樣本下,DLWCB取得優勢,具備一定應用價值。

然而,各個故障類別訓練樣本數量平衡,在實踐中,可能會遇到類別不平衡的數據集。針對少量不平衡數據集,在未來的工作中,考慮從生成對抗網絡、時間卷積網絡、注意力機制等拓撲結構或者元學習、小樣本學習、遷移學習、集成學習等機器學習技術方面入手擴展該方法的應用,探究它們的性能表現。

猜你喜歡
故障診斷特征信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
抓住特征巧觀察
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 新SSS无码手机在线观看| 国产成人91精品| 人妻熟妇日韩AV在线播放| 欧美一区二区人人喊爽| 国产亚洲欧美日韩在线一区| 极品av一区二区| 91成人在线观看视频| 国产黄色视频综合| 好紧好深好大乳无码中文字幕| 亚洲va欧美ⅴa国产va影院| 亚洲婷婷六月| 一本一本大道香蕉久在线播放| 88国产经典欧美一区二区三区| 亚洲av色吊丝无码| 亚洲男人天堂2018| 无码啪啪精品天堂浪潮av| 91精品国产福利| 欧美精品xx| 久久亚洲精少妇毛片午夜无码| 午夜少妇精品视频小电影| 久久精品免费国产大片| 日韩中文无码av超清| 成人福利一区二区视频在线| 日韩国产高清无码| 国产一二三区在线| 国产在线拍偷自揄观看视频网站| 亚洲男人的天堂网| 一级一级特黄女人精品毛片| 成人一级免费视频| 亚洲人成网18禁| 久久天天躁夜夜躁狠狠| 韩日无码在线不卡| 精品久久久久久中文字幕女| 国产91无码福利在线| 四虎影视无码永久免费观看| 99在线小视频| 2048国产精品原创综合在线| 国产乱子伦手机在线| 美臀人妻中出中文字幕在线| 青草视频网站在线观看| 亚洲综合专区| 97影院午夜在线观看视频| 国产精品手机视频一区二区| 亚洲精品视频免费看| 国产69精品久久久久妇女| 国内精品一区二区在线观看| 999精品在线视频| 国产午夜人做人免费视频中文 | 久久人人爽人人爽人人片aV东京热| 一本久道热中字伊人| 久久国产精品无码hdav| 亚洲熟女偷拍| 东京热av无码电影一区二区| 广东一级毛片| 欧美三级视频在线播放| 欧美激情福利| 久久久久88色偷偷| 一级高清毛片免费a级高清毛片| 国产99精品久久| 色综合久久88| 一级毛片高清| 人妻丝袜无码视频| 国产在线麻豆波多野结衣| 国产日韩欧美精品区性色| 欧美yw精品日本国产精品| 五月婷婷欧美| 国产菊爆视频在线观看| 日韩天堂视频| 色婷婷成人网| 亚洲欧美不卡| 成人福利在线视频免费观看| 一区二区午夜| 无码国产伊人| 91麻豆国产视频| 精品国产99久久| 人妻中文字幕无码久久一区| 国产精品偷伦在线观看| 国产成在线观看免费视频| 中文字幕久久精品波多野结| 亚洲免费黄色网| 97精品国产高清久久久久蜜芽| 国产精品99久久久久久董美香|