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

滾動(dòng)軸承振動(dòng)的非線性超混沌特性研究

2016-07-25 08:18:28李兆飛任小洪黃臣程
軸承 2016年7期
關(guān)鍵詞:振動(dòng)故障信號(hào)

李兆飛,任小洪,黃臣程

(1.四川理工學(xué)院 自動(dòng)化與電子信息學(xué)院,四川 自貢 643000;2.人工智能四川省重點(diǎn)實(shí)驗(yàn)室,四川 自貢 643000)

研究滾動(dòng)軸承振動(dòng)的非線性演化動(dòng)態(tài)特性有助于認(rèn)識(shí)其內(nèi)在變化規(guī)律,并利用這些規(guī)律對(duì)其進(jìn)行診斷和監(jiān)測,具有一定的普遍意義。隨著對(duì)軸承振動(dòng)非線性機(jī)理的認(rèn)識(shí)逐漸加深,采用關(guān)聯(lián)維數(shù)[1]、最大Lyapunov指數(shù)[2-3]、Kolmogorov熵及相關(guān)復(fù)雜性測度[4-6]等方法對(duì)振動(dòng)的非線性混沌特性進(jìn)行判斷,很好地解釋了軸承振動(dòng)現(xiàn)象的復(fù)雜運(yùn)動(dòng)現(xiàn)象,判定軸承振動(dòng)是一種在有限相空間軌道永不重復(fù)、形態(tài)復(fù)雜的運(yùn)動(dòng)。然而,這些特征量只能表現(xiàn)軸承振動(dòng)的某方面特性,而且受信號(hào)長度及噪聲影響,這些特征量的計(jì)算結(jié)果會(huì)有所偏差,導(dǎo)致判斷錯(cuò)誤。因此,嘗試采用遞歸圖、CLY方法和功率譜3種方法對(duì)軸承不同振動(dòng)狀態(tài)實(shí)測信號(hào)的混沌特性進(jìn)行試驗(yàn)分析。

如前所述,目前對(duì)軸承振動(dòng)的混沌特性研究多是基于混沌的辨別,但對(duì)其振動(dòng)超混沌特性的分析還未見報(bào)道,分析軸承振動(dòng)的超混沌特性有助于研究振動(dòng)在相空間多個(gè)不同方向的復(fù)雜演化規(guī)律。通常通過計(jì)算實(shí)測滾動(dòng)軸承振動(dòng)信號(hào)的Lyapunov指數(shù)譜判斷其超混沌特性。Lyapunov指數(shù)譜的計(jì)算往往要通過重構(gòu)相空間,之后再判定未知系統(tǒng)方程狀態(tài)。判斷未知系統(tǒng)方程通常采用人工神經(jīng)網(wǎng)絡(luò)[7-9],存在訓(xùn)練時(shí)間較長、易陷入局部極小點(diǎn)及過擬合,且確定網(wǎng)絡(luò)結(jié)構(gòu)困難的問題[10]。因此,采用訓(xùn)練耗時(shí)較小且較少出現(xiàn)過擬合現(xiàn)象的最小二乘支持向量機(jī)[11](Least Squares Support Vector Maehine,LS-SVM)判別系統(tǒng)方程計(jì)算信號(hào)Lyapunov指數(shù)譜的方法,分析軸承振動(dòng)的超混沌特性。

1 軸承振動(dòng)信號(hào)采集及相空間重構(gòu)

1.1 軸承振動(dòng)信號(hào)的采集

采用Case West Reserve University滾動(dòng)軸承實(shí)驗(yàn)室實(shí)測的內(nèi)、外圈和滾動(dòng)體的不同故障程度與轉(zhuǎn)速振動(dòng)數(shù)據(jù)庫,該數(shù)據(jù)庫采用深溝球軸承6205-2RS JEM SKF,電動(dòng)機(jī)空載,采樣頻率fs=12 kHz,轉(zhuǎn)速為1 797 r/min,軸承故障由電火花加工機(jī)在球軸承內(nèi)、外圈及鋼球上模擬損傷性故障,故障直徑分別為0.18,0.36,0.54和0.72 mm(文中選取內(nèi)圈模擬點(diǎn)蝕故障并在后續(xù)分析中分別對(duì)應(yīng)簡化表示為f1~f4),故障深度為0.28 mm。正常狀態(tài)和內(nèi)圈故障振動(dòng)信號(hào)的時(shí)域波形如圖1所示,并分別取1 024個(gè)點(diǎn)研究不同故障類型及不同程度故障軸承振動(dòng)的混沌特性。

(b)內(nèi)圈故障

1.2 相空間重構(gòu)技術(shù)

相空間重構(gòu)技術(shù)是混沌特性研究的前提。混沌吸引子是混沌系統(tǒng)的特征之一,體現(xiàn)了混沌系統(tǒng)的運(yùn)動(dòng)規(guī)律,相空間重構(gòu)目的就是為了在高維空間中恢復(fù)混沌吸引子。系統(tǒng)運(yùn)動(dòng)是有n個(gè)變量的動(dòng)力系統(tǒng),即一組n個(gè)變量的一階微分方程

(1)

由連續(xù)變量坐標(biāo)x(t)及其(n-1)階導(dǎo)數(shù)構(gòu)成的n維相空間就表示了系統(tǒng)隨時(shí)間演變的狀態(tài)空間,即

X(t)=[x(t),x1(t),…,xn-1(t)]T。

(2)

Ruelle用離散序列x(t)及其(n-1)時(shí)滯位移將其代替,構(gòu)成一個(gè)新的n維嵌入相空間,即

X(t)= [x(t),x(t+τ),…,x(t+(m-1)τ)]T

,(3)

式中:τ為時(shí)延;m為重構(gòu)相空間的維數(shù),m≥2d+1(d為狀態(tài)空間的關(guān)聯(lián)維數(shù),也稱拓?fù)渚S)。對(duì)xn(n=1,2,…,N),取定τ,其重構(gòu)相空間為

Xl= [xl,xl+τ,…,xl+(m-1)τ]T=

(4)

式中:l=N-(m-1)τ。求取τ使原序列時(shí)延后能作為獨(dú)立坐標(biāo),可根據(jù)平均位移[14]、自相關(guān)[13]、互信息[12]及(去偏)復(fù)自相關(guān)技術(shù)等計(jì)算。好的重構(gòu)相空間是使重構(gòu)后的吸引子與系統(tǒng)真正的吸引子拓?fù)涞葍r(jià)。m是嵌入維,可根據(jù)預(yù)測誤差最小[16]、幾何不變量[14]、虛假鄰近點(diǎn)[17]及其改進(jìn)后的Cao方法[18]計(jì)算。確定時(shí)延與嵌入維有2種爭議:一種認(rèn)為可先求得時(shí)延再選擇嵌入維;另一種則認(rèn)為時(shí)延與嵌入維是相關(guān)的。研究者提出較常用的C-C[19]及其改進(jìn)法[20-21]能同時(shí)計(jì)算時(shí)延與嵌入窗。因?yàn)椴荒艿玫交煦缦到y(tǒng)先驗(yàn)知識(shí),這些方法均有一定的主觀性[22]。

分析中采用自相關(guān)法[13]求τ,根據(jù)xn與xn+τ有些獨(dú)立又不完全不相關(guān)的特性,為使其能在重構(gòu)相空間中作為獨(dú)立坐標(biāo)處理,計(jì)算自相關(guān)函數(shù)為

(5)

則C(τ)首次下降到初值的1-1/e時(shí)的延時(shí)為最優(yōu)τ,不同嵌入維下正常及內(nèi)圈故障振動(dòng)信號(hào)的自相關(guān)函數(shù)分別如圖2所示,從圖中可以看出,正常狀態(tài)最佳延時(shí)為3,不同內(nèi)圈故障狀態(tài)的最佳延時(shí)均為2。

(b)內(nèi)圈故障

再用虛假鄰近點(diǎn)法[17]計(jì)算m,先選定m,求得Xl的最近鄰點(diǎn)Xη(l),即

‖Xl-Xη(l)‖m=min{‖Xl-Xj‖:j=1,2,…,l,l≠j},

(6)

將m增加到m+1,當(dāng)滿足

(7)

Xη(l)即為Xl的虛假鄰近點(diǎn)。(7)式說明當(dāng)相空間從m維演變到m+1維時(shí),2個(gè)相點(diǎn)相差很大距離是因?yàn)楦呔S混沌吸引子中2個(gè)不相鄰的點(diǎn)投影到低維軌道上時(shí)變成相鄰的2個(gè)臨近虛假點(diǎn)造成的。求取最佳嵌入維數(shù)m時(shí)一般取RT=10[17],統(tǒng)計(jì)每個(gè)m時(shí)其虛假鄰近點(diǎn)占所有重構(gòu)向量數(shù)的比例P(m),直到P(m)不再隨m的增加而減小或小于某個(gè)值(如l%或5%)時(shí),吸引子被認(rèn)為完全展開。

在不同嵌入維下正常及故障狀態(tài)統(tǒng)計(jì)的P(m)如圖3所示。從圖中可以看出,當(dāng)m=5和m=4時(shí),正常狀態(tài)及內(nèi)圈故障狀態(tài)的P(m)不再明顯減少,因此正常狀態(tài)最佳嵌入維取5,故障狀態(tài)最佳嵌入維均取4。

(a)正常狀態(tài)

根據(jù)嵌入維和延時(shí),采用相空間重構(gòu)技術(shù)繪制正常和內(nèi)圈故障振動(dòng)信號(hào)的混沌吸引子如圖4所示。從圖中可以看出,混沌吸引子在有限的空間不斷的纏繞,呈現(xiàn)一定自相似的狀態(tài)。

圖4 滾動(dòng)軸承振動(dòng)信號(hào)混沌吸引子

2 滾動(dòng)軸承振動(dòng)的混沌特性分析

首先采用遞歸圖對(duì)實(shí)測軸承不同振動(dòng)狀態(tài)的相關(guān)性和確定性進(jìn)行檢驗(yàn),然后使用CLY方法判定軸承振動(dòng)的混沌特性,最后結(jié)合功率譜對(duì)振動(dòng)進(jìn)行頻譜分析。

2.1 遞歸圖

為分析序列周期及非周期成分的相對(duì)大小,揭示動(dòng)力系統(tǒng)內(nèi)部相似性的先驗(yàn)知識(shí),文獻(xiàn)[23]提出了遞歸圖(Recurrence Plot,RP)方法。周期系統(tǒng)信號(hào)被噪聲污染時(shí)可能表現(xiàn)出非周期性,這表明系統(tǒng)不會(huì)準(zhǔn)確地回到以前的狀態(tài)變量,但可能返回很接近以前某時(shí)刻的狀態(tài),因此可用來表現(xiàn)周期及非周期的程度。由(4)式采用嵌入技術(shù)把信號(hào)重構(gòu)為向量Xl(l=1,2,…,N-(m-1)τ)。然后求Xl的任意2個(gè)時(shí)刻i,j(i,j=1,2,…,N-(m-1)τ)的距離,即

dij=‖Xi-Xj‖ 。

(8)

如果dij≤r,以(i,j)為縱坐標(biāo)做一點(diǎn),即繪制出遞歸圖,表現(xiàn)重構(gòu)的軌線重復(fù)或遞歸其自身的信息,表示了系統(tǒng)的時(shí)間關(guān)聯(lián)情況。試驗(yàn)中,計(jì)算正常狀態(tài)及內(nèi)圈故障狀態(tài)信號(hào)時(shí)分別取r=2和r=1,得到的遞歸圖如圖5所示。

圖5 滾動(dòng)軸承振動(dòng)信號(hào)的遞歸圖

與文獻(xiàn)[24]中周期序列、混沌序列及隨機(jī)序列的遞歸圖進(jìn)行對(duì)比分析發(fā)現(xiàn):軸承正常狀態(tài)的遞歸圖存在較明顯的平行于對(duì)角線的線段,但帶狀線不是連續(xù)的,說明正常振動(dòng)為非嚴(yán)格周期狀態(tài),存在隨機(jī)性因素(部分為噪聲干擾);內(nèi)圈故障狀態(tài)遞歸圖也存在較明顯的平行于對(duì)角線的線段,與混沌系統(tǒng)類似,表明其也具有周期成分,但介于周期與隨機(jī)之間,振動(dòng)為混沌特征,且故障狀態(tài)越大,遞歸圖與白噪聲信號(hào)越相似;但軸承正常和內(nèi)圈故障狀態(tài)周期性均不明顯,由于繪制遞歸圖受噪聲和r的影響,其周期和非周期性大小只可作為定性的判斷,還需要參考其他方法進(jìn)行綜合的分析判定。

2.2 CLY方法

為判別信號(hào)的混沌特性,以虛假鄰域方法為基礎(chǔ),文獻(xiàn)[25]提出了CLY算法:由(4)式將xn(n=1,2,…,N)嵌入到m維的相空間,對(duì)選定的m,由 (6) 式計(jì)算每個(gè)Xl(m)(l=1,2,…,N-(m-1)τ)的最近鄰點(diǎn)Xη(l)(m),得到平均預(yù)測誤差

。(9)

定義函數(shù)E(m)=E*(m+1)/E*(m)。對(duì)選擇的任意m,通常隨機(jī)信號(hào)的E(m)=1或在其上下波動(dòng);而對(duì)混沌系統(tǒng),隨m的增大穩(wěn)定在1附近。故根據(jù)m-E(m)曲線可簡單的判別信號(hào)的混沌特性。正常及內(nèi)圈故障狀態(tài)振動(dòng)信號(hào)的m-E(m)曲線如圖6所示,從圖中可以看出,當(dāng)m≥5時(shí)E(m)趨于1,可簡單判別出不同振動(dòng)狀態(tài)都存在混沌特性。

2.3 功率譜

對(duì)xn(n=1,2,…,N)加上周期條件xn=xn+j,求其自相關(guān)函數(shù)(離散卷積)

(10)

然后對(duì)Cn進(jìn)行離散Fourier變換,計(jì)算Fourier變換系數(shù)[26]

(11)

圖6 滾動(dòng)軸承振動(dòng)信號(hào)的m-E(m)曲線

圖7 不同信號(hào)的功率譜

正常和內(nèi)圈故障狀態(tài)振動(dòng)信號(hào)的功率譜如圖8所示,從圖中可以看出:正常狀態(tài)振動(dòng)信號(hào)有較明顯的獨(dú)立尖峰譜線,表明出較強(qiáng)的周期性,而存在連續(xù)的譜線,又表明其也具有混沌性,部分譜線分布類似于隨機(jī)信號(hào)譜線,則對(duì)應(yīng)于其中由噪聲引起的隨機(jī)成分;內(nèi)圈故障狀態(tài)振動(dòng)信號(hào)的功率譜類似于Lorenz序列的功率譜圖,其功率譜線不同于周期和隨機(jī)信號(hào)的離散譜,表現(xiàn)為寬帶連續(xù)譜,還表現(xiàn)出在高頻段隨頻率指數(shù)衰減的特性,因此故障狀態(tài)存在很強(qiáng)的混沌性,但周期性并不明顯。綜合分析可知,軸承不同狀態(tài)振動(dòng)信號(hào)的頻率分布特性存在顯著的差異,正常狀態(tài)具有明顯的周期性,而內(nèi)圈故障狀態(tài)具有明顯的混沌特性,與遞歸圖判定結(jié)論基本一致。

圖8 滾動(dòng)軸承振動(dòng)信號(hào)的功率譜

3 滾動(dòng)軸承振動(dòng)的超混沌特性分析

Lyapunov指數(shù)譜λi(i=1,2,…,m-1)是表示相空間各方向軌道旁的膨脹及收縮的平均度量,由系統(tǒng)長時(shí)間演變決定,反映了初始時(shí)刻2個(gè)無限靠近的點(diǎn)隨時(shí)間分離的情況,每個(gè)λi都是相空間不同方向上相對(duì)運(yùn)動(dòng)局部變形的平均,不同特征的運(yùn)動(dòng)與λi的關(guān)系見表1。

混沌吸引子的局部不穩(wěn)定和某方向上的指數(shù)發(fā)散或膨脹相對(duì)應(yīng),有2個(gè)或2個(gè)以上正的λi系統(tǒng)就為超混沌,否則系統(tǒng)為非混沌。現(xiàn)采用LS-SVM判別系統(tǒng)方程方法得到軸承振動(dòng)的Lyapunov指數(shù)譜[27-28],研究滾動(dòng)軸承不同狀態(tài)振動(dòng)的超混沌特性。

表1 不同運(yùn)動(dòng)與Lyapunov指數(shù)的關(guān)系

3.1 LS-SVM算法

由(4)式重構(gòu)xn(n=1,2,…,N)在m維歐式空間的一條軌道Xl(l=1,2,…,N-(m-1)τ),則Xl能再現(xiàn)原系統(tǒng)的動(dòng)力學(xué)特性,設(shè)Xl到Xl+τ的映射關(guān)系為G[6],即

(12)

顯然Xi+τ的前m-1個(gè)數(shù)值即為Xi的后m-1個(gè)數(shù)值,只要確定Xl→Xl+mτ的關(guān)系F即可確定G,設(shè)

xl+mτ=F(Xl) 。

(13)

1)LS-SVM辨識(shí)未知系統(tǒng)

采用LS-SVM辨識(shí)F。先由原始信號(hào)xn(n=1,2,…,N)構(gòu)造訓(xùn)練數(shù)據(jù)集(Xl,xj+τ),其中l(wèi)=1,2,…,N-mτ。利用高維特征空間的線性函數(shù)擬合樣本集[29]

xl+mτ=ωTφ(Xl)+b,

(14)

則原空間中的約束條件和優(yōu)化問題可表達(dá)為

s.t.xj+mτ=ωTφ(Xl)+b+el,

(15)

式中:el為誤差,其組成的誤差向量為e,偏差為b,權(quán)重為C。用Lagrange乘子al有

b+ej-xj+mτ}。

(16)

由KKT條件并消除el和ω,得到線性方程組

(17)

式中:Y=[x1+mτ,x2+mτ,…,xN]T;α=[α1,α2,…,αN-mτ]T;I=[1,1,…,1]T;E為N-mτ階單位陣;l=(ki,j)(N-mτ)×(N-mτ),采用徑向基函數(shù)得到K的每個(gè)元素

(18)

解方程(17)得到α和b,則LS-SVM回歸函數(shù)逼近 (12) 式中的F為

(19)

2)計(jì)算Lyapunov指數(shù)譜

由 (12) 式得到

ΔXl+τ=DG(Xl)·ΔXl,

(20)

(21)

式中:DG(Xn)為映射G的Jacobian矩陣;ΔXl表示重構(gòu)的向量Xl在l時(shí)刻的微小變化。通過DG(Xl)矩陣,在l+τ時(shí)刻重構(gòu)相空間的向量Xl+τ值的微小變化中將被表現(xiàn)出來。繼續(xù)下去,該作用會(huì)累積到l+rτ(r∈N)時(shí)刻重構(gòu)相空間的向量Xl+rτ取值的變化,其關(guān)系為

ΔXl+rτ=DG(Xl+(r-1)τ)…DG(Xl)·ΔXl。(22)

(23)

按大小排列λi(i=1,2,…,m)。因?yàn)?(20) 式中矩陣會(huì)產(chǎn)生分?jǐn)?shù)冪和指數(shù),所以,ΔXl+τ往往是病態(tài)陣,不能求得其全部精確特征值。可利用長乘積矩陣的分解方法對(duì) (23) 式求解。首先定義

Ar=DG(Xl+(r-1)τ)…DG(Xl)=A(r)A(r-1)…A(1)

,(24)

對(duì)此長乘積矩陣,遞歸計(jì)算

A(i)Q(i-1)=Q(i)R(i),i=1,2…r,

(25)

式中:Q(i)為正交矩陣;R(i)為上三角矩陣;Q(0)為m階單位陣。按 (25) 式方法將QR分解r次,則Ar的QR分解為Ar=Q(r)R(r)R(r-1) …R(1),故 (22) 式的本征值可以通過下式求出

(26)

進(jìn)而求得系統(tǒng)的Lyapunov指數(shù)譜為

(27)

3.2 計(jì)算結(jié)果分析

采用前文給出的實(shí)測滾動(dòng)軸承正常和內(nèi)圈故障狀態(tài)的振動(dòng)信號(hào),分析其振動(dòng)的超混沌特性。計(jì)算同樣采用相關(guān)法和虛假鄰近點(diǎn)法分別確定延時(shí)τ和嵌入維m。LS-SVM計(jì)算參數(shù)取為:C=1 000,a=1,計(jì)算出的Lyapunov指數(shù)譜見表2。

表2 不同狀態(tài)振動(dòng)的Lyapunov指數(shù)譜

由表可知,軸承正常及內(nèi)圈故障狀態(tài)振動(dòng)均為超混沌,但正常狀態(tài)振動(dòng)的超混沌特征不明顯。當(dāng)故障增大時(shí),軸承超混沌特性程度也增強(qiáng),并且表現(xiàn)的更加明顯。

此外,對(duì)軸承外圈、滾動(dòng)體及保持架不同故障狀態(tài)進(jìn)行了相同的混沌及超混沌試驗(yàn)分析,均得到同樣的結(jié)果。因此軸承振動(dòng)的混沌程度與故障大小存在一定的對(duì)應(yīng)關(guān)系,正常狀態(tài)振動(dòng)比較平穩(wěn)有序,較大故障的振動(dòng)狀態(tài)會(huì)增加混亂程度。

4 結(jié)束語

應(yīng)用4種方法,從多個(gè)方面對(duì)實(shí)測軸承正常及故障狀態(tài)振動(dòng)的非線性混沌特性進(jìn)行了判定和分析,研究表明:軸承正常及故障狀態(tài)振動(dòng)的混沌特征具有明顯的差異性,研究軸承振動(dòng)的非線性混沌特性要結(jié)合不同的振動(dòng)狀態(tài)。上述研究結(jié)果有助于理解軸承振動(dòng)的復(fù)雜演化規(guī)律,為今后使用其振動(dòng)信號(hào)混沌特性提取故障特征及進(jìn)行故障診斷奠定一定的基礎(chǔ)。

猜你喜歡
振動(dòng)故障信號(hào)
振動(dòng)的思考
信號(hào)
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
振動(dòng)與頻率
故障一點(diǎn)通
基于FPGA的多功能信號(hào)發(fā)生器的設(shè)計(jì)
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動(dòng)性
奔馳R320車ABS、ESP故障燈異常點(diǎn)亮
基于LabVIEW的力加載信號(hào)采集與PID控制
故障一點(diǎn)通
主站蜘蛛池模板: 亚洲成年网站在线观看| 91精品专区国产盗摄| 亚洲国产综合自在线另类| 香蕉eeww99国产在线观看| 国产精彩视频在线观看| 国内嫩模私拍精品视频| 成人国产一区二区三区| 亚洲男人天堂2018| 亚洲日韩久久综合中文字幕| 国产95在线 | 亚洲二三区| 亚洲欧美一区二区三区麻豆| 成人va亚洲va欧美天堂| 婷婷午夜天| 麻豆精品久久久久久久99蜜桃| 久久女人网| 亚洲成人在线免费| 毛片免费试看| 中文字幕免费视频| 精品久久久久久久久久久| 亚洲精品午夜无码电影网| 久久99国产精品成人欧美| 欧美成人午夜影院| 成人第一页| 国产午夜福利片在线观看| 亚洲福利片无码最新在线播放| 夜色爽爽影院18禁妓女影院| 亚洲 欧美 偷自乱 图片| 亚洲伊人久久精品影院| 国产在线精彩视频论坛| 亚洲精品爱草草视频在线| 精品国产www| 国产一级无码不卡视频| 日韩成人午夜| 国产色网站| 思思热在线视频精品| 性色在线视频精品| 国产美女精品一区二区| 97精品久久久大香线焦| 国产精品林美惠子在线观看| 亚亚洲乱码一二三四区| 99久久精品无码专区免费| A级全黄试看30分钟小视频| 无码又爽又刺激的高潮视频| 欧美一区二区三区不卡免费| 五月婷婷亚洲综合| 国产综合精品日本亚洲777| 久久精品66| 国产午夜福利亚洲第一| 国产精品亚洲一区二区在线观看| 99久久性生片| 无码精品一区二区久久久| 成人福利在线观看| 欧美另类图片视频无弹跳第一页| 91精品国产91久久久久久三级| 在线中文字幕网| 中国一级特黄大片在线观看| 农村乱人伦一区二区| 午夜丁香婷婷| 色综合中文综合网| 亚洲欧美另类色图| 国产在线91在线电影| 欧美精品三级在线| 97人人做人人爽香蕉精品| 亚洲三级a| 久青草网站| 国产亚洲欧美在线中文bt天堂| 69视频国产| 伊人国产无码高清视频| 亚洲成人精品久久| 国产欧美日韩综合一区在线播放| 国产另类视频| 午夜电影在线观看国产1区| 国内精品久久人妻无码大片高| 97se亚洲综合不卡| 激情综合激情| 中文字幕亚洲精品2页| 四虎国产成人免费观看| 91精品最新国内在线播放| 日韩国产一区二区三区无码| 国产精品亚洲欧美日韩久久| 天天摸夜夜操|