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

基于ASTFA的廣義解調方法及應用

2015-12-30 05:20:14李寶慶,程軍圣,彭延峰
中國機械工程 2015年19期

基于ASTFA的廣義解調方法及應用

李寶慶程軍圣彭延峰楊宇

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

摘要:結合自適應最稀疏時頻分析(ASTFA)和廣義解調的優點提出了基于ASTFA的廣義解調方法。該方法首先采用ASTFA對原始信號進行分解,得到分量信號及其相位函數;然后,提取該相位函數的二次項及高次項,獲得解調相位函數;之后利用解調相位函數對分量信號進行廣義解調;最后對廣義解調后的信號進行頻域分析,提取特征信息。仿真和實驗分析結果表明,基于ASTFA的廣義解調方法非常適用于處理多分量頻率調制信號,能夠有效提取滾動軸承在變速工況下的故障特征信息。

關鍵詞:自適應最稀疏時頻分析;廣義解調;相位函數;頻率調制;滾動軸承

中圖分類號:TH113.1;TH911.7

收稿日期:2015-03-23

基金項目:國家自然科學基金資助項目(51375152)

作者簡介:李寶慶,男,1984年生。湖南大學機械與運載工程學院博士研究生。主要研究方向為機械設備故障診斷、動態信號分析與處理。發表論文2篇。程軍圣,男,1968年生。湖南大學機械與運載工程學院教授、博士研究生導師。彭延峰,男,1988年生。湖南大學機械與運載工程學院博士研究生。楊宇,女,1971年生。湖南大學機械與運載工程學院教授、博士研究生導師。

Generalized Demodulation Method Based on ASTFA Method and Its Applications

Li BaoqingCheng JunshengPeng YanfengYang Yu

State key Laboratory of Advanced Design and Manufacture for Vehicle Body,

Hunan University, Changsha, 410082

Abstract:Combining with the advantages of the ASTFA and the generalized demodulation, a generalized demodulation method was proposed based on the ASTFA. Firstly, the original signals were decomposed via ASTFA,and the component and its corresponding phase function were obtained. Secondly, the demodulation phase function was formatted by extracting the secondary and higher order terms of the phase function. Then, the generalized demodulation method was applied to the component using the demodulation phase demodulation function. Finally, the frequency domain analyses were applied to the demodulated component to extract feature information. Simulation and experimental analyses show that the method is very suitable for processing multi-component frequency modulated signals, and can effectively extract the damage characteristic informations of the rolling bearing damages in shifting conditions.

Key words: adaptive and sparsest time-frequency analysis(ASTFA); generalized demodulation; phase function; frequency modulation; rolling bearing

0引言

當滾動軸承出現故障時,會產生以軸承元件固有頻率為載波頻率的多頻率調制信號(表現為不斷重復的脈沖衰減信號)。對信號進行頻域分析,獲取脈沖出現的頻率成為軸承損傷判別的關鍵[1]。滾動軸承在升降速過程中的振動信號包含了豐富的狀態信息,一些在平穩運行時不易反映的特征信息可能會被充分地表現出來,對升降速過程中的振動信號進行分析具有重要意義。轉速變化時,軸承的特征信息表現為頻率曲線變化的非平穩信號,傅里葉變換等傳統的故障診斷方法不再適用[2]。工程上一般采用階次分析方法來實現非平穩信號的平穩化,它可以將變速過程中產生的與轉速有關的振動信號有效地分離出來[3]。但是,階次分析方法首先需要獲取轉速信號,采用傳感器等硬件獲取轉速需要較高的成本,而且在一些場合無法通過硬件獲取轉速[4]。有學者基于原始振動信號的瞬時頻率估計來提取轉速信號[4-5],但是通過這種方式得到的轉速信號準確度欠佳。另外,階次分析中的角域重采樣階次不容易確定。文獻[6]提出了一種濾波定階方法,但在選擇濾波器截止頻率時存在一定的主觀性。Olhede等[7]提出的廣義解調方法同樣能夠實現非平穩信號的平穩化,它把時頻分布是傾斜、非線性或曲線的非平穩信號變換為時頻分布是線性和平行于時間坐標軸的平穩信號。廣義解調方法的關鍵是解調相位函數的選擇[8-9]。目前都是基于先驗知識來估計解調相位函數,文獻[9-10]先采用最大離散小波包變換獲得原始信號的近似時頻分布,再估計解調相位函數。文獻[11]利用基于多尺度線調頻基的稀疏信號分解(multi-scale chirplet path pursuit,MCPP)方法先將頻率曲線變化的多分量信號分解為單分量信號,獲得各分量的瞬時頻率曲線,然后對各分量的瞬時頻率進行積分估計解調相位函數。無論是基于近似時頻分布,還是基于瞬時頻率的積分來估計解調相位函數都存在誤差,嚴重時會造成錯誤的廣義解調結果。

近期, Hou等[12-13]提出了一種自適應最稀疏時頻分析(adaptive and sparsest time-frequency analysis,ASTFA)方法。該方法在信號分解的同時可以直接獲得各分量的相位函數,無需先驗知識。提取該相位函數的二次項及高次項就可以獲得解調相位函數,利用該解調相位函數對分量信號進行廣義解調后,分量信號的時頻分布為平行于時間坐標軸的直線,實現了非平穩信號的平穩化。利用ASTFA方法獲得的解調相位函數具有更好的準確度。

本文結合ASTFA和廣義解調的優點,提出了基于ASTFA的廣義解調方法。首先使用ASTFA對原始信號進行分解,得到分量信號并獲得相位函數,然后提取該相位函數的二次項及高次項作為解調相位函數用于廣義解調,最后對廣義解調后的信號進行頻域分析提取特征信息。本文利用該方法對變速工況下存在故障的滾動軸承振動信號進行了分析,實際應用表明該方法能夠有效提取相應的特征信息。

1基于ASTFA的廣義解調方法

1.1廣義解調

廣義解調方法能把時頻分布是傾斜、非線性或曲線的非平穩信號變換為時頻分布是線性和平行于時間軸的平穩信號,本質上是廣義傅里葉變換(generalized Fourier transform,GFT)。信號x(t)的GFT定義為

(1)

其中,s0(t)為隨時間變化的實值函數。GFT實際上是對x(t)exp(-j2πs0(t))作標準傅里葉變換。設信號x(t)=exp(-j2π(f0t+s0(t))),如果找到近似s0(t)的相位函數s(t),對信號x(t)進行GFT:

exp(-j2π(ft+s(t)))dt=

(2)

1.2ASTFA方法

多尺度信號在時頻分布上具有稀疏分布的特點,基于此,ASTFA在建立包含所有自適應基的過完備字典庫的基礎上尋找信號的最稀疏表達。ASTFA首先建立合適的包含所有自適應基的過完備字典庫,然后在過完備字典庫中搜索對信號的匹配性最好的自適應基。

(1)過完備字典庫:

(3)

k=0,1,…,λn;l=1,2,…,λn}

(4)

(2)最優化問題:

(5)

上述優化問題是一個L0優化問題,即NP-Hard問題。根據逼近論,采用對L0范數不同的逼近方式可以將上述問題轉化為其他的優化問題,如L1或L2優化問題[14]。ASTFA方法將上述L0優化問題轉化為L2優化問題,并采用高斯牛頓迭代方法解決該問題。主要的計算步驟如下。

(1)令初始殘差r0=f(t);

(2)求解非線性最小二乘問題:

(6)

式中,ui(t)為分解得到的第i個內稟模態函數分量;ai(t)為分量的瞬時幅值;θi(t)為分量的相位函數。

(3)更新殘差:

ri=ri-1-ai(t)cosθi(t)

(7)

1.3廣義解調相位函數的選擇

為區分廣義解調的相位函數與ASTFA方法中的相位函數,本文將廣義解調的相位函數稱為解調相位函數,用s(t)表示;將ASTFA方法中的相位函數稱為相位函數,用θ(t)表示。

(8)

那么

(9)

式(9)表明解調相位函數s(t)為時間t的二次項及高次項之和。

按ASTFA方法原理及式(9)可得相位函數θ(t)與解調相位函數s(t)之間的關系:

s(t)=Θ(t)/(2π)

(10)

其中,Θ(t)為相位函數θ(t)的二次項及高次項之和。提取ASTFA獲得的相位函數θ(t)的二次項及高次項,再除以系數2π就可以獲得廣義解調的解調相位函數s(t),無需獲得原始信號的先驗時頻分布。

1.4基于ASTFA的廣義解調方法

基于ASTFA的廣義解調方法步驟如下:

(1)采用ASTFA方法對原始信號進行分析,同時獲得分量信號及其相位函數θi(t)。

(2)采用式(9)對相位函數θi(t)進行處理,獲得解調相位函數si(t)。

(3)使用解調相位函數si(t)對各分量信號進行廣義解調。

(4)對廣義解調后的信號進行頻域分析,可采用FFT變換產生幅值譜。

2仿真信號分析

為驗證基于自適應最稀疏時頻分析的廣義解調方法對多分量調頻信號分析的有效性,考察兩個幅值不等的調頻信號組成的仿真信號:

x(t)=x1(t)+x2(t)

(11)

x1(t)=3cos(120t3-80t2+500t)

x2(t)=cos(120t3-80t2+800t)

其中,t∈[0,1],單位:s。

仿真信號的時域波形如圖1所示,分解結果如圖2所示。從圖2可以看出,兩個調頻信號已經完全分離。

圖1 仿真信號x(t)的時域波形

圖2 仿真信號x(t)的ASTFA分解結果

在信號分解的過程中可以直接獲得相位函數θ(t),函數曲線如圖3所示。按照式(10)獲得解調相位函數s(t),得到的各分量的解調相位函數:

s1(t)=19.161t3-12.822t2

(12)

s2(t)=19.296t3-13.022t2

(13)

同時給出解調相位函數的理論值:

(14)

圖3 ASTFA方法得到的相位函數

方便起見,將各分量的解調相位函數與理論值曲線在同一坐標系下進行對比。從圖4可以看出,3條曲線基本完全重合,說明由ASTFA方法得到的解調相位函數完全正確。

圖4 解調相位函數曲線對比圖

利用解調相位函數s1(t)、s2(t)對分量信號進行廣義解調,將廣義解調后的信號疊加混合并進行頻譜分析,得到的幅值譜如圖5所示。式(11)給定的仿真信號x1(t)、x2(t)的載波頻率為79.6Hz和127.4Hz。圖5中峰值頻率為80Hz、128Hz,與載波頻率接近,說明基于ASTFA的廣義解調方法能有效地對多分量調頻信號進行解調。

圖5 廣義解調后的分量信號的幅值譜

為進一步說明基于ASTFA的廣義解調方法對頻率調制信號解調的有效性,我們對仿真信號的頻率進行分析。給出仿真信號的理論瞬時頻率曲線(圖6)、由ASTFA方法分解直接得到的瞬時頻率曲線(圖7)。通過頻率對比分析可知,由ASTFA方法獲得的瞬時頻率與理論值基本相同,表明ASTFA能有效分解多分量頻率調制信號。由廣義解調后的疊加混合信號的Hilbert譜(圖8)可知,經過廣義解調后的信號的時頻分布為平行于時間軸的直線,實現了信號的平穩化。以上分析表明,基于ASTFA的廣義解調方法可對多分量的頻率調制信號實現有效分解并將分量信號變為平穩信號。

圖6 仿真信號x(t)的理論瞬時頻率

圖7 ASTFA方法得到的瞬時頻率

3應用實例

3.1滾動軸承的故障特征

假設滾動體與內外圈之間為純滾動接觸,軸承外圈故障特征頻率為

(15)

內圈故障特征頻率為[15]

(16)

式中,f為軸的轉頻;Z為滾動體個數;d為滾動體直徑;D為軸承節徑;α為接觸角。

從式(15)、式(16)可以看出,特征頻率與軸的轉頻間存在固定的比例關系,比例由軸承本身的特性決定。當轉速非平穩時,滾動軸承的故障特征頻率與軸的轉頻都是非平穩的。因此,當轉速非平穩時,可以采用廣義解調方法對信號進行平穩化處理。但是,當滾動軸承存在故障時,其振動信號還表現為多頻率調制特性,而且在廣義解調之前還需要確定解調相位函數。因此,基于ASTFA的廣義解調方法適用于變轉速工況下的滾動軸承非平穩信號處理。

3.2實驗信號分析

如圖9所示,在實驗臺上進行具有外圈故障、內圈故障的滾動軸承瞬態實驗。軸承為6307型,利用激光分別在外圈和內圈上切割寬0.15mm、深0.13mm的槽來模擬故障。在軸承座上采集振動加速度信號,為了便于觀察與驗證,采用光電式轉速傳感器采集轉速信號。

3.2.1外圈故障分析

外圈故障實驗的主軸轉速如圖10所示,轉速由809r/min升為1281r/min,轉速非平穩。振動信號如圖11所示,圖中存在明顯的沖擊脈沖,但無法判定是否為外圈故障沖擊。

圖10 外圈故障主軸轉速信號

圖11 滾動軸承外圈故障振動信號

對圖10所示的轉速信號進行廣義解調,可以得到廣義解調后的轉頻13.07Hz,根據式(15)得到滾動軸承的外圈故障特征頻率39.99Hz。采用基于ASTFA的廣義解調方法對圖11所示的振動加速度信號進行分析,得到的分量信號如圖12所示,得到的相位函數如圖13所示,根據式(9)計算得到的解調相位函數如圖14所示。采用圖14所示的解調相位函數s1(t)、s2(t)、s3(t)對相對應的3個分量信號u1(t)、u2(t)、u3(t)進行廣義解調,并對廣義解調后的合成信號進行FFT分析得到頻譜,如圖15所示。從圖15可以看出,在滾動軸承外圈故障特征頻率及其2倍頻處具有明顯的峰值,表明基于ASTFA的廣義解調方法能夠有效提取滾動軸承外圈故障特征。

圖12 滾動軸承外圈故障振動信號ASTFA分解結果

圖13 ASTFA方法得到的外圈故障相位函數

圖14 外圈故障解調相位函數

圖15 外圈故障分量信號廣義解調后的頻譜(ASTFA)

圖16 基于MCPP的廣義解調方法的外圈故障頻譜

采用文獻[11]基于MCPP的廣義解調方法對圖11所示的信號進行分析,結果如圖16所示。對比圖15、圖16可知,基于ASTFA的廣義解調方法得到的故障特征頻率及其2倍頻更接近于理論值,準確度更高。

3.2.2內圈故障分析

內圈故障實驗的主軸轉速如圖17所示,轉速由714r/min降為532r/min。振動加速信號如圖18所示。對圖17所示的轉速信號進行廣義解調,可以得到廣義解調后的轉頻11.9Hz,根據式(16)得到滾動軸承的內圈故障特征頻率58.8Hz。采用基于ASTFA的廣義解調方法對圖18所示的振動加速度信號進行分析,得到的分量信號如圖19所示,得到的相位函數如圖20所示,根據式(10)計算得到的解調相位函數曲線如圖21所示。

圖17 內圈故障主軸轉速信號

圖18 滾動軸承內圈故障振動信號

圖19 滾動軸承內圈故障振動信號ASTFA分解結果

圖20 ASTFA方法得到的內圈故障相位函數

圖21 內圈故障解調相位函數

圖22 內圈故障分量信號廣義解調后的頻譜(ASTFA)

采用圖21所示的解調相位函數s1(t)、s2(t)、s3(t)對相對應的三個分量信號u1(t)、u2(t)、u3(t)進行廣義解調,并對廣義解調后的合成信號進行FFT分析得到頻譜,如圖22所示。從圖22可以看出,在滾動軸承內圈故障特征頻率處具有明顯的峰值,表明基于ASTFA的廣義解調方法能夠有效提取滾動軸承內圈故障特征。

類似于外圈故障,同樣采用文獻[11]基于MCPP的廣義解調方法對圖18所示的信號進行分析,結果如圖23所示。從圖23可以看出,圖中存在明顯的故障特征峰值頻率55Hz,但是圖示的故障特征頻率與理論值相比偏差較大。因此,基于ASTFA的廣義解調方法得到的故障特征頻率更接近于理論值,準確度更高。

圖23 基于MCPP的廣義解調方法的內圈故障頻譜

對比圖22、圖23還可以看出,圖23中存在明顯的2倍頻及近3倍頻的峰值,但是圖22中并無明顯的轉頻峰值。出現這種情況的主要的原因是,ASTFA方法基于最優化方法實現信號的分解,當滿足收斂條件時才能夠實現信號的分離。對于滾動軸承內圈故障來講,滾動軸承故障振動信號與主軸旋轉振動信號混雜在一起。采用ASTFA方法進行分解時,轉頻信號與故障信號的收斂條件并不一致,因此ASTFA可以準確地提取到滾動軸承的故障振動信號。多尺度線調頻基稀疏信號分解是基于多尺度線調頻原子的匹配追蹤方法,本質上是利用多個基函數對信號進行分段表示,從而實現信號的分離。因此,其在對信號進行分段描述時無法將比較明顯的轉頻信號完全分離。

4結語

ASTFA方法作為一種新的自適應時頻分析方法能有效地對多分量頻率調制信號進行解調,并且在實現信號分解的過程中直接獲得相位函數,提取其二次項及高次項就可得到廣義解調方法所需的解調相位函數。

相比于其他的解調相位函數估計方法,由ASTFA得到的解調相位函數更加準確,更徹底地解決了廣義解調時頻分析方法中相位函數的選擇問題,提高了廣義解調方法的應用性。實驗分析結果表明,基于ASTFA的廣義解調方法能夠有效提取變速工況下的滾動軸承故障特征信息,具有一定的應用價值。

但是,ASTFA方法在計算效率上存在一定的不足,由于其基于最優化方法實現信號的分解,因此計算效率較低。若原始信號為N維向量,則ASTFA需要同時對N個參數進行優化。當原始信號維數較大時,使用ASTFA進行信號分解會產生很大的計算量。在后續的研究中,需在ASTFA方法的快速算法上進行重點研究。

參考文獻:

[1]王冬云,張文志.基于小波變換的滾動軸承故障診斷[J].中國機械工程,2012,23(3):295-298.

WangDongyun,ZhangWenzhi.FaultDiagnosisStudyofBallBearingBasedonWaveletPacketTransform[J].ChinaMechanicalEngineering,2012,23(3):295-298.

[2]徐亞軍,于德介,孫云嵩,等.滾動軸承故障診斷的階比多尺度形態學解調方法[J]. 振動工程學報,2013,26(2):252-259.

XuYajun,YuDejie,SunYunsong,etal.RollerBearingFaultDiagnosisUsingOrderMulti-scaleMorphologyDemodulation[J].JournalofVibrationEngineering,2013,26(2):252-259.

[3]康海英,欒軍英,鄭海起,等.基于階次跟蹤和經驗模態分解的滾動軸承包絡解調分析[J].機械工程學報,2007,43(8):119-122.

KangYajun,LuanJunying,ZhengHaiqi,etal.EnvelopeDemodulationAnalysisofBearingBasedonOrderTrackingandEmpiricalModeDecomposition[J].JournalofMechanicalEngineering,2007,43(8):119-122.

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

GuoYu,QinShuren,TangBaoping,etal.OrderTrackingofRotatingMachineryBasedonInstantaneousFrequencyEstimation[J].JournalofMechanicalEngineering,2003,39(3):32-36.

[5]李蓉,于德介,陳向民,等.基于階次分析與循環平穩解調的齒輪箱復合故障診斷方法[J]. 中國機械工程,2013,24(5):1320-1327.

LiRong,YuDejie,ChengXiangmin,etal.ACompoundFaultDiagnosisMethodforGearboxBasedonOrderTrackingandCyclostationaryDemodulation[J].ChinaMechanicalEngineering,2013,24(5):1320-1327.

[6]叢華,吳廣平,饒國強,等.計算階次分析中避免階次混疊的濾波定階方法及其應用[J]. 振動與沖擊,2012,31(12):42-44.

CongHua,WuGuangping,RaoGuoqiang,etal.AMethodtoAvoidOrderAliasinginCOTBasedonFiltering[J].JournalofVibrationandShock,2012,31(12):42-44.

[7]OlhedeS,WaldenAT.TheHilbertSpectrumviaWaveletProjections[J].Proc.R.Soc.London.A,2004,460(2044):955-975.

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

ChengJunsheng,YangYu,YuDejie.AMulti-componentSignalDecompositionMethodBasedontheGeneralizedDemodulationTime-frequencyAnalysis[J].JournalofVibrationEngineering,2007,20(6) :563-569.

[9]楊宇,程軍圣,于德介. 廣義解調時頻分析方法中的若干問題探討[J]. 振動與沖擊,2008,27(2):19-24.

YangYu,ChengJunsheng,YuDejie.StudyonSomeProblemsintheGeneralizedDemodulationTime-frequencyAnalysisMethod[J].JournalofVibrationandShock,2008,27(2):19-24.

[10]張曉菲,劉振興,陳棟.基于廣義解調時頻分析的多分量復雜信號分解方法[J].數據采集與處理,2012,27(5) :630-634.

ZhangXiaofei,LiuZhenxing,ChenDong.Multi-componentandComplicatedSignalDecompositionMethodBasedonGeneralizedDemodulationTime-frequencyAnalysis[J].JournalofDataAcquisition&Processing,2012,27(5) :630-634.

[11]任凌志,于德介,彭富強.基于多尺度線調頻基稀疏信號分解的廣義解調方法及其在滾動軸承故障診斷中的應用[J]. 中國電機工程學報,2010,30(11):102-108.

RenLingzhi,YuDejie,PengFuqiang.GeneralizedDemodulationMethodBasedonMulti-scaleChirpletandSparseSignalDecompositionandItsApplicationtoRollerBearingFaultDiagnosis[J].ProceedingsoftheCSEE,2010,30(11):102-108.

[12]HouTY,SHIZQ.AdaptiveDataAnalysisviaSparseTime-frequencyRepresentation[J].AdvancesinAdaptiveDataAnalysis,2011,3(1/2):1-28.

[13]HouTY,ShiZQ.Data-drivenTime-frequencyAnalysis[J].AppliedandComputationalHarmonicAnalysis,2013,35(2):284-308.

[14]MascarenasD,CattaneoA,TheilereJ,etal.CompressedSensingTechniquesforDetectingDamageinStructures[J].StructuralHealthMonitoring,2013,12(4):325-338.

[15]鐘秉林,黃仁.機械故障診斷學[M].北京:機械工業出版社,2006.

(編輯張洋)

主站蜘蛛池模板: 欧美不卡视频在线| 在线免费观看AV| 国产成熟女人性满足视频| 成人福利在线观看| 精品剧情v国产在线观看| 亚洲美女视频一区| 亚洲天堂网2014| 日韩色图在线观看| 精品无码人妻一区二区| 成人免费午间影院在线观看| 亚洲午夜福利在线| www.亚洲国产| 亚洲伊人久久精品影院| 日韩免费毛片| 四虎在线高清无码| 久久精品91麻豆| 亚洲视频黄| 无码 在线 在线| 日韩毛片免费视频| 九九热这里只有国产精品| 亚洲综合亚洲国产尤物| 91久久国产热精品免费| 免费国产不卡午夜福在线观看| 日韩久久精品无码aV| 国产97公开成人免费视频| 国产熟女一级毛片| 色九九视频| 91福利国产成人精品导航| 成人福利在线观看| 亚洲成人77777| 亚洲性日韩精品一区二区| 国产精品综合色区在线观看| 国产9191精品免费观看| 亚洲精品成人福利在线电影| 亚洲乱强伦| 免费看黄片一区二区三区| 亚洲啪啪网| 亚洲欧美日韩动漫| 久久综合亚洲鲁鲁九月天| 亚洲精品在线影院| 国产精品无码在线看| 无码一区18禁| 99无码熟妇丰满人妻啪啪| 亚洲精品无码久久毛片波多野吉| 国产精品一区二区在线播放| 国产成人精品免费视频大全五级 | 中文字幕66页| 国产精品久久久久久久伊一| 美女无遮挡被啪啪到高潮免费| 一级毛片免费的| 97国产精品视频自在拍| 国产91全国探花系列在线播放| 2022精品国偷自产免费观看| 久久人妻xunleige无码| 国产福利在线观看精品| 久久亚洲美女精品国产精品| 日本免费福利视频| 精品视频一区在线观看| 亚洲高清在线天堂精品| 免费不卡视频| 国产特级毛片| 青青草国产在线视频| 狠狠色狠狠色综合久久第一次| 亚洲国产天堂久久综合226114| 亚洲最大福利网站| 漂亮人妻被中出中文字幕久久| 性色一区| 亚洲欧美另类中文字幕| 中国国产一级毛片| 国产女主播一区| 日本精品影院| 欧美午夜一区| 日本午夜三级| 成人亚洲视频| 亚洲欧洲日韩久久狠狠爱 | 欧美综合区自拍亚洲综合绿色| 69国产精品视频免费| 中文字幕第4页| 亚洲av无码人妻| 免费在线国产一区二区三区精品| 亚洲国产综合精品中文第一| 亚洲色图欧美视频|