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

改進(jìn)LSTM神經(jīng)網(wǎng)絡(luò)在電機(jī)故障診斷中的應(yīng)用

2021-05-07 09:41:44林谷燁
關(guān)鍵詞:故障診斷故障模型

張 凱,林谷燁,羅 權(quán)

(華南理工大學(xué) 電力學(xué)院,廣州 510000)

0 引言

伴隨著人工智能領(lǐng)域的發(fā)展,越來(lái)越多的傳統(tǒng)領(lǐng)域開始嘗試采用神經(jīng)網(wǎng)絡(luò)理論來(lái)解決工程實(shí)踐中難以精確建立數(shù)學(xué)模型的問題。異步電動(dòng)機(jī)由于具有故障種類較多,受外界干擾大等特點(diǎn),而傳統(tǒng)基于建模分析的方法難以建立適應(yīng)各個(gè)復(fù)雜場(chǎng)景下的三相異步電機(jī)運(yùn)行時(shí)的數(shù)學(xué)模型,所以使用神經(jīng)網(wǎng)絡(luò)進(jìn)行故障診斷分析相比于傳統(tǒng)的電機(jī)建模分析等方法有著得天獨(dú)厚的優(yōu)勢(shì)。目前國(guó)內(nèi)外也已經(jīng)有一些研究采用神經(jīng)網(wǎng)絡(luò)進(jìn)行電機(jī)故障診斷的研究[1-3]。但是目前采用的方法仍然存在一些不足:首先,基于模型的方法雖然精度高,但計(jì)算量大,且難以得到精確的模型;其次,對(duì)于傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)模型存在無(wú)法排除擾動(dòng),收斂精度不足或者無(wú)法保證收斂等問題。

在這樣的背景下,本文基于最新提出的RAdam優(yōu)化器構(gòu)建LSTM神經(jīng)網(wǎng)絡(luò),相比于傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)而言:1)LSTM神經(jīng)網(wǎng)絡(luò)的網(wǎng)絡(luò)結(jié)構(gòu)更適合進(jìn)行時(shí)序自回歸預(yù)測(cè)[4];2)使用RAdam優(yōu)化器解決了傳統(tǒng)神經(jīng)網(wǎng)絡(luò)模型在訓(xùn)練的初期學(xué)習(xí)率的方差較大、收斂速度較慢、收斂精度較低等問題;3)針對(duì)輸出精度較低以及模型易受干擾等問題[5-6],采用雙峰插值法[7]進(jìn)行滑窗濾波,并采用時(shí)間序列進(jìn)行輸入以提高精度。

1 雙峰譜線插值與滑窗法提取諧波

在實(shí)際工業(yè)場(chǎng)景下,由于電網(wǎng)電壓存在波動(dòng)等情況,導(dǎo)致三相異步電機(jī)在運(yùn)行過(guò)程中無(wú)法針對(duì)一個(gè)長(zhǎng)的時(shí)間序列進(jìn)行高精度的諧波分析,所以在實(shí)際應(yīng)用中采用滑窗法對(duì)實(shí)際運(yùn)行的波形進(jìn)行分段分析[8]。

1.1 雙峰譜線插值算法

假設(shè)含有H次諧波分量的信號(hào)x(t)經(jīng)過(guò)采樣頻率為f(s)的均勻采樣后得到如下形式的離散時(shí)間信號(hào):

(1)

式中,H為諧波次數(shù);當(dāng)i=1時(shí),A1、f1和θ1分別為基波的幅值、頻率和相位;當(dāng)i≠1時(shí),Ai、fi和θi分別為第i次諧波的幅值、頻率和相位。

對(duì)式(1)的信號(hào)加窗處理后的離散傅立葉變換表達(dá)式為:

(2)

式中,Δf=fs/N為采樣頻率分辨率;k為采樣序號(hào);N為數(shù)據(jù)截?cái)嚅L(zhǎng)度。

忽略其他諧波對(duì)第i次諧波的泄露影響以及負(fù)頻點(diǎn)-f1處頻峰的旁瓣影響,并有rif1=kioΔf,假設(shè)要求解的是第i次諧波分量,則有:

(3)

由于實(shí)際采樣時(shí)會(huì)出現(xiàn)采樣周期非整數(shù)次以及采樣非同步而導(dǎo)致頻譜泄露,所以ki0一般不為整數(shù)。設(shè)峰值點(diǎn)左右兩側(cè)的譜線為第ki1和ki2,則有ki1≤ki0≤ki2=ki1+1。引入?yún)?shù)α=ki0-≤ki1-0.5且α∈[-0.5,0.5],并設(shè)β=(y2-y1)/(y2+y1),并由式(3)可得:

f(α)

(4)

采用余弦組合窗函數(shù),并計(jì)算可得:

(5)

α=f-1(β)≈

a1·β+a3·β3+…+a2/l+1·β2/l+1

(6)

由α可求得第i次諧波的頻率修正公式、諧波相位修正公式分別為:

fi=ki0Δf=(α+ki1+0.5)Δf

(7)

針對(duì)t=1,2兩種情況的結(jié)果進(jìn)行加權(quán)平均最終可得到相位修正公式:

θi=θi1·(α+0.5)+θi2·(-α+0.5)

(8)

第i次諧波的幅值修正是對(duì)第ki1和ki2根譜線取加權(quán)平均,其計(jì)算公式為:

Ai=

(9)

1.2 滑窗法提取基波

在雙峰譜線插值算法的基礎(chǔ)上,針對(duì)時(shí)間跨度較大的波形采用分段滑窗求取整個(gè)時(shí)間域上精確的諧波。

圖1為利用雙峰插值和滑窗法提取基波的示意圖。

圖1 滑窗法提取基波原理圖

對(duì)于相鄰的兩個(gè)窗之間需要進(jìn)行插值運(yùn)算,設(shè)第i個(gè)插值點(diǎn)左側(cè)的基波幅值、頻率以及相角為A1、f1、θ1,插值點(diǎn)右側(cè)的基波幅值、頻率以及相角為A2、f2、θ2,且插值點(diǎn)所在時(shí)刻為ti,則根據(jù)左右兩側(cè)極限的平均值可得插值點(diǎn)的值為:

(10)

由式(7)~(10)推廣至整個(gè)采樣時(shí)間域,并設(shè)采樣起始時(shí)刻為t0,采樣終止時(shí)刻為tn窗口大小為ΔT,則經(jīng)過(guò)雙峰插值與滑窗法提取后的采樣波形的時(shí)間域函數(shù)φ(t)表達(dá)式為:

(11)

其中:

(12)

1.3 實(shí)際數(shù)據(jù)處理效果

以實(shí)驗(yàn)時(shí)A相電壓為例,采樣頻率為8 000 Hz,采樣時(shí)間為10 s,則使用滑窗法和不使用滑窗法提取并去除基波后的諧波波形與頻譜如圖2所示。

圖2 滑窗法與不使用滑窗法在時(shí)域與頻域的對(duì)比圖

從圖2可以看出,由于電網(wǎng)電壓的基頻是隨時(shí)間波動(dòng)的,所以如果不使用滑窗法從一個(gè)小窗口提取基波將會(huì)造成很大的誤差。在頻譜圖中也說(shuō)明了使用滑窗法提取并去除基波時(shí)可以很好地動(dòng)態(tài)跟蹤電網(wǎng)基頻的變化,而不使用滑窗法則會(huì)受到電網(wǎng)電壓基頻波動(dòng)的影響。

2 改進(jìn)LSTM神經(jīng)網(wǎng)絡(luò)的時(shí)序預(yù)測(cè)模型

2.1 LSTM神經(jīng)網(wǎng)絡(luò)與RAdam優(yōu)化器

LSTM神經(jīng)網(wǎng)絡(luò)是一個(gè)深度學(xué)習(xí)網(wǎng)絡(luò),被設(shè)計(jì)為學(xué)習(xí)長(zhǎng)期依賴關(guān)系。通過(guò)引入特殊類型的門來(lái)實(shí)現(xiàn)長(zhǎng)期保存信息。遺忘門可以丟棄冗余信息;輸入門可以選擇要存儲(chǔ)在內(nèi)部狀態(tài)的關(guān)鍵信息;而輸出門用于識(shí)別輸出信息。

通常,LSTM神經(jīng)網(wǎng)絡(luò)訓(xùn)練采用Adam優(yōu)化器,是由Kingma和Lei Ba兩位學(xué)者提出的結(jié)合AdaGrad和RMSProp兩種優(yōu)化算法的優(yōu)點(diǎn)的優(yōu)化器。對(duì)梯度的一階矩估計(jì)(first moment estimation,梯度的均值)和二階矩估計(jì)(second moment estimation,梯度的未中心化的方差)進(jìn)行綜合考慮,計(jì)算出更新步長(zhǎng)。

雖然Adam優(yōu)化器中加入了一階矩估計(jì)和二階矩估計(jì)可以產(chǎn)生自適應(yīng)式的學(xué)習(xí)率從而加速了收斂過(guò)程并提升了model的泛化能力。但是由于在訓(xùn)練的初期,訓(xùn)練數(shù)據(jù)較少,此時(shí)學(xué)習(xí)率的方差大,不確定性因素較多,所以會(huì)使得訓(xùn)練初期不太穩(wěn)定,進(jìn)而導(dǎo)致模型訓(xùn)練時(shí)陷入局部最優(yōu)解。常用解決Adam收斂問題的方法有:采用自適應(yīng)啟動(dòng)方法;指數(shù)加權(quán)平均的啟發(fā)式warmup算法等。

而RAdam(Rectified Adam)優(yōu)化器在其導(dǎo)數(shù)的基礎(chǔ)上,通過(guò)估算自由度來(lái)實(shí)現(xiàn)自適應(yīng)的學(xué)習(xí)率,并進(jìn)行顯式地修正。RAdam算法的輸入有:步長(zhǎng)αt;衰減率{β1,β2}用于計(jì)算移動(dòng)平均值和它的二階矩。

首先,初始化移動(dòng)量的一階矩m1和二階矩m2并估算簡(jiǎn)單移動(dòng)平均值(SMA)的最大長(zhǎng)度ρ∞=2/(1-β2)-1;

(13)

移動(dòng)偏差經(jīng)過(guò)計(jì)算之后再針對(duì)ρt判斷步長(zhǎng)的穩(wěn)定度:

如果ρt>4則進(jìn)行自適應(yīng)動(dòng)量來(lái)更新參數(shù)。移動(dòng)二階矩的修正值和方差修正值為:

(14)

此時(shí)的模型參數(shù)迭代規(guī)則如下:

(15)

如果ρt≤4,則使用非自適應(yīng)動(dòng)量更新參數(shù):

(16)

輸出為θt表示第t步后的模型參數(shù)。

2.3 LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行時(shí)序預(yù)測(cè)的原理

時(shí)序預(yù)測(cè)建模是使用現(xiàn)有的數(shù)據(jù)來(lái)建立未來(lái)數(shù)據(jù)的模型。建模的前提是所要預(yù)測(cè)的時(shí)間序列是一個(gè)穩(wěn)定的、可預(yù)測(cè)的序列,通常也稱為自回歸(AR,autoregressive)預(yù)測(cè),其表達(dá)式如下:

(17)

式中,c是常數(shù)項(xiàng),φi是建模參數(shù),εt是隨機(jī)誤差值。則基于此模型擴(kuò)展的神經(jīng)網(wǎng)絡(luò)自適應(yīng)學(xué)習(xí)算法的公式如下所示:

xt=f(xt-1,xt-2,…,xt-n)t∈(n+1,k)

(18)

本文采用LSTM神經(jīng)網(wǎng)絡(luò)建立k-n數(shù)據(jù)集的訓(xùn)練模型,即得到網(wǎng)絡(luò)f。與傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)不同,LSTM神經(jīng)網(wǎng)絡(luò)在訓(xùn)練階段不是獨(dú)立的,網(wǎng)絡(luò)不是直接從當(dāng)前時(shí)刻t開始訓(xùn)練,而是從以前的訓(xùn)練結(jié)果0到t-1繼續(xù)訓(xùn)練,充分利用了時(shí)間序列的時(shí)間連續(xù)性。LSTM網(wǎng)絡(luò)訓(xùn)練模型時(shí),以xt-1為輸入,xt為回歸數(shù)據(jù)的輸出。

最后,在進(jìn)行訓(xùn)練之前對(duì)LSTM神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行構(gòu)建。神經(jīng)網(wǎng)絡(luò)的層結(jié)構(gòu)如表1所示。整個(gè)神經(jīng)網(wǎng)絡(luò)由一個(gè)多維向量的輸入層,多個(gè)LSTM和Dropout組成的中間層和一個(gè)回歸輸出層組成。其中加入Dropout層是避免模型出現(xiàn)過(guò)擬合問題。

表1 時(shí)間序列的LSTM模型

3 利用LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行電機(jī)故障診斷

3.1 電機(jī)故障診斷的原理

(19)

其中:f是2.3小節(jié)所述的時(shí)間序列模型,ε為允許誤差范圍,Φ(t)為校驗(yàn)函數(shù)。

在實(shí)際應(yīng)用中,采用電網(wǎng)三相電壓以及運(yùn)行時(shí)三相電流作為整個(gè)LSTM網(wǎng)絡(luò)的輸入。針對(duì)電網(wǎng)中基波數(shù)值較大導(dǎo)致在訓(xùn)練時(shí)無(wú)法正確捕捉電網(wǎng)諧波,從而無(wú)法正確預(yù)測(cè)在特定電網(wǎng)諧波下電機(jī)實(shí)際運(yùn)行特性的問題[10-12],文章提出采用1.2小節(jié)中所述的方法,訓(xùn)練時(shí)采用滑窗法處理后的電網(wǎng)三相電壓的諧波以及電機(jī)實(shí)際運(yùn)行時(shí)的三相電流的基波作為輸入向量,并且以三相電流的實(shí)際波形作為輸出,并且將預(yù)測(cè)理論值的頻譜與實(shí)際采樣時(shí)的電流頻譜作對(duì)比,并利用特征電流法理論對(duì)比兩電流的實(shí)際頻譜,便可得出實(shí)際電機(jī)有無(wú)故障以及故障類型。

同時(shí),傳統(tǒng)的訓(xùn)練方法采用單個(gè)輸入?yún)?shù)進(jìn)行訓(xùn)練,魯棒性較差。而在本文中對(duì)此方法進(jìn)行了改進(jìn),采用6×n的向量形式作為層內(nèi)數(shù)據(jù),以提取數(shù)據(jù)特性,如圖3所示。

圖3 LSTM神經(jīng)網(wǎng)絡(luò)滑窗算法

3.2 故障診斷的算法流程

電機(jī)故障診斷算法分為模型訓(xùn)練和模型預(yù)測(cè)與故障診斷兩個(gè)部分。整個(gè)算法的流程如圖4所示。

圖4 故障識(shí)別算法流程圖

通過(guò)訓(xùn)練數(shù)據(jù)集,得到了一個(gè)LSTM模型。當(dāng)模型遇到正常工作數(shù)據(jù)時(shí),它能在下一時(shí)刻精確地預(yù)測(cè)參數(shù)。當(dāng)出現(xiàn)異常狀態(tài)時(shí),預(yù)測(cè)結(jié)果會(huì)偏離實(shí)際值,從而判斷故障的發(fā)生。并通過(guò)將故障頻譜與正常運(yùn)行時(shí)頻譜進(jìn)行比較,基于特征電流法可以判斷出故障類型。

4 實(shí)際算例分析

4.1 故障診斷系統(tǒng)

為了獲得待測(cè)電機(jī)分別在各種故障情況下的定子電流,需要搭建電機(jī)故障監(jiān)測(cè)實(shí)驗(yàn)平臺(tái)。

實(shí)驗(yàn)平臺(tái)的結(jié)構(gòu)如圖5(a)所示。圖中M為待測(cè)異步電機(jī),G為發(fā)電機(jī),R為水阻箱,G和R組成待測(cè)電機(jī)的負(fù)載。利用電機(jī)信號(hào)采集儀,監(jiān)測(cè)并采集定子電信號(hào),將采集到的信號(hào)上傳至上位機(jī)中,由上位機(jī)內(nèi)置的故障診斷系統(tǒng)進(jìn)行數(shù)據(jù)分析與故障識(shí)別。圖5(b)為感應(yīng)電機(jī)信息采集分析系統(tǒng)示意圖。

圖5 故障診斷系統(tǒng)實(shí)驗(yàn)示意圖

考慮到實(shí)驗(yàn)電機(jī)的選型要滿足易發(fā)生高頻故障,又考慮到價(jià)格需要在實(shí)驗(yàn)成本之內(nèi)。本次實(shí)驗(yàn)選擇一臺(tái)低壓電機(jī)作為實(shí)驗(yàn)的待測(cè)電機(jī),電機(jī)銘牌參數(shù)如表2所示。

表2 實(shí)驗(yàn)電機(jī)銘牌參數(shù)

對(duì)于信號(hào)采集單元,本文實(shí)際監(jiān)測(cè)設(shè)備接線如圖6所示。

圖6 監(jiān)測(cè)設(shè)備接線示意圖

采樣模塊通過(guò)網(wǎng)線連接至局域網(wǎng)的電腦,即可實(shí)時(shí)獲得測(cè)量的數(shù)據(jù)。采樣模塊擁有7個(gè)通道接入測(cè)量數(shù)據(jù),本文采用其三相電壓、電流作為實(shí)驗(yàn)測(cè)試。其中電壓變比kU=5 928 256*515,電流變比kI=5 928 256*10.31。采樣頻率為8 kHz,即每個(gè)周波160個(gè)點(diǎn)。采樣設(shè)備內(nèi)部基于UDP通信協(xié)議并以Modbus隧道的通信方式將電壓、電流波形數(shù)據(jù)向外發(fā)送。

本文針對(duì)異步電機(jī)運(yùn)行時(shí)最經(jīng)常出現(xiàn)的偏心故障進(jìn)行實(shí)驗(yàn),分別測(cè)取各型號(hào)電機(jī)在正常運(yùn)行以及人為制造偏心故障時(shí)的運(yùn)行數(shù)據(jù)。其中,人為制造偏心故障的方法是:打磨電機(jī)端蓋一側(cè),并在另一側(cè)塞入銅箔,使電機(jī)轉(zhuǎn)子斜偏心,偏心約10絲(1絲=0.01 mm),處理結(jié)果如圖7所示。

圖7 氣隙偏心的實(shí)現(xiàn)

經(jīng)過(guò)試驗(yàn)后,針對(duì)電流數(shù)據(jù)過(guò)濾出整數(shù)次諧波,并做傅里葉分析,可得如圖8所示的偏心故障局部特征頻譜圖。

圖8 22 kW偏心故障局部特征頻譜圖

由圖8可見,當(dāng)異步電機(jī)發(fā)生靜態(tài)偏心故障時(shí),定子電流出現(xiàn)頻率為99.5 Hz的邊頻帶。通過(guò)傳感器測(cè)速得電機(jī)發(fā)生斷條故障時(shí)的轉(zhuǎn)速n=2 970 r/min,轉(zhuǎn)差率s=0.01。根據(jù)特征電流法理論分析[13],定子電流的邊頻帶頻率應(yīng)為:

(2-s)f1=(2-0.01)*50=99.5 Hz

(20)

實(shí)驗(yàn)結(jié)果與理論分析完全一致。

4.2 與傳統(tǒng)神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)對(duì)比

為驗(yàn)證所提出的基于RAdam優(yōu)化器的LSTM算法的有效性,選擇了經(jīng)典的BP網(wǎng)絡(luò)、標(biāo)準(zhǔn)的LSTM神經(jīng)網(wǎng)絡(luò)、采用RAdam優(yōu)化器的LSTM神經(jīng)網(wǎng)絡(luò)以及經(jīng)過(guò)數(shù)據(jù)預(yù)熱的RAdam優(yōu)化器的LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行比較。其中BP網(wǎng)絡(luò)結(jié)構(gòu),以及LSTM相關(guān)網(wǎng)絡(luò)的結(jié)構(gòu)如表3所示。

表3 神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)

圖9 4種網(wǎng)絡(luò)結(jié)構(gòu)在驗(yàn)證集上的表現(xiàn)

訓(xùn)練數(shù)據(jù)長(zhǎng)度為100 s,采樣率8 000 Hz,即800 000個(gè)點(diǎn)。使用均方誤差(MSE,mean square error)來(lái)衡量訓(xùn)練準(zhǔn)確性。4種神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)在驗(yàn)證集的比較結(jié)果如圖9所示。從驗(yàn)證集的結(jié)果可以看出,BP網(wǎng)絡(luò)相較于LSTM網(wǎng)絡(luò)在頻譜預(yù)測(cè)上存在較大的誤差,并且在300 Hz與400 Hz處均存在錯(cuò)誤的預(yù)測(cè)結(jié)果。而LSTM網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果明顯更優(yōu),并且使用RAdam優(yōu)化器進(jìn)行優(yōu)化的網(wǎng)絡(luò)具有更好的測(cè)試數(shù)據(jù)集擬合能力。使用Warm-up進(jìn)行數(shù)據(jù)預(yù)熱的結(jié)果與使用RAdam優(yōu)化器的結(jié)果極為相似,但是收斂速度更快。

4.3 實(shí)際故障診斷效果

故障檢測(cè)的原理是檢測(cè)模型中的數(shù)據(jù)與實(shí)際數(shù)據(jù)之間的差距。通過(guò)4.2小節(jié)的比較,可以看出RAdam-LSTM方法的時(shí)間序列建模更好。因此,本節(jié)針對(duì)這個(gè)方法的故障檢測(cè)效果進(jìn)行實(shí)驗(yàn),實(shí)驗(yàn)預(yù)測(cè)結(jié)果如圖10所示。

圖10 基于RAdam-LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行故障識(shí)別

在圖10中,展示了基于RAdam-LSTM時(shí)間序列模型在故障條件下的實(shí)測(cè)數(shù)據(jù)和預(yù)測(cè)數(shù)據(jù)。該圖顯示了電機(jī)在發(fā)生4.1小節(jié)所述的偏心運(yùn)行情況下的故障識(shí)別。當(dāng)發(fā)生故障時(shí),RAdam-LSTM方法在預(yù)測(cè)值和測(cè)量值之間有較大的差距,使得故障得以確定。

5 結(jié)束語(yǔ)

本文創(chuàng)新性地將LSTM神經(jīng)網(wǎng)絡(luò)應(yīng)用于三相異步電機(jī)的故障檢測(cè),通過(guò)滑動(dòng)窗口剔除基波,并建立預(yù)測(cè)模型。研究結(jié)果表明,使用LSTM神經(jīng)網(wǎng)絡(luò)進(jìn)行時(shí)間序列建模時(shí),基于RAdam優(yōu)化器的LSTM神經(jīng)網(wǎng)絡(luò)的性能優(yōu)于傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)模型。此外,當(dāng)該方法預(yù)測(cè)故障實(shí)例時(shí),在LSTM網(wǎng)絡(luò)上應(yīng)用滑動(dòng)窗口比傳統(tǒng)網(wǎng)絡(luò)方法提供了更好的分類精度,可以實(shí)現(xiàn)更少的訓(xùn)練迭代和更高的精度。最后通過(guò)實(shí)際數(shù)據(jù)實(shí)驗(yàn)表明此方法具有實(shí)際可行性,為之后電機(jī)故障診斷提供了新的思路。

猜你喜歡
故障診斷故障模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
故障一點(diǎn)通
3D打印中的模型分割與打包
奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
因果圖定性分析法及其在故障診斷中的應(yīng)用
故障一點(diǎn)通
江淮車故障3例
基于LCD和排列熵的滾動(dòng)軸承故障診斷
主站蜘蛛池模板: 精品久久久久无码| 狠狠色综合久久狠狠色综合| 国产一级一级毛片永久| www.精品国产| 毛片基地视频| 亚洲三级色| 伊人成色综合网| 欧美日韩国产在线人成app| 色网在线视频| 亚洲无码A视频在线| 伊伊人成亚洲综合人网7777 | 欧美 国产 人人视频| 亚洲国产天堂久久综合226114| 人人艹人人爽| 亚洲女同欧美在线| 伊人国产无码高清视频| 国产精品私拍在线爆乳| 成人年鲁鲁在线观看视频| 婷婷伊人五月| 午夜精品福利影院| 国产在线视频导航| 波多野结衣中文字幕一区二区| 国产老女人精品免费视频| 五月婷婷中文字幕| 黄色成年视频| 老司国产精品视频91| 无码aaa视频| 国产va免费精品观看| 思思热在线视频精品| 久久香蕉国产线看精品| 最新精品久久精品| 色悠久久综合| 中文字幕自拍偷拍| 国产成人a在线观看视频| 亚洲制服丝袜第一页| 精品91视频| 亚洲欧美日韩另类在线一| 99er精品视频| 欧美国产日韩一区二区三区精品影视| 亚洲二三区| 国产精品无码制服丝袜| 91精品情国产情侣高潮对白蜜| 亚洲欧美综合精品久久成人网| 日韩精品成人在线| 精品福利网| 久久青草热| 亚洲综合激情另类专区| 夜夜拍夜夜爽| 久久永久视频| 国产精品一区在线观看你懂的| 亚洲欧美自拍视频| 中国精品久久| 国产一区亚洲一区| 视频一区视频二区中文精品| 成人第一页| 中国特黄美女一级视频| 国产欧美在线观看精品一区污| 国产欧美日韩免费| 欧美一区二区丝袜高跟鞋| 亚洲人成网站18禁动漫无码| 国产伦片中文免费观看| 亚洲综合色婷婷中文字幕| 日韩久草视频| 国产无码制服丝袜| 国产精品无码AV中文| 狼友av永久网站免费观看| 黄色网在线| 国产成人无码AV在线播放动漫 | 日本午夜精品一本在线观看| 国产精品流白浆在线观看| 黄色一及毛片| 成年看免费观看视频拍拍| 亚洲国产成人麻豆精品| 在线国产三级| 中日无码在线观看| 亚洲精品国产自在现线最新| 国产精品国产三级国产专业不| 情侣午夜国产在线一区无码| 五月六月伊人狠狠丁香网| 亚洲视频免费在线| 伊伊人成亚洲综合人网7777| 国产日本欧美亚洲精品视|