• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    用于滾動(dòng)軸承轉(zhuǎn)頻估計(jì)的EMD軟閾值降噪算法

    2016-08-01 10:00:34程衛(wèi)東趙德尊
    關(guān)鍵詞:振動(dòng)故障信號(hào)

    程衛(wèi)東,趙德尊

    (北京交通大學(xué) 機(jī)械與電子控制工程學(xué)院,北京 100044)

    ?

    用于滾動(dòng)軸承轉(zhuǎn)頻估計(jì)的EMD軟閾值降噪算法

    程衛(wèi)東,趙德尊

    (北京交通大學(xué) 機(jī)械與電子控制工程學(xué)院,北京 100044)

    摘要:針對(duì)經(jīng)驗(yàn)?zāi)J椒纸?EMD)互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法與小波閾值降噪方法的不足,提出EMD與小波軟閾值降噪相結(jié)合的降噪方法.該方法主要包括以下4部分:1)對(duì)原始信號(hào)進(jìn)行EMD分解得到固有模態(tài)函數(shù)(IMF)的集合;2)計(jì)算各個(gè)IMF與原始信號(hào)的互相關(guān)系數(shù)以及各IMF的峭度值;3)利用互相關(guān)系數(shù)-峭度準(zhǔn)則選擇需要降噪的IMF以及需要剔除的IMF;4)對(duì)選定的IMF進(jìn)行閾值降噪后與剩余的IMF相加重構(gòu)信號(hào).利用仿真和實(shí)測(cè)的故障軸承信號(hào)對(duì)所提出算法以及EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法進(jìn)行對(duì)比驗(yàn)證.結(jié)果表明:采用EMD軟閾值降噪方法比采用EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法對(duì)信號(hào)進(jìn)行預(yù)處理,更能確保軸承振動(dòng)信號(hào)的完整性,突出信號(hào)的故障特征,降低瞬時(shí)轉(zhuǎn)頻估計(jì)的誤差.

    關(guān)鍵詞:經(jīng)驗(yàn)?zāi)J椒纸?EMD);軟閾值降噪;滾動(dòng)軸承;信號(hào)預(yù)處理;瞬時(shí)轉(zhuǎn)頻(IRF)估計(jì)

    為了克服相關(guān)硬件設(shè)備安裝空間及安裝成本對(duì)傳統(tǒng)階比分析技術(shù)的限制[1-3],階比分析技術(shù)框架下不依靠轉(zhuǎn)速計(jì)直接根據(jù)時(shí)域振動(dòng)信號(hào)對(duì)瞬時(shí)轉(zhuǎn)頻進(jìn)行估計(jì)的研究引起了廣泛關(guān)注[4-5].Boashash等[6-7]詳細(xì)討論了瞬時(shí)頻率這一概念并從時(shí)頻分布的角度對(duì)其進(jìn)行估計(jì).郭瑜等[8]利用峰值搜索算法對(duì)電動(dòng)機(jī)變轉(zhuǎn)速振動(dòng)信號(hào)時(shí)頻圖中的瞬時(shí)頻率進(jìn)行提取以估計(jì)其瞬時(shí)轉(zhuǎn)頻.趙曉平等[9]則提出了短時(shí)傅里葉與Viterbi擬合算法相結(jié)合的方法對(duì)與臥旋離心機(jī)瞬時(shí)轉(zhuǎn)頻相關(guān)的瞬時(shí)頻率進(jìn)行提取.然而,對(duì)于變轉(zhuǎn)速滾動(dòng)軸承振動(dòng)信號(hào),基于時(shí)頻譜的瞬時(shí)轉(zhuǎn)頻估計(jì)方法無(wú)能為力.為此,Wang等[10]提出了基于瞬時(shí)故障特征頻率(instantaneous fault characteristic frequency, IFCF)的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻(instantaneous rotational frequency, IRF)估計(jì)方法.該方法對(duì)于滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻的提取有較好的效果.但實(shí)際工況中滾動(dòng)軸承低轉(zhuǎn)速運(yùn)行[11-12]和早期故障[13-15]普遍存在.此時(shí),故障沖擊信息易被背景噪音覆蓋,較低信號(hào)的信噪比.基于IFCF的IRF估計(jì)方法由于受到噪聲干擾,估計(jì)的瞬時(shí)轉(zhuǎn)頻將會(huì)產(chǎn)生較大的誤差,其中根據(jù)實(shí)驗(yàn)數(shù)據(jù)可知低信噪比情況下誤差可達(dá)到25%,甚至更多,這將嚴(yán)重影響后續(xù)滾動(dòng)軸承的故障診斷.因此,需要對(duì)原始信號(hào)進(jìn)行降噪預(yù)處理以降低滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)的誤差.

    隨著經(jīng)驗(yàn)?zāi)B(tài)分解(emprical mode decomposition, EMD)方法的提出,近年來(lái)很多學(xué)者嘗試將該技術(shù)及其擴(kuò)展算法用于振動(dòng)信號(hào)的降噪預(yù)處理.蘇文勝等[16]提出了基于互相關(guān)系數(shù)-峭度準(zhǔn)則的EMD 降噪方法并運(yùn)用于滾動(dòng)軸承的早期故障診斷.彭暢等[17]則提出了EEMD(ensemble EMD)度量因子(相似性度量和基于距離的度量因子)的降噪方法.唐貴基等[18]將互相關(guān)系數(shù)-峭度準(zhǔn)則與EEMD相結(jié)合選擇最優(yōu)固有模態(tài)函數(shù)(intrinsic mode function, IMF)分量重構(gòu)信號(hào)以降低滾動(dòng)軸承振動(dòng)信號(hào)中的噪聲成分.該類算法的主要思想是根據(jù)一定準(zhǔn)則篩選出含有有用成分的IMF分量對(duì)信號(hào)進(jìn)行重構(gòu),其他分量則被認(rèn)為是偽分量或者噪聲予以剔除.然而,在低信噪比條件下,比如:早期故障或者低轉(zhuǎn)速運(yùn)行,滾動(dòng)軸承振動(dòng)信號(hào)中含有強(qiáng)烈的背景噪音,往往噪聲起主導(dǎo)作用的IMF分量中也含有故障沖擊成分.因此,單純根據(jù)相關(guān)準(zhǔn)則選擇某幾個(gè)IMF分量重構(gòu)信號(hào)容易丟失原有信號(hào)的有用信息,對(duì)后續(xù)的瞬時(shí)轉(zhuǎn)頻估計(jì)以及階比分析的結(jié)果造成影響.小波變換[19-20]由于基函數(shù)固定、多分辨率恒定等因素的限制缺乏自適應(yīng)性.對(duì)于滾動(dòng)軸承變轉(zhuǎn)速振動(dòng)信號(hào)而言,小波閾值去噪過(guò)程中最優(yōu)分解層次以及小波基函數(shù)不易確定,降噪步驟較為繁瑣.

    為此,本文在上述方法的基礎(chǔ)上提出EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪與小波軟閾值降噪相結(jié)合的方法對(duì)滾動(dòng)軸承變轉(zhuǎn)速信號(hào)進(jìn)行預(yù)處理,然后根據(jù)基于瞬時(shí)故障特征頻率的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法對(duì)轉(zhuǎn)頻進(jìn)行估計(jì).EMD軟閾值降噪方法通過(guò)對(duì)原始信號(hào)進(jìn)行EMD分解得到本征模態(tài)函數(shù);設(shè)定互相關(guān)系數(shù)和譜峭度值閾值,根據(jù)互相關(guān)系數(shù)-峭度準(zhǔn)則剔除偽分量、選出需要進(jìn)行降噪的IMF分量,對(duì)選定的IMF分量進(jìn)行軟閾值降噪;最后將降噪后的IMF分量與剩余的IMF分量相加重構(gòu)信號(hào).

    1EMD降噪算法在滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)中的問(wèn)題

    EMD是由Huang等[21]提出的一種自適應(yīng)的信號(hào)分解方法,基于信號(hào)的局部特征,無(wú)須選擇基函數(shù)就可將原始信號(hào)分解為許多窄帶分量,即IMF.每個(gè)IMF分量滿足以下2個(gè)條件:1)在整個(gè)數(shù)據(jù)段內(nèi),極值點(diǎn)的個(gè)數(shù)和過(guò)零點(diǎn)的個(gè)數(shù)必須相等或最多相差一個(gè);2)在任意時(shí)刻,由局部極大值點(diǎn)形成的上包絡(luò)線和由局部極小值點(diǎn)形成的下包絡(luò)線的平均值為零.IMF分量按頻率從高到低依次分布,每個(gè)IMF分量中的頻率成分隨信號(hào)本身變化而變化.不同于傳統(tǒng)的濾波器,EMD適用于非線性或者非平穩(wěn)信號(hào)[21].另外,EMD方法不受Heisenberg 測(cè)不準(zhǔn)原理[21]的限制,能夠獲得較高的分辨率.因此,EMD方法被廣泛應(yīng)用于機(jī)械振動(dòng)和地震分析等領(lǐng)域.

    基于互相關(guān)系數(shù)-峭度準(zhǔn)則的EMD降噪方法中,互相關(guān)系數(shù)代表IMF分量與原始信號(hào)的相關(guān)度,其值越大,說(shuō)明IMF分量與原始信號(hào)的相關(guān)度越高.因此,可以根據(jù)互相關(guān)系數(shù)辨別EMD分解過(guò)程中由于插值誤差.端點(diǎn)效應(yīng)以及過(guò)分解等諸多因素影響而產(chǎn)生的虛假分量.互相關(guān)系數(shù)的計(jì)算公式如下:

    (1)

    式中:E表示求數(shù)學(xué)期望;μx和μy分別為原始信號(hào)x和y的均值;σx和σy分別為原始信號(hào)x和y的標(biāo)準(zhǔn)差.

    峭度是描述波形尖峰度的一個(gè)無(wú)量綱參數(shù),其對(duì)信號(hào)的沖擊成分非常敏感,即沖擊成分比重越大,峭度值就越大.健康滾動(dòng)軸承的振動(dòng)信號(hào)峭度值分布接近正態(tài)分布,其值大約為3[16].當(dāng)滾動(dòng)軸承出現(xiàn)故障時(shí),峭度值會(huì)明顯增加,因此可以用峭度值衡量IMF分量中故障沖擊成分的多少.峭度值可由下式計(jì)算:

    (2)

    EMD 互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法根據(jù)互相關(guān)系數(shù)-峭度準(zhǔn)則,篩選出相關(guān)系數(shù)較大并且峭度值大于3的IMF分量對(duì)信號(hào)進(jìn)行重構(gòu),其他分量則被認(rèn)為是偽分量或者噪聲予以剔除.然而,噪聲往往表現(xiàn)為高頻信號(hào),滾動(dòng)軸承的高頻共振信號(hào)會(huì)集中分布于前幾個(gè)IMF分量.當(dāng)信號(hào)的信噪比較低時(shí),噪聲起主導(dǎo)作用的IMF分量中也存在故障沖擊成分.因此,單純的根據(jù)互相關(guān)系數(shù)-峭度準(zhǔn)則剔除峭度值小于3而互相關(guān)系數(shù)較大的IMF分量會(huì)造成重構(gòu)信號(hào)丟失原有的故障沖擊成分.

    本文對(duì)故障滾動(dòng)軸承變轉(zhuǎn)速振動(dòng)信號(hào)進(jìn)行仿真(仿真信號(hào)的具體參數(shù)見(jiàn)3.1節(jié)),人為添加高斯白噪聲和信噪比rSNB=-8 dB.采用基于瞬時(shí)故障特征頻率的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法對(duì)轉(zhuǎn)頻進(jìn)行估計(jì)(見(jiàn)圖1),圖1圓圈表示未經(jīng)過(guò)降噪獲取的IFCF趨勢(shì)線,虛線代表通過(guò)瞬時(shí)故障特征頻率趨勢(shì)線和故障特征系數(shù)(仿真信號(hào)故障特征系數(shù)C=2)計(jì)算得到的瞬時(shí)轉(zhuǎn)頻,實(shí)線表示理論瞬時(shí)轉(zhuǎn)頻.其中t代表時(shí)間,f為頻率.雖然估計(jì)的瞬時(shí)轉(zhuǎn)頻在一定程度上反映出理論轉(zhuǎn)頻的變化趨勢(shì),但是由于噪聲干擾,估計(jì)的瞬時(shí)轉(zhuǎn)頻出現(xiàn)誤差,其中最大差值為14.77 Hz,平均差值為9.36 Hz,嚴(yán)重影響了基于瞬時(shí)故障特征頻率的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法的適用性.圖2表示經(jīng)過(guò)EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪后獲取的IFCF趨勢(shì)線、計(jì)算的瞬時(shí)轉(zhuǎn)頻以及理論瞬時(shí)轉(zhuǎn)頻,計(jì)算的瞬時(shí)轉(zhuǎn)頻(虛線)隨時(shí)間的增加反而減小,與理論的瞬時(shí)轉(zhuǎn)頻(實(shí)線)相比,變化趨勢(shì)完全相反.經(jīng)過(guò)上述討論可知:在低信噪比條件下,基于IFCF的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法受噪聲干擾嚴(yán)重,需要對(duì)原始信號(hào)進(jìn)行降噪預(yù)處理,而EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法容易造成滾動(dòng)軸承振動(dòng)信號(hào)中的故障沖擊成分缺失,出現(xiàn)瞬時(shí)轉(zhuǎn)頻估計(jì)不準(zhǔn)的現(xiàn)象.

    圖1 未經(jīng)過(guò)降噪獲取的瞬時(shí)轉(zhuǎn)頻和理論瞬時(shí)轉(zhuǎn)頻Fig.1 Calculated IRF without denoising and theoretical IRF

    圖2 基于EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法得到的瞬時(shí)轉(zhuǎn)頻和理論瞬時(shí)轉(zhuǎn)頻Fig.2 Calculated IRF with EMD cross-correlation coefficient and kurtosis criterion denoising and theoretical IRF

    2EMD軟閾值降噪

    小波閾值降噪是由Donohon等[22-24]提出的一種降噪算法.其原理是原始振動(dòng)信號(hào)經(jīng)過(guò)小波變換之后,有用信號(hào)對(duì)應(yīng)的小波系數(shù)較大,而噪聲信號(hào)一般對(duì)應(yīng)較小的小波系數(shù),根據(jù)這一思想,通過(guò)選取合理的閾值對(duì)小波系數(shù)進(jìn)行量化處理,以達(dá)到降噪的目的.小波閾值降噪主要包括以下3個(gè)基本步驟:1)選擇小波基以及分解層次j,對(duì)信號(hào)進(jìn)行分解;2)對(duì)從第1層到第j層的每一層高頻系數(shù),選擇一個(gè)閾值進(jìn)行閾值量化處理;3)根據(jù)小波分解的第j層低頻系數(shù)和經(jīng)過(guò)量化處理后的第1層到第j層的高頻系數(shù),進(jìn)行小波重構(gòu).

    小波閾值函數(shù)包括軟閾值函數(shù)和硬閾值函數(shù),其中軟閾值處理信號(hào)相對(duì)平滑,軟閾值函數(shù)[23]為

    (3)

    其中,N為某個(gè)尺度j的小波系數(shù)個(gè)數(shù),J是小波分解的最大層次數(shù).

    然而,小波閾值去噪過(guò)程中最優(yōu)分解層數(shù)不易確定,分解以及重構(gòu)小波基函數(shù)選擇不當(dāng)同樣會(huì)造成降噪效果不理想.滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法中運(yùn)用小波閾值降噪,步驟比較繁瑣,不易得到最優(yōu)降噪效果.

    文獻(xiàn)[25]提出將EMD與小波閾值降噪相結(jié)合用于齒輪故障模式識(shí)別與診斷,原理是采用小波閾值降噪處理高頻IMF分量后與低頻的IMF分量相加重構(gòu)信號(hào),最后通過(guò)時(shí)頻分析完成齒輪的故障診斷,并取得了較好的效果.然而,僅根據(jù)文獻(xiàn)[25]無(wú)法確定如何選擇需要降噪的IMF分量以及選擇幾個(gè)高頻IMF分量做降噪處理,小波閾值降噪函數(shù)也沒(méi)有具體指出.另外,由于EMD算法本身的缺陷,容易產(chǎn)生偽分量,將閾值去噪后的IMF分量與其余IMF直接相加也會(huì)忽略低頻噪聲信號(hào)和偽分量對(duì)重構(gòu)信號(hào)的影響.為此,本文將EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪與小波軟閾值降噪相結(jié)合,設(shè)定互相關(guān)系數(shù)閾值θ=0.1,剔除互相關(guān)系數(shù)ρxy<0.1的IMF分量,對(duì)互相關(guān)系數(shù)ρxy≥0.1而峭度值小于3的IMF分量進(jìn)行軟閾值降噪,將降噪后IMF分量與剩余的IMF分量相加得到重構(gòu)信號(hào).對(duì)于預(yù)設(shè)的互相關(guān)系數(shù)閾值,本文經(jīng)過(guò)多次試驗(yàn)對(duì)比確定其值為0.1較為合理,選取互相關(guān)系數(shù)大于0.1的IMF分量進(jìn)行后續(xù)處理,既能去除EMD分解過(guò)程中由于插值誤差和過(guò)分解等諸多因素而產(chǎn)生的虛假分量,也能避免含有軸承信息的IMF分量的缺失.本文選取峭度值3作為分界點(diǎn)(第一章已經(jīng)指出,健康軸承的峭度值成正態(tài)分布,其值大約為3,當(dāng)滾動(dòng)軸承出現(xiàn)故障時(shí),峭度值會(huì)明顯增加).因此峭度值大于3的IMF分量軸承故障沖擊較為明顯,無(wú)須進(jìn)一步降噪處理;而峭度值小于3的IMF分量由于噪聲信號(hào)的干擾軸承的故障特征不明顯,需要進(jìn)行降噪處理以突出故障信息.

    根據(jù)式(3)給出本文所用到的軟閾值函數(shù)公式:

    (4)

    式中:ηj表示第j個(gè)IMF分量的函數(shù)式,γj為第j個(gè)IMF分量的閾值,計(jì)算公式如下:

    EMD軟閾值降噪算法的步驟如下,具體的流程圖如圖3所示.1)對(duì)原始信號(hào)進(jìn)行EMD分解得到IMF分量;2)計(jì)算各IMF分量與原始信號(hào)的互相關(guān)系數(shù)以及各IMF分量的峭度值;3)根據(jù)互相關(guān)系數(shù)-峭度準(zhǔn)則選出需要進(jìn)行閾值降噪的IMF分量以及剔除偽IMF分量;4)對(duì)選定的IMF分量進(jìn)行閾值降噪后與剩余的IMF分量相加得到重構(gòu)信號(hào).

    圖3 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)軟閾值降噪流程圖Fig.3 Flowchart of empirical mode decomposition (EMD) soft-thresholding denoising

    EMD軟閾值降噪方法的優(yōu)點(diǎn)主要包括:1)具備EMD的自適應(yīng)性、多分辨率、完備性等優(yōu)勢(shì);2)確保噪聲起主導(dǎo)作用并且含故障沖擊成分的IMF分量不被剔除;3)避免了小波最優(yōu)分解層數(shù)的選擇,同時(shí)也避免了小波分解過(guò)程中小波基函數(shù)選擇的繁瑣性,從而更容易對(duì)滾動(dòng)軸承非平穩(wěn)振動(dòng)信號(hào)進(jìn)行降噪處理.

    3仿真與實(shí)測(cè)信號(hào)驗(yàn)證

    3.1仿真信號(hào)分析

    為了驗(yàn)證本文降噪方法在低信噪比條件下滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)中的有效性,首先對(duì)故障軸承變轉(zhuǎn)速振動(dòng)信號(hào)進(jìn)行仿真.故障軸承非平穩(wěn)信號(hào)構(gòu)造公式[26]如下:

    sin[wr(t-tm)]}.

    (5)

    式中:Am表示由故障引起的第m個(gè)沖擊的幅值;β為結(jié)構(gòu)衰減系數(shù);wr為由軸承故障引起的共振頻率;u(t)為單位階躍函數(shù);tm為第m個(gè)沖擊發(fā)生的時(shí)間,可由遞推下式確定:

    t1=(1+μ)·1/f(t0)/n,

    tm=(1+μ)·1/f(tm-1)/n,m=2,3,…,M.

    (6)

    式中:f(t)是軸承轉(zhuǎn)頻隨時(shí)間的變化規(guī)律;t0=0;n為每轉(zhuǎn)出現(xiàn)的故障沖擊數(shù).

    本文根據(jù)式(6)對(duì)故障軸承低信噪比非平穩(wěn)信號(hào)進(jìn)行仿真,如下式:

    gb=g(t)+n(t).

    (7)

    式中:Am與時(shí)間成線性變化,軸承瞬時(shí)轉(zhuǎn)頻隨時(shí)間的變化規(guī)律設(shè)定為f(t)=10(t)+30,系統(tǒng)的共振頻率wr=3 500Hz,衰減系數(shù)β=8,采樣頻率Fs=5 000 Hz,采樣點(diǎn)數(shù)為20 000,故障特征系數(shù)C=2;n(t)為高斯白噪聲,信噪比rSNB=-8 dB.仿真信號(hào)的時(shí)域波形圖如圖4所示,U為電壓幅值.

    圖4 仿真信號(hào)的時(shí)域波形圖Fig.4 Simulated vibration signal in time domain

    圖5 前6個(gè)IMF分量(IMF1-6)Fig.5 Intrinsic mode functions (IMF1-6)

    根據(jù)重構(gòu)信號(hào)對(duì)滾動(dòng)軸承的瞬時(shí)轉(zhuǎn)頻進(jìn)行估計(jì),具體步驟如下:1)通過(guò)重構(gòu)信號(hào)的快速譜峭度圖,確定帶通濾波器參數(shù)進(jìn)行帶通濾波;2)對(duì)濾波信號(hào)進(jìn)行Hilbert變換以及STFT變換獲取包絡(luò)時(shí)頻圖,利用基于幅值累加的峰值搜索算法從包絡(luò)時(shí)頻圖中提取瞬時(shí)故障特征頻率并對(duì)其擬合得到瞬時(shí)故障特征頻率趨勢(shì)線;3)瞬時(shí)故障特征頻率趨勢(shì)線與故障特征系數(shù)相除即可得到瞬時(shí)轉(zhuǎn)頻.圖6(a)給出了獲取的瞬時(shí)故障特征頻率趨勢(shì)線、計(jì)算的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻與理論瞬時(shí)轉(zhuǎn)頻的對(duì)比,理論瞬時(shí)轉(zhuǎn)頻f(t)=10(t)+30,估計(jì)的瞬時(shí)轉(zhuǎn)頻與其相吻合,其中兩者的最大差值為2.00 Hz,平均差值為1.03 Hz.分別采用本文降噪方法、EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法以及未采用降噪方法估計(jì)的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻如圖6(b)所示.其中方塊表示經(jīng)EMD軟閾值降噪后估計(jì)的瞬時(shí)轉(zhuǎn)頻,虛線表示經(jīng)過(guò)EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪后獲取的瞬時(shí)轉(zhuǎn)頻,圓圈表示未降噪直接從原始信號(hào)中估計(jì)的瞬時(shí)轉(zhuǎn)頻,實(shí)線為理論瞬時(shí)轉(zhuǎn)頻.通過(guò)綜合對(duì)比可知,單純的采用EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法對(duì)原始信號(hào)進(jìn)行預(yù)處理容易造成信號(hào)中的故障沖擊成分的缺失,出現(xiàn)瞬時(shí)轉(zhuǎn)頻估計(jì)不準(zhǔn)現(xiàn)象;EMD軟閾值降噪優(yōu)于EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪方法,保證信號(hào)完整性的同時(shí),降低了噪聲干擾,突出了信號(hào)的故障信息,進(jìn)而瞬時(shí)轉(zhuǎn)頻的估計(jì)結(jié)果更接近理論值.

    表1IMF1-IMF7與原始信號(hào)的互相關(guān)系數(shù)及其各自峭度值和IMF′的峭度值

    Tab.1Cross-correlation coefficients and kurtosis values of IMF1-IMF7 and IMF′ kurtosis values

    IMF序號(hào)互相關(guān)系數(shù)KIMFIMF'IMF10.72442.54313.8331IMF20.47663.3816—IMF30.30892.96164.6204IMF40.21633.0943—IMF50.15142.87054.4870IMF60.11192.80393.0255IMF70.07283.1745—

    圖6 本文算法獲取的轉(zhuǎn)頻及其與傳統(tǒng)算法的綜合比對(duì)Fig.6 Calculated IRF by proposed algorithm and comparison with traditional algorithms

    3.2實(shí)測(cè)信號(hào)分析

    為進(jìn)一步驗(yàn)證EMD軟閾值降噪在實(shí)際工況低信噪比非平穩(wěn)條件下估計(jì)滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻方法中的實(shí)用性,利用北京交通大學(xué)機(jī)電學(xué)院現(xiàn)代制造技術(shù)綜合實(shí)驗(yàn)中心ST-5000A型多功能柔性轉(zhuǎn)子試驗(yàn)臺(tái)(見(jiàn)圖7)進(jìn)行實(shí)測(cè)驗(yàn)證.該試驗(yàn)臺(tái)由速度調(diào)節(jié)器控制電機(jī)轉(zhuǎn)速,加速度傳感器安裝在離軸承較近的位置以準(zhǔn)確的測(cè)量其振動(dòng)信號(hào),轉(zhuǎn)速計(jì)安裝在最右側(cè)用于測(cè)量轉(zhuǎn)速,采集裝置為YE6231采集卡及其配套的采集軟件.采集卡以相同采樣率同時(shí)采集兩路信號(hào),一路為滾動(dòng)軸承故障信號(hào),另一路為故障信號(hào)相對(duì)應(yīng)的轉(zhuǎn)速脈沖信號(hào).利用電火花切割凹槽模擬軸承外圈早期裂紋故障,實(shí)驗(yàn)軸承型號(hào)為6 000.該軸承的滾珠個(gè)數(shù)為7,滾動(dòng)體直徑為4.8 mm,節(jié)圓直徑為17.65 mm,接觸角為0.根據(jù)其幾何參數(shù)計(jì)算出外圈故障特征系數(shù)Co=2.548.實(shí)測(cè)信號(hào)的時(shí)域波形圖如圖8所示,采樣率為24 000 Hz,采樣時(shí)長(zhǎng)為2 s,瞬時(shí)轉(zhuǎn)頻隨時(shí)間增加而逐步增加.根據(jù)圖8可知,由于滾動(dòng)軸承的故障輕微,故障沖擊淹沒(méi)在背景噪音中,表現(xiàn)不明顯.

    圖7 滾動(dòng)軸承試驗(yàn)臺(tái)布局Fig.7 Experimental setup of rolling element bearing

    圖8 實(shí)測(cè)信號(hào)的時(shí)域波形圖Fig.8 Raw vibration signal in time domain

    圖9 本文算法獲取的轉(zhuǎn)頻及其與傳統(tǒng)算法的綜合比對(duì)Fig.9 Calculated IRF by proposed algorithm and comparison with traditional algorithms

    根據(jù)基于瞬時(shí)故障特征頻率的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法,從未降噪的滾動(dòng)軸承實(shí)測(cè)信號(hào)中獲取的瞬時(shí)故障特征頻率趨勢(shì)線,使其與故障特征系數(shù)Co=2.548相除得到的滾動(dòng)軸承的瞬時(shí)轉(zhuǎn)頻,以及通過(guò)轉(zhuǎn)速計(jì)實(shí)測(cè)的瞬時(shí)轉(zhuǎn)頻分別如圖9(a)所示.其中,實(shí)測(cè)的瞬時(shí)轉(zhuǎn)頻首先根據(jù)轉(zhuǎn)速計(jì)測(cè)取的脈沖信號(hào)計(jì)算出離散的轉(zhuǎn)速點(diǎn),進(jìn)而通過(guò)曲線擬合所得到.實(shí)測(cè)的瞬時(shí)轉(zhuǎn)頻直接由硬件測(cè)量計(jì)算而來(lái),與振動(dòng)信號(hào)具有同步性,因此理論上可以作為驗(yàn)證本文算法的標(biāo)準(zhǔn).根據(jù)圖9(a)易知,估計(jì)的瞬時(shí)轉(zhuǎn)頻有明顯的誤差,其中最大差值約為46.10 Hz,平均差值為28.05 Hz.圖9(b)表示根據(jù)互相關(guān)系數(shù)-峭度準(zhǔn)則進(jìn)行EMD降噪后得到的瞬時(shí)故障特征頻率趨勢(shì)線、計(jì)算得到的滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻以及實(shí)測(cè)的瞬時(shí)轉(zhuǎn)頻. 可知,估計(jì)的瞬時(shí)轉(zhuǎn)頻與實(shí)測(cè)值的變化趨勢(shì)不符,出現(xiàn)了明顯失真現(xiàn)象,沒(méi)有參考價(jià)值.利用本文方法對(duì)滾動(dòng)軸承實(shí)測(cè)信號(hào)進(jìn)行降噪后獲取的瞬時(shí)故障特征頻率趨勢(shì)線、計(jì)算得到的瞬時(shí)轉(zhuǎn)頻以及實(shí)測(cè)的瞬時(shí)轉(zhuǎn)頻分別如圖9(c)所示,估計(jì)的瞬時(shí)轉(zhuǎn)頻準(zhǔn)確的反應(yīng)了實(shí)測(cè)瞬時(shí)轉(zhuǎn)頻的變化趨勢(shì),數(shù)值與實(shí)測(cè)值十分接近,最大差值約為6.30 Hz,平均差值約為2.60 Hz.圖9(d)表示實(shí)測(cè)瞬時(shí)轉(zhuǎn)頻以及分別利用3種方法估計(jì)得到的瞬時(shí)轉(zhuǎn)頻,其中,方塊表示經(jīng)EMD軟閾值降噪后估計(jì)的瞬時(shí)轉(zhuǎn)頻,虛線表示經(jīng)過(guò)EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪后獲取的瞬時(shí)轉(zhuǎn)頻,圓圈表示未降噪直接從原始信號(hào)中估計(jì)的瞬時(shí)轉(zhuǎn)頻,實(shí)線為實(shí)測(cè)瞬時(shí)轉(zhuǎn)頻.通過(guò)對(duì)滾動(dòng)軸承實(shí)測(cè)振動(dòng)信號(hào)分析進(jìn)一步說(shuō)明:EMD軟閾值降噪算法優(yōu)于EMD互相關(guān)系數(shù)-峭度準(zhǔn)則降噪算法,避免了重構(gòu)信號(hào)的失真,具有較好的降噪效果,適用于變轉(zhuǎn)速條件下故障軸承的瞬時(shí)轉(zhuǎn)頻估計(jì).

    4結(jié)語(yǔ)

    對(duì)仿真信號(hào)分析的結(jié)果表明,該方法能夠保證信號(hào)完整性的同時(shí),突出了故障信息,降低了噪聲干擾,進(jìn)而提高了瞬時(shí)轉(zhuǎn)頻估計(jì)的準(zhǔn)確性.本文對(duì)早期故障滾動(dòng)軸承變轉(zhuǎn)速條件下的實(shí)測(cè)信號(hào)進(jìn)行了分析,采用EMD軟閾值降噪后估計(jì)的瞬時(shí)轉(zhuǎn)頻的誤差明顯減小,進(jìn)一步驗(yàn)證了本文方法在滾動(dòng)軸承瞬時(shí)轉(zhuǎn)頻估計(jì)方法中信號(hào)預(yù)處理的有效性.因此,該方法在低信噪比條件下的滾動(dòng)軸承轉(zhuǎn)頻估計(jì)方面具有較好的應(yīng)用前景.

    參考文獻(xiàn)(References):

    [1] MING Z, JING L, WANG X F, et al. A tacho-less order tracking technique for large speed variations[J]. Mechanical Systems and Signal Processing, 2013, 40(1): 76-90.

    [2] BONNARDOT F,BADAOUI E M, RANDALL R B, et al. Use of the acceleration signal of a gearbox in order to perform angular resampling (with limited speed fluctuation)[J]. Mechanical Systems and Signal Processing, 2005, 19 (4): 766-785.

    [3] URBANEK J, BARSZCZ T, ANTONI J. A two-step procedure for estimation of instantaneous rotational speed with large fluctuations [J]. Mechanical System and Signal Processing, 2013, 38:96-102.

    [4] COMBET F, ZIMROZ R. A new method for the estimation of the instantaneous speed relative fluctuation in a vibration signal based on the short time scale transform[J].Mechanical System and Signal Processing, 2009,23: 1382-1397.

    [5] ZHAO D Z, LI J Y, CHENG W D. Feature extraction of faulty rolling element bearing under variable rotational speed and gear interferences conditions[J]. Shock and Vibration, 2015,2015:1-9.

    [6] BOASHASH B. Interpreting and estimating the instantaneous frequency of a signal-Part1: fundamentals [J]. Proceedings of the IEEE, 1992, 80(4):520-538.

    [7] BOASHASH B. Interpreting and estimating the instantaneous frequency of a signal-Part 2: algorithms and applications [J]. Proceedings of the IEEE, 1992, 80(4):540-568. [8] 郭瑜, 秦樹(shù)人, 湯寶平, 等.基于瞬時(shí)頻率估計(jì)的旋轉(zhuǎn)機(jī)械階比跟蹤[J]. 機(jī)械工程學(xué)報(bào), 2003,39(3):32-36.

    GUO Yu, QIN Shu-ren, TANG Bao-ping, et al. Order tracking of rotating machinery based on instantaneous frequency estimation [J].Chinese Journal of Mechanical Engineering, 2003, 39(3):32-36.

    [9] 趙曉平,趙秀莉,侯榮濤,等.一種新的旋轉(zhuǎn)機(jī)械升降速階段振動(dòng)信號(hào)的瞬時(shí)頻率估計(jì)算法[J].機(jī)械工程學(xué)報(bào),2011,47(7):103-108.

    ZHAO Xiao-ping, ZHAO Xiu-li, HOU Rong-tao, et al.A new method for instantaneous frequency estimation of run-up or run-down vibration signal for rotating machinery[J]. Chinese Journal of Mechanical Engineering, 2011, 47(7):103-108.

    [10] WANG T Y, LIANG M, LI J Y, et al. Rolling element bearing fault diagnosis via fault characteristic order (FCO)analysis[J].Mechanical System and Signal Processing, 2014, 45(1): 139-153.[11] DONG G M, CHEN J, ZHAO F G, A frequency-shifted bispectrum for rolling element bearing diagnosis [J]. Journal of Sound and Vibration, 2015, 339:396-418.

    [12] WANG J, HE Q B, KONG F R. An improved multiscale noise tuning of stochastic resonance for identifying multiple transient faults in rolling element bearings[J]. Journal of Sound and Vibration, 2014, 333(26):7401-7421.

    [13] 馮輔周,司愛(ài)威,饒國(guó)強(qiáng),等. 基于小波相關(guān)排列熵的軸承早期故障診斷技術(shù)[J]. 機(jī)械工程學(xué)報(bào),2012,48(13): 73-79.

    FENG Fu-zhou, SI Ai-wei, RAO Guo-qiang, et al. Early fault diagnosis technology for bearing based on wavelet correlation permutation entropy [J]. Journal of Mechanical Engineering, 2012, 48(13): 73-79.

    [14] 王仲生, 姜洪開(kāi), 徐一艷. 發(fā)動(dòng)機(jī)轉(zhuǎn)子系統(tǒng)早期故障智能診斷[J]. 航空學(xué)報(bào),2009,30(2):242-246.

    WANG Zhong-sheng, JIANG Hong-kai, XU Yi-yan. Early fault intelligent diagnosis of aero-engine rotor systems [J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(2):242-246.

    [15] 胥永剛, 馬海龍, 付勝, 等. 機(jī)電設(shè)備早期故障微弱信號(hào)的非線性檢測(cè)方法及工程應(yīng)用[J]. 振動(dòng)工程學(xué)報(bào),2011,24(5):529-538.

    XU Yong-gang, MA Hai-long, FU Sheng, et al. Theory and applications of weak signal non-linear detection method for incipient fault diagnosis of mechanical equipments [J]. Journal of Vibration Engineering, 2011, 24(5):529-538.

    [16] 蘇文勝,王奉濤,張志新,等.EMD降噪和譜峭度法在滾動(dòng)軸承早期故障診斷中的應(yīng)用[J].振動(dòng)與沖擊,2010,29(3): 18-21.SU Wen-sheng, WANG Feng-tao, ZHANG Zhi-xin, et al. Application of EMD denoising and spectral kurtosis in early fault diagnosis of rolling element bearing [J]. Journal of Vibration and Shock, 2010, 29(3):18-21.

    [17] 彭暢,柏林,謝小亮. 基于EEMD、度量因子和快速譜峭度圖的滾動(dòng)軸承故障診斷方法[J].振動(dòng)與沖擊2012, 31(20):144-146.

    PENG Chang, BO Lin, XIE Xiao-liang. Fault diagnosis method of rolling element bearing based on EEMD, measure-factor and fast kurtogram[J].Journal of Vibration and Shock, 2012, 31(20): 144-146. [18] 唐貴基,王曉龍.基于EEMD降噪和1.5維能量譜的滾動(dòng)軸承故障診斷研究[J].振動(dòng)與沖擊, 2014,33(1):6-10.

    TANG Gui-ji, WANG Xiao-long. Fault diagnosis for roller bearings based on EEMD de-noising and 1.5-dimensional energy spectrum [J]. Journal of Vibration and Shock, 2014, 33(1):6-10.

    [19] 程正興,楊守志,馮曉霞.小波分析的理論、算法、進(jìn)展與應(yīng)用[M].北京:國(guó)防工業(yè)出版社,2007:332-335.

    [20] BADOUR A F, SUNAR M, CHEDED L. Vibration analysis of rotating machinery using time-frequency analysis and wavelet techniques [J].Mechanical Systems and Signal Processing,2011,25(6): 2083-2101.

    [21] 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.

    [22] DONOHO D L, JOHNSTONE I M. Threshold selection for wavelet shrinkage of noisy data [J].Engineering in Medicine and Biology Society, 1994, 1:A24-A25.

    [23] DONOHO D L. Denoising by soft-thresholding[J]. IEEE Transaction on Information Theory, 1995, 41(3):613- 627.

    [24] DONOHO D L, JOHNSTONE I. Wavelet shrinkage asymptotic [J].Journal of Royal Statistical Society, 1995, 57(2):301-369.

    [25] 邵忍平,曹精明,李永龍.基于EMD小波閾值去噪和時(shí)頻分析的齒輪故障模式識(shí)別與診斷[J].振動(dòng)與沖擊,2012,31(8):96-101.

    SHAO Ren-ping, CAO Jing-ming, LI Yong-long. Gear fault pattern identification and diagnosis using time- frequency analysis and wavelet threshold de-noising based on EMD [J].Journal of Vibration and Shock, 2012, 31(8):96-101.

    [26] 王天楊, 李建勇, 程衛(wèi)東.基于改進(jìn)的自適應(yīng)噪聲消除和故障特征階比譜的齒輪噪源干擾下變轉(zhuǎn)速滾動(dòng)軸承故障診斷[J].振動(dòng)與沖擊, 2014,33(18):7-13.

    WANG Tian-yang, LI Jian-yong, CHENG Wei-dong. Fault diagnosis of rolling bearing under a variable rotational speed and gear vibration noise based on revised ANC algorithm and FCO spectrum [J]. Journal of Vibration and Shock, 2014, 33 (18):7-13.

    DOI:10.3785/j.issn.1008-973X.2016.03.005

    收稿日期:2015-04-10.

    基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(51275030).

    作者簡(jiǎn)介:程衛(wèi)東(1967-), 男,副教授,從事制造裝備智能測(cè)控與故障診斷教學(xué)科研工作.ORCID:0000-0001-5085-4758. E-mail: wdcheng@bjtu.edu.cn

    中圖分類號(hào):TH 113.1

    文獻(xiàn)標(biāo)志碼:A

    文章編號(hào):1008-973X(2016)03-08-0428

    EMD soft-thresholding denoising algorithm for rolling element bearing rotational frequency estimation

    CHENG Wei-dong, ZHAO De-zun

    (SchoolofMechanicalElectronicandControlEngineering,BeijngJiaotongUniversity,Beijing100044,China)

    Abstract:A method based on the empirical mode decomposition (EMD) and wavelet shrinkage was proposed due to the shortcomings of the EMD cross-correlation coefficient and kurtosis criterion denoising and wavelet shrinkage. The method consists of four main steps: (i) IMFs were obtained by decomposing raw signal, (ii) The cross-correlation coefficient between IMFs and the raw signal, and the kurtosis values of the IMFs were calculated, (iii) IMFs with the noise were selected and the false IMFs were removed, (vi) the noise of the selected IMFs was removed by the soft-thresholding denoising method, and then the signal with the rest of the IMFs was reconstructed. The proposed method was tested based on both the simulated and experimental bearing vibration signals. Results show that, compared with EMD cross-correlation coefficient and kurtosis criterion denoising, the method of EMD soft-thresholding denoising can ensure the integrity of the signal, highlight fault features and reduce the error of the instantaneous rotational frequency (IRF).

    Key words:empirical mode decomposition (EMD); soft-thresholding denoising; rolling element bearing; signal preprocessing; instantaneous rotational frequency (IRF) estimation

    猜你喜歡
    振動(dòng)故障信號(hào)
    振動(dòng)的思考
    信號(hào)
    鴨綠江(2021年35期)2021-04-19 12:24:18
    完形填空二則
    振動(dòng)與頻率
    故障一點(diǎn)通
    基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
    電子制作(2018年11期)2018-08-04 03:25:42
    中立型Emden-Fowler微分方程的振動(dòng)性
    奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
    基于LabVIEW的力加載信號(hào)采集與PID控制
    故障一點(diǎn)通
    亚洲国产欧美在线一区| 国产成人精品在线电影| 永久网站在线| 成年女人毛片免费观看观看9 | 精品亚洲成国产av| 性高湖久久久久久久久免费观看| 久久久久精品久久久久真实原创| 国产成人av激情在线播放| 2021少妇久久久久久久久久久| 妹子高潮喷水视频| 日韩,欧美,国产一区二区三区| 欧美精品人与动牲交sv欧美| av又黄又爽大尺度在线免费看| av国产久精品久网站免费入址| 激情视频va一区二区三区| 久久这里有精品视频免费| 亚洲欧美精品自产自拍| 国产欧美日韩综合在线一区二区| 日日摸夜夜添夜夜爱| 国产无遮挡羞羞视频在线观看| 另类亚洲欧美激情| 大话2 男鬼变身卡| 免费不卡的大黄色大毛片视频在线观看| 午夜福利视频精品| 乱人伦中国视频| 亚洲中文av在线| 免费观看a级毛片全部| 中国国产av一级| 国产黄色视频一区二区在线观看| 女人被躁到高潮嗷嗷叫费观| 好男人视频免费观看在线| 国产又色又爽无遮挡免| 美国免费a级毛片| 久久午夜福利片| 毛片一级片免费看久久久久| 9热在线视频观看99| av不卡在线播放| 高清欧美精品videossex| av在线播放精品| 男人添女人高潮全过程视频| 91成人精品电影| 99re6热这里在线精品视频| 久久精品国产亚洲av高清一级| 看非洲黑人一级黄片| 多毛熟女@视频| 久久人人爽人人片av| 国产成人精品婷婷| 日日爽夜夜爽网站| 久久国产亚洲av麻豆专区| 好男人视频免费观看在线| 美女国产高潮福利片在线看| 精品少妇一区二区三区视频日本电影 | 久久久久国产一级毛片高清牌| 亚洲一码二码三码区别大吗| 电影成人av| a级毛片在线看网站| 中文欧美无线码| 亚洲成国产人片在线观看| av一本久久久久| 天堂中文最新版在线下载| 国产精品二区激情视频| 午夜免费男女啪啪视频观看| 青春草视频在线免费观看| 久久精品夜色国产| 亚洲国产日韩一区二区| 一个人免费看片子| 十分钟在线观看高清视频www| 搡女人真爽免费视频火全软件| 精品一区二区免费观看| 18禁裸乳无遮挡动漫免费视频| 免费在线观看黄色视频的| 日韩三级伦理在线观看| 熟女少妇亚洲综合色aaa.| 大码成人一级视频| 久久久久久久久免费视频了| 久久久久人妻精品一区果冻| 午夜日韩欧美国产| 亚洲欧美精品综合一区二区三区 | 午夜精品国产一区二区电影| 女人被躁到高潮嗷嗷叫费观| 国产成人欧美| 亚洲国产看品久久| kizo精华| 国产xxxxx性猛交| 国产黄色视频一区二区在线观看| 毛片一级片免费看久久久久| 91在线精品国自产拍蜜月| 综合色丁香网| 久久 成人 亚洲| 啦啦啦视频在线资源免费观看| av又黄又爽大尺度在线免费看| 2022亚洲国产成人精品| 欧美xxⅹ黑人| 男女免费视频国产| 91aial.com中文字幕在线观看| 国产亚洲午夜精品一区二区久久| 中文字幕精品免费在线观看视频| 久久久久久免费高清国产稀缺| 亚洲美女搞黄在线观看| tube8黄色片| 九草在线视频观看| 高清黄色对白视频在线免费看| 最黄视频免费看| 精品久久蜜臀av无| 丝袜美足系列| 国精品久久久久久国模美| 久久99热这里只频精品6学生| 亚洲精品av麻豆狂野| 国语对白做爰xxxⅹ性视频网站| 秋霞在线观看毛片| 18禁观看日本| 精品一品国产午夜福利视频| 90打野战视频偷拍视频| 永久免费av网站大全| 国产熟女午夜一区二区三区| 大码成人一级视频| 亚洲色图综合在线观看| 日本欧美国产在线视频| 女人被躁到高潮嗷嗷叫费观| 精品人妻在线不人妻| 色94色欧美一区二区| a级毛片在线看网站| 日韩中文字幕视频在线看片| 男男h啪啪无遮挡| 一二三四在线观看免费中文在| 免费少妇av软件| 纵有疾风起免费观看全集完整版| 久久影院123| 国产成人欧美| 午夜福利在线观看免费完整高清在| 国产激情久久老熟女| 一二三四在线观看免费中文在| 午夜福利,免费看| 性高湖久久久久久久久免费观看| 亚洲av福利一区| 亚洲欧美一区二区三区国产| 国产精品一区二区在线观看99| 亚洲一级一片aⅴ在线观看| 黄色视频在线播放观看不卡| 亚洲国产精品999| 哪个播放器可以免费观看大片| 美女福利国产在线| 久久狼人影院| 99香蕉大伊视频| 久久午夜福利片| 日韩电影二区| 91aial.com中文字幕在线观看| 一二三四中文在线观看免费高清| 亚洲三级黄色毛片| 五月伊人婷婷丁香| 色婷婷久久久亚洲欧美| 人妻一区二区av| 26uuu在线亚洲综合色| 一级毛片我不卡| 男女无遮挡免费网站观看| 超色免费av| 高清不卡的av网站| 秋霞伦理黄片| 日本av免费视频播放| 一二三四中文在线观看免费高清| 久久狼人影院| 亚洲av成人精品一二三区| 又黄又粗又硬又大视频| 1024香蕉在线观看| 亚洲av.av天堂| 美女国产高潮福利片在线看| videosex国产| 男女边吃奶边做爰视频| 亚洲精品国产av蜜桃| 在线天堂最新版资源| www.熟女人妻精品国产| 高清在线视频一区二区三区| 国产精品国产三级专区第一集| 一二三四中文在线观看免费高清| 少妇的逼水好多| 国产欧美日韩综合在线一区二区| 少妇 在线观看| 国产在线免费精品| 日本猛色少妇xxxxx猛交久久| 国产精品久久久久久精品古装| 国产高清不卡午夜福利| 精品第一国产精品| 久久久久久久亚洲中文字幕| 亚洲天堂av无毛| 亚洲欧美清纯卡通| 99国产精品免费福利视频| 色婷婷av一区二区三区视频| 亚洲精品中文字幕在线视频| 亚洲成人av在线免费| 国产亚洲午夜精品一区二区久久| 涩涩av久久男人的天堂| 香蕉国产在线看| 夜夜骑夜夜射夜夜干| 1024香蕉在线观看| 成年女人毛片免费观看观看9 | 国产熟女午夜一区二区三区| 国产精品免费视频内射| 狠狠婷婷综合久久久久久88av| 久久ye,这里只有精品| 亚洲av综合色区一区| 精品国产露脸久久av麻豆| 久久久精品94久久精品| 国产亚洲精品第一综合不卡| av免费观看日本| 春色校园在线视频观看| 精品人妻熟女毛片av久久网站| 日韩三级伦理在线观看| 另类精品久久| 男女啪啪激烈高潮av片| 美国免费a级毛片| 免费在线观看完整版高清| 成人国产av品久久久| 久久久久久人妻| 大话2 男鬼变身卡| 最近中文字幕高清免费大全6| 少妇精品久久久久久久| 亚洲精品,欧美精品| 夜夜骑夜夜射夜夜干| 人妻一区二区av| 久久久久视频综合| 亚洲图色成人| 亚洲av国产av综合av卡| 欧美日韩国产mv在线观看视频| 国产成人午夜福利电影在线观看| 性少妇av在线| 999精品在线视频| 久久人妻熟女aⅴ| 亚洲情色 制服丝袜| 激情视频va一区二区三区| 亚洲国产看品久久| 久久久亚洲精品成人影院| av卡一久久| 国产淫语在线视频| 精品国产一区二区久久| 国产 精品1| 性高湖久久久久久久久免费观看| 在线观看美女被高潮喷水网站| 波多野结衣av一区二区av| 如日韩欧美国产精品一区二区三区| 国产精品嫩草影院av在线观看| 久久久久精品人妻al黑| 国产高清国产精品国产三级| 精品国产一区二区三区四区第35| 日本vs欧美在线观看视频| 老司机亚洲免费影院| 亚洲三级黄色毛片| 欧美亚洲 丝袜 人妻 在线| 欧美亚洲 丝袜 人妻 在线| 成年人免费黄色播放视频| 国产成人91sexporn| 中文乱码字字幕精品一区二区三区| 国产精品久久久av美女十八| 91午夜精品亚洲一区二区三区| 国产成人免费无遮挡视频| 老汉色∧v一级毛片| 伦理电影大哥的女人| 国产精品秋霞免费鲁丝片| av国产久精品久网站免费入址| 精品亚洲成国产av| 国产有黄有色有爽视频| 日韩中文字幕视频在线看片| 777久久人妻少妇嫩草av网站| 少妇被粗大的猛进出69影院| 黄片无遮挡物在线观看| 美女xxoo啪啪120秒动态图| 99热全是精品| 女性生殖器流出的白浆| 国产男人的电影天堂91| 日韩大片免费观看网站| 免费人妻精品一区二区三区视频| 青青草视频在线视频观看| 国产精品久久久久久精品古装| 女人高潮潮喷娇喘18禁视频| 熟女电影av网| 亚洲美女视频黄频| 色播在线永久视频| 久久精品国产a三级三级三级| 色94色欧美一区二区| 一本色道久久久久久精品综合| 国产成人a∨麻豆精品| 国产爽快片一区二区三区| 日本免费在线观看一区| 国产精品成人在线| 中文字幕最新亚洲高清| 日韩欧美精品免费久久| 狠狠精品人妻久久久久久综合| 午夜福利网站1000一区二区三区| 一区二区三区乱码不卡18| 国产探花极品一区二区| av一本久久久久| 国产色婷婷99| av片东京热男人的天堂| 99久久人妻综合| 国产精品一国产av| 免费看av在线观看网站| 欧美精品一区二区免费开放| 久久人人爽人人片av| 成人毛片60女人毛片免费| 丰满少妇做爰视频| 国产高清不卡午夜福利| 王馨瑶露胸无遮挡在线观看| 精品酒店卫生间| 亚洲欧洲日产国产| 亚洲三级黄色毛片| 午夜激情久久久久久久| 巨乳人妻的诱惑在线观看| 亚洲av成人精品一二三区| 黄色 视频免费看| 日产精品乱码卡一卡2卡三| 日本欧美国产在线视频| 午夜福利在线观看免费完整高清在| 丝袜美足系列| 午夜福利在线免费观看网站| 9191精品国产免费久久| 国产一区二区激情短视频 | 涩涩av久久男人的天堂| 性少妇av在线| 亚洲少妇的诱惑av| 少妇人妻久久综合中文| 女人高潮潮喷娇喘18禁视频| 最新中文字幕久久久久| 国产欧美日韩一区二区三区在线| av免费观看日本| 性色av一级| 久久久久久久久久久久大奶| 亚洲国产精品成人久久小说| 满18在线观看网站| 亚洲国产色片| 一级片免费观看大全| 老司机影院成人| 叶爱在线成人免费视频播放| 亚洲国产最新在线播放| 久久久久久久久久久免费av| 色播在线永久视频| 人体艺术视频欧美日本| 日韩免费高清中文字幕av| 国产成人精品在线电影| 亚洲欧美一区二区三区久久| 成人亚洲欧美一区二区av| 一二三四中文在线观看免费高清| 伦理电影大哥的女人| 91在线精品国自产拍蜜月| 欧美人与善性xxx| 久久热在线av| 精品少妇黑人巨大在线播放| 天堂8中文在线网| 日韩在线高清观看一区二区三区| 黑人欧美特级aaaaaa片| 亚洲男人天堂网一区| 亚洲国产精品一区三区| 尾随美女入室| 青青草视频在线视频观看| 尾随美女入室| 看免费av毛片| 欧美激情极品国产一区二区三区| 91精品伊人久久大香线蕉| 国产在线视频一区二区| 亚洲,一卡二卡三卡| 不卡视频在线观看欧美| 久久人人爽av亚洲精品天堂| 精品人妻熟女毛片av久久网站| 午夜免费鲁丝| 捣出白浆h1v1| 国产免费视频播放在线视频| 大香蕉久久网| 亚洲国产最新在线播放| 欧美xxⅹ黑人| 午夜激情av网站| 国产女主播在线喷水免费视频网站| 五月天丁香电影| 久久99一区二区三区| 欧美国产精品va在线观看不卡| 成年美女黄网站色视频大全免费| 夜夜骑夜夜射夜夜干| 中文字幕人妻丝袜制服| 99久国产av精品国产电影| 日本-黄色视频高清免费观看| 我要看黄色一级片免费的| 叶爱在线成人免费视频播放| www日本在线高清视频| 日本-黄色视频高清免费观看| av不卡在线播放| 最近最新中文字幕大全免费视频 | 久久97久久精品| 黄色一级大片看看| 国产精品久久久久久av不卡| 免费高清在线观看日韩| 2021少妇久久久久久久久久久| 大话2 男鬼变身卡| 久久精品夜色国产| 日本av免费视频播放| 国产成人欧美| 国产麻豆69| 极品少妇高潮喷水抽搐| 成人国产av品久久久| 欧美老熟妇乱子伦牲交| 国产精品女同一区二区软件| 成年人免费黄色播放视频| 高清不卡的av网站| 99国产综合亚洲精品| 国产在视频线精品| 国产免费福利视频在线观看| 一级a爱视频在线免费观看| 一级,二级,三级黄色视频| 最近2019中文字幕mv第一页| 在线观看免费日韩欧美大片| 国产成人欧美| 菩萨蛮人人尽说江南好唐韦庄| 高清在线视频一区二区三区| 啦啦啦在线观看免费高清www| 亚洲三级黄色毛片| 久久午夜综合久久蜜桃| 国产欧美亚洲国产| 午夜91福利影院| 99国产综合亚洲精品| 亚洲四区av| 母亲3免费完整高清在线观看 | 最近最新中文字幕大全免费视频 | 欧美日韩国产mv在线观看视频| 你懂的网址亚洲精品在线观看| 老司机亚洲免费影院| 天天躁夜夜躁狠狠久久av| 国产乱来视频区| 两性夫妻黄色片| 国产一级毛片在线| 少妇人妻久久综合中文| 亚洲国产看品久久| 丰满迷人的少妇在线观看| 久久精品国产亚洲av涩爱| 看免费av毛片| 精品国产一区二区三区久久久樱花| 免费在线观看视频国产中文字幕亚洲 | 国产亚洲最大av| 欧美人与善性xxx| 少妇猛男粗大的猛烈进出视频| 国产男女超爽视频在线观看| 亚洲av免费高清在线观看| 天天躁夜夜躁狠狠久久av| 国产精品一二三区在线看| 69精品国产乱码久久久| 久久久国产欧美日韩av| 多毛熟女@视频| 国产毛片在线视频| 亚洲婷婷狠狠爱综合网| 亚洲欧美一区二区三区黑人 | 亚洲熟女精品中文字幕| 国产乱人偷精品视频| av在线观看视频网站免费| 美女视频免费永久观看网站| 国产成人91sexporn| 欧美在线黄色| 久久午夜综合久久蜜桃| 伦理电影大哥的女人| 亚洲精品国产av成人精品| 中文天堂在线官网| 在线观看www视频免费| 最新中文字幕久久久久| 巨乳人妻的诱惑在线观看| 伦理电影免费视频| 亚洲精品美女久久av网站| 丰满少妇做爰视频| 男女啪啪激烈高潮av片| 18禁动态无遮挡网站| 日韩伦理黄色片| 韩国高清视频一区二区三区| 欧美bdsm另类| 巨乳人妻的诱惑在线观看| 国产视频首页在线观看| 国产精品.久久久| 超碰成人久久| 一区二区三区精品91| 久久久久久人妻| 成年女人毛片免费观看观看9 | 免费黄色在线免费观看| 男人操女人黄网站| 一个人免费看片子| 久久久久人妻精品一区果冻| 久久久国产精品麻豆| 五月天丁香电影| 大香蕉久久成人网| 十八禁高潮呻吟视频| 国产精品三级大全| 欧美日本中文国产一区发布| 一本久久精品| 亚洲人成电影观看| 亚洲国产精品国产精品| 有码 亚洲区| 午夜福利,免费看| 9191精品国产免费久久| 日韩av不卡免费在线播放| 另类精品久久| 一本—道久久a久久精品蜜桃钙片| 日韩一区二区视频免费看| 天堂8中文在线网| 大陆偷拍与自拍| 亚洲国产色片| 69精品国产乱码久久久| 国产成人精品一,二区| 丰满迷人的少妇在线观看| 夫妻性生交免费视频一级片| 欧美日韩精品网址| 五月伊人婷婷丁香| 精品一区二区三卡| 天美传媒精品一区二区| 熟女电影av网| 一区二区三区激情视频| 99久久中文字幕三级久久日本| 久久精品aⅴ一区二区三区四区 | a级片在线免费高清观看视频| 久久 成人 亚洲| 国产成人精品久久二区二区91 | 新久久久久国产一级毛片| 欧美另类一区| 色94色欧美一区二区| a 毛片基地| 美女主播在线视频| 男女无遮挡免费网站观看| 精品少妇久久久久久888优播| av.在线天堂| 一本久久精品| 免费观看在线日韩| 日本欧美国产在线视频| 国产精品 国内视频| 一区二区三区激情视频| av免费在线看不卡| 国产成人午夜福利电影在线观看| 王馨瑶露胸无遮挡在线观看| 午夜日本视频在线| 欧美成人午夜免费资源| 亚洲欧美清纯卡通| 亚洲av日韩在线播放| 国产成人欧美| 国产成人精品婷婷| 欧美激情极品国产一区二区三区| 亚洲欧美中文字幕日韩二区| 欧美日韩一级在线毛片| 午夜免费观看性视频| 综合色丁香网| 国产福利在线免费观看视频| 中文字幕制服av| 一边亲一边摸免费视频| 五月开心婷婷网| 亚洲国产色片| 日韩不卡一区二区三区视频在线| 成年女人在线观看亚洲视频| 国产免费视频播放在线视频| 男女啪啪激烈高潮av片| 久久精品国产鲁丝片午夜精品| 亚洲欧洲精品一区二区精品久久久 | 久久精品aⅴ一区二区三区四区 | 亚洲精品美女久久久久99蜜臀 | 久久韩国三级中文字幕| 日韩人妻精品一区2区三区| 日韩制服丝袜自拍偷拍| 97在线视频观看| 欧美最新免费一区二区三区| 国产精品国产三级国产专区5o| 一级毛片 在线播放| a级片在线免费高清观看视频| 亚洲欧美清纯卡通| 成人黄色视频免费在线看| 免费黄频网站在线观看国产| 亚洲精品国产色婷婷电影| 亚洲成色77777| 久久人人爽人人片av| 七月丁香在线播放| 99国产精品免费福利视频| 五月天丁香电影| 日韩伦理黄色片| 少妇的逼水好多| 一区二区三区乱码不卡18| 欧美精品一区二区免费开放| 日韩视频在线欧美| 一二三四中文在线观看免费高清| 麻豆av在线久日| 精品亚洲成国产av| 日本爱情动作片www.在线观看| 色哟哟·www| 国产精品.久久久| 国产精品熟女久久久久浪| 欧美精品人与动牲交sv欧美| 超碰97精品在线观看| 久久久久精品久久久久真实原创| 在线观看人妻少妇| 国产男女内射视频| 老汉色av国产亚洲站长工具| 人妻少妇偷人精品九色| 亚洲av成人精品一二三区| 国产精品不卡视频一区二区| 一本大道久久a久久精品| 五月伊人婷婷丁香| 少妇被粗大的猛进出69影院| 亚洲欧美清纯卡通| 日韩中字成人| 免费看不卡的av| 中文字幕av电影在线播放| 天天躁夜夜躁狠狠久久av| av在线app专区| 我的亚洲天堂| 精品人妻偷拍中文字幕| 久久精品国产a三级三级三级| 国产精品女同一区二区软件| 纵有疾风起免费观看全集完整版| 99久久人妻综合| 国产综合精华液| 久久99一区二区三区| av线在线观看网站| 久久人人爽av亚洲精品天堂| 人成视频在线观看免费观看| 十八禁高潮呻吟视频| 校园人妻丝袜中文字幕| 日韩av不卡免费在线播放| 天天躁狠狠躁夜夜躁狠狠躁| 久久青草综合色| 91午夜精品亚洲一区二区三区|