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

能量加權時間特征用于淺海聲源深度類型判別?

2020-09-29 05:56:56于夢梟周士弘
應用聲學 2020年5期
關鍵詞:深度

于夢梟 周士弘

(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

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

0 引言

聲吶基陣檢測到的海上聲源目標復雜多樣,通常水面聲源目標吃水深度在十幾米以內,而水下聲源目標完全浸水或正常工作深度均在十幾米或幾十米以下,因此,從多樣化的目標類型和復雜變化的特征中脫離出來,解決水面聲源和水下聲源的深度判別問題,對分類和識別水下聲源目標十分重要。

水面及水下聲源判別問題一直是水聲學領域中研究的難題之一。長期以來,聲吶目標識別技術是解決該問題的主要途徑[1?2],該途徑依賴于對目標數據的獲取和訓練,由于目標及其工況的多樣性、海洋環境復雜性、系統自身條件限制以及用于分類識別的目標特征可分性等多方面因素的制約,聲吶目標識別的應用尚存在諸多問題需要解決。從水聲物理角度出發,在特定的海洋環境條件下,尋求利用聲信號場傳播物理特征來判別水面和水下聲源值得探索。

用于水面和水下聲源判別的聲信號場傳播物理特征主要包括聲源深度或者具有深度差異性的能量、相位或干涉結構特征等,在淺海波導中多從簡正模及其干涉結構特征出發進行研究。Shang[3]基于模態濾波技術,使用垂直陣接收聲場信號,將包含有限階簡正模的信號分離并轉化為帶有各階模態本征函數及聲源距離所對應相位信息的向量數據,利用函數正交性及相關算子處理估計聲源深度;Yang[4]同樣使用垂直陣分解簡正模,當實際聲場與拷貝場中各階簡正模幅值相關系數達到最大時,其對應深度即為聲源深度;Premus 等[5?6]在聲線向下折射的淺海波導中使用水平線列陣作為接收陣列,利用多基地聲吶回聲探測技術及模態濾波器,分析簡正模低階子空間的能量差異,從而實現深度類型判別;郭曉樂等[7]利用消頻散變換分離各階簡正模,基于水平波數差與波導不變量的關系實現聲源深度估計;李鵬等[8]使用水平陣接收并在已知目標方位的條件下,在模態域進行波束形成,對各階簡正模的模態強度與理論聲場進行匹配估計聲源深度,其陣元數越多,頻帶選擇越寬,測量精度越高;Conan 等[9]使用水平陣對于簡正模進行模態濾波,將低階簡正模與高階簡正模的能量比作為評判量,判別聲源深度類型,但方法受陣列孔徑影響較大;劉志韜等[10]利用信號自相關函數warping 變換,分離簡正模相干項,通過不同深度聲源激發的聲場能量中,占主導的簡正模相干項特征頻率不同的性質實現水面和水下聲源判別,當頻率較高或海深變深或水文環境變得較為復雜后,方法存在一定的困難。以上研究中大多需要對簡正模進行分解后再進行深度估計或判別。

水平直線陣作為接收陣列,對頻域波束形成后的聲場信號做傅里葉逆變換可以獲得波束強度-時間分布圖。Lee 等[11]基于此提出基陣不變量概念,波束輸出最大值對應波束遷移線(Beamformer migration line)的倒數隨時間沿固定斜率線性偏移,在信號方位固定且海水聲速已知時,該偏移斜率僅與聲源距離有關,即基陣不變量,利用基陣不變量方法可以實現聲源被動測距。由于各階簡正模存在頻散現象,在波束強度-時間分布圖中,接收到的各階簡正模分布在不同的時間段內,而聲源深度不同時,各階簡正模能量的時間分布也將有所不同。基于以上物理思想,本文將研究一種基于水平線列陣接收以及利用能量加權到達時間特征作為特征量的脈沖聲源深度類別判別方法,相比于上述基于簡正模或干涉結構特征過濾的方法而言,力求在水平短陣情況下無需進行簡正模分解即可進行聲源深度類型判別。

1 基本原理

根據簡正模理論,淺海波導中遠程傳播的聲壓場可表示為

式(1)中,S(f)為聲源的頻譜,um(·)為第m階簡正模本征函數(關于深度的實函數),km是第m階簡正模水平波數,M為簡正模數目,zs為聲源深度,zr為接收深度,ρ為海水介質密度,f為頻率。

圖1 為水平直線陣信號接收示意圖,對于N元等間距水平陣接收的信號,經常規波束形成(Conventional beamforming,CBF)后,聲壓可表示為與掃描角θ有關的函數[11],

式(2)中,θ為信號入射方向與陣列法向夾角,θs為聲源真實方位,B(x)=sin(kxNd/2)/sin(kxd/2)為波束指向性函數,s= sinθ,sm= sin?msinθs,?m為第m階簡正模俯仰角,sin?m=krm/k,k= 2πf/c。τn= (n ?1)dsinθ/c,d為陣元間距,c為接收深度處海水聲速,rs為聲源距離。

圖1 水平線列陣信號接收示意圖Fig.1 Schematic diagram of horizontal line array signal reception

對式(2)進行逆傅里葉變換,得到時域表達式:

式(3)中,Re{·}代表實部值,PB+(s,t)為聲壓的復數表達形式。

波束輸出的波束強度-時間分布由式(3)來表示:

波束強度-時間分布圖的橫軸為掃描波束s,縱軸為時間t。

1.1 方法原理

下面通過仿真示例來具體說明聲源深度類型判別方法的基本原理。仿真條件:等聲速剖面,海水聲速為1500 m/s,海深為70 m;單層半無限海底,海底聲速為1700 m/s,密度為1.7 g/cm3,吸收系數為0.2 dB/λ;聲源類型為信道沖擊響應,信號入射方向為60?,聲源與水平線列陣參考陣元距離為5 km。陣列孔徑為441 m,陣元數為64 個,陣元間距為7 m,接收深度為25 m。采用常規波束形成方法獲得陣列輸出。圖2 為聲源頻率為50~100 Hz、聲源深度為5 m 及40 m 時的波束強度-時間分布圖(單位:dB,下同)。

圖2 聲源頻率為50 ~100 Hz,不同聲源深度時的波束強度-時間分布圖Fig.2 Beam-time image with different source depth at 50–100 Hz

圖2中脈沖寬度皆為0.36 s,脈沖寬度確定方法將在1.3 節中作介紹。聲場中,第m階簡正模到達時間為

其中,vgm為第m階簡正模的群速度。而群速度與俯仰角可通過式(6)關聯[11]:

當聲源頻帶足夠寬且聲源方位固定時,各階簡正模的波束遷移線可以連接成一條曲線,即

利用式(7)可解釋圖2 中各時間點能量最大值對應波束的偏移情況。同時,在脈沖寬度內,各個時間點的陣列輸出能量也不同。將脈沖寬度內時間記為ti,波束輸出的聲強最大值為該時刻簡正模貢獻的最大能量,記作Ai,Ai=LB(s,i)。由于水面聲源激發聲場中高階簡正模能量占優,水下聲源激發聲場中低階簡正模能量占優,兩類聲源隨時間的Ai分布有所不同。以Ai作為權值,對時間ti進行加權平均并構造函數如下:

其中,t1為脈沖到達接收陣的初始時刻。脈沖寬度及脈沖到達接收陣的初始時刻的確定將在1.3 節中作介紹。

該函數利用信號到達時間的加權平均值與脈沖到達接收陣初始時刻的差值來表示不同深度聲源的Ai分布,取時間相對量可以減小因聲源距離或聲速誤差帶來的影響。通過水平基陣的基陣不變量被動測距方法[11]或其他被動測距方法[12?17]獲得聲源距離,并在拷貝聲場中對函數H在全海深范圍內進行統計。一般來講,水下聲源激發的聲場中低階簡正模先到達接收器且能量占優,H(zs)較小;水面聲源及近海底聲源激發的聲場中高階簡正模較遲到達接收器且能量占比增大,H(zs)較大。考慮船只吃水深度等因素,在海深不小于幾十米的淺海海域,本文定義深度大致在15 m以上的聲源為水面聲源,15 m 以下的聲源為水下聲源。

1.2 工作頻率選取

對于理想波導,簡正模垂直波數[18]為

式(9)中,m= 1,2,···。各階簡正模本征函數可表示成

研究發現,水面和水下聲源深度類型判別閾值深度與簡正模本征函數在水面下的第一個節點深度有關,并可以此確定聲場中最高階簡正模的階數,從而根據截止頻率公式計算最高工作頻率。由式(9)可知,kzmz?=π.因此,判別閾值對應的深度可近似表示為

反過來,若要以z?≈15 m 作為水面和水下聲源判別的基本依據,那么應滿足

根據截止頻率計算公式,可得到用于判別處理的最高頻率

以上是理想波導中最高工作頻率的計算過程。針對一般波導,已知聲場環境參數后,在拷貝場中計算各階簡正模的本征函數隨深度變化示意圖,可選擇第一節點深度最接近15 m 的簡正模階數作為M,在此也可將海深D帶入式(12)中計算得到M,從而計算工作頻段的最高頻率。此時最高頻率計算公式為

以1.1 節的環境參數為例,此時聲場中的前6階簡正模本征函數歸一化幅度示意圖如圖3 所示。由圖6 可知各階簡正模第一節點深度最接近15 m的是第5 階簡正模,節點深度為15.22 m。將海深D帶入式(12)同樣求得M= 5。根據式(14)的計算f?約為102.46 Hz,故本文選取的工作頻段為50~100 Hz。值得注意的是,在選取工作頻率時不宜過低,使得聲場中僅存1階簡正模,從而導致判別方法性能下降,甚至失效。

考慮到實際應用情況,本文忽略近海底深度范圍(這里選取為海底以上深度15 m 范圍),選擇聲源深度為15 m的H值作為判別閾值,即

忽略的深度范圍可根據實際海上環境進行適當調整。一般情況下,Q值的確定需要考慮聲源距離rs和方位θs以及基陣接收深度zr的影響。最終由實際聲場中計算得到的與拷貝場中的閾值進行對比,判斷目標為水面聲源或水下聲源,即

Q值的確定應在一定頻帶范圍要求下進行,以確保能夠利用波導中低階、高階簡正波的時間擴展及其能量分布特征。

圖3 前6 階簡正模本征函數歸一化幅度示意圖Fig.3 Normalized amplitude schematic diagram of first six order normal mode eigenfunctions

圖4 為聲源距離5 km、 工作頻段分別為50~100 Hz 及70~120 Hz 時,函數H在全海深范圍內隨聲源深度變化情況。在忽略近海底附近范圍后,由圖4(a)可知,深度15 m 聲源H值為0.09984 s,即選取閾值Q= 0.09984 s 后,可通過H值與閾值的比較判別聲源深度;而圖4(b)中,最高頻率超出理論計算值。此時深度15 m 聲源的H為0.0822 s,即閾值Q=0.0822 s,但深度27 m 聲源的H為0.088 s,大于閾值Q,會被判定為水面聲源,即當最高工作頻率高于102.46 Hz 時,會使得方法出現誤判的情況。在此情況下,若因客觀因素必須使用該頻段時,則可以通過選擇水下聲源H的最大值作為閾值解決該問題,閾值深度根據選定閾值映射至15 m附近,此時判別深度類型時出現一定的誤判深度區間,圖4(b)中若選擇0.088 s 作為閾值,其閾值深度映射至14.4 m 附近,出現0.6 m 的誤判深度區間,工作最高頻率越高,誤判區間越大。故應盡量滿足工作最大頻率小于理論計算值以保證方法的有效性。

圖4 不同聲源頻率時,H 值隨聲源深度變化情況Fig.4 Variation of H versus source depths at different frequency

1.3 脈沖寬度及脈沖到達接收陣初始時刻的確定

若脈沖寬度及脈沖到達接收陣的初始時刻誤差較大,計算得到的H值將受到較大影響,所以確定以上兩物理量十分重要。在脈沖寬度內,根據式(7),聲強最大值對應波束的倒數隨時間的偏移線斜率固定,可以根據該性質確定脈沖寬度和脈沖到達接收陣的初始時刻。圖5 為聲源深度5 m 時聲強最大值對應的波束倒數隨時間變化示意圖。

圖5 波束遷移線的倒數隨時間變化示意圖Fig.5 Schematic diagram of the reciprocal of beamformer migration line versus time

圖5 中t1為脈沖到達接收陣的初始時刻,te為脈沖結束時刻。在時間為3.33~3.69 s的區間內,聲強最大值對應波束的倒數產生小幅波動,但可近似認為隨時間沿一定斜率變化。而當時間小于3.33 s及大于3.69 s 的區間內,不再存在這種規律,由此確定t1= 3.33 s,脈沖寬度為0.36 s,時間區間為3.33~3.69 s。

2 性能分析

下面考察在不同聲源距離、方位、海深以及接收深度和陣列孔徑的條件下方法的性能。仿真參數同第1節,聲源頻率為50~100 Hz。

2.1 聲源距離

圖6 為聲源距離分別為5 km、10 km 及20 km時,H函數值隨聲源深度的變化情況。

圖6 不同聲源距離條件下,H 值隨聲源深度變化情況Fig.6 Variation of H versus source depths at different source range

由圖6 可知,不同的聲源距離對應于不同的H曲線,閾值也有所不同。距離變遠時,脈沖寬度增加,各階簡正模到達時間差異變大,且相比于低階簡正模,高階簡正模的能量衰減要大,因而閾值變大。可見,式(14)中Q值的選取應根據聲源距離的變化而改變,例如,圖6 中5 km、10 km 和20 km 時的閾值分別為0.09627 s、0.1712 s、0.2339 s。當給定了聲源距離所對應的合適的Q值,即可實現聲源深度類型判別。

2.2 聲源方位

圖7 為聲源從不同方位入射時,H函數值隨聲源深度的變化情況。

圖7 不同聲源方位條件,H 值隨聲源深度變化情況Fig.7 Variation of H versus source depths at different source azimuth

由圖7可知,對于不同方向入射的信號,H值隨聲源深度變化情況接近,表明閾值Q及閾值深度與聲源方位近似無關。需要指出的是,如果信號從正橫及其鄰近方向入射,入射角的正弦值為0,此時波束強度-時間分布圖中可能無法產生波束遷移線,因此實際應用中無法通過1.3 節方法確定脈沖寬度及脈沖到達接收陣初始時刻,可能導致H值計算出現誤差,造成誤判的情況。而從端射方向入射時也可能會對高階簡正模不敏感,應盡量避免使信號從陣列正橫及端射方向入射。

2.3 海深

驗證在不同的海深時,根據1.2 節的頻率選定標準,方法是否可行。海深為100 m 和200 m 時,本征函數第一節點深度最接近15 m 的簡正模階數分別為第7 階和第14 階,計算最大工作頻率分別為103.59 Hz及107.58 Hz,頻率選擇50~100 Hz。圖8為海深為100 m 及200 m 時,H函數值隨聲源深度的變化情況。

圖8 不同海深條件下,H 值隨聲源深度變化情況Fig.8 Variation of H versus source depths at different sea depth

圖9 陣深70 m、海深200 m 時,H 值隨聲源深度變化情況Fig.9 Variation of H versus source depths at horizontal line array depth of 70 m and sea depth of 200 m

由圖8 中標記黑點為聲源深度15 m 時的H值,在海深為100 m 及200 m 的條件下,在不改變接收器深度時,通過拷貝聲場中H值曲線確定閾值后,可能會導致聲源深度類型誤判,如圖8(b)所示,在深度119 m 附近的聲源會被判定為水面聲源。而將接收陣深度改為70 m 時,如圖9 所示,此時誤判區間消失。這是由于不同的接收深度對各階簡正模的敏感性不同,具體內容會在2.4節中分析。海深的改變不會影響方法的性能。

2.4 水平陣接收深度

在布置水平陣時,要考慮不同接收深度對于各階簡正模敏感性的問題。圖10 為水平陣分別布置于25 m、30 m 和69 m 深時,H函數值隨聲源深度的變化情況。根據圖3 中前5 階簡正模本征函數隨深度變化情況,深度25 m 位于簡正模第3 階節點位置附近,深度30 m非任意階簡正模節點,深度69 m低階簡正模幅度較小,高階簡正模幅度較大。

圖10 不同的水平陣深度,H 函數值隨聲源深度的變化情況Fig.10 Variation of Hversus source depths at different horizontal line array depth

由圖9 可知,當水平陣布置于69 m 時,水下聲源的H值與陣深25 m及30 m時差異較大。由于海底附近位置的低階簡正模本征函數幅值較小,接收陣對于低階簡正模敏感性變弱,第1 階簡正模尤為明顯。水下聲源第1階簡正模的能量占比大,因此對其H值改變明顯。但在已知接收深度和忽略近海底附近范圍的前提下,依舊可以根據Q值,對聲源進行深度類型判別。不同的條件使用不同的閾值,對于水面與水下聲源深度類型的判斷性能影響無影響。方法在陣任意深度皆可使用,但值得一提的是,在布置水平陣的接收深度時,應盡量保證低階簡正模良好的敏感性以獲得水面、水下聲源區別較大的函數H曲線。同理,接收陣對于聲場中激發的最高階簡正模也應保證良好的敏感性,保證脈沖寬度判斷的準確性,有利于獲得較好的聲源深度類型判別性能。這也解釋了2.3 節中海深為200 m、接收深度為25 m時,出現了誤判深度區間的原因。

2.5 陣列孔徑

不改變陣元間距,圖11為32 個陣元,即陣列孔徑為217 m時,聲源深度為5 m及40 m 時的波束強度-時間分布圖。圖12 為32 陣元與64 陣元條件下H函數值隨聲源深度的變化情況的對比。

圖11 陣元數為32 個時,不同聲源深度波束強度- 時間分布圖Fig.11 Beam-time image with different source depth when the number of array elements is 32

將圖11(孔徑217 m)與圖2(孔徑441 m)進行對比可知,當陣列孔徑為217 m 時,波束輸出方位譜的主瓣寬度明顯增加。兩種陣列孔徑條件下,水面聲源對應的H函數值略有差別,而水下聲源對應的H函數值則近似不變;由圖12 可知,陣列孔徑為217 m 條件下,H函數值隨聲源深度的變化情況與陣列孔徑為441 m的H函數值差異很小。這說明盡管接收陣孔徑影響接收陣波束分辨力和空間增益,但對由波束輸出聲強提取的能量加權時間特征隨深度的分布關系而言影響較小。

圖12 不同的陣列孔徑,H 函數值隨聲源深度的變化情況Fig.12 Variation of H versus source depths at different array aperture

2.6 水文類型

下面在負梯度水文環境及溫躍層水文環境中驗證方法的可行性。聲速剖面如圖13所示,其中溫躍層的躍層深度為15~25 m。圖14為兩種水文環境下H函數值隨聲源深度的變化情況。可見,在負梯度及溫躍層水文條件下,方法依然保持了對水面聲源和水下聲源的可判別性。

圖13 負梯度及溫躍層水文環境的聲速剖面Fig.13 Sound speed profile of negative gradient velocity and thermocline

3 失配分析

下面考察聲源距離、海底縱波聲速、海深參數失配時方法的魯棒性。仿真參數均第1 節的模型為基礎,聲源頻率選取50~100 Hz。

3.1 聲源距離失配

聲源距離是本文方法進一步判別聲源深度的前提條件,因為聲源距離會影響到拷貝場計算值的準確性,但實際應用中,聲源距離的估計值會存在誤差,因此有必要討論方法對于聲源距離失配的魯棒性。圖15為聲源距離分別設置為4 km、5 km、6 km時H函數值隨聲源深度的變化情況。

由圖15 可知,不同聲源距離的H值隨深度的變化趨勢接近。隨著聲源距離的增加,低階簡正模與高階簡正模到達接收陣的時間差增加,同時水面聲源中到達接收陣列較晚的高階簡正模能量占優,從而使得水面位置H值隨著距離增加而變大。

距離失配導致實際聲場與拷貝聲場各階簡正模到達時間及脈沖寬度有所差異,但差異量很小,相鄰兩階簡正模的到達時間差在10?2秒級以下,由于H為時間相對量,受小范圍變化的距離影響也非常小,因此較小的聲源距離失配對于判別閾值影響很小,方法對于聲源距離失配具有較好的魯棒性。

圖15 聲源距離小范圍變化時,H 值隨聲源深度變化情況Fig.15 Variation of H versus source depths due to the mismatch of source range

3.2 海底縱波聲速失配

圖16 為海底縱波聲速分別為1650、1700、1750 m/s時H函數值隨聲源深度的變化情況。

圖16 不同海底縱波聲速時H 值隨聲源深度變化情況Fig.16 Variation of H versus source depths due to the mismatch of bottom longitudinal wave velocity

海底縱波聲速的失配將影響聲場中簡正模激發階數及本征函數在深度上的分布。圖17為聲源深度30 m、不同海底縱波聲速條件下的時頻分析圖,由圖17可知,海底縱波聲速為1650 m/s時,在25 m的接收深度下(第3 階簡正模本征函數節點位置附近),聲場中激發兩階簡正模(1 階、2 階),第4 階可以忽略不計,與其他海底縱波聲速環境(1階、2階、4階)有所不同。

圖17 海底縱波聲速為1650 m/s、1700m/s 和1750m/s 時信號時頻分析圖Fig.17 Signal time-frequency analysis image corresponding to bottom longitudinal wave velocity of 1650 m/s,1700 m/s and 1750 m/s respectively

當簡正模總階數變少時,脈沖寬度變短,H值變小。根據圖16中H值隨聲源深度變化情況及2.3節關于海深變化的性能分析可知,聲場中簡正模的階數雖然對方法影響較小,但不同階數情況下,水面與水下聲源的判別閾值卻有所差異,所以,當海底縱波聲速失配導致聲場中出現簡正模階數失配時,方法性能則可能下降。因此在實際應用中應通過底質采樣、反演以及資料查詢等手段盡量保證海底縱波聲速的準確性。

3.3 海深失配

海深的失配同樣會影響本征函數隨深度的分布情況。3.2 節中已提到,若拷貝聲場與實際聲場激發的簡正模階數不同,H值差異較大,會導致方法性能下降,本節僅考慮激發簡正模階數相同時,方法對于海深略有失配時的魯棒性分析。圖18 為海深分別為70 m、72 m 及75 m 時,H函數值隨聲源深度的變化情況,海深70 m。由圖18 可知,在不影響聲場中簡正模階數的情況下,海深略有失配對方法產生的影響可忽略不計。

圖18 不同海深,H 函數值隨聲源深度的變化情況Fig.18 Variation of H versus source depths due to the mismatch of water depth

4 結論

根據簡正模理論,各階簡正模到達接收器的時間有所不同,簡正模階數越低,群速度越大,到達時間越早。利用水平直線陣作為接收陣列,可獲得波束強度-時間分布圖,在脈沖寬度內對每一時刻的陣列輸出最大聲強值(即該時刻信號的能量)進行搜索,以其作為權值,計算各階簡正模到達時間的加權平均值。構造函數H,將信號到達時間加權平均值與脈沖到達接收陣的初始時刻相減,以表示簡正模高低階能量的分布情況。在環境參數及聲源距離已知的前提下,可以獲得拷貝聲場中水面與水下聲源判別閾值,將實際聲場中值與拷貝聲場中判別閾值進行比較,實現寬帶脈沖聲源的深度類型判別。

在具有一般性意義的等聲速淺海波導條件下進行數值仿真表明,以第一節點深度最接近15 m的某階簡正模的截止頻率為上限,給出工作頻率選取準則。無參數失配時,方法對于定義的水面與水下聲源界限方法判別效果好。聲源方位及聲源距離、海深、陣接收深度和陣列孔徑均不影響方法性能。有參數失配時,除海底縱波聲速失配對方法影響較大外,方法對于其他參數(如聲源距離、海深等)的失配魯棒性好,同時方法可在負梯度及溫躍層水文環境下使用。

猜你喜歡
深度
深度理解不等關系
四增四減 深度推進
深度理解一元一次方程
深度觀察
深度觀察
深度觀察
深度觀察
芻議深度報道的深度與“文”度
新聞傳播(2016年10期)2016-09-26 12:14:59
提升深度報道量與質
新聞傳播(2015年10期)2015-07-18 11:05:40
微小提議 深度思考
主站蜘蛛池模板: 日韩成人午夜| 91综合色区亚洲熟妇p| 国产美女一级毛片| 亚洲天堂首页| 久久精品视频亚洲| 日韩无码视频专区| 超碰91免费人妻| 高清不卡一区二区三区香蕉| 亚洲91在线精品| 国产原创演绎剧情有字幕的| 91福利免费视频| 亚洲香蕉在线| 国产日本欧美亚洲精品视| 成人亚洲国产| 亚洲中文久久精品无玛| 在线永久免费观看的毛片| 在线观看国产精品日本不卡网| 伊人久久婷婷| 欧美yw精品日本国产精品| 欧美一级99在线观看国产| 亚洲国产日韩在线观看| 一级毛片不卡片免费观看| 人人爽人人爽人人片| 国产真实乱子伦精品视手机观看| 狠狠综合久久| 亚洲国产理论片在线播放| 怡春院欧美一区二区三区免费| 午夜视频在线观看免费网站| 性69交片免费看| 久久综合成人| 亚洲中文精品人人永久免费| 天堂亚洲网| 亚洲无码91视频| 综合天天色| 欧美在线精品一区二区三区| 人妻无码中文字幕一区二区三区| 久久黄色毛片| 亚洲av无码牛牛影视在线二区| 99视频在线免费| 久久国产精品嫖妓| 99热这里只有精品国产99| 欧美国产精品拍自| 一级毛片网| 精品国产乱码久久久久久一区二区| 992Tv视频国产精品| 国产不卡网| 国产亚洲视频免费播放| 中文字幕免费播放| 影音先锋亚洲无码| 色综合网址| 国产微拍一区| 在线观看无码a∨| 中文字幕精品一区二区三区视频| 无码一区18禁| 亚洲第一成人在线| 无码中字出轨中文人妻中文中| 2020极品精品国产 | 女人av社区男人的天堂| 综合色88| 欧美自慰一级看片免费| 在线国产毛片手机小视频| 久久免费视频播放| 日韩无码白| 久草国产在线观看| 日本精品视频| 久久亚洲天堂| 青青草一区| 在线播放真实国产乱子伦| 成人福利视频网| 成人第一页| 色综合久久久久8天国| 日韩成人午夜| 午夜精品久久久久久久无码软件 | 国产美女人喷水在线观看| 91精品国产一区| 久久久久无码国产精品不卡| 狠狠亚洲五月天| 黄色污网站在线观看| 四虎影视8848永久精品| 国产欧美日韩资源在线观看 | 国产v精品成人免费视频71pao| 五月天福利视频|