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

基于小波包和拉普拉斯特征值映射的柱塞泵健康評估方法

2017-11-30 05:49:33王浩任黃亦翔劉成良王雙園張大慶
振動與沖擊 2017年22期
關鍵詞:特征信號方法

王浩任, 黃亦翔, 趙 帥, 劉成良, 王雙園, 張大慶

(1.上海交通大學 機械系統與振動國家重點實驗室,上海 200240; 2.山河智能裝備集團,長沙 410100)

基于小波包和拉普拉斯特征值映射的柱塞泵健康評估方法

王浩任1, 黃亦翔1, 趙 帥1, 劉成良1, 王雙園1, 張大慶2

(1.上海交通大學 機械系統與振動國家重點實驗室,上海 200240; 2.山河智能裝備集團,長沙 410100)

柱塞泵是液壓系統的關鍵部件之一,監測其健康狀態對液壓系統的可靠運行具有重要意義。提出一種基于小波包和流形學習的方法,用于分析柱塞泵出口振動信號,從而對其進行健康評估;該方法利用小波包對原始信號進行分解,從中提取用于描述柱塞泵健康狀態的有效特征群;把提取的高維特征群作為輸入,利用并比較多種流形學習方法進行特征降維,選取狀態識別準確率最高的拉普拉斯特征映射方法,建立起的特征向量到健康狀態之間的對應關系,實現液壓泵健康狀態監測的分類要求。實驗結果表明,采用小波包和拉普拉斯特征映射相結合的方法可以有效提高柱塞泵狀態評估的準確性。

小波包分析;流形學習;柱塞泵;拉普拉斯特征映射;健康狀態評估

液壓泵廣泛應用于各種工業場合,但液壓泵健康評估機理復雜,缺乏理論研究模型,且檢測信號所包含的噪聲較多,所以難以對健康狀態識別。當液壓泵處于故障初期,一般表現為振動、沖擊、噪聲增加,并制約生產效率提高;隨著故障不斷加劇,液壓泵往往會因為故障引起壓力下降,最終導致液壓泵不能正常工作,甚至造成嚴重的安全事故。因此,通過對液壓泵的健康狀態的評估,判斷液壓泵的運行狀態,從而對可能產生的故障進行預防,具有極其重要的作用。小波分解(Wavelet Decomposition, WD)[1-2]是一種時頻域分解的方法,適用于非平穩信號的分解。小波包分解(Wavelet Package Decomposition, WPD)[3]利用一對高通濾波器和低通濾波器,能對信號的高頻和低頻成分作同樣的精細的分解。運用小波包的分析方法,彌補了小波分解對高頻信號分辨率的不足。對于具有非平穩特性的實際柱塞泵振動信號,需要對高頻和低頻信號成分進行較為全面的分析。使用小波包分解,能夠更好地提取信號的特征,獲得信號的潛在信息。

流形學習是一種用于高維空間數據降維的算法,作為一種維數約簡的手段,能挖掘高維非線性數據的幾何特征,從而在復雜的非線性空間中獲取信息。自Science雜志首次發表與流形學習相關算法后,相繼出現了局部線性嵌入算法(Local Linear Embedding, LLE)[4]、等距特征映射(Isometric Mapping, ISOMAP)[5]、拉普拉斯特征值映射(Laplacian Eigenmaps, LE)[6]、局部保持投影映射算法(Locality Preserving Projections, LPP)[7-8]等。其中流形學習在一些常見機械設備的故障診斷中已有應用:Jiang等[9]提出了監督拉普拉斯映射算法,并成功應用到了齒輪箱的故障診斷中;Yu等[10]提出了LPP的改進算法,并運用在軸承的故障檢測和狀態評估中;王紅軍等[11]提出了局部切空間排列算法(Local Tangent Space Alignment,LTSA)在數控機床上的應用,并成功用于機床主軸的狀態評估。雖然流形學習方法在一些常見的機械設備已經有了一定的研究,但對柱塞泵等液壓元件尚未應用。

針對目前國內外研究現狀,本文提出了一種針對柱塞泵健康評估研究方法。采用小波包分解對信號預處理,并提取時域特征作為高維特征向量。通過拉普拉斯特征映射等流形學習方法對高維特征向量進行降維處理,并使用K最近鄰方法(K-Nearest Neighbors, KNN)驗證分類效果。實驗結果表明,所提出的健康評估的方法能夠提高辨識精度,對柱塞泵的健康狀態識別具有重要的作用。

1 小波包系數分解

本文利用小波包對非平穩隨機的柱塞泵振動信號進行分解,能夠獲得信號在某一時間點的頻率特性。同時,小波包能夠更加精細地對原始信號進行分解和重構,細化不同健康狀態對應的信號頻率范圍,有利于準確評估健康狀態,提高模型對柱塞泵可靠性。小波包通過分解應用一對相關聯的低通濾波器和高通濾波器,將信號序列分解為某一尺度下的低頻和高頻兩部分。在改變尺度對已分解的低頻部分和高頻部分再次進一步分解,獲取更為細化的頻率成分。這種分解可以進行多次,以達到所需要的頻率分辨率。圖1為小波包三層分解示意圖。

圖1 小波包三層分解示意圖Fig.1 Wavelet package decomposition of three layers schemes

(1)

(2)

(3)

2 拉普拉斯特征值映射方法

拉普拉斯特征值映射(LE)[12]是一種非常有效的流形學習的算法。拉普拉斯特征值映射的基本思想是平均意義上保持數據點的鄰近信息,即通過特征映射,原本在高維空間上距離接近的點在低維空間中的映射后也應該距離接近。以兩個數據點之間的加權距離作為罰函數,利用圖拉普拉斯算子進行求解。拉普拉斯特征值映射具有良好的算法保持性,具有對噪聲不敏感,使用局部距離不易短路等優點。

LE算法可表述為以下三個步驟:

步驟1 構建近鄰圖。如果Xi和Xj是近鄰點,可以在結點i和j之間置一條邊,目前有ε-和K近鄰兩種方法。本文選用k=5的K近鄰方法構建鄰接圖。

步驟2 確定邊的權值Wij。邊的權值的確定一般有兩種方法:① 熱核方法(Heat Kernel),若第i個節點和第j個結點之間的是連接的,則定義邊的權值:Wij=exp(-‖xi-xj‖2/σ2),否則Wij=0;② 直接法,如果第i個和第j個節點之間有邊連接,則定義邊的權值為Wij=1,否則Wij=0。

步驟3 特征映射。假設前面所建立的鄰接圖是連通的,尋找低維的嵌入的問題實際上是對廣義特征向量的求解

Ly=λDy

(4)

(5)

經過特征值分解得到的第2項~第d+1項的特征值對應的特征向量也就是算法的d維輸出坐標。LE算法把問題轉化為了特征值的求解,不需要迭代值的計算,所以只需要很少的計算量。

3 故障特征的提取方法

3.1 高維原始特征向量的形成

柱塞泵運行時的健康狀態的識別的關鍵是通過對柱塞泵運行狀態特征的提取。由于實驗的設備復雜,所測得的振動信號中包含大量噪聲數據和冗余信息,所以難以直接通過振動信號進行評估。借助于特征提取的方法,把原始信號變換到高維的特征空間中,進而通過對特征的分析掌握設備運行的健康狀態。通過原始信號提取的時域,頻域和時頻混合域的參數已經被廣泛使用,但是對于不同的設備對應采集信號的性質各不相同,甄選出具有較高分辨率和規律性強的特征參數是狀態評估的關鍵。

本文采用常見的小波系數能量為計算特征向量的方法[13],同時根據柱塞泵振動信號的特點,選取最大絕對均值,最小絕對均值,絕對均值,方差,峰度和偏度六個參數作為信號的時頻信息,采用db6小波進行8層分解,得到全頻帶均勻劃分的256個特征作為子頻帶濾波信號,將各子頻帶的能量作為頻域的統計特征。表1所列的時域和頻域的262個特征構成數據表示為X∈RN×m,其中N為樣本數目,m為原始特征數,m=262。

表1 時域與頻域特征參數表

3.2 特征向量的歸一化處理

為了避免某些節點系數能量的絕對值過大導致在分析中造成偏差,一般需要對所得小波系數能量均進行規范化處理,使得各維度取值在一定范圍內避免較大偏差,且方便連續算法處理。采用規范化方法z-score算法,其公式為

(6)

式中:μtrain為從某一訓練集計算的均值;σtrain則是從同一訓練集中計算的標準差。

3.3 基于拉普拉斯特征映射對特征的提取

柱塞泵系統的結構復雜,要獲得精確的數學模型并不容易,因此模式識別是對柱塞泵的健康狀態進行評估的有效方法。

在對故障樣本進行特征提取后,再利用拉普拉斯特征向量良好的映射能力,對故障數據樣本的特征進行相關,將故障樣本特征映射到特征空間蘊涵的幾何關系作為分類特征,進而識別故障的類別。本文通過k=5的K最近鄰法(KNN)構建鄰接圖,使用熱核法確定權值,進而使用局部線性投影和低維特征提取。

如圖2所示,基于拉普拉斯特征向量柱塞泵診斷方法主要步驟為:

步驟1 從柱塞泵測試系統采集到的原始數據作為振動信號的樣本空間;

步驟2 對信號進行預處理,分別提取振動信號時域、時頻域的特征,共計262個特征參數,組成特征空間;

步驟3 選取特征空間的高維矩陣,對特征向量矩陣進行降維處理,選取的主要方法有LLE、LE、LPP、PCA(Principal Component Analysis)[14]、KPCA(Kernel Principal Component Analysis)[15-16]、Isomap[17-19]等;

步驟4 通過降維后的特征向量使用K最近鄰算法進行分類,計算步驟4不同的降維方法獲得的分類精度;

步驟5 評價不同降維方法的分類效果。

圖2 健康評估方案流程Fig.2 Health assessment program flow

本文方法的最大優點在于通過對特征集的時域、頻域和時頻域特征的選擇,得到了一組能夠通過表征系統故障特性的特征向量,通過降維之后的低維矩陣進行健康狀態分類,降低了分類的難度。

4 健康狀態分類實驗驗證

本文中的柱塞泵測試系統參考國標“液壓傳動電控液壓泵性能試驗方法:GB/T 23253—2009 ”中的標準進行搭建。實驗的對象為川崎K3V系列,是一種具有9個柱塞的斜盤式軸向柱塞泵。柱塞泵的壓力脈動頻率ω和柱塞泵的轉速n成正比,ω=18πn/60 rad/s。實驗選取1#~5#柱塞泵,在柱塞泵轉速為2 200 r/min工況下運行800 s,泵的出口處的脈動頻率為2 073.5 Hz。一般情況下,在柱塞泵保持運行狀態時,很少會出現兩種故障同時發生的情況,每次運轉柱塞泵的故障所造成的壓力沖擊只有一次,于是柱塞泵的沖擊頻率為ω1=2πn/60 rad/s。當柱塞泵的轉速為2 200 r/min時,故障特征頻率的集中在230.4 Hz。柱塞泵的振動信號通過JX61G系列振動傳感器測得。實驗使用的振動傳感器有3個,分別安裝在柱塞后泵油泵出口處的左上、右上和正下方,任意兩個傳感器之間的夾角為120°。傳感器的采樣頻率為50 kHz,滿足奈奎斯特采樣定理的要求。實驗臺搭建如圖3所示。

圖3 實驗臺系統Fig.3 Test-bed system

4.1 實驗Ⅰ柱塞泵磨損情況驗證

實驗用的柱塞泵在前400 h是柱塞泵的磨合期,柱塞泵振動信號的性質變化較大。為了細化柱塞泵在前400 h的工作情況,我們分別選取1#,2#和4#作為柱塞泵的磨合期測試實驗泵,選取5#作為正常工作(健康狀態)的測試柱塞泵,選取3#作為即將報廢的實驗泵。表2對應柱塞泵的實際運行時間和簡記代號。

實驗中選取了工作時間最長的3#柱塞泵進行拆解,對應部件的磨損情況如圖4所示。經過實際測量,可測得斜盤支撐座的最大磨損量為1.24 mm,最小磨損量為0.22 mm,而柱塞磨損量為0.06 mm。從圖4(c)可以看出,斜盤支撐座的高壓區的銅鍍層已經磨盡,而低壓區的銅鍍層磨損量相對較少,表明即使是同一部件,由于工況不同造成不同位置的磨損量不同。

表2 柱塞泵代號

4.2 實驗Ⅱ 不同柱塞泵健康狀態評估實驗

我們分別從P1、P2、P3、P4和P5柱塞泵選取了共500個數據樣本點,其中每種狀態樣本各選100個數據樣本點,特征空間維數為262,特征空間可以表示為X500×262。分別以LE、LPP、LLE、Isomap、PCA和KPCA算法對數據進行維數約簡處理,轉化到低維空間進行特征提取以進行模式識別分析。設參數的鄰近因子k=5,使用最大似然估計計算嵌入維數d′=23。

圖4 3#柱塞泵不同部件磨損情況Fig.4 3# wearing condition of different piston parts

如圖5(a)所示,是拉普拉斯故障信號明顯被分成了5類,類內距小,聚合明顯,容易區分,其中P1和P2類間距較近,說明柱塞泵在前100 h和100~250 h的特征區別不明顯,可以作為同一類型的健康狀態處理。P3、P4和P5類間距大,類內距小,健康狀態區別明顯,P3和其他幾類間距更大,區別明顯。如圖5(b)和圖5(e)所示分別采用LPP和PCA方法進行分類,類間距大,樣本分散,聚類性差。如圖5(f)所示,采用KPCA方法進行分類,雖然使用了核函數方法使類內距變小,但是類間距小,分類效果差。如圖5(c)和圖5(d)所示,分別采用LLE方法和Isomap方法對數據進行降維,和LE方法同屬非線性降維方法,特點是降維后類內距小,聚類效果優于線性降維方法,但是P4和P5類出現了不同程度的混疊,聚類效果也不如LE方法。

圖5 柱塞泵1#~5#特征向量降維結果Fig.5 Feature vectors reduction results of piston 1#~5#

為了比較不同嵌入維數對分類精度的影響,我們計算了使用LE方法降維后的數據在不同嵌入維數的下的降維結果,然后計算出使用KNN作為分類器的辨識精度,實驗結果如圖6所示。實驗結果表明,當嵌入維數23lt;d′≤262,此時特征維數較高,原始特征矩陣包含大量噪聲等冗余信息,最小特征值對應的特征向量往往和噪聲有關,降維能起到降噪的效果;當嵌入維數0lt;d′lt;23,此時特征維數較低,一定程度上損失了大量有用信息,此時的辨識率隨著維數的不斷下降而減少;使用最大似然估計計算嵌入維數d′=23獲得最佳分類精度98.649%。

圖6 不同嵌入維數下的分類精度Fig.6 Classification accuracy of different embedding dimension

為了評價不同降維方法對應的分類效果,我們選用K最近鄰分類法驗證模型分類精度。設置近鄰因子k=5,把實驗樣本按照7∶3分為訓練樣本和測試樣本兩類,訓練樣本用于訓練模型,利用測試樣本數據來驗證模型的準確性。不同降維方法對應模型的分類精度按10次平均計算,并由此獲得最佳的效果。如表3所示,使用LE、LLE、KPCA等不同降維方法獲得的降維結果經過K最近鄰方法驗證,比較不同驗證模型的準確性,LE方法具有最高分類精度97.973%,表明使用LE方法對于柱塞泵的運行狀態具有較高的區分性,評估最為準確。

表3 柱塞泵1#~5#降維后的數據辨識結果

5 結 論

本文針對柱塞泵的健康評估問題,提出了一種基于小波包和拉普拉斯特征值映射的分析方法。該方法利用小波包對高頻信號的分辨率高的特點,選取時域特征和小波包系數的能量作為特征,使用流形學習方法降維處理,并做出狀態評估。分析實驗結果表明:

(1) 對于柱塞泵柱塞泵健康狀態評估方法,通過LE方法降維后的數據都具有更好的聚類性,狀態識別的準確率比LPP、LLE、Isomap、PCA和KPCA等方法更高。通過對柱塞泵狀態的準確評估,能對有可能產生的故障進行預防,具有重要的意義。

(2) 利用小波包對信號進行分解,能同時對高頻和低頻信號進行精細劃分,再加上能表征信號特性的時域特征,實驗表明,所選取的特征能較好地反映柱塞泵的狀態信息,對柱塞泵的狀態的評估具有重要作用。

(3) 在今后的工作中,將會采用不同型號的柱塞泵,分析驗證該模型的準確性,同時使用在工程機械的實際應用場合,結合周圍元件對柱塞泵的影響,驗證實際工況下不同條件對液壓泵健康狀態的影響。

[ 1 ] 何正嘉,訾艷陽,孟慶豐.機械設備非平穩信號的故障診斷原理及應用[M].北京:北京高等教育出版社,2001.

[ 2 ] DAUBECHIES I, BATES B J. Ten lectures on wavelets[J]. The Journal of the Acoustical Society of America, 1993, 93(3): 1671.

[ 3 ] 高英杰,孔祥東. 基于小波包分析的液壓泵狀態監測方法[J]. 機械工程學報, 2009, 45(8): 80-88.

GAO Yingjie, KONG Xiangdong. Wavelet packets analysis based method for hydraulic pump condition monitoring [J]. Journal of Mechanical Engineering, 2009, 45(8): 80-88.

[ 4 ] ROWEIS S T, SAUL L K. Nonlinear dimensionality reduction by locally linear embedding[J]. Science, 2000, 290: 2323-2326.

[ 5 ] TENENBAUM J B, SILVA V D, LANGFORD J C. A global geometric framework for nonlinear dimensionality reduction[J]. Science, 2000, 290: 2319-2323.

[ 6 ] 蔣全勝,賈民平,胡建中,等. 基于拉普拉斯特征映射的故障模式識別方法[J].系統仿真學報,2008,20(20): 5710-5713.

JIANG Quansheng,JIA Minping,HU Jianzhong,et al. Method of fault pattern recognition based on Laplacian eigenmaps[J]. Journal of System Simulation,2008,20(20):5710-5713.

[ 7 ] HE X,NIYOGI P. Locality preserving projections[M]. Cambridge: MIT Press,2003:153-160.

[ 8 ] 徐金梧,呂勇,王海峰. 局部投影算法及其在非線性時間序列分析中的應用[J]. 機械工程學報,2003,39(9):146-150.

XU Jinwu, Lü Yong, WANG Haifeng. Local projective method andit’s application onnonlinear time series[J]. Journal of Mechanical Engineering,2003,39(9): 146-150.

[ 9 ] JIANG Quansheng,JIA Minping,HU Jianzhong,et al. Machinery fault diagnosis using supervised manifold learning[J]. Mechanical Systems and Signal Processing,2009,23(7): 2301-2311.

[10] YU Jianbo. Local and non-local preserving projection for bearing defect classification and performance assessment[J]. IEEE Trans. Ind. Electron,2012,59(5): 2363-2376.

[11] 王紅軍, 徐小力,萬鵬. 基于軸心軌跡流形拓撲空間的轉子系統故障診斷[J]. 機械工程學報, 2014, 50(5): 95-101.

WANG Hongjun, XU Xiaoli, WAN Peng. Rotor system fault diagnosis based on orbit manifold topological[J]. Journal of Mechanical Engineering, 2014, 50(5): 95-101.

[12] 黃宏臣,韓振南,張倩倩,等. 基于拉普拉斯特征映射的滾動軸承故障識別[J]. 振動與沖擊, 2015, 34(5): 128-134.

HUANG Hongchen,HAN Zhennan,ZHANG Qianqian, et al. A method for fault diagnosis of rolling bearings based on laplacian eigenmap[J]. Journal of Vibration and Shock, 2015, 34(5): 128-134.

[13] 丁曉喜,何清波. 基于WPD和LPP的設備故障診斷方法研究[J].振動與沖擊, 2014, 33(3): 89-93.

DING Xiaoxi, HE Qingbo. Machine fault diagnosis based on WPD and LPP[J]. Journal of Vibration and Shock, 2014, 33(3): 89-93.

[14] 高東,吳重光,張貝克,等. 基于PCA和SDG的傳感器故障診斷方法研究及應用[J].系統仿真學報, 2011, 23(3): 567-573.

GAO Dong, WU Chongguang, ZHANG Beike, et al. Method combining PCA and SDG for fault diagnosisof sensors and its application[J]. Journal of System Simulation, 2011, 23(3): 567-573.

[15] XU Yong, ZHANG D, SONG F, et al. A method for speeding up feature extraction based on KPCA[J]. Neurocomputing, 2007, 70(4/5/6): 1056-1061.

[16] 楊先勇,周曉軍,張文斌,等. 基于局域波法和 KPCA-LSSVM的滾動軸承故障診斷[J]. 浙江大學學報(工學版), 2010, 44(8):1519-1524.

YANG Xianyong, ZHOU Xiaojun, ZHANG Wenbin, et al. Rolling bearing fault diagnosis based on local wave methodand KPCA-LSSVM[J]. Journal of Zhejiang University (Engineering Science), 2010, 44(8): 1519-1524.

[17] BALASUBRAMANIAN M, SCHWARTZ E L. The isomap algorithm and topological stability[J]. Science, 2002, 295: 7.

[18] 孟德宇,徐晨,徐宗本. 基于Isomap的流形結構重建方法[J]. 計算機學報, 2010, 33(3): 545-555.

MENG Deyu, XU Chen, XU Zongben. A new manifold reconstruction method based on Isomap[J]. Chinese Journal of Computers, 2010, 33(3): 545-555.

[19] 劉愛萍, 王洪元, 程起才, 等. 基于KISOMAP-LDA-KNN算法TE過程故障診斷研究[J]. 計算機與數字工程, 2010, 38(11): 34-37.

LIU Aiping, WANG Hongyuan, CHENG Qicai, et al. An algorithm of KISOMAP-LDA-KNN for the research of the fault diagnosis of TE process[J]. Computer and Digital Engineering, 2010, 38(11): 34-37.

HealthassessmentforapistonpumpbasedonWPDandLE

WANG Haoren1, HUANG Yixiang1, ZHAO Shuai1, LIU Chengliang1, WANG Shuangyuan1, ZHANG Daqing2

(1.State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai 200240, China;2. Sunward Intelligent Equipment Group, Changsha 410100, China)

A piston pump is one of the key components in a hydraulic system, monitoring its health condition is of significant importance for reliable performance of the hydraulic system. Therefore, a health state evaluation method was proposed based on the wavelet packet decomposition (WPD) and the manifold learning by analyzing the vibration signals at the piston outlet. The wavelet packet method was used to decompose original signals and effectively extract health state features group from them. The high dimensional feature group was set as an input and multiple manifold learning methods were conducted and compared for dimensional reduction. The Laplacian eigenmaps (LE) method of the highest accuracy was used to establish a relationship between feature vectors and health states, which achieved the aim of health state classification. It is shown that the combination of the wavelet packet decomposition and the manifold learning method improve the accuracy of piston pump health state evaluation.

wavelet packet analysis; manifold learning; piston pump; Laplacian eigenmaps; health state evaluation

國家科技支撐計劃(2014BAA04B01);國家自然科學基金(51305258);上海市科委項目(1411104600)

2015-12-10 修改稿收到日期: 2016-09-14

王浩任 男,博士生,1991年生

黃亦翔 男,博士,講師,1980年生

TP306

A

10.13465/j.cnki.jvs.2017.22.008

猜你喜歡
特征信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
抓住特征巧觀察
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 中国国产A一级毛片| 波多野结衣一区二区三区四区视频 | 97视频在线精品国自产拍| AV色爱天堂网| 国产熟女一级毛片| 亚洲综合婷婷激情| 狠狠v日韩v欧美v| 片在线无码观看| 国产精品主播| 国产午夜福利亚洲第一| 久久久久国色AV免费观看性色| 国产一级精品毛片基地| 奇米影视狠狠精品7777| 狠狠v日韩v欧美v| 成人91在线| 成人一级免费视频| 国产亚洲成AⅤ人片在线观看| 91精品小视频| 国产第一页第二页| 91香蕉国产亚洲一二三区| 波多野衣结在线精品二区| 青青青亚洲精品国产| 亚洲人成电影在线播放| 久久精品国产精品青草app| 伊人久久大香线蕉影院| 亚洲午夜福利在线| 日本不卡在线视频| 亚洲高清中文字幕在线看不卡| 日韩不卡免费视频| 欧美亚洲香蕉| 无码日韩精品91超碰| 69国产精品视频免费| 亚洲欧洲免费视频| 制服丝袜在线视频香蕉| 欧美a在线看| …亚洲 欧洲 另类 春色| 国产AV毛片| 亚洲色无码专线精品观看| a级免费视频| 99热国产这里只有精品9九| 伊人久综合| 国产成人h在线观看网站站| 精品欧美日韩国产日漫一区不卡| 91色爱欧美精品www| 国产精品成人第一区| 欧美日本在线观看| 91口爆吞精国产对白第三集| 免费观看欧美性一级| 91在线精品免费免费播放| 国产美女视频黄a视频全免费网站| 日韩欧美高清视频| 一本大道香蕉久中文在线播放| 自拍偷拍欧美日韩| 超碰aⅴ人人做人人爽欧美| 一边摸一边做爽的视频17国产| 国产午夜精品鲁丝片| 欧美在线三级| 欧美日韩综合网| 国产美女无遮挡免费视频| 亚卅精品无码久久毛片乌克兰 | 色婷婷成人网| 欧美亚洲国产视频| 国产一级做美女做受视频| 国产精女同一区二区三区久| 日韩无码视频播放| 黄色免费在线网址| 国产在线视频导航| 性网站在线观看| 国产呦视频免费视频在线观看| 波多野结衣一二三| 青青草欧美| 欧美一级在线| 国产网站免费| 欧美亚洲国产一区| 亚洲va欧美ⅴa国产va影院| 国产一级在线播放| 欧美在线伊人| 91网址在线播放| 99伊人精品| AV色爱天堂网| 亚洲欧美成人网| h网址在线观看|