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

數字孿生驅動的小樣本旋轉機械剩余壽命預測

2023-12-18 03:21:36張誠馬梓瑋劉斌孫錚徐俊
西安交通大學學報 2023年12期
關鍵詞:機械信號模型

張誠, 馬梓瑋, 劉斌, 孫錚, 徐俊

(1. 巴斯夫新材料有限公司研發部, 200137, 上海; 2. 西安交通大學未來技術學院, 710049, 西安;3. 西安交通大學機械制造系統工程國家重點實驗室, 710049, 西安)

旋轉機械廣泛應用于各種生產生活場景,并且隨著生產需求與制造水平的提高,旋轉機械的集成化、精密化程度也日益增長。旋轉機械應用場景復雜,容易在可變重載、高溫摩擦、周期性沖擊等惡劣工況下產生退化乃至失效,輕則降低設備精度,影響加工質量,重則設備損壞,生產中斷,造成巨大的經濟損失和安全威脅。因此,研究旋轉機械的剩余使用壽命預測方法有助于指導旋轉機械的使用與維護,對保證旋轉機械工作質量、提高生產效率有重大意義。

旋轉機械剩余使用壽命預測方法主要包括兩種:基于機理模型的預測方法與數據驅動的預測方法。前者利用失效機理、概率統計等方法分析旋轉機械退化過程,通過理論分析、試驗驗證等方法構建退化機理模型并進行剩余壽命預測。例如,Cerrada等[1]提取傳感器信號的時域、頻域、時頻域特征,輸入隨機森林模型進行齒輪故障診斷。李乃鵬等[2]融合多傳感器數據建立Wiener過程模型,可有效預測銑刀剩余使用壽命。此外,小波變換[3-4]、排列熵[5-6]、變分模態分解(variational mode decomposition, VMD)[7-8]等方法也被廣泛應用于剩余壽命預測。數據驅動的剩余壽命預測則是基于各類預測算法構建傳感器信號與旋轉機械剩余使用壽命之間的端到端模型,在大數據支持下無需專家經驗與領域知識即可進行旋轉機械剩余壽命預測。宋亞等[9]使用自編碼器提取信號高維特征,將提取到的特征,輸入雙向長短時記憶神經網絡進行預測得到剩余使用壽命。Kong等[10]和Yu等[11]都結合CNN與LSTM網絡構建時空融合網絡來進行剩余使用壽命預測。康濤等[12]和鄧飛躍等[13]通過向卷積神經網絡引入多注意力機制和多尺度特征級聯結構,在軸承剩余壽命預測問題上取得了比傳統卷積網絡更好的抗噪性能和預測精度。部分學者也嘗試構建機理模型與數據驅動模型結合的預測方法,如鄭小霞等[14]利用改進VMD方法構建高維故障特征向量并輸入深度置信網絡進行故障診斷。張弘斌等[15]使用連續小波變換提取信號時頻域特征構造多通道樣本,使用多通道樣本訓練深度卷積神經網絡進行剩余壽命預測。Zhao等[16]將不同頻率的小波包系數輸入深度殘差神經網絡,達到向神經網絡引入領域知識的目的,可有效提高預測準確率。

受到應用場景、使用工況、工作時間等因素的影響,旋轉機械的剩余壽命通常是動態演化的,上述方法只能預測旋轉機械在設計階段或者使用過程中某一特定時刻的靜態平均壽命,無法描述旋轉機械的動態退化過程,難以實現剩余壽命的實時在線預測。

數字孿生概念由Michael Grieves在2002年提出,其核心在于構建實體對象幾何、物理、行為、規則等多維度、多領域特征的孿生模型,反應實體對象的運作規律。隨著數字孿生概念在工業場景中應用的深入,用戶已不滿足只通過幾何模型反應實體設備運動這類表象型應用,而是希望通過實時數據交互挖掘數據潛在特征,實現孿生模型對物理實體全生命周期的精準映射[17-18]。使用數字孿生技術可以模擬機械設備的動態退化過程,為剩余壽命精準實時預測提供新方法。Ren等[19]提出一種機器學習與數字孿生結合的復雜設備全生命周期數字孿生管理系統框架。許敏俊等[20]融合鉆削毛刺機理模型與循環神經網絡,實現對鉆削過程的實時仿真監測。付洋等[21]分析航空發動機渦輪盤損傷機理,使用動態貝葉斯網絡與粒子濾波算法建立渦輪盤數字孿生模型,并基于實時采集傳感器數據對孿生模型進行動態更新,解決了渦輪盤剩余壽命的在線預測問題。Wang等[22]提出一種基于參數靈敏度分析的旋轉機械數字孿生模型修正方案并應用于轉子不平衡故障預測問題,預測誤差低于5%。

旋轉機械通常具有退化周期長、數據密度大、采集成本高的特點,難以采集到大量全生命周期退化信號去構建和訓練機理預測模型或數據驅動預測模型,因此傳統模型難以在小樣本場景下準確預測旋轉機械剩余使用壽命。張西寧等[23]基于遷移學習技術使用數據量充足的源域樣本訓練神經網絡,將神經網絡結構和參數遷移到樣本不足的目標域進行預測,在小樣本下軸承故障診斷任務中取得了較好效果,但是本質上仍然需要源域上大量數據樣本的支撐。然而,使用數字孿生技術可以構建旋轉機械退化模型以監測退化趨勢并預測剩余壽命,同時根據采集到的實時信號迭代更新退化模型,通過退化模型與旋轉機械的交互實現旋轉機械剩余壽命的動態預測,是解決小樣本旋轉機械壽命預測難題的有效方法。

基于上述分析,本文提出一種Weibull可靠性理論與卷積自編碼器(convolutional autoencoder, CAE)結合的數字孿生驅動的旋轉機械剩余壽命預測方法Weibull-CAE。將旋轉機械視作實體模型,選擇Weibull可靠度函數作為旋轉機械的退化行為模型。在訓練階段僅需要使用旋轉機械早期退化信號訓練卷積自編碼器。預測階段將從旋轉機械實體模型中獲取的實時數據經重構誤差計算與映射后得到健康因子,根據健康因子對退化行為模型進行實時迭代更新以反映旋轉機械實體模型的實時退化趨勢,通過實體模型與退化行為模型的虛實交互實現旋轉機械的剩余使用壽命預測。公開數據集上的驗證結果與實際應用表明,本文方法具有在小樣本旋轉機械剩余使用壽命預測問題中的可行性與有效性,為數字孿生技術在旋轉機械設備監測運維領域的應用提供了新思路。

1 卷積自編碼器

卷積自編碼器是一種無監督神經網絡算法,由編碼器和解碼器兩個子模塊構成:編碼器對輸入信號進行編碼獲得輸入信號的高階特征信號,解碼器接收高階特征信號并對其進行解碼得到重構信號。卷積自編碼器網絡結構如圖1所示。

圖1 一維卷積自編碼器結構

(1)

(2)

(3)

CAE的訓練目標是最小化訓練信號的重構誤差。訓練后重構誤差較小說明CAE學習到了訓練信號的特征模式,可以在不需要額外信息的前提下提取訓練信號的高階特征。

基于旋轉機械早期退化信號訓練的CAE可以學習到早期退化信號模式,對早期退化信號的重構誤差接近0。隨著工作時間的增長和退化程度的加劇,旋轉機械的退化信號模式也會逐漸發生變化。由于CAE只學習了早期退化信號模式,所以CAE對其他退化時期信號的重構誤差會隨著退化程度的增大而增大,可以利用CAE對輸入信號的重構誤差來判斷輸入信號對應的旋轉機械退化程度。

2 Weibull退化行為模型

數字孿生技術驅動旋轉機械剩余壽命預測的關鍵在于構建能夠準確描述旋轉機械退化規律的退化行為模型,并根據旋轉機械實時信號對退化行為模型進行迭代更新,使其準確實時地反映旋轉機械的當前健康狀態、剩余使用壽命以及未來退化趨勢。Weibull分布是可靠性理論中的經典分布,具有靈活性高、可解釋性強的優點,廣泛應用于各類機械零部件的可靠性建模[24],適合作為旋轉機械的退化行為模型。

典型三參數Weibull概率密度函數表達式如下

(4)

式中:t為設備已工作時間;η為尺度參數;β為形狀參數;γ為位置參數。

實際應用中,通常不考慮位置參數,即默認位置參數為0,此時三參數Weibull分布簡化為兩參數Weibull分布。對兩參數Weibull概率密度函數積分可得兩參數Weibull累積分布函數,也稱Weibull不可靠度,表達式如下

(5)

式中:F表征設備的不可靠度,F(t)=1表示設備已經失效,F(t)=0表示設備還未開始退化。考慮計算使用的方便性,實際應用中更常用到的是Weibull可靠度函數

(6)

式(6)描述了隨著工作時間t的增加,設備可靠度的下降趨勢。給定 [0,1]區間上的可靠度r0,代入式(6),即可求得設備從初始時刻退化到可靠度等于r0時的理論工作時間T0,表達式如下

T0=η(-lnr0)1/β

(7)

若給定失效閾值可靠度rf,則可根據式(7)計算設備的理論壽命Tf,進而計算出當設備可靠度為r0時的剩余使用壽命Trul

Trul=η[(-lnrf)1/β-(-lnr0)1/β]

(8)

旋轉機械的退化行為模型應當根據實時信號進行模型更新,以保持與旋轉機械實體的狀態同步,因此需要基于實時數據對Weibull可靠度函數進行參數更新。Weibull可靠度函數是非線性函數,通常使用梯度下降法進行參數更新。

考慮t={t1,t2,…,tN}時刻根據實時數據計算得到的可靠度r={r1,r2,…,rN},選擇平方差作為損失函數。此時,Weibull可靠度函數的擬合誤差為

(9)

擬合誤差對形狀參數與尺度參數的偏導數分別為

(10)

(11)

給定學習率a,對形狀參數與尺度參數進行更新

(12)

(13)

3 Weibull-CAE剩余壽命預測方法

3.1 健康因子映射

旋轉機械有多種退化失效模式,不同失效情況對應的CAE重構誤差范圍不同,為統一不同失效情況,本文提出映射函數將 [0,+∞)區間上的重構誤差映射為 [0,1]區間上的健康因子(HI)。HI用于定量描述旋轉機械的健康狀態,HI為1表示設備還沒有開始退化,HI為0表示設備已經完全退化失效。根據重構誤差、健康因子的定義與旋轉機械的退化規律可知,映射函數應當滿足以下要求:

(1)定義域為 [0,+∞),值域為 [0,1];

(2)在定義域上單調遞減且平滑可導,方便計算;

(3)具有飽和性,在自變量趨近于0時函數值趨近于1,自變量趨近于+∞時函數值趨近于0。

神經網絡中常用的Sigmoid函數具有單調性與飽和性

(14)

但是,Sigmoid函數是單調遞增函數,且定義域為(-∞,+∞),值域為(0,1)。為了使Sigmoid函數滿足上述要求,改進Sigmoid函數得到指數映射函數

(15)

式中:h表示重構誤差與HI的指數映射函數;k為形狀系數;b為偏置常數偏置常數b保證重構誤差較小時HI趨近于1,形狀系數k用于調節映射函數的遞減速率。

健康因子與可靠度定義相同,均為表征旋轉機械健康狀態的指標,可以互換使用。此外由式(6)可知,可靠度與健康因子的分布區間一致,都為[0,1]。因此,使用式(15)中計算得到的健康因子直接替換可靠度,代入Weibull可靠度函數進行相應的擬合構建和參數更新。

3.2 Weibull-CAE整體結構

本文提出的Weibull-CAE剩余壽命預測模型基于CAE計算信號重構誤差,并通過映射函數將重構誤差映射為健康因子,把健康因子作為可靠度代入Weibull退化行為模型進行參數擬合與更新,進而利用式(8)預測剩余使用壽命。數字孿生驅動的Weibull-CAE剩余壽命預測模型整體框架如圖2所示。

圖2 數字孿生驅動的Weibull-CAE預測方法架構

方法架構包括CAE與Weibull可靠度函數兩個部分。CAE的作用是計算輸入信號的重構誤差用于后續的健康因子構建,在構建、訓練CAE時的目標是使CAE較好地學習到旋轉機械早期退化信號模式。CAE模型的網絡結構參數決定CAE的模型復雜度,復雜度過低時CAE無法完整反映早期退化信號模式,復雜度過高時CAE會記住包含噪聲和冗余信息在內的早期退化信號所有特征,兩者都會導致CAE無法有效學習旋轉機械早期退化信號模式,不能準確重構早期退化信號。

Weibull可靠度函數部分接收CAE計算的實時信號重構誤差,將重構誤差經指數映射函數映射為健康因子后代入Weibull可靠度函數進行參數更新,進而預測當前剩余使用壽命。總體預測流程如圖3所示。

圖3 數字孿生驅動的Weibull-CAE方法預測流程

圖3中剩余壽命預測分為模型訓練和實時預測兩個階段。訓練階段具體步驟如下:

(1)采集旋轉機械早期退化信號,進行降噪、有效信號截取、歸一化等預處理操作;

(2)將預處理后的早期退化信號輸入CAE,計算原始信號與重構信號之間的重構誤差;

(3)反向傳播重構誤差,更新CAE網絡參數;

(4)重復步驟(2)、(3),直至CAE網絡收斂。

實時預測階段具體步驟如下:

(1)采集旋轉機械實時振動信號,進行降噪、有效信號截取、歸一化等預處理操作;

(2)將預處理后的振動信號輸入CAE,計算原始信號與重構信號之間的重構誤差;

(3)使用式(15)將重構誤差映射為健康因子;

(4)根據健康因子,使用式(9)~(13)更新Weibull可靠度函數參數;

(5)使用式(8)更新當前旋轉機械剩余可使用壽命。

4 軸承試驗臺數據驗證

4.1 數據集

驗證數據選擇法國勃艮第-弗朗什孔泰大學FEMTO研究所基于PRO-NOSTIA軸承加速退化試驗臺制作的PHM2012軸承數據集[25]。試驗臺結構如圖4所示。研究人員在PRO-NOSTIA試驗臺上通過加速壽命試驗使滾動軸承在短時間內快速退化失效,在軸承水平和垂直兩個方向上布置加速度傳感器采集軸承振動信號,采樣頻率為25 600 Hz,每隔10 s采集一次信號,每次采集0.1 s信號,每組采樣數據包含2 560個數據點。

圖4 PRO-NOSTIA軸承加速退化試驗臺

PHM2012數據集包括3種工況共17組滾動軸承的全生命周期退化信號,每種工況都有2組軸承數據作為訓練集,其余均為測試集。PHM2012數據集詳細信息如表1所示。

表1 PHM2012數據集工況信息

由于試驗中僅給軸承施加水平徑向載荷,因此使用水平方向振動信號進行后續分析預測。3種工況下軸承水平方向振動時域信號如圖5所示。

(a)工況1

4.2 訓練參數設置

原始信號中每組數據包含2 560個樣本點,為減小CAE網絡復雜度與計算量,考慮到每組數據都包含至少兩個旋轉周期的信號,因此只選取前1 280個樣本點作為輸入信號。將每個訓練集軸承樣本全生命周期前15%時間內信號作為早期退化數據。

CAE模型中:卷積層、轉置卷積層的卷積核大小設置為5,步長設置為2;池化層與反池化層的池化大小與步長均設置為4;訓練網絡時批尺寸大小為64,訓練輪數為50,學習率設置為0.001。根據指數映射函數的定義,考慮到偏置常數的作用是使重構誤差較小時的HI近似為1,將偏置常數設為5,此時h(0)=0.993,滿足定義。CAE重構誤差為1時說明輸入信號與早期退化信號有一定差異,旋轉機械已處于退化中期或晚期,即h(1)<0.5。結合式(15)可得k>5,因此將形狀參數k設置為5。Weibull失效可靠度閾值設定為0.05,即軸承健康因子小于0.05時認為軸承失效。

4.3 模型訓練

采用均方根誤差(RMSE)與平均絕對百分比誤差(MAPE)作為評價指標來衡量模型預測結果,公式分別為

(16)

(17)

預測方法的目標是盡可能減小剩余使用壽命真實值與預測值之間的誤差,使得ERMSE、EMAPE趨近于0。

使用3.2小節中的參數訓練CAE網絡,訓練集工況1中第一個軸承樣本(以“1_1”的形式表示)的第一組與最后一組振動信號的重構結果如圖6所示。可以看出,使用早期退化數據訓練的CAE能夠以較高精度重構早期退化信號,但是無法重構退化末期信號。

(a)退化早期信號

使用訓練好的CAE網絡對訓練集中軸承樣本1_1、2_1、3_2進行重構并分別計算重構誤差,結果如圖7所示。可以看出,3個軸承樣本的重構誤差均是前期緩慢上升或保持穩定,后期劇烈上升,符合軸承實際退化規律。

圖7 訓練集軸承信號重構誤差

將圖7中的振動信號重構誤差映射為健康因子,隨后進行Weibull可靠度函數參數擬合,得到3個軸承健康狀態變化趨勢及對應的Weibull可靠度曲線,如圖8所示。可以看出,各個軸承樣本在退化早期、中期的Weibull可靠度函數擬合誤差偏大,而在退化晚期時的擬合誤差則較小。這是因為軸承工作時間越長,對應的健康因子數量越多,因此可靠度函數曲線的擬合效果越好,越能更加精確地反映軸承退化趨勢。上述結果也驗證了使用Weibull可靠度函數描述軸承退化趨勢的可行性。

圖8 訓練集軸承健康因子及Weibull可靠度函數曲線

4.4 結果對比與分析

使用同樣方法對測試集軸承樣本進行信號重構與重構誤差計算、健康因子映射、Weibull可靠度函數參數更新。以測試集軸承樣本1_3為例,其重構誤差如圖9所示。

圖9 測試集軸承1_3重構誤差

從圖9可知,軸承樣本1_3前12 500 s內的重構誤差在0.3~0.4范圍內保持穩定,從13 000 s開始重構誤差開始隨工作時間增加而快速上升,說明此時軸承1_3已經進入快速退化階段。

軸承樣本1_3測試集包含軸承全生命周期內前18 010 s信號,選擇前2 700 s信號作為CAE訓練數據,此后每隔1 000 s更新一次Weibull可靠度函數并預測當前軸承剩余壽命,結果如圖10所示。可以看出,隨著實時數據的累積與Weibull可靠度函數參數的迭代更新,Weibull-CAE預測方法的預測退化趨勢與旋轉機械實際退化趨勢的誤差逐漸減小,剩余壽命預測精度也逐漸提高,EMAPE從60%逐漸下降到20%。

圖10 測試集軸承1_3剩余壽命實時迭代預測結果

圖11為使用Weibull-CAE方法預測軸承1_3樣本剩余壽命時對應的健康因子及Weibull可靠度函數。可以看出,此時實時信號的重構誤差為0.75,對應可靠度為0.845。迭代更新得到的Weibull可靠度函數尺度參數η=21 865,形狀參數β=8.86,根據式(8)可得此時剩余壽命為6 876 s,實際剩余壽命為5 730 s,預測EMAPE為20%。

圖11 測試集軸承1_3剩余壽命預測結果

對其余各測試集軸承樣本進行相同操作預測剩余使用壽命,計算ERMSE、EMAPE,與其他常用傳統方法進行對比,結果如表2所示。可以看出,本文Weibull-CAE預測方法的ERMSE、EMAPE指標都好于其余兩種方法,EMAPE指標相比于RNN-HI、SOM-HI方法分別提升4.15%、10.62%,ERMSE指標則分別提升了30.1%、34.4%,并且本文方法的剩余壽命預測偏差方差較小,預測結果更穩定。這是因為本文Weibull-CAE預測方法引入Weibull可靠度函數描述軸承退化趨勢,并基于數字孿生技術使用實時數據對Weibull可靠度函數進行迭代更新,能夠實時準確反映軸承退化趨勢,從而使得預測結果具有更高精度。同時,本文數字孿生驅動的Weibull-CAE預測方法僅需少量早期退化數據用于模型訓練,能夠在小樣本情況下預測旋轉機械剩余使用壽命并通過與實時數據的交互進行模型迭代更新,使預測精度隨著設備工作時間的增加和數據量的累積而提高。

表2 PHM2012數據集不同方法預測結果對比

5 工程應用

本文數字孿生驅動的Weibull-CAE剩余壽命預測方法目前已在上海巴斯夫新材料有限公司某實驗室自動化測試系統數字孿生預測性維護平臺上部署。該數字孿生預測性維護平臺旨在融合自動化測試系統的幾何、物理、行為、規則多維模型構建測試系統孿生模型,方便實驗人員與運維人員通過孿生模型進行測試過程的遠程監測與運維。Weibull-CAE剩余壽命預測方法以圖12所示自動化測試平臺中某試驗機構聯軸器為應用對象,構建聯軸器退化行為模型以預測其健康狀態與剩余使用壽命,并將預測結果實時展示在數字孿生預測性維護平臺中。

圖12 自動化測試系統某試驗機構三維模型

聯軸器工作時的力矩信號可以反映聯軸器的健康狀態,健康狀態下的力矩信號如圖13所示。

圖13 自動化測試系統聯軸器對應力矩信號

由于該自動化測試系統搭建完成后投入使用時間較短且聯軸器的工作壽命較長,因此缺少聯軸器的全生命周期退化數據,傳統剩余壽命預測方法無法用于預測聯軸器剩余使用壽命,使用本文方法則可以解決這一問題。根據歷史情況及專家經驗給定一個合適的聯軸器初始Weibull可靠度函數,既可以在早期退化數據不足時作為對退化趨勢的事先預估為現場人員提供參考,也可以使Weibull可靠度函數更快更準確地逼近旋轉機械真實退化趨勢。通過截取聯軸器工作時的力矩信號,使用本文數字孿生驅動的Weibull-CAE方法對給定的初始Weibull可靠度函數進行基于實時數據的參數迭代更新,實時預測聯軸器當前的剩余使用壽命以及健康狀態變化趨勢,最終預測結果如圖14所示。

圖14 聯軸器剩余壽命預測結果

6 結 論

針對旋轉機械剩余使用壽命中存在的退化數據少、失效模式多等問題,本文提出一種數字孿生驅動的卷積自編碼器與Weibull分布結合的Weibull-CAE預測方法來建立旋轉機械退化行為模型,本文主要結論如下。

(1)本文所提方法根據旋轉機械的早期退化信號,在無全生命周期訓練數據集的情況下使用卷積自編碼器評估旋轉機械當前信號與早期信號的差異,基于該差異構建健康因子并擬合Weibull可靠度函數以描述旋轉機械動態退化過程并實時預測當前設備剩余使用壽命。試驗結果表明,該方法僅需要早期退化數據即可預測剩余壽命,并且隨著實時數據的增多,預測精度也會逐漸提高,具有較高的實用性。

(2)在PHM2012軸承公開數據集上的剩余壽命預測對比試驗中,對于無全生命周期數據測試集,本文方法的EMAPE、ERMSE兩個指標分別為42.62%、938.8 s,優于其他現有預測方法,具有更好的預測精度。考慮到實際應用中通常會進行提前設備維護以避免對生產運行造成影響,本文方法能夠反映設備大體退化趨勢,其預測結果也能夠在一定程度上作為維護更換時間的輔助決策依據。

(3)在缺乏全生命周期退化數據的情況下對某自動化測試系統中的聯軸器部件進行退化趨勢與剩余使用壽命預測,驗證了本文方法在小樣本情況下對旋轉機械剩余壽命預測的可行性與有效性。

(4)下一步工作是在已有方法基礎上研究不同失效模式下CAE重構誤差的變化規律,以期建立更加準確的重構誤差與健康因子映射關系。

猜你喜歡
機械信號模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
調試機械臂
當代工人(2020年8期)2020-05-25 09:07:38
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
簡單機械
3D打印中的模型分割與打包
機械班長
主站蜘蛛池模板: 亚洲AV无码不卡无码| 国产一级毛片yw| 精品一区二区无码av| 香蕉在线视频网站| 亚洲三级电影在线播放| 亚洲精品在线91| 国产欧美视频在线| 18禁黄无遮挡免费动漫网站| 国产亚洲精久久久久久久91| 美女潮喷出白浆在线观看视频| jijzzizz老师出水喷水喷出| 亚洲日本中文字幕天堂网| 国产原创第一页在线观看| 日韩在线影院| 亚洲成肉网| 亚洲婷婷丁香| 一级毛片无毒不卡直接观看| 国产91小视频在线观看| 欧美成人亚洲综合精品欧美激情| 91精品国产自产91精品资源| 久久久久中文字幕精品视频| 国产日韩精品一区在线不卡| 高清色本在线www| 亚洲天堂.com| 中国精品久久| 无码福利视频| 国产精鲁鲁网在线视频| 99re66精品视频在线观看| 日韩在线视频网| 亚洲人免费视频| 欧美精品影院| 久久99国产综合精品1| 狠狠色香婷婷久久亚洲精品| 人妻中文字幕无码久久一区| 综合社区亚洲熟妇p| 亚洲一区网站| 久久99热66这里只有精品一| 国产主播福利在线观看| 中文字幕乱妇无码AV在线| 操美女免费网站| 欧美亚洲日韩中文| 99视频在线精品免费观看6| 欧美成人精品欧美一级乱黄| 国产国产人成免费视频77777| 国产日本欧美亚洲精品视| 五月天综合网亚洲综合天堂网| 91福利国产成人精品导航| AV天堂资源福利在线观看| 狠狠v日韩v欧美v| 97久久免费视频| 国产美女在线免费观看| 伊人久久婷婷| 欧美午夜在线播放| 午夜精品久久久久久久99热下载| 午夜不卡福利| 欧美成人一级| 久久伊人操| 亚洲视频色图| 久青草国产高清在线视频| 小蝌蚪亚洲精品国产| 亚洲国产天堂在线观看| 亚洲视频一区| 国产成人精品免费视频大全五级| 狂欢视频在线观看不卡| 日韩精品久久久久久久电影蜜臀 | 国产成人一区二区| 国内熟女少妇一线天| 亚洲天堂精品视频| 久久免费观看视频| 国产精品永久在线| 久久综合国产乱子免费| 久久综合色播五月男人的天堂| 十八禁美女裸体网站| 欧美色视频日本| 制服丝袜一区| 欧美国产综合色视频| 在线免费看黄的网站| 黄色网站不卡无码| 国产高清毛片| 激情午夜婷婷| 日本免费福利视频| 强奷白丝美女在线观看|