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

基于層次診斷的水下結(jié)構(gòu)振動噪聲源分離量化

2014-09-05 07:14:26時(shí)勝國于樹華
振動與沖擊 2014年9期
關(guān)鍵詞:振動

時(shí)勝國,于樹華,韓 闖,時(shí) 潔

(1.哈爾濱工程大學(xué) 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室,哈爾濱 150001;2.哈爾濱工程大學(xué) 水聲工程學(xué)院,哈爾濱 150001)

聲隱身化程度是評估現(xiàn)代潛艇總體性能先進(jìn)與否的重要標(biāo)志。潛艇的水下輻射噪聲可被分為機(jī)械結(jié)構(gòu)噪聲、推進(jìn)系統(tǒng)噪聲和流體動力噪聲[1]。潛艇在中低速巡航過程中的主要噪聲源是機(jī)械結(jié)構(gòu)噪聲,在振源處開展減振降噪工作從而控制潛艇的機(jī)械結(jié)構(gòu)噪聲是潛艇聲學(xué)優(yōu)化設(shè)計(jì)過程中的一個(gè)有效環(huán)節(jié)。為了對潛艇機(jī)械設(shè)備噪聲指標(biāo)的分解提供科學(xué)依據(jù),確定各振動噪聲源的貢獻(xiàn)比例是很有必要的。目前,對于相干分析和偏相干分析等傳統(tǒng)噪聲源識別方法已經(jīng)開展了大量的研究工作[2-6]。復(fù)雜噪聲源識別需要處理的是對多種噪聲源識別方法和多個(gè)特征線譜的分離量化結(jié)果進(jìn)行融合與集成,但是由于缺乏一種將相干分析和偏相干分析結(jié)果與定量因素統(tǒng)一處理的方法,傳統(tǒng)噪聲源識別方法對此無法得到合理的結(jié)果。將層次分析法[7]與偏相干分析相結(jié)合的層次診斷應(yīng)用到復(fù)雜機(jī)械系統(tǒng)振動噪聲源分離量化問題中[8-12],能夠?qū)?fù)雜的振動噪聲源分離量化問題表示為有序的遞階層次結(jié)構(gòu),根據(jù)偏相干分析得到的結(jié)果在各特征頻帶內(nèi)對振動噪聲源進(jìn)行排序,從而通過融合計(jì)算實(shí)現(xiàn)頻域內(nèi)的振動噪聲源分離量化。在層次診斷中,通常是根據(jù)振動信號與評價(jià)點(diǎn)信號的相干函數(shù)或偏相干函數(shù)的大小借助于標(biāo)度建立判斷矩陣,但是相干函數(shù)和偏相干函數(shù)不具有能量的物理意義,而且通過兩兩比較對各要素進(jìn)行測度的過程帶有一定的主觀色彩。針對以上問題,本文研究了一種改進(jìn)的水下結(jié)構(gòu)振動噪聲源分離量化方法,該方法以多種噪聲源識別方法為準(zhǔn)則,并借助于具有實(shí)際能量物理意義的標(biāo)度對判斷矩陣元素進(jìn)行量化,從而在線譜上和頻帶內(nèi)實(shí)現(xiàn)了振動噪聲源的分離量化。對各環(huán)節(jié)的計(jì)算結(jié)果采用信息融合技術(shù)進(jìn)行融合得到艙段內(nèi)各機(jī)械設(shè)備在水聲場評價(jià)點(diǎn)處的能量貢獻(xiàn)比例。

1 層次診斷基本原理

1.1 遞階層次結(jié)構(gòu)的建立

根據(jù)層次分析理論,結(jié)合振動噪聲源分離量化的特點(diǎn)建立具有四個(gè)層次的遞階層次結(jié)構(gòu),如圖1所示。

圖1 遞階層次結(jié)構(gòu)模型

其中,目標(biāo)層為振動噪聲源的分離量化結(jié)果,用A表示;準(zhǔn)則層是采用的包括頻帶能量分布分析和偏相干輸出譜分析的特征提取方法,用M表示;頻率層為噪聲評價(jià)點(diǎn)處信號的特征線譜或特征頻帶,用F表示;聲源層為各振動噪聲源,用S表示。

1.2 層次分析法中標(biāo)度的改進(jìn)

在建立遞階層次結(jié)構(gòu)之后,層次分析法借助于合適的標(biāo)度建立具有滿意一致性的判斷矩陣。在水聲工程領(lǐng)域的問題中,測度對象的屬性大多具有明確的物理意義。為了使層次分析法中的標(biāo)度適用于特定的具有物理意義的決策問題,就要對標(biāo)度進(jìn)行調(diào)整使其能夠反映測度對象的物理特性。

假設(shè)偏相干輸出功率計(jì)算得到某頻帶內(nèi)的兩個(gè)振動噪聲源在輸出譜中形成的線性部分分別為L1dB和L2dB,并假設(shè)L1≥L2,那么差值為ΔL=L1dB-L2dB,在這個(gè)頻帶內(nèi)的輸出譜為L=L1dB+L+dB,其中L+為輸出譜的增量。ΔL與L+的關(guān)系如圖2所示。

圖2 分貝數(shù)相加曲線

采用層次診斷對振動噪聲源進(jìn)行分離量化時(shí),需要對各評價(jià)對象的頻帶能量或偏相干輸出譜計(jì)算結(jié)果進(jìn)行兩兩比較,并根據(jù)標(biāo)度建立頻率層與聲源層之間的判斷矩陣。此時(shí),被比較的對象具有明確的物理意義,而層次分析法中的標(biāo)度不具有可以測定這種屬性的物理意義。對于多輸入/輸出模型,在各個(gè)頻率上輸出譜可以寫成各個(gè)輸入的偏相干輸出譜相加的形式,所以將標(biāo)度與分貝數(shù)相加曲線進(jìn)行結(jié)合,得到適用于振動噪聲源分離量化問題的標(biāo)度方法。具體實(shí)現(xiàn)方法為:在某個(gè)頻段內(nèi),如果ΔL=0,則認(rèn)為這兩個(gè)振動噪聲源是同等重要的,其量化值為1;如果ΔL=10 dB,則認(rèn)為與L2相比L1是極端重要的,其量化值為9;對ΔL在0 dB與10 dB之間的輸出譜增量進(jìn)行線性劃分,并在實(shí)數(shù)1~9范圍內(nèi)進(jìn)行量化取值,結(jié)果如表1所示。

表1 1~9標(biāo)度在層次診斷中的應(yīng)用

1.3 D-S證據(jù)理論

設(shè)Ω為識別框架,M條證據(jù)源得到了如下證據(jù):

mi=(mi1mi2…miN),i=1,2,…,M

(1)

這M條證據(jù)由于重要程度不同,將被賦予不同的權(quán)重系數(shù)μi,它們滿足的條件如下:

(2)

在M條證據(jù)中選取兩條進(jìn)行融合計(jì)算,首先需要根據(jù)式(3)將這兩條證據(jù)的全局權(quán)重轉(zhuǎn)化為局部權(quán)重。

(3)

(4)

其中:n=1,2,…,N。為了保留原始證據(jù)的分布形態(tài),在對兩條證據(jù)進(jìn)行調(diào)整時(shí),權(quán)重系數(shù)較大的證據(jù)將被直接保留下來,并對權(quán)重較小的證據(jù)進(jìn)行調(diào)整[13]。

(5)

(6)

對這兩條新證據(jù)進(jìn)行融合后得到的證據(jù)的權(quán)重為μi+μj,然后再進(jìn)一步與其它的證據(jù)進(jìn)行類似的調(diào)整與融合得到最終的結(jié)果。

2 仿真研究

仿真模型如圖3所示。輸入信號x1(t)、x2(t)和x3(t)通過均為一的頻率響應(yīng)函數(shù)H1(f)、H2(f)和H3(f),n(t)為加性高斯白噪聲,輸出信號為y(t)。采樣頻率為8 192 Hz,時(shí)間長度為4 s。

圖3 三輸入單輸出模型

x1(t)=7.5×sin(2πf1t+φ1)+1.5×sin(2πf2t+φ1)+

1.2×sin(2πf3t+φ3)+0.7×randn(1,length(t))

x2(t)=2.5×sin(2πf1t+φ4)+4.5×sin(2πf2t+φ3)+

1.7×sin(2πf3t+φ5)+0.5×randn(1,length(t))

x3(t)=1.6×sin(2πf1t+φ3)+1.2×sin(2πf2t+φ5)+

3.8×sin(2πf3t+φ6)+0.3×randn(1,length(t))

y(t)=x1(t)+x2(t)+x3(t)+n(t)

其中:f1=40 Hz,f2=50 Hz和f3=350 Hz,φ1、φ2、φ3、φ4、φ5和φ6是以隨機(jī)相位的形式出現(xiàn)的。由仿真條件可知,x1(t)對y(t)的貢獻(xiàn)最大,為最主要的振動噪聲源,x2(t)次之,x3(t)的貢獻(xiàn)最小。根據(jù)層次診斷步驟,建立診斷模型如圖4所示。

目標(biāo)層與準(zhǔn)則層之間判斷矩陣建立的依據(jù)是功率譜分析和偏相干輸出譜分析對于振動噪聲源診斷能力的強(qiáng)弱,并按照1~9標(biāo)度確定,如表2所示。

根據(jù)輸出譜的特征線譜可以得到頻率層的元素,并根據(jù)其能量的大小確定準(zhǔn)則層與頻率層之間的判斷矩陣,如表3所示。

圖4 振動噪聲源層次診斷模型

表2 目標(biāo)層與準(zhǔn)則層之間的判斷矩陣

表3 準(zhǔn)則層與頻率層之間的判斷矩陣

對輸入信號進(jìn)行功率譜分析得到其在特征線譜上的能量大小,并進(jìn)行兩兩比較得到以功率譜分析為準(zhǔn)則的頻率層與聲源層之間的判斷矩陣,如表4~表6所示。

表4 在40 Hz特征線譜上頻率層與聲源層之間的判斷矩陣

表5 在50 Hz特征線譜上頻率層與聲源層之間的判斷矩陣

表6 在350 Hz特征線譜上頻率層與聲源層之間的判斷矩陣

在特征線譜上根據(jù)各輸入信號對輸出信號的偏相干輸出譜的大小通過進(jìn)行兩兩比較得到以偏相干輸出譜分析為準(zhǔn)則的頻率層與聲源層之間的判斷矩陣,如表7~表9所示。

表7 在40 Hz特征線譜上頻率層與聲源層之間的判斷矩陣

根據(jù)各層判斷矩陣求出的單準(zhǔn)則下的相對權(quán)重,計(jì)算聲源層各元素對于目標(biāo)層的總排序權(quán)重。

將單準(zhǔn)則下的權(quán)重向量帶入可得:

α=[0.646 8 0.218 8 0.134 4]

3 試驗(yàn)研究

為了驗(yàn)證層次診斷在振動噪聲源分離量化應(yīng)用中的可行性,開展艙段模型振動與噪聲測試試驗(yàn)研究。整個(gè)艙段包括兩個(gè)艙體,分別為Ⅱ艙和Ⅲ艙如圖5所示。激振機(jī)通過基座直接剛性安裝在Ⅱ艙中部的艙底;海水泵通過一套雙層隔振裝置安裝在Ⅲ艙內(nèi)的Ⅱ艙和Ⅲ艙之間的艙壁上;在Ⅲ艙內(nèi)靠近艉部的位置主疏水泵通過隔振器側(cè)掛在浮筏裝置上。為了提取各主要振動噪聲源的特征信息,開展了艙段內(nèi)機(jī)電設(shè)備組合單機(jī)激勵試驗(yàn),主要測量艙段內(nèi)機(jī)械設(shè)備、管路系統(tǒng)、殼體振動信號和水下輻射噪聲。

圖5 艙段模型示意圖

3.1 殼體測點(diǎn)對水聲場評價(jià)點(diǎn)貢獻(xiàn)量計(jì)算

在設(shè)備全開工況下,以殼體測點(diǎn)作為系統(tǒng)輸入,水聲場測點(diǎn)作為評價(jià)點(diǎn)進(jìn)行層次診斷。對艙段模型水下輻射噪聲信號進(jìn)行功率譜分析,結(jié)果如圖6所示。

圖6 水下輻射噪聲頻譜圖

圖7 水下輻射噪聲層次診斷模型

選取水下輻射噪聲的特征線譜或特征頻帶作為頻率層的元素,建立水下輻射噪聲層次診斷模型如圖7所示。

對于目標(biāo)層與準(zhǔn)則層之間判斷矩陣的建立與仿真過程中的方法相同,如表2所示。根據(jù)水下輻射噪聲在特征線譜或特征頻帶內(nèi)的聲壓級的大小,建立準(zhǔn)則層與頻率層的判斷矩陣,如表10所示。

表10 準(zhǔn)則層與頻率層之間的判斷矩陣

根據(jù)S1~S6振動信號的能量大小以及水下輻射噪聲層次診斷模型,得到以頻帶能量分布分析為準(zhǔn)則的頻率層與聲源層之間的判斷矩陣。這里僅列出以F1元素為例的情況,如表11所示。其余在F2~F6元素支配下構(gòu)建判斷矩陣的思路與此類似,不再贅述。

表11 在F1頻帶內(nèi)頻率層與聲源層之間的判斷矩陣

根據(jù)S1~S6振動信號的偏相干輸出譜大小以及水下輻射噪聲層次診斷模型,得到以偏相干輸出譜分析為準(zhǔn)則的頻率層與聲源層之間的判斷矩陣,同樣僅給出以F1元素為例的情況。

表12 在F1頻帶內(nèi)頻率層與聲源層之間的判斷矩陣

以頻帶能量分布分析和偏相干輸出譜分析為準(zhǔn)則的各殼體測點(diǎn)排序權(quán)重如表13所示。

表13 殼體測點(diǎn)排序結(jié)果

根據(jù)單準(zhǔn)則下的殼體測點(diǎn)排序權(quán)重向量,可以得到分別以頻帶能量分布分析和偏相干輸出譜分析為準(zhǔn)則的殼體測點(diǎn)排序結(jié)果β和γ。

根據(jù)單準(zhǔn)則下的權(quán)重向量可以得到各殼體測點(diǎn)對水聲場聲壓信號貢獻(xiàn)的總排序權(quán)重為α。

[0.398 6 0.183 0 0.040 7 0.107 2 0.113 0 0.157 5]

3.2 設(shè)備測點(diǎn)對殼體評價(jià)點(diǎn)貢獻(xiàn)量計(jì)算

同樣是在設(shè)備全開工況下,以各設(shè)備測點(diǎn)作為系統(tǒng)輸入,主疏水泵管路系統(tǒng)進(jìn)水口測點(diǎn)S1作為評價(jià)點(diǎn)為例,根據(jù)層次診斷步驟,首先對主疏水泵管路系統(tǒng)進(jìn)水口測點(diǎn)S1處振動信號進(jìn)行功率譜分析,得到的結(jié)果如圖8所示。

同樣選取特征線譜或特征頻帶作為頻率層元素,建立的殼體振動信號層次診斷模型如圖9所示。

同樣,目標(biāo)層與準(zhǔn)則層之間判斷矩陣如表2所示。根據(jù)主疏水泵管路系統(tǒng)進(jìn)水口振動信號在特征線譜或特征頻帶內(nèi)能量的大小建立準(zhǔn)則層與頻率層之間的判斷矩陣如表14所示。

圖8 主疏水泵管路系統(tǒng)進(jìn)水口振動信號頻譜圖

圖9 殼體振動信號層次診斷模型

表14 準(zhǔn)則層與頻率層之間的判斷矩陣

分別以頻帶能量分布分析和偏相干輸出譜分析為準(zhǔn)則建立頻率層與聲源層之間的判斷矩陣。同樣只列出以F1元素為例的情況,這兩種準(zhǔn)則支配下的判斷矩陣的元素均相同,如表15所示。

表15 在F1頻帶內(nèi)頻率層與聲源層之間的判斷矩陣

以頻帶能量分布分析和偏相干輸出譜分析為準(zhǔn)則的主疏水泵、海水泵和激振機(jī)的排序權(quán)重如表16所示。

表16 設(shè)備測點(diǎn)排序結(jié)果

根據(jù)單準(zhǔn)則下的各機(jī)械設(shè)備排序權(quán)重向量,可以得到分別以頻帶能量分布分析和偏相干輸出譜分析為準(zhǔn)則的機(jī)械設(shè)備排序結(jié)果β和γ。

根據(jù)單準(zhǔn)則下的權(quán)重向量可以得到各機(jī)械設(shè)備對主疏水泵管路系統(tǒng)進(jìn)水口貢獻(xiàn)的總排序權(quán)重為α1。

[0.115 1 0.160 2 0.724 7]

各設(shè)備測點(diǎn)對殼體評價(jià)點(diǎn)的貢獻(xiàn)比例如圖10所示。

3.3 設(shè)備測點(diǎn)對水聲場評價(jià)點(diǎn)貢獻(xiàn)量計(jì)算

將各設(shè)備測點(diǎn)對殼體測點(diǎn)的貢獻(xiàn)比例以及各殼體測點(diǎn)對水聲場評價(jià)點(diǎn)的貢獻(xiàn)比例分別采用加權(quán)平均方法和D-S證據(jù)理論進(jìn)行融合計(jì)算得到各設(shè)備測點(diǎn)對水聲場評價(jià)點(diǎn)的貢獻(xiàn)比例。

首先根據(jù)加權(quán)平均方法對各設(shè)備測點(diǎn)對殼體測點(diǎn)的貢獻(xiàn)比例計(jì)算結(jié)果和各殼體測點(diǎn)對水聲場評價(jià)點(diǎn)的貢獻(xiàn)比例計(jì)算結(jié)果進(jìn)行融合計(jì)算,結(jié)果如圖11所示。

根據(jù)圖10中艙段內(nèi)各設(shè)備測點(diǎn)對主疏水泵管路系統(tǒng)進(jìn)水口測點(diǎn)及其出水口測點(diǎn)的貢獻(xiàn)比例這兩條證據(jù),采用D-S證據(jù)理論進(jìn)行融合。首先根據(jù)式(3)、式(4)和式(5)對這兩條證據(jù)進(jìn)行調(diào)整得到的新證據(jù)為

根據(jù)式(6)得到的這兩條證據(jù)的融合結(jié)果為:

α102=[0.047 3 0.132 6 0.820 1]

該證據(jù)的權(quán)重系數(shù)為0.581 6。然后采用類似的調(diào)整和融合方法可以得到圖10中的六條證據(jù)的推理結(jié)果為激振機(jī)是最主要的振動噪聲源,這與加權(quán)平均方法得到的結(jié)論是一致的。

圖10 各設(shè)備對殼體測點(diǎn)貢獻(xiàn)比例

圖11 各設(shè)備對水聲場評價(jià)點(diǎn)貢獻(xiàn)比例

4 結(jié) 論

根據(jù)層次診斷的基本原理,給出了將頻帶能量分布分析與偏相干輸出譜分析相結(jié)合在特征線譜上和特征頻帶內(nèi)進(jìn)行振動噪聲源分離量化的實(shí)現(xiàn)方法。針對層次分析法中存在的標(biāo)度不能對水聲工程領(lǐng)域中的物理量進(jìn)行定量測度的問題,對判斷矩陣元素的量化方法進(jìn)行了改進(jìn)。在對艙段模型試驗(yàn)數(shù)據(jù)進(jìn)行分析的過程中,分別計(jì)算了各設(shè)備測點(diǎn)對殼體評價(jià)點(diǎn)的貢獻(xiàn)比例以及各殼體測點(diǎn)對水聲場評價(jià)點(diǎn)的貢獻(xiàn)比例,并通過加權(quán)平均方法和D-S證據(jù)理論對這兩個(gè)環(huán)節(jié)的計(jì)算結(jié)果進(jìn)行了融合,分析結(jié)果表明將基于頻帶能量分布分析和偏相干輸出譜分析的層次診斷應(yīng)用到艙段模型振動噪聲源分離量化中是可行的。

[1]王京齊,楊衛(wèi)東,朱春景,等.潛艇低噪聲安靜操縱控制技術(shù)研究[J].武漢理工大學(xué)學(xué)報(bào)(交通科學(xué)與工程版),2011,35(2):257-260.

WANG Jing-qi,YANG Wei-dong,ZHU Chun-jing,et al.Study on submarine's low noise quiet maneuver control technology[J].Journal of Wuhan University of Technology(Transportation Science & Engineering),2011,35(2):257-260.

[2]Trethewey M W,Evensen H A.Development and application of multiple input models for structural noise source identification of forge hammers,part I: development[J].Journal of Acoustical Society of America,1984,75(4): 1092-1098.

[3]Trethewey M W,Evensen H A.Development and application of multiple input models for structural noise source identification of forge hammers,part II: application[J].Journal of Acoustical Society of America,1984,75(4): 1099-1104.

[4]章林柯,江涌,何琳.喬列斯基分解應(yīng)用于偏相干問題的可行性研究[J].船舶力學(xué),2011,15(5):556-562.

ZHANG Lin-ke,JIANG Yong,HE Lin.Study on the feasibility of applying cholesky decomposition to partial coherence[J].Journal of Ship Mechanics,2011,15(5): 556-562.

[5]Davis I,Bennett G J.Experimental investigations of coherence based noise source identification techniques for turbomachinery applications-class and novel techniques[C]// 17th AIAA/CEAS Aeroacoustics Conference(32nd AIAA Aeroacoustics Conference).Portland,Oregon,USA: American Institute of Aeronautics and Astronautics,2011.

[6]Kim H C,Cho M G,Kim J,et al.Coherence technique for noise reduction in rotary compressor[J].Journal of Mechanical Science and Technology,2012,26(7): 2073-2076.

[7]孫宏才,田平,王蓮芬.網(wǎng)絡(luò)層次分析法與決策科學(xué)[M].北京:國防工業(yè)出版社,2011.

[8]黃其柏.復(fù)雜噪聲源層次診斷方法及其在風(fēng)機(jī)中的應(yīng)用研究[J].風(fēng)機(jī)技術(shù),1997(6): 13-16.

HUANG Qi-bai.Complex noise source hierarchial diagnosis method and its application research on fans[J].Compressor Blower & Fan Technology,1997(6): 13-16.

[9]Shu G Q,Liang X Y.Identification of complex diesel engine noise sources based on coherent power spetrum analysis[J].Mechanical System and Signal Processing,2007(21): 405-416.

[10]余桐奎,時(shí)勝國,熊草根.層次分析法在復(fù)雜噪聲源識別中的應(yīng)用[C]//第十三屆船舶水下噪聲學(xué)術(shù)討論會.中國,江西,鷹潭:中國造船工程學(xué)會船舶力學(xué)學(xué)術(shù)委員會、《船舶力學(xué)》編輯部,2011:403-409.

[11]余桐奎.改進(jìn)的復(fù)雜噪聲源識別方法[J].振動與沖擊,2012,31(14):152-156.

YU Tong-kui.Improved method for noise source identification[J].Journal of Vibration and Shock,2012,31(14):152-156.

[12]Wu H P,Lou J J,Liu W W.Identification of noise source based on partial coherence analysis[J].Applied Mechanics and Materials,2012(239-240): 482-486.

[13]孫亮,于雷,黃文卿,等.改進(jìn)加權(quán)D-S證據(jù)理論在目標(biāo)意圖預(yù)測中的應(yīng)用[J].空軍工程大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,10(1):17-22.

SUN Liang,YU Lei,HUANG Wen-qing,et al.Application of Dempster-Shafer evidence theory to target intention prediction[J].Journal of Air Force Engineering University(Natural Science Edition),2009,10(1):17-22.

猜你喜歡
振動
振動的思考
某調(diào)相機(jī)振動異常診斷分析與處理
振動與頻率
This “Singing Highway”plays music
具非線性中立項(xiàng)的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應(yīng)分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發(fā)態(tài)分子的振動-振動馳豫
帶有強(qiáng)迫項(xiàng)的高階差分方程解的振動性
主站蜘蛛池模板: 日韩欧美国产成人| 国产网站黄| 一级成人欧美一区在线观看| 亚洲成人一区在线| 精久久久久无码区中文字幕| 国产免费久久精品99re不卡| 亚洲无码视频图片| 亚洲区视频在线观看| 欧美成人精品欧美一级乱黄| 亚洲视频欧美不卡| 呦女亚洲一区精品| 久久精品丝袜高跟鞋| 国产亚洲精品资源在线26u| 99国产在线视频| 91福利片| 国产日韩欧美精品区性色| 国产毛片不卡| 久久男人视频| 欧美色99| 日韩在线成年视频人网站观看| Aⅴ无码专区在线观看| 男人天堂亚洲天堂| 国产成人免费高清AⅤ| 亚洲天堂视频网站| 精品伊人久久久久7777人| 91伊人国产| 欧美日韩国产高清一区二区三区| 亚洲浓毛av| 色丁丁毛片在线观看| www精品久久| 一区二区三区四区在线| 色135综合网| 成人免费视频一区二区三区 | 日韩小视频在线观看| 久久青青草原亚洲av无码| 国产丝袜一区二区三区视频免下载| 中文字幕无码av专区久久| 日韩精品成人在线| 99久久精品免费观看国产| 亚洲另类第一页| 亚洲免费三区| 妇女自拍偷自拍亚洲精品| 免费jizz在线播放| 国产人成在线视频| 欧美区日韩区| 久久午夜夜伦鲁鲁片无码免费| 国产成人免费高清AⅤ| 91麻豆国产视频| 国产成人一区二区| 国产va免费精品观看| 人妻少妇乱子伦精品无码专区毛片| 欧美成人精品高清在线下载| 国产精品亚欧美一区二区三区 | 91精品免费久久久| 真人高潮娇喘嗯啊在线观看| 又污又黄又无遮挡网站| 女人18毛片一级毛片在线 | 亚洲另类色| 国产一区在线观看无码| 在线观看精品自拍视频| 成人福利视频网| 丰满的熟女一区二区三区l| 精品国产网站| 99久久国产精品无码| 久久9966精品国产免费| 精品伊人久久久久7777人| 九九九九热精品视频| 91久久偷偷做嫩草影院电| 四虎成人精品在永久免费| 久久久无码人妻精品无码| 88av在线播放| 高清无码手机在线观看| 欧美曰批视频免费播放免费| 黄片一区二区三区| 亚洲国产日韩在线成人蜜芽| 亚洲成A人V欧美综合天堂| 亚洲欧美日韩天堂| 72种姿势欧美久久久大黄蕉| 欧美一级视频免费| 国产H片无码不卡在线视频| 无码人妻免费| 国产一二三区在线|