亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種模態(tài)參數(shù)識別的虛假模態(tài)剔除技術(shù)*

        2012-01-05 08:14:48王樹青王典鶴林裕裕王巍巍
        關(guān)鍵詞:模態(tài)信號模型

        王樹青,王典鶴,林裕裕,王巍巍,李 慧

        (1.中國海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室,山東青島266100;2.海洋石油工程股份有限公司,天津300451)

        一種模態(tài)參數(shù)識別的虛假模態(tài)剔除技術(shù)*

        王樹青1,王典鶴1,林裕裕1,王巍巍2,李 慧2

        (1.中國海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室,山東青島266100;2.海洋石油工程股份有限公司,天津300451)

        針對背景噪聲下結(jié)構(gòu)的模態(tài)參數(shù)識別結(jié)果存在虛假模態(tài)的情況,分析虛假模態(tài)產(chǎn)生的原因,通過理論推導(dǎo)提出辨別虛假模態(tài)的判斷條件,根據(jù)判斷條件剔除識別模態(tài)參數(shù)中包含的虛假模態(tài)。以比例阻尼海洋平臺結(jié)構(gòu)數(shù)值模型為例,在模擬噪聲水平5%的工況下,獲得結(jié)構(gòu)脈沖響應(yīng)信號,對基于奇異值分解定階消噪后的信號用復(fù)指數(shù)法進(jìn)行模態(tài)參數(shù)識別,對識別的模態(tài)頻率和阻尼比進(jìn)行虛假模態(tài)剔除。結(jié)果表明:根據(jù)判斷條件可有效剔除所識別參數(shù)中的虛假模態(tài)。

        模型定階;信號消噪;復(fù)指數(shù)法;模態(tài)參數(shù)識別;虛假模態(tài)

        模態(tài)參數(shù)識別技術(shù)廣泛應(yīng)用于工程領(lǐng)域,成為解決各類工程問題的重要手段。傳統(tǒng)的模態(tài)參數(shù)識別方法建立在已知系統(tǒng)輸入激勵和輸出響應(yīng)的基礎(chǔ)上,可以根據(jù)輸入輸出信息利用實(shí)驗(yàn)?zāi)B(tài)分析得到系統(tǒng)模態(tài)參數(shù)。但工程實(shí)際中,激勵信號很難獲得,往往只能獲得含噪聲的響應(yīng)信號,因此基于輸出響應(yīng)的模態(tài)參數(shù)識別技術(shù)受到重視。

        近年來,基于輸出響應(yīng)的模態(tài)參數(shù)識別技術(shù)已得到一定發(fā)展。受測量噪聲影響,識別結(jié)果通常包含虛假模態(tài),如何辨識并剔除虛假模態(tài)成為參數(shù)識別過程的關(guān)鍵問題,直接影響了參數(shù)識別技術(shù)在工程實(shí)際中應(yīng)用的效果。目前虛假模態(tài)的剔除主要是從模型定階和信號消噪方面入手,傳統(tǒng)的模型定階方法有利用模態(tài)置信因子定階的方法、曲線擬合法、穩(wěn)定圖法等。Guid De Roeck[1]證明了模態(tài)置信因子有時會由于噪聲的影響而失效;曲線擬合法根據(jù)實(shí)測信號產(chǎn)生數(shù)值模型,進(jìn)而生成模擬信號,由于實(shí)測信號本身存在誤差導(dǎo)致這一方法的使用受到了限制;穩(wěn)定圖表示的是模態(tài)參數(shù)與模型階次的關(guān)系,理論上隨模型階次的增加,真實(shí)的模態(tài)參數(shù)會趨于穩(wěn)定而虛假模態(tài)卻不會,但實(shí)際上識別結(jié)果中真實(shí)模態(tài)參數(shù)也會存在不穩(wěn)定性,導(dǎo)致穩(wěn)定圖法有時也會失效。

        近年來發(fā)展了一些基于奇異值分解的定階消噪技術(shù),易偉建[2]等提出了根據(jù)殘差期望比進(jìn)行模型定階的方法;周幫友[3]等提出根據(jù)奇異值差值進(jìn)行模型定階的方法;針對密集模態(tài)模型,黃應(yīng)來[4]等提出加逆衰減指數(shù)窗與帶通濾波結(jié)合的密集模態(tài)分離方法。王樹青[5-6]等采用歸一化出現(xiàn)次數(shù)的方法來確定模型階次。對于較簡單的結(jié)構(gòu),通過前面的模型定階和信號消噪方法可以有效剔除虛假模態(tài),得到較理想的模態(tài)參數(shù)識別結(jié)果。但當(dāng)結(jié)構(gòu)變得復(fù)雜,通常計(jì)算采用的模型階次高于真實(shí)的模型階次,考慮了噪聲模態(tài)的影響[7],僅通過模型定階和信號消噪過程不能完全剔除虛假模態(tài)。本文提出對模態(tài)參數(shù)識別后進(jìn)行虛假模態(tài)剔除的方法,該方法以比例阻尼結(jié)構(gòu)為模型,基本思想是:比例阻尼模型中模態(tài)頻率與阻尼比之間滿足一定的比例關(guān)系,據(jù)此提出2個判別條件:條件一,各階模態(tài)頻率和阻尼比的乘積呈遞增趨勢;條件二,各階阻尼比與模態(tài)頻率的比值近似相同。通過此判別條件可以判斷所識別模態(tài)頻率和阻尼比的真假。

        1 模型階次的確定和信號消噪

        1.1 Hankel矩陣秩估計(jì)[8]

        用{hl},l=1,…,s表示N自由度系統(tǒng)的脈沖響應(yīng)序列,利用該序列構(gòu)建Hankel矩陣Hm×n

        對矩陣進(jìn)行奇異值分解,有

        式(2)中:矩陣U和V為正交矩陣;∑為對角矩陣,其對角元素為降序排列的奇異值。理論上,信號不受噪聲影響時,超出矩陣秩的奇異值為0,若矩陣的秩為r,則有

        對于受噪聲影響的實(shí)測信號,超出矩陣秩的奇異值并不等于0,但會趨于1個較小值ε,使式(4)成立,可據(jù)此確定模型的階次。

        1.2 基于迭代平均的信號消噪技術(shù)[9]

        對脈沖響應(yīng)信號的降噪處理屬于線性數(shù)學(xué)中的低秩逼近(Low rank approximation)范疇,即對線性時不變系統(tǒng)的脈沖響應(yīng)序列逼近,以達(dá)到Frobenius范數(shù)優(yōu)化準(zhǔn)則。對給定的Hankel矩陣H,要得到滿足Frobenius范數(shù)優(yōu)化的低秩逼近矩陣,可采用截斷奇異值分解法,但得到的結(jié)果矩陣不滿足Hankel矩陣的形式,為維持Hankel矩陣形式,可采用迭代平均的方法,概述如下:

        (1)選定適當(dāng)?shù)闹萺,采用截斷奇異值分解法得到對原Hankel矩陣H低秩逼近的矩陣。此時,不滿足Hankel矩陣的形式。

        將步驟(1)和(2)交替迭代使用,直到滿足收斂標(biāo)準(zhǔn)。

        1.3 模態(tài)參數(shù)識別方法

        模態(tài)參數(shù)識別方法主要分兩類:頻域法和時域法。時域法僅利用輸出的響應(yīng)信號識別系統(tǒng)的模態(tài)參數(shù),因而受到專家學(xué)者的關(guān)注。發(fā)展相對成熟的時域識別方法有復(fù)指數(shù)法、特征系統(tǒng)實(shí)現(xiàn)算法(ERA法)、隨機(jī)子空間法等。本文中采用單參考點(diǎn)復(fù)指數(shù)法進(jìn)行模態(tài)參數(shù)識別,詳細(xì)理論見參考文獻(xiàn)[10]。

        2 虛假模態(tài)的剔除

        虛假模態(tài)產(chǎn)生的原因主要有兩方面:一方面是由于計(jì)算過程中不合理的模型階次估計(jì)導(dǎo)致的;另一方面是由于實(shí)際應(yīng)用中輸入信號不滿足白噪聲的假定和/或輸出信號受到環(huán)境的干擾而導(dǎo)致的[11]?;贖ankel矩陣奇異值分解的模型定階消噪技術(shù)多數(shù)情況下可以實(shí)現(xiàn)合理的模型定階和信號消噪,但在結(jié)構(gòu)較復(fù)雜或噪聲嚴(yán)重的情況下識別結(jié)果中仍可能包含虛假(或準(zhǔn)確度低)的模態(tài)。對于比例阻尼系統(tǒng),通過下面的推導(dǎo)尋找系統(tǒng)模態(tài)參數(shù)(模態(tài)頻率和阻尼比)之間的關(guān)系,并提出剔除虛假模態(tài)的判別條件,根據(jù)判別條件對識別結(jié)果進(jìn)行虛假模態(tài)剔除,從而獲得真實(shí)可信的模態(tài)參數(shù)。具體推導(dǎo)過程如下:

        1個N自由度的動力系統(tǒng)可以用1個二階線性微分方程表示:

        其中:M、C和K∈RN×N分別為質(zhì)量、阻尼和剛度矩陣;分別為位移、速度和加速度向量。系統(tǒng)的第k階特征值表示如下:

        式中:ωk為系統(tǒng)第k階模態(tài)頻率;ζk為第k階模態(tài)阻尼比。

        假定結(jié)構(gòu)滿足Rayleigh阻尼形式,即阻尼矩陣[C]=α[M]+β[K],比例系數(shù)α、β可以用任意已知的兩階振型的阻尼比來確定,如式(7):

        對第k階模態(tài),有

        定義參數(shù)Ak,令

        式中:α、β均為常數(shù),且β>0,對于固有頻率ωk,有

        因此,Ak滿足下式條件

        定義參數(shù)γk,令

        比例系數(shù)α、β∈(0,1),當(dāng)模態(tài)頻率ωk較小時,<α,γk大于1且接近1;隨著模態(tài)頻率值的增大>α,此時可近似認(rèn)為

        此時,有

        令ηk=ζk/ωk,可定義1個合理的系數(shù)區(qū)間(a,b),使

        式(15)中系數(shù)區(qū)間(a,b)的選取可根據(jù)要求的識別精度而定,例如可取區(qū)間為(0.9,1.1),a,b的乘積約為1,如果要求精度較高,則需取1個范圍較小的系數(shù)區(qū)間,相反,擴(kuò)大區(qū)間范圍也就相應(yīng)降低了識別精度,因?yàn)檎`差較大的識別參數(shù)會被作為真實(shí)結(jié)果而輸出。由上述推導(dǎo)可以看出,條件一為模態(tài)頻率和阻尼比的乘積,為遞增關(guān)系;而條件為阻尼比與頻率的比,

        其值對各階模態(tài)近似相等。利用這2個條件就可以判斷識別模態(tài)的虛假。同時,也應(yīng)指出系數(shù)區(qū)間(a,b)的選取依賴于技術(shù)人員的經(jīng)驗(yàn)。剔除虛假模態(tài)的流程如圖1。

        圖1 剔除虛假模態(tài)流程圖Fig.1 Flowchart of false mode elimination

        3 數(shù)值模擬

        3.1 海洋平臺數(shù)值模型

        本文采用海洋平臺結(jié)構(gòu)模型(見圖2)。平臺位于渤海海域,為四樁腿導(dǎo)管架結(jié)構(gòu),就位水深為11.1 m。平臺甲板面積為21 m×16 m,甲板采用板、梁結(jié)構(gòu),Y向設(shè)4根主梁,X向設(shè)4根主梁,主梁為H700工字鋼。導(dǎo)管架設(shè)3層水平橫撐,分別位于-11.1、-4.0和+4.0 m處,在-11.1~-4.0 m處設(shè)豎向斜撐。工作點(diǎn)(17~20點(diǎn))高程為6.6 m。導(dǎo)管架頂(高程EL+5.0m)為13 m×8 m,導(dǎo)管架底(高程EL-12.6 m)為16.84 m×11.84 m,導(dǎo)管尺寸為Φ1 350 m×24 m,斜度為1/10,樁尺寸為Φ1 200 m×26 m,固結(jié)點(diǎn)(45~48)高程取為-21.8 m。水平外圍橫撐尺寸為Φ700 mm×22 mm,水平內(nèi)圍橫撐及斜撐尺寸均為Φ500 mm×18 mm。鋼材密度取7 850 kg/m3,彈性模量為2.1×1011Pa,泊松比0.3。結(jié)構(gòu)采用Rayleigh阻尼(C=αM+βK),考慮到斜向激勵,因此采用第1階和2階頻率來計(jì)算系數(shù)α、β,阻尼比取為0.02,從而得出α、β分別為0.038和0.002。結(jié)構(gòu)固有頻率(前6階):2.199 7、2.393 2、2.614 0、4.079 7、4.366 6、5.480 4 Hz;阻尼比(前6階):0.015 2、0.016 3、0.017 6、0.026 4、0.028 1、0.035 0。結(jié)構(gòu)第2階和第4階模態(tài)分別為Y向第1階和第2階模態(tài),結(jié)果列于表1中。

        圖2 Matlab建立平臺結(jié)構(gòu)數(shù)值模型Fig.2 Platform model produced by MATLAB

        表1 平臺結(jié)構(gòu)Y向前2階模態(tài)頻率(Hz)和阻尼比Table 1 The first two modal frequencies(Hz)and damping ratios of platform in Y direction

        利用Matlab編程建立平臺結(jié)構(gòu)數(shù)值模型,模擬產(chǎn)生單位脈沖激勵,同時作用在29和31結(jié)點(diǎn)Y方向上,獲得平臺結(jié)構(gòu)的脈沖響應(yīng)信號,分別取14、19、20、24結(jié)點(diǎn)的信號進(jìn)行模態(tài)參數(shù)識別。模擬白噪聲信號疊加到脈沖響應(yīng)序列中,噪聲水平通過一個百分比來定量描述,該百分比定義為白噪聲的標(biāo)準(zhǔn)差和精確信號的標(biāo)準(zhǔn)差之比。算例中采樣頻率為200 Hz,信號長度為1 024個點(diǎn)。

        3.2 結(jié)構(gòu)模態(tài)參數(shù)識別及虛假模態(tài)剔除

        分別獲得14、19、20、24結(jié)點(diǎn)Y向含5%噪聲的脈沖響應(yīng)信號。根據(jù)響應(yīng)信號分別構(gòu)建Hankel矩陣得到的奇異值歸一化曲線如圖3,按曲線趨于平穩(wěn)處確定模型階次,則各結(jié)點(diǎn)可確定的模型階次均為4階,可識別的模態(tài)階數(shù)均為2。為避免遺漏模態(tài),并驗(yàn)證所提出方法的有效性,下面按2種情況確定模型階次:第一種情況將各點(diǎn)模型階次確定為6階,則識別結(jié)果中包含3階模態(tài);第二種情況將各點(diǎn)模型階次確定為10階,則識別結(jié)果中包含5階模態(tài)。

        圖3 14、19、20、24結(jié)點(diǎn)信號奇異值歸一化曲線Fig.3 Normalized singular value of node signal at 14,19,20 and 24

        分別按2種定階結(jié)果進(jìn)行信號消噪,用復(fù)指數(shù)法參數(shù)識別得到模態(tài)頻率和阻尼比,按模態(tài)階數(shù)N=3識別的結(jié)果見表2,按模態(tài)階數(shù)N=5識別的結(jié)果見表3。

        當(dāng)按模態(tài)階數(shù)N=3進(jìn)行識別時,以14結(jié)點(diǎn)為例,對Y向響應(yīng)信號的識別結(jié)果進(jìn)行虛假模態(tài)剔除如下:ζ1=0.016 6∈(0.01,0.1),第1階模態(tài)為真實(shí)模態(tài),予以保留;,因?yàn)锳2>A1,滿足條件一,(0.9,1.1)η1,不滿足條件二,所以第2階模態(tài)為虛假模態(tài),予以剔除10-3,A3>A1且η3=0.94η1∈(0.91,1.1)η1,同時滿足條件一和條件二,第3階模態(tài)為真實(shí)模態(tài),予以保留。

        表2中斜體部分的模態(tài)參數(shù)均被作為虛假模態(tài)而剔除,與表1中結(jié)構(gòu)的真實(shí)模態(tài)頻率和阻尼比相比,可以看出剔除虛假模態(tài)后得到了準(zhǔn)確的結(jié)果。

        表2 14、19、20、24結(jié)點(diǎn)信號模態(tài)頻率(Hz)和阻尼比的識別結(jié)果(N=3)Table 2 Identification of modal frequencies(Hz)and damping ratios using response signal of nodes14,19,20 and 24(N=3)

        表3 14、19、20、24結(jié)點(diǎn)信號模態(tài)頻率(Hz)和阻尼比的識別結(jié)果(N=5)Table 3 Identification of modal frequencies(Hz)and damping ratios using response signal of nodes14,19,20 and 24(N=5)

        當(dāng)按模態(tài)階數(shù)N=5進(jìn)行識別時,以14結(jié)點(diǎn)為例,對Y向響應(yīng)信號的識別結(jié)果進(jìn)行虛假模態(tài)剔除如下:ζ1=0.016 7∈(0.01,0.1),第1階模態(tài)為真實(shí)模態(tài),予以保留;同時滿足條件一和條件二,第2階模態(tài)為真實(shí)模態(tài),予以保留;不滿足條件二,第3階模態(tài)為虛假模態(tài),予以剔除;,滿足條件一但不滿足條件二,第4階模態(tài)為虛假模態(tài),予以剔除;不滿足條件一,第4階模態(tài)為虛假模態(tài),予以剔除。

        值得注意的是按N=5進(jìn)行參數(shù)識別時,24結(jié)點(diǎn)信號識別的第2階模態(tài)參數(shù)值被作為虛假模態(tài)剔除,這是因?yàn)闂l件二選取的系數(shù)區(qū)間(0.9,1.1)范圍較小,對識別結(jié)果的精度要求較高,與其它3個結(jié)點(diǎn)的第2階識別結(jié)果相比,可明顯看出24結(jié)點(diǎn)的識別誤差較大,所以被作為虛假模態(tài)而剔除;對24結(jié)點(diǎn)識別的第3階模態(tài)是結(jié)構(gòu)的高階(第17階)模態(tài),該階模態(tài)參數(shù)真實(shí)值為:頻率12.486 1Hz,阻尼比0.078 7。所以按模態(tài)階數(shù)N=5識別,剔除虛假模態(tài)后同樣得到了準(zhǔn)確的結(jié)果。

        4 結(jié)語

        對復(fù)雜的比例阻尼海洋平臺結(jié)構(gòu),根據(jù)結(jié)構(gòu)的模態(tài)頻率與阻尼比之間存在的關(guān)系推導(dǎo)出判別虛假模態(tài)的條件,利用復(fù)指數(shù)法對定階消噪后的脈沖響應(yīng)信號(噪聲水平5%)進(jìn)行模態(tài)參數(shù)識別,然后根據(jù)判別條件對識別結(jié)果進(jìn)行虛假模態(tài)剔除,結(jié)果表明該方法可以有效剔除識別結(jié)果中包含的虛假模態(tài)參數(shù)值。文中僅對比例阻尼系統(tǒng)的單參考點(diǎn)復(fù)指數(shù)法模態(tài)參數(shù)識別進(jìn)行了分析,對采用其它方法識別的參數(shù)進(jìn)行虛假模態(tài)剔除效果還需進(jìn)一步驗(yàn)證。當(dāng)結(jié)構(gòu)阻尼形式發(fā)生改變時,模態(tài)參數(shù)之間的關(guān)系也隨之發(fā)生相應(yīng)改變,但仍可根據(jù)模態(tài)參數(shù)間存在一定關(guān)系的思想,推導(dǎo)出剔除虛假模態(tài)的判斷條件,用于定階消噪?yún)?shù)識別后的虛假模態(tài)剔除。

        [1] Guido De Roeck,Bart Peeters,Wei-Xin Ren.Benchmark study on system identification through ambient vibration measurements[C].San Antonio:Proceedings of 18thIMAC,2000.

        [2] 易偉建,劉翔.動力系統(tǒng)模型階次的確定[J].振動與沖擊,2008,27(11):12-16.

        [3] 周幫友,胡紹全,杜強(qiáng).特征系統(tǒng)實(shí)現(xiàn)算法中的模型定階方法研究[J].科學(xué)技術(shù)與工程,2009,9(10):2715-2722.

        [4] 黃應(yīng)來,董大偉,閆兵.密集模態(tài)分離及其參數(shù)識別方法研究[J].機(jī)械強(qiáng)度,2009,31(1):8-13.

        [5] Wang Shuqing,Li Huajun,Takayama T.Modal identification of offshore platform using statistical method based on ERA[J].China Ocean Engineering,2005,19(2):175-184.

        [6] Wang Shuqing,Liu Fushun.New accuracy indicator to quantify the true and false modes for eigensystem realization algorithm[J].Structural Engineering and Mechanics,2010,34(5):625-634.

        [7] Allenmang R J,Brown D L.A unified matrix polynomial approach to modal identification[J].Journal of Sound and Vibration,1998,211(3):301-322.

        [8] Hu S L J,Bao Xingxian,Li Huajun.Model order determination and noise removal for modal parameter estimation[J].Mechanical System and Signal Processing,2010,24(6):1605-1620.

        [9] Tufts D,Shah A.Estimation of a signal waveform from noisy data using low-rank approximation to a data matrix[J].IEEE Trans Signal Process,1993,41(4):1716-1721.

        [10] 傅志方,華宏星.模態(tài)分析理論與應(yīng)用[M].上海:上海交通大學(xué)出版社,2000.

        [11] 常軍,張啟偉,孫利民.隨機(jī)子空間產(chǎn)生虛假模態(tài)及模態(tài)遺漏的原因分析[J].工程力學(xué),2007,24(11):57-62.

        A Technique of False Mode Elimination in Modal Parameter Identification

        WANG Shu-Qing1,WANG Dian-He1,LIN Yu-Yu1,WANG Wei-Wei2,LI Hui2
        (1.The Key Laboratory of Ocean Engineering of Shandong Province,Ocean University of China,Qingdao 266100,China;2.Offshore Oil Engineering CO,Ltd,Tianjin 300451,China)

        Under noisy background,structural modal parameter identification includes false modes.The causes of false mode in the identified modes were analyzed and two criteria of discriminating false modes from system modes were derived in this paper for eliminating the false modes.Numerical study was conducted of an offshore platform with Rayleigh damping system.First the structural impulse response signal was acquired and a level of 5%noise was added to the simulated signal.Then singular value decomposition was used for noise elimination and Least Square Complex Exponential method(LSCE)was taken to identify the modal parameters.Finally the two discrimination criteria were used for false modes elimination.The results show that the criteria are efficient in eliminating false modes.

        mode order determination;signal de-noising;complex exponential method;modal parameter identification;false mode

        P752

        A

        1672-5174(2012)09-097-05

        國家自然科學(xué)基金項(xiàng)目(50909088;51010009);國家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2008AA092701);教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃項(xiàng)目(NCET-10-0762)資助

        2011-06-23;

        2011-10-12

        王樹青(1975-),男,教授。E-mail:shuqing@ouc.edu.cn

        責(zé)任編輯 陳呈超

        猜你喜歡
        模態(tài)信號模型
        一半模型
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        重要模型『一線三等角』
        完形填空二則
        重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        基于FPGA的多功能信號發(fā)生器的設(shè)計(jì)
        電子制作(2018年11期)2018-08-04 03:25:42
        3D打印中的模型分割與打包
        基于LabVIEW的力加載信號采集與PID控制
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
        国产山东熟女48嗷嗷叫| 亚洲国产综合久久天堂| 亚洲国产精彩中文乱码av| 国产又滑又嫩又白| 国产精品综合久久久久久久免费 | 亚洲中字永久一区二区三区| 人妻久久久一区二区三区蜜臀| 国产精品无码一本二本三本色| 欧美成人三级网站在线观看| 日韩熟女一区二区三区| 久久久中文字幕日韩精品| 亚洲一区二区三区香蕉| 国产又爽又黄的激情精品视频| 国产亚洲一区二区三区夜夜骚| 久久精品国产亚洲av成人文字| 又大又紧又粉嫩18p少妇| 午夜一级韩国欧美日本国产| 一区二区三区四区在线观看视频| 日本一区二区视频在线| 欧美成人片在线观看| 欧美黑人性色黄在线视频| 中文字幕精品亚洲一区二区三区| 激情五月婷婷一区二区| 麻豆果冻传媒在线观看| 亚洲V在线激情| 日本在线一区二区三区视频| 吃奶摸下高潮60分钟免费视频| 国产色综合天天综合网| 中文无码免费在线| 青青草手机在线免费观看视频 | 久久人人妻人人做人人爽| 玩弄人妻奶水无码AV在线| 成人大片免费在线观看视频| 亚洲妇女自偷自偷图片| 免费a级毛片在线观看| 国产av一区二区内射| 无码一区二区三区| 亚洲欧美日韩国产综合一区二区| 国产精品国产三级国产在线观| 久久久天堂国产精品女人| 亚洲码国产精品高潮在线|