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

基于OHF Elman-AdaBoost算法的滾動軸承故障多時期診斷方法

2021-03-31 06:32:12卓鵬程夏唐斌鄭美妹奚立峰
振動與沖擊 2021年6期
關鍵詞:故障診斷故障信號

卓鵬程,夏唐斌,2,鄭美妹,鄭 宇,奚立峰,2

(1.上海交通大學 機械與動力工程學院,上海 200240;2.上海交通大學 機械系統與振動國家重點實驗室,上海 200240)

隨著我國工業水平的不斷發展,工業領域的設備在不斷趨向大型化、復雜化。大型復雜裝備的平穩安全運行關系到生產安全的關鍵核心[1]。滾動軸承作為大型復雜設備的重要組成部分,其性能退化或失效將影響整機性能甚至導致設備非計劃停機,造成嚴重經濟損失甚至重大人員傷亡。據統計,約45%~55%的旋轉機械故障是由于滾動軸承的損傷所造成的[2]。然而當滾動軸承處于故障初期時,其故障信息往往表現比較模糊,振動信號異常不明顯。同時因為缺乏對軸承故障程度的掌控,往往當滾動軸承處于故障晚期時才有明顯表征,診斷滯后會造成嚴重的失效風險。因此,對于滾動軸承故障中早期的診斷則十分必要,可以有效地預先設計維修與更換方案。在當前制造業中,對滾動軸承建立退化模型或進行剩余壽命預測需要大量全生命周期狀態數據,且預測模型與實際退化狀態之間的誤差較大。因此,實現滾動軸承實時精確的故障初期、中期及晚期診斷,可以有效降低事故發生率,對于提高企業的制造運維和健康管理水平具有重要意義。

滾動軸承故障診斷技術目前主要包括兩個方面:故障信號提取與故障模式識別。在故障信號提取方面,目前常用的方法有Hilbert-Huang變換、小波變換(wavelet transform,WT)、經驗模態分解(empirical mode decomposition,EMD)與集合經驗模態分解(ensemble empirical mode decomposition,EEMD)等[3-6]。EEMD對于滾動軸承振動信號有著較好的去噪效果,可以有效地提取出滾動軸承早期故障信息。在故障模式識別方面,常用的方法包括概率統計方法、支持向量機(support vector machine,SVM)和人工神經網絡(artificial neural networks,ANN)等[7-12]。滾動軸承故障早期數據具有非線性、強干擾與數據源復雜等特性,概率統計方法與支持向量機均存在一定的局限性,目前神經網絡由于其自適用能力與容錯性,成為最熱門的滾動軸承故障診斷方法。Shi等[13]基于Elman神經網絡提出OHF Elman (output hidden feedback Elman) 網絡,并證明了OHF Elman相較Elman神經網絡的優越性。但是,單個OHF Elman神經網絡無法保證對所有測試樣本都有較好的適用性。為了實現更為精確的滾動軸承全樣本故障診斷方法,Singh等[14]通過Bagging (bootstrap aggregating) 算法將樸素貝葉斯(Na?ve Bayes)與決策樹結合,提高了對多類分類的準確度。Li等[15]提出了一種AdaBoost-BP檢測器,實現了對液壓執行器內部泄露故障的自動評估。鑒于滾動軸承不同故障時期特征信號的相似性,本文提出一種基于OHF Elman神經網絡與AdaBoost的新型集成算法,以適用早期故障信號的多樣性與魯棒性。

為了對滾動軸承不同時期的故障進行有效的診斷,本文提出整合EEMD與OHF Elman-AdaBoost算法的整體診斷框架,重點將集成算法引入故障診斷領域。鑒于Bagging算法是通過選擇不同的樣本數據訓練得到不同的弱學習器,弱學習器的個數與樣本數量有關,而AdaBoost算法則是調整樣本權重得到不同的弱學習器。相比較而言,AdaBoost算法原理更易理解且更具有適用性,因此,本文選擇以OHF Elman神經網絡作為弱回歸器,利用AdaBoost集成算法構造一個強回歸器,從而實現對滾動軸承故障診斷誤差的顯著降低。

1 問題描述

以滾動軸承振動信號作為原始數據進行故障多時期診斷。在故障診斷中,不僅需要診斷出滾動軸承所處故障時期,同時也需要準確剖析出故障源。本研究全面分析了滾動體、內圈與外圈這三種主要故障部件,因此所涵蓋的滾動軸承故障多時期診斷可轉化為多分類問題,如圖1所示。

圖1 滾動軸承故障分類圖Fig.1 Classification diagram of rolling bearing faults

由于滾動軸承振動信號的非平穩性、非線性,且摻雜不同程度的噪聲,選用EEMD進行模態分解與信號重構,從振動信號中提取有效的特征向量。在此基礎上,選擇OHF Elman神經網絡作為故障模式識別的方法,并結合AdaBoost算法集成得到強回歸器,研究對滾動軸承進行更為準確地故障診斷。所構建的技術方法框架如圖2所示。

圖2 滾動軸承故障診斷技術框架圖Fig.2 Technical framework of fault diagnosis for rolling bearing

該框架主要由故障信號提取與故障模式識別兩部分組成。故障信號提取方面,采用峭度與相關系數作為EEMD模態分解得到基本模態分量(intrinsic mode function,IMF)篩選的判斷準則。故障模式識別方面,通過Elman神經網絡的診斷誤差效果確定OHF Elman的隱含層數,從而保證兩種神經網絡的參數一致性,進而證明OHF Elman神經網絡相較于Elman神經網絡的優越性。在與AdaBoost算法集成時,采用樣本誤差閾值作為樣本權重變化的判斷準則。利用所構建的機制,最終形成的OHF Elman-AdaBoost強回歸器在全樣本數據上可實現更為準確的故障診斷水平。

OHF Elman基于Elman神經網絡增加輸出層對隱含層的反饋實現處理動態數據的能力。鑒于單個OHF Elman神經網絡無法保證對所有測試樣本都有較好的適用性,本研究為了提供更為精確地診斷滾動軸承全樣本故障,將AdaBoost集成算法引入故障診斷領域。所建立的OHF Elman-AdaBoost算法的優勢在于:

(1) 選擇以OHF Elman神經網絡作為診斷算法核心,以其獨特的反饋結構,擁有對數據的動態處理能力,特別適用于滾動軸承的故障診斷。

(2) AdaBoost算法通過集成多個弱學習器得到一個強學習器,從而擁有對樣本數據更強的學習能力。

(3) 以OHF Elman神經網絡作為弱回歸器,結合AdaBoost算法集成出新的強回歸器,從而實現對滾動軸承全樣本數據更有效的故障診斷。

2 診斷方法建模

2.1 OHF Elman神經網絡模型

人工神經網絡在理論上可以任意精度逼近任意非線性映射,無需建立精確的數學模型,一旦輸入輸出數據確定,神經網絡即可挖掘輸入輸出之間的連接關系。與BP神經網絡相比,Elman神經網絡在結構上增加了一層隱含層對輸入層的反饋,具有更高的動態性與準確性。而OHF Elman神經網絡則是在Elman神經網絡基礎上增加了輸出與隱含層的反饋,進一步提高了網絡對歷史數據的記憶功能。OHF Elman結構如圖3所示。

圖3 OHF Elman網絡結構示意圖Fig.3 Illustration of OHF Elman network structure

OHF Elman神經網絡的數學模型

x(k)=f(wI1xc(k)+wI2u(k-1))

(1)

xc(k)=αxc(k-1)+x(k-1)

(2)

yc(k)=γyc(k-1)+y(k-1)

(3)

y(k)=g(wI3x(k)+wI4yc(k))

(4)

式中:g(·)為線性函數;α,γ均為反饋增益因子,其取值范圍為(0,1);f(t)通常取為sigmoid函數,即

(5)

誤差函數為

(6)

式中:k為系統迭代步數;u為r維輸入向量;x,xc分別為n維隱含層與承接層輸出;y,yc分別為m維輸出層與承接層2的輸出;wI1,wI2,wI3,wI4分別為n×n,n×r,m×n,m×m維權值矩陣;yd為實際輸出。

2.2 原AdaBoost算法模型

AdaBoost是一種常用于提高數據訓練性能的機器學習集成算法。一般AdaBoost算法在使用的時候會選擇一種弱學習器,目前常用的有分類樹、支持向量機與神經網絡等。AdaBoost主要是通過在下一次弱學習器訓練之前適當調整訓練數據分布中的權重,來強調先前學習效果錯誤或者學習誤差較大的樣本。因此,通過集成每一步的弱學習器可以更好地對測試數據進行預測診斷[16-17]。

AdaBoost算法主要解決分類與回歸問題,分類算法目前主要有AdaBoost基本算法、AdaBoost M1與AdaBoost M2算法,分別適用于兩類分類、多類離散分類與多類連續分類問題。回歸算法主要有AdaBoost R2與AdaBoost RT算法。AdaBoost算法構建強學習器時比較靈活且結果不易發生過擬合。AdaBoost R2模型的流程如下。

弱回歸器個數:K

弱回歸器:Gj(X),j=1,…,K

輸出:強回歸器G(X)

(1) 訓練集樣本權重初始化:w1(i)=1/N,i=1,…,N

(2)j=1,…,K循環訓練

① 對樣本進行訓練并得到弱回歸器Gj(X),計算弱回歸器Gj(X)在樣本數據上的最大誤差Ej

Ej=maxi|Oji-Yji|,i=1,…,N

(7)

式中,Oji為弱回歸器Gj(X)對訓練集Xi的預測結果。

(8)

(9)

(10)

上面分別為線性誤差、平方誤差與指數誤差情況下的相對誤差計算公式。

③ 計算弱回歸器Gj(X)的誤差率ej

(11)

④ 計算弱回歸器Gj(X)的權重系數αj

(12)

⑤ 訓練集樣本權重更新

(13)

(14)

(3) 強回歸器G(X)獲得

(15)

2.3 OHF Elman-AdaBoost算法模型

為了進一步提高對滾動軸承故障多時期診斷的準確性,本研究構建了一種OHF Elman-AdaBoost算法:選擇以OHF Elman神經網絡作為弱回歸器,反復訓練OHF Elman神經網絡預測輸出,通過AdaBoost算法集成得到一個由多個OHF Elman神經網絡組成的強回歸器。

OHF Elman-AdaBoost模型流程如下。

弱回歸器個數:K

弱回歸器:OHF Elman神經網絡Cj(X),j=1,…,K

輸出:強回歸器H(X)

(1) 訓練集樣本權重初始化:w1(i)=1/N,i=1,…,N

(2)j=1,…,K循環訓練

① 初始化ej=0,對樣本Xi與wj(i)(i=1,…,N)進行弱回歸器Cj(X)訓練,并計算弱回歸器Cj(X)在各個樣本數據上的絕對誤差εi

(16)

② 計算弱回歸器Cj(X)的誤差率ej

ej=∑wj(i),εi>τ

(17)

式中,τ為預先設定的閾值,高于τ即為預測不準確。

③ 計算弱回歸器Cj(X)的權重系數αj

(18)

④ 訓練集樣本權重更新

(19)

(20)

(3) 弱回歸器Cj(X)權重系數歸一化

(21)

(4) 強回歸器H(X)獲得

(22)

OHF Elman-AdaBoost模型與AdaBoost R2模型的區別:①在計算弱回歸器的誤差率部分,OHF Elman-AdaBoost模型使用樣本的絕對預測誤差閾值進行樣本劃分,因為當閾值設置合理時,小于閾值的樣本預測誤差已經可以達到要求,只考慮預測誤差較大樣本會加快模型迭代;②在計算弱回歸器的權重系數部分,OHF Elman-AdaBoost模型借鑒了AdaBoost分類模型,通過對數函數減慢弱分類器權重系數的變化,從而得到更好的強回歸器;③在訓練樣本權重更新部分,OHF Elman-AdaBoost模型選擇使用不同的常數系數來更新,這樣會加快樣本權重的變化,也是對第二點的呼應,中和弱分類器權重變化速度,從而加快強回歸器的集成。

3 案例分析

3.1 數據集說明

本研究采用自凱斯西儲大學軸承中心的滾動軸承振動數據,即采用SKF軸承進行的0.017 78 cm,0.035 56 cm和0.053 34 cm故障直徑的實驗。分析電機轉速為1 797 r/min,采樣頻率為12 kHz的內圈、滾動體與外圈故障數據(包括正常模式),共有十種故障模式。每種故障狀態選擇120 000個樣本點,正常狀態選擇240 000個樣本點。該實驗是采用電火花加工將單點故障引入試驗軸承,0.017 78 cm,0.035 56 cm和0.053 34 cm均為故障點的直徑,現實設備的滾動軸承一旦發生單點故障,故障點的直徑會隨著故障時間的發展而增大。因此,本文選擇故障直徑為0.017 78 cm的故障狀態為故障初期,故障直徑為0.035 56 cm的故障狀態為故障中期,故障直徑為0.053 34 cm的故障狀態為故障晚期。下文中b017,b035,b053代表直徑為0.017 78 cm,0.035 56 cm和0.053 34 cm下的滾動體故障;IR017,IR035,IR053代表直徑為0.017 78 cm,0.035 56 cm和0.053 34 cm下的內圈故障;OR017,OR035,OR053代表直徑為0.017 78 cm,0.035 56 cm和0.053 34 cm下的外圈故障。

3.2 數據處理

軸承故障振動信號通常表現為周期性瞬態脈沖,為了有效地顯示時頻域內的特征,采樣信號應至少覆蓋兩個或三個周期。考慮采樣頻率和特征頻率,本文選擇2 000個樣本點代表故障振動信號的特征。因此,按時序每2 000個樣本點設定為一組輸入數據,即每種故障狀態(涵蓋滾動體、內圈、外圈的各自初期、中期、晚期故障共九種)各有60組樣本數據,正常狀態有120組樣本數據。

以2 000個樣本點正常信號的EEMD分解波形圖為例,其結果如圖4所示。從上至下,10個波形圖分別為2 000樣本點分解出的9個IMF分量和一個殘余分量R。從圖4中可見,各個IMF分量之間基本消除了模態混疊且殘余分量R單調遞減。

圖4 正常信號EEMD分解波形圖Fig.4 Decomposition waveform of normal signal

3.2.1 信號重構

由于噪聲與原時序信號完全無關與無量綱參數:峭度K對信號的沖擊成分比較敏感,因此可將IMF分量與原時序信號的相關系數與分量的峭度值作為IMF分量選擇的指標。本文選取相關系數0.1作為篩選閾值,相關系數高于0.1的IMF分量保留下來,相關系數低于0.1的IMF分量移除。不同故障模式保留的IMF分量不盡相同,正常模式選擇保留IMF1~IMF5,而OR017故障模式選擇保留IMF1與IMF2。通過計算保留下來的IMF分量峭度值發現,內圈與外圈故障保留的IMF分量峭度值較高,且在故障直徑上存在較大差異。選取峭度值較高的前兩個IMF分量進行信號重構:正常狀態選擇IMF3與IMF5;b017與b053故障狀態選擇IMF1與IMF3;b035故障狀態選擇IMF1與IMF4;IR017,IR035,IR053,OR017與OR053故障狀態選擇IMF1和IMF2;OR035故障狀態選擇IMF3和IMF4。

3.2.2 故障特征提取

選擇合適的故障特征是進行準確故障診斷的前提,通過EEMD的重構信號分析可以計算信號的特征參數,以此來判斷軸承是否發生故障與發生故障的類型。在時域內分別選擇三種有量綱參數:均值、標準差、均方根值與三種無量綱參數:偏度、峭度、裕度,見表1。

表1 三種有量綱與三種無量綱參數Tab.1 Three dimensional and three dimensionless parameters

3.3 參數設計

在模型參數設計之前,首先明確樣本數據。樣本數據的輸入由故障特征提取的六種時域參數組成,每種故障狀態下選擇60組樣本數據,正常狀態下選擇120組樣本數據。樣本數據的輸出則由[0,…,0,1,0,…,0]T向量表示(1對應的位置i表示為第i種故障模式,其余為0),正常、b017,b035,b053,IR017,IR035,IR053,OR017,OR035,OR053分別代表第1~第10種故障模式。同時將樣本數據按照5∶1的比例進行隨機分成訓練數據與測試數據。

3.3.1 OHF Elman神經網絡參數設計

OHF Elman神經網絡算法的參數主要有精度指標與隱含層神經元數。參數確定過程與結果如下:

(1) 精度指標,選取均方誤差MSE(mean square error)

(23)

式中:Ti為真實輸出數據;Oi為網絡輸出數據;s為數據維度。

(2) 神經網絡隱含層神經元數n的確定

根據式(24)首先確定隱含層神經元數在4~14。通過Elman神經網絡在不同隱含層神經元數下進行訓練,選擇測試誤差MSE較小的神經元數作為隱含層神經元數。如圖5所示,隱含層數n確定為10。

圖5 神經元個數不同的Elman神經網絡預測誤差圖Fig.5 Elman neural network prediction error graph with different number of neurons

(24)

式中:n為隱含層神經元數;r為輸入層神經元數;m為輸出層神經元數;a為1~10的常數。

3.3.2 OHF Elman-AdaBoost模型參數設計

根據OHF Elman-AdaBoost模型流程,OHF Elman-AdaBoost模型參數主要是閾值τ與弱回歸器數量K,參數確定過程與結果如下:

(1) 閾值τ的確定

閾值τ的確定直接決定了樣本的權重變化,從而決定了強回歸器對樣本的預測誤差水平。為了確定閾值τ,首先選擇弱回歸器OHF Elman神經網絡進行樣本數據訓練與測試,結果為:平均絕對誤差為0.051 1,最大絕對誤差為2,最小絕對誤差為8.427 2×10-19,因此τ確定為0.01。

(2) 弱回歸器數量K的確定

由于強回歸器的診斷效果與弱回歸器數量之間的相關性并不確定,弱回歸器數量的增加甚至可能造成過擬合。考慮到集成學習的迭代時間和有效性,首先選擇8~15作為弱回歸器的數量范圍,并分別進行OHF Elman-AdaBoost算法的測試訓練,試驗結果如圖6所示。

為了選擇弱回歸器數量,以所有測試樣本的MSE誤差之和作為判別準則,如圖6中的縱軸所示。圖6說明了弱回歸器數量對強回歸器的影響是非線性的,并且在弱回歸器數量為13時預測誤差之和較小,因而確定弱回歸器的數量K=13。

圖6 不同數量弱回歸器下的強回歸器預測誤差Fig.6 Prediction error of strong regressor under different numbers of weak regressors

3.4 模型訓練

時小虎等學者在理論上證明了OHF Elman神經網絡的收斂性與相較于Elman神經網絡的優越性,但是這一結論還是需要進行實例驗證。因此,本文在進行OHF Elman-AdaBoost模型的測試訓練之前,首先進行Elman神經網絡與OHF Elman神經網絡的對比訓練,其結果如圖7所示。

圖7 Elman與OHF Elman神經網絡預測誤差圖Fig.7 Prediction error of Elman and OHF Elman neural network

因為一組測試樣本會得到一個10×1的絕對誤差矩陣,將0.2作為誤差數據的上限,10個誤差結果均低于0.2為診斷正確,否則為診斷錯誤。根據診斷結果,Elman神經網絡的診斷準確率為90.91%,OHF Elman神經網絡的診斷準確率為91.82%。表2為OHF Elman神經網絡的部分測試樣本診斷結果。

表2 OHF Elman神經網絡的部分測試樣本診斷結果Tab.2 Diagnostic results of some test samples under OHF Elman neural network

圖7是兩種神經網絡在各個測試樣本上的診斷效果。可以看出,除去一組樣本數據的預測誤差上兩種神經網絡基本相同,在其他樣本數據上的預測誤差,OHF Elman神經網絡皆小于Elman神經網絡,證明OHF Elman神經網絡的預測性能確實優于Elman神經網絡。

在證明了OHF Elman神經網絡更適合滾動軸承的故障診斷后,進行OHF Elman-AdaBoost的訓練與測試,測試結果如圖8與圖9所示。

圖8展示了每個弱回歸器與強回歸器對滾動軸承樣本數據的平均診斷效果。如圖8所示,強回歸器的預測誤差要遠小于每一個弱回歸器的預測誤差。各個弱回歸器的平均預測誤差具體數據如表3所示。

表3 各個弱回歸器的樣本平均預測誤差MSETab.3 Average of sample prediction error MSE for each weak regressor

圖8 強回歸器與各個弱回歸器的預測誤差對比圖Fig.8 Comparison of prediction error between strong regressor and each weak regressor

圖9展示了OHF Elman-AdaBoost算法得到的強回歸器與OHF Elman神經網絡弱回歸器在各個測試樣本上的預測誤差效果。可以看出,除去在一組樣本上強回歸器與弱回歸器預測誤差幾乎相等,強回歸器基本上在所有測試樣本上都取得了更好的預測效果。

圖9 強、弱回歸器在測試樣本上的預測誤差MSE對比圖Fig.9 Comparison of prediction errors between strong and weak regressors on test samples

表4計算了表3中弱回歸器的預測誤差均值,并與強回歸器預測誤差進行了對比。結果可見,通過OHF Elman-AdaBoost算法得到的強回歸器在弱回歸器OHF Elman神經網絡的基礎上降低了61.7%的預測誤差。

表4 強、弱回歸器預測誤差對比表Tab.4 Comparison of prediction error between strong and weak regressors

表5對比了所有弱回歸器對于不同故障時期的平均預測誤差與強回歸器對于不同故障時期的預測誤差,強回歸器在不同故障時期的預測效果均優于弱回歸器,其中在故障晚期預測效果提升最為明顯,同時也在故障初期與故障中期上降低了約25%與66.5%的預測誤差。在不同的故障部位上,強回歸器也擁有較好的診斷效果,其中,對于內圈故障診斷效果最好,外圈故障次之,滾動體故障診斷誤差最大。

表5 強、弱回歸器下不同故障時期診斷預測誤差對比圖Tab.5 Comparison of prediction errors between strong regressor and weak regressor in different fault periods

表6則進一步細分對比了所有弱回歸器對于不同故障不同時期的平均預測誤差與強回歸器對于不同故障不同時期的預測誤差。如表所示,強回歸器對于故障初期與故障中期的診斷誤差主要集中于滾動體故障診斷,對于故障晚期診斷均具有非常好的效果。對于內圈故障的三種時期均具有較低的診斷誤差,對于外圈的故障診斷誤差主要來源于故障中期,而對于滾動體的故障初期與中期診斷效果較差。與弱回歸器診斷結果相比,強回歸器同時提高了在三種故障的不同時期診斷精度。由此可以驗證,所提出的OHF Elman-AdaBoost算法在滾動軸承故障多時期診斷的精確性方面有著很好的性能優勢。

表6 強、弱回歸器下三種故障的不同時期診斷預測誤差對比圖Tab.6 Comparison of prediction errors for different periods of three faults under strong regressor and weak regressor

4 結 論

針對滾動軸承不同故障時期故障信息的相似性與振動信號的非平穩性和非線性,本文提出了一種基于EEMD與OHF Elman-AdaBoost算法的研究框架。通過EEMD進行信號分解與重構,提取了有效的故障特征。然后選擇OHF Elman神經網絡作為滾動軸承故障診斷的主要方法,通過其獨特的反饋結構,提高了對不同故障時期信息的動態處理功能。然后結合AdaBoost算法設計了一種新穎的OHF Elman-AdaBoost算法模型。經算例分析,OHF Elman-AdaBoost算法可以有效地對滾動軸承進行故障多時期診斷,同時提高了對全樣本數據的預測精度。未來將進一步考慮優化神經網絡結構與AdaBoost算法參數設計,以期在滾動軸承故障多時期診斷上提供更為精確的解決方案。

猜你喜歡
故障診斷故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
故障一點通
江淮車故障3例
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 国产精品成人啪精品视频| 精品中文字幕一区在线| 亚洲成av人无码综合在线观看| 日韩无码真实干出血视频| 国产欧美中文字幕| 国产亚洲日韩av在线| 日韩欧美中文亚洲高清在线| 亚洲色图综合在线| 国产最爽的乱婬视频国语对白| 国产亚洲精品自在线| 国产99精品久久| 狠狠躁天天躁夜夜躁婷婷| 91娇喘视频| 欧美天堂在线| 国产午夜精品一区二区三区软件| 久久国产精品影院| 免费A级毛片无码无遮挡| 国产精品视频导航| 欧美怡红院视频一区二区三区| 波多野结衣无码视频在线观看| 国产原创自拍不卡第一页| 亚洲国产欧美国产综合久久| 亚洲av无码久久无遮挡| 97在线观看视频免费| 国产成人综合亚洲网址| 最新痴汉在线无码AV| 成人av专区精品无码国产| 久久天天躁狠狠躁夜夜2020一| 人人妻人人澡人人爽欧美一区| 久久黄色一级视频| 国产91全国探花系列在线播放| 国产成人在线无码免费视频| 一区二区三区成人| 999精品色在线观看| 久久精品丝袜高跟鞋| 精品久久久无码专区中文字幕| 亚洲高清在线天堂精品| 亚洲国产精品人久久电影| 国产欧美日韩在线一区| 欧美日韩国产在线人| 国产精品免费久久久久影院无码| 国产小视频在线高清播放 | 国产精品9| 国产九九精品视频| 97人人做人人爽香蕉精品| 国产成人av一区二区三区| 久久久久青草大香线综合精品| 免费国产黄线在线观看| julia中文字幕久久亚洲| 99精品在线看| 国产性生大片免费观看性欧美| 国产二级毛片| 精品伊人久久久大香线蕉欧美| 亚瑟天堂久久一区二区影院| 999精品在线视频| 欧美成在线视频| 尤物精品视频一区二区三区| 亚洲男人的天堂久久香蕉网| 亚洲AV无码久久天堂| 动漫精品中文字幕无码| 视频国产精品丝袜第一页| 国产69精品久久久久孕妇大杂乱| 国产在线第二页| 亚洲第一在线播放| 欧美另类视频一区二区三区| 欧美一级一级做性视频| 久久精品亚洲热综合一区二区| 天天躁日日躁狠狠躁中文字幕| 强乱中文字幕在线播放不卡| 国产导航在线| 欧美一级片在线| 又大又硬又爽免费视频| 久久国产亚洲偷自| 国产sm重味一区二区三区| 亚洲天堂久久新| 99热在线只有精品| 55夜色66夜色国产精品视频| 又黄又爽视频好爽视频| 欧美日韩中文字幕二区三区| 久久免费精品琪琪| 美女内射视频WWW网站午夜| 凹凸精品免费精品视频|