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

基于變量概率信息的因子分析監控方法

2017-07-18 11:43:33胡婷婷王帆侍洪波
化工學報 2017年7期
關鍵詞:故障信息模型

胡婷婷,王帆,侍洪波

(華東理工大學化工過程先進控制和優化技術教育部重點實驗室,上海 200237)

基于變量概率信息的因子分析監控方法

胡婷婷,王帆,侍洪波

(華東理工大學化工過程先進控制和優化技術教育部重點實驗室,上海 200237)

因子分析(factor analysis,FA)將噪聲因素加入到建模過程中,可通過最大期望(expectation maximum,EM)算法建立模型。傳統的 FA(ST)指標僅利用了變量的期望信息而忽略了更能代表不確定性的方差信息,這可能會導致故障的漏報。通過對過程變量的概率分析,從本質上揭示了FA(ST)的這一缺陷。建模過程中的另一個重要因素是確定因子個數,使得在降維的同時能最大程度地保留對過程有用的信息。針對傳統監控指標信息不足的問題,提出的負對數似然概率(negative log likelihood probability,NLLP)指標整合了更全面的概率信息;針對因子個數給定的問題,提出了一種整體-局部因子數確定法,使得因子和變量對于過程的信息解釋率都達到收斂。最后通過數值例子和Tennessee Eastman(TE)過程驗證了所提方法的有效性和優越性。

因子分析;參數估計;過程控制;信息解釋率;統計分析

引 言

過程監控對于提高工業生產過程的安全性和穩定性具有重要的意義[1]。分散控制系統(distributed control system,DCS)技術的發展促進了生產過程數據的采集和存儲[2-3],多元統計過程監控(multivariate statistical process monitoring,MSPM)方法可以從具有高維度、強耦合等特性的大量數據中挖掘出有利于過程監控的信息[3-5]。目前應用最廣泛的MSPM方法如主元分析(principal component analysis,PCA)方法[6-10]是在主元空間和殘差空間分別建立統計量,而二者的量度不同,往往會造成監控結果不一致的問題[11-13]。另一方面,環境噪聲等非確定性因素是客觀存在的,忽略噪聲因素而建立的模型是不準確的。概率潛變量生成模型可以同時解決這兩個問題。概率主元分析(probabilistic principle component analysis,PPCA)[13]和因子分析[14]是兩種噪聲條件假設下的概率潛變量模型,其中FA釋放了PPCA中噪聲各向同性的約束,更符合實際工況,應用于過程監控領域可取得較 PPCA更好的效果。由于生成模型的特點可以使用一個綜合指標來進行過程監控,從而避免了監控結果不一致的缺陷[15]。Ge等[16]將ICA與FA結合來處理含有非高斯成分變量的過程,介紹了一種具有新型相似因子的ICA-FA監控方法。文獻[17]在FA的基礎上,提出了極大似然混合因子分析模型(MLMFA)以解決噪聲因素、非高斯成分和多模態問題。文獻[18]在混合因子分析(MFA)的基礎上,提出了一種直接極大似然估計方法。Zhao等[19]在FA框架的基礎上,從過程數據的內在特性和數據丟失現象出發,說明了方差信息對于過程監控的必要性。然而以上列舉方法中的監控統計量僅僅利用了變量的期望信息,而更能代表非確定性的方差信息卻被忽略了。

用EM算法建立FA模型時,需給定因子數,目前確定因子數的方法[11,20-23]中應用最多的是基于信息解釋率的方法,但是這些方法只能從整體或者局部上保證方差貢獻度達到收斂。MSPM過程包含了各種量綱不一致的變量,且這些變量之間存在著某種聯系,使得整體解釋率和單個變量解釋率之間存在不對等的關系,然而整體的信息解釋率并不能理解為單個變量解釋率的簡單疊加,這種聯系必然與整個系統的結構和功能等有關。為確保整體和變量的信息解釋率都達到收斂,本文介紹了一種基于FA模型的可自動確定因子數的整體-局部因子數確定法。此外,本文從本質上證明了傳統的FA(ST)指標只包含變量的期望信息而忽略了更能代表不確定性的方差信息,為了克服這一不足,本文利用增加了變量方差信息的 NLLP[24]來構造監控統計量。由于NLLP是過程變量概率分布的函數,因此其包含了更全面的概率信息,且只用一個監控量可以避免多監控圖監控結果不一致的問題,其監控控制限用核密度估計(kernel density estimation, KDE)得到。最后,通過兩個仿真驗證了本文所提算法的有效性。

1 基于FA(NLLP)模型的過程監控

設經過歸一化預處理后的過程觀測數據x ∈?d×n,其中d為過程變量個數,n為樣本觀測數,生成模型將x看成是由不相關的因子z的線性組合外加服從一定高斯分布的噪聲變量e組合而成[12],該生成模型認為因子z是輸入,觀測變量x是輸出,即

其中,A∈?d×m為載荷矩陣,m<d為主元個數,z∈?m,e ∈?d,同時假設z~N(0,I),e ~N(0,Ψ),則x~N(0,C),C=AAT+Ψ,其中Ψ=diag(λ1,λ2,…,λd)為噪聲協方差矩陣。若λ1=λ2=…=λd,則 FA退化成 PPCA;當λ1=λ2=…=λd≈0時,PPCA 就退化成 PCA[19],因此FA的限制條件更少,更能反映過程數據的本質特性。

1.1 因子數的確定

為了建立精確的FA模型,需合理地選擇因子數,利用因子數增加到一定程度時信息解釋率收斂的性質[11],本文提出整體-局部的因子數選取法,保證整體解釋率和單個變量解釋率都達到收斂。

因子數選取步驟如下。

(1)根據文獻[25]可確定因子數m的上限m0

(2)計算因子數l取不同值時的整體解釋率per1(l)

其中,Al和Ψl為不同l值下的由EM算法估計出的參數值,為因子數為l時的所有變量的總信息,為當前因子數所包含的信息,tr()為求跡算子。

(3)計算整體解釋率的增長值rep1(l)

(4)判斷整體解釋率的收斂性,確定整體解釋率達到收斂的因子數m1

其中,為足夠小的數。判斷式(5)是否成立,選擇式(5)成立時的最小因子數l為m1。

(5)計算因子數l取不同值時對變量i的解釋率per2(l,i)

式中,diag()i表示取矩陣對角線上第i 個元素。

(6)計算每個變量解釋率的增長值rep2(l,i);

(7)判斷每個變量解釋率的收斂性,確定每個變量解釋率都達到收斂的因子數m2

式中,li表示第 i個變量解釋率達到收斂時的最小因子數。

(8)確定綜合因子數m

需要說明的是,步驟(2)~步驟(4)為計算整體信息解釋率過程,步驟(5)~步驟(7)為計算局部信息解釋率過程,這兩部分沒有先后順序之分,實際操作時同時進行。

1.2 FA概率生成模型的建立

由式(1)可知,給定因子數m后,該FA模型就由A和Ψ決定。對這些參數的估計可通過最大化測量值的似然函數來確定,但是直接進行極大似然計算存在困難,因為A和Ψ都與因子z有關,而z又是未知的,能得到的僅是觀測數據 x,EM 算法將因子z視為遺失變量,利用完整數據集(x, z)來進行參數估計,以期實現參數的極大似然估計[14,18-19]。EM算法一般分為E步(求期望)和M步(最大化),具體計算公式見文獻[14]。

1.3 監控指標的確定

1.3.1 過程變量的概率分析 根據概率生成模型,可以推導出過程變量x在第i個采樣時刻的期望和方差

其中,β=ATC-1,E[·]和分別表示·的期望和方差。由式(11)和式(12)可以看出,過程變量的方差在任意采樣時刻都保持恒定,而期望與采樣時刻有關,是隨著觀測值變化的。

1.3.2 基于過程變量的綜合監控指標 由歷史正常數據離線建立好FA模型后,可以直接用一個綜合指標ST實現對過程的在線監測。ST統計量的定義如下[11]

其中,ST服從自由度為d的χ2分布,控制限(d)表示χ2分布χ2(d)的上側分位數,為置信度。從式(13)可知 ST是過程變量對其空間中心馬氏距離的檢驗,但公式中僅僅利用了過程變量觀測值的在線期望,而不隨觀測值變化的方差信息卻被忽略了。在線監測時,即使實際過程變量的方差因為故障的發生而導致變化了,由于 ST本身的缺陷使得故障數據的特征仍符合正常模型,這可能會造成更多的漏報。

1.3.3 基于NLLP的監控指標 基于以上的討論,為了獲得更好的監控性能,除了利用變量的期望信息外,還需要將方差信息也加入到監控指標的構建中,本文使用NLLP作為監控統計量。對于一個給定的FA模型,NLLP衡量了新測試樣本與FA模型之間概率分布的匹配程度。當故障發生時,故障樣本的NLLP值會明顯高于正常點的NLLP值,NLLP定義為

從NLLP的定義式(14)可以看出,NLLP是關于觀測變量x的概率分布的函數,其必然同時包含期望和方差信息,可檢測出導致變量期望和方差變化的故障。相較于ST,NLLP包含更全面的過程變量概率信息,能更有效地檢測故障的發生。

1.4 算法步驟

基于 FA(NLLP)算法的過程監控方法具體步驟如下。

(1)對正常樣本數據集進行歸一化預處理;

(2)由式(2)~式(10)確定因子數m;

(3)用EM算法建立FA模型;

(4)計算正常樣本的NLLP并用KDE估算控制限limNLLP(α);

(5)對新的測試數據集進行歸一化預處理;

(6)計算測試數據的 NLLP值,判斷其是否超過控制限。

其中,步驟(1)~步驟(4)為離線建模過程,步驟(5)~步驟(6)為在線監測過程。

2 仿真實驗

本節將進行以下兩個例子的仿真以驗證本文所提算法的有效性,一個是數值例子,另一個是Tennessee Eastman(TE)過程。為了比較結果的公平性,因子數的選取方法均采用本文提出的整體-局部確定法,將FA(ST)和FA(NLLP)仿真結果進行對比。

2.1 數值仿真

為了說明方差信息對于過程監控的必要性,Zhao等[19]提出了一個數據結構模型,本節將采用這一模型來驗證本文所提方法的有效性。具體模型結構如下

其中,z ~N(0,0.4),過程噪聲e1和e2都服從N(0,0.1)。很顯然,這兩個變量的過程僅含有一個因子,由結構式(16)產生的1000個正常訓練樣本用于建立FA模型。故障數據的產生通過在過程變量上施加如下故障

其中,h和Ξ分別表示故障幅度和方向。假設Ξ是A的第1列,h=1.8。故障在第51個采樣點處加入,共采集100個測試樣本,前50個樣本為正常樣本,后50個樣本為故障樣本。分別用ST和NLLP作監控統計量,其控制限分別由χ2(d)分布和 KDE確定,取99%置信度。表1列出了這兩種方法的誤報率和漏報率,監控圖如圖1所示。

顯然該故障導致了變量方差的變化,由圖1可以看出,雖然 FA(ST)能迅速地檢測到故障的發生,但是在故障持續期間仍有很多采樣點的 ST值落到了控制限以下,造成了高的漏報,從表1的數據可以看出FA(ST)的漏報率高達46%。相比而言,FA(NLLP)的漏報率降低到4%,從圖1的部分放大圖可以看出,只有第57和第78個點的NLLP值落到控制限以下,這說明 FA(NLLP)能更有效地檢測出故障。

表1 數值例子的誤報率和漏報率Table 1 Fault alarm rate and miss alarm ratein numerical example

圖1 FA(ST)和FA(NLLP)對數值例子的監控結果Fig 1 Monitoring result of FA(ST) and FA(NLLP) methods in numerical example

2.2 TE過程仿真

TE過程[28-29]是一個復雜的多變量化工反應過程,已經被廣泛地用于評價MSPM方法的有效性。TE過程共有41個測量變量(22個連續過程測量變量和19個非連續過程測量變量)和12個控制變量。在仿真過程中,控制變量5、9和12保持不變,因此本節選取剩下的9個控制變量和22個連續過程測量變量進行仿真實驗。TE過程運行一次的時間為48h,采樣間隔 3min,因此運行一次能采集到 960組數據。TE過程預設21種故障,每一種故障都在480min后開始加入,即故障數據從第161次采樣開始產生,這樣就得到了21種故障類型的960組測試數據,其中故障3、9、15的測試數據在均值、方差上沒有明顯的變化,所以統計方法的故障監測率都很低[6],本節不對這3種故障進行分析。

本節的TE過程中n=960,d=31,根據式(2)計算因子數的上限m0=23。用EM算法分別估算出因子數從1~23的FA模型參數 A和Ψ,根據式(3)計算出不同因子數下的信息解釋率,如圖 2所示。

根據式(4)和式(5)確定m1=14。由圖2也可以看出,因子數達到14及以上時,整體解釋率收斂并達到了 95%以上。根據式(6)計算出不同因子數下每個變量的信息解釋率,如圖3所示。

根據式(7)~式(9)確定出m2=18。由圖3可以看出,因子數達到18及以上時,每個變量的解釋率都收斂到70%以上。最終選取的因子數確定為m=18。需要說明的是,計算不同因子數下的 FA模型參數時都要通過EM算法,由于EM算法是局部最優算法,使得解釋率曲線會出現波動的情況,但總體趨勢是上升收斂的[30]。另外每次得到的m1和m2的值會在一定小范圍內波動且相對大小也會不同,即因子數m會有不同,但總能保證在該因子數下的總體和局部信息解釋率都達到收斂。

圖2 不同因子數下的整體信息解釋率Fig.2 Global information explanation ratio under different numbers of factors

圖3 不同因子數下的每個變量的信息解釋率Fig.3 Variable information explanation ratio under different numbers of factors

確定好因子數m后,利用EM算法獲得該因子數下的FA模型參數的估計值,用NLLP建立統計量,取0.99置信度,用KDE確定控制限。測試數據由48 h采集到的960個采樣點組成,其中前160個為正常樣本,后800個為故障樣本點,這樣的數據大小和規格被普遍應用于比較算法性能,從而保證了算法比較的公平性。表2列出了PPCA(NLLP)、FA(ST)和FA(NLLP)這3種方法對于TE過程的 18個故障的漏報率,黑體部分為監測結果的最佳值。

表2 TE過程故障數據的漏報率Table 2 Miss alarm rate of fault database in TE process

從表2的監控效果數據可以看出,基于FA模型的監控性能要優于基于PPCA的監控性能,而FA(NLLP)的故障監測性能最佳。對比FA(NLLP)與PPCA(NLLP),故障5和10的漏報率分別減少了47.875%和43%,這兩種方法對故障5的監測圖如圖4所示,其他故障類型(故障11、16、17、19、20和 21)的故障檢測率也有明顯的提高,這說明FA模型的約束條件更少,更符合實際過程的本質特性,監控性能更好。

對比FA(NLLP)與FA(ST),對于故障1、4、5、6、7和14,這些故障的發生造成了過程數據期望的明顯變化,所以這兩種方法的故障檢測率都達到了100%。此外,這兩種方法對故障2、8、12的漏報率都在2%以內,且FA(NLLP)的故障漏報率都略好于FA(ST)。對于其他的故障類型,如故障10、11、16、19、21,FA(NLLP)的監控性能較FA(ST)提高得較明顯,其中故障16的FA(NLLP)性能提高了10.375%,其他故障的FA(NLLP)性能都提高了5%左右。從整體上來看,FA(NLLP)的故障監測效果最佳,這是因為NLLP統計量作為概率分布的函數,包含了變量的期望和方差信息,監測性能要優于比只包含了期望信息的ST統計量。

圖4 故障5的監測結果圖Fig 4 Monitoring result of PPCA(NLLP) and FA(NLLP) for fault 5

3 結 論

在FA模型的基礎上,本文提出了一種可自動確定因子數的方法,同時建立了包含更全面的變量概率信息的統計量 NLLP。從因子和變量的信息解釋率角度出發,本文提出的因子數確定法能保證因子和變量對于過程的信息解釋率都達到收斂。相較于傳統的FA(ST)指標,本文提出的NLLP統計量加入了更全面的統計信息,并通過數值例子和TE過程仿真說明了加入方差信息后的NLLP的優越性。

[1] ZHOU D H, HU Y Y. Fault diagnosis techniques for dynamic systems[J]. Acta Automatica Sinica, 2009, 35 (6): 748-758.

[2] SONG B, TAN S, SHI H B. Process monitoring via enhanced neighborhood preserving embedding [J]. Control Engineering Practice, 2016, 50: 48-56.

[3] GE Z Q, SONG Z H, GAO F R. Review of recent research on data-based process monitoring [J]. Industrial & Engineering Chemistry Research, 2013, 52 (10): 3543-3562.

[4] QIN S J. Survey on data-driven industrial process monitoring and diagnosis [J]. Annual Reviews in Control, 2012, 36 (2): 220-234.

[5] 王帆, 楊雅偉, 譚帥,等. 基于稀疏性非負矩陣分解的故障監測方法 [J]. 化工學報, 2015, 66 (5): 1798-1805.WANG F, YANG Y W, TAN S, et al. Fault detection method based on sparse non-negative matrix factorization [J]. CIESC Journal, 2015, 66(5): 1798-1805.

[6] GE Z Q, SONG Z H. Distributed PCA model for plant-wide process monitoring [J]. Industrial & Engineering Chemistry Research, 2013,52 (5): 1947-1957.

[7] JIN H D, LEE Y H, LEE G, et al. Robust recursive principal component analysis modeling for adaptive monitoring [J]. Industrial& Engineering Chemistry Research, 2006, 45 (2): 696-703.

[8] DENG X G, TIAN X M. Multiple component analysis and its application in process monitoring with prior fault data [J].IFAC-PapersOnLine, 2015, 48 (21): 1383-1388.

[9] KRUGER U, ZHOU Y Q, IRWIN G W. Improved principal component monitoring of large-scale processes [J]. Journal of Process Control, 2004, 14 (8): 879-888.

[10] WOLD S, ESBENSEN K, GELADI P. Principal component analysis[J]. Chemometrics & Intelligent Laboratory Systems, 1987, 2 (1/2/3):37-52.

[11] KIM D, LEE I B. Process monitoring based on probabilistic PCA [J].Chemometrics & Intelligent Laboratory Systems, 2003, 67 (2):109-123.

[12] TIPPING M E, BISHOP C M. Mixtures of probabilistic principal component analyzers. [J]. Neural Computation, 1997, 11 (2): 443-482.

[13] TIPPING M E, BISHOP C M. Probabilistic principal component analysis [J]. Journal of the Royal Statistical Society, 1999, 61 (3):611-622.

[14] GHAHRAMANI Z, HINTON G E. The EM algorithm for mixtures of factor analyzers [J]. University of Toronto Technical Report, 1997.

[15] GE Z Q, SONG Z H. Mixture Bayesian regularization method of PPCA for multimode process monitoring [J]. AIChE Journal, 2010,56 (11): 2838-2849.

[16] GE Z Q, LEI X, SONG Z H. A novel statistical-based monitoring approach for complex multivariate processes [J]. Industrial &Engineering Chemistry Research, 2009, 48 (10): 4892-4898.

[17] GE Z Q, SONG Z H. Maximum-likelihood mixture factor analysis model and its application for process monitoring [J]. Chemometrics &Intelligent Laboratory Systems, 2010, 102 (1): 53-61.

[18] MONTANARI A, VIROLI C. Maximum likelihood estimation of mixtures of factor analyzers [J]. Computational Statistics & Data Analysis, 2011, 55 (9): 2712-2723.

[19] ZHAO Z G, LI Q H, HUANG B, et al. Process monitoring based on factor analysis: probabilistic analysis of monitoring statistics in presence of both complete and incomplete measurements [J].Chemometrics & Intelligent Laboratory Systems, 2015, 142: 18-27.

[20] QIN S J, DUNIA R. Determining the number of principal components for best reconstruction [J]. Journal of Process Control, 2000, 10 (2/3):245-250.

[21] TAMURA M, TSUJITA S. A study on the number of principal components and sensitivity of fault detection using PCA [J]. Computers& Chemical Engineering, 2007, 31 (9): 1035-1046.

[22] VALLE S, LI W H, QIN S J. Selection of the number of principal components:? the variance of the reconstruction error criterion with a comparison to other methods [J]. Industrial & Engineering Chemistry Research, 1999, 38 (11): 653-658.

[23] VIROLI C. Choosing the number of factors in independent factoranalysis model [J]. Metodolo?ki Zvezki, 2005, 2 (2): 219-229.

[24] WANG F, TAN S, YANG Y W, et al. Hidden Markov model-based fault detection approach for a multimode process [J]. Industrial &Engineering Chemistry Research, 2016, 55 (16): 4613-4621.

[25] FEITAL T, KRUGER U, XIE L, et al. A unified statistical framework for monitoring multivariate systems with unknown source and error signals [J]. Chemometrics & Intelligent Laboratory Systems, 2010,104 (2): 223-232.

[26] CHEN Q, WYNNE R J, GOULDING P, et al. The application of principal component analysis and kernel density estimation to enhance process monitoring [J]. Control Engineering Practice, 2000,8 (5): 531-543.

[27] CHEN Q, KRUGER U, LEUNG A T. Regularised kernel density estimation for clustered process data [J]. Control Engineering Practice,2004, 12 (3): 267-274.

[28] DOWNS J J, VOGEL E F. A plant-wide industrial process control problem [J]. Computers & Chemical Engineering, 1993, 17 (3):245-255.

[29] CHIANG L H, RUSSELL E L, BRAATZ R D. Fault Detection and Diagnosis in Industrial Systems [M]. London: Springer, 2001:103-109.

[30] 趙忠蓋, 劉飛. 因子分析及其在過程監控中的應用 [J]. 化工學報,2007, 58 (4): 970-974.ZHAO Z G, LIU F. Factor analysis and its application to process monitoring [J]. Journal of Chemical Industry and Engineering (China),2007, 58 (4): 970-974.

Factor analysis process monitoring method based on probabilistic information of variables

HU Tingting, WANG Fan, SHI Hongbo
(Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of Education, East China University of Science and Technology, Shanghai 200237, China)

Factor analysis (FA), which noise factors are taken into consideration, can establish probabilistic generative model by the expectation maximum (EM) algorithm. However, traditional FA (ST) index may lead to missed fault alarms by utilizing only expectation information of variables and ignoring variance information of variables that is more representative of uncertainty. The drawback of FA (ST) index was revealed by probabilistic analysis of process variables. Another important part in the modeling process was to determine number of factors,which could preserve most useful process information in the meantime of dimension reducing. A negative log likelihood probability (NLLP) index, which integrated more comprehensive probabilistic information, was proposed to overcome dilemma of insufficient information of traditional monitoring index. For the determination of number of factors, a novel global-local method was introduced so that information explanation ratios of global factors and variables over process information reached convergence simultaneously. Numerical simulation and Tennessee Eastman (TE) process study illustrated effectiveness and superiority of the proposed method.

factor analysis; parameter estimation; process control; information explanation ratio; statistical analysis

date: 2017-01-13.

Prof. SHI Hongbo, hbshi@ecust.edu.cn

supported by the National Natural Science Foundation of China (61374140, 61673173) and Fundamental Research Funds for the Central Universities (222201714031, 222201717006).

TP 277

A

0438—1157(2017)07—2844—07

10.11949/j.issn.0438-1157.20170057

2017-01-13收到初稿,2017-03-23收到修改稿。

聯系人:侍洪波。

胡婷婷(1991—),女,碩士研究生。

國家自然科學基金項目(61374140,61673173);中央高校基本科研業務費專項資金(222201714031);中央高校基本科研業務費重點科研基地創新基金項目(222201717006)。

猜你喜歡
故障信息模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
故障一點通
訂閱信息
中華手工(2017年2期)2017-06-06 23:00:31
3D打印中的模型分割與打包
奔馳R320車ABS、ESP故障燈異常點亮
故障一點通
江淮車故障3例
展會信息
中外會展(2014年4期)2014-11-27 07:46:46
主站蜘蛛池模板: 日韩美一区二区| 伊人久久大线影院首页| 老司机午夜精品视频你懂的| 中文字幕欧美日韩高清| 欧美日韩免费在线视频| 久久中文字幕av不卡一区二区| 丝袜久久剧情精品国产| 中国国产A一级毛片| 久久黄色视频影| 青青久久91| jijzzizz老师出水喷水喷出| 精品无码一区二区三区在线视频| 免费毛片网站在线观看| 99久久精品美女高潮喷水| 中文字幕在线欧美| 国产成人精品男人的天堂下载| 91在线一9|永久视频在线| 人妻一区二区三区无码精品一区| 女人18毛片一级毛片在线 | 国产精品毛片一区视频播| 欧美日韩另类在线| 在线观看国产精品日本不卡网| 久久免费视频6| 在线国产毛片| 91小视频在线| AⅤ色综合久久天堂AV色综合| 精品一区二区无码av| 国产精品私拍在线爆乳| 综合亚洲网| 露脸国产精品自产在线播| 成人福利在线视频免费观看| 精品亚洲麻豆1区2区3区| 色综合天天娱乐综合网| 啪啪永久免费av| 91精品专区国产盗摄| 99青青青精品视频在线| 亚洲色图在线观看| 毛片网站免费在线观看| 亚洲欧洲自拍拍偷午夜色无码| AV无码一区二区三区四区| 国模沟沟一区二区三区| 极品性荡少妇一区二区色欲 | 亚洲国产黄色| 精品五夜婷香蕉国产线看观看| 超清无码一区二区三区| 宅男噜噜噜66国产在线观看| 国产免费福利网站| 欧美色综合网站| 国产成人免费视频精品一区二区| 久久永久精品免费视频| 日本不卡在线| 中文字幕无码电影| 免费一看一级毛片| 日韩精品一区二区三区swag| 一级爆乳无码av| 蜜臀AV在线播放| 中文字幕久久波多野结衣| 99久久国产综合精品女同| 福利在线一区| 国产精品区视频中文字幕 | 久久伊人色| 国产原创自拍不卡第一页| 波多野结衣久久高清免费| 91精品啪在线观看国产60岁 | 国产精品免费福利久久播放| 国产va免费精品观看| 国产网站免费| 无码免费的亚洲视频| 黄片一区二区三区| 亚洲欧美激情小说另类| 国产精品永久在线| 中文字幕亚洲综久久2021| 国产成人久视频免费| 国产精品无码久久久久AV| 日韩小视频在线观看| 重口调教一区二区视频| 精品久久久无码专区中文字幕| 日本精品αv中文字幕| 国产精品va免费视频| 久久精品无码国产一区二区三区| 国模私拍一区二区 | 午夜精品久久久久久久99热下载|