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

基于拉普拉斯特征映射的滾動軸承故障識別

2016-01-12 10:37:45黃宏臣,韓振南,張倩倩
振動與沖擊 2015年5期
關鍵詞:模式識別特征提取

第一作者黃宏臣男,碩士生,1986年10月生

通信作者韓振南男,教授,博士生導師,1958年2月生

基于拉普拉斯特征映射的滾動軸承故障識別

黃宏臣,韓振南,張倩倩,李月仙,張志偉

(太原理工大學機械工程學院, 太原030024)

摘要:用傳統的線性方法對非平穩和非線性運行狀態的滾動軸承進行故障診斷時,效果欠佳。為了及時、準確地監測軸承的運行狀態,提出了將拉普拉斯特征映射算法(Laplacian Eigenmap LE)應用到滾動軸承的故障識別中。在振動信號構建的時域和頻域高維特征空間矩陣中,充分利用LE算法在非線性特征提取和降維的優點,進行學習,提取表征軸承狀態的特征量,并以可視化的聚類結果進行表示。實驗模擬了軸承的4種不同類型故障以及滾動體的4種不同受損程度,采用模式識別中聚類性的類內距和類間距兩個參數作為衡量指標。與PCA和KPCA兩種方法對比,LE不僅明顯識別出四種故障類型和有效的區分出滾動體的不同受損程度,而且識別率大大提高。并通過測試樣本組驗證了LE方法的有效性。

關鍵詞:滾動軸承故障;流形學習;模式識別;拉普拉斯特征映射;特征空間的構建;特征提取;測試樣本驗證

基金項目:國家自然科學基金資助項目(50775157); 山西省高等學校留學回國人員科研資助項目(2011-12); 山西省基礎研究項目(2012011012-1)

收稿日期:2013-10-09修改稿收到日期:2014-07-07

中圖分類號:TH133.3文獻標志碼:A

Method of fault diagnosis for rolling bearings based on Laplacian eigenmap

HUANGHong-chen,HANZhen-nan,ZHANGQian-qian,LIYue-xian,ZHANGZhi-wei(College of Mechanical Engineering, Taiyuan University of Technology, Taiyuan 030024, China)

Abstract:The traditional linear methods for rolling bearing fault diagnosis under non-stationary and nonlinear running status are not effective. In order to monitor the rolling bearing status accurately and timely, a new diagnosis method was put forward by applying the algorithm of Laplacian eigenmap (LE) to the fault diagnosis of rolling bearings. With the method, the advantage of LE algorithm was fully utilized for extracting nonlinear features and reducing dimensions of characteristic space matrix constructed with vibration signals in time domain and frequency domain, the features of running status were extracted and the clustering results were visualized. Two parameters (between-class scatter and within-distance in pattern recognition) were used in experiments as the measurable indicators to identify four different faults of bearings and four different levels of ball damages in bearings. Compared with PCA & KPCA, it was shown that the LE algorithm can more clearly identify the four different faults and the different levels of ball damages, and its recognition efficiency rises greatly. The effectiveness of LE was verified with test samples.

Key words:rolling bearing fault; manifold learning; pattern recognition; Laplacian eigenmap; construction of characteristic space; feature extraction; validation by using test samples

滾動軸承是各類旋轉機械中最常用、最易損壞的通用零部件。據統計,旋轉機械中30%的故障是由軸承缺陷引起,軸承運行的好壞對機器的工作狀況影響極大,它的缺陷會導致機器劇烈振動和產生噪聲,甚至會導致設備的損壞而引起巨大的經濟損失[1]。因此,對軸承進行工況監測和故障診斷是非常必要的。基于滾動軸承故障沖擊特征的振動信號提取方法的研究不斷深入,新方法不斷被提出。滾動軸承的所測得振動信號是典型的非線性和非平穩的,在出現故障時振動表現得更加明顯。為了準確的監測設備運行狀態,通常采用多個相關傳感器進行多方面的監測,數據之間存在著高度相關性,但也造成了大量的冗余信息形成了海量的高維數據。如何在海量、高維、非線性數據中提取出有用的信息是機器學習領域的熱點也是目前設備故障診斷所面臨的最大問題之一。

主元分析法(Principal Components Analysis PCA)和獨立主元分析(Independent Components Analysis ICA)是傳統的兩種線性分析方法,對線性結構的數據集具有較好的效果,但由于這些方法的線性本質使得無法揭示復雜的非線性流形結構[1]。核方法的出現為非線性學習提供了一條解決途徑,核函數的引入將數據映射到可以發現線性關系的高維特征空間中,在特征空間進行線性學習,實現了非線性問題的高效求解,同時避免了復雜的非線性映射的求取。核主元分析法(Kernel Principal Component Analysis KPCA)就是其中的典型方法。KPCA巧妙運用核函數在希爾伯特特征空間執行線性PCA,從而實現非線性特征的提取。KPCA的非線性特征提取在故障診斷應用上明顯優越于PCA[2-4]。但是,核函數及參數在選取時是需要預先給定的,且在使用核函數時不需要知道具體的特征空間,使得核方法缺乏物理直觀性。并且KPCA中的核函數對于故障診斷中的類可分性測度影響很大,選擇不當會導致類可分性測度很小,使得不同類別的樣本點在特征空間中出現相互重疊覆蓋。而流形學習的出現給機械故障診斷提供了一種新的思路。

提出使用流形學習方法從滾動軸承故障非線性數據中提取出我們想要的信息,針對高維非線性數據的處理困難,采用拉普拉斯特征映射方法(Laplacian Eigenmap LE)對由軸承振動信號提取出的時域和頻域所構成的原始高維特征空間進行學習,提取出隱藏在高維特征數據中的內在低維流形特征,以可視化的方式識別出軸承的異常狀態以及受損程度,通過與PCA和KPCA兩種方法對比來驗證它的可行性和有效性,為軸承的運行狀態分析提供了一條新的途徑。

1流形學習方法

流形學習是Science[5-7]在同期發表了三篇有關流形學習的論文后而被人們所認知,由于流形算法在非線性、高維數據處理方面具有強大的優勢,一經提出便迅速成為眾多領域的研究熱點。并相繼提出多種經典方法,如LLE、ISOMAP、LE和LTSA等,并在不斷改進。較之于傳統方法,流形學習方法能夠快速、有效地發現非線性高維數據中的本質結構并用低維表示,實現了維數約簡和非線性數據分析。

流形學習方法是一種基于微分幾何與拓撲學的非線性高維數據處理方法,目的就是從高維特征空間中尋找到最能表征滾動軸承運行狀態的低維特征數據,為準確診斷滾動軸承的健康狀態提供依據。

2拉普拉斯特征映射算法(LE)

拉普拉斯特征映射(LE)是一種典型的流形學習方法,它是運用圖拉普拉斯概念計算高維特征集得到低維流形表示。LE算法降維的實質就是尋找到一個平均意義上保持數據點局部鄰近信息,即在原來高維特征空間中是近鄰的點在低維表示中也應該是近鄰的。LE算法所產生的映射可以看作是對幾何流形的一種連續離散逼近的映射,用數據點的鄰域圖來近似表示流形,并用Laplace-Beltrami算子近似表示鄰域圖的權值矩陣,實現高維流形的最優嵌入。這樣就將問題簡單轉化為對拉普拉斯算子的特征函數的求解。因為LE算法有良好的保持局部特征的性質,所以使得它具有對噪聲不敏感的優點,即使只使用局部距離,也不易出現短路等問題。

LE算法可表述為以下三步:

(1)構建近鄰圖。如果Xi和Xj為近鄰點,就在節點i和j之間置一條邊,目前有ε-鄰域法和K近鄰域法兩種方法,根據參考文獻[8]在本文中將采用K=8的K近鄰域法構建近鄰圖。

(2)確定邊的權值Wij。同樣的,邊的權值確定也有兩種方法:①熱核法(Heat Kernel)。若第i個和第j個節點之間是連接的,則邊的權值為:Wij=exp(-‖xi-xj‖2/σ2);否則,Wij= 0。在本文中將采用熱核法確定權值。②簡單法。這種方法較簡單,如果第i個和第j個節點之間有邊連接,則定義邊的權值Wij= 1,否則Wij= 0 。③特征映射。假設前面所建立的近鄰圖是連通的,那么尋找低維嵌入的問題就歸結為對廣義特征向量問題的求解:

Ly=λDy

3基于LE算法的故障特征提取方法

滾動軸承故障診斷的過程是典型的模式識別的過程,對于流形學習算法這種無監督模式識別系統的構成如圖1所示。特征提取與選擇和分類識別是模式識別的核心問題。本文中將對振動信號先作預處理提取出多個時域和頻域特征量構建高維特征空間,繼而利用LE算法的非線性和降維的優點在特征空間中進行學習,進行特征融合,提取出最能準確、敏感表征軸承的運行狀態的低維流形特征量以便做出最有利的決策。特征提取是模式識別中最重要和最困難的環節,而特征空間構建的好壞將決定降維提取特征質量。

圖1 模式識別的典型過程 Fig.1 The typical process of pattern recognition

3.1構建原始高維特征空間

滾動軸承運行狀態識別和分類的關鍵是對表征狀態特征的提取[9]。由于現在的設備越來越復雜,所測得的振動信號所包含的信息量大,因而設備的運行狀況很難直接利用振動信號進行評估。通過特征提取技術把原始信號變換到特征空間,進而提取出準確描述設備運行狀態的信息。從原始振動信號中提取幅值、時域、頻域以及時頻域的特征參數已被廣泛使用,但是由于提取的每個特征參數對設備健康狀態的規律性、敏感性和聚類性各不相同且表征規律不一,所以提取出規律性強、敏感性好的狀態特征參數一直是設備故

障診斷領域所追求的,也是研究的熱點和難點。

為了提高對運行狀態和故障樣本的識別率,本文將從時域和頻域綜合進行分析。根據滾動軸承故障信號特點,將選擇最大峰值、絕對均值、均方根等13個能反映軸承故障的時域特征指標。針對時域特征量的局限性,在頻域特征中除了選擇平均頻率外還采用db4小波包函數進行3層正交小波包分解,得到均勻劃分的8個子頻帶的濾波信號,將得到的8個子頻帶對總頻率的相對能量比作為頻域統計特征參數[10],共計22種特征量作為滾動軸承的運行狀態原始特征量,如表1所示,詳細的計算公式和定義見文獻[11]。由于上述22個特征參數的規律性、敏感性以及聚類性各不相同,很難用其中一個或某幾個特征參數來準確地表征滾動軸承的運行狀態。因此,本文將以上述特征量作為訓練樣本的維度,使用降維方法對特征矩陣進行學習和特征融合,提取出敏感和穩定的軸承振動特征參數,有效、準確地描述滾動軸承運行狀態。

表1 時域和頻域特征參數表

3.2LE算法對原始特征空間進行學習

參考文獻根據LE算法的描述,局部鄰域的選取合適與否對低維流形提取效果至關重要。根據[8],在本文中選用K=8的K近鄰域法(KNN)構建鄰域圖,使用熱核法確定權值,進而進行局部線性投影和低維特征的提取。由于LE算法是以保持高維流形的局部近鄰點的信息為目標,即原空間中是近鄰點在低維空間中也應該是近鄰點,采用模式識別中聚類方法進行聚類識別。這樣LE方法就將問題轉化成對矩陣特征值的求解,過程就簡單化,不需要進行迭代計算,和其他幾種流形學習算法相比,計算量和運算時間大大減少。表2為三種方法在配置為2 G內存、CPU為2 GHz T5870的Thinkpad筆記本上運行Matlab R2008軟件實驗平臺上進行兩組模擬實驗所需的時間。

表2 計算時間

從上表可知,PCA所用時間最短,LE所花時間比KPCA少,這是由于LE算法計算過程簡單、執行速度快,而KPCA則是以付出大量計算為代價,所以比PCA和LE用時多。

4故障分類實驗驗證

本文中所采用的實驗數據均來自美國凱斯西儲大學(CWRU)實驗室所進行的滾動軸承實驗[12],實驗對象為SKF6205-2RS滾動軸承,將選擇滾動軸承不同故障類型和滾動體上不同受損程度進行兩組模擬實驗驗證。所有實驗數據是在下列條件下測得:電機負載功率為2 hp,轉速為1 750 r/min,數據采樣頻率為48 K。在故障類型識別試驗中,軸承內圈、外圈和滾動體上損傷小坑的尺寸為:直徑0.053 cm、深度0.028 cm。而在滾動體故障受損程度實驗中,損傷尺寸的直徑分別為0.018 cm、0.036 cm和0.053 cm,深度均為0.028 cm。在兩組對比實驗中每種狀態都選取50組樣本,構成N=50×4=200,m=22,N×m=200×22的特征矩陣。

4.1實驗I:故障類型實驗

在實驗中將選取滾動軸承的內圈、外圈和滾動體3種典型不同的故障位置及健康無故障狀態進行分析。

對4種不同狀態檢測都選擇同一條件下驅動端的軸承加速度值進行分析,目標維數選取為d=3。

從圖2和表3中可知,圖2(a)為PCA方法所得到的結果。除無故障類型的樣本很好的識別,聚類性比較好以外其他三種故障類型樣本之間出現了相互混疊,尤其外圈故障分布比較離散且識別率最低。總體來說識別率和聚類效果不理想。

圖2(b)為通過KPCA提取的前三個主元K1-K3的三維圖。從圖中可知,由于KPCA采用了高斯核函數,四種狀態樣本基本區分開,無故障的50個樣本數據實現了很好的識別和聚集,其它三種故障樣本的聚類性與PCA的聚類性相比有了較大的提高。但是,三種故障樣本的聚類性仍不太理想,出現了相互混疊現象,聚集性仍有待提高。

圖2(c)是采用流形學習中的LE算法對高維特征空間進行特征融合學習提取的前三個主分量L1-L3的結果。滾動軸承的四種樣本狀態在三維空間中呈現了很好的聚類性,較好地聚集在了同一點,并且類間距有明顯的分離便于區分。

圖2(d)和圖2(c)相輔相成,圖2(d)是由LE方法提取的前三個低維嵌入主分量L1-L3確定的4種狀態的表示。可知內圈故障樣本在圖2(c)中的中心坐標為(0.1543,0,0)、滾動體故障坐標為(0,0,0.1601)、外圈故障為(0,0,0)及無故障樣本坐標為(0,-0.1622,0)。與PCA和KPCA兩種方法相比,LE方法能對高維樣本空間提取出不同故障狀態下的有效信息,實現了對軸承故障狀態的有效識別。

圖2 不同的故障位置分類識別 Fig.2 Different fault location classification

4.2實驗II:滾動體故障不同受損程度實驗

本實驗選取深度均為0.028 cm,直徑分別為0.018 cm、0.036 cm和0.053 cm的滾動體故障,以及無故障滾動軸承四種類型進行對比驗證。

圖3中的(a)和(b)分別為滾動軸承的不同受損程度下的時域和頻域圖。在圖(a)中,至上而下分別是無故障、故障直徑為0.018 cm、0.036 cm和0.053 cm的時域圖,從圖中只能區分出無故障數據和有故障數據。結合故障頻譜圖(b),發現滾動體故障的頻率為137.5 Hz,但并不能實現對故障受損程度的識別。

圖3 不同的故障損傷 Fig.3 Different fault damage

表3為利用三種方法對多特征矩陣進行特征融合優化,采用K近鄰分類測試得到的樣本識別率,圖4為三種方法得到的降維結果。同樣,對于每一種類型也都選取50組樣本特征數據構造出N×m=200×22矩陣。目標維數同樣選取d=3。

表3 三種方法識別準確率

由圖4(a)可知,四種狀態沒有很明顯的區分,只有無故障樣本呈現了較好的聚類性,其它三種的聚類性表現并不理想,故障直徑為0.018 cm的樣本最分散、聚類性最差,且無故障、0.036 cm和0.053 cm故障的三種樣本之間的現對距離較近,出現了邊界和樣本混疊。

圖4(b)是KPCA方法降維的結果:無故障樣本呈現了很好的聚集;直徑為0.053 cm的樣本只有少部分被錯誤識別,也基本實現了識別并聚集在了一起;0.018 cm和0.036 cm的樣本之間沒有很明顯的界限出現了大量的混疊面。總體而言識別率和聚類性比PCA有了提高,仍有待提高。

圖4(c)是用LE方法對特征空間學習提取的前三個主分量得到的結果。滾動軸承的四種樣本狀態在三維空間中呈現了很好的聚類性,幾乎聚集在了同一點,并且類間距也非常明顯地分離開來。圖4(d)與圖4(c)是相輔相成,圖4(d)是經LE方法提取出來的前三個低

圖4 滾動體不同的受損程度分類識別 Fig.4 Different fault location classification

維嵌入主分量L1-L3確定的4種狀態的表示,可知從圖4(c)和圖4(d)兩幅圖中得知四種狀態在圖4(c)中的三維位置:0.018 cm的坐標為(0,-0.1414,0)、0.014的坐標為(0,0,0.1401)、0.053 cm的坐標為(0,0,0)以及無故障樣本坐標為(-0.1814,0,0)。較PCA和KPCA兩種方法,LE方法對高維特征空間進行學習,提取出了不同故障狀態下的有效信息,實現了對滾動體故障不同的受損程度的有效識別和分離。

4.3參數驗證

故障診斷是典型的模式識別問題,因此將采用類間距Sb和類內距Sw來衡量提取特征的聚類效果[13-15]。假定特征向量為{f1、f2、…fd},d是特征向量的目標維數,本文選取d=3。Sb和Sw兩個參數的描述如下:

由表中數據可知:LE方法的效果明顯好于PCA和KPCA方法。具體而言,PCA的效果最差;由于KPCA采用了核函數,聚類性比PCA有了很大提高,尤其是類內聚;LE的效果最好,Sb/ Sw的比值遠遠大于其他兩種方法的結果。

表4 類間距和類內聚的值

5測試樣本驗證

上一部分的實驗已經證明了LE算法對滾動軸承不同的故障類型以及滾動體不同受損程度有很好的診斷和識別。為進一步證明其有效性以便能應用到工業生產中,本節將選取另外50組實測特征樣本和作為訓練樣本原50組特征數據進行驗證實驗。

由于LE算法采用批量處理的方式處理數據,對新增樣本進行特征提取時,不具備增量處理能力,不得不重新建立低維流形結構。因此,在本節中將采用增量LE算法進行故障樣本的特征提取。根據參考文獻[16]對增量LE算法進行如下描述:

(1)先用LE算法對原始樣本集X進行學習,得到對應的低維嵌入坐標Y。

增量LE算法較好的保持了對新增樣本的領域局部映射關系,以鄰域內樣本的權值矩陣求得新樣本對應低維空間的位置坐標。當數據集實時更新的情況下,可以保證提取的特征信息能夠隨新樣本數據的增加而動態更新,從而有利于更新故障樣本的診斷分析。

5.1驗證實驗Ⅰ

此驗證實驗中將對滾動軸承的內圈、外圈、滾動體故障以及無故障四種狀態測試樣本與訓練樣本進行驗證,得到圖5和表5所示結果。可得到如下結論:PCA的線性本質使得各種狀態之間不能很好的分離,四種狀態間出現了混疊,測試樣本和訓練樣本也出現大量混疊,識別率不高。圖5(b)由于KPCA采用了核函數,各種類型間的類間距比較明顯,基本區分出了不同類型故障;滾動體故障和無故障的測試和訓練樣本基本重疊在一點,其他兩種也實現了較好的重疊,識別率方面有了很大的提高。LE的聚類性和識別率很好,幾乎聚集在了一點,很好的實現了測試和訓練樣本的重疊。

5.2驗證實驗Ⅱ

采用同樣的方法對軸承無故障、滾動體故障直徑三種不同受損程度下的樣本進行驗證識別。訓練樣本和測試樣本驗證得到的結果圖6和表5,與驗證實驗Ⅰ的結論類似,進一步證明了LE算法在軸承故障診斷中的有效性。

表5 三種方法識別準確率

圖5 軸承不同故障類型分類驗證 Fig.5 Different fault location classification

圖6 滾動體不同受損程度分類驗證 Fig.6 Different fault location classification

6結論

本文在流形學習的基礎上提出了基于LE算法的滾動軸承故障診斷模型。將LE的非線性降維能力應用于高維復雜樣本特征空間中,這種空間是由多個時域特征量和經小波包變換得到多個頻域特征量構成。LE算法進行自學習,進行特征融合,提取出了最能表征滾動軸承狀態的敏感特征量,實現了從高維復雜特征數據中的低維提取并用低維空間表示。LE算法不僅識別和診斷出滾動軸承的不同故障類型而且實現了對滾動體上不同故障受損程度的識別。使用類內距Sw和類間距Sb個參數來衡量聚類性的質量,以及使用樣本識別率證明了LE方法比PCA和KPCA優越。并通過采用測試樣本和訓練樣本進一步驗證其可使用性,為以后將LE算法以及流形學習中其他算法推廣到工業生產中進行故障診斷奠定了基礎。

參考文獻

[1]鐘秉林,黃仁.機械故障診斷學[M].北京:機械工業出版社,2006:298-299;171-180.

[2]Lee J M, Yoo C K. Nonlinear process monitoring using kernel principal component analysis [J], Chemical Engineering Science,2004(59):223-234.

[3]鄧曉剛,田學民.一種基于KPCA 的非線性故障診斷方法[J].山東大學學報,2005,35(3):103-106.

DENG Xiao-gang, TIAN Xue-min.Nonlinear process fault diagnosis method using kernel principal component analysis[J]. Journal of Shandong University,2005,35(3):103-106.

[4]蔣玲莉.基于核方法的旋轉機械故障診斷技術與模式分析方法研究[D].長沙:中南大學,2010.

[5]Seung H S, Daniel D L. The manifold ways of perception [J]. Science, 2000, 290 (5500): 2268-2269.

[6]Roweis S, Saul L. Nonlinear dimensionality reduction by locally linear embedding[J]. Science, 2000, 290(5500): 2323-2326.

[7]Tenenbaum J, Silva D D, Langford J. A global geometric framework for nonlinear dimensionality reduction [J]. Science,2000,290 (5500): 2319-2323.

[8]王澤杰,胡浩民.流形學習算法中的參數選擇問題研究[J].計算機應用與軟件,2010,27(6):84-85.

WANG Ze-jie, HU Hao-min. On parameter selection in manifold learning algorithm [J]. Computer Applications and Software,2010,27(6):84-85.

[9]何清波.多元統計分析在設備狀態監測診斷中的應用研究[D].合肥:中國科學技術大學,2007.

[10]賈茂林,王孫安,梁霖.利用非線性流形學習的軸承早期故障特征提取方法[J].西安交通大學學報,2010,44(5):45-49.

JIA Mao-lin,WANG Sun-an,LIANG Lin. Feature extraction for incipient fault diagnosis of rolling bearings based on nonlinear manifold learning[J].Journal of Xi’an Jiaotong University,2010,44(5):45-49.

[11]盛兆順,尹琦嶺.設備狀態監測與故障診斷技術及應用[M]. 北京:化學工業出版社, 1991:40-41.

[12]The Case Western Reserve University Bearing Data Center.http://csegroups.case.edu/bearingdatacenter/home.

[13]張學工.模式識別[M].北京:清華大學出版社, 2010,146.

[14]He Qing-bo. Vibration signal classification by wavelet packet energy flow manifold learning[J]. Journal of Sound and Vibration,2013,332:1881-1894.

[15]He Qing-bo. Time-frequency manifold for nonlinear feature extraction in machinery fault diagnosis[J]. Mechanical Systems and Signal Processing,2013,35:200-218.

[16]蔣全勝,李華榮,黃鵬. 一種基于非線性流形學習的故障特征提取模型[J]. 振動與沖擊,2012,31 (23):132-136.

JIANG Quan-sheng,LI Hua-rong,HUANG Peng. A fault feature extraction model based on nonlinear manifold learning[J]. Journal of Vibration and Shock,2012,31 (23):132-136.

猜你喜歡
模式識別特征提取
特征提取和最小二乘支持向量機的水下目標識別
基于Gazebo仿真環境的ORB特征提取與比對的研究
電子制作(2019年15期)2019-08-27 01:12:00
基于Daubechies(dbN)的飛行器音頻特征提取
電子制作(2018年19期)2018-11-14 02:37:08
紫地榆HPLC指紋圖譜建立及模式識別
中成藥(2018年2期)2018-05-09 07:19:52
淺談模式識別在圖像識別中的應用
電子測試(2017年23期)2017-04-04 05:06:50
Bagging RCSP腦電特征提取算法
第四屆亞洲模式識別會議
可拓模式識別算法中經典域的確定方法
第3屆亞洲模式識別會議
基于MED和循環域解調的多故障特征提取
主站蜘蛛池模板: www.狠狠| 亚洲欧美一区二区三区蜜芽| 欧美在线一二区| 亚洲第一页在线观看| 亚洲无码37.| 免费无码网站| 人人91人人澡人人妻人人爽 | 亚洲综合九九| 国产成人在线无码免费视频| 黄色福利在线| 国产视频大全| 在线一级毛片| 久久综合干| 在线观看亚洲国产| 看国产毛片| 乱人伦中文视频在线观看免费| 婷婷午夜天| 99久久国产综合精品2023| 亚洲va在线观看| 亚洲色图另类| 激情无码视频在线看| 国产精品不卡永久免费| 热九九精品| 亚洲aaa视频| 免费一级毛片在线播放傲雪网| 亚洲精品无码在线播放网站| 久久久久亚洲av成人网人人软件| 全部免费特黄特色大片视频| 狂欢视频在线观看不卡| 国产a v无码专区亚洲av| 高清精品美女在线播放| 国产玖玖玖精品视频| 热re99久久精品国99热| 亚洲精品爱草草视频在线| 乱系列中文字幕在线视频| 国产亚洲精久久久久久久91| 91无码人妻精品一区| 日本一区二区不卡视频| 亚洲国产精品无码AV| 伊人成人在线| 狼友视频一区二区三区| 亚洲一级无毛片无码在线免费视频 | 午夜一区二区三区| 精品成人一区二区| 国产丝袜一区二区三区视频免下载| 亚洲无码不卡网| 国产免费网址| 久草国产在线观看| 日本五区在线不卡精品| 久久国产精品麻豆系列| 一级毛片无毒不卡直接观看| 2024av在线无码中文最新| 精品视频一区在线观看| 日韩AV无码一区| 91精品国产自产91精品资源| a色毛片免费视频| 亚洲毛片在线看| 国产成人综合欧美精品久久| 日韩精品成人网页视频在线| 97超级碰碰碰碰精品| 欧美一区福利| 亚洲精品波多野结衣| 日本中文字幕久久网站| 色哟哟国产成人精品| 国产精品女人呻吟在线观看| 亚洲精品图区| 丁香婷婷激情综合激情| 免费中文字幕一级毛片| 亚洲精品日产AⅤ| 一区二区自拍| 亚洲成人精品久久| 国产精品大白天新婚身材| 免费观看精品视频999| 人妻精品全国免费视频| 久久天天躁狠狠躁夜夜躁| 中国一级毛片免费观看| 亚洲人成色77777在线观看| 成人伊人色一区二区三区| 777国产精品永久免费观看| 婷婷色在线视频| 美女被操黄色视频网站| 午夜毛片福利|