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

基于信源估計和頻域反卷積的滾動軸承故障特征分離與辨識

2017-02-10 07:19:36張建宇胥永剛北京工業(yè)大學先進制造技術(shù)北京市重點實驗室北京004臥龍電氣章丘海爾電機有限公司章丘5000北京市精密測控技術(shù)與儀器工程技術(shù)研究中心北京004
中國機械工程 2017年1期
關(guān)鍵詞:故障信號

張建宇 孟 浩 胥永剛. 北京工業(yè)大學先進制造技術(shù)北京市重點實驗室,北京,004.臥龍電氣章丘海爾電機有限公司,章丘,5000. 北京市精密測控技術(shù)與儀器工程技術(shù)研究中心,北京,004

基于信源估計和頻域反卷積的滾動軸承故障特征分離與辨識

張建宇1孟 浩2胥永剛3
1. 北京工業(yè)大學先進制造技術(shù)北京市重點實驗室,北京,1001242.臥龍電氣章丘海爾電機有限公司,章丘,2502003. 北京市精密測控技術(shù)與儀器工程技術(shù)研究中心,北京,100124

針對軸承故障的振動特征由于受到強振源的抑制作用而增加了故障分離與辨識難度的問題,建立了基于信源估計和頻域反卷積的故障診斷方法。利用小波包分解將信號分離成多個子帶信號,并和奇異值分解相結(jié)合,解決欠定條件下的信號源數(shù)估計問題;根據(jù)估計的源數(shù),選取相應(yīng)維數(shù)的觀測信號,通過短時傅里葉變換、復(fù)數(shù)域獨立分量分析、相關(guān)排序、短時傅里葉逆變換,完成頻域反卷積的分析過程,實現(xiàn)故障特征的分離與提取。仿真信號和實驗數(shù)據(jù)均驗證了該方法在故障特征分離與微弱特征辨識中的有效性。

小波包分解;奇異值分解;短時傅里葉變換;復(fù)數(shù)域獨立分量分析;頻域反卷積

0 引言

滾動軸承是機械系統(tǒng)中最重要的零部件之一,其運行狀態(tài)既是監(jiān)測診斷的重點,也是一大難點。軸承作為支撐部件,其結(jié)構(gòu)振動受到系統(tǒng)內(nèi)部多個振源的影響,如何實現(xiàn)故障源的有效解耦和準確辨識是一個非常棘手的問題。

盲信號處理可根據(jù)若干觀測信號恢復(fù)出無法直接觀測的各個源信號,這種優(yōu)勢有助于多源的故障分離與識別。獨立分量分析(independent component analysis,ICA)是其中的典型代表,作為快速算法,F(xiàn)ast-ICA自從被HYVRINEN等[1]學者提出,已經(jīng)在旋轉(zhuǎn)機械[2]、感應(yīng)電機[3]以及滾動軸承[4]的故障診斷中獲得了應(yīng)用。

然而,ICA只適用于滿足瞬時混合模型的機械振動信號,其要求傳感器的安裝位置和振源位置非常接近,以避免信號在傳遞中產(chǎn)生時延。實際上,從設(shè)備內(nèi)部振源(包括故障源)到傳感器的傳遞路徑復(fù)雜多變,信號中存在明顯的時間延遲,觀測信號多表現(xiàn)為故障源與傳輸通道的卷積混合。針對卷積混合模型,TORKKOLA[5]提出了基于信息最大化的反卷積方法,而SERVIERE等[6]和PELED等[7]學者在旋轉(zhuǎn)機械和滾動軸承的故障診斷中,均研究了反卷積的應(yīng)用效果。劉婷婷等[8]在軸承轉(zhuǎn)子的裂紋試驗中,采用反卷積算法提取故障特征。葉紅仙等[9]和潘楠等[10]分別討論了盲反卷積的時域和頻域算法。

從目前來看,盲反卷積的時域和頻域算法均是基于觀測信號數(shù)目不小于源信號數(shù)目的假設(shè),但在實際的監(jiān)測診斷中,信號源數(shù)目大多未知,如果觀測信號少于源信號,則無法完成分離或者只能分離出部分信息。因此,信源數(shù)量的正確估計是特征辨識的必要步驟,目前已有學者對振動信號的源數(shù)估計問題進行了系統(tǒng)研究[11-12]。本文把源數(shù)估計方法和頻域反卷積結(jié)合,以實現(xiàn)多源故障的特征分離與微弱故障的準確辨識。

1 基于小波包分解的源數(shù)估計方法

特征值方法是一種經(jīng)典的源數(shù)估計方法。該算法首先對多維觀測信號的相關(guān)矩陣進行奇異值分解,得到不同的特征值,然后根據(jù)特征值的分布規(guī)律來確定源信號的數(shù)目。但是,這種方法只適用于信號源數(shù)目不大于觀測信號數(shù)目的情況。由于機械設(shè)備中振動源數(shù)量多且未知,而有時在狹窄的設(shè)備空間中安裝的傳感器數(shù)量有限,可能導(dǎo)致觀測信號數(shù)量小于振動源數(shù)。小波包分解可以將單路信號分解為不同頻段內(nèi)的子帶信號,和小波分解相比,不僅可以分解低頻子帶信號,而且可以對高頻子帶信號作進一步劃分,精細化程度更高。因此,本文將通過小波包分解將單路機械振動信號拓展為多路信號,以解決觀測信號小于信號源數(shù)的問題。該方法只需要一路觀測信號就可以完成信號的源數(shù)估計,因此具有較強的適用性。本文源數(shù)估計的步驟如下:

(1)提取一路觀測信號x(t),對該信號進行小波包分解,分解層數(shù)為l,從而得到2l個小波包系數(shù)。

(2)對每個小波包的節(jié)點系數(shù)進行單支重構(gòu),組成一個2l維的矩陣xn。

(3)計算矩陣xn的協(xié)方差矩陣Rxx=E(xnxnT)

(5)通常占優(yōu)特征值(有用信號)與非占優(yōu)特征值(噪聲)有較大差別,因此兩者比值K較大。本文通過相鄰特征值的比值分布曲線確定占優(yōu)特征值的個數(shù),提取曲線中極值點所對應(yīng)的比值

K=max(λ1/λ2,…,λl/λl+1,…,λτ-1/λτ)

(1)

K對應(yīng)的分子即為最小占優(yōu)特征值,該分子的下標即為占優(yōu)特征值的個數(shù),即信號源數(shù)。

2 基于STFT的頻域反卷積方法

2.1 卷積混合的基本模型

頻域反卷積(frequency-domain blind deconvolution,F(xiàn)DBD)算法是針對符合卷積混合模型的信號而提出的一種恢復(fù)源信號的算法,卷積混合模型表達式為

(2)

其中,s(t)為源信號;X(t)為觀測信號;A(p)是時間延遲為p時的混合濾波器矩陣,表達式為

(3)

式中,每一個元素Aij都表示一個時延濾波器。

卷積混合過程可以用圖1表示。

圖1 卷積模型示意圖Fig.1 Schematic diagram of convolution model

卷積混合就是源信號si(t)和存在不同時延的傳遞路徑Aij的卷積,由于時域上卷積作傅里葉變換后會轉(zhuǎn)變?yōu)轭l域上的線性乘積,由此可實現(xiàn)信號從時域卷積混合到頻域瞬時混合的過程

x(t)*h(t)?X(f)·H(f)

(4)

2.2 時域卷積信號的頻域線性化

機械故障信號可以看作一系列短時平穩(wěn)信號的疊加。由于傅里葉變換更適合于平穩(wěn)信號,因此本文首先將非平穩(wěn)的故障信號分割為多個短時平穩(wěn)信號,再利用傅里葉變換將長時卷積混合信號轉(zhuǎn)變?yōu)轭l域多段瞬時混合信號。短時傅里葉變換公式為

(5)

通過移動窗口?(τ-t)將原始信號x(τ)分成多段,對每一段信號做傅里葉變換,從而得到該信號的時頻信息。

2.3 復(fù)數(shù)域Fast-ICA算法

STFT處理后的多段頻域信號滿足線性混合關(guān)系,因此適用于瞬時盲分離算法。考慮到STFT處理后的復(fù)信號特征,本文將引入HYVRINEN提出的復(fù)數(shù)Fast-ICA算法[13]。該方法和實數(shù)域Fast-ICA類似,都是選取負熵最大化作為表征信號之間獨立性的目標函數(shù),通過牛頓梯度法優(yōu)化解混矩陣w,從而使目標函數(shù)達到極大值。

復(fù)數(shù)域負熵可表示為

(6)

復(fù)數(shù)信息熵表達式為

(7)

其中,p表示概率密度函數(shù)。負熵Ng的計算需要已知概率密度函數(shù),這在工程應(yīng)用中顯然是不實際的。HYVRINEN給出負熵的近似公式:

Ng(w)≈E([G(wx)]2)

(8)

其中,E為數(shù)學期望。G(·)的三種近似函數(shù)及其一階導(dǎo)數(shù)分別為

(9)

分離矩陣w的求解迭代算法為

w+=E(x(wx)g(|wx|2))-E(g(|wx|2)+|wx|2g′(|wx|2))w

(10)

此外還需要對分離矩陣進行歸一化處理以解決輸出向量幅值不確定問題,每一步分離矩陣迭代公式為

(11)

2.4 源信號的復(fù)原過程

復(fù)數(shù)Fast-ICA存在輸出信號次序不確定問題,相鄰頻率段輸出信號次序不一致會引起連接混亂,進而導(dǎo)致分離失敗。由于屬于同一源信號頻率段相關(guān)性能良好,因此可以利用相鄰頻率段的相關(guān)性來調(diào)整信號次序。兩個復(fù)值信號之間的相關(guān)系數(shù)可表示為

ρ(xi(w,k),xj(w,v))=

(12)

式中,w表示時間信息的幀序號;k、v表示頻段序列號;cov(·)表示協(xié)方差。

頻域反卷積的基本流程如圖2所示。

圖2 頻域反卷積基本流程圖Fig.2 Basic procedure of deconvolution in frequency domain

具體算法流程如下:

(1)對每個輸入信號xi(t)進行短時傅里葉變換得到xi(w,k),k=1,2,…,L。

(2)將具有相同頻序號k的頻段做復(fù)數(shù)域ICA處理,共得到L組輸出信號。

(3)固定第一組輸出,然后通過復(fù)數(shù)域相關(guān)系數(shù)法調(diào)整后面每組輸出的次序,使得屬于同一源信號的輸出信號在同一個位置。

3 仿真信號算法驗證

為了驗證以上方法的有效性,構(gòu)造三組仿真信號進行卷積混合,函數(shù)形式如下所示:

(13)

仿真信號的采樣點數(shù)均為7680點,采樣頻率為15 360 Hz。三組信號均符合軸承沖擊缺陷的理論模型,δ(t)表示沖擊函數(shù),沖擊間隔分別為0.0078 s、0.0156 s、0.0039 s。因此,各組信號模擬的軸承故障特征頻率分別為127.5 Hz、63.75 Hz和255 Hz。源信號的時域波形和包絡(luò)頻譜如圖3、圖4所示。

(a)仿真信號x1

(b)仿真信號x2

(c)仿真信號x3圖3 源信號時域波形Fig.3 Waveform of source signal

(a)仿真信號x1

(b)仿真信號x2

(c)仿真信號x3圖4 源信號包絡(luò)解調(diào)譜Fig.4 Envelop demodulation spectrum of source signals

構(gòu)造如下式所示的混合矩陣:

(14)

其中,Aij為5階隨機濾波器。將式(13)三組信號進行卷積混合,混合后觀測信號的包絡(luò)頻譜如圖5所示。

從圖5可以看出,每一路混合信號中都能找到各種源信號的特征頻率,但由于各組源信號的特征頻率之間存在倍數(shù)關(guān)系,所以無法根據(jù)包絡(luò)譜判斷出混合信號中包含哪些源。

為了判定混合信號中包含幾個源,本文選取一路混合信號進行源數(shù)估計。采用db10小波進行3層小波包分解,得到8個小波包節(jié)點系數(shù)。然后對系數(shù)進行重構(gòu),計算8組重構(gòu)信號的協(xié)方差矩陣,再對協(xié)方差矩陣進行SVD處理,獲得8個特征值,計算非零特征值的下降速比分布如圖6所示。

(a)混合信號1

(b)混合信號2

(c)混合信號3圖5 混合信號包絡(luò)解調(diào)譜Fig.5 Envelop demodulation of mixed signals

圖6 混合信號特征值下降速比曲線Fig.6 Descending ratio curve of eigenvalue of mixed signal

由圖6可知,當N=3時下降速比最大,表明占優(yōu)特征值有三個,即信號的源數(shù)為3,與預(yù)設(shè)情況一致,因此輸入信號的維數(shù)應(yīng)選為3。為了與頻域反卷積方法作對比,先用經(jīng)典的時域Fast-ICA算法嘗試,得到的分離信號解調(diào)譜如圖7所示。

從圖7可以看出,三種模擬故障特征依然混雜在一起,無法區(qū)分。同時也證實了對于卷積混合信號,時域Fast-ICA算法并不適用。

應(yīng)用本文FDBD算法進行處理,算法步驟如圖2所示,STFT采用矩形窗,窗寬和傅里葉變換的長度均設(shè)為960,步長為240,復(fù)數(shù)域Fast-ICA算法采用CFPA算法,通過復(fù)值相關(guān)系數(shù)法進行頻域重排序。圖8為采用頻域反卷積算法得到的解調(diào)譜圖,比較圖8和圖4中混合信號包絡(luò)譜圖可以看出,頻域反卷積算法完全將三種不同且存在諧波關(guān)系的信號分離開來,只是出現(xiàn)次序有所不同。考慮到在實際的故障診斷中,分離次序?qū)τ诠收献R別沒有影響,因此該結(jié)果表明本文所建立的方法對故障特征的分離是有效的。

(a)分離信號1

(b)分離信號2

(c)分離信號3圖7 Fast-ICA分離信號包絡(luò)解調(diào)譜Fig.7 Envelop demodulation spectrum of separated signals by Fast-ICA

(a)分離信號1

(b)分離信號2

(c)分離信號3圖8 FDBD分離信號的包絡(luò)解調(diào)譜Fig.8 Envelop demodulation spectrum of separated signals by FDBD

4 微弱故障特征的分離與識別

在復(fù)雜的齒輪箱系統(tǒng)中,滾動軸承的故障沖擊多會被背景噪聲(如齒輪嚙合振動)所淹沒,給故障診斷帶來很大困擾。本文選取聲振論壇提供的試驗數(shù)據(jù)加以分析。

圖9 齒輪箱試驗臺示意圖Fig.9 Schematic map of gearbox test rig

試驗系統(tǒng)包括齒輪箱試驗臺、MFD310數(shù)據(jù)采集儀、計算機等,其結(jié)構(gòu)簡圖見圖9,齒輪箱的支撐軸承為滾動軸承。軸Ⅱ軸承存在外圈點蝕故障,在箱體表面采集振動加速度信號,采樣點數(shù)為4096,采樣頻率為2048 Hz,其中齒輪嚙合頻率為164.5 Hz,軸Ⅰ轉(zhuǎn)頻為10 Hz,軸Ⅱ轉(zhuǎn)頻為19 Hz。實驗軸承的外圈故障特征頻率為63.78 Hz,內(nèi)圈故障特征頻率為97.38 Hz。

選取試驗臺的一組振動數(shù)據(jù),經(jīng)過本文的源數(shù)估計方法,獲得特征值下降速比的分布如圖10所示,可以判定該信號包含2個振源成分。

圖10 實驗信號特征值下降速比曲線Fig.10 Descending ratio curve of eigenvalue of experimental signal

圖11、圖12所示為試驗臺振動信號的波形及其解調(diào)譜,可見,齒輪嚙合振動是原始信號最突出的一個振源,圖中164.5 Hz即為齒輪的嚙合頻率。與齒輪嚙合振動相比,軸承故障引起的沖擊比較微弱,無論時域還是頻域觀察均無法找到對應(yīng)特征。

圖11 實驗信號時域波形Fig.11 Waveform of experimental signal

圖12 實驗信號包絡(luò)解調(diào)譜Fig.12 Envelop demodulation spectrum of experimental signal

本文利用經(jīng)典的Fast-ICA算法對該信號進行故障識別,處理結(jié)果如圖13所示,圖中依然無法找到和軸承故障相關(guān)的特征。

(a)分離信號1

(b)分離信號2圖13 Fast-ICA分離信號包絡(luò)解調(diào)譜Fig.13 Envelop demodulation spectrum of separated signals by Fast-ICA

圖14 FDBD分離信號1的包絡(luò)解調(diào)譜Fig.14 Envelop demodulation spectrum of separated signal 1 by FDBD

圖14是應(yīng)用本文FDBD算法分離后得到的一路信號包絡(luò)譜圖。從圖中可以明顯看出64 Hz的譜線成分,這和軸承外圈故障特征頻率63.78 Hz非常接近,說明已經(jīng)把軸承故障特征成功提取出來。本例的實踐表明,頻域反卷積技術(shù)對于微弱故障特征的分離和識別也具有較高的研究價值。

5 結(jié)論

(1)針對振動系統(tǒng)中的多振源卷積耦合問題,本文建立了基于信源估計和頻域反卷積的特征分離方法。首先通過小波包分解和SVD實現(xiàn)單通道信號的源數(shù)估計,然后根據(jù)源數(shù)選取相應(yīng)維數(shù)的混合信號,輸入頻域反卷積算法,實現(xiàn)故障特征的分離與識別。

(2)通過本文方法實現(xiàn)了三個呈諧波關(guān)系的仿真沖擊信號的特征分離,并和經(jīng)典的時域Fast-ICA方法進行對比,結(jié)果證明了本文方法的有效性。

(3)針對被齒輪嚙合振動所淹沒的微弱軸承故障特征,本文提供的方法成功實現(xiàn)了微弱特征辨識。同時,與時域Fast-ICA分離結(jié)果的對比表明,齒輪箱振動信號更加符合卷積混合模型。

[2] YPMA A, PAJUNEN P. Rotating Machine Vibration Analysis with Second-order Independent Compon-ent Analysis[C]// Proceedings of the First International Workshop on Independent Component Analysis and Signal Separation Ica. Aussois,1998:37-42.

[4] 張海軍, 屈梁生. 一種提高診斷信息質(zhì)量的方法[J]. 西安交通大學學報,2002, 13(2):295-299. ZHANG Haijun, QU Liangsheng. Method to Improve the Quality of Diagnostic Information[J]. Journal of Xi’an Jiaotong University, 2002,13(2):295-299.

[5] TORKKOLA K. Blind Separation of Delayed Sources Based on Information Maximization[C]// IEEE Proceedings of International Conference on Acoustics, Speech, and Signal Processing. Atlanta, 1996:3510-3513.

[6] SERVIERE C, FABRY P. Blind Source Separation of Noisy Harmonic Signals for Rotating Machine Diagnosis[J]. Journal of Sound and Vibration, 2004, 272:317-339.

[7] PELED R, BRAUN S, Zacksenhouse M. A Blind Deconvolution Separation of Multiple Sources, with Application to Bearing Diagnostics[J]. Mechanical Systems and Signal Processing, 2005, 19:1181-1195.

[8] 劉婷婷, 任興民, 楊永鋒,等. 盲解卷積的機械振動信號分離技術(shù)[J]. 振動、測試與診斷,2009,29(4):419-423. LIU Tingting, REN Xingmin, YANG Yongfeng, et al. The Technology of Mechanical Vibration Signals Separation Based on Blind Deconvolution[J].Journal of Vibration, Measurement and Diagnosis,2009,29(4):419-423.

[9] 葉紅仙, 楊世錫, 楊將新. 多振源卷積混合的時域盲源分離算法[J]. 機械工程學報,2009, 45(1):190-199. YE Hongxian, YANG Shixi, YANG Jiangxin. Temporal Blind Source Separation Algorithm for Convolution Mixtures with Muti Vibration Sources[J]. Chinese Journal of Mechanical Engineering, 2009, 45(1):190-199.

[10] 潘楠, 伍星, 遲毅林, 等. 基于頻域盲解卷積的機械設(shè)備狀態(tài)監(jiān)測與故障診斷[J]. 振動與沖擊,2012, 31(12):34-41. PAN Nan, WU Xing, CHI Yilin, et al. Mechanical Equipment Condition Monitoring and Fault Diagnosis Based on Frequency-domain Blind Deconvolution[J]. Journal of Vibration and Shock, 2012, 31(12):34-41.

[11] 焦衛(wèi)東, 楊世錫, 吳昭同. 基于源數(shù)估計的旋轉(zhuǎn)機械源盲分離[J]. 中國機械工程, 2003, 14(14):1184-1188. JIAO Weidong, YANG Shixi, WU Zhaotong. A Method of Blind Source Saparation for Rotating Machinery Based on Estimation of the Number of Sources[J]. China Mechanical Engineering, 2003, 14(14):1184-1188.

[12] 葉紅仙, 楊世錫, 楊將新. 基于EMD-SVD-BIC的機械振動源數(shù)估計方法[J]. 振動、測試與診斷, 2010,30(6):330-335. YE Hongxian, YANG Shixi, YANG Jiangxin. Mechanical Vibration Source Number Estimation Based on EMD-SVD-BIC[J]. Journal of Vibration, Measurement &Diagnosis, 2010, 30(6):330-335.

[13] BINGHAM E,HYVRINEN A. A Fast Fixed-point Algorithm for Independent Component Analysis of Complex-valued Signals[J]. International Journal of Neural Systems, 2000, 10(1):1-8.

(編輯 王旻玥)

Fault Separation and Recognition of Rolling Bearings Based on Source Number Estimation and FDBD

ZHANG Jianyu1MENG Hao2XU Yonggang3

1.Key Laboratory of Advanced Manufacturing Technology, Beijing University of Technology, Beijing, 100124 2.Zhangqiu Haier Motor Co. Ltd., Wolong Electrics, Zhangqiu, Shandong,2502003. Beijing Engineering Research Center of Precision Measurement Technology and Instrument, Beijing, 100124

According to the problems that the vibration features of bearing faults were hard to separate and recognize in strong vibration source inhibition, a diagnosis method was established based on source number estimation and FDBD algorithm. Wavelet packet decomposition was used to divide the signals into multiple sub band signals, and SVD was selected to estimate the signal source numbers in underdetermined conditions. The multiple dimension signals were constructed based on the source number estimation. The FDBD algorithm, which included STFT, fast-ICA in complex domain, relevance ranking and inverse STFT, was finally applied on fault feature separation and extraction. The effectiveness of the method was validated in fault feature separation and weak feature recognition by the simulation signals and experimental data of rolling bearing faults.

wavelet packet decomposition; singular value decomposition(SVD); short time Fourier transform (STFT); fast-ICA in complex domain; frequency domain blind deconvolution (FDBD)

2016-01-21

北京市教委科技計劃資助項目(KM201410005027)

TH16

10.3969/j.issn.1004-132X.2017.01.008

張建宇,男,1975年生。北京工業(yè)大學機電學院副教授。主要研究方向為機電設(shè)備故障診斷。發(fā)表論文20余篇。E-mail:zhiy-1999@bjut.edu.cn。孟 浩,男,1989年生。臥龍電氣章丘海爾電機有限公司工程師。胥永剛,男,1975年生。北京工業(yè)大學機電學院副教授。

猜你喜歡
故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
孩子停止長個的信號
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
故障一點通
故障一點通
故障一點通
主站蜘蛛池模板: 777午夜精品电影免费看| 在线播放国产99re| 国产91在线免费视频| 久久无码av三级| 亚洲精品无码日韩国产不卡| 日本一本在线视频| 欧美性天天| 香蕉蕉亚亚洲aav综合| 国产乱子伦视频在线播放| 亚洲AV永久无码精品古装片| 激情影院内射美女| 高潮毛片无遮挡高清视频播放| 亚洲国产在一区二区三区| 免费看黄片一区二区三区| 国产成人精品视频一区视频二区| 久久久久88色偷偷| 成人字幕网视频在线观看| 熟女视频91| 国产永久无码观看在线| 好久久免费视频高清| 久久亚洲综合伊人| 欧美色图久久| 一级毛片不卡片免费观看| 一区二区三区毛片无码| 欧美一级99在线观看国产| 自拍偷拍欧美| 久久99久久无码毛片一区二区| 亚洲人成电影在线播放| 熟妇无码人妻| 最新国产在线| 色综合手机在线| 午夜精品福利影院| 毛片大全免费观看| 依依成人精品无v国产| 国产精品理论片| 97在线碰| 久久这里只精品国产99热8| 久草中文网| 国产在线啪| 日韩一区二区三免费高清| 亚洲毛片网站| 免费看美女自慰的网站| 亚洲欧美激情小说另类| 夜夜操狠狠操| 97青草最新免费精品视频| 国产性精品| 国产黄在线免费观看| 日韩第八页| 亚洲永久视频| 免费观看欧美性一级| 91美女在线| 67194亚洲无码| 亚洲国产系列| 国产视频自拍一区| 高清无码一本到东京热| 亚洲第一区欧美国产综合| 国产精品人人做人人爽人人添| 2021国产精品自产拍在线| 亚洲品质国产精品无码| 一级毛片中文字幕| 国产jizzjizz视频| 欧美日韩一区二区在线播放 | 99九九成人免费视频精品| 精品99在线观看| 人人看人人鲁狠狠高清| 国产精品成人观看视频国产| 国产精品女主播| 久久国产拍爱| 亚洲中文字幕久久精品无码一区 | 手机看片1024久久精品你懂的| 国产精品亚洲va在线观看| 亚洲黄色激情网站| 女人18一级毛片免费观看| 99热这里只有精品在线观看| 毛片视频网| 无码内射在线| 色老二精品视频在线观看| 亚洲国内精品自在自线官| 刘亦菲一区二区在线观看| 国产综合网站| 青青操视频免费观看| 中文纯内无码H|