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

基于改進變分模態分解的旋轉機械故障時頻分析方法①

2016-02-09 11:23:07劉尚坤唐貴基王曉龍
振動工程學報 2016年6期
關鍵詞:模態故障信號

劉尚坤, 唐貴基, 王曉龍

(華北電力大學機械工程系, 河北 保定 071003)

基于改進變分模態分解的旋轉機械故障時頻分析方法①

劉尚坤, 唐貴基, 王曉龍

(華北電力大學機械工程系, 河北 保定 071003)

變分模態分解(variational mode decomposition, VMD)算法作為一種新的多分量信號分解方法,采用頻域迭代的求解方式,具有明確的理論基礎和高的分解精度,為了自適應確定其分解分量的個數,以互信息為判據對原方法進行了迭代停止條件的改進,并結合Teager能量算子具有對單分量信號解調速度快、精度高的優點,提出了Teager-VMD時頻分析新方法。仿真信號分析表明該方法能將含有變頻、頻率突變和頻率相近的多分量信號進行有效的分離;實測轉子系統局部碰摩、油膜振蕩信號分析表明,該方法能夠精確提取各分量的清晰時頻特征,能準確地辨別碰摩故障的嚴重程度和診斷油膜振蕩故障,較HHT方法精確有效,具有一定的工程應用價值。

故障診斷;時頻分析; 碰摩;Teager能量算子; 油膜振蕩

引 言

汽輪機等旋轉機械由于結構復雜、工況多變,故障振動信號常表現為多分量非平穩特點。由于時頻分析方法能同時提供非平穩信號在時域和頻域的局部化信息而得到了廣泛的應用[1-6],而對多分量信號分析的關鍵是找到一種有效的信號分解方法[1]。針對這一問題學者們進行了不斷的探索,文獻[2]利用希爾伯特-黃變換(Hilbert-Huang Transform,HHT)分析和監測裂紋轉子的瞬態響應,文獻[3]利用建立在EMD基礎上的Teager-Huang對轉子碰摩故障進行了診斷,但由于EMD存在模式混淆、過包絡、欠包絡等問題[4]影響分析效果;文獻[5]利用經驗小波變換(Empirical Wavelet Transform, EWT)分析了碰摩故障的頻率分布以判別碰摩故障的嚴重程度;文獻[6]采用希爾伯特振動分解(Hilbert Vibration Decomposition,HVD)診斷了轉子系統油膜渦動故障。

Dragomiretskiy等提出了一種新的多分量信號自適應分解方法——變分模態分解(VMD)[7],該方法采用頻域非遞歸的迭代求解方式搜尋變分模型最優解來確定每個調幅調頻分量的中心頻率及帶寬,最終自適應的實現信號的頻域剖分及各分量的分離,分解精度高且能有效避免模式混淆問題[8]。文獻[9]利用VMD分解得到的中心頻率作為聚類中心并結合標準模糊C均值聚類診斷了滾動軸承故障。本文在介紹VMD算法的基礎上,提出基于互信息準則的自適應確定變分模態分解分量個數的改進算法,并與Teager能量算子解調相結合,提出了Teager-VMD時頻分析新方法,仿真信號和實測轉子局部碰摩、油膜失穩信號分析效果說明了該方法的有效性。

1 Teager-VMD時頻分析方法原理

1.1 變分模態分解

變分模態分解將本征模態函數(IMF)定義為一個調幅調頻信號,即

uk(t)=Ak(t)cos[φk(t)]

(1)

假設原信號f為多分量信號,由K個有限帶寬的本征模態函數(IMF)分量uk組成,且各IMF的中心頻率為ωk,VMD算法建立的約束變分模型為

(2)

式中 uk={u1, … ,uK},ωk={ω1, …,ωK}。

該模型中,通過Hilbert變換得到uk(t)的解析信號進而得到單邊譜,乘以指數函數e-jωkt是將所估計的uk(t)的中心頻帶調整到基頻帶上,為求上述約束變分問題的最優解,需將其轉換為非約束變分問題,為此引入如下形式的增廣Lagrange函數

(3)

式中 α為二次項懲罰參數;λ為Lagrange乘子。利用乘子交替方向算法不斷更新各IMF及其中心頻率,最終所求式(3)的鞍點即為原問題的最優解,而所有的IMF可從頻域中通過下式獲得

(4)

(5)

1.2 VMD分量個數的自適應確定

由VMD算法原理可知,其處理信號時需要預先設定分解IMF分量的個數K,然而受不同設備不同工況的限制,K值通常難以準確設定。為解決這一問題,本文利用互信息為判據來控制VMD算法的迭代過程,提出自適應確定分量個數K的改進VMD算法。

互信息(mutual information,MI)由信息論中熵的概念引申而來,用兩個隨機變量間不確定度的差值表示,能夠表明其統計相關性,比相關系數法更準確[10],能更好地辨別相關程度。互信息表達式如下

MI(X,Y)=H(Y)-H(Y|X)

(6)

式中H(Y)為Y的熵,H(Y|X)為已知X時Y的條件熵。

上式表明,X與Y的相關性越強,條件熵值H(Y|X)越小,則互信息MI(X,Y)越大[11],本文將其由下式

δi=MIi/max(MIi)

(7)

進行歸一化后,用于判斷分解余量與原信號的相關程度,當歸一化互信息值低于閾值(參考文獻[10],本文取閾值δ=0.02)時,認為分解余量不再含有重要信息,原信號已被完全分解,結束整個運算過程,從而能夠根據實際分析信號自動地確定IMF分量的個數,具體實現步驟如下:

(1) 初始化K為1;

(2)K=K+1,執行外層循環;

(4) n=n+1,開始內層循環;

1.3 Teager能量算子解調

本研究針對當前應用型心理學課程教學存在的弊端,在心理衛生學教學中,整合合作學習和探究性學習兩種學習方式,以學生學習小組為單位,結合相關教學內容,通過探究性學習,研究和解決實際心理問題,深化學生對心理衛生學理論的理解,培養學生的心理衛生學應用能力。

要表達K個分量的時頻分布信息,還需對得到的每個IMF分量進行解調分析,Teager能量算子是一種非線性差分算子,解調能力比Hilbert變換解調能力強[12],同時具有時間分辨率高、解調速度快的優點,能夠對單分量的調幅調頻信號解調[13],近些年被廣泛用于信號解調分析[14-15]。

對于調幅調頻的連續時間信號x(t),其Teager能量算子定義為

(8)

對應離散信號x(n)的能量算子定義為

ψ[x(n)]=x2(n)-x(n+1)x(n-1)

(9)

由式(9)可知每一時刻算子的計算只需要3個樣本數據,計算量較小,由式(10)和(11)實現單分量調幅調頻信號的瞬時頻率與瞬時幅值的求解。

(10)

(11)

相對于Hilbert變換,Teager能量算子在求信號的瞬時頻率時,無需進行復數計算,具有快速響應能力,可以迅速跟蹤瞬時變化信號的幅值和頻率,能更好地揭示瞬時頻率的突變情況。

1.4 Teager-VMD時頻分析方法

本文提出的Teager-VMD時頻分析方法,首先利用改進VMD算法自適應將一多分量信號分解為K個表征信號特征的IMF分量(同時參考文獻[7],本文設定懲罰參數α=300),再通過Teager能量算子計算每個IMF分量的瞬時幅值和瞬時頻率,得到原信號的時頻分布信息,最后通過分析時頻分布圖中的故障特征頻率診斷故障類型。基于改進變分模態分解的機械故障時頻分析方法的分析流程如圖1所示。

圖1 Teager-VMD時頻分析流程圖Fig.1 Teager-VMD time-frequency analysis flow chart

具體實現步驟如下:

(2)通過Teager能量算子解調計算每一個IMF分量的瞬時幅值和瞬時頻率。

(3)繪制三維時頻分布圖,分析其中是否存在相應的故障特征頻率,進而判斷故障類型,實現故障診斷。

2 仿真信號分析

利用仿真信號來驗證本文方法的時頻分析能力,考查一含有噪聲的多分量信號x(t),其中分量x1(t)為調幅調頻信號,x2(t)是在兩個不同時間段頻率單一的頻率突變信號,x3(t)是頻率接近的兩個余弦信號,x4(t)是噪聲模擬信號。具體表達式如下式

(12)

采樣頻率為2048 Hz,分析點數為1024點。圖2為多分量信號x(t)的時域波形。首先,利用改進VMD算法對x(t)進行分解,得到5個IMF分量,如圖3所示,圖中各個分量被較好地分解出來,各個IMF與原信號的歸一化互信息如表1所示,其中前4個分量大于閾值0.02。圖4是對前4個IMF分量進行能量算子解調得到的時頻分布圖,由圖可知,VMD能較好的把調幅調頻分量x1(t)、頻率突變分量x2(t)以及x3(t) 中頻率接近的200與240 Hz分量進行有效分離,清晰地表示了各分量的時頻信息。

為了對比分析,圖5是采用HHT分析該仿真信號得到的Hilbert譜,由圖可知,各分量頻率調制現象明顯,分解得到的調幅調頻分量x1(t)不夠明確,頻率突變分量x2(t)中的130 Hz分量尚能基本識別,而160 Hz分量不能識別,x3(t)中兩個頻率接近的分量不能有效分離,效果較本文方法差。

圖2 仿真信號的時域波形Fig.2 Time domain waveform of simulation signal

圖3 仿真信號的自適應變分模態分解Fig.3 Adaptive variational mode decomposition of simulation signal

表1 各IMF與仿真信號的歸一化互信息值

Tab.1 Normalized mutual information value between IMFS and simulation signal

IMF1IMF2IMF3IMF4IMF50.90130.74901.00000.88720.0003

圖4 仿真信號的Teager-VMD時頻分布圖Fig.4 Teager-VMD time-frequency distribution of simulation signal

圖5 仿真信號HHT分析的Hilbert譜Fig.5 Simulation signal′s Hilbert spectrum by HHT

3 轉子實驗振動信號分析

實驗裝置采用Bently RK-4轉子系統故障模擬實驗臺(如圖6所示),配有信號前置適配器、轉速控制調節裝置和做油膜失穩所需的軸承、油泵系統。在單盤轉子兩側安裝電渦流傳感器測量轉軸徑向振動位移,由美國Iotech 公司生產的 ZonicBook/618E設備采集轉子在碰摩過程中以及在油膜失穩狀態下的振動信號,采樣頻率為1280 Hz,本文取轉子各實驗狀態的1280個采樣點加以分析。

圖6 轉子實驗臺Fig.6 Rotor experiment platform

3.1 碰摩故障分析

在距離轉軸小間隙處由支架固定一銅質頂針,轉軸升速過程中與銅質頂針發生局部碰摩,隨著轉速的升高局部碰摩程度由輕微發展到嚴重,采集振動信號進行分析。

圖7是轉速為1 560 r/min時采集的輕微碰摩信號的時域波形和頻譜圖,從時域波形上基本上看不出發生了碰摩故障,頻譜圖中也只是轉頻26 Hz幅值明顯,難以判斷已經發生了碰摩故障。

圖7 輕微碰摩信號的時域波形和頻譜圖Fig.7 Time domain waveform and spectrum of slight rub-impact signal

圖8 輕微碰摩信號的自適應變分模態分解Fig.8 Adaptive variational mode decomposition of slight rub-impact signal

圖9 輕微碰摩信號的Teager-VMD時頻分布圖Fig.9 Teager-VMD time-frequency distribution of slight rub-impact signal

圖8是采用改進VMD分解得到的8個IMF分量,周期性明顯的分量中,IMF7對應轉頻,IMF4對應二倍頻,IMF3對應三倍頻。對大于閾值的前7個IMF分量進行Teager能量算子解調得到圖9所示的輕微碰摩信號的Teager-VMD時頻分布圖,圖中轉頻幅值最大,同時存在幅值較小的二倍頻、三倍頻成分,也存在幅值很小且頻率波動的三分頻及高次諧波,通過分析,結合碰摩故障特征[16]可以準確地確定轉子已經發生了碰摩且程度輕微,實現了早期轉子碰摩故障的識別。

為了對比本文方法的分析效果,圖10是將輕微碰摩信號進行HHT分析得到的Hilbert譜,圖中只能分析出轉頻具有頻率調制特點,難以準確地判斷一定發生了碰摩故障,相比本文分析方法較差。

圖10 輕微碰摩信號HHT分析的Hilbert譜Fig.10 Slight rub-impact signal′s Hilbert spectrum by HHT

圖11是轉速為1740 r/min時采集的嚴重碰摩信號時域波形和頻譜圖,圖中幅值變大說明碰摩程度及能量變大,時域波形中已經發生了較明顯的變形,頻譜圖中仍是轉頻29 Hz幅值突出,但能夠觀察到小幅值的二倍頻、三倍頻成分。該信號的自適應變分模態分解如圖12所示,共得到11個IMF分量,分量間存在著倍頻特征,對大于閾值的前10個IMF分量做能量算子解調得到的時頻分布圖如圖13所示,圖中各分量的時頻譜清晰可見,其中2~7倍頻成分被準確檢測出來,且2~4倍頻幅值相對較大,5~7倍頻幅值相對較小,同時圖中還有10倍頻的調制頻率成分出現,此頻譜特征充分表明轉子發生了嚴重的碰摩故障。

圖11 嚴重碰摩信號的時域波形和頻譜圖Fig.11 Time domain waveform and spectrum of serious rub-impact signal

圖12 嚴重碰摩信號的自適應變分模態分解Fig.12 Adaptive variational mode decomposition of serious rub-impact signal

圖13 嚴重碰摩信號的Teager-VMD時頻分布圖Fig.13 Teager-VMD time-frequency distribution of serious rub-impact signal

為了對比分析,圖14是將嚴重碰摩信號進行HHT分析得到的Hilbert譜,圖中雖能顯示出被調制的轉頻29 Hz及高頻成分,但不能詳細分析各頻率特征,對碰摩的嚴重程度不易準確辨別,分析效果明顯不如本文分析效果。

圖14 嚴重碰摩信號HHT分析的Hilbert譜Fig.14 Serious rub-impact signal′s Hilbert spectrum by HHT

3.2 油膜失穩故障分析

轉子系統油膜失穩常由半速油膜渦動發展為油膜振蕩[17],由于油膜振蕩的頻譜結構較油膜渦動復雜,下面分析實驗臺轉子1階油膜振蕩時的頻譜特征。

圖15 油膜振蕩信號的時域波形和頻譜圖Fig.15 Time domain waveform and spectrum of oil-whip

圖16 油膜振蕩信號的自適應變分模態分解Fig.16 Adaptive variational mode decomposition of oil-whip

圖17 油膜振蕩信號的Teager-VMD時頻分布圖Fig.17 Teager-VMD time-frequency distribution of oil-whip

為了進行分析效果對比,圖18是采用HHT分析油膜振蕩信號得到的Hilbert譜,圖中800~1000采樣點時段出現了明顯的頻率混淆現象,雖然分析出振蕩頻率31 Hz及其二分頻成分,但轉頻76 Hz出現明顯的調制現象,其他組合頻率也不能分解出來,比本文方法分析效果差。

圖18 油膜振蕩信號HHT分析的Hilbert譜Fig.18 Oil-whip signal′s Hilbert spectrum by HHT

4 結 論

(1)結合變分模態分解與Teager能量算子解調各自的優點,本文提出了一種Teager-VMD時頻分析方法。為了在信號分解過程中自適應地確定VMD分量的個數,采用互信息準則對VMD的迭代停止條件進行了改進,實現了分量個數的自適應確定。對比分析表明該方法效果明顯優于HHT方法,具有一定的工程應用價值。

(2)將Teager-VMD時頻分析方法應用于轉子實驗信號,精確地分析出轉子輕微碰摩和嚴重碰摩時的頻率成分以及油膜振蕩時的轉頻、振蕩頻率及其和差組合頻率成分,準確地診斷了故障,也表明該方法提取信號時頻特征的精確性和準確性。

(3)本文僅對單盤轉子的局部碰摩及油膜振蕩現象進行了時頻分析與故障診斷,取得了較好的效果,所提方法在其他復雜設備、復雜工況及其他行業領域中的應用有待進一步研究。

[1] 程軍圣, 楊宇, 于德介. 基于廣義解調時頻分析的多分量信號分解方法[J]. 振動工程學報, 2007, 20(6):563—569.

CHENG Junsheng, YANG Yu, YU Dejie. A multi component signal decomposition method based on the generalized demodulation time frequency analysis[J]. Journal of Vibration Engineering, 2007, 20(6):563—569.

[2] Ramesh Babu T, Srikanth S, Sekhar A S. Hilbert-Huang transform for detection and monitoring of crack in a transient rotor [J]. Mech. Syst. Signal Process, 2008, 22 (4): 905—914.

[3] Liu Shangkun,Tang Guiji,Pang Bin. Rotor local rubbing fault feature analysis based on Teager-Huang transform [C]. 4th International Conference on Materials Science and Information Technology, Tianjin, China, 2014:3244—3247.

[4] 向玲,鄢小安. 汽輪機轉子故障診斷中LMD法和EMD法的性能對比研究[J]. 動力工程學報, 2014, 34(12): 945—951.

XIANG Ling, YAN Xiaoan. Performance contrast between LMD and EMD in fault diagnosis of turbine Rotors [J]. Journal of Chinese Society of Power Engineering, 2014, 34(12): 945—951.

[5] 李志農,朱明,褚福磊,等. 基于經驗小波變換的機械故障診斷方法研究[J]. 儀器儀表學報,2014, 35(11): 2423—2432.

LI Zhinong, ZHU Ming, CHU Fulei, et al. Mechanical fault diagnosis method based on empirical wavelet transform [J]. Chinese Journal of Scientific Instrument, 2014, 35(11): 2423—2432.

[6] 唐貴基,龐彬. 基于改進的希爾伯特振動分解的機械故障診斷方法研究[J]. 振動與沖擊,2015, 34(3): 167—171.

Tang Guiji, Pang Bin. Research for a mechanical fault diagnosis method based on improved Hilbert vibration decomposition [J]. Journal of Vibration and Shock, 2015, 34(3): 167—171.

[7] Dragomiretskiy Konstantin, Zosso Dominique. Variational mode decomposition [J]. IEEE Transactions on Signal Processing,2014,62(3): 531—544.

[8] 唐貴基,王曉龍. 參數優化變分模態分解方法在滾動軸承早期故障診斷中的應用[J]. 西安交通大學學報,2015,49(5): 73—81.

TANG Guiji, WANG Xiaolong. Parameter optimized variational mode decomposition method with application to incipient fault diagnosis of rolling bearing [J]. Journal of Xi′an Jiaotong University, 2015, 49(5):73—81.

[9] 劉長良, 武英杰, 甄成剛. 基于變分模態分解和模糊C均值聚類的滾動軸承故障診斷[J]. 中國電機工程學報, 2015, 35(13): 3358—3365.

LIU Changliang, WU Yingjie, ZHEN Chenggang. Rolling bearing fault diagnosis based on variational mode decomposition and fuzzy c means clustering [J]. Proceedings of the CSEE, 2015, 35(13): 3358—3365.

[10]胡愛軍. Hilbert-Huang 變換在旋轉機械振動信號分析中的應用研究[D].保定:華北電力大學, 2008.

Hu Aijun. Research on the application of Hilbert-Huang Transform in vibration signal analysis of rotating machinery [D].Baoding, North China Electric Power University, 2008.

[11]張志剛,石曉輝,施全,等. 基于改進EMD和譜峭度法滾動軸承故障特征提取[J]. 振動、測試與診斷, 2013, 33(3): 478—530.

Zhang Zhigang, Shi Xiaohui, Shi Quan, et al. Fault feature extraction of rolling element bearing based on improved emd and spectral kurtosis [J].Journal of Vibration,Measurement & Diagnosis,2013, 33(3): 478—530.

[12]鐘先友.旋轉機械故障診斷的時頻分析方法及其應用研究[D].武漢:武漢科技大學,2014.

Zhong Xianyou. Research on time-frequencyanalysis methods and its applications to rotating machinery fault diagnosis[D].Wuhan:Wuhan University of Science and Technology,2014.

[13]Marago S P,Kaiser J F,Quatieri T F.Energy separation in signal modulations with applications to speech analysis [J].IEEE Transactions on Signal Processing,1993,41(10):3024—3051.

[14]孟宗, 李姍姍, 季艷. 基于對稱差分能量算子解調的局部均值分解端點效應抑制方法[J]. 機械工程學報, 2014, 50(13): 80—87.

MENG Zong, LI Shanshan, JI Yan. Restraining method for end effect of local mean decomposition based on energy operator demodulation of symmetrical differencing[J]. Chinese Journal of Mechanical Engineering , 2014, 50(13): 80—87.

[15]張文義,于德介,陳向民.齒輪箱復合故障診斷的信號共振分量能量算子解調方法[J].振動工程學報,2015,28(1):148—155.

Zhang Wenyi,YU Dejie,Chen Xiangmin.Energy operator demodulating of signal′s resonance components for the compound fault diagnosis of gearbox [J]. Journal of Vibration Engineering, 2015,28(1):148—155.

[16]馬輝, 楊健, 宋溶澤, 等. 轉子系統碰摩故障實驗研究進展與展望[J]. 振動與沖擊, 2014, 33(6): 1—12.

MA Hui, YANG Jian, SONG Rongze, et al. Review and prospect on the research of rub-impact experiment of rotor systems [J]. Journal of Vibration and Shock, 2014, 33(6): 1—12.

[17]唐貴基,向玲,朱永利. 基于 HHT的旋轉機械油膜渦動和油膜振蕩故障特征分析[J]. 中國電機工程學報,2008, 28(2):77—81.

TANG Guiji, XIANG Ling, ZHU Yongli. Fault analysis of oil whirl and oil whip based on Hilbert Huang transform for rotor system [J]. Proceedings of the Chinese Society for Electrical Engineering, 2008, 28(2): 77—81.

[18]馬輝,李輝,唐玉生,等. 兩種不同載荷形式下轉子系統油膜失穩的數值研究[J]. 振動工程學報,2013, 26(1): 105—111.

MA Hui, LI Hui, TANG Yusheng, et al. Numerical research on oil-film instability in a rotor system under two types of load conditions [J]. Journal of Vibration Engineering, 2013, 26(1): 105—111.

Time frequency analysis method for rotary mechanical fault based
on improved variational mode decomposition

LIUShang-kun,TANGGui-ji,WANGXiao-long

(School of Mechanical Engineering,North China Electric Power University, Baoding 071003, China)

As a new method for multi-component signal's decomposition, variational mode decomposition (VMD) has explicit theoretical basis and high decomposition accuracy. In order to adaptively determine the number of its decomposition components, the criterion of mutual information is used to improve the iterative stopping conditions to the original method. Combining with the advantages of fast and high accuracy of Teager energy operator demodulation, a new method named Teager-VMD is proposed for time frequency analysis. The simulation signal analysis results show that the proposed method can effectively separate the multi-component signal with frequency conversion, frequency mutation and frequency similar signal. The analysis results of the local rub-impact and oil whip of the experiment rotor system show that the proposed method can accurately extract the clear time-frequency characteristics and can accurately identify the severity of the rub-impact fault and diagnosis the oil whip fault. The proposed method is more accurate and effective than the HHT method and has certain engineering application value.

fault diagnosis; time frequency analysis; rub-impact; Teager energy operator; oil whip

2015-10-21;

2016-02-28

國家自然科學基金資助項目(51307058,51475164);河北省自然科學基金資助項目(E2014502052, E2015502013);中央高校基本科研業務費專項資金資助項目(2014MS156,2015XS120)

TH165+.3;TN911.7

1004-4523(2016)06-1119-08

10.16385/j.cnki.issn.1004-4523.2016.06.022

劉尚坤(1979—),男,博士研究生,講師。電話:(0312)7523442; E-mail:lsk1213@163.com

唐貴基(1962—) ,男, 教授, 博士生導師。電話:(0312)7525028; E-mail:tanggjlk@ncepubd.edu.cn

猜你喜歡
模態故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
故障一點通
基于HHT和Prony算法的電力系統低頻振蕩模態識別
江淮車故障3例
主站蜘蛛池模板: 动漫精品中文字幕无码| 欧美啪啪一区| 国产门事件在线| 中国特黄美女一级视频| 精品国产一区91在线| 国产老女人精品免费视频| 国产超碰在线观看| 国产精品极品美女自在线| 国产超碰在线观看| 欲色天天综合网| 国产精品入口麻豆| 国产精品久久精品| 国产精品久久国产精麻豆99网站| 亚洲精品视频免费看| 成人在线不卡| 国产成人精品一区二区三区| A级全黄试看30分钟小视频| 日本a级免费| 91最新精品视频发布页| 伊人91视频| a级免费视频| www中文字幕在线观看| 亚洲精品国产自在现线最新| 亚洲香蕉在线| 欧美国产日韩一区二区三区精品影视| 日韩欧美国产精品| 国产自在自线午夜精品视频| 国产成人精品亚洲日本对白优播| 日本不卡视频在线| 人妻无码一区二区视频| 在线视频一区二区三区不卡| 欧美丝袜高跟鞋一区二区| 一级毛片免费高清视频| 亚洲国产午夜精华无码福利| 伊人91在线| 日韩高清无码免费| 国产久草视频| 在线观看国产小视频| 午夜日韩久久影院| 超碰精品无码一区二区| 波多野结衣无码中文字幕在线观看一区二区 | 97无码免费人妻超级碰碰碰| 亚洲欧美综合另类图片小说区| 国模极品一区二区三区| 嫩草国产在线| 99热线精品大全在线观看| 亚洲69视频| 欧美在线中文字幕| 怡春院欧美一区二区三区免费| 伊人婷婷色香五月综合缴缴情| 亚洲国产精品无码AV| 国产香蕉97碰碰视频VA碰碰看| 成人亚洲天堂| 欧美亚洲激情| 欧美成人午夜视频免看| 欧美午夜在线观看| 日韩免费视频播播| 免费在线观看av| 午夜性刺激在线观看免费| 亚洲婷婷丁香| 日韩乱码免费一区二区三区| 亚洲日本精品一区二区| 欧美一区国产| 国产精品亚洲精品爽爽| 久久99精品久久久久纯品| 欧洲高清无码在线| 国产午夜精品鲁丝片| 伊人久久婷婷| 免费观看亚洲人成网站| 国产成人免费视频精品一区二区 | 亚洲性日韩精品一区二区| 中文字幕天无码久久精品视频免费| 精品国产福利在线| 亚洲精品无码抽插日韩| 亚洲国产日韩在线观看| 亚洲第七页| 国产精品精品视频| 国产呦精品一区二区三区下载| 亚洲欧美激情另类| 秋霞一区二区三区| 女人爽到高潮免费视频大全| 亚洲欧美激情另类|