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

warping變換提取單模態反演海底衰減系數?

2017-11-12 17:07:56李佳蔚鹿力成郭圣明馬力
物理學報 2017年20期
關鍵詞:模態信號

李佳蔚 鹿力成 郭圣明 馬力

1)(中國科學院聲學研究所,水聲環境特性重點實驗室,北京 100190)

2)(中國科學院大學,北京 100049)

warping變換提取單模態反演海底衰減系數?

李佳蔚1)2)鹿力成1)?郭圣明1)馬力1)

1)(中國科學院聲學研究所,水聲環境特性重點實驗室,北京 100190)

2)(中國科學院大學,北京 100049)

warping變換,模態提取,曲線擬合,海底衰減

為了獲得淺海海底地聲模型參數,利用warping變換方法分離出單模態簡正波.對于接收深度固定、定深爆炸聲源情況,以簡正波理論為基礎定義了距離歸一化的簡正波傳播損失,并且其隨傳播的距離呈線性關系,故可通過此變化規律得到聲壓值實部的衰減因子,進而可求得海底地聲模型參數:海底衰減系數.為驗證此方法的有效性,仿真了warping變換提取單模態簡正波的過程,同時將warping變換提取的單模態簡正波與數值計算的結果進行比較驗證;并針對某次黃海試驗數據進行了處理,得到在150—550 Hz頻帶范圍內海底衰減隨頻率的變化規律為α=0.581f1.86k(dB/m).通過與其他學者在相同海域試驗結果的對比驗證,變化規律基本相同.此外不同模態間反演相同頻點的衰減系數接近也較好地支撐了結果.

1 引 言

海底是海洋環境的一個重要組成部分,是聲波反射和散射的重要邊界,也是影響水下聲傳播的一個主要因素.由于海底邊界的主要影響作用,國內外很多水聲工作者都致力于海底聲學參數的獲取工作,其中反演方法是主要手段之一,其基本思想是對通過攜帶海底信息量的聲場物理量進行最優匹配獲得模型參數.由于對匹配場的不敏感性,海底衰減系數在多參數反演中誤差較大,且衰減系數為頻率的函數,這也對反演的可靠性造成了較大影響.故常在分步反演獲得其他參數的情況下,通過傳播損失擬合[1]、水平縱向相關、測量海底反射系數[2]、模式幅值比等[3,4]來獲得衰減系數.在海底衰減與頻率關系的研究中,在較寬頻帶范圍內,Hamilton總結出海底衰減隨頻率的變化基本呈一次關系;而根據Biot理論,在低頻段較窄的頻帶范圍內,海底衰減隨頻率的變化會呈現出非線性.Holmes總結了獲得海底衰減的主要方法,以及多位學者在100—1000 Hz范圍內砂質海底為主的海底衰減系數隨距離的變化呈非線性關系[5],所得結果與Biot理論更加符合.

簡正波理論可以很好地表征水下聲場,通過有效的消頻散手段能夠提取頻散曲線、分離單模態簡正波等,這些在水聲學研究中具有重要的意義.Zhou[6]利用濾波可分離的模態方法提取第一模式簡正波幅度來反演海底衰減系數.李整林等[7]利用簡正波過濾技術提取出簡正波系數,獲得單模態簡正波反演海底衰減.warping變換是一種處理復雜信號的重要手段,Baraniuk和Jones[8]首次將其應用于信號處理,之后被引入水聲信號處理中.Bonnel等利用warping變換補償信號的頻散特征[9]、通過單聲源進行地聲參數反演[10].Zeng等[11]利用warping變換提取模態簡正波并利用模態簡正波之間的關系來反演海底參數.鹿力成和馬力[12]結合波導不變量,給出了一般意義的warping變換算子.Duan等[13]在獲得聲源信息的情況下,利用warping變換對美國新澤西州海岸的地聲參數進行了反演.此外warping變換也被用于定位[14]、測距等[15,16]諸多方面.

本文結合warping變換分離出單模態簡正波和單模態簡正波聲壓幅值(單位為dB,下同)與距離的關系來獲得海底衰減系數.第一部分推導了淺海環境下以簡正波理論為基礎定義了第m階簡正波的距離歸一化傳播損失,介紹了理想波導的warping變換算子對信號變換的物理本質,為反演衰減系數提供了理論依據;第二部分針對淺海Pekeris波導環境下利用warping變換算子對單模態簡正波的分離進行了仿真,并將分離結果與數值計算結果進行比較;同時通過對2002年某次黃海海試數據的處理得到海底衰減隨頻率的變化關系,驗證了該方法的有效性;最后給出全文的總結.

2 簡正波理論以及warping變換

2.1 簡正波理論

在水平均勻的環境中,單頻點聲源激發聲場可以用一系列簡正波疊加的形式表示

其中

其中,時間因子取 eiωt形式;km,ψm(z)和βm分別為簡正波本征值(或水平波數)、本征函數和衰減系數,一般統稱為簡正波參數,通過求解給定聲速剖面以及海底和海面邊界條件下的簡正波本征方程得到,它們原則上都是頻率f的函數;r是聲源與接收水聽器的水平距離;zs和zr分別是聲源和接收深度;ρ(zs)為聲源深度上的介質密度;M是遠距離聲場有顯著貢獻簡正波的最高階數.定義第m階簡正波的距離歸一化傳播損失

其中

由(2)式可以知道,距離歸一化的簡正波傳播損失Qm(f,r)是隨距離線性變化的,其斜率Km(f)=?8.686βm(f)與簡正波衰減系數直接對應.注意到簡正波衰減系數的微擾近似計算公式為[17]

對于圖1所示的環境模型,認為海底是均勻半無限均勻介質,有

其中,ρw(z),cw(z)和αw(f,z)分別是海水介質密度、聲速和吸收系數剖面;ρb,cb和αb(f)分別是海底介質密度、聲速和衰減系數;D是海深.在頻率比較低的情況下(1 kHz以下),海水聲吸收系數比海底要小幾個數量級,其貢獻βm1可以忽略不計,近似得到

圖1 淺海地聲模型Fig.1.The shallow water acoustic model.

2.2 warping變換

由簡正波傳播損失來提取簡正波衰減系數,進一步再換算得到海底聲衰減系數,關鍵性的工作就是如何從接收到的不同階簡正波彼此疊加的信號中分離出單個簡正波成分.最容易想到的就是對接收信號進行窄帶(例如中心頻率附近1/3倍頻程帶寬)濾波,利用不同階簡正波群速度差異形成的時域波形分離效應,直接提取得到不同階簡正波的時域波形,進一步得到對應的簡正波傳播損失,該方法原則上單通道接收即可,但要求水平距離足夠遠、簡正波群速度差異足夠大,不同簡正波在時域的接收信號波形能夠彼此分開.其次就是利用垂直陣接收,結合簡正波過濾技術來提取不同簡正波分量,該方法盡管不要求簡正波在時域波形上彼此分開,但需要用垂直陣接收,實際應用中容易受到陣形誤差影響,本文利用warping變換來實現不同簡正波成分的分離,進一步提取得到簡正波傳播損失,用于海底衰減系數反演.

warping變換是一種時間軸的非線性變換,對于接收到的多個簡正波疊加的信號波形

其中

s(t)是聲源發射的寬帶脈沖信號,S(f)是對應的信號頻譜,采用時間變換

將接收信號p(t,r)變換成為

以及

式中˙w(t′)=t′/w(t′),t0=r/c0,r是聲源到接收點的水平距離,c0是參考聲速,一般取水中(最小)聲速,t′稱為warping時間.對于接收信號p(t,r)經過warping變換之后得到的信號g(t′,r),不同簡正波在新的頻率域(warping頻率)上是彼此分離的,可以通過濾波得到單個簡正波成分gm(t′,r),再利用warping反變換得到所需的pm(t,r).

2.3 模態提取仿真

為了說明warping變換分離簡正波成分的效果,選取具有兩層均勻分層介質的Pekeris模型進行數值計算,其中海深為30 m,海水和海底介質聲速分別為1500和1600 m/s,密度比為1/1.6;聲源深度為7 m,接收深度為29 m,水平距離為10 km;海底介質聲衰減系數為αb(f)=0.3f1.9kdB/m,fk=f/1000是以kHz為單位的信號頻率;發射信號s(t)選取為頻帶100—400 Hz的寬帶高斯信號.

圖2給出了接收信號p(t)對應的時域波形、頻譜結構和時頻分析結果及其warping變換信號g(t)對應的時域波形、頻譜結構和時頻分析結果.從圖中可以看出原始的接收信號時域波形,能夠大致分辨出對應前三階簡正波的三部分信號波形(100 Hz時,應該有三階波導簡正波,400 Hz時應該有六階波導簡正波,更高階的簡正波由于聲源級和激發強度等方面的原因,并沒有明顯觀察到).在三階簡正波對應的三部分波形結構中,可以看到頻散效應導致的不同的頻率成分是不同時間到達的(高頻先到,低頻后到,對應于上排時頻圖中三條亮條紋結構).由于簡正波之間的群速度差異和水平距離不甚合適、簡正波頻散引起的時域擴展等原因,很難從時域波形結構上(圖2(a))將三階簡正波分開,也無法從頻域結構上(圖2(b))進行分離,但warping變換后的信號則有明顯改善.在時域波形結構上,盡管三階簡正波對應的三部分信號波形在時域上也沒有完全分開(圖2(d)),但變成三個單頻成分,分別對應前三階簡正波在理想波導的截止頻率[5](0.5mcw/D):25,50,75 Hz,頻譜上完全分離開了(圖2(e)). 這樣就可以對warping信號g(t′,r)進行頻率濾波處理,得到不同的簡正波對應的信號波形gm(t′,r),然后再反變換得到其原始的接收信號波形pm(t,r)以及對應的信號頻譜Pm(f,r).

圖3給出了利用warping變換提取得到的簡正波接收信號波形與直接數值計算得到的結果比較.圖3(a)—圖3(c)依次為前三階簡正波在warping頻率域濾波以后的分離結果;圖3(d)—圖3(f)依次為前三階簡正波通過warping變換分離得到的信號時域波形(虛線)及其數值計算結果(實線).從圖中可以看到前三階簡正波被正確地提取出來,得到的結果與直接數值計算的結果完全一致.注意圖中第三階簡正波的接收信號波形,由于它在圖1所示的Pekeris波導中的截止頻率約為173.8 Hz,故低于100 Hz的頻段頻譜被波導截止了.

圖2 (網刊彩色)(a)原始接收信號;(b)原始接收信號的頻譜圖;(c)原始接收信號的時頻圖;(d)warping變換的信號;(e)warping變換信號的頻譜圖;(f)warping變換信號的時頻圖Fig.2.(color online)(a)Original signal;(b)the spectrum of original signal;(c)spectrogram of the original signal;(d)warped signal;(e)the spectrum of the warped signal;(f)spectrogram of the warped signal.

圖3 (網刊彩色)warping變換通過窄帶濾波提取模態簡正波以及提取結果與數值計算在頻譜的對比圖Fig.3.(color online)Normal modes extracted from the warped signal by narrow band fi lter and the contrast between extracted result and numerical calculation in spectrum.

3 試驗與結果分析

海上實驗數據來自于2002年冬季黃海北部進行的一次聲傳播試驗.試驗海區海深約為30 m,海底平坦,沉積層厚度超過10 m以上,對于數據分析的100—1000 Hz頻段范圍內,可以等效近似為均勻半無限流體介質海底;接收船拋錨,垂直接收陣船舷吊放入水,水中共29個陣元,陣元間距1 m,分布在1—29 m深度范圍;接收信號經過濾波放大后用(RSR512)磁帶錄音機記錄,采樣頻率10 kHz;發射船沿設定航線航行,航速約12節,沿途投放7 m定深、38 g裝藥量的爆炸聲源,前50枚(距離約11 km以內)間隔30 s投放一枚,對應間隔約為200 m,后面50枚間隔1 min投放一枚,對應距離間隔約400 m,這樣最遠距離約31 km.試驗期間海況良好.圖4給出了測量得到的聲速剖面以及發射、接收布置示意圖.從圖中可以看到,聲速剖面上下差異不大,表層6 m聲速約為1477.3 m/s,16 m以下聲速約為1475.1 m/s,6—16 m深度范圍是一個弱躍層,聲速差約為2.2 m/s.

圖4 發射接收布置示意圖以及測量聲速剖面Fig.4.Diagram of the transmitter and receiver system and the sound velocity pro fi le.

圖5 (網刊彩色)(a)原始接收信號;(b)原始接收信號的頻譜圖;(c)原始接收信號的時頻圖;(d)warping變換的信號;(e)warping變換信號的頻譜圖;(f)warping變換信號的時頻圖;(g)第一階簡正波頻譜圖;(h)第二階簡正波頻譜圖;(i)第三階簡正波頻譜圖Fig.5.(color online)(a)Original signal;(b)the spectrum of original signal;(c)spectrogram of the original signal;(d)warped signal;(e)the spectrum of the warped signal;(f)spectrogram of the warped signal;(g)the spectrum of the fi rst mode wave;(h)the spectrum of the second mode wave;(i)the spectrum of the third mode wave.

選取接收陣最下面一個陣元接收的7個距離點上的爆炸聲信號進行數據處理,這樣接收深度約為29 m,聲源深度約為7 m,距離點分別為1.735,3,4.16,5.31,7.78,9.1,10.25 km;圖5給出了4.16 km距離上的爆炸聲接收信號的分析處理結果.可以看到,在800 Hz頻段范圍內,大致有5階明顯的簡正波存在,已經很難將不同階的簡正波進行分離,在經過濾波以及反變換以后,得到了所需要的單個簡正波成分Pm(f,r)(前三階簡正波結果見圖5(g)—圖5(i)).

獲得單個簡正波成分以后進一步計算出(2)式對應的距離歸一化的簡正波傳播損失;對7個距離點的處理結果進行線性擬合,得到對應的斜率Km(f)以及簡正波衰減系數βm(f),結果如圖6—圖8所示,用于后續的海底衰減系數反演.從圖中結果可以看到,簡正波衰減系數(或者簡正波傳播損失斜率),隨著簡正波階數和頻率是逐漸增大的,這與“海底聲衰減隨頻率增加以及高階簡正波海底反射損失增加”的物理機理是符合的.

圖6 模態一擬合所得斜率Fig.6.Slope of the fi rst normal mode.

表1給出了根據(5)式、由前三階簡正波衰減系數βm(f)換算得到的海底衰減系數αb(f)結果,其中海底介質密度和聲速采用文獻[6]給出的數據.從表中可以看到,海底衰減系數隨頻率的增加而增大,并且不同簡正波得到的結果也相當一致.

為了進一步分析海底衰減系數αb(f)的頻率變化關系,對得到的數據結果按照關系式αb(f)=α0dB/m,fk=f/1000進行了擬合,在150—550 Hz頻段范圍內,擬合結果為α0=0.581,n=1.86,如圖9所示.

圖7 模態二擬合所得斜率Fig.7.Slope of the second normal mode.

圖8 模態三擬合所得斜率Fig.8.Slope of the third normal mode.

表1 海底衰減反演結果Table 1.Invert results of seabed attenuation.

從處理數據的過程分析,試驗的誤差來自于兩方面:海洋波導參數的偏差、海洋噪聲對信號的干擾,以及原始信號的預處理等對信號造成的影響.爆炸聲接收信號包含直達波和氣泡脈動,而直達波與一次脈動的時間間隔較短,使得選擇的信號長度有限.且隨著信號傳播距離增加,有限長度的直達波由于頻散效應其低頻部分會有部分能量流失.故利用warping變換提取單模態簡正波可使用的距離有限,而簡正波距離歸一化的傳播損失其斜率對距離是敏感的,可用數據在距離上太短會造成數據處理的誤差.鑒于這些因素,單模態簡正波中信噪比高、較為穩健的頻帶較窄.

圖9 海底衰減擬合圖Fig.9.Fitting of seabed attenuation.

4 結 論

在海底地聲反演問題中,通過建立海底地聲模型反演所得參數的本質為在頻帶范圍內真實海底的等效模型參數.在本研究中,利用warping變換較好地分離出前三階簡正波,由簡正波幅值與距離的關系獲得半無限海底地聲模型下的海底衰減.通過數據處理和分析,給出在150—550 Hz范圍內的衰減隨頻率的變化關系為α=在頻點500,550 Hz處,由第二模態和第三模態所得海底衰減系數相近.圖10為本文反演結果與文獻[6]海域A、海域B以及文獻[11]中反演結果的對比圖.試驗海區均為黃海海域,其變化規律基本一致.

圖10 (網刊彩色)海底衰減隨頻率的變化關系Fig.10.(color online)Seabed attenuation as a function of frequency.

[1]Peng Z H,Zhou J X 2004IEEE J.Oceanic Engineer.29 4

[2]Jiang Y M,Chapman N R 2009J.Acoust.Soc.Am.125 4

[3]Tindle C T 1982J.Acoust.Soc.Am.71 5

[4]Potty G R,Miller J H,Lynch J F 2003J.Acoust.Soc.Am.114 4

[5]Holmes J D,Carey W M,Dediu S M,Siegmann W L 2007J.Acoust.Soc.Am.121 5

[6]Zhou J X 1985J.Acoust.Soc.Am.78 3

[7]Li Z L,Yan J,Li F H,Guo L H 2002Acta Acoust.27 6(in Chinese)[李整林,鄢錦,李風華,郭良浩2002聲學學報27 6]

[8]Baraniuk R G,Jones D L 1995IEEE Trans.Sign.Proc.43 2269

[9]Bonnel J,Le Touzé G,Nicolas B,Mars J I 2013IEEE Signal Proc.Magazine30 120

[10]Bonnel J,Barbara N 2010J.Acoust.Soc.Am.128 719

[11]Zeng J,Chapman N R,Bonnel J 2013J.Acoust.Soc.Am.134 EL394

[12]Lu L C,Ma L 2015Acta Phys.Sin.64 024305(in Chinese)[鹿力成,馬力 2015物理學報 64 024305]

[13]Duan R,Chapman N R,Yang K D,Ma Y L 2016J.Acoust.Soc.Am.139 70

[14]Yao M J,Lu L C,Ma L,Guo S M 2016Acta Acoust.41 1(in Chinese)[姚美娟,鹿力成,馬力,郭圣明2016聲學學報41 1]

[15]Qi Y B,Zhou S H,Zhang R H,Zhang B,Ren Y 2014Acta Phys.Sin.63 044303(in Chinese)[戚聿波 2014物理學報63 044303]

[16]Wang D,Guo L H,Liu J J,Qi Y B 2016Acta Phys.Sin.65 104302(in Chinese)[王冬,郭良浩,劉建軍,戚聿波2016物理學報65 104302]

[17]Jensen F B,Kuperman W A,Porter M B,Schmidt H 1992Computational Ocean Acoustics(New York:Springer)pp385–389

Inversion of seabed attenuation by using single mode extracted by warping transform?

Li Jia-Wei1)2)Lu Li-Cheng1)?Guo Sheng-Ming1)Ma Li1)

1)(Key Laboratory of Underwater Acoustic Environment,Institute of Acoustics,Chinese Acdemy of Sciences,Beijing 100190,China)
2)(University of Chinese Academy of Sciences,Beijing 100049,China)

6 April 2017;revised manuscript

24 June 2017)

Seabed is an important part of the marine environment and it has a signi fi cant in fl uence on sound propagation.Considering the fact that geoacoustic parameters are directly acquired with difficulty and complexity,a lot of researchers have focused on the inversion of them.The seabed attenuation coefficient is insensitivity to the matching fi eld.However it has great e ff ects on the transmission loss,mode amplitude ratios,etc.It can be inverted from measurements of these quantities.In this paper,we present an inversion scheme based on warping transform technique for estimating the seabed attenuation coefficient.It utilizes an equivalent seabed model which is constructed by using a prior and posterior knowledge.The dispersion characteristics of normal modes can be observed using the time-frequency analysis of the explosive signal recorded.The dispersion curve can be used to invert the seabed sound speed and density.The results presented by other scholars in the same circle are cited in this paper that focuses on how to obtain the seabed attenuation.Warping transform technique is used to separate and extract the normal modes.The main advantage of warping transform is that it can transform the time-frequency spectrogram into linear relationship which makes it easier to extract the normal modes.The feature of this paper lies in determining the distance normalized normal mode transmission loss.If the depths of receiving hydrophone and the explosion source are constant,the plot of normalized normal mode transmission loss versus distance is a straight line from the normal modes theory,which can be used to obtain the attenuation factor of real part of pressure.Then the seabed attenuation coefficient of the shallow water acoustic model can be calculated.In order to verify the e ff ectiveness of this method,the warping transformation technology is used to separate and extract the fi rst three modes from the simulated Gaussian pulse signal which is obtained in a simulated environment which is similar to the real marine environment.The extracted results are completely consistent with the numerical results.After that,the impulsive signal data collected in the Yellow Sea are analyzed according to the scheme process,and the relationship between the seabed attenuation and frequency isα=0.581f1.86k(dB/m)in a range from 150 Hz to 550 Hz.The results are in good agreement with those obtained by other scholars in the same circle.On the other hand,the inversion results of seabed attenuation from different modes can be used for comparison at the same frequency,which can be a good support for the result.

warping transform,extract the mode,matched curve,seabed attenuation

(2017年4月6日收到;2017年6月24日收到修改稿)

10.7498/aps.66.204301

?國家自然科學基金(批準號:11004214,11274338)資助的課題.

?通信作者.E-mail:luce_1983@sina.com

?2017中國物理學會Chinese Physical Society

http://wulixb.iphy.ac.cn

PACS:43.30.Bp,43.30.Pc,43.60.–c,43.60.PtDOI:10.7498/aps.66.204301

*Project supported by the National Natural Science Foundation of China(Grant Nos.11004214,11274338).

?Corresponding author.E-mail:luce_1983@sina.com

猜你喜歡
模態信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
孩子停止長個的信號
車輛CAE分析中自由模態和約束模態的應用與對比
基于LabVIEW的力加載信號采集與PID控制
國內多模態教學研究回顧與展望
一種基于極大似然估計的信號盲抽取算法
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 国产资源站| 2020国产精品视频| 日韩免费中文字幕| 三级欧美在线| 精品丝袜美腿国产一区| 色网站在线免费观看| 91精品啪在线观看国产60岁| www.狠狠| 在线国产综合一区二区三区 | 国产成人综合日韩精品无码首页 | 婷婷久久综合九色综合88| 色妞永久免费视频| 91亚洲精选| 亚洲成在线观看| 国产精品第一区| 国产精品刺激对白在线| 欧美伦理一区| 一级高清毛片免费a级高清毛片| 亚洲综合一区国产精品| 老熟妇喷水一区二区三区| 亚洲女同一区二区| 免费99精品国产自在现线| 露脸国产精品自产在线播| 无码精品国产dvd在线观看9久| 国产69精品久久久久孕妇大杂乱 | 亚洲aⅴ天堂| 国产在线视频导航| 国产欧美视频综合二区| 老司机精品99在线播放| 久久a级片| 欧美精品v欧洲精品| 国产亚洲欧美另类一区二区| 国产 在线视频无码| 色综合五月婷婷| 国产中文在线亚洲精品官网| 99手机在线视频| 中文字幕亚洲电影| 特级欧美视频aaaaaa| 国产精品尹人在线观看| 亚国产欧美在线人成| 日韩欧美国产三级| 日韩 欧美 国产 精品 综合| 欧美一区中文字幕| 国产日本欧美亚洲精品视| 欧美成人手机在线观看网址| 国产色婷婷| 国产欧美视频一区二区三区| 72种姿势欧美久久久大黄蕉| 婷婷成人综合| 国产精品一区在线观看你懂的| 好久久免费视频高清| AV在线天堂进入| 国产三级成人| 国产一级片网址| 午夜视频免费一区二区在线看| 伊人91在线| 伊人成人在线视频| 日韩久草视频| 中文字幕在线播放不卡| 中文一区二区视频| 日韩午夜福利在线观看| 久996视频精品免费观看| 亚洲视频黄| 一本色道久久88综合日韩精品| 国产精品入口麻豆| 97精品伊人久久大香线蕉| 久久综合国产乱子免费| 毛片免费网址| 久久黄色小视频| 免费国产高清视频| 高清无码手机在线观看| 亚洲开心婷婷中文字幕| 国产xx在线观看| 精品国产一区二区三区在线观看| 国产成人久视频免费| 四虎在线观看视频高清无码| 老色鬼欧美精品| 日本免费福利视频| 精品欧美视频| 四虎永久免费地址| 国产成本人片免费a∨短片| 精品国产中文一级毛片在线看|