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

基于隨機減量法的非平穩激勵下模態參數識別

2015-05-24 16:14:04黃宗明
振動與沖擊 2015年21期
關鍵詞:模態信號結構

羅 鈞,劉 綱,2,黃宗明,2

(1.重慶大學土木工程學院,重慶 400045;2.山地城鎮建設與新技術教育部重點實驗室,重慶 400045)

基于隨機減量法的非平穩激勵下模態參數識別

羅 鈞1,劉 綱1,2,黃宗明1,2

(1.重慶大學土木工程學院,重慶 400045;2.山地城鎮建設與新技術教育部重點實驗室,重慶 400045)

針對現有結構模態參數識別方法平穩激勵假定的不足,提出了一種非平穩激勵下的結構模態參數識別方法。以白噪聲調幅激勵、調幅相關激勵和頻率幅值調制共三類非平穩隨機激勵為對象,采用演化譜理論推導了隨機激勵分段疊加后的均值和方差表達式,然后在此基礎上對傳統的隨機減量法進行拓展,將非平穩激勵下模態參數識別問題轉化為結構自由振動曲線的模態參數識別,最后利用特征實現算法識別結構的模態參數。通過一個懸臂梁數值模型和兩層實驗室鋼框架進行驗證,識別結果表明運用該算法可以準確識別出非平穩激勵下系統的固有頻率和振型,且固有頻率的識別誤差在2%以內,具有較高的識別精度。

隨機減量法,非平穩環境激勵,模態參數識別,自由振動

橋梁、輸電塔和高層建筑等大型實際結構模態參數的識別,是進行有限元模型修正、損傷識別和優化設計等工作的基礎。傳統的模態參數識別方法多采用人工在結構上施加激勵,是在結構輸入、輸出均已知的情況下進行結構模態參數的識別。但對于大型土木結構而言,結構構件多、幾何尺寸大,很難甚至無法采用人工激勵的方式進行激振,因此,依靠車輛、脈動風、地脈動和地震動等隨機激勵方式(通常無法量測)識別出結構模態參數具有重要的理論意義和現實工程價值,現已發展為土木工程領域內一個重要的研究方向——工作模態分析。

經過30多年的研究,業界已提出了很多的工作模態分析方法,大致可分為基于頻域、時域和時頻域共三類方法。在頻域內,基于白噪聲激勵下結構體系響應的互功率譜估計和奇異值分解,Brincker等[1]提出了頻率域分解法;Guillaume等[2]在激勵和響應均可測量的前提下,基于權矩陣分數模型提出了PolyMax方法;Peeters[3]將PolyMax的應用范圍擴展到白噪聲激勵的情況下,推導了利用輸出功率譜進行模態參數識別的PolyMax方法。在時域內,Van Overschee等[4]將LQ分解和SVD分解引入子空間識別,提出了適用于白噪聲激勵的隨機子空間識別法。在時頻域內,Sun等[5]提出了基于協方差驅動的小波識別方法,由于涉及協方差計算,該方法并不適用于非平穩激勵情況。

為便于理論推導,以上方法均假設隨機激勵是具有時間不相關性的高斯白噪聲,即空間上各點激勵為互不相關的平穩隨機過程。但這一假設與真實激勵并不相符,例如胡強等[6]通過對大地脈動實測數據的分析,表明大地脈動含有豐富的低頻周期分量,具有隨機性和非平穩性。實測脈動風也有一定的主要頻率和方向,且在時域具有明顯的間歇性和非高斯性[7]?,F有工作模態分析方法理論假設的不足,導致識別精度遠小于已知激勵的模態參數識別方法,限制了其在大型土木結構中的應用。

大量研究表明,地脈動、脈動風荷載和地震動等環境激勵都可認為是均值為0,方差隨時間改變的非平穩隨機信號[7-8]。本文以土木結構中以上常用的非平穩隨機激勵為主要研究對象,從演化譜的角度,通過理論推導得出其統計特性,然后對傳統的隨機減量法進行擴展,并利用擴展后的隨機減量法將結構非平穩響應轉化為結構的自由衰減信號,再利用特征系統實現算法(ERA)識別結構的模態參數。通過對白噪聲調幅、相關調幅和頻率幅值調制這三類非平穩激勵下結構模態識別算例表明,所提方法能成功實現非平穩隨機激勵下結構模態參數的識別,從而為大型土木結構模態參數的獲取提供了一種有效的新方法。

1 非平穩隨機荷載的統計特性

地脈動等隨機荷載往往無法準確測量。但在理論分析上,演化譜作為平穩隨機信號譜分解的推廣,廣泛應用于地震動、風荷載等非平穩信號的時頻特性分析和數值模擬[9-10],故以演化譜為工具,研究非平穩隨機荷載的統計特性。

1.1 非平穩隨機荷載的均值和方差

從信號的角度,非平穩隨機荷載可視為非平穩隨機信號。對演化譜而言,該非平穩隨機信號可認為是白噪聲信號通過某線性時變濾波系統所得的輸出[11],例如地脈動可理解為在場地土底部輸入的白噪聲通過場地土后的輸出。當均值為0、方差為的白噪聲e(t)通過一個線性時變濾波系統的輸出為:

式中:f(k)表示k×dt時刻的非平穩隨機荷載;he(k,τ)表示k×dt時刻線性濾波系統脈沖響應函數的第τ個時延值,e(k-τ)表示第(k-τ)×dt時刻的白噪音信號;p為整數,是時延的最大長度值;dt為離散采樣的時間間隔。

如假定線性時變濾波系統為確定性系統,則根據白噪聲信號在時間上具有相互獨立的特性,可得非平穩隨機荷載f(k)的均值和方差為:

1.2 非平穩隨機荷載分段疊加的統計特性

將非平穩隨機荷載f(k)分段疊加后求平均,即將f(k)分為長度為P的n段,然后將各段信號進行平均,即:

式中,kl為第l段數據的起始時間。考察平均后所得信號F的均值和方差:

1.3 信號F在不同荷載形式下的方差特性

從式(4)、(5)可知,將非平穩隨機荷載f(k)分段平均所形成信號F的方差是時間的變量,下面就線性濾波系統的脈沖響應函數he(k,τ)是否與時間相關來分析該方差的特性。

當脈沖響應函數he(k,τ)為時間的獨立函數時,即脈沖響應函數滿足:

式中,A(kl)為調幅函數。此時,非平穩隨機荷載f(k)為白噪聲調幅信號。將式(7)代入式(4)、(5),可得疊加信號F的方差為:

當脈沖響應函數he(k,τ)為時間相關函數時,則所得的非平穩隨機荷載f(k)在各瞬時也是相關的,例如地脈動、脈動風荷載和地震動等激勵均符合這一特性。此時若假設非平穩隨機荷載f(k)分段時,各段信號之間的間隔大于線性濾波系統自由衰減振動的時間常數,則可得:

此時,公式(6)中的第2項將比第1項小得多,可以忽略,則疊加信號F的方差可簡化為:

從式(4)、(9)和(10)可知,無論零均值非平穩隨機荷載為白噪聲調幅信號,還是時間相關信號,經過分段疊加后的信號F的均值為0,方差隨著疊加次數n的增加將逐步減小。

2 基于擴展隨機減量法的模態參數識別

2.1 隨機減量法的擴展

對單自由度結構體系,隨機減量法首先選取某一閾值A,去截取隨機激勵下體系的響應信號,得到起點時刻為tk(k=1,2,…,n),長度為s的n段時間序列信號,如圖1所示。然后再對這n段信號疊加后求平均,從而得出體系的自由振動信號,

當隨機激勵滿足零均值的高斯分布時,理論推導表明該振動信號為結構系統在初始位移情況下的自由衰減振動信號[12]。

下面分析零均值非平穩隨機荷載作用下,隨機減量法所得多自由度體系的自由振動信號δ(s)。對多自由度線性土木結構體系,其運動微分方程可表述為:

式中:M,C和K分別為體系的質量矩陣、阻尼矩陣和剛度矩陣;x為位移向量;f(t)為零均值非平穩的外荷載向量。

圖1 隨機減量法流程示意圖Fig.1 Procedure of the Random decretion technique

由結構動力學基本原理,該線性結構體系第i自由度的位移響應可表示為多個模態位移響應的疊加,即:

式中,φir為第r階振型向量在第i個自由度上的數值;m為模態的階數;Yr為第r階模態位移響應。

按式(11)的方式截取信號進行平均,有:

由式(14)可見,通過線性結構的模態分解,可將求第i自由度上的自由振動信號δi(s)轉化為求各階模態位移響應分段疊加再取平均。

因單位脈沖響應函數均具有衰減性,則可認為gr(t-τ)是有限長度的,同時考慮到模態荷載和模態脈沖響應函數具有單邊性質,則有:

式中,q為模態脈沖響應函數的長度。

將式(15)和式(16)中的t換為tk+s,推導可得Yr(tk+s)的表達式:

式中,Ar,tk,Br,tk分別為:

從式(18)可知,在零均值非平穩荷載作用下,隨機減量法獲得的信號可分為2項的疊加,第1項僅與結構自身及其初始模態位移和模態速度相關,第2項與非平穩外荷載相關。由第1節的推導,F的均值為0,方差隨著隨機減量法平均次數n的增加而趨近于0,故當n較大時,隨機減量法信號δi(s)中第2項的值將趨近于0。

值得注意的是,當外荷載為平穩激勵時,式(18)中的第2項中F的均值為0,方差隨n增大趨近于0,故第2項的值在0附近極小范圍波動;當外荷載為非平穩激勵時,式(18)中的第2項中F均值為0,同時考慮到F的方差隨n增大而逐漸減小,故可將第2項視為噪聲項,則非平穩激勵下采用隨機減量法所得信號為包含噪聲的自由振動衰減信號。本文將隨機減量法從處理白噪音信號擴展到處理非平穩信號,稱為擴展隨機減量法。

2.2 模態參數的識別

為提高非平穩激勵下擴展隨機減量法所得信號的信噪比,本文首先采用奇異值消噪[14]對隨機減量信號進行預處理,再利用文獻[15]提出的改進特征系統實現算法進行土木結構體系模態參數的識別。進行非平穩激勵下模態參數識別的流程如圖2所示。

圖2 非平穩激勵下模態參數識別的流程圖Fig.2 Procedure of themodal parameter identification under non-stationary excitation

3 算例及實驗驗證

3.1 數值算例

以一個懸臂梁模型說明擴展隨機減量法在非平穩激勵下的模態參數識別效果,如圖3所示。梁總長為0.7 m,截面形狀為矩形,寬取為30 mm,高取為3 mm。梁等間距劃分為7段,模型采用模態阻尼假定,各階阻尼比均為0.01。

在固定端分別施加白噪聲平穩激勵、白噪聲調幅激勵、相關調幅激勵。白噪聲調幅激勵和相關調幅激勵采用白噪聲輸入一個線性時變濾波系統的輸出進行模擬。白噪聲調幅激勵模擬時采用的線性時變濾波系統脈沖響應函數為he1(t,τ),相關調幅激勵模擬時采用的線性時變濾波系統脈沖響應函數為he2(t,τ),其表達式為:

圖3 懸臂梁計算模型Fig.3 The cantilever beam model

取圖3中傳感器位置處的加速度響應,采樣頻率為200 Hz,利用極值觸發的隨機減量法提取隨機減量曲線。圖4為參考測點為第2測點時,白噪聲調幅激勵和相關調幅激勵下提取的第4測點隨機減量曲線。圖中右下角的小圖為對應的激勵曲線,橫坐標為時間,單位為s,縱坐標為幅值,單位為m/s2。

值得注意的是,圖4中圖(a)和圖(b)響應衰減規律有差別,這主要是由于自由衰減信號是多頻率成分的組合信號,而在不同激勵下將激起結構不同頻率的信號,則自由振動的衰減信號主要表現為所激振頻率的衰減規律。在相關調幅激勵下,算例結構的自由振動響應以第一階頻率為主;而在白噪聲調幅激勵下,算例結構的響應在各階模態處相差不大,加之一階模態的衰減速率較其余模態階次衰減速率明顯偏小,因此頻率成分較均衡的白噪聲調幅激勵下結構自由振動響應的衰減速率比以一階頻率為主的相關調幅激勵下結構自由振動響應的衰減速率大,衰減快。

采用奇異值消噪技術處理后,截取前面信噪比較高的部分信號,利用特征系統實現算法識別懸臂梁的模態參數,頻率、阻尼比、模態振型、模態能量容差分別取0.01、0.05、0.02、0.1。得到各激勵形式下的識別結果如表1和圖5所示。從識別結果可以看出,本文方法識別的頻率和振型具有較高的精度,頻率的最大相對誤差僅為0.33%,振型的MAC值均為1。

圖4 不同激勵形式下第4測點的隨機減量曲線Fig.4 Random decrement curve of the fourth sensor under various loads

表1 各激勵形式下懸臂梁的模態參數Tab.1 M odal parameters of cantilever beam under various loads

3.2 試驗驗證

采用寬65 mm,厚4 mm的鋼板組成框架的梁和柱,并通過節點板和螺栓連接,每個節點板安裝4顆螺栓,2顆與柱相連,2顆與梁或剛性基座相連,如圖6所示。

圖5 各激勵形式下識別的頻率相對誤差Fig.5 Relative error of frequency identified under various loads

圖6 兩層鋼框架模型Fig.6 2-story steel frame

在梁柱節點處布置兩個加速度傳感器,下部的為測點1,上部的為測點2。為了采集較多的數據點,結構響應的采樣頻率為800 Hz。采用KDJ-50型電磁激振器在第一層右柱下側輸入激振信號。激勵信號為白噪聲調幅激勵、相關信號調幅激勵和地震動激勵。白噪聲調幅激勵和相關信號調幅激勵的構造方式與數值模擬相同,時間間隔取0.001 25 s。地震波選用1976 年7月27日我國唐山地震中密云水庫臺站地震記錄,由于實際地震波的時間間隔為0.01 s,因此本文采用樣條插值對地震波進行處理,使其時間間隔為0.001 25 s。

圖7 不同激勵形式下第2測點的隨機減量曲線Fig.7 Random decrement curve of the 2th sensor under various loads

對采集的原始數據首先通過低通濾波,再利用極值觸發的隨機減量法提取隨機減量曲線。圖7為參考測點為第2測點時,白噪聲調幅激勵、相關調幅激勵和地震激勵下提取的第2測點隨機減量曲線。由于試驗中未測量激振力的大小,激勵的大小由結構響應加以控制,故本圖中右上角給出第2測點的響應信號,橫坐標為時間,單位為s,縱坐標為幅值,單位為m/s2。

采用本文方法識別結構的模態參數,識別結果如表2所示。頻率、阻尼比、模態振型、模態能量容差分別取0.01、0.2、0.02、0.1。由于模態參數的精確值未知,本文將脈沖激勵下識別的模態參數作為精確值來衡量其余激勵形式下的識別精度。從表2可以看出,本文方法可以較準確的識別出結構的頻率和振型。然而阻尼比的識別誤差較大,主要原因在于模態參數識別理論上的瑞雷阻尼假定、各階模態成分的互相干擾、提取隨機減量曲線的平均次數和噪聲的影響。

表2 各激勵形式下識別的模態參數Tab.2 M odal parameters of the steel frame under various loads

4 結 論

本文從演化譜的角度對隨機非平穩激勵荷載的非平穩特性進行分析,然后對傳統的隨機減量法進行擴展,將非平穩激勵下模態參數識別問題轉化為基于包含非平穩噪聲自由振動曲線的模態參數識別,并利用特征系統實現算法(ERA)識別結構的模態參數。最后利用一個懸臂梁數值模型和兩層鋼框架實驗討論了激勵為白噪聲調幅、相關調幅和頻率幅值調制下的模態識別精度,結果表明:

(1)作用于結構的隨機激勵可視為白噪聲通過線性時不變系統的輸出,通過假定線性時不變系統具有的不同特征,得到白噪聲調幅、相關調幅和頻率幅值調制三類具有零均值的非平穩激勵(信號)。這三類信號分段疊加所得信號具有均值為0,方差隨疊加次數增加而線性減小的特征。

(2)通過模態疊加法推導了非平穩激勵下隨機減量法的理論計算公式,結果表明在零均值非平穩隨機激勵下,隨機減量所得結果為包含噪聲的自由振動衰減信號。

(3)本文方法適用于零均值非平穩激勵下工程結構的模態參數識別。數值算例和實驗室模型的結構模態參數識別結果表明懸臂梁數值算例的固有頻率的識別誤差在1%以內,振型的MAC值為1.0;兩層實驗室鋼框架結構的固有頻率識別結果在2%以內,振型的MAC值為1.0。

(4)本文提出方法在實際高層建筑、大跨橋梁結構中的實際應用有待進一步研究。對高層建筑物和橋梁荷載輸入及響應進行實測,研究輸入、輸出信號的特性并進行模態參數識別是本文后續研究的方向。

[1]Brincker R,Zhang L,Andersen P.Modal identification from ambient responses using frequency domain decomposition [C]//Proc.of,2000,18:625-630.

[2]Guillaume P,Verboven P,Vanlanduit S,et al.A polyreference implementation of the least-squares complex frequency-domain estimator[C]//Proceedings of IMAC,2003,21:183-192.

[3]Peeters B,Vanhollebeke F,Van der Auweraer H.Operational PolyMAX for estimating the dynamic properties of a stadium structure during a football game[C]//Proceedings of the IMAC,2005,23.

[4]Van Overschee P,De Moor B.Subspace identification for linear systems:theory,implementation,applications[C]//status:published,1996.

[5]Sun Z,Chang C C.Covariance-driven wavelet technique for structural damage assessment[J].Smart Structures and Systems,2006,2(2):127-140.

[6]胡強,程耀東,齊津.大地脈動數據的分析及建模[J].浙江大學學報,1997,31(6):767-773.

HU Qiang,CHENG Yao-dong,QIJin.Analysis and modeling of earth pulsation data[J].Journal of Zhejiang University,1997,31(6):767-773.

[7]伊廷華,李宏男,王國新.基于小波變換的高層建筑脈動風速模擬與實測研究[J].振動與沖擊,2007,26(3):13-18.

YI Ting-hua,LI Hong-nan,WANG Guo-xin.Expermental study and numerical simulation of fluctuating wind speed on top of a tall building based on wavelet transform[J].Journal of Vibration and Shock,2007,26(3):13-18.

[8]杜修力,陳厚群.地震動隨機模擬及其參數確定方法[J].地震工程與工程振動,1994,14(4):1-5.

DU Xiu-li,CHEN Hou-qun.Random simulation and its parameter determination method of earthquake ground motion [J].Earthquake Engineering and Engineering Vibration,1994,14(4):1-5.

[9]Failla G,Pappatico M,Cundari G A.A wavelet-based spectrum for non-stationary processes[J].Mechanics Research Communications,2011,38(5):361-367.

[10]丁幼亮,周廣東,李萬恒等.基于演化譜理論的橋梁風致響應非平穩性分析[J].中國公路學報,2013,26(5):54-61.

DING You-liang,ZHOU Guang-dong,LIWan-heng,et al.Nonstationary analysis of wind-induced responses of bridges based on evolutionary spectrum theory[J].China Journal of Highway and Transport,2013,26(5):54-61.

[11]王宏禹,邱天爽,陳喆.非平穩隨機信號分析與處理[M].北京:國防工業出版社,2008.

[12]李德葆,陸秋海.實驗模態分析及其應用[M].北京:科學出版社,2001.

[13]克拉夫RW,彭津J著.結構動力學[M].2版.王光遠等,譯.北京:高等教育出版社,2006:177.

[14]孫鑫暉,張令彌,王彤.基于奇異值分解的頻響函數降噪方法[J].振動、測試與診斷,2009,29(3):325-328.

SUN Xin-hui,ZHANG Ling-mi,WANG Tong.Noise reduction of frequency response function using singular value decomposition[J].Journal of Vibration Measurement&Diagnosis,2009,29(3):325-328.

[15]章國穩.環境激勵下結構模態參數自動識別與算法優化[D].重慶:重慶大學,2012.

M odal parametric identification under non-stationary excitation based on random decrementmethod

LUO Jun1,LIU Gang1,2,HUANG Zong-ming1,2

(1.College of Civil Engineering,Chongqing University,Chongqing 400045,China;2.The Key Laboratory of New Technology for Construction of Cities in Mountain Area of the Ministry of Education,Chongqing University,Chongqing 400045,China)

Aiming at the shortcoming of the existingmodal parametric identification method with the assumption of stationary excitations,amodal parametric identification method under non-stationary excitation was presented here based on the radom decrementmethod.Firstly,the variance and the mean of the excitations were derived on the basis of the evolutionary spectral amplitude modulation theory for white noise amplitude modulation excitation,the related excitation with amplitude modulation and the non-stationary excitation with frequency and amplitude modulation.Then,the traditional random decrementmethod was expanded,and themodal parametric identification problem under non-stationary excitation was converted into a modal parametric identification problem based on free vibration response curves of structures.Finally,themodal parameters were identified using the eigen-system realization algorithms.The method was verified through a numerical model of a cantilevered beam and a two-story steel frame tests.The identification results showed that the proposed method can accurately identify natural frequencies and vibration modal shapes of the two structures under non-stationary excitation and the identification error for natural frequencies is less than 2%,it has a higher identification accuracy.

random decrementmethod;non-stationary environment excitation;modal parametric identification;free vibration response

TN911.7

A

10.13465/j.cnki.jvs.2015.21.004

中央高?;究蒲袠I務費(CDJXS11200007)資助

2014-07-14 修改稿收到日期:2014-09-18

羅鈞男,博士生,1986年5月生

黃宗明男,博士,教授,1957年5月生

猜你喜歡
模態信號結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
論《日出》的結構
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
主站蜘蛛池模板: 91青青视频| 色久综合在线| 激情综合激情| 国产精品男人的天堂| 免费在线一区| 一级毛片免费观看不卡视频| 激情六月丁香婷婷四房播| 一级毛片在线播放免费| 国产亚洲高清视频| 黄色一级视频欧美| 在线观看国产网址你懂的| 亚洲国产成人久久精品软件 | 国产午夜在线观看视频| 中文字幕中文字字幕码一二区| 国产精品综合色区在线观看| 色欲色欲久久综合网| 蜜芽国产尤物av尤物在线看| 亚洲第一色网站| 欧美日韩v| 国产xxxxx免费视频| 国产精品一区二区不卡的视频| 国产精品毛片一区视频播| 亚洲欧美国产高清va在线播放| 高清不卡毛片| 国产99视频精品免费观看9e| 国产成人精品一区二区秒拍1o| 国产91无码福利在线 | 国产色图在线观看| 欧美在线一二区| 丁香五月激情图片| 久久久精品无码一二三区| 最新国产精品第1页| 色精品视频| 国产一区免费在线观看| 国产chinese男男gay视频网| 在线播放国产99re| 免费国产高清视频| 免费观看三级毛片| h网址在线观看| 91最新精品视频发布页| 高清无码手机在线观看 | 国产视频a| 国产91蝌蚪窝| 女人18毛片一级毛片在线 | 91精品国产情侣高潮露脸| 亚洲视频二| yjizz视频最新网站在线| 国产一区二区三区日韩精品| 亚洲国产天堂久久综合226114| 亚洲精品自产拍在线观看APP| 无码福利视频| 国产女人18水真多毛片18精品 | 亚洲aaa视频| 亚洲AV无码一区二区三区牲色| 中文字幕欧美日韩高清| 丁香六月综合网| 国产亚洲精久久久久久无码AV| 嫩草国产在线| 狠狠做深爱婷婷综合一区| 日韩最新中文字幕| 亚洲第一成年人网站| 国产午夜一级淫片| 亚洲久悠悠色悠在线播放| 日韩欧美中文| 中文字幕第4页| 99这里只有精品免费视频| 成人日韩视频| 天堂岛国av无码免费无禁网站| 国产精品夜夜嗨视频免费视频| 日韩一区二区三免费高清| 日韩欧美一区在线观看| 毛片免费在线视频| 亚洲欧美精品在线| 五月激情婷婷综合| 72种姿势欧美久久久久大黄蕉| 亚洲熟女偷拍| 国产日韩精品欧美一区灰| 欧美日韩免费观看| 午夜性爽视频男人的天堂| 91亚洲精品国产自在现线| 欧美日本不卡| 日本不卡在线播放|