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

分數階時頻局部二值模式譜在齒輪故障診斷中的應用

2016-01-15 02:23:57張云強張培林
振動與沖擊 2015年11期
關鍵詞:特征提取故障診斷

張云強,張培林,李 兵

(軍械工程學院車輛與電氣工程系,石家莊 050003)

第一作者張云強男,博士生,1987年生

分數階時頻局部二值模式譜在齒輪故障診斷中的應用

張云強,張培林,李兵

(軍械工程學院車輛與電氣工程系,石家莊050003)

摘要:針對齒輪故障信號分析,提出采用分數階時頻LBP譜表達齒輪故障信號的時頻特征。首先為克服S變換對高頻信號的時頻分辨性能差的不足,基于分數階Fourier變換良好的時頻旋轉特性設計了一種分數階S變換,用于獲取齒輪信號的二維時頻表示;然后引入局部二值模式(LBP)算子,將LBP算子作用于分數階S變換時頻圖,提取分數階時頻LBP譜;最后結合“uniform”模式LBP的概念和類內類間距準則,對分數階時頻局部二值模式譜進行特征優選,用于表達齒輪故障特征。對5種不同狀態的齒輪信號進行了分析,結果表明優選后的分數階時頻LBP譜具有較強的特征描述能力,是齒輪故障信號的一類新的有效特征參數。

關鍵詞:故障診斷;齒輪;特征提取;分數階S變換;局部二值模式

基金項目:國家自然科學基金資助項目(E51205405,51305454)

收稿日期:2014-03-05修改稿收到日期:2014-06-06

中圖分類號:TH165.3;TH113

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2015.11.022

Abstract:A novel scheme utilizing fractional order time-frequency local binary pattern(LBP) spectra was proposed to express time-frequency characteristics of gear fault signals. Since the time-frequency resolution of S transformation was low to high-frequency signals, the fractional Fourier transformation with excellent time-frequency rotating character was integrated into S transformation and as a result, a fractional S transformation was designed. The fractional S transformation was employed to obtain 2D time-frequency representations of gear fault signals. Then, LBP operator was introduced and imposed on time-frequency images of the fractional S transformation to extract fractional order time-frequency LBP spectram. With the concept of uniform pattern of LBP and the separability evaluation criterion, the fractional order time-frequency LBP spectra were optimally selected for succinct description of gear fault features. Gear vibration signals under five different states were analyzed. The results indicated that the optimally selected fractional order time-frequency LBP spectra have an outstanding feature description capacity; they are a kind of new and effective feature parameters for gear fault signals.

Application of fractional order time-frequency local binary pattern spectra in gear fault diagnosis

ZHANGYun-qiang,ZHANGPei-lin,LIBing(Department of Vehicles and Electrical Engineering, Ordnance Engineering College, Shijiazhuang 050003, China)

Key words:fault diagnosis; gear; feature extraction; fractional order S transformation; local binary pattern (LBP)

齒輪是在機械設備中廣泛應用的傳動部件之一,當發生故障時,其振動信號屬于典型的非線性、非平穩信號[1]。在非平穩信號處理領域,時頻分析技術采用時域和頻域的二維聯合表示,可以對非平穩信號的局部特性精確描述,具有時域和頻域等傳統技術無法比擬的優勢[2-3]。

常用的時頻分析技術有小波變換[4]、S變換[5]和經驗模式分解[6]等。這些技術雖然在機械信號處理中具有廣泛應用,但是都存在一些不足。在小波分析中,母小波的選擇對信號分析結果影響很大,而選擇合適的母小波并非易事,并且小波尺度與信號頻率沒有良好的對應關系。在短時Fourier變換和小波變換的基礎上提出的S變換雖然克服了小波變換的一些不足,但是其高斯窗函數的標準差固定為頻率的倒數,導致對高頻信號的時頻分辨力較差。經驗模式分解對噪聲敏感,存在模式混疊現象,容易產生虛假分量。

分數階Fourier變換作為Fourier變換的推廣,是一種統一的時頻分析方法,可以看成時頻面內的旋轉因子,在不同階次的分數階頻域具有不同的時頻聚集性[7-8]。為了提高信號時頻分析的靈活性和S變換的時頻分辨性能,本文基于分數階Fourier變換良好的時頻旋轉特性,提出一種新的時頻分析技術,分數階S變換,并應用于齒輪故障信號處理。

分數階S變換二維時頻圖雖能有效描述信號的時頻局部特性,但其維數巨大,不能直接作為信號的時頻特征參數,還需要進一步提取低維特征。鑒于時頻圖本質上是一種圖像,不同狀態齒輪信號的時頻圖會呈現出不同的紋理,引入常用于描述圖像紋理特征的局部二值模式[9-10]理論,在分數階S變換對齒輪信號處理的基礎上,提出一種分數階時頻局部二值模式譜,并結合“uniform”局部二值模式[11]概念和類內類間距準則[12],對分數階時頻局部二值模式譜進行特征優選,用于表達齒輪故障信號特征。

1分數階S變換

1.1S變換和分數階Fourier變換

S變換結合了短時Fourier變換和小波變換的優點,是一種較新的時頻分析方法。對于非平穩時間信號x(t),其S變換如下[13]

ST(τ,f)=

(1)

式中:w(t)為高斯窗函數,表達式為

(2)

由式(2)可知,S變換的高斯窗函數的標準差σ為頻率f的倒數,即σ=1/|f|,其時窗寬度隨著頻率f的增大而減小。由此可知,S變換對低頻信號具有較高的頻率分辨力,而對高頻信號具有較高的時間分辨力。

分數階Fourier變換是Fourier變換的廣義形式,能將時域信號變換到分數階頻域,并且得到信號新的表示。信號x(t)的分數階Fourier變換定義為

(3)

式中:Ka(t,u)為分數階Fourier變換核。

(4)

分數階Fourier變換可以看成二維時頻面內的旋轉因子,具有一定的時頻旋轉特性。a階分數階Fourier變換可以將信號從時間軸逆時針旋轉任意角度α到分數階頻率軸,并且獲得其在分數階頻域中新的表示。當a=1時,分數階Fourier變換即為傳統的Fourier變換。

1.2分數階S變換

S變換雖然具有比小波變換更好的時頻特性,但是對高頻信號的時頻分辨力不理想。為了提高S變換的時頻分辨力,本文利用分數階Fourier變換良好的時頻旋轉特性,定義分數階S變換如下

FrST(τ,u)=

(5)

式中:高斯窗函數w(t,u)是時間t和分數階頻率u的函數,如式(6)所示。

(6)

式中:p為窗函數調整參數,p∈(0,1]。

由式(5)和式(6)可知,隨著a和p取值的變化,分數階S變換能將信號變換到不同的分數階頻域分析,具有更大的靈活性。同時,分數階S變換的時窗寬度隨著分數階頻率的增大而減小,所以它對分數階低頻信號具有較高的頻率分辨力,而對分數階高頻信號具有較高的時間分辨率,繼承了S變換良好的時頻特性。通過引入調整參數p,可以有效調節信號頻率對時窗寬度的影響,使低頻和高頻信號的分辨力達到較好的折中,從而提高信號的整體時頻分辨性能。需要注意的是,當a=1且p=1時,分數階頻域即傳統的Fourier頻域,此時分數階S變換退化為S變換。

1.3分數階S變換快速算法

為了快速實現分數階S變換,根據卷積的定義,首先將分數階S變換的定義式改寫為

FrST(τ,u)=[x(t)Ka(t,u)]*w(t,u)

(7)

令F(x(t))和F-1(x(t))分別表示x(t)的Fourier變換及其逆變換,根據時域卷積定理有

F[FrST(τ,u)]=F[x(t)Ka(t,u)]F[w(t,u)]

(8)

對上式兩邊同時進行Fourier逆變換得

FrST(τ,u)=

F-1{F[x(t)Ka(t,u)]F[w(t,u)]}

(9)

由式(9)可知,信號x(t)的分數階S變換可以通過Fourier變換及其逆變換得到。由于Fourier變換及其逆變換已有成熟的快速算法,因此可以借助它們實現分數階S變換的快速計算,具體步驟不再贅述。

1.4仿真信號分析

為了驗證分數階S變換的有效性,構造仿真信號x(t)如下

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

x1(t)=sin(30πt)

x2(t)=sin(30+50πt·4t)

x3(t)=sin[200πt+200πt2+20πtcos(2πt)]

(10)

由式(10)可知,x(t)由1個正弦分量和2個非線性調頻分量組成。設置采樣頻率和采樣時間分別為1 024 Hz和1 s,其時域波形見圖1。

圖1 仿真信號Fig.1 The simulated signal

對仿真信號x(t)分別進行短時Fourier變換、Morlet連續小波變換、S變換和分數階S變換(見圖2)。

對比圖2中子圖可以看出,由于短時Fourier變換的時窗寬度固定,其低頻和高頻的時頻分辨率都很差;連續小波變換由于尺度的大小與信號的頻率沒有較好的對應關系,導致時頻圖的頻率分辨力較差,不能很好地表達信號各分量的頻率隨時間的變化情況;S變換雖然對低頻信號具有較好的時頻分辨力,但在高頻的時頻分辨力較差;分數階S變換由于具有類似于分數階Fourier變換的時頻旋轉特性,并引入了窗函數調整參數,時頻分辨性能明顯提高,尤其是在信號高頻部分。因此,分數階S變換具有更好的時頻特性,本文選其獲取齒輪信號的二維時頻圖。

圖2 仿真信號二維時頻圖Fig.2 2D time-frequency images of the simulated signal

2分數階時頻局部二值模式譜特征

2.1局部二值模式算子

局部二值模式是由Ojala等為了度量圖像局部對比度而提出的,其基本思想是根據圖像局部區域的中心像素與鄰域像素的灰度差異進行二進制編碼,從而對圖像紋理進行描述[14]。

LBP算子最初定義在3×3的矩形鄰域上,編碼過程見圖3,編碼規則為:當fi≥fc時,對應位置編碼為1,否則編碼為零,然后按順時針方向構成8位二進制數,即為中心像素fc的LBP。由于定義在3×3矩形鄰域上的LBP存在不能提取大尺寸紋理特征的局限性,Ojala等又將LBP算子的鄰域從3×3的矩形擴展到了任意半徑R和任意鄰域點數P的圓環形,對應的算子記為LBPP,R。圖4為三種常見的圓環形鄰域。

圖3 LBP編碼過程Fig.3 The coding process of LBP

圖4 三種常見的圓環形鄰域Fig.4 Three circular neighborhoods for LBP

LBPP,R的計算公式如下[11]

(11)

(12)

2.2分數階時頻局部二值模式譜

齒輪信號的時頻圖經過局部二值模式算子處理后仍然是一個維數巨大的矩陣。為了有效描述齒輪信號時頻圖的紋理特征,在分數階S變換時頻圖經局部二值模式算子處理的基礎上,定義分數階時頻局部二值模式譜H(h),數學描述如下

(13)

式中,n=2P;h=0,1,…,n-1;I是一個語句判斷函數,當LBPP,R(x,y)=h時,I{LBPP,R(x,y)=h}=1,否則I{LBPP,R(x,y)=h}=0;M、N分別為時頻圖矩陣的行數和列數。

2.3分數階局部二值模式譜特征優選

通過研究發現,LBP只有少部分模式屬于描述圖像紋理的重要模式,其出現概率達到90%以上,這些模式稱為“uniform”模式。如果把二進制串連成一個圓,“uniform”模式就是串中從0到1和1到0的轉換次數不超過2次的模式。以P=8,R=1為例,由式(11)可知LBPP,R的取值可達256種之多,而“uniform”模式LBPP,R的模式數目最多只有59種。雖然“uniform”模式僅占LBP的小部分,但是它反映了絕大部分紋理信息。另外, “uniform”模式LBP譜中不同參數的可分性差異很大,因此對齒輪信號特征描述和分類的貢獻差異也很大。

鑒于特征優選的目的是尋求類間散度最大而類內散度最小,本文首先利用“uniform”模式的概念對分數階時頻局部二值模式譜進行初選,即將分數階時頻局部二值模式譜中非“uniform”模式對應的數值舍棄,僅保留“uniform”模式對應的數值。在此基礎上,進一步通過類內類間判據對特征進一步優選,其定量評價指標為

(14)

式中:Sb和Sw分別為類內散度和類間散度。

3齒輪故障信號分析

實驗齒輪故障信號選自某型二級傳動齒輪箱振動試驗。該試驗系統中的兩對齒輪副的齒數分別為25/50和18/91,輸入軸由電機帶動,其轉速為1 491 r/min,輸出軸連接一個磁性阻尼器作為負載。試驗中模擬了齒輪在正常、中間軸齒根裂紋和齒面磨損、輸出軸齒根裂紋和齒面磨損5種狀態。選擇的采樣頻率和采樣點數為6 400 Hz和1 024。圖5為5種不同狀態的齒輪故障信號。

圖5 5種不同狀態的齒輪故障信號Fig.5 Gear fault signals from five different states

3.1齒輪故障信號的分數階時頻表示

采用分數階S變換分別對上述5種不同狀態下的齒輪故障信號進行處理,獲取其二維時頻圖。對于同一種狀態下的齒輪信號,設定參數a和p的取值相同,并通過多次實驗進行確定。最終,5種信號的時頻表示見圖6。

圖6 齒輪信號的分數階S變換時頻圖Fig.6 Time-frequency images of gear signals from fractional S transform

由圖6可知,不同狀態的齒輪信號的時頻圖具有明顯差異,呈現出不同的紋理特征。因此,齒輪故障信號的分數階S變換二維時頻圖就紋理角度而言具有可分性。

3.2分數階局部二值模式譜特征提取

在分數階S變換的基礎上,分別采用圖4中三種常見的環形鄰域提取分數階局部二值模式譜。

圖7為不同信號的分數階時頻LBP8,1譜、LBP8,2譜和LBP16,2譜,其中每種信號包含5個樣本,維數分別為256、256和65 536。基于“uniform”模式的概念對圖7中分數階局部二值模式譜特征進行初選的結果見圖8,其維數分別降為59、59和243。

對比圖7和圖8可知,由于“uniform”模式是局部二值模式中少量且重要的模式,初選后的“uniform”模式分數階時頻LBP譜特征維數大大降低,并且包含了原分數階LBP譜的主要信息。

根據類內類間判據準則易知,當式(14)中定量評價指標較大時,對應的特征具有較好的可分性,但是評價指標閾值不能選擇太大,太大會造成選出的特征數量太少,不能很好地描述信號。因此,通過多次實驗,最終選擇定量評價指標閾值為0.9,對上述初選后的分數階時頻LBP譜特征進一步優選。優選結果見圖9,特征維數分別為13、15和97。由圖9中優選結果可以看出,雖然“uniform”模式是局部二值模式中比較重要的模式,但是就可分性而言,它們仍存在明顯差異,表現出不同的類內聚合性和類間分散性。

圖7~圖9中H(h)和h的含義同式(13),且均無量綱。

綜上可知本文特征優選方法具有合理性和有效性。

圖7 齒輪信號的分數階時頻局部二值模式譜Fig.7 Fractional time-frequency LBP spectrums of gear signals

圖8 “uniform”模式分數階時頻局部二值模式譜Fig.8 Uniform pattern fractional time-frequency LBP spectrums of gear signals

圖9 優選后的分數階時頻局部二值模式譜Fig.9 Optimally selected fractional time-frequency LBP spectrums of gear signals

3.3齒輪故障診斷

將分數階時頻局部二值模式譜應用于齒輪故障診斷。從每類狀態的齒輪故障信號中分別隨機抽取40個樣本,其中20個組成訓練樣本,其余20個組成測試樣本,選擇樸素貝葉斯分類器進行識別。同時,引入基于S變換的時頻局部二值模式譜(時頻局部二值模式譜)作為對比實驗。實驗一共重復10次,平均結果見圖10。

由圖10柱狀圖可知,由于分數階S變換比S變換更靈活、時頻分辨性能更好,表達齒輪故障信號時頻特性更加精確,基于分數階S變換提取的特征分類效果優于基于S變換提取的特征,并且隨著特征優選過程的進行,這種優勢更加明顯。

無論基于S變換還是分數階S變換,所提取的局部二值模式譜特征在未進行特征優選前的分類精度較低,說明原始局部二值模式譜中包含的大量不重要信息不但對齒輪故障診斷沒有貢獻,反而嚴重影響了診斷精度。通過“uniform”模式概念優選后的特征由于舍棄了不重要模式,特征維數大幅度降低,上述不重要信息對故障診斷的影響減小,因此診斷精度有所提高。經過兩次優選后的局部二值模式譜特征通過類內類間判據準則篩選出的特征參數具有較大的類間散度和較小的類內散度,使對齒輪故障信號的描述能力和區分能力增強,因而樸素貝葉斯分類器獲得了更好的分類精度,進一步提高了齒輪故障診斷精度。另外需要指出的是,隨著特征優選過程的深入,特征維數迅速減少,故優選后的特征具有更快的診斷速度。

總之,分數階S變換改善了S變換的時頻特性,分析齒輪故障信號的靈活性更大、精度更高。結合“uniform”局部二值模的概念和類內類間距準則優選后的分數階時頻局部二值模式譜具有特征維數低、描述能力強和鑒別性能好的優點,適用于齒輪故障診斷。

圖10 齒輪故障診斷精度Fig.10 Gear fault diagnosis accuracy

4結論

針對齒輪故障信號時頻特征提取問題,本文借鑒分數階Fourier變換的時頻旋轉特性,改進了S變換。通過在S變換中引入分數階Fourier變換核和增加窗函數調整參數,提出了一種分數階S變換。仿真信號分析表明,分數階S變換具有更大的靈活性和更好的時頻分辨力。基于分數階S變換,定義了分數階時頻局部二值模式譜,并結合“uniform”局部二值模式的概念和類內類間距準則對分數階時頻局部二值模式譜特征優選方法進行了研究。通過齒輪故障信號分析,驗證了分數階時頻局部二值模式譜用于表達齒輪故障信號時頻特征的可行性和有效性。與優選前相比,優選后的分數階時頻局部二值模式譜具有特征維數低、可分性好的優點,可以有效提高齒輪故障診斷精度。

參考文獻

[1]李兵, 米雙山, 劉鵬遠, 等. 二維非負矩陣分解在齒輪故障診斷中的應用[J]. 振動、測試與診斷, 2012, 32(5): 836-840.

LI Bing, MI Shuang-shan, LIU Peng-yuan, et al.Application of two-dimensional non-negative factorization for gear fault diagnosis[J]. Journal of Vibration, Measurement & Diagnosis, 2012, 32(5): 836-840.

[2]姜鳴, 陳進, 汪慰軍. 幾種Cohen 類時頻分布的比較及應用 [J]. 機械工程學報, 2003, 39(8): 129-134.

JIANG Ming, CHEN Jin, WANG Wei-jun. Comparison and application of some time-frequency distributions belonging to Cohen class[J]. Chinese Journal of Mechanical Engineering, 2003, 39(8): 129-134.

[3]Meltzer G, Ivanov Y Y. Fault detection in gear drives with non-stationary rotational speed—part I: the time-frequency approach[J]. Mechanical Systems and Signal Processing, 2003, 17(5): 1033-1047.

[4]馬倫, 康建設, 孟妍, 等. 基于Morlet小波變換的滾動軸承早期故障特征提取研究[J]. 儀器儀表學報, 2013, 34(4): 920-926.

MA Lun, KANG Jian-she, MENG Yan, et al. Research on feature extraction of rolling bearing incipient fault based on Morlet wavelet transform[J]. Chinese Journal of Scientific Instrument, 2013, 34(4): 920-926.

[5]趙淑紅, 朱光明. S變換時頻濾波去噪方法[J]. 石油地球物理勘探, 2007, 42(4): 402-408.

ZHAO Shu-hong, ZHU Guang-ming. Time-frequency filtering to denoise by S transform[J]. Oil Geophysical Prospecting, 2007, 42(4): 402-408.

[6]邵忍平, 曹精明, 李永龍. 基于EMD 小波閾值去噪和時頻分析的齒輪故障模式識別與診斷[J]. 振動與沖擊, 2012, 31(8): 96-101.

SHAO Ren-ping, CAO Jing-ming, LI Yong-long. Gear fault pattern identification and diagnosis using time-frequency analysis and wavelet threshold de-noising based on EMD[J]. Journal of Vibration and Shock, 2012, 31(8): 96-101.

[7]陶然, 鄧兵, 王越. 分數階傅里葉變換及其應用[M]. 北京:清華大學出版社, 2009.

[8]Ervin S, Igor D, Ljubisa S. Fractional fourier transform as a signal processing tool: An overview of recent development[J]. Signal Processing, 2011, 91: 1351-1369.

[9]Loris N, Alessandra L, Sheryl B. Survey on LBP based texture descriptors for image classification[J]. Expert Systems with Applications, 2012, 39: 3634-3641.

[10]SONG Tie-cheng, LI Hong-liang. WaveLBP based hierarchical features for image classification[J]. Pattern Recognition Letters 2013, 34: 1323-1328.

[11]Ojala T, Pietikainen M, Maeenpaa T. Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(7): 971-987.

[12]周曉英, 李巍華, 丁康.基于關聯向量機的齒輪故障檢測方法研究[J]. 振動與沖擊, 2008, 27(6): 51-54.

ZHOU Xiao-ying, LI Wei-hua, DING Kang. Gear fault detection based on relevance vector machine[J]. Journal of Vibration and Shock, 2008, 27(6): 51-54.

[13]Djurovi I, Sejdi E, Jiang J. Frequency-based window width optimization for S-transform[J]. AEU-International Journal of Electronics and Communications, 2008, 62(4): 245-250.

[14]白航, 趙擁軍, 胡德秀. 時頻圖像局部二值模式特征在雷達信號分類識別中的應用[J]. 宇航學報, 2013, 34(1): 139-146.

BAI Hang, ZHAO Yong-jun, HU De-xiu. Radar signal recognition based on the local binary pattern feature of time-frequency image[J]. Journal of Astronautics, 2013, 34(1): 139-146.

猜你喜歡
特征提取故障診斷
特征提取和最小二乘支持向量機的水下目標識別
凍干機常見故障診斷與維修
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
基于Daubechies(dbN)的飛行器音頻特征提取
電子制作(2018年19期)2018-11-14 02:37:08
Bagging RCSP腦電特征提取算法
基于量子萬有引力搜索的SVM自駕故障診斷
因果圖定性分析法及其在故障診斷中的應用
基于MED和循環域解調的多故障特征提取
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 特级欧美视频aaaaaa| 女人毛片a级大学毛片免费| 潮喷在线无码白浆| 欧美日本在线观看| 亚洲男人天堂网址| 国产精品福利社| 亚洲人成色在线观看| 亚洲看片网| 中文字幕无码电影| 国产91成人| 伊伊人成亚洲综合人网7777 | 亚洲色图另类| www.99在线观看| 狠狠ⅴ日韩v欧美v天堂| 丰满的少妇人妻无码区| 亚卅精品无码久久毛片乌克兰 | 国产毛片片精品天天看视频| julia中文字幕久久亚洲| 亚洲一区二区三区国产精华液| 亚洲人成网站18禁动漫无码| 国产欧美视频综合二区| 久久99热66这里只有精品一| 国产手机在线小视频免费观看| 久久人与动人物A级毛片| 久久青草免费91观看| 国产成人精品在线| 91区国产福利在线观看午夜| 国产精品三区四区| 亚洲欧美国产视频| 亚洲第一成年网| 成人av专区精品无码国产| 草草影院国产第一页| 国产va在线观看免费| 亚洲黄色激情网站| 野花国产精品入口| 亚洲精品手机在线| 国产在线小视频| 福利在线不卡| 国产高清在线丝袜精品一区| 欧美日韩在线观看一区二区三区| 国产精品网址你懂的| 久久精品66| 专干老肥熟女视频网站| 999精品色在线观看| 日韩精品成人在线| 青青青视频免费一区二区| 日韩国产亚洲一区二区在线观看| 中文字幕亚洲乱码熟女1区2区| 亚洲无线观看| 国产精品成人一区二区不卡| 无码精品国产dvd在线观看9久| 色婷婷亚洲综合五月| 国产粉嫩粉嫩的18在线播放91| 国产九九精品视频| 91成人试看福利体验区| 中文字幕在线观| 国产欧美日韩综合一区在线播放| 久久综合五月婷婷| 亚洲综合精品香蕉久久网| 亚洲 欧美 偷自乱 图片 | 国产一级在线播放| 亚洲一区无码在线| 亚洲AV一二三区无码AV蜜桃| 亚洲一区二区约美女探花| 国产一区二区三区在线观看视频 | 中日无码在线观看| 五月六月伊人狠狠丁香网| 久久精品无码一区二区国产区| 风韵丰满熟妇啪啪区老熟熟女| 久久中文无码精品| 91精品亚洲| 国产一级精品毛片基地| 国产精品va免费视频| 国产成人精品一区二区免费看京| 在线视频亚洲欧美| 国产激爽大片高清在线观看| 国产精品成人一区二区不卡| 国产永久无码观看在线| 日韩在线网址| 久久亚洲高清国产| 精品三级网站| 欧美一区二区自偷自拍视频|