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

基于經驗耦合函數的動態運行參數異常檢測

2023-12-15 05:27:26宋柯錢唐江武彬陳勇旭鐘婷周帆
科學技術與工程 2023年33期
關鍵詞:檢測方法

宋柯, 錢唐江, 武彬, 陳勇旭, 鐘婷, 周帆*

(1.國家能源集團國能大渡河流域水電開發有限公司, 成都 610016; 2. 電子科技大學信息與軟件工程學院, 成都 610054)

在當今工業生產系統中,對核心設備的運維大多通過智能管控系統完成,其中對關鍵設備的智能異常檢測是很重要的環節,會直接影響到工業生產系統的工作效率。設備的維護不足可能會造成巨大的經濟損失,而過度維護又會浪費不必要的人力物力資源,因此需要通過智能管控系統及時識別和處理關鍵設備故障,才能避免生產受到影響而帶來經濟損失。

動態異常檢測是工業制造領域的一個重要研究領域。國內外學者以及研究人員在這一領域都進行了廣泛而深入的研究。在國外,近年來隨著深度學習技術的發展,基于深度學習的異常檢測方法逐漸成為研究熱點[1-5]。在國內,目前學者對動態異常檢測的研究主要集中于對時間序列建模的研究,主要方法包括時間序列降維、隱空間特征提取以及時序深度學習模型等[6-13]。

然而,現有的工業異常檢測方法仍存在一些不足。首先,由于工業生產過程的復雜性和多變性,工業設備的動態運行參數大多為高維復雜時序數據,現有的異常檢測算法往往不能滿足異常檢測任務實時、可靠的需求。其次,異常檢測算法的可解釋性也需要進一步提高,以便能識別并解釋異常問題,從而為工業設備的實際維護工作提供可靠依據。

因此,針對工業設備運行動態參數數據量大、數據維度較高的特點,現提出一種基于耦合函數的動態參數異常檢測方法,利用聯合分布的耦合函數簡單高效地完成對整體時序數據的建模。在具體算法設計中,采用無參數的經驗估計法來構建耦合函數,進一步提升算法的整體效率。同時,在此基礎上構建概率模型,使得模型具備一定的可解釋性。為了驗證該方法的有效性,對大渡河流域水電站排水系統的廠房滲漏排水泵動態參數進行異常檢測。對算法進行可解釋性分析,驗證方法自身的可靠性,也為排水泵類設備后續的維護提供重要的參考依據。

1 相關工作

在近年來,許多專家學者針對不同種類的異常提出了不同的異常檢測方法,這些方法的原理不盡相同[9]。為了找出適合排水泵動態參數異常檢測這一問題的方法,首先對異常的類型進行了探究,之后研究了現有方法的基本原理及其它們的優勢與不足。

1.1 異常概述

在數據挖掘和統計學領域中,異常也被稱為偏差、離群值。異常大多是由于系統故障而產生,但也有人為因素故意造成的異常,即欺詐。在數據層面,與正常值相比異常值數量更為稀少,但往往能反映出系統的漏洞,更具有研究價值。因此異常檢測是一項重要的課題,目前,在數據科學以及應用領域已經受到廣泛且深入的研究。

異常可以被大致分成三類:點異常、條件異常和集合異常[10]。點異常表示隨機發生的偏差。條件異常表示在特定的時間或空間中表現出的偏差,不同于點異常,它需要一定環境條件的約束,通常出現在時序數據中。集合異常表現為由正常點組成的集合會在一組數據中表現為異常,通常出現在交易欺詐中。

對于動態運行參數來說,它屬于時序數據,其中發生的異常可以看作是條件異常。但由于條件異常問題的處理較為復雜,通常通過數據處理,將其轉換為點異常問題進行分析。

1.2 異常檢測方法概述

1.2.1 經典異常檢測方法

經典異常檢測方法主要利用了機器學習,通過歷史數據來檢測異常,主要分為以下四類。

(1)基于距離的方法。這種方法主要通過計算數據點之間的距離來判斷異常,最常用的是K-近鄰(k-nearest neighbor,KNN)算法[11]。

(2)基于密度的方法。這類方法也是最早的異常檢測方法之一,其核心思想認為正常值通常位于高密度區域,而異常值位于低密度區域。局部異常因子算法(local outlier factor,LOF)是最早被提出的這類算法,它在KNN的基礎上加入了局部密度理論。此后,一些基于LOF的改進算法也被相繼提出,例如連接異常因子(connectivity outlier factor,COF)分析算法,局部相關積分(local correlation integral,LOCI)算法、被動式異常因子分析(influenced outlierness,INFLO)算法等。

(3)基于統計的方法。這類方法依賴于數據分布建模,通過模型檢測數據點間的關系從而發現異常。具體包括核密度估計(kernel density estimation,KDE)為代表的無參數方法以及為混合高斯模型(Gaussian mixture model,GMM)為代表的有參數方法。

(4)基于集成的方法。這類方法使用決策樹,通過組合不同的子模型,從而提高整體模型的泛化能力和健壯性。主要包括極端梯度提升異常檢測(extreme gradient boosting,XGBoost)算法[14]和孤立森林(isolation forest,iForest)算法。

以上經典的異常檢測算法雖然在工業領域的動態異常檢測上應用廣泛,但在面對高維數據時容易陷入高維災難,具體表現為計算復雜度提高,算法效率以及準確率下降。

1.2.2 基于學習的異常檢測方法

這種方法基于深度學習或主動學習,使用這些學習方法去構建不同的分類模型,以達到異常檢測的目的。近年來,在各種領域中,深度學習都取得了很高的關注度,在異常檢測領域也不例外。深度學習的方法適用于異常檢測的各類情形,尤其適用于大規模數據。深度學習方法一般分為有監督的、半監督的和無監督的方法。其中半監督的和無監督的方法受到廣泛使用,主要原因在于異常標簽數據的缺乏和數據集類別的不平衡性,這些因素會導致有監督的方法表現欠佳。在無監督方法的模型中,自動編碼器(auto-encoder)最為常見。此外,在針對時間序列類型數據的情形下,傳統的長短期記憶(long short-term memory,LSTM)模型以及循環神經網絡(regression neural network, RNN)同樣能夠發揮重要作用[10]。

基于深度學習的方法可以提取出數據集深層次特征,能夠更好地適應大規模數據。不過,這類方法通常具有較高的計算復雜度,并且大多數深度學習的方法缺乏可解釋性。

1.3 具體問題定義

水泵廠房場景如圖1所示。排水泵類設備具有用途廣、分散性強的特點,在水電站的發電水輪機和檢修系統中均有分布,排水泵出現異常會直接影響到發電系統的工作效率。通過多次對大崗山水電站的排水系統進行實地考察及數據調研,發現排水泵是否異常與其運行的各項相關參數(如:電壓、電流等)有緊密聯系。具體來說,排水泵在工作期間會因各種外界因素或自身因素出現無法正常工作或工作效率下降的情況,此時水泵的動態運行參數很可能出現了異常。如果短時間內各項動態運行參數出現多次異常,可以表明排水泵整體出現了異常。

圖1 大渡河流域水電站排水系統水泵廠房圖Fig.1 The photo of drainage system water pump plant in the hydropower station of Dadu River Basin

2 算法設計

2.1 異常檢測方法設計

通過對排水泵的實時數據采集,得到電氣屬性參數、出口壓力參數等多種動態運行參數。這些動態參數可以反映水泵運行情況,如果某一項或某幾項動態參數出現異常,那么排水泵也會有很大概率出現異常。

由此,設計了一種從實時和整體兩個角度檢測排水泵異常的方法。

實時運行狀況檢測方法:在固定的時間間隔對動態運行參數進行采樣,對每次采集的各項參數進行異常檢測。

整體運行狀況檢測方法:一次性采集排水泵一個工作周期內的動態運行參數,計算出這段周期內每項參數的均值方差等統計特征,對各項統計特征進行異常檢測。

實時檢測可以及時發現排水泵的異常,而整體檢測可以評估排水泵工作的穩定性,為后續設備維護提供參考依據。這兩個部分的檢測結果將共同作為判斷排水泵異常的依據。

具體地,實時運行狀況檢測方法每間隔20 s采集排水泵工作期間的10維動態運行參數,分別為三相電流(Ia,Ib,Ic)、三相電壓(Ua,Ub,Uc)、三相線間電壓(Uab,Ubc,Uca)和出口壓力(fo)。將這些數據作為樣本X(X∈Rn×d,d=10),使用異常檢測方法進行檢測后得到異常分數S(S∈Rn×1)。通過設定閾值,對異常分數進行二分類,分數超過閾值的樣本判定為異常。類似地,整體運行狀況檢測方法只是更換了樣本數據,采集排水泵一個工作周期內的進出水速度(si,so)、電流的均值和方差(Iam,Icm,Iavar,Icvar)、電壓的均值和方差(Uam,Ubm,Ucm,Uavar,Ubvar,Ucvar)、出口壓力的均值和方差(fom,fovar),共14維數據。

2.2 算法基本原理

2.2.1 基于聯合分布的異常檢測思想

基于統計的異常檢測使用概率分布對觀測數據進行建模,然后可以根據數據分布模型的關系將某些符合條件的觀測點標記為異常。

在一維的情況下,對于任意一個數據點xi,如果xi是異常點,它會處于其分布的概率密度函數的高密度區域之外,也就是尾端概率P(X≤xi)或P(X≥xi)的數值會非常小。那么在其累積分布函數FX中,FX(xi)和1-FX(xi)的數值就會非常小。因此,累積分布函數中FX(xi)和1-FX(xi)的數值可以作為判斷異常的重要依據。類似地,對于多維數據,也可以使用這樣的方法來判斷異常。然而在多維情況下僅僅依靠單變量的邊緣分布函數是不夠的,這忽略了變量之間的依賴關系,這種依賴關系在水泵的異常檢測中是極其重要的。例如,由于啟停水泵的操作,三相電流與出口壓力的變動往往是協同的,當兩者的數據變動差異較大時,則可能出現異常。為此,在多維數據的情況下,通過多維變量的聯合分布函數進行建模計算尾端概率。

通常而言,對數據的聯合分布進行建模是一個復雜的問題,基于耦合函數(Copula)的方法可以解決這一問題,并在此基礎上進行異常檢測。

2.2.2 基于Copula的聯合分布建模

Copula是統計學中一種用來處理隨機變量相關性問題的工具,它可以將一個多維隨機變量的聯合分布分解為它們各自的邊緣分布和一個Copula函數的形式[12]。具體地,在聯合分布建模的過程中,排水泵的動態參數數據矩陣可以寫為(X1,X2,…,Xi,…,Xd),其中Xi表示每一類動態參數的數據向量,它們的邊緣分布Fi(x)=P(Xi≤x)連續,因此可以對它們使用反變換法,可得

(1)

通過式(1)的變換,可以得到隨機向量U=[F1(X1),F2(X2),…,Fd(Xd)]。這樣數據矩陣(X1,X2,…,Xd)的Copula函數就表示為隨機向量(U1,U2,…,Ud)的聯合分布函數,即

C(u1,u2,…,ud)=P(U1≤u1,U2≤u2,…,

Ud≤ud)

(2)

根據Sklar定理,任意的聯合概率分布函數可以被寫為單變量邊緣分布函數和一個描述它們之間依賴結構的Copula函數。使用Copula函數來捕捉多個單變量之間的依賴關系,并構建聯合分布為

F(X)=C[F1(x1),F2(x2),…,Fn(xn)]

(3)

式(3)中:F(X)為聯合分布函數;C(·)為Copula函數。

在根據動態運行參數檢測排水泵異常的問題中,使用Copula函數將多個不同的動態運行參數的邊緣分布進行建模,構造一個整體的多維數據的聯合分布。聯合分布的相關性性質僅由Copula函數決定,這使得在對排水泵動態運行參數這種高維度數據建模時更加靈活高效。具體說來,對每一個維度分別建模,得到每個維度的邊緣分布,最后通過Copula函數將它們聯系在一起,得到關于多種動態運行參數的多維隨機變量的聯合分布。這個聯合分布能夠考慮到各個參數分布之間的相互關聯,構建更準確的后驗分布。

2.2.3 經驗Copula估計方法

通過以上Copula函數的性質可知,對數據的聯合分布進行建模需要獲得數據的Copula函數的具體形式,因此首先需要找到合適的Copula函數估計方法。

Copula函數有多種表達形式,其中有些具有固定的函數形式,例如橢圓Copula函數簇、阿基米德Copula函數簇等,這些有固定形式的Copula函數通常帶有一個或多個參數。還有的Copula函數沒有固定的表達形式,它們可以由觀測數據構建,這稱為經驗Copula函數[15]。

對于Copula函數的估計,基于現有的理論有參數法、半參數法和無參數法[16]。參數法和半參數法通常使用形式固定的Copula函數,通過極大似然法計算其中的參數。然而,當數據維度過高時,這兩種方法的優化過程會變得十分復雜,容易遭受高維災難。相對地,無參數法可以構建經驗Copula函數,它在高維數據的情況下表現更好,具有更強的泛用性。因此,采用基于無參數的經驗Copula函數估計方法實現異常檢測。

(4)

式(4)中:I(·)為指示函數。

(5)

2.3 算法描述

在上述基于經驗Copula的排水泵動態運行參數異常檢測算法的基礎上,根據實際問題的需要,進行了如下修改,最終的算法流程如算法1所示。

算法1基于經驗Copula的動態運行參數異常檢測算法

輸入:由排水泵動態運行參數組成的多維觀測值矩陣X∈Rn×d,經驗分布函數計算函數eCDF(·),偏度計算函數Skew(·)

輸出:排水泵的異常分數向量S∈Rn×1

(1)forj←1 toddo

(2)end for

(3)fori←1 tondo

ifbd<0

else

(4)end for

(5)returnS=(s1,s2,…,sn)

在實際應用中,數據的概率分布通常具有不對稱性。分布的主體不一定位于中央,其概率密度函數的長尾部可能在左側,也可能在右側,如圖2所示。此時,如果使用單一的尾端概率作為異常檢測的標準,就可能出現誤判的情況。因此,需要首先了解數據分布不對稱性的情況,從而選擇合適的尾端概率作為異常檢測的參考。式(6)用于計算數據分布的偏度。當偏度值為負時,分布左偏,意味著左尾部更長,此時應該參考左尾端概率進行異常檢測。同樣地,當偏度為正時,分布右偏,意味著右尾部更長,此時應該參考右尾端概率進行異常檢測。因此在異常檢測算法中,先分別計算了左尾端概率分布和右尾端概率分布,然后根據偏度的正負選擇合適的分布作為異常分數計算的參考依據。

圖2 分布函數左偏和右偏示意圖Fig.2 The schematic diagram of positive and negative skewed distribution function

(6)

在最終結果的表現上,通過計算尾端概率可以將該觀測值的異常程度量化,概率越小說明該值越可能為異常。然而這樣獲得的異常值是一個概率值,其取值區間為[0,1],取值的分布區間過于集中,容易引發數值計算上的問題。由此,采取了將概率值取負對數的方法,將異常值區間轉換到[0,+∞),并將該值作為異常分數。概率值越小,異常的可能性越大,相應地,異常分數也會越大。

3 實驗及結果分析

3.1 數據集介紹

首先采集了2018—2020年大渡河流域水電站中三臺滲漏排水泵的歷史動態運行參數。針對實時和整體運行狀況,進行了不同的數據預處理,得到了數據集(Ⅰ)和數據集(Ⅱ)。

根據2.1節的實時運行狀況檢測方法,把在三臺排水泵上采集的電壓、電流等數據處理后作為數據集(Ⅰ)。具體信息如表1所示。

表1 數據集(Ⅰ)基本信息Table 1 The basic information of dataset(Ⅰ)

根據2.1節的整體運行狀況檢測方法,把在三臺排水泵上采集的電壓、電流等數據的統計特征進行處理,得到數據集(Ⅱ)。具體信息如表2所示。

表2 數據集(Ⅱ)基本信息Table 2 The basic information of dataset(Ⅱ)

3.2 評價標準

異常檢測問題可以看作是一個二分類問題,即分為異常和正常兩類,可以使用二分類問題中的評價標準。不同的是,檢測到異常為正例,而檢測到正常為負例,通常正例數目要遠小于負例數目。在最終結果中,記錄算法判定為異常且實際為異常的數目(true positive,TP),判定為異常而實際為正常的數目(false positive,FP),判定為正常且實際為正常的數目(true negative,TN),判定為正常而實際為異常的數目(false negative,FN)。由于異常檢測問題樣本的不平衡性,需要計算精確率(precision,P)和召回率(recall,R),如式(7)和式(8)所示。為了綜合考慮精確率和召回率帶來的影響,選擇了平均精確率(average precision,AP)作為算法的評價指標。

(7)

(8)

除此之外,還計算了曲線下的面積(area under curve,AUC)值。對于分類器的接受者操作特征(receiver operating characteristics,ROC),平面的橫坐標是FP率,縱坐標是TP率。可以根據算法在測試樣本上的表現得到一個TP率和FP率點對,這樣該分類器就可以映射成ROC平面上的一個點。調整這個分類器進行分類時使用的閾值,就可以得到一條ROC曲線。AUC值取值在0.5~1,其數值越接近1表示算法表現越好。

為了衡量算法運行效率,還選擇了運行時間作為評價指標,在相同計算資源、相同實驗數據集的條件下,運行時間越短表示算法效率越高。

3.3 對比算法

選取了以下4種常用的,并且基于不同原理的異常檢測算法和本文算法進行比較。

(1)K-近鄰(KNN)算法,一種基于距離的異常檢測算法[17]。

(2)連接的異常因子分析(COF)算法[18],一種基于密度的異常檢測算法。

(3)基于角度的異常檢測(ABOD)算法[19],一種基于統計的異常檢測算法。

(4)孤立森林(iForest)算法[20],一種基于集成的異常檢測算法。

3.4 實驗結果分析

對比算法的ROC值、AP值如圖3和圖4所示,在數據集(Ⅰ)上的運行時間如表3所示。

表3 各算法在數據集(Ⅰ)上的運行時間Table 3 The running time of each algorithm on dataset (Ⅰ)

圖3 各算法在數據集(Ⅰ)上的表現Fig.3 The performance of each algorithm on dataset (Ⅰ)

圖4 各算法在數據集(Ⅱ)上的表現Fig.4 The performance of each algorithm on dataset (Ⅱ)

從實驗結果可以看出,對于在數據集(Ⅰ)上進行的實驗,基于Copula的算法總體表現最好:其中AUC值為0.987,比表現排名第二的KNN算法高1.96%;AP值為0.873,和表現最高的iForest算法相比僅低8.49%;運行時間僅為0.781 s,在所有算法中最短。

具體分析數據集(Ⅰ)的結果可以發現,從AUC和AP兩個指標來看,iForest算法和基于Copula的算法表現較佳,這是因為iForest算法可以通過控制二叉樹的數量,使算法表現變得穩定,而基于Copula的異常檢測算法能夠建模聯合分布,因此算法性能也較好。而COF算法本身用于處理低密度下的異常值,隨著數據量以及數據維度的增加COF算法表現顯著下降。KNN算法則無法找出局部異常點,同樣不能達到最優表現。從算法效率來看,由于高數據維度的影響,計算數據點的距離或密度函數的計算復雜度變大,基于距離的算法KNN和基于密度的算法COF運行時間顯著變長,基于集成的算法iForest具有線性復雜度,運行時間相對較短。因為采用了無參數的經驗估計法,基于Copula的異常檢測算法效率明顯優于其他算法。

對于在數據集(Ⅱ)上進行的實驗,與在數據集(Ⅰ)上的實驗不同的是,數據集(Ⅱ)的數據維度更高,高維度的影響使得所有算法的表現均有所下降,其中基于Copula的算法總體表現下降最少,AUC值下降了2.83%,AP下降了0.23%。這體現了基于Copula的算法對高維度數據具有更強的適應性。在數據集(Ⅱ)上并沒有比較運行時間,這是由于數據集(Ⅱ)數據量較少,各算法在運行時間上相差無幾。

3.5 可解釋性分析

在實時運行狀況檢測中,算法的最終結果是一個整體的異常分數,但仍然可以在計算過程中獲得各個維度單獨的異常分數,對于10維度的數據點xi,可以得到10個維度的異常分數si=(s1,s2,…,s10)。它可以描述在異常發生時,哪個維度的數據對整體異常的貢獻最大,這為更加精確地定位異常提供了重要的依據。由此,任意選取實驗中兩個標為正常和異常的數據點,得到異常分數對比圖。如圖5所示,對這一可解釋性結果進行描述,其中95%為數據集中正常數據的比例。

a,b,c為三相線;i為對應相的電流;u為電壓;uab,ubc,uca為相應的線電壓;of表示出水口壓力圖5 異常分數對比圖Fig.5 Anomaly score comparison diagram

可以發現,正常點46號的所有特征的異常分數全部在95%臨界線以下,而異常點45號的特征ua的異常分數已經超過99%臨界線,特征of的異常分數超過了95%臨界線。這表明在數據點45號采集時,排水泵的A相電壓和出口壓力很可能出現了異常,這些異常很可能會在后續導致排水泵整體工作異常,這使得工作人員的故障排查能夠更具有目的性。

4 結論

提出了一種基于經驗估計耦合函數的動態參數異常檢測方法。首先對大渡河流域水電站排水泵開展實地數據調研,具體分析了排水泵動態運行參數,通過對時序數據采樣,將問題轉化為點異常檢測問題進行研究。之后根據基于聯合分布的異常檢測方法,使用了經驗Copula函數這一數學工具對數據分布進行建模,完成了異常檢測任務。對比實驗顯示,算法具有最優的算法性能并較大地提升了算法效率。最后,對異常檢測結果進行可解釋性分析,證實了本文方法的可靠性,且該方法能夠為故障排查和后續維護提供重要依據。

研究仍存在數據豐富度不足的問題。在數據調研與采集的過程中,只能從水電站數據中心獲得部分類別的水泵動態運行數據,其中大部分是電氣屬性相關的數據。根據泵類設備研究經驗,對水泵運行影響較大的參數包括電機轉子振動或溫度等參數。采集這些數據并進行分析可以更加全面地檢測排水泵故障,從而提高方法整體的可靠性。針對以上的不足,將在以后的研究工作中進行進一步改進與優化。

猜你喜歡
檢測方法
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
“幾何圖形”檢測題
“角”檢測題
學習方法
小波變換在PCB缺陷檢測中的應用
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
主站蜘蛛池模板: 五月天在线网站| 91蝌蚪视频在线观看| 另类重口100页在线播放| 97se亚洲| www.99精品视频在线播放| 97se亚洲综合在线韩国专区福利| 一区二区三区在线不卡免费| 久久中文电影| 国产精品黑色丝袜的老师| 高清视频一区| 2020亚洲精品无码| 国产一区在线视频观看| 国产最爽的乱婬视频国语对白| 欧美在线免费| 嫩草影院在线观看精品视频| 免费在线a视频| 亚洲黄色视频在线观看一区| 国产精品va免费视频| 久久精品无码国产一区二区三区| 国产福利免费在线观看| 毛片免费视频| 国产精品自在自线免费观看| 国产一区成人| 华人在线亚洲欧美精品| AⅤ色综合久久天堂AV色综合| 91精品小视频| 波多野结衣一区二区三区AV| 精品丝袜美腿国产一区| 伊人色天堂| 一区二区理伦视频| 国产黄色免费看| 成人亚洲视频| 免费a级毛片视频| 精品国产一区91在线| 99热在线只有精品| 丝袜亚洲综合| jizz亚洲高清在线观看| 欧美午夜理伦三级在线观看| 成人免费黄色小视频| 尤物成AV人片在线观看| 久久精品嫩草研究院| 久久午夜夜伦鲁鲁片不卡| 秋霞午夜国产精品成人片| 中文字幕va| 免费看美女自慰的网站| 久久久久亚洲AV成人网站软件| 亚洲AV无码一二区三区在线播放| 国产精品永久久久久| 国产亚洲精品91| 国产成人8x视频一区二区| 91丨九色丨首页在线播放| 丁香综合在线| 亚洲成aⅴ人在线观看| 18禁影院亚洲专区| 91热爆在线| 一区二区三区精品视频在线观看| 99在线免费播放| 国产九九精品视频| 久久永久视频| 国产视频欧美| 亚洲人成网站18禁动漫无码| 男人天堂亚洲天堂| 精品亚洲麻豆1区2区3区 | 国产精品自在在线午夜| 五月天在线网站| 91年精品国产福利线观看久久| 久久人与动人物A级毛片| 尤物成AV人片在线观看| 久草网视频在线| 内射人妻无套中出无码| 国产男人天堂| 国产男人的天堂| 久操线在视频在线观看| 亚洲精品少妇熟女| 国产男女免费视频| 成年A级毛片| 91成人在线观看| 国产自无码视频在线观看| 亚洲无码视频一区二区三区| 99在线观看精品视频| 久久青草视频| 在线播放精品一区二区啪视频|