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

基于改進的多元離群檢測方法的風機齒輪箱早期故障診斷

2016-08-09 06:35:47顧煜炯賈子文任玉亭
中國機械工程 2016年14期

顧煜炯 賈子文 王 瑞 任玉亭

1.華北電力大學新能源電力系統國家重點實驗室,北京,1022062.國華能源投資有限公司,北京,100007

?

基于改進的多元離群檢測方法的風機齒輪箱早期故障診斷

顧煜炯1賈子文1王瑞1任玉亭2

1.華北電力大學新能源電力系統國家重點實驗室,北京,1022062.國華能源投資有限公司,北京,100007

摘要:針對風電機組運行工況波動性以及機組早期故障特征不易提取的特點,提出一種基于改進的多元離群監測方法來實現風機齒輪箱故障的早期診斷。運用階比重采樣方法對原始振動信號進行預處理,并對處理結果進行量綱一因子分析;通過馬氏距離建立風電齒輪箱的早期故障識別模型;利用多元線性回歸改進多元離群檢測算法進行實際數據的分析計算。結果表明,該方法較原始方法能夠更早地察覺出風電齒輪箱早期故障。

關鍵詞:階比重采樣;量綱一因子分析;多元線性回歸;多元離群檢測

0引言

風電場通常建設在地處偏遠、交通閉塞、環境因素變化劇烈的區域,這使得風電機組的健康運行遭受極大挑戰[1]。風機齒輪箱作為整個系統的傳動機構,內部結構緊湊,部件之間耦合性較強,在運行過程中長期受到交變載荷與沖擊載荷作用,容易造成齒輪點蝕、磨損等故障。同時,齒輪箱故障維修過程較為復雜,維修時間長,長時間的停機維修給風場業主帶來很大的經濟損失[2]。因此,進行風電機組齒輪箱早期故障診斷的研究,尋找故障早期信號特征,確定故障模式,在故障還沒有發展到嚴重程度時及時排除安全隱患,對保證機組正常運行和提高風場經濟效益具有重要意義。

目前,針對齒輪箱故障診斷的方法有很多,如經驗模態分解(EMD)與隱馬爾科夫模型[3]、小波分析[4]、神經網絡技術[5]等方法。這些方法雖然能夠實現齒輪箱故障的診斷,但都主要面對設備的中晚期故障問題,對齒輪箱早期故障特征提取和分析的過程較少。風電機組齒輪箱早期故障診斷具有如下特性:①風電機組運行工況具有波動性、間歇性的特點,造成現場采集的振動信號具有明顯的非線性和非平穩特征,需要有對應的方法進行處理,否則在很大程度上會影響振動信號特征提取的效果;②機組齒輪箱早期故障信號特征表現不明顯,需要找到對早期故障數據敏感的參數并將其作為征兆,實現機組齒輪箱的早期診斷。針對這一現狀,本文提出一種改進的多元離群檢測方法對風電齒輪箱進行故障的早期診斷。通過對齒輪箱故障特征參數進行隸屬度劃分,建立多元統計數理模型,運用多元離群檢測方法實現風電機組齒輪箱故障的在線診斷。

1階比重采樣角域信號特征的獲取

有量綱幅域參數受到工作負載、轉速等運行條件的影響,使得故障趨勢判斷存在誤差。峰值、翹度、波形等量綱一幅域參數指標對幅值能量變化不敏感,即與設備運行條件關系不大,但對設備故障較為敏感,并且計算過程簡單,容易應用到實際工程中[6]。郭厚明等[7]運用量綱一因子實現了礦用低速重載齒輪的故障診斷;岑少起等[8]通過對監測參數進行量綱一因子轉換,實現了滑動軸承的動力學分析。目前,量綱一因子在旋轉機械早期故障診斷領域的應用較少,所以,筆者利用量綱一因子的特性,通過提取齒輪箱振動數據信息,計算齒輪箱振動數據的因子值,并將其作為故障早期診斷的特征。

受風電機組工況影響,齒輪箱振動信號存在較為明顯的非線性特征。所以本文先通過階比重采樣技術實現非平穩時域信號向平穩角域信號的轉化,建立了角域信號量綱一因子變化趨勢指標。

1.1階比重采樣技術

風電機組運行工況的特殊性導致機組振動信號具有明顯的非線性、非平穩性特征,運用傳統的信號分析方法很難提取出反映設備運行狀況的振動特征。階比重采樣通過對振動信號和轉速信號進行同步采集,并結合插值計算的方法保證設備每一轉采樣的點數相等,實現非平穩時域信號向平穩角域信號的等價轉換[9]。

假設機組角度編碼器碼盤孔數為n,則孔間角度Δθ=360/n(°),角度編碼器計數脈沖每秒脈沖數序列記為{m1,m2,…,mk},k為時間序列序號??梢缘贸鲲L電機組第i秒的平均轉速(r/min):

Ri=60miΔθ

(1)

則可依據每秒平均轉速得出累計轉角曲線。假設等時間間隔{t1,t2,…,tk}采樣下的秒級時域序列為{x1,x2,…,xk},其中,xk為時域振動幅值。則時域序列轉化成角域序列的公式如下:

(2)

將等時間時域信號{x1,x2,…,xk}轉化為等時間間隔采樣下的角域信號{y1,y2,…,yk},因為機組在此時間不一定是恒定轉速運行的,所以此角域信號非嚴格意義上的等角度間隔角域信號序列,為解決此問題,采用數值插值方法對其進行進一步計算。

1.2齒輪箱早期故障特征提取

新的量綱一因子值趨勢分析技術主要解決風電機組齒輪箱早期故障振動信號特征值提取困難的問題。通過對齒輪箱等角域信號量綱一因子值的計算,分析量綱一數值變化規律來反映風電齒輪箱的故障發展趨勢[10]:機組運行正常時,各個量綱一因子值沒有明顯的變化趨勢;機組出現早期故障時,各個量綱一因子的值會出現不同的單調性變化,通過這些因子值的變化趨勢以及數值變化程度,實現機組早期故障模式的甄別。

筆者為克服傳統有量綱幅域參數與能量有關,或與能量無關但是只能定性分析特性的缺點,基于時域波動統計分析,引入對能量不敏感的量綱一幅域參數,將其作為故障特征參數。各因子的具體描述如下:

(1)翹度Kf。Kf與旋轉部件的尺寸、設計參數、運行工況等無關,對沖擊信號反應敏感,其數值隨故障發展程度加深而增大,且變化較為明顯。Kf的計算公式為

(3)

式中,xi為數據數列中的第i個數據的值;σ為數據標準差;N為數據個數。

(2)波形裕度CL與偏態因子SK。這兩個因子對振動信號形狀和趨勢的細微波動變化敏感,大量實驗證明,這兩個因子可以作為判斷旋轉機械早期故障的指標,其計算公式分別為

(4)

(5)

式中,Xmax為數據數列中的最大值。

(3)重復性因子Rf。該指標適合對波形的重復性進行定量分析,隨故障的發展其波形重復性變差,因子值發生變化。將數據按整周期截取成s段,每段t個數,則可得數據{x11,x12,…,x1t;x21,x22,…,x2t;…;xs1,xs2,…,xst},每個時間段的數據可以看作一個角域序列,則對應的平均重復性波形定義為

(6)

計算重復波形平均差分值,運用鏈碼技術對波形差分符號進行編輯,規則如圖1所示,圖中,0表示波形下降,差分值為負值;1表示波形不變,差分值為0;2表示波形上升,差分值為正值。對平均重復波形編碼進行比較,得到編碼中不相同點的個數k,則重復性因子為

Rf=k/t

(7)

圖1 波形差分編碼規則

(4)相似性因子Ff。此因子涉及到分形理論內容,它運用盒維數的概念反映故障信號早期變化情況。首先對角域信號進行符合其運算要求的標準化處理:

(8)

式中,xγ(ti)為ti時刻的幅值;max|xγ(j)|表示記錄數據中幅值絕對值最大的值;K為比例放大因子,一般取整數,需根據經驗與實際情況而定。

分形盒維數以單元方格的方式對信號波形進行覆蓋,要求覆蓋單元有較強相似性,且整個信號波形有嚴格的自相似性。設F是實數集合Rn中任意非空有限子集,記N(F,φ)為最大邊長φ能覆蓋F集合的最小數字,則F的盒維數計算表達式為

(9)

相似性因子為

Ff=dimBF

(10)

(5)跳躍性因子Jf。此因子主要反映振動波形的幅度調制,通過式(8)對原始角域信號進行標準化處理,對處理后的數據進行整周期截取,此過程類似重復性因子對數據的分段處理,計算方差值:

(11)

1.3實際案例分析

以滄州某風電場1.5MW風電機組為例,經調查,2013年2月初的第3次全場巡檢發現某一風電機組一級行星太陽輪發生較嚴重的磨損。從風場中繼室數據庫調取該機組整個2月內齒輪箱一級行星輪位置測點全部振動數據,對原始振動數據進行階比重采樣計算,將時域信號轉化成平穩的等角度角域信號后,以天為單位對2月份29天的歷史數據進行傳統量綱一因子和新量綱一因子的分析計算。

由圖2可以看出:在機組故障初期,新的量綱一因子數值隨故障發展有較為明顯的變化。對比圖2、圖3可以看出:新量綱一因子較傳統量綱一因子的變化趨勢更為明顯,數據波動小,說明新的量綱一因子對機組早期故障信號的微弱波動變化敏感。所以通過量綱一因子對角域信號的特征提取,能夠準確分析出風機齒輪箱早期故障。

圖2 新量鋼一翹度Kf和波形裕度CL變化趨勢圖

圖3 傳統量綱一方差值D和峰峰值Xmax變化趨勢圖

2多元離群檢測方法的改進

多元故障特征離群檢測可以有效地將多元時序特征融合成綜合評價指標,通過綜合指標的差異反映設備故障的嚴重程度,實現對設備早期故障的等級劃分。本文選取基于距離的離群檢測方法來解決融合多元故障特征的風電機組早期故障模式預警問題。基于距離的離群點檢測算法,一方面避免了數據分布模型不確定時,檢測精確度較低的問題;另一方面,當處理屬性較多的空間數據集時,基于距離的離群點檢測算法比基于密度的離群點檢測方法效率更高,這對于將數據異常檢測用于旋轉機械設備的在線監測與數據分析至關重要。

離群檢測中,距離求解方法很多,這里采用馬氏(Mahalanobis)距離實現對機組齒輪箱運行狀態的監測。傳統的馬氏距離方法在進行樣本數據與被檢測數據距離比較的時候,沒有對影響故障模式的各個特征參數進行隸屬度的劃分。因為故障數據對各個故障特征波動的影響不同,同時,每個故障特征對數據變化的敏感程度也有差異,所以,找到各個故障模式下故障特征之間的隸屬關系,可以提高對故障模式變化顯著的因子的貢獻度。

為解決計算馬氏距離時各影響因素隸屬度的問題,這里采用多元線性回歸的方法。

2.1多元線性回歸方程

多元線性回歸[11]方程是描述因變量Y的平均值或期望值如何依賴于自變量X1、X2、…、Xp的方程。多元線性回歸方程形式為

E(Y)=b0+b1X1+b2X2+…+bpXp

(12)

其中,bi(i=0,1,…,p)為待定參數,表示假定其他變量不變,僅Xi每變動一個單位時,Y的平均變量。

總體參數b0、b1、…、bp是未知的,必須利用樣本觀測值去估計它們。估計的多元線性回歸方程為

(13)

2.2參數最小二乘估計

(14)

式中,yi為在第i個時刻實際觀察的因變量數值。

計算出各個自變量的待定參數后,要對回歸方程擬合優度、回歸方程顯著性和回歸系數顯著性進行檢驗。

2.3改進馬氏距離的計算

馬氏距離[12]的計算公式為

(15)

在實際應用中,Xi、Xj均可量化處理出對應的數值,例如信號能量、平均溫度等,作為多元線性回歸方程的因變量。由式(14)計算得出Xi、Xj向量中各個數據的參數,假設有p個特征參數,則記參數隸屬度序列為{b0、b1、b2、…、bp}。這樣就對參與計算的向量Xi、Xj中每一個數據完成隸屬度劃分,則改進的馬氏距離為

(16)

3齒輪箱早期故障模型建立

3.1故障特征隸屬度計算

首先要確定風機齒輪箱故障的分析周期。假定檢驗周期數為m,故障模式對應的故障特征元數目為n,則風電機組某故障模式的待檢測多元故障特征矩陣經過標準化處理后可表示為

3.2齒輪箱故障特征值的確定

已知風電機組多元故障特征指標中的傳動鏈振動數據為6種量綱一幅域參數,分別對風機齒輪箱正常工況和不同故障情況下的數據進行收集,選取某機組不同時間段內的歷史數據。為保證計算結果隨數值的波動性變化不大,對不同工況、不同階段的故障特征值取均值,確定機組齒輪箱故障特征值:

式中,x0ρ為融合隸屬度標準參考樣本的特征矩陣;xiρ為考慮了特征隸屬度的待測樣本矩陣。

3.3基于馬氏距離的多元離群檢測的故障模式確定

基于馬氏距離的多元故障特征離群檢測公式為

4案例分析

依舊以某風電場1.5 MW風電機組一級行星太陽輪磨損為例,選取2013年1月份的31天為研究時間段,以每天同一時間記錄下來的10 s振動數據為分析樣本,來挖掘該時間段內機組齒輪故障發生過程中的多元故障特征的演變情況。其中,參與分析的數據包括:信號能量E(單位為J)、翹度Kf、波形裕度CL、偏態因子SK、重復性因子Rf、相似性因子Ff、跳躍性因子Jf、溫度異常率φ。

齒輪箱故障發生初期,機組信號數據故障特征并不明顯,各量綱一因子的值變化較為緩慢,而且每個因子值的變化趨勢因故障模式不同的會略有差異,不易統一劃定出機組故障的判斷特征。故障發生時,油溫綜合了各故障模式的機理特性,即故障引起的齒輪箱零件間振動與摩擦程度加大,導致油溫發生改變。故將齒輪箱油溫異常率作為進行診斷的閾值,同時將各個量綱一因子值變化趨勢作為故障模式特征,實現風電機組齒輪箱的早期故障診斷,這里將溫度異常率的閾值設為0.1。

通過對數據進行多元線性回歸分析計算,獲得的風電機組一級太陽輪磨損故障與多元故障特征指標間的量化隸屬關系式為

E=19 382.43+178.12Kf+235.65CL+453.07SK+231.54Rf+603.11Ff+112.45Jf

(17)

對式(17)進行擬合優度評定、回歸方程顯著性檢驗(F檢驗)和回歸系數顯著性檢驗(t檢驗)。

計算擬合優度評定系數為

R2=SSR/SST=0.9640

式中,SSR為總離差平方和;SST為回歸平方和。

調整后的擬合優度判定系數為

n=31p=6

4.1F檢驗

(i=1,2,…,6)

由表1可知,F=107.18,顯著水平α=0.05,查表得F0.05(6,24)=2.51,明顯小于F。所以拒絕H0,接受H1,即故障特征隸屬度值不全為零,模型線性關系在95%置信水平下顯著成立。

表1 方差分析表

4.2t檢驗

因為模型自由度為6,這里只對第一個自由度t1進行假設檢驗:

H0∶b1=0;H1∶b1≠0。

給定一個顯著性水平α=0.05,查t分布表,得到臨界值Ta=0.5000。因為t1=2.144>Ta,所以拒絕H0,接受備擇假設,即回歸系統b1≠0。對于其他回歸系數bi(i=2,3,…,6),用上述同樣方法可得出各回歸系數是顯著不為0的。

將同一臺風機2011年5月~9月的正常數據作為樣本數據,得出正常工況下6個量綱一因子閾值集合和隸屬度集合:

{Kf,CL,SK,Rf,Ff,Jf}=

{0.525,5.211,1.103,0.017,12.342,5.178}

{b1,b2,b3,b4,b5,b6}=

{162.76,213.77,412.32,221.21,569.01,100.99}

對機組1月份31天數據進行馬氏距離計算,結果如表2所示。

表2 1月份故障數據與樣本數據馬氏距離

由表2可得到,該機組2013年1月實際數據與標準樣本數據的馬氏距離從第5天開始,直到月末26天均屬于離群點,因此可判斷該機組從1月5日開始逐漸檢測出一級太陽輪發生磨損故障。

最后,計算多元離群檢測因子:

ODF=W/A

式中,W為待測對象的離群點數目;A為所有待測數據總數。

從5日到31日,因子值由1%增長到24%,說明機組故障程度正在不斷加深。

表3 傳統方法計算馬氏距離

通過對比可以發現,傳統算法在第11天發現機組齒輪箱出現一級太陽輪磨損故障,比改進方法晚了5天,說明改進方法能夠更快速有效地識別風電機組齒輪箱早期的故障。

5結論

(1)應用軟件算法實現了階比重采樣計算,將非平穩時域信號轉化成平穩角域信號,不僅保證了后續計算結果的準確性,同時在實際工程應用中,節省了硬件費用。

(2)利用量綱一因子對設備監測數據波動變化敏感的特點,將因子趨勢變化作為風機齒輪箱早期故障的數據基礎,為后續的故障模式識別打下堅實基礎。

(3)改進的多元離群檢測算法對影響機組故障的各個特征參數進行了隸屬度分析計算,使得基于距離的多元離群檢測更有說服力,計算結果更為可靠,有很好的工程應用前景。

參考文獻:

[1]劉德順,戴巨川,胡燕平,等.現代大型風電機組現狀與發展趨勢[J].中國機械工程,2013,24(1):125-134.

LiuDeshun,DaiJuchuan,HuYanping,etal.StatusandDevelopmentTrendsofModernLarge-scaleWindTurbines[J].ChinaMechanicalEngineering, 2013,24(1):125-134.

[2]張亮.風力發電機組齒輪箱早期故障診斷方法研究[D]. 大連:大連理工大學,2010.

[3]何邵燦,高宏力,許明恒.基于隱馬爾科夫模型的機床部件故障預警技術[J].機械設計與制造,2012 (8):159-161.

HeShaocan,GaoHongli,XuMingheng.ComponentsofMachineTools’FailureWarningBasedonHiddenMarkovModel[J].MachineryDesignandManufacture, 2012 (8):159-161.

[4]李蓉,于德介,陳向民,等.基于階次分析與循環平穩解調的齒輪箱復合故障診斷方法[J]. 中國機械工程,2013,24(10):1320-1327.

LiRong,YuDejie,ChenXiangmin,etal.ACompoundFaultDiagnosisMethodforGearboxBasedonOrderTrackingandCyclostationaryDemodulation[J].ChinaMechanicalEngineering, 2013,24(10):1320-1327.

[5]祁麗婉,梁庚,童國煒.基于果蠅算法優化BP神經網絡的齒輪箱故障診斷[J]. 電網與清潔能源,2014,30(9):31-42.

QiLiwan,LiangGeng,TongGuowei.AGearboxDiagnosisMethodBasedonFruitFlyOptimizationAlgorithmtoOptimizetheBPNeuralNetwork[J].PowerSystemandCleanEnergy, 2014,30(9):31-42.

[6]馬懷祥,徐明新,金玨.應用油液分析法和振動法診斷往復機械故障[J].石家莊鐵道學院學報,1998,11(2):39-43.

MaHuaixiang,XuMingxin,JinJue.BreakdownDiagnosisofReciprocatingMachinerybyMeansofOilAnalysisandVibrationAnalysis[J].JournalofShijiazhuangRailwayInstitute, 1998,11(2):39-43.

[7]郭厚明,行志剛. 量綱一參數在礦用低速重載齒輪故障診斷中的應用[J]. 煤炭科學技術,2006,34(8):28-31.

GuoHouming,XingZhigang.DimensionlessParametersAppliedtoFaultDiagnosisofMineLowSpeedHeavyLoadedGear[J].CoalScienceandTechnology, 2006,34(8):28-31.

[8]岑少起,于衛東,李瑞珍. 動靜壓滑動軸承基本方程和特性系數量綱一化研究[J].機械傳動,2008,32(3):4-7.

CenShaoqi,YuWeidong,LiRuizhen.TheStudyofDimensionlessBasicEquationandCharacteristicsCoefficientsofHydrodynamicHydrostaticBearings[J].JournalofMechanicalTransmission, 2008,32(3):4-7.

[9]HochmanD,SadokM.TheoryofSynchronousAveraging[C]//Proceedingsofthe2004IEEEAerospaceConference.BigSky,Montana,USA,2004:3636-3653.

[10]WonlenK,MansonG,FiellerNJ.DamageDetectionUsingOutlierAnalysis[J].JournalofSoundandVibration, 2000, 229(3):647-667.

[11]扈靜,劉明周,龔任波,等. 基于主成分分析的汽車操縱力舒適性分析與評價[J]. 中國機械工程,2011,22(20):2456-2459.

HuJing,LiuMingzhou,GongRenbo,etal.EvaluationofVehicleOperatingForceComfortBasedonPrincipalComponentAnalysis[J].ChinaMechanicalEngineering, 2011,22(20):2456-2459.

[12]申志剛,何寧,李亮. 高速硬銑削加工刀具磨損監測研究[J].中國機械工程,2009,20(13):1582-1586.

ShenZhigang,HeNing,LiLiang.MonitoringofToolWearinHardMillingProcess[J].ChinaMechanicalEngineering, 2009,20(13):1582-1586.

(編輯張洋)

收稿日期:2015-09-10

基金項目:神華集團科技創新項目(GTKJ-12-02);華能集團科學技術項目(HNKJ-H27)

中圖分類號:TK83

DOI:10.3969/j.issn.1004-132X.2016.14.011

作者簡介:顧煜炯,男,1968年生。華北電力大學能源動力與機械工程學院教授、博士研究生導師。主要研究方向為電站大型旋轉機械振動監測與故障診斷。發表論文80余篇。賈子文,男,1986年生。華北電力大學能源動力與機械工程學院博士研究生。王瑞,男,1989年生。華北電力大學能源動力與機械工程學院碩士研究生。任玉亭,男,1966年生。國華能源投資有限公司副總工程師。

Early Fault Diagnosis for Wind Turbine Gearbox Based on Improved Multivariate Outlier Detection

Gu Yujiong1Jia Ziwen1Wang Rui1Ren Yuting2

1.State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources,North China Electric Power University,Beijing,102206 2.Guohua Energy Investment Limited Company,Beijing,100007

Abstract:An improved method of multivariate outlier detection was used in early fault diagnosis for wind turbine gearboxes, which might extract the early fault features under fluctuation working conditions. First, the primitive vibration signals were preprocessed by order resampling, and the processed results were analyzed by dimensionless parameter analysis. Second, a model of early fault diagnosis was created based on Mahalanobis distance, which was used for turbine gearboxes. Finally, the actual data was analyzed by the method of multivariate outlier detection, which was improved by multiple linear regression. The results show that the new method may detect the gearbox faults earlier than original one.

Key words:order resampling; dimensionless parameter analysis; multiple linear regression(MLR); multivariate outlier detection

主站蜘蛛池模板: 国产精品亚洲а∨天堂免下载| 国产亚洲精久久久久久久91| 一本久道久久综合多人| 国产免费怡红院视频| 日本午夜影院| 99久久精品国产麻豆婷婷| 91av成人日本不卡三区| 亚洲人成网线在线播放va| 色悠久久久久久久综合网伊人| 无码内射在线| 亚洲综合二区| 亚洲欧美不卡视频| 亚洲爱婷婷色69堂| 一本大道香蕉久中文在线播放| 一级毛片免费高清视频| 国产精品视频a| 国产日韩精品一区在线不卡| 精品国产乱码久久久久久一区二区| 国产日韩欧美精品区性色| 播五月综合| 成人精品免费视频| 亚洲国产欧美中日韩成人综合视频| 亚洲国产精品国自产拍A| 国产伦精品一区二区三区视频优播| 在线观看国产精品一区| 日韩欧美国产综合| 全色黄大色大片免费久久老太| www.狠狠| 亚洲浓毛av| 一区二区三区四区在线| 91国语视频| 亚洲中文字幕久久无码精品A| 99伊人精品| 狠狠色综合网| 久久精品无码一区二区日韩免费| 91精品情国产情侣高潮对白蜜| AV色爱天堂网| 欧美三级日韩三级| 国产女人在线| 久草热视频在线| 国产丝袜无码一区二区视频| 日韩免费毛片视频| 日本手机在线视频| 97色婷婷成人综合在线观看| 欧美中出一区二区| 人妻丝袜无码视频| 2020国产在线视精品在| 综合人妻久久一区二区精品 | 中文字幕精品一区二区三区视频| 亚洲成a人片在线观看88| 日韩一二三区视频精品| 国产呦视频免费视频在线观看 | 免费视频在线2021入口| 亚洲视频欧美不卡| 精品无码国产一区二区三区AV| 18禁色诱爆乳网站| 国产真实二区一区在线亚洲| 91小视频在线观看免费版高清| 欧美综合在线观看| 国产一区二区色淫影院| 欧美激情福利| 午夜国产精品视频| 99国产精品免费观看视频| 久久精品无码专区免费| 1024国产在线| 老色鬼久久亚洲AV综合| 人妻丝袜无码视频| 亚洲国产欧美中日韩成人综合视频| 91极品美女高潮叫床在线观看| 国产精品自在线天天看片| 亚洲精品波多野结衣| 国产一级无码不卡视频| 亚洲男人的天堂在线| 欧美专区在线观看| 蝴蝶伊人久久中文娱乐网| 欧美亚洲香蕉| 日韩AV无码免费一二三区| 亚洲精品色AV无码看| 免费一级毛片| 大香伊人久久| 久久先锋资源| 久久天天躁狠狠躁夜夜躁|