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

基于集合經驗模態分解敏感固有模態函數選擇算法的滾動軸承狀態識別方法

2014-05-29 08:42:08王玉靜康守強姜義成Mikulovich
電子與信息學報 2014年3期
關鍵詞:特征提取模態振動

王玉靜 康守強 張 云 劉 學 姜義成 Mikulovich V I

?

基于集合經驗模態分解敏感固有模態函數選擇算法的滾動軸承狀態識別方法

王玉靜①②康守強②張 云①劉 學②姜義成*①Mikulovich V I③

①(哈爾濱工業大學電子與信息工程學院 哈爾濱 150001)②(哈爾濱理工大學電氣與電子工程學院 哈爾濱 150080)③(白俄羅斯國立大學 明斯克 220030)

為了更有效地提取滾動軸承各狀態振動信號的特征,該文提出了一種基于集合經驗模態分解(EEMD)的敏感固有模態函數(IMF)選擇算法。該算法對振動信號經EEMD分解后得到的固有模態函數采用峭度值、相關系數相結合的方法自動提取其敏感分量,以此獲得振動信號的初始特征。再運用奇異值分解和自回歸(AR)模型方法得到滾動軸承各狀態振動信號的特征向量,并將其輸入到改進的超球多類支持向量機中進行智能識別,從而實現滾動軸承的正常狀態,不同故障類型及不同性能退化程度的各狀態識別。實驗結果表明,相比基于經驗模態分解結合自回歸模型或奇異值分解的特征提取方法,該方法可更有效地提取滾動軸承故障特征信息,且識別精度更高。

信號處理;狀態識別;非平穩信號;集合經驗模態分解(EEMD);敏感固有模態函數(IMF)

1 引言

滾動軸承是支撐旋轉軸的關鍵部件,廣泛應用于各種機械設備之中。現有的診斷方法大多集中在滾動軸承故障位置的診斷,而不同故障程度的識別是最近研究的新方向[1]。這也對提取滾動軸承振動信號的有效特征提出了新的挑戰。

經驗模態分解(Empirical Mode Decomposition, EMD)方法適合分析處理非平穩、非線性信號[2]。基于EMD的特征提取方法在各領域得到了廣泛的應用,如文獻[3]采用基于EMD的希爾伯特能量譜作為語音和非語音的鑒別特征,結合順序統計濾波器提高了對語音信號的檢測率;文獻[4]提出一種基于EMD的彈道目標平動補償與微多普勒特征提取方法,并通過仿真實驗驗證了其有效性和魯棒性;文獻[5]采用EMD方法來提取齒輪故障振動信號的特征;文獻[6,7]采用EMD結合AR模型分別對滾動軸承故障類型、故障類型及不同性能退化程度振動信號進行特征提取,獲得了較好效果。

EMD本身也存在一些不足,為了抑制EMD的模態混疊現象,文獻[8]提出了一種集合經驗模態分解(Ensemble EMD, EEMD)方法。文獻[9]對被測信號序列進行等間隔二次采樣,之后再運用EEMD的思想對信號進一步處理,應用到穿墻雷達人的運動微多普勒特性分析中,獲得了較好效果。但EEMD分解結果取決于人為設定的總體平均次數與加入噪聲的幅值大小這兩個參數,還受篩選條件的限制,導致EEMD分解出的固有模態函數(IMF)仍可能存在模態混疊。再者,在EEMD分解過程中,分解出IMF的個數受原始信號點數的制約,這就可能造成EEMD的過分解[8]。從另一個角度,滾動軸承振動信號經EEMD,得到的IMF分量中,通常只有部分IMF包含大量故障信息,對故障敏感。因此,非常有必要研究一種方法來選擇有效的,且對故障敏感的IMF,剔除噪聲干擾成份或者與故障無關的偽分量,凸顯故障特征。文獻[10]提出了利用各IMF峭度最大來選取滾動軸承振動信號EEMD分解后的敏感IMF的方法。各IMF峭度值可能有多個大于3,對利用IMF進一步提取特征來說,取最大的峭度值對應的一個IMF會丟失部分故障信息。目前廣泛應用的敏感IMF選擇算法是相關系數法,計算各IMF與原信號的歸一化相關系數,但該方法需要人為確定一個閾值來選擇敏感IMF分量,文獻[11, 12]中,設定閾值為最大相關系數的1/10。對于不同的信號,該閾值同樣需要人為調整。

另一方面,在智能分類方法中,支持向量機(SVM)有其獨特的優勢并被廣泛應用。文獻[13]在基于EMD的提取特征基礎上,采用SVM對常規低分辨雷達體制下空中飛機目標進行分類,獲得了較好的識別效果。對于多分類問題,基于超球結構支持向量機[14]及改進的分類規則[15],文獻[7]又進一步進行了完善,并應用到滾動軸承的多類故障分類中,獲得了較好效果。

因此,本文在分析EEMD方法理論的基礎上,提出了峭度值與相關系數相結合的自動提取滾動軸承振動信號敏感IMF算法,之后結合奇異值分解(SVD),自回歸(AR)模型參數估計進一步提取特征,并引入改進的超球多類支持向量機來進行多狀態識別,最終實現滾動軸承的正常狀態、故障類型(內、外環)及不同性能退化程度的各狀態識別。

2 振動信號特征提取

2.1 集合經驗模態分解

EEMD方法的本質是一種疊加高斯白噪聲的多次經驗模式分解,利用了高斯白噪聲具有頻率均勻分布的統計特性,通過每次加入同等幅值的不同白噪聲來改變信號的極值點特性,之后對多次EMD得到的相應IMF進行總體平均來抵消加入的白噪聲,從而有效減少模式混淆問題[8]。EEMD分解步驟如下:

(1)初始化總體平均次數;

(4)對所得含噪聲的信號x()分別進行EMD分解,得到如式(2)的各自IMF和的形式:

式中c,j()為第次加入白噪聲后分解得到的第個IMF,r,j()是殘余函數,代表信號的平均趨勢,是IMF的數量;

(5)重復步驟(3)和步驟(4)進行次,每次分解加入幅值不同的白噪聲信號得到IMF的集合為

(6)利用不相關序列的統計平均值為零的原理,將上述對應的IMF進行集合平均運算,得到EEMD后的最終IMF,即

(7)EEMD分解的最終結果為

式中()為殘余分量。

2.2 自動提取敏感IMF分量算法

峭度系數是一個常用的時域統計指標,正常狀態下,高斯分布的峭度系數約為3,當峭度系數大于3時就預示著故障的出現,根據不同頻帶的峭度變化趨勢可以更有效地選擇對故障敏感頻帶[16]。

峭度指標定義:

峭度系數是無量綱的時域指標,它波動較大,不適合強噪聲環境下的狀態監測問題。對最初始退化時微弱故障反應不十分敏感,在故障后期(故障程度較大)退化過程中峭度值呈高-低-高的變化[17]。因此,單獨根據峭度值不能有效衡量故障程度,也不能準確判斷出故障的類型。

綜上,為了不丟失故障信息,并且達到初步特征提取,獲取敏感IMF的目的,本文采用自動相關系數法和IMF峭度系數法(峭度系數大于3)分別選取IMF,并以二者的并集作為敏感IMF。為進一步進行特征提取做準備。

2.3 構造特征向量矩陣及特征提取流程

為了進一步提取振動信號的特征,在上述提出敏感的IMF分量的基礎上,進一步采用AR模型和奇異值分解的方法來構造特征向量。

由于AR模型參數、殘差方差和奇異值在數值單位上差距較大,為了使后續分類不出現權值分配混亂現象,對所有特征值采用線性函數進行歸一化處理。因此,構造每個振動信號歸一化后的特征向量矩陣

式中max為AR模型的最大階數,小于max的模型參數補零,使所有參數個數一致;為敏感IMF的數量。

特征提取流程如圖1所示。

圖1 振動信號特征提取流程圖

對于滾動軸承多狀態分類,由于學習信號和測試信號不止一個,而多個信號分別經過EEMD分解提取敏感IMF分量的個數可能不同,設最大個數為,不足的用零向量進行補充。之后對每個振動信號的個IMF分別求AR模型參數、殘差方差,同時將個IMF組成特征向量矩陣求取個奇異值,構造歸一化特征向量矩陣為

3 滾動軸承各狀態識別

由于改進分類規則的超球結構多類支持向量機處理多分類問題時,在補充新類樣本、各樣本間數據不均衡、分類精度方面具有一定的優勢[7]。因此,本文選用該分類器對滾動軸承多狀態進行分類。根據文獻[7]改進的分類規則重寫如下:

測試特征向量不包含在關鍵區域,則

測試特征向量包含在關鍵區域,該區域中的訓練樣本集合為空時,則

測試特征向量包含在關鍵區域,該區域中的訓練樣本集合不為空時,則

式中所有參數的含義參見文獻[7]。

4 滾動軸承各狀態識別方法流程

基于集合經驗模態分解,自動提取敏感IMF分量算法,奇異值分解和改進超球多類支持向量機的滾動軸承各狀態識別方法的步驟如下:

(1)采集大量的滾動軸承振動信號,包括正常狀態、不同故障類型及不同性能退化程度振動信號;

(2)按8倍交叉驗證法將信號分為訓練信號和測試信號;

(3)對各狀態的每個訓練信號進行特征提取,按照式(9)構造特征向量矩陣;

(4)各狀態的所有學習信號構造特征向量矩陣

(5)求得各狀態超球的參數,確定初始狀態超球;

(6)采用相同的方法求得測試信號的特征,按照式(9)構造特征向量矩陣;

(7)基于特征向量矩陣,根據式(10),式(11)和式(12)判斷測試信號對應的滾動軸承的狀態。當診斷正確率最高時,確定超球的參數;

(8)實際現場振動信號,同樣提取出特征向量,根據式(10),式(11)和式(12)評判現場滾動軸承的健康狀態。

5 應用與分析

基于美國西儲大學軸承數據中心的型號為6205-2RS滾動軸承振動數據[18],本實驗利用該滾動軸承正常工作情況下、內外環故障情況下的5種狀態數據進行研究。其中滾動軸承內、外環不同故障程度的數據包括:內、外環損傷直徑0.18mm的故障數據和內、外環損傷直徑0.53mm的故障數據。5種狀態滾動軸承數據,每個樣本數據長度1024點,信號的采樣頻率為12 kHz。

表1給出了部分滾動軸承在正常狀態和故障狀態下,對振動信號進行EEMD分解得到的IMF的峭度值和振動信號與各IMF之間的相關系數。

由表1可以看出,滾動軸承故障狀態振動信號敏感IMF的選取取決于峭度值法和相關系數法,為了獲取全面的故障信息,取二者的并集。比如內環故障,大于3峭度值的IMF有IMF1, IMF2, IMF4,相關系數法選取小于閾值0.2973的敏感IMF有IMF1, IMF2, IMF3,事實上通過求取IMF3和IMF4的包絡譜可以得到其包含有故障頻率164 Hz的信息,因此取二者的并集作為敏感IMF。對于滾動軸承正常狀態震動信號,敏感IMF的選擇取決于相關系數法。另外,如果采用最大相關系數的1/10作為相關系數法的閾值(分別為0.05436, 0.08892, 0.09835),則選取故障振動信號的敏感IMF均多出2~3個IMF,多出的IMF包絡譜上故障頻率處的幅值約為零。

表1 IMF的峭度值及與振動信號的相關系數

圖2 內環故障信號及其EEMD分解

從對表1中給出的滾動軸承3個振動信號的分析可以得出,取峭度值法和相關系數法的并集確定敏感IMF數量,正常的為前5個,內環故障為前4個,外環故障為前4個。通過對實驗中大量滾動軸承振動信號進行敏感IMF提存,得知敏感IMF的最大個數=5,因此對于敏感IMF個數小于5的,用零向量進行補充。按照式(9)構造一個信號的特征向量矩陣,滾動軸承第狀態的N個振動信號按照式(13)構造特征向量矩陣。將輸入到改進的超球支持向量機中進行訓練,獲得的分類器的參數如表2所示。表2給出了基于其它特征提取方法的分類器參數。在該參數下獲得的平均識別率和一個振動信號識別的平均時間如圖3所示。

表2故障診斷方法比較

特征提取方法改進的超球支持向量機參數 Cs? EMD+Yule-Walker0.232.10.965 EMD+Ulrych-Clayton0.132.10.940 EMD+SVD0.585.10.975 EEMD+敏感IMF+(SVD結合Ulrych-Clayton)0.142.10.955

從圖3可以看出,本文方法比其它方法在識別一個振動信號狀態時所消耗的平均時間多。原因為敏感IMF的提取算法計算量小,并且該算法選出來的敏感IMF由原來經驗選取的6個降到5個,使后續的特征提取和分類耗時減小。但由于EEMD多次加入高斯白噪聲進行多次分解,其耗時比EMD耗時多,AR模型的Ulrych-Clayton[7]參數估計較Yule- Walker[7]或SVD耗時多;上述原因綜合使得總的識別時間增加。但同時換取了更高的平均識別率,比已有的基于EMD結合Ulrych-Clayton的特征提取方法平均識別率高出2%,可見本文所提出的特征提取方法更有效。

6 結論

本文提出了一種基于集合經驗模態分解的峭度值結合相關系數的敏感固有模態函數選擇算法,相關系數法的閾值采用所有相關系數的標準差。為了獲取更有效更全面的故障特征,敏感IMF由峭度值法和相關系數法的敏感IMF的并集構成。該方法減小了人為選取敏感IMF的主觀性。基于此,提出EEMD敏感IMF, SVD, AR模型結合改進的超球多類支持向量機方法實現滾動軸承的正常狀態、不同故障位置及性能退化程度的各狀態識別。

滾動軸承振動實驗結果表明,敏感IMF選擇算法可以全面有效地提取包含故障信息的IMF。基于敏感IMF的EEMD, SVD, AR模型的Ulrych- Clayton參數估計的特征提取方法,再利用改進的超球多類支持向量機來識別滾動軸承各狀態的方法是有效的。與基于EMD結合AR模型參數或SVD的方法相比所耗時間雖然長,但換取了較高的平均識別率。

[1] Pan Yu-na, Chen Jin, and Li Xing-lin. Bearing performance degradation assessment based on lifting wavelet packet decomposition and fuzzy c-means[J]., 2010, 24(2): 559-566.

[2] Huang N E, Shen Z, Long S R,.. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time A series analysis[J]., 1998, 454(1971): 903-995.

[3] 盧志茂, 金輝, 張春祥, 等. 基于HHT和OSF的復雜環境語音端點檢測[J]. 電子與信息學報, 2012, 34(1): 213-217.

Lu Zhi-mao, Jin Hui, Zhang Chun-xiang,.. Voice activity detection in complex environment based on Hilbert-Huang transform and order statistics filter[J].&, 2012, 34(1): 213-217.

[4] 羅迎, 柏又青, 張群, 等. 彈道目標平動補償與微多普勒特征提取方法[J]. 電子與信息學報, 2012, 34(3): 602-608.

Luo Ying, Bai You-qing, Zhang Qun,.. Translational motion compensation and micro-doppler feature extraction of ballistic targets[J].&, 2012, 34(3): 602-608.

[5] 程軍圣, 于德介, 楊宇, 等. 基于EMD的齒輪故障識別研究[J]. 電子與信息學報, 2004, 26(5): 825-829.

Cheng Jun-sheng, Yu De-jie, Yang Yu,.. Research on gear fault diagnosis based on EMD[J].&, 2004, 26(5): 825-829.

[6] Cheng Ju-sheng, Yu De-jie, and Yang Yu. A fault diagnosis approach for roller bearings based on EMD method and AR model[J]., 2006, 20(2): 350-362.

[7] 康守強, 王玉靜, 楊廣學, 等. 基于經驗模態分解和超球多類支持向量機的滾動軸承故障診斷方法[J]. 中國電機工程學報, 2011, 31(14): 96-102.

Kang Shou-qiang, Wang Yu-jing, Yang Guang-xue,.. Rolling bearing fault diagnosis method using empirical mode decomposition and hypersphere multiclass support vector machine[J]., 2011, 31(14): 96-102.

[8] Wu Zhao-hua and Huang Norden E. Ensemble empirical mode decomposition: a noise assisted data analysis method[J]., 2009, 1(1): 1-41.

[9] 王宏, Narayanan R M, 周正歐, 等. 基于改進EEMD的穿墻雷達動目標微多普勒特性分析[J]. 電子與信息學報, 2010, 32(6): 1355-1360.

Wang Hong, Narayanan R M, Zhou Zheng-ou,.. Micro- doppler character analysis of moving objects using through- wall radar based on improved EEMD[J].&, 2010, 32(6): 1355-1360.

[10] 胡愛軍, 馬萬里, 唐貴基. 基于集成經驗模態分解和峭度準則的滾動軸承故障特征提取方法[J]. 中國電機工程學報, 2012, 32(11): 106-111.

Hu Ai-jun, Ma Wan-li, and Tang Gui-ji. Rolling bearing fault feature extraction method based on ensemble empirical mode decomposition and kurtosis criterion[J].,2012, 32(11): 106-111.

[11] 林麗, 余輪. 基于相關系數的EMD改進算法[J]. 計算機與數字工程, 2008, 36(12): 28-29,38.

Lin Li and Yu Lun. Improvement on empirical mode decomposition based on correlation coefficient[J].&, 2008, 36(12): 28-29,38.

[12] Peng Z K, Tse P W, and Chu F L. A comparison study of improved Hilbert-Huang transform and wavelet transform: application to fault diagnosis for rolling bearing[J]., 2005, 19(5): 974-988.

[13] 王寶帥, 杜蘭, 劉宏偉, 等. 基于經驗模態分解的空中飛機目標分類[J]. 電子與信息學報, 2012, 34(9): 2116-2121.

Wang Bao-shuai, Du Lan, Liu Hong-wei,.. Aircraft classification based on empirical mode decomposition[J].&, 2012, 34(9): 2116-2121.

[14] 朱美琳, 劉向東, 陳世福. 用球結構的支持向量機解決多分類問題[J]. 南京大學學報 (自然科學版), 2003, 39(2): 153-158.

Zhu Mei-lin, Liu Xiang-dong, and Chen Shi-fu. Salving the problem of Multi-class pattern recognition with sphere- structured support vector machines[J].(), 2003, 39(2): 153-158.

[15] 故磊, 吳慧中, 肖亮. 基于新的決策規則的球形支持向量機分類算法[J]. 系統仿真學報, 2008, 20(11): 2901-2904.

Gu Lei, Wu Hui-zhong, and Xiao Liang. Sphere-structured support vector machines classification algorithm based on new decision rule[J]., 2008, 20(11): 2901-2904.

[16] 楊國安. 機械設備故障診斷實用技術[M]. 北京: 中國石化出版社, 2007: 244-245.

[17] 潘玉娜, 陳進. 小波包-支持向量數據描述在軸承性能退化評估中的應用研究[J]. 振動與沖擊, 2009, 28(4): 164-167.

Pan Yu-na and Chen Jin. Wavelet package-support vector data description applied in bearing performance degradation assessment[J]., 2009, 28(4): 164-167.

[18] Loparo K A. Bearing data center [EB/OL]. http://www. eecs.case.edu/laboratory/bearing/welcome_overview.htm, 2013.

[19] 陳略, 訾艷陽, 何正嘉, 等. 總體平均經驗模式分解與1.5維譜方法的研究[J]. 西安交通大學學報, 2009, 43(5): 94-98.

Chen Lue, Zi Yan-yang, He Zheng-jia,.. Research and application of ensemble empirical mode decomposition principle and 1.5 dimension spectrum method[J]., 2009, 43(5): 94-98.

王玉靜: 女,1983年生,博士生,研究方向為信號與信息處理.

康守強: 男,1980年生,副教授,研究方向為信號與信息處理.

姜義成: 男,1964年生,教授,研究方向為信號與信息處理.

Condition Recognition Method of Rolling Bearing Based on Ensemble Empirical Mode Decomposition Sensitive Intrinsic Mode Function Selection Algorithm

Wang Yu-jing①②Kang Shou-qiang②Zhang Yun①Liu Xue②Jiang Yi-cheng①Mikulovich V I③

①(,,150001,)②(,,150080,)③(,220030,)

In order to extract effectively the characteristics of each condition vibration signal for rolling bearing, a sensitive Intrinsic Mode Function (IMF) selection algorithm which based on Ensemble Empirical Mode Decomposition (EEMD) is proposed. First, for obtaining the initial characteristics of the vibration signal, the vibration signal is decomposed by using EEMD, and the sensitive components of obtained IMFs are extracted automatically by using kurtosis combined with correlation coefficient. Then, the feature vectors of each condition vibration signal of rolling bearing are obtained by using Singular Value Decomposition (SVD) and AutoRegressive (AR) model. The obtained feature vectors are regarded as the input of the improved hyper-sphere multi-class Support Vector Machine (SVM) for intelligent recognition. Thereby, the condition recognition of normal state, different fault types and different degrees of performance degradation of rolling bearing can be achieved. The experimental results show that, the proposed method can effectively extract fault characteristics information of rolling bearing more than EMD combined with AR model and EMD combined with SVD method, and the recognition rate is higher.

Signal processing; Condition recognition; Nonstationary signal; Ensemble Empirical Mode Decomposition (EEMD); Sensitive Intrinsic Mode Function (IMF)

TN911.7; TP18

A

1009-5896(2014)03-0595-06

10.3724/SP.J.1146.2013.00434

2013-04-02收到,2013-10-31改回

國家自然科學基金(51305109),高等學校博士學科點專項科研基金(20122303120010),留學人員科技活動項目擇優資助和哈爾濱市科技創新人才專項基金(留學回國人員)(2013RFLXJ019)資助課題

姜義成 jiangyc@hit.edu.cn

猜你喜歡
特征提取模態振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
中立型Emden-Fowler微分方程的振動性
一種基于LBP 特征提取和稀疏表示的肝病識別算法
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
基于MED和循環域解調的多故障特征提取
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 亚洲国产91人成在线| 日本成人一区| 55夜色66夜色国产精品视频| 国产精品亚洲欧美日韩久久| 亚洲欧美精品在线| 国产免费久久精品44| 中文天堂在线视频| 亚洲欧美日韩另类在线一| 亚洲自拍另类| 亚洲欧美成人综合| 99re热精品视频国产免费| 幺女国产一级毛片| 欧美日韩午夜| 欧美激情福利| 伊人天堂网| 日本黄色a视频| 久久精品无码国产一区二区三区| 成人日韩精品| 美女一区二区在线观看| 国产丝袜无码一区二区视频| 免费女人18毛片a级毛片视频| 国产丝袜丝视频在线观看| 99国产精品一区二区| 在线国产91| 国产精品内射视频| 亚洲免费三区| 亚洲成综合人影院在院播放| 欧美激情网址| 99在线观看国产| 99久久国产精品无码| 国产尤物视频网址导航| 亚洲人成网站18禁动漫无码| 毛片一级在线| 日韩精品无码免费专网站| 精品99在线观看| 亚洲精品视频免费看| 欧美精品亚洲精品日韩专区| 18禁黄无遮挡免费动漫网站| 91久久精品国产| 99视频在线免费观看| 欧美不卡在线视频| 色婷婷成人| 伊人久久久久久久久久| 亚洲无码一区在线观看| 久久福利网| 色妞永久免费视频| 久草网视频在线| 国产幂在线无码精品| 亚洲制服丝袜第一页| 欧美日本视频在线观看| 久久中文字幕不卡一二区| 中文字幕va| 在线精品视频成人网| 色一情一乱一伦一区二区三区小说| 日韩小视频在线观看| 亚洲高清在线天堂精品| 国产精品一区二区不卡的视频| 久久精品无码国产一区二区三区| av天堂最新版在线| 好久久免费视频高清| 一级片一区| 亚洲欧美日韩成人高清在线一区| 一边摸一边做爽的视频17国产| 国产 日韩 欧美 第二页| 日本三级黄在线观看| 国产aⅴ无码专区亚洲av综合网| 久久精品国产精品青草app| 亚洲免费三区| 亚洲啪啪网| 国产成年女人特黄特色毛片免 | 婷婷色丁香综合激情| 亚洲av片在线免费观看| 亚洲一级毛片在线观| 欧美成人精品一级在线观看| 粉嫩国产白浆在线观看| 97免费在线观看视频| 精品一区二区三区视频免费观看| 欧美日韩激情在线| 26uuu国产精品视频| 国产精品白浆无码流出在线看| 亚洲男女天堂| 综合色区亚洲熟妇在线|