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

基于生成對抗網(wǎng)絡(luò)的離心泵時序數(shù)據(jù)異常檢測*

2023-12-20 12:39:12李思漢張鑫宇李云鵬
機電工程 2023年12期
關(guān)鍵詞:檢測方法模型

李思漢,黃 倩,付 強*,張鑫宇,李云鵬

(1.江蘇大學(xué) 流體機械工程技術(shù)研究中心,江蘇 鎮(zhèn)江 212001;2.中國核電工程有限公司,北京 100840; 3.核電泵及裝置智能診斷運維聯(lián)合實驗室,江蘇 鎮(zhèn)江 212013)

0 引 言

隨著人工智能的發(fā)展,機器學(xué)習(xí)進(jìn)一步提升了旋轉(zhuǎn)機械故障診斷的精度。機器學(xué)習(xí)五要素分別為:數(shù)據(jù)、模型、學(xué)習(xí)準(zhǔn)則、優(yōu)化算法和評價指標(biāo)[1]。其中,數(shù)據(jù)是機器學(xué)習(xí)的基礎(chǔ),是進(jìn)行模型訓(xùn)練和測試的要素。采用異常數(shù)據(jù)少和覆蓋范圍廣的數(shù)據(jù)通常能夠訓(xùn)練出良好的模型。

異常數(shù)據(jù)主要有兩種:一種是由于設(shè)備故障或設(shè)備安裝錯誤而產(chǎn)生的異常數(shù)據(jù),另一種則是不符合正常數(shù)據(jù)分布的潛在相關(guān)性的數(shù)據(jù)[3]。

在采集離心泵運行數(shù)據(jù)時,難免會采集到上述兩種情況的異常數(shù)據(jù)。因采集到的正常數(shù)據(jù)和異常數(shù)據(jù)是無標(biāo)簽數(shù)據(jù),所以筆者選用更適用于無標(biāo)簽數(shù)據(jù)的無監(jiān)督學(xué)習(xí)模型進(jìn)行異常數(shù)據(jù)檢測。

無監(jiān)督學(xué)習(xí)模型的異常檢測方法可分為以下4類:基于線性模型的方法、基于距離的方法、基于概率和密度估計的方法以及基于機器學(xué)習(xí)的方法。

許曉東等人[4]采用無監(jiān)督學(xué)習(xí)K-Means聚類方法,對網(wǎng)絡(luò)流量異常進(jìn)行了檢測,解決了異常數(shù)據(jù)檢測中準(zhǔn)確率低和誤報警率高的問題;但k值的大小對K-Means聚類方法的影響較大,且對數(shù)據(jù)量較大的數(shù)據(jù)集處理能力較差。王晨曦等人[5]提出了一種基于閾值的聚類方法;但旋轉(zhuǎn)機械數(shù)據(jù)呈多變性和突變性,故閾值的聚類方法不適用于進(jìn)行離心泵的異常數(shù)據(jù)檢測。

劉峰麟等人[6]采用了基于密度估計的空間聚類算法(density-based spatial clustering of applications with noise,DBSCAN),該算法可以根據(jù)數(shù)據(jù)特點來確定聚類半徑和k值,解決了數(shù)據(jù)的異常檢測問題和閾值選擇問題;但密度估計的方法并不能用于處理數(shù)據(jù)量較大的時序數(shù)據(jù)。唐立等人[7]基于無人機產(chǎn)生的異常軌跡數(shù)據(jù),提出了改進(jìn)孤立森林的無監(jiān)督學(xué)習(xí)算法,該算法可標(biāo)明無人機異常數(shù)據(jù)類型和異常節(jié)點;但孤立森林算法隨著數(shù)據(jù)量的增加,分割次數(shù)也隨之增多,導(dǎo)致孤立森林的參數(shù)難以調(diào)整和固定[8]。

KO J U等人[9]對某火力發(fā)電廠的異常數(shù)據(jù)進(jìn)行了研究,發(fā)現(xiàn)自編碼器(AE)可以有效地對多維數(shù)據(jù)進(jìn)行異常檢測;但自編碼器的重構(gòu)能力較差。AN J等人[10]在自編碼器的基礎(chǔ)上,提出了變分自編碼器(variational auto-encoders,VAE);但當(dāng)變分自編碼器處理未訓(xùn)練過的數(shù)據(jù)類型時,其檢測能力較差。

生成對抗網(wǎng)絡(luò)作為目前深度學(xué)習(xí)領(lǐng)域的新興模型,受到了廣泛的關(guān)注。它在生成圖像和圖像分類方面的成功應(yīng)用,已經(jīng)證明了該模型的杰出性能。此外,最近的研究發(fā)現(xiàn),生成對抗網(wǎng)絡(luò)可以處理序列形式的數(shù)據(jù),這表明了生成對抗網(wǎng)絡(luò)對時序數(shù)據(jù)進(jìn)行異常檢測方面的潛力。

綜上所述,筆者提出采用生成對抗網(wǎng)絡(luò)對離心泵時序數(shù)據(jù)進(jìn)行異常檢測的方法(該方法可以優(yōu)化生成對抗網(wǎng)絡(luò),解決梯度消失問題),隨后搭建離心泵異常數(shù)據(jù)檢測試驗臺,基于試驗數(shù)據(jù)訓(xùn)練生成器和判別器,將重構(gòu)損失和判別損失構(gòu)建成檢測閾值,完成對離心泵異常數(shù)據(jù)的檢測。

1 異常檢測模型

1.1 GAN模型

與變分自編碼器、深度信念網(wǎng)絡(luò)等算法不同,生成對抗網(wǎng)絡(luò)采用了隱式密度模型。它充分利用了神經(jīng)網(wǎng)絡(luò)的擬合能力,估算了數(shù)據(jù)分布的隱式密度函數(shù)[11-12]。

生成對抗網(wǎng)絡(luò)異常數(shù)據(jù)檢測訓(xùn)練模型與檢測模型總體流程圖如圖1所示。

在訓(xùn)練模型時,向GAN模型中的生成器輸入高斯分布白噪聲,生成器會根據(jù)輸入的白噪聲隨機生成時序數(shù)據(jù),生成器在判別器的監(jiān)督下不斷提高生成水平,使得判別器對生成器所產(chǎn)生的數(shù)據(jù)判別為“真”。同時,進(jìn)行訓(xùn)練的方式是將正常數(shù)據(jù)輸入判別器中,提高判別器的判別能力。

判別器與生成器交替迭代訓(xùn)練時,如果生成器生成數(shù)據(jù)的能力過強,則會傾向生成容易“欺騙”判別器的數(shù)據(jù)。因此,當(dāng)生成器與判別器達(dá)到納什平衡時即認(rèn)為訓(xùn)練完成。

筆者根據(jù)訓(xùn)練完成時生成器產(chǎn)生的重構(gòu)損失和判別器產(chǎn)生的判別損失來構(gòu)建閾值,利用構(gòu)建的閾值來檢測數(shù)據(jù)是否異常。

如圖1右半部分所示,筆者使用傳感器和振動采集儀,采集的信號會上傳至上機位;采集的數(shù)據(jù)輸入至訓(xùn)練好的生成器和判別器中,進(jìn)而計算得到重構(gòu)損失和判別損失;結(jié)合重構(gòu)損失和判別損失構(gòu)建異常損失值,若異常損失值超過閾值,則該數(shù)據(jù)為異常數(shù)據(jù)。

1.2 基礎(chǔ)網(wǎng)絡(luò)LSTM

LSTM是由ft∈[0,1]D、it∈[0,1]D、ot∈[0,1]D這3個門控制數(shù)據(jù)傳遞信息的。其中“0”代表“門”為關(guān)閉狀態(tài),“1”代表“門”為全開狀態(tài)。LSTM網(wǎng)絡(luò)中的“門”為“軟門”,其值可在0到1中調(diào)節(jié)。輸入門(it)調(diào)節(jié)當(dāng)前時刻記憶單元(Ct′)中的信息保存數(shù)量;遺忘門(ft)調(diào)節(jié)上一時刻內(nèi)部狀態(tài)(Ct-1)中的信息遺忘數(shù)量;輸出門(ot)調(diào)節(jié)當(dāng)前時刻內(nèi)部狀態(tài)(Ct)輸出給外部狀態(tài)(ht)的信息數(shù)量[13]。

LSTM 3個門的計算公式表示如下:

it=σ(Wixt+Uiht-1+bi)

(1)

ft=σ(Wfxt+Ufht-1+bf)

(2)

ot=σ(Woxt+Uoht-1+bo)

(3)

式中:σ(·)為Logistic函數(shù);W,U為系數(shù)組成的矩陣;b為偏執(zhí)向量;xt為當(dāng)前時刻的輸入狀態(tài);ht-1為上一時刻的外部狀態(tài)[14-16]。

首先,LSTM會結(jié)合上個時刻外部狀態(tài)(ht-1)和當(dāng)前輸入狀態(tài)(xt)來計算3個門控;其次結(jié)合遺忘門(ft)和輸入門(it)來調(diào)控記憶單元(Ct′);最后由輸出門與內(nèi)部狀態(tài)來確定外部狀態(tài)(ht)。

ht的計算公式表示如下:

ht=ot⊙tanh(Ct),

(4)

式中:⊙為向量元素乘。

1.3 基于GAN的異常檢測

離心泵數(shù)據(jù)采集為八通道同步采集,數(shù)據(jù)采集頻率范圍為8 Hz~204 kHz。數(shù)據(jù)異常檢測即檢測每個通道的數(shù)據(jù)是否處于異常狀態(tài)[17-18]。訓(xùn)練GAN模型時的訓(xùn)練數(shù)據(jù)皆為正常數(shù)據(jù),訓(xùn)練模型的目的是將正常數(shù)據(jù)標(biāo)記為0,異常數(shù)據(jù)標(biāo)記為1。

模型訓(xùn)練時,生成器和判別器會被看成一個整體。同時,整個生成對抗網(wǎng)絡(luò)的目標(biāo)函數(shù)被視為大小博弈的過程。

博弈公式表示如下:

(5)

式中:Pr為真實數(shù)據(jù)分布;P(z)為低維空間數(shù)據(jù)分布;θ,φ為生成器與判別器的參數(shù);D(x),G(z)為生成器和判別器的映射函數(shù)[19-21]。

1.4 梯度消失問題解決

當(dāng)真實數(shù)據(jù)分布(Pr)與生成數(shù)據(jù)分布(Pθ)不相交時,生成對抗網(wǎng)絡(luò)會出現(xiàn)梯度消失的問題[22]。筆者提出了使用Wasserstein距離方法來解決上述問題。

當(dāng)真實數(shù)據(jù)分布和生成數(shù)據(jù)分布不相交或相交點較少時,Wasserstein距離方法可用于對Pr與Pθ之間的距離進(jìn)行計算。Wasserstein距離方法計算公式表示如下:

W(pr,pθ)=infγ~∏(pr,pθ)E(x,y)~γ[‖x-y‖]

(6)

式中:∏(Pr,Pθ)為Pr和Pθ組合起來的聯(lián)合分布的集合;[║x-y║]為Pr和Pθ間的距離;γ為聯(lián)合分布;E(x,y)~γ[║x-y║]為Pr和Pθ間的距離在聯(lián)合分布下的期望值。

由于Wasserstein距離公式中的Infγ~∏(Pr,Pθ)無法直接求出,因此根據(jù)Kantorovich-Rubinstein對偶定理,筆者將該模型中的真實數(shù)據(jù)和生成數(shù)據(jù)分布之間的1st-Wasserstein距離轉(zhuǎn)為1-Lipschitz連續(xù)函數(shù)期望差值的上界,有效解決了Infγ~∏(Pr,Pθ)無法直接求出的問題。

為了進(jìn)一步優(yōu)化,筆者提出了評價網(wǎng)絡(luò)(f(x;φ))。評價網(wǎng)絡(luò)同時滿足K-Lipschitz連續(xù)函數(shù)的特點,使上界的計算變得更簡潔,其上界公式表示如下:

(7)

1.5 異常檢測方式

生成對抗網(wǎng)絡(luò)對時序數(shù)據(jù)進(jìn)行異常檢測時,筆者采用了雙重判別的方法。時序數(shù)據(jù)分別輸入訓(xùn)練好的生成器和判別器,判別器通過判別損失來判斷數(shù)據(jù)是否異常,生成器通過重構(gòu)損失來判斷數(shù)據(jù)是否異常。

在對異常數(shù)據(jù)進(jìn)行檢測時,測試數(shù)據(jù)和生成器的輸出樣本之間會產(chǎn)生一定的殘差,殘差的大小被作為判斷數(shù)據(jù)是否異常的依據(jù)。因受到兩時序數(shù)列起點、終點以及連續(xù)性等條件的約束,所以筆者使用適應(yīng)性較強的動態(tài)時間扭曲方法(dynamic time warping,TW)來計算兩時序數(shù)列的最小距離。其距離公式表示如下:

(8)

式中:ωk為初始兩數(shù)列的度量。

再結(jié)合判別器的判別損失,筆者將兩個判別標(biāo)準(zhǔn)結(jié)合,得到了異常檢測損失值(δ),計算公式表示如下:

δ=λδG+(1-λ)δD

(9)

式中:δG為生成器的重構(gòu)損失;δD為判別損失;λ為參數(shù)。

2 異常數(shù)據(jù)檢測試驗

2.1 試驗裝置

此處所用異常數(shù)據(jù)檢測試驗數(shù)據(jù)集由江蘇省某泵站數(shù)據(jù)集和江蘇大學(xué)流體機械工程技術(shù)研究中心試驗臺數(shù)據(jù)集組成。

臥式離心泵異常數(shù)據(jù)檢測試驗臺如圖2所示。

采集儀器詳細(xì)參數(shù)如表1所示。

圖2 臥式離心泵異常數(shù)據(jù)檢測試驗臺Fig.2 Experiment bench for abnormal data detection of horizontal centrifugalpump 1為15 kW的電機;2為扭矩儀,設(shè)備型號為JCIA;3為振動傳感器;4為力傳感器;5為試驗水泵,水泵為臥式離心泵,額定流量為100 m3/h,額定揚程為32 m,額定效率77.9%,轉(zhuǎn)速2 900 r/min;6為進(jìn)出口壓力變送器;7為故障診斷信號采集儀,八通道同步采集且最大采集率為204 kHz,采集范圍為±10 V電壓信號模擬量,其內(nèi)置恒流適配器可直接采集振動信號;8為上位機。

表1 采集設(shè)備參數(shù)表

計算機使用Windows 10 64位操作系統(tǒng),LabVIEW的版本為2020版,Python版本為3.8。由LabVIEW采集和處理數(shù)據(jù),并調(diào)用由Python編寫的生成對抗網(wǎng)絡(luò)代碼來進(jìn)行異常數(shù)據(jù)檢測。

2.2 數(shù)據(jù)預(yù)處理

泵站共有4臺立式離心泵,每2臺水泵為一組,執(zhí)行周輪換制。離心泵的參數(shù)為:揚程H=16 m,流量Q=2 700 m3/h,功率P=160 kW。此處以離心泵連續(xù)運行7天產(chǎn)生的數(shù)據(jù)作為泵站的數(shù)據(jù)集。

筆者提取泵站的正常數(shù)據(jù)的80%和試驗臺的正常數(shù)據(jù)的80%作為生成對抗網(wǎng)絡(luò)模型的訓(xùn)練集。生成對抗網(wǎng)絡(luò)模型的驗證集包括泵站正常數(shù)據(jù)的20%、試驗臺正常數(shù)據(jù)的20%以及泵站和試驗臺數(shù)據(jù)中的異常數(shù)據(jù)。

筆者將采樣頻率設(shè)置為10 000 Hz,采樣數(shù)設(shè)置為5 000。LabVIEW能夠用于顯示采集的數(shù)據(jù)并進(jìn)行傅里葉變換。

數(shù)據(jù)異常原因及數(shù)據(jù)量如表2所示。

表2 數(shù)據(jù)異常原因及數(shù)據(jù)量

數(shù)據(jù)的時域?qū)Ρ葓D如圖3所示。

圖3 數(shù)據(jù)時域?qū)Ρ葓DFig.3 Time domain comparison of data

圖3(c)中被圈出的數(shù)據(jù)是由于傳感器故障而產(chǎn)生的數(shù)據(jù),可見其幅值瞬間減小并持續(xù)到采樣結(jié)束。采集后的數(shù)據(jù)在LabVIEW程序中同時進(jìn)行時頻域轉(zhuǎn)換。

數(shù)據(jù)的頻域?qū)Ρ葓D如圖4所示。

圖4 數(shù)據(jù)頻域?qū)Ρ葓DFig.4 Comparison of data frequency domain

圖4中的1和4區(qū)域分別為正常數(shù)據(jù)50 Hz和150 Hz中的頻域幅值,特征較為明顯。因試驗數(shù)據(jù)100 Hz(二倍頻)不夠明顯突出,所以不作為異常數(shù)據(jù)的判斷依據(jù)。圖4中的2和5區(qū)域為異常數(shù)據(jù)的頻域幅值,其幅值遠(yuǎn)低于正常數(shù)據(jù)的幅值。圖4中的3區(qū)域的基頻在20 Hz,明顯不同于正常數(shù)據(jù)的基頻50 Hz。

2.3 評價指標(biāo)

為驗證生成對抗網(wǎng)絡(luò)的異常數(shù)據(jù)檢測能力,筆者將生成對抗網(wǎng)絡(luò)與孤立森林、AE、K-Means算法和一類支持向量機(one-class support vector machine,OC-SVM)進(jìn)行對比,并應(yīng)用精確率(precision)、召回率(recall)和F1值來評估生成對抗網(wǎng)絡(luò)的異常檢測效果。評價指標(biāo)公式表示如下:

(10)

(11)

(12)

式中:TP為實際數(shù)據(jù)正常,檢測結(jié)果也為正常,即正樣本被正確識別的數(shù)量;TN為實際數(shù)據(jù)異常,檢測結(jié)果也為異常,即負(fù)樣本被正確識別的數(shù)量;FP為實際數(shù)據(jù)異常,檢測結(jié)果為正常,即誤報的負(fù)樣本數(shù)量;FN為實際數(shù)據(jù)正常,檢測結(jié)果為異常,即漏報的正樣本數(shù)量;Pre,Rec為精確率與召回率的縮寫。

數(shù)據(jù)的優(yōu)劣會影響旋轉(zhuǎn)機械故障診斷模型的性能。如果正常數(shù)據(jù)被當(dāng)成異常數(shù)據(jù)處理,數(shù)據(jù)將不夠完整并丟失重要特征。因此召回率是評價異常數(shù)據(jù)檢測模型性能的重要指標(biāo)。

2.4 時序數(shù)據(jù)樣本生成

在GAN模型訓(xùn)練過程中,生成器會根據(jù)采集到的數(shù)據(jù)生成逼真的數(shù)據(jù)樣本。為了使生成器生成良好的數(shù)據(jù),筆者在時域特征中選取1個特征,在頻域特征中選取2個特征,觀察訓(xùn)練過程中生成數(shù)據(jù)的特征分布,確定最佳訓(xùn)練輪次。

試驗結(jié)果特征分布對比如圖5所示。

圖5 試驗結(jié)果特征分布對比圖Fig.5 Comparison of experimental results feature distribution

由圖5可知:通過多次試驗對比發(fā)現(xiàn),生成器訓(xùn)練50輪時可達(dá)到最優(yōu)生成結(jié)果。

2.5 異常類型分類

在采集信號的過程中,傳感器的穩(wěn)定性十分重要。當(dāng)傳感器脫落或損壞時,傳感器上傳的數(shù)據(jù)將無法正確體現(xiàn)設(shè)備狀態(tài),并產(chǎn)生大量的異常數(shù)據(jù)。因此需要對發(fā)生該情形時所產(chǎn)生的異常數(shù)據(jù)進(jìn)行報警。

筆者在試驗中發(fā)現(xiàn),當(dāng)試驗臺發(fā)生傳感器故障時,生成對抗網(wǎng)絡(luò)平均重構(gòu)損失值約為277。泵站發(fā)生傳感器故障時,平均重構(gòu)損失值約為271。為了防止出現(xiàn)傳感器發(fā)生故障時,生成對抗網(wǎng)絡(luò)未進(jìn)行報警的情況,故筆者選用最低值271作為閾值。當(dāng)重構(gòu)損失超過閾值,以報警的方式提醒監(jiān)測人員,設(shè)備可能出現(xiàn)故障。

3 結(jié)果分析

因泵站和試驗臺兩組數(shù)據(jù)集略有差異,故筆者將兩組數(shù)據(jù)集的結(jié)果取平均值并進(jìn)行了分析。對于精確率這一指標(biāo),生成對抗網(wǎng)絡(luò)約為89.5%,比AE平均值高2.5%,比孤立森林平均值高7.5%,高出K-Means平均值和OC-SVM平均值約20%;對于召回率這一指標(biāo),生成對抗網(wǎng)絡(luò)為69.5%,比孤立森林平均值高9.5%,比K-Means平均值高18%,高出AE平均值和OC-SVM平均值約30%;對于F1值這一指標(biāo),生成對抗網(wǎng)絡(luò)為78%,比孤立森林平均值高10%,比K-Means平均值高21%,同時也遠(yuǎn)高于AE和OC-SVM算法的平均值。

綜合上述3個指標(biāo)可以看出,生成對抗網(wǎng)絡(luò)優(yōu)于其他4種算法。盡管應(yīng)用在試驗臺上比應(yīng)用在泵站上的精確度略有下降,但從整體表現(xiàn)來看,生成對抗網(wǎng)絡(luò)的精確率可達(dá)到90%左右,且召回率這一指標(biāo)相比其他算法表現(xiàn)優(yōu)異,表明生成對抗網(wǎng)絡(luò)可以用于對離心泵的異常數(shù)據(jù)進(jìn)行檢測。

精確率對比結(jié)果如圖6所示。

圖6 精確率對比圖Fig.6 Accuracy comparison chart

召回率對比結(jié)果如圖7所示。

圖7 召回率對比圖Fig.7 Comparison of recall rates

F1值對比結(jié)果如圖8所示。

圖8 F1值對比圖Fig.8 Comparison of F1 scores

4 結(jié)束語

為了實現(xiàn)離心泵在線異常數(shù)據(jù)檢測的目的,筆者提出了一種基于生成對抗網(wǎng)絡(luò)的離心泵時序數(shù)據(jù)異常檢測方法(該方法可以優(yōu)化生成對抗網(wǎng)絡(luò),解決梯度消失問題)。

該方法通過在生成對抗網(wǎng)絡(luò)中引入基礎(chǔ)網(wǎng)絡(luò)LSTM,以適應(yīng)時序數(shù)據(jù)的特性,使用生成器和判別器對采集數(shù)據(jù)進(jìn)行了雙重檢測,提高了異常數(shù)據(jù)檢測的精度。

研究結(jié)論如下:

1)利用LabVIEW程序采集數(shù)據(jù)并進(jìn)行了數(shù)據(jù)的預(yù)處理,利用頻域中的倍頻信息進(jìn)行了初步的異常數(shù)據(jù)檢測;

2)采用Wasserstein距離方法解決了生成對抗網(wǎng)絡(luò)梯度消失的問題,增強了網(wǎng)絡(luò)在現(xiàn)實中的可用性。生成對抗網(wǎng)絡(luò)使用雙重判別方法,增強了異常數(shù)據(jù)檢測的能力;

3)試驗結(jié)果表明,生成對抗網(wǎng)絡(luò)在精確率、召回率和F1值3個評價指標(biāo)中的表現(xiàn)良好,能夠有效檢測出異常數(shù)據(jù)。

在后續(xù)工作中,筆者將進(jìn)一步優(yōu)化算法,由單維數(shù)據(jù)異常檢測提升為多維數(shù)據(jù)異常檢測,以縮短異常數(shù)據(jù)檢測的時間。

猜你喜歡
檢測方法模型
一半模型
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
小波變換在PCB缺陷檢測中的應(yīng)用
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 欧美日韩激情在线| 欧美日韩中文字幕二区三区| 欧美精品v| 特级做a爰片毛片免费69| 天天激情综合| 免费人成视频在线观看网站| 欧美激情视频一区二区三区免费| 国产精品v欧美| 爱做久久久久久| 人妻精品久久无码区| 国产美女在线免费观看| 国产H片无码不卡在线视频| 国产亚洲精| 亚洲手机在线| 国产自产视频一区二区三区| 国产在线一区视频| 中文国产成人精品久久| 精品久久香蕉国产线看观看gif| 日韩欧美国产区| 国产欧美精品一区二区| 福利片91| 国产成人永久免费视频| 亚洲精品动漫| 99热这里都是国产精品| 国产尤物视频网址导航| 欧美色亚洲| 国产亚洲现在一区二区中文| 2020极品精品国产| 国产精品手机视频| 亚洲欧美不卡中文字幕| 国产成人精品一区二区三区| 亚洲Aⅴ无码专区在线观看q| 国产精品三区四区| 最新国产高清在线| 亚洲A∨无码精品午夜在线观看| 思思热精品在线8| 制服丝袜一区| 凹凸精品免费精品视频| 国产乱人伦AV在线A| 精品无码日韩国产不卡av| 99久久99这里只有免费的精品| 欧美日本视频在线观看| 亚洲成在线观看| 欧美特级AAAAAA视频免费观看| 亚洲无码不卡网| 国产成人一区| 福利国产在线| 在线看片中文字幕| 五月综合色婷婷| www.亚洲色图.com| 无码人妻免费| 一区二区日韩国产精久久| 高清免费毛片| 五月天丁香婷婷综合久久| 久青草国产高清在线视频| 亚洲成人黄色在线观看| 干中文字幕| 色播五月婷婷| 在线精品亚洲一区二区古装| 亚洲日本韩在线观看| 日韩大乳视频中文字幕| 久青草免费在线视频| 国产成人亚洲无码淙合青草| 亚洲欧美精品日韩欧美| 九九九久久国产精品| 欧美乱妇高清无乱码免费| 欧美人人干| 特级毛片8级毛片免费观看| 欧美在线中文字幕| 国产jizz| 57pao国产成视频免费播放| 天天摸夜夜操| 午夜激情福利视频| 久久国产乱子| 91网址在线播放| 美女免费黄网站| 午夜欧美理论2019理论| 91在线高清视频| 五月激情婷婷综合| 999精品视频在线| 久久久精品国产亚洲AV日韩| 欧美无遮挡国产欧美另类|