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

基于變分模態(tài)分解和魯棒性獨立成分分析的內(nèi)燃機缸蓋振動信號分離

2018-05-02 07:22:30姚家馳錢思沖張冠軍
中國機械工程 2018年8期
關(guān)鍵詞:模態(tài)振動信號

姚家馳 向 陽 錢思沖 張冠軍

1.武漢理工大學(xué)能源與動力工程學(xué)院,武漢,4300632.武漢理工大學(xué)船舶動力系統(tǒng)運用技術(shù)交通行業(yè)重點實驗室,武漢,430063

0 引言

內(nèi)燃機的缸蓋振動信號中蘊含著內(nèi)燃機工作狀態(tài)的重要信息,對其進行研究可以判斷內(nèi)燃機的燃燒異常、氣門間隙異常、活塞敲擊異常等故障[1-4],因此一直是學(xué)者們研究的熱點。內(nèi)燃機缸蓋振動信號由許多激勵源產(chǎn)生的信號混疊而成,利用時頻分析方法不能很好地區(qū)分各獨立源信號[5]。近年來,盲源分離方法廣泛地應(yīng)用于生物醫(yī)學(xué)信號處理、語音分離、圖像處理等領(lǐng)域,許多學(xué)者用盲源分離方法對內(nèi)燃機振動信號進行分離。LIU等[6-7]用盲最小均方差算法從機體振動信號中分離出燃燒信號和活塞敲擊信號。康斌等[8]基于Fixed-point ICA算法從缸蓋振動信號中分離出燃燒壓力激勵信號和氣門落座激勵信號。韓春楊等[9]用多通道盲最小均方差與縮減盲源方法從內(nèi)燃機表面混合振動信號中分離出燃燒信號和活塞撞擊缸體信號。但盲源分離方法要求觀測信號數(shù)目與源信號數(shù)目相同,需要測量內(nèi)燃機多個通道的振動信號,在實際工程測試中受造價和安裝條件的限制,往往只能使用較少的傳感器且難以確定源信號的數(shù)目,該方法在實際工程應(yīng)用中受到限制。后來,又有學(xué)者提出了內(nèi)燃機振動信號的單通道分離方法。DU等[10]利用經(jīng)驗?zāi)B(tài)分解和獨立分量分析相結(jié)合的方法從內(nèi)燃機單通道機體振動信號中分離識別出燃燒信號和活塞敲擊信號,但是經(jīng)驗?zāi)B(tài)分解方法缺乏嚴密的數(shù)學(xué)推導(dǎo),產(chǎn)生的端點效應(yīng)和模態(tài)混疊問題會嚴重影響振動信號的分離效果[11-13]。目前,在缸蓋振動信號中,燃燒信號成分和活塞敲擊信號成分都發(fā)生在上止點附近,在時域和頻域上混疊嚴重,對其進行分離和識別一直是研究的熱點和難點。

針對這些問題,本文提出采用變分模態(tài)分解(variational mode decomposition, VMD)和魯棒性獨立成分分析(robust independent component analysis, RobustICA)相結(jié)合的方法從單通道缸蓋振動信號中分離出燃燒信號成分和活塞敲擊信號成分,并通過對內(nèi)燃機不同試驗工況下的單通道缸蓋振動信號進行分離和識別,驗證該方法的有效性。

1 VMD和RobustICA算法的基本原理

1.1 VMD基本原理

變分模態(tài)分解算法[14]是一種信號分解算法,采用非遞歸的分解模式,能夠很好地對多分量信號進行分解。假設(shè)將一個信號通過變分模態(tài)分解算法分解為K個變分模態(tài)分量,則其相應(yīng)的約束變分問題為

(1)

{uk}={u1,u2,…,uK}

{ωk}={ω1,ω2,…,ωK}

在對該約束變分問題進行求解時,需要引入二次懲罰因子α和拉格朗日乘法算子λ將其轉(zhuǎn)變?yōu)闊o約束問題,擴展的拉格朗日表達式:

(2)

(2)n←n+1;

(3)對k=1,2,…,K,執(zhí)行

(3)

(4)

(5)

式中,τ為噪聲容限參數(shù)。

(4)判斷是否滿足收斂條件:

(6)

式中,c為預(yù)先設(shè)置的迭代停止值。

若滿足條件,則停止程序運行;否則繼續(xù)執(zhí)行步驟(2)。

式(3)可寫為

(7)

采用Parseval/Plancherel傅里葉等距變換將式(7)變換到頻域進行求解:

(8)

對其進行求解后可得到

(9)

同理,可以將式(4)等價為

(10)

對其進行求解后可得

(11)

1.2 RobustICA基本原理

獨立分量分析算法的實質(zhì)是在假設(shè)源信號統(tǒng)計獨立的基礎(chǔ)上,在不知道源信號及混合矩陣任何信息的情況下,試圖將一組隨機變量表示成統(tǒng)計獨立的變量的線性組合。

設(shè)X(t)是一組由若干隱含變量產(chǎn)生的p維觀測信號,其中i=1,2,…,p表示觀測信號的序數(shù),獨立分量分析要從p維觀測信號X(t)中找出隱含在其中的變量Y(t),通常為解決多變量數(shù)據(jù)分析問題,將其簡化為一個線性問題進行處理,即

(12)

i=1,2,…,pj=1,2,…,m

寫成矩陣形式為Y=WX,即

FastICA算法是應(yīng)用比較廣泛的一種獨立分量分析算法[15],隨后ZARZOSO等[16]學(xué)者提出了魯棒性更好和收斂速度更快的RobustICA算法,RobustICA算法是一種基于峭度和最優(yōu)步長的迭代算法,將一組信號經(jīng)過RobustICA算法處理后,可以得到各個相互獨立的信號成分。

1.3 VMD-RobustICA方法的計算流程

本文提出的VMD-RobustICA方法可從單通道混合信號中恢復(fù)出各獨立源信號,其計算流程如圖1所示。首先用VMD算法對采集的單通道信號進行分解,然后將分解得到的各個窄帶模態(tài)分量和采集的單通道信號組成一個新的信號組,然后用RobustICA算法對其進行解耦,并用組合模態(tài)函數(shù)法(combined mode function,CMF)對解耦得到的在時域和頻域上相似性較高的分量進行組合[17],最后結(jié)合頻譜分析、連續(xù)小波變換(continuous wavelet transform,CWT)、相干函數(shù)法及倒拖試驗對得到的各獨立源信號進行識別,其中相干函數(shù)定義為

圖1 VMD-RobustICA方法的計算流程Fig.1 The calculation process of VMD-RobustICA method

(13)

式中,pf(ω)和px(ω)分別為輸入函數(shù)f(t)和輸出函數(shù)x(t)的自功率譜密度;pfx(ω)為輸入輸出函數(shù)的互功率譜密度,相干函數(shù)值在0~1之間。

2 內(nèi)燃機試驗

試驗對象為WP10-240國Ⅲ型六缸四沖程水冷直列式內(nèi)燃機,點火順序為1-5-3-6-2-4。試驗臺架由WP10-240型內(nèi)燃機、德國Siemens 1PL6交流電機、倒拖控制臺及其控制配件等組成。在試驗過程中,以WP10-240型內(nèi)燃機第6號缸為研究對象,用LC0158T型加速度傳感器測量內(nèi)燃機缸蓋的單通道振動信號,用缸壓傳感器Kistler7013C和單通道電荷放大器5018A1000采集內(nèi)燃機的缸壓信號,用SM-12-100型磁電式傳感器采集上止點信號,在進行倒拖試驗時,將加速度傳感器布置在第6號缸活塞敲擊處,測量活塞敲擊振動信號,試驗測量系統(tǒng)如圖2所示,試驗工況如表1所示。

圖2 試驗測量系統(tǒng)Fig.2 Test measurement system

序號試驗工況類型轉(zhuǎn)速(r/min)負荷百分數(shù)(%)1正常運轉(zhuǎn)100002倒拖100003正常運轉(zhuǎn)2100254倒拖210025

3 缸蓋振動信號的分離與識別

在低轉(zhuǎn)速工況下,內(nèi)燃機的各振動激勵源信號相對容易被分離,在高轉(zhuǎn)速高負荷工況下,內(nèi)燃機的各振動激勵源信號中會包含有更多其他干擾成分。本文首先采用VMD-RobustICA方法對內(nèi)燃機在1 000 r/min空載工況下的缸蓋振動信號進行分離和識別,然后對內(nèi)燃機在額定轉(zhuǎn)速高轉(zhuǎn)速2 100 r/min、25%負荷工況下的缸蓋振動信號進行分離和識別。

3.1 內(nèi)燃機在1 000 r/min空載工況下缸蓋振動信號的分離與識別

內(nèi)燃機試驗工況為1 000 r/min空載,采樣頻率為25.6 kHz,通過試驗測得內(nèi)燃機一個工作循環(huán)的缸壓pc和缸蓋振動加速度信號a如圖3所示。

圖3 一個工作循環(huán)的缸壓和缸蓋振動信號Fig.3 A working cycle of cylinder pressure and cylinder head vibration signals

為了更好地對缸蓋振動信號進行分離,減少在測量信號的過程中產(chǎn)生的隨機誤差成分,需要對采集的缸蓋振動信號進行消除趨勢項及滑動平均等預(yù)處理,經(jīng)過預(yù)處理后的信號如圖4所示。

首先對預(yù)處理后的單通道缸蓋振動信號進行變分模態(tài)分解,在進行變分模態(tài)分解之前,需要通過觀察對比各個模態(tài)的中心頻率來確定最佳的模態(tài)數(shù)K值,通過計算得到的結(jié)果如表2所示。

表2 模態(tài)數(shù)K和中心頻率f

從表2中可看出,當模態(tài)數(shù)K取8時,有兩個變分模態(tài)分量的中心頻率分別為5 014 Hz和5 594 Hz,相距較近,因此可認為出現(xiàn)了過分解,故最佳模態(tài)數(shù)K應(yīng)取7。

將預(yù)處理后的缸蓋振動信號經(jīng)過變分模態(tài)分解處理后可得到7個變分模態(tài)分量,由于得到的7個變分模態(tài)分量之間不總是相互獨立的,因此需要進一步地采用RobustICA算法提取其獨立成分。將這7個變分模態(tài)分量和預(yù)處理后的缸蓋振動信號組成一個新的信號組,用RobustICA算法提取獨立成分,通過計算得到的結(jié)果如圖5所示。

圖5 RobustICA計算結(jié)果Fig.5 The calculation results of RobustICA

通過對圖5中的計算結(jié)果進行分析后可知,分量IC1和分量IC2在時域波形上比較相似,并進一步地分析其頻譜后發(fā)現(xiàn),它們在頻域上也比較相似,因此采用組合模態(tài)函數(shù)法將分量IC1和分量IC2組合為CMF12,SCMF12=SIC1+SIC2,組合分量CMF12的時域波形、頻譜和時頻圖如圖6所示。同時,通過分析發(fā)現(xiàn)分量IC3可能為活塞敲擊信號,對其進行FFT分析和連續(xù)小波時頻分析,結(jié)果如圖7所示。

從圖6中可知,組合分量CMF12時域波形的幅值在140°CA和380°CA左右變化較大,根據(jù)內(nèi)燃機的先驗知識,內(nèi)燃機的發(fā)火順序為1-5-3-6-2-4,第5號缸和第6號缸的發(fā)火角度分別在140°CA和380°CA左右,兩缸工作間隔為240°CA,同時從頻譜圖中可以看出,組合分量CMF12的頻率成分主要集中在4 350 Hz,結(jié)合缸壓和缸蓋振動的相干函數(shù)(圖8),在該頻率成分附近,缸壓和缸蓋振動的相干性很好,缸壓變化主要由缸內(nèi)燃燒引起,通過缸壓可以計算出燃燒信號[18],這里主要對分離得到的分量進行定性判斷,并進一步地結(jié)合時頻圖6c可知,在4 350 Hz附近,140°CA和380°CA左右的頻率能量值較大,并且在380°CA左右的頻率能量值要大于140°CA左右的頻率能量值,這是因為測量的是第6號缸的缸蓋振動信號,第5號缸是第6號缸的相鄰缸,燃燒產(chǎn)生的振動信號會傳遞到第6號缸,但是要比第6號缸由于燃燒引起的缸蓋振動信號小,因此可以判斷組合分量CMF12為燃燒信號。

(a)時域波形

(b)FFT

(c)CWT圖6 CMF12的時域波形、FFT和CWT圖Fig.6 The time domain waveform, FFT and CWT of CMF12

(a)時域波形

(b)FFT

(c)CWT圖7 IC3的時域波形、FFT和CWT圖Fig.7 The time domain waveform, FFT and CWT of IC3

圖8 1 000 r/min空載工況下的缸壓和缸蓋振動的相干函數(shù)Fig.8 The coherence function of cylinder pressure and cylinder head vibration of 1 000 r/min no-load condition

根據(jù)圖7可知,分量IC3的時域波形的幅值在380°CA左右變化較大,同時在140°CA左右也有一定的變化,并且分量IC3的頻率成分主要集中在1 150 Hz,根據(jù)缸壓和缸蓋振動的相干函數(shù)(圖8),在該頻率成分附近,缸壓和缸蓋振動的相干性不好,同時由圖7c可看出,在1 150 Hz附近,380°CA左右的頻率能量值很大,且在140°CA左右也有一定的頻率能量,內(nèi)燃機的發(fā)火順序1-5-3-6-2-4,與活塞敲擊缸壁的時刻一致,將倒拖試驗測得的活塞敲擊振動信號頻譜與分量IC3的頻譜進行對比,如圖9所示,可以發(fā)現(xiàn)活塞敲擊振動頻譜與分量IC3的頻譜基本吻合,因此可以判斷分量IC3主要為活塞敲擊信號。但是在2 000 Hz附近,活塞敲擊振動頻譜較大,這可能是由于在倒拖試驗中內(nèi)燃機的其他運動部件產(chǎn)生的振動所致,還有待深入研究。

圖9 1 000 r/min倒拖工況下的活塞敲擊振動和分量IC3的頻譜Fig.9 The spectrum of IC3 and piston slap vibration of 1 000 r/min drag condition

3.2 內(nèi)燃機在2 100 r/min、25%負荷工況下缸蓋振動信號的分離與識別

內(nèi)燃機在2 100 r/min、25%負荷工況下的一個工作循環(huán)的缸壓和缸蓋振動信號如圖10所示。

圖10 缸壓和缸蓋振動信號Fig.10 Cylinder pressure and cylinder head vibration signals

對采集的單通道缸蓋振動信號進行消除趨勢項及滑動平均等預(yù)處理,然后對預(yù)處理后的信號進行分解,通過觀察對比各個模態(tài)的中心頻率后發(fā)現(xiàn),當模態(tài)數(shù)K取為11時,有兩個模態(tài)的中心頻率為4 018 Hz和4 720 Hz,相距較近且出現(xiàn)了過分解,因此最佳的模態(tài)數(shù)取為10。

將通過變分模態(tài)分解算法計算得到的10個變分模態(tài)分量和預(yù)處理后的缸蓋振動信號組成一個新的信號組,用RobustICA算法提取其獨立成分,對得到的各個獨立分量進行分析,并用組合模態(tài)函數(shù)法將時域和頻域相似性較高的分量進行組合。通過分析后發(fā)現(xiàn)分量IC1和組合分量CMF23(SCMF23=SIC2+SIC3)可能為內(nèi)燃機的燃燒信號和活塞敲擊信號,下面結(jié)合頻譜分析、連續(xù)小波時頻分析、相干函數(shù)法、倒拖工況下的活塞敲擊信號以及內(nèi)燃機的先驗知識對分離得到的結(jié)果進行進一步識別。分量IC1和組合分量CMF23的時域波形、FFT和CWT圖分別見圖11和圖12。

從圖11可知,分量IC1的時域波形在380°CA處變化較大,經(jīng)過FFT頻譜分析后發(fā)現(xiàn)其頻率主要集中在4 250 Hz,由圖13可知,在該頻率成分附近缸壓和缸蓋振動的相干性較好,缸壓變化主要由缸內(nèi)燃燒引起,內(nèi)燃機的發(fā)火順序為1-5-3-6-2-4,第6號缸的發(fā)火角度在380°CA處,與時頻圖上在該位置處的能量幅值較大一致,因此可判斷分量IC1為燃燒信號。

根據(jù)圖12可知,組合分量CMF23的頻率成分主要集中在1 500 Hz,從CWT時頻圖中可看出在380°CA處的能量幅值要大于在140°CA處的能量幅值,與內(nèi)燃機的發(fā)火順序1-5-3-6-2-4相對應(yīng),即5號缸和6號缸的發(fā)火角度分別在140°CA和380°CA處,其工作間隔為240°CA,并進一步地將倒拖工況下測得的活塞敲擊振動信號頻譜和組合分量CMF23的頻譜進行對比,從圖14中可看出,活塞敲擊振動信號頻譜和組合分量CMF23的頻譜基本吻合,因此可判斷組合分量CMF23主要為活塞敲擊信號。但是活塞敲擊振動信號頻譜在500 Hz和2 500 Hz附近仍有其他頻率成分,這可能是由于在倒拖工況時其他零部件運動時產(chǎn)生的,還有待進一步研究。

(a)時域波形

(b)FFT

(c)CWT圖11 IC1的時域波形、FFT和CWT圖Fig.11 The time domain waveform, FFT and CWT of IC1

(a)時域波形

(b)FFT

(c)CWT圖12 CMF23的時域波形、FFT和CWT圖Fig.12 The time domain waveform, FFT and CWT of CMF23

圖13 2 100 r/min、25%負荷下的缸壓和缸蓋振動的相干函數(shù)Fig.13 The coherence function of cylinder pressure and cylinder head vibration of 2 100 r/min and 25% load

圖14 2 100 r/min倒拖工況下的活塞敲擊振動和分量CMF23的頻譜Fig.14 The spectrum of CMF23 and piston slap vibration of 2 100 r/min and 25% load

4 結(jié)論

(1)通過試驗測量內(nèi)燃機的單通道缸蓋振動信號后,用VMD算法、RobustICA算法及組合模態(tài)函數(shù)法,并結(jié)合頻譜分析、連續(xù)小波時頻分析、相干函數(shù)法及倒拖試驗準確有效地分離識別出了內(nèi)燃機的燃燒信號和活塞敲擊信號。

(2)在內(nèi)燃機1 000 r/min空載工況下,燃燒信號和活塞敲擊信號的頻率成分分別集中在4 350 Hz和1 150 Hz;在2 100 r/min、25%負荷工況下,燃燒信號和活塞敲擊信號的頻率成分分別集中在4 250 Hz和1 500 Hz。由此可看出燃燒信號的頻率成分都集中在4 300 Hz附近,這主要是由于內(nèi)燃機燃燒階段的高頻振蕩所引起;但是活塞敲擊信號的頻率在不同工況下有一定的差別,主要是因為內(nèi)燃機在轉(zhuǎn)速較高負荷較大的工況下,活塞撞擊缸壁的頻率較高,撞擊力較大,從而導(dǎo)致活塞敲擊信號的頻率較高。

參考文獻:

[1] CHEN J, RANDALL R B, PEETERS B. Advanced Diagnostic System for Piston Slap Faults in IC Engines, Based on the Non-stationary Characteristics of the Vibration Signals[J]. Mechanical Systems & Signal Processing, 2016, 75:434-454.

[2] CHEN J, RANDALL R B. Intelligent Diagnosis of Bearing Knock Faults in Internal Combustion Engines Using Vibration Simulation[J]. Mechanism & Machine Theory, 2016, 104:161-176.

[3] 吳虎勝,呂建新,吳廬山,等.基于EMD和SVM的柴油機氣閥機構(gòu)故障診斷[J].中國機械工程, 2010,21(22):2710-2714.

WU Husheng, LYU Jianxin, WU Lushan, et al. Fault Diagnosis for Diesel Valve Train Based on SVM and EMD[J]. China Mechanical Engineering, 2010,21(22):2710-2714.

[4] 丁巖,邵毅敏.基于角域振動信號特征統(tǒng)計量的發(fā)動機故障分類方法[J]. 中國機械工程, 2014, 25(10):1374-1380.

DING Yan, SHAO Yimin. Engine Fault Classification Based on Characteristic Statistics of Vibration Signals in Angle Domain[J]. China Mechanical Engineering, 2014, 25(10):1374-1380.

[5] 劉建敏,李曉磊,喬新勇,等.基于EMD和STFT柴油機缸蓋振動信號時頻分析[J]. 噪聲與振動控制, 2013(2):133-137.

LIU Jianmin, LI Xiaolei, QIAO Xinyong, et al. Time-Frequency Analysis of Vibration Signal from Cylinder Head of a Diesel Engine Based on EMD and STFT[J]. Noise and Vibration Control, 2013(2):133-137.

[6] LIU X, RANDALL R B. Blind Source Separation of Internal Combustion Engine Piston Slap from Other Measured Vibration Signals[J]. Mechanical Systems & Signal Processing, 2005, 19(6):1196-1208.

[7] LIU X, RANDALL R B, ANTONI J. Blind Separation of Internal Combustion Engine Vibration Signals by a Deflation Method[J]. Mechanical Systems & Signal Processing, 2008, 22(5):1082-1091.

[8] 康斌,向陽.獨立分量分析在柴油機缸蓋振動信號分離中的應(yīng)用[J]. 武漢理工大學(xué)學(xué)報(交通科學(xué)與工程版), 2012, 36(4):753-756.

KANG Bin, XIANG Yang. Application of Independent Component Analysis to Decomposition of Engine Vibration Signal[J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2012, 36(4):753-756.

[9] 韓春楊, 姚國鳳, 趙建,等.柴油機振動信號盲分離組合算法[J]. 振動與沖擊, 2014, 33(6):44-47.

HAN Chunyang, YAO Guofeng, ZHAO Jian, et al. Combination Algorithm for Blind Separation of Diesel Engine Vibration Signal [J]. Journal of Vibration and Shock, 2014, 33(6):44-47.

[10] DU X, LI Z, BI F, et al. Source Separation of Diesel Engine Vibration Based on the Empirical Mode Decomposition and Independent Component Analysis[J]. Chinese Journal of Mechanical Engineering, 2012, 25(3):557-563.

[11] HUANG N E, SHEN Z, LONG S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J]. Proceedings of the Royal Society: A Mathematical Physical & Engineering Sciences, 1998, 454:903-995.

[12] ?VOKELJ M, ZUPAN S, PREBIL I. EEMD-based Multiscale ICA Method for Slewing Bearing Fault Detection and Diagnosis [J]. Journal of Sound and Vibration, 2016, 370: 394-423.

[13] ZHANG M, JIANG Z, FENG K. Research on Variational Mode Decomposition in Rolling Bearings Fault Diagnosis of the Multistage Centrifugal Pump[J]. Mechanical Systems & Signal Processing, 2017, 93:460-493.

[14] DRAGOMIRETSKIY K, ZOSSO D. Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing, 2014, 62(3):531-544.

[15] HYVARINEN A. Fast and Robust Fixed-point Algorithms for Independent Component Analysis[J]. IEEE Transactions on Neural Networks, 1999, 10(3):626-634.

[16] ZARZOSO V, COMON P. Robust Independent Component Analysis by Iterative Maximization of the Kurtosis Contrast with Algebraic Optimal Step Size [J]. IEEE Transactions on Neural Networks, 2010, 21(2):248-261.

[18] ANTONI J, DUCLEAUX N, NGHIEM G, et al. Separation of Combustion Noise in IC Engines under Cyclo-non-stationary Regime[J]. Mechanical Systems & Signal Processing, 2013, 38(1):223-236.

猜你喜歡
模態(tài)振動信號
振動的思考
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
基于FPGA的多功能信號發(fā)生器的設(shè)計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
基于LabVIEW的力加載信號采集與PID控制
國內(nèi)多模態(tài)教學(xué)研究回顧與展望
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 亚洲国产精品日韩av专区| 国产成人福利在线视老湿机| 亚洲an第二区国产精品| 国产欧美日韩另类精彩视频| 免费人欧美成又黄又爽的视频| 九九热免费在线视频| 国产福利不卡视频| 人与鲁专区| 丁香婷婷久久| 日韩精品亚洲精品第一页| 亚洲综合婷婷激情| 亚洲日韩国产精品无码专区| 国产亚洲男人的天堂在线观看 | 高清久久精品亚洲日韩Av| 国产日韩欧美在线播放| 婷婷色婷婷| 国产精品9| 青青草一区| 成人亚洲国产| 日韩毛片在线播放| 日本爱爱精品一区二区| 亚洲人成网站18禁动漫无码| 久久一级电影| 色成人亚洲| 国产成人av一区二区三区| 亚洲人成人无码www| 欧美日韩免费在线视频| 亚洲视频四区| 人人91人人澡人人妻人人爽| 中字无码av在线电影| 日韩无码真实干出血视频| 女人一级毛片| 国产三级精品三级在线观看| 久久国产精品电影| 一区二区三区四区日韩| 久久无码av三级| 久久亚洲日本不卡一区二区| 99久久精品免费看国产电影| 亚洲第一网站男人都懂| 乱人伦视频中文字幕在线| 国产高清在线观看| 夜夜爽免费视频| www.91在线播放| 国产精品漂亮美女在线观看| 亚洲一区二区三区麻豆| 亚洲欧美成人在线视频| 亚洲第一黄色网| 国产福利免费观看| Aⅴ无码专区在线观看| 2048国产精品原创综合在线| 99999久久久久久亚洲| 久久精品国产精品青草app| 亚洲男人天堂久久| 国产三级国产精品国产普男人| 久久久久久久久亚洲精品| 亚洲综合久久成人AV| 国产亚卅精品无码| 亚洲国产精品VA在线看黑人| 2018日日摸夜夜添狠狠躁| 2020国产免费久久精品99| 国产成人一区免费观看| 亚洲日韩第九十九页| 四虎影视永久在线精品| 国产免费a级片| 欧美日韩午夜| 欧美成人A视频| 三区在线视频| 五月综合色婷婷| 91午夜福利在线观看| 亚洲高清无码久久久| 国产精品3p视频| 欧美五月婷婷| 国产精品19p| 亚洲最猛黑人xxxx黑人猛交 | av在线5g无码天天| 亚洲精品在线观看91| 91久久国产综合精品女同我| 国产精品自拍合集| 欧美成人午夜视频免看| www.99精品视频在线播放| 午夜欧美理论2019理论| 欧美视频在线播放观看免费福利资源 |