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

一種新穎的深度因果圖建模及其故障診斷方法

2022-07-03 02:11:50彭開香
自動化學報 2022年6期
關鍵詞:深度故障檢測

唐 鵬 彭開香 董 潔

現代工業過程是一個高度復雜的系統.由于實際的物理連接、控制回路的作用,工業過程中的設備、部件、過程變量相互耦合,構成了復雜的互連網絡.這種互聯耦合關系使得系統某一部位一旦發生異常,將會隨著系統網絡傳播并演變演化,進而導致更加嚴重的故障.采用先進的故障檢測和診斷技術是保證工業過程安全有效穩定運行的重要手段.傳統的基于知識或者模型的方法很難構建大規模系統變量間的復雜關系.隨著工業自動化、信息化的快速發展,工業過程收集了越來越多的傳感器數據,這為數據驅動方法提供強有力的支撐.數據驅動的故障檢測和診斷方法也因此受到了廣泛關注和研究,并大量應用到化工、半導體制造等過程,尤其是多元統計過程監測 (Multivariate statistical process monitoring,MSPM) 方法[1].

MSPM 利用多元投影技術將高維觀測數據投影到低維主元子空間和殘差子空間,并設計相應的多元統計量(如T2,SPE) 及其控制限來監測數據是否超出正常工作范圍,常用的多元統計方法有主元分析 (Principal component analysis,PCA)[2]、偏最小二乘 (Partial least squares,PLS)[3]、獨立主元分析 (Independent component analysis,ICA)[4]和典型相關分析(Canonical correlation analysis,CCA)[5]等.一旦檢測出故障,需要采用貢獻圖[6]、重構貢獻圖[7]等故障隔離方法來辨識故障相關變量.這些方法由于拖尾效應[8]的影響并不能準確地辨識出所有的故障變量.隨后,格蘭杰因果分析[9]、傳遞熵[10]等因果分析方法被用來進行故障根源診斷和傳播路徑辨識,然而這些算法很難達到預期效果,且需要較長的分析時間.深度學習能夠處理大規模數據,自動地提取深層次特征,在工業領域也取得了較成功的應用.近些年來,自動編碼器 (Autoencoder,AE)[11]、深度信念網(Deep belief network,DBN)[12]、變分自編碼器(Variational autoencoder,VAE)[13]等深度學習技術越來越多地應用到過程監測中,但是由于神經網絡的黑箱特性,變量貢獻率難以計算,限制了其在故障隔離、根源診斷和傳播路徑辨識方面的應用.

圖論技術利用由若干節點和連接節點的線構成的圖模型來定性或者定量地表征變量之間關系,能夠較好地描述工業過程的復雜網絡結構[14-15].其中貝葉斯網絡 (Bayesian network,BN) 作為一種概率有向圖模型,已經應用在風險分析、可靠性可維護性分析等領域[16],在故障檢測和診斷領域也取得了較成功的應用.Mehranbod 等使用BN 進行穩定和過渡階段的故障檢測和隔離[17-18].Azhdari 和Mehranbod 驗證了BN 在田納西-伊斯曼 (Tennessee Eastman,TE) 過程中故障檢測與診斷應用的有效性[19].Gonzalez 等將BN 應用在過程監測中的維數化簡過程[20].Chen 和Ge提出了一個分層BN (Hierachical Bayesian network,HBN)建??蚣躘21],通過對工業過程進行分解,構建局部單元以及單元間的貝葉斯網絡,實現大規模工業過程的故障檢測與診斷.隨后,Chen 和Ge 又針對過程監測中低質量數據建模問題,提出了魯棒BN (Robust Bayesian networks,RBN)[22].針對工業過程中存在的動態特性,Yu 和Rashid 采用動態貝葉斯網絡 (Dynamic Bayesian networks,DBN),用異常似然指標 (Abnormality likelihood index,ALI) 和動態貝葉斯概率指標 (Dynamic Bayesian probability index,DBPI)分別進行故障檢測和根源診斷[23].Zhang 和Dong 將高斯混合模型(Gaussian mixture model,GMM)和三時間切片DBN 整合來解決過程監測中存在的數據缺失和非高斯問題[24].

BN 為過程變量的因果關系提供了條件概率表示,但是,該條件概率關系一般是線性的,無法描述過程中存在的非線性特性.同時,大多數基于BN 的方法需要通過過程知識得到BN 的網絡結構,這在很多復雜工業過程中是很難確定的.BN 作為一種有向無環圖模型,也較難表示系統中存在的閉環結構.為解決上述問題,并考慮系統中存在的動態特性,本文提出了一種深度因果圖 (Deep causality graph,DCG) 建模方法,利用多層感知器 (Multilayer perceptron,MLP) 和門控循環單元 (Gate recurrent unit,GRU)對每一個過程變量建立概率預測模型.在模型訓練過程中,引入組稀疏懲罰項,自動地檢測變量間因果關系,從而得到過程變量的因果有向圖結構以及定量的條件概率表征;然后基于DCG 模型的條件后驗概率分布建立單變量監測統計指標,并通過貝葉斯推理融合,構建綜合的監測統計量,實現工業過程的整體監測.進一步通過計算變量貢獻度指標,隔離出故障相關變量;最后根據深度因果圖模型獲得的有向圖結構,診斷出故障根源,并辨識故障的傳播路徑.

論文的結構如下:第 1 節詳細介紹了深度因果圖推導和建模過程;第 2 節提出了基于深度因果圖模型的故障檢測和診斷框架;隨后,在第 3 節利用TE 過程數據對所提算法進行了驗證,并在最后一節中進行了總結.

1 深度因果圖建模方法

1.1 圖結構已知的因果關系建模方法

考慮一個有向圖模型,其網絡結構為G=〈V,E〉,V表示節點的集合,E表示連接節點的有向邊集合.在工業過程中,節點可以指過程變量特征,有向邊對應著連接節點之間的因果關系,有向邊的首和尾連接的節點分別表示為父節點和子節點.

對一個有n維觀測變量的工業過程,t時刻觀測變量表示為xt=[x1,t,x2,t,···,xn,t], 則節點i在t時刻觀測變量為xi,t, 其父節點集在t時刻的狀態用xpa(i),t表示.考慮系統的動態特性,觀測變量xi,t不僅與當前時刻父節點狀態有關,而且與歷史時刻自節點和父節點的狀態有關.t時刻觀測變量xi,t用節點i的歷史觀測數據xi,t-T:t-1和其父節點集的觀測數據xpa(i),t-T:t非線性表示:

式中,εi,t為隨機噪聲項.利用概率形式來表示這種變量間的依賴關系:

式中,xc,t-T:t為xpa(i),t-T:t和xi,t-T:t-1組合的觀測值向量;分別為xi,t的后驗概率分布的均值和方差,均為非線性函數,其非線性關系可以用神經網絡表示.

通過神經網絡構建節點i預測模型,其條件概率對數似然表示為lnp(xi,t|xc,t-T:t).設預測模型參數集為 Θi,模型參數Θi學習的目標是最大化對數似然期望值,目標函數表示為:

由上述討論可知觀測變量xi,t由節點i與其父節點的歷史狀態以及父節點的t時刻狀態共同決定.將節點i及其父節點的歷史狀態中與變量xi,t相關的動態特征信息表示為zi,t-1, 父節點當前時刻的相關特征信息表示為hi,t,式 (1) 中的非線性預測模型于是可以轉變為以下形式:

當前時刻相關特征hi,t的信息提取可以利用多層感知器實現,歷史時刻的動態特征信息zi,t-1采用GRU 進行特征提取.GRU 作為循環神經網絡 (Recurrent neural network,RNN) 的一種變體能夠提取時間序列的動態特征,同時有效地解決標準RNN 算法存在的長期記憶和反向傳播中的梯度問題.GRU 的詳細介紹可以參考文獻[25].

1.2 圖結構未知的因果圖建模方法

上述節點預測模型是建立在因果有向圖結構已知的情況.而在很多復雜工業過程中,變量之間的因果關系是未知的.為了尋找每個變量的因變量,本文首先將除自變量外的其他變量都默認為因變量,然后在模型訓練過程中引入Group Lasso 懲罰項,使得節點預測模型中與輸入變量相關的連接稀疏化,使用盡可能少的輸入獲取盡可能準確的預測,進而實現變量之間的因果關系自動檢測,其單一節點i的后驗概率預測模型如圖1 所示.節點的概率分布預測模型由兩個動態特征提取單元和一個預測輸出單元構成.在第一個動態特征提取單元中,輸入為t時刻觀測變量值xt={xi,t,xi-,t},xi-,t表示除節點i之外其他節點的觀測值.為在t時刻從輸入xi-,t提取的m1維隱藏特征,其過程表示為:

圖1 深度因果圖的單節點預測模型網絡結構Fig.1 The network structure of single node prediction model for deep causality graph

式中,σ為激活函數;W1表示第一層全連接層的連接權重,是一個m1×(n-1)大小的權重矩陣.將W1展開為則式(5) 重寫為:

為了獲取深層次的特征,增大節點特征提取的視野域,以便用較少的變量連接獲取準確的預測結果,將除自節點外的其他節點提取的動態特征合并作為第二個動態特征提取單元的輸入,再進一步地提取節點的動態特征.其網絡結構與第一個動態特征提取單元一致,第二層全連接層連接權重為其中

預測輸出單元將t-1時刻的動態特征以及t時刻隱藏特征合并作為輸入,通過多層感知器輸出節點i在t時刻觀測變量預測概率分布的均值μi,t和方差對數

組最小絕對收縮和選擇算子(Group least absolute shrinkage and selection operator,Group Lasso) 懲罰項被選擇為稀疏懲罰項加入到單節點的概率預測模型目標函數中.Group Lasso 是Yuan 在 2006 年在Lasso 稀疏約束方法上到組上的推廣[26].它根據輸入變量對權重矩陣進行分組,然后在目標函數中懲罰每一組的L2 范數,這樣就可以將部分組的全部系數同時消成零,即抹掉整組的變量,這種方法叫做Group Lasso 分組最小角回歸算法.對節點i的預測模型網絡權重矩陣W1和W2,根據連接的過程變量節點劃分成n-1組.與過程變量節點j對應的權重范數為加入Group Lasso 懲罰項后,節點i的預測模型通過最小化目標函數(7) 來迭代更新網絡參數:

超參數λ控制著連接權的稀疏性,λ越大,越多組的權重會被消零,使其節點的因果關系連接呈現為稀疏集.

由于在單一節點預測過程中需要利用其他節點提出的動態特征信息,因此在模型訓練過程中需將n個過程變量節點的預測模型作為整體進行訓練,故深度因果圖整體模型的目標函數是n個目標函數的和形式:

在深度因果圖模型訓練過程中,首先要用大小為T+1的滑動窗口將標準化后的時間序列訓練數據集規范為Xtrain∈RNtrain×(T+1)×n.由于目標函數中存在Group Lasso 懲罰項,其在權重范數為 0 處不可微.為解決目標函數存在不可微的問題,采用了近端梯度下降 (Proximal gradient descent,PGD) 算法進行模型優化訓練.對包含可微和不可微函數的目標函數f(θ),優化問題分解為:

式中,g(θ)是可微函數,h(θ) 是不可微函數.

定義一個近端映射函數:

參數θ的迭代過程可以表示為:

對給定的一個批次的訓練數據,PGD 的一次迭代算法可以理解為:給定優化參數起點θ(k-1), 可微函數g沿著起點的負梯度方向,做步長為tk的梯度下降得到一個預更新值,然后使用近端映射尋找一個z,使得不可微函數h(θ)足夠小,并接近預更新值,用z作為本次迭代的更新值θ(k).

在模型訓練中,選擇軟閾值算子作為PGD 的近端算子.該優化算法可以使非因果變量的連接權重精確地收斂到 0,以便能夠用權重值來解釋變量之間的因果性.

2 基于深度因果圖模型的故障檢測與診斷

2.1 故障檢測

完成深度因果圖模型構建后,可以利用該圖模型得到的條件后驗概率分布設計過程監測統計量.首先,對每一個變量節點設計單變量統計指標.考慮觀測變量為n維的動態連續過程,t時刻觀測變量表示為xxx1:n,t={x1,t,x2,t,···,xn,t}.節點i有實際觀測值xi,t.根據深度因果圖模型,計算其預測值概率分布為其平方馬氏距離服從自由度為 1 的χ2分布:

為了獲取一個綜合評價指標來監測整個過程,本文利用貝葉斯融合對各節點檢測信息進行融合.根據貝葉斯推理,變量節點i的統計量為故障的概率表示為:

調節因子ηp可以調整概率值對故障狀態的靈敏度,在本文中ηp設置為 1.

根據上述推導過程,可以將單節點的監測統計指標轉換為概率形式表示.當1-α時,節點i發生異常.然后對單變量概率監測指標加權平均,來計算綜合監測指標PD2:

對收集的歷史正常操作數據,計算每個時刻的綜合監測指標.利用核密度估計方法,可以估計PD2的控制限為在線故障檢測階段,可以直接通過判斷PD2是否超出控制限來判斷該過程是否發生故障.

2.2 故障診斷

根據深度因果圖模型的稀疏化的連接權重矩陣確定過程變量之間的定性因果關系,構建描述過程變量因果關性的因果矩陣以及相對應的因果有向圖.利用有向圖可以有效地分析故障的傳播路徑,辨識故障根源.但是,深度圖學習模型的節點預測考慮到了其他節點當前時刻的輸入,同一時刻不同變量之間可能存在強相關性,使得兩個變量互為因果.對深度圖學習模型中存在的互為因果的變量進行格蘭杰因果分析,可以去除部分錯誤的連接.

一旦檢測出故障,需要設計合理的變量貢獻圖指標來辨識出故障相關變量.設置指標τi,t為一個二進制數,表示t時刻節點i的概率檢測指標是否超出設定的閾值αc,如下式所示:

在一設定的時間周期T1內節點i的概率檢測指標超出閾值的次數累計和為:

然后,變量貢獻度指標(Variable contribution index,VCI)設計為如下形式:

由于系統噪聲的影響,即使在正常運行狀態下單變量檢測指標也會有超出設定閾值的可能性.設定一個 VCI 閾值為1/n,去除掉異常頻度較小的變量,將 VCI大于閾值的節點被辨識為故障相關變量.然后,根據深度因果圖模型得到的因果關系矩陣,確定由故障相關變量構成的局部因果有向圖網絡.根據分析網絡的因果連接,定性地確定出故障的根源和傳播路徑.在進行根源診斷時,故障根源一般為沒有其他異常變量影響的變量,即在由故障相關變量構成的局部因果有向圖中,將局部有向圖中的根節點變量辨識為故障的根源.

實際模型中,由于岔路的遍布,很多故障在追溯源頭的過程中存在多條傳播路徑,在利用局部有向圖進行根源診斷時也會存在多個根節點或者根節點不存在的情況.在使用深度因果圖模型進行故障檢測過程中,已經提供了每一個變量的故障概率對辨識出來的故障相關變量,可以通過對故障發生后T1時間段內的故障概率均值,以此來衡量每一個故障相關變量的故障程度.然后,再結合故障相關變量的局部有向圖,以概率均值為依據,確定最有可能的故障傳播路徑,進而在存在多個根節點或者根節點不存在的情況確定其根源.

2.3 故障檢測和診斷流程

對深度因果圖模型獲得后驗概率分布進行統計分析,并利用模型得到的定性因果關系,所提的方法可以確定故障發生的時間,并隔離故障相關的變量,以及辨識出故障的根源和傳播路徑.基于深度因果圖模型的故障檢測和診斷框架包含離線建模和在線監測兩部分:

離線階段:

1) 獲取歷史正常操作時間序列數據并標準化,滑動窗口時間片處理;

2) 近端梯度下降法迭代訓練深度因果圖模型直至收斂;

3) 計算歷史數據集的綜合監測指標PD2;

在線階段:

1) 獲取t-T到t時刻在線觀測數據,標準化處理;

2) 通過深度因果圖模型計算單變量節點后驗概率分布,計算單變量監測指標

3) 計算全過程綜合監測指標PD2;

5) 計算變量貢獻度指標VCI,確定故障相關變量;

6) 使用故障相關變量構建局部因果有向圖模型,辨識故障根源和傳播路徑.所提框架的算法流程圖如圖2 所示.

圖2 基于深度因果圖模型的故障檢測和診斷框架Fig.2 The fault detection and diagnosis based on deep causality graph model

3 案例研究

3.1 TE 過程描述

本文用TE 過程來驗證算法有效性.TE 過程是由美國化工公司的Downs 和Vogel 提出的一個仿真過程模型,被廣泛應用于過程控制和監測方法的性能評估中,其工藝流程圖如圖3 所示,由五個生產單元構成,分別是冷凝器、反應器、分離器、汽提塔和壓縮機.TE 過程仿真包括了41個測量變量和 12個控制變量,其測量變量中 22個是連續過程變量, 19個是成分測量值.本文采用Bathelt 提供的仿真模型進行數據生成,模型的采樣周期為 36秒,選擇測量變量XMEAS1~XMEAS22, 控制變量XMV1~XMV4、XMV6、XMV7、XMV9、XMV10 共31個變量作為監測變量.在正常模態下仿真模型運行1000個小時,生成 105個正常操作數據樣本,作為訓練集.對 21種故障模型,分別運行 10個小時,并在第 3 小時引入故障,即每個故障樣本量為 1000個,故障發生在第301個樣本及之后,故障檢測置信度設置為α=95%.

圖3 TE 過程工藝流程圖Fig.3 The flowchart of TE process

3.2 實驗結果分析

深度因果圖學習模型參數設置如下:1) 模型特征層維數設置為30;2) 圖4 展示了超參數λ在[0.3,0.8]范圍以0.1為間隔選取下的變量連接數和預測誤差的關系,隨著λ的增大,預測誤差跟隨增大,變量連接數量隨之減小.當超參數λ為 0.5 時近似為曲線拐點,能較好保證預測誤差盡量小的同時盡可能減少變量連接數量;3) 滑動時間窗大小設置為10.硬件平臺使用了Core i7 4.2 GHz 四核CPU,NVIDIA GTX 1070 6 GB 顯存,16 GB RAM,采用Pytorch進行算法仿真.

圖4 變量連接數和預測誤差的關系曲線Fig.4 The relation curve between prediction error and the number of variable connections

利用標準化的歷史正常操作數據集對深度因果學習模型進行迭代訓練.根據稀疏化連接權重可以確定變量之間的因果關系,從而確定TE 過程的有向圖結構.格蘭杰因果分析進一步對互為因果的變量進行分析,去除部分虛假連接.TE 過程的過程變量因果關系構成的因果矩陣如表1所示,表的行名對應果變量,列名對應因變量.在單變量統計指標控制限設定中,核密度估計的帶寬為 10-3,并根據正常操作歷史數據,通過核密度估計確定控制限

在線故障檢測中, 21 個故障集的故障檢測率(Fault detection rate,FDR)結果如表2 所示,平均故障誤報率(False alarm rate,FAR)為2.7%,單個樣本的平均計算時間為10 ms.可以發現基于深度因果圖學習模型的故障檢測方法可以快速有效地檢測出大部分類型的故障.由于所提出的故障檢測方法是基于實際觀測值在預測概率分布中的偏離程度,當故障程度較小時,變量之間的因果依賴關系可能沒有被破壞,實際的預測值仍然能夠跟蹤變量變化,所以所提故障檢測方法不能有效地檢測出故障5、9、15、21 這幾種微小故障.

表2 21個故障類型的FDRs (%)Table 2 The FDRs of 21 faults (%)

故障 4 是反應器冷卻水入口溫度的一個階躍變化故障,其直接導致影響反應器內部溫度(v9),然后反應器冷卻水出口溫度(v21)也會上升.但是,由于反饋控制,反應器冷凝水流量增大,補償反應器冷卻水入口溫度的上升造成的影響.反應器溫度會逐漸恢復正常,冷卻水出口溫度也趨于正常.采用深度因果圖模型進行故障檢測的結果如圖5所示,由圖可以看出檢測指標PD2在第 301 個采樣點及時檢測出過程異常.利用故障發生后的 50 個連續樣本點,計算所有變量的VCI,結果如圖6 所示.VCI 大于 1/n的變量為故障相關變量,由圖可知v9、v11、v16、v21、v27、v30共6個變量被診斷為故障相關變量.根據深度因果圖模型提供的因果矩陣,利用這 6 個故障相關變量構建局部因果有向圖網絡,如圖7 所示,節點的數值表示 50 個連續樣本點的平均故障概率,反映節點的故障程度.根據故障程度的大小,確定主要故障傳播路線為根節點是變量v9,即反應器內部溫度異常為故障 4 發生的根源,這與故障 4 的機理分析結果一致.

圖5 故障4 的故障檢測結果Fig.5 The fault detection result for Fault 4

圖6 故障4 的VCI 圖Fig.6 The plot of VCI for Fault 4

圖7 故障4 的故障相關變量的因果關系Fig.7 The causalities among fault-related variables for Fault 4

故障 8 是A、B、C 進料量的一個隨機變化故障.由圖8的檢測曲線可知過程異常在第 354 個采樣點被檢測出來.根據圖9 所示的VCI,確定故障相關變量為v6、v7、v10、v16、v20、v25、v27、v29共 8 個變量.根據表1 的因果矩陣,用故障相關變量構建局部因果有向圖網絡,如圖10所示.分析構建的有向圖網絡,故障的主要傳播路徑為v25→v16→v20→{v6,v29},v25辨識為故障根源變量.v25是A 進料控制量,這和故障導致的原因一致.根據TE 過程的故障描述,故障8會導致反應器壓力 (v7) 和汽提塔壓力 (v16) 發生明顯變化,由于閉環反饋調節作用,排放閥 (v27)、壓縮機工作功率(v20)也跟隨變化,從而導致反應器進料率(v6) 發生異常.圖10 可以反映出這一故障傳播路徑.

圖8 故障8 的故障檢測結果Fig.8 The fault detection result for Fault 8

圖9 故障8 的VCI 圖Fig.9 The plot of VCI for Fault 8

圖10 故障8 的故障相關變量的因果關系Fig.10 The causalities among fault-related variables for Fault 8

表1 TE 過程的因果矩陣Table 1 The causality matrix of TE process

通過對兩個故障案例的故障檢測和診斷結果分析,發現所提出的深度因果圖模型能夠及時有效地檢測出故障,并且在沒有經驗知識來確定過程變量的有向圖結構的情況下,仍然能夠提供較為精準的故障隔離和故障根源診斷結果.

4 結論

本文提出了一種新穎的深度因果圖建模方法,通過在模型訓練目標函數中引進Group Lasso 稀疏懲罰項,可自動地獲取變量之間的因果關系.不同于貝葉斯網絡,該深度因果圖無需先驗過程知識即可得到過程變量的網絡結構及其條件概率分布.然后,基于該深度因果圖模型計算每個變量的檢測統計指標,通過融合得到綜合檢測指標,從而實現整體過程的故障檢測.一旦檢測到故障,分析變量貢獻指標和變量因果關系,來隔離故障相關變量,診斷故障根源,并辨識故障傳播路徑.將所提出的方法在TE 過程驗證.實驗結果表明所提方法能夠有效地檢測出故障,辨識故障傳播路徑,并定位故障的根源.深度因果圖模型的因果有向圖結構受訓練模型中稀疏懲罰項的超參數影響,合理的超參數能夠減小網絡結構的復雜性,提升故障根源診斷的精準度.由于沒有利用先驗過程知識,確定的因果有向圖結構可能與實際網絡結構有差別.在后期的研究中,考慮將先驗知識引入到深度因果圖模型建模中,以獲得更加精準的網絡結構,提升故障診斷性能.

猜你喜歡
深度故障檢測
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
深度理解一元一次方程
故障一點通
深度觀察
深度觀察
深度觀察
奔馳R320車ABS、ESP故障燈異常點亮
小波變換在PCB缺陷檢測中的應用
主站蜘蛛池模板: 国产91无码福利在线| 亚洲无码视频喷水| 一本大道东京热无码av| 成人自拍视频在线观看| 亚洲第一黄色网址| 色婷婷成人| 国内精自视频品线一二区| 黄色片中文字幕| 香蕉色综合| 成人综合网址| 一本大道香蕉久中文在线播放 | 人妻丰满熟妇αv无码| 国产小视频网站| 福利一区三区| 欧美亚洲综合免费精品高清在线观看| 欧美一区二区三区欧美日韩亚洲| 亚洲AV无码乱码在线观看代蜜桃| 午夜精品福利影院| av无码久久精品| 中文字幕人妻无码系列第三区| 午夜福利无码一区二区| 少妇精品在线| 国产激情国语对白普通话| 午夜欧美理论2019理论| 日本91视频| 风韵丰满熟妇啪啪区老熟熟女| 性视频一区| 人妻出轨无码中文一区二区| 99尹人香蕉国产免费天天拍| 色哟哟国产精品一区二区| 高清无码一本到东京热| 国产 日韩 欧美 第二页| 91成人试看福利体验区| 亚洲欧美h| 国产精品内射视频| 在线看片免费人成视久网下载| 午夜精品福利影院| 在线观看网站国产| 国产一二三区视频| 国模极品一区二区三区| 亚洲精品色AV无码看| 91精品国产91久无码网站| 国产精品毛片一区| 欧美高清日韩| 欧美性天天| 精品无码一区二区三区在线视频| 制服丝袜无码每日更新| 老司机午夜精品网站在线观看| 色妞www精品视频一级下载| 久久网欧美| 高潮毛片免费观看| 亚洲综合国产一区二区三区| 伊人狠狠丁香婷婷综合色| 白浆视频在线观看| 九九精品在线观看| 欧美一级高清片欧美国产欧美| 在线观看国产小视频| 日本黄色a视频| 国产精品亚欧美一区二区| 国产精鲁鲁网在线视频| 国产激情第一页| 国产精品女熟高潮视频| 69视频国产| 久久国产香蕉| 久久美女精品| 免费 国产 无码久久久| 毛片国产精品完整版| 国模沟沟一区二区三区| 亚瑟天堂久久一区二区影院| 亚洲精品亚洲人成在线| 日韩精品专区免费无码aⅴ| 国内老司机精品视频在线播出| 亚洲小视频网站| 欧美日韩中文国产va另类| 91欧美亚洲国产五月天| 国产精品xxx| 99久久免费精品特色大片| 亚洲国产成人超福利久久精品| 91在线视频福利| 国产高清在线丝袜精品一区| 精品少妇人妻无码久久| 亚洲欧美人成电影在线观看|