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

基于多尺度關聯維數和流形學習的自動機故障診斷

2017-11-07 02:34:04房立清齊子元
中國測試 2017年10期
關鍵詞:故障診斷關聯振動

杜 偉,房立清,呂 巖,齊子元

(軍械工程學院火炮工程系,河北 石家莊 050003)

基于多尺度關聯維數和流形學習的自動機故障診斷

杜 偉,房立清,呂 巖,齊子元

(軍械工程學院火炮工程系,河北 石家莊 050003)

針對自動機振動信號非平穩、非線性的特點,提出基于多尺度關聯維數和線性局部切空間排列(linear local tangent space alignment,LLTSA)相結合的自動機故障診斷方法。首先,利用局部特征尺度分解(local characteristicscale decomposition,LCD)將自動機振動信號分解為不同尺度下的內稟尺度分量(intrinsic scale component),提取出反映狀態信息的主要分量并計算各分量的關聯維數。然后,利用線性局部切空間排列算法挖掘出可區分度更高的特征子集。最后,將得到的低維特征輸入支持向量機進行識別,自動機故障診斷實驗表明,所提方法具備較高的診斷準確率。此外,將LCD與經驗模態分解(empirical mode decomposition,EMD)和局部均值分解(local mean decomposition,LMD)方法的診斷結果進行比較,驗證所提方法的優勢。

故障診斷;多尺度關聯維數;流形學習;自動機

0 引言

自動機是高炮武器系統故障率最高的部分之一,零部件較多且大部分安裝在狹小的空間內。自動機在高溫、高壓的火藥氣體作用下完成開鎖和閉鎖、開閂和關閂、抽筒和拋筒、擊發等動作,其往復運動過程經常伴隨著激烈的撞擊、摩擦、振動和跳動等[1],再加上惡劣工作環境的影響,其組件容易出現磨損、失效等故障狀態。

傳統的自動機維護檢測方式如“看、摸、聽”和開箱解體,維修周期長、成本高。利用性能參數、頻譜分析、鐵譜分析等方法對故障進行定位需要復雜的濾波和頻譜分析設備。戴涌等[2]采用故障樹診斷方法分析了某自動機異常發射的故障模式,但利用故障樹進行故障診斷需要以大量的經驗知識和數據為基礎,并進行反復試驗和逐一排查。隨著信號處理技術和模式識別技術的發展,基于振動信號測量與分析的故障診斷方法得到了廣泛應用。自動機工作時產生的非線性、非平穩的振動信號含有豐富的工況信息,并在一定尺度范圍內具有分形特性。許多學者應用分形理論在非線性行為的定量描述中作了許多探索和研究,如蘭海龍等[3]將局域波分解和分形理論相結合,通過求解各模式分量的關聯維數,有效提取了自動機的狀態特征。Zhang等[4]利用提升小波分解對振動信號進行分解和重構,然后以時頻域能量的形態學分形維數為特征,準確地區分了故障類別。因此將自動機振動信號量化為關聯維數進行故障診斷是可行的。

經查閱文獻可知,噪聲的存在對關聯維數計算結果影響較大[5]。當信號的非線性特性突出時,可通過對信號進行多尺度自適應分解并剔除冗余的單分量來達到降噪的目的。楊宇等[6]用內稟時間尺度分解改進算法對振動信號進行分解,選取包含主要信息的分量重新組合,以此作為降噪后的信號并計算關聯維數,故障診斷結果表明降噪后的關聯維數區分性更高。李琳等[7]對振動信號進行經驗模態分解(EMD)[8],根據固有模態分量與原信號的相關性選取主要分量,取得了較好的降噪效果。局部特征尺度分解(LCD)算法[9]是一種新的自適應信號分解方法,能夠將非平穩信號分解成一系列內稟尺度分量(ISC),表征信號在不同尺度下隨時間變化的局部特征。程軍圣等[10]對振動信號的ISC分量做包絡譜分析,有效提取了故障信息,并用仿真信號驗證了LCD在運算速度、抑制端點效應和模態混疊方面較EMD和局部均值分解(LMD)[11-12]方法性能更優。由于LCD中提取到的特征集中依然存在一定數量的混疊信息,為了進一步消除非敏感特征,可采用流形學習進行二次特征提取。線性局部切空間排列(LLTSA)是一種新型流形學習算法[13],能夠充分挖掘高維非線性特征集的本質結構,進一步提高特征的可辨識性。張良等[14]提取出LCD各分量的模糊熵,然后采用LLTSA得到了敏感度更高的低維特征,輸入支持向量機進行故障分類,取得了較高的識別精度。以上研究都為自動機故障診斷提供了新的思路。

基于以上分析,提出了基于多尺度關聯維數和流形學習的故障診斷方法。結合LCD多尺度分析能力和關聯維數對非線性程度反應敏感的特點,提取自動機振動信號多個尺度的故障特征,而后采用 LLTSA得到維數低、聚類性好的故障特征集,并采用支持向量機判斷故障類型。自動機的故障實驗表明,所提方法能提高故障診斷的準確率。

1 多尺度關聯維數與LLTSA

1.1 分形與關聯維數

相比于其他分形維數如盒子維數、信息維數和容積維數等,關聯維數的計算最為簡單,能夠敏感反映系統吸引子的不均勻度,較好地識別出系統的工作狀態。

基于延時嵌入相空間重構思想的G-P算法[15],得到一維時間序列在嵌入維數為m時的關聯函數Cm(r)。Cm(r)反映了重構的相空間內吸引子上距離小于r(r>0)的點對所占的比例,在r的無標度區內滿足 Cm(r)=rD2(m,r),當 r→0 時,可以得到關聯維數:

畫出雙對數曲線lnCm(r)-lnr,采用點間斜率遠離均值剔除算法[16]得到無標度區,對其進行最小二乘線性擬合,計算其斜率即可得到該時間序列的關聯維數。

嵌入維數m和延遲時間 的選取對關聯維數計算影響較大。參數 確定方法有自相關函數法、C-C法和互信息法[17]等。參數m的確定方法有關聯積分法、觀察法和Cao法[18]等,這里分別選用互信息法和Cao法來確定參數 和m。

1.2 局部特征尺度分解

假設一個多分量信號x(t)可被分解成有限個瞬時頻率具有物理意義的內稟尺度分量(ISC)之和,且任意兩個ISC分量相互獨立,根據所定義的ISC分量,對信號x(t)進行LCD分解,其分解過程如下[9]:

1)確定 x(t)的所有極值點(k,Xk),k=1,2,3,…,M,其中M為極值點的個數。連續相鄰的兩個極值點可以將x(t)分割成若干個區間,在任意兩個相鄰極值點間對x(t)進行線性變換:

式中Lt(k)表示由原信號的第k個區間的線性變換所得到的基線信號段,將相鄰區間內的基線信號段首尾相連即得基線信號Lt,式(2)中:

其中,參數a一般取值為0.5。

2)將基線信號Lt從原始信號x(t)中分離出來,得到剩余信號P1(t)。如果P1(t)滿足ISC分量判據,則令 ISC1(t)=P1(t)。 否則將 P1(t)作為原始信號重復步驟1)、步驟2),循環k次,直到得到內稟尺度分量Pk,即 ISC1(t)。

3)將 ISC1(t)從原始信號 x(t)中分離出來,重復步驟1)、步驟2),重復循環n次,得到n個滿足ISC分量判據的分量,直到殘余分量rn為一單調函數或者小于預設閾值為止。于是:

LCD方法從原始信號中逐步分離出ISC分量,實現信號數據的自適應多尺度化,得到信號中不同層次的信息。通過剔除殘余分量,可達到降噪的目的。

1.3 線性局部切空間排列算法

LLTSA是一種非線性維數約簡方法,通過構建樣本點鄰域的低維切空間并進行全局排列,得到樣本點的低維全局坐標。即尋找一個轉換矩陣A,將RD空間中具有N個點的含噪數據集XORG(故障樣本集)映射為 Rd空間數據集 Y=[y1,…,yN],即:

其中HN=I-eeT/N為中心矩陣,I為單位矩陣,e為k維全1向量。Y為XORG潛在的d維非線性流行。包含以下3個步驟[13]:

1)構建鄰域。采用K-近鄰法(KNN)得到每個數據樣本點 xi(i=1,…,N)的鄰域 Xi=[xi1,…,xik],k 為鄰近點個數。

2)獲取局部信息。尋找一組正交基,提取Xi的局部低維坐標Θi,正交基的求取過程相當于在Xi上進行主成分分析(principle component analysis,PCA)。

3)局部切空間全局排列。局部切空間全局排列的目的是重構數據集的本征結構,使得將所有樣本點xi的局部切空間映射到全局低維坐標的誤差之和最小,即如下目標函數:

式中:Yi為Xi的全局低維坐標,Li為映射矩陣,且當Li=YiHkΘi+時重構誤差最小,是 Θi的 Moore-Penrose廣義逆。式(6)的最小化問題可以進一步轉化為式(7)的廣義特征值求解:

式中:B=SWWTST為全局排列矩陣,S=[S1,S2,…,SN],Si為 0~1 之間的選擇向量,W=diag(W1,W2,…,WN),其中。式(7)中最小的 d 個非零廣義特征值對應的特征向量 α1,α2,…,αd組成的矩陣ALLTSA=[α1,α2,…,αd]即為從高維空間到低維流形子空間的映射矩陣,最后求得低維空間坐標。

2 基于多尺度關聯維數和流形學習的診斷方法

該方法利用LCD在多個尺度內提取信號的關聯維數,結合流形學習算法對非線性特征集的約簡能力,獲得聚類性好的低維特征子集,并通過支持向量機判斷故障類型,實現故障診斷。該故障診斷方法的主要流程如圖1所示,具體步驟如下:

1)對原始振動信號進行LCD分解,獲得一系列ISC分量。

2)不同數據樣本被分解為若干ISC分量和1個殘余分量,剔除殘余分量,計算其余D個有用分量的關聯維數,組成k維的故障特征向量。

3)將訓練樣本和測試樣本同時輸入LLTSA進行維數約簡,經多次試驗確定目標維數d和鄰域參數k,其中 1≤d<D。

圖1 基于多尺度關聯維數和流形學習的故障診斷流程

4)將訓練樣本的d維特征集用于SVM訓練,使用已訓練的SVM對測試樣本進行分類識別,確定故障類型。

3 實驗分析

3.1 實驗數據采集

實測振動信號來自某型高炮導氣式自動機,自動機藥筒與身管以及炮閂組件之間的相互作用構成了自動機振動的主要激勵源。閉鎖機構在射擊過程中用于實現閉鎖炮膛,閉鎖塊在高溫、高壓和高沖擊的惡劣工作環境下易出現磨損、點蝕和裂紋等故障。輸彈機構中輸彈簧的斷裂或疲勞也會對自動機的工作狀態產生較大影響。根據自動機常見故障模式,設置閉鎖塊磨損(F1)、輸彈簧失效(F2)和正常(N)3種工作狀態。其中,閉鎖塊上設置面積約為7mm×10mm的磨損故障,如圖2所示;輸彈簧失效方式為塑性變形。考慮實際應用的便捷性和實驗操作的安全性,采用手氣開閂的方式使自動機完成關閂、閉鎖、開鎖、開閂等過程。

圖2 閉鎖塊磨損

圖3 傳感器安裝位置

分別測取3種狀態工作過程的振動信號,加速度傳感器安裝在相對空曠、光滑的箱體表面,如圖3所示。采樣頻率為20kHz,每種狀態測取40組數據,3種狀態共120組。為突出故障特征,減小實驗操作中無關信號的干擾,采用分時段分析的方法,即只截取閉鎖塊工作過程內的信號進行分析,從閉鎖完畢時閉鎖塊與支撐塊的碰撞開始,按時間順序以1600個采樣點為一組數據樣本,自動機3種運行狀態的振動信號如圖4所示。從圖中可以看出,各個狀態振動信號的時域波形差異并不明顯,因此,僅憑時域波形判斷是不準確的,需要進一步分析識別。

3.2 LCD關聯維數特征提取

采用G-P算法直接計算各狀態原始振動信號的關聯維數,其中,分別由互信息法和Cao法確定正常狀態和閉鎖塊磨損狀態的時間延遲 =3,嵌入維數m=10,輸彈簧失效狀態的時間延遲 =6,嵌入維數m=10。限于篇幅,表1給出了每種狀態前5組數據的計算結果,以及全部40組關聯維數的均值。可以看出正常狀態下的關聯維數較小,相對于故障狀態區分較為明顯;閉鎖塊磨損和輸彈簧失效故障的關聯維數較大,且差別較小,不易區分。原因可能為閉鎖塊的磨損和輸彈簧的失效導致運動產生故障性沖擊成分,系統的非線性響應增強,反映非線性系統復雜度的關聯維數相應增大,且振動信號的形態具有一定的相似性,故關聯維數比較接近。因此,僅在1個尺度計算關聯維數值的大小難以判別故障類型,需進一步挖掘信號中的深層次信息。

表1 原始信號關聯維數

圖4 自動機3種狀態時域波形

圖5 閉鎖塊磨損故障振動信號分解

對數據樣本進行LCD分解,多數數據樣本被分解為7個ISC分量和1個殘余分量。圖5為閉鎖塊磨損故障振動信號中一組樣本數據的分解結果。計算7個ISC分量的關聯維數,共得到3個40×7維的特征矩陣。圖6為自動機3種狀態各個ISC分量的關聯維數誤差線,計算結果為全部數據的統計。可以看出,部分分量所計算出的關聯維數具有一定的可分性,但在正常狀態和閉鎖塊磨損狀態的第一個ISC分量處、兩種故障狀態的第3個ISC分量處,誤差線出現重合,數據區分度有所下降;不同狀態曲線的排序趨勢在第4個分量之后發生變化,并且第5和第6個ISC分量的誤差線重合較嚴重,數據點的區分度下降,說明存在一定的冗余信息。據此,選取全部分量用于關聯維數計算,提高特征提取的全面性。

圖6 ISC分量關聯維數

圖7 SVM-LLTSA識別結果

3.3 LLTSA降維

將特征集輸入LLTSA中進行維數約簡,作為比較,分別將LMD和EMD分解得到關聯維數特征同時利用LLTSA進行降維處理,并將得到的低維特征集輸入SVM進行識別。隨機抽取20組數據作為訓練樣本,其余20組作為測試樣本。其中SVM的核函數選用性能較優的徑向基核函數(RBF),懲罰參數C=1,核函數參數g=1。圖7所示為不同分解方法的識別率隨目標維數d(范圍[1,6])和鄰域參數k(范圍[1,38])的變化趨勢。

圖7(a)為LCD關聯維數的識別結果,當目標維數和鄰域參數都小于2時,識別率均不足52%,隨著目標維數從2至5逐漸增大,故障識別率逐漸升高,由60%增至80%以上,并在d=5、k=31處達到了最大識別率88.33%,當d增至6時,識別率有所下降;從圖7(b)中可以看出,LMD關聯維數的識別率絕大部分在65%以上,在d=5處,不同鄰域參數值下的識別率相差較小,且在k=36處有最大識別率78.33%;圖7(c)為EMD方法的識別結果,當目標維數為1至3時,整體識別率較低,當目標維數增至5時,識別率隨不同鄰域參數值的變化波動較大,在k=9處達到最大識別率78.33%,繼續增大目標維數,識別率呈下降趨勢。比較圖7可以看出,當維數降為5維時,提取到的特征集聚類效果最佳,且EMD特征參量的最高識別率較低,LMD方法的識別率較高且變化范圍小,LCD關聯維數的最高識別率較前兩種方法提升了10%,識別率達到最高。

為進一步突出所提方法的優勢,將原始信號的關聯維數(None)、未經約簡的EMD、LMD和LCD關聯維數以及采用PCA約簡的LCD關聯維數特征輸入SVM進行分類識別,結果如表2所示。

表2 支持向量機識別結果

由表可以看出,輸彈簧故障振動信號相對于另外兩種狀態的識別效果較好,識別率達到了100%;由于原始信號關聯維數僅在1個尺度進行計算,忽略了信號中的深層次信息,使分析具有很大的局限性,因此故障識別率偏低,而利用信號自適應分解并提取多尺度特征參量的方法使故障識別率提高5%以上,也進一步證明了進行多尺度分析的必要性;LCD方法的自適應分解性能較EMD和LMD更優,所以提取的特征具有更高的識別率;由于LCD方法提取的特征中包含一定的冗余信息,對SVM的識別產生了影響,所以經維數約簡后的識別率比未約簡的識別率有明顯的提高;由于PCA屬于線性降維方法,忽略了樣本數據的非線性流形本質結構,降維效果不理想,因此用LLTSA方法約簡的結果比用PCA提高約8.33%。

4 結束語

1)針對非線性、非平穩振動信號提出了基于多尺度關聯維數和流形學習的故障診斷方法。實測信號分析表明,將振動信號的關聯維數組成特征集可對狀態識別提供良好的依據,而僅從1個尺度或部分尺度提取特征進行診斷存在一定的局限性。

2)利用LCD的自適應分解性能,提取到的多尺度關聯維數能夠從各個尺度反映信號的復雜性,將其輸入SVM進行故障診斷具有較高的識別率,且由LCD與EMD、LMD方法的對比實驗表明,LCD方法得到的特征集具有更優的識別效果。

3)利用LLTSA對非線性振動信號進行維數約簡,通過多次調整目標維數和鄰域參數得出最大識別率,實驗結果表明,所提方法比結合PCA降維的故障診斷方法具有更高的識別率,可有效用于自動機故障診斷。

[1]潘宏俠,都衡,馬春茂.局域波信息熵在高速自動機故障診斷中的應用[J].振動、測試與診斷,2015(6):1159-1164,1205.

[2]戴涌,張國平,王茂林,等.某轉膛自動機異常發射故障分析[J].火炮發射與控制學報,2013(3):67-71.

[3]蘭海龍,潘宏俠,龔明.局域波分形方法在柴油機故障診斷中的應用[J].機床與液壓,2013(5):180-183.

[4]ZHANG Z, WU J, MA J, et al.Fault diagnosis for rolling bearing based on lifting wavelet and morphological fractal dimension[C]∥Control and Decision Conference.IEEE,2015.

[5]ARGYRIS J, ANDREADIS I, PAVLOS G, et al.The influence ofnoise on the correlation dimension of chaotic attractors[J].Chaos Solitons and Fractals,1998,9(3):343-361.

[6]楊宇,王歡歡,喻鎮濤,等.基于ITD改進算法和關聯維數的轉子故障診斷方法[J].振動與沖擊,2012(23):67-70,76.

[7] 李琳,張永祥,明廷濤.EMD降噪的關聯維數在齒輪故障診斷中的應用研究[J].振動與沖擊,2009,28(4):145-148.

[8]HUANG N E, SHEN Z, LONG S R.A new view of nonlinear water waves:the Hilbert spectrum[J].Annual Review of Fluid Mechanics,1999(31):417-457.

[9]程軍圣,鄭近德,楊宇.一種新的非平穩信號分析方法:局部特征尺度分解法[J].振動工程學報,2012(2):215-220.

[10]程軍圣,楊怡,楊宇.局部特征尺度分解方法及其在齒輪故障診斷中的應用[J].機械工程學報,2012(9):64-71.

[11]SMITH J S.The local mean decomposition and its application to EEG perception data[J].Journal of the Royal Society Interface,2005,2(5):443-454.

[12]LEI Y G, HE Z J, ZI Y Y.Application of the EEMD method to rot or fault diagnosis of rotating machinery[J].Mechanical Systems and Signal Processing,2009,23(4):1327-1338.

[13]李鋒,湯寶平,陳法法.基于線性局部切空間排列維數化簡的故障診斷[J].振動與沖擊,2012(13):36-40.

[14]張良,張前圖.基于LCD模糊熵和流行學習的故障特征提取方法[J].機械強度,2016(2):225-230.

[15]GRASSBERGER P,PROCACCIA I.Characterization of strange attractors[J].Physical Review Letters,1983,50(5):346-349.

[16]李奎為,胡瑾秋,張來斌,等.關聯維數無標度區確定方法及診斷應用[J].石油機械,2007(4):43-45.

[17]PENG J E, QIAO H, XU Z B.A new approach to stability of neural networks with time-varying delays[J].Neural Networks,2002,15(1):95-103.

[18]CAO L Y.Practical method for determining the minimum embedding dimension of a scalar time series[J].Physica.Section D: Nonlinear Phenomena,1997,110(1/2):43-50.

Fault diagnosis of automaton based on multiscale correlation dimension and manifold learning

DU Wei, FANG Liqing, Lü Yan, QI Ziyuan
(Department of Artillery Engineering,Ordnance Engineering College,Shijiazhuang 050003,China)

Aiming at the non-stationary and non-linear characteristics of automaton vibration signal,a fault diagnosis method of automaton based on multiscale correlation dimension and linear local tangent space alignment (LLTSA)was proposed.Firstly, the vibration signal of automaton was decomposed with local characteristic-scale decomposition (LCD) to obtain intrinsic scale componentsin differentscales, and the correlation dimension ofeach principalcomponent reflected the state information was calculated.Then,the mining performance of the feature subset with higher distinguishability was furtherimplemented by using linear localtangentspace alignment.Finally,low-dimensional feature was put into SVM to recognize the state types.The results of automaton fault test indicate that the proposed method is of high accuracy.In addition,the diagnostic resultscalculated with LCD, empiricalmode decomposition and localmean decomposition were compared,verifying the advantage of the method.

fault diagnosis; multiscale correlation dimension; manifold learning; automaton

A

1674-5124(2017)10-0102-07

10.11857/j.issn.1674-5124.2017.10.020

2017-01-12;

2017-02-08

河北省自然科學基金資助項目(E2016506003)

杜 偉(1992-),男,山東臨沂市人,碩士研究生,專業方向為兵器性能檢測與故障診斷。

(編輯:劉楊)

猜你喜歡
故障診斷關聯振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
奇趣搭配
中立型Emden-Fowler微分方程的振動性
智趣
讀者(2017年5期)2017-02-15 18:04:18
因果圖定性分析法及其在故障診斷中的應用
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
基于LCD和排列熵的滾動軸承故障診斷
基于WPD-HHT的滾動軸承故障診斷
機械與電子(2014年1期)2014-02-28 02:07:31
主站蜘蛛池模板: 一本色道久久88| 国产尹人香蕉综合在线电影| 亚洲区视频在线观看| 精品少妇人妻av无码久久| 欧美日韩国产系列在线观看| 欧美综合一区二区三区| 99伊人精品| 亚洲国产成人超福利久久精品| 91无码人妻精品一区二区蜜桃| 国产成人AV男人的天堂| 19国产精品麻豆免费观看| 在线国产毛片手机小视频| 五月综合色婷婷| 999国内精品久久免费视频| 国内毛片视频| 欧美、日韩、国产综合一区| 国产精品一区在线麻豆| 欧美日韩一区二区三| 五月婷婷欧美| 狠狠ⅴ日韩v欧美v天堂| 91偷拍一区| 国产黄网永久免费| 午夜国产精品视频黄| 成人精品区| 影音先锋亚洲无码| 妇女自拍偷自拍亚洲精品| 亚洲天堂免费观看| 狠狠操夜夜爽| 亚洲精品无码在线播放网站| 高清精品美女在线播放| 丰满人妻一区二区三区视频| 欧美精品啪啪一区二区三区| 91九色国产porny| 最新国产精品鲁鲁免费视频| 精品国产污污免费网站| 2019年国产精品自拍不卡| 一级成人欧美一区在线观看| 国产精品人成在线播放| 中文毛片无遮挡播放免费| 中文国产成人精品久久一| 天堂岛国av无码免费无禁网站| 亚洲av无码人妻| 久久福利网| 日韩在线第三页| 亚洲日本在线免费观看| 99福利视频导航| 99久久99这里只有免费的精品| www中文字幕在线观看| 国产日韩精品一区在线不卡| 精品视频一区在线观看| 伊伊人成亚洲综合人网7777| 人妻21p大胆| 成年人午夜免费视频| 精品人妻一区无码视频| 欧美另类视频一区二区三区| 天天做天天爱天天爽综合区| 午夜一级做a爰片久久毛片| 日韩东京热无码人妻| аv天堂最新中文在线| 成人亚洲视频| 国产噜噜在线视频观看| 9cao视频精品| 亚洲国产高清精品线久久| 国产亚洲精久久久久久久91| 亚洲天堂精品在线| 亚洲人成色77777在线观看| 色欲综合久久中文字幕网| 午夜毛片福利| 9999在线视频| 精品人妻系列无码专区久久| 国产成人高清精品免费软件| 免费啪啪网址| 午夜国产不卡在线观看视频| 91破解版在线亚洲| 97亚洲色综久久精品| 成人国产三级在线播放| 中文字幕va| 97亚洲色综久久精品| 色妺妺在线视频喷水| 91小视频版在线观看www| 亚洲成人一区二区| 91美女在线|