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

高階聲場傳感器的噪聲空間相關性及其陣列超指向性方法

2020-12-15 02:29:58朱少豪汪勇楊益新
哈爾濱工程大學學報 2020年10期
關鍵詞:方向

朱少豪,汪勇,楊益新

(西北工業大學 航海學院,陜西 西安 710072)

隨著聲隱身技術的發展,潛艇等水下航行器的噪聲越來越小,并且主要集中在低頻,因此對水下安靜型目標的探測成為水聲領域迫切需要解決的問題[1-2]。傳統的陣列信號處理通常采用加大陣列孔徑的做法來保證低頻的波束寬度和陣增益,但是大孔徑往往會降低工程實用性。如何利用小孔徑陣列實現高增益是長久以來的難題,可以同時從單個陣元的選擇和陣列波束形成方法2個方面入手:高階聲場傳感器和超指向性波束形成方法,兩者結合將為解決這類難題提供一種思路。高階聲場傳感器陣列在消除柵瓣、降低旁瓣、抑制左右舷模糊等方面有很大的潛力,而超指向性波束形成方法則為進一步提高陣列探測性能奠定基礎[3-6]。

高階聲場傳感器的概念最早于1994年由D’Spain[7]提出,他將聲場通過泰勒級數展開,指出了各階項與聲壓以及質點振速之間的關系,并提出利用聲壓二階導數提高陣列的指向性。Cray[8]通過泰勒級數展開的各分量之間的關系,引入了二階聲傳感器,可以提供65°的波束寬度和9.5 dB的最大陣增益,而一般的矢量傳感器只能同時提供105°的波束寬度和6 dB的最大陣增益。可見,高階量的獲取對提高信號處理的能力至關重要。二階矢量傳感器能同步共點測量聲壓、振速以及振速梯度3種聲場分量[9-11],但矢量傳感器的研究最高停留在二階的振速梯度分量上,并且二階分量對非聲學噪聲尤其敏感,穩健性很差,在實際中暫時難以應用。

利用泰勒級數展開并利用有限差分來估算高階量的方法,在理論上是一種近似的求解過程,并不是一個解析解。本文從另一個全新的角度—高階超指向性波束來闡釋高階傳感器以及高階分量的物理含義。Ma等[3, 5]提出了圓環陣特征波束分解與綜合(eigen-beam decomposition and synthesis, EBDS)超指向性模型,能夠提供超指向性的精確閉合解,將最優波束圖表示成了有限階子分量的疊加和形式,并且階數越高,能提供的指向性因子也越大。圓環陣能夠提供圓周360°較為均勻的波束,且沒有左右舷模糊的影響,通過EBDS模型實現高階波束分量的空間指向性,因此可選取圓環陣作為高階傳感器陣元,各階波束分別代表了高階傳感器的各階分量。相比于矢量傳感器目前能提供的二階分量,圓環陣能夠提供更高的有效階數,并且觀測方向靈活可調,而將EBDS各階分量合成后的波束將具有更窄的波束和更高的指向性因子,這相當于是高階傳感器多通道的聯合處理。

在高階傳感器陣列的波束形成研究方面,以往的研究主要針對常規波束形成(conventional beamforming,CBF)分析陣列性能,而對超指向性方法研究較少,仍然停留在傳統的最小方差無失真響應[12](minimum variance distortionless response,MVDR)波束形成。Wang等[13]提出的Gram-Schmidt(GS)模態波束分解與綜合超指向性波束形成適用于任意聲壓陣,但無法處理矢量傳感器或高階聲場傳感器陣列,這是因為GS方法是基于噪聲互譜矩陣全部為實數的假設下求解的,而高階傳感器的噪聲空間相關性往往存在虛部,是一個復數,因此GS方法不再適用。

本文選取圓環陣作為高階傳感器陣元,從特征波束分解與綜合的角度闡釋其高階分量的物理意義;通過將各階分量相加進行的聯合處理,得到高階傳感器的空間接收響應,求解了空間均勻噪聲場中高階傳感器的噪聲空間相關性;最后利用改進的GS(improved gram-schmidt,IGS)模型得到了高階傳感器陣列的超指向性波束、指向性因子以及穩健性水平等,而IGS模型成為真正的適用于帶任意指向性的陣元組成的任意陣型的超指向性方法,填補了具有空間指向性的陣元組成的陣列沒有穩健超指向性方法的空白。

1 高階聲場傳感器

1.1 高階聲場傳感器的構成

傳統的高階聲場傳感器一般指矢量傳感器,其階數最高到2階。本文采用多個無指向性聲壓傳感器組成一個陣列來體現高階傳感器。為方便能夠靈活在周向360°范圍內都能有相對均勻的空間接收響應,可以設計一個圓環形的高階傳感器,見圖1。相對于矢量水聽器中的振速分量的觀測方向固定的缺點,采用圓環陣作為高階傳感器陣元具有觀測方向靈活可控、穩健性可調的優點,并且各階分量在觀測方向上的到達信號都具有無失真的空間接收響應。

圖1 圓環形高階聲場傳感器坐標示意Fig.1 Coordinates of high-order circular array

1.2 高階聲場傳感器的指向性和物理意義

圓環陣EBDS超指向性模型[3]能夠提供超指向性的精確閉合解,將最優波束圖表示成了有限階子分量的疊加和形式,在較寬的頻帶內具有恒定束寬的特性。基于圓環陣的以上特點,在本文中定義:圓環陣為高階傳感器陣元,EBDS子分量波束分別為高階傳感器的各階分量。對于高階聲場傳感器而言,所采用的階數越高,能提供的指向性因子越大,但同時對誤差也越敏感,穩健性越差;如果階數較低,其指向性因子也會降低,因此一個穩健的超指向性波束需要在階數和指向性之間進行權衡。

下面用EBDS超指向性模型具體分析高階聲場傳感器。假設圓環陣均勻分布有M個陣元,如圖2所示,θ和φ分別表示入射平面波的俯仰角和水平方位角,EBDS最優加權向量可表示關于特征值和特征向量的疊加和的形式[3],有:

(1)

vm=M-1/2[1,exp(imβ),…,exp(i(M-1)mβ)]T

(2)

(3)

式中:β=2π/M是相鄰陣元之間的夾角;ρs=sinc(k·Δds);波數k=2π/λ;λ為波長。Δds是第m個陣元和第m′個陣元之間的距離,s=|m-m′|。

圓環陣的最優波束響應為[3]:

(4)

(5)

式中:M為偶數,而Bm在本文中是高階傳感器的第m階分量。假設高階傳感器(圓環陣)由12個無指向性聲壓傳感器組成,半徑為1.5 m,信號頻率為50 Hz。當觀測方向為(θ0,φ0)=(90°,180°)時,圖2給出了第1~4階分量(理論上最高階為第6階)的空間指向性(或空間接收響應)圖,可見第m階分量與第m階的柱面貝塞爾函數具有類似的幾何形狀,而與矢量傳感器cosmθ是不一樣的。

最優波束響應可以由所有特征波束來合成,即:

(6)

這可以看作是高階傳感器的空間響應是各階分量空間響應的疊加和表示。在實際應用中,由于高階波束在低頻時不穩健,因此可以使用降秩技術,舍棄不穩健的高階項而保留低階項來保持穩健性,第0~m階的合成波束可表示為:

(7)

圖2 高階傳感器第1~4階分量的空間指向性(即第1~4階EBDS波束),(θ0,φ0)=(90°,180°)Fig.2 The 1st-4th components of the high-order sensor,(θ0,φ0)=(90°,180°)

圖3 高階傳感器第1~4階合成分量的空間指向性(即第1~4階EBDS合成波束),(θ0,φ0)=(90°,180°)Fig.3 Synthesized parts of the 1st-4th components of the high-order sensor

1.3 高階聲場傳感器的空間噪聲相關性

高階聲場傳感器具有空間指向性,其組陣后第l個陣元的空間接收響應可表示為:

(8)

三維各向同性空間均勻噪聲場可認為是從空間所有方向到達的具有均勻譜密度的獨立平面波的疊加,根據此原理,第l個和第l′個高階傳感器間的噪聲相關性可表示為:

(9)

pl(θ,φ)pl′(θ,φ)*sinθdθdφ

(10)

式(10)較難求解得到解析表達式,但是可以通過數值積分得到一個近似值。

2 IGS超指向性波束形成方法

高階聲場傳感器的空間接收響應具有明顯的指向性,這導致其組陣后的噪聲互譜矩陣的元素(陣元間的噪聲空間相關性)在某些觀測方向上存在虛部。Wang等[13]提出的超指向性方法值適用于噪聲互譜矩陣中的元素全部為實數的陣列,并不適用于矩陣元素存在虛部的情況。針對這一問題,需要對GS超指向性方法的正交變換矩陣重新推導以改進波束形成的性能,以便對有指向性的陣元組成陣列的穩健超指向性波束形成奠定基礎。

對于由Q個陣元組成的具有任意陣型的聲壓陣,基于GS模態分解與綜合的超指向性方法中的噪聲互譜矩陣ρ的逆矩陣可表示為[13]:

ρ-1=CTD-2C

(11)

式中:D為一個Q×Q對角矩陣;D-2用來對輸入噪聲進行歸一化;C是一個正交變換矩陣,為:

C=[c0c1…cq…cQ-1]T

(12)

式中:cq=[cq0cq1…cqq0 … 0]T,當噪聲互譜矩陣ρ中的元素全部是實數時,cqi為:

(13)

(14)

式(14)與式(13)在表達上略有區別,但式(14)也適用于噪聲互譜矩陣中元素全部為實數的情形。IGS方法最優加權向量為:

wopt=αCTD-2C·P(θ0,φ0)

(15)

IGS超指向性波束響應可進一步寫為:

α[C·P(θ0,φ0)]HD-2[C·P(θ,φ)]=

(16)

式中:Bq為第q(q=0,1,…,Q-1)階IGS模態波束,即:

(17)

DF=P(θ0,φ0)Hρ-1P(θ0,φ0)=

[C·P(θ0,φ0)]HD-2[C·P(θ0,φ0)]=

(18)

式中:

Dq=|Eq(θ0,φ0)|2/ηq

(19)

是第q階IGS模態波束的指向性因子。由于式(17)~(19)中的分母均有ηq,并且與誤差敏感度函數有關,是誤差敏感度函數中的分母[13],因此其值的大小對波束形成穩健性有直接的影響,ηq值越大代表第q階波束越穩健,可以將ηq作為對穩健性表征的因子。指向性指數(directivity index, DI)是指向性因子(directivity factors, DF)的分貝表示,有

DI=10lg(DF)

(20)

3 高階聲場傳感器陣列波束形成具體實施步驟

本文中單個高階傳感器是基于EBDS超指向性模型得到的,而高階傳感器所組成陣列的超指向性波束形成是基于IGS模型進行計算的。下面以高階傳感器組成直線陣為例,分析超指向性波束形成具體實施步驟。

假設直線陣由均勻分布在z軸的L個高階傳感器組成,陣元間距為d,如圖4所示;每個高階傳感器半徑為r,由M個無指向性聲壓傳感器組成。高階傳感器陣列波束形成主要分為以下幾步:

圖4 高階聲場傳感器直線陣坐標示意Fig.4 High-order sensor line array

1)得到單個高階傳感器的接收空間響應。

現有一遠場窄帶平面波信號從角度(θ,φ)入射到該直線陣,其中θ為垂直俯仰角,φ為水平方位角,那么該信號傳播方向的單位向量可定義為:

u=-[sinθcosφsinθsinφcosθ]T

(21)

如果高階傳感器陣列在y-z坐標平面上,第1個高階傳感器(圓環陣)的中心位于原點,那么其第m個陣元的位置坐標為:

Lm=[0rcosθmrsinθm]T

(22)

式中:θm=2π(m-1)/M。第m個陣元的接收響應可表示為:

pm(θ,φ)=exp(-ikuTLm)=

exp[ikr(sinθsinφcosθm-cosθsinθm)]

(23)

那么,單個高階傳感器的陣列流形向量為:

Pc(θ,φ)=[p0(θ,φ)p1(θ,φ) …pM-1(θ,φ)]T

(24)

令w為加權向量,那么其波束響應為:

Bc(θ,φ,θ0,φ0)=wc(θ0,φ0)HPc(θ,φ)

(25)

此波束響應具體可由將式(4)和式(7)得到EBDS第m階合成波束,即得到高階傳感器第m階合成分量,從而得到其空間接收響應。

2)確定直線陣陣列流形向量。

假設一個具有L個陣元的直線陣位于z軸上,如圖4所示,相鄰陣元間距為d,如果第1個陣元位于原點,那么第l個陣元的位置坐標:

Ll=[0 0 (l-1)d]T

(26)

對于由無指向性陣元組成的聲壓直線陣,第l個陣元的接收響應可表示:

pl(θ,φ)=exp(-ikuTLl)=exp[ik(l-1)dcosθ]

(27)

當直線陣的陣元為高階聲場傳感器時,每個陣元都有指向性,當觀測方向為(θ0,φ0)時,第l個陣元(圓環陣)的接收響應變:

(28)

這相當于波束圖乘積定理。此時,高階傳感器直線陣的陣列流形向量為:

(29)

3)求高階傳感器陣列噪聲互譜矩陣。

高階傳感器直線陣歸一化的噪聲互譜矩陣為

(30)

式中:矩陣ρL中第l行、第l′列的元素代表了第l個和第l′個陣元間的空間相關性,可由式(10)計算。

4)高階傳感器直線陣波束形成。

得到陣列噪聲互譜矩陣后,依據IGS超指向性方法進行波束形成,通過式(17)、(19)等得到陣列超指向性波束圖、指向性指數(陣增益)、穩健性因子等性能指標。

4 仿真分析

當圓環陣的陣元數較少時,其無法在周向360°范圍內提供均勻的波束;陣元數過多會增加經濟成本,因此應選擇適當的陣元數目以使利益最大化。仿真中選擇由M=12個無指向性聲壓傳感器組成一個高階聲場傳感器。為了使陣列在可探測低頻目標的同時保持較小的孔徑,選取高階傳感器的半徑為r=1.5 m,直線陣由L=7個陣元組成,如圖4所示,陣元間距為d=10 m,水下聲速設為1 500 m/s,在d/λ=0.1時,理論上可探測的目標信號頻率低至15 Hz。仿真中選取直線陣端射和舷側2個典型方向為例進行分析。

4.1 高階聲場傳感器的噪聲空間相關性

噪聲的空間相關性是組成噪聲互譜矩陣的元素,當觀測方向沿z軸的端射方向(即θ0=0°)時,圖5給出了直線陣第1個和第2個高階傳感器陣元間的第1~4階合成分量隨d/λ變化的噪聲相關性曲線。從圖5(a)中可以看到,階數越高,噪聲相關性的實部在d/λ→0時越小;由于在低頻時高階波束對計算誤差較為敏感的原因,第4階在d/λ→0時出現了偏差。噪聲相關性的虛部(見圖5(b))在d/λ→0時接近于0,同時階數越高,曲線初始的起伏振蕩的幅度越小,但隨著頻率的升高,收斂速度也越慢。

當觀測方向為舷側方向(即(θ0,φ0)=(90°,90°))時,圖6給出了直線陣第1個和第2個高階傳感器陣元間的第1~4階合成分量隨d/λ變化的噪聲相關性曲線。舷側方向的噪聲相關性實部的曲線(見圖6(a))變化規律與端射方向(見圖5(a))明顯不同,各階曲線的振蕩起伏更加平緩,相關性實部第1次到達0點的d/λ值遠高于0.5,這提醒我們若使用常規波束形成方法,陣元間距的設置應大于半波長,以使陣元間的噪聲相關性盡量趨于0。

圖5 端射方向高階聲場傳感器的噪聲空間相關性Fig.5 Spatial correlation of the noise in the endfire look direction

圖6 舷側方向高階聲場傳感器的噪聲空間相關性Fig.6 Spatial correlation of the noise in the broadside look direction

4.2 端射方向的陣列性能

高階聲場傳感器由M=12個無指向性聲壓傳感器組成,理論上其最高階分量可取到第6階,但由于波束在不同頻率對誤差的敏感程度不同,可選取適當的階數調整穩健性,這樣既能在低頻時提供穩健的恒定束寬波束,也能在高頻是提供較高的指向性指數,根據圖5和圖6的噪聲相關性曲線,本文中高階傳感器EBDS模型合成分量中的最高階數選取到第4階。

直線陣采用IGS超指向性模型,7元陣的IGS最高階數為第6階。根據式(16)、(17),圖7給出了傳統聲壓直線陣和高階傳感器直線陣在端射方向隨d/λ值變化的IGS方法的波束圖,為節省篇幅,仿真中的波束響應直接采用第0~6階的合成波束,不再分析單階波束響應。從圖7(a)中可以看到,聲壓陣在d/λ約為0.5處會出現柵瓣;而高階傳感器陣(見圖7(b))在較寬的d/λ值范圍內沒有出現柵瓣,表現出較強的抑制柵瓣的能力,同時高階傳感器陣的波束主瓣更窄、旁瓣更低,說明高階傳感器陣能有效抑制來自非觀測方向的干擾和噪聲。

圖7 觀測方向為端射方向時隨d/λ變化的IGS波束圖Fig.7 Beampatterns vs d/λ in the endfire look direction

根據式(19),圖8給出了隨d/λ值變化的IGS第0~6階波束的指向性因子。從圖8(a)中可以看出,在d/λ< 0.5時,聲壓陣階數越高,指向性因子越大;高階傳感器陣則無此規律,反而第0階在較大的頻率范圍內的指向性因子保持最大,而且各階的變化沒有明顯的規律;除第0階外,第1~6階指向性因子隨著頻率的升高逐漸接近。

圖8 端射方向隨d/λ變化的IGS各階指向性因子Fig.8 Different-order directivity factors vs d/λ in the endfire look direction

圖9中IGS的指向性指數是圖8中IGS第0~6階波束的指向性因子相加后的分貝表示,并與CBF進行了對比。作為超指向性方法,2種陣列的IGS的指向性指數在d/λ<0.5時遠遠大于CBF方法。當d/λ→0時,高階傳感器陣的IGS指向性指數比聲壓陣高約3 dB;當d/λ>0.5時,高階傳感器陣在較寬的頻帶內比聲壓陣高約10 dB。聲壓陣在d/λ為0.5時,CBF和IGS方法的指向性指數的曲線基本已經重合;在圖9所示的頻帶內,高階傳感器陣IGS方法的指向性指數始終大于CBF。可見,對于分析高階傳感器陣列的性能不能照搬聲壓陣的規律。

由于ηm與各階IGS模態波束的穩健性有關,可以表征IGS方法的穩健性水平,其值越大,代表第m階波束的穩健性越好。圖10給出了2種陣列隨d/λ值變化的IGS第0~6階波束的穩健性水平。對比圖10(a)和(b),高階傳感器陣的各階穩健性水平均比聲壓陣要低,對誤差更敏感,因此穩健的超指向性波束形成對高階傳感器陣列尤為重要。與EBDS模型類似,在有位置誤差、相位誤差等存在的情況下,可以舍棄不穩健的高階波束,只選取較為穩健的低階波束合成所需的穩健超指向性波束。

圖9 端射方向隨d/λ變化的CBF和IGS方法的指向性指數Fig.9 Directivity index of CBF and IGS vs d/λ in the endfire look direction

圖10 端射方向隨d/λ變化的IGS穩健性水平Fig.10 Robustness levels vs d/λ in the endfire look direction

4.3 舷側方向的陣列性能

根據式(16)和(17),圖11給出了傳統聲壓直線陣和高階傳感器直線陣在舷側方向(θ0,φ0)=(90°,90°)隨d/λ值變化的IGS方法的第0~6階波束的合成圖。從圖11(a)中可以看到,聲壓陣在舷側方向會出現左右舷模糊的現象,并且在d/λ約為1處會出現柵瓣;而高階傳感器陣(見圖11(b))不僅成功消除了左右舷模糊,并且在較寬的頻率范圍內沒有柵瓣。

圖11 舷側方向隨d/λ變化的IGS波束圖Fig.11 Beampatterns vs d/λ in the broadside look direction

根據式(19),圖12給出了舷側方向隨d/λ值變化的IGS第0~6階波束的指向性因子。與圖8(a)對比可以發現,聲壓陣在d/λ<0.5時并沒有出現階數越高,指向性因子越大的現象,除了第0階外,其他階數的指向性因子在整個頻帶內都很接近;高階傳感器陣偶數階(第0,2,4,6階)比奇數階(第1,3,5階)的指向性因子要大,奇數階的指向性因子在低頻時比較接近。除第0階外,第1~6階指向性因子的數值隨著頻率的升高逐漸趨于一致。

圖12 舷側方向隨d/λ變化的IGS各階指向性因子Fig.12 Different-order directivity factors vs d/λ in the broadside look direction

圖13中IGS的指向性指數是圖12中IGS第0~6階波束的指向性因子相加后的分貝表示。通過圖9與圖13對比可以發現,直線陣在端射方向的超指向性比舷側方向更加明顯,這符合以往的理論常識。2種陣列在d/λ> 0.5后,IGS方法的性能都迅速退化為CBF,從而失去了超指向性的效果。總體而言,高階傳感器陣在舷側方向的指向性指數仍然高于聲壓陣約8 dB。

圖14給出了舷側方向2種陣列隨d/λ值變化的IGS第0~6階波束的穩健性水平。與圖10(a)相比,聲壓陣的穩健性水平變化較小(見圖14(a)); 與圖10(b)相比,高階傳感器陣的各階穩健性水平(見圖14(b))有所提高,說明高階傳感器陣在舷側方向的超指向性波束更穩健。

圖13 端射方向隨d/λ變化的CBF和IGS方法的指向性指數Fig.13 Directivity index of CBF and IGS vs d/λ in the endfire look direction

圖14 舷側方向隨d/λ變化的IGS穩健性水平Fig.14 Robustness levels vs d/λ in the broadside look direction

5 結論

1)本文通過EBDS高階超指向性模型給出了以圓環陣作為高階聲場傳感器的物理意義,然后研究了高階傳感器陣元的噪聲空間相關性及其陣列IGS超指向性波束形成方法,并舉例給出了波束形成具體實施步驟。

2)仿真中采用了12元均勻圓環陣作為高階傳感器陣元,并選取EBDS第0~4階的合成分量給出了高階傳感器陣元的空間接收響應以及噪聲相關性,然后將高階傳感器陣元組成直線陣采用IGS超指向性模型進行性能分析。仿真結果表明,具有空間指向的高階傳感器陣列具有較強的消除柵瓣和抑制左右舷模糊的能力,7個4階傳感器組成的直線陣比傳統聲壓陣的指向性指數在較寬的頻帶內高約8 dB以上。

3)高階聲場傳感器在水聲應用中可組成端射陣和拖曳線列陣等,這對于探測水下低頻安靜型潛艇等目標具有重要的實際意義。

猜你喜歡
方向
2023年組稿方向
計算機應用(2023年1期)2023-02-03 03:09:28
方向
青年運動的方向(節選)
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
如何確定位置與方向
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
大自然中的方向
主站蜘蛛池模板: 国产精品无码AⅤ在线观看播放| 日韩精品免费一线在线观看 | 88av在线看| 国产精品 欧美激情 在线播放| 在线不卡免费视频| 青草精品视频| 又粗又硬又大又爽免费视频播放| 亚洲精品你懂的| 日韩欧美视频第一区在线观看| 一区二区影院| 99久久免费精品特色大片| 色吊丝av中文字幕| 人妻丰满熟妇av五码区| 四虎在线观看视频高清无码 | av天堂最新版在线| 五月天久久婷婷| 美女国内精品自产拍在线播放| 日本精品视频| 青青操国产视频| a级毛片免费网站| 欧美一级99在线观看国产| 久久国产热| 欧美怡红院视频一区二区三区| 久青草免费在线视频| 欧美五月婷婷| 看你懂的巨臀中文字幕一区二区| 91人妻日韩人妻无码专区精品| 日本不卡在线视频| 91精品专区国产盗摄| 国产va在线观看| 国产精品视频999| 日韩毛片免费观看| 男人天堂亚洲天堂| 亚洲娇小与黑人巨大交| 精品国产91爱| 99re视频在线| 热99精品视频| 99久久精品无码专区免费| 精品乱码久久久久久久| 日韩视频免费| 国产成人亚洲毛片| 成人日韩视频| 国内自拍久第一页| 久久亚洲黄色视频| 免费a在线观看播放| 一级做a爰片久久免费| 国产99热| 午夜少妇精品视频小电影| 国产精品13页| 久热这里只有精品6| 国产情侣一区| 国产精品欧美在线观看| 国语少妇高潮| 一级高清毛片免费a级高清毛片| 欧美丝袜高跟鞋一区二区| 波多野结衣二区| 一级香蕉人体视频| 九色在线观看视频| 中文字幕欧美日韩| 国产精品一区二区不卡的视频| 亚洲色图另类| 青青国产在线| 亚洲有无码中文网| 青草视频久久| 18禁高潮出水呻吟娇喘蜜芽| 四虎影院国产| 色偷偷综合网| 亚洲天堂区| 毛片免费高清免费| 老司国产精品视频91| 国产一级毛片网站| 亚洲色婷婷一区二区| 国产精品入口麻豆| 热re99久久精品国99热| 国产亚洲精品yxsp| 欧美日韩激情| 99视频精品全国免费品| 免费毛片全部不收费的| 国产精品人成在线播放| 色婷婷狠狠干| 毛片视频网址| 国产精品播放|