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

    質(zhì)心動(dòng)態(tài)路徑規(guī)劃的變轉(zhuǎn)速齒輪箱時(shí)頻脊線索引算法

    2024-12-03 00:00:00張伯麟,萬書亭,趙曉艷,張雄,顧曉輝
    振動(dòng)工程學(xué)報(bào) 2024年6期

    摘要: 針對(duì)強(qiáng)噪聲下難以估計(jì)信號(hào)瞬時(shí)頻率的問題,提出了基于質(zhì)心動(dòng)態(tài)路徑規(guī)劃(DPPB)的變轉(zhuǎn)速齒輪箱時(shí)頻脊線索引算法。該算法在剖析多路徑匹配追蹤(MMP)脊線索引算法及其在強(qiáng)噪聲下失效原因的基礎(chǔ)上,通過對(duì)MMP算法得到的脊線集加窗求質(zhì)心,構(gòu)建信號(hào)的脊線質(zhì)心稀疏矩陣,并針對(duì)質(zhì)心稀疏矩陣設(shè)計(jì)動(dòng)態(tài)路徑規(guī)劃函數(shù)索引脊線上的質(zhì)心點(diǎn),根據(jù)脊線代價(jià)函數(shù)值計(jì)算最優(yōu)時(shí)頻脊線。將相似度系數(shù)Ra和置信度σRa作為脊線提取效果的衡量指標(biāo),通過仿真和試驗(yàn)驗(yàn)證了DPPB算法可有效提取強(qiáng)噪聲信號(hào)的時(shí)頻脊線,且不同程度噪聲下的可靠性和魯棒性均優(yōu)于峰值索引算法和MMP脊線索引算法。

    關(guān)鍵詞: 故障診斷; 齒輪箱; 質(zhì)心動(dòng)態(tài)路徑規(guī)劃(DPPB); 時(shí)頻脊線索引; 多路徑匹配追蹤(MMP); 時(shí)頻分析

    中圖分類號(hào): TH165+.3; TH132.41 文獻(xiàn)標(biāo)志碼: A 文章編號(hào): 1004-4523(2024)06-1077-12

    DOI:10.16385/j.cnki.issn.1004-4523.2024.06.018

    引 言

    齒輪箱在變轉(zhuǎn)速工況下會(huì)產(chǎn)生與轉(zhuǎn)速相關(guān)的非平穩(wěn)振動(dòng)信號(hào),使用傳統(tǒng)的平穩(wěn)信號(hào)分析方法難以診斷齒輪箱故障[1]。瞬時(shí)轉(zhuǎn)速不僅是反映齒輪箱運(yùn)行狀態(tài)的重要參數(shù),也是階次跟蹤算法、廣義解調(diào)分析等非平穩(wěn)信號(hào)分析方法的必需參量[2?4]。所以準(zhǔn)確獲得變轉(zhuǎn)速齒輪箱瞬時(shí)轉(zhuǎn)速對(duì)齒輪箱健康狀態(tài)監(jiān)測(cè)和故障診斷具有重要意義。

    利用高精度轉(zhuǎn)速測(cè)量設(shè)備可獲得振動(dòng)信號(hào)對(duì)應(yīng)的轉(zhuǎn)速,但是由于實(shí)際工作環(huán)境的制約,大部分齒輪箱不具備安裝轉(zhuǎn)速測(cè)量設(shè)備的條件。齒輪箱振動(dòng)信號(hào)中的特征頻率通常與轉(zhuǎn)速相關(guān),所以通過估計(jì)瞬時(shí)頻率(Instantaneous Frequency,IF)可達(dá)到求取瞬時(shí)轉(zhuǎn)速的目的[5]。郭瑜等[6]針對(duì)變轉(zhuǎn)速旋轉(zhuǎn)機(jī)械,利用峰值索引算法在振動(dòng)信號(hào)時(shí)頻分布中獲得一階轉(zhuǎn)速對(duì)應(yīng)的IF時(shí)頻脊線,實(shí)現(xiàn)了無轉(zhuǎn)速計(jì)條件下估計(jì)振動(dòng)信號(hào)的轉(zhuǎn)頻信息。峰值索引算法計(jì)算簡(jiǎn)單但對(duì)噪聲較為敏感,在實(shí)際工程中往往難以得到滿足要求的時(shí)頻曲線。為了提高估計(jì)算法對(duì)噪聲的魯棒性, latsenko等[7]以頻率平均值和標(biāo)準(zhǔn)差構(gòu)造自適應(yīng)代價(jià)函數(shù)代替時(shí)頻峰值進(jìn)行脊線索引,引入跳頻懲罰保證脊線的連續(xù)性,一定程度上降低了噪聲對(duì)脊線索引的影響。江星星等[8]提出了一種雙向搜索時(shí)頻脊融合算法,運(yùn)用代價(jià)函數(shù)脊線索引方法從正反兩個(gè)方希搜索同一IF分量的時(shí)頻脊線,通過脊線融合算法合并最優(yōu)脊線段從而消除噪聲引起的局部誤差。Ding等[9]則根據(jù)軸承振動(dòng)信號(hào)的頻率特征,運(yùn)用概率密度函數(shù)將基頻與其對(duì)應(yīng)的倍頻進(jìn)行融合,從而達(dá)到修正時(shí)頻脊線、提高轉(zhuǎn)速估計(jì)精度的目的。由于受到Heisenberg不確定原理的影響,時(shí)頻分布中頻率發(fā)散、多分量頻率信號(hào)混疊等現(xiàn)象嚴(yán)重影響時(shí)頻脊線的索引結(jié)果。張炎等[10]運(yùn)用同步壓縮小波變換(Synchrosqueezing Wavelet Transform,SWT)對(duì)軸承振動(dòng)信號(hào)的包絡(luò)譜進(jìn)行時(shí)頻重排,通過提高時(shí)頻包絡(luò)譜的分辨率和聚集性,減小頻譜發(fā)散對(duì)脊線提取的影響,從而準(zhǔn)確提取軸承的特征頻率曲線。上述算法一定程度上提高了時(shí)頻脊線索引的精度,但對(duì)索引初始點(diǎn)的依賴性很高,若選擇錯(cuò)誤的索引初始點(diǎn)便無法提取正確的時(shí)頻脊線。多路徑匹配追蹤(Multipath Matching Pursuit,MMP)[11]通過建立基于樹的路徑搜索模型,運(yùn)用貪心算法搜尋最優(yōu)路徑,降低初始搜索點(diǎn)選擇對(duì)搜索結(jié)果的影響。相繼有學(xué)者利用同步壓縮變換(Synchro?Squeezing Transform, SST)[12]、Multisynchro?squeezing Transform(MSST)[13]、同步提取變換(Synchro?Extracting Transform,SET)[14]、Synchro?Reassigning Transform(SRT)[15]算法獲得變轉(zhuǎn)速齒輪箱振動(dòng)信號(hào)的重排時(shí)頻分布,運(yùn)用MMP算法提取到了較為準(zhǔn)確的瞬時(shí)頻率。然而對(duì)于強(qiáng)噪聲信號(hào),重排時(shí)頻分布會(huì)出現(xiàn)大量噪聲分量,MMP算法在干擾下依然難以識(shí)別正確的時(shí)頻脊線。

    目前脊線索引算法均是以時(shí)頻分布上的離散時(shí)頻點(diǎn)作為索引對(duì)象,強(qiáng)噪聲產(chǎn)生的局部幅值極值會(huì)對(duì)脊線索引造成很強(qiáng)的干擾;同時(shí)現(xiàn)有算法將局部最優(yōu)點(diǎn)的集合作為最終的脊線,所以其結(jié)果不一定滿足全局最優(yōu)。針對(duì)現(xiàn)有算法的不足,本文提出一種質(zhì)心動(dòng)態(tài)路徑規(guī)劃(Dynamic Path Planning for Barycenter, DPPB)時(shí)頻脊線索引算法。算法在MMP算法的基礎(chǔ)上構(gòu)建脊線質(zhì)心稀疏矩陣,將索引的對(duì)象由時(shí)頻離散點(diǎn)轉(zhuǎn)換為質(zhì)心點(diǎn),減少噪聲對(duì)脊線索引的干擾并降低時(shí)頻分布的維度;運(yùn)用動(dòng)態(tài)路徑規(guī)劃算法索引時(shí)頻脊線,盡可能使索引結(jié)果達(dá)到全局最優(yōu)。

    1 MMP脊線索引算法

    1.1 MMP算法原理及實(shí)現(xiàn)

    對(duì)于信號(hào),窗函數(shù),通過滑動(dòng)窗口可得到不同時(shí)間段信號(hào)的頻譜圖,從而獲得信號(hào)的時(shí)頻分布如下式所示:

    (1)

    式中 和表示時(shí)間;表示頻率;表示g的共軛。

    根據(jù)時(shí)頻曲線附近局部能量突出且具有一定連續(xù)性的特點(diǎn),建立最優(yōu)脊線路徑代價(jià)函數(shù)[16]:

    (2)

    式中 表示幅值對(duì)代價(jià)函數(shù)的權(quán)重系數(shù);表示第k階導(dǎo)數(shù)對(duì)代價(jià)函數(shù)的權(quán)重系數(shù);k為求導(dǎo)階數(shù)。

    對(duì)式(2)進(jìn)行離散變換,忽略高階項(xiàng)的影響,同時(shí)設(shè)置脊線搜索頻率半徑為,得到如下式所示的脊線索引代價(jià)函數(shù):

    (3)

    式中 為權(quán)重系數(shù),用于確定幅值(時(shí)頻分布)和頻率變化量對(duì)代價(jià)函數(shù)的貢獻(xiàn)率。由于噪聲的影響,時(shí)頻分布中可能出現(xiàn)噪聲時(shí)頻點(diǎn)幅值高于脊線上的幅值,若索引初始點(diǎn)落在噪聲點(diǎn)上則無法索引到正確脊線。

    MMP脊線索引算法通過選擇多個(gè)索引初始點(diǎn)的方法消除索引結(jié)果過于依賴初始點(diǎn)的問題。對(duì)于時(shí)頻分布構(gòu)建基于樹的脊線索引模型,選取n個(gè)起始點(diǎn)作為路徑搜索樹的支集起點(diǎn),以代價(jià)函數(shù)C更新支集原子,通過貪心算法搜索每個(gè)支集上的原子,得到個(gè)支集作為時(shí)頻脊線,然后根據(jù)下式確定最優(yōu)時(shí)頻脊線:

    (4)

    式中 為離散信號(hào)總長(zhǎng)度。

    MMP脊線索引算法具體步驟如下:

    (1) 對(duì)于時(shí)頻矩陣,取個(gè)索引初始點(diǎn)將時(shí)間平均分為n段,每一段的時(shí)間間隔為,每個(gè)初始點(diǎn)對(duì)應(yīng)的時(shí)刻為,初始化;

    (2) 令,找到第個(gè)初始時(shí)刻的最大幅值點(diǎn)作為脊線搜索的起始點(diǎn);

    (3) 正向搜索:令,根據(jù)式(3)以半徑搜索處的脊線點(diǎn);

    (4) 重復(fù)步驟(3),直到停止迭代,;

    (5) 負(fù)向搜索:令,根據(jù)式(3)以半徑搜索處的脊線所在位置,;

    (6) 重復(fù)步驟(5) ,直到停止迭代,得到時(shí)頻脊線;

    (7) 重復(fù)步驟(2)~(6) ,直到停止迭代;根據(jù)式(4)得到最終的脊線。

    1.2 強(qiáng)噪聲下MMP算法失效機(jī)理分析

    圖1(a)為信噪比為-10 dB的強(qiáng)噪聲下利用MMP算法得到的時(shí)頻脊線,其中黑色虛線代表理論時(shí)頻曲線,紅色實(shí)線為MMP算法提取的結(jié)果。從圖1(a)中可以看出,時(shí)間在0~7.5 s時(shí)與基本重合,提取效果較好;在7.5 s時(shí)與發(fā)生偏離,然后便沿著偏離后的結(jié)果繼續(xù)搜尋,致使7.5 s后未能提取到有效的時(shí)頻脊線。

    圖1(b)為圖1(a)中位置A的局部放大圖,紅色虛線表示式(3)中的搜索頻率半徑。在時(shí)刻,時(shí)頻點(diǎn)a雖然在附近,但在噪聲的影響下幅值并不是搜索半徑內(nèi)的最大值,存在一個(gè)遠(yuǎn)離理論曲線的時(shí)頻點(diǎn)b,其幅值大于,所以在時(shí)刻會(huì)將b作為上的點(diǎn)。當(dāng)脊線索引至?xí)r刻,附近的時(shí)頻點(diǎn)c的幅值高于遠(yuǎn)離的d處幅值,然而此刻c點(diǎn)已經(jīng)超出的范圍,算法只能在的范圍內(nèi)選擇最大幅值點(diǎn)d作為f上的點(diǎn)而無法索引到正確的時(shí)頻點(diǎn)c,使逐步偏離正確的索引方向。

    由于MMP算法在索引處的脊線點(diǎn)時(shí)僅根據(jù)確定的索引范圍,得到關(guān)于的局部最優(yōu)值;若在噪聲的影響下偏離,的索引范圍也隨之發(fā)生偏離;當(dāng)?shù)钠屏吭谒饕秶鷥?nèi),索引結(jié)果僅在附近波動(dòng);當(dāng)?shù)睦塾?jì)偏移量超出索引范圍,偏離了并且偏離量可能越來越大,得到如圖1(a)中的結(jié)果。如果根據(jù)已搜索到的所有脊線點(diǎn)確定的索引范圍,可以有效抑制局部最優(yōu)導(dǎo)致的脊線偏移,但是在高采樣頻率下數(shù)據(jù)點(diǎn)多,影響計(jì)算效率。所以本文通過構(gòu)造脊線質(zhì)心稀疏矩陣重構(gòu)時(shí)頻分布,減少噪聲干擾并降低脊線索引的數(shù)據(jù)量,進(jìn)而利用動(dòng)態(tài)路徑規(guī)劃算法解決脊線索引時(shí)容易陷入局部最優(yōu)的問題。

    2 DPPB時(shí)頻脊線索引算法

    2.1 構(gòu)建脊線質(zhì)心稀疏矩陣

    由第1節(jié)可知,MMP算法得到最終脊線之前會(huì)生成一個(gè)脊線集,如圖2所示。圖2中彩色細(xì)實(shí)線為算法索引的所有脊線,紅色粗實(shí)線為最終的提取結(jié)果,紫色虛線區(qū)域?yàn)槔碚摃r(shí)頻曲線所在的搜索范圍。

    雖然噪聲導(dǎo)致脊線集中不一定存在一條脊線能夠較好地描述信號(hào)真實(shí)的時(shí)頻分布,但是對(duì)于,存在和集合,使得:

    (5)

    即時(shí)頻脊線上存在某一段正好落在內(nèi),所有滿足式(5)的脊線點(diǎn)集合便構(gòu)成目標(biāo)脊線。

    質(zhì)心可以反映質(zhì)點(diǎn)集的整體分布情況,將式(2)中的視為脊線的質(zhì)量,設(shè)窗函數(shù)為,時(shí)窗半徑為,窗口平移步長(zhǎng)為λ,那么第k條脊線的質(zhì)心可由下式求得:

    (6)

    針對(duì)離散的時(shí)頻分布,將式(6)改寫為離散形式:

    (7)

    式中 和分別為質(zhì)心對(duì)應(yīng)的時(shí)間和頻率;為質(zhì)心質(zhì)量;m為時(shí)間序列;為質(zhì)心序列,其中,表示向上取整運(yùn)算,為采樣頻率。

    截取一段脊線加窗求質(zhì)心的結(jié)果如圖3所示,其中不同顏色代表加窗后的脊線段,“○”為該段脊線的質(zhì)心,藍(lán)色虛線表示理論曲線。從圖3中可以看出,加窗質(zhì)心運(yùn)算弱化了異常脊線點(diǎn)對(duì)時(shí)頻分布的影響,從而達(dá)到消除噪聲干擾的目的。同時(shí),質(zhì)心點(diǎn)反映了窗口內(nèi)脊線點(diǎn)的整體分布,所以通過質(zhì)心點(diǎn)描述時(shí)頻脊線可以降低時(shí)頻脊線的維度。

    利用式(7)對(duì)圖2中的n段脊線加窗求質(zhì)心,得到一個(gè)由質(zhì)心點(diǎn)構(gòu)成的稀疏矩陣V如圖4所示,圖4(b)為4(a)中區(qū)域B的局部放大圖。圖中“·”為質(zhì)心點(diǎn)的位置,顏色深淺代表該質(zhì)心質(zhì)量的大小。與圖1(a)對(duì)比,圖4中的時(shí)頻脊線更加清晰。

    采用二維數(shù)組存儲(chǔ)稀疏矩陣會(huì)浪費(fèi)存儲(chǔ)單元存放零元素,在運(yùn)算中需要花費(fèi)大量時(shí)間對(duì)零元素進(jìn)行無效運(yùn)算。三元組存儲(chǔ)格式(Coordinate,COO)是一種直觀、簡(jiǎn)單的稀疏矩陣存儲(chǔ)格式,分別將二維數(shù)組中非零元素的行、列和數(shù)值存在三個(gè)一維矩陣中,極大地壓縮了原始數(shù)據(jù)量[17]。但是COO存儲(chǔ)數(shù)據(jù)不考慮存儲(chǔ)順序,在數(shù)據(jù)讀取時(shí)需要遍歷整個(gè)矩陣索引需要的數(shù)據(jù),一定程度上降低運(yùn)算效率。針對(duì)時(shí)頻脊線索引的特點(diǎn),按照時(shí)間順序運(yùn)用COO方法存儲(chǔ)稀疏矩陣,減少索引時(shí)的數(shù)據(jù)讀取量。首先通過COO方法將b時(shí)刻的質(zhì)心儲(chǔ)存在三個(gè)一維矩陣中,即

    (8)

    然后將所有時(shí)刻的存儲(chǔ)矩陣組成質(zhì)心稀疏矩陣V:

    (9)

    通過構(gòu)建質(zhì)心稀疏矩陣將時(shí)頻脊線索引對(duì)象由時(shí)頻分布轉(zhuǎn)化為稀疏矩陣V,保留了時(shí)頻分布中的關(guān)鍵信息,極大地壓縮了數(shù)據(jù)量;在索引某一時(shí)刻的脊線點(diǎn)時(shí),只需在該時(shí)刻的COO矩陣中查找目標(biāo)元素,從而避免遍歷稀疏矩陣的全部數(shù)據(jù),達(dá)到降低運(yùn)算量的目的。

    2.2 構(gòu)造動(dòng)態(tài)路徑規(guī)劃脊線索引函數(shù)

    路徑規(guī)劃可通過現(xiàn)有的數(shù)據(jù)集預(yù)測(cè)接下來一段時(shí)間的數(shù)據(jù)分布,多項(xiàng)式擬合是常用的路徑規(guī)劃方法。但是數(shù)據(jù)量和擬合階數(shù)對(duì)多項(xiàng)式擬合的結(jié)果影響較大,決定了擬合結(jié)果的精度和泛化能力。為了提高脊線點(diǎn)預(yù)測(cè)的魯棒性且盡可能保留脊線局部的細(xì)節(jié)信息,本文采用了動(dòng)態(tài)路徑規(guī)劃的脊線點(diǎn)索引算法。

    設(shè)置閾值確定擬合點(diǎn)數(shù)量和擬合階數(shù),當(dāng)已搜索到脊線的質(zhì)心個(gè)數(shù)時(shí)選用低階數(shù)對(duì)序列的質(zhì)心集進(jìn)行擬合;當(dāng)時(shí)選用高階數(shù)對(duì)序列的質(zhì)心集進(jìn)行擬合;然后計(jì)算序列為的預(yù)測(cè)值,如下式所示:

    (10)

    式中 ,分別為預(yù)測(cè)點(diǎn)的時(shí)間集和頻率集;;。

    定義新的代價(jià)函數(shù)如下式所示:

    (11)

    式中 為時(shí)刻的預(yù)測(cè)頻率;表示距離。代價(jià)函數(shù)是一個(gè)與距離和質(zhì)量有關(guān)的函數(shù)。

    將式(11)得到的代價(jià)函數(shù)值取最大值時(shí)的質(zhì)心點(diǎn)作為的脊線點(diǎn),并把此時(shí)的作為脊線點(diǎn)的幅值,如下式所示:

    (12)

    將每條脊線各點(diǎn)的幅值合作為脊線的判別參數(shù),選擇最大幅值合的脊線作為最終脊線,如下式所示:

    (13)

    2.3 DPPB脊線索引算法實(shí)現(xiàn)

    對(duì)于一段時(shí)間長(zhǎng)度為的信號(hào),運(yùn)用DPPB算法搜索時(shí)頻脊線的具體步驟如下:

    (1) 根據(jù)第1節(jié)的MMP脊線索引算法搜索時(shí)頻分布的時(shí)頻脊線,得到脊線集合。

    (2) 設(shè)置窗函數(shù),時(shí)窗半徑,窗口平移步長(zhǎng)λ,對(duì)脊線集加窗求質(zhì)心,構(gòu)建質(zhì)心稀疏矩陣,計(jì)算初始化,確定擬合階數(shù)閾值m。

    (3) 令,搜索中的最大值所對(duì)應(yīng)的質(zhì)點(diǎn)作為第i條脊線的搜尋起始點(diǎn),即。

    (4) 正向搜索:,若,執(zhí)行步驟(5);否則,令,執(zhí)行步驟(6)。

    (5) 由式(10)確定預(yù)測(cè)值,由式(11)計(jì)算時(shí)刻的代價(jià)函數(shù)集,由式(12)搜索脊線點(diǎn);若,重復(fù)步驟(4)。

    (6) 反向搜索:,若,執(zhí)行步驟(7);否則,執(zhí)行步驟(8)。

    (7) 由式(10)確定預(yù)測(cè)值,由式(11)計(jì)算時(shí)刻的代價(jià)函數(shù)集,由式(12)搜索脊線點(diǎn),重復(fù)步驟(6)。

    (8) 令,重復(fù)步驟(3)~(6),直到停止。由式(13)得到最終的時(shí)頻脊線。

    DPPB算法流程圖如圖5所示。

    2.4 關(guān)于算法的討論

    2.4.1 參數(shù)設(shè)置對(duì)結(jié)果的影響

    在DPPB算法中,搜索半徑、窗寬和步長(zhǎng)是算法的三個(gè)關(guān)鍵參數(shù)。搜索半徑的大小影響MMP算法對(duì)噪聲的魯棒性。減小可以降低索引區(qū)域內(nèi)的噪聲點(diǎn)數(shù)量,避免脊線波動(dòng),但同時(shí)也降低了索引到正確脊線點(diǎn)的概率,容易在噪聲的干擾下偏離正確索引路徑;增大可以增加索引到正確脊線點(diǎn)的概率,當(dāng)索引路徑偏離理論脊線時(shí)有更大的機(jī)會(huì)重新返回正確的索引路徑,但是也增加了噪聲點(diǎn)的數(shù)量,容易使提取結(jié)果產(chǎn)生波動(dòng)。加窗求質(zhì)心的目的是降低噪聲對(duì)脊線提取的干擾,使質(zhì)心點(diǎn)能夠準(zhǔn)確描述理論脊線。窗口大小和步長(zhǎng)會(huì)影響質(zhì)心計(jì)算的準(zhǔn)確性、抗噪聲干擾能力和復(fù)雜度。隨著增大,質(zhì)心計(jì)算的抗噪聲干擾能力增強(qiáng),而準(zhǔn)確性則先增高后降低;增大時(shí),質(zhì)心計(jì)算的準(zhǔn)確性提高,但是提升效果會(huì)隨的減小而降低。同時(shí),和影響算法的復(fù)雜度,增大和會(huì)增加運(yùn)算量。

    目前尚未有合適的理論確定數(shù)據(jù)的最優(yōu)搜索半徑、窗寬和步長(zhǎng),本文僅根據(jù)經(jīng)驗(yàn)和試驗(yàn)分析獲得各參數(shù)的較優(yōu)取值范圍。對(duì)于,一般取10~35 Hz效果較好。對(duì)于,如果待提取的脊線變化較為劇烈,通常選擇0.1~0.5 s的較小窗口;如果變化較為緩慢,通常選擇0.5~1 s的較大窗口。對(duì)于,原則上越小越好,但是考慮到運(yùn)算量,一般選擇。

    2.4.2 算法的計(jì)算成本

    本文所提算法的計(jì)算成本主要包括MMP索引、質(zhì)心求解和動(dòng)態(tài)路徑規(guī)劃三部分。若信號(hào)時(shí)頻分布的數(shù)據(jù)量為×,為離散信號(hào)總長(zhǎng)度,為頻率維度,算法的計(jì)算成本可以表示為:

    (14)

    式中 ,,和分別為算法總成本、MMP索引成本、質(zhì)心計(jì)算成本和動(dòng)態(tài)路徑規(guī)劃成本;為脊線段的數(shù)量;為個(gè)點(diǎn)的多項(xiàng)式擬合計(jì)算成本。

    所提算法以較高的計(jì)算成本來換取較好的提取性能,所以計(jì)算成本較大。特別注意的是,的系數(shù)對(duì)計(jì)算成本產(chǎn)生較大的影響,所以在保證索引精度的前提下,合理減小和、增大可以有效降低算法的計(jì)算成本。

    3 仿真分析

    3.1 構(gòu)建仿真模型

    設(shè)仿真信號(hào)x及瞬時(shí)頻率曲線f如下式所示:

    (15)

    設(shè)置采樣頻率fs=800 Hz,采樣時(shí)間t=10 s,信號(hào)的時(shí)域圖和STFT變換后的時(shí)頻分布如圖6所示。

    對(duì)信號(hào)加高斯白噪聲后,取,運(yùn)用第1節(jié)MMP脊線提取算法得到脊線集如圖7所示。

    由式(7)構(gòu)建質(zhì)心稀疏矩陣V,其中總采樣點(diǎn)數(shù)N=fs·t=8000,窗寬度g=0.5 s,步長(zhǎng)λ=0.1 s,時(shí)間序列總長(zhǎng)度B=t/λ=100。

    令式(10)中的,,,,,則仿真模型的動(dòng)態(tài)路徑規(guī)劃函數(shù)可表示為:

    (16)

    根據(jù)式(11)構(gòu)造仿真信號(hào)的代價(jià)函數(shù)為:

    (17)

    式中 為歸一化質(zhì)心質(zhì)量;為距離權(quán)重,用來平衡歸一化質(zhì)心質(zhì)量和距離對(duì)代價(jià)函數(shù)的貢獻(xiàn)度,通常的取值范圍為。

    3.2 時(shí)頻脊線提取效果衡量指標(biāo)

    為了更直觀地分析不同算法的優(yōu)劣,此處構(gòu)造兩個(gè)量化衡量指標(biāo):相似度系數(shù)和置信度。

    絕對(duì)值倒數(shù)法是確定兩個(gè)向量間相似度的一種方法,對(duì)于提取的第j條脊線上任意一點(diǎn)的頻率與對(duì)應(yīng)的理論值之間的相關(guān)程度可以由下式確定:

    (18)

    式中 c為誤差允許值,表示接受與理論值的距離不大于c的提取結(jié)果。

    對(duì)于某一算法O提取的第j條脊線與理論曲線的相似度,計(jì)算公式如下:

    (19)

    式中 Ra∈[0,1],Ra越接近1,說明算法提取的脊線與理論曲線相似度越高,脊線提取效果越好;反之Ra越接近0,則說明相似度越低,提取效果越差。

    圖8所示為不同相似度Ra時(shí)的時(shí)頻脊線對(duì)比圖,圖中誤差上、下限之間為允許誤差c的范圍。當(dāng)Ra≥0.9時(shí),脊線與理論時(shí)頻曲線基本吻合,提取結(jié)果誤差較??;當(dāng)0.9≤Ra≤0.8時(shí),脊線在理論曲線附近波動(dòng)變大,且存在少量脊線點(diǎn)超出了誤差范圍,但是脊線基本落在誤差范圍以內(nèi);當(dāng)Ra≤0.7時(shí),脊線與理論曲線保持相同的變化趨勢(shì),但是有部分脊線嚴(yán)重超出誤差范圍,提取效果欠佳。

    運(yùn)用算法O對(duì)M個(gè)信號(hào)進(jìn)行脊線提取,獲得M條脊線,根據(jù)式(19)計(jì)算得到所有脊線與理論曲線的相關(guān)系數(shù),算法O的平均相似度系數(shù)為:

    (20)

    對(duì)于任意,存在k條脊線滿足,則置信水平在上的置信度為:

    (21)

    平均相似系數(shù)表示提取結(jié)果與理論結(jié)果偏差的統(tǒng)計(jì)平均值,可以反映算法的穩(wěn)定程度。置信度表示將與理論結(jié)果的相似度在的提取結(jié)果視為合格結(jié)果,則在整個(gè)試驗(yàn)結(jié)果中合格的提取結(jié)果所占的比率為,值越大,代表算法的提取效果越好。

    3.3 仿真結(jié)果及分析

    脊線提取結(jié)果及其與理論曲線對(duì)比結(jié)果如圖9所示。圖9(a)中紅色曲線為提取到的時(shí)頻脊線,脊線索引到了理論脊線附近所有的高幅值時(shí)頻點(diǎn);從圖9(b)中可以看到,提取到的脊線與理論曲線基本重合,表明本算法可以提取得到較為準(zhǔn)確的時(shí)頻脊線。

    為了分析噪聲對(duì)脊線提取結(jié)果的影響,運(yùn)用峰值索引、MMP脊線索引、DPPB脊線索引三種算法對(duì)信噪比的信號(hào)進(jìn)行脊線提取,結(jié)果如圖10所示。

    在圖10(a)中,對(duì)于信噪比SNR=-5 dB的低噪聲信號(hào),3種算法均能較好地提取時(shí)頻脊線。在SNR=-8和-10 dB的較強(qiáng)噪聲時(shí)(圖10(b)和(c)所示),峰值索引算法得到的脊線已經(jīng)完全偏離理論曲線,MMP算法得到的脊線與理論曲線有相同的變化趨勢(shì),但存在局部誤差而偏離理論曲線,DPPB算法提取的脊線與理論曲線基本重合。在SNR=-12 dB的強(qiáng)噪聲下(圖10(d)所示),其他兩種算法已經(jīng)完全偏離理論曲線,DPPB算法可得到與理論曲線變化趨勢(shì)相同的脊線,僅存在局部誤差。由此可以初步得出,峰值索引算法和MMP算法容易受噪聲的影響,強(qiáng)噪聲下會(huì)失去脊線提取能力,而DPPB算法可以在強(qiáng)噪聲下有效提取時(shí)頻脊線,提取效果和對(duì)噪聲的魯棒性優(yōu)于其他兩種算法。

    重復(fù)進(jìn)行圖10中的仿真試驗(yàn),試驗(yàn)次數(shù)M=10000,根據(jù)式(20)計(jì)算三種算法的平均相關(guān)系數(shù)如表1所示,其中取 Hz。

    由表1可知,當(dāng)SNR=-5 dB時(shí),三種算法的平均相關(guān)系數(shù)均大于0.9,提取結(jié)果比較穩(wěn)定。峰值索引算法在SNR=-8,-10和-12 dB的強(qiáng)噪聲下值下降明顯;MMP算法在三種強(qiáng)噪聲下的值雖然高于峰值索引算法,但是也有較大幅度的下降;DPPB算法在三種強(qiáng)噪聲下的值降幅較小且均在0.8以上。隨著噪聲的增強(qiáng),三種算法的值均有不同程度的下降,但是DPPB算法的值始終高于其他兩種算法,且噪聲越強(qiáng)越明顯,說明對(duì)噪聲的魯棒性優(yōu)于其他算法。

    根據(jù)圖8的仿真結(jié)果,取,由式(21)計(jì)算置信水平在[0.9,1]上的置信度,結(jié)果如表2所示。

    由表2可知,在SNR=-5 dB時(shí),MMP脊線索引和DPPB算法的置信度分別為98.01%和100%,可以準(zhǔn)確地提取到符合要求的時(shí)頻脊線;峰值索引算法的也達(dá)到80%,具備較好的提取效果。當(dāng)噪聲強(qiáng)度高于-8 dB時(shí),峰值索引算法的均為0,已經(jīng)無法提取到符合要求的脊線。MMP算法在SNR=-8 dB時(shí)為54.75%,接近一半的提取結(jié)果不符合要求,在SNR=-10和-12 dB時(shí)的分別為13.95%和2.87%,已經(jīng)很難提取到合格的脊線。DPPB算法在SNR=-8和-10 dB的較強(qiáng)噪聲條件下分別為99.91%和93.51%,合格脊線所占的比率均在90%以上,在SNR=-12 dB的強(qiáng)噪聲下為57.12%,依然有一半以上的提取結(jié)果符合要求。雖然噪聲的增強(qiáng)降低了各算法的提取效果,但DPPB算法的提取效果始終優(yōu)于其他兩種算法,且針對(duì)強(qiáng)噪聲具有良好的脊線提取能力。

    通過以上仿真結(jié)果和分析可以得出,DPPB算法可以抵抗信號(hào)中強(qiáng)噪聲的干擾,相比于峰值索引算法和MMP算法更容易從強(qiáng)噪聲信號(hào)中提取到有效的時(shí)頻脊線,對(duì)噪聲有更好的魯棒性。

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

    4.1 試驗(yàn)說明

    齒輪箱振動(dòng)信號(hào)中的嚙合頻率分量是分析故障的重要特征頻率之一,本文利用SpectraQuest公司研發(fā)的動(dòng)力傳動(dòng)故障診斷綜合試驗(yàn)臺(tái)(DDS)獲取不同轉(zhuǎn)速下齒輪箱的振動(dòng)信號(hào),通過提取不同轉(zhuǎn)速下的嚙合頻率時(shí)頻脊線驗(yàn)證所提算法的有效性,試驗(yàn)臺(tái)如圖11所示。

    其中齒輪箱5為二級(jí)減速器減速齒輪箱,由4組直齒輪組合而成,主要參數(shù)如表3所示。設(shè)置加速度振動(dòng)傳感器6的采樣頻率為25.6 kHz,采集齒輪箱徑向的振動(dòng)信號(hào)。

    由于DDS試驗(yàn)臺(tái)不易輸出非線性變化的轉(zhuǎn)速,所以試驗(yàn)?zāi)M三組線性變化的輸入轉(zhuǎn)速,分別為震蕩波動(dòng)曲線、三角形單波谷波動(dòng)曲線、三角形三波谷波動(dòng)曲線,如圖12所示。為了與時(shí)頻脊線對(duì)應(yīng),下文均用“轉(zhuǎn)頻”代替“轉(zhuǎn)速”。

    齒輪箱振動(dòng)信號(hào)中包含與轉(zhuǎn)頻相關(guān)的第i級(jí)傳動(dòng)齒輪嚙合頻率及其倍頻,與的關(guān)系如下式所示:

    (22)

    因此,獲得齒輪箱振動(dòng)信號(hào)的第i級(jí)傳統(tǒng)齒輪嚙合頻率的n倍頻,便可計(jì)算得到齒輪箱的輸入轉(zhuǎn)頻。

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

    為了計(jì)算方便,試驗(yàn)通過提取1級(jí)傳動(dòng)齒輪嚙合頻率的基頻計(jì)算輸入轉(zhuǎn)頻。三種不同輸入轉(zhuǎn)速下齒輪箱的振動(dòng)時(shí)域圖和STFT變換后附近的時(shí)頻分布如圖13所示。

    運(yùn)用DPPB算法提取不同轉(zhuǎn)頻曲線下嚙合頻率的時(shí)頻脊線,提取結(jié)果如圖14所示。

    在圖13(b)和(d)的時(shí)頻分布中,信號(hào)的SNR較高,噪聲對(duì)瞬時(shí)頻率的影響較小,可以觀察到較為明顯的時(shí)頻脊線;在圖13(f)中,頻率波動(dòng)部分受噪聲污染嚴(yán)重,SNR較低,難以觀察到時(shí)頻脊線的變化趨勢(shì)。通過圖14的提取結(jié)果可以看出,雖然三種轉(zhuǎn)頻下的信噪比差異較大,但是DPPB算法均能提取到較為準(zhǔn)確的時(shí)頻脊線,表明該算法對(duì)強(qiáng)噪聲有良好的魯棒性。

    為了驗(yàn)證DPPB算法提取轉(zhuǎn)頻的準(zhǔn)確性,將圖14得到的嚙合頻率和表3中齒輪箱參數(shù)及代入式(22)可計(jì)算三組信號(hào)的轉(zhuǎn)頻。同時(shí),運(yùn)用MMP算法提取的時(shí)頻脊線,計(jì)算轉(zhuǎn)頻;用試驗(yàn)臺(tái)中的角度編碼器測(cè)得齒輪箱輸出軸轉(zhuǎn)頻,根據(jù)傳動(dòng)比計(jì)算得到輸入軸實(shí)測(cè)轉(zhuǎn)頻。將,作為對(duì)照組與進(jìn)行對(duì)比,結(jié)果如圖15所示。

    由圖15可以看出,在震蕩波動(dòng)轉(zhuǎn)頻和三角形單波谷轉(zhuǎn)頻下,MMP脊線索引和DPPB兩種算法提取到的時(shí)頻脊線基本與實(shí)測(cè)轉(zhuǎn)頻重合(圖15(a)和(b)所示),但由MMP脊線索引算法提取的三角形單波谷波時(shí)頻脊線存在較小的誤差。在三波谷轉(zhuǎn)頻波動(dòng)下,MMP脊線索引算法提取的脊線在恒頻部分與重合較好,但波谷處與存在較大誤差;DPPB算法得到的脊線無論在恒頻部分還是波谷處均與重合較好,只有局部存在較小的誤差(如圖15(c)所示)。由此表明DPPB算法可以在強(qiáng)噪聲下較為準(zhǔn)確地估計(jì)信號(hào)的轉(zhuǎn)頻曲線。

    將實(shí)測(cè)轉(zhuǎn)頻作為理論輸入轉(zhuǎn)頻,根據(jù)式(19)計(jì)算三組信號(hào)在不同算法下提取結(jié)果的相似度Ra衡量MMP算法和DPPB算法的提取效果,如表4所示。

    由表4可以看出,MMP算法和DPPB算法得到的震蕩波動(dòng)信號(hào)轉(zhuǎn)頻的Ra值均為1,三角形單波谷波動(dòng)信號(hào)轉(zhuǎn)頻的Ra值分別為0.9939和1,均有較好的提取效果。對(duì)于三角形三波谷波動(dòng)信號(hào),MMP算法得到的轉(zhuǎn)頻的Ra值為0.7678,與實(shí)際轉(zhuǎn)頻的偏差較大,而DPPB算法得到的轉(zhuǎn)頻的Ra值為0.9669,遠(yuǎn)高于MMP算法,進(jìn)一步說明DPPB算法的有效性和對(duì)噪聲的魯棒性。

    5 結(jié) 論

    mhqytu3B96llJVmdX9+jYw==本文提出了一種改進(jìn)的時(shí)頻脊線索引算法來提高時(shí)頻估計(jì)的準(zhǔn)確性,主要結(jié)論如下:

    (1) 信號(hào)時(shí)頻分布在強(qiáng)噪聲下會(huì)產(chǎn)生多個(gè)幅值極值,MMP算法索引時(shí)頻脊線時(shí)容易受到噪聲極值點(diǎn)的干擾陷入局部最優(yōu),導(dǎo)致算法提取到錯(cuò)誤的時(shí)頻脊線。

    (2) 通過構(gòu)建脊線質(zhì)心稀疏矩陣,將脊線索引對(duì)象由時(shí)頻離散點(diǎn)轉(zhuǎn)化為質(zhì)心點(diǎn),降低了噪聲對(duì)脊線索引的干擾,同時(shí)壓縮了索引數(shù)據(jù)量,提高了運(yùn)算速率。

    (3) 提出了一種基于質(zhì)心動(dòng)態(tài)路徑規(guī)劃(DPPB)的時(shí)頻脊線索引算法,設(shè)計(jì)動(dòng)態(tài)路徑優(yōu)化函數(shù)和脊線索引代價(jià)函數(shù)識(shí)別質(zhì)心稀疏矩陣中的時(shí)頻脊線點(diǎn),從而避免了脊線索引過程中陷入局部最優(yōu),提高了算法對(duì)噪聲的魯棒性。通過仿真和試驗(yàn)分析,驗(yàn)證了DPPB算法索引時(shí)頻脊線的精度和魯棒性優(yōu)于傳統(tǒng)的脊線索引算法。所提算法可用于實(shí)際工程中變轉(zhuǎn)速齒輪箱的轉(zhuǎn)速估計(jì)和故障診斷,具有一定的工程應(yīng)用價(jià)值。

    參考文獻(xiàn):

    [1] 陳雪峰, 郭艷婕, 許才彬, 等. 風(fēng)電裝備故障診斷與健康監(jiān)測(cè)研究綜述[J]. 中國(guó)機(jī)械工程, 2020, 31(2): 175-189.

    Cheng Xuefeng, Guo Yanjie, Xu Caibin, et al. Review of fault diagnosis and health monitoring for wind power equipment[J]. China Mechanical Engineering, 2020, 31(2): 175-189.

    [2] 孫鑫暉, 董翔文, 郝木明. 基于小波脊提取的扭轉(zhuǎn)振動(dòng)測(cè)試方法研究[J]. 振動(dòng)、測(cè)試與診斷, 2020, 40(1): 135-139.

    Sun Xinhui, Dong Xiangwen, Hao Muming. Research on testing method of torsional vibration based on wavelet ridge extraction[J]. Journal of Vibration,Measurement & Diagnosis, 2020, 40(1): 135-139.

    [3] 趙德尊, 王天楊, 褚福磊. 基于自適應(yīng)廣義解調(diào)變換的滾動(dòng)軸承時(shí)變非平穩(wěn)故障特征提?。跩]. 機(jī)械工程學(xué)報(bào), 2020, 56(3): 80-87.

    Zhao Dezun, Wang Tianyang, Chu Fulei. Adaptive generalized demodulation transform based rolling bearing time-varying nonstationary fault feature extraction[J]. Journal of Mechanical Engineering, 2020, 56(3): 80-87.

    [4] 王曉龍, 閆曉麗, 何玉靈. 變速工況下基于IEWT能量階次譜的風(fēng)電機(jī)組軸承故障診斷[J]. 太陽能學(xué)報(bào), 2021, 42(4): 479-486.

    Wang Xiaolong, Yan Xiaoli, He Yuling. Fault diagnosis wind turbine bearing based on IEWT energy order spec-trum under variable speed condition[J]. Acta Energiae Solaris Sinica, 2021, 42(4): 479-486.

    [5] 陳是扦, 彭志科, 周鵬. 信號(hào)分解及其在機(jī)械故障診斷中的應(yīng)用研究綜述[J]. 機(jī)械工程學(xué)報(bào), 2020, 56(17): 91-107.

    Chen Shiqian, Peng Zhike, Zhou Peng. Review of signal decomposition theory and its applications in machine fault diagnosis[J]. Journal of Mechanicla Egnineering, 2020, 56(17): 91-107.

    [6] 郭瑜, 秦樹人, 湯寶平, 等. 基于瞬時(shí)頻率估計(jì)的旋轉(zhuǎn)機(jī)械階比跟蹤[J]. 機(jī)械工程學(xué)報(bào), 2003, 39(3): 32-36.

    Guo Yu, Qin Shuren, Tang Baoping, et al. Order tracking of rotating machinery based on instantaneous frequency estimation[J]. Journal of Mechanical Engineering, 2003, 39(3): 32-36.

    [7] Iatsenko D, McClintock P V E, Stefanovska A. Extrac-tion of instantaneous frequencies from ridges in time–frequency representations of signals[J]. Signal Processing, 2016, 125: 290-303.

    [8] 江星星, 李舜酩, 周東旺, 等. 時(shí)頻脊融合方法及時(shí)變工況行星齒輪箱故障識(shí)別[J]. 振動(dòng)工程學(xué)報(bào), 2017, 30(1): 127-134.

    Jiang Xingxing, Li Shunming, Zhou Dongwang, et al. Time-frequency ridge fusion method and defective identification of planetary gearbox running on time-varying condition[J]. Journal of Vibration Engineering, 2017, 30(1): 127-134.

    [9] Ding Rongmei, Shi Juanjuan, Jiang Xingxing, et al. Multiple instantaneous frequency ridge based integration strategy for bearing fault diagnosis under variable speed operations[J]. Measurement Science and Technology, 2018, 29(11): 115002.

    [10] 張焱, 何姝鋇, 王平, 等. 無轉(zhuǎn)速計(jì)下變工況滾動(dòng)軸承故障特征量化表征提?。跩]. 儀器儀表學(xué)報(bào), 2021, 42(8): 104-114.

    Zhang Yan, He Shubei, Wang Ping, et al. Tacholess quanti-tative characterization of rolling bearing fault feature under varying conditions[J]. Chinese Journal of Scientific Instrument, 2021, 42(8): 104-114.

    [11] Wu M, Wu F, Yang K, et al. A multipath matching pursuit algorithm based on improved-inner product matching criterion[C].2020 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC). IEEE, 2020: 1-5.

    [12] Fourer D, Auger F, Flandrin P. Recursive versions of the Levenberg-Marquardt reassigned spectrogram and of the synchrosqueezed STFT[C].2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2016: 4880-4884.

    [13] Yu G, Wang Z H, Zhao P. Multisynchrosqueezing transform[J]. IEEE Transactions on Industrial Electronics, 2018, 66(7): 5441-5455.

    [14] Yu G, Yu M J, Xu C Y. Synchroextracting transform[J]. IEEE Transactions on Industrial Electronics, 2017, 64(10): 8042-8054.

    [15] Li M F, Wang T Y, Kong Y, et al. Synchro-reassigning transform for instantaneous frequency estimation and signal reconstruction[J]. IEEE Transactions on Industrial Electronics, 2021, 69(7): 7263-7274.

    [16] 江星星. 齒輪箱關(guān)鍵部件非平穩(wěn)振動(dòng)信號(hào)分析及診斷方法研究[D]. 南京:南京航空航天大學(xué), 2016.

    Jiang Xingxing. Research on the nonstationary vibration signal of key parts of gearbox and its fault diagnosis method[D]. Nanjing:Nanjing University of Aeronautics and Astronautics, 2016.

    [17] 楊世偉, 蔣國(guó)平, 宋玉蓉, 等. 基于GPU的稀疏矩陣存儲(chǔ)格式優(yōu)化研究[J]. 計(jì)算機(jī)工程, 2019, 45(9): 23-31.

    Yang Shiwei, Jiang Guoping, Song Yurong, et al. Research on storage format optimization of sparse matrix based on GPU[J]. Computer Engineering, 2019, 45(9): 23-31.

    Time-frequency ridge index algorithm for gearbox under variable speed based on dynamic path planning of barycenter

    ZHANG Bo-lin1, WAN Shu-ting1, ZHAO Xiao-yan1, ZHANG Xiong1, GU Xiao-hui2

    (1.Hebei Key Laboratory of Electric Machinery Health Maintenance & Failure Prevention,North China Electric Power University, Baoding 071003, China;2.State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures,Shijiazhuang Tiedao University, Shijiazhuang 050043, China)

    Abstract: The paper proposes a time-frequency ridge index algorithm for gearboxes under variable speed conditions, based on Dynamic Path Planning of Barycenter (DPPB). This algorithm addresses the challenge of estimating the instantaneous frequency of signals in a high-noise environment. The algorithm builds upon the analysis of the Multi-Path Matching Pursuit (MMP) ridge index algorithm and its limitations under high noise. By adding windows to the ridge set obtained by the MMP algorithm, a ridge barycenter sparse matrix of the signal is constructed. A dynamic path planning function is then designed for the barycenter sparse matrix to index the barycenters on the ridge line. The optimal time-frequency ridge line is calculated based on the values of the ridge line cost function. The similarity coefficient Ra and confidence σRa are used as measures of the ridge extraction effect. Simulations and experiments indicate that the DPPB algorithm can effectively extract the time-frequency ridge of signals in high-noise environments, and it is more reliable and robust than the peak index algorithm and the MMP algorithm under various noise intensities.

    Key words: fault diagnosis; gearbox;dynamic path planning barycenter (DPPB); time-frequency ridge index; multipath matching pursuit (MMP); time-frequency analysis

    作者簡(jiǎn)介: 張伯麟(1993—),男,博士研究生。電話:(0312)7525428; E-mail:393696838@qq.com。

    通訊作者: 萬書亭(1970—),男,博士,教授,博士生導(dǎo)師。電話:(0312)7525455; E-mail:52450809@ncepu.edu.cn。

    国产精品二区激情视频| 免费女性裸体啪啪无遮挡网站| 视频在线观看一区二区三区| 一级毛片电影观看| 精品国产国语对白av| 999精品在线视频| 午夜福利一区二区在线看| 视频区欧美日本亚洲| 色精品久久人妻99蜜桃| 成人三级做爰电影| 国产在视频线精品| 国产精品久久久av美女十八| www.999成人在线观看| 大香蕉久久网| 1024视频免费在线观看| 国产男人的电影天堂91| 脱女人内裤的视频| 999久久久精品免费观看国产| 国产av一区二区精品久久| 久久精品国产亚洲av高清一级| 国产一区二区在线观看av| 欧美日韩福利视频一区二区| 青青草视频在线视频观看| 国产精品久久久久久精品电影小说| 国产亚洲精品第一综合不卡| 色94色欧美一区二区| 亚洲一卡2卡3卡4卡5卡精品中文| 一级毛片精品| 亚洲国产欧美日韩在线播放| 亚洲中文av在线| 国产av又大| 男女国产视频网站| 国产亚洲一区二区精品| 精品少妇黑人巨大在线播放| 欧美日韩福利视频一区二区| 久久精品国产综合久久久| 欧美 日韩 精品 国产| 一本一本久久a久久精品综合妖精| 亚洲av国产av综合av卡| 国产亚洲欧美精品永久| 久久久久国产精品人妻一区二区| 亚洲成人国产一区在线观看| 香蕉国产在线看| 制服人妻中文乱码| 999久久久精品免费观看国产| 欧美激情 高清一区二区三区| 男女边摸边吃奶| 久久久国产一区二区| 亚洲va日本ⅴa欧美va伊人久久 | 人人妻人人澡人人看| 精品少妇久久久久久888优播| 亚洲成人免费电影在线观看| 18禁国产床啪视频网站| 一区在线观看完整版| 国产精品一区二区免费欧美 | 99久久国产精品久久久| 窝窝影院91人妻| 国产成人欧美| 久久女婷五月综合色啪小说| 一本色道久久久久久精品综合| 9191精品国产免费久久| 国产成人系列免费观看| 久久久久久久久久久久大奶| 免费av中文字幕在线| 日韩 亚洲 欧美在线| 在线观看免费视频网站a站| 在线看a的网站| 亚洲第一欧美日韩一区二区三区 | 亚洲专区国产一区二区| 亚洲成av片中文字幕在线观看| 黄片大片在线免费观看| 国产av国产精品国产| 好男人电影高清在线观看| videosex国产| 精品久久蜜臀av无| 狂野欧美激情性bbbbbb| 91精品三级在线观看| 制服诱惑二区| 国产免费一区二区三区四区乱码| 国产在线观看jvid| 亚洲精品自拍成人| 脱女人内裤的视频| 一个人免费看片子| a在线观看视频网站| 91大片在线观看| 亚洲中文字幕日韩| 人人妻人人澡人人爽人人夜夜| 久久久精品国产亚洲av高清涩受| 欧美国产精品va在线观看不卡| 在线观看一区二区三区激情| 永久免费av网站大全| 高清av免费在线| 手机成人av网站| 国产精品香港三级国产av潘金莲| 亚洲一区二区三区欧美精品| 脱女人内裤的视频| 亚洲成人免费av在线播放| 天堂中文最新版在线下载| 久久久久国产一级毛片高清牌| 国产精品熟女久久久久浪| 午夜免费观看性视频| 一级,二级,三级黄色视频| 正在播放国产对白刺激| 伊人久久大香线蕉亚洲五| 国产精品一区二区在线观看99| 日本vs欧美在线观看视频| av在线app专区| 十八禁网站网址无遮挡| 亚洲伊人久久精品综合| a级毛片黄视频| 亚洲视频免费观看视频| 一进一出抽搐动态| 欧美性长视频在线观看| 宅男免费午夜| 国产成人一区二区三区免费视频网站| 欧美性长视频在线观看| 18禁国产床啪视频网站| 日韩一卡2卡3卡4卡2021年| 亚洲欧美清纯卡通| 美女午夜性视频免费| 电影成人av| 久久久久久久久免费视频了| 天天躁夜夜躁狠狠躁躁| 亚洲精品中文字幕一二三四区 | 手机成人av网站| 国产又色又爽无遮挡免| 嫁个100分男人电影在线观看| 黄色a级毛片大全视频| 国产成人欧美| 最近最新中文字幕大全免费视频| 国产高清国产精品国产三级| 纯流量卡能插随身wifi吗| 午夜福利在线免费观看网站| 女人精品久久久久毛片| 曰老女人黄片| 久久亚洲国产成人精品v| 久9热在线精品视频| 日日爽夜夜爽网站| 啦啦啦免费观看视频1| 一级片免费观看大全| 久久国产精品人妻蜜桃| 中文字幕制服av| 天天添夜夜摸| 亚洲国产中文字幕在线视频| 成年av动漫网址| 性色av乱码一区二区三区2| 日本猛色少妇xxxxx猛交久久| 99热全是精品| 精品久久久久久久毛片微露脸 | 成年av动漫网址| 天堂中文最新版在线下载| 亚洲欧洲精品一区二区精品久久久| 另类精品久久| 国内毛片毛片毛片毛片毛片| 午夜久久久在线观看| 久久久精品区二区三区| 成人亚洲精品一区在线观看| 一级,二级,三级黄色视频| 久久国产精品人妻蜜桃| 久久精品国产a三级三级三级| 满18在线观看网站| 久久久久久久精品精品| 青春草视频在线免费观看| 最新的欧美精品一区二区| 亚洲色图综合在线观看| 国产精品香港三级国产av潘金莲| 18禁黄网站禁片午夜丰满| 99国产精品一区二区三区| 操出白浆在线播放| 久久久久网色| 人妻一区二区av| bbb黄色大片| 亚洲自偷自拍图片 自拍| 91麻豆精品激情在线观看国产 | tube8黄色片| 亚洲精品自拍成人| 一本大道久久a久久精品| kizo精华| 亚洲成国产人片在线观看| 一边摸一边抽搐一进一出视频| a级片在线免费高清观看视频| 欧美日韩一级在线毛片| 欧美日韩av久久| 啦啦啦 在线观看视频| 最新在线观看一区二区三区| 热re99久久精品国产66热6| 午夜福利免费观看在线| 久久ye,这里只有精品| 午夜精品国产一区二区电影| av超薄肉色丝袜交足视频| 青青草视频在线视频观看| 涩涩av久久男人的天堂| 午夜精品久久久久久毛片777| 女人爽到高潮嗷嗷叫在线视频| 在线观看免费视频网站a站| 午夜精品久久久久久毛片777| 欧美国产精品va在线观看不卡| 欧美黑人欧美精品刺激| 老司机福利观看| 不卡一级毛片| 高清黄色对白视频在线免费看| 黄色a级毛片大全视频| 天天躁狠狠躁夜夜躁狠狠躁| 久久精品亚洲av国产电影网| a 毛片基地| 捣出白浆h1v1| 国产精品欧美亚洲77777| 老熟妇仑乱视频hdxx| 狠狠狠狠99中文字幕| 久久久久久久久免费视频了| h视频一区二区三区| 午夜精品久久久久久毛片777| 女性被躁到高潮视频| 男女床上黄色一级片免费看| 欧美国产精品一级二级三级| 涩涩av久久男人的天堂| 在线观看免费午夜福利视频| 色婷婷久久久亚洲欧美| 色94色欧美一区二区| 精品一区在线观看国产| 久久国产精品男人的天堂亚洲| 日韩大码丰满熟妇| 9热在线视频观看99| 女人久久www免费人成看片| 国产亚洲午夜精品一区二区久久| 国产日韩欧美亚洲二区| 欧美亚洲 丝袜 人妻 在线| 91精品伊人久久大香线蕉| 欧美在线黄色| 国产一区二区在线观看av| 国产亚洲av高清不卡| 国产亚洲欧美在线一区二区| 国产人伦9x9x在线观看| 一区福利在线观看| 成年女人毛片免费观看观看9 | 精品国产乱子伦一区二区三区 | 亚洲欧美日韩高清在线视频 | 满18在线观看网站| 精品久久久精品久久久| 日韩熟女老妇一区二区性免费视频| 国产精品欧美亚洲77777| 97在线人人人人妻| 久久久国产成人免费| 老司机影院成人| 成年人免费黄色播放视频| www.精华液| 日韩欧美国产一区二区入口| 最新的欧美精品一区二区| 又黄又粗又硬又大视频| 亚洲av日韩精品久久久久久密| 午夜成年电影在线免费观看| 欧美av亚洲av综合av国产av| 性少妇av在线| 久久久久久人人人人人| 十八禁人妻一区二区| av电影中文网址| 亚洲精品美女久久久久99蜜臀| 可以免费在线观看a视频的电影网站| 精品国产超薄肉色丝袜足j| 国产成人精品无人区| 国精品久久久久久国模美| 精品一品国产午夜福利视频| 国产精品亚洲av一区麻豆| 久久久久久人人人人人| 后天国语完整版免费观看| 在线十欧美十亚洲十日本专区| 国产精品.久久久| 亚洲欧洲精品一区二区精品久久久| 无遮挡黄片免费观看| av福利片在线| 不卡一级毛片| 91国产中文字幕| 不卡av一区二区三区| 亚洲av日韩精品久久久久久密| 久久热在线av| 考比视频在线观看| 老司机福利观看| 99久久人妻综合| 国产黄频视频在线观看| 亚洲精品一卡2卡三卡4卡5卡 | 午夜福利视频精品| 另类亚洲欧美激情| 亚洲综合色网址| 99久久人妻综合| 国产av又大| 少妇 在线观看| 中文字幕人妻丝袜制服| 在线观看免费高清a一片| 91麻豆精品激情在线观看国产 | 国产欧美日韩精品亚洲av| 亚洲国产中文字幕在线视频| 男人添女人高潮全过程视频| 一本—道久久a久久精品蜜桃钙片| 黄色 视频免费看| 日韩大码丰满熟妇| 一级黄色大片毛片| 午夜福利在线观看吧| 亚洲欧美清纯卡通| 欧美激情高清一区二区三区| 亚洲国产欧美网| 成年人免费黄色播放视频| 久久影院123| 国产主播在线观看一区二区| 中文精品一卡2卡3卡4更新| 免费在线观看影片大全网站| 日本a在线网址| 欧美少妇被猛烈插入视频| 丝瓜视频免费看黄片| 欧美xxⅹ黑人| 国产成人精品久久二区二区免费| 国产免费现黄频在线看| 国产1区2区3区精品| 人妻一区二区av| 国产高清videossex| 欧美精品高潮呻吟av久久| 久久久国产成人免费| 又紧又爽又黄一区二区| 色精品久久人妻99蜜桃| 一级a爱视频在线免费观看| 久久人妻熟女aⅴ| 国产成人免费无遮挡视频| 宅男免费午夜| 亚洲情色 制服丝袜| 在线观看免费视频网站a站| av又黄又爽大尺度在线免费看| 精品人妻一区二区三区麻豆| 女人精品久久久久毛片| 波多野结衣av一区二区av| 国产一区二区在线观看av| 亚洲国产成人一精品久久久| 久9热在线精品视频| 亚洲精品在线美女| 欧美精品一区二区免费开放| 国产一级毛片在线| 久久香蕉激情| 18禁观看日本| 一区在线观看完整版| a 毛片基地| 国产成人精品久久二区二区91| 我要看黄色一级片免费的| 久久精品人人爽人人爽视色| 国产成人a∨麻豆精品| 国产欧美日韩精品亚洲av| 久久综合国产亚洲精品| 十八禁网站网址无遮挡| 亚洲精品av麻豆狂野| 老汉色av国产亚洲站长工具| 一级毛片精品| 午夜福利在线免费观看网站| 人妻人人澡人人爽人人| 91成年电影在线观看| 精品国内亚洲2022精品成人 | 一区二区三区激情视频| 超色免费av| 午夜两性在线视频| 香蕉国产在线看| 久久久久久免费高清国产稀缺| 中亚洲国语对白在线视频| av不卡在线播放| 性色av乱码一区二区三区2| 夜夜骑夜夜射夜夜干| 国产一区二区三区在线臀色熟女 | 狠狠狠狠99中文字幕| 免费在线观看视频国产中文字幕亚洲 | 一级片免费观看大全| 欧美日韩av久久| 欧美精品人与动牲交sv欧美| 秋霞在线观看毛片| 欧美亚洲 丝袜 人妻 在线| 国产日韩欧美视频二区| 欧美成狂野欧美在线观看| 亚洲三区欧美一区| 青草久久国产| 在线亚洲精品国产二区图片欧美| 日韩有码中文字幕| 国产深夜福利视频在线观看| 国产精品亚洲av一区麻豆| tube8黄色片| 老司机午夜福利在线观看视频 | 男女床上黄色一级片免费看| 亚洲 欧美一区二区三区| 亚洲成av片中文字幕在线观看| 50天的宝宝边吃奶边哭怎么回事| 捣出白浆h1v1| 亚洲国产欧美网| 国产精品av久久久久免费| 亚洲国产欧美网| 俄罗斯特黄特色一大片| 两性午夜刺激爽爽歪歪视频在线观看 | 一级毛片女人18水好多| 美女高潮到喷水免费观看| 中文字幕高清在线视频| 肉色欧美久久久久久久蜜桃| 色婷婷久久久亚洲欧美| 亚洲一区二区三区欧美精品| 欧美日韩黄片免| 美女主播在线视频| 十八禁高潮呻吟视频| 在线av久久热| 午夜免费观看性视频| 精品一区在线观看国产| 亚洲成av片中文字幕在线观看| 精品国内亚洲2022精品成人 | 国产精品影院久久| 中文字幕色久视频| 午夜福利一区二区在线看| 国产精品一区二区在线不卡| 成年人免费黄色播放视频| 国产成人a∨麻豆精品| 99精国产麻豆久久婷婷| 亚洲一区中文字幕在线| 久久天躁狠狠躁夜夜2o2o| 一级a爱视频在线免费观看| 超色免费av| 国产成人欧美| 一区二区av电影网| 搡老岳熟女国产| 自拍欧美九色日韩亚洲蝌蚪91| 男女无遮挡免费网站观看| svipshipincom国产片| 亚洲成人国产一区在线观看| 脱女人内裤的视频| 久久国产精品大桥未久av| 亚洲自偷自拍图片 自拍| 精品福利永久在线观看| 一级黄色大片毛片| 欧美人与性动交α欧美精品济南到| 成年av动漫网址| 国产区一区二久久| 精品人妻在线不人妻| 成人18禁高潮啪啪吃奶动态图| 精品乱码久久久久久99久播| 少妇猛男粗大的猛烈进出视频| 午夜福利在线免费观看网站| 亚洲一区中文字幕在线| 久久性视频一级片| 午夜精品国产一区二区电影| 99热网站在线观看| 成年动漫av网址| 成人18禁高潮啪啪吃奶动态图| 日本黄色日本黄色录像| 国产亚洲精品一区二区www | 91成年电影在线观看| 一本大道久久a久久精品| 欧美精品人与动牲交sv欧美| 亚洲精品国产区一区二| 一区二区av电影网| 男人添女人高潮全过程视频| 黑人欧美特级aaaaaa片| 欧美精品啪啪一区二区三区 | 大片免费播放器 马上看| 免费观看人在逋| a级毛片黄视频| 十八禁人妻一区二区| 亚洲国产成人一精品久久久| 少妇的丰满在线观看| 日韩视频在线欧美| 久久午夜综合久久蜜桃| 19禁男女啪啪无遮挡网站| 精品卡一卡二卡四卡免费| 日韩欧美免费精品| 亚洲精品成人av观看孕妇| 国产一区二区三区av在线| 久久精品国产a三级三级三级| 自拍欧美九色日韩亚洲蝌蚪91| 蜜桃在线观看..| 一区二区三区四区激情视频| 少妇裸体淫交视频免费看高清 | 满18在线观看网站| 亚洲精品在线美女| 操美女的视频在线观看| 成年动漫av网址| 91老司机精品| 亚洲av电影在线进入| 欧美激情久久久久久爽电影 | 久久99热这里只频精品6学生| 成在线人永久免费视频| 97人妻天天添夜夜摸| 欧美日韩视频精品一区| 国产野战对白在线观看| 国产亚洲午夜精品一区二区久久| 欧美日本中文国产一区发布| av不卡在线播放| 亚洲精品久久久久久婷婷小说| 悠悠久久av| 首页视频小说图片口味搜索| 久久亚洲国产成人精品v| 男男h啪啪无遮挡| 少妇裸体淫交视频免费看高清 | 国产成人系列免费观看| 黑人巨大精品欧美一区二区蜜桃| 久久久久久久久久久久大奶| 两个人看的免费小视频| 亚洲欧美色中文字幕在线| bbb黄色大片| 日韩 亚洲 欧美在线| 亚洲性夜色夜夜综合| 亚洲色图综合在线观看| 另类亚洲欧美激情| 亚洲免费av在线视频| 国产亚洲一区二区精品| 亚洲国产精品一区二区三区在线| 久久国产亚洲av麻豆专区| 在线永久观看黄色视频| 各种免费的搞黄视频| 国产深夜福利视频在线观看| 成人黄色视频免费在线看| 老司机在亚洲福利影院| 丝袜美腿诱惑在线| 91精品国产国语对白视频| 国产不卡av网站在线观看| 十分钟在线观看高清视频www| 老熟女久久久| 亚洲国产av新网站| 12—13女人毛片做爰片一| 色精品久久人妻99蜜桃| 精品卡一卡二卡四卡免费| 日本黄色日本黄色录像| 午夜成年电影在线免费观看| 国产亚洲精品第一综合不卡| 精品国产一区二区久久| 黄色 视频免费看| 亚洲美女黄色视频免费看| 手机成人av网站| 人成视频在线观看免费观看| 两性午夜刺激爽爽歪歪视频在线观看 | 可以免费在线观看a视频的电影网站| 三上悠亚av全集在线观看| 在线观看免费日韩欧美大片| 少妇 在线观看| 久久99一区二区三区| 老司机深夜福利视频在线观看 | 一级毛片女人18水好多| 亚洲精品一二三| 国产高清videossex| 久久久精品免费免费高清| 日韩大片免费观看网站| 国产在线一区二区三区精| 91国产中文字幕| 黄网站色视频无遮挡免费观看| 亚洲av国产av综合av卡| 99热全是精品| 嫁个100分男人电影在线观看| 麻豆乱淫一区二区| 大型av网站在线播放| 欧美 亚洲 国产 日韩一| 亚洲九九香蕉| 中文字幕高清在线视频| 女人久久www免费人成看片| 久久青草综合色| 精品少妇黑人巨大在线播放| 久久综合国产亚洲精品| 亚洲男人天堂网一区| 午夜久久久在线观看| 亚洲第一av免费看| 日本撒尿小便嘘嘘汇集6| 欧美日韩亚洲高清精品| 波多野结衣一区麻豆| 久久人人爽人人片av| 亚洲一码二码三码区别大吗| 99国产精品一区二区三区| 香蕉丝袜av| 一本大道久久a久久精品| 叶爱在线成人免费视频播放| 色播在线永久视频| 精品人妻熟女毛片av久久网站| 青青草视频在线视频观看| 久久亚洲国产成人精品v| 久久午夜综合久久蜜桃| 久久亚洲国产成人精品v| 国产片内射在线| 亚洲成人国产一区在线观看| 久久ye,这里只有精品| 国产欧美日韩一区二区三 | 亚洲国产精品成人久久小说| 黄频高清免费视频| 咕卡用的链子| 国产精品自产拍在线观看55亚洲 | 欧美黑人欧美精品刺激| 91大片在线观看| 欧美黄色片欧美黄色片| 久久ye,这里只有精品| 亚洲视频免费观看视频| 正在播放国产对白刺激| 考比视频在线观看| 国产片内射在线| 人人妻人人爽人人添夜夜欢视频| 国产精品一区二区在线不卡| 欧美激情 高清一区二区三区| 中文字幕最新亚洲高清| 亚洲成人国产一区在线观看| 国产精品一区二区免费欧美 | 夫妻午夜视频| 国产欧美日韩一区二区三 | 午夜两性在线视频| 国产色视频综合| 亚洲精品国产区一区二| 亚洲熟女精品中文字幕| 亚洲精品国产av蜜桃| 精品一区二区三区四区五区乱码| 亚洲少妇的诱惑av| 美女主播在线视频| 老熟女久久久| 亚洲精品乱久久久久久| 天天影视国产精品| 亚洲avbb在线观看| 久久天堂一区二区三区四区| 可以免费在线观看a视频的电影网站| 一个人免费看片子| 18禁裸乳无遮挡动漫免费视频| 国内毛片毛片毛片毛片毛片| 亚洲欧美清纯卡通| 国产精品久久久久久人妻精品电影 | 久久久国产一区二区| 青春草视频在线免费观看|