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

    基于MSET重構(gòu)模型整體優(yōu)化的軸承性能退化評(píng)估方法

    2023-09-05 02:11:10劉楊遠(yuǎn)吳榮真承志恒顏秋宏
    振動(dòng)與沖擊 2023年16期
    關(guān)鍵詞:軸承向量觀測(cè)

    張 龍, 劉楊遠(yuǎn), 吳榮真, 王 良, 承志恒, 顏秋宏

    (華東交通大學(xué) 載運(yùn)工具與裝備教育部重點(diǎn)實(shí)驗(yàn)室,南昌 330013)

    軸承的性能從健康到失效是一個(gè)漸變的過程,性能退化評(píng)估不僅是對(duì)該過程實(shí)現(xiàn)定量評(píng)估的有效手段,而且是實(shí)現(xiàn)軸承故障預(yù)診斷的前提與基礎(chǔ),更是發(fā)揮狀態(tài)預(yù)測(cè)維修潛力的重要一環(huán)[1]。

    軸承性能退化評(píng)估的一般流程為特征提取、特征選擇、建立劣化指標(biāo)、比較實(shí)測(cè)樣本與無故障基準(zhǔn)之間相似性[2]等幾個(gè)方面。其中比較度量相似性方法又可分為基于邊界的檢測(cè)方法、基于概率密度估計(jì)和基于模型重構(gòu)的檢測(cè)方法3個(gè)方面。典型的邊界檢測(cè)方法有支持向量機(jī)[3]和支持向量數(shù)據(jù)描述[4],該類方法通過無故障訓(xùn)練數(shù)據(jù)集來創(chuàng)建邊界,使用待評(píng)估數(shù)據(jù)與邊界的距離確定其相對(duì)于無故障狀態(tài)的類隸屬。但是基于邊界的異常檢測(cè)方法問題在于距離指標(biāo)的上限難以界定,從而難于從退化指標(biāo)中定量推斷性能退化程度?;诟怕拭芏裙烙?jì)的方法一般使用高斯混合模型[5]和隱馬爾科夫模型[6],該類方法首先對(duì)訓(xùn)練期間的無故障數(shù)據(jù)進(jìn)行概率密度建模,之后將待測(cè)數(shù)據(jù)輸入到訓(xùn)練好的模型中,樣本若落入低密度區(qū)域則不屬于無故障樣本。但是,高斯混合模型假設(shè)無故障數(shù)據(jù)的概率密度分布可由多個(gè)高斯分布組合得到,其獨(dú)立高斯函數(shù)數(shù)量難于精確估計(jì),不能保證與實(shí)際情況完全相符;而對(duì)于復(fù)雜的系統(tǒng)而言,隱馬爾科夫模型利用指數(shù)表達(dá)形式模擬健康狀態(tài)區(qū)間,該策略與現(xiàn)實(shí)情況相違背,往往導(dǎo)致精度較低。

    基于模型重構(gòu)的方法通過映射目標(biāo)數(shù)據(jù)集,將待測(cè)數(shù)據(jù)輸入到經(jīng)過無故障數(shù)據(jù)訓(xùn)練的模型進(jìn)行映射,其重建誤差大小(例如歐幾里德距離)與退化性能程度密切相關(guān)。多元狀態(tài)估計(jì)(multivariate state estimation technique, MSET)就是一個(gè)典型的基于模型重構(gòu)的性能退化評(píng)估方法,本質(zhì)上是一種非線性的多元預(yù)測(cè)診斷技術(shù),也可稱為一種非參數(shù)建模方法。MSET通過對(duì)健康狀態(tài)基準(zhǔn)進(jìn)行相似性建模,然后輸入待測(cè)的觀測(cè)向量進(jìn)行模型重構(gòu),比較度量前后差異性從而判斷目標(biāo)的性能退化程度。例如Wang等[7]將最大條件互信息(conditional mutual information, CMI)特征選擇算法應(yīng)用于MSET的訓(xùn)練數(shù)據(jù),用K近鄰算法(K-nearest neighbor,KNN)動(dòng)態(tài)存儲(chǔ)歷史記憶矩陣,成功實(shí)現(xiàn)了風(fēng)電機(jī)組性能狀態(tài)檢測(cè)。MSET重構(gòu)模型相比于基于邊界和概率密度估計(jì)的性能退化方法優(yōu)勢(shì)明顯,其模型簡(jiǎn)單且物理意義明確,能及時(shí)準(zhǔn)確地反映設(shè)備性能退化趨勢(shì)。

    在性能退化特征提取和選擇時(shí),常用的方法有時(shí)域、頻域特征提取、小波特征提取和時(shí)序模型分析等等。崔凱等[8]將加速度有效值、峭度值兩個(gè)時(shí)域指標(biāo)結(jié)合在一起聯(lián)合判斷風(fēng)機(jī)軸承運(yùn)行狀態(tài)。Zhou等[9]采用軸承信號(hào)小波包能量熵作為原始特征,然后將全壽命數(shù)據(jù)輸入到利用健康狀態(tài)特征建立好的RBF網(wǎng)絡(luò)模型中,從而得到性能退化指標(biāo)。Cong等[10]利用基于時(shí)序模型中的自回歸模型(AR模型)和Kolmogorov-Smirnov檢驗(yàn)相似概率值來反映滾動(dòng)軸承性能劣化的趨勢(shì)。但上述傳統(tǒng)的單個(gè)域指標(biāo)特征提取效果單一,對(duì)狀態(tài)變化感知不明顯,無法有效囊括軸承性能狀態(tài)且不利于早期故障的發(fā)現(xiàn)。

    許多研究將多域特征指標(biāo)融合進(jìn)行性能退化評(píng)估,例如Liu等[11]從軸承振動(dòng)信號(hào)中提取多域統(tǒng)計(jì)特征,然后利用核聯(lián)合特征矩陣近似對(duì)角化算法將其融合為敏感特征,最后計(jì)算健康數(shù)據(jù)和當(dāng)前數(shù)據(jù)的類間和類內(nèi)離散度作為性能退化指標(biāo)。對(duì)于將多域特征組合的方法,雖然相較于單域指標(biāo)能更細(xì)致地描繪軸承退化過程,但是高維特征存在信息冗余且計(jì)算量大的問題,從而降低性能退化評(píng)估模型的性能。Dong等[12]將時(shí)域、頻域、時(shí)頻域特征利用主成分分析(principal component analysis,PCA)降維后作為性能退化指標(biāo)。叢華等[13]采用遺傳算法(genetic algorithm,GA)優(yōu)化提取特征,然后計(jì)算特征與支持向量數(shù)據(jù)描述超球體模型間的距離從而評(píng)估軸承的性能退化程度。然而,①PCA降維與其他優(yōu)化算法不同,并非直接從高維特征中挑選出具有代表性的低維特征,而是利用空間轉(zhuǎn)換將高維空間映射成低維空間,導(dǎo)致破壞了原有數(shù)據(jù)的結(jié)構(gòu),改變了原有數(shù)據(jù)的物理意義;②利用PCA或者GA單純對(duì)輸入特征進(jìn)行降維或優(yōu)化雖然一定程度上減少了冗余信息,但傳統(tǒng)的時(shí)頻域特征降維后的指標(biāo)不一定適配所需模型,匹配度有待進(jìn)一步考量。而目前做聯(lián)合降維性能退化評(píng)估的長(zhǎng)短期記憶網(wǎng)絡(luò)[14-15](long short-time memory, LSTM)等參數(shù)過于復(fù)雜,時(shí)間跨度大且網(wǎng)絡(luò)較深,導(dǎo)致計(jì)算量大且耗時(shí)費(fèi)力。

    綜上,本文提出一種基于MSET重構(gòu)模型整體優(yōu)化的軸承性能退化評(píng)估方法。首先,基于軸承振動(dòng)信號(hào)提取m組時(shí)域和頻域特征、n階AR模型系數(shù)和k維三層小波包Renyi 熵歸一化后組成m+n+k維多域特征向量指標(biāo),并將得到健康狀態(tài)的高維特征向量構(gòu)建MSET模型的歷史觀測(cè)矩陣;然后,利用遺傳算法對(duì)軸承高維特征向量和多元狀態(tài)估計(jì)中的歷史記憶矩陣的觀測(cè)向量個(gè)數(shù)和維度進(jìn)行同步聯(lián)合優(yōu)化,從而得到低維特征向量和歷史記憶矩陣觀測(cè)向量個(gè)數(shù)與該低維特征向量最匹配的MSET 模型,實(shí)現(xiàn)了特征優(yōu)選和滾動(dòng)軸承性能評(píng)估模型的整體自適應(yīng)優(yōu)化;最后,構(gòu)建余弦相似度劣化指標(biāo)實(shí)現(xiàn)軸承性能退化評(píng)估。

    1 理論基礎(chǔ)

    1.1 小波包Renyi熵

    小波包相較于小波而言,在信號(hào)高頻部分分解得更加精細(xì),能更好地刻畫信號(hào)的突變性和非平穩(wěn)性。其本質(zhì)就是集合數(shù)個(gè)高低通濾波器從而將頻譜劃分多個(gè)層次,得到相應(yīng)的高低頻成分,更好地表達(dá)信號(hào)。其分解所得的頻帶之間沒有疏漏和冗余,相互獨(dú)立。

    小波包分解表達(dá)式[16]如下

    (1)

    小波包重構(gòu)表達(dá)式如下

    (2)

    小波包變換實(shí)際上就是小波基函數(shù)與信號(hào)進(jìn)行卷積的過程。卷積在時(shí)域中的表現(xiàn)形式為移動(dòng)的內(nèi)積,從而使基函數(shù)和信號(hào)之間相似部分放大,反之則被抑制,因此小波基函數(shù)的選擇將直接影響最后結(jié)果。信號(hào)進(jìn)行j層小波包分解后將分成2j個(gè)子帶信號(hào)。Nikolaou等[17]提出選擇3層或4層小波包分解時(shí),對(duì)所得到的信號(hào)分析結(jié)果影響不大。因此,為提高計(jì)算效率,本文采用db5小波基函數(shù),對(duì)信號(hào)進(jìn)行3層小波包分解。

    當(dāng)設(shè)備發(fā)生故障后,信號(hào)的不確定信息及能量分布將較正常運(yùn)轉(zhuǎn)時(shí)發(fā)生改變,不同的故障類型以及不同程度的故障其能量分布更是不同。信息熵具有表征信號(hào)復(fù)雜程度的能力,信號(hào)隨機(jī)性越高、越復(fù)雜時(shí)信息熵越大,反之信號(hào)越規(guī)則信息熵越小。

    信息熵計(jì)算公式如下

    (3)

    Renyi熵作為信息熵概念的擴(kuò)展,具有表征信號(hào)在時(shí)頻分布特征的能力[18]。Renyi 熵較信息熵而言對(duì)信號(hào)細(xì)微變化更加敏感,尤其是在端點(diǎn)處,能夠捕捉到隨機(jī)變量在概率分布中的差異并進(jìn)行放大。Renyi熵計(jì)算公式如下

    (4)

    式中,α≠1且α≥0。當(dāng)α=0時(shí)Hα(x)取最大值,Markel等[19]對(duì)Renyi熵的α進(jìn)行了研究,表明α=0.5為結(jié)果最優(yōu)值,故本文α選為0.5。

    1.2 AR模型

    自回歸模型(autoregressive model, AR)是一種隨機(jī)信號(hào)參數(shù)化建模方法,常用于處理時(shí)間序列。AR模型利用有限的知識(shí)和自身的參數(shù)對(duì)前面多個(gè)觀測(cè)值進(jìn)行學(xué)習(xí),從而對(duì)下一時(shí)刻的輸出進(jìn)行預(yù)測(cè),其對(duì)系統(tǒng)狀態(tài)變化極其敏感,且能表征系統(tǒng)狀態(tài)特征。

    取時(shí)間序列y(t),AR模型的分析階數(shù)為p,則關(guān)于時(shí)間序列的p階AR模型可以表示為

    (5)

    式中:L(t)為AR模型的殘差;αj為第j項(xiàng)的系數(shù)。

    殘差L(t)為真實(shí)值和預(yù)測(cè)值兩者間的誤差,對(duì)于平穩(wěn)信號(hào)L(t)主要表征白噪聲。所以AR模型實(shí)際上就是調(diào)整模型參數(shù)使殘差更接近高斯白噪聲。本文選用最常用的最小二乘法來進(jìn)行AR譜估計(jì)以及模型參數(shù)的選擇,根據(jù)貝葉斯信息準(zhǔn)則(Bayesian information criterion,BIC)確定模型階數(shù)??纱笾路譃槿?

    步驟1確定分析階數(shù)p的范圍繼而調(diào)整AR模型,本文p分別取1,2,…,300;

    步驟2通過最小二乘法分別求得各階次下的自回歸參數(shù)αj(j=1,2,…,300),構(gòu)造式(5)所示的AR模型,進(jìn)而得到殘差L(t);

    步驟3根據(jù)各階次殘差L(t)的BIC值最小原則,確定最優(yōu)階數(shù)。

    1.3 多元狀態(tài)估計(jì)

    多元狀態(tài)估計(jì)是一種非線性的非參數(shù)建模方法。MSET通過比較新輸入數(shù)據(jù)與無故障時(shí)的健康數(shù)據(jù),對(duì)各變量數(shù)據(jù)之間的內(nèi)在關(guān)聯(lián)做出判斷,根據(jù)歷史數(shù)據(jù)所獲得的權(quán)值向量等知識(shí)對(duì)實(shí)際運(yùn)行狀態(tài)做出估計(jì),通過殘差比較估計(jì)向量與新輸入的觀測(cè)向量的相似性,從而對(duì)設(shè)備進(jìn)行診斷。

    假設(shè)在某一過程或系統(tǒng)中某個(gè)時(shí)刻tj有n個(gè)互相關(guān)聯(lián)的變量,將其記為觀測(cè)向量X(tj),如式(6)所示

    X(tj)=[x1(tj)x2(tj) …xn(tj)]T

    (6)

    式中,X(tj)為tj時(shí)刻狀態(tài)變量xi的觀測(cè)值。

    首先,通過一定的方法合理地挑選m列正常運(yùn)行狀態(tài)下的觀測(cè)變量即可組成歷史記憶矩陣Dn×m,如式(7)所示,從而完整地表示某一過程或系統(tǒng)正常運(yùn)行狀態(tài)下的整個(gè)動(dòng)態(tài)過程,換言之D即是對(duì)設(shè)備正常運(yùn)行數(shù)據(jù)的特征提取、選擇、學(xué)習(xí)和記憶的過程。

    (7)

    式中:n為在某一時(shí)刻互相關(guān)聯(lián)的狀態(tài)變量有n個(gè);m為取m個(gè)觀測(cè)向量X(tj),即取m個(gè)時(shí)刻。

    然后,假設(shè)在某一過程或系統(tǒng)中某個(gè)時(shí)刻模型的輸入為觀測(cè)向量Xobs,將其輸入歷史記憶矩陣D之后的輸出為估計(jì)向量Xest。Xest是對(duì)輸入Xobs的預(yù)測(cè)值,即對(duì)當(dāng)前狀態(tài)的估計(jì)。對(duì)于任一輸入觀測(cè)向量Xobs,MSET均對(duì)應(yīng)生成一個(gè)m維的權(quán)值向量W

    (8)

    根據(jù)權(quán)值向量W對(duì)D中m個(gè)歷史觀測(cè)向量進(jìn)行線性組合即可得到估計(jì)向量Xest,如下所示

    w1X(t1)+w2X(t2)+…+wmX(tm)

    (9)

    ε=Xest-Xobs

    (10)

    W=(DTD)-1DTXobs

    (11)

    在式(11)中由于Xobs的個(gè)數(shù)遠(yuǎn)大于監(jiān)測(cè)系統(tǒng)參數(shù)個(gè)數(shù),導(dǎo)致DTD不可逆,因此可采用非線性運(yùn)算符代替普通矩陣乘積運(yùn)算,即可得權(quán)重向量式(12)

    W=(DT?D)-1·(DT?Xobs)

    (12)

    式中,?為非線性運(yùn)算符,避免因變量之間的相關(guān)性造成矩陣點(diǎn)乘之后的不可逆現(xiàn)象,從而擴(kuò)大關(guān)系式(12)的通用性。本文MSET重構(gòu)模型采用兩向量之間的歐式距離(Euclidean distance)作為非線性運(yùn)算符

    (13)

    式中,xi和yi均為變量。

    將式(12)代入式(9)可得估計(jì)向量的表達(dá)式

    Xest=D·(DT?D)-1·(DT?Xobs)

    (14)

    式(10)、式(13)和式(14)的物理意義為:當(dāng)新輸入的觀測(cè)向量Xobs仍然處于設(shè)備正常工作狀態(tài)區(qū)間范圍時(shí),其與歷史記憶矩陣D中某些觀測(cè)向量的組合會(huì)很相似,歐氏距離會(huì)很小,導(dǎo)致輸出的估計(jì)向量Xest預(yù)測(cè)值精度高,重構(gòu)精度高,殘差ε小;當(dāng)新輸入的觀測(cè)向量Xobs為設(shè)備出現(xiàn)故障時(shí),勢(shì)必導(dǎo)致其動(dòng)態(tài)特性發(fā)生變化,使Xobs偏離原始正常工作范圍,從而無法利用D中某些觀測(cè)向量的組合進(jìn)行重構(gòu),兩者相似度低,歐氏距離大,導(dǎo)致輸出的估計(jì)向量Xest預(yù)測(cè)值精度低,殘差ε大。由于Xobs和Xest兩向量間的異同性包含豐富的運(yùn)行狀態(tài)信息,因此常利用二者間的殘差進(jìn)行故障模式識(shí)別或性能退化評(píng)估。

    2 整體優(yōu)化的性能退化評(píng)估模型

    由上可知MSET重構(gòu)模型能有效地獲取故障信息并對(duì)性能退化程度進(jìn)行評(píng)估。但應(yīng)進(jìn)一步考慮以下情況:①定量分析模型重構(gòu)前后的殘差從而對(duì)退化趨勢(shì)進(jìn)行量化描述;②確保高維多域特征經(jīng)降維后得到的低維特征與MSET模型歷史記憶矩陣中觀測(cè)向量個(gè)數(shù)是匹配的;③設(shè)置合理的預(yù)警閾值從而對(duì)設(shè)備故障進(jìn)行及時(shí)預(yù)警。這三點(diǎn)都是本文方法的關(guān)鍵所在,后續(xù)將進(jìn)一步進(jìn)行研究和分析。

    2.1 試驗(yàn)數(shù)據(jù)分析

    軸承在實(shí)際運(yùn)行中,其故障的產(chǎn)生和發(fā)展是一個(gè)突發(fā)的、連續(xù)變化的過程,為了真實(shí)地反映故障連續(xù)變化的過程和驗(yàn)證本文所提方法有效性和可靠性,利用西安交大-昇陽科技聯(lián)合實(shí)驗(yàn)室(XJTU-SY)所提供的全壽命加速疲勞軸承試驗(yàn)數(shù)據(jù)[20]進(jìn)行分析。試驗(yàn)臺(tái)如圖1所示,包括驅(qū)動(dòng)部分、測(cè)試部分和液壓加載部分,具體器件如圖中所示。DT9837采集器通過兩個(gè)PCB 352C33單向加速度傳感器對(duì)型號(hào)為L(zhǎng)DK UER204的試驗(yàn)軸承進(jìn)行信號(hào)采集,轉(zhuǎn)軸轉(zhuǎn)頻為40 Hz,采樣頻率為25.6 kHz。試驗(yàn)設(shè)計(jì)了三種工況,每種工況有 5組測(cè)試軸承,工況類別見表1。

    表1 加速疲勞試驗(yàn)工況Tab.1 Accelerated fatigue test condition

    圖1 XJTU-SY試驗(yàn)臺(tái)Fig.1 The XJTU-SY test rig

    本文以工況3第一個(gè)軸承Bearing 3_1為例進(jìn)行分析,該軸承從開始運(yùn)行到出現(xiàn)外圈故障直至失效累計(jì)時(shí)長(zhǎng)42 h 18 min,試驗(yàn)中每隔1 min采樣1.28 s數(shù)據(jù),共計(jì)2 538組數(shù)據(jù),每組數(shù)據(jù)32 768個(gè)點(diǎn),圖2為該軸承全壽命時(shí)域波形圖。利用軸承故障特征頻率計(jì)算公式所得的軸承各元件故障頻率如表2所示。

    表2 故障特征頻率Tab.2 The fault characteristic frequency 單位:Hz

    圖2 軸承3_1的時(shí)域振動(dòng)信號(hào)Fig.2 The temporal vibration signal of bearing 3_1

    由前言可知,不同域的特征包含滾動(dòng)軸承振動(dòng)信號(hào)在不同狀態(tài)空間的信息,例如表3所示的11個(gè)時(shí)域特征中絕對(duì)平均值、標(biāo)準(zhǔn)差、均方根值和峰值4個(gè)指標(biāo)反映時(shí)域信號(hào)的能量大小,其余7個(gè)指標(biāo)反映信號(hào)在時(shí)域的分布特征。而為盡可能地反映軸承信號(hào)在頻域中所包涵的信息,取表4中13個(gè)常見的頻域統(tǒng)計(jì)特征,其中序號(hào)1的特征反映信號(hào)在頻域中振動(dòng)能量大小,序號(hào)2~9的特征反映頻譜的集中程度,序號(hào)10~13的特征反映主頻帶在頻譜中的位置變化。

    表3 11個(gè)時(shí)域特征指標(biāo)Tab.3 11 time domain characteristic indexes

    表4 13個(gè)頻域特征指標(biāo)Tab.4 13 frequency domain characteristic indexes

    對(duì)各樣本提取AR 系數(shù)特征時(shí),AR 模型階數(shù)的確定至關(guān)重要,階數(shù)過低則無法準(zhǔn)確對(duì)模型進(jìn)行描述,反之則會(huì)導(dǎo)致過擬合,本文利用常見的BIC 準(zhǔn)則對(duì)階數(shù)進(jìn)行確定。Bearing 3_1軸承全壽命周期2 538 min,本文認(rèn)為試驗(yàn)剛開始時(shí)的前60 min為設(shè)備開機(jī)后的磨合階段,試驗(yàn)數(shù)據(jù)不準(zhǔn)確,于是取第61~第80共 20組測(cè)試軸承正常運(yùn)行下的壽命信號(hào)分別計(jì)算其AR 最優(yōu)階數(shù),由圖3可知,最優(yōu)AR 模型階數(shù)為26階,所以后續(xù)使用該軸承數(shù)據(jù)做故障診斷提取AR 系數(shù)特征時(shí),一律取26階AR 模型系數(shù)。

    圖3 20個(gè)無故障信號(hào)的最優(yōu)階數(shù)Fig.3 Optimal orders of 20 bearing signals without faults

    由1.1節(jié)可知,信號(hào)數(shù)據(jù)進(jìn)行j層小波包分解后將分成2j個(gè)子帶信號(hào),所以三層小波包分解將產(chǎn)生8個(gè)子帶信號(hào),依次對(duì)每個(gè)子帶信號(hào)計(jì)算其Renyi 熵,即可組成一個(gè)8維三層小波包Renyi 熵特征向量。

    綜上,本文采用上述24組時(shí)域和頻域特征、26階AR模型系數(shù)和8維三層小波包Renyi 熵歸一化處理后共同組成58維多域特征向量指標(biāo)。

    2.2 性能退化評(píng)估模型

    本文基于MSET重構(gòu)模型整體優(yōu)化的軸承性能退化評(píng)估模型如圖4所示,具體步驟如下:

    圖4 性能退化評(píng)估模型Fig.4 Evaluation model of performance degradation

    步驟1基于滾動(dòng)軸承振動(dòng)信號(hào)對(duì)原始時(shí)域數(shù)據(jù)提取24組時(shí)域和頻域特征、26階AR模型系數(shù)和8維三層小波包Renyi 熵歸一化后組成58維多域特征向量指標(biāo);

    步驟2利用滾動(dòng)軸承正常運(yùn)行狀態(tài)下健康數(shù)據(jù)特征構(gòu)建多元狀態(tài)估計(jì)模型的歷史記憶矩陣;

    步驟3利用遺傳算法對(duì)訓(xùn)練集高維特征向量和MSET中歷史記憶矩陣的觀測(cè)向量個(gè)數(shù)和維度進(jìn)行同步聯(lián)合優(yōu)化,實(shí)現(xiàn)特征降維和個(gè)數(shù)優(yōu)選;

    步驟4對(duì)測(cè)試集數(shù)據(jù)提取低維特征,維度與GA優(yōu)化后一致,將得到的測(cè)試樣本低維特征向量作為觀測(cè)向量輸入MSET模型中;

    步驟5輸出估計(jì)向量,計(jì)算觀測(cè)向量與估計(jì)向量之間相似性。

    經(jīng)MSET模型重構(gòu)后的Xest為單重構(gòu)的低維向量,并不能直接反映劣化程度,本文選用重構(gòu)前后兩向量之間的余弦相似度(cosine similarity, CS)作為性能退化指標(biāo),從兩向量空間距離的角度評(píng)估滾動(dòng)軸承性能退化程度,CS表達(dá)式如下

    (15)

    GA相較于粒子群優(yōu)化[21]和模擬退火[22]等優(yōu)化算法而言,對(duì)初始條件的設(shè)置要求較低,只需設(shè)定相應(yīng)的適應(yīng)度函數(shù)即可,具有較強(qiáng)的魯棒性,能更快地搜索到全局最優(yōu)解。上述步驟3中利用GA對(duì)模型進(jìn)行同步優(yōu)化的具體操作如下:

    步驟1采用二進(jìn)制編碼方式進(jìn)行編碼,對(duì)測(cè)試軸承前n個(gè)樣本和58維特征向量組成染色體長(zhǎng)度為n+58的基因組,其排列方式如圖5所示;

    圖5 編碼排列方式Fig.5 Coding arrangement mode

    步驟2初始化GA參數(shù),設(shè)置初始種群大小為50,遺傳代數(shù)為100,采用多點(diǎn)交叉方式,交叉率為0.7,變異率為0.05,代溝值為0.9;

    步驟3設(shè)置訓(xùn)練樣本集合,主要由無故障階段和失效階段兩部分信號(hào)構(gòu)成;并采用所有正常樣本CS均值與所有故障樣本CS均值的比值設(shè)計(jì)適應(yīng)度函數(shù),從而計(jì)算個(gè)體適應(yīng)度,適應(yīng)度越高,個(gè)體越優(yōu);

    (16)

    式中:m為正常樣本數(shù)量;n為故障樣本數(shù)量;CSN為正常樣本的余弦相似度;CSF為故障樣本的余弦相似度。

    步驟4根據(jù)步驟2中所確定的交叉和變異遺傳算子對(duì)當(dāng)前種群進(jìn)行處理,并判別其種群個(gè)體適應(yīng)度是否達(dá)到要求。若達(dá)到要求則進(jìn)行解碼,得到最優(yōu)低維向量的維數(shù)及其索引以及最優(yōu)歷史記憶矩陣觀測(cè)向量個(gè)數(shù)和維度,執(zhí)行步驟6,否則繼續(xù)下一步;

    步驟5進(jìn)行選擇、交叉、變異過程,返回步驟3;

    步驟6獲得最優(yōu)參數(shù)建立MSET單重構(gòu)模型,利用模型進(jìn)行性能退化評(píng)估。

    2.3 自適應(yīng)預(yù)警閾值

    余弦相似度CS是一條連續(xù)變化的曲線,將其作為劣化指標(biāo)能準(zhǔn)確地量化表征軸承目前運(yùn)行狀態(tài)與無故障運(yùn)行時(shí)的偏離程度。進(jìn)一步還需要對(duì)該曲線標(biāo)定一個(gè)恰當(dāng)?shù)念A(yù)警值從而及時(shí)對(duì)故障進(jìn)行示警。

    由統(tǒng)計(jì)學(xué)知識(shí)可知,數(shù)學(xué)期望為μ,方差為σ2的正態(tài)分布(μ,σ2)數(shù)據(jù),樣本數(shù)值落在(μ-3σ,μ+3σ)區(qū)間內(nèi)的概率約為0.997 3,因此判定在樣本故障程度相近時(shí),CS值也符合正態(tài)分布要求。進(jìn)一步當(dāng)連續(xù)多個(gè)樣本的CS值超過由歷史樣本CS值所確定的(μ-3σ,μ+3σ)區(qū)間范圍時(shí),判定該時(shí)刻故障程度于前一刻已經(jīng)發(fā)生較大變化。CS值是連續(xù)單調(diào)遞減的,故只計(jì)算下限μ-3σ即可,如式(17)所示。本文利用3σ準(zhǔn)則對(duì)CS值進(jìn)行預(yù)警值標(biāo)定,未超過該值的視為正常運(yùn)行狀態(tài),同時(shí)利用當(dāng)前所有無故障樣本更新CS值的3σ區(qū)間,達(dá)到自適應(yīng)示警的目的。

    式中:{D(i),i=1,2,…,n}為t=[1,n]時(shí)刻的CS值;T(t)為t時(shí)刻預(yù)警值。當(dāng)首次連續(xù)多個(gè)CS值小于te時(shí)刻的預(yù)警值T(te)時(shí),定義該時(shí)刻為早期故障時(shí)刻或故障程度加深時(shí)刻。為防止外界偶然因素的影響,本文早期故障時(shí)刻要求連續(xù)5個(gè)測(cè)試數(shù)據(jù)的CS值小于或等于te時(shí)刻的預(yù)警閾值。

    3 試驗(yàn)結(jié)果分析

    3.1 聯(lián)合優(yōu)化特征降維

    對(duì)試驗(yàn)數(shù)據(jù)按照步驟1~步驟5進(jìn)行整體優(yōu)化性能退化評(píng)估操作,具體如下:首先從健康運(yùn)行軸承的采集信號(hào)中選擇300個(gè)樣本提取其58維特征作為觀測(cè)向量用來構(gòu)造歷史記憶矩陣D。然后額外選取100個(gè)正常運(yùn)行狀態(tài)的樣本和100個(gè)故障樣本提取其高維特征作為訓(xùn)練數(shù)據(jù),利用遺傳算法對(duì)訓(xùn)練集高維特征向量和MSET 模型中歷史記憶矩陣的觀測(cè)向量個(gè)數(shù)和維度進(jìn)行同步聯(lián)合優(yōu)化。最后將Bearing 3_1 軸承全壽命周期2 538個(gè)樣本作為測(cè)試數(shù)據(jù),提取其低維特征,維度與GA優(yōu)化后一致,將得到的測(cè)試樣本低維特征向量作為新輸入觀測(cè)向量輸入MSET模型中進(jìn)行后續(xù)分析。

    對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行樣本劃分后按上述過程進(jìn)行整體優(yōu)化性能退化評(píng)估,優(yōu)化結(jié)果如圖6所示。從圖6(a)中可以看出,經(jīng)優(yōu)化后歷史記憶矩陣所包含的觀測(cè)向量個(gè)數(shù)由原始的300個(gè)變?yōu)?33個(gè),且分布較為均勻,有效地剔除了部分冗余觀測(cè)向量,減少了后續(xù)性能退化評(píng)估的計(jì)算量;從圖6(b)中可以看出,在步驟3同步聯(lián)合優(yōu)化中,58維特征優(yōu)選成了9維特征,由小波Renyi 熵后三項(xiàng)和AR模型中6個(gè)系數(shù)組成,其中24組時(shí)域和頻域特征降維成了0維,表明時(shí)、頻域指標(biāo)對(duì)西安交大軸承數(shù)據(jù)中的Bearing 3_1試驗(yàn)數(shù)據(jù)分析結(jié)果較差。

    圖6 優(yōu)化結(jié)果Fig.6 Optimization results

    為進(jìn)一步分析該試驗(yàn)數(shù)據(jù)特征優(yōu)選后時(shí)域、頻域特征降維后變成了0維的原因,對(duì)3個(gè)時(shí)域指標(biāo)和3個(gè)頻域指標(biāo)進(jìn)行作圖,如圖7所示。從圖7中可以看出,脈沖因子和峭度值與軸承故障發(fā)展趨勢(shì)基本一致,但后期均出現(xiàn)先減小后增大的情況,與故障實(shí)際發(fā)展情況不相符,并且脈沖因子指標(biāo)前期波動(dòng)過大,不利于早期故障的發(fā)現(xiàn);加速度有效值波形后期振蕩過于嚴(yán)重,與故障實(shí)際發(fā)展情況不相符;其余3個(gè)頻域特征均不能準(zhǔn)確表征故障程度和性能退化趨勢(shì)。

    圖7 軸承時(shí)域指標(biāo)Fig.7 Time domain index of bearing

    3.2 全壽命曲線

    對(duì)Bearing 3_1軸承全壽命2 538個(gè)樣本進(jìn)行分析,性能退化評(píng)估CS 值曲線經(jīng)過5點(diǎn)平滑處理后如圖8所示。由于故障程度具有漸變性,平滑處理不僅不會(huì)影響故障預(yù)警,更能削弱噪聲干擾,更好地表達(dá)信號(hào)。圖中虛線為自適應(yīng)預(yù)警值,實(shí)線為故障指標(biāo)CS 值,曲線整體呈現(xiàn)下降趨勢(shì),表明CS 值能準(zhǔn)確、一致地表征軸承故障嚴(yán)重度趨勢(shì),且相較于圖7的時(shí)域指標(biāo)對(duì)軸承早期故障更加敏感。

    圖8 軸承全壽命CS值Fig.8 CS value of bearing life

    圖8中性能退化評(píng)估曲線大致可劃分成無故障時(shí)期、初始故障時(shí)期、中度故障時(shí)期和惡化失效時(shí)期4個(gè)階段。由于該軸承的CS 值在534時(shí)刻起連續(xù)5個(gè)值小于自適應(yīng)預(yù)警值,故0~533時(shí)刻為無故障時(shí)期,534時(shí)刻為初始故障發(fā)生時(shí)刻。534~1 119時(shí)刻內(nèi)的CS 值呈明顯下降趨勢(shì),但在1 119時(shí)刻曲線出現(xiàn)階躍現(xiàn)象,CS 值急劇下降,這是因?yàn)樵摃r(shí)刻軸承故障程度持續(xù)加重,導(dǎo)致新輸入觀測(cè)向量與歷史記憶矩陣的偏離程度持續(xù)加深,最終在1 119時(shí)刻達(dá)到中度故障時(shí)期。1 119~1 378時(shí)刻內(nèi)的CS 值振蕩降低后垂直上升,這是因?yàn)楣收系陌l(fā)展是一個(gè)漸變的過程,從微小故障逐漸擴(kuò)大,然后磨平,隨著時(shí)間增加故障再度擴(kuò)大,所以導(dǎo)致CS 值出現(xiàn)了大范圍地反復(fù)波動(dòng)。1 378時(shí)刻,原始故障點(diǎn)已被磨平,因故障產(chǎn)生的突變沖擊減小,新輸入觀測(cè)向量與歷史記憶矩陣的偏離程度減小,CS 值突變?cè)龃蟆? 378~2 538時(shí)刻,軸承故障程度不斷加重,且CS 值波動(dòng)程度比之前各階段更大,波動(dòng)時(shí)間更短,表明故障點(diǎn)的“磨平- 惡化”情況更加明顯,且故障程度快速加深直至軸承失效。從圖中還可以看出,CS 值曲線整體毛刺較多,局部波動(dòng)性大,這是由于軸承在整體運(yùn)轉(zhuǎn)過程中微小的“故障磨平、再度劣化”現(xiàn)象一直存在。

    為進(jìn)一步分析CS 值與故障程度之間的關(guān)聯(lián)性,選取無故障階段(0~533)、 輕度故障階段(534~1 118)、 中度故障階段(1 119~1 377)、重度故障階段(1 378~2 538) 四種狀態(tài)中各3個(gè)樣本,計(jì)算其低維特征后輸入聯(lián)合優(yōu)化模型,可得到輸出估計(jì)向量與觀測(cè)向量之間的余弦相似度差異,不同故障程度下CS 值的變化如圖9所示。從圖中可以看出,無故障樣本的CS 值最大,3個(gè)樣本的均值約為0.985,接近于1,這是因?yàn)闊o故障樣本提取的特征與歷史記憶矩陣中各向量相似程度高,輸出估計(jì)向量與觀測(cè)向量的偏斜度小,兩空間向量之間的夾角小,余弦相似度大。而隨著外圈故障尺寸的變大,估計(jì)向量與觀測(cè)向量之間的CS值越來越小,這是由于隨著故障嚴(yán)重度加深,軸承當(dāng)前狀態(tài)越來越偏離無故障運(yùn)行狀態(tài),輸入特征越來越偏離依靠健康樣本所建立的歷史記憶矩陣的觀測(cè)向量,兩者之間的偏斜度越來越大,在CS 值上的表現(xiàn)則呈現(xiàn)出遞減趨勢(shì),從而很好地區(qū)分軸承的故障程度。

    圖9 不同程度外圈故障CS值Fig.9 CS value of outer ring faults in different degrees

    4 試驗(yàn)驗(yàn)證及對(duì)比

    4.1 試驗(yàn)驗(yàn)證

    為驗(yàn)證聯(lián)合優(yōu)化多域特征-MSET 重構(gòu)模型方法的有效性,對(duì)No.534和No.533樣本信號(hào)進(jìn)行以包絡(luò)譜譜峰因子為優(yōu)化指標(biāo)自適應(yīng)共振解調(diào)[23]處理,共振解調(diào)后對(duì)軸承信號(hào)進(jìn)行傅里葉變換,如圖10所示。

    圖10 故障樣本包絡(luò)解調(diào)Fig.10 Envelope demodulation of fault samples

    圖10(a) No.534樣本共振解調(diào)結(jié)果中出現(xiàn)明顯脈沖現(xiàn)象,且在頻譜的包絡(luò)譜中出現(xiàn)了與外圈故障頻率BPFO(123.3 Hz)相近的一倍頻124 Hz、二倍頻246 Hz和三倍頻370 Hz,此外轉(zhuǎn)軸轉(zhuǎn)頻40 Hz也清晰可見。但由于該時(shí)刻為早期故障點(diǎn),軸承外圈剛剛出現(xiàn)微小故障,導(dǎo)致二倍頻和三倍頻幾乎泯沒在噪聲當(dāng)中。而在圖10(b)No.533樣本包絡(luò)解調(diào)中沖激現(xiàn)象和123.3 Hz的外圈故障頻率成分并不明顯,且在534時(shí)刻的最優(yōu)濾波頻帶相較于533時(shí)刻的最優(yōu)濾波頻帶發(fā)生了明顯改變。因此,證明534時(shí)刻為軸承早期故障點(diǎn),驗(yàn)證了方法的有效性和可靠性。

    4.2 對(duì)比試驗(yàn)分析

    為驗(yàn)證本文所提方法的優(yōu)越性,首先采用僅進(jìn)行特征維度優(yōu)化的MSET軸承性能退化模型進(jìn)行對(duì)比分析。利用GA對(duì)58維特征向量進(jìn)行降維,如圖11所示,從圖中可以看出僅對(duì)特征進(jìn)行優(yōu)化時(shí),時(shí)域、頻域的特征向量仍有所保留,與圖6的分析不相符,導(dǎo)致輸入將含有許多冗余信息。進(jìn)一步分析,將降維后的34維特征向量輸入未經(jīng)優(yōu)化的MSET模型,其性能退化曲線和自適應(yīng)預(yù)警值如圖12(a)所示,其早期故障點(diǎn)為No.306樣本,遠(yuǎn)遠(yuǎn)早于本文方法所判斷的早期故障點(diǎn),原因是因?yàn)槲唇?jīng)優(yōu)化的MSET模型歷史記憶中的觀測(cè)向量個(gè)數(shù)過多,而輸入的特征向量維度也過多,導(dǎo)致含有許多冗余信息從而提前誤判故障點(diǎn)。為驗(yàn)證該時(shí)刻軸承并未發(fā)生故障,對(duì)No.306樣本信號(hào)進(jìn)行包絡(luò)解調(diào),如圖12(b)所示,其頻譜中并無外圈故障頻率成分,證明此刻軸承仍處于正常運(yùn)行狀態(tài),并未發(fā)生故障,早期故障點(diǎn)判斷失誤。

    圖11 多域特征向量?jī)?yōu)化結(jié)果Fig.11 Multi-domain eigenvector optimization result

    圖12 對(duì)比試驗(yàn)一分析Fig.12 Analysis of the comparative experiment one

    其次,采用僅優(yōu)化MSET中歷史記憶矩陣的軸承性能退化模型進(jìn)行對(duì)比分析。利用GA對(duì)觀測(cè)向量的個(gè)數(shù)進(jìn)行優(yōu)選后,如圖13所示。從圖13中可以看出,歷史記憶矩陣觀測(cè)向量個(gè)數(shù)由原來的300個(gè)變成了156個(gè),分布較為均勻,有效剔除了冗余觀測(cè)向量。進(jìn)一步分析,將提取的高維特征向量直接輸入優(yōu)化后的MSET模型,其性能退化曲線和自適應(yīng)預(yù)警值如圖14所示。從圖14中可以看出,在No.320樣本處CS值曲線與自適應(yīng)預(yù)警值曲線已有切點(diǎn),預(yù)示軸承已經(jīng)開始出現(xiàn)故障,但從圖14(b)的包絡(luò)解調(diào)結(jié)果可以看到該樣本僅包含轉(zhuǎn)頻40 Hz,因此該預(yù)測(cè)結(jié)果不準(zhǔn)確。而由于未對(duì)輸入高維特征進(jìn)行優(yōu)化,導(dǎo)致性能退化曲線并不能很好地區(qū)分軸承的故障程度,自適應(yīng)預(yù)警值也不能準(zhǔn)確有效地反映出早期故障發(fā)生時(shí)刻,致使長(zhǎng)時(shí)間后在No.2296樣本處兩曲線再次出現(xiàn)交點(diǎn),而該時(shí)刻雖已發(fā)生故障,但遠(yuǎn)遠(yuǎn)晚于本文方法所判斷的早期故障點(diǎn),時(shí)效性不強(qiáng),不具有使用意義。

    圖13 觀測(cè)向量?jī)?yōu)化結(jié)果Fig.13 Optimization result of observation vector

    圖14 對(duì)比試驗(yàn)二分析Fig.14 Analysis of the comparative experiment two

    由以上兩個(gè)對(duì)比試驗(yàn)綜合分析可知,相比于本文所提出的聯(lián)合優(yōu)化多域特征-MSET模型,單優(yōu)化輸入特征維度和單優(yōu)化歷史記憶矩陣的模型效果均一般,與軸承的實(shí)際發(fā)展趨勢(shì)不一致,不能及時(shí)地、準(zhǔn)確地反映軸承早期故障,也無法很好地區(qū)分故障程度,由此說明本文所提的整體優(yōu)化MSET單重構(gòu)模型方法的準(zhǔn)確有效性和優(yōu)越性。

    進(jìn)一步為驗(yàn)證所提模型的泛化性和優(yōu)越性,采用多個(gè)軸承全壽命數(shù)據(jù)進(jìn)行性能退化評(píng)估,并與其他文獻(xiàn)中結(jié)果進(jìn)行對(duì)比。從表5中可以看到,所用方法在不同軸承數(shù)據(jù)集下均能更準(zhǔn)確地判別軸承性能退化趨勢(shì),模型具有一定的泛化性。

    表5 不同軸承數(shù)據(jù)間的試驗(yàn)結(jié)果Tab.5 The results between different bearing data

    5 結(jié) 論

    針對(duì)傳統(tǒng)單域特征指標(biāo)無法充分表征軸承狀態(tài)信息,而目前性能退化評(píng)估中基于多域高維特征向量的重構(gòu)評(píng)估模型存在信息冗余且易受到不一致優(yōu)化目標(biāo)的影響,不能在低維空間中保留必要的信息,而導(dǎo)致模型次優(yōu)性能。提出一種基于MSET重構(gòu)模型整體優(yōu)化的軸承性能退化評(píng)估方法,通過同步聯(lián)合優(yōu)化多域特征和MSET 模型歷史記憶矩陣從而實(shí)現(xiàn)特征優(yōu)選和性能評(píng)估模型的整體自適應(yīng)優(yōu)化。

    (1)XJTU-SY 軸承全壽命數(shù)據(jù)分析表明經(jīng)聯(lián)合優(yōu)化多域特征-MSET 重構(gòu)模型后的所得CS 值能很好地表達(dá)軸承故障嚴(yán)重度變化,且隨著程度加深具有單調(diào)遞減趨勢(shì);并通過對(duì)傳統(tǒng)時(shí)域指標(biāo)進(jìn)行比較,表明本文所提方法對(duì)軸承早期故障更加敏感,且與故障發(fā)展趨勢(shì)的一致性更好。

    (2)通過對(duì)兩組單優(yōu)化對(duì)比試驗(yàn)和不同方法進(jìn)行比較分析,并利用復(fù)小波共振解調(diào)法對(duì)早期故障點(diǎn)進(jìn)行驗(yàn)證,表明本文基于MSET重構(gòu)模型整體優(yōu)化的軸承性能退化評(píng)估方法更具有效性、可靠性和優(yōu)越性。

    猜你喜歡
    軸承向量觀測(cè)
    觀測(cè)到恒星死亡瞬間
    軍事文摘(2023年18期)2023-11-03 09:45:42
    軸承知識(shí)
    哈爾濱軸承(2022年2期)2022-07-22 06:39:32
    軸承知識(shí)
    哈爾濱軸承(2022年1期)2022-05-23 13:13:24
    向量的分解
    軸承知識(shí)
    哈爾濱軸承(2021年2期)2021-08-12 06:11:46
    軸承知識(shí)
    哈爾濱軸承(2021年1期)2021-07-21 05:43:16
    聚焦“向量與三角”創(chuàng)新題
    天測(cè)與測(cè)地VLBI 測(cè)地站周圍地形觀測(cè)遮掩的討論
    可觀測(cè)宇宙
    太空探索(2016年7期)2016-07-10 12:10:15
    向量垂直在解析幾何中的應(yīng)用
    久久午夜亚洲精品久久| 特大巨黑吊av在线直播| 成人综合一区亚洲| 亚洲在线观看片| 狂野欧美白嫩少妇大欣赏| 成熟少妇高潮喷水视频| 国产免费男女视频| 日韩一本色道免费dvd| www日本黄色视频网| 国产一区二区激情短视频| 男女视频在线观看网站免费| 国产日本99.免费观看| 久久精品国产鲁丝片午夜精品| 看片在线看免费视频| .国产精品久久| 自拍偷自拍亚洲精品老妇| 亚州av有码| 久久久精品欧美日韩精品| 亚洲四区av| 久久精品国产99精品国产亚洲性色| 可以在线观看毛片的网站| 在线观看66精品国产| 最好的美女福利视频网| 亚洲最大成人av| 在现免费观看毛片| 精品久久久久久久久久免费视频| 美女脱内裤让男人舔精品视频 | 色哟哟·www| 性插视频无遮挡在线免费观看| 桃色一区二区三区在线观看| 波多野结衣巨乳人妻| 国产精品.久久久| 精品久久久久久久久久免费视频| 亚洲成人精品中文字幕电影| 黄色一级大片看看| 国产乱人偷精品视频| 欧美3d第一页| 国内久久婷婷六月综合欲色啪| 免费观看精品视频网站| a级毛片免费高清观看在线播放| 高清毛片免费观看视频网站| 欧美极品一区二区三区四区| 伦理电影大哥的女人| 亚洲精品自拍成人| 色综合站精品国产| 久久久久性生活片| 国产午夜精品一二区理论片| 久久精品人妻少妇| 激情 狠狠 欧美| 亚洲国产精品成人久久小说 | 精品国产三级普通话版| 一边摸一边抽搐一进一小说| 国产精品麻豆人妻色哟哟久久 | 在线国产一区二区在线| 22中文网久久字幕| 亚洲欧洲日产国产| 久久人人爽人人片av| 日韩一区二区视频免费看| 欧美区成人在线视频| 久久这里只有精品中国| 日韩欧美精品免费久久| 亚洲人成网站高清观看| 成人毛片60女人毛片免费| 免费电影在线观看免费观看| 男人的好看免费观看在线视频| 亚洲精品乱码久久久v下载方式| 久久久久国产网址| 久久这里只有精品中国| 午夜福利在线观看吧| 九九爱精品视频在线观看| 国产精品国产高清国产av| 哪个播放器可以免费观看大片| 18禁在线播放成人免费| 欧美成人免费av一区二区三区| 欧美日韩国产亚洲二区| 少妇丰满av| 亚洲av一区综合| 国产精品一二三区在线看| 亚洲精华国产精华液的使用体验 | 91精品国产九色| 国产伦一二天堂av在线观看| 男人舔女人下体高潮全视频| 一级毛片我不卡| 精品久久久久久久久久久久久| 黄色配什么色好看| 丰满的人妻完整版| 午夜精品在线福利| 国产色爽女视频免费观看| 中文亚洲av片在线观看爽| 日日撸夜夜添| 国产淫片久久久久久久久| 欧美成人免费av一区二区三区| 午夜久久久久精精品| 久久热精品热| 国产爱豆传媒在线观看| 日韩人妻高清精品专区| 免费观看a级毛片全部| 在线免费观看的www视频| 偷拍熟女少妇极品色| 性色avwww在线观看| 日本黄色片子视频| 色综合色国产| 国产精品人妻久久久影院| 国产老妇伦熟女老妇高清| 日韩视频在线欧美| 日韩精品青青久久久久久| 国产一区二区三区在线臀色熟女| 一卡2卡三卡四卡精品乱码亚洲| 国内精品一区二区在线观看| 日日撸夜夜添| 狂野欧美白嫩少妇大欣赏| 亚洲精华国产精华液的使用体验 | 简卡轻食公司| 久久精品国产亚洲av涩爱 | 亚洲av中文字字幕乱码综合| 变态另类成人亚洲欧美熟女| 日韩欧美一区二区三区在线观看| 亚洲精品国产av成人精品| 国产色爽女视频免费观看| 精品一区二区三区人妻视频| 成人漫画全彩无遮挡| 成人二区视频| 欧美最新免费一区二区三区| 丰满人妻一区二区三区视频av| videossex国产| 男女那种视频在线观看| 国产成人精品久久久久久| 国产av不卡久久| 欧美xxxx性猛交bbbb| 国产精品人妻久久久久久| 国产三级在线视频| 午夜福利在线观看吧| 亚洲人成网站高清观看| 欧美日本亚洲视频在线播放| 亚洲精品影视一区二区三区av| 一本精品99久久精品77| 国产高潮美女av| 欧美精品国产亚洲| 久久99热6这里只有精品| 免费看光身美女| 18禁黄网站禁片免费观看直播| 五月玫瑰六月丁香| 中国国产av一级| 一本一本综合久久| 嘟嘟电影网在线观看| 国产亚洲av片在线观看秒播厂 | 国产不卡一卡二| 一级二级三级毛片免费看| 精品欧美国产一区二区三| 亚洲美女搞黄在线观看| 色综合色国产| 97超视频在线观看视频| 久久99热6这里只有精品| 99久久九九国产精品国产免费| 日韩欧美一区二区三区在线观看| 男人和女人高潮做爰伦理| 可以在线观看毛片的网站| 2021天堂中文幕一二区在线观| 精品免费久久久久久久清纯| 成人美女网站在线观看视频| 亚洲av中文av极速乱| 观看美女的网站| 九九爱精品视频在线观看| 免费观看在线日韩| av在线天堂中文字幕| 欧美成人免费av一区二区三区| 九九在线视频观看精品| 中国国产av一级| 99热这里只有精品一区| 国产一级毛片七仙女欲春2| 成人特级av手机在线观看| 国产单亲对白刺激| 三级国产精品欧美在线观看| 观看美女的网站| 国内久久婷婷六月综合欲色啪| 永久网站在线| 亚洲七黄色美女视频| 午夜久久久久精精品| 日韩一本色道免费dvd| 最近中文字幕高清免费大全6| 亚洲第一电影网av| 2022亚洲国产成人精品| 春色校园在线视频观看| 亚洲成av人片在线播放无| 黄片wwwwww| 美女cb高潮喷水在线观看| 国产激情偷乱视频一区二区| 草草在线视频免费看| 国产淫片久久久久久久久| 成人永久免费在线观看视频| 精品一区二区免费观看| 非洲黑人性xxxx精品又粗又长| 精品少妇黑人巨大在线播放 | 97热精品久久久久久| 免费人成在线观看视频色| 国产午夜精品久久久久久一区二区三区| 日本免费a在线| 又爽又黄a免费视频| ponron亚洲| 91在线精品国自产拍蜜月| 男女那种视频在线观看| 少妇被粗大猛烈的视频| 最近2019中文字幕mv第一页| av福利片在线观看| 国产成人福利小说| 久久草成人影院| 国产亚洲精品久久久久久毛片| 亚洲国产欧美人成| 亚洲天堂国产精品一区在线| 成年女人看的毛片在线观看| 久久婷婷人人爽人人干人人爱| 亚洲精品国产成人久久av| 天天躁日日操中文字幕| 晚上一个人看的免费电影| 欧美变态另类bdsm刘玥| 边亲边吃奶的免费视频| 久久久久久久久久黄片| 欧美xxxx黑人xx丫x性爽| 老司机福利观看| 波野结衣二区三区在线| 欧美性猛交黑人性爽| 国产91av在线免费观看| 国产乱人视频| 久久久久九九精品影院| 少妇熟女aⅴ在线视频| 91麻豆精品激情在线观看国产| 亚洲精品国产av成人精品| 亚洲图色成人| kizo精华| 人妻夜夜爽99麻豆av| 99热只有精品国产| 男插女下体视频免费在线播放| 久久久久久久久久久丰满| 久久久久久久久久久免费av| 亚洲欧洲国产日韩| 两个人视频免费观看高清| 精品99又大又爽又粗少妇毛片| 九草在线视频观看| 国产精品精品国产色婷婷| 级片在线观看| 18禁黄网站禁片免费观看直播| 在现免费观看毛片| 亚洲成人av在线免费| 国产视频首页在线观看| 高清在线视频一区二区三区 | 国产亚洲5aaaaa淫片| 97热精品久久久久久| 在线播放无遮挡| av免费在线看不卡| 亚洲av电影不卡..在线观看| 亚洲真实伦在线观看| 日韩高清综合在线| 免费看光身美女| 一边亲一边摸免费视频| 亚洲精品自拍成人| 欧美在线一区亚洲| 国产伦理片在线播放av一区 | 一卡2卡三卡四卡精品乱码亚洲| 一区二区三区高清视频在线| 日韩高清综合在线| 久久久久久久久久久免费av| 日韩成人av中文字幕在线观看| 午夜精品一区二区三区免费看| 免费观看a级毛片全部| 亚洲国产精品成人久久小说 | 97热精品久久久久久| 久久久久免费精品人妻一区二区| 精品久久久久久久久av| 亚洲人成网站高清观看| 色综合亚洲欧美另类图片| 少妇熟女欧美另类| 日本黄色视频三级网站网址| 午夜福利在线观看吧| 色5月婷婷丁香| АⅤ资源中文在线天堂| 日本一本二区三区精品| 九九久久精品国产亚洲av麻豆| 日韩一区二区三区影片| 人体艺术视频欧美日本| 99在线人妻在线中文字幕| 高清毛片免费观看视频网站| 天堂影院成人在线观看| 伦理电影大哥的女人| 婷婷六月久久综合丁香| 精品久久久久久久末码| 国产高清三级在线| 成人毛片60女人毛片免费| 久久久精品94久久精品| 欧美色欧美亚洲另类二区| 伦精品一区二区三区| av国产免费在线观看| 久久午夜福利片| 国产伦一二天堂av在线观看| 国产精品爽爽va在线观看网站| 国产大屁股一区二区在线视频| 精品人妻视频免费看| 亚洲成人久久爱视频| 欧美日韩综合久久久久久| 免费大片18禁| 看黄色毛片网站| 亚洲一级一片aⅴ在线观看| 女的被弄到高潮叫床怎么办| 亚洲人成网站在线播放欧美日韩| 自拍偷自拍亚洲精品老妇| 少妇熟女aⅴ在线视频| 午夜福利高清视频| 又粗又硬又长又爽又黄的视频 | 99久久精品国产国产毛片| 乱码一卡2卡4卡精品| 亚洲欧美精品自产自拍| 国产国拍精品亚洲av在线观看| 国产成人aa在线观看| 内射极品少妇av片p| 国产成人91sexporn| 可以在线观看的亚洲视频| 麻豆国产97在线/欧美| 一本久久中文字幕| 天堂网av新在线| 亚洲精品日韩在线中文字幕 | 欧洲精品卡2卡3卡4卡5卡区| 三级国产精品欧美在线观看| 亚洲国产欧美人成| 日本黄色片子视频| 亚洲色图av天堂| 波多野结衣高清作品| 少妇人妻精品综合一区二区 | 自拍偷自拍亚洲精品老妇| av视频在线观看入口| 中文资源天堂在线| 最近视频中文字幕2019在线8| 亚洲在线观看片| 99热网站在线观看| 久久久久久久久久久丰满| 国产不卡一卡二| 又爽又黄无遮挡网站| 午夜福利在线在线| 人妻系列 视频| 99久久久亚洲精品蜜臀av| 午夜免费激情av| .国产精品久久| 69av精品久久久久久| 国产午夜精品论理片| 久久久精品94久久精品| 麻豆一二三区av精品| 又粗又硬又长又爽又黄的视频 | 亚洲精品成人久久久久久| 99热这里只有是精品在线观看| 身体一侧抽搐| 1024手机看黄色片| 国产成人精品久久久久久| 免费观看人在逋| 亚洲人成网站在线播| 国产黄片美女视频| 精品一区二区免费观看| 日本一二三区视频观看| 成人一区二区视频在线观看| 国产精品无大码| 99久国产av精品国产电影| 九九在线视频观看精品| 亚洲精品国产av成人精品| 三级男女做爰猛烈吃奶摸视频| av在线播放精品| 内射极品少妇av片p| 国产精品女同一区二区软件| 久久久久久久久久成人| 精品免费久久久久久久清纯| 欧美激情久久久久久爽电影| 99久久无色码亚洲精品果冻| 人人妻人人澡人人爽人人夜夜 | 国产毛片a区久久久久| 免费不卡的大黄色大毛片视频在线观看 | 卡戴珊不雅视频在线播放| 黄色欧美视频在线观看| 国产真实伦视频高清在线观看| 一本久久精品| 国产精品不卡视频一区二区| www.色视频.com| 高清毛片免费看| 久久亚洲精品不卡| 国产精品永久免费网站| 国产精品一及| 欧美一级a爱片免费观看看| 久久久久久久久久黄片| 国产探花在线观看一区二区| 久久久久久久午夜电影| 亚洲成av人片在线播放无| 国产精品爽爽va在线观看网站| 午夜免费男女啪啪视频观看| 99九九线精品视频在线观看视频| 亚洲在久久综合| 青青草视频在线视频观看| 久久久久久久久中文| 中文字幕熟女人妻在线| 国产午夜精品论理片| 黄色视频,在线免费观看| 日韩成人伦理影院| 天美传媒精品一区二区| 老熟妇乱子伦视频在线观看| 亚洲精品久久国产高清桃花| 久久精品91蜜桃| 美女cb高潮喷水在线观看| 三级男女做爰猛烈吃奶摸视频| 97热精品久久久久久| 熟女电影av网| 国产亚洲精品久久久com| 两性午夜刺激爽爽歪歪视频在线观看| 国国产精品蜜臀av免费| 亚洲成av人片在线播放无| av视频在线观看入口| 成人高潮视频无遮挡免费网站| 日本av手机在线免费观看| 亚洲最大成人中文| 亚洲婷婷狠狠爱综合网| 美女脱内裤让男人舔精品视频 | 国产精品乱码一区二三区的特点| 午夜亚洲福利在线播放| 网址你懂的国产日韩在线| 亚洲五月天丁香| 九九在线视频观看精品| 久久久国产成人免费| 久久热精品热| 啦啦啦啦在线视频资源| 看非洲黑人一级黄片| 亚洲久久久久久中文字幕| 日韩欧美一区二区三区在线观看| 国产在线精品亚洲第一网站| 九色成人免费人妻av| 小蜜桃在线观看免费完整版高清| 亚洲av一区综合| 久久久久国产网址| 最后的刺客免费高清国语| 你懂的网址亚洲精品在线观看 | 啦啦啦观看免费观看视频高清| 精品久久国产蜜桃| 91狼人影院| 国产精品三级大全| 国产伦一二天堂av在线观看| 少妇猛男粗大的猛烈进出视频 | 亚洲国产色片| 亚洲精品久久久久久婷婷小说 | 国产乱人偷精品视频| 国产不卡一卡二| 国产av一区在线观看免费| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 99久久人妻综合| 国产亚洲精品久久久久久毛片| 亚洲,欧美,日韩| 久久久成人免费电影| 可以在线观看的亚洲视频| 国产免费一级a男人的天堂| 亚洲va在线va天堂va国产| 婷婷六月久久综合丁香| 成人三级黄色视频| 国产真实乱freesex| 国产女主播在线喷水免费视频网站 | 国产成人精品久久久久久| 亚洲成人久久爱视频| 免费看a级黄色片| 成熟少妇高潮喷水视频| 91麻豆精品激情在线观看国产| 午夜精品一区二区三区免费看| 亚洲国产欧美人成| 欧美+日韩+精品| 99国产精品一区二区蜜桃av| 欧美日韩综合久久久久久| 好男人在线观看高清免费视频| 国产三级中文精品| 又爽又黄无遮挡网站| 国产高清激情床上av| 国产亚洲av片在线观看秒播厂 | 悠悠久久av| 日韩一区二区三区影片| 18禁黄网站禁片免费观看直播| 国产乱人偷精品视频| 天堂√8在线中文| 欧美日本视频| 免费人成视频x8x8入口观看| 国产一级毛片七仙女欲春2| 日韩 亚洲 欧美在线| 2021天堂中文幕一二区在线观| 国产91av在线免费观看| 国产成年人精品一区二区| 午夜视频国产福利| 国产黄片视频在线免费观看| 亚洲成人av在线免费| 成人性生交大片免费视频hd| 精品人妻视频免费看| 久久精品国产鲁丝片午夜精品| 又粗又硬又长又爽又黄的视频 | 亚洲av成人av| 国产伦在线观看视频一区| 99久久无色码亚洲精品果冻| 国产精品一区二区性色av| 99热6这里只有精品| 91av网一区二区| 一本精品99久久精品77| 成人特级av手机在线观看| 亚洲最大成人中文| 免费观看的影片在线观看| 国产乱人视频| 亚洲av电影不卡..在线观看| 亚洲国产欧洲综合997久久,| 成人永久免费在线观看视频| 一级毛片我不卡| 午夜福利高清视频| 老司机福利观看| 亚洲自偷自拍三级| 国产伦理片在线播放av一区 | 波多野结衣巨乳人妻| 久久久精品94久久精品| 校园春色视频在线观看| 国产91av在线免费观看| 欧美又色又爽又黄视频| 国产精品一区www在线观看| 欧美人与善性xxx| 我要看日韩黄色一级片| 91狼人影院| av视频在线观看入口| 欧美最新免费一区二区三区| 国产高清有码在线观看视频| 可以在线观看毛片的网站| 少妇丰满av| 如何舔出高潮| 国产不卡一卡二| 色尼玛亚洲综合影院| 亚洲综合色惰| 国产一级毛片在线| 中文欧美无线码| 精品日产1卡2卡| 亚洲自拍偷在线| 成人av在线播放网站| 日本成人三级电影网站| 亚洲av熟女| 九草在线视频观看| 九九热线精品视视频播放| av专区在线播放| 91久久精品国产一区二区三区| kizo精华| 综合色丁香网| 成人高潮视频无遮挡免费网站| av视频在线观看入口| 婷婷亚洲欧美| 男女视频在线观看网站免费| 三级男女做爰猛烈吃奶摸视频| 国模一区二区三区四区视频| 变态另类成人亚洲欧美熟女| 亚洲av熟女| 成年版毛片免费区| videossex国产| 欧美日韩综合久久久久久| 亚洲最大成人手机在线| 99久久久亚洲精品蜜臀av| 99热这里只有是精品在线观看| av免费观看日本| 一级av片app| 久久国产乱子免费精品| 99热这里只有精品一区| 亚洲av男天堂| 亚洲在线观看片| 青春草亚洲视频在线观看| 国产精品乱码一区二三区的特点| 国模一区二区三区四区视频| 日韩欧美 国产精品| 麻豆国产av国片精品| 22中文网久久字幕| 欧美色欧美亚洲另类二区| 男人舔女人下体高潮全视频| 亚洲人成网站高清观看| 九九热线精品视视频播放| 看免费成人av毛片| 最后的刺客免费高清国语| ponron亚洲| 欧美色欧美亚洲另类二区| 成人综合一区亚洲| 国产精品人妻久久久影院| 欧美变态另类bdsm刘玥| 欧美一区二区亚洲| 国语自产精品视频在线第100页| 日本欧美国产在线视频| 国产高潮美女av| 精品人妻熟女av久视频| 蜜臀久久99精品久久宅男| 国产v大片淫在线免费观看| 18禁黄网站禁片免费观看直播| 内地一区二区视频在线| 久久中文看片网| 国产av麻豆久久久久久久| 美女xxoo啪啪120秒动态图| 岛国毛片在线播放| 亚洲精品影视一区二区三区av| 国产白丝娇喘喷水9色精品| 欧美xxxx黑人xx丫x性爽| 国产极品精品免费视频能看的| 人妻少妇偷人精品九色| 日日啪夜夜撸| 久久草成人影院| 日韩欧美一区二区三区在线观看| 午夜精品一区二区三区免费看| a级毛片a级免费在线| 两性午夜刺激爽爽歪歪视频在线观看| 国产精品久久久久久av不卡| 成年女人看的毛片在线观看| 久久人人爽人人片av| 97超视频在线观看视频| 国内精品一区二区在线观看| 久久精品夜色国产| 国产成人91sexporn| 搞女人的毛片| 国产成人aa在线观看| av免费在线看不卡| 天天一区二区日本电影三级| 一进一出抽搐动态| 99久久中文字幕三级久久日本| 久久精品国产亚洲av涩爱 | 成年女人看的毛片在线观看|