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

基于SMO-SVM的單點金剛筆鈍化監測

2015-10-29 02:14:03李郝林遲玉倫
中國機械工程 2015年20期
關鍵詞:分類信號模型

岳 泰 李郝林 遲玉倫

上海理工大學,上海,200093

基于SMO-SVM的單點金剛筆鈍化監測

岳泰李郝林遲玉倫

上海理工大學,上海,200093

針對單點金剛筆在砂輪修整過程中易于鈍化且難以檢測的問題,使用支持向量機建立智能模型。為了得到建立模型所需的樣本庫,使用小波包分析等方法在線提取修整時聲發射信號中的特征信息,并引入鈍化平臺直徑定義鈍化臨界值。模型本身選用基于串行優化算法的支持向量分類機,使用交叉驗證法搭配遺傳算法以達到優化模型參數的目的。實驗結果表明,該模型在分類精度和計算時間上均優于一般的智能模型,可以有效地監測金剛筆的鈍化。

單點金剛筆;支持向量分類機;聲發射信號;串行優化算法;鈍化平臺直徑

0 引言

一直以來,磨削都是精密加工的主要實現方法。砂輪作為磨削加工的主角,其表面形貌直接決定著被加工工件的質量[1]。單點金剛筆是應用最為廣泛的砂輪修整工具之一,它具有原理簡單、易于操作和成本低廉的優點,但也常常因為自身的鈍化而影響砂輪修整的質量。

判斷金剛筆是否鈍化的傳統方法是離線的直接測量法,這種方法需要額外的人工操作以及停機拆裝,過程繁瑣且可靠性差[2];在線的光學監測法[3]可靠性好且原理簡單,但卻因為其光學系統的本身限制而不適合用于砂輪修整的惡劣環境中;利用聲發射信號反映金剛筆的鈍化信息,依靠類似于專家系統的模型實現對金剛筆鈍化的智能判斷,具有成本低、效率高的優點,受到國內外學者的廣泛關注[2,4]。對比相關文獻使用的監測方法發現,一般的統計過程法或實驗回歸法對訓練樣本本身分類較準確,對未知樣本則誤差較大[4-5];文獻[6]中的隱馬爾科夫模型法分類準確率較高,但模型較為復雜且魯棒性差;文獻[7]應用的貝葉斯網絡模型易于構造和使用,但分類精度相對較低;文獻[2]中人工神經網絡在處理狀態監測的非線性問題時效果較好,但由于其算法本身研究的是樣本數趨于無窮大時的漸進理論,在處理樣本數有限的實際問題時,存在分類精度差、易陷于局部最優等諸多問題。因此,分類方法的選擇歸結為尋找一種模型參數簡單、分類準確率高、運算簡潔且適合小樣本的模式識別模型。

相比其他的智能模型,支持向量機[8](support vector machines,SVM)有著樣本需求量少、模型參數簡單、全局尋優求解等優勢,但也存在著樣本較多、訓練時間長且模型參數難以選擇等缺點。針對此問題,本文選擇串行優化算法(sequential minimal optimization,SMO)來簡化支持向量機訓練中的二次規劃問題,并將其應用于金剛筆的鈍化監測。在模型輸入樣本的提取中,通過能夠實時反映細微變化的聲發射信號表征金剛筆狀態,并使用濾波分段及小波包分析等方法提取特征信息。對于模型的輸出部分,通過掃描電子顯微鏡觀察金剛筆的形貌特征,并引入鈍化平臺等相關理論對金剛筆的鈍化標準進行量化的評價。最后,將訓練好的模型分別與標準支持向量機、最小二乘支持向量機以及神經網絡進行對比。

1 單點金剛筆鈍化監測模型的總體結構和流程

1.1硬件系統

硬件系統由砂輪修整、聲發射信號采集以及金剛筆測量三部分組成。砂輪修整基于STUDER K-C33精密外圓磨床,使用L1102型金剛筆修整粒度為60目的平形鉻剛玉砂輪(400 mm×50 mm×203 mm);在AE信號采集方面,選擇AE傳感器以及高速數據采集卡,配以計算機和輔助器件;金剛筆的測量則通過3D輪廓儀以及掃描電子顯微鏡實現。

圖1 砂輪修整實驗臺

圖1所示為金剛筆修整砂輪的實驗現場,在整個過程中,選用相同的修整工藝以保證結果的準確性(表1)。考慮到金剛筆的修整深度對后文中的鈍化判別和特征選擇有較大影響,選取多個修整深度以增加實驗的泛化性能。AE傳感器依靠磁力固定在工作臺上,距離修整區域100 mm。采集卡的采樣頻率設置為1 MHz,從而避免高頻電氣噪聲,并盡可能多地采集有用的信號頻段。此外,傳感器和金剛筆的安裝位置都被嚴格固定,以保證信號的穩定和實驗的可靠性。

表1 修整工藝參數

1.2軟件系統

如圖2所示,軟件系統有模型建立和模型使用兩個步驟,模型的使用需要建立在模型建立的基礎上,且每個步驟中各個模塊的作用方式都不盡相同。整個軟件系統包括輸入模塊、輸出模塊以及支持向量機模塊。

圖2 金剛筆鈍化監測系統結構和流程

2 模型輸入:聲發射信號的處理和提取

2.1信號預處理

(a)噪聲信號

(b)修整信號圖3 環境噪聲和修整信號的頻譜對比

在實驗中,為了充分采集金剛筆狀態信息,設置較高的采樣頻率,并且不可避免地采集了修整周期中的非修整部分。因此,對AE信號進行適當的預處理十分必要。通過觀察發現,在修整周期的駛入過程中,AE信號只受環境噪聲影響。因此,基于此時的AE信號對環境噪聲以及修整信號進行頻譜分析(圖3)。可見,環境噪聲集中在100 kHz以下,而修整信號則大多分布在100 kHz以上。根據此結論,設計了IIR高通數字濾波器,設置截止頻率為100 kHz。由于AE信號有很強的突發性和衰減性,故需要濾波器有較強的衰減能力。因此,本文使用高階的巴特沃斯型IIR濾波器來實現濾波[9]。

觀察圖4發現,經過濾波后的AE信號數據過于龐大,而且駛入、駛出、修整開始和修整結束這四個階段的AE信號具有很強的不確定性,明顯不能用于之后的特征提取。出于避免無用信號及減少計算量的目的,編寫數控系統的內部PLC程序,使系統在周期開始時給出一個高電平信號,在周期結束時再恢復至低電平,通過此方波信號的上升沿和下降沿得到修整周期的時間信息。在修整周期確定的基礎上,提取出修整周期中第10~15 s的穩定修整的信號并將其按秒分段,從而得到了數據量小且代表性強的AE信號。

(a)時域圖

(b)頻域圖圖4 修整周期濾波后的時域和頻域

2.2信號特征提取

AE信號的特征提取方法通常有時域法和頻域法兩種。本文綜合以上兩種方法,先對AE信號進行頻域分段,再將特征明顯頻段的信號進行時域處理。相比單獨的頻域法或時域法,該方法可以更深層次地挖掘AE信號中的特征信息。

通常,小波分析和小波包分析都可以完成信號頻域和時域的雙尺度分解。一般的小波分解只是將低頻的近似信號進行多次細化處理,未能提高高頻細節信號的頻率分辨率[10]。小波包算法將一對鏡像濾波器組同時作用在近似信號和細節信號上,使它們同時進行分解,是小波算法的進一步深化。就金剛筆的狀態監測來說,AE信號的低頻部分大多是環境噪聲,而高頻才包含了真正有用的信息。因此,只對低頻信號進行分解處理的小波分析顯然不適用。

基于以上討論,選用可以高頻細分的小波包分析對信號進行處理。假設對信號進行n層的小波包分解,將原始信號按頻率平均分解至2n個頻段上,則可以得到同樣數量的小波包系數(圖5)。設Si,j表示第i層分解中第j個小波包系數,fs表示原始信號的采樣頻率,則Si,j包含著原始信號中fs(j-1)/2i+1到fsj/2i+1的時域信息。一般情況下,小波包系數相當于一個獨立的時域信號。

圖5 小波包分解原理

由小波包的定義可知,不同的小波包系數對應完全相同的頻帶寬度。這使得小波包系數的能量分布可以較好地反映原始信號的頻譜特性。同時,在金剛筆的狀態監測中,AE信號的頻譜特性同金剛筆的鈍化狀態具有很強的聯系。因此,利用小波包系數能量、頻譜特性、金剛筆鈍化狀態三者之間的遞推關系,將小波包系數的能量分布用于特征提取中。設第n層分解的總能量為E,第j個小波包系數Sn,j的能量為Ej,則有

(1)

(2)

其中,z表示小波包系數的離散點幅值,m表示采樣點數,P為能量比,利用參數j和P可以得到各個小波包系數的能量關系。

基于以上流程,對經過預處理的AE信號進行小波包分析。多次實驗結果表明,使用4層的小波包分析可以較好地識別頻率特征,又能避免頻段的過分細分。而在小波包基的選取上,針對AE信號的瞬態突發性以及在時域波形上的沖擊特性,選擇具有正交性、時域緊支性、頻域快速衰減性,且時域波形和AE信號類似的多貝西一類小波(db1)。在參數確定的基礎上,分別計算信號在金剛筆鈍化和銳利兩種狀態下的能量分布,再通過對比絕對值的方法獲取鈍化和銳利之間的能量變化。

通過分析不同修整深度下的AE信號發現,隨修整深度的減小,各個頻段的能量變化值相應變小,但3、6、7、14所代表的頻段一直保持變化最大的趨勢(圖6)。由此可知,上述四個頻段即為所需的特征頻段。根據小波包系數可以對應到時域的特點,分別計算它們的峭度作為特征樣本。為了方便支持向量機的使用,對特征樣本進行必要的歸一化處理和主成分分析處理,設定歸一化范圍為-1~1,主成分分析閾值為95%。經過處理后,樣本的維數從4降為3。

圖6 金剛筆在不同修整深度下的能量差異對比圖

3 模型輸出:金剛筆鈍化標準的確立和測量

3.1鈍化標準的確立

在評價金剛筆是否鈍化前,需要確定合理的鈍化評判標準。本文根據分析金剛筆鈍化對修整的影響來確定金剛筆的鈍化標準,以此評價金剛筆的狀態。

單點砂輪修整的直接評價指標主要有磨粒出刃高度和磨粒分布間隔兩種。圖7為銳利金剛筆和鈍化金剛筆修整砂輪時的對比示意圖。其中,hg是當前被修整砂輪的磨粒出刃高度,代表磨粒的縱向尺寸,其值越大,磨粒所能承受的磨削時間越長。fd為磨粒分布間隔,即砂輪轉過一周時金剛筆的進給量,其值越大代表磨粒間隔越小,所對應的被磨削工件的表面粗糙度越低。磨粒分布間隔在數值上等于金剛筆進給速度vd和砂輪轉速nw的比值。

圖7 不同狀態下的金剛筆修整過程對比示意圖

文獻[11]提出,在不考慮金剛筆鈍化的情況下,可將其刃端截面等效為二次拋物線輪廓。本文在此基礎上考慮了鈍化的影響,通過使用掃描電子顯微鏡對金剛筆進行輪廓評估(圖8),發現金剛筆的鈍化程度主要體現為刃端的圓形鈍化平臺大小。用dw表示鈍化平臺的等效直徑,結合圖7不難看出,在工藝參數不變的情況下,若dw的值和磨粒分布間隔fd相當,由于鈍化平臺的作用,砂輪在理論上將不會有任何磨粒產生,此時,金剛筆即達到鈍化臨界狀態。使用表1中的修整工藝參數求得fd=0.0698 mm,所以當dw大于此值時,表示金剛筆已鈍化。

圖8 掃描電鏡下的鈍化金剛筆刃端

3.2實際鈍化信息的測量

在砂輪修整實驗中,即使是對金剛筆位置或角度的細微調整,都會對其鈍化狀態產生極大的影響。為了保證實驗的可靠性和連續性,選擇使用復印法配合3D輪廓儀來對金剛筆的鈍化狀態進行在位地測量。設置測量間隔為10個修整周期,在測量時,停止砂輪和金剛筆的運動,在不拆下金剛筆的情況下對其進行復印并評估。當測量到金剛筆鈍化平臺的等效直徑約為0.0698 mm時,繼續修整十個周期,再通過掃描電鏡觀察其細微狀態,以確認其鈍化。值得注意的是,由于每次測量后都需要重新對刀,而對刀后的前幾個周期可能會出現修整深度過淺或過深的問題,所以需要確定信號穩定后,再進行AE信號的采集。

結合前文中AE信號的特征向量,設達到鈍化標準的金剛筆所對應的樣本為-1類,未鈍化的為1類。至此,模型的樣本庫已建立。

4 支持向量機模型的建立和驗證

4.1支持向量機的選擇

本文的金剛筆狀態監測是一個典型的非線性二值分類問題,考慮到模型的實用性要求,選擇使用C型支持向量分類機進行分類。

假設給定兩類訓練輸入樣本xi和輸出樣本yi,代入C型支持向量分類機中可得到模型函數:

(3)

使用函數φ將輸入向量映射至高維的特征空間以滿足非線性分類的需要。其中,w為向量機決策超平面的法向量(對應為神經網絡的權值),b為超平面的截距(對應為神經網絡的截距),C表示對于錯誤分類的懲罰系數,ξi代表用于提高模型容錯率的松弛變量。在模型函數的基礎上,引入核函數的定義:

K(xi,xj)≡kT(xi)φ(xj)

(4)

顧名思義,核函數是支持向量機的核心,它賦予模型從低維向高維映射的能力以解決非線性問題。由于在高維空間中,樣本之間只進行內積運算,所以只要滿足Mercer條件,核函數就可以在原空間中映射高維空間的任何內積。在此過程中,可調參數并沒有增加,從而避免因升維導致的維數災難。常用的核函數有線性核函數、徑向基(radial basis function,RBF)核函數、多項式核函數、多層感知核函數(sigmoid)四種。根據文獻[12]的論述,線性核函數是RBF核函數的特殊情況;多項式核函數在樣本維數較高時,會導致數值的畸形(趨于零或無窮);而Sigmoid核函數雖然性能和RBF核函數相當,但只有在其核矩陣對稱半正定時才能發揮作用。因此,本文選擇使用結構簡單、性能優越的RBF函數作為模型的核函數。

在此基礎上,引入拉格朗日函數對模型進行處理,得到了模型的對偶問題:

(5)

4.2串行優化算法的使用

串行優化算法是分解算法的特殊形式,即僅選擇兩個數據點作為工作集,考察其中不滿足KKT條件的點,并啟發式地將數據換入換出,反復迭代直到所有的數據都滿足KKT條件[13]。此算法的優勢在于可以通過兩個數據點得到解析解,從而大大節省計算時間和儲存空間。在訓練SVM模型時,SMO算法每一次只選取兩個拉格朗日乘子進行優化,再更新模型。首先確定根據式(5)得到當前的分類判別函數:

(6)

其中,b可以通過KKT條件中的互補松弛性求出。設第i個點的誤差為Ei,則有

Ei=f(xi)-yi

(7)

在此基礎上,可以求解兩個乘子的二次規劃。分別設α1和α2對應兩個乘子 ,基于以上超線性約束和超區域約束,可以得到兩個乘子在類別相同時和不同時的關系,進一步確定相關目標函數并將偏導數置零,可得到兩個乘子的解:

(8)

η=K(x1,x1)+K(x2,x2)-2K(x1,x2)

此外,需要簡化KKT條件并得到SMO算法的停機準則:

(9)

當pi為零時,即表示所有樣本滿足KKT條件,算法停止。拉格朗日乘子的選擇基于啟發式算法,分別針對兩個乘子設計兩個嵌套的循環以保證模型的決策函數在優化過程中有明顯的更新。初始化α=b=0,α1所對應的外部循環只考慮違反KKT條件的樣本,首先在非邊界樣本中選擇,若沒有找到,則在所有樣本中串行選擇;若仍沒有找到,則退出程序。在找到α1后,進入算法的內部循環,從違反KKT條件的非邊界樣本中找出使|E1-E2|最大的α2,重復進行兩個循環直到滿足停機準則。

4.3結果驗證和對比

前文已確定了模型的類型結構以及所需的樣本庫,現只需訓練模型參數即可完成SVM模型的構建。為了提高模型泛化能力,不失一般性地從樣本庫中隨機選取100組樣本作為訓練集,100組樣本作為測試集。在訓練時,針對懲罰系數和徑向基寬度系數難以選擇的問題,應用具有較好收斂性的遺傳算法對模型參數進行全局優化,同時配合交叉驗證法以充分利用訓練集。設置種群數量為20,進化代數為100,指定分類精度為適應度函數,將訓練集分為5組進行交叉驗證。

圖9和圖10分別為訓練樣本在三維空間中的分布圖以及測試樣本經過SVM分類后的混淆矩陣。不難看出,訓練樣本在分布上具有明顯的非線性,單單從三維空間中觀察很難找到規律。

圖9 支持向量機訓練樣本分布

圖10 分類結果的混淆矩陣

圖11 SMO-SVM模型性能

圖12 LS-SVM模型性能

圖13 O-SVM模型性能

圖14 BPANN模型性能

為了驗證模型的分類性能,使用修整深度為5 μm的樣本對串行優化支持向量機(SMO-SVM)、最小二乘支持向量機(LS-SVM)、原始支持向量機(O-SVM)以及BP神經網絡(BPANN)進行訓練和驗證,并都使用遺傳算法進行交叉尋優,得到了圖11~圖14所示的性能比較結果。可見,在訓練時間上,SMO-SVM模型和LS-SVM模型相當,BPANN模型和O-SVM模型有較大劣勢,且BPANN模型相對于其他模型有較強的不穩定性。而在分類精度上,SMO-SVM模型和O-SVM模型相當,明顯優于LS-SVM模型和BPANN模型。綜合評價后得到,本文所使用的SMO-SVM模型在整體性能上具有明顯的優越性。

分析對比結果發現,模型性能的差異主要是因為算法本身的不同。BPANN模型在訓練過程中參變量過多,每個神經元的各個參數都需要依據梯度最大方向更新,從而使得訓練時間變長。神經網絡的分類精度較低的問題,則歸結于其經驗誤差最小化原理和基于無限大樣本的特點。由于修整過程中的樣本數目較少,使得神經網絡可不能有效地對未知樣本進行判斷。LS-SVM模型、SMO-SVM模型以及O-SVM模型的差別在于它們針對訓練中的二次規劃問題使用了不同的方法。LS-SVM模型將二次規劃問題簡化為一個線性方程,失去了凸優化的最優解特性,犧牲了分類精度卻大大加快了訓練速度。O-SVM模型的分類性能較好,但是由于二次規劃的計算量較大而導致訓練時間較長。SMO-SVM模型利用串行優化算法巧妙地簡化了二次規劃問題,并采用啟發式的算子搜索,所以能在較短的時間內訓練出分類精度較高的模型。通過計算可知,對于本文的SMO-SVM模型,分類精度和金剛筆的修整深度成正比。在特征最不明顯的5 μm修整深度的樣本中,模型的平均訓練時間僅為3.3893 s,平均分類精度可達到95.2571%,這充分證明該模型可以有效地監測金剛筆的鈍化狀態。

5 結論

(1)鈍化平臺直徑可以較好地量化評價金剛筆的鈍化程度。

(2)小波包分解可以有效地從高頻聲發射信號中提取特征信息。

(3)基于串行優化算法的支持向量機在金剛筆鈍化分類上具有比原始支持向量機、最小二乘支持向量機以及神經網絡更優越的性能。

[1]Marinescu M,Hitchiner E,Uhlmann,et al.Handbook of Machining with Grinding Wheels[M].Boca Raton:CRC/Taylor & Francis Group,2007.

[2]Martins C H R,Aguiar P R,Frech A,et al.Tool Condition Monitoring of Single-point Dresser Using Acoustic Emission and Neural Networks Models[J].Instrumentation and Measurement,2014,63(3):667-679.

[3]Habrat W,Batsch A,Porzycki J.Monitoring of the Single Point Diamond Dresser Wear[J].Archives Civil Mech.,2005,5(1):13-18.

[4]Shi D,Axinte D A,Gindy N N.Online Machining Process Monitoring Using Wavelet Transform and SPC[C]//Proceedings of the IMTC.Sorrento,2006:2081-2086.

[5]Karpuschewski B,Wehmeier M,Inasaki I.Grinding Monitoring System Based on Power and Acoustic Emission Sensors[J].Annals of the CIRP,2000,49(1):235-240.

[6]Liao T W,Hua G G,Qu J,et al.Grinding Wheel Condition Monitoring with Hidden Markov Model-based Clustering Methods[J].Machining Science and Technology,2006,10(4):511-538.

[7]林峰,焦慧鋒,傅建中.基于貝葉斯網絡的平面磨削狀態智能監測技術研究[J].中國機械工程,2011,22(11):1269-1273.

Lin Feng,Jiao Huifeng,Fu Jianzhong.Research on Intelligent Monitoring Technique of Machining State for Surface Grinder Based on Bayesian Network[J].China Mechanical Engineering,2011,22(11):1269-1273.

[8]張學工.關于統計學習理論與支持向量機[J].自動化學報,2000,26(1):32-42.

Zhang Xuegong.Introduction to Statistical Learning Theory and Support Vector Machines[J].Acta Automatica Sinica,2000,26(1):32-42.

[9]Yang Z S,Yu Z H.Grinding Wheel Wear Monitoring Based on Wavelet Analysis and Support Vector Machine[J].The International Journal of Advanced Manufacturing Technology,2012,62(1/4):107-121.

[10]崔錦泰.小波分析導論[M].西安:西安交通大學出版社,1997.

[11]Doman D A,Warkentin A,Bauer R.A Survey of Recent Grinding Wheel Topography Models[J].International Journal of Machine Tools and Manufacture,2006,46(3/4):343-352.

[12]Keerthi S S,Lin C J.Asymptotic Behaviors of Support Vector Machines with Gaussian Kernel[J].Neural Computation,2003,15(7):1667-1689.

[13]Keerthi S S,Shevade S K,Bhattacharyya C,et al. Improvements to Platt’s SMO Algorithm for SVM Classifier Design[J].Neural Computation,2001,13(3):637-649.

(編輯陳勇)

Tool Wear Monitoring of Diamond Single-point Dresser Based on SMO-SVM

Yue TaiLi HaolinChi Yulun

University of Shanghai for Science and Technology,Shanghai,200093

An intelligent monitoring model was proposed based on support vector machine to solve the problem of identifying the wear of diamond single-point dresser in the dressing process of grinding wheel.To obtain the required samples for modeling,wavelet packet analysis was used to extract the feature informations from acoustic emission signals during the dressing process,and the diameter of wear platform was employed to define the threshold of dresser wear. Besides, for improving the practicability of the model, a SOM method was applied to train the support vector classifier,the parameters of the model were selected by using genetic algorithm as well as cross validation method. Experimental results show that the model has higher performance than general intelligent model, and can monitor the wear of the dresser effectively.

diamond single-point dresser;support vector classifier;acoustic emission signal;sequential minimal optimization(SMO) method;diameter of wear platform

2015-04-20

國家科技重大專項(2013ZX04008-011)

TH117.1DOI:10.3969/j.issn.1004-132X.2015.20.007

岳泰,男,1992年生。上海理工大學機械學院碩士研究生。主要研究方向為精密檢測技術。李郝林,男,1961年生。上海理工大學機械工程學院教授、博士研究生導師。遲玉倫,男,1983年生。上海理工大學機械工程學院講師。

猜你喜歡
分類信號模型
一半模型
分類算一算
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
分類討論求坐標
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
數據分析中的分類討論
教你一招:數的分類
主站蜘蛛池模板: 日韩av在线直播| www.亚洲天堂| 国语少妇高潮| 美女视频黄频a免费高清不卡| 在线中文字幕网| 国产人人干| 国产精品第一区在线观看| www.亚洲一区| 一级做a爰片久久毛片毛片| 91小视频在线| 五月婷婷欧美| 成人免费一区二区三区| 免费在线播放毛片| 亚洲欧美激情另类| 97无码免费人妻超级碰碰碰| 久久人人97超碰人人澡爱香蕉| 精品国产一区91在线| 欧美一级在线| 内射人妻无码色AV天堂| 亚洲国产综合自在线另类| 青青草原国产| 久久黄色视频影| 亚亚洲乱码一二三四区| 國產尤物AV尤物在線觀看| 欧美精品v欧洲精品| 在线国产三级| 国产熟女一级毛片| 99re在线视频观看| 二级特黄绝大片免费视频大片| 国产精品香蕉在线| 精品少妇人妻无码久久| 久久人妻系列无码一区| a级毛片毛片免费观看久潮| 99视频免费观看| 欧美天天干| 久久96热在精品国产高清| 国产美女久久久久不卡| 91精选国产大片| 国产第一页第二页| 中国一级特黄视频| 亚洲成aⅴ人在线观看| 狠狠色噜噜狠狠狠狠色综合久 | 国产99视频精品免费观看9e| 欧美区一区| 国产拍在线| 国产精品免费露脸视频| 免费jjzz在在线播放国产| 日本午夜影院| 久久久久久久久18禁秘| 午夜视频在线观看免费网站| 成年网址网站在线观看| 国产一级毛片在线| 国产肉感大码AV无码| 亚洲日本在线免费观看| 亚洲成AV人手机在线观看网站| 欧美精品xx| 第九色区aⅴ天堂久久香| 无码精品一区二区久久久| 日韩欧美91| 欧美有码在线| 国产剧情国内精品原创| 中文字幕无码电影| 99精品高清在线播放| 亚洲国产成熟视频在线多多 | 夜精品a一区二区三区| 成年女人a毛片免费视频| 亚洲综合欧美在线一区在线播放| 亚洲中文字幕国产av| 91麻豆国产在线| 国产又大又粗又猛又爽的视频| 亚洲欧美综合另类图片小说区| 91免费在线看| 国产经典免费播放视频| 色综合中文| 亚洲黄色成人| 99视频在线免费观看| 亚洲精品第一页不卡| 成年人福利视频| 亚洲欧美在线看片AI| 精品乱码久久久久久久| 国产白浆视频| 欧美激情视频一区|