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

基于EWT-OPRCMDE-ELM的風電機組齒輪箱故障診斷研究

2022-01-17 08:28:10羅興琦
自動化儀表 2021年11期
關鍵詞:故障診斷特征故障

李 輝,李 宣,賈 嶸,羅興琦,白 亮

(1.西安理工大學電氣工程學院,陜西 西安 710048;2.西安理工大學水利水電學院,陜西 西安 710048)

0 引言

惡劣的工作環(huán)境導致風機故障頻發(fā)[1]。齒輪箱作為風電機組傳動系統(tǒng)的主要設備,更是故障的高發(fā)部位。據(jù)統(tǒng)計,齒輪箱故障占機組總故障60%以上[2],而一臺風力機故障停機時間的20%是由齒輪箱故障引起的。一旦發(fā)生故障,機組將面臨長時間停機和昂貴的維修費用,損失巨大[3]。因此,準確、高效地對風機齒輪箱進行故障診斷,對于保障機組的安全、穩(wěn)定運行具有重要意義。

由于齒輪箱故障振動信號受背景噪聲、交變載荷和傳輸路徑等影響,呈現(xiàn)非平穩(wěn)、非線性的特點,通常用時頻分析方法對其處理,如小波變換、經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)、變分模態(tài)分解(variational mode decomposition,VMD)等。但是小波變換中小波基的選擇不當會對處理效果產(chǎn)生不良影響。EMD方法缺乏理論基礎且存在模態(tài)混疊、端點效應等問題,所以無法有效分析故障信號。VMD方法依賴于分解層數(shù)和懲罰因子的設定,如參數(shù)設置不當,將極大地影響地處理效果。針對傳統(tǒng)信號處理方法的不足,Gilles[4]結(jié)合經(jīng)驗模態(tài)分解的自適應性和小波分析的框架理論,提出了經(jīng)驗小波變換(empirical wavelet transform,EWT)方法,有效地改善了傳統(tǒng)時頻分析方法的弊端。文獻[5]將EWT應用于旋轉(zhuǎn)機械故障診斷中,并證明EWT處理效果優(yōu)于傳統(tǒng)EMD方法。

高效地提取區(qū)分度大的穩(wěn)定特征是風電機組齒輪箱故障診斷的重要環(huán)節(jié)。學者們通常提取故障信號的時域、頻域或時頻域特征,進行故障診斷。文獻[6]提取軸承振動信號的均值、峭度等時域特征,進行故障診斷,取得較好的診斷效果。文獻[7]提取軸承振動信號的時域、頻域和時頻域的多個混合特征參量,利用XGBoost進行風機軸承故障診斷。文獻[8]融合風機齒輪箱振動信號的時域、頻域和時頻域特征,特征降維后利用隨機森林算法進行故障模式識別。但是,由于振動信號受強烈噪聲干擾、多振動源耦合和傳輸路徑等影響,導致基于傳統(tǒng)時域、頻域和時頻域特征提取的故障診斷方法效果不佳,故障診斷正確率和診斷時間受特征參量的維數(shù)影響較大。精細復合多尺度散布熵[9](refined composite multi-scale dispersion entropy,RCMDE)作為一種衡量時間序列復雜度的非線性動力學參數(shù),已被成功應用于心電信號和腦電信號的分析中,并被證明在計算速度和穩(wěn)定性方面優(yōu)于樣本熵和排列熵。其特征提取效果優(yōu)于傳統(tǒng)時域、頻域特征。

極限學習機(extreme learning machine,ELM)[10]是一種泛化的單隱層前饋神經(jīng)網(wǎng)絡的機器學習算法。其輸入層和隱含層之間的連接權(quán)值和隱含層閾值隨機產(chǎn)生,訓練過程無需迭代調(diào)整,與傳統(tǒng)反向傳播(back propagation,BP)神經(jīng)網(wǎng)絡、支持向量機等智能分類算法相比,輸入?yún)?shù)少、運算速度更快。

為提高風電機組齒輪箱故障診斷效果,本文提出一種基于EWT重構(gòu)、最優(yōu)參數(shù)精細復合多尺度散布熵和ELM結(jié)合的故障診斷方法,并將其應用于風電機組齒輪箱。

1 經(jīng)驗小波變換

EWT是通過對信號的傅里葉頻譜進行自適應分割,根據(jù)分割邊界,構(gòu)建正交小波濾波器組以提取包含故障信息的調(diào)幅調(diào)頻分量。步驟如下。

①對故障信號進行傅里葉變換,將信號的傅里葉頻譜定義在[0,π]范圍內(nèi),并將[0,π]分成N個帶寬不等的頻帶,則每個頻帶可以表示為Λn=[ωn-1,ωn]。以ωn為中心,定義Tn=2τn為過渡帶。傅里葉軸分割如圖1所示。

圖1 傅里葉軸分割Fig.1 Partitioning of the Fourier axis

②經(jīng)驗小波的尺度函數(shù)φn(ω)和小波函數(shù)ψn(ω)在頻域里的計算式如下:

(1)

(2)

③根據(jù)經(jīng)典小波變換的方法構(gòu)造經(jīng)驗小波變換,細節(jié)系數(shù)由經(jīng)驗小波函數(shù)內(nèi)積計算得出,近似系數(shù)由經(jīng)驗尺度函數(shù)內(nèi)積計算得出:

(3)

(4)

原信號重構(gòu)表達式如下:

(5)

式中:∧表示傅里葉變換;∨表示傅里葉逆變換;—表示復共軛;?表示卷積。

經(jīng)驗模態(tài)fk(t)定義如下:

(6)

2 精細復合多尺度散布熵

散布熵計算步驟如下。

①x={xj|j=1,2,…,N}為一時間序列,利用正態(tài)分布函數(shù)將時間序列x映射到y(tǒng)={yj|j=1,2,…,N},yj∈(0,1),即:

(7)

式中:μ、σ分別為期望和標準差。

②通過線性變換將y映射到{1,2,…,c}范圍內(nèi):

(8)

(9)

式中:m、d分別為嵌入維數(shù)和時延。

④計算散布模式πv0v1…vm-1,v=1,2,…,c。

(10)

⑤計算每種散布模式的概率πv0v1…vm-1:

(11)

⑥根據(jù)香農(nóng)熵定義,有:

(12)

散步熵(dispersion entropy,DE)的值越大,信號時間序列不規(guī)則程度越高;反之,不規(guī)則程度越低。

RCMDE計算步驟如下。

RCMDE是對原始數(shù)據(jù)精細化和多尺度化,其作用是減小初始點位置對計算結(jié)果的影響和充分利用數(shù)據(jù)序列的信息。計算過程如下。當尺度因子為S時:首先,將原始信號按初始點連續(xù)分割成長度為S的小段,求取每一段的平均值,再將這些平均值按順序排列成S個粗粒化序列;其次,計算每個粗粒化序列散布模式π的概率,并計算這些散布模式概率的平均值;最后,根據(jù)香農(nóng)熵的定義計算RCMDE的值。

(13)

②對于每個尺度S,RCMDE值定義如下:

(14)

3 Relief-F與ELM

提取齒輪箱故障特征后,由于存在特征冗余,會使診斷時長和精度受影響。故采用Relief-F算法[11]計算分類權(quán)重,剔除冗余。Relief-F算法是從Relief算法擴展而來的,可以處理多類別分類權(quán)重問題。該算法計算過程如下。

①從所有樣本中,隨機取出一個樣本a。

②在與樣本a相同分類的樣本組內(nèi),取出k個最近鄰樣本。

③在所有其他與樣本a不同分類的樣本組內(nèi), 也分別取出k個最近鄰樣本。

④計算每個特征的權(quán)重。

極限學習機是一種單隱層前饋神經(jīng)網(wǎng)絡,由輸入層、隱含層和輸出層構(gòu)成。與傳統(tǒng)神經(jīng)網(wǎng)絡不同的是,極限學習機的輸入層與隱含層之間的連接權(quán)值和隱含層神經(jīng)元的閾值為隨機設定,在開始創(chuàng)建網(wǎng)絡時隨機產(chǎn)生,且設置完成后無需調(diào)整。隱含層與輸出層之間的連接權(quán)值由方程組依次解出,無需迭代計算。故極限學習機可在保證學習精度的前提下極大地加快學習速度。

極限學習機訓練過程如下。假設有N個任意樣本(Xi,ti)。其中,Xi=[xi1xi2…xin]T∈Rn,ti=[ti1ti2…tin]∈Rm,n為輸入層的維度,m為輸出層的維度。

當隱層節(jié)點的個數(shù)為L時,單隱層前饋神經(jīng)網(wǎng)絡輸出可以表示為:

(15)

式中:i=1,2,…,N;Wj=[wj1wj2…wjn]T為輸入權(quán)值;βj=[βj1βj2…βjm]T為輸出權(quán)重;bj是第j個隱層單元的偏置;WjXj為Wj與Xj的內(nèi)積;g()為激活函數(shù)。

單隱層前饋神經(jīng)網(wǎng)絡學習的目標是在最小的誤差下逼近N個樣本,可以表示為:

(16)

即存在Wj、βj、bj,使得:

(17)

式(17)可簡化為Hβ=T。其中:H為隱層節(jié)點輸出矩陣;β為輸出權(quán)重;T為期望輸出。

(18)

1,2,…,L

(19)

當輸入權(quán)重Wj和隱含層偏置bj隨機確定后,隱層的輸出矩陣H就會被確定。此時,單隱層神經(jīng)網(wǎng)絡求解的問題被轉(zhuǎn)換成線性系統(tǒng)Hβ=T的求解,則輸出權(quán)重β可由式(20)確定。

(20)

極限學習機診斷過程如下。

①設定激活函數(shù)g()和隱層神經(jīng)元個數(shù)N。

②輸入訓練集進行網(wǎng)絡訓練。

③利用訓練好的ELM,輸入測試集進行故障診斷。

4 故障診斷流程

故障診斷流程如下。

①利用振動傳感器采集風電機組齒輪箱軸承故障振動信號。

②利用EWT分解原始振動信號,將得到若干子模態(tài)分量EWF。

③計算各子模態(tài)分量與原始信號的相關系數(shù)Rn。設定相關系數(shù)閾值為0.3,選取大于閾值的子模態(tài)分量EWF進行信號重構(gòu),提取故障主要信息。

相關系數(shù)是一種度量2個反映之間相關程度的系數(shù),其值介于1和-1之間。其中:1表示變量完全正相關;0表示無關;-1表示完全負相關。皮爾遜相關系數(shù)定義為兩個變量之間的協(xié)方差和標準差的商。函數(shù)表達式為:

(21)

④以重構(gòu)信號RCMDE偏度值的平方函數(shù)為適應度函數(shù),計算重構(gòu)信號最優(yōu)RCMDE參數(shù),提取重構(gòu)信號最優(yōu)參數(shù)RCMDE(optimal parameter RCMDE,OPRCMDE)構(gòu)成故障特征向量。

⑤利用Relief-F計算OPRCMDE各尺度特征分類權(quán)重,并進行特征參量篩選,剔除冗余特征。

⑥將篩選得到的特征向量作為輸入,利用ELM進行故障診斷。

5 信號仿真驗證

為驗證EWT分解在處理復雜信號時的優(yōu)越性。采用仿真信號進行驗證。仿真信號f由指數(shù)信號f1、調(diào)頻信號f2、噪聲信號f3(強度為0.2 dB·W)構(gòu)成。仿真信號時域如圖2所示。

圖2 仿真信號時域圖Fig.2 Time domain diagram of simulation signal

利用EWT分解信號f,同時與EMD方法對比。EWT分解結(jié)果時域如圖3所示。EWT分解結(jié)果的Hilbert譜如圖4所示。

圖3 EWT分解結(jié)果時域圖Fig.3 Time domain diagram of EWT decomposition results

圖4 EWT分解結(jié)果的Hilbert譜圖Fig.4 Hilbert spectrum of EWT decomposition results

(22)

根據(jù)圖3和圖4,由 EWT分解結(jié)果可得,EWT將模擬信號f分解為4個子模態(tài)分量。 其中,EWF1在[0,1]上單調(diào)增加,最大值為5,對應于f1的變化趨勢。 EWF2的周期約為0.1 s,其波形和周期對應于f2中的sin(20πt)。而EWF3的周期約為0.05 s,對應于f2中cos(40πt)。 EWF4沒有明顯的特征,并且在時域中很雜亂,因此被判斷為噪聲信號。 同時,在Hilbert頻譜中可以觀察到清晰、穩(wěn)定的10 Hz和20 Hz對應于f2的特征頻率。EMD分解結(jié)果時域如圖5所示。

圖5 EMD分解結(jié)果時域圖Fig.5 Time domain diagram of EMD decomposition results

EMD分解結(jié)果的Hilbert譜圖如圖6所示。觀察圖5和圖6,EMD將模擬信號分解為8個子模態(tài)分量,但只有Res分量具有與其相同的變化趨勢,并且很難清楚地觀察到其余IMF分量中與模擬信號分量相對應的分量。 而且,由于噪聲干擾,Hilbert頻譜中的特征頻率不規(guī)律且雜亂混疊,并且難以觀察到穩(wěn)定的模擬信號中包含的特征頻率。

圖6 EMD分解結(jié)果的Hilbert譜圖Fig.6 Hilbert spectrum of EMD decomposition results

上述仿真試驗中,通過比較EWT和EMD的分解結(jié)果,可以發(fā)現(xiàn)EMD方法會過度分解信號,并且存在模態(tài)混疊現(xiàn)象,不能有效地提取故障特征分量。 EWT可以在嘈雜的環(huán)境中有效地提取信號的主要成分,以獲得更好的特征提取結(jié)果。

6 試驗分析

本節(jié)結(jié)合了EWT信號重構(gòu)和OPRCMDE,并使用極限學習機進行故障診斷。采用基于千鵬 QPZZ-Ⅱ 旋轉(zhuǎn)機械振動及故障模擬試驗平臺的齒輪箱振動數(shù)據(jù)進行分析。試驗平臺由變速驅(qū)動電機(0.75 kW)、軸承、齒輪箱、軸、偏重轉(zhuǎn)盤、調(diào)速器等組成。通過調(diào)節(jié)配重、調(diào)節(jié)部分的安裝位置以及組件的有機組合,快速模擬各種齒輪箱故障。

齒輪箱參數(shù)如下。

①齒輪。大齒輪:模數(shù)2、齒數(shù)75、材質(zhì)S45C(3只,其中2只備件齒輪);小齒輪:模數(shù)2、齒數(shù)55、材質(zhì)S45C(2只,其中1只備件齒輪)。

②潤滑:浸油式。

③軸承:滾動軸承(N205)。

設定轉(zhuǎn)速為880 r/min,采樣頻率為5 120 Hz,通過安裝在齒輪箱輸出軸負載側(cè)軸承X方向的加速度傳感器,采集齒輪箱正常、點蝕、磨損和斷齒這4類故障狀態(tài)的振動信號各25組樣本。每個樣本長度為2 048。

6.1 EWT分解與重構(gòu)

以齒輪箱齒輪磨損故障為例。EWT將磨損故障信號分解為10個EWF。根據(jù)本文所提到的相關系數(shù)閾值,選取EWF2、EWF6為最優(yōu)分量重構(gòu)故障信號并進行下一步的故障特征提取及模式識別。

EWF相關系數(shù)如圖7所示。

圖7 EWF相關系數(shù)Fig.7 Correlation coefficient of EWF

6.2 OPRCMDE參數(shù)計算及特征提取

選取合適的RCMDE參數(shù),能夠有效提取區(qū)分度較大的故障特征向量,同時減少計算時間。針對不同的故障類型選取對應的最佳參數(shù),可獲取最佳的故障特征區(qū)分度。

在精細復合多尺度散布熵的計算中,涉及4個參數(shù)的選取,即尺度因子S、嵌入維數(shù)m、類別個數(shù)C和時延d。為充分利用軸承故障數(shù)據(jù),設置尺度因子S=15,時延d對算法影響較小,故設置為1。嵌入維數(shù)m、類別個數(shù)C的設定對特征提取效果影響重大。m過大,將導致RCMDE運算時間加長;m過小,則算法檢測信號突變性減弱。C過大,將導致RCMDE對噪聲十分敏感;C過小,則彼此相距較遠的兩個振幅可能被分到相同的類別內(nèi)。考慮到m和C之間的交互作用,使用網(wǎng)格搜索算法對范圍內(nèi)的最優(yōu)m值和C值進行同步搜索。

為分析多尺度熵值的集中趨勢,可計算其均值。但僅憑均值不足以表征多尺度熵值的總體趨勢。本文引入精細復合多尺度散布熵的偏度。其絕對值越大,表明均值的效能越有問題;其絕對值越小,表明均值越可信賴。本文以EWT重構(gòu)信號精細復合多尺度散布熵偏度的平方函數(shù)[12]作為適應度函數(shù),求取其最小值。偏度計算過程如下:

信號序列X={xi},i=1,2,…,N。15個尺度下的精細復合多尺度散布熵值組成序列:

BP(X)={BP(1),BP(2),…,BP(15)}

(23)

偏度值為:

(24)

適應度函數(shù)取為:

(25)

設定m搜索范圍為[2,6]、C搜索范圍為[3,8]、設定步長均為1。通過網(wǎng)格搜索算法,得到4類故障OPRCMDE參數(shù),如表1所示。

表1 4類故障OPRCMDE參數(shù)

提取4類故障狀態(tài)EWT重構(gòu)信號的OPRCMDE值,得到4組25×15故障特征矩陣。原始信號推薦參數(shù)RCMDE如圖8所示(4類故障狀態(tài)均設定m=3、C=6)。EWT重構(gòu)信號OPRCMDE如圖9所示。

圖8 原始信號推薦參數(shù)RCMDEFig.8 RCMDE of recommended parameters of original signal

圖9 EWT重構(gòu)信號OPRCMDEFig.9 OPRCMDE of the EWT reconstructed signal

對比圖8、圖9,可觀察到本文方法所提取的4種故障狀態(tài)特征區(qū)別明顯,有利于后期的故障特征提取和故障診斷。

由于高維OPRCMDE特征向量存在冗余,利用Relief-F算法(參數(shù)設置迭代次數(shù)M=20,近鄰值KL=5)計算每個尺度因子下OPRCMDE值的特征權(quán)重,并進行特征向量二次篩選和約簡。各尺度因子分類特征權(quán)重如圖10所示。

圖10 各尺度因子分類特征權(quán)重Fig.10 Feature weights at each scale

根據(jù)權(quán)重大小,選取第1、第2、第6共3個尺度因子下的OPRCMDE值作為故障特征向量,約簡后得到4組25×3的故障特征矩陣。Relief-F約簡特征散點圖如圖11所示。觀察圖11可發(fā)現(xiàn):各類故障間聚集緊密從而可將4種不同故障狀態(tài)明顯地區(qū)分開。

圖11 Relief-F約簡特征散點圖Fig.11 Relief-F reduced feature scatter

6.3 極限學習機故障診斷

為4故障狀態(tài)貼上標簽,正常、點蝕、磨損、斷齒的標簽分別為1、2、3、4。選取60%樣本作為訓練樣本,剩余40%樣本作為測試樣本。測試集40個樣本分類結(jié)果與實際類別一致,故障診斷正確率100%,對應類別未出現(xiàn)錯誤分類。ELM診斷結(jié)果如圖12所示。

圖12 ELM診斷結(jié)果Fig.12 Results of ELM diagnosis

根據(jù)試驗結(jié)果,證明本文提出的方法能夠提取風機齒輪箱區(qū)分度較大的故障特征,實現(xiàn)了對齒輪箱不同故障的準確診斷。為進一步說明本文方法的優(yōu)越性,將本文方法與多尺度排列熵(multi-scale permutation entropy,MPE)-ELM(MPE參數(shù)設置為m=5、S=12、T=1)、MPE-ReliefF-ELM、RCMDE-ELM(RCMDE參數(shù)推薦設置為m=3、C=6)、時頻域融合特征-ReliefF-ELM(時頻域融合特征共11個,其中8個時域特征分別為平均值、有效值、峰值、峰值因子、峭度、波形因子、脈沖因子、裕度因子,3個頻域特征分別為重心頻率、均方頻率、頻率方差)進行對比。每種方法運行10次,取診斷正確率平均值。各診斷模型診斷效果對比如表2所示。

表2 各診斷模型診斷效果對比

由表2中結(jié)果可知,當本文方法所提取的齒輪箱故障特征應用于齒輪箱故障診斷時,診斷準確率更高且更加穩(wěn)定。

7 結(jié)論

本文針對傳統(tǒng)風電機組齒輪箱故障診斷方法的不足,提出一種基于EWT重構(gòu)、最優(yōu)參數(shù)精細復合多尺度散布熵和極限學習機結(jié)合的故障診斷方法。具體結(jié)論如下。

①EWT解決了傳統(tǒng)時頻處理方法存在的模態(tài)混疊、端點效應、參數(shù)設置等問題。通過相關系數(shù)選取子模態(tài)分量進行信號重構(gòu),提取了故障的主要信息。

②為提高精細多尺度散布熵算法的特征提取性能,以其偏度值的平方函數(shù)作為適應度函數(shù),通過網(wǎng)格搜索算法同步搜索計算2個關鍵參數(shù)(即m和C)。通過試驗對比,證明EWT重構(gòu)信號最優(yōu)參數(shù)精細復合多尺度散布熵在提取各類故障特征時區(qū)分度更好,診斷結(jié)果更穩(wěn)定準確。

③針對特征向量冗余和一般分類算法參數(shù)多且參數(shù)設定影響分類準確率的問題,采用Relief-F算法計算特征向量的分類權(quán)重,選擇權(quán)重大者構(gòu)成最終的特征向量,剔除了冗余特征。最后,再利用運算速度快、參數(shù)設置少的ELM進行故障診斷。試驗結(jié)果證明,該方法能夠有效運用于風電機組齒輪箱故障診斷中。

猜你喜歡
故障診斷特征故障
故障一點通
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
奔馳R320車ABS、ESP故障燈異常點亮
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: P尤物久久99国产综合精品| 午夜毛片福利| 国产精品视频第一专区| 亚洲国产日韩在线成人蜜芽| 人妻精品久久无码区| 在线观看亚洲精品福利片| 91福利国产成人精品导航| 免费在线成人网| 97综合久久| 国产成人一区| 国产精品综合久久久 | 亚洲人成网18禁| 欧美精品二区| 日本成人不卡视频| 精品一区二区无码av| 欧美啪啪一区| 久久精品一品道久久精品| 久热中文字幕在线| 婷婷激情五月网| 很黄的网站在线观看| 亚洲一级毛片免费观看| 粉嫩国产白浆在线观看| 免费人成黄页在线观看国产| 成人午夜在线播放| 成年看免费观看视频拍拍| 精品午夜国产福利观看| 亚洲欧洲一区二区三区| 白浆免费视频国产精品视频| 欧美在线伊人| 呦视频在线一区二区三区| 色婷婷狠狠干| 亚洲欧美一级一级a| 日本人妻丰满熟妇区| 欧美精品H在线播放| 午夜限制老子影院888| 精品国产美女福到在线直播| 欧美国产视频| 五月综合色婷婷| 不卡午夜视频| 日韩国产高清无码| 日韩精品无码不卡无码| 国产真实自在自线免费精品| 蜜臀AV在线播放| 国产精品九九视频| 国产区人妖精品人妖精品视频| 欧美成人怡春院在线激情| 永久免费无码成人网站| 欧美性爱精品一区二区三区| 九九视频在线免费观看| 99久久精彩视频| 丁香婷婷激情综合激情| 制服丝袜亚洲| 日韩麻豆小视频| 一级毛片免费观看久| 九九久久99精品| 中文字幕亚洲另类天堂| 欧美一级黄片一区2区| jizz亚洲高清在线观看| 成人在线第一页| 性69交片免费看| A级全黄试看30分钟小视频| 97精品久久久大香线焦| 色视频国产| 伊人久综合| 人妖无码第一页| 99精品免费在线| 色九九视频| 亚洲精品大秀视频| 久久久精品久久久久三级| 中文字幕在线日本| 毛片在线播放a| www成人国产在线观看网站| 婷婷六月天激情| 欧日韩在线不卡视频| 欧美、日韩、国产综合一区| 尤物精品视频一区二区三区 | 高清免费毛片| 亚洲成a人片| 精品三级网站| 一区二区偷拍美女撒尿视频| 国产精品永久免费嫩草研究院| 色爽网免费视频|