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

GRU-BP在數字化車間關鍵部件壽命預測中的研究

2020-05-12 09:09:46郝俊虎崔寧寧韓豐羽徐崇良
小型微型計算機系統 2020年3期
關鍵詞:特征方法模型

郝俊虎,胡 毅,崔寧寧,韓豐羽,徐崇良

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

2(中國科學院大學 計算機與控制學院,北京 101408)

3(沈陽高精數控智能技術股份有限公司, 沈陽 110168)

E-mail :ekkohao@126.com

1 引 言

大數據的發展正驅動著傳統數控車間的數字化和智能化變革,其研究的本質在于探測事務無法被直接觀測的內在規律,從而對事務未來的發展方向進行預測.工業4.0可以通過應用大數據處理方法,數字化系統,以及更多面向未來的智能技術等使數字化車間變得智能化.然而,數字化車間的高復雜性、高自動化和靈活性給數控設備的可靠性和安全性帶來了新的挑戰.

預測性維護作為數字化車間制造維護的重要組成部分,在調度、維護管理和質量改進中起著至關重要的作用[1].數字化車間安全保證中最重要的環節是故障的排查和維修,但是目前故障的排查和維修都是需要在停機的狀態下完成,這不僅可能會影響生產加工的進度,還可能會增加人工干預的風險.人工排查和維修還有一個最難解決的問題是檢查周期的選擇,周期過短則會影響機器的正常運行,周期過長又會增加生產的風險.如果能夠對數字化車間關鍵設備的關鍵部位的剩余使用壽命(RUL)進行預測,那么將會極大的改善這個問題.

數字化車間的部件或系統的RUL定義為從當前時間到該部件或系統使用壽命結束時的長度[2].RUL預測可以看作是數字化車間安全運行的基礎,自提出以來受到了大量的研究學者對其預測方法進行研究.Camci[3]和Liao[4]等人將剩余使用壽命預測分為了基于數學和物理模型的、基于經驗的和基于數據的三種方法.基于數學和物理模型的的方法需要構建退化的物理故障模型,例如裂縫,磨損,腐蝕等.物理模型對于在沒有足夠數據可用的情況下解決RUL問題非常有用.在這種情況下,通過物理模型結合實際數據學習出RUL估計模型,可以很大程度上減少故障的發生.但是,這種物理故障模型非常復雜且難以構建,并且許多部件可能也并不存在物理模型,所以這種方法很難應用于復雜的系統.基于經驗的包括專家系統和模糊邏輯理論,這種方法非常依賴經驗的準確程度以及經驗知識的多少,需要大量的經驗積累和總結.基于數據的方法則不需要探索和理解繁瑣的機械工作原理,只需要采集足夠的數據用于模型訓練和預測即可.即對于具有足夠退化數據的設備,使用傳感器數據和操作條件數據來估計RUL的數據驅動方法是優選的且有效的.在數字化車間中,這些數據也相對容易收集,為數據驅動方法提供準確的RUL估計提供了堅實的支持.

多種預測剩余使用壽命的方法近年來被需多研究學者提出并驗證.陳雄姿等人[5]提出了使用貝葉斯最小二乘支持向量回歸的時間序列預測模型,通過NASA的公開數據集做對比實驗,證明了這種基于貝葉斯的改良方法在預測準確度是十分高的,但是由于模型對計算效率的要求,這種方法無法處理大數據條件下的預測.豆金昌等人[6]提出了使用ARIMA自回歸模型和粒子濾波(PF)預測剩余使用壽命的框架,將粒子濾波應用于更長時間的預測,使得該框架在長期和短期預測中都有較高的正確率,但該方法也并沒有考慮時間序列的輸入順序對預測結果的影響.

源于人工神經網絡的深度學習具有很強的學習能力[7],使得深度學習在很多方面都極大地改進了現有技術.張國輝[8]提出了一種深度置信網絡的時間序列預測方法,驗證了深度置信網絡在時間序列預測上的可行性,但是并沒考慮實際應用中模型訓練需要動態更新的問題;高育林[9]提出了一種基于深度學習的多模態剩余壽命預測方法,研究了多模態在描述準確性中的提高方法,并著重解決了傳統深度學習的精度不足的問題,成功將多模態下的深度學習網絡進行有效訓練并進行預測,但是多模態模型需要準確的對數據模態進行分組歸類,面對沒有明確分組的數據該模型無法使用.Yan等人[10]為了減少專家經驗和人類決策對預測的影響,提出了一種設備心電圖(DECG)的概念和一種基于深度去噪自動編碼器(DDA)和回歸操作的算法來預測工業設備的剩余使用壽命.

2 RUL預測系統架構設計

本文將數字化車間關鍵部位的RUL預測分為兩個階段,分別為采集階段和分析預測階段.采集階段從加工中的數控車床采集數據,如刀具和軸承可采集的信號有振動、溫度、空間位移、切削力以及主軸電流等.采集得到的數據將用于分析預測階段,該階段首先對數據進行預處理,需要處理的情況有缺失值、異常值以及信號中的噪音.然后通過特征提取和特征選擇得到與RUL預測最相關的特征.最后輸入進選定的模型進行訓練和預測.

圖1描述了用于數字化車間關鍵部件RUL預測的系統架構.如前所述,分析預測階段包含了五個步驟.第一步為對信號的預處理,一般是對信號的異常值進行檢查,但對于設備采集的信號,如振動等信號,一般還要包括降噪的過程,本文使用了小波降噪的方法進行降噪.第二步為特征提取,一般采集到的信號都是時域的,很多頻域和時頻特征不能直接拿到,也無法觀測,所以還要經過頻域和時頻域轉換來挖掘更多特征,本文使用了傅里葉變換和小波包變換來完成此項工作.第三步為特征選擇,經過第二階段得到的特征可能是幾十維甚至上百維,所以要通過相關分析或主成分分析等一些可以降維的方法來選擇與目標結果最相關的特征.第四步為模型訓練,經過前面幾步的處理,得到的數據已經滿足訓練的要求,把這些數據輸入進合適的模型來訓練,這個過程最常見的問題之一就是過擬合,改善的方法就是優化結構風險和增加訓練數據.最后一個階段為預測階段,即應用模型來完成預測,并根據預測結果做出相應的建議.

圖1 RUL預測系統架構

3 理論基礎

3.1 BP神經網絡

傳統的BP神經網絡是經典的3層模型.第一層是輸入層,它用來接收自變量,其神經元個數由具體的預測輸入維度決定.第二層為隱藏層,該層的神經元個數可任意.最后一層為輸出層,為最終的預測結果.BP網絡的訓練過程首先由輸入層計算出前向傳播的輸出,然后根據輸出層的誤差反向調整各層連接的權重.

3.2 門限循環單元

門控循環單元(GRU)是由Cho和Chung等人提出的一種深度學習網絡[11].GRU是對長短期記憶網絡(LSTM)在計算效率上的改進算法,它們的共同基礎是循環神經網絡(RNN).RNN改進了傳統的神經網絡沒有考慮輸入時序相關性的問題.RNN網絡隱狀態的值不僅僅取決于當前計算時序的輸入,還與前一個計算時序的隱狀態有關.隱狀態和輸出值的計算公式如式(1-2)所示.

Ht=σ(U*Xt+W*Ht-1)

(1)

yt=Softmax(V*Ht)

(2)

式(1-2)中σ為sigmoid激活函數,U為輸入層到隱藏層權重向量,V為隱藏層到輸出層權重向量,W為時間序列計算向量.RNN最大的改進是把前一次計算的信息帶到當前的計算過程,但是,當兩次計算的時間幀相距過遠時,需要十分苛刻的設置RNN訓練參數才能使長距離的時間信息被計算到,這種苛刻的條件實際計算時很難達到,一般將這種問題稱為長期依賴丟失問題.

Hochreiter和Schmidhuber提出的LSTM[12]旨解決RNN存在的長期依賴丟失的問題,它可以認為是一種能夠學習計算時序中的長期依賴性的特殊的RNN.RNN在時序計算中僅簡單的將上一步的輸入作為下一步的輸入,LSTM對此做出了改進.LSTM的核心是設計了一個記憶單元的用來存儲長期記憶,在時間序列的訓練中,記憶單元的值緩慢更新.每個LSTM單元接收三個輸入,一個是當前時序的輸入,另外兩個輸入是由前一時序的計算結果而來,一個結果保存了長期記憶,它在每次計算中僅作比較小的更新,另一個結果保存了短期記憶,完全由上一時序的輸入決定.LSTM網絡通過將新的信息遞增地添加到單個存儲器槽中來處理可變長度序列x=(x1,x2,...,xn),使用門來控制信應該記住的新信息,應該刪除的舊信息,以及應該暴露輸出的當前信息.LSTM的內部結構如圖2所示.

圖2 LSTM內部結構

圖2描述了在計算時序t,記憶狀態和隱藏狀態的計算過程.圖中從下往上第3層為3個控制門f、i和o,它們的值范圍都是0到1,0代表著完全刪除信息,1代表完全保留該信息.f為忘記門控,用來從上一步的記憶單元忘記不需要的信息.i為輸入門控,用來從輸入中選擇需要更新的信息,其右邊的tanh激活函數則是用來創建新的記憶候選項.o為輸出門控,用來從tanh激活的記憶單元中輸出當前需要的信息.記憶狀態和隱藏狀態的具體公式如式(3-5)所示.

(3)

(4)

ht=ot⊙tanh(Ct)

(5)

式(3-5)中,i、f和o是控制門,⊙表示哈達瑪乘積.結合內部結構和上述公式就可以理解LSTM與經典的RNN網絡的不同,LSTM保持了額外的記憶狀態更新,將記憶狀態與預測時與時序環境相互作用的隱藏狀態分開.但也可以看到,LSTM的計算是較為復雜,導致其計算速度慢,訓練耗時較長.針對這個問題,GRU設計了重置門和更新門來簡化LSTM的結構.

圖3 GRU內部結構

(6)

Zt=σ(xtWxz+ht-1Whz+bz)

(7)

式(6-7)中,W為權重參數,b為偏差參數.更新門和重置門也是通過sigmoid激活得到,所以它們的取值范圍也是從0到1,0同樣代表著完全刪除信息,1也同樣代表完全保留該信息.重置門的作用在于刪除預測無關信息.GRU在計算隱狀態之前,最重要的一部是首先計算出一個候選隱狀態,它是從兩部分得到,一部分是上一時序的隱狀態與重置門的哈達瑪乘積,另一部分是當前時序的輸入,這兩部分作為輸入經過tanh激活得到候選隱隱狀態.總結其計算公式如式(8)所示.

(8)

從式(8)中可以看出,重置門是用來控制上一計算時許隱狀態的值對當前輸入的影響,上一時序的隱狀態包含了長期記憶和短期記憶(上一時序的輸出)兩部分,經過重置門控制后可以將這兩部分中的垃圾信息(以后預測不再使用的信息)刪除.然后與當前時序的輸入Xt一起求和并經過tanh激活就得到了候選隱狀態,從名字上就可以看出,候選隱狀態的計算是為下一步得到真正的隱狀態做準備.

最終當前時序的隱狀態由兩部分求和得到,一部分是經過更新門更新的上一時許的隱狀態,另一部分是更新門更新的候選隱狀態.其計算公式如式(9)所示.

(9)

綜上,可以總結GRU單元的計算分為4步.第一步,使用重置門刪除上一時序的隱狀態中的與之后預測無關的信息,第二步,第一部分得到的結果同當前時序的輸入一起經過tanh激活得到候選隱狀態.第三步,使用更新門的帶經過遺忘后的上一時序隱狀態,第四步,將第三步的結果結合候選隱狀態最終得到當前時序的隱狀態.

4 RUL預測模型

本文的RUL預測模型輸入為時序多維特征,輸出為單維預測序列,是典型的多對一結構.為充分使用輸入提供的時間序列信息,設計了一種多對一的雙GRU神經網絡模型,其結構如圖4所示.

圖4中的神經網絡模型共包括6層,各層的輸入輸出參數中,B為Batch個數,T為時序數,I為輸入的特征數,O為輸出向量長度,U為隱藏層個數.RUL預測模型是個多對一的模型,所以輸出長度O為1.第1層為GRU層,輸出全部時間序列的結果;第2層為TimeDistributed層,是對第一層結果在每個時間序列上進行多對多的映射,用來提高單個時序內的學習能力;第3層為第二個GRU層,用來增強時間序列信息,但只輸出最終的時序結果;第4層為Dropout層,用來丟棄某些輸入,防止過擬合;第5層為全鏈接層,選用relu激活函數,其作用為保證學習率的情況下降低隱藏層個數;第六層為全鏈接層,也是最終輸出層,輸出為1維.下一節將使用此模型訓練神經網絡,然后使用訓練好的神經網絡預測刀具和軸承的RUL序列.

圖4 基于GRU-BP神經網絡的RUL預測模型

5 仿真實驗

5.1 數據來源

1)銑削刀磨損數據.本為使用的銑削刀數據來源于PHM 2010 Data Challenge.該數據集包含6個6mm球鼻碳化鎢鋼刀的工作數據記錄,分別用c1,c2,…,c6表示,取前4組作為訓練數據,后2組作為測試數據.每個刀具反復用于切割同一類工件,主軸轉速為10400RPM,進給率為1555 mm/min,橫向切深為0.125 mm,縱向切深為0.2 mm,采樣率為50 KHz進行實驗.每次切割0.001 mm后測量刀具三個凹槽的磨損.此外,使用測力計測量X,Y,Z三個軸方向的切削力.三個加速度計用于測量X,Y,Z方向上的切割過程的振動.聲發射(AE)傳感器用于測量切割過程中工件的聲學特征(AE-RMS).每組刀具數據包含315個切割文件,包含X軸切削力、Y軸切削力、Z軸切削力、X軸振動、Y軸振動、Z軸振動、聲音信號RMS以及每次切割后產生的三個凹槽的磨損數據.

2)軸承退化數據.本文使用的軸承數據來源于法國國家應用力學實驗室FEMTO-ST為IEEE PHM 2012 Prognostic challenge 提供的比賽數據[13].該數據由PRONOSTIA平臺完成實驗,共提供了17個軸承全部生命周期的溫度和加速度數據(包括徑向水平和垂直兩個方向).每組數據包含約7萬行溫度數據和約180萬行加速度數據.每個測量值都以100Hz的頻率測量.振動傳感器由兩個微型加速度計組成,彼此夾角成90°,分別位于垂直軸和水平軸上.兩個加速度計徑向放置在軸承的外部通道上.加速度的采樣頻率為25.6kHz,溫度采樣頻率為0.1Hz.17個軸承被分成三組負載進行實驗,其中一組有7個實驗軸承,負載為1650rpm轉速和4200N徑向力,分別將它們編號為b1,b2,…,b7,本文選取此組軸承的b1到b5作為訓練組,b6和b7作為測試組數據進行實驗.

5.2 小波降噪

刀具和軸承工作于數字化車間的關鍵部位,在高精度加工過程中,傳感器獲取到的信號往往含有噪聲.因此,在正式分析數據之前,降噪是一個非常必要的過程.本文對刀具3個軸的振動信號和軸承2個方向共5組振動信號進行小波降噪,結合采集數據的實際環境,振動信號的噪聲頻帶和有效頻帶未知,經過比較最終選用小波系數閾值法完成降噪.小波基函數選取了Daubechie、Coiflets和Symlets三種來對比效果,分解層數為5.閾值的選擇采用啟發式閾值選擇方法,考慮到軟閾值和硬閾值都有欠缺之處,轉而采用了一種軟硬閾值改良折衷法[14]來調整閾值.該改良方法在閾值附近接近硬閾值法,在向無窮大增長的過程中接近軟閾值法.選用信噪比(SNR)衡量去噪效果,其計算方法如公式(10)所示.

SNR=10lg(ps/pn)

(10)

式(10)中,ps為有效信號的功率,pn為純噪聲信號的功率.表1為刀具3個軸振動信號在小波Daubechie、Coiflets和Symlets各個階數時的平均信噪比.

表1 Daubechie、Coiflets和Symlets各個階數的信噪比

表1可以看出,使用db6小波降噪時的平均信噪比最高.圖5顯示了該小波在c1刀具x軸上的降噪效果.

圖5 刀具c1在x軸方向降噪前后比較

5.3 特征提取和選擇

完成信號降噪后,接下來就是分析并提取與刀具和軸承的RUL結果相關的特征.其步驟是首先提取特征,提取的維度包括時域特征、頻域特征和時頻域特征三個方面,然后再選擇與預測結果最相關的特征集合[15].

1.3 評價指標 ①急診效率指標:包括就診至首次球囊擴張(DtoB)的時間和介入治療前谷草轉氨酶(AST)、乳酸脫氫酶(LDH)、肌酸激酶(CK)、肌酸激酶同工酶(CK-MB),均采用全自動生化分析儀器檢測;②治療7 d后采用心尖四腔切面彩色多普勒超聲檢測患者的心功能E峰值、左心室舒張末期直徑(LVEDD)、E/A值及左心室射血分數(LEVF);③治療7 d后測定氨基末端腦鈉素前體(NT-proBNP)的水平;④治療7 d后采用問卷調查評估患者的滿意度,包括心理護理得分、生理護理得分、護理服務得分,滿分為100分。

時域特征本文選取了均值、方差、標準差、偏度、峰度、均方根、峰-峰值、峰值因子、脈沖因子和裕度因子等十個維度.頻域特征使用功率譜密度(PSD)來提取.獲取PSD的直接方法是求傅里葉變換在一較大的時間間隔內幅度的平均值,但更簡單的方法是直接通過信號的自相關方程進行傅里葉變換得到,本文使用的就是后一種方法.

時頻域特征是特征提取時一個非常重要的考量方面,其對象主要是非平穩信號.小波包變換是提取時頻域特征的一個非常合適的選擇,它是連續小波變換和離散小波變換的折衷方法,計算效率和分辨率精度都滿足要求.如圖6為刀具c1的5級小波包分解圖:

圖6 5級小波包分解各節點的能量變化

圖6中可以看到左上和右下兩個圖內的節點能量變化和剩余使用壽命的負相關性最強.圖中還可以看到某些子節點的能量較高,說明這個節點對應的頻段為能量的集中區域.小波包分解總共分解的到25個頻段節點,將這32個節點的節點系數的歸一化能量作為提取的時頻域特征.

通過前面的特征提取,每個軸向通道分別分別獲取了10個時域特征,7個頻域特征以及32個時頻域特征,共49個特征.刀具的振動信號有xyz三個軸向通道,軸承有徑向水平和垂直兩個軸向通道,所以刀具和軸承分別共有147和98個特征.這個特征數是非常龐大的,需要進一步選擇出與RUL最相關的特征用于預測模型.本文使用Pearson相關系數來進行特征選擇.分別計算刀具147個特征和軸承98個特征與剩余壽命相關性.刀具的每個候選特征與剩余壽命的相關性取4組用于訓練的數據的平均值作為相關系數參考值.同樣的,軸承的每個候選特征與剩余壽命的相關性取5組訓練數據的平均值作為參考值.最終選取相關系數絕對值大于0.95為輸入特征.

最終選擇出的刀具輸入特征有,x軸的小波包節點6、節點11、節點16和節點17,y軸的小波包節點7、節點10和節點11,z軸的小波包節點1和節點8,以及z軸的方差、標準差和均方根,總共12個特征.

最終選擇軸承的輸入特征有,x軸的小波包節點7、節點12、節點14和節點18,y軸的小波節點20和節點21,總共6個特征.

5.4 模型訓練與預測

經過前面降噪和特征分析得到了用于模型訓練的特征.將6組刀具的數據進行分區以便于提取特征和訓練模型,最終選擇用于訓練的數據有25200行,用于訓練的數據有12600行;同樣對軸承的數據進行分區特征分析,得到6150行訓練數據,1500行測試數據.將訓練數據輸入進模型進行訓練.模型訓練使用的損失函數為均方根誤差,優化器算法采用自適應矩估計優化法.經過網格搜索超參數法,得到刀具和軸承RUL預測模型的參數如表2所示.

表2 RUL預測模型參數

將前面訓練好的模型分別對刀具和軸承的測試組數據進行RUL預測,最終的結果如圖7所示,其上半部分為刀具c5和c6的RUL預測結果,下半部分為軸承b6和b7的RUL預測結果.四個預測曲線圖中,前20-50的樣本區域誤差相對其他區域稍大(推測可能的原因是此時部件工作處于磨合期),其余區域預測值和真實值的誤差基本在50s以內.如果將剩余10%壽命作為維修報警的閾值,刀具和軸承分別在1250、1400、580和650樣本區報警.整體上預測結果達到預期效果.

圖7 RUL預測結果

本次實驗除了使用常用的均方根誤差(RMSE)評價方法外,還使用了PHM2012比賽針對RUL設計的打分方法RUL Score,該打分方法考慮了RUL預測產生的正誤差(預測值大于目標值)和負誤差(預測值小于目標值)對結果的影響不同.Score的計算公式如式(11)所示.

(11)

為了方便衡量本文的實驗結果,將本文提出的方法同未考慮輸入時序因素的梯度增強算法XGBoost和傳統的BP網絡預測結果相比較,比較結果如表3所示.通過表3的數據對比可以得到,兩種評分方法都證明了考慮了輸入時序因素的GRU-BP神經網絡在RUL預測精度上要優于XGBoost和BP神經網絡.

表3 RUL預測結果比較

6 結 語

本文主要針對數字化車間關鍵部件的RUL預測提出了一種多對一雙GRU層神經網絡的RUL預測方法.首先設計了RUL預測的基本系統架構,闡述了GRU對于RNN和LSTM的改進以及GRU的基本原理,然后搭建了基于GRU-BP的神經網絡模型,在信號預處理時使用了改進的軟硬閾值適中法進行小波降噪,并充分的考慮了輸入時序因素進行RUL預測.最后通過實驗驗證了本文方法的可行性和正確性,并將其同XGBoost和BP神經網絡方法做對比,進一步顯示了GRU應用于RUL預測的準確性.

猜你喜歡
特征方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 精品撒尿视频一区二区三区| 美女毛片在线| 国产乱人免费视频| 欧美福利在线观看| 五月激激激综合网色播免费| 亚洲综合激情另类专区| 亚洲成在线观看| 日韩乱码免费一区二区三区| 亚洲成年人网| 国产网友愉拍精品视频| 国产一级精品毛片基地| 欧美一级黄色影院| 欧美午夜精品| 9966国产精品视频| yjizz国产在线视频网| 一级爱做片免费观看久久| 亚洲欧美自拍视频| 国产精品99久久久| 四虎影视无码永久免费观看| 国产毛片网站| 欧洲一区二区三区无码| 国产黄网永久免费| 无码中文AⅤ在线观看| 美女被操91视频| 91精品国产丝袜| 成年人免费国产视频| 欧美精品v日韩精品v国产精品| 久久精品国产免费观看频道| 国产在线啪| 国产成人禁片在线观看| 久久久久久久久久国产精品| 国产成人资源| 国产h视频免费观看| 亚洲色图欧美激情| 婷婷六月在线| 美女毛片在线| 亚洲成人一区二区| 成人毛片在线播放| 国产在线精品美女观看| 色婷婷亚洲综合五月| 一级毛片免费高清视频| 欧美午夜久久| 国产SUV精品一区二区| 老熟妇喷水一区二区三区| 尤物精品国产福利网站| 欧洲一区二区三区无码| 久久精品国产精品一区二区| a色毛片免费视频| 亚洲天堂成人在线观看| 欧美精品三级在线| 久久久噜噜噜久久中文字幕色伊伊| 成人精品区| 一本一道波多野结衣一区二区| 凹凸国产分类在线观看| 啪啪永久免费av| 午夜成人在线视频| 毛片一级在线| 日韩精品无码免费专网站| 国产啪在线| 丰满少妇αⅴ无码区| 久久这里只精品热免费99| 欧美a网站| 久久鸭综合久久国产| Jizz国产色系免费| 免费看a毛片| 97se亚洲| 亚洲一区二区日韩欧美gif| 国产成人1024精品| 色成人综合| 色亚洲激情综合精品无码视频| 少妇被粗大的猛烈进出免费视频| 免费看一级毛片波多结衣| 久久中文字幕不卡一二区| 伊人中文网| 91精品伊人久久大香线蕉| 国产视频大全| 欧美一区精品| 91九色视频网| 国产欧美专区在线观看| 亚洲高清资源| 91精品伊人久久大香线蕉| 国产微拍一区二区三区四区|