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

基于極點對稱模態分解和概率神經網絡的軸承故障診斷

2017-03-09 06:49:27張淑清徐劍濤姜安琦李軍鋒宿新爽姜萬錄
中國機械工程 2017年4期
關鍵詞:故障診斷模態故障

張淑清 徐劍濤 姜安琦 李軍鋒 宿新爽 姜萬錄

1.燕山大學電氣工程學院,秦皇島,0660042.中南大學信息工程學院,長沙,4100063.燕山大學機械工程學院,秦皇島,066004

基于極點對稱模態分解和概率神經網絡的軸承故障診斷

張淑清1徐劍濤1姜安琦2李軍鋒1宿新爽1姜萬錄3

1.燕山大學電氣工程學院,秦皇島,0660042.中南大學信息工程學院,長沙,4100063.燕山大學機械工程學院,秦皇島,066004

針對復雜非線性的滾動軸承系統,提出了極點對稱模態分解(ESMD)和概率神經網絡(PNN)相結合的滾動軸承故障診斷方法。ESMD將固有模態函數的定義進行擴充,采用內部極點對稱直接插值的方法替代外部包絡線插值,引入最優的自適應全局曲線(AGM)的概念優化分解的趨勢線,并由此確定最佳的模態分解次數。PNN是一種基于核函數逼近的神經網絡分類器,將指數函數引入神經網絡用來替代S型激活函數并進行重新構造,突出體現了梯度最速下降法的概念,減少實際和預測的輸出函數之間的誤差。通過對經驗模態分解(EMD)、屏蔽經驗模態分解(MEMD)和ESMD方法進行信號仿真分解對比,以及采用ESMD和PNN對故障數據進行處理,結果表明,該方法能夠更加有效地對故障信號進行識別。

滾動軸承;極點對稱模態分解;概率神經網絡;故障診斷

0 引言

滾動軸承是旋轉機械的重要組成部分,其振動信號是各種噪聲和軸承的振動信號的疊加,具有非線性和非平穩性特征。處理方法包括傅里葉變換、小波變換以及自適應信號處理方法經驗模態分解(empirical mode decomposition,EMD)等[1-3]。傅里葉變換和小波變換等時頻信號分析方法都是由線性疊加原理發展而來的,對振動信號的非線性非平穩信號不適用[4]。EMD將采集的滾動軸承的振動信號進行包絡分解,獲得包含各個主分量的固有模態函數,從而體現出振動信號的特征分量信號[5]。然而其在信號分解過程中出現的模態混疊以及不能自適應地解決欠包絡和過包絡問題,限制了EMD處理信號的能力[6]。自適應信號處理方法屏蔽經驗模態分解(masking empirical mode decomposition,MEMD) 使EMD法的混疊問題得到有效的解決[7]。然而,由于MEMD在處理信號的過程中添加了掩蔽信號,并且掩蔽信號的選擇和處理具有不確定性,在信號的處理過程中不能完全綜合所添加的信號,造成了其性能缺陷。

文獻[8]研究極點對稱模態分解算法,借鑒經驗模態分解,采用內部極點對稱直接插值的方法替代外部包絡線插值,引入最優的自適應全局曲線(adaptive global mean,AGM)的概念來優化分解的趨勢線,并由此確定最佳模態分解次數,取得了較好效果。

概率神經網絡(probabilistic neural network,PNN)與BP(back propagation)神經網絡相比具有一個相對較快的訓練過程,這是因為預先設定的代表訓練集的規模保證了最大化和訓練樣本可以隨時添加或刪除而無需廣泛地再訓[8]。

概率神經網絡是一種基于核函數逼近的神經網絡分類器,將指數函數引入神經網絡用來替代S型激活函數并進行重新構造,突出體現了梯度最速下降法的概念,減小實際和預測的輸出函數之間的誤差。PNN在處理傳統意義上復雜的、非線性的和不精確的問題時具有獨特的優勢。本文采用PNN神經網絡故障數據進行處理,更加有效地對故障信號進行識別。

1 信號的ESMD分解

1.1 經驗模態分解

經驗模態分解的原理可以簡單表述為

(1)

其中,R為殘余量;cj為信號經過EMD分解所產生的固有模態函數(intrinsic mode function,IMF)[8],n為信號經EMD分解所產生的IMF的個數,j表示第j個IMF分量。

1.2 極點對稱模態分解

經驗模態分解中,所分解出的IMF分量不是完全正交的,如果信號中的極點分布不均勻或者極點數過于少,EMD分解就不能進行并且會出現模態混疊現象。極點對稱模態分解的方法根據信號的特點,引入最優自適應全局曲線概念優化分解的趨勢線,并由此確定最佳模態分解次數[9-11]。

ESMD重新定義了IMF:①局部的極大值點和極小值點是不同的,相鄰且相等的極值點作為一個極值點加入到信號分解中,極大值必須是正數,極小值必須是負數;②廣義上,IMF分量應該幾乎是包絡對稱或者極點對稱的。

ESMD并不將信號分解到最后剩余最多的一個極值點,相反,它允許殘留成分包含一定數量的極值點。這一過程的優勢一方面在于分解的殘余分量可以反映整體數據的變化趨勢,可以理解為一種最優的自適應全局曲線(adaptive global mean,AGM);另一方面在于可以優化最小二乘法意義上的AGM曲線來確定最優的篩選次數。它提供了一個很好的自適應方法的數據擬合,優于一般的最小二乘法。

不同于EMD中構建兩個不同的外部包絡線,ESMD采用多次對極值點連線的重點進行內部曲線插值。根據以往實驗,一般在數據處理中,兩次內部曲線插值即可達到標準[9-11]。數據處理過程中,篩選次數的確定參考均值曲線的原理,采用最佳篩選次數的分解條件。

1.3 極點模態分解算法實現過程

(1) 用Ei(i=1,2,…,n)分別將信號X所有的極大值點和極小值點進行標記。

(2)用線段相互連接所有的極值點,并標記這些線段的中點為Fi(i=1,2,…,n-1)。

(3)通過一定的方法分別在左右兩個邊界添加F0和Fn。

(4)根據上述1+n個點構造p個內插值曲線L1,L2,…,Lp(p=1,2,…),并且計算它們的均值:

L*=(L1+L2+…+Lp)/p

(2)

(5)令X-L*,并重復上述4個步驟,直到

|L*|≤ε

(3)

式中,ε為允許的誤差。

或者篩選次數達到最大值K。此時,得到第一個模態函數M1。

(6)對于剩余部分X-M1重復上述5個步驟,可以依次得到M2,M3,…,直到最后的剩余分量R含有不超過約定的極點個數。

(7)在有限整數區間[Kmin,Kmax]改變最大篩選次數K的取值,重復上述6個步驟。然后計算X-R的方差σ2并作關于σ/σ0和K的曲線圖,其中σ0為信號X的標準差。

(8)找到區間[Kmin,Kmax]中對應于σ/σ0最小值時K的最小值K0,用K0重復前6個步驟,得到最佳的ESMD分解出的各IMF分量。與此同時,得到的R即為最優的AGM曲線。

通常選擇ε=0.001σ0,并且選擇比率υ=σ/σ0去反映AGM曲線的最優程度。

此時原始數據和AGM曲線可分別用X={xi}和R={ri}表示,i=1,2,…,n。

定義數據的總體均值為

(4)

數據的方差為

(5)

AGM曲線的方差為

(6)

經過分解, 原始的時間序列X可表示為

(7)

利用ESMD將時間序列X分解成一系列IMF和一個趨勢項。

1.4 ESMD與EMD、MEMD的對比

EMD將信號分解成包含不同主頻分量的IMF分量,然而,當信號存在異常干擾或信號發生間斷跳躍時,這些突變信號的頻率往往要高于原信號的頻率,出現模態混疊,無法將不同的主頻分量有效地分離到不同的IMF分量中,使得一個固有模態函數分量中出現多個相近的主頻分量,所以信號內部的特征無法真正通過EMD分解出的IMF分量表現出來[12]。模態混疊的出現導致EMD對信號的分解不能進行有效分析。屏蔽經驗模態分解(MEMD)法由于掩蔽信號的不確定性造成其對信號的分解也會出現不確定性。

取兩個等幅值的單頻信號的疊加形式如下:

x1(t)=sin2πt+sin1.2πt

(8)

圖1 x1(t)的EMD分解Fig.1 EMD of x1(t)

對該信號x1(t)進行EMD分解,如圖1所示,信號x1(t)的頻譜如圖1a所示,f1=1 Hz,f2=0.6 Hz;圖1b是EMD分解后的IMF1的頻譜圖。從圖1b中可以看出,IMF1分量并不只包含一個主頻率分量,它不僅包含主要的1 Hz高頻成分而且含有部分0.6 Hz的低頻成分,充分體現了經驗模態分解過程中出現的模態混疊問題。在實際應用中,信號形式十分復雜,相鄰的頻率成分大小比較接近,難以區分,這就使得EMD分解不可避免地存在模態混疊問題。

采用ESMD對信號x1(t)進行處理后的頻譜圖見圖2。通過對比圖2和圖1,可以看出ESMD分解后的IMF1分量中只包含單一的主頻分量,不存在ESMD分解所存在的模態混疊問題。

圖2 x1(t)的ESMD分解Fig.2 ESMD of x1(t)

為了綜合對比EMD、MEMD和ESMD的效果,取信號:

x2(t)=sin(2πf1fsn)+sin(2πf2fsn)

(9)

其中,fs=1024 Hz,f1=177.6 Hz,f2=100 Hz。信號經EMD、MEMD和ESMD處理結果如圖3所示。

圖3 x2(t)的EMD、MEMD、ESMD對比Fig.3 Comparison among EMD, MEMD and ESMDof x2(t)

由圖3三種方法的對比可以看出,MEMD由于掩蔽信號的加入,在分解后會出現信號的失真,造成物理意義上的分解錯誤,顯然ESMD分解方法很好地解決了模態混疊問題。

2 概率神經網絡

典型的PNN神經網絡如圖4所示。它一般包括4個組成部分,即輸入層、模式層、疊加層和輸出層。輸入層是第一層神經元。每個輸入神經元在訓練/測試數據集代表一個單獨的屬性,輸入的數目與數據集中的屬性數相等。輸入數據的值乘以適當的權重,并被傳送到模式層。輸出層通常只包含一個類的最后一層,因為通常只有一個輸出被要求。在訓練階段,目標是確定最準確的權重分配給連接線。在這個階段,輸出被反復計算,并與訓練/測試數據集所產生的優選輸出的結果進行比較[13-15]。

圖4 概率神經網絡圖Fig.4 Schematic of probabilistic neural network

(10)

(11)

PNN具有如下優良特性: ①網絡學習訓練的過程簡單,能較快地達到訓練標準;②PNN神經網絡對數據的分類準確,只要訓練樣本足夠大,可以達到非常準確的分類效果,并且噪聲對其分類效果影響低;③PNN可以用線性計算方法計算非線性數據,并能夠以任意精度逼近[16]。

3 采用ESMD和PNN的軸承故障診斷流程

(1)對振動信號進行ESMD分解。

(2)分別將得到的IMF分量與原信號的相關系數進行計算并排序,取前n個IMF分量用來表征振動信號(n的取值根據具體的振動信號分解而定)。相關表達式如下:

D(X)=E(X-E(X))2

(12)

cov(X,Y)=E((X-E(X))(Y-E(Y)))

(13)

(14)

式中,D(X)為信號X的標準差;E(X)為信號X的均值;cov(X,Y)為信號X和信號Y的協方差;ρXY為信號X和信號Y的相關性系數。

(3)為了便于對數據進行PNN訓練和分類,將上述所得相關性大的IMF分量進行能量的計算[17]:

(15)

i=1,2,…,n

構造出向量T,并對所得向量T進行整體歸一化處理:

T=(E1,E2,…,En)

(16)

(17)

(18)

向量T′即為歸一化后的能量特征向量。

(4)將T′作為特征向量輸入到PNN中進行訓練。

4 故障診斷實例

本實驗采用美國凱斯西儲大學(Case Western Reserve University)的滾動軸承故障數據。實驗裝置如圖5所示,主要以SKF6205型深溝球軸承數據為仿真數據。

圖5 軸承故障診斷實驗儀器Fig.5 Experimental instrument of bearingfault diagnosis

實驗中,分別隨機選取風扇端以采樣頻率12 kHz采集滾動軸承的正常、滾動體故障、外圈故障和內圈故障4種類型的故障信號共240組,每組選取1500個點。240組數據中,200組數據作為訓練樣本,40組作為測試樣本。

首先對上述4種軸承振動數據進行ESMD分解,并計算分解后的包含各個單主頻的固有模態函數與信號的相關系數,根據相關系數選取了7組經ESMD處理過的內蘊模態函數作為PNN神經網絡訓練的特征分量。然后對這些特征分量進行整體歸一化,分別標記滾動軸承正常、滾動體故障、外圈故障和內圈故障4種類型的故障數據為1、2、3、4,然后輸入到PNN神經網絡中進行分類器的分類。ESMD分解振動信號如圖6所示。

圖6 ESMD分解振動信號Fig.6 The experimental results after ESMD

PNN神經網絡訓練中初始的SPREAD系數ηS=1。隨著SPREAD系數的逐漸減小,訓練效果會逐步變好。當ηS=0.1時,PNN神經網絡訓練結果見圖7。當ηS=0.08時,PNN神經網絡訓練結果見圖8。從圖中可知,實驗結果達到了預期效果。

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖7 ηS=0.1時ESMD分解后PNN訓練結果Fig.7 The result of using ESMD and PNN to diagnosis the rolling bearing fault under ηS=0.1

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖8 ηS=0.08時ESMD分解后PNN訓練結果Fig.8 The result of using ESMD and PNN to diagnosis the rolling bearing fault under ηS=0.08

為了說明本文方法的優越性,同樣信號采用EMD分解、MEMD分解和PNN訓練,PNN訓練結果如圖9~圖12所示。

圖9和圖10所示均為在信號分解的方法中采用了EMD后的效果。與圖7、圖8對比可以明顯看出,ESMD方法處理過的數據相對于EMD處理過的數據在初始條件和變化條件下,數據的收斂狀態與速度都要好很多。

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖9 ηS=0.1時EMD分解后PNN訓練結果Fig.9 The result of using EMD and PNN to diagnosis the rolling bearing fault under ηS=0.1

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖10 ηS=0.08時EMD分解后PNN訓練結果Fig.10 The result of using EMD and PNN to diagnosis the rolling bearing fault under ηS=0.08

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖11 ηS=0.1時MEMD分解后PNN訓練結果Fig.11 The result of using MEMD and PNN to diagnosis the rolling bearing fault under ηS=0.1

(a)PNN訓練后的效果 (b)PNN訓練后的誤差圖12 ηS=0.08時MEMD分解后PNN訓練結果Fig.12 The result of using MEMD and PNN todiagnosis the rolling bearing fault under ηS=0.08

圖11與圖12所示均為在信號分解的方法中采用了MEMD后的效果。與圖9、圖10對比發現,MEMD處理信號后經PNN訓練效果要好于EMD方法,在信號的初始條件以及變化條件下,其收斂速度都要高于EMD方法。但是與圖7、圖8相比,ESMD對信號的處理優于MEMD方法。

由上述分析可知,隨著SPREAD系數的逐漸減小,經ESMD分解的信號收斂較快,分類結果準確、誤差小,證明了ESMD方法的優越性。

5 結束語

本文通過極點模態分解將機械振動信號分解成多個IMF函數,以IMF分量作為信號的特征向量,采用概率神經網絡對樣本數據進行訓練,對測試數據進行辨識,構成了旋轉機械故障診斷的新方法。通過實驗證明該方法可以有效地區分滾動軸承的正常、內圈、外圈、滾動體故障等故障狀態,與傳統的經驗模態分解等方法相比,在提取信號的特征向量方面具有分解徹底、殘余量可以表征信號振動趨勢等優點,便于在故障初期發現故障并及時消除。實驗結果表明,采用ESMD和PNN相結合的方法能夠更有效地對不同故障類型信號進行診斷與識別。

[1]LIYongbo,XUMinqiang,WANGRixin,etal.AFaultDiagnosisSchemeforRollingBearingBasedonLocalMeanDecompositionandImprovedMultiscaleFuzzyEntropy[J].JournalofSoundandVibration, 2016,360:277-299.

[2] 趙志宏. 基于振動信號的機械故障特征提取與診斷研究[D]. 北京:北京交通大學,2012.ZHAOZhihong.ResearchonVibrationSignalBasedMachineryFaultFeatureExtractionandDiagnosis[D].Beijing:BeijingJiaotongUniversity, 2012.

[3] 鐘先友. 旋轉機械故障診斷的時頻分析方法及其應用研究[D]. 武漢:武漢科技大學, 2014.ZHONGXianyou.ResearchonTime-frequencyAnalysisMethodsandItsApplicationstoRotatingMachineryFaultDiagnosis[D].Wuhan:WuhanUniversityofScienceandTechnology, 2014.

[4]GILLESJ.EmpiricalWaveletTransform[J].IEEETransactionsonSignalProcessing, 2013, 61(16): 3999-4010.

[5] 黃建,胡曉光,鞏玉楠.基于經驗模態分解的高壓斷路器機械故障診斷方法[J]. 中國電機工程學報,2011,31(12):108-113.HUANGJian,HUXiaoguang,GONGYunan.MachineryFaultDiagnosisofHighVoltageCircuitBreakerBasedonEmpiricalModeDecomposition[J].ProceedingsoftheCSEE, 2011,31(12):108-113.

[6] 時培明, 李庚, 韓東穎. 基于改進EMD的旋轉機械耦合故障診斷方法研究[J].中國機械工程, 2013, 24(17): 2367-2372.SHIPeiming,LIGeng,HANDongying.StudyonCouplingFaultsofRotaryMachineryDiagnosisMethodBasedonImprovedEMD[J].ChinaMechanicalEngineering, 2013, 24(17): 2367-2372.

[7]DEERINGR,KAISERJF.TheUseofaMaskingSignaltoImproveEmpiricalModeDecomposition[C]//Acoustics,Speech,andSignalProcessing, 2005.Pennsylvania, 2005: 485-488.

[8] 全學海, 丁宣浩, 蔣英春. 基于EMD和概率神經網絡的說話人識別[J]. 桂林電子科技大學學報, 2012, 30(2): 108-112.QUANXuehai,DINGXuanhao,JIANGYingchun.SpeakerRecognitionBasedonEMDandProbabilisticNeuralNetworks[J].JournalofGuilinUniversityofElectronicTechnology, 2012,30(2): 108-112.

[9]WANGJinliang,LIZongjun.Extreme-pointSymmetricModeDecompositionMethodforDataAnalysis[J].AdvancesinAdaptiveDataAnalysis, 2013,5(3):1350015.

[10]LIHuifeng,WANGJinliang,LIZongjun.ApplicationofESMDMethodtoAir-seaFluxInvestigation[J].InternationalJournalofGeosciences, 2013, 4(5): 8-11.

[11] 王金良, 李宗軍. 可用于氣候數據分析的ESMD方法[J]. 氣候變化研究快報, 2014(3): 1-5.WANGJinliang,LIZongjun.TheESMDMethodforClimateDataAnalysis[J].ClimateChangeResearchLetters, 2014(3): 1-5.

[12] 包紅燕. 基于MEMD和條件熵相空間重構的滾動軸承故障診斷[D]. 秦皇島:燕山大學, 2014.BAOHongyan.RollingBearingFaultDiagnosisBasedonMaskingEmpiricalModeDecompositionandPhaseSpaceReconstructionofConditionalEntroy[D].Qinhuangdao:YanshanUniversity, 2014.

[13]CHENXianyue,ZHOUJianzhong,XIAOHan.FaultDiagnosisBasedonComprehensiveGeometricCharacteristicandProbabilityNeuralNetwork[J].AppliedMathematicsandComputation, 2014, 230(3) : 542-554.

[14]WANGChangqing,ZHOUJianzhong,QINGHui.FaultDiagnosisBasedonPulseCoupledNeuralNetworkandProbabilityNeuralNetwork[J].ExpertSystemswithApplications, 2011, 38 (11):14207-14313.

[15] 劉鳳龍, 宋藝. 基于EMD與PNN的機械故障檢測[J]. 計算機應用與軟件, 2010, 27(9): 237-239.LIUFenglong,SONGYi.MachineryFaultDiagnosisBasedonEMDandPNN[J].ComputerApplicationsandSoftware, 2010,27(9): 237-239.

[16] 孟宗,胡猛,谷偉明,等. 基于LMD多尺度熵和概率神經網絡的滾動軸承故障診斷方法[J]. 中國機械工程, 2016, 27(4): 433-437.MENGZong,HUMeng,GUWeiming,etal.RollingBearingFaultDiagnosisMethodBasedonLMDMulti-scaleEntropyandProbabilisticNeuralNetwork[J].ChinaMechanicalEngineering, 2016,27(4): 433-437.

[17] 肖韜,袁興中,唐清華,等. 基于概率神經網絡的城市湖泊生態系統健康評價研究[J]. 環境科學學報, 2013, 33(11): 3166-3172.XIAOTao,YUANXingzhong,TANGQinghua,etal.InvestigationofHealthAssessmentforUrbanLakesSystemBasedonProbabilisticNeuralNetworks[J].ActaScientiaeCircumstantiae, 2013, 33(11): 3166-3172.

(編輯 王旻玥)

Fault Diagnosis of Bearings Based on Extreme-point Symmetric Mode Decomposition and Probabilistic Neural Network

ZHANG Shuqing1XU Jiantao1JIANG Anqi2LI Junfeng1SU Xinshuang1JIANG Wanlu3

1.Institute of Electrical Engineering,Yanshan University,Qinhuangdao,Hebei,066004 2.School of Information Engineering,Central South University,Changsha,410006 3.College of Mechanical Engineering,Yanshan University,Qinhuangdao,Hebei,066004

Aimed at the complex non-linear rolling bearing systems, a new method combining ESMD and PNN was introduced for bearing fault diagnosis. ESMD expanded the definition of the intrinsic mode function, and changed the external envelope interpolation to internal pole symmetric interpolation. An idea of adaptive global mean(AGM) was used to optimize the last remaining modal, thus to determine the optimal number of decomposition. PNN was a neural network classifier based on kernel function approximation. The exponential functions were introduced to the neural network to replace the S type activation functions and to reconstruct the functions, representing the notion of gradient steepest descent method prominently, and reducing the errors between the actual and predicted output functions. Through the decomposing comparison to the simulation signals among empirical mode decomposition(EMD), making empirical mode decomposition(MEMD) and ESMD, and the diagnosis to the bearing vibration data by ESMD and PNN shows that the new method introduced may diagnose the bearing faults more effectively.

rolling bearing; extreme-point symmetric mode decomposition(ESMD); probabilistic neural network(PNN); fault diagnosis

2016-04-01

國家自然科學基金資助項目(51475405,61077071);河北省自然科學基金資助項目(F2016203496,F2015203413)

TH17

10.3969/j.issn.1004-132X.2017.04.009

張淑清,女,1966年生。燕山大學電氣工程學院教授、博士研究生導師。主要研究方向為弱信號檢測、智能信號處理、故障診斷。發表論文50余篇。E-mail:zhshq-yd@163.com。 徐劍濤,男,1992年生。燕山大學電氣工程學院碩士研究生。姜安琦,女,1995年生。中南大學信息科學與工程學院本科生。李軍峰,男,1994年生。燕山大學電氣工程學院碩士研究生。宿新爽,女,1993年生。燕山大學電氣工程學院碩士研究生。姜萬錄,男,1964年生。燕山大學機械工程學院教授、博士研究生導師。

猜你喜歡
故障診斷模態故障
故障一點通
奔馳R320車ABS、ESP故障燈異常點亮
國內多模態教學研究回顧與展望
因果圖定性分析法及其在故障診斷中的應用
故障一點通
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 欧美a级完整在线观看| 波多野结衣二区| 国产精品丝袜在线| 日韩黄色精品| 国产亚洲欧美在线中文bt天堂| 国产乱子伦一区二区=| 久久精品亚洲热综合一区二区| 亚洲人成电影在线播放| 欧美第一页在线| 国产性生大片免费观看性欧美| 亚洲日韩精品欧美中文字幕 | 日韩精品一区二区三区大桥未久| 久久综合九色综合97网| 91在线中文| 久久无码av三级| 国产精品美女自慰喷水| 亚洲国产成人久久精品软件| 国产91全国探花系列在线播放 | 亚洲中文字幕23页在线| 免费人成在线观看成人片| 国产精品无码久久久久久| 欧美一区中文字幕| 国产性猛交XXXX免费看| 午夜精品影院| 亚洲国产系列| 99在线免费播放| 亚洲综合第一区| 在线观看精品自拍视频| 精品免费在线视频| 干中文字幕| 亚洲aaa视频| 久久香蕉国产线| 在线观看国产精品一区| 国产欧美精品午夜在线播放| 日本日韩欧美| 国产成人精品一区二区不卡 | 久久亚洲国产视频| 亚洲VA中文字幕| 国产精品成人啪精品视频| 国产精品太粉嫩高中在线观看| 欧美日韩第三页| 内射人妻无套中出无码| 在线观看91精品国产剧情免费| 九九视频免费看| 日韩在线第三页| 综合色天天| 色亚洲激情综合精品无码视频| 国产精品99在线观看| 四虎成人在线视频| 亚洲精品成人片在线观看| 99精品一区二区免费视频| 国产精品综合色区在线观看| 亚洲男人天堂网址| 国产精品人莉莉成在线播放| 青青操国产| 国产91视频观看| 奇米精品一区二区三区在线观看| 91精品免费高清在线| 亚洲国产精品一区二区第一页免| 欧美精品1区| 国产黄在线观看| 亚洲欧美在线综合一区二区三区| 国产精品熟女亚洲AV麻豆| 欧美日韩在线国产| 国产在线专区| 国产午夜无码专区喷水| 经典三级久久| 亚洲天堂.com| 亚洲男人天堂2020| 国产精品亚洲а∨天堂免下载| 最新国产在线| 免费在线视频a| 国产三级国产精品国产普男人| 国产18在线| 亚洲国产成人无码AV在线影院L| 亚瑟天堂久久一区二区影院| 9cao视频精品| 青草视频在线观看国产| 国产亚洲精久久久久久久91| 国产97区一区二区三区无码| 玖玖精品视频在线观看| 精品人妻系列无码专区久久|