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

基于EMD-WVD與LNMF的內燃機故障診斷

2017-01-10 08:14:42牟偉杰石林鎖蔡艷平金廣智
振動與沖擊 2016年23期
關鍵詞:故障診斷振動信號

牟偉杰, 石林鎖, 蔡艷平, 劉 浩, 金廣智

(第二炮兵工程大學 五系,陜西 西安 710025)

基于EMD-WVD與LNMF的內燃機故障診斷

牟偉杰, 石林鎖, 蔡艷平, 劉 浩, 金廣智

(第二炮兵工程大學 五系,陜西 西安 710025)

內燃機的振動信號是復雜非平穩信號,準確提取內燃機振動信號中的特征信息進行模式識別,是對振動信號進行故障診斷的關鍵。基于經驗模態分解的維格納時頻分析方法,不但保留了維格納分布的所有優良特,而且還避免了交叉項的干擾,能夠有效地提取內燃機振動信號的特征信息;在此基礎之上,針對傳統非負矩陣分解非正交的基矩陣導致數據冗余性較大、影響后續故障分類準確率提高的問題,提出采用局部非負矩陣分解的方法,直接對EMD-WVD時頻圖像的矩陣進行分解,計算用于內燃機故障診斷的特征參數,并利用特征參數進行故障分類。對內燃機4種不同工況的振動信號進行實驗,證明基于EMD-WVD與局部非負矩陣分解的方法對內燃機氣門間隙的故障診斷的有效性。

內燃機;故障診斷;時頻分布;特征提取;局部非負矩陣分解

在內燃機故障診斷中,振動信號往往含有豐富的故障信息,因此內燃機的振動診斷是其故障診斷的重要研究領域。但傳統的振動診斷方法存在許多不足,時域分析方法所能提供的信息量非常有限,因此基本用于機械設備的簡易診斷;基于傅里葉變換的頻譜分析建立在信號平穩性的假設條件之上,且分析的結果只有頻域信息,喪失了時域特征。而大多數機械故障振動信號是非平穩和非線性的,對這些非平穩信號,由于傅里葉變換的本質缺陷,使得提取的故障特征信息有缺陷,從而降低了故障診斷的準確率。內燃機的振動信號屬于非平穩信號,因此時頻分析方法成為內燃機故障診斷的主要方法,如短時傅里葉變換[1-2]、維格納時頻分析[3-4]、小波變換[5-6]以及希爾伯特-黃變換等方法,這些時頻分析方法為處理內燃機振動信號中突變、不連續、非平穩信號提供了必要的手段。在用這些方法對內燃機振動信號進行分析時, 短時傅里葉變換時窗大小及形狀都固定不變,不能同時滿足時頻分辨率的要求,小波分解會明顯出現多余信號,而希爾伯特-黃變換由于經驗模態分解時的模態混疊現象的存在,使得到的時頻圖像雜亂不堪。用維格納時頻分析方法對內燃機振動信號進行分析時,對信號的瞬時頻率和局域化都有很好的描述,但其交叉項成為應用時的瓶頸,交叉項的存在很難將有多個頻率成分的信號表示清楚,本文采用基于經驗模態分解的維格納時頻分析方法(Empirical Mode Decomposition-Wigner-Ville Distribution,EMD-WVD)對內燃機振動信號進行分析,不但避免了交叉項的干擾,而且還保留了WVD分布的所有優良特性,能夠有效地提取內燃機缸蓋振動特征信息。

近年來,通過對振動信號的時頻圖像提取特征參數,然后進行故障分類的方法,逐漸應用于對內燃機故障的智能診斷。一些學者對這方面也進行了大量研究,但是這些方法的研究重點基本在于對時頻圖像的特征提取上[7-11],而且這些方法大多特征指標選擇方法模糊,計算復雜。

非負矩陣分解(Nonnegative Matrix Factorization,,NMF)算法是由LEE等[12]提出的特征提取方法。該方法可以針對不同的圖像自適應地計算圖像的特征參數,避免了不同圖像特征指標的選擇,而且計算簡單,具有明確的物理意義, 已廣泛用于圖像工程、計算機視覺和模式識別等研究領域[13-15]。有學者將非負矩陣分解技術與時頻分析相結合應用于故障診斷[16-21],取得了較好的效果。但是傳統NMF中存在非正交的基矩陣,導致數據冗余性較大,影響識別率的提高。

本文提出采用基于EMD-WVD與局部非負矩陣分解(Local Nonngeative Matrix Factorization,LNMF)的內燃機故障診斷方法,對內燃機4種工況下的振動信號進行EMD-WVD時頻分析,得到EMD-WVD的時頻分布圖像, 再分別采用非負矩陣分解和局部非負分解方法對EMD-WVD時頻分布圖進行特征提取并進行分類,對比分析兩種方法在不同的特征維數下的分類性能。

1 基于經驗模態分解的維格納時頻分析[22]

維格納分布是一種二次時頻分布,根據二次疊加原理,如果信號存在兩個以上的分量,就會存在交叉干擾項,這也成了維格納分布在實際應用中受到限制的主要原因。為消除多分量信號維格枘分析時的交叉干擾項,可以將一個多分量信號分解成幾個單分量信號,然后把每個單分量信號進行維格納分析,并將時頻分析后的結果相加組合成原信號的維格納時頻分布。而基于經驗模態分解的維格納時頻分析方法,就是利用這一思想,利用EMD方法可以將多分量信號分解成多個單分量信號的性質,既有效地消除了WVD的交叉干擾項,又保留了原來WVD方法的優良性質[22]。

信號x(t)的EMD-WVD時頻分布定義為

(1)

其算法的基本流程是:先對多分量信號進行EMD分解,生成內稟模態分量C1,C2,…,Cn,然后對各內稟分量分別進行WVD計算,最后將各WVD結果進行線性疊加,所得結果即為信號的EMD-WVD時頻分析結果,EMD-WVD算法的流程見圖1。

圖1 EMD-WVD方法的算法流程

2 局部非負矩陣分解

2.1 非負矩陣分解

非負矩陣分解的本質是一種矩陣分解和投影技術,其基本原理如下:

對于非負矩陣V,存在W≥0、H≥0,滿足:

Vn×m=Wn×rHr×m

(2)

式中:n為數據樣本的維數;m為集合中數據樣本的個數;r為特征維數,一般情況下,r的選擇滿足nm>(n+m)r,所以得到的W和H小于原始矩陣V,從而實現數據的降維。稱W為基矩陣,H為系數矩陣。

為了描述V≈WH的近似效果,一般通過分解前后V與WH兩個非負矩陣之間一些距離來定義目標函數。一種方法是利用矩陣V與WH間的歐氏距離,對應的目標函數為

(3)

另一種方法是利用矩陣V與WH間的廣義K-L散度,對應的目標函數

(4)

式(3)和式(4) 對應的優化問題為

minE(V‖WH),并且W≥0,H≥0;

minD(V‖WH),并且W≥0,H≥0。

對目標函數式(3),相應的迭代規則為

(5)

(6)

對目標函數式(4),相應的迭代規則為

Wik←Wik∑jHkjVij/(WH)ij

(7)

(8)

Hkj←Hkj∑iWikVij/(WH)ij

(9)

2.2 局部非負矩陣分解

局部非負矩陣分解方法是LI等[23]提出的,其根本目的是為了使分解后的基圖像能夠得到更為局部的特征,并將該方法成功用在人臉識別中。

LNMF具體算法是通過對K-L散度目標函數(式(4))施加基的列正交性約束,以減少基向量之間的冗余。

LNMF的目標函數如下所示:

(10)

其迭代式為:

(11)

(12)

(13)

3 內燃機故障診斷實例

基于EMD-WVD與局部非負矩陣分解的故障診斷方法的步驟見圖2,首先對采集到的信號進行EMD-WVD時頻分析得到時頻圖像,然后對時頻圖像進行局部非負矩陣分解得到圖像的特征參數,最后用分類器對特征參數進行分類完成對內燃機的故障診斷。

圖2 基于EMD-WVD與LNMF的故障診斷方法的步驟

3.1 內燃機試驗工況

選取道依茨BF4L1011F型柴油機為研究對象,采樣頻率為25 kHz,轉速為3 000 r/min,測試過程中,內燃機空載運行,共設置了4種氣門間隙狀況,分別對應氣門間隙正常與不當,具體情況如表1 。其中,正常進、排氣門間隙0.3 mm、0.4 mm,氣門間隙過小為0.06 mm,氣門間隙過大為0.6 mm,嚴重漏氣為開口4 mm(長)×1 mm(寬)。

試驗中采集內燃機氣門4種故障狀態下各70種振動信號樣本,共計280個,圖6給出了4種工況下振動信號的波形。每種工況的前35種作為訓練樣本,其余作為測試樣本。

表1 四種實驗工況設置

3.2 內燃機缸蓋振動信號的EMD-WVD時頻分析

圖3~圖5分別給出了氣門間隙正常時缸蓋振動信號的時域圖、功率譜圖、WVD時頻分布圖和EMD-WVD時頻分布圖。從圖3可以看出,在62~66 ms之間、以及74 ms處,振動信號有較大的振動幅度,并且振動信號的頻譜主要集中在6~10 kHz。比較圖4和圖5可知,用WVD方法對缸蓋振動信號進行分析時存在嚴重的交叉干擾項,而EMD-WVD方法可以清晰地分辨出62~66 ms時振動信號的頻率為6~10 kHz,在74 ms時振動信號的頻率為6. 5~ 9. 5 kHz,有效地抑制了WVD方法的交叉項。

圖3 振動信號的時域和功率譜圖

圖4 振動信號的WVD時頻圖

圖5 振動信號的EMD-WVD時頻圖

圖6為采用EMD-WVD時頻分析方法對四種典型工況下整循環缸蓋振動信號進行時頻分析而得到的時頻分布圖,為使圖像更清晰, 便于觀察, 圖6 給出的時頻圖像是用等高線表示的, 而在實際的計算及分類識別研究中, 使用的都是灰度圖像。在圖像中,橫坐標表示時間,縱坐標表示頻率,圖像的灰度代表幅值,灰度越深幅值越大。從圖6中可知,缸蓋振動信號的頻率主要出現在5~12 kHz之間的高頻區,在氣門間隙正常時,信號的頻率主要出現在6.5~9.5 kHz之間,其它幾種故障狀態下,頻率主要出現在6~12 kHz之間,從頻率上可以分別出故障狀態。另外,比較圖6(b)、圖6(c)和圖6(d),雖然3種故障狀態的頻率主要出現在6~12 kHz的高頻區,但在時頻圖上出現和消失時間并不相同,幅值大小不同,而且它們頻率組成更不相同, 例如隨著氣門間隙的增大,振動信號在62~74 ms時的頻率有向高頻移動和集中的趨勢。對應于灰度圖上,反映不同工況灰度圖像的形狀(頻率組成)、出現的位置(時間)和灰度值深淺(幅值)具有較大差異,可以對四種工況進行分類。

圖6 四種典型工況缸蓋振動信號的EMD-WVD時頻分析Fig.6 EMD-WVD time-frequency image for four states

3.3 對EMD-WVD時頻圖的非負矩陣分解

由于EMD-WVD時頻分布圖的矩陣維數高,計算量大,不利于進行特征參數的提取,因此在保留時頻圖像主要特征的前提下,將EMD-WVD時頻圖轉換為175×55的灰度圖像,降低了矩陣維數,大大減少的計算量。

在進行非負矩陣分解時要先通過矩陣重排,將時頻分布灰度圖像的矩陣向量化,將每一個灰度圖像矩陣由175×55維變形為9 625×1維的列向量,并對其進行歸一化處理。

圖7給出了特征維數為40時,EMD-WVD時頻圖像測試集對應的特征系數,為直觀表示,把每種工況的特征系數用編碼矩陣的形式表示,每一個樣本的特征系數對應一個8×5的編碼矩陣,編碼矩陣圖中每個像素的灰度值與該樣本系數一一對應。每種工況下選取5個樣本,圖中每一行代表一種內燃機工況,從上到下依次為氣門間隙正常、過小、過大和漏氣。從圖中可以看出,每種工況編碼矩陣圖中各像素灰度相似,說明每種工況的系數相近,便于后續分類計算。

圖7 LNMF提取的測試集特征系數

3.4 內燃機故障分類

由圖8可知,對于基于EMD-WVD時頻圖和局部非負矩陣分解的方法可以達到很好的分類效果,而且在相同特征維數的情況下,LNMF特征提取方法的分類精度始終優于NMF方法。其根本原因在于基圖像的特征維數相同時,LNMF方法在對時頻圖像矩陣進行分解時,對基矩陣各向量之間的正交性施加了約束,使分解后得到的圖像的基矩陣各向量之間盡量正交,減少基矩陣的冗余,使重構圖像的信息更精確,提高了分類精度。

但通過圖8也可知,當基向量的特征維數為40時,所提取的特征參數取得了最優的分類精度,隨著特征維數的增加,分類精度反而有下降的趨勢。因為,特征維數的值反映了用于表達原始圖像的基圖像的個數,對于用相同方法進行分非負矩陣分解的具體圖像,特征維數的值是確定的。當特征維數的取值與具體圖像一致時,得到的基圖像剛好可以重構所有原始圖像;當維數的取值小時,降維效果顯著,但圖像特征提取的信息不足,難以滿足后續分類的要求;當維數太大時,有可能帶入信息冗余或噪聲,影響識別率。目前,對于如何確定特征維數的值,以獲得最佳的特征參數的問題還沒有很好的計算方法。目前一般采用的方法是選取多個特征維數的值分別進行實驗,取得最佳識別效果的值即為特征維數的最終取值。

圖8 LNMF和NMF特征提取方法的識別率

為驗證結果的通用性,用另外兩種常見的分類器,樸素貝葉斯分類器(NBC)和支持向量機分類器(SVM)重復以上實驗30次取平均值。表2為特征維數為40時三種分類器分類精度的對比,由表可以看出對于不同的分類器,基于EMD-WVD時頻圖和局部非負矩陣分解的方法有很好的分類效果,并且LNMF方法的精度優于NMF方法,說明該方法具有很好的通用性。

4 結 論

(1)EMD-WVD時頻分析方法可以有效抑制WVD方法的交叉干擾項,明顯區分不同氣門間隙工況的時頻分布特點,能夠用于內燃機氣門間隙的故障診斷。

表2 內燃機故障診斷精度

(2)用局部非負矩陣分解方法對EMD-WVD時頻圖像直接進行特征參數提取,并進行故障診斷的方法有較高的分類精度,而且LNMF方法比NMF方法的分類精度更高,可以用于對EMD-WVD時頻圖像分類。

[ 1 ] KIM Y H. Fault detection in a ball bearing system using a moving window [J].Mechanical Systems and Signal Processing,1991,5(6):461-473.

[ 2 ] SAMIN B, RIZZONI G. Engine knock analysis and detection using time-frequency analysis[C]//SAE960618,1996.

[ 3 ] XIAO J, FLANDFIN P. Muhitaper time-frequency reassignment for nonstationary spectrum estimation and chirp enhancement[J].IEEE Transaction on Signal Processing,2007,55(6):2851-2860.

[ 4 ] 毛永芳,秦樹.重分配譜圖和多窗譜在機械故障診斷中的應用[J].振動與沖擊,2009,28(1):161-165. MAO Yongfang, QIN Shu. Re-allocation of spectrum and multi-spectral windows and its application in machinery fault diagnosis [J].Journal of Vibration and Shock,2009,28(1):161-165.

[ 5 ] BO L, QIN S R, LIU X F. Theory and application of Wavelet analysis instrument library[J]. Chinese Journal of Mechanical Engineering: English Edition,2007,19(3):464-467.

[ 6 ] CARY S, AKUJUOBI C M. An approach to vibration analysis using wavelets in an application of aircraft health monitoring[J].Mechanical Systems and Signal Processing,2007,21:1255-1272.

[ 7 ] 張金玉,張優云,謝友柏. 機械圖像的識別與重構[J].西安交通大學學報,2000,34(2):80-84. ZHANG Jinyu, ZHANG Youyun, XIE Youbai. Recognition and reconstruction of mechanical graphics[J]. Jouranl of Xi’an Jiaotong University,2000,34(2), 80-84.

[ 8 ] 夏勇,商斌梁,張振仁, 等. 基于小波包與圖像處理的內燃機故障診斷研究[J].內燃機學報,2001,19(1):62-68. XIA Yong, SHANG Binliang, ZHANG Zhenren, et al. Research on fault diagnosis for internal combustion engines using wavelet packet and image processing[J].Transaction of CSICE, 2001,19(1):62-68.

[ 9 ] 沈虹,趙紅東,梅檢民,等. 基于高階累積量圖像特征的柴油機故障診斷研究[J].振動與沖擊,2015,34(11):133-138. SHEN Hong, ZHAO Hongdong, MEI Jianmin, et al. Diesel engine fault diagnosis based on high-order cumulant image features[J]. Journal of Vibration and Shock, 2015,34(11):133-138.

[10] 王成棟,張優云,夏勇. 基S變換的柴油機氣閥機構故障診斷研究[J].內燃機學報,2003,21(4):271-275. WANG Chengdong, ZHANG Youyun, XIA Yong. Fault diagnosis for diesel valve train based on S Transform[J]. Transaction of CSICE, 2003,21(4):271-275.

[11] 王成棟,張優云,夏勇. 模糊函數圖像在柴油機氣閥故障診斷中的應用研究[J]. 內燃機學報,2004,22(4):162-168. WANG Chengdong, ZHANG Youyun, XIA Yong. Study on the use of ambiguity function images in fault diagnosis for diesel valve train[J]. Transaction of CSICE, 2004,22(4):162-168.

[12] LEE D D, SEUNG H S. Learning the parts of objects by no n-negative matrix factorization [J].Nature,1999,401(6755):788-791.

[13] LIU W, ZHENG N. Non-negative matrix factorization based methods for object recognition [J].Pattern Recognition Letters,2004,25(8):893-897.

[14] ZHANG T, FANG B, TANG Y Y, et al. Topology preserving non-negative matrix factorization for face recognition[J].IEEE Transactions on Image Processing,2008, 17(4):574-584.

[15] BUCAK S S, GUNSEL B. Incremental subspace learning via non-negative matrix factorization [ J ] . Pattern Recognition, 2009, 42(5): 788-797.

[16] QING H W, YUN Z Y, LEI C, et al. Fault diagnosis for diesel valve trains based on non-negative matrix factorization and neural network ensemble[J].Mechanical Systems and Signal Processing,2009,23(5):1683-1695.

[17] 蔡蕾, 朱永生. 基于稀疏性非負矩陣分解和支持向量機的時頻圖像識別[J].自動化學報, 2009, 35(10):1272-1277. CAI Lei, ZHU Yongsheng. Time-frequency spectra recognition based on sparse non-negative matrix factorization and support vector machine[J]. Acta Automatica Sinica, 2009, 35(10):1272-1277.

[18] 李兵,高敏,張旭光,等. 用形態梯度法與非負矩陣分解的齒輪故障診斷[J].振動、測試與診斷,2014,34(2):295-300. LI Bing, GAO Min, ZHANG Xuguang, et al. Feature extraction for engine fault diagnosis by utilizing adaptive multi-scale morphological gradient and non-negative matrix factorization[J].Journal of Vibration, Measurement and Diagnosis, 2014,34(2):295-300.

[19] 李兵,米雙山,劉鵬遠,等. 二維非負矩陣分解在齒輪故障診斷中的應用[J]. 振動、測試與診斷,2012,32(5):836-840. LI Bing, MI Shuangshan, LIU Pengyuan, et al. Application of two-dimensional non-negative factorization for gear fault diagnosis[J]. Journal of Vibration, Measurement and Diagnosis, 2014,34(2):295-300.

[20] 李兵,徐榕,賈春寧,等. 基于自適應形態提升小波與改進非負矩陣分解的發動機故障診斷方法[J].兵工學報,2013,34(3):353-360. LI Bing, XU Rong, JIA Chunning, et al. Engine fault diagnosis Utilizing adaptive morphological lifting wavelet and improved non-negative matrix factorization[J].Acta Armamentarii, 2013,34(3):353-360.

[21] 劉東升,李元,楊博文,等. 基于自適應多尺度形態梯度與非負矩陣分解的軸承故障診斷[J].振動與沖擊,2013,32(19):106-110. LIU Dongsheng, LI Yuan, YANG Bowen, et al. Bearing fault diagnosis based on adaptive multi-scale morphological gradient andnon-negative matrix factorization[J].Journal of Vibration and Shock, 2013,32(19):106-110.

[22] 蔡艷平,李艾華,王濤,等.基于EMD-Wigner-Ville的內燃機振動時頻分析[J].振動工程學報,2010,23(4):430-437. CAI Yanping, LI Aihua, WANG Tao, et al. Engine vibration time-frequency analysis based on EMD-Wigner-Ville[J].Journal of Vibration Engineering, 2010,23(4):430-437.

[23] LI S Z, HOU X W, ZHANG H J, et al. Learning spatially localized,parts-based representation[J]. Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition,2001(1):207-212.

IC engine fault diagnosis method based EMD-WVD and LNMF

MU Weijie, SHI Linsuo, CAI Yanping, LIU Hao, JIN Guangzhi

(Department NO.5, the Second Artillery Engineering University, Xi’an 710025, China)

IC engine vibration signals are complex and non-stationary, therefore, how to accurately extract the feature information of IC engine vibration signals for pattern recognition is very important for fault diagnosis of an IC engine. Here, Wigner-Ville time-frequency distribution analysis method based on empirical mode decomposition (EMD-WVD) was used for IC engine vibration fault diagnosis. The method of EMD-WVD could not only avoid interferences of cross-terms, but also retain all the excellent characteristics of Wigner-Ville distribution to effectively extract fault features of IC vibration signals. A new method of locals non-negative matrix factorization (LNMF) was directly used for feature extraction due to the traditional NMF non-orthogonal basis matrix with larger data redundancy, and being not conducive to subsequently improving the correct rate of faut classification. The application of LNMF in practical IC engine fault diagnosis showed that LNMF algorithm’s fault classification accuracy is better than that of the traditional NMF algorithm. The test results for IC engine vibration signals under 4 different operational conditions verified the effectiveness of EMD-WVD and LNMF methods in IC engine fault diagnosis.

IC engine; fault diagnosis; time-frequency distribution; feature extraction; local nonnegative matrix factorization (LNMF)

國家自然科學基金青年基金項目(51405498);陜西省自然科學基金項目(2013JQ8023)

2015-09-01 修改稿收到日期:2015-11-10

牟偉杰 男,博士生,1984年生

石林鎖 男,教授,博士生導師,1958年生

TK428

猜你喜歡
故障診斷振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
基于LabVIEW的力加載信號采集與PID控制
因果圖定性分析法及其在故障診斷中的應用
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
基于LCD和排列熵的滾動軸承故障診斷
主站蜘蛛池模板: 免费国产不卡午夜福在线观看| 欧美不卡在线视频| 91在线播放国产| 亚洲视屏在线观看| 青青草国产免费国产| 欧美丝袜高跟鞋一区二区| 亚洲国产精品久久久久秋霞影院| 亚洲最大福利视频网| 91娇喘视频| 国产自无码视频在线观看| 亚洲性日韩精品一区二区| 99热最新网址| 中文字幕在线看| 国产探花在线视频| 欧美区一区| 99久久国产精品无码| 伊人五月丁香综合AⅤ| 伊人久久婷婷五月综合97色| 国产无码制服丝袜| 中文字幕 91| 看看一级毛片| 2021天堂在线亚洲精品专区| 国产成人免费高清AⅤ| 少妇露出福利视频| 国产成人无码AV在线播放动漫 | 亚洲专区一区二区在线观看| 毛片三级在线观看| 四虎国产永久在线观看| 在线国产三级| 国产精品亚洲天堂| 日韩一级二级三级| 9啪在线视频| a级毛片免费看| 国产第一页屁屁影院| 2021国产在线视频| 亚洲无码高清视频在线观看| 老熟妇喷水一区二区三区| 青青草原国产av福利网站| 无码中文字幕乱码免费2| 久久成人国产精品免费软件 | a级毛片免费网站| 一级毛片中文字幕| 午夜国产大片免费观看| 亚洲男人天堂2020| 欧美成人第一页| 中文一区二区视频| 亚洲国产欧美国产综合久久| 91福利免费| 亚洲国产91人成在线| 久热中文字幕在线观看| 青草视频在线观看国产| 久久精品这里只有精99品| 99热这里只有精品免费国产| 久久综合九色综合97网| 制服丝袜在线视频香蕉| 久久网欧美| 91久久夜色精品| www中文字幕在线观看| 91青青草视频| 免费毛片全部不收费的| 日本欧美精品| 国产美女视频黄a视频全免费网站| 一本大道香蕉中文日本不卡高清二区| 亚洲日韩欧美在线观看| 毛片最新网址| 中文字幕永久视频| 亚洲精品成人7777在线观看| 九九视频免费在线观看| 欧美成人手机在线观看网址| 欧美日韩成人| 亚洲欧美日韩另类| 1769国产精品免费视频| 国产精品污视频| 98超碰在线观看| 三级欧美在线| 性色生活片在线观看| 日韩欧美国产综合| 夜色爽爽影院18禁妓女影院| 亚洲欧洲日韩久久狠狠爱| 日韩成人午夜| 最新国产成人剧情在线播放| 最新国产在线|