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

基于EMD-SVM的鈦合金銑削過程刀具磨損監(jiān)測*

2022-11-04 08:37:32謝振龍岳彩旭劉獻(xiàn)禮嚴(yán)復(fù)鋼劉智博穆殿方梁越昇
振動、測試與診斷 2022年5期
關(guān)鍵詞:信號模型

謝振龍,岳彩旭,劉獻(xiàn)禮,嚴(yán)復(fù)鋼,劉智博,穆殿方,梁越昇

(1.哈爾濱理工大學(xué)先進(jìn)制造智能化技術(shù)教育部重點實驗室 哈爾濱,150080)

(2.佐治亞理工學(xué)院喬治·W·伍德拉夫機(jī)械工程學(xué)院 亞特蘭大,30332)

引言

智能制造是包括原材料、加工及測量檢測等一系列環(huán)節(jié)的生產(chǎn)過程,刀具是其中非常關(guān)鍵的環(huán)節(jié)[1]。刀具磨損作為最主要的刀具失效形式,關(guān)乎著制造的精度及產(chǎn)品的表面質(zhì)量,如何精確識別刀具磨損已成為一個研究熱點。

近些年,專家們提出了許多刀具磨損識別的方法,主要分為直接觀察法和間接觀察法。直接觀察法就是測量刀具磨損量來衡量刀具磨損程度,包括監(jiān)測磨損寬度、磨損面大小等。常用的方法有接觸法、輻射法和光學(xué)檢測法。它存在2個主要缺陷:①不能對刀具磨損狀態(tài)進(jìn)行實時監(jiān)測;②停機(jī)檢測極大降低工作效率。間接觀察法就是通過檢測刀具銑削時產(chǎn)生的信號,構(gòu)建基于信號的刀具磨損模型,從而監(jiān)測刀具磨損的方法。此方法不對加工過程造成干擾,且可以連續(xù)監(jiān)測加工過程,更適合于在線監(jiān)測。常用的方法有切削力監(jiān)測[2]、聲發(fā)射監(jiān)測[3]、振動監(jiān)測[4]、電流與功率監(jiān)測[5]、超聲波監(jiān)測[6]和溫度監(jiān)測[7]。然而單純采集銑削時的信號不能準(zhǔn)確反應(yīng)刀具狀態(tài),機(jī)器學(xué)習(xí)及信號處理技術(shù)為刀具磨損識別提供了新的思路。Rizal等[8]提出了一種結(jié)合多傳感器信號和馬氏田口系統(tǒng)的檢測刀具磨損的方法,并提出馬氏距離作為識別刀具磨損的評價指標(biāo)。Bhuiyan等[9]對切削時的聲發(fā)射及振動信號進(jìn)行研究,探究刀具磨損與信號之間的關(guān)系。Wu等[10]通過多信號融合探究刀具剩余壽命,揭示信號與刀具磨損之間的關(guān)系。Chen等[11]通過深度置信網(wǎng)絡(luò)學(xué)習(xí)模型融合刀具力信號、振動信號及聲發(fā)射信號對刀具磨損進(jìn)行識別。謝峰云等[12]提出的廣義BP神經(jīng)網(wǎng)絡(luò)識別模型,提取了振動信號廣義均方根、廣義功率譜密度均方根及小波包系數(shù)均方根為特征值,將特征值輸入BP神經(jīng)網(wǎng)絡(luò)中識別鋁合金銑削時的顫振。甘梓舜等[13]通過對機(jī)床的各種信號進(jìn)行監(jiān)測從而監(jiān)測刀具磨損。陳剛等[14]通過BP神經(jīng)網(wǎng)絡(luò)對振動及切削力信號進(jìn)行學(xué)習(xí)從而識別刀具磨損。但是由于深度學(xué)習(xí)及其他學(xué)習(xí)模型所需數(shù)據(jù)量大、運算速度慢,BP神經(jīng)網(wǎng)絡(luò)準(zhǔn)確性相對不穩(wěn)定,而SVM由于其算法簡單以及具有優(yōu)秀的“學(xué)習(xí)”能力[15],所以筆者采用SVM對數(shù)據(jù)進(jìn)行學(xué)習(xí),從而達(dá)到刀具磨損識別的效果。

利用信號對刀具磨損狀態(tài)識別會面臨一個問題,即采取的信號往往摻雜著噪聲,所以對信號進(jìn)行處理至關(guān)重要。Plaza等[16]通過對比不同的信號處理方法,選擇小波包與神經(jīng)網(wǎng)絡(luò)結(jié)合的方式對信號進(jìn)行處理監(jiān)測表面質(zhì)量,但是由于小波包基有選擇困難的局限,選用小波包基會花費較長的時間且監(jiān)測的信號會出現(xiàn)非線性、非平穩(wěn)的現(xiàn)象。為了解決這個問題,Huang等[17]提出了經(jīng)驗?zāi)B(tài)分解,該方法對于分析信號非線性與非平穩(wěn)具有明顯的優(yōu)越性,同時克服了小波包基選擇困難的特點。在處理信號時選擇EMD將會簡化數(shù)據(jù)的計算且不會失去數(shù)據(jù)原有的特性,因此將EMD引用到信號監(jiān)測中對信號進(jìn)行處理。

由于采取的信號數(shù)據(jù)量大,直接輸入信號進(jìn)行學(xué)習(xí)及判斷會增加模型運算時間,所以應(yīng)當(dāng)提取信號特征值來對數(shù)據(jù)進(jìn)行簡化。刀具磨損時,加工過程中的時、頻域和能量分布將發(fā)生變化。為了有效地監(jiān)測信號的變化,選取的特征值應(yīng)當(dāng)囊括以上特征。標(biāo)準(zhǔn)差可以反映信號能量的變化,它的值隨信號幅度的增大而增大。功率譜熵[18]是一個無量綱指標(biāo),其可以反映不同頻率在頻帶內(nèi)的分布。當(dāng)頻率分量廣泛分布在頻帶上時,分布的不確定度最大,其所對應(yīng)的功率譜熵最大;相反,當(dāng)頻率分量分布在一定頻帶上時,頻率分布的不確定性最小,其所對應(yīng)的功率譜熵最小。I-kazTM是一種多分辨率分析方法,能準(zhǔn)確得出信號的變化[19]。因此,將標(biāo)準(zhǔn)差、功率譜熵和I-kazTM作為信號的特征值。

基于以上分析,筆者針對鈦合金切削過程中的刀具磨損識別問題,提出了一種基于EMD及SVM刀具磨損階段識別方法。首先,將原始加速度信號及力信號分解為一系列IMF,選擇有用的IMF來組合一個新的信號;其次,計算新信號的特征值,將得到的特征值矩陣作為SVM的輸入;最后,得到了刀具磨損識別模型,能對刀具磨損階段進(jìn)行準(zhǔn)確識別。

1 信號處理及特征值的選取

1.1 EMD分解思想

EMD方法假設(shè)任何信號都由不同IMF組成,每個IMF可以是線性的,也可以是非線性的。IMF分量必須滿足2個條件:①其極值點個數(shù)和過零點數(shù)相同或最多相差1個;②其上下包絡(luò)關(guān)于時間軸局部對稱。

EMD分解過程基于以下假設(shè):

1)信號最少有2個極值,即1個極大值和1個極小值;

2)時域特性由極值間隔決定;

3)如果數(shù)據(jù)序列完全缺乏極值但僅包含拐點,那么它也可以通過1次或多次求導(dǎo)來表示極值點,而最終結(jié)果可以由這些成分求積分來獲得。

具體方法是由一個“篩選”過程完成。

首先,找出信號s(t)所有的極大值點并將其用3次樣條函數(shù)擬合成原數(shù)據(jù)序列上的包絡(luò)線,再找出所有的極小值點并將其用3次樣條函數(shù)擬合成原數(shù)據(jù)[20]。

其次,計算上下包絡(luò)線的均值,記為m1(t),并將原數(shù)據(jù)序列s(t)減去該均值即可得到1個去掉低頻的新數(shù)據(jù)序列h1

因為h1(t)一般不是1個IMF分量序列,為此需要對它重復(fù)進(jìn)行上述處理過程k次,直到h1(t)符合IMF的定義要求:所得到的均值趨于零。這樣就得到了第1個IMF分量c1(t),它代表信號s(t)中最高頻率的分量

然后,將c1(t)從s(t)中分離出來,即得到1個去掉高頻分量的差值信號r1(t)。將r1(t)作為原始數(shù)據(jù),重復(fù)以上步驟,得到第2個IMF分量c2(t),重復(fù)n次,得到n個IMF分量[21]

最后,當(dāng)cn(t)或rn(t)滿足給定的終止條件(通常rn(t)成為一個單調(diào)函數(shù))時,循環(huán)結(jié)束。由式(3)和式(4)得到

1.2 信號特征值的選取

1.2.1 I-kazTM

I-kazTM是由I-kaz指數(shù)演化出來的,其主要思想是將動態(tài)信號分解為3個頻率范圍,其中:低頻(LF)范圍為0~0.25fmax;高頻(HF)范圍為0.25fmax~0.5fmax;極高頻(VF)范圍為大于0.5fmax。為了測量數(shù)據(jù)分布的散射,計算了每個頻帶的方差σ2

由于I-kazTM方法是基于數(shù)據(jù)質(zhì)心的數(shù)據(jù)發(fā)散概念發(fā)展起來的,所以I-kazTM的系數(shù)可以用Z∞來表示

將式(6)~(8)代入式(9),得到

式(10)可以用峰度和標(biāo)準(zhǔn)差表示出來,峰度K的公式為

其中:N為數(shù)據(jù)的數(shù)量;s為標(biāo)準(zhǔn)偏差。

因此,根據(jù)峰度K和標(biāo)準(zhǔn)偏差s可得

其中:KL,KH和KV分別為LF,HF和VF范圍中的信號的峰度;sL,sH,sV分別為LF,HF和VF范圍內(nèi)的信號的標(biāo)準(zhǔn)差。

根據(jù)I-kazTM來提取適應(yīng)EMD分解的方法,其不同于I-kazTM,即信號不需要分解為三差分。EMD處理完信號后,選擇有效的模態(tài)分量,分別計算其峰度及標(biāo)準(zhǔn)偏差,最終計算I-kazTM系數(shù)。因此,I-kazTM系數(shù)可以寫成如下形式

其中:N為有效分解模態(tài)的數(shù)目;Ki為分解后模態(tài)的峰度;si為第i個分解信號的標(biāo)準(zhǔn)偏差。

1.2.2 功率譜熵

功率譜熵(power spectral entropy,簡稱PSE)是“信息熵”在頻域中的擴(kuò)展,其與頻率分量的分布有關(guān)。具體計算步驟如下。

1)根據(jù)式(14)獲得信號x(t)的功率譜

其中:N為數(shù)據(jù)長度;X(w)為X(t)的傅里葉變換。

2)功率譜的概率密度函數(shù)可以通過對所有頻率分量的歸一化來計算

其中:s(fi)為頻率分量fi的光譜能量;pi為相應(yīng)的概率密度;N為快速傅里葉變換中頻率分量的總數(shù)。

3)相應(yīng)的功率譜熵定義為

為了比較不同的工作條件,結(jié)果由因子logN標(biāo)準(zhǔn)化,得

功率譜熵E是在[0,1]的范圍內(nèi)的無量綱指示符,其中1表示頻率分量分布不確定度最大,0表示分布不確定度最小。

1.2.3 標(biāo)準(zhǔn)差

標(biāo)準(zhǔn)差(standard deviation,簡稱SD)能很客觀準(zhǔn)確地反映一組數(shù)據(jù)的離散程度,其計算公式為

2 基于EMD-SVM刀具磨損識別模式

SVM是一種二分類模型,其基本模型是定義在特征空間上的間隔最大的線性分類器[22]。SVM學(xué)習(xí)的基本想法是求解能夠正確劃分訓(xùn)練數(shù)據(jù)集并且?guī)缀伍g隔最大的分離超平面,其將空間的點進(jìn)行劃分并根據(jù)劃分的直線對點進(jìn)行識別,如圖1所示。對于線性可分的數(shù)據(jù)來說,這樣的超平面有無窮多個(即感知機(jī)),但是幾何間隔最大的分離超平面卻是唯一的。由于信號特征值相對于信號整體來說,數(shù)據(jù)量小,易于學(xué)習(xí),所以采用SVM對信號進(jìn)行學(xué)習(xí)。

圖1 支持向量機(jī)Fig.1 Support vector machine

在選擇核函數(shù)時,若對給出的數(shù)據(jù)沒有先驗知識,徑向基函數(shù)(radial basis function,簡稱RBF)就是最好的選擇,而且RBF核的支持向量機(jī)可以獲得非常平滑的估計[23],所以筆者采用高斯核函數(shù)作為SVM的核函數(shù)。

為了消除噪聲和其他無關(guān)的信號,首先,利用EMD將加速度信號分解為多個IMF,將分解完的信號重構(gòu)為新信號;其次,將新信號劃分成一系列段,分別計算每個分段的I-kazTM、功率譜熵及均方根,得到銑削狀態(tài)的特征向量矩陣;最后,以特征向量矩陣為輸入特征,建立以SVM為基礎(chǔ)的刀具磨損識別模型。根據(jù)該模型可以判斷刀具磨損階段。刀具磨損階段識別的流程如圖2所示。

圖2 基于EMD-SVM的刀具磨損階段識別流程Fig.2 Process of tool wear stage based on EMD-SVM

3 鈦合金銑削實驗及信號分析

3.1 實驗設(shè)備

實驗在三軸機(jī)床VDL-1000E進(jìn)行,實驗工件材料為Ti-6Al-4V(TC4),實驗刀具為無涂層硬質(zhì)合金四齒銑刀,刀具螺旋角、前角、后角分別為75°,8°和9°。采用PCB加速度傳感器及Kister9171A旋轉(zhuǎn)式測力儀對加速度信號及銑削力信號進(jìn)行采集。采用東華DHDAS動態(tài)信號采集分析系統(tǒng)對加速度信號及力信號進(jìn)行處理,采樣頻率為1 kHz。圖3所示為銑削加工實驗現(xiàn)場,旋轉(zhuǎn)測力儀安裝在刀柄處對加工實驗的3向力信號進(jìn)行測量,將加速度傳感器連接在工件上對加速度信號進(jìn)行測量,選取有效的信號點作為2種信號的共同起始點。

圖3 銑削加工實驗現(xiàn)場Fig.3 The milling experiment site

3.2 實驗方案及數(shù)據(jù)分析

3.2.1 實驗方案

為了驗證該方法,實驗準(zhǔn)備獲取在不同刀具磨損情況下的力及加速度信號,實驗參數(shù)如表1所示。

表1 實驗參數(shù)Tab.1 Experimental parameters

如果要對刀具磨損階段進(jìn)行劃分,就要先確定刀具磨損曲線。在本研究實驗參數(shù)下進(jìn)行實驗,取3次走刀的平均磨損值作為1次磨損取值,建立如圖4所示的刀具磨損曲線。表2為刀具磨損狀態(tài)分類表,根據(jù)磨損率變化情況將刀具磨損過程分為3個階段,即初期磨損階段[0 mm,0.12 mm)、正常磨損階 段[0.12 mm,0.17 mm)和 急 劇 磨 損 階 段[0.17 mm,0.30 mm)。預(yù)先將3個磨損階段標(biāo)記為1,2,3以輸入SVM中。

圖4 刀具磨損過程Fig.4 Tool wear process

表2 刀具磨損狀態(tài)分類表Tab.2 Classification of tool wear status

3.2.2 特征值提取及數(shù)據(jù)分析

在實驗中很難根據(jù)原始的加速度及力信號來判斷銑削狀態(tài)。特征提取的目的是為了減少原始加工信號的尺寸,同時保持所提取的特征中具有刀具狀態(tài)的相關(guān)信息。本節(jié)主要以加速度信號為例進(jìn)行數(shù)據(jù)分析。

選取1組銑削時檢測的初期振動信號,經(jīng)過頻域變換處理后的加速度信號及頻域分布見圖5。

圖5 加速度信號及其頻域Fig.5 Acceleration signal and its frequency domain

圖5所示的加速度信號經(jīng)EMD分解得到的時域和頻域如圖6所示。圖6(a)中,從左到右依次為IMF1,IMF2,···,IMF6。從圖6(b)可以看出,第4~6個模態(tài)固有函數(shù)的頻域極其相近且主頻過小,此時的模態(tài)固有函數(shù)為虛假分量,沒有實際意義,可以忽略不計,所以計算時值采取IMF1~I(xiàn)MF3。圖中為無綱量單位。

圖6 EMD分解時域及頻域圖Fig.6 EMD decomposition time domain and frequency domain diagram

圖7為3個刀具磨損階段的力信號,刀具磨損量分別為0.05,0.15及0.25 mm。由圖可知,當(dāng)后刀面磨損處于正常磨損時切向力及徑向力分別為270和160 N,而當(dāng)后刀面磨損處于急劇磨損階段時,2種力信號的幅值急劇增加,比例分別為41%及27%,這是因為刀具后刀面磨損后,作用在后面的法向力及摩擦力都增大,故切向力及徑向力增加,而軸向力卻未有較大的變化。

圖7 不同階段力變化圖Fig.7 Force variations at different stages

圖8為刀具在初期磨損、中期磨損及急劇磨損的特征值隨時間的變化圖,其中特征值為無量綱單位。由圖可以看出,在不同的刀具磨損情況下,刀具磨損的特征值趨于一個穩(wěn)定值,且3個階段的特征值不同,由此可以得出選取的特征值對信號變化敏感,可以用于模型的學(xué)習(xí)。

圖8 不同階段加速度信號及其特征值變化圖Fig.8 Acceleration signal and its eigenvalue change in different stages

通過對初期、中期及急劇磨損階段的加速度信號進(jìn)行對比發(fā)現(xiàn):急劇磨損階段的加速度信號相較于初期及中期磨損階段來說,均方根及功率譜熵指數(shù)增大而I-kazTM減小,即刀具磨損時信號能量增加,信號的分散性變高。

為了驗證上述結(jié)論的有效性,重建加速度信號的形態(tài)特征、時域能量和頻率分量,對信號進(jìn)行了深入分析。圖9為在不同刀具磨損情況下加速度信號的形態(tài)特征及頻域分布圖,其中頻域分布圖中的頻率分量是由刀齒通過頻率(270 Hz)和諧波組成。對比不同刀具磨損情況下的頻域分布發(fā)現(xiàn):隨著刀具磨損的增加,信號的頻域由低頻向高頻部分轉(zhuǎn)移,高次諧波增加,這是因為切削力隨刀具磨損的增加而增加進(jìn)而引起信號的變化;信號呈現(xiàn)周期性波動,這是刀具切入切出的接觸特性發(fā)生了變化引起的。

圖9 不同刀具磨損情況下的加速度信號及其頻域Fig.9 Acceleration signal and frequency domain under different tool wear conditions

監(jiān)測的特征值信號以加速度信號特征值為主,建立加速度信號與力信號特征值特征矩[I-L,E-L,s-L,I-Z,E-Z,s-Z]作為線性分類器的輸入,其中:I,E,s分別為I-kazTM、功率譜熵和標(biāo)準(zhǔn)差;L,Z分別為力信號和加速度信號。圖10所示為加速度信號構(gòu)成的特征值空間,由于PSE為接近0的數(shù),為了可以直觀地看出特征值在空間的分布,采用1-PSE作為z軸,圖中為無量綱單位。由圖可見,不同刀具磨損情況下的特征值相互獨立且相同刀具磨損階段下的特征值集中在一起。

圖10 特征值空間Fig.10 The feature space

將每個刀具磨損階段的160組數(shù)據(jù)作為訓(xùn)練集,將其余40組數(shù)據(jù)作為測試集。最后,得到1個大小為160的訓(xùn)練集和1個大小為40的測試集的數(shù)據(jù)集對EMD-SVM模式進(jìn)行訓(xùn)練及測試。

為進(jìn)一步驗證建立模型的準(zhǔn)確度及運算速度,使用同1組數(shù)據(jù)作為輸入,建立SVM,BP神經(jīng)網(wǎng)絡(luò)、EMD-SVM及小波包-SVM這4種模型進(jìn)行對比。

對于SVM模型及BP神經(jīng)網(wǎng)絡(luò)模型,將實驗所獲取的信號直接輸入學(xué)習(xí)模型中進(jìn)行訓(xùn)練,然后將測試信號輸入到學(xué)習(xí)好的模型中進(jìn)行判斷。對于小波包-SVM,先將信號進(jìn)行小波處理后輸入SVM構(gòu)建刀具磨損模型,然后將檢測信號輸入模型中進(jìn)行判別。通過對4種模型進(jìn)行訓(xùn)練后得到了4種模型的混淆矩陣,其準(zhǔn)確率對比如圖11所示。其中:混淆矩陣縱坐標(biāo)為預(yù)期判斷的刀具磨損階段;橫坐標(biāo)為實際判斷出的刀具磨損階段;兩者相交的黑色對角線為預(yù)測正確的測試數(shù)據(jù)。表3為4種模型輸入訓(xùn)練信號進(jìn)行訓(xùn)練及輸入檢測信號進(jìn)行識別的運算時間及精度。

圖11 4種模型準(zhǔn)確率對比圖Fig.11 The confusion matrix of three models

表3 不同模型運算時間表Tab.3 Operating time of different models

由于BP網(wǎng)絡(luò)預(yù)測精度不穩(wěn)定,SVM的泛化能力不穩(wěn)定,容易出現(xiàn)過擬合,小波包分解由于挑選小波包基困難和其對非平穩(wěn)及非線性數(shù)據(jù)上的劣勢,訓(xùn)練時間較長;而EMD-SVM由于對信號進(jìn)行預(yù)處理,提取的特征值能更好地表達(dá)樣本特征,不容易出現(xiàn)過擬合現(xiàn)象且計算效率相對于其他模型高,所以模型準(zhǔn)確度較高,泛化能力強,學(xué)習(xí)所用時間較其他機(jī)器學(xué)習(xí)模型短。對比后發(fā)現(xiàn)EMD-SVM的精度最高,運算速度最快。相對于BP神經(jīng)網(wǎng)絡(luò),同1組數(shù)據(jù)的運算時間節(jié)省60%,精度提高了15.12%,對刀具磨損狀態(tài)有很高的靈敏性,能監(jiān)測刀具磨損的發(fā)生。

4 結(jié)論

1)刀具發(fā)生磨損時,加工信號的時頻域發(fā)生變化,頻域會由低頻轉(zhuǎn)向高頻,同時信號的時域部分能量會增大。刀具后刀面磨損后,切向力及徑向力增加,而軸向力卻未有較大的變化。

2)設(shè)計了EMD-SVM的刀具識別模式,采用EMD方法對信號進(jìn)行重構(gòu),對重構(gòu)后的信號進(jìn)行特征提取并組成特征值矩陣,將不同刀具磨損階段特征值矩陣輸入SVM中,能對刀具磨損進(jìn)行識別。

3)選用160組每個刀具磨損階段的特征值矩陣對本研究提出的刀具磨損識別模式進(jìn)行訓(xùn)練,選用40組特征值矩陣進(jìn)行驗證,對刀具磨損階段的識別成功率達(dá)到了99.17%。實驗結(jié)果表明,該模式具有良好的識別效果,可以準(zhǔn)確實現(xiàn)刀具磨損階段的識別。

4)通過對比BP神經(jīng)網(wǎng)絡(luò)模型,本研究提出的EMD-SVM模式運算時間提高60%以上,識別精度提高15.12%,對鈦合金銑削過程中刀具磨損階段具有較高的敏感度,能準(zhǔn)確識別刀具磨損的發(fā)生。

猜你喜歡
信號模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
孩子停止長個的信號
3D打印中的模型分割與打包
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国产在线91在线电影| 免费高清自慰一区二区三区| 国产一在线观看| 天天摸天天操免费播放小视频| 日韩亚洲综合在线| 中文字幕亚洲另类天堂| 欧美不卡二区| 国产在线97| 国产区免费精品视频| 青青草原偷拍视频| 91无码视频在线观看| 亚洲精品色AV无码看| 国产在线观看第二页| 中字无码av在线电影| 影音先锋丝袜制服| 国产av一码二码三码无码| Jizz国产色系免费| 国产又色又刺激高潮免费看| 日本伊人色综合网| 日本三区视频| 亚洲第一福利视频导航| 女人18一级毛片免费观看| 亚洲精品手机在线| 欧美成人亚洲综合精品欧美激情| 97一区二区在线播放| 久久semm亚洲国产| 视频在线观看一区二区| 国产激爽大片高清在线观看| 98精品全国免费观看视频| 日本精品视频| 自拍偷拍欧美日韩| 免费国产小视频在线观看| 欧美 国产 人人视频| 国产精品九九视频| 色视频久久| 极品国产一区二区三区| 在线观看国产小视频| 国产成人无码AV在线播放动漫| 精品国产成人三级在线观看| 亚洲av无码片一区二区三区| 亚洲AⅤ永久无码精品毛片| 国产精品美女免费视频大全| 8090成人午夜精品| 白丝美女办公室高潮喷水视频| 亚洲视频黄| 久久精品丝袜| 99久久无色码中文字幕| 日韩少妇激情一区二区| 欧美日韩午夜视频在线观看| 亚洲无限乱码一二三四区| 精品一区二区久久久久网站| 亚洲91在线精品| 久久综合丝袜长腿丝袜| 欧美精品不卡| 另类综合视频| 成人在线综合| 国产无码网站在线观看| 一区二区影院| 欧美日韩午夜| 久久精品无码中文字幕| 中国一级特黄大片在线观看| 亚洲Av综合日韩精品久久久| av天堂最新版在线| 国产成人综合亚洲欧美在| 久久久久久久97| 99久久亚洲精品影院| 91热爆在线| 无码精品国产dvd在线观看9久 | 园内精品自拍视频在线播放| 国产亚洲欧美在线中文bt天堂| 亚洲色图在线观看| 99视频在线观看免费| 亚洲欧洲自拍拍偷午夜色无码| 日韩 欧美 小说 综合网 另类| 亚洲精选无码久久久| 尤物视频一区| 久久综合激情网| 欧美亚洲一区二区三区导航| 久久综合九色综合97网| 亚洲成人一区二区三区| 国内精品视频| 亚洲欧美日韩精品专区|