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

多分辨奇異值分解在滾動軸承振動信號解調分析中的應用

2019-03-12 07:49:29羅潔思張紹輝李葉妮
振動工程學報 2019年6期
關鍵詞:故障診斷

羅潔思 張紹輝 李葉妮

摘要:針對滾動軸承在自身諧振干擾及強背景噪聲影響下,滾動軸承損傷時引起調制現象難以檢測的問題,提出基于多分辨奇異值分解(Multi-resolution Singular Value Decomposition,MRSVD)的包絡解調方法。該方法首先采用MRSVD方法將振動信號逐層分解獲得具有不同分辨率的近似信號和細節信號,經理論分析得到的第一個細節信號主要成分為噪聲,且最后一個近似信號主要成分為諧波干擾。進一步結合峭度指標從其他細節信號(第一個細節信號除外)中提取其中隱藏的周期性沖擊信號,根據周期性沖擊信號的包絡解調譜進行軸承故障的診斷。仿真分析和應用實例證明了該方法的有效性。

關鍵詞:故障診斷;滾動軸承;多分辨奇異值分解;信噪比;解調分析

中圖分類號:TH165+.3;TH13 3.3;TN911.7

文獻標志碼:A

文章編號:1004-4523 (2019) 06-1114-07

DOI:10. 16 385/j. cnki. issn. 1004-4523. 2019. 06. 021

引言

當滾動軸承的內圈、外圈或滾動體有損傷時,隨著軸承的周期性旋轉,損傷表面與其他元件表面在接觸過程中會發生周期性脈動沖擊,激起內、外圈的固有頻率振動[1-2],其振動信號中往往出現周期性的瞬態沖擊信號,形成調制現象,頻譜上表現為固有頻率兩側出現等間隔的調制邊頻帶[1,3]。因此,對軸承故障振動信號中的周期性沖擊成分進行提取和解調,根據解調譜的強度和頻次判斷軸承損傷程度和部位是軸承故障診斷廣泛使用的一種方法[3-4]。

各種信號處理方法已被用于從軸承故障振動信號中提取故障特征信息,如EMD[5]、小波分析(WT)[6]、及它們的各種改進算法[7-10]、稀疏分解[11-12]、流行學習[13]、譜峭度[14]、包絡分析[15-19]等。上述方法從時域、頻域、或時一頻聯合域檢測由軸承故障引起的周期性沖擊特征表現出的調制現象。然而,在故障診斷的初期階段,故障特征較微弱,加之環境噪聲和系統固有諧波振動的影響,使得共振帶內的周期性沖擊振動并不明顯,上述方法的有效性也大打折扣。因此,有些學者研究采用信號預處理的方法從軸承故障振動信號中去除噪聲和諧波振動干擾,提取共振帶內與故障有關的周期性沖擊特征[19-23]。目前,大多采用基于頻率的信號預處理方法,即試圖將噪聲、諧波振動、周期性的沖擊振動分解到不同的頻帶。事實上,噪聲分布頻帶較廣,共振頻帶也存在噪聲干擾,諧波振動與周期性沖擊振動也并不一定是頻域可分的。

鑒于噪聲、諧波振動、周期性沖擊振動在奇異值分解后有不同的奇異值特征,本文將新近提出的一種信號處理方法——多分辨率奇異值分解方法(MRSVD) [24]用于軸承故障振動信號的預處理,從奇異值角度將噪聲與諧波振動干擾從軸承故障振動信號中去除,可解決上述預處理方法遇到的頻率不可分問題。本文第二部分的理論分析表明:若對包含噪聲、諧波振動、周期性沖擊振動的軸承故障振動信號進行多分辨率奇異值分解,分解結果的第一個細節信號主要成分為噪聲,最后一個近似信號的主要成分為諧波振動。

因此,本文采用MRSVD方法對滾動軸承振動信號進行預處理,擬去除其中的背景噪聲信號和諧波干擾信號,并結合峭度( Kurtosis)指標從其他(除第一個細節信號外的)細節信號中提取振動信號所包含與故障有關的的周期性沖擊成分,最后根據周期性沖擊成分的Hilbert解調譜進行滾動軸承的故障診斷。仿真分析和應用實例表明該方法能在強背景噪聲及諧波振動干擾的影響下有效提取軸承故障引起的周期性沖擊成分,凸顯故障特征。

1 多分辨奇異值分解

MRSVD是在奇異值分解的基礎上,結合矩陣二分遞推結構原理,參考小波多分辨率分析思想,將信號分解到不同層次子空間的一種信號分解方法。該算法的具體步驟如下:

(1)對待分析信號A0=(x1,x2,…,xN),構造行數為 2 的 Hankel 矩 陣 H0 =

(2)對矩陣H0進行奇異值分解,可以得到且只能得到兩個奇異值,分別記為σ1,σ2,其中σ1>σ2。從σ1,σ2分別可得到兩個信號分量,記為A1,D1。這兩個分量信號A1,D1對原始信號的貢獻量是有輕重之分的,其中A1較D1對原信號的貢獻量大,是信號的主要成分,它反映了此次分解時從原始信號獲取的主要概貌,類似于小波分析中的近似信號,稱其為SVD近似信號,D1則反映此次分解時從原信號獲取的細枝末節,這類似于小波分析中的細節信號,稱其為SVD細節信號。

保留細節信號D1,對近似信號A1繼續構造行數為2的Hankel矩陣Hi,進行下一層次的奇異值分解,如此逐層遞推就可以得到一系列SVD細節信號和近似信號。

信號的MRSVD分解過程如圖1所示。

2 MRSVD在滾動軸承故障診斷中的應用研究

2.1 MRSVD的消噪與抗諧波干擾原理分析

滾動軸承振動信號通常包含三種信號分量,即由故障引起的周期性沖擊信號,由設備自身振動產生的諧振信號以及背景噪聲信號。假設采集得到的滾動軸承振動信號為x(i)。x(i)可以表達為

x(i)=S(i)+n(i)+h(i),i=1,2,…,N(1)式中 s(i)為周期性沖擊信號,n(i)為背景噪聲信號,h(i)為諧振信號,N為信號長度。

根據式(1),信號x(i)構造的Hankerl矩陣H可以表示為

H=Hs+Hn+Hh

(2)式中 Hs,Hn,Hh分別為信號S(i),n(i),h(i)構造的Hankerl矩陣。在MRSVD中,Hankerl矩陣的行數為2,其特點為:下一行矢量比上一行矢量僅僅滯后一個數據點。下面分別對這三類信號的Han-kerl矩陣奇異值分解進行分析:

(1)對于諧振信號h(i)而言,它所構造的Han-kerl矩陣兩行將高度相關,矩陣的秩為1,若Hh奇異值分解后得到兩個奇異值為σ=(Hs)=(σs,ξ),則有ξ為一很小的正數,且σs》ξ。即諧振信號h(i)的能量主要被分配至第一個奇異值中,只有很微小的一部分被分配到第二個奇異值中。

(2)對于噪聲信號n(i)而言,其自相關函數為Rn(τ)=d2δ(τ),其中d為n(i)的標準差,δ(T)為單位脈沖函數。對于噪聲信號n(i)所構造的Hankerl矩陣而言,盡管相鄰兩個行矢量同樣只滯后一位,但是它們根本不會相關,矩陣的秩為2,它的兩個奇異值接近相等,可表示為σ(Hn)=(σn,σn)。

(3)對周期性沖擊信號S(i)而言,它所構造的Hankerl矩陣兩行具有一定的相關性,但相關程度次于諧振信號h(i),記Hh奇異值分解后得到的兩個奇異值為σ(Hh)=(σhb,σhs),其中σhb為兩個奇異值中的大者,σhs為兩個奇異值中的小者。

綜合以上分析,對信號x(i)所構造的奇異值矩陣H,其奇異值矢量滿足下式

σ(H)≤σ(Hs)+σ(Hn)+σ(Hh)=

(σs+σn+ σhb,ξ十σn+σhs)

(3)

實際上,由于奇異值矩陣H的行數很小,近似有σ(H)≈(σs+σn+σhb,ξ十σn+σhs)。近似信號和細節信號的能量與其奇異值的平方呈正比,從式(3)可知,MRSVD第一次分解獲得的細節信號Di對應的奇異值為ξ十σn+σhs,其成分主要為噪聲信號和小部分的周期性沖擊信號,且D1包含的噪聲信號的總能量占x(i)所包含的噪聲總能量的一半。近似信號A1對應的奇異值為σs+σn+σhb,與信號x(i) 一樣,近似信號A1包含了噪聲成分、諧振成分和周期性沖擊成分,對其做進一步的分解,得到下一層次的細節信號D2,D2同樣包含了噪聲信號(D2包含的噪聲信號的總能量占Ai所包含的噪聲總能量的一半)和小部分的周期性沖擊信號,逐層分解得到的下一層細節信號的噪聲能量將比上一層細節信號的噪聲能量減少一半,即Dj較Dj-1具有更高的信噪比,其周期性的沖擊成分將慢慢凸顯。由于每次分解諧振信號都只有非常微小的能量(ξ2/(σs十ξ)2→0)被分解到細節信號中,因此對諧振信號h(i)而言,它的大部分能量將被保留在最后一層近似信號中。

綜合上述分析,對強背景噪聲影響下的滾動軸承信號而言,其分解得到的第一個細節信號D1主要為噪聲信號,其峭度系數大約為3,隨著分解層次的深入,細節信號所包含的噪聲成分逐層減小,周期性沖擊成分逐層增加,因此峭度系數會隨著層數的增加而增大,當峭度系數達到某個最大值后,又將隨著分解層數的增加而減少,這是因為殘余的近似信號中包含的周期性沖擊成分也在逐層減小,這將在仿真分析中得到證實。

本文結合峭度指標,從MRSVD分解得到的細節信號中選取峭度值最大的幾個進行和運算,則和信號的成分主要是與故障相關的周期性沖擊振動信號,若對該和信號進行Hilbert解調,即可從解調譜上對軸承故障進行診斷。

2.2MRSVD的分解層數的確定

由2.1節的分析可知:MRSVD的分解層數太少,得到的細節信號將包含有較多的噪聲干擾。然而分解層數也不是越多越好,分解層數多計算量也將增大,與滾動軸承故障診斷的實時性要求不符;且分解層數過多,分解結果將出現偽信號。因此本文中MRSVD的分解層數在分解過程中自適應獲取,當細節信號的峭度值達到最大,即當前近似信號若繼續分解得到的下一層細節信號的峭度值出現下降趨勢時,分解層數則確定為當前分解層數加5即可。

3 仿真分析

為了驗證基于MRSVD的包絡解調方法的有效性,下式所示的數學模型將被用于模擬滾動軸承故障振動信號[11]。式中 n(t)為噪聲信號,采用Matlab函數awgn.m進行添加;h(t)為諧波干擾信號,其包含N個諧波信號,且第j個諧波信號的幅值和頻率分別為Bj,fj。由故障引起的周期性沖擊信號共包含2M+1個沖擊,式(4)中Am為第r個沖擊的幅值,β為由結構阻尼引起的衰減系數,Tp為與故障特征頻率相對應的時間周期,即故障特征頻率fe= l/Tp,u(t)為單位沖擊函數,ωr為激勵頻率,τi為第i個均值為零,方差介于(0. Ol Tp,0.02Tp)的隨機變量,此變量是由滾子隨機滑移引起的。表1和2列出了式(4)中信號s(t),h(t)所涉及參數的取值。

設置采樣頻率為20000 Hz,采樣點數為20000,對其進行數字化采樣,添加噪聲使得信噪比為-6dB,得到離散信號的時域波形及其Hilbert包絡解調譜分別如圖2(a),(b)所示。由于噪聲和諧振干擾的影響從圖2(a)無法判斷滾動軸承故障,圖2(b)所示的Hilbert解調譜將不包含故障信息的諧波信號以其各諧振分量的頻率之差作為解調信號而解出,而實際故障特征頻率在圖2(b)中并未凸顯。

對圖2所示的仿真信號及其包含的周期性沖擊分量、噪聲分量、諧振分量分別構造行數為2的Hankerl矩陣,記為H,Hs,Hn,Hh,分別對H,Hs,Hn,Hh進行奇異值分解,各自得到的兩個奇異值如表3所示。表3所示的數值結果與本文2.1節中對3類信號(周期性沖擊信號、噪聲信號、諧振信號)的奇異值分析是相吻合的,即諧波信號的分解得到的兩個奇異值相差較大,噪聲信號分解得到的兩個奇異值數值大小接近相等,周期性沖擊沖擊信號分解得到的奇異值存在大小之分,但兩者差值既不會太大,也不會接近于零。

對圖2所示的仿真信號進行MRSVD分解,設置分解層數為15,得到的各層細節信號對應的奇異值及其峭度系數分別如圖3(a),(b)所示。由圖3(a)可見,前幾個細節信號的奇異值下降速度較快,這是因為開始時噪聲的去除比較迅速,之后隨著層數的增加細節信號的奇異值下降變得越來越平緩,這表明噪聲的去除隨著分解層數的增加而變得緩慢。如圖3(b)所示,第一個細節信號由于包含了大量的噪聲,其峭度系數接近等于3。之后隨著分解層數的增加,峭度值也在增加,這是因為細節信號所包含的噪聲成分逐層減小,周期性沖擊成分逐層增加。在分解層數為5時,細節信號的峭度值達到最大,之后峭度值又將隨著分解層數的增加而減少,這是因為殘余的近似信號中包含的周期性沖擊成分在逐層減小。因此,本文2.1節中對細節信號的峭度值分析結論在圖3(b)中得到驗證。

圖4為仿真信號MRSVD分解得到的第一層細節信號的時域、頻域波形圖,圖4信號主要成分為噪聲信號。因此,頻譜圖上看不出任何凸顯的頻率。圖5為仿真信號MRSVD分解得到的最后一層近似信號的時域、頻域波形圖。正如本文2.1節所言“諧振信號的大部分能量將被保留在最后一層近似信號中”,因此圖5(b)中凸顯的頻率為諧波頻率fi,f2,f3。

從仿真信號MRSVD分解得到細節信號中選取峭度值最大的5個進行和運算,其和信號的時域波形及其Hilbert包絡解調譜如圖6(a),(b)所示。圖6(a)所示信號主要成分為與故障有關的周期性沖擊振動信號,對其進行Hilbert解調,可凸顯故障特征頻率,因此從圖6(b)可明顯地看到故障特征頻率及其倍頻。

結合2.2節及圖3(b),MRSVD分解層數只需設置為1 0即可。

4 實測信號分析

為了驗證MRSVD在實際應用中的有效性,從如圖7所示的實驗臺上采集外圈故障滾動軸承振動信號,其中軸承的相關參數如表4所示,表4中BPFO (ballpass frequency,outer race)和fr分別表示外圈故障特征頻率、轉軸旋轉頻率。由于實際軸承振動信號通常混有噪聲和諧波干擾,為了模擬更真實的工況,采用圖7所示轉盤來引入諧振成分,即利用轉盤的不平衡引起頻率為轉軸旋轉頻率fr的諧振信號。實驗測試時,斷開傳感器信號調節器的連接,使得測試信號未經過濾波與放大,含有更多的噪聲成分,并且測試信號還包含了頻率為電力線頻率(60 Hz)的諧波干擾信號。

設置采樣頻率為20 kHz,采用點數為2 0000,轉軸旋轉速度為1440 r/min,則有fr=24 Hz,BPFO一85. 728 Hz,采集得到的含有噪聲及諧振干擾的振動信號的時域波形及Hilbert包絡解調譜分別如圖8(a),(b)所示。對圖8(b)縮小其坐標軸的顯示區間,得到圖9所示的細節圖。從圖9中能清楚地看到由于轉子不平衡引起的諧振頻率24 Hz及電力線頻率60 Hz,但卻無法看到軸承外圈故障特征頻率。

對圖8(a)所示的實測信號進行MRSVD分解,取其分解結果中峭度值最大的5個細節信號求和信號,得到如圖10 (a)所示的去噪聲和諧波干擾后的周期性沖擊振動,其Hibert包絡解調譜如圖10(b)所示,進一步縮小圖10(b)中頻譜的顯示區間,得到圖11所示的頻譜圖。從圖11所示的頻譜圖上可以明顯地看到軸承外圈故障特征頻率及其倍頻,進而判斷軸承發生了外圈故障。

5 結 論

本文提出了一種基于MRSVD的包絡解調方法,并將該方法應用于強背景噪聲及諧波干擾影響下的滾動軸承振動信號的解調。仿真分析和應用實例均表明:該方法可從含有噪聲和諧振成分的滾動軸承振動信號中提取與故障有關的周期性沖擊振動,進而從解調譜上凸顯故障特征頻率。

所提方法存在以下優點:(1)該方法根據噪聲、諧波振動、周期性沖擊振動在奇異值分解后表現出不同的奇異值特征去除滾動軸承振動信號中混有的噪聲及諧波干擾,從而提取與故障有關的周期性沖擊振動信號。因此,該方法與傳統的基于頻帶濾波的共振解調方法不同,無需預先知道故障引起的共振頻帶信息,并可解決噪聲、諧波振動、周期性沖擊振動的頻率不可分離問題;(2)該方法計算效率高,除了計算行數為2的矩陣奇異值分解外,該方法僅涉及矩陣的加、減運算。

參考文獻:

[1] Randall R B,Antoni J. Rolling element bearing diag-nostics-A tutorial[J]. Mechanical Systems&SignalProcessing, 2011, 25(2):485-520.

[2] Jaskaran S,Ashish D,Satinder P S.Rolling elementbearing fault diagnosis based on over-complete rationaldilation wavelet transform and auto-correlation of ana-lytic energy operator[J]. Mechanical Systems &Sig-nal Processing, 2018, 100:662-693.

[3] Yaqub M F,Gondal I,Kamruzzaman J.Inchoate faultdetection framework: Adaptive selection of waveletnodes and cumulant orders[J]. IEEE Transactions onInstrumentation and Measurement, 2012,61(3):685-695.

[4] 張 進,馮志鵬,褚福磊,滾動軸承故障特征的時間——小波能量譜提取方法[J].機械工程學報,2011, 47(17):44-49.

Zhang Jin, Feng Zhipeng, Chu Fulei. Extraction ofrolling bearing fault feature based on time-wavelet en-ergy spectrum[J]. Journal of Mechanical Engineering,2011, 47(17):44-49.

[5] Wang Yanxue, He Zhengjia, Zi Yanyang.A compara-tive study on the local mean decomposition and empiri-cal mode decomposition and their applications to rota-ting machinery health diagnosis[J]. Journal of Vibra-tion and Acoustics, 2010,132(2): 613-624.

[6] Chen Ruqiang, Robert X Gao, Chen Xuefeng. Wave-lets for fault diagnosis of rotary machines: A reviewwith applications[J]. Signal Processing, 2014, 96(5):1-15.

[7]Wang Hongchao, Chen Jin, Dong Guangming. Fea-ture extraction of rolling bearing's early weak faultbased on EEMD and tunable Q-factor wavelet trans-form[J]. Mechanical Systems&Signal Processing,2014, 48 (1-2):103-119.

[8]Deng Linfeng, Zhao Rongzhen. Fault feature extrac-tion of a rotor system based on local mean decomposi-tion and Teager energy kurtosis[J]. Journal of Me-chanical Science&Technology, 2014, 28 (4):1161-1169.

[9] Yi Cancan, Lu Yong, Xiao Han, et al.Matching syn-chrosqueezing wavelet transform and application toaeroengine vibration monitoring[J]. IEEE Transac-tions on Instrumentation & Measurement, 2017, 66(2) :360-372.

[10] Li Chuan, Liang Ming. Time-frequency signal analysisfor gearbox fault diagnosis using a generalized syn-chrosqueezing transform[J]. Mechanical Systems&Signal Processing, 2012, 26(1) :205-217.

[11] Shazali Osman, Wilson Q Wang. A leakage-free reso-nance sparse decomposition technique for bearing faultdetection in gearboxes[J]. Measurement Science andTechnology ,2017 ,29 (3) :1-12.

[12] Cui Lingli, Gong Xiangyang, Zhang Jianyu, et al.Double-dictionary matching pursuit for fault extent e-valuation of rolling bearing based on the Lempel-Zivcomplexity[J]. Journal of Sound& Vibration, 2016,385 :372-388.

[13] Ding Xiaoxi, He Qingbo. Time-frequency manifoldsparse reconstruction: A novel method for bearingfault feature extraction[J]. Mechanical Systems& Sig-nal Processing, 2016 ,80 :392-413.

[14] Wang Yanxue, Xiang Jiawei, Richard Markert, et al.Spectral kurtosis for fault detection, diagnosis andprognostics of rotating machines: A review with appli-cations[J]. Mechanical Systems& Signal Processing,2016 ,66-67 :679-698.

[15] Liang Ming, Bozchalooi I S. An energy operator ap-proach to joint application of amplitude and frequency-demodulations for bearing fault detection[J]. Mechan-ical Systems and Signal Processing, 2010, 24 (5):1473-14 9 4.

[16] Rodriguez P H, Alonso J B, Ferrer M A. Applicationof the Teager-Kaiser energy operator in bearing faultdiagnosis[J]. Isa Transactions, 2013,52(2) :278-284.

[17] Wang Hui, Gao Guige, Zeng Xianwen. Wind turbinefearbox gault diagnosis based on wavelet theory andHilbert demodulation spectrum[J]. Applied Mechan-ics& Materials, 2015,703:390-393.

[18] Liu Xiaofeng, Lin Bo, Peng Chang. Application of or-der cyclostationary demodulation to damage detectionin a direct-driven wind turbine bearing[J]. Measure-ment Science & Technology, 2014, 25 (2): 25004-25011.

[19] Wang Chun, Jiang Deping, He Zhengchun. Researchon fault diagnoses of rolling bearing based on the ener-gy operator demodulation approach[J]. Advanced Ma-terials Research,2012 ,588-589 :134-140.

[20] Sheen Y T. On the study of applying Morlet waveletto the Hilbert transform for the envelope detection ofbearing vibrations[J]. Mechanical Systems& SignalProcessing, 2009 ,23 (5) :1518-1527.

[21] Deng Lifeng, Zhao Rongzhen. Fault feature extractionof a rotor system based on local mean decompositionand Teager energy kurtosis[J]. Journal of MechanicalScience& Technology, 2014,28(4) :1161-1169.

[22] Han Guodong, Wan Shuting, Lu Zhanjie. The analy-SiS of gearbox fault diagnosis research based on theEMD and Hilbert Envelope Demodulation[J]. Ad-vanced Materials Research, 2014,926-930:1800-1805.

[23] Luo Jiesi, Yu Dejie, Liang Ming. A kurtosis-guided a-daptive demodulation technique for bearing fault detec-tion based on tunable-Q wavelet transform[J]. Meas-urement Science& Technology, 2013,24(5) :1-11.

[24] Zhao Xuezhi, Ye Bangyan, Chen Tongjian. Theory ofmulti-resolution singular value decomposition and its ap-plication to signal processing and fault diagnosis[J]. Jour-nal of Mechanical Engineering,2010,46 (20) : 64-75.

猜你喜歡
故障診斷
基于包絡解調原理的低轉速滾動軸承故障診斷
一重技術(2021年5期)2022-01-18 05:42:10
ILWT-EEMD數據處理的ELM滾動軸承故障診斷
水泵技術(2021年3期)2021-08-14 02:09:20
凍干機常見故障診斷與維修
基于EWT-SVDP的旋轉機械故障診斷
數控機床電氣系統的故障診斷與維修
電子制作(2018年10期)2018-08-04 03:24:46
基于改進的G-SVS LMS 與冗余提升小波的滾動軸承故障診斷
因果圖定性分析法及其在故障診斷中的應用
改進的奇異值分解在軸承故障診斷中的應用
基于LCD和排列熵的滾動軸承故障診斷
基于KPCA和PSOSVM的異步電機故障診斷
主站蜘蛛池模板: 91精品国产91久久久久久三级| 欧美一级高清免费a| 三上悠亚一区二区| 亚洲精品va| 伊人蕉久影院| 国产在线精彩视频二区| 欧美成人精品一区二区| 成人年鲁鲁在线观看视频| 日韩色图在线观看| 色噜噜狠狠色综合网图区| 无码国内精品人妻少妇蜜桃视频| 色哟哟国产精品| 精品一区二区三区水蜜桃| 91青青视频| a级毛片在线免费观看| 2020精品极品国产色在线观看| 久久视精品| 试看120秒男女啪啪免费| 国产精品亚洲五月天高清| 1024你懂的国产精品| 99爱在线| 无码有码中文字幕| 免费黄色国产视频| 亚洲综合国产一区二区三区| 免费在线a视频| 91美女视频在线观看| 亚州AV秘 一区二区三区| 青青草国产免费国产| 国产99精品久久| 日韩国产精品无码一区二区三区| a色毛片免费视频| 亚洲永久色| 亚洲黄网在线| 国产91麻豆免费观看| 日本欧美中文字幕精品亚洲| 日韩东京热无码人妻| 伊人色综合久久天天| 久久国产免费观看| 在线国产欧美| 国产特级毛片aaaaaa| 九九九久久国产精品| 免费无码网站| 在线观看网站国产| 色呦呦手机在线精品| 精品视频一区二区观看| 丰满人妻中出白浆| 丁香婷婷久久| 91蝌蚪视频在线观看| 99热这里只有精品在线观看| 在线观看国产网址你懂的| 黑人巨大精品欧美一区二区区| 日本伊人色综合网| 综合色88| 香蕉综合在线视频91| 欧美午夜在线视频| 91娇喘视频| 国产成人精品在线1区| 国产成人a毛片在线| 免费国产小视频在线观看| 国产在线观看人成激情视频| 久久综合一个色综合网| 青青草一区二区免费精品| 亚洲一区二区三区麻豆| 深夜福利视频一区二区| 国产一区二区影院| 99在线视频免费| 久久综合亚洲鲁鲁九月天| 熟妇无码人妻| 精品久久久久久中文字幕女| 欧美成a人片在线观看| 四虎永久免费在线| 91精品专区国产盗摄| 激情在线网| 欧美在线中文字幕| 美女无遮挡被啪啪到高潮免费| 亚洲国产高清精品线久久| 免费高清毛片| 日韩高清成人| 亚洲综合二区| P尤物久久99国产综合精品| 在线观看无码av免费不卡网站| 青草精品视频|