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

ConceFT在滾動軸承瞬時轉(zhuǎn)頻估計中的應(yīng)用

2019-10-21 04:11:16馬增強阮婉瑩陳明義
振動與沖擊 2019年19期
關(guān)鍵詞:信號方法

馬增強, 阮婉瑩, 陳明義

(1. 省部共建交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國家重點實驗室, 石家莊 050043;2. 石家莊鐵道大學(xué) 電氣與電子工程學(xué)院, 石家莊 050043)

基于瞬時轉(zhuǎn)頻估計的階比分析方法無需安裝鍵相裝置,節(jié)省了空間及成本,是變轉(zhuǎn)速階比分析[1]的研究熱點,而瞬時轉(zhuǎn)頻的準確估計是其成功的關(guān)鍵和前提。現(xiàn)有瞬時轉(zhuǎn)頻估計方法很多元,第一大類為基于相位解調(diào)法,例如:Combet等[2]采用短時尺度變換法實現(xiàn)回轉(zhuǎn)軸瞬時轉(zhuǎn)頻的提取;Coats等[3]利用多級迭代法實現(xiàn)相位解調(diào),提取瞬時轉(zhuǎn)頻信息。該類方法在轉(zhuǎn)速大范圍變化時,會在相鄰階次諧波之間產(chǎn)生交叉項,使得該方法應(yīng)用受限。第二大類為基于時頻分析法,例如:郭瑜等[4]將短時傅里葉變換(Short-time Fourier Transform, STFT)與峰值搜索結(jié)合,對STFT時頻譜進行峰值搜索,進而對變轉(zhuǎn)速電機瞬時轉(zhuǎn)頻進行估計,取得較好的測試效果,但由于STFT的固有缺陷使得其抗噪性較差;趙曉平等[5]提出STFT結(jié)合Viterbi算法估計變轉(zhuǎn)速振動信號瞬時轉(zhuǎn)頻,成功應(yīng)用于渦流離心機升速狀態(tài)的轉(zhuǎn)頻估計,但Viterbi算法復(fù)雜度高,計算效率低;Urbanek等[6]對比幅值和相位的時頻信息,尋找對應(yīng)關(guān)系成功提取齒輪箱的瞬時轉(zhuǎn)頻。此外,還有很多其他轉(zhuǎn)頻估計算法,如:羅潔思等[7-8]、程衛(wèi)東等[9]將多尺度線調(diào)頻路徑追蹤應(yīng)用于旋轉(zhuǎn)機械瞬時轉(zhuǎn)頻提取,但該算法最大缺點是計算速度太慢,難以處理大量實驗數(shù)據(jù),很難應(yīng)用于實際;Wang等[10]利用經(jīng)驗?zāi)B(tài)分解處理變轉(zhuǎn)速信號,進而估計瞬時轉(zhuǎn)頻,但該方法存在模態(tài)混疊、端點效應(yīng)等問題,處理復(fù)雜多分量信號時會產(chǎn)生較大誤差。

綜上所述,基于時頻分析的瞬時轉(zhuǎn)頻估計方法原理簡單,不受轉(zhuǎn)速變化影響,適用范圍廣,更具實際應(yīng)用價值,是現(xiàn)階段研究瞬時轉(zhuǎn)頻估計的熱門方法。該類算法成功的關(guān)鍵在于所用時頻分析方法是否具備較高的抗噪性和時頻分辨率,傳統(tǒng)的時頻分析方法(STFT、小波變換和Wigner-Ville分布)在時頻分辨率、交叉項及抗噪性上均存在不同程度的問題,不便于直接應(yīng)用于信號分析,需要不同的優(yōu)化算法進行改進。現(xiàn)有的基于時頻分析的瞬時轉(zhuǎn)頻估計方法中,通常所用的都是上述傳統(tǒng)時頻分析方法,很難取得滿意的效果,這大大限制了此類方法的應(yīng)用。筆者由此入手,將一種新型的時頻分析方法——時頻聚集(Concentration of Frequency and Time,ConceFT)應(yīng)用其中。

ConceFT是Daubechies等[11]在同步壓縮變換[12](Synchrosqueezing Wavelet Transform, SST)基礎(chǔ)上,于2016年提出的全新的時頻分析方法,該方法集合了多正交窗和SST的優(yōu)點,不但保證了高時頻分辨率,而且在抗噪性方面有極大改善,恰好彌補現(xiàn)有時頻分析方法的不足,其能夠勝任多分量復(fù)雜信號的處理,有著廣闊的發(fā)展空間和應(yīng)用前景。這一嶄新的時頻分析方法自被提出以來還尚未被廣泛應(yīng)用,本文首次提出將其聯(lián)合峰值搜索算法用于滾動軸承的瞬時轉(zhuǎn)頻估計,著重研究算法的抗噪性能及瞬時轉(zhuǎn)頻估計精度,結(jié)果驗證了所提方法的有效性,ConceFT聯(lián)合峰值搜索能夠在較低信噪比下得到瞬時轉(zhuǎn)頻的精確估計。

1 理論基礎(chǔ)

1.1 時頻聚集方法原理

1.1.1 同步壓縮變換

同步壓縮變換是從連續(xù)小波變換出發(fā)的,首先對信號進行連續(xù)小波變換,再進行時頻重排,對小波系數(shù)進行壓縮變換,使得時頻分布更清晰、時頻分辨率更高。

Daubechies以一個幅值恒定的諧波信號s(t)=Acos(ωt) 為例進行闡述。

首先求出信號s(t)的連續(xù)小波變換為

(1)

然后計算信號s(t)的瞬時頻率

(2)

式中:?bWs(a,b)為Ws(a,b)對b的一階偏導(dǎo)。

通過估計出的瞬時頻率ωs(a,b)可以將原來的時間-尺度平面轉(zhuǎn)換到時間-頻率平面:(a,b)→[ωs(a,b),b]。小波系數(shù)變?yōu)閃s(ωs(a,b),b) ,通過壓縮時間-頻率平面內(nèi)的小波系數(shù),將其壓縮至中心頻率ωl的一個鄰域內(nèi):[ωl-1/2Δω,ωl+1/2Δω] ,從而得到同步壓縮系數(shù)Ts(ωl,b) 。考慮計算機在計算過程中,a,b,ω均需離散化,設(shè)ai-ai-1=(Δa)i,ωi-ωi-1=Δω,則同步壓縮系數(shù)可表示為

(3)

同步壓縮變換相比于重排算法的優(yōu)勢在于其支持信號重構(gòu),根據(jù)同步壓縮系數(shù)Ts(ωl,b)可以重構(gòu)信號s(b),表達式如下

(4)

對于離散參數(shù),s(b)近似為

(5)

SST將時頻分布沿尺度/頻率方向進行壓縮,從而有效提高了時頻分辨率。但是,由于雜散噪聲的存在使真實信號的時頻分布變得模糊,時頻分辨率大大降低,為此,引入多正交窗來降低噪聲的干擾。

1.1.2 多正交窗同步壓縮變換

在同步壓縮變換的時頻分布中,對于不同母小波函數(shù),時頻分布不受影響,然而噪聲分布卻隨母小波函數(shù)的不同而不同。因為小波變換本質(zhì)上就是信號與母小波函數(shù)的卷積,因此基于不同母小波函數(shù)的小波變換相互獨立,恒等分布,這一特性引出了多正交窗同步壓縮變換的想法。

對于信號x(t),給定I個標準正交母小波φi(t),i=1,…,I,根據(jù)式(2)、式(3)計算出每個母小波對應(yīng)的瞬時頻率和同步壓縮變換,多正交窗同步壓縮變換即定義為各母小波同步壓縮變換的平均

(6)

對大量標準正交母小波進行平均能夠抵消噪聲引起的時頻模糊。由此認為標準正交母小波數(shù)量越多抑制噪聲干擾能力越強,但是隨著標準正交母小波的增多,時頻分布中模糊區(qū)域增多,因此,標準正交小波數(shù)目需要設(shè)置一個平衡,這在一定程度上限制了其使用。

1.1.3 ConceFT時頻聚集方法

為了克服多正交窗同步壓縮變換上述缺陷,Daubechies根據(jù)同步壓縮變換的非線性特點,拓展了多正交窗方法提出ConceFT方法。

步驟1選擇I個標準正交母小波φi(t),i=1,…,I,使其滿足良好的時頻聚集性。

步驟2選擇N個單位隨機向量rn,n=1,…,N。

步驟5對隨機向量rn做平均,得ConceFT的時頻表達:

(7)

實際應(yīng)用中,為了提高時頻分辨率,標準正交母小波可選則2個,而為了抑制噪聲干擾,權(quán)重向量N可根據(jù)需要盡可能取大[13]。

1.2 峰值搜索原理

對于一個高精度、高時頻聚集性的時頻分析方法,應(yīng)用峰值搜索算法即可實現(xiàn)從時頻圖中準確提取瞬時頻率,且峰值搜索原理簡單,效率較高,故本文利用峰值搜索法實現(xiàn)轉(zhuǎn)頻曲線的提取。下面介紹峰值搜索法是如何從時頻圖中進行瞬時頻率估計的步驟:

步驟1時頻分析得時頻圖。

步驟2選取搜索起始點。在時頻圖中被跟蹤分量峰值突出的區(qū)域內(nèi)選擇一點作為起始點。選定起始點之后,對時頻圖進行峰值搜索,按照下式進行

(8)

式中:ni=n1±1,n1±2,…,ni∈(0,M-1);ki∈(0,N-1);M時頻網(wǎng)格中時間線數(shù);N為時頻網(wǎng)格中頻率線數(shù);IFE為峰值搜索函數(shù);argmax為目標函數(shù)取最大值時所取參數(shù);SPEC為對應(yīng)的時頻圖;(n1,k1) 為以(n1,k0) 為起始點進行峰值搜索所得的第一個瞬時頻率坐標;p為峰值搜索的范圍;(ni,ki) 為經(jīng)過峰值搜索所得的各時刻對應(yīng)的瞬時頻率坐標。

步驟3計算各點瞬時轉(zhuǎn)頻。按照下式進行

(9)

式中:fq(ni)為峰值搜索所得各點瞬時頻率;q為瞬時頻率對應(yīng)的故障特征階次;ni為相應(yīng)的時間點。

步驟4轉(zhuǎn)頻曲線擬合。對上述所得離散瞬時轉(zhuǎn)頻進行最小二乘擬合。根據(jù)各點瞬時轉(zhuǎn)頻變化趨勢,選擇多項式次數(shù)。通常情況下,轉(zhuǎn)速不會發(fā)生突變,可選擇低次多項式擬合,以二次為例,擬合公式如下

(10)

平方誤差如下

(11)

2 ConceFT估計滾動軸承瞬時轉(zhuǎn)頻實現(xiàn)過程

實際滾動軸承振動信號常伴有強烈噪聲干擾,鑒于ConceFT具有強抗噪性和高時頻分辨率的優(yōu)點,提出了將其與峰值搜索結(jié)合用于滾動軸承的瞬時轉(zhuǎn)頻估計,以解決現(xiàn)有方法的抗噪性差及估計精度不足的問題。對于實際的滾動軸承振動信號,為了提高其瞬時轉(zhuǎn)頻估計精度,在進行ConceFT之前需要對振動信號進行預(yù)處理[14-16],步驟如下:

步驟1低通濾波,避免頻率混疊。低通濾波的截止頻率選為降采樣后頻率的1/2。

步驟2降采樣,凸出時頻圖中的轉(zhuǎn)頻分量,提高估計精度。按照式(12)設(shè)置采樣倍數(shù)。

(12)

式中:d為降采樣倍數(shù);int為取整函數(shù);fs為實際采樣頻率;fRmax為最大轉(zhuǎn)頻;q為搜索分量的階次。

步驟3去趨勢項,提高振動信號的信噪比。

步驟4Hilbert包絡(luò)解調(diào),得轉(zhuǎn)頻相關(guān)信息。

對振動信號預(yù)處理后再進行ConceFT分析,進一步峰值搜索即可提取瞬時轉(zhuǎn)頻曲線。本文方法估計滾動軸承瞬時轉(zhuǎn)頻的總體流程圖,如圖1所示。

3 仿真信號驗證

變轉(zhuǎn)速下滾動軸承振動信號為復(fù)雜的多分量信號,在不同工況會有不同的頻率調(diào)制現(xiàn)象、不同類型的轉(zhuǎn)頻曲線,且伴隨的強烈背景噪聲會極大地影響轉(zhuǎn)頻曲線提取精度。本文分別設(shè)計一個線調(diào)頻多分量信號和一個正弦調(diào)頻多分量信號來模擬兩種工況下滾動軸承的振動信號,以更好地驗證ConceFT方法在估計瞬時轉(zhuǎn)頻曲線問題上的優(yōu)勢。文獻[13]已證明母小波的選擇對ConceFT分析效果幾乎沒有影響,故在仿真信號分析中,采用2個標準正交Morse小波,權(quán)重向量N取5。

圖1 ConceFT估計滾動軸承瞬時轉(zhuǎn)頻流程圖

3.1 線調(diào)頻信號驗證

模擬線調(diào)頻的多分量振動信號,設(shè)信號瞬時角頻率為

ω(t)=2.5πt+30π

(13)

則瞬時轉(zhuǎn)頻為

f(t)=1.25t+15

(14)

多分量線調(diào)頻信號為

(15)

式中:η(t)為高斯白噪聲;0≤t≤20 s。

該信號信噪比(Signal Noise Ratio,SNR)取為-10 dB,采樣頻率為100 Hz。信號波形如圖2,由ConceFT所得時頻圖如圖3所示。由圖3可知,信號的一倍頻、0.66倍頻和0.5倍頻的頻率曲線,基本不受強噪聲的影響。為了凸顯ConceFT方法的優(yōu)勢,我們與經(jīng)典的瞬時轉(zhuǎn)頻估計方法——STFT結(jié)合峰值搜索進行對比,圖4為信號的STFT時頻圖。由圖4可知,其受噪聲干擾嚴重,各頻率成分被淹沒在強噪聲中,難以識別。利用峰值搜索法分別對ConceFT和STFT的時頻分布圖進行瞬時轉(zhuǎn)頻提取,結(jié)果如圖5所示。可見本文方法提取的瞬時轉(zhuǎn)頻曲線與真實轉(zhuǎn)頻曲線基本吻合,而基于STFT的轉(zhuǎn)頻估計結(jié)果偏離真實值太大。

為了定量說明兩種方法的轉(zhuǎn)頻估計精度,利用式(16)計算估計誤差的百分比值。通過計算可得,本文方法估計誤差為2.56%,STFT結(jié)合峰值索搜的估計誤差為52.68%。

(16)

圖3 線調(diào)頻信號波形圖

圖3 線調(diào)頻信號ConceFT時頻圖

圖4 線調(diào)頻信號STFT時頻圖

圖5 線調(diào)頻信號ConceFT與STFT轉(zhuǎn)頻估計對比圖

3.2 正弦調(diào)頻信號驗證

模擬正弦調(diào)頻的多分量振動信號,設(shè)信號角頻率為

ω(t)=2π(7.5-2.5cos(t))

(17)

則瞬時轉(zhuǎn)頻為

f(t)=7.5-2.5cos(t)

(18)

多分量線調(diào)頻信號為

(19)

式中:η(t)為高斯白噪聲;0≤t≤20 s。

信號信噪比取為-8 dB,采樣頻率為100 Hz,信號波形如圖6所示。ConceFT時頻圖如圖7所示,圖7中三個頻率分量清晰可見,而STFT的時頻結(jié)果如圖8所示。由圖8可知,在強烈噪聲干擾下,時頻能量發(fā)散嚴重,不能識別出真實信號分量。分別對兩個時頻圖進行峰值索搜,提取出的瞬時轉(zhuǎn)頻曲線如圖9所示。由圖9可知,STFT峰值搜索所得曲線較理論曲線波動較大,原因是信號信噪比低,部分噪聲能量高于信號能量,導(dǎo)致峰值索搜時出現(xiàn)誤搜索,導(dǎo)致最后擬合曲線出現(xiàn)較大偏差。由于本文方法有較強的抗噪性與高時頻分辨率,能夠弱化噪聲影響,使頻率成分能量集中,可以對其時頻圖進行精確的峰值搜索,所得擬合的曲線與理論曲線基本吻合。經(jīng)過計算,本文方法估計誤差為1.26%,STFT結(jié)合峰值索搜的估計誤差為42.85%。

圖6 正弦調(diào)頻信號波形圖

圖7 正弦調(diào)頻信號ConceFT時頻圖

圖8 正弦調(diào)頻信號STFT時頻圖

圖9 正弦調(diào)頻信號ConceFT與STFT轉(zhuǎn)頻估計對比圖

Fig.9 Contrast figure between ConceFT and STFT for rotational frequency estimation of sinusoidal frequency modulation signal

3.3 抗噪性及估計精度分析討論

在此,我們進一步分析討論ConceFT方法的抗噪性及用于瞬時轉(zhuǎn)頻估計的精度問題。為了更具說服力,筆者采用大量數(shù)據(jù)進行測試。分別對上述線調(diào)頻信號和正弦調(diào)頻信號添加白噪聲,使信號信噪比從-20~10 dB之間變化,對各狀態(tài)下的信號進行基于ConceFT的轉(zhuǎn)頻估計工作,計算不同信噪比下的轉(zhuǎn)頻估計誤差。

選4組典型信噪比下的信號,做ConceFT時頻圖,圖10對應(yīng)線調(diào)頻信號。圖11對應(yīng)正弦調(diào)頻信號。從兩幅圖中可知,在信噪比很低時,時頻分辨率仍然很高,仍能清晰顯示信號的瞬時頻率曲線。

為了定量說明該方法在低信噪比下仍能保持高時頻分辨率,采用Rényi熵作為評價指標,Rényi熵值越小表示時頻分辨率越高,其定義式如式(20)所示

(a) 信噪比為10 dB

(b) 信噪比為0

(c) 信噪比為-10 dB

(d) 信噪比為-20 dB

(a) 信噪比為10 dB

(b) 信噪比為0

(c) 信噪比為-10 dB

(d) 信噪比為-20 dB

(20)

式中,R為Rényi熵,q≥0,q≠1,(p1,p2,…,pn)為任意離散變量的概率分布。

就圖10和圖11所用的4種典型信噪比下的兩種信號分別求ConceFT和STFT的Rényi熵,如表1和表2所示。從表1和表2可知,兩個表中ConceFT方法的Rényi熵遠小于STFT,且隨著信噪比的降低,ConceFT的Rényi熵變化并不大,說明其時頻分辨率受信噪比影響不大,說明該方法具有強抗噪性,低信噪比時仍具有高時頻分辨率,這是精確提取瞬時轉(zhuǎn)頻曲線的前提。

表1 線調(diào)頻信號不同信噪比下兩種方法的Rényi熵

表2 正弦調(diào)頻信號不同信噪比下兩種方法的Rényi熵

為了驗證本文方法的瞬時轉(zhuǎn)頻估計精度,再取16組不同信噪比的信號進行分析,計算轉(zhuǎn)頻估計誤差,結(jié)果如圖12所示,從圖12可知,信噪比在-20 dB以上時,估計誤差都保持在3%以內(nèi),足夠滿足精度要求。

圖12 ConceFT對于兩種信號轉(zhuǎn)頻估計的信噪比-估計誤差曲線圖

4 實測信號驗證

利用實測滾動軸承變轉(zhuǎn)速故障振動信號來進一步驗證本文方法的實用性。模擬故障實驗臺如圖13所示。從圖13可知,“1”為CA-YD-188型加速度傳感器,“2”為ICP激光轉(zhuǎn)速計,試驗軸承為外圈有輕微裂紋故障,型號為NU205EM。振動信號采樣頻率為25 600 Hz,激光轉(zhuǎn)速計采樣頻率為1 kHz,據(jù)此利用五點公式法擬合轉(zhuǎn)頻曲線,作為理論轉(zhuǎn)頻曲線。本文利用兩組數(shù)據(jù)進行充分驗證,分別取轉(zhuǎn)速上升階段和轉(zhuǎn)速復(fù)雜波動變化階段的軸承故障振動信號。ConceFT算法中選2個標準正交Morse小波,N取為10。

4.1 升速工況下瞬時轉(zhuǎn)頻估計精度驗證

振動信號波形如圖14所示。為強背景噪聲干擾下的升速工況,轉(zhuǎn)頻由11.4 Hz上升至24.6 Hz。首先對信號進行預(yù)處理,低通濾波截止頻率設(shè)為1 500 Hz,降采樣倍數(shù)設(shè)為5,所得降采樣后的采樣頻率為5 120 Hz,對降采樣后的信號去趨勢項,再取Hilbert包絡(luò),預(yù)處理結(jié)束。對預(yù)處理后的信號進行ConceFT得時頻圖如圖15所示。從圖15可知,一倍轉(zhuǎn)頻能量最高,最突出,將其作為峰值索搜的目標得轉(zhuǎn)頻曲線,經(jīng)過計算轉(zhuǎn)頻估計誤差為0.59%,足以滿足精度要求。

圖13 QPZZ-Ⅱ旋轉(zhuǎn)機械故障試驗臺

圖14 升速工況振動信號波形

圖15 升速工況振動信號ConceFT時頻圖

4.2 轉(zhuǎn)速劇烈波動工況下瞬時轉(zhuǎn)頻估計精度驗證

振動信號波形如圖17所示。從圖17可知,轉(zhuǎn)速升降復(fù)雜變化的工況,可見信號中夾雜著大量噪聲成分,該信號轉(zhuǎn)頻在6~20 Hz之間升降往復(fù)變化,轉(zhuǎn)頻曲線的提取難度大大增加。首先仍然先對信號進行預(yù)處理,同“4.1”所述,對預(yù)處理后的信號進行ConceFT得時頻圖,如圖18所示。由圖18能清晰識別一倍轉(zhuǎn)頻分量,將其作為峰值索搜的目標,得轉(zhuǎn)頻曲線如圖19中虛線所示。圖19中實線為理論轉(zhuǎn)頻,可見二者相似度極高,基本重合,經(jīng)過計算該轉(zhuǎn)頻估計誤差僅為0.51%。可見,在復(fù)雜工況下,該方法仍能保持非常高的估計精度。

圖16 升速工況ConceFT轉(zhuǎn)頻曲線估計結(jié)果與理論轉(zhuǎn)頻對比圖

圖17 轉(zhuǎn)速波動工況振動信號波形

圖18 轉(zhuǎn)速波動工況振動信號ConceFT時頻圖

5 結(jié) 論

(1) 分析表明ConceFT方法是一種具有強抗噪性、高分辨率的時頻分析方法,將其與峰值索搜結(jié)合用于滾動軸承瞬時轉(zhuǎn)頻估計,通過仿真信號與兩種復(fù)雜工況下的變轉(zhuǎn)速故障信號分析,結(jié)果驗證了該方法在復(fù)雜噪聲干擾下仍能準確提取信號瞬時轉(zhuǎn)頻,準確度高達99%以上,滿足高精度的要求,可以用于無轉(zhuǎn)速計情況下的轉(zhuǎn)頻估計與轉(zhuǎn)速曲線提取工作,能夠解決現(xiàn)有方法抗噪性差及估計精度不足的問題,是一種有實際應(yīng)用價值的轉(zhuǎn)頻估計方法。

圖19 轉(zhuǎn)速波動工況ConceFT轉(zhuǎn)頻曲線估計結(jié)果與理論轉(zhuǎn)頻對比圖

(2) ConceFT作為一種新興的時頻分析方法,由于其具有強抗噪性及較高的時頻分辨率而有廣闊的應(yīng)用空間。如,基于時頻分析的變轉(zhuǎn)速軸承故障診斷方法較階比分析簡單直接,但目前的時頻分析方法在抗噪性及時頻分辨率上通常不能滿足要求,利用ConceFT可以很好的解決此問題;齒輪等其他機械設(shè)備的故障診斷問題同樣可以借助ConceFT來解決;其他領(lǐng)域,如地震信號、油氣檢測等凡需時頻分析的領(lǐng)域均可嘗試用ConceFT來處理。

猜你喜歡
信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
學(xué)習(xí)方法
孩子停止長個的信號
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 97精品久久久大香线焦| 在线欧美日韩| 精品国产成人av免费| 色亚洲激情综合精品无码视频| 毛片在线播放网址| 国产乱人伦偷精品视频AAA| 国产网站在线看| 欧美a级完整在线观看| 亚洲天堂免费| 国产99精品久久| 91免费国产在线观看尤物| 午夜不卡福利| 最新国语自产精品视频在| 青青操国产| 欧美一级夜夜爽www| 成人午夜亚洲影视在线观看| 欧美 亚洲 日韩 国产| 国产成人精品午夜视频'| 成人年鲁鲁在线观看视频| 精品视频第一页| 国产sm重味一区二区三区| 国产精品私拍99pans大尺度| 天天爽免费视频| 色成人亚洲| www欧美在线观看| 亚洲丝袜中文字幕| 欧洲精品视频在线观看| 88国产经典欧美一区二区三区| 国产精品视频公开费视频| 99视频在线免费观看| 国产av色站网站| 久草中文网| 婷婷亚洲天堂| 91欧美亚洲国产五月天| 黄色一级视频欧美| 亚洲精品无码高潮喷水A| 国产成人一区免费观看| 试看120秒男女啪啪免费| 激情视频综合网| 精品视频免费在线| 国产草草影院18成年视频| a毛片免费在线观看| 成人字幕网视频在线观看| 成人一级免费视频| 国产精品亚洲专区一区| 色综合天天综合| 超碰免费91| 国产资源免费观看| 国产福利一区二区在线观看| 99激情网| 国产激情在线视频| 欧美一级在线看| 丰满的熟女一区二区三区l| 婷婷综合在线观看丁香| 国产另类视频| 91欧美在线| 亚洲人成网站18禁动漫无码| jizz在线免费播放| 性欧美久久| 亚洲性影院| 日韩区欧美区| 午夜成人在线视频| 亚洲国产成人久久77| 国产噜噜在线视频观看| 亚洲精品福利网站| 无码免费的亚洲视频| 久久久久人妻一区精品色奶水| 国产精品视频免费网站| 久草青青在线视频| 青青久久91| 秘书高跟黑色丝袜国产91在线| 久久99热这里只有精品免费看| 日日噜噜夜夜狠狠视频| 亚洲乱伦视频| 狠狠色成人综合首页| 91在线视频福利| 亚洲第一天堂无码专区| 亚洲成人黄色网址| 久久精品最新免费国产成人| 少妇精品久久久一区二区三区| 日韩小视频网站hq| 亚洲欧美一区二区三区图片|