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

基于自動搜峰瞬時頻率估計的自適應多階比分析技術研究

2016-01-15 05:24:18王棟,丁雪娟
振動與沖擊 2015年18期
關鍵詞:故障診斷

第一作者王棟男,碩士,中級工程師,1982年生

通信作者丁雪娟女,碩士,1986年生

基于自動搜峰瞬時頻率估計的自適應多階比分析技術研究

王棟,丁雪娟

(中能電力科技開發有限公司,北京100034)

摘要:提出一種利用自動搜峰的瞬時頻率估計技術來實現旋轉機械自適應多階比分析(Adaptive Multiple Order Tracking, AMOT)的新方法。首先,通過時頻分析得到振動信號的時頻分布,根據頻率峰值坐標自動選取搜峰起始點,自適應搜索出不同階比分量的時頻峰值。其次,利用最小二乘法將不同頻率分量進行擬合實現瞬時頻率估計,然后根據參考分量計算出重采樣的鑒相時標對原始信號進行重采樣,最后通過FFT變換實現階比分析。該方法通過瞬時頻率估計能夠自動識別出所有階比分量,實現優中選優,避免了傳統算法中人為直觀選取一個分量進行遮掩濾波提取分量的方法,減少了人為選取分量及起始點造成的誤差,具有自適應性。并且無需同步采集轉速信號,大大簡化了應用條件,同時減少了人為因素,提高了分析精度,為旋轉機械故障診斷提供了新方法。仿真實驗和應用實例驗證了方法的有效性。

關鍵詞:故障診斷;自適應多階比分析;瞬時頻率;角域重采樣

基金項目:龍源電力集團股份有限公司科技項目

收稿日期:2014-05-21修改稿收到日期:2014-08-19

中圖分類號:TN911.6; TH115文獻標志碼:A

Tacholess adaptive multi-order tracking technology based on instantaneous frequency estimation with automatic peak search algorithm and its application

WANGDong,DINGXue-juan(Zhong Neng Power-Tech. Development Co. Ltd., Beijing 100034, China)

Abstract:A new method of adaptive multi-order tracking (AMOT) in rotating machinery was proposed based on the instantaneous frequency estimation with an adaptive peak search algorithm. The time-frequency distribution was obtained by using a certain time-frequency analysis method. According to the frequency peak coordinate, the starting point of peak search was selected automatically and time-frequency peaks of different orders were searched adaptively. The different frequency components were fitted to achieve the instantaneous frequency estimation by using the least square polynomial, and then the resampling time was calculated according to the reference component to resample the original signal. Finally multi-order tracking spectrums were obtained by FFT transform. By virtue of the method of instantaneous frequency estimation, the proposed method can automatically identify all-order components and choose the best one, avoiding the human error owing to choosing reference component and starting point subjectively in traditional algorithms. The method can identify the components of different orders without sampling synchronously the speed signal, greatly simplify the application conditions and improve the accuracy of analysis. It provides a new method for fault diagnosis of rotating machinery. The method was proved to be effective by simulation examples and practical applications.

Key words:fault diagnosis; adaptive multiple-order tracking; instantaneous frequency; angle domain resampling

旋轉機械是工業生產中最重要的一類機械設備,廣泛應用于電力、能源、礦山、冶金、國防等各個領域。為了確保旋轉機械的安全平穩運行,狀態監測與故障診斷技術得到了應用和發展。目前,基于振動信號的狀態監測與故障診斷技術是最主要的方式。實際的振動信號一般都是非平穩信號,另外,旋轉機械在啟、停機和升降速的過程中表現為間歇性和非平穩性,包含豐富的故障信息。對于非平穩振動信號的處理,常規的信號處理方法往往失效,無法實現準確的故障診斷。

階比分析是20世紀90年代發展起來的一種針對非平穩旋轉機械振動信號的變頻采樣技術[1]。由于階比分析在處理非平穩振動信號,尤其是處理轉速變化時的振動信號方面的優勢,對其研究得到了廣泛的關注。目前,階比分析技術主要分為兩大類:一類是采用鑒相裝置和轉速計的硬件式階比跟蹤技術[2-3],另一類是計算階比跟蹤技術[4-5]。對于硬件階比技術,受到附加硬件及現場安裝條件的制約,應用受到了限制。雖然計算階比跟蹤技術降低了對硬件的依賴,但仍然需要提供轉速脈沖。

近年來提出的基于瞬時頻率估計的階比分析技術,根據頻率和轉速曲線的線性對應關系,完全靠軟件來獲得轉速信息。通過時頻分析方法獲得信號的轉速曲線,進而求得信號重采樣的鑒相時刻,擺脫了階比分析過程中對硬件的依賴,簡化了階比采樣的實現要求[6-7]。文獻[8]提出STFT-Viterbi擬合法進行瞬時頻率估計,取得了一定效果。該方法需將頻率范圍分段,并給定參考軸轉速。文獻[9]提出改進的峰值搜索法,通過小波時頻分布圖進行峰值搜索,利用導數差值搜索峰值,避免虛假峰值,提高了瞬時頻率的估計精度。但無法實現多階比的自適應搜索。

本文提出一種基于自動搜峰瞬時頻率估計的自適應多階比分析(Adaptive Multiple Order Tracking,AMOT)新方法。通過自適應搜峰算法在振動信號的時頻分布上得到多階比的時頻峰值頻段,避免了人為因素的影響,能夠自適應的選取起始點并根據判據條件自適應的提取出多個階比分量。經最小二乘多項式擬合實現參考分量瞬時頻率曲線擬合,計算信號重采樣的鑒相時標,對原始信號進行角度重采樣,最后進行頻譜變換實現多階比分析。仿真信號和實例信號結果均表明,該方法能夠自動識別出不同階比分量,具有自適應性,同時無需同步采集轉速信號,簡化了應用條件,可操作性強,為無轉速計的階比分析技術的在線監測應用提供了新方法。

1基于時頻分析的自適應搜峰瞬時頻率估計

階比分析是將時域非平穩振動信號通過等角度重采樣轉換為角域信號獲得階比譜的一種信號分析方法。階比能夠很好地反映與轉速相關的振動特征。對于非平穩振動信號,利用時頻分析技術可以進行瞬時頻率估計,進而得到鑒相信號,實現基于時頻分析獲取瞬時頻率的無轉速計階比跟蹤技術。

1.1瞬時頻率

瞬時頻率的概念具有一定的爭議。原因有兩個:一個是人們深受Fourier變換的影響,認為頻率只能是用諧波函數在整個信號內定義。這樣,瞬時頻率的定義就和正余弦函數聯系起來,當短于一個正余弦波形周期的頻率是沒有意義。另一個是瞬時頻率沒有統一的定義方法。直到應用Hilbert變換來構造解析信號時,解決了瞬時頻率的統一定義問題。

給定信號x(t),其Hilbert變換可以定義為函數x(t)與1/πt的卷積

(1)

式中:p.v.表示取積分的主值

(2)

式中:

(3)

那么a(t)是z(t)的瞬時幅值,反映了信號能量隨時間的變化情況。θ(t)是z(t)的瞬時相位。

理論上,可以有無數多種方法定義z(t)的虛部。而Hilbert變換用一種特別的方法來定義,通過這種方法來定義的虛部,z(t)能滿足作為x(t)的解析信號的條件。由式(1)可以看出,Hilbert變換實際上是信號x(t)和1/t的卷積。因此,Hilbert變換可以表示x(t)的局部特性。式(2)的幅值和相位隨時間變化,很好地表示了x(t)的局部特征。通過Hilbert變換,瞬時頻率定義如下:

(4)

式(4)表明,瞬時頻率是t的函數,任意時刻的瞬時頻率只有一個值,只表示一個單分量頻率。

1.2自適應峰值搜索算法

基于瞬時頻率估計的階比分析中,需利用原始信號估計參考軸轉速。文獻[8-9]利用峰值搜索法獲得轉速,并計算原始信號的時頻分布圖,對時頻分布圖一階轉速附近進行峰值搜索,將峰值時刻對應的頻率視為該時刻轉子的瞬時頻率。但沒有實現自適應多階比分析。本文提出一種自適應多階比分析算法,采用局部自動搜峰法,從最低頻率開始自動搜索第一個階比分量的時頻分布峰值帶寬,繼續向上搜索,按照從低到高的順序,依次搜索出所有階比分量,具體算法為:在時頻面上,設時頻面的時間長度為N,頻率長度為M,具體步驟描述如下:

(2)對nNB時刻(一般取中間時刻)這一列頻率進行自動分段,即從最低頻率開始,當出現連續f(nNB,mi)≥f0時,搜索出Q個頻段,即,

Pj(nNB,hkj)=f(nNB,mi)

(5)

式中:Pj(nNB,hkj)為第j個階比分量所包含的頻率數組;hkj=1,…,kj;kj為階比分量瞬時頻率寬度,j=1,…,Q;Q為階比個數;i=1,…,l1,l1+1,…,l1+k1-1,…,l2,l2+1,…,l2+k2-1,…,lQ,lQ+1,…,lQ+kQ-1。

Pj(nNB,hkj)=f(nNB,mi)具體表達式為:

(6)

(3)依次對每一個階比帶寬頻率Pj(nNB,hkj),j=1,…,Q,分別進行峰值搜索,得到相應的最大頻率值,具體算法為:

(7)

其中:arg max表示取最大值;Pj(nNB,mMj)為最大頻率峰值。

(4)從nNB時刻,改變時間ng,g=1,NB-1,NB+1,…,N。分別向兩側重復“(2)”逐列進行自動搜峰,得到所有時間列的階比頻率帶寬。Pj(nNg,hkj),為第g個時間列中,第j個階比分量所包含的頻率數組;hkj=1,…,kj;kj為階比分量瞬時頻率寬度,g=1,NB-1,NB+1,…,N;j=1,…,Q。

(5)對每一個階比帶寬頻率Pj(nNg,hkj),g=1,NB-1,NB+1,…,N,j=1,…,Q,分布進行峰值搜索,具體算法為:

(8)

式中:Pj(nk,mMj)為最大頻率峰值,k=1,NB-1,NB+1,…,N;j=1,…,Q。

(6)最后得到Q個階比分量峰值Pj(ns,mMj),s=1,NB-1,NB,NB+1,…,N;j=1,…,Q。

自適應多階比搜峰示意圖見圖1。

圖1 自適應多階比搜峰示意圖 Fig.1 The sketch map of searching peaks adaptively for AMOT

2基于瞬時頻率估計的自適應多階比分析(AMOT)算法

基于瞬時頻率估計的自適應多階比分析(AMOT)算法流程圖見圖2。

AMOT算法具體步驟為:

(1)采用振動加速度傳感器在被測齒輪箱上采集振動信號(無需同步采集轉速信號);

(2)對采集的振動信號進行加窗濾波處理;

(3)對初步處理后的振動信號進行時頻分析,得到時頻譜;

(4)采用自適應多階比搜峰法,從最低頻率開始自動搜索第一個階比分量的時頻分布峰值帶寬,確定搜峰的起始點坐標,繼續向上搜索,按照從低到高的順序,依次搜索出所有階比分量對應帶寬及起始點坐標;

圖2 基于瞬時頻率估計的AMOT算法流程圖 Fig.2 The flow chart of the AMOT algorithm based on instantaneous frequency estimation

(5)采用最小二乘方法擬合搜索到的所有階比分量曲線;

對于給定的數據點(xi,yi),1≤i≤N,可用下面的n階多項式進行擬合,即

(9)

為了使擬合出的近似曲線能盡量反映所給數據的變化趨勢,要求在所有數據點上的殘差都較小。

|δi|=|f(xi)-yi|

(10)

為達到上述目標,可以令上述偏差的平方和最小,即

(11)

這種方法稱為最小二乘原則,利用這一原則確定擬合多項式f(x)的方法即為最小二乘法多項式擬合。

(6)根據瞬時頻率的擬合曲線產生相應的鑒相時標Tn。

(12)

對于旋轉機械升、降速階段,若fi(t)滿足光滑條件,在小范圍內用擬合多項式或樣條方程實現高精度逼近是可行的,因此若fi(t)=at2+bt+c。

則上式可以表述為

(13)

求解此方程的有效值,即可求得等角度采樣的鑒相時標Tn。然后根據等角度階比重采樣時刻Tn,對初始等時間采樣得到的信號進行三次樣條插值以實現等角度階比重采樣,得到的角域周期循環平穩信號x(Tn)。

根據鑒相時標Tn對原始采樣數據進行插值以實現信號的等角度插值重采樣,進而得到所有階比采樣的角域離散轉角序列x(Tn);

(7)對離散轉角序列x(Tn)進行FFT,得到全部的階比分量。

3仿真實驗

3.1仿真信號1

圖3 x(t)的時域圖Fig.3Thetimedomainwaveformofx(t)圖4 x(t)的頻域圖Fig.4Thefrequencyspectrumofx(t)

圖5為信號x(t)的STFT時頻譜圖,從圖5可知信號所包含的頻率分量以及其隨時間變化的趨勢。圖6為采用本文提出的自適應搜峰算法的結果圖,可以看出,仿真信號包含的兩個分量均可以提取出來,圖6(a)為低階分量的瞬時頻率曲線所對應的搜峰結果,圖6(b)為高階分量信號的瞬時頻率曲線所對應的搜峰結果。對比圖5和圖6,直觀上可以看出自動搜峰方法能夠實現信號的搜峰過程,具有自適應性。圖7(a)為對低階分量的最小二乘的擬合過程,圖7(b)為對高階分量的最小二乘的擬合過程,從圖7可知,最小二乘的擬合結果比較光滑,能夠表現信號峰值的變化過程。通過對比理想曲線和擬合曲線可知,二者基本重合且具有較高的精度。

圖5 信號的STFT時頻譜圖 Fig.5 The time-frequency spectrum of signal by STFT

圖6 AMOT自動搜峰結果 Fig.6 Theresults of AMOT automatic peak search

圖7 分量的最小二乘擬合結果 Fig.7 The least square polynomial result of the component

以不同分量為參考分量進行階比分析,圖8為按照不同參考分量計算鑒相時標重采樣后的角域圖,可以看出信號比較平穩。

圖8 重采樣后的角域信號 Fig.8 Angle domain signal after resampling

圖9為以兩個分量分別為參考分量計算的階比譜圖,通過計算各個階比分量的幅值關系,幅值較大的階比譜圖作為最后的顯示分析結果,實現優中選優,避免了傳統算法中人為直觀選取一個分量進行遮掩濾波提取分量的方法,減少了人為選取分量及起始點造成的誤差,具有自適應性。圖9(a)中低階分量為基頻1階,對應的高階分量對應為4階,兩個分量的幅值與理想值非常接近。圖9(b)中高階分量為1階,則相應的低階分量為0.25階,由于環節誤差及能量泄露導致低頻分量的幅值與理想值具有一定偏差。

從圖9階比譜圖中能夠清晰的顯示出與基頻分量相關的頻率分量,與FFT頻譜圖形成了鮮明的對比,可以得出,階比分析技術比較適合于分析變頻振動信號,特別是頻率隨轉速變化的變頻振動信號。

圖9 仿真信號1的階比譜 Fig.9 The order spectrum of the first simulation signal

3.2升降速仿真信號2

為進一步說明算法的有效性,用旋轉機械升速、平穩和降速的仿真實驗信號進行階比分析。仿真信號表達式如下:

y(t)=sin4πf(t)t+sin8πf(t)t+

sin12πf(t)t+sin2πf0t

(14)

式中:f0=6 Hz,t∈[2,7.913](s),用以仿真旋轉機械的低頻共振。調頻函數f(t)分為三段:第一段仿真旋轉機械的升速過程,第二段仿真旋轉機械的平穩運轉過程,第三段仿真旋轉機械的降速過程。其瞬時頻率fi(t)定義如下:

(15)

圖10為信號y(t)的時域圖,用來仿真信號的升速、平穩、降速的過程。圖11為信號的STFT時頻譜圖,很明顯的可以看出該信號中有1階,2階,3階分量,同時還有一個頻率很低的直線。從圖11可知,某一瞬時所對應的頻率成分及其分布情況。

圖10 y(t)的時域圖Fig.10Timedomainwaveformofy(t)圖11 y(t)的STFT時頻圖Fig.11Thetime-frequencyspectrumofy(t)bySTFT

圖12~14圖為自適應多階比分析系統最后的顯示結果,圖12為自動搜索出來的第一階分量,圖13為重采樣信號的角域信號,圖14為按第一階分量分析的階比譜圖,各階分量的幅值譜都比較清晰,驗證了方法的有效性。

圖12 第一階分量的峰值Fig.12Thepeaksresultofthefirstordercomponent圖13 重采樣角域信號Fig.13Theangledomainsignalafterresampling

圖14 仿真信號2的階比譜 Fig.14 The order spectrum of the second simulation signal

4應用實例

某850 kW風電機組,齒輪箱結構由一級行星輪傳動和二級平行軸傳動組成。傳感器測點布置位置見圖15,齒輪箱高速端傳感器布置實物見圖16。

圖15 傳感器測點布置圖 Fig.15 Measuring points layout of sensor

圖16 高速端測點傳感器布置圖 Fig.16 The sensor layout of the high speed shaft

對機組高速端采集的振動信號見圖17,采樣率為2 500 Hz,采樣時間為3.28 s,時域圖中存在明顯的沖擊,但是不能判斷故障類型及位置。通過本文提出的自適應多階比方法分析得到時頻圖見圖18,其中,自適應階比分析方法展示的結果是以高速端轉頻的3倍頻(圖18中圈和箭頭所示)作為參考頻率進行分析,搜峰的結果見圖19,圖中X坐標為時間,Y坐標為頻率,Z坐標為瞬時頻率強度。經過最小二乘擬合得到圖20,最后經過階比分析得到階比譜圖21及其局部放大圖22。從階比圖圖21可知,高速端轉頻倍頻現象明顯,圖中在3.574階及其2倍階附近均存在間隔為高速端轉頻的諧波分量,且高速端小齒輪與中速端小齒輪的嚙合頻率(8.965階,約為高速端轉頻的25倍)較突出,其邊頻調制現象也較明顯,邊頻帶寬為高速端轉頻。圖22將階比譜放大(0~5階之間,圖21圈位所示),觀察低頻段的頻譜特征,存在高速端轉頻的1/3、2/3等倍頻成分,根據齒輪的傳動比關系,判斷為中速軸的轉頻,其幅值較高,進而判斷齒輪箱中速軸齒輪及高速軸齒輪存在損傷。圖23為原始信號的頻譜圖,從圖可知,圖中存在諧波分量,但由于底部噪聲較大,信噪比較低,信號特征頻率的幅值較低,同時如圖24,原始頻譜圖的放大效果顯示(0~200 Hz之間,圖23圈位所示),在低頻段的特征信息不明顯。對比分析結果,可以看出本文提出的方法能夠有效地實現齒輪故障的診斷。

圖17 高速端振動信號Fig.17Vibrationsignalofthehighspeedshaft圖18 高速端振動信號時頻圖Fig.18Thetime-frequencyspectrumofthehighspeedshaftsignal

圖19 階比搜峰結果Fig.19Theresultsoforderpeaksearch圖20 分量擬合結果Fig.20Fittingresultsofcomponent

圖21 高速端振動信號的階比譜 Fig.21 The order spectrum of vibration signal of high speed shaft

圖22 高速端振動信號的階比譜低頻段局部放大 Fig.22 The order spectrum local amplification of vibration signal of high speed shaft

圖23 高速端振動信號的原始頻譜圖 Fig.23 The original spectrum of vibration signal of high speed shaft

圖24 高速端振動信號的原始頻譜低頻段放大圖 Fig.24 The original spectrum local amplification of vibration signal of high speed shaft

經現場工程師登塔對齒輪箱進行內窺鏡檢查,如圖25所示,高速端小齒輪存在明顯的磨損和劃痕(圖25中左側兩個圈位所示),中速端齒輪也存在磨損(圖25中右側圈位所示),驗證了方法的有效性。

圖25 實際故障照片 Fig.25 The actual fault photos

5結論

提出一種用于提取旋轉機械瞬時頻率的自適應頻率峰值搜索算法,在此基礎上提出基于瞬時頻率估計的無轉速計自適應多階比分析新方法。給出了算法的

步驟和實現過程。該方法無需同步采集轉速信號,也不需要指定具體搜索階比,能夠自動識別出所有階比分量。仿真信號和實例分析結果驗證了方法的有效性。實例分析表明,所提出的自適應峰值搜索法能夠有效提取出不同階比分量,準確識別出風電機組高速級齒輪故障。

參考文獻

[1]Potter R. A new order tracking method for rotating machinery[J]. Sound and Vibration, 1990, 24(9):30-34.

[2]Fyfe K R, Munck E D S. Analysis of computed order tracking[J]. Mechanical Systems and Signal Processing,1997, 11(2):187-205.

[3]郭瑜, 秦樹人, 梁玉前. 時頻分析階比跟蹤技術[J]. 重慶大學學報, 2002, 25(5): 17-20.

GUO Yu, QIN Shu-ren, LIANG Yu-qian. Order tracking method based on time-frequency analysis[J]. Journal of Zhongqing University, 2002, 25(5):17-20.

[4]趙曉平, 張令彌, 郭勤濤. 基于瞬時頻率估計的自適應Vold-Kalman 階比跟蹤研究[J]. 振動與沖擊, 2008, 27(12):32-36.

ZHAO Xiao-ping, ZHANG Ling-mi, GUO Qin-tao. Adaptive Vold-Kalman order tracking based on instantaneous frequency estimation[J]. Journal of Vibration and Shock, 2008, 27(12):112-116.

[5]楊志堅, 丁康, 楊茜. 基于頻譜校正理論的階比跟蹤分析[J]. 機械工程學報, 2009, 45(12):41-45.

YANG Zhi-jian, DING Kang, YANG Xi. Novel method of order tracking analysis based on spectrum correction[J]. Journal of Mechanical Engineering, 2009,45(12):41-45.

[6]Blough J R. Development analysis of time variant discrete Fourier transform order tracking[J]. Mechanical Systems and Signal Processing, 2003, 17(6): 1185-1199.

[7]郭瑜, 秦樹人,湯寶平,等. 基于瞬時頻率估計的旋轉機械階比跟蹤[J]. 機械工程學報, 2003, 40(9): 1648-1651.

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

[8]趙曉平,趙秀莉,侯榮濤,等.一種新的旋轉機械升降速階段振動信號瞬時頻率的估計算法[J]. 機械工程學報, 2011, 47(7):104-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]. Journal of Mechanical Engineering, 2011, 47(7):104-108.

[9]胡愛軍,朱瑜. 基于改進峰值搜索法的旋轉機械瞬時頻率估計, 振動與沖擊, 2013, 32(7):113-117.

HU Ai-jun, ZHU Yu. Instantaneous frequency estimation of a rotating machinery based on an improved peak searchmethod[J]. Journal of Vibration and Shock, 2013, 32(7):113-117.

猜你喜歡
故障診斷
基于包絡解調原理的低轉速滾動軸承故障診斷
一重技術(2021年5期)2022-01-18 05:42:10
ILWT-EEMD數據處理的ELM滾動軸承故障診斷
水泵技術(2021年3期)2021-08-14 02:09:20
凍干機常見故障診斷與維修
基于EWT-SVDP的旋轉機械故障診斷
數控機床電氣系統的故障診斷與維修
電子制作(2018年10期)2018-08-04 03:24:46
基于改進的G-SVS LMS 與冗余提升小波的滾動軸承故障診斷
因果圖定性分析法及其在故障診斷中的應用
改進的奇異值分解在軸承故障診斷中的應用
基于LCD和排列熵的滾動軸承故障診斷
基于KPCA和PSOSVM的異步電機故障診斷
主站蜘蛛池模板: 在线看AV天堂| 99在线视频免费| a色毛片免费视频| 亚洲精品午夜天堂网页| 亚洲欧美激情另类| 亚洲bt欧美bt精品| 欧美激情第一区| 精品无码人妻一区二区| 97无码免费人妻超级碰碰碰| 久久6免费视频| 喷潮白浆直流在线播放| 欧美亚洲香蕉| 久久国产精品国产自线拍| 中文字幕佐山爱一区二区免费| 蜜桃视频一区二区三区| 欧美色视频在线| 中文字幕亚洲专区第19页| 丁香婷婷久久| 久草青青在线视频| 欧美日韩国产精品va| 精品人妻无码中字系列| 亚洲精品无码专区在线观看| 国产肉感大码AV无码| 国产成a人片在线播放| 国产视频 第一页| 亚洲制服丝袜第一页| 久久国产乱子伦视频无卡顿| 狠狠色综合久久狠狠色综合| 久久国产黑丝袜视频| 在线a网站| 久久女人网| 狠狠干欧美| 嫩草影院在线观看精品视频| 宅男噜噜噜66国产在线观看| 国产毛片高清一级国语 | 一边摸一边做爽的视频17国产| 91无码网站| 青青草一区| 久久久久国产一区二区| 中文字幕天无码久久精品视频免费 | 91在线国内在线播放老师| 欧美a在线| 久草视频精品| 亚洲精选无码久久久| 国产精品偷伦在线观看| 国产欧美在线观看一区| v天堂中文在线| 欧美中文字幕在线二区| 久久无码av一区二区三区| 免费在线一区| 日韩少妇激情一区二区| 久久黄色影院| 无码免费的亚洲视频| 91久久国产热精品免费| 五月激情婷婷综合| 亚洲成人高清无码| 免费在线国产一区二区三区精品| 午夜国产小视频| 高清亚洲欧美在线看| 四虎成人精品| 欧美精品亚洲精品日韩专区va| 台湾AV国片精品女同性| 国产精品美女在线| 国产精品99r8在线观看| 亚洲大学生视频在线播放| 亚洲九九视频| 一区二区三区国产精品视频| 波多野结衣无码AV在线| 欧美日韩另类在线| 九九免费观看全部免费视频| 韩日无码在线不卡| 久久综合色88| 曰韩免费无码AV一区二区| 波多野结衣无码AV在线| 2020国产在线视精品在| 91视频精品| 99久久精品免费观看国产| 午夜视频免费一区二区在线看| 成人综合网址| 无码中文字幕精品推荐| 亚洲中文在线视频| 国产成人av大片在线播放|