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

橋梁動力測試信號的自適應分解與重構

2015-12-30 03:42:56單德山,李喬,黃珍
振動與沖擊 2015年3期
關鍵詞:橋梁

第一作者單德山男,博士,教授,博士生導師,1969年生

橋梁動力測試信號的自適應分解與重構

單德山,李喬,黃珍

(西南交通大學土木工程學院橋梁工程系,成都610031)

摘要:針對橋梁結構動力測試信號噪聲水平高、難以分離結構有效信號的特點,在總體平均經驗模態分解方法和主成分分析的基礎上,建立了自適應分解與重構方法。對經驗模態分解結果的模態混疊現象進行深入分析,利用白噪聲概率密度函數的均勻性對模態混疊模式一進行了改進,基于相關性分析改進了模態混疊模式二,改進后的分解方法在計算效率和分解精度上均有較大提升;隨后對所有分解獲得的固有模態函數進行多尺度主成分分析,實現降噪和選擇并重構測試信號。分別用模擬信號和實際橋梁測試信號對所提方法的有效性進行了驗證。結果表明:改進后的信號自適應分解和重構方法能在降噪的同時,有效地提取橋梁結構信息,可用于實際橋梁結構的動力測試分析中。

關鍵詞:橋梁;動力測試;EEMD;信號分解;信號重構

基金項目:國家自然科學基金(51078316);國家重點基礎研究發展計劃(2013CB036300-2);四川省科技計劃項目(2011JY0032);鐵路科技研究開發計劃項目(2011G026-E、2012G013-C)

收稿日期:2013-12-05

中圖分類號:U448.25;U448.27文獻標志碼:A

基金項目:國家自然科學

Adaptive decomposition and reconstruction for bridge structural dynamic testing signals

SHANDe-shan,LIQiao,HUANGZhen(Bridge Engineering Department, Southwest Jiaotong University, Chengdu 610031, China)

Abstract:In order to extract structural information from bridge structural dynamic signals with high noise level, a novel adaptive decomposition and reconstruction method was proposed by combining the ensemble empirical mode decomposition (EEMD) method and the principal component analysis (PCA) method for the specific characteristics of bridge structural dynamic signals. Based on the in-depth analysis of mode mixing in results of empirical mode decomposition, the uniformity of probability density function of white noise was adopted to improve the pattern I of mode mixing, and the correlation analysis was used to ameliorate the pattern II of mode mixing, then the calculation efficiency and decomposition accuracy were upgraded greatly for the improved EEMD. The multi-scale principal components analysis was implemented for all of the intrinsic mode functions (IMFs) obtained with the improved EEMD to reduce noise and select IMFs. Moreover, the dynamic signals were reconstructed. The effectiveness of the proposed method was verified with both the simulated signals and testing signals from real bridge structures. The results showed that the proposed method can be used to decompose adaptively and denoise effectively the bridge dynamic signals with high noise, and extract accurately the structural information from the testing signals, furthermore, it is applicable for the dynamic testing analysis of real bridge structures.

Key words:bridge; dynamic test; ensemble empirical mode decomposition (EEMD); signal decomposition; signal reconstruction

信號分解技術廣泛應用于橋梁結構的模態參數識別中,如Fourier變換、小波變換以及HHT變換[1]。因HHT變換中的經驗模式分解(Empirical Mode Decomposition, EMD)能對非線性、非平穩信號進行自適應分解而受到廣泛關注。與Fourier變換和小波變換不同的是,EMD無需基函數,屬于數據驅動的自適應信號分解[1]。然而,EMD也存在一些問題,如分解獲得的同一固有模態(Intrinsic Mode Function, IMF)中存在振幅完全不同的振動模式(本文定義模態混疊模式一),或不同的IMF中存在相似的振動模式(本文定義模態混疊模式二)[2]。

針對模態混疊問題,國內外學者提出了很多改進方法:Huang等[3]以降低EMD自適應性為代價,提出間歇檢測方法;譚善文等[4]通過約束特征尺度大小避免EMD的模態混疊,其方法也犧牲了EMD自適應性。近年來,隨著白噪聲統計特性研究的進展,基于EMD二進自適應濾波的特點[5],建立了噪聲輔助數據分析(Noise-Assisted Data Analysis, NADA)[6],在此基礎上Wu等[7]提出了總體平均經驗模態分解方法(Ensemble Empirical Mode Decomposition, EEMD),該方法在EMD中加入Gaussian白噪聲對模態混疊進行改善;研究表明白噪聲添加次數越多,效果越顯著[7]。

橋梁結構測試信號分解和重構的目的是從分解后的成分中選擇有效成分以表征結構的真實響應。對于EMD/EEMD方法的分解結果來說,有效成分的選擇實際上就是IMF的選擇。陳仁祥等[8]提出了自動選擇IMF分量的算法,并應用于滾動軸承振動信號的降噪;陳雋等[9]通過去除前若干階IMF分量和余量,實現IMF的選擇,進而對疲勞應變信號的降噪;基于IMF中噪聲信號的統計分析,Terrien等[10]對IMF1進行平穩性檢驗,確定其中的有效成分,進而實現IMF的選擇,并應用于機械和生物信號的處理中。

然而,EEMD必會增加新的模態混疊,即重構信號中必然包含殘余的噪聲,且由于白噪聲的隨機性使得不同噪聲的疊加將獲得不同的模態數目,進而影響總體平均的結果[7];另外EEMD方法是通過犧牲計算效率來解決模態混疊的。既有IMF的選擇方法基本上均依據白噪聲EMD得到的IMF成分的能量密度與其平均周期的乘積為常量的特點[7]進行的。實際橋梁的動力測試多采用環境激勵的方法進行,其動力測試信號具有噪聲水平高、難以分離結構有效信號的特點[11]。針對橋梁結構動力測試信號特點,本文對EEMD算法進行了改進,并用多尺度主成分分析對IMF進行降噪和選擇,達到精確提取結構有效信號的目的。

1EEMD的改進

1.1EEMD

(1)

EEMD整個過程如圖 1所示,由該圖可知,確保不同xi(t)EMD得到的IMF數量相同是關鍵所在,在不同噪聲情況,要使得每一xi(t)分解得到相同數量IMF難以保證[7]。

圖1 EEMD基本流程 Fig.1 The Flowchart of EEMD

1.2改進EEMD

EMD中存在兩種不同性質的模態混疊情況[2],本文從信號處理的角度出發,分別進行改進。

1.2.1模態混疊模式一的改進

(1)在信號x(t)中分別添加N個滿足N(0,1)分布的白噪聲,即:

xi(t)=x(t)+wi(t)(i=1,2,…,N)

(2)

(3)

殘量為:

(4)

(3)用EMD分解噪聲wi(t),獲得各噪聲的IMF分量以及殘量

(5)

(6)

(7)

(8)

(9)

(7)重復第(6)步,直到所得的殘量不再可分。測試信號x(t)的分解殘量R(t)為

(10)

測試信號x(t)可表示為

(11)

1.2.2模態混疊模式二的改進

若不同IMF中存在相同的振動模式,即不同IMF的形狀相似。為此,本文從信號相關性[12]的角度對分解獲得IMF進行相關性檢驗,進而將相同振動模式的不同IMF合并,進而改進該類模態混疊模式。

設分解后的IMFi=si(t)+ni(t)(i=1,2,…,K),其中si(t)為IMF中的特征信號,ni(t)為噪聲。不同IMF間的互相關函數為

Rij(τ)=Rsisj(τ)+Rsinj(τ)+

Rsjni(τ)+Rninj(τ)(i≠j)

(12)

因噪聲互不相關,即Rsinj(τ)=Rsjni(τ)=Rninj(τ)=0。式(12)可寫為

Rij(τ)=Rsisj(τ)(i≠j)

(13)

(14)

2IMF的降噪與選擇

在IMF選擇中,Flandrin等[13]認為IMF1均為噪聲,獲得IMF1的噪聲能量后,估計其他IMF的噪聲水平和置信區間,并據此對IMF進行降噪和選擇。近年的研究表明IMF1中仍可能含有結構信息[10];且各IMF噪聲能量的關系是在某些假設條件下建立的,對橋梁結構動力測試信號的適應性不好。為克服該問題,采用多尺度主成分分析(Multi-scale Principal Component Analysis, MPCA)[14]處理每一IMF,實現IMF的降噪和選擇。

2.1主成分分析

主成分分析[12](Principal Component Analysis, PCA)在方差分析的基礎上,將數據投影到方差最大的正交主成分上,使得多維數據的互相關最小,進而實現多維數據的降維。

C=E(XXT)=Um×mΛm×mUTm×m

(15)

式中:Λ為協方差矩陣特征值的對角矩陣,Λ=diag(λ1,λ2,…,λm),且λ1≥λ2≥…≥λm;U為特征向量組成的正交矩陣。

(16)

pi對應的特征值λi定義為該主成分的貢獻率φi

(17)

前l個主成分得累計貢獻率ψi定義為:

(18)

實際應用時,確定累計貢獻率的閾值后,即可確定主成分的選擇數量。

2.2多尺度主成分分析

多尺度主成分分析將主成分分析的正交相關性能力與小波多尺度分解[15]能力相結合[14],利用噪聲的互不相關特點對測試信號進行降噪,從橋梁動力信號中提取結構特征。實現過程如下:

(1)對IMFi(i=1,2,…,K)進行多尺度小波分解;

(2)對每一尺度進行主成分分析,計算各自小波系數的協方差矩陣、主成分分量;確定主成分數量并計算閾值,獲得大于或等于閾值的小波系數;

(3)對檢測到的顯著事件尺度進行組合,依據所選分值和閾值重構IMF,獲得降噪后的IMFi,并計算其PCA;

(4)對降噪后IMFi進行相關性檢驗,合并具有相同振動模式的IMF,用PCA方法選擇并獲得最終的分解結果。

3測試與驗證

利用改進EEMD分解橋梁結構測試信號后,應用MPCA對每一分解成分進行降噪,實現結構動力測試信號的處理,分別用模擬信號和實際橋梁健康監測系統的動力測試信號對本文所提方法進行驗證,以說明本文改進方法的效果。

3.1模擬信號

模擬信號由兩個1 Hz和5 Hz的余弦信號疊加噪聲水平約15%的隨機噪聲組成:

分解結果如圖2所示,圖中(a)為EEMD分解結果,(b)為改進EEMD的分解結果,(c)為MPCA降噪后的結構信號。比較圖2(a)和(b)可知:EEMD獲得10個IMF,其中IMF3為1 Hz的分量,IMF4和IMF5為5 Hz的分量,即不同的IMF包含相同的頻率成分,仍然

存在模態混疊;本文改進EEMD獲得9個IMF,其中IMF4和IMF5為主要頻率成分,不存在明顯的模態混疊;另外,EEMD方法得到IMF4和IMF5的幅值變化較大,特別是IMF5的幅值變化十分明顯,與原信號的該成分差異顯著。分析認為,本文方法的每次分解均是從原信號的殘量開始,且僅取每次分解的第1階固有模態,從而確保了模態分量的一致性,避免了EEMD中不同分解過程中模態分量不一致的問題。由圖2(b)和(c)可知改進EEMD的每一IMF中均包含結構振動信息,若直接采用圖2(b)中的IMF4和IMF5重構結構信息必然造成結構信息的丟失;MPCA處理后的各IMF噪聲水平分別99.81%,98.49%,97.48%,0.61%,0.09%,96.04%,99.98%,99.01%,其中IMF4和IMF5的噪聲水平最低,分別為0.61%和0.09%,其他IMF分量的噪聲水平很高,均大于96%;圖2(c)的結果還說明了盡管IMF1的噪聲水平很高,但不全是噪聲。

圖2 模擬信號EMD及MPCA結果 Fig.2 EMD results of simulation signal

將分解所得的所有IMF疊加后重構信號,并與原信號進行比較,得到分解誤差,兩種方法的重構誤差如圖 3所示,由圖可知本文方法的誤差在10-14量級,EEMD方法的誤差在10-1量級,顯然本文方法的分解精度遠大于原EEMD方法的分解精度。

圖4示出了兩種方法EMD執行次數箱線圖,本文方法總的迭代次數為45 918,最大迭代次數出現IMF3的計算過程中,為73次;EEMD方法總的次數為59 940,最大迭代次數出現IMF5的計算過程中,為309次。本文方法的計算減少約25%,即本文方法的計算效率更好。

兩種分解方法獲得的瞬時頻率如圖 5所示,為凸現主要頻率成分,圖中瞬時頻率顏色的深淺與該頻率成分能量大小相關,能量越大顏色越深。由圖中可知,兩種EEMD方法均能凸顯模擬信號的頻率成分,但本文改進EEMD方法獲得的瞬時頻率更為連續和清晰,特別是1 Hz瞬時頻率,由圖 2亦可發現該結論。圖5中還示出了MPCA降噪后的瞬時頻率,由圖5(c)可知,降噪后的主要頻率成分比改進EEMD方法的更為清晰和連續。

圖3 模擬信號重構誤差 Fig.3 Reconstruction error of simulation signal

圖4 模擬信號分解迭代次數 Fig.4 Decomposition iteration count of simulation signal

為評判降噪后的重構信號與不含噪聲的模擬信號的差異,計算了MPCA降噪前后的相關性(MPCA降噪前的相關性計算只選擇了改進EEMD得到的IMF4和IMF5),降噪前后的相關性分別為0.986 7和0.986 2。雖然降噪后的相關性略有降低,但其瞬時頻率連續且清晰(如圖5(c)所示)。

圖5 模擬信號分解后的Hilbert-Huang譜 Fig.5 Hilbert-Huang spectrum of simulation signal

圖6 主梁加速度傳感器布置(單位:m) Fig.6 Layout of accelerometer on main girder(Unit: m)

圖7 13#加速度傳感器測試信號分解后的Hilbert-Huang譜 Fig.7 Hilbert-Huang spectrum of testing signal for 13# accelerometer

圖8 數據驅動隨機子空間參數識別穩定圖 Fig.8 Stable graph of data-driven stochastic subspace system identification

3.2橋梁實測動力信號

采用某主跨1 088 m斜拉橋的主梁加速度測試數據對本文方法進行檢驗。在主梁主跨1/6截面和次邊跨1/2截面的上游和下游布置豎向加速度傳感器,共布置14個加速度計,如圖6所示。

用本文方法對所有加速度信號進行分解和重構方法進行驗證。信號采樣頻率為20 Hz,測試時長700 s。分解和重構的情況與模擬信號分解與重構情況類似,改進EEMD的分解誤差與次數均遠低于EEMD分解的相應數據。圖 7示出了13#加速度信號的Hilbert-Huang譜,該圖的情況與模擬信號的Hilbert-Huang譜基本一致;因該橋跨度大,結構柔,模態密集,因此該圖中頻率不像圖5一樣能清晰表征主要成分。為說明該問題,采用數據驅動隨機子空間方法[16]對重構后的信號進行參數識別,識別結果如圖8所示,圖8(a)為原始測試數據的穩定圖,圖8(b)為本文方法處理后的穩定圖,對比這兩個穩定圖可知,用本文方法處理后數據識別的頻率值較多,且穩定軸更為清晰,即本文方法提取的結構信息更為豐富,圖中還標識出了具體的頻率值,與文獻[17]的豎向頻率理論計算值吻合得很好。

4結論

經模擬信號與實橋實測數據的驗證,可得如下結論:

(1)改進后的EEMD方法具有較好的計算效率和分解精度,分解誤差比EEMD小很多,且其Hilbert-Huang譜更為清晰和連續;

(2)多尺度主成分分析能有效地實現固有模態函數的降噪和選擇;降噪后Hilbert-Huang譜清晰和連續;

(3)本文所提方法能對橋梁結構的動力測試信號進行有效的分解和降噪,提取的結構信息更為豐富準確;

(4)實橋測試數據的處理結果表明,本文方法能應用于實際橋梁的動力測試分析中。

參考文獻

[1]Huang N, Attoh-Okine N O. The Hilbert-Huang transform in engineering[M]. Boca Raton, USA:CRC Prese,2005.

[2]Feldman M. Hilbert transform applications in mechanical vibration [M]. John Wiley & Sons, Ltd., 2011.

[3]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.

[4]譚善文,秦樹人,湯寶平. Hilbert-Huang變換的濾波特性及其應用[J]. 重慶大學學報,2004.27(2): 9-12

TAN Shan-wen, QIN Shu-ren, TANG Bao-ping. The filtering character of Hilbert-Huang transform and Its application[J]. Journal of Chongqing University,2004.27(2): 9-12.

[5]Flandrin P, Rilling G, Gonclaves P. Empirical mode decomposition as a filter bank [J]. IEEESig. Pros, 2004(11):112-114.

[6]Gledhill R J. Methods for investigating conformational change in bimolecular simulations[D]. University of Southampton, 2004:201.

[7]Wu Zhao-hua,Huang N E. Ensemble empirical mode decomposition: a noise assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1):1-41.

[8]陳仁祥, 湯寶平, 馬婧華.基于EEMD的振動信號自適應降噪方法[J].振動與沖擊, 2012, 31(15): 82-86.

CHEN Ren-xiang, TANG Bao-ping, MA Jing-hua. Adaptive de-noising method based on ensemble empirical mode decomposition for vibration signal [J]. Journal of Vibration and Shock, 2012, 31(15): 82-86.

[9]陳雋, 李想.運用總體經驗模態分解的疲勞信號降噪方法[J].振動、測試與診斷, 2011, 31(1):15-20.

CHEN Jun, LI Xiang. Application of ensemble empirical mode decomposition to noise reduction of fatigue signal [J]. Journal of Vibrat ion, Measurement & Diagnosis, 2011, 31(1):15-20.

[10]Terrien J, Marque C, Karlsson B. Automatic detection of mode mixing in empirical mode decomposition using non-stationarity detection: application to selecting IMFs of interest and denoising[J]. EURASIP Journal on Advances in Signal Processing, 2011, 2011: 37, doi:10.1186/1687-6180-2011-37.

[11]單德山, 李喬, 付春雨, 等.橋梁健康監測與損傷評估[M].北京:人民交通出版社,2011:25-78.

[12]朱建平. 應用多元統計分析[M].北京: 科學出版社,2006:93-100.

[13]Flandrin P, Gon?alvés P, Rilling G. Detrending and denosing with empirical mode decompositions[C]. Proceedings of the European Signal Processing Conference (EUSIPCO ’04), Aalborg,Denmark ,2004.

[14]Aminghafari M, Cheze N, Poggi J M. Multivariate de-noising using wavelets and principal component analysis [J]. Computational Statistics & Data Analysis, 2006(50)2381-2398.

[15]侯遵澤, 楊文采. 小波多尺度分析應用[M]. 北京:科學出版社, 2012.

[16]單德山,徐敏. 數據驅動隨機子空間算法的橋梁運營模態分析[J]. 橋梁建設,2011(06): 16-21.

SHAN Des-han, XU Min. Operational modal analysis of bridge structure based on data-driven stochastic subspace algorithm [J]. Bridge Construction, 2011(06): 16-21.

[17]陳文元.考慮樁土水耦合的大跨度斜拉橋地震響應與可靠度研究[D].成都:西南交通大學,2013: 30-39. P, Rilling G, Gonclaves P. Empirical mode decomposition as a filter bank [J]. IEEESig. Pros, 2004(11):112-114.

[6]Gledhill R J. Methods for investigating conformational change in bimolecular simulations[D]. University of Southampton, 2004:201.

[7]Wu Zhao-hua,Huang N E. Ensemble empirical mode decomposition: a noise assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1):1-41.

[8]陳仁祥, 湯寶平, 馬婧華.基于EEMD的振動信號自適應降噪方法[J].振動與沖擊, 2012, 31(15): 82-86.

CHEN Ren-xiang, TANG Bao-ping, MA Jing-hua. Adaptive de-noising method based on ensemble empirical mode decomposition for vibration signal [J]. Journal of Vibration and Shock, 2012, 31(15): 82-86.

[9]陳雋, 李想.運用總體經驗模態分解的疲勞信號降噪方法[J].振動、測試與診斷, 2011, 31(1):15-20.

CHEN Jun, LI Xiang. Application of ensemble empirical mode decomposition to noise reduction of fatigue signal [J]. Journal of Vibrat ion, Measurement & Diagnosis, 2011, 31(1):15-20.

[10]Terrien J, Marque C, Karlsson B. Automatic detection of mode mixing in empirical mode decomposition using non-stationarity detection: application to selecting IMFs of interest and denoising[J]. EURASIP Journal on Advances in Signal Processing, 2011, 2011: 37, doi:10.1186/1687-6180-2011-37.

[11]單德山, 李喬, 付春雨, 等.橋梁健康監測與損傷評估[M].北京:人民交通出版社,2011:25-78.

[12]朱建平. 應用多元統計分析[M].北京: 科學出版社,2006:93-100.

[13]Flandrin P, Gon?alvés P, Rilling G. Detrending and denosing with empirical mode decompositions[C]. Proceedings of the European Signal Processing Conference (EUSIPCO ’04), Aalborg,Denmark ,2004.

[14]Aminghafari M, Cheze N, Poggi J M. Multivariate de-noising using wavelets and principal component analysis [J]. Computational Statistics & Data Analysis, 2006(50)2381-2398.

[15]侯遵澤, 楊文采. 小波多尺度分析應用[M]. 北京:科學出版社, 2012.

[16]單德山,徐敏. 數據驅動隨機子空間算法的橋梁運營模態分析[J]. 橋梁建設,2011(06): 16-21.

SHAN Des-han, XU Min. Operational modal analysis of bridge structure based on data-driven stochastic subspace algorithm [J]. Bridge Construction, 2011(06): 16-21.

[17]陳文元.考慮樁土水耦合的大跨度斜拉橋地震響應與可靠度研究[D].成都:西南交通大學,2013: 30-39.

猜你喜歡
橋梁
一種橋梁伸縮縫防滲水裝置
工程與建設(2019年4期)2019-10-10 01:45:56
手拉手 共搭愛的橋梁
句子也需要橋梁
加固技術創新,為橋梁健康保駕護航
中國公路(2017年11期)2017-07-31 17:56:30
無人機在橋梁檢測中的應用
中國公路(2017年10期)2017-07-21 14:02:37
高性能砼在橋梁中的應用
現代鋼橋制造對橋梁鋼的更高要求
焊接(2016年8期)2016-02-27 13:05:15
城鄉建設一體化要注重橋梁的建筑設計
南昌54座橋梁進行兩個月的夏季體檢
橋梁伸縮縫損壞因素與加固
主站蜘蛛池模板: 毛片三级在线观看| 国产成人综合网| 欧美97欧美综合色伦图| 男女男免费视频网站国产| 免费无遮挡AV| www.精品国产| 国产成人无码综合亚洲日韩不卡| 国产十八禁在线观看免费| 美女国产在线| 久久国产精品影院| 中文字幕欧美日韩| 国产女人在线视频| 国产第一福利影院| 91热爆在线| 欧美成人综合在线| 亚洲水蜜桃久久综合网站 | 另类欧美日韩| 一个色综合久久| 国内精品视频| 真实国产乱子伦高清| 青青久久91| 色婷婷狠狠干| 成人在线不卡| 国产美女丝袜高潮| 在线欧美日韩| 欧美在线黄| 国产欧美高清| 91在线一9|永久视频在线| 国产成人高清精品免费软件 | 美女国产在线| 日韩成人在线网站| 国产精品自拍露脸视频| 国产三区二区| 一级成人a毛片免费播放| 国产欧美另类| 亚洲欧美另类中文字幕| 久久香蕉国产线| 精品国产成人国产在线| 亚洲男人天堂2020| 亚洲国产欧美国产综合久久| www.亚洲一区二区三区| 国产99视频精品免费观看9e| 国产精品高清国产三级囯产AV| 依依成人精品无v国产| 国产 在线视频无码| 在线观看91精品国产剧情免费| 国产精品永久在线| 亚洲一级色| 成人在线天堂| 国产欧美高清| 亚洲天堂视频在线免费观看| 天天色天天操综合网| 午夜福利在线观看入口| 国产精品三级专区| 国产高清在线丝袜精品一区| 欧美日韩精品在线播放| 国产精品天干天干在线观看| 视频在线观看一区二区| 国产三级韩国三级理| 久久免费视频6| 强奷白丝美女在线观看| 欧美成一级| 全裸无码专区| 国产毛片不卡| 欧美成人日韩| 国产精品自在自线免费观看| 欧美 亚洲 日韩 国产| 久久精品国产在热久久2019| 9久久伊人精品综合| 少妇精品网站| 欧美一区福利| 成人亚洲视频| 久久综合AV免费观看| 伦伦影院精品一区| 操国产美女| 国产精品妖精视频| av在线5g无码天天| 国产精品妖精视频| 丁香五月激情图片| 精品人妻无码中字系列| 最新精品久久精品| 91无码国产视频|