亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于MCKD和VMD的滾動(dòng)軸承微弱故障特征提取

        2017-11-04 01:27:34夏均忠白云川于明奇汪治安
        振動(dòng)與沖擊 2017年20期
        關(guān)鍵詞:模態(tài)故障信號(hào)

        夏均忠, 趙 磊, 白云川, 于明奇, 汪治安

        (軍事交通學(xué)院 軍用車輛工程技術(shù)研究中心,天津 300161)

        基于MCKD和VMD的滾動(dòng)軸承微弱故障特征提取

        夏均忠, 趙 磊, 白云川, 于明奇, 汪治安

        (軍事交通學(xué)院 軍用車輛工程技術(shù)研究中心,天津 300161)

        針對(duì)滾動(dòng)軸承早期故障特征非常微弱,易受隨機(jī)噪聲和其他信號(hào)干擾而難以提取等現(xiàn)象,提出了用最大相關(guān)峭度解卷積(Maximum Correlated Kurtosis Deconvolution,MCKD)和變分模態(tài)分解(Variational Mode Decomposition,VMD)相結(jié)合的方法提取滾動(dòng)軸承故障特征。首先用MCKD進(jìn)行信號(hào)增強(qiáng),然后利用VMD得到一系列模態(tài),應(yīng)用互相關(guān)系數(shù)和峭度準(zhǔn)則篩選包含故障信息較為豐富的模態(tài)進(jìn)行重構(gòu)降噪,最后對(duì)重構(gòu)信號(hào)進(jìn)行包絡(luò)解調(diào)提取故障特征。通過(guò)仿真分析和軸承故障模擬實(shí)驗(yàn)驗(yàn)證了該方法的有效性,可以精確地分離軸承故障振動(dòng)信號(hào)的不同頻率成分。

        滾動(dòng)軸承;最大相關(guān)峭度解卷積;變分模態(tài)分解;互相關(guān)系數(shù);峭度準(zhǔn)則

        滾動(dòng)軸承頻發(fā)的各種故障大部分以局部缺陷的形式發(fā)生在軸承工作周期的早期,而且多數(shù)是潛在故障,具有特征信息微弱、易被噪聲淹沒(méi)及信噪比低等特點(diǎn),使得產(chǎn)生的周期性脈沖往往淹沒(méi)在背景噪聲當(dāng)中不易識(shí)別和提取。為了準(zhǔn)確提取軸承早期微弱故障特征,排除噪聲對(duì)信號(hào)的干擾,需要增強(qiáng)周期沖擊信號(hào)、降低噪聲來(lái)提高信噪比。

        最大相關(guān)峭度解卷積(Maximum Correlated Kurtosis Deconvolution,MCKD)是McDonald等[1]在最小熵解卷積的基礎(chǔ)上提出的一種新的增強(qiáng)信號(hào)周期性沖擊成分的解卷積技術(shù),并應(yīng)用于齒輪的故障診斷。該方法以相關(guān)峭度為評(píng)價(jià)指標(biāo),通過(guò)迭代過(guò)程實(shí)現(xiàn)解卷積,進(jìn)而突出信號(hào)中被噪聲淹沒(méi)的連續(xù)脈沖序列。文獻(xiàn)[2]應(yīng)用MCKD預(yù)處理信號(hào),再使用1.5維譜分析判斷軸承故障特征頻率成分,有效地提取了早期故障特征信息。文獻(xiàn)[3]針對(duì)機(jī)械信號(hào)中存在的噪聲會(huì)降低重分配小波尺度譜的時(shí)頻分布可讀性,提出基于MCKD和重分配小波尺度譜的方法對(duì)旋轉(zhuǎn)機(jī)械復(fù)合故障進(jìn)行診斷。

        經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是Huang在上個(gè)世紀(jì)末期提出的一種分析非線性、非平穩(wěn)信號(hào)的遞歸式模態(tài)分解方法[4]。由于缺乏良好的數(shù)學(xué)理論以及在分解時(shí)容易受到噪聲和采樣率的限制,EMD存在一些問(wèn)題,如模態(tài)混疊、端點(diǎn)效應(yīng)以及過(guò)包絡(luò)、欠包絡(luò)等。鑒于此,Dragomiretskiy K和Zosso D于2014年提出了一種新的非遞歸式自適應(yīng)模態(tài)分解算法——變分模態(tài)分解(Variational Mode Decomposition,VMD)[5]。文獻(xiàn)[6]應(yīng)用VMD和消除趨勢(shì)波動(dòng)分析對(duì)信號(hào)進(jìn)行降噪處理。

        為了增強(qiáng)軸承故障振動(dòng)信號(hào)中的周期性沖擊成分,降低噪聲的干擾,更加精確地提取軸承故障特征,提出了基于MCKD和VMD的滾動(dòng)軸承微弱故障特征提取方法。首先用MCKD對(duì)軸承振動(dòng)信號(hào)進(jìn)行增強(qiáng),然后利用VMD得到一系列模態(tài),并結(jié)合互相關(guān)系數(shù)和峭度準(zhǔn)則篩選包含故障信息較為豐富的模態(tài)進(jìn)行重構(gòu)降噪,最后對(duì)重構(gòu)信號(hào)進(jìn)行包絡(luò)解調(diào),精確有效地提取出滾動(dòng)軸承故障特征頻率。

        1 最大相關(guān)峭度解卷積(MCKD)

        MCKD的本質(zhì)就是尋找一個(gè)FIR濾波器f(l)(l為濾波器的長(zhǎng)度),使原始沖擊序列的相關(guān)峭度最大,以此恢復(fù)其所具有的特性,達(dá)到增強(qiáng)信號(hào)的目的。

        相關(guān)峭度的定義為

        (1)

        為選取一個(gè)最優(yōu)濾波器f(l),使CKM(T)最大,令

        (2)

        式中:T為沖擊信號(hào)的周期;M為位移數(shù);f為濾波器向量,f=[f1,f2,…,fL]T;L為FIR濾波器的長(zhǎng)度。

        上述優(yōu)化求解問(wèn)題等價(jià)于求解方程

        (3)

        用矩陣的形式來(lái)表示最終的濾波器系數(shù)

        (4)

        其中:

        r=[0T2T…mT]

        通過(guò)迭代方式求濾波器f參數(shù)的過(guò)程如下:

        (1) 確定周期T、移位數(shù)M和濾波器的長(zhǎng)度L。

        (3) 計(jì)算濾波后的輸出信號(hào)y。

        (4) 根據(jù)y計(jì)算αm和β。

        (5) 計(jì)算新濾波器的系數(shù)f。

        (6) 如果濾波前后信號(hào)的ΔCKM(T)>ε時(shí),跳轉(zhuǎn)到步驟(3)繼續(xù)循環(huán),否則停止遞歸。ε是用來(lái)控制迭代終止的較小正數(shù)。

        2 變分模態(tài)分解(VMD)

        2.1 基本原理

        VMD是一種完全非遞歸的信號(hào)分解方法。它可以將任意信號(hào)f(t)分解成許多圍繞在中心頻率ωk周圍的模態(tài)分量信號(hào)。具體步驟如下:

        (1) 首先用Hilbert變換計(jì)算每個(gè)模態(tài)uk的相關(guān)解析信號(hào)以獲得一個(gè)單邊頻譜,然后通過(guò)加入一個(gè)指數(shù)項(xiàng)來(lái)調(diào)整各自估計(jì)的中心頻率,把每個(gè)模態(tài)的頻譜轉(zhuǎn)移到基帶上,最后通過(guò)解調(diào)信號(hào)的高斯平滑來(lái)估計(jì)帶寬,即梯度的二范數(shù)平方。由此產(chǎn)生了一個(gè)由變分問(wèn)題組成的目標(biāo)函數(shù)

        (5)

        (2) 通過(guò)引入拉格朗日乘子λ和二次懲罰項(xiàng)將上述約束性變分問(wèn)題轉(zhuǎn)化為非約束性變分問(wèn)題。增廣拉格朗日表達(dá)式如下

        (6)

        式中:α為懲罰因子;λ(t)是加強(qiáng)約束。二次懲罰項(xiàng)的作用是提高收斂性。

        (3) 通過(guò)尋找迭代子優(yōu)化序列中增廣拉格朗日的鞍點(diǎn),即應(yīng)用交替方向乘子法(優(yōu)化算法)求解式(5)的最小化問(wèn)題。求解式(6)的迭代方法為

        (7)

        (8)

        (9)

        (10)

        式(10)為收斂條件,ε為精度(ε>0),n為迭代次數(shù)。

        利用L2范數(shù)下Parseval/Plancherel傅里葉等距變換在頻域?qū)κ?7)~式(9)進(jìn)行求解

        (11)

        (12)

        (13)

        因此,完整的VMD算法實(shí)現(xiàn)過(guò)程如下:

        一團(tuán)玫瑰的甜香隨即旋轉(zhuǎn)到夏天的身邊,他知道是葉曉曉剛剛洗完澡,稍稍定了定神,說(shuō):“曉曉,我渴了,能給我倒杯水嗎?”

        (3) 對(duì)于所有的ω≥0:

        (4) 判斷是否滿足式(10)的收斂條件,若滿足則停止迭代,否則返回到步驟(2)。

        2.2VMD和EMD的對(duì)比分析

        設(shè)包含噪聲的仿真信號(hào)為

        (14)

        式中:f1=2,f2=24,f3=288,n(t)為隨機(jī)噪聲,采樣頻率為1 000 Hz,采樣點(diǎn)數(shù)為1 000個(gè)。仿真信號(hào)的波形如圖1所示。分別對(duì)其進(jìn)行EMD和VMD分解,得到的分解結(jié)果如圖2、3所示。EMD分解結(jié)果中,IMF1、IMF2包含著高頻諧波和噪聲,IMF3和IMF4出現(xiàn)了模態(tài)混疊現(xiàn)象,IMF4最接近于中頻諧波,IMF6包含有低頻諧波但出現(xiàn)了嚴(yán)重扭曲。VMD分解結(jié)果中,模態(tài)U1、U2包含信號(hào)的低頻、中頻成分,模態(tài)U3中的高頻成分由于噪聲的影響較難分辨,但是通過(guò)VMD算法可以計(jì)算出其中心頻率。所以,VMD比EMD分解誤差小,而且克服了模態(tài)混疊現(xiàn)象,物理意義更明確。

        圖1 仿真信號(hào)的波形圖Fig.1 The waveform of simulated signal

        圖2 EMD分解結(jié)果Fig.2 The result of EMD

        圖3 VMD分解結(jié)果Fig.3 The result of VMD

        3 仿真分析

        為驗(yàn)證本文方法的有效性,利用故障模型模擬軸承故障產(chǎn)生的周期性沖擊信號(hào),并向其中添加白噪聲來(lái)模擬軸承早期故障信號(hào)[7]。仿真信號(hào)如下:

        (15)

        式中:載波頻率fn=3 000 Hz,位移常數(shù)x0=5,阻尼系數(shù)ξ=0.1,沖擊故障的周期T=0.01 s,采樣頻率fs=20 kHz,采樣點(diǎn)數(shù)N=4 096,t為采樣時(shí)刻,n(t)為白噪聲。

        仿真信號(hào)的時(shí)域波形和包絡(luò)譜如圖4所示,軸承故障沖擊信號(hào)已經(jīng)完全被噪聲淹沒(méi),無(wú)法識(shí)別故障信息。對(duì)仿真信號(hào)進(jìn)行MCKD信號(hào)增強(qiáng),得到的時(shí)域波形和包絡(luò)譜如圖5所示,可以看到時(shí)域波形里的沖擊信號(hào)更加明顯,而頻譜圖中的沖擊頻率受到噪聲影響還不能清晰精確地從中識(shí)別出來(lái)。

        圖4 仿真信號(hào)的時(shí)域波形和包絡(luò)譜Fig.4 Time-domain and envelope spectrum of the simulated signal

        對(duì)MCKD處理后的信號(hào)進(jìn)行VMD,根據(jù)互相關(guān)系數(shù)和峭度準(zhǔn)則篩選包含故障信息較為豐富的模態(tài)分量進(jìn)行重構(gòu)降噪,然后對(duì)重構(gòu)信號(hào)進(jìn)行包絡(luò)解調(diào),得到的包絡(luò)譜圖如圖6所示。從圖中可以明顯看出軸承故障周期沖擊頻率(100 Hz)及其倍頻(200 Hz、300 Hz、400 Hz、500 Hz等)。仿真分析證明了該方法的有效性,適用于強(qiáng)噪聲背景下軸承早期故障特征提取。

        圖5 MCKD處理后的時(shí)域波形和包絡(luò)譜Fig.5 Time-domain and envelope spectrum based on MCKD

        圖6 VMD重構(gòu)后的包絡(luò)譜Fig.6 Envelope spectrum based on VMD reconstruction

        4 滾動(dòng)軸承故障特征提取

        4.1 算法流程

        如圖7所示為基于MCKD和VMD的滾動(dòng)軸承故障特征提取流程圖,具體步驟如下:

        (1) 對(duì)采集的振動(dòng)信號(hào)用MCKD增強(qiáng)周期沖擊成分;

        (2) 對(duì)處理后的信號(hào)進(jìn)行VMD分解得到各模態(tài)分量;

        (3) 取互相關(guān)系數(shù)和峭度值較大的模態(tài)分量進(jìn)行重構(gòu)降噪,并求出重構(gòu)信號(hào)的包絡(luò)譜;

        (4) 將滾動(dòng)軸承故障特征頻率與包絡(luò)譜峰值較大處的頻率進(jìn)行比較,從而精確有效地提取出滾動(dòng)軸承故障特征。

        4.2 應(yīng)用實(shí)例

        實(shí)驗(yàn)裝置由驅(qū)動(dòng)電機(jī)、振動(dòng)加速度傳感器、扭矩解碼/編碼器、聯(lián)軸器和功率計(jì)等組成,如圖8所示[8]。試驗(yàn)軸承為SKF 6205-2RS深溝球軸承,其技術(shù)參數(shù)見(jiàn)表1。使用電火花在軸承內(nèi)、外圈加工直徑均為0.18 mm(深度均為0.28 mm),模擬內(nèi)、外圈故障。電機(jī)轉(zhuǎn)速為1 750 r/min,采樣頻率為12 kHz。

        圖7 滾動(dòng)軸承故障特征提取流程圖Fig.7 Flow chart of rolling element bearing fault feature extraction

        圖8 實(shí)驗(yàn)裝置示意圖Fig.8 The schematic diagram of experimental device表1 試驗(yàn)軸承技術(shù)參數(shù)Tab.1 The specification of tested bearing

        滾動(dòng)體直徑d/mm節(jié)圓直徑D/mm內(nèi)徑dm/mm外徑Do/mm滾動(dòng)體數(shù)Z接觸角α/(°)839255290

        如圖9所示,為滾動(dòng)軸承內(nèi)、外圈故障振動(dòng)信號(hào)的時(shí)域波形及頻譜。從圖中可以發(fā)現(xiàn)內(nèi)、外圈故障信號(hào)時(shí)域波形由于噪聲的影響已無(wú)法識(shí)別故障特征,且頻譜圖中出現(xiàn)了一些高頻成分。

        (a) 內(nèi)圈故障

        (b) 外圈故障圖9 軸承振動(dòng)信號(hào)的時(shí)域波形及頻譜Fig.9 Time-domain waveform and spectrum of bearing vibration signal

        對(duì)軸承振動(dòng)信號(hào)用MCKD進(jìn)行增強(qiáng),信號(hào)的時(shí)域波形和包絡(luò)譜如圖10所示。內(nèi)圈故障振動(dòng)信號(hào)的峭度值由5.50增加到10.82,外圈故障振動(dòng)信號(hào)的峭度值由7.84增加到23.05,增加幅度較為明顯;時(shí)域波形中出現(xiàn)了明顯的沖擊成分,且頻譜中的低頻沖擊成分較為增強(qiáng),但是還存在噪聲的影響,仍無(wú)法清晰識(shí)別軸承故障特征頻率。

        (a) 內(nèi)圈故障

        (b) 外圈故障圖10 MCKD處理后的時(shí)域波形和包絡(luò)譜Fig.10 Time-domain and envelope spectrum based on MCKD

        再對(duì)MCKD增強(qiáng)后的信號(hào)進(jìn)行VMD分解,結(jié)果如圖11所示。分別計(jì)算各模態(tài)與原信號(hào)的相關(guān)系數(shù)及其峭度值,見(jiàn)表2。從表中可以看出模態(tài)分量U2、U3、U4與原信號(hào)的相關(guān)系數(shù)、峭度值較大,于是選擇上述3個(gè)模態(tài)分量進(jìn)行信號(hào)重構(gòu),再對(duì)重構(gòu)后的信號(hào)進(jìn)行包絡(luò)解調(diào),得到其包絡(luò)譜如圖12所示。

        (a) 內(nèi)圈故障

        (b) 外圈故障圖11 VMD分解結(jié)果Fig.11 The result of VMD表2 內(nèi)、外圈故障的相關(guān)系數(shù)和峭度值Tab.2 Correlation coefficient and kurtosis of innerand outer race fault

        特征參數(shù)模態(tài)UU1U2U3U4內(nèi)圈峭度值相關(guān)系數(shù)3.380.446.490.456.820.534.990.55外圈峭度值相關(guān)系數(shù)7.850.467.340.558.590.5711.540.53

        從圖12(a)可以看出,軸承內(nèi)圈故障時(shí)振動(dòng)信號(hào)的包絡(luò)譜包括回轉(zhuǎn)頻率,內(nèi)圈故障特征頻率及其二次諧波(2BPFI=316.26 Hz)、三次諧波(3BPFI=474.39 Hz),調(diào)制邊帶(BPFI-fr=128.97 Hz、BPFI+fr=187.29 Hz)。從圖12(b)可以看出,軸承外圈故障時(shí)的振動(dòng)信號(hào)包括回轉(zhuǎn)頻率,外圈故障特征頻率及其二次諧波(2BPFO=208.60 Hz)、三次諧波(3BPFO=312.90 Hz)、四次諧波(4BPFO=417.20 Hz)、五次諧波(5BPFO=521.50 Hz)。綜合分析可以精確判斷出滾動(dòng)軸承內(nèi)、外圈故障,分析結(jié)果與實(shí)際相符。因此,基于MCKD和VMD的方法可以精確地提取滾動(dòng)軸承微弱故障特征頻率。

        (a) 內(nèi)圈故障

        (b) 外圈故障圖12 VMD重構(gòu)后的包絡(luò)譜Fig.12 Envelope spectrum based on VMD reconstruction

        5 結(jié) 論

        (1) 滾動(dòng)軸承故障早期階段沖擊特征非常微弱,而且容易受到其他信號(hào)的干擾,MCKD以相關(guān)峭度最大化為優(yōu)化目標(biāo),通過(guò)迭代過(guò)程實(shí)現(xiàn)信號(hào)的解卷積,增強(qiáng)信號(hào)中被強(qiáng)噪聲干擾的脈沖沖擊,提高了峭度值和信噪比,非常適用于增強(qiáng)軸承早期故障信號(hào)的沖擊特性。

        (2) 故障軸承產(chǎn)生的振動(dòng)信號(hào)中往往含有較多的調(diào)幅調(diào)頻分量,解調(diào)后的信號(hào)中包含了軸承回轉(zhuǎn)頻率及其倍頻,故障特征頻率及其倍頻等調(diào)制成分。VMD是一種新的非遞歸式自適應(yīng)模態(tài)分解方法,它能自適應(yīng)地確定相關(guān)頻帶,同時(shí)估計(jì)對(duì)應(yīng)的模態(tài)??朔藗鹘y(tǒng)的遞歸式模態(tài)分解(如EMD、LMD)因沖擊成分和噪

        聲引起的模態(tài)混疊,可以更加精確地提取不同故障類型的特征信息。

        [1] MCDONALD G L,ZHAO Q,ZUO M J.Maximum correlated kurtosis deconvolution and application on gear tooth chip fault detection[J].Mechanical Systems and Signal Processing,2012,33: 237-255.

        [2] 唐貴基,王曉龍.最大相關(guān)峭度解卷積結(jié)合1.5維譜的滾動(dòng)軸承早期故障特征提取方法[J].振動(dòng)與沖擊,2015,34(12): 79-84.

        TANG Guiji,WANG Xiaolong.Feature extraction for rolling bearing incipient fault based on maximum correlated kurtosis deconvolution and 1.5 dimension spectrum[J]. Journal of Vibration and Shock,2015,34(12): 79-84.

        [3] 鐘先友,趙春華,陳保家,等.基于MCKD和重分配小波尺度譜的旋轉(zhuǎn)機(jī)械復(fù)合故障診斷研究[J].振動(dòng)與沖擊,2015,34(7): 156-161.

        ZHONG Xianyou,ZHAO Chunhua,CHEN Baojia,et al.Rotating machinery fault diagnosis based on maximum correlation kurtosis deconvolution and reassigned wavelet scalogram[J].Journal of Vibration and Shock,2015,34(7): 156-161.

        [4] 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 of London,1998,454(1971): 903-995.

        [5] DRAGOMIRETSKIY K,ZOSSO D.Variational mode decomposition[J].IEEE Transaction on Signal Processing,2014,62(3): 531-544.

        [6] LIU Y Y,YANG G L,LI M,et al.Variational mode decomposition denoising combined the detrended fluctuation analysis[J]. Signal Processing,2016,125: 349-364.

        [7] WANG X D,ZI Y Y,HE Z J.Multiwavelet construction via an adaptive symmetric lifting scheme and its applications for rotating machinery fault diagnosis[J]. Measurement Science and Technology,2009,20(4): 1-17.

        [8] SMITH W A,RANDALL R B.Rolling element bearing diagnostics using the Case Western Reserve University data:a benchmark study[J].Mechanical Systems and Signal Processing,2015,64: 100-131.

        FeatureextractionforrollingelementbearingweakfaultbasedonMCKDandVMD

        XIA Junzhong, ZHAO Lei, BAI Yunchuan, YU Mingqi, WANG Zhi’an

        (Research Center of Military Vehicle Engineering & Technology,Academy of Military Transportation,Tianjin 300161,China)

        The fault feature of rolling element bearings in early failure period is so weak and susceptible to random noise and other signal interferences that it is very difficult to be extracted. To solve this problem, the maximum correlated kurtosis deconvolution was combined with the variational mode decomposition for extracting rolling element bearing fault feature. Firstly the signal was enhanced by MCKD and decomposed into several modes by VMD. The signal then was reconstructed and noise was reduced with the mode,which was selected by the comparative correlation coefficient and the kurtosis criterion. Fnally the envelope demodulation was used to extract fault feature. The simulating signal analysis and bearing fault simulator show the validity of the method. This method can accurately separate different frequency components of bearing fault vibration signals.

        rolling element bearing;maximum correlated kurtosis deconvolution(MCKD);variational mode decomposition(VMD);correlation coefficient;kurtosis criterion

        2016-04-07 修改稿收到日期: 2016-08-21

        夏均忠 男,博士,教授,1967年生

        趙磊 男,研究生,1991年生

        E-mail:1016091436@qq.com

        TH133.33

        A

        10.13465/j.cnki.jvs.2017.20.013

        猜你喜歡
        模態(tài)故障信號(hào)
        信號(hào)
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        故障一點(diǎn)通
        基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
        電子制作(2018年11期)2018-08-04 03:25:42
        奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
        基于LabVIEW的力加載信號(hào)采集與PID控制
        國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
        故障一點(diǎn)通
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
        江淮車故障3例
        亚洲av日韩精品久久久久久a| 国产一级片内射在线视频| 丰满人妻一区二区三区52| 一区二区三区天堂在线| 性无码免费一区二区三区在线| 青青草国产成人99久久| 亚洲日韩AV无码美腿丝袜 | 久久精品国产亚洲AV无码不| 久久无人码人妻一区二区三区| 极品尤物人妻堕落沉沦| 中文乱码字慕人妻熟女人妻| 一本大道香蕉视频在线观看| 国产一区二区三区色区| 亚洲国产国语在线对白观看| 天堂资源中文最新版在线一区| 婷婷综合五月| 中文在线最新版天堂av| 青青草成人在线免费视频| 99精品欧美一区二区三区| 国产精品久久久久尤物| 日本一道高清在线一区二区| 无码熟妇人妻av在线网站| 国产精品毛片久久久久久久| 久久免费视亚洲无码视频| 免费在线观看草逼视频| 国产免费久久精品99久久| 免费a级毛片出奶水| 在线观看精品国产福利片87| 精品少妇人妻av一区二区蜜桃| 久久精品亚洲一区二区三区浴池| 国产精品成人va| 中文字幕精品乱码一区| 欧美性白人极品1819hd| 亚洲有码转帖| 久久精品韩国日本国产| 白白色发布的在线视频| 成年无码av片在线| 亚洲国产精品线观看不卡| av免费在线播放观看| 成人影院yy111111在线| 中文字幕日韩高清|