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

第三屆全國隨機動力學學術會議優(yōu)秀論文基于EEMD與FastICA的損傷異常識別與定位

2016-07-26 02:21:48姜紹飛陳志剛沈清華吳銘昊麻勝蘭
振動與沖擊 2016年1期

姜紹飛, 陳志剛, 沈清華, 吳銘昊, 麻勝蘭

(福州大學 土木工程學院,福州 350108)

?

第三屆全國隨機動力學學術會議優(yōu)秀論文基于EEMD與FastICA的損傷異常識別與定位

姜紹飛, 陳志剛, 沈清華, 吳銘昊, 麻勝蘭

(福州大學 土木工程學院,福州350108)

摘要:為了準確地提取結構損傷異常信息,消除小波奇異值分解時存在需要特定的小波基和分解層數(shù)以及經(jīng)驗模態(tài)分解(EMD)方法存在諸如虛假模態(tài)混疊等問題,提出一種基于改進的總體平均經(jīng)驗模態(tài)分解(EEMD)與快速獨立分量分析(FastICA)相結合的提取結構損傷特征并進行識別與定位的新方法。首先,通過EEMD對結構動力響應信號進行預處理并用 FastICA提取出包含損傷信息的特征分量對結構響應異常進行識別和初步定位;然后,計算歸一化的源分布向量(NSDV)的最大值,并根據(jù)該最大值精確定位結構損傷。最后,通過框架數(shù)值算例和試驗進行了所提方法的驗證,結果表明該算法能夠較好地進行結構損傷異常的識別與定位。

關鍵詞:總體平均經(jīng)驗模態(tài)分解;快速獨立分量分析;損傷定位;源分布向量

當結構發(fā)生損傷后,造成其固有頻率和剛度的改變,進而使得結構的動力響應發(fā)生變化,通過現(xiàn)代的信號處理方法可有效地從動力響應中提取結構早期損傷的特征,從而提高損傷預測的準確性和可靠性[1]。因此,選擇恰當?shù)男盘柼幚矸椒ǎ崛〕鼋Y構早期的損傷異常信息,是進行結構健康診斷與預后分析的前提[2]。

結構損傷往往對應著信號的奇異點,小波分析[3]在時域和頻域分析中都具有表征信號局部特征的能力。但是,小波分析在工程應用中存在需要人為選擇小波基及其分析結果易受假定閾值的影響而丟失有用信息等缺陷[4]。1998年出現(xiàn)的EMD方法在處理非平穩(wěn)信號方面較小波變換等其他傳統(tǒng)方法更具靈活性[5],可根據(jù)信號自身的內在特性進行自適應分解,但也存在信號模態(tài)混疊和端部效應等問題。為了解決EMD存在的問題,Huang等[6]提出了改進EMD的方法,即EEMD,通過加噪聲輔助分析的方法,使信號在表征不同尺度上具有連續(xù)性,克服了存在的模態(tài)混疊問題,同時EEMD分解是隨信號本身的變化而變化的,能更有效地提取出信號本身的特征。

近些年,盲源分離在各個領域得到了廣泛的應用,其基本思想是將多個觀測信號按照統(tǒng)計獨立的原則通過優(yōu)化算法分解為若干個獨立成分[7],從而實現(xiàn)信號的增強和分析。FastICA 算法[8]是盲信號處理中一種較成熟的算法,它能消除各個輸入量之間的互信息和信息冗余,分離出信息之間隱藏的內部相互獨立的成分。鑒于此,本文提出了EEMD-FastICA相結合的結構損傷特征提取的新方法,利用前者能夠自適應地分解測量響應,提取含有結構損傷信息的高階IMF成分,后者能夠分離出IMF隱藏的損傷異常信息并以此初步定位損傷,之后計算歸一化的源分布向量(NSDV)的最大值精確定位損傷異常,來解決損傷檢測中測量響應損傷異常奇異點不敏感的問題[9]。并通過框架結構在地震波激勵下的結構損傷驗證了所提方法有效性。

1相關原理

1.1總體平均經(jīng)驗模態(tài)分解(EEMD)

EEMD方法本質是一種改進的EMD算法,通過加入正態(tài)的白噪聲,使信號在不同尺度上具有連續(xù)性[10],從而達到在EMD分解中避免模式混疊的目的。

EEMD方法步驟[11]如下:

步驟1初始信號x(t)中多次加入高斯白噪聲ni(t)序列,即:

si(t)=s(t)+ni(t)

(1)

式中:si(t)為第i次加入高斯白噪聲后的信號,ni(t)為第i次加入的高斯白噪聲。

步驟2對所得的含高斯白噪聲的信號si(t)按照EMD算法分解,得到各自的IMF記為cij(t),與一個余項記為ri(t)。

步驟3重復步驟(1)和(2),每次加入的應是不同幅值的白噪聲序列。

步驟4將上述得到的各個IMF的均值進行總體平均運算,得到EEMD分解后最終的IMF,即:

(2)

式中:N為經(jīng)驗模態(tài)分解的聚合次數(shù)。

1.2快速獨立分量分析(FastICA)

獨立分量分析[12]是把觀測信號線性分解為幾個相互獨立成分的盲信號處理方法,具體描述為:

X=AS

(3)

其中:X=[x1,x2, …xm]T為m維觀測信號矩陣;S=[s1,s2,…sn]T為n維相互獨立的源信號分量矩陣;A為m×n維混合矩陣。

ICA所要進行的工作就是找到混合矩陣A的逆矩陣,然后通過式(4)就可以估計出S:

S≈Y=WX

(4)

快速獨立分量分析[13](FastICA)是基于極大化非高斯性的目標函數(shù)以負熵作為衡量各個分量之間獨立性的目標函數(shù),采用牛頓迭代算法對觀測信號的大量采樣點進行批處理,每次從觀測信號中分離出一個獨立分量,直到分離出所有的獨立分量。負熵定義式:

J(yi)={E[G(yi)-E{G(v)}]}2

(5)

式(5)中yi是具有零均值和單位方差的輸出變量,yi=WTZ, Z為對X進行中心化和白化處理后的數(shù)據(jù);v是具有零均值和單位方差的高斯隨機變量,G是一個非二次非線性的函數(shù)。在W正交的約束條件下,求取式(6)的極大值可得到W的迭代式:

(6)

k為迭代次數(shù),g是G函數(shù)的導數(shù)。重復此過程,便可逐次提取多個獨立源信號。

通過一個數(shù)值算例來說明FastICA算法,假設有兩個信號源,如圖1(a)所示,一個是含有突變特征的信號(S1(t)=3,t=5 s;S1(t)=0,其他時刻);另一個是高斯白噪聲信號(S2(t))。二者的采樣頻率均取100 Hz,通過混合矩陣A,可以得到兩個觀測信號X1(t)和X2(t),如圖1(b)所示。從圖1(b)中可以看出,源信號中的突變特征完全被噪聲所淹沒。此時,若對觀測信號進行快速獨立分量分析,則可以重新提取出信號中的突變特征,如圖1(c)中的IC1所示。

圖1 快速獨立分量分析Fig.1 Fast independent component analysis

1.3歸一化源分布向量(NSDV)

假設X=[x1,x2, …xm]T為包含突變信息的m維觀測信號矩IC=[IC1,IC2, …,ICn]T為利用FastICA得到的n維獨立分量矩陣,則

X=A×IC

(7)

式中:A為m×n維混合矩陣,式(7)展開為:

Xi(t)=ai1·IC1(t)+…+aij·ICj(t)+…+

ain·ICn(t),i=(1,2…m)

(8)

其中:aij表示獨立分量ICj(t)在觀察信號X(t)中的比重大小,稱為源分布因子。

假設第j個獨立分量ICj(t)中包含突變信息,則稱之為特征獨立分量(Feature Independent Component,F(xiàn)IC),令SDV=[a1j,a2j, …amj]T代表FIC在觀測信號X(t)中的源分布向量(Source Distribution Vector,SDV),進一步歸一化得到:

NSDV=[NSDF1,NSDF2,…,NSDFi,NSDFm]T

其中:

(9)

這里,稱NSDV為歸一化的源分布向量。

2基于EEMD-FastICA損傷定位方法

(1) EEMD預處理:采用EEMD對結構動力響應進行處理:首先,將m個傳感器采集到的結構動力響應數(shù)據(jù)X(i)(i=1,2,…,m)集合成信號矩陣X(m×l),l為數(shù)據(jù)長度;其次,對信號矩陣X的每一行按照EEMD處理;第三,高階IMF成分能夠反映信號的真實信息[14],分解后按照式(2)提取IMF中的高階成分EIMF。

(2) FastICA特征分量提取:運用FastICA對信號提取的EIMF分量進行特征分量提取:首先,對EIMF分量信號矩陣進行FastICA得到m個獨立分量IC(j) (j=1,2,…,m);其次,從m個獨立分IC(j)(j=1,2,…,m)中提取出具有信號突變點的特征分FIC(k)(k=1,2, …,m),m為特征分量的個數(shù)。

(3) 初步定位損傷:根據(jù)每個特征分量FIC(k)中信號突變點在時間軸上的位置確定相應的損傷發(fā)生時刻Tdk(k=1,2,…,m)以及FIC(k) 信號突變所在的層數(shù)初步定位損傷。

(4) 精確定位損傷:提取出每個特征分量FIC(k)所對應的SDV,按照式(9)進行歸一化,從歸一化的NSDV圖中找出最大值點,根據(jù)最大值點對應的位置精確定位損傷發(fā)生位置。

3數(shù)值算例

3.1結構模型

圖2 框架結構簡化模型Fig.2 Simplified model of frame structure

考慮一個結構底層受地震激勵作用的三層剪切結構,質量Mi= 125.53 kg,采用Rayleigh阻尼,阻尼比ζi=0.04,剛度Ki=24.2 kN/m,i=1,2,3,激勵采用EL-Centro地震波(PGA=50 gal),采樣時間為30 s,采樣頻率為50 Hz,假設有兩種工況,見表1。將該結構模型簡化為一個三自由度集中質量模型進行分析,如圖2所示。結構動力加速度響應計算采用采New mark逐步積分法,其中參數(shù)γ=0.5,β=0.25。

表1 結構損傷工況

3.2損傷模型建立

(1) EEMD預處理:利用EEMD分別對各層加速度數(shù)據(jù)Acc(i)(i=1,2,3)進行預處理,然后提取分解得到的IMF的高階成分組成3維信號矩陣EAcc;

(2) FastICA特征分量提取:根據(jù)FastICA對高階信號矩陣EAcc進行獨立分量分析得到3個獨立分量IC(j)(j=1,2,3),并從中提取出特征分量FIC(k);

(3) 損傷時刻確定和定位:特征分量FIC(k)中信號突變點在時間軸上的位置確定損傷發(fā)生時刻及根據(jù)FIC(k)所在的層數(shù)初步確定損傷位置,同時得到與之相對應的歸一化的NSDV,精確定位損傷。

3.3損傷異常判別

3.3.1工況一

無背景噪聲環(huán)境,在地震激勵作用下的工況一仿真得到加速度響應輸出如圖3(a)所示。之后,對得到的加速度響應進行EEMD分解,圖3(b)是EEMD分解之后的高階IMF成分,其能夠較好的反應原始響應的特性,但無法直接從結構加速度對應的IMF成分中識別出明顯的損傷突變特征;按3.2節(jié)步驟進行FastICA特征獨立分量提取,得到如圖3(c)所示的結果,對于損傷的框架結構,在第一層得到一個FIC,即IC1,圖中顯而易見,信號在t=10 s時存在一個明顯的突變值,因此可以認為結構可能于t=10 s時在第一層發(fā)生損傷異常;為了進一步確定損傷發(fā)生位置,按照式(9)計算圖中FIC對應的歸一化源向量(NSDV)如圖3(d)所示,可以看出NSDF的最大值出現(xiàn)在對應的第一層,由此可以確定,在t=10 s時,結構第一層發(fā)生了損傷異常,其識別結果與假設情況相符。

圖3 工況一損傷異常判別Fig.3 Abnormal damage identification in Case 1

圖4 工況二損傷異常判別Fig.4 Abnormal damage identification in Case 1

3.3.2工況二

工況二相對來說是多損傷,地震作用下的加速度響應如圖4(a)所示;采用EEMD對各層加速度響應進行分解,取其高階信號IMF成分如圖4(b)所示,圖中并沒有明顯的現(xiàn)象;結合FastICA對IMF進行特征獨立分量的提取,圖4(c)可以明顯看出在第一層IC1,t=10 s處和第三層IC3,t=20 s處,信號都存在明顯的突變點,于是可知結構的損傷可能發(fā)生在t=10 s時結構第一層,和t=20 s時結構第三層處;在此基礎上分析工況二的NSDF的值域分布如圖4(d),IC1的NSDF最大值位于框架結構第一層,IC3的NSDF在框架的第三處值最大。綜合以上可以判定,結構第一層于t=10s的時候發(fā)生損傷,t=20 s時,在第三層發(fā)生損傷異常,這與實際假定的情況相符合。

數(shù)值分析結果表明,本文所提的方法能夠直接根據(jù)獨立分量所在的位置快速地對結構損傷發(fā)生的位置進行初步確定,其效果優(yōu)于文獻[4]中所提及的小波分析法。

3.3.3不同噪聲水平下的討論比較

為了研究所提方法在不同噪聲水平的實用性,按照式(10)向各層動加速度數(shù)據(jù)和地震波數(shù)據(jù)中加入噪聲信噪比為40 dB,35 dB,30 dB,25 dB的高斯白噪聲來探討所提方法的容噪性能力。

SNR=20lg(Asignal/Anoise)

(10)

式中:A代表信號幅值。

為了方便比較,將工況一在四種信噪比下提取的IC1集成,如圖5(a);工況二在不同信噪比下提取的IC1,IC3集成,如圖5(b),黑色虛線代表IC1,黑色實線代表IC3,兩種工況不同噪聲水平下提取得到的NSDV見圖6。SNR=25 dB時,在兩種工況的特征獨立分量中仍能夠明顯識別出損傷發(fā)生時刻t=10 s,t=20 s時的突變點,但隨著噪聲水平的提高,突變點越來越不明顯,并且隨著信噪比的降低,噪聲干擾作用越來越明顯;此外,對于兩種工況在t=10 s時發(fā)生的損傷,根據(jù)對應的NSDV圖,NSDF最大值均位于結構第一層。對于工況二在t=20 s時發(fā)生的損傷對應的NSDV圖的最大值分布于第三層,不同噪聲水平下的兩種工況均符合實際情況。由此可見,所提方法在不同噪聲情況下也能夠較好的提取損傷特征量與定位損傷,具有較好的容噪性。

圖5 不同信噪比下的特征分量Fig.5 The characteristic independent components with different SNR

圖6 不同信噪比下的NSDV圖Fig.6 The NSDV with different SNR

4試驗驗證

4.1結構模型

試驗框架模型是一個7層、2跨×1跨的鋼框架縮尺模型[15],該模型平面尺寸為0.4 m0.2 m,高1.412 5 m,梁間距和柱間距均為200 mm,柱子是薄鋼板,梁是空鋼管,梁柱截面特性如表2所示,模型如圖7所示。在模型頂層跨中處進行激振,每層記錄加速度,結構損傷是通過切割柱子模擬的,通過在柱子上下端切割出4個7.5×2.5 mm2大小的缺口來模擬小損傷;通過將柱子從中間完全切斷來模擬大損傷。根據(jù)不同損傷程度和損傷位置及損傷時刻的組合,考慮兩種損傷工況,見表3。

表2 結構構件特性

損傷工況損傷情況工況一T=0.4s,第4層大損傷、第6層小損傷工況二T=0.4s,第4層大損傷、第6層小損傷,T=0.6s,第3層小損傷

4.2損傷異常判別

用加速度傳感器采集兩種工況下七層框架的各層加速度響應,利用第3.2節(jié)中所提出的EEMD-FastICA的方法,分別對兩種工況所輸出的加速度響應進行特征量提取。對于工況一,如圖8所示,可以看出,損傷突變點在T=0.4 s只出現(xiàn)在獨立分量IC1中。而工況二則分別T=0.4 s在IC1和T=0.6 s在IC2中出現(xiàn)了信號突變點,突變的時間均與試驗工況損傷發(fā)生時刻相對應,結果如圖9所示。試驗的結果初步表明,所提及的方法可以對損傷的異常點進行判斷。

圖8 工況一對應的獨立分量 Fig.8 Feature component of case l

圖9 工況二對應的獨立分量Fig.9 Feature component of case 2

圖10 工況一對應的NSDV圖Fig.10 NSDV of case 1

在此基礎上利用式(9)分析兩種工況的NSDF的值域分布。如圖10所示,對于工況一,IC1的NSDF最大值位于框架結構第四層和第六層;對于工況二,IC1的NSDF最大值在框架的第四層和第六層處,IC2最大值位于框架的第三層處,如圖11所示,這與實際損傷發(fā)生的位置相同。試驗的結果表明,利用本文所提的方法能夠較好地進行結構損傷異常的判定。

圖11 工況二對應的NSDV圖Fig.11 NSDV of case 2

5結論

(1) 本文提出的基于EEMD-FastICA相結合的損傷異常檢測與定位方法,通過EEMD可以自適應的對結構響應進行分解而不丟失有用信息,結合FastICA對高階IMF分量進行損傷敏感特征量的提取,可以有效地檢測出測量響應中存在的損傷異常值,從而確定損傷發(fā)生的時刻及初步定位損傷。

(2) 通過計算特征分量的NSDF值,并進行歸一化, NSDF值的最大值分布情況可以更加精確地定位損傷發(fā)生的位置。

(3) 計算了不同信噪比下提取得到的特征獨立分量及相應的NSDV值,當信噪比SNR大于25 dB本文所提的方法能夠較好地檢測出損傷異常信息并精確定位損傷。

可見,基于外部激勵和結構動力響應信號,運用本文提出的方法可以提取結構的損傷異常并進行有效地識別和定位。同時,對于外部激勵未知、測量數(shù)據(jù)不完備,更強背景噪聲環(huán)境下的損傷策略還需要進一步研究和探討。

參 考 文 獻

[1] 姜紹飛.結構健康監(jiān)測-智能信息處理及應用[J].工程力學,2009,26(增刊2):184-212.

JIANG Shao-fei.Structural health monitoring intelligent information processing and application[J]. Journal of Engineering Mechanics,2009, 26(Sup2):1-7.

[2] 高維成,劉偉,鄒經(jīng)湘.基于結構振動參數(shù)變化的損傷探測方法綜述[J].振動與沖擊, 2004, 23(4):1-7.

GAO Wei-cheng, LIU Wei, ZOU Jing-xiang. Damage detection methods through changes of vibration parameters: a summary review [J]. Journal of Vibration and Shock, 2004, 23(4): 1-7.

[3] 李宏男,孫鴻敏.小波分析在土木工程領域中的應用[J].世界地震工程,2003,19(2): 16-22.

LIHong-nan, SUN Hong-min. Application of wavelet analytical method to civil engineering[J]. World Earthquake Engineering, 2003,19(2): 16-22.

[4] 劉濤,李愛群,丁幼亮.小波分析在結構損傷識別中的應用[J].地震工程與工程振動,2008, 28(2):29-35.

LIU Tao, LI Ai-qun, DING You-liang. Application of wavelet analysis to structural damage identification[J]. Journal of Earthquake Engineering and Engineeing Vibration, 2008,24(2):29-35.

[5] Huang N E, Shen Z, Long S R. A new view of nonlinear water waves: The Hilbert Spectrum1[J]. Annual Review of Fluid Mechanics, 1999, 31(1):417-457.

[6] Wu Z H, Huang N E. Ensemble empirical mode decomposition: a noise-assisted data analysis method [J].Advances in Adaptive Data Analysis, 2009, 1: 1-41.

[7] Hyv?rinen A, Oja E. Independent component analysis: a tutorial[J].Neural Networks, 2000, 13(45): 411-430.

[8] Hyvarinen A. Fast and robust fixed-point for independent component Analysis [J]. IEEE Transactions on Neural Networks, 1999,10(3):626-634.

[9] 李富強,劉國華,吳志根.基于雙譜和奇異值分解的結構損傷試驗[J].浙江大學學報(工版),2012,46(10):1872-1879.

LI Fu-qiang,LIU Guo-hua,WU Zhi-gen.Experimental study of structural damage based on bispectral analysis and singular value decomposition[J]. Journal of Zhejiang University, 2012, 46 (10): 1872-1879.

[10] 鄭近德,程軍圣,楊宇.改進的 EEMD 算法及其應用研究[J].振動與沖擊,2013,32(21): 21-26.

ZHENG Jin-de,CHENG Jun-sheng,YANG Yu. Modified EEMD algorithm and its applications [J].Journal of Vibration and Shock,2013,32(21): 21-26.

[11] Wu Z. Huang N E.Ensemble empirical mode decomposition: a noise assiste data analysis method[J]. Advances in Adaptive Data Analysis,2009,1:1-41.

[12] Shimizu S, Hyvarinen A, Kano Y, et al. Discovery of non-gaussian linear causal models using ICA[C]//21st Conference on Uncertainty in Artificial Intelligence. Edinburgh,Scotland:UAI,2005.

[13] 路亮,龍源,鐘明壽,等. FastICA算法在低信噪比爆破振動信號信噪分離中的應用研究[J].振動與沖擊,2012,31(17):32-37.

LU Liang,LONG Yuan,ZHONG Ming-shou, et al.Separating noise from ablasting vibration signal based on fastICA[J]. Journal of Vibration and Shock, 2012,31(17):32-37.

[14] Lei Y, He Z, Zi Y. EEMD method and WNN for fault diagnosis of locomotive roller bearings[J].Expert Systems with Applications, 2011, 38(6): 7334-7341.

[15] Jiang S F,Wu S Y,Dong L Q.A time-domain structural damage detection based on improved multi-particle swarm coevolution optimization algorithm[J]. Mathematical Problems in Engineering,2014:232763, 11papes.

基金項目:國家自然科學基金(51278127,50878057);國家十二五科技支撐計劃(2012BAJ14B05)

收稿日期:2014-12-01修改稿收到日期:2015-03-06

中圖分類號:TU317

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.01.032

Damage detection and locating based on EEMD-FastICA

JIANG Shao-fei, CHEN Zhi-gang, SHEN Qing-hua, WU Ming-hao, MA Sheng-lan

(College of Civil Engineering, Fuzhou University, Fuzhou 350108, China)

Abstract:It is known that the wavelet decomposition requires specific wavelet basis functions and decomposition layers. Meanwhile, there exist some problems in the empirical mode decomposition (EMD), such as, false modes. To avoid the disadvantages above, a method of structural damage detection and locating based on the ensemble empirical mode decomposition (EEMD) and the fast independent componentan analysis (FastICA) was presented to accurately extract the structural damage novelty information. At first, the measured structural dynamic responses were preprocessed with EEMD,and then FastICA algorithm was used to extract the feature components involving the damage information so as to detect the structural response anomalies and preliminarily locate damage. After ward, the maximum value of the normalized source distribution vector (NSDV) was computed to accurately locate the structural damage. Finally, a frame numerical example and test were conducted, the results showed that the proposed method can successfully detect damages and locate them.

Key words:EEMD; FastICA; damage locating; NSDV

第一作者 姜紹飛 男,博士,教授,博士生導師,1969年生

主站蜘蛛池模板: 午夜无码一区二区三区| 亚洲免费播放| 99伊人精品| 欧美日韩一区二区在线免费观看| 欧美一级在线播放| 亚洲人成在线精品| 青青草久久伊人| 亚洲 成人国产| 国产成人a毛片在线| 国产91高清视频| 欧美亚洲国产日韩电影在线| 成人在线亚洲| 毛片手机在线看| 成人免费一区二区三区| 伊人五月丁香综合AⅤ| 亚洲综合天堂网| 人人爽人人爽人人片| 91人妻日韩人妻无码专区精品| 亚洲日韩精品欧美中文字幕| 亚洲精品老司机| 日韩美毛片| 国产又粗又猛又爽视频| 最新国产高清在线| 免费A∨中文乱码专区| 亚洲成人一区在线| 久久一本精品久久久ー99| 97青青青国产在线播放| 日韩一区二区三免费高清| 1级黄色毛片| 欧美a√在线| 国产又爽又黄无遮挡免费观看| a级毛片一区二区免费视频| 亚洲乱伦视频| 日韩欧美网址| 国产精品一区二区不卡的视频| 亚洲精品成人7777在线观看| 国产成人无码久久久久毛片| 国产男女免费视频| 欧美成人午夜影院| 亚洲成aⅴ人在线观看| 免费女人18毛片a级毛片视频| 激情综合网激情综合| 欧美成人怡春院在线激情| 最新精品久久精品| 国产青青草视频| 欧美.成人.综合在线| 99手机在线视频| 亚洲AV一二三区无码AV蜜桃| 99久久精品无码专区免费| 五月婷婷亚洲综合| jijzzizz老师出水喷水喷出| a级毛片一区二区免费视频| 国产成人在线无码免费视频| 欧美精品黑人粗大| 91久久国产成人免费观看| 午夜日本永久乱码免费播放片| 亚洲精品成人片在线观看| 在线精品亚洲国产| 欧美黄色网站在线看| 欧美午夜精品| 欧美精品1区| 亚洲中文无码av永久伊人| 亚洲国产日韩视频观看| 免费看美女毛片| 成人免费午间影院在线观看| 久久香蕉国产线| 视频一本大道香蕉久在线播放| 亚洲av无码久久无遮挡| 中文字幕亚洲无线码一区女同| 国产日本视频91| 日本一区二区不卡视频| 欧美成人精品在线| 色丁丁毛片在线观看| 日本福利视频网站| 成人国产免费| 亚洲日韩精品综合在线一区二区 | 国产在线观看一区二区三区| 国产日韩欧美一区二区三区在线| 少妇精品在线| 欧美自慰一级看片免费| 91精品国产91久久久久久三级| 亚洲成人在线免费|