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

基于高斯模型的工業過程數據的故障預測

2018-01-22 06:02:49楊為惠陳彥萍溫福喜
太原理工大學學報 2018年1期
關鍵詞:故障方法模型

楊為惠,陳彥萍,溫福喜,高 聰

(西安郵電大學 a.計算機學院,b.大數據處理研究中心,西安 710121)

隨著信息技術和自動化技術的快速發展,現代工業系統的集成度和復雜度越來越高。各部分之間的相互影響也越來越復雜,導致系統發生故障和功能失效的概率逐漸加大,且故障一旦發生,危害影響極大,嚴重的會導致系統整個失效和癱瘓。故障診斷就是在這種背景下,在故障已經發生后,確定故障的性質及其對未發生故障部分的影響,給出補救策略。

但如果在故障發生的早期,即在其還未對系統造成任何損害的情況下及時地檢測出故障的發生,就可以提前實施可靠的維修策略排除故障,從而避免產品損壞、系統癱瘓以及災難性事故的發生。基于這一背景,故障預測與健康管理(prognosis and health management,PHM)技術應運而生。

現代工業系統已經可以采集并且積累了海量的過程運行數據,基于數據驅動(data-driven)的故障預測方法[2]逐漸獲得重視并取得快速發展。如何有效利用這些數據,依據一些智能算法來建立系統的故障預測模型,是保障復雜工程系統安全高效運行亟待解決的問題。而工業過程數據的非線性、小樣本及自回歸特性,以及故障預測中的不確定性問題等,為預測建模提升了難度[3]。

高斯過程[4]回歸是近年發展起來的一種機器學習回歸方法,它有著嚴格的統計學習理論基礎, 對處理高維數、小樣本、非線性等復雜的問題具有很好的適應性, 且泛化能力強。與神經網絡、支持向量機相比, 具有容易實現、超參數自適應獲取、非參數推斷靈活以及輸出具有概率意義等優點。

1 相關工作與背景知識

在故障預測中,預測方法是最基礎,也是最核心的部分。

數據驅動故障預測方法主要集中在機器學習、計算智能或其他統計模型[6-7]。應用較為廣泛和相對成熟的方法有自回歸(AR)模型[8],支持向量機(SVM)[9-10],相關向量機(RVM)和神經網絡(NN,neural network)[12-13]以及隨機過程模型如馬爾可夫模型、高斯過程回歸模型等。

以上方法中,AR是線性模型,只適用于線性時間序列建模;其余方法為非線性模型,適用于工業過程中的非線性數據建模。文獻[14]對NN建模技術進行了多步預測的預測性能研究,文獻[15]提出先對時間序列進行EMD分解得到若干個分量,然后用神經網絡模型分別對時間序列的EMD 分量進行多步預測,最后將該若干個分量作為輸入,用神經網絡模型進行時間序列的多步預測。文獻[16]基于改進的小波神經網絡算法對機床軸承故障進行預測。然而,神經網絡的結構設計尚沒有統一的理論依據,比如網絡的層數、神經元的個數、傳遞函數的選擇等,因此網絡構建上有一定困難。神經網絡的訓練需要大量的數據,而在工業過程中往往只有少量的故障數據樣本,樣本缺乏致使網絡性能下降。人工神經網絡收斂速度慢,訓練過程中存在過度擬合或擬合不足等問題,并且輸出無概率意義,無法確定預測結果的置信度。

相關向量機、高斯過程回歸模型這兩種方法具備不確定性的管理能力[17]。較之NN模型,高斯過程具有完備的貝葉斯理論基礎,在給出預測結果的同時還會給出一個置信區間,用作故障預測不確定性的度量。

文獻[18]基于高斯過程建立TE仿真數據的單步預測模型,并將之運用于在線的故障檢測,但是忽略了過程數據的時間序列特性。文獻[19]提出用高斯過程模型對混沌時間序列進行單步與多步預測,并提出用混合的高斯過程核函數進行建模。盡管如此,以上文獻提出的方法(包括基于NN模型的方法),都是根據數據的統計特征建立模型,沒有將數據具有的物理意義作為先驗信息加入到模型中,而不同的工業過程、同一工業過程不同的模態下產生的數據都具有不同的物理特性。文獻[20]提出了一種根據EMD分解環境監測數據,根據數據的物理特性構建核函數的方法。本文基于高斯過程回歸(GPR,gaussian process regression)建立工業過程中時間序列數據的模型。通過針對特定數據構建高斯過程核函數來更好的描述數據特性,從而提升模型性能;并且在基于單步預測基礎上,實現了不同種方法的多步預測。實驗數據為在TE平臺[5]上采集到復雜工業過程的模擬數據,通過大量的預測精度性能對比實驗表明,構建的高斯過程預測方法較之對比的預測方法具有較高的平均預測精度。

2 基于高斯過程回歸的故障預測

2.1 單步預測模型

2.1.1 高斯過程回歸模型

假設觀測目標值yi由一個未知函數f決定,并且被一個獨立同分布的高斯噪聲ωi腐蝕,即

yi=f(xi)+ωi.

(1)

式中:ωi為獨立的隨機變量,符合高斯分布,均值為0,方差為σ2(假設訓練和測試數據的噪聲具有相同的分布,即ωi∶N(0,K).

在高斯過程回歸方法中,通常賦予f一個均值為0的高斯過程先驗概率分布:

f|x,θ∶N|(0,K) .

(2)

式中:K是一個N×N維的對稱正定的協方差矩陣,Ki,j=k(xi,yj),k(·,·)是協方差函數,也稱核函數。它是一個以超參數θ為自變量的正定函數,因此,訓練輸出y的概率分布為:

y|x,θ∶N(0,K+σ2I) .

(3)

給定一個測試輸入數據x*和對應的測試輸出數據y*以及相應的先驗概率分布抽樣得到f(x*),f(x)與f(x*)的聯合概率分布也是一個均值為0的多維高斯過程:

(4)

其中,K*=[k(x*,x1),…,k(x*,xN)]T,K**=k(x*,x*).

基于式(1)中的高斯噪聲假設,訓練輸出y和測試輸出y*的聯合概率分布為:

(5)

根據條件高斯法則,可以推導出

y*|y,x,θ,σ2∶N(m(x*),v(x*)) .

(6)

其中,預測值y*的均值m(x*)和標準差v(x*)為:

(7)

(8)

由此可以看出,已知協方差函數和它的超參數,就可以對未來時刻的數據做出預測。

2.1.2 模型訓練

模型的訓練就是獲取協方差函數的參數最優值。目前GPR模型大多采用極大似然估計(maximum likelihood estimation,MLE)獲得超參數的最優取值,即通過建立訓練樣本條件概率的對數似然函數,并對超參數求偏導,再采用共軛梯度優化方法搜索出超參數的最優解,具體求解過程可見參考文獻[4].

其中對數似然函數的形式為:

(9)

對上式求偏導的結果為:

(10)

2.1.3 基于EMD的核函數構建方法

高斯過程的協方差函數為模型的訓練增加了很大靈活性,但是怎樣構建協方差函數就變成了一個難點。運用不同的核函數訓練的模型性能會有較大的差距。以下為平方指數核函數(squared exponential,SE),周期核函數(periodic,Per),線性核函數(linear,Lin)和二次有理數核函數(rational quadratic,RQ) 的表達式:

(11)

(12)

(13)

(14)

Merce原理指出,一個矩陣只要滿足半正定就可以是一個協方差矩陣。通過將兩個核函數相加和相乘的操作組合依然是一個有效的核函數。DAVID[21]在研究中總結了基本核函數所表述的數據特性和在組合后他們的特性變化規律。

在以上的核函數中,SE描述數據局部變化的特征,Per描述數據周期性變化的特征,并且是一個全局性的周期核函數,Lin描述數據的長期變化趨勢,RQ描述數據的不規則變動;核函數的相加與相乘都會有相應的特性變化,例如,Lin核加Per核具有趨勢周期特性,SE乘Per核函數具有局部周期的特性。

HUANG et al提出的經驗模態分解(empirical mode decomposition,EMD)是一種基于信號局部特征的自適應分解方法[22],可以有效地將原始時域信號分解為具有多尺度時頻特性的本征模態函數(intrinsic mode functions,IMFs).該方法在旋轉機械故障診斷領域、高層建筑和橋梁健康狀態監測得到了廣泛應用。工業過程中傳感器監測到的過程數據分解得到的不同IMFs往往具有不同的物理含義[23]。

根據以上的規律,就可以依據特定數據的物理特性,構建具有針對性的高斯過程核函數,使模型具有更加好的預測性能。

2.2 基于時間序列的多步預測

2.2.1 預測策略

為了掌握更多的未來信息,需要對高斯過程模型進行多步預測[19]。多步預測是建立在上文單步的模型基礎上的。預測策略是多步預測的重要研究內容。

在預測策略方面,現階段應用的預測策略主要有迭代策略和直接策略。多步預測策略在TAIEB的文獻[24]中有比較全面的闡述和研究。

2.2.1.1迭代策略

迭代預測策略中,首先需要建立單步預測模型,

yt=f(xt)+ω.

(15)

式中,xt={(xt,…,xt-d+1)},yt,d為輸入維數,即時延步數,ω是高斯噪聲。

在接下來的迭代過程中,將yt當做其中一個輸入變量,同時剔除時延最久的輸入變量,從而在輸入變量個數不變和模型不變的情況下實現下一時刻的yt+1預測,以此類推,直到預測得到yt+H-1.

2.2.1.2直接策略

直接多步預測策略,即在相同的輸入數據的基礎上對每一步預測建立相應的預測模型。每一步的預測輸出值如下式:

[yt+h-1,yt+h-2,…,yt]=fh(xt,…,xt-d+1)+ω.

(16)

其中,yt,d為輸入維數,ω是高斯噪聲。

式中沒有任何的預測輸出作為下一時刻的輸入,因此,多步預測誤差并不會受到前面預測值的影響而不斷的放大。然而,因為直接的多步預測之間關聯性較弱,可能導致預測輸出值偏差較大。

2.2.2 預測原理

高斯過程回歸模型的預測方差不僅能夠解釋模型參數不確定性的影響,還可以描述過程的不確定性[17]。通過預測值和方差,可確定實際值可能的范圍,以2倍方差來設置預測值的置信區間[yt+h-1-2σ,yt+h-1+σ],其中,yt+h-1表示提前h步的預測值。從而確定可能出現故障的最早或者最遲時間。

3 實驗驗證

3.1 TE模擬平臺

由于在實際情況中,對于任何大型工業過程而言,數據是無法公開得到的。為了驗證本文提出的算法在時序數據異常檢測中的有效性,選用了TE(Tennessee Eastman)仿真平臺對本文算法進行驗證。

TE化工過程是Eastman一個實際工藝流程的原型。整個過程包括5個主要操作單元,即反應器、冷凝器、循環壓縮機、氣液分離器和產品解吸塔。TE過程共有41個測量變量和12個操控變量,并設定21種故障類型。由于攪拌速度XMV(12)在反應過程中保持恒定不變,其對系統的建模沒有任何影響,因而選擇前52個變量[XMEAS(1),…,XMEAS(41),XMV(1),…,XMV(11)]T作為某一特定時刻的原始觀測信號向量。

圖1 TE工業過程原理圖Fig.1 A process ftowsheet of the TE

在仿真實驗中,設定TE過程總反應時間為48 h,采樣間隔為3 min.所有的異常都是在仿真開始8 h后引進。因此此類故障的仿真可采集到48×20=960個數據樣本,每個樣本數據包含52維特征,每一類異常數據樣本是一個960×52的樣本矩陣。

接下來就以異常一為例進行實驗。異常一的定義是A/C進料比發生了變化,主要影響的特征變量為XMSE(1)和XMSE(4).

異常一的XMSE(1)特征變量作為數據集a1,XMSE(4)特征變量作為數據集a2.每個數據集均包含了480個訓練樣本和960個測試樣本。

3.2 對比實驗

本節將提出的高斯過程模型與BP神經網絡和應用最廣泛的RBF核函數高斯過程模型進行對比,并與兩種預測策略相結合,共形成6種多步預測方法:高斯過程迭代策略模型(GPI)、高斯過程直接策略模型(GPD)、神經網絡迭代策略模型(NNI)、神經網絡直接策略模型(NND)、高斯過程RBF迭代策略模型(GPrI)和高斯過程RBF直接策略模型(GPrD).

高斯過程模型的先驗分布為自由度為4的Student's T分布。

在應用BP神經網絡模型時,根據輸入要求,設置輸入節點個數為4,設置隱藏層的神經元個數為10.BP網絡的訓練函數為trainlm函數。

3.3 評價指標

本節采用均方根誤差(RMSE)評價多步預測方法的性能。方根誤差ERMS計算表達式為:

(17)

3.4 實驗核函數構建

圖2為原始數據通過經驗模態分解得到的4個IMF分量。IMF1表達了數據有不規則變化與局部變化的物理特性,IMF2,IMF3與IMF4表示了數據有局部周期特性與周期特性。所以在核函數的構建中通過加入SE核與RQ核描述數據的局部變化與不規則變化,通過加入Per核函與SE核相乘來描述局部周期變化。

圖2 原始數據與其IMF分量Fig.2 Raw data and its IMF component

最終構建的核函數為:

K=SE×Per+SE+RQ .

(18)

3.5 實驗結果分析

首先,基于a1和a2數據集,比較GPI,GPD,GPrI,GPrD,NNI,NND 6種方法的五步預測結果。圖3—圖7為6種方法在a1數據集上五步預測的預測值與真實值的對比。可以明顯看出,對神經網絡而言,基于高斯過程模型的預測方法得到的預測值曲線與真實值曲線擬合得更好。

圖3 a1數據集上GPD與GPI的五步預測結果與真實值對比Fig.3 Five-step ahead prediction results of GPD and GPI on a1 data set compare with the real value

6種預測方法對2種數據集的5步預測性能(通過RMSE計算得到)的統計結果如表1和表2所示,由于神經網絡的方法需要交叉驗證,故神經網絡的RMSE是5次實驗求平均值得到。

由表1和表2可以看出,本文提出的預測方法在兩個數據集上的5步預測都有最小的RMSE.這說明提出的基于高斯過程的預測方法得到的預測值與其他方法相比更接近于真實的傳感器測量值。而基于高斯過程的方法相較于神經網絡都有較小的RMSE,表明高斯過程回歸方法處理具有小樣本特性的工業過程故障數據時更有優勢,而神經網絡模型的性能可能會因為訓練樣本不夠充分而受限。

固定回歸階數P=4,預測步長依次為[1,5,10,20,50],相應的預測結果的RMSE的變化情況如表3,表4所示。可以看出,隨著預測步長的增大,各種方法的預測結果的RMSE呈逐漸增大的趨勢。當步長低于20步時,一般迭代預測方法的RMSE低于直接預測方法,而當步長大于20時,往往直接預測取得更低的RMSE,這是因為迭代的過程中會累積誤差,步長越大,迭代次數越多,誤差就越大。特別地,GPD預測方法的RMSE絕大多數低于其他方法,表明在各種預測步長之下,提出的基于高斯過程的預測方法能夠對工業過程數據進行更準確的多步預測。

圖4 a1數據集上GPrD與GPrI的五步預測結果與真實值對比Fig.4 Five-step ahead prediction results of GPrD and GPrI on a1 data set compare with the real value

圖5 a1數據集上NND與NNI的5步預測結果與真實值對比Fig.5 Five-step ahead prediction results of NND and NNI on a1 data set compare with the real value

RMSEGPGPrNNI0.04100.04180.0578D0.04010.04130.0514

表2 a2數據集上5步預測的性能統計結果Table 2 The performance of 5-step ahead prediction on a2 data set

圖6 a1數據集上GPI與GPrI的10步預測結果與真實值對比Fig.6 Ten-step ahead prediction results of GPI and GPrI on a1 data set compare with the real value

4 結束語

本文基于高斯過程模型對工業過程中的歷史數據建立預測模型,通過對過程中監測變量的多步預測,實現了故障的提前預知。針對不同的工業過程,不同的工況,數據會呈現出不同的物理特性。通過分析特定數據的物理特性,構建出可以描述該特性的核函數;同時結合多種多步預測策略和建模方法,在TE模擬數據集上進行了一系列對比實驗。實驗結果表明,本文提出的方法有更好的預測精度。下一步將TE模型的更多種類故障數據上、在多模態情況下驗證該方法的性能。

圖7 a1數據集上GPD與NND的20步預測結果與真實值對比Fig.7 Twenty-step ahead prediction results of GPD and NND on a1 data set compare with the real value

表3 a1數據集上不同步長預測性能比較Table 3 The performance of multi-step ahead prediction based on different method on a1 data set

表4 a2數據集上不同步長預測性能比較Table 4 The performance of multi-step ahead prediction based on different method on a2 data set

[1] 曾聲奎,MICHAEL G PECHT,吳際.故障預測與健康管理(PHM)技術的現狀與發展[J].航空學報,2005,26(5):627-632.

ZENG S K,MIEHAEL G PECHT,WU J.Status and perspectives of prognostics and health management technologies[J].Acta Aeronautica et Astronautica Sinica,2005,26(5):627-632.

[2] 彭宇,劉大同.數據驅動故障預測和健康管理綜述[J].儀器儀表學報,2014,35(3):481-495.

PENG Y,LIU D T.Data-driven prognostics and health management:a review of recent advances[J].Chinese Journal of Scientific Instrument,2014,35(3):481-495.

[3] 劉強,秦泗釗.過程工業大數據建模研究展望[J].自動化學報,2016,42(2):161-171.

LIU Q,QIN S Joe.Perspectives on big data modeling of process industries[J].Acta Automatica Sinica,2016,42(2):161-171.

[4] RASMUSSEN C E,WILLIAMS C K I.Gaussian processes for machine learning[M].London:The MIT Press,2006.

[5] ANDREAS BATHELT,LAWRENCE RICKER N,MOHIEDDINE JELALI.Revision of the tennessee eastman process model[J].International Symposium on Advanced Control of Chemical Processes,2014(9):

[6] HENG A,ZHANG S,TAN ACC,et al.Rotating machineryprognostics:State of the art,challenges and opportunities[J].Mechanical Systems and Signal Processing,2009,23(3):724-739.

[7] PENG Y,DONG M,ZUO M J.Current status of machine prognostics in condition-based maintenance:a review[J].The International Journal of Advanced Manufacturing Technology,2010,50(1/4):297-313.

[8] LIU D,LUO Y,PENG Y,et al.Lithium-ion battery remaining useful life estimation based on nonlinear AR model combined with degradation feature[C]∥Annual Conference of the Prognostics and Health Management Society 2012,Minneapolis,Minnesota,USA.2012:1-7.

[9] QU J,ZUO M J.An LSSVR-based algorithm for online system condition prognostics[J].Expert Systems with Applications,2012,39(5):6089-6102.

[10] WIDODO A,SHIM M C,CAESARENDRA W,et al.Intelligent prognostics for battery health monitoring based on sample entropy[J].Expert Systems with Applications,2011,38(9):11763-11769.

[11] 朱永利,尹金良.組合核相關向量機在電力變壓器故障診斷中的應用研究[J].中國電機工程學報,2013,33(22):68-74.

ZHU Y L,YIN J L.Study on application of multi-kernel learning relevance vector machines in fault diagnosis of power transformers[J].Proceedings of the CSEE,2013,33(22):68-74.

[12] LIU J,SAXENA A,GOEBEL K,et al.An adaptive recurrent neural network for remaining useful life prediction of lthium-ion batteries[C]∥Annual Conference of thePrognostics and Health Management Society 2010.2010.

[13] PEEL L,GOLD I.Data driven prognostics using a Kalman filter ensemble of neural network models[C]∥International Conference on Prognostics and Health Management 2008.2008:1-6.

[14] GUO Z H,ZHAO H Y,LU H Y,et al.Multi-step forecasting for wind speed using a modified EMD-based artificial neural network model[J].Renew Energy,2012,37(1):241-249.

[15] 謝景新,程春田,周桂紅.基于經驗模式分解與混沌分析的直接多步預測模型[J].自動化學報,2008,34(6):684-689.

XIE J X,CHENG C T,ZHOU G H.A new nulti-step ahead prediction model on EMD and chaos analysis[J].Acta Automatica Sinica,2008,34(6):684-689.

[16] LIU Q,LIU M,ZHANG F,et al.A fault prediction method based on modified genetic algorithm using BP neural network algorithm[C]∥2016 IEEE International Conference on Systems,Man,and Cybernetics.2016.

[17] 孫強,岳繼光.基于不確定性的故障預測方法綜述[J].控制與決策,2015,29(5):769-778.

SUN Q,YUE J G.Review on fault prognostic methods based on uncertainty[J].Control and Decision,2015,29(5):769-778.

[18] 于冰潔,夏戰國,王久龍.基于高斯過程模型的異常檢測算法[J].計算機工程與設計,2016,37(4):914-920.

YU B J,XIA Z G,WANG J L.Anomaly detection algorithm based on Gaussian process model[J].Computer Engineering and Design,2016,37(4):914-920.

[19] 李軍,張友鵬.基于高斯過程的混沌時間序列單步與多步預測[J].物理學報,2011,60(7):143-152.

LI J,ZHANG Y P.Single-step and multiple-step prediction of chaotic time series using Gaussian process model[J].Acta Phys Sin,2011,60(7):143-152.

[20] 陳艷,王子健,趙澤,等.傳感器網絡環境監測時間序列數據的高斯過程建模與多步預測[J].通信學報,2015,36(10):252-262.

CHEN Y,WANG Z J,ZHAO Z,et al.Gaussian process modeling and multi-step prediction for time series data in wireless sensor network environmental monitoring[J].Journal on Communications,2015,36(10):252-262.

[21] DUVENAUD D,LLOYD J R,GROSSE R,et al.Structure discovery in nonparametric regression through compositional kernel search[C]∥Proc of the 30th International Conference on Machine Learning,Atlanta GA,USA.2013.1166-1174.

[22] HUANG N E,LONG S R,SHEN Z.The mechanism for frequency downshift in nonlinear wave evolution[J].Advances in Applied Mechanics,1996,32:59.

[23] CHEN J.Application of empirical mode decomposition in structural health monitoring:some experience[J].Advances in Adaptive Data Analysis,2009,1(4):601:621.

[24] TAIEB SB,SORJAMAA A,BONTEMPI G.Multiple-output modeling for multi-step-ahead time series forecasting[J].Neurocomputing,2010,73(10/12):1950-1957.

猜你喜歡
故障方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
故障一點通
3D打印中的模型分割與打包
奔馳R320車ABS、ESP故障燈異常點亮
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
故障一點通
主站蜘蛛池模板: 特黄日韩免费一区二区三区| 亚洲国产精品VA在线看黑人| 日韩欧美成人高清在线观看| 亚洲一区毛片| 四虎成人精品| 国产精品99r8在线观看| 风韵丰满熟妇啪啪区老熟熟女| 毛片网站在线播放| 98精品全国免费观看视频| 嫩草影院在线观看精品视频| 国产精品亚洲天堂| 19国产精品麻豆免费观看| 亚洲动漫h| 国产在线观看人成激情视频| 日韩国产精品无码一区二区三区| 欧美色伊人| 国产在线一区视频| 一级香蕉人体视频| 在线亚洲天堂| 国产成人精品男人的天堂| 成人精品午夜福利在线播放| 国产激情在线视频| 欧美不卡二区| 在线观看亚洲成人| 国产三级毛片| 99无码熟妇丰满人妻啪啪| 男人天堂伊人网| 亚洲综合久久成人AV| 欧美综合区自拍亚洲综合绿色| 成人免费午夜视频| AV无码一区二区三区四区| 亚洲一区二区约美女探花| 国产精品无码一区二区桃花视频| 青青热久免费精品视频6| 亚洲男女天堂| 九色综合伊人久久富二代| 无码电影在线观看| 亚洲国产高清精品线久久| 欧美日韩精品在线播放| 久久久久久久久久国产精品| 一级毛片在线播放免费| jizz在线观看| 国产欧美日韩18| 黄色福利在线| 91精品免费久久久| 国产精品国产三级国产专业不| 国产极品美女在线| 2022精品国偷自产免费观看| 日韩国产高清无码| 波多野结衣视频一区二区| 91精选国产大片| 欧美午夜精品| 亚洲AⅤ波多系列中文字幕| 综合亚洲色图| 国产婬乱a一级毛片多女| 亚洲成人动漫在线观看| 中文字幕av无码不卡免费| 波多野结衣视频网站| 亚洲第一黄色网| 精品自窥自偷在线看| 久久精品国产在热久久2019| 国产又黄又硬又粗| 欧美成a人片在线观看| 亚洲欧美极品| 天堂网亚洲系列亚洲系列| 久久国产香蕉| 成人福利在线视频免费观看| 欧美成人精品一级在线观看| 亚洲码一区二区三区| 久久伊人色| 狠狠综合久久| 国产好痛疼轻点好爽的视频| 蜜桃臀无码内射一区二区三区| 久草中文网| 国产在线观看第二页| 国产一级毛片在线| 天堂亚洲网| 香蕉视频国产精品人| 亚洲成a人片在线观看88| 性喷潮久久久久久久久| 少妇被粗大的猛烈进出免费视频| 精品无码一区二区三区在线视频|