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

基于回歸樹和AdaBoost方法的刀具磨損評估①

2018-01-08 03:12:04陶耀東曾廣圣
計算機系統應用 2017年12期
關鍵詞:深度特征方法

陶耀東,曾廣圣,李 寧

1(中國科學院 沈陽計算技術研究所,沈陽 110168)

2(中國科學院大學,北京 100049)

基于回歸樹和AdaBoost方法的刀具磨損評估①

陶耀東1,2,曾廣圣1,李 寧2

1(中國科學院 沈陽計算技術研究所,沈陽 110168)

2(中國科學院大學,北京 100049)

本文利用高速數控銑刀銑削中不同側面方向的切削力和振動信號以及聲發射信號均方根值,以數據驅動的形式對刀具磨損進行了擬合評估. 在本次研究中,分別從時域、頻域和時頻聯合域上探索與刀具磨損相關的敏感特征,具體特征提取方法包括時域統計分析、頻域上的快速傅里葉變換(FFT)和時頻聯合分析的小波變換(WT). 本文中,決策樹被用于回歸問題而非分類問題,用于評估刀具磨損值. 同時,引入AdaBoost算法對回歸樹模型進行提升,并從模型的準確性、穩定性和適用性三個方面上綜合對比了提升的決策樹回歸模型和原模型的性能. 研究表明,AdaBoost算法提升的回歸決策樹模型在預測的準確性和穩定性上都有一定程度上提高,并且在面向全新刀具磨損預測的適用性上也取得了不錯的提升效果.

刀具磨損; PHM; 統計分析; FFT; 小波變換; 回歸樹; AdaBoost

在高速的切削工件的過程中,刀具的磨損對于加工的工件表面光滑程度和尺寸精度的影響至關重要.因此,刀具磨損的狀態評估成為了一個重要的研究課題,但更為重要的是如何精確和穩定地擬合或預測一個全新刀具的磨損值. 因為,過高地預測磨損值,將可能導致刀具材料的浪費,而過低地預測磨損值將會增加工件產品的次品率,甚至出現廢品或發生機械安全事件. 同時,準確但不穩定地預測也將沒有足夠工程應用上的可信度.

目前,在對刀具磨損評估方面的研究包括: S.Huang、X. Li等人提出了一種基于支持向量回歸(Support vector regression,SVR)的刀具磨損評估模型,并通過實驗驗證了該方法的靈活性和有效性[1]; T.Benkedouh等人用SVR模型預測了刀具的剩余使用壽命 (Remainning useful life,RUL)[2]; 作者韓玉輝將粒子群優化算法引入SVR模型,并通過仿真實驗表明,該方法具有較強的推廣能力和較高的預測精度[3]; X.Li、B.S. Lim等人對于干銑過程中刀具磨損建立了模糊神經網絡 (Fuzzy neural network,FNN)模型,并闡述了FNN相對于其他傳統神經網絡模型在預測準確率和學習速度上的優越性[4]; D.A. Tobon-Mejia 等人論述了使用動態貝葉斯網絡(Dynamic Bayesian networks,DBN)來對刀具磨損狀態進行評估,并預測其RUL[5];Huimin Chen提出了一種基于多模型聯合的刀具磨損預測的方法,具有很強的針對性和擬合能力[6]. 從上述有關刀具磨損的故障診斷和健康管理(Prognostic and health management,PHM)方面的研究中可知,越來越趨向于: (1)基于模型和數據驅動式的智能診斷; (2)基于統計學和機器學習兩大理論方法的結合; (3)屬于預測性的維護,而非基于故障發生后的條件狀態數據的后發式維護. 這類基于數據驅動的智能診斷和預測性維護的方法,能夠在故障發生前幾個周期即可評估出故障發生的概率,從而給出預警信號. 除此之外,該類方法還能夠預測工件的RUL,以達到合理利用原材料和提高工作效率的目的. 但是,上述模型中也存在幾個缺點: (1)大部分模型過于復雜或者可解釋性不強;(2)單模型的回歸擬合不夠穩定,突出地體現在對于類似刀具這樣磨損模式不太固定的預測問題上.

本文在此背景下,針對數控銑刀刀具的磨損,提出基于回歸決策樹 (Decision tree regression,DTR)和集成方法AdaBoost相互結合的方式,用于擬合和預測刀具的磨損. DTR模型復雜性容易控制,并且可解釋性強;AdaBoost方法是一種多模型集成的方法,“集成”使得預測的穩定性得以提高. 在本文研究中,通過對比和實驗分析,分別從擬合的準確度、預測穩定性和方法的適用性三個方面上驗證基于AdaBoost提升的回歸決策樹 (即提升回歸樹模型,Adative boosted decision tree regression,DTR-Ada)對于刀具磨損預測的有效性.

1 分析流程

本次研究分析的主要流程包括: 數據預處理、特征提取、特征選擇、建模擬合和性能對比分析. 如圖1所示.

圖1 分析流程

2 數據分析與特征工程

2.1 數據采集過程

本次研究的數據采用2010年PHM Society所提供的數據競賽的數據集[7]. 數據集共6份,其中c1、c4和 c6為訓練集,c2、c3和 c5為測試集. 每份數據集對應315次銑削,每次銑削采集的信號包括X、Y、Z三個軸向的銑削力和振動信號,以及銑削過程中聲發射信號的均方根值,共計7個通道數據. 每次銑削過程的磨損值分為flute1、flute2和flute3三個槽向的值.整個數據采集過程中,采集頻率為 50 KHz/通道,銑削主軸速度為 10400 RPM(轉/分),進給速度為 1555 mm/min,Y 軸徑向銑削深度為 0.125 mm,Z 軸軸向銑削深度為 0.2 mm. 有關具體的數據采集系統[4],如圖2 所示.

本文研究中,只采用了c1、c4和c6三個具有真實磨損值的數據集作為實驗數據,共計3組實驗數據,每組數據空間為315×7.

圖2 實驗數據采集系統

2.2 數據預處理和特征工程

2.2.1 數據預處理

該過程的作用是清洗數據,剔除無效和異常值數據點,采用的方法是截斷法、采樣法和濾波法.“截斷”法用于剔除進刀和退刀過程中形成的無效數據. 以c1數據集X軸銑削力的第200次銑削的時序信號為例,如圖3所示,圖中顯示了尾部無效數據(箭頭范圍內虛線部分).

圖3 c1_200 剔除尾部無效數據的前后對比

“截斷”方法具體步驟為: 首先,求一次銑削過程中原始數據的上四分位值Q3,作為首尾無效數據限界點;然后,從原始數據的開頭向后和末尾向前,剔除所有小于Q3數值,直至碰到第1個不小于Q3的數據,則停止截斷.

剔除異常值之后,對數據進行了降采樣處理,降采樣比例1/10. 因此,得到的降采樣后的數據頻率為5KHz,每次銑削過程數據量大約仍有2萬個.

最后,對降采樣后的數據進行了hampel濾波處理,目的是剔除中間異常值數據. hampel濾波原理中采用了滑動窗口的機制和中值濾波的原理[8]. 由于篇幅原因,不在贅述. 經過 hampel濾波前后的信號對比,如圖4所示. 圖4中顯示了c1數據集中Y軸銑削力的第001次銑削的信號剔除中間異常值前后的對比結果(圓圈內為異常數據點). 可見,hample濾波能夠有效剔除中間數據中的異常值.

圖4 c1_001 剔除中間異常數據的前后對比

2.2.2 特征工程——特征提取

由上文所述,每次銑削過程都具有一個真實磨損值,而每次銑削的數據經過數據預處理之后也有上萬條數據,與磨損值并非一一對應關系,不能直接作為特征使用,需要“歸約”數據中的信息,提取出有效的特征.

時域統計提取特征: 本文中對信號進行時域上的統計特征提取過程中,參考S. Huang、陳侃等作者有關這方面的研究[1,9]. 時域的統計分析值包括絕對平均幅值、平均值、標準差、均方根、方根幅值、峰度、偏斜度、裕度因子、波形因子等.

時頻變換和頻譜分析提取特征: 本文對銑削力、振動信號采用的時頻變換分析方法為快速傅里葉變換(Fast Fourier transform,FFT). 頻域中提取的特征包括了幅值譜、功率譜(密度)和相位譜等.

時頻聯合能量分析提取特征: 對于時頻聯合的特征提取,本文采用小波變換 (Wavelet transform,WT)[10,11]對銑削力和振動信號進行特征提取. 其中對銑削力和振動信號進行小波變換分解的層數分別為2和4,采用的小波基函數為‘db3’. 提取的特征包括在分解出的近似系數和細節系數相應頻段上能量比.

特征提取的最終結果為,時域統計特征計40個,頻域特征計22個,時頻聯合特征計13個. 共計總特征為 75個. 樣本數量為 3組 (c1、c4和 c6),315個樣本/組. 因此,最終的特征空間大小為 945×75.

2.2.3 特征工程——特征選擇

經過特征提取過程,數據轉化為特征,每組數據集(c1、c4和c6)分別對應一組315×75大小的特征空間.但并非有所特征都與磨損值相關或有較強關系,擇優選取特征,既能夠減少特征空間大小,降低模型訓練/測試時間,又能夠提升模型擬合性能.

本文中采用的特征評估和選擇方法包括互信息(Mutual information,MI)評估方法和 F-test檢驗方法. MI的基本原理是一種基于k-最近鄰居距離的熵估計的非參數方法[12,13]. F-test檢驗用于特征評估原理描述如下:

假設特征向量為X,i表示特征序號,目標值為y,mean()、std()和size()分別表示求向量均值、標準差和長度,centered表示是否去均值中心化,則有:

相關系數Corr:

自由度n:

最終,F-test計算公式:

F-test檢驗可以評估特征與目標值具有的線性關系; MI可用于評估特征與目標值任何的非線性相關關系. 特征選擇的標準: MI評估值>0.5、F-test檢驗值>0.5. 如圖5所示,為特征選擇過程的解析.

圖5 最優特征選擇過程

通過最優特征選擇過程,可以得到flute1磨損值較優特征,共計 34 個,如表1 所示 (其中,A4 表示第4層近似(Approximate)系數).

表1 最優特征 (MI>0.5,F-test>0.5)

從表1可知,通過特征選擇出來的特征數量仍不少,但在本文中并未進行特征降維處理,主要出于以下兩點考慮:

第一,不同刀具的磨損模式不太固定,即使在同一個銑削槽向上(例如,flute1)的磨損模式也存在一定差異性,這無形增加了擬合預測的難度. 在特征數量上保持一定的規模,能夠一定程度上保證預測的準確度. 第二,決策樹模型在結點分支過程中,本身就具有一定的特征選擇作用. 因此,保持足夠的特征數量能夠在分支特征選擇過程中給予更優的決策.

3 基于回歸樹與提升回歸樹模型擬合分析與性能對比

決策樹用于回歸問題的基本原理是: 使用最小化剩余方差的方法來決定結點的最優劃分. 劃分的目的是使得子樹誤差方差最小[14,15]. 而本文中引入的AdaBoost在回歸決策樹上應用的基本原理是[16]: 首先,回歸問題的損失函數為平方損失函數. 通過擬合當前模型的殘差的方式學習一個回歸樹,并通過前向分步算法和經驗風險最小化的方式確定下一顆決策樹的參數. 最后,通過加法模型得到最終的提升回歸樹(DTR-Ada)模型.

下面將以flute1銑削槽上的磨損值擬合為例進行分析. 如上文所述,flute1提取出來的最終特征空間為945個樣本(其中第1~315來自c1數據集,第316~630來自 c4數據集,第631~945來自 c6數據集),34個維度. 主要的研究內容包括決策樹回歸的準確度、AdaBoost應用于DTR提升的準確度和穩定性分析以及AdaBoost的適用性分析.

3.1 單數據集的獨立實驗分析

單數據集獨立實驗分析的目的在于比較DTR和DTR-Ada模型對于同組內未知數據的擬合準確度和預測穩定性.

對3個數據集c1、c4和c6采取相同的數據劃分:訓練集比例 2/3,計 210 個樣本; 測試集比例 1/3,計105個樣本. DTR模型決策樹最大深度(max_depth)選擇范圍為2~10. 通過100次重復訓練和測試過程(劃分比例固定,樣本選擇隨機),取每個最大深度的100次試驗的平均決定系數R2和平均均方誤差MSE作為評估值. 確定決策樹最大深度的原則:

(1) 平均 R2 最大,平均 MSE 最小;

(2) 若多個最大深度的R2和MSE相差不大,盡量選擇較小最大深度,以簡化模型復雜度.

如圖6和圖7所示,為對c4和c6特征數據100次重復DTR和DTR-Ada訓練/測試,在不同最大決策樹深度下的結果. 從圖中可知,決策樹最大深度分別達到 6 和 4 時,R2 停止增加,MSE 停止減少. 因此,flute1c4、flute1c6所選擇決策樹最大深度分別為6和4.

在確定了決策樹最大深度的前提下,對DTR和DTR-Ada模型測試的準確度和穩定性進行分析. 如表2~表4所示,為100次重復訓練/測試過程中得到的3個獨立測試集的平均結果. 通過比較可知,DTR-Ada得到的兩項表征準確度的評估指標均優于DTR.

如圖8和圖9所示,為100次重復訓練/測試過程中分別得到的c4和c6的R2和MSE變化情況. 從圖中可知,在相同的最大決策樹深度前提下,DTR-Ada的R2和MSE變化波動比DTR對應指標小的多,即DTR-Ada面向同組未知數據的預測擬合更加穩定.

圖6 c4不同決策樹最大深度下的平均性能

圖7 c6不同決策樹最大深度下的平均性能

表2 c1 數據集兩種模型測試結果對比

3.2 多數據集聯合的實驗分析(“K-Fold”法)

多數據集聯合(“K-Fold”交叉驗證)的實驗分析的主要目的在于驗證DTR-Ada模型對于預測不同組的未知數據中的穩定性和適用性.

表3 c4 數據集兩種模型測試結果對比

表4 c6 數據集兩種模型測試結果對比

圖8 c4獨立實驗中模型測試的穩定性對比

在本次實驗中,數據集劃分的方法為: 從每組特征數據中得到訓練集和測試集的方法是“K-Fold”方法,設置K=5,即訓練集和測試集比例為4:1. 如圖10所示,為本次多數據集聯合的特征數據劃分和實驗過程. 重復20次訓練/測試過程,但訓練集和測試集仍是隨機產生.

確定最大決策樹深度的方法如同3.1節所述,通過20次重復5-Fold交叉試驗,確定6為最合理的決策樹最大深度.

如表5所示,為20次重復5-Fold訓練/測試過程中得到的聯合測試集擬合的平均R2和平均MSE. 通過比較可知,DTR-Ada得到的兩項表征準確度的評估指標仍然都優于DTR. 如圖11所示,為本次實驗中20次5-Fold重復測試得到的平均擬合效果圖,可知DTR-Ada對于測試集的擬合效果優于DTR. 說明了,DTR-Ada對于不同組特征集的擬合預測具有一定的適用性.

圖9 c6獨立實驗中模型測試的穩定性對比

圖10 c1c4c6 多數據集聯合實驗過程 (K-Fold)

表5 聯合特征數據下測試結果對比 (“5-Fold”法)

圖11 20 次“5-Fold”重復試驗的預測結果對比

如圖12所示,為20次重復訓練/測試過程中得到的R2和MSE變化情況. 從圖中可知,在最大決策樹深度為6時,模型測試結果為: DTR-Ada模型對于測試集的擬合預測仍要比DTR模型更加穩定.

3.3 多數據集聯合的實驗分析(“留一組”法)

多數據集聯合(“留一組”交叉驗證)的實驗分析的主要目的驗證DTR-Ada方法對于未知的一整組數據集的預測準確性和適用性.

本次實驗中,從每組特征數據中得到訓練集和測試集的方法是“留一組”方法: c1、c4和c6三組特征集中,兩組作為訓練集,一組作為測試集,進行 10 次重復試驗. 由于訓練集和測試集生成過程是固定的,所以10次重復試驗將呈現周期性的結果,并且模型的穩定性也無從驗證.

圖12 多數據集聯合實驗中測試的穩定性對比

確定最大決策樹深度的方法如同3.1和3.2節節所述. 在訓練集為 630×34,測試集為 315×34 情況下,通過10次重復“留一組”試驗對DTR和DTR-Ada模型進行訓練/測試,可確定決策樹最大深度為3. 下面將在最大深度為3的前提下,對DTR和DTR-Ada模型預測的準確性和模型的適用性進行分析.

如表6和圖13所示,分析過程可以參考3.1和3.2節,可知: DTR-Ada模型對于一整組未知的磨損值的擬合預測結果要比DTR更加準確,適用性也比DTR更強.

表6 聯合特征數據下測試結果對比 (“留一組”法)

4 總結

本文研究中,利用銑削過程中采集的原始信號數據,通過統計學方法、信號處理相關知識和特征工程理論的結合,對數據和特征進行了分析,形成了數據驅動式的分析流程. 同時,引入回歸決策樹模型和集成方法AdaBoost,對刀具磨損進行了擬合、預測和模型性能的對比分析. 研究表明,單獨地運用DTR模型來擬合/預測未知磨損值的效果很不理想. 通過引入AdaBoost算法對DTR模型進行提升,在擬合準確度、預測穩定性以及模型適用性三個方面上均能夠得到有效的提升.本文研究中,雖然在準確度、穩定性和適用性三方面上取得了提升效果,但仍存在有待提升的空間.

圖13 “留一組”交叉驗證重復試驗的預測對比

特征工程方面: 1)特征提取時缺乏足夠的專家知識指導,未能預先人工地剔除一些無關特征,以致特征空間過于龐大. 2)本文只是泛性地應用FFT和WT方法去提取每個信號每個軸向上的特征,對于不同信號沒有針對性地選擇時域/頻域上的分析方法. 3)在特征選擇后的數量上,在能夠保證模型足夠性能的前提下,可以通過調整特征評估的閥值,例如提高F-test閥值至6.7165,或者通過特征降維的方法,來減少輸入模型的特征空間大小.

DTR-Ada模型方面: 1) 通過第3節表中結果可知,DTR-Ada模型帶來了準確度、穩定性和適用性提升的同時,其訓練和預測時間也增加了,主要由于AdaBoost算法本身“串行計算”的特點造成的. 這方面可以通過適當地減少決策樹數量、應用類似Bagging這樣的具有“并行計算”特點的集成方法或利用集群計算來解決,后兩者也是未來工業大數據處理發展的方向之一. 2)DTR-Ada模型只設置了固定參數(學習率0.5,DTR模型個數300等),沒有深入優化,其性能方面有待提升.

1Huang S,Li X,Gan OP. Tool wear estimation using support vector machines in ball-nose end milling. Proc. of the 2010 Annual Conference of the Prognostics and Health Management Society. Portland,Oregon,USA. 2010.

2Benkedjouh T,Medjaher K,Zerhouni N,et al. Health assessment and life prediction of cutting tools based on support vector regression. Journal of Intelligent Manufacturing,2015,26(2): 213–223. [doi: 10.1007/s10845-013-0774-6]

3韓玉輝. 基于粒子群優化SVM的刀具磨損量預測. 工具技術,2016,50(11): 109–112. [doi: 10.3969/j.issn.1000-7008.2016.11.031]

4Li X,Lim BS,Zhou JH,et al. Fuzzy neural network modelling for tool wear estimation in dry milling operation.Proc. of 2009 International Conference of Prognostics and Health Management. San Diego,USA. 2009.

5Tobon-Mejia DA,Medjaher K,Zerhouni N. CNC machine tool’s wear diagnostic and prognostic by using dynamic Bayesian networks. Mechanical Systems and Signal Processing,2012,28: 167–182. [doi: 10.1016/j.ymssp.2011.10.018]

6Chen HM. A multiple model prediction algorithm for cnc machine wear pHM. International Journal of Prognostics and Health Management,2011,2(2): 27.

7PHM Society. 2010 phm society conference data challenge.https://www.phmsociety.org/competition/phm/10. [2010].

8MathWorks. Signal processing ToolBox. http://cn.mathworks.com/help/signal/ref/hampel.html. [2016].

9陳侃. 基于多模型決策融合的刀具磨損狀態監測系統關鍵技術研究[博士學位論文]. 成都: 西南交通大學,2012.

10MathWorks. Wavelet ToolBox. http://cn.mathworks.com/help/wavelet/ref/wavedec.html. [2016].

11Mallat SG. A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans. on Pattern Analysis and Machine Intelligence,1989,11(7): 674–693.[doi: 10.1109/34.192463]

12Kraskov A,St?gbauer H,Grassberger P. Estimating mutual information. Physical Review E,2004,69(6 Pt 2): 066138.

13Ross BC. Mutual information between discrete and continuous data sets. PLoS One,2014,9(2): e87357. [doi: 10.1371/journal.pone.0087357]

14Wikipedia. Decision tree learning. https://en.wikipedia.org/wiki/Decision_tree_learning#Variance_reduction.

15Breiman L,Friedman J,Stone CJ,et al. Classification and regression trees. Boca Raton: Chapman and Hall/CRC,1984:358.

16李航. 統計學習方法. 北京: 清華大學出版社,2012: 137–151.

Tool Wear Evaluation Based on Decision Tree Regression and AdaBoost Algorithm

TAO Yao-Dong1,2,ZENG Guang-Sheng1,LI Ning2

1(Shenyang Institute of Computing Technology,Chinese Academy of Sciences,Shenyang 110168,China)
2(University of Chinese Academy of Sciences,Beijing 100049,China)

In this paper,the cutting force and vibration signals in different axial directions and the RMS of the acoustic emission signal in the milling of the high speed CNC cutters are fully utilized to evaluate the tool wear in the data-driven method. In this study,the sensitive features related to tool wear are explored from three aspects: time-domain,frequencydomain and joint time-frequency domain,and the feature extraction methods include time-domain statistical analysis,fast Fourier transform (FFT) between time-domain and frequency-domain,and wavelet transform (WT) in time-frequency domain. In this paper,the decision tree will be used for regression problems,rather than classification issues,to assess the tool wear value. And then,the AdaBoost algorithm is introduced to improve the performance of the decision tree regression (DTR),and the performance of the adaptive boosted decision tree regression (DTR-Ada) model and the original model are compared at the aspects of the accuracy,steadiness and applicability. The result shows the DTR-Ada model can improve the accuracy and stability of the fitting and prediction,and it also achieves a good effect on the applicability of the new tool wears prediction.

tool wear; prognostic and health management; statistical analysis; fast Fourier transform; wavelet transform; decision tree regression; AdaBoost

陶耀東,曾廣圣,李寧.基于回歸樹和 AdaBoost方法的刀具磨損評估.計算機系統應用,2017,26(12):212–219. http://www.c-sa.org.cn/1003-3254/6117.html

沈陽市2014年科學計劃項目(F14-056-7-00)

2017-03-20; 修改時間: 2017-04-10; 采用時間: 2017-04-13

猜你喜歡
深度特征方法
深度理解一元一次方程
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
深度觀察
深度觀察
深度觀察
抓住特征巧觀察
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 色综合a怡红院怡红院首页| 亚洲AⅤ综合在线欧美一区| 成人夜夜嗨| 日本午夜影院| 又粗又大又爽又紧免费视频| 97成人在线观看| 丁香婷婷久久| 91毛片网| 超清人妻系列无码专区| 亚洲成人一区在线| 精品無碼一區在線觀看 | 国产成人AV男人的天堂| 免费在线色| 欧美精品黑人粗大| 日韩免费毛片视频| 一级成人a毛片免费播放| 国产精品天干天干在线观看 | 成年午夜精品久久精品| 国产福利拍拍拍| 蜜桃视频一区| 欧美日韩在线国产| av在线无码浏览| 夜夜操国产| 亚洲国产91人成在线| 青青青伊人色综合久久| 狠狠亚洲五月天| 国产精品成人啪精品视频| 麻豆国产在线观看一区二区 | 中文字幕亚洲精品2页| 久久香蕉国产线看精品| 国产成人禁片在线观看| 中国一级毛片免费观看| 亚洲 欧美 日韩综合一区| 香蕉网久久| 亚洲国产综合精品中文第一| 不卡无码网| 人妻一区二区三区无码精品一区| 国产大片喷水在线在线视频| 精品欧美日韩国产日漫一区不卡| 亚洲伊人天堂| 久久精品66| 99久久99视频| 少妇人妻无码首页| 日韩精品成人在线| 2021国产在线视频| 国产精品xxx| 国产对白刺激真实精品91| 在线精品欧美日韩| 国产剧情国内精品原创| a级毛片毛片免费观看久潮| 国产女人18水真多毛片18精品| 1769国产精品免费视频| 久久一色本道亚洲| 久久男人视频| 91偷拍一区| 精品免费在线视频| 欧美国产三级| 日韩大片免费观看视频播放| a天堂视频| 永久天堂网Av| 精品国产成人a在线观看| 亚洲色中色| 欧美19综合中文字幕| 国产精品视频白浆免费视频| 在线精品亚洲一区二区古装| 国产后式a一视频| 亚洲综合色吧| 精品久久国产综合精麻豆 | 国产亚洲欧美在线专区| 国产好痛疼轻点好爽的视频| 欧美亚洲欧美区| 国产无码高清视频不卡| 国产精品自在在线午夜| 在线欧美a| аⅴ资源中文在线天堂| 91丝袜在线观看| 18禁色诱爆乳网站| 综合天天色| 亚洲欧美激情小说另类| 国产黄色爱视频| 国产91精品久久| 中文字幕波多野不卡一区|