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

基于諧波小波和去趨勢波動分析的摩擦振動信號研究

2017-08-30 12:22:29李精明魏海軍魏立隊楊智遠
振動與沖擊 2017年15期
關鍵詞:趨勢振動信號

李精明, 魏海軍, 魏立隊, 楊智遠, 劉 沖, 劉 竑

(1. 上海海事大學 商船學院, 上海 201306; 2. 大連海事大學 輪機工程學院, 大連 116026)

基于諧波小波和去趨勢波動分析的摩擦振動信號研究

李精明1,2, 魏海軍1, 魏立隊1, 楊智遠1, 劉 沖1, 劉 竑1

(1. 上海海事大學 商船學院, 上海 201306; 2. 大連海事大學 輪機工程學院, 大連 116026)

為實現摩擦振動信號的降噪和摩擦振動信號特征提取,在往復式摩擦磨損試驗機上進行了摩擦副摩合磨損試驗。應用諧波小波對獲得的非線性、非平穩的摩擦振動信號進行分解,實現摩擦振動信號的降噪。應用去趨勢波動分析算法對摩擦振動信號進行分析,獲得不同階數下的Hurst指數,判別數據序列的屬性及其趨勢增強的程度。研究結果表明,隨著磨合磨損試驗的進行,摩擦振動信號的標度指數呈現逐漸增大的變化趨勢,去趨勢波動分析算法能夠實現摩擦振動信號特征提取,摩擦振動信號的標度指數變化能夠用于摩擦副的磨合磨損狀態監測和識別。

諧波小波; 去趨勢波動分析; 標度指數; 摩擦振動

摩擦振動是機械運動摩擦副在摩擦磨損過程中產生的現象,蘊含著許多反映系統摩擦學特征和摩擦狀態的信息[1]。摩擦學系統的輸出信息包括摩擦振動、摩擦力矩、摩擦因數、磨粒形貌、磨損表面形貌[2]等。許多學者從這些摩擦學參數來研究摩擦副摩合磨損狀態以及摩擦副相關的機械故障診斷,但通過摩擦力矩、摩擦因數、磨損表面形貌來提取摩擦學特征存在一定的困難[3]。而通過磨粒形貌來分析則程序繁瑣和低效,且分析結論與研究人員的經驗有密切關系[4]。機械設備的摩擦振動信號可以通過加速度傳感器獲得,可在機械設備運行的情況下實時在線采集,因此可以通過提取機械設備的摩擦振動信號的特征來實時監測機械設備摩擦副的狀態。亦可通過分析摩擦副相關機械各狀態摩擦振動信號的特征,實現摩擦副相關機械的故障診斷。

諧波小波(Harmonic Wavelet, HW)理論[5]是Newland于1993年提出的信號處理方法。諧波小波是一種復小波,在頻域緊支,有明確的函數表達式,諧波小波函數具有頻域盒形緊支譜特性及良好的相位定位能力,在信號分解過程中數據信息量不變,頻域成分不相交,在頻域具有良好的細化能力,可以將信號細分到任意頻段,廣泛應用于信號去噪[6],機械故障診斷[7]地震波分析[8]等領域。基于分形自相似理論發展起來的去趨勢波動分析(Detrended Fluctuation Analysis, DFA)是Peng等[9]生物物理學家于1994年探測DNA內部分子鏈的相關可能性程度時首先提出的。通過對信號累積和的分段去趨勢運算來突出信號中的弱相關成分,能夠檢測非平穩信號中的長程相關信息,被成功應用到信號去噪[10]、故障診斷[11]、醫學數據分析[12]等領域,取得很好的效果。

通過摩擦振動信號來進行摩擦副狀態識別和摩擦磨損行為研究的核心是摩擦振動信號的特征提取。有學者應用時頻譜圖[13]、諧波小波、奇異值分解[14]、混沌吸引子[15]、多重分形[16]等提取摩擦振動信號的特征,取得了較好的效果。本文針對以船用柴油機缸套和活塞環為材質的摩擦副摩擦磨損試驗過程中獲得的摩擦振動信號進行諧波小波分解,重構獲得降噪的摩擦振動信號。應用去趨勢波動分析算法對降噪后的摩擦振動信號進行分析,得到去趨勢波動分析譜圖及不同階數下Hurst指數,實現特征參數對摩擦副摩擦磨損狀態的表征,為基于摩擦振動信號的摩擦副狀態識別和摩擦磨損行為研究提供了新的方法。

1 分析方法

1.1 諧波小波算法

摩擦振動信號振幅小、能量微弱,常規的短時傅里葉變換、魏格納維爾分布會因為摩擦振動信號的譜峰微弱而不易分辨。諧波小波(Harmonic Wavelet, HW)是Newland提出的信號處理方法,具有頻域盒形緊支譜特性及良好的相位定位能力,能有效提取信號中的奇異成分。圖1為信號采用諧波小波變換后的頻域分布。作為諧波小波變換的拓展,諧波小波包變換可以將信號分解到任意頻段,獲得信號所有的低頻分量和高頻分量,準確地顯現微弱信號。

在諧波小波分析過程中,設信號的最高分析頻率為fh,則第j(j=1, 2, 3,…)層頻段的帶寬B可以表示為:

B=2-jfh

(1)

s為頻段數,分析頻帶的上限m和下限n分別為:

m=sB

n= (s+1)Bs=0,1,2,…,2j-1

(2)

根據諧波小波包變換原理,摩擦振動信號的諧波小波包變換過程如下:

(1) 根據摩擦振動信號的特點確定諧波小波包的分解層數,確定頻帶的上下限m,n;

(2) 計算信號的離散時間傅里葉變換;

(3) 計算m,n確定尺度下的頻域;

(4) 求離散傅里葉逆變換FFT,對所需頻段的信號進行時域分析。

1.2 去趨勢波動分析理論及算法

DFA是一種計算長程相關性的方法,能定量地描述時間序列標度不變性,可用于非靜態、非平穩數據分析。與重標極差法(Rescaled range analysis)相比[17],DFA 分析法消除了序列的局部趨勢,避免了將時間序列的短程相關、非平穩性虛假地檢測為長程相關性。

圖1 諧波小波包分解的頻域分布圖

對于一個非平穩時間序列{xi}(i=1,2,…,N),應用DFA方法分析的主要步驟如下:

(1) 計算時間序列{xi}(i=1,2,…,N)的離差序列Y(i):

(3)

(2) 把Y(i)等分成Ns個不重疊的等時間長度s的區間,其中Ns=[N/s](即取整數)。由于序列長度并不總是時間長度s的倍數,使得有小部分序列數據剩余. 因此,對Y(i)的進行逆序的同樣操作,這樣共有2Ns個等長度的區間。

(3) 對每個區間v,用最小二乘法擬合數據,得到局部趨勢,擬合多項式的階數可以是一次的、二次的、更高次的,擬合的階數反映了“趨勢”被消除的程度,階數越高,“趨勢”消除的效果越好,但計算時間也相應增加。濾去該趨勢后的時間序列記為Ys(i),表示原序列與擬合值之差,即:

Ys(i)=Y(i)-Pv(i)i=1,2,…,N

式中,Pv(i)是第v區間的擬合多項式。擬合多項式采用線性、二次、三次,可以是更高階m的多項式,則分別記為DFA1,DFA2,DFA3,…,DFAk等。顯然k階的DFA 濾去了累積離差中的m階趨勢成分以及原始序列中的k-1階趨勢成分。

(4) 計算每個區間濾去趨勢后的方差,以區間(v=1,2,…,2Ns)為例,進行k階多項式擬合:

yv(i)=a1ik+a2ik-1+…+ak+1k=1,2,…

(4)

對于區間(v=1,2,…,Ns):

(5)

對于區間(v=N+1,N+2,…,2Ns)

(6)

(5) 確定數據序列的q階波動函數:

(7)

(6) 通過分析雙對數坐標圖s~Fq(s)的關系,確定波動函數的標度指數h(q),即存在冪率關系:

Fq(s)∝sh(q)

(8)

不同階數q下的Hurst指數即標度指數h(q),可考察數據DFA波動函數的標度行為,通過標度不變性來刻畫時間序列的長程相關特性,即可通過數據序列不同階數q下的Hurst指數數值來判別數據序列的屬性及其趨勢增強的程度。

1.3 標度指數的含義

標度指數h(q)存在于一定標度區間,可用于表征數據序列的相關性,能將時間序列區分為隨機序列與非隨機序列。q為2時,對于一個時間序列[18]:

當0

當h(q)=0.5時,標志著所研究的數據序列是一個隨機序列,即過去的增量與未來的增量不相關,序列具有標度不變性。

當0.50.5的程度。

2 實驗部分

2.1 試驗方法與材料

摩擦副摩擦磨損試驗的設備采用CFT-I型摩擦磨損試驗機,圖2為該試驗機的原理。

圖2 CFT-I型摩擦磨損試驗機原理圖

銷試樣和盤試樣為采用線切割方法從船用柴油機的活塞環和缸套上截取。銷試樣固定在主軸下端部,盤試樣固定在臺架上,電機驅動的偏心輪連桿機構使臺架和盤試樣作往復運動。載荷由加載機構經銷試樣施加于盤試樣,摩擦振動信號采用美國壓電公司生產的356A16型三軸加速度傳感器來測量,靈敏度為100 mV/g,量程為±50 g,檢測位置在盤試樣下方。銷試樣材質為合金鑄鐵,主要成分為Fe、C、Si、Mo、P等,硬度HV600~680,原始表面粗糙度Ra=0.687 μm,截面尺寸3 mm×4 mm。盤試樣接觸表面為弧面,材質為合金鑄鐵,主要成分為Fe、C、Si、Mn、P等,硬度HV320~420,原始表面粗糙度Ra=1.641 μm,尺寸φ30 mm×10 mm。施加的試驗力為50 N,驅動電機轉速為600 r/min,銷試樣-盤試樣相對運動的平均速度為0.1 m/s,試驗時間為480 min。潤滑方式為滴油潤滑,潤滑介質為Mobilgard 412船用潤滑油,黏度141 cSt(40 ℃)。

2.2 摩擦振動信號的采集

摩擦振動信號的采集裝置為比利時LMS公司生產的SCADAS(Supervisory Control and Data Acquisition System)型前端數據采集系統,采樣頻率25 600 Hz,采樣點數4 096。采集的數據每0.16 s自動生成一個文本文件,存入計算機。

圖3 摩擦振動信號的時域波形

圖3為摩擦磨損試驗過程中,傳感器在初期、中期、末期采集的摩擦振動信號時域波形,即對于試驗獲得的數據,在第1 min、240 min、480 min附近各選取一個0.16 s的數據。從圖3可以看出,所獲得的摩擦振動信號是非線性、非平穩信號,信號波動復雜,無規律可循,微弱的摩擦振動信號堙沒于背景噪聲之中,需要去噪才能正確提取摩擦振動信號特征。

3 試驗數據分析

3.1 摩擦振動信號的諧波小波降噪

應用諧波小波實現摩擦振動信號的降噪,文獻[5]指出摩擦振動信號的特點是振幅小、頻率高,其中蘊含著有價值的信息。文獻[14]應用諧波小波包變換對摩擦振動信號進行降噪,提取能夠反映系統摩擦學特征和摩擦狀態的摩擦振動信號,指出低頻信號不含有反映摩擦磨損狀態的摩擦振動信號,驗證5 000~6 000 Hz的摩擦振動信號能夠反映摩擦副摩擦磨損狀態。

本文應用諧波小波,將采集的摩擦振動信號分解成7層64個頻帶,帶寬度為100 Hz,重構第51~60頻帶,降噪后的摩擦振動信號時域波形見圖4。從圖4可以看出,降噪后摩擦振動信號幅值顯著減小,周期成分減弱,摩擦振動的沖擊信息清晰出現。隨著磨合磨損試驗的進行,摩擦振動信號的振幅呈現減小的變化趨勢。

圖4 重構的摩擦振動信號時域波形

3.2 摩擦振動信號的去趨勢波動分析

應用去趨勢波動分析算法分析去噪后的摩擦振動信號,對試驗初期的摩擦振動信號,當尺度s取16,階數q分別取-3,-1,1,3時,分析得到的波動函數如圖5所示。從圖中可以看出,不同階數對信號的小振幅和大振幅有不同的分辨能力。階數為負時,波動函數敏感地反映摩擦振動信號的小振幅變化,階數為正時,波動函數敏感地反映摩擦振動信號的大振幅變化。可以通過不同階數下的波動函數值來考察摩擦振動信號的振幅分布情況。

圖5 摩擦振動信號的q階波動函數

圖6(a)、6(b)、6(c)所示為摩擦磨損試驗初期、中期、末期摩擦振動信號應用去趨勢波動分析,在不同尺度s、不同階數q下的波動函數值,在某一階數q下,波動函數的回歸線斜率即標度指數h(q)。從圖6(a)、6(b)、6(c)可以看出,不同階數q下的波動函數值在小尺度較大尺度更為明顯。由于小尺度能夠更好地嵌入數據的局部周期,因此小尺度能夠更明顯地區分數據局部周期的大波動和小波動。隨著磨合磨損試驗的進行,圖6(a)、6(b)、6(c)呈現出漸進的變化規律,表現在波動函數的回歸線斜率的變化。圖6(d)為摩擦磨損試驗初期、中期、末期摩擦振動信號在不同階數q下的Hurst指數即標度指數h(q),隨著磨合磨損試驗的進行,標度指數h(q)呈現增大的趨勢。

(a)(b)

(c)(d)

圖6 摩擦振動信號的冪率關系譜圖

Fig.6 Power-law relationship of frictional vibration signals

表1為摩擦振動特征信號的不同階數q下的Hurst指數即標度指數h(q)。從表1可以看出,隨著磨合磨損試驗的進行,標度指數h(q)呈現增大的趨勢變化,表明摩擦磨損試驗獲得的摩擦振動信號是反持久性數據序列,摩擦振動未來的趨勢表現為與之前的相反。摩擦振動信號回歸線的斜率呈現逐漸增大的變化,表明反持續性逐漸較小,數據序列的突變跳躍逆轉性趨于緩和,表征摩擦副的摩擦磨損逐漸趨于穩定。

表1 摩擦振動特征信號的q階Hurst指數

3.3 機理分析

摩擦振動信號標度指數的變化與摩擦副磨合磨損過程中磨損表面變化密切相關。摩合磨損試驗摩擦副磨合初期,摩擦副表面粗糙度較大,摩擦副往復運動過程中產生的能量較大,摩擦振動信號劇烈,摩擦振動的主峰明顯,平均振幅較大,摩擦副磨損過程處于不穩定的狀態。隨著磨合磨損試驗的進行,摩擦副表面粗糙度逐漸減小。磨合磨損試驗初期的摩擦振動信號標度指數數值最小,反持續性最強。隨著磨合磨損試驗的進行,摩擦副表面的粗糙度逐漸下降,激發的摩擦振動也減弱,磨合過程摩擦表面消耗的能量減弱,摩擦磨損逐漸趨于穩定,標度指數數值增大。第240 min時,盤試樣粗糙度從最初的1.641 μm降至1.263 μm,激發的摩擦振動強度減小,振幅明顯減小,主峰減少,摩擦振動信號的變化表明磨合過程損耗的能量減小,摩擦磨損趨于穩定。在磨合磨損試驗的中后期,標度指數進一步增大,摩擦磨損逐漸進入穩定磨損階段,盤試樣表面粗糙度進一步減小為1.093 μm,摩擦副表面粗糙度減小的幅度減弱,表面磨損輕微,激發的摩擦振動減弱。摩擦副進入穩定磨損狀態,摩擦振動主峰不明顯,平均振幅減小,摩擦振動標度指數數值增大,反持續性減弱。

上述分析表明,標度指數能夠參數化摩擦振動信號的特征。摩擦振動信號去趨勢波動分析冪率關系譜圖及其參數能體現摩擦振動的特征,反映摩擦副摩合磨損過程中所處的摩擦振動狀態。

4 結 論

摩擦振動蘊含著反映系統摩擦學特征和摩擦狀態的信息,本文利用諧波小波和去趨勢波動分析研究摩擦副摩合磨損過程中的摩擦振動信號,結論如下:

(1) 應用去趨勢分析算法能有效地分析摩擦振動信號的非線性特征,反映摩擦副摩合磨損過程中摩擦振動漸變過程,可以用摩擦振動信號去趨勢波動分析冪率關系譜圖及其參數表征摩擦振動信號的特征。

(2) 去趨勢分析標度指數能夠參數化摩擦振動信號的特征,隨著磨合磨損試驗的進行,摩擦振動信號的標度指數呈現逐漸增大的變化趨勢。

[1] 葛世榮,朱華. 摩擦學復雜系統及其問題的量化研究方法[J]. 摩擦學學報,2002, 22(5): 405-408.

GE Shirong, ZHU Hua. Complicate tribological system and quantitative study methods of their problems[J]. Tribology, 2002, 22(5): 405-408.

[2] ZHU Hua, GE Shirong, CAO Xucheng, et al. The Changes of fractal dimensions of frictional signals in the running-in wear process[J]. Wear, 2007, 263(7): 1502-1507.

[3] SUN Di, LI Guobin, WEI Haijun, et al. Investigation on frictional vibration behavior of tribological pairs under different wear states[J]. Journal of Tribology, 2015, 137(4): 74-81.

[4] MACIAN V, PAYRI R, TORMOS B, et al. Applying analytical ferrography as a technique to detect failures in diesel engine fuel injection systems[J]. Wear, 2005, 260(2):562-566.

[5] 李國賓,任宗英,王宏志,等. 摩擦振動信號諧波小波包特征提取[J]. 摩擦學學報,2011, 31(5): 452-456.

LI Guobin, REN Zongying, WANG Hongzhi, et al. Characteristics extraction of friction vibration singal using harmonic wavelet packet transforms[J]. Tribology, 2011, 31(5): 452-456.

[6] 李舜酩,徐慶余.微弱振動信號的諧波小波頻域提取[J].西安交通大學學報,2004, 38(2):51-55.

LI Shunming, XU Qingyu. Harmonic wavelet extraction for weak vibration signal in frequency domain[J]. Journal of Xi’an Jiaotong University, 2004, 38(2): 51-55.

[7] 李方,李友榮,王志剛.諧波小波時頻圖在齒輪故障診斷中的應用[J]. 振動與沖擊,2007, 26(3): 128-134.

LI Fang, LI Yourong, WANG Zhigang, Harmonic wavelet time-spectrum plot with applications in gear fault diagnosis[J]. Journal of Vibration and Shock, 2007, 26(3): 128-134.

[8] CECINI D, PALMERI A. Spectrum-compatible accel-erograms with harmonic wavelets[J]. Computers and Structures, 2014, 147(10): 26-35.

[9] ECHEVERRIA C, RODRIGUEZ E, AGUILAR-CORNEJO M, et al. Linear combination of power-law functions for detecting multiscaling using detrended fluctuation analysis[J]. Physica A: Statistical Mechanics and its Applications, 2016, 460(15): 283-293.

[10] 何文平,吳瓊,成海英,等.不同濾波方法在去趨勢波動分析中去噪的應用比較[J]. 物理學報,2011, 60(2): 45-53.

HE Wenping, WU Qiong, CHENG Haiying, et al. Comparison of applications of different filter methods for de-noising detrended fluctuation analysis[J]. Acta Phys. Sin., 2011, 60(2): 45-53.

[11] 李力,彭中笑,彭書志,等.去趨勢波動分析在齒輪故障診斷中的應用研究[J].中國機械工程,2009, 20(19): 2311-2314.

LI Li, PENG Zhongxiao, PENG Shuzhi. Detrended fluctuation analysis for gear fault diagnosis[J]. China Mechanical Engineering, 2009, 20(19): 2311-2314.

[12] LIM J H, KHANG E J, LEE T H, et al. Detrended fluctuation analysis and Kolmogorov-Sinai entropy of electroencephalogram signals[J]. Physics Letters A, 2013, 377(38): 2542-2545.

[13] 黃朝明,于洪亮,關德林,等. 摩擦振動時頻圖像特征提取[J]. 振動與沖擊,2012, 31(7): 45-49.

HUANG Chaoming, YU Hongliang, GUAN Delin, et al. Feature extraction of frictional vibration based on time-frequency image[J]. Journal of Vibration and Shock, 2012, 31(7): 45-49.

[14] 孫迪,李國賓,魏海軍,等. 磨合磨損過程中摩擦振動變化規律研究[J]. 哈爾濱工程大學學報,2015,36(2): 166-170.

SUN Di, LI Guobin, WEI Haijun, et al. Study on variation rules of friction vibration in the process of friction and wear[J]. Journal of Harbin Engineering University, 2015, 36(2):166-170.

[15] SUN Di, LI Guobin, WEI Haijun, et al. Experimental study on the chaotic attractor evolvement of the friction vibration in a running-in process[J]. Tribology International, 2015, 88(3): 290-297.

[16] 李精明,魏海軍,魏立隊,等. 摩擦振動信號的經驗模式分解和多重分形研究[J]. 振動與沖擊,2016, 35(3): 198-203.

LI Jingming, WEI Haijun, WEI Lidui, et al. Empirical mode decomposition and multifractal of frictional vibration signal[J]. Journal of Vibration and Shock, 2016, 35(3): 198-203.

[17] 吳建軍,徐尚義,孫會君. 混合交通流時間序列的去趨勢波動分析[J].物理學報,2011, 60(1):78-85.

WU Jianjun, XU Shangyi, SUN Huijun. Detrended fluctuation analysis of time series in mixed traffic flow[J]. Acta Phys. Sin., 2011, 60(1): 78-85.

[18] KANTELHARDT J W, ZSCHIEGNER S A, KOSCIELN Y, et al. Multifractal detrended fluctuation analysis of nonstationary time series[J]. Physical A, 2002, 316(2): 87-114.

Frictional vibration signals based on harmonic wavelet and detrended fluctuation analysis

LI Jingming1,2, WEI Haijun1, WEI Lidui1, YANG Zhiyuan1, LIU Chong1, LIU Hong1

(1. Merchant Marine College, Shanghai Maritime University, Shanghai 201306, China;2. Marine Engineering College, Dalian Maritime University, Dalian 116026, China)

For the purpose of de-noising and characteristic-extracting of frictional vibration signals,the friction and wear tests of friction pairs were conducted on a reciprocating pin-on-disc tester. The nonlinear and non-stationary frictional vibration signals were decomposed with the harmonic wavelet to realize their de-noising. Then the frictional vibration signals were analyzed with the detrended fluctuation analysis (DFA) to get Hurst exponents under different orders and judge the attribute of data series and their trend enhancement level. The results showed that with friction and wear tests going on, the scaling exponents of the frictional vibration signals reveal a gradual increase trend; the characteristics of the frictional vibration signals can be extracted with the detrended fluctuation analysis; the changing of the scaling exponents of the frictional vibration signals can be used to monitor and identify friction and wear states of friction pairs.

harmonic wavelet; detrended fluctuation analysis; scaling exponent; frictional vibration

國家高技術研究發展計劃(863計劃)(2013AA040203)

2016-01-29 修改稿收到日期:2016-06-12

李精明 男,博士生,講師,1981年5月生

魏海軍 男,博士,教授,博士生導師,1971年7月生

TH117.2

A

10.13465/j.cnki.jvs.2017.15.035

猜你喜歡
趨勢振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
趨勢
第一財經(2021年6期)2021-06-10 13:19:08
信號
鴨綠江(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微分方程的振動性
初秋唇妝趨勢
Coco薇(2017年9期)2017-09-07 21:23:49
SPINEXPO?2017春夏流行趨勢
基于LabVIEW的力加載信號采集與PID控制
主站蜘蛛池模板: 亚洲乱伦视频| 国产麻豆精品久久一二三| 日本黄色a视频| 久久久久88色偷偷| 在线观看免费人成视频色快速| 999国产精品| 91成人在线免费视频| 中文字幕在线一区二区在线| 国产啪在线91| 日韩精品成人网页视频在线| 97视频免费在线观看| 国产一级在线观看www色| 香蕉久人久人青草青草| 中文国产成人久久精品小说| 久久久精品无码一区二区三区| 亚洲成人免费在线| 天天色天天综合| 国产va在线观看| 91久久偷偷做嫩草影院免费看| 国产欧美日韩视频一区二区三区| 免费看a级毛片| 精品一区二区三区视频免费观看| 欧美日韩国产高清一区二区三区| 成人在线天堂| 青青草原偷拍视频| 1769国产精品视频免费观看| 精品国产免费观看一区| 亚洲婷婷六月| 99成人在线观看| 99精品免费在线| 欧美曰批视频免费播放免费| 亚洲成a人在线观看| 无码又爽又刺激的高潮视频| 欧美色视频日本| 一区二区三区成人| 欧美国产综合色视频| 亚洲自偷自拍另类小说| 亚洲AⅤ波多系列中文字幕| 国产凹凸视频在线观看| av大片在线无码免费| 日韩AV手机在线观看蜜芽| 波多野结衣在线一区二区| 青青草91视频| 中国一级特黄视频| 国产成+人+综合+亚洲欧美| 色香蕉影院| 国产精品福利导航| 国产高清在线观看91精品| 全部免费毛片免费播放| 亚洲一级毛片在线观播放| 女同久久精品国产99国| 亚洲激情区| 九色在线观看视频| 国产欧美日韩精品综合在线| 一本一道波多野结衣一区二区| 蝴蝶伊人久久中文娱乐网| av色爱 天堂网| 99在线观看视频免费| 亚洲天堂网站在线| 国产成人精品一区二区| 欧洲熟妇精品视频| 亚洲性视频网站| 五月婷婷激情四射| A级毛片无码久久精品免费| 国产免费久久精品99re丫丫一| 欧美亚洲日韩不卡在线在线观看| 精品99在线观看| 激情亚洲天堂| 她的性爱视频| 久久久久久久久18禁秘| 久久婷婷人人澡人人爱91| 黄片一区二区三区| 国产精品第页| 一本一本大道香蕉久在线播放| 国产激情国语对白普通话| 伊人无码视屏| 国产欧美精品专区一区二区| 久久亚洲国产一区二区| 中文字幕乱码二三区免费| 亚洲IV视频免费在线光看| 无码精品国产VA在线观看DVD| 1级黄色毛片|