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

基于多通道時頻域信號的卷積神經(jīng)網(wǎng)絡(luò)智能故障診斷技術(shù)

2021-06-26 04:06:04孫仕鑫杜勁松
科學(xué)技術(shù)與工程 2021年15期
關(guān)鍵詞:故障診斷特征信號

孫仕鑫, 高 潔, 王 偉, 杜勁松, 楊 旭

(1.中國科學(xué)院沈陽自動化研究所, 沈陽 110016; 2.中國科學(xué)院機(jī)器人與智能制造創(chuàng)新研究院, 沈陽 110169; 3.中國科學(xué)院大學(xué), 北京 100049; 4.遼寧省智能檢測與裝備技術(shù)重點(diǎn)實(shí)驗(yàn)室, 沈陽 110179)

滾動軸承作為旋轉(zhuǎn)機(jī)械如齒輪箱的關(guān)鍵部件[1],長時間工作在嘈雜的環(huán)境中,不可避免地會出現(xiàn)各種故障。為防范可能存在的隱患,最常見的方法是通過分析其振動信號來進(jìn)行早期的故障診斷。傳統(tǒng)智能故障診斷算法主要包括特征設(shè)計(jì)、特征提取和選擇、狀態(tài)識別[2]。而特征設(shè)計(jì)需要專家知識,并且設(shè)計(jì)的特征具有專一性,即在此數(shù)據(jù)集上設(shè)計(jì)的特征不適用于另一數(shù)據(jù)集,因此傳統(tǒng)智能故障診斷算法具有特征設(shè)計(jì)復(fù)雜、不利于推廣等缺點(diǎn)。近年來,深度學(xué)習(xí)的崛起為克服上述缺點(diǎn)提供了新的方法[3]。深度學(xué)習(xí)是一種端到端的學(xué)習(xí)方法,不需要專家知識,原始數(shù)據(jù)直接作為輸入,分類結(jié)果作為輸出,省去了繁雜的特征設(shè)計(jì),且利于推廣,因此在滾動軸承故障診斷領(lǐng)域中廣泛使用。

BN是批標(biāo)準(zhǔn)化層圖1 卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.1 Convolutional neural network structure

然而,在實(shí)際工業(yè)環(huán)境中,算法不能學(xué)習(xí)到所有負(fù)載下的健康狀態(tài)特征,例如在極端情況下,算法只能在1 hp(馬力,horsepower)負(fù)載下訓(xùn)練,在實(shí)際測試時,卻需要在3 hp負(fù)載下有效識別健康狀態(tài),所以就要求算法有很強(qiáng)的負(fù)載域適應(yīng)能力。因此,中外專家學(xué)者相繼提出了在變負(fù)載下的滾動軸承故障診斷算法。Jiao等[4]提出了雙級對抗域自適應(yīng)網(wǎng)絡(luò)(double-level adversarial domain adaptation network,DL-ADAN),該模型由基于深度卷積網(wǎng)絡(luò)的特征提取器、域識別器和兩個標(biāo)簽分類器組成,用于識別跨領(lǐng)域的故障狀態(tài)。Qiao等[5]提出了基于卷積神經(jīng)網(wǎng)絡(luò)和長短期記憶神經(jīng)網(wǎng)絡(luò)的時頻雙輸入模型(time-frequency dual-input model based on a convolutional neural network and a long-short term memory neural network,F(xiàn)T-WConvLSTM),此模型的卷積層和長短期記憶層分別提取時頻圖像的空間特征和時域信號的序列特征,最后的密集層被用于故障分類。Li等[6]提出了深度平衡域自適應(yīng)神經(jīng)網(wǎng)絡(luò)(deep balanced domain adaptation neural network,DBDANN),多個卷積層提取特征后,采用多層平衡域自適應(yīng)方法減小不同域之間邊際概率分布和條件概率分布的差異。祝道強(qiáng)等[7]提出了一種基于一維卷積神經(jīng)網(wǎng)絡(luò)的變負(fù)載適應(yīng)軸承故障診斷模型,并在兩種近鄰負(fù)載條件的軸承數(shù)據(jù)構(gòu)成的變負(fù)載數(shù)據(jù)集上取得了較好的診斷效果。

上述算法在變負(fù)載下的滾動軸承故障診斷方面都取得了不錯的效果,但由于時域信號的特征規(guī)律不明顯,以及網(wǎng)絡(luò)丟失重要特征的緣故,使得診斷精度仍有提升的空間。針對上述問題,提出了基于多通道時頻域信號的卷積神經(jīng)網(wǎng)絡(luò)算法(convolutional neural network based on multi-channel time-frequency signals, CNN-MCTFS),用于診斷滾動軸承在變負(fù)載下的健康狀態(tài)。該算法只在源域下訓(xùn)練,無需目標(biāo)域數(shù)據(jù)的參與,即可在目標(biāo)域下得到良好的診斷效果。此算法的貢獻(xiàn)如下:

(1)卷積層采用空洞卷積,也就是向原卷積核中填充空洞,以擴(kuò)大卷積核的感受野。

(2)全局最大池化替代最大池化,在全局范圍內(nèi)提取卷積層的最大激活,而非有限寬度的滑動窗口中的最大激活。

(3)使用不同小波的離散小波變換將時域信號變換為多通道時頻域信號,其中不同的小波提取不同的特征,可以為卷積神經(jīng)網(wǎng)絡(luò)提供更豐富的健康狀態(tài)信息。

1 卷積神經(jīng)網(wǎng)絡(luò)的搭建

卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network, CNN)由輸入層、卷積層、批標(biāo)準(zhǔn)化層、池化層、分類層組成,如圖1所示。

1.1 卷積層(空洞卷積)

卷積層(convolutional layer,Conv)的每個卷積核使用相同的數(shù)值來提取輸入的局部特征[8],如果感受野越大,也就是卷積核越寬那么卷積的輸出包含的信息就越多[9]。VGGNet模型中,三層3×1Conv的感受野相當(dāng)于一層7×1 Conv,因此可以通過加深網(wǎng)路獲取更大的感受野,但CNN層數(shù)過深在訓(xùn)練時容易出現(xiàn)退化現(xiàn)象[10],降低了網(wǎng)絡(luò)的分類精度。如果在Conv中采用7×1寬度的卷積核,雖然擴(kuò)大了卷積核的感受野,但也增加了計(jì)算量。為此,可以使用空洞卷積,也就是向原卷積核中添加空洞以擴(kuò)大卷積核的寬度,從而擴(kuò)大卷積核的感受野。如圖2所示,空洞卷積核的寬度由空洞率(dilation rate)控制,則空洞卷積核的寬度為

Kernel+(Kernel-1)(Dilationrate-1)

(1)

式(1)中:Kernel是原卷積核的寬度。

1.2 池化層(全局最大池化)

池化層的作用就是在提取那些被激活成較大值的重要特征,用于削減特征尺寸,以及保持平移不變性。然而,最大池化只會提取被固定寬度的滑動窗口內(nèi)的最大激活,但如果滑動窗口中都是最大激活,它們只能選出其中之一,而其他值則被拋棄,如圖3(a)所示。因此,在每一Conv后,采用全局最大池化[圖3(b)]代替最大池化,用于提取全局范圍內(nèi)的最大激活,即無需考慮位置,也可保留更多的重要特征。

圖3 兩種池化操作Fig.3 Two pooling operations

(2)

式(2)中:activation和Activation分別是池化前的激活和池化后的激活;GlobalMaxPool(·)是全局最大池化;C和L分別是通道數(shù)和長度。

w和a分別是卷積的可學(xué)習(xí)參數(shù)和輸入值圖2 空洞卷積Fig.2 Dilated convolution

1.3 卷積神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)

在每一Conv之后,激活函數(shù)之前放置批標(biāo)準(zhǔn)化層,以加速網(wǎng)絡(luò)收斂速度,并且起到正則化的作用,防止網(wǎng)絡(luò)過擬合。接著放置激活函數(shù),文中激活函數(shù)選擇了ReLU函數(shù),相比Sigmoid和Tanh等函數(shù),此函數(shù)不僅克服了梯度消失等問題,并且加快了網(wǎng)絡(luò)的訓(xùn)練速度。再放置全局最大池化用于提取全局最大激活,以及削減特征維度。選擇4層Conv,每一Conv的卷積核都相當(dāng)于是7×1的感受野,這樣不僅可以充分學(xué)習(xí)到信號的健康狀態(tài)特征,而且也無需耗費(fèi)太多計(jì)算量。最后全局平均池化將每一通道擠壓為每一點(diǎn),用于最終的分類,并且無論輸入信號的長度或者分解層數(shù)是多少,全局最大池化和全局平均池化的使用可以令網(wǎng)絡(luò)按照表1所示參數(shù)搭建。分類層采用Softmax函數(shù)進(jìn)行分類,以判斷信號的健康狀態(tài)。

表1 卷積神經(jīng)網(wǎng)絡(luò)參數(shù)

2 多通道時頻域信號的構(gòu)造

振動信號在時域下呈波形,特征規(guī)律不明顯,為了凸顯滾動軸承的健康狀態(tài)特征,可以通過傅里葉變換或離散小波變換將時域信號轉(zhuǎn)換到頻域或時頻域。其中傅里葉變換以三角函數(shù)為基擬合時域信號,考察的是信號的整體頻域特征,對于平穩(wěn)信號具有良好效果,而在非平穩(wěn)信號上的效果較差。離散小波變換利用小波,通過伸縮、平移等運(yùn)算對信號進(jìn)行多尺度細(xì)化分析,可以很好地分析非平穩(wěn)信號,從時域信號中提取時頻信息。并且不同的小波從時域信號中提取的時頻信息不同,多種小波可以提供更多樣的時頻信息,以便CNN可以抽取更抽象的健康狀態(tài)特征。

設(shè)f(x)是時域信號,ψ(x)為小波,則連續(xù)小波變換(continuous wavelet transform, CWT)為

(3)

式(3)中:s是尺度因子;τ平移因子; *表示復(fù)共軛操作。離散小波變換是把尺度因子s與平移因子τ離散化,即在離散點(diǎn)取值[11]為

(4)

式(4)中:s≠1;τ是常數(shù)。相應(yīng)的離散小波變換(discrete wavelet transform,DWT)為

(5)

尺寸為(C=1,L)的時域信號首先經(jīng)過level級DWT后得到近似值A(chǔ)level(level表示分解層數(shù))和細(xì)節(jié)系數(shù)(Dlevel,Dlevel-1,…,D1),再將它們串聯(lián)在一起得到尺寸為(C=1,L′)的時頻域信號[Alevel,Dlevel,Dlevel-1,…,D1]。

然而,有多種可用于離散小波變換的小波,例如daubechies(db)小波、symlets(sym)小波、coiflets(coif)小波、haar小波、dmeyer(dmey)小波、biorthogonal(bior)小波和reverse biorthogonal(rbio)小波。因此,在進(jìn)行離散小波變換時,為了保證變換前后信息的保留,避免在變換后引入誤差,應(yīng)選用重構(gòu)誤差低的小波[12-13]。重構(gòu)誤差為

error=sum(|f(x)-iDWT{DWT[f(x),ψ(x)]}|)

(6)

式(6)中:DWT和iDWT分別是離散小波變換和逆離散小波變換;sum(·)表示元素求和操作;|·|表示求絕對值操作。挑選4.1節(jié)數(shù)據(jù)集A的10 000個樣本用于求取重構(gòu)誤差,則各小波的重構(gòu)誤差如表2所示(只列出有限個小波的重構(gòu)誤差,并且level=1)。

表2 各小波的重構(gòu)誤差

從表2可以看到,rbio2.2、bior2.2、coif1和db3小波的重構(gòu)誤差較低且相近,由它們得到的多通道時頻域信號,比只使用rbio2.2小波得到的時頻域信號,蘊(yùn)含著更加多樣化的健康狀態(tài)信息。雖然db5小波的重構(gòu)誤差也和這些小波相似,但由db5小波變換得到的時頻域信號尺寸L′與它們不同,如果通過伸縮操作使得尺寸L′強(qiáng)行一致,可能會缺失信息或引入錯誤信息,故放棄使用。所以,rbio2.2、bior2.2、coif1和db3小波分別對時域信號進(jìn)行l(wèi)evel級(最佳分解層數(shù)之后由實(shí)驗(yàn)選擇)離散小波變換,得到4種不同的時頻域信號,再將這4種不同的時頻域信號在維度C上組合成尺寸為(C=4,L′)的四通道時頻域信號,如圖4所示。具體步驟如下。

圖4 多通道時頻域信號Fig.4 Multi-channel time-frequency signals

3 滾動軸承智能故障診斷算法流程

基于上述闡述,建立起滾動軸承智能故障診斷流程圖,如圖5所示。具體步驟如下。

圖5 滾動軸承智能故障診斷算法流程圖Fig.5 Flow chart of intelligent fault diagnosis algorithm for rolling bearing

Step1重采樣。從加速度傳感器采集原始時域信號,通過重采樣增加每一類的樣本數(shù)量。

Step2構(gòu)造多通道時頻域信號。重構(gòu)誤差選取合適的小波,時域信號被多個不同的小波轉(zhuǎn)換成時頻域信號,從而組成多通道時頻域信號。

Step3確定超參數(shù)。使用驗(yàn)證集確定網(wǎng)絡(luò)的各個超參數(shù),以搭建CNN。本文所用的損失函數(shù)為

(7)

式(7)中:第一項(xiàng)是交叉熵?fù)p失函數(shù);θ表示各個可學(xué)習(xí)參數(shù);p(θ)表示標(biāo)簽的正確概率;q(θ)表示預(yù)測概率。第二項(xiàng)是L2正則化項(xiàng),λ是正則項(xiàng)權(quán)重。使用Adam算法最小化損失函數(shù),并不斷更新可學(xué)習(xí)參數(shù),表示為

(8)

式(8)中:J表示損失函數(shù)值;η表示學(xué)習(xí)率;t表示第t次更新參數(shù)。

Step4訓(xùn)練網(wǎng)絡(luò)。直至訓(xùn)練epoch次后,得到一個訓(xùn)練過的CNN-MCTFS。

Step5測試網(wǎng)絡(luò)。用測試集測試CNN-MCTFS在變負(fù)載下的故障診斷能力。

4 實(shí)驗(yàn)方案及分析

4.1 數(shù)據(jù)集描述

所用到的滾動軸承時域振動信號公共數(shù)據(jù)來自凱斯西儲大學(xué)(Case Western Reserve University,CWRU)軸承數(shù)據(jù)中心。采樣頻率為12 kHz,負(fù)載為1~3 hp,待采集振動信號的滾動軸承型號是SKF6205,用電火花在軸承的內(nèi)圈(inner race,IR)、滾動體(ball)和外圈(outer race,OR)加工損傷直徑分別為0.177 8、0.355 6、0.533 4 mm的單點(diǎn)損傷,再加上正常的類別,即共有10個類別,如表3所示。

從加速度傳感器采集到的原始振動信號是一維時域信號,為了滿足CNN對數(shù)據(jù)量的要求,采用重采樣實(shí)現(xiàn)數(shù)據(jù)集擴(kuò)容。本文的訓(xùn)練集制作過程,就是從各段信號第1個數(shù)據(jù)點(diǎn)開始,每移動32個數(shù)據(jù)點(diǎn)作為一個樣本采集起始點(diǎn),采集連續(xù)的256個數(shù)據(jù)點(diǎn)作為一個樣本,這樣既可以采集到每一類的大量數(shù)據(jù),又可以充分學(xué)習(xí)到類別的特征,如圖6所示。在負(fù)載1~3 hp下,每一類采集到2 100個樣本,其中2 000個屬于訓(xùn)練集,另外的100個屬于驗(yàn)證集。而測試集不可與訓(xùn)練集重疊,所以在每一類別中,測試集需要在第2 100個樣本的最后一個數(shù)據(jù)點(diǎn)之后采集,每256個點(diǎn)采集為一個測試集樣本,從每個類別里各采集100個樣本作為測試集,如表3所示。

4.2 實(shí)驗(yàn)環(huán)境及模型參數(shù)

實(shí)驗(yàn)環(huán)境是Ubuntu 16.04 LTS,Intel Core i7-67000 CPU@3.40 GHz,16 GB RAM,GeForce GTX 1080,CUDA 11.0,Python 3.6.9,PyTorch 0.4.1。根據(jù)上述滾動軸承智能故障診斷流程,在驗(yàn)證集上經(jīng)過多次試驗(yàn),得到本文所用算法超參數(shù)如表4所示。

表4 超參數(shù)

4.3 實(shí)驗(yàn)結(jié)果及分析

4.3.1 分解層數(shù)的選取

level不同得到的多通道時頻域信號不同,對最終的分類結(jié)果也會造成不同的影響。最大分解層數(shù)由信號長度和小波長度共同決定,當(dāng)信號長度分解到小于小波長度時,就達(dá)到了最大分解層數(shù)。在本文信號長度下,rbio2.2、bior2.2、coif1和db3小波的最大分解層數(shù)是5,因此level=1~5的分類精度如圖7所示。從圖7中可知,所有l(wèi)evel的平均精度(average accuracy, AVG)均達(dá)到了0.950以上,由此證明該算法在任何level下都達(dá)到了較好的分類效果。但level=5的效果最好,AVG達(dá)到了0.975,相比其他level至少高出了0.012。這是由于level越大,輸入到CNN的多通道時頻域信號蘊(yùn)含的多尺度信息越多,使得網(wǎng)絡(luò)的分類精度就越高。因此,在接下來的實(shí)驗(yàn)中,均是對時域信號進(jìn)行l(wèi)evel=5的DWT。

圖7 不同level的精度Fig.7 Different levels of accuracy

4.3.2 與其他算法的比較

在此實(shí)驗(yàn)中,將CNN-MCTFS與SVM(support vector machines)、MLP(multi-Layer perceptron)、文獻(xiàn)[14]、WDCNN[15]和TICNN[8]進(jìn)行對比實(shí)驗(yàn),實(shí)驗(yàn)結(jié)果如圖8所示。SVM的AVG只達(dá)到了0.614,也從側(cè)面驗(yàn)證了傳統(tǒng)智能故障診斷算法在變負(fù)載下的效果較差。當(dāng)由MLP進(jìn)行實(shí)驗(yàn)時,AVG明顯提升了0.069,說明神經(jīng)網(wǎng)絡(luò)的擬合輸入的效果強(qiáng)于SVM等傳統(tǒng)機(jī)器學(xué)習(xí)算法。之后基于CNN的算法的AVG都達(dá)到了0.910以上,證明CNN的泛化性能較強(qiáng)。并且CNN-MCTFS的AVG達(dá)到了0.975,相比TICNN提高了0.043,相比WDCNN提高了0.010,因此該算法在變負(fù)載下可以有效識別滾動軸承的健康狀態(tài)。

圖8 不同算法之間的比較Fig.8 Comparison between different algorithms

而且,在CNN-MCTFS中,A-B的0.992要高于A-C的0.986,C-B的0.987要高于C-A的0.920,可以從中發(fā)現(xiàn)算法在相鄰負(fù)載的效果要好于負(fù)載相差較大的,因此要盡可能收集多種負(fù)載下的數(shù)據(jù),縮小訓(xùn)練集負(fù)載與測試集負(fù)載的差距。并且算法在中低負(fù)載下訓(xùn)練的效果要強(qiáng)于在高負(fù)載下的,例如A-C的0.986要高于C-A的0.920,B-C的0.996要高于C-B的0.987,因此也要盡可能地多收集中低負(fù)載下的信號作為訓(xùn)練集。

4.3.3 消融實(shí)驗(yàn)

將CNN-MCTFS與最原始的CNN(MaxPool, no-dilation)-T,將多通道時頻域信號替換成時域信號、頻域信號、單通道時頻域信號的CNN-T、CNN-F、CNN-rbio2.2,將全局最大池化替換成最大池化的CNN(MaxPool)-MCTFS,將空洞卷積替換成3×1卷積的CNN(no-dilation)-MCTFS相比較,實(shí)驗(yàn)結(jié)果如圖9所示。由于時域信號的特征不明顯,時域信號在淺層網(wǎng)絡(luò)中的AVG最高達(dá)到了0.916,而如果將時域信號變化到頻域或時頻域中,AVG則達(dá)到了0.948和0.939,而CNN-MCTFS的AVG相比它們至少多出了0.027,說明多個小波的共同使用可以挖掘出時域信號更豐富的特征,可以讓CNN提取到更抽象的特征,從而提高分類精度。在將最大池化替換成全局最大池化,提取了每一Conv的最大激活后,CNN-MCTFS的AVG相比CNN(MaxPool)-MCTFS高出了0.024,而提高每一個Conv的卷積核的感受野之后,CNN-MCTFS的AVG相比CNN(no-dilation)-MCTFS提升了0.036。因此,CNN-MCTFS在變負(fù)載下取得了較高的診斷精度。

圖9 各貢獻(xiàn)點(diǎn)之間的比較Fig.9 Comparison between contributions

4.3.4 可視化各層特征

為進(jìn)一步說明CNN-MCTFS的負(fù)載域適應(yīng)能力,在C-A中,利用t-SNE(t-distributed stochastic neighbor embedding)算法[16]可視化數(shù)據(jù)集A在每一層的分布情況,如圖10所示。在輸入層中,所有樣本聚成一簇。經(jīng)過第1層Conv后,正常類樣本逐漸聚集在一起,并與其他類樣本分離,但每一類樣本已經(jīng)有了聚類的趨勢。經(jīng)過第2層Conv之后,正常類樣本首先聚集在一起,并游離在外。經(jīng)過第3、4層Conv之后,所有樣本被分成10類,但還有少量樣本散落在其他簇中。在分類層,基本上所有樣本都被正確分類,但只有少數(shù)樣本被錯分。如圖11所示,主要是IR0.355 6 mm被錯分為OR0.355 6 mm和OR0.533 4 mm,ball0.533 4 mm被錯分為ball0.177 8 mm,因此之后的研究要強(qiáng)化對這兩類樣本的正確分類。

圖10 C-A中各層的特征分布Fig.10 In C-A, the feature distribution of each layer

圖11 混淆矩陣Fig.11 Confusion matrix

5 結(jié)論

在滾動軸承故障診斷中,為有效診斷滾動軸承在變負(fù)載下的健康狀態(tài),需要算法具有較強(qiáng)的負(fù)載域適應(yīng)能力。提出了基于多通道時頻域信號的卷積神經(jīng)網(wǎng)絡(luò)算法(CNN-MCTFS),得出如下結(jié)論。

(1)空洞卷積擴(kuò)大了每一Conv的卷積核寬度,使得CNN在無需加深網(wǎng)絡(luò)以及計(jì)算量的前提下,獲得了較大的感受野。

(2)最大池化被替換成了全局最大池化,從全局范圍內(nèi)選取最大激活,從而最大程度地保留了許多重要特征。

(3)與時域信號、頻域信號、單通道時頻域信號相比,多個小波的共同使用為CNN提供了更豐富的特征。

(4)最終,CNN-MCTFS在變負(fù)載下的AVG達(dá)到了0.975。上述實(shí)驗(yàn)過程表明,在沒有目標(biāo)域數(shù)據(jù)的參與下,CNN-MCTFS也可以有效識別出滾動軸承的健康狀態(tài),相比于其他算法具有更強(qiáng)的負(fù)載域適應(yīng)能力,在軸承類設(shè)備的故障診斷中具有廣闊的應(yīng)用前景。

猜你喜歡
故障診斷特征信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
如何表達(dá)“特征”
不忠誠的四個特征
基于FPGA的多功能信號發(fā)生器的設(shè)計(jì)
電子制作(2018年11期)2018-08-04 03:25:42
抓住特征巧觀察
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應(yīng)用
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
主站蜘蛛池模板: 91麻豆精品国产高清在线| 欧美全免费aaaaaa特黄在线| 国产福利在线免费| 澳门av无码| 人人澡人人爽欧美一区| 国产欧美视频在线观看| 国产人成午夜免费看| jizz国产视频| 不卡色老大久久综合网| 日韩欧美中文| 国产尤物在线播放| 国产杨幂丝袜av在线播放| 精品無碼一區在線觀看 | 亚洲精品无码专区在线观看| 国产精品成人AⅤ在线一二三四| 91视频区| 人妻丰满熟妇AV无码区| 在线视频亚洲欧美| 精品国产乱码久久久久久一区二区| 国产日韩欧美视频| 999福利激情视频| 亚洲Av激情网五月天| 日韩成人午夜| 91香蕉国产亚洲一二三区 | 欧美α片免费观看| 免费一级全黄少妇性色生活片| 亚洲成人播放| 国产成人调教在线视频| 亚洲av中文无码乱人伦在线r| 无码免费试看| 久久久久国产精品免费免费不卡| 国产精品jizz在线观看软件| 九九视频免费在线观看| 中文字幕欧美日韩| 国产乱人乱偷精品视频a人人澡| 国产自产视频一区二区三区| 成人福利在线免费观看| 欧美一级大片在线观看| 视频一区视频二区中文精品| 91在线激情在线观看| 欧美一级大片在线观看| 日本人真淫视频一区二区三区| 国产swag在线观看| 日本在线国产| 亚洲一区二区日韩欧美gif| 亚洲中文字幕久久无码精品A| 久久久黄色片| 精品少妇人妻av无码久久| 在线观看免费黄色网址| 色妞永久免费视频| 女人18一级毛片免费观看| 国产原创第一页在线观看| 四虎影视永久在线精品| 日本妇乱子伦视频| 欧美日韩在线亚洲国产人| 亚洲精品自拍区在线观看| 国产尤物jk自慰制服喷水| 超碰精品无码一区二区| 一级毛片在线播放| 国产资源站| 全免费a级毛片免费看不卡| 国产欧美日韩91| 国产精品久久久久久搜索| 国内精品久久人妻无码大片高| 日韩美毛片| 欧美色亚洲| 亚洲伊人电影| 欧美日韩国产高清一区二区三区| 精品一区二区三区视频免费观看| 欧美亚洲欧美区| 久久久久久久久18禁秘| 四虎国产永久在线观看| 亚洲,国产,日韩,综合一区 | 青青草国产一区二区三区| 福利一区在线| 亚洲中文字幕23页在线| 成人在线综合| 国产激爽大片高清在线观看| 日韩欧美高清视频| 婷五月综合| 日本免费福利视频| 欧美成人A视频|