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

高階頻率加權能量算子在列車軸箱軸承故障診斷中的應用

2019-12-21 03:05:54劉澤潮張兵易彩吳文逸黃晨光
西安交通大學學報 2019年12期
關鍵詞:振動故障信號

劉澤潮,張兵,易彩,吳文逸,黃晨光

(西南交通大學牽引動力國家重點實驗室,610031,成都)

軸箱軸承作為軌道交通車輛走行部中的關鍵旋轉部件,不但承受著車輛全部的垂向載荷,同時遭受著強烈的輪對沖擊。惡劣的服役環境使其容易發生突發性失效,從而威脅列車的運行安全[1],因此需要對軸承的服役狀態進行實時監控,保證列車的運行安全。復雜的服役環境導致采集到的振動信號成分更為復雜,對故障沖擊的識別提出了更高的要求[2]。

根據循環平穩理論可知,軸承的故障信號具有典型的二階循環平穩特性,即信號的時域信號是非周期的,但是由故障沖擊引起的瞬時能量波動具有典型的周期特性[3]。因此,可通過識別軸承振動信號中的瞬時能量變化,實現軸承故障診斷與識別。

根據振動方程可知,在計算故障沖擊引起的能量波動時不僅應當包含瞬時幅值部分,而且應當包含瞬時頻率部分的權重。希爾伯特變換(HT)通過對故障沖擊包絡信息的計算,實現信號中沖擊信息的識別[4-5]。但HT在對信號進行解調時,計算的僅僅是信號的瞬時幅值,不包含瞬時頻率部分的權重。當信號信噪比較低且存在其他干擾頻率時,解調效果明顯下降[6]。特別是當干擾頻率的幅值較大時,在HT解調譜中的主要譜線是干擾頻率及其調制后的頻率成分[7]。為更好地提取信號中瞬時能量的波動,Teager提出了Teager能量算子(TEO)[8]。TEO通過原始信號、信號一階導數和二階導數的非線性組合,對信號中的瞬時能量波動進行了有效估計,實現了瞬態沖擊特征的提取[9-13],但是TEO的本質是對信號中瞬時能量的估計,原則上其不應當存在無意義的負能量。為避免像TEO一樣存在無意義的負值,且又同時在解調時保留信號的瞬時頻率部分的權重,O′Toole提出了一種非負的能量算子——頻率加權能量算子(FWEO)[14]。FWEO通過計算信號導數的HT,得到信號的瞬時能量估計,故而其不存在負能量部分[15]。眾多學者通過對FWEO的研究表明,FWEO具有比HT和TEO更好的抗干擾特性[16]。但是,當信號中的噪聲與干擾頻率的能量較大時,FWEO依然無法實現信號瞬時沖擊特征的識別。

FWEO的核心是通過求導加入瞬時頻率的權重,在此基礎上本文通過計算信號的高階導數,加大瞬時頻率部分的權重,提高解調的可靠性。將所提出的解調方法稱為高階頻率加權能量算子(HFWEO)。同時,提出相關峭度準則,確定合適的HFWEO階次,保證解調的可靠性,并通過仿真信號驗證了所提算法在強干擾環境下依然具有較好的魯棒性。

1 基礎理論

1.1 振動方程

根據牛頓力學方程,無阻尼單自由度線性系統自由振動方程可表示為[17]

(1)

x(t)=Acos(ωt+φ)

(2)

式中:A為振動的瞬時幅值;ω為系統的共振頻率;φ為初始相位。

系統振動的總能量E為

(3)

將式(2)代入式(3),則系統振動的總能量E可近似表示為

(4)

1.2 頻率加權能量算子

FWEO通過對信號求導,在計算信號的瞬時能量時加入瞬時頻率部分的權重,故而FWEO又叫做導數包絡能量算子[12]。

頻率加權能量算子的表達式為

(5)

式中:S[·]表示信號的復包絡;H[·]表示希爾伯特變換;φ[x(t)]為頻率加權能量算子。

根據式(2)可知

(6)

(7)

將式(6)(7)代入式(5),可得

φ[x(t)]=(-Aωsin(ωt+φ))2+

(Aωcos(ωt+φ))2=A2ω2

(8)

式(8)表明,FWEO與TEO一樣計算的是信號的瞬時能量,而非HT的瞬時幅值平方[15]。由于FWEO本質上是計算信號導數的包絡,所以沒有像TEO一樣存在無意義負能量部分。

2 高階頻率加權能量算子

2.1 基本原理

FWEO通過對信號求導,在計算瞬時能量時引入瞬時頻率部分的權重,從而使得FWEO具有同Teager能量算子相同的性質。鑒于此,本文提出通過求取信號的高階導數,加大瞬時頻率部分的權重,進而提高能量算子在干擾情況下解調的魯棒性。因此,將所提的方法稱為高階頻率加權能量算子(HFWEO)。

HFWEO的表達式為

ξ[x(t),m]=S[xm(t)]=|xm(t)+

jH[xm(t)]|2=xm(t)2+H[xm(t)]2

(9)

式中:xm(t)表示x(t)的m階導數;ξ[x(t),m]為m階頻率加權能量算子。

根據式(2),可得

xm(t)=Aωmcos(ωt+φ+mπ/2)

(10)

H[xm(t)]=Aωmsin(ωt+φ+mπ/2)

(11)

將式(10)和(11)代入式(9),可得

ξ[x(t),m]=(xm(t))2+(H[xm(t)])2=

(Aωmcos(ωt+φ+mπ/2))2+

(Aωmsin(ωt+φ+mπ/2))2=A2ω2m

(12)

根據式(12)可知,當階數m等于1時,HFWEO即為FWEO。隨著階數m的增加,瞬時頻率ω部分的權重呈指數級增長,故HFWEO可以更好地追蹤信號中瞬時能量的微弱變化,從而實現信號中故障沖擊的分離。

(13)

根據遞推公式,離散信號x(n)的m階導數xm(n)可表示為

xm(n)=xm-1(n+1)-xm-1(n)

(14)

將式(14)代入式(9),可得HFWEO的離散形式

ξ[x(n),m]=S[xm(n)]=

|xm(n)+jH[xm(n)]|2=

(xm-1(n+1)-xm-1(n))2+

H[xm-1(n+1)-xm-1(n)]2

(15)

式(2)的離散形式為

x(n)=Acos(ωn+φ)

(16)

將式(16)代入式(15),可得

ξ[x(n),m]=22mA2sin2m(ω/2)

(17)

當ω/2<π/4時sin(ω/2)≈ω/20,式(17)可近似為

(18)

此時,離散HFWEO與連續HFWEO解調結果近似相等。通過上述分析可知,離散化誤差是由sin(ω/2)≈ω/2導致,且sin2m(ω/2)與(ω/2)2m的差值會隨著階次m的增大而增大。連續的差分運算會導致計算誤差的增大,故應當選取合適的階次m。

2.2 解調特性分析

通常單個軸承的故障信號x(t)可以簡化為指數信號與正弦信號的調幅信號[3],表達式為

x(t)=Ae-β tcos(tω+φ)

(19)

式中:β為阻尼系數。

諧波干擾信號為

v(t)=Lhcosωht

(20)

式中:Lh為諧波干擾信號的幅值;ωh為諧波干擾信號的頻率。

信號r(t)=x(t)+v(t)中包含軸承的故障信號x(t)和諧波干擾信號v(t)。信號r(t)的FWEO解調結果為

(21)

式中

λ(t)=-2AβLhωhe-β tsin(-tωh+tω+φ)+

2AωLhωhe-β tcos(-tωh+tω+φ)

(22)

從式(21)可以看出,第1部分A2(β2+ω2)e-2β t為軸承故障信號的平方包絡A2e-2β t再乘以加權系數(β2+ω2),此部分對應的是軸承故障特征頻率的低頻部分。第2部分為干擾信號與軸承信號的互相調制部分,該部分位于高頻。第3部分與干擾信號相關。FWEO的信號干擾比σSIR可以通過幅值調制部分(式(21)的第1部分)和干擾部分(式(21)的第3部分)進行計算[19]

σSIR(φ(r(t)))=

(23)

式中:Tp為軸承故障沖擊間隔。

信號r(t)的二階HFWEO解調結果為

式中

二階HFWEO的σSIR為

(24)

二階HFWEO與FWEO的σSIR之比為

(25)

信號r(t)的m階HFWEO的σSIR為

(26)

(27)

通過式(27)可知,隨著HFWEO階次m的增大,HFWEO的σSIR逐漸提高。階次m越大,HFWEO對信號中的諧波干擾抑制效果越好。

通過對比不同HFWEO解調結果的信噪比σSNR,評估HFWEO階次對噪聲的影響。HFWEO解調結果的σSNR定義為

(28)

式中:x(t)為軸承的故障信號;n(t)為白噪聲。

原始信號的σSNR為5 dB,圖1為根據式(28)計算的1至10階HFWEO解調結果的σSNR,其中一階HFWEO即FWEO。從圖中可以看出,隨著階次的增大HFWEO的σSNR逐漸降低,這是由于連續的差分運算,導致解調結果中包含更多的高頻噪聲。圖2為σSNR=5 dB軸承故障仿真信號FWEO和10階HFWEO的解調結果,顯然FWEO結果中可以發現明顯的故障沖擊,但是10階HFWEO結果中的故障沖擊很微弱。

圖1 不同階次HFWEO解調結果σSNR

(a)FWEO結果(σSNR=6.22 dB)

(b)10階HFWE結果(σSNR=0.463 7 dB)圖2 兩種方法解調結果的對比

綜上可知,隨著階次的增大,HFWEO對諧波干擾信號的抑制效果越來越好,但是隨著連續的差分運算,導致HFWEO解調結果中引入了更多的高頻噪聲。故應當設置合適的階次m,保證HFWEO在對諧波干擾具有較好的抑制效果的同時又避免引入更多的噪聲。

2.3 階次確定

通過上述分析可知,階次對HFWEO的解調結果有著至關重要的影響,故在此提出相關峭度準則(CK)KCK來確定合適的階次m。

KCK指標不僅可以對序列中的沖擊特性進行表征,同時通過平移周期T,可以對序列中的周期性進行表征[20]。序列yn的KCK表達式為

(29)

式中:H為平移周期的個數;T為平移周期長度;N為序列長度。

在此,yn代表不同階次下的HFWEO解調結果ξ(yn,m)。通過式(29)可以得到不同階次HFWEO解調結果的KCK值,則最優的HFWEO階次mbest為

mbest=arg maxm(KCK(ξ(yn,m),T,H))

(30)

3 仿真驗證

使用仿真信號對所提算法進行驗證,仿真信號的表達式為

u(t-mTp-τm)

(31)

式中:Am是沖擊的幅值,Am=1;L是沖擊個數;Tp是故障沖擊的間隔,特征頻率fc=1/Tp=104.5 Hz;ωr是故障引起共振的頻率,ωr=4 000 Hz;τm為隨機滑移系數,通常取0.01Tp~0.02Tp;u(t)是單位階躍函數;β是阻尼系數,β=1 500 N·s/m。

信號的采樣頻率fs=12 000 Hz,信號時間長度為2 s。在信號中加入不同信號噪聲比的白噪聲與不同信號干擾比的正弦干擾信號[17],以研究所提方法在不同噪聲與干擾情況下的解調效果。

首先,對比幾種不同的解調方法它們在高σSNR與高σSIR情況下的解調效果。在信號中加入σSNR=0 dB的白噪聲和σSIR=-5 dB的正弦干擾信號,干擾信號的頻率分別為50、270、2 043和2 810 Hz。加噪后信號的時域波形與頻譜如圖3a和3b所示。雖然從仿真信號的頻譜在圖3b中可以發現很明顯的干擾頻率及噪聲,但由于信號的σSNR與σSIR都較高,依然可以很明顯地觀察到共振頻帶。

(a)時域信號

(b)傅里葉譜圖3 σSNR=0 dB和σSIR=-5 dB時仿真信號

(a)包絡譜

(b)Teager譜

(c)FWEO譜圖4 σSNR=0 dB和σSIR=-5 dB時3種解調方法結果的對比

使用HT、TEO、FWEO與HFWEO 4種方法對仿真信號進行解調,HT、TEO、FWEO解調結果如圖4所示??梢园l現,在高σSNR與σSIR情況下,HT、TEO、FWEO 3種方法的解調譜中都可以發現故障特征頻率fc及其倍頻,解調譜中存在干擾頻率。這說明3種解調方法對干擾信號都較敏感,只是由于干擾頻率的能量較低,所以從解調譜中依然可以發現故障特征頻率及其倍頻。

圖5 HFWEO相關峭度

圖6 σSNR=0 dB和σSIR=-5 dB時的仿真信號HFWEO譜

圖5為1至10階的HFWEO相關峭度,第2階HFWEO的相關峭度值最大,當階次繼續增大時,由于連續差分運算引入了更多的噪聲,所以相關峭度值逐漸降低。圖6為歸一化后1至10階的HFWEO解調譜。從不同階次的解調譜中都可以識別到故障特征頻率fc及其倍頻,同時可以發現解調譜中不存在干擾頻率。但隨著階次的增大,解調譜中的背景噪聲也逐漸增大,故障特征頻率的倍頻數量減少。綜上所述,HFWEO通過增加瞬時頻率部分的權重,保證在解調時對干擾頻率具有更好的抑制效果。同時,通過相關峭度指標可以確定合適的HFWEO階次,確保HFWEO在對諧波干擾進行抑制的同時引入更少的高頻噪聲。

保持信號的σSNR=0 dB不變,將σSIR降低為-20 dB,對比4種方法對高σSNR、低σSIR信號的解調效果。加噪信號的時域波形與頻譜如圖7a和7b所示。由于信號的σSIR較低,此時諧波頻率的幅值較大,故從信號的頻譜中很難識別出共振頻帶。

(a)時域信號

(b)傅里葉譜圖7 σSNR=0 dB和σSIR=-20 dB時的仿真信號

(a)包絡譜

(b)Teager譜

(c)FWEO譜圖8 σSNR=0 dB和σSIR=-20 dB時的3種解調方法結果對比

HT、TEO和FWEO的解調譜如圖8所示。由于信號的σSIR較低,所以從這3種方法的解調譜中都只能識別出干擾頻率,而無法識別出故障特征頻率fc及其倍頻。使用HFWEO對仿真信號進行解調,1至10階HFWEO相關峭度如圖9所示,歸一化后1至10階的HFWEO譜如圖10所示。在前4階HFWEO解調譜中只能識別到干擾頻率,此時相關峭度值也較低;到第5階HFWEO解調譜中發現故障特征頻率fc及其倍頻,但是依然可以發現干擾頻率,此時相關峭度值依然較小。隨著階次的增大,在7階以后HFWEO解調譜中只能識別到故障特征頻率fc及其倍頻,而不再存在干擾頻率。此時,相關峭度值較大,且在第9階時HFWEO的相關峭度值達到最大。同時,隨著階次的增大,高頻噪聲也在逐漸增大,第10階HFWEO的相關峭度值反而降低。這說明,HT、TEO與FWEO的抗諧波干擾性較差,而本文所提的HFWEO通過連續的導數運算,加大了瞬時頻率部分的權重,使得其在較強干擾時仍然具有較好的魯棒性。同時,相關峭度準則可以準確地得到合適的HFWEO階次,保證HFWEO的可靠性。

圖9 HFWEO相關峭度

圖10 σSNR=0 dB和σSIR=-20 dB時的仿真信號HFWEO譜

(a)時域信號

為進一步驗證4種解調算法對低σSNR和低σSIR信號的解調效果,在仿真信號中加入σSNR=-5 dB的白噪聲與σSIR=-20 dB的干擾信號。加噪仿真信號的時域波形與頻譜如圖11a與11b所示。由于信號的σSNR與σSIR都很低,所以從頻譜中同樣很難識別出共振頻帶。

(b)傅里葉譜圖11 σSNR=-5 dB和σSIR=-20 dB時的仿真信號

分別使用HT、TEO和FWEO對加噪后仿真信號進行解調,解調譜如12所示??梢钥闯?解調譜中的主要頻率成分均是干擾頻率,無法識別出故障特征頻率fc及其倍頻。

(a)包絡譜

(b)Teager譜

(c)FWEO譜圖12 σSNR=-5 dB和σSIR=-20 dB時的3種解調方法結果對比

圖13和圖14所示分別為1至10階HFWEO的相關峭度值和歸一化后的解調譜,顯然當階次較低時,HFWEO同樣無法對信號中的瞬態沖擊進行提取,解調譜中只能發現干擾頻率,故此時相關峭度較低。隨著階次的增加,干擾頻率的幅值逐漸減小,而故障特征頻率fc及其倍頻的幅值逐漸增加。當階次大于7以后,HFWEO解調譜中不存在干擾頻率,但隨著階次的增加,解調譜中的噪聲也在增大,從而影響對故障特征頻率的觀測效果,所以第7階時的HFWEO相關峭度值最大。

圖13 HFWEO相關峭度

圖14 σSNR=-5 dB和σSIR=-20 dB時的仿真信號HFWEO譜

綜上所述,本文所提的HFWEO能量算子對干擾具有很好的抑制作用,可以在低σSIR情況下具有良好的解調效果,同時所提的相關峭度準則可以準確地確定最佳的HFWEO階次。

4 試驗驗證

4.1 鐵路貨車輪對跑合試驗臺

為驗證HFWEO對實際振動信號的解調效果,本文對貨車輪對跑合試驗臺的振動數據進行分析,關于試驗臺的詳細介紹見文獻[18]。測試軸承型號為197726,在實際運行過程中,軸承形成內圈損傷,內圈滾道出現明顯剝落。輪對的旋轉速度為465 r/min,軸承的旋轉頻率為fr=7.75 Hz,采樣頻率為fs=12.8 kHz。根據軸承參數,計算得到內圈故障頻率fBPFI=88.25 Hz?,F對1 s的振動數據進行分析,振動信號的時域波形及其頻譜如圖15a和15b所示。當軸承存在內圈損傷時,振動信號的解調譜中應當出現轉頻fr、內圈故障頻率fBPFI及其倍頻和以內圈故障頻率fBPFI為中心且以轉頻fr為間隔的邊頻帶[3]。

使用HT對信號進行解調,包絡譜如16a所示。包絡譜中故障特征頻率fBPFI的譜線并不明顯,并且沒有發現以fBPFI為中心、以轉頻fr為間隔的邊頻帶。然后,使用TEO與FWEO對信號進行解調,解調譜如圖16b和16c所示。同包絡譜相同,從解調譜中同樣只能隱約觀測到fBPFI,但依然無法識別到以fBPFI為中心且以轉頻為fr間隔的邊頻帶。這說明HT、TEO和FWEO都無法對信號中的軸承故障沖擊進行提取,因此無法從解調譜中發現與內圈故障對應的頻譜特征。

(a)時域信號

(b)傅里葉譜圖15 鐵路貨車輪對跑合試驗臺振動信號及頻譜

(a)包絡譜

(b)Teager譜

(c)FWEO譜圖16 3種解調方法結果鐵路貨車輪對跑合試驗臺振動數據解調對比

使用本文所提的HFWEO對信號進行解調,并計算各階次HFWEO相關峭度值,其結果如圖17a所示。第7階HFWEO的相關峭度值最大,其解調譜如17b所示。從解調譜中不但可以識別到轉頻fr、內圈故障特征頻率fBPFI,而且還可以識別到以fBPFI為中心且以轉頻fr為間隔的邊頻帶。這與內圈故障對應的頻譜特征相吻合,說明HFWEO可以很好地從信號中分離出故障沖擊信息。

(a)HFWEO相關峭度

(b)第6階HFWEO譜圖17 鐵路貨車輪對跑合試驗臺振動數據HFWEO解調結果

在貨車輪對跑合試驗臺的振動信號中噪聲與干擾較大,軸承故障沖擊的能量較微弱,因此傳統解調方法無法對故障沖擊進行提取。根據式(27)可知,隨著階數m的增加,HFWEO的抗干擾性越來越好,故從高階HFWEO譜中可以識別到故障特征頻率。同時,相關峭度準則可以準確地確定HFWEO的階次。

4.2 高速列車輪對跑合試驗臺

為進一步驗證本文所提算法在復雜環境下的有效性,使用高速列車輪對跑合試驗臺的振動數據對算法進行驗證。圖18所示為高速列車輪對跑合試驗臺,試驗臺包括電機、驅動輪對、高速動車組輪對、高速動車組軸箱和載荷加載裝置。電機帶動驅動輪對驅動高速動車組輪對,進而帶動軸箱內的軸承旋轉,載荷加載裝置可以對軸箱施加垂向載荷,振動傳感器位于軸箱上。

圖18 高速列車輪對跑合試驗臺

在軸承的外滾道表面以120°夾角設置3處人工損傷,3處損傷的深度為1 mm,寬度分別為1、3和5 mm,損傷部位如圖19所示。

圖19 軸承外圈損傷示意圖

試驗時,在軸箱垂向上施加50 kN的垂向載荷,采樣頻率為10 kHz。驅動輪對的轉速為100 km/h,則軸承的旋轉頻率fr=10.28 Hz,因此根據軸承的幾何參數和轉頻可以計算得到外圈故障特征頻率fBPFO=83.29 Hz。

圖20a和20b所示為1 s的高速列車輪對跑合試驗臺振動信號及其頻譜。顯然,由于試驗臺的結構復雜,因而信號中的頻率成分更加豐富?,F場惡劣的采集環境導致信號中存在50、100 Hz的諧波干擾頻率。故高速列車輪對跑合試驗臺的振動信號比貨車輪對跑合試驗臺的振動信號更加復雜,對算法的魯棒性要求更高。

(a)時域信號

(b)傅里葉譜圖20 高速列車輪對跑合試驗臺振動信號及頻譜

使用HT、TEO和FWEO對信號進行解調,解調譜如圖21所示。由于諧波干擾和噪聲的影響,從它們的解調譜中都無法識別到故障特征頻率。在包絡譜中可以發現很明顯的干擾頻率,這說明HT對諧波干擾較敏感。TEO和FWEO計算的是信號的瞬時能量,對諧波干擾的抑制效果優于HT。但是,振動信號中還存在輪對沖擊的成分,且其能量大于軸承故障沖擊能量,因此從TEO與FWEO解調譜使用本文所提HFWEO對振動信號進行解調,1至10階HFWEO相關峭度如22a所示,第10階相關峭度最大,其對應的HFWEO譜如22b所示。從HFWEO中可以發現很明顯的故障特征頻率fBPFO及其二倍頻,說明本文算法成功診斷出軸承外圈故障。這是因為,輪對沖擊的能量雖然較大,但是軸承故障沖擊的瞬時頻率比輪對沖擊的瞬時頻率更高[1],HFWEO通過不斷加大瞬時頻率ω的權重,從而對瞬時頻率較高的沖擊提取效果更好,從HFWEO譜中可以更明顯識別出軸承外圈故障特征頻率fBPFO及其倍頻。

(a)包絡譜

(b)Teager譜

中僅可以識別到轉頻fr及其二倍頻。

(c)FWEO譜圖21 高速列車車輪對跑合試驗臺振動數據4種解調方法結果對比

(a)HFWEO相關峭度

(b)第10階HFWEO譜圖22 高速列車車輪對跑合試驗臺振動數據HFWEO解調結果

5 結 論

通過增加瞬時頻率部分的權重,提高了HFWEO抗干擾的能力,保證在干擾均較大時,依然具有良好的解調效果。但是,連續的差分運算導致解調時引入了高頻噪聲,因此需要確定合適的階數。通過相關峭度準則可以確定合適的HFWEO階次,保證算法在對干擾信號具有較好魯棒性的同時引入更少的高頻噪聲。仿真信號與試驗信號分析表明,本文算法可以在信號中存在嚴重干擾情況下依然可以實現軸承故障沖擊的解調,從解調譜中識別故障特征頻率。

猜你喜歡
振動故障信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
故障一點通
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
故障一點通
主站蜘蛛池模板: 亚洲欧美综合在线观看| 九九热精品视频在线| 久久国产精品电影| 欧美曰批视频免费播放免费| Jizz国产色系免费| 美女扒开下面流白浆在线试听| 女人毛片a级大学毛片免费| 99久久这里只精品麻豆| 国产91特黄特色A级毛片| 国产精品综合久久久| 原味小视频在线www国产| 五月丁香在线视频| 亚洲欧洲天堂色AV| 国产一级视频在线观看网站| AV网站中文| 99国产精品一区二区| 久久精品免费看一| 最新国产网站| 麻豆国产精品视频| 婷婷伊人久久| 91精品国产综合久久不国产大片| 免费Aⅴ片在线观看蜜芽Tⅴ| 又大又硬又爽免费视频| 欧美性精品| 亚洲精品老司机| 尤物午夜福利视频| 国产又大又粗又猛又爽的视频| 国产精品香蕉| 波多野结衣中文字幕一区二区| 国产中文在线亚洲精品官网| 亚洲a免费| 国产91视频免费| 欧美视频免费一区二区三区| 久久这里只精品国产99热8| 亚洲精品日产精品乱码不卡| 亚洲人成亚洲精品| 精品人妻系列无码专区久久| 91高清在线视频| 精品少妇人妻无码久久| 99热亚洲精品6码| 无码一区二区波多野结衣播放搜索| 激情国产精品一区| 精品福利网| 久久不卡国产精品无码| 亚洲天堂免费| 91福利片| 国产h视频在线观看视频| 精品国产成人高清在线| 亚洲AV无码不卡无码 | 国产成人精品午夜视频'| 91久久国产综合精品| 熟妇无码人妻| 国产在线视频导航| 91探花在线观看国产最新| 青青青草国产| 在线观看无码a∨| 中文字幕亚洲精品2页| 无码一区18禁| 日韩午夜伦| 玖玖精品视频在线观看| 不卡色老大久久综合网| 日韩国产综合精选| 国产肉感大码AV无码| 九九热视频精品在线| 亚洲精品天堂自在久久77| 中文天堂在线视频| 91国内在线观看| 国产成人久视频免费| 久久精品娱乐亚洲领先| 国产精品免费电影| 一级一级特黄女人精品毛片| 欧美第九页| 伊人成色综合网| 激情午夜婷婷| 91在线日韩在线播放| 五月综合色婷婷| 久久综合一个色综合网| 国产网站一区二区三区| 亚洲欧美h| 国产真实乱子伦视频播放| 先锋资源久久| 无码中文AⅤ在线观看|