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

近場寬帶LFM信號被動測向和測距方法

2016-08-15 08:09:52林秋華楊秀庭康春玉
系統工程與電子技術 2016年8期
關鍵詞:信號方法

李 軍, 林秋華, 楊秀庭, 康春玉

(1. 大連理工大學電子信息與電氣工程學部, 遼寧 大連 116024;2. 海軍大連艦艇學院信息作戰系, 遼寧 大連 116018)

?

近場寬帶LFM信號被動測向和測距方法

李軍1,2, 林秋華1, 楊秀庭2, 康春玉2

(1. 大連理工大學電子信息與電氣工程學部, 遼寧 大連116024;2. 海軍大連艦艇學院信息作戰系, 遼寧 大連 116018)

提出了一種采用分數階傅里葉變換的聚焦波束形成被動定位方法,實現了水聲近場寬帶線性調頻(linear frequency modulated, LFM)信號的被動測向和測距。建立了基于球面波模型的近場寬帶LFM信號接收數據模型,應用分數階傅里葉變換(fractional Fourier transform, FRFT)將LFM信號的時變陣列流形矩陣變換為固定陣列流形矩陣,結合近場聲源的聚集波束形成技術,利用多重信號分類算法實現了對多個寬帶LFM信號的方位與距離聯合估計。數值仿真驗證了該方法對水聲目標方位和距離估計的有效性,并仿真分析信噪比、聲源距離、聲源個數等對該算法性能的影響。

線性調頻信號; 被動測距; 分數階傅里葉變換; 多重信號分類法

0 引 言

隱身技術的發展使現代潛艇變得越來越安靜,大大降低了被動聲吶的作用距離;傳統主動聲吶由于收發裝置合置,在探測目標的同時也容易暴露自己。近年來,采用收發裝置分置的多基地聲吶系統由于具有隱蔽性能好、抗干擾能力強、定位精度高等優點正成為人們研究熱點[1-3]。為了使聲吶基陣具備良好的隱蔽性、更佳的環境適用性和更大的基陣孔徑,發射基地一般遠離目標發射寬帶線性調頻(linear frequency modulated,LFM)信號,接收基地接近目標應用被動拖曳線列陣聲吶接收回波信號來實現目標探測,因此針對近場寬帶LFM信號的被動測向和測距問題研究具有重要意義。

當目標位于近場環境時,影響各陣元接收信號時延的主要因素不但包括目標的方位信息,而且還包括目標的距離信息,此時應采用球面波模型進行建模[4]。文獻[5-6]結合近場球面波傳播模型,利用常規聚焦波束形成方法(conventional focused beamformer, C-FB),實現了水下目標方位和距離的估計。文獻[7-8]提出了一種最小方差無失真響應聚焦波束形成方法(minimum variance distortionless response focused beamformer, MVDR-FB),進一步提高了定位的精度。文獻[9-10]改進了二階多重信號分類方法(multiple signal classification,MUSIC),實現了近場和遠場目標的聯合定位。針對窄帶信號,這些方法能夠取得較好的定位效果;但是針對寬帶信號,這些方法多數采用計算量比較大的頻域分子帶技術[11-13]完成水下目標的定位,并沒有較好地利用寬帶LFM信號自身的特點進行定位。

根據寬帶LFM信號在分數階傅里葉域具有極高的聚焦性,文獻[14-16]提出了基于分數階傅里葉變換(fractional Fourier transform, FRFT)的LFM信號測向方法。文獻[17]修正了該方法的中心頻率估計問題,進一步提高了測向精度。文獻[18]提出了基于FRFT的水聲信道參數估計方法,為水聲多途環境下目標定位提供了理論支持。然而,這些方法主要是基于遠場環境平面波模型進行目標方位的估計,對近場環境中的目標方位估計則存在較大誤差,而且不能對目標距離實現有效估計。

針對上述問題,本文提出了一種采用FRFT的聚焦波束形成被動定位方法,簡稱FRFT-MUSIC-FB。針對近場寬帶LFM信號,構建了基于球面波模型的信號接收數據模型,應用FRFT將傳統的陣列數據模型轉化成分數階傅里葉域陣列數據模型,同時將LFM信號的時變陣列流形矩陣變換為固定陣列流形矩陣。結合近場聲源的聚集波束形成技術,通過對分數階傅里葉域的相關矩陣進行特征值分解,估計信號子空間和噪聲子空間,利用MUSIC算法實現了對多個LFM信號的方位與距離聯合估計。通過仿真實驗,對該方法的有效性進行了驗證,重點研究了信噪比、聲源距離、聲源個數對該算法性能的影響。

1 近場寬帶LFM信號接收模型

均勻線列陣近場信號接收模型[4,10]如圖1所示。

設有M個陣元,間距為d,陣列的第一個陣元為參考陣元,并標記為1點。設有Q個LFM信源,第i個LFM信號模型[14]為

(1)

式中,fi為第i個信號的初始頻率;μi為第i個信號的調頻率。

設θi、ri分別為第i個信號到參考陣元的入射角和距離,c為聲速,rim為第i個信號與第m個陣元的距離,第i個信號到達第m個陣元時相對于參考陣元的時延τim為[5-6]

(2)

則第m個陣元的接收信號時域表示形式為

(3)

式中,si(t)表示接收的第i個信號;nm(t)是附加的高斯白噪聲。將所有陣元輸出寫成向量形式,并將式(1)、式(2)代入式(3),得到觀測信號為

(4)

式中,N(t)=[n1(t),n2(t),…,nM(t)]T是高斯白噪聲,空間協方差為σ2I;ai(θi,ri,t)為第i個信號在M線陣上的近場導向矢量矩陣,可表示為

(5)

式中,時延τim可由式(2)求得。由式(5)可以看出ai(θi,ri,t)是由信源初始頻率fi、調頻率μi、信源相對于參考陣元的方位θi和距離ri等未知參數所組成的函數構成。

2 分數階傅里葉域近場寬帶LFM信號定位方法

根據球面波模型建立均勻線列陣近場寬帶LFM信號接收數據模型,利用分數階傅里葉域數據模型代替傳統的陣列數據模型,將近場寬帶LFM信號的時變陣列流形矩陣變換為固定陣列流形矩陣,結合聚集波束形成技術,利用MUSIC算法實現信號方位和距離的聯合估計。

2.1分數階傅里葉域近場LFM信號接收模型

設si1(t)為參考陣元(1號陣元)上接收的第i個LFM信號,則有

(6)

以采樣率fs采樣得到N點si1(n),有

(7)

式中,N為快拍數。根據Ozaktas提出的快速離散FRFT算法[19],離散信號si1(n)關于α角的N點離散FRFT(digitalFRFT,DFRFT)信號Si1(α,v)可表示為

(8)

通過對Si1(α,v)在(α,v)平面進行峰值搜索可得到與LFM信號對應的峰值坐標(αi1,vi1),與信源初始頻率fi、調頻率μi的對應關系為

(9)

設sim(t)為第m個陣元上接收到的第i個LFM信號,其與原始信號si(t)、參考陣元接收信號si1(t)的關系為

(10)

設sim(t)的分數階傅里葉變換為Sim(α,v),其出現峰值的坐標(αim,vim)滿足[14]

(11)

(12)

(13)

式中

(14)

由Ai(θi,ri,τim)可構成第i個信號在分數階傅里葉域的方向向量為

(15)

Ai(θi,ri)取值僅與時延τim有關,即只與第i個信號的入射角θi和距離ri有關,是時不變的。由Ai(θi,ri)可組成分數階傅里葉域的陣列流形矩陣A(θ,r),表示為

(16)

Q個LFM信號的疊加會在分數階傅里葉域出現Q個峰值,分別選擇分數階傅里葉域上Q個峰值點上的數據作為該陣元的觀測數據,則第m個陣元上的分數階傅里葉域輸出為

(17)

所有陣元上的分數階傅里葉域輸出為

(18)

對式(3)兩端采樣并進行離散分數階傅里葉變換,并利用式(16)~式(18)可以得到近場寬帶LFM信號在分數階傅里葉域的接收數據模型為

(19)

(20)

2.2基于聚焦波束形成的LFM信號定位方法

設RXX為輸出數據在分數階傅里葉域的相關矩陣,根據MUSIC算法的基本原理,對RXX進行特征分解

(21)

式中,US是由大特征值對應的特征矢量張成的信號子空間;UN是由小特征值對應的特征矢量張成的噪聲子空間。利用分數階傅里葉域的相關矩陣代替傳統的陣列相關矩陣,根據MUSIC算法可得到第i個信號的FRFT-MUSIC-FB空間譜為

(22)

對式(22)進行關于(θ,r)的二維搜索,Pi(θ,r)的最大值對應的坐標(θi,ri)即為第i個信號的入射角θi和距離ri。

FRFT-MUSIC-FB算法具體計算步驟如下:

步驟 1對1號陣元(參考陣元)接收數據進行DFRFT得到X1(α,v),對(α,v)平面的峰值作二維搜索,分別得到Q個LFM信號對應的(αi1,vi1)。

步驟 2根據式(11),各個LFM信號對應的αim相等,由式(8)分別計算各個陣元上接收數據關于αi1角的DFRFT得到Xm(αim,v),作一維峰值搜索得到第i個信號對應的vim。

步驟 3由各陣元數據Xm(αim,vim)構造第i個LFM信號在分數階傅里葉域對應的接收數據矩陣為

步驟 4計算分數階相關矩陣RXX=X(αi,vi)X(αi,vi)H,對RXX進行特征分解,估計噪聲子空間UN。

步驟 5明確方位θ和距離r的二維搜索范圍及步長,由式(2)、式(14)、式(15)計算每一組(θ,r)對應第i個LFM信號的方向向量Ai(θi,ri)。

步驟 6由式(22)得到第i個LFM信號對應的FRFT-MUSIC-FB空間譜Pi(θ,r),極大值對應的方位和距離(θmax,rmax)就是第i個LFM信號的方位和距離估計。

步驟 7存在多個LFM信號,重復步驟2~步驟6,分別得到各個信號的方位與距離估計。

3 仿真數據驗證與分析

為了驗證本文方法的有效性,下面將進行3種仿真實驗,分別測試本文方法在單一LFM信號、雙LFM信號條件下,信噪比、目標方位及距離等因素對定位性能的影響,并與頻域分子帶C-FB、MVDR-FB等方法定位性能進行比較分析。

采用蒙特卡羅仿真分析算法性能時,一般應用均方根誤差(rootmean-squareerror,RMSE)來評價系統方位和距離估計的性能

(23)

(24)

(25)

仿真試驗中,陣元數M=32,陣元間距d=1m,聲速c=1 500m/s,距離掃描范圍1~1 000m,步長1m,方位掃描范圍(-60°,60°),步長1°,信號的快拍數K=4 800,系統采樣頻率fs=48kHz,噪聲為加性高斯白噪聲。

仿真分析 1單一LFM信號定位分析

仿真產生LFM信號初始頻率為2kHz,調頻斜率為10kHz/s,脈沖寬度為0.1s,頻率變化范圍為2 ~3kHz。當目標的方位為10°,距離分別為200m、500m、800m時,應用本文方法進行目標定位得到的FRFT-MUSIC-FB譜如圖2所示,可以看出,當目標距離小于500m時,能夠對目標方位和距離里進行較精確的估計。當目標距離逐漸增大時,FRFT-MUSIC-FB譜由細變粗,聚焦性能下降,此時方位估計依然比較準確,但距離估計誤差逐漸增大。

圖2 單LFM信號被動定位的FRFT-MUSIC-FB譜圖

仿真分析 2信噪比、目標方位及距離對定位性能的影響

分別研究信噪比對頻域分子帶C-FB、MVDR-FB和本文FRFT-MUSIC-FB方法定位性能的影響。仿真中只考慮存在一個LFM信號的情況,信號形式同仿真1,蒙特卡羅仿真次數為100次。當目標位于(20°,500 m)時,分別應用3種方法進行定位,圖3畫出了3種方法在不同信噪比下目標方位估計的RMSE,可以看出當信噪比大于-12 dB時,3種方法對目標方位估計的RMSE開始收斂,隨著信噪比的增大,頻域分子帶MVDR-FB和本文FRFT-MUSIC-FB方法的測向誤差將趨近于0,性能要優于頻域分子帶C-FB方法。圖4畫出了3種方法在不同信噪比下目標距離估計的R-RMSE,可以看出當信噪比大于-10 dB時,3種方法對目標距離估計的R-RMSE收斂于不同的較小值,本文方法的抗噪性和估計精度要優于其他兩種方法。

圖3 不同信噪比下3種方法對目標方位估計的RMSE

圖4 不同信噪比下3種方法對目標距離估計的R-RMSE

為了更好地分析本文FRFT-MUSIC-FB方法在不同目標距離和方位的條件下對目標進行定位的性能。圖5畫出了當信噪比為0 dB,目標距離從100 m變化到500 m,目標方位從-40°變化到40°時,應用本文方法進行目標距離估計的R-RMSE。可以看出,目標靠近正橫方位(0°)時本文方法距離估計比較準確,距離估計的R-RMSE隨著目標方位角的增大而逐漸增大。

為了增加定位精度,本文根據拖曳線列陣聲吶自身的特點選擇距離母船最近的1號陣元作為參考陣元,建立近場球面波數據模型,如圖1所示。當目標i位于(θi,ri)(θi>0)時,目標與2至M號陣元的實際距離ri2,ri3,…,riM必然小于目標位于(-θi,ri)時與對應陣元的實際距離,線列陣所得接收信號也優于目標位于(-θi,ri)時的接收信號。因此,對位于(θi,ri)的目標距離估計精度要優于與其方位對稱的(-θi,ri)點。圖5中的距離誤差數據存在明顯的方位不對稱性也進一步驗證了上述理論。

圖5 不同目標位置對本文方法距離估計的R-RMSE

仿真分析 3海洋多徑干擾對定位性能的影響

分別研究海洋多徑干擾對3種方法定位性能的影響,仿真基本條件同仿真分析2。仿真采用的水聲多徑信道為典型的3徑水聲信道[18,20],每條路徑的衰減系數和相對時延在表1中給出。

表1 水聲信道多路徑參數

當目標位于(20°,500 m),信噪比為0 dB,多徑1的衰減系數從0變化到1時,圖6和圖7分別畫出了3種方法在不同多徑參數下目標方位估計的RMSE和目標距離估計的R-RMSE。

圖6 不同多徑參數下3種方法對目標方位估計的RMSE

圖7 不同多徑參數下3種方法對目標距離估計的R-RMSE

可以看出,當多徑1的衰減系數小于0.7時,多徑干擾對算法性能的影響有限,但是當衰減系數大于0.8時,多徑干擾對算法性能的影響逐步增大,本文方法的抗多徑干擾性能要優于其他兩種方法。

仿真分析 4雙LFM信號定位分析

設存在兩個LFM信號時,基本參數如表2所示。

表2 雙LFM信號基本參數

當信噪比為0 dB時,分別應用頻域分子帶C-FB、MVDR-FB和本文FRFT-MUSIC-FB方法進行目標定位,所得到的空間譜如圖8、圖9和圖10所示。可以看出,頻域分子帶C-FB方法的偽峰數量明顯增多,真假目標的判斷較難,波束圖的聚焦性能也較差。頻域分子帶MVDR-FB方法的各種性能有所提升,FRFT-MUSIC-FB方法的聚焦性能和定位精度要優于其他兩種方法。

圖8 頻域分子帶C-FB譜(目標位置(0°,200 m)、(10°,100 m))

圖9 頻域分子帶MVDR-FB譜(目標位置 (0°,200 m)、(10°,100 m))

圖10 FRFT-MUSIC-FB譜(目標位置(0°,200 m)、(10°,100 m))

4 結束語

針對發射信號采用寬帶LFM信號的多基地聲吶系統目標定位問題,提出了一種基于分數階傅里葉變換的聚焦波束形成被動定位方法,將傳統的陣列數據模型轉化為分數階傅里葉域陣列數據模型,應用MUSIC聚焦波束形成技術,實現了近場寬帶LFM信號的被動測向和測距。通過仿真實驗,對該方法的有效性進行了驗證,該方法對目標方位的估計性能受外界條件的影響較小,對目標距離的估計精度受目標原始方位距離的影響較大,距離估計的誤差隨著目標方位角的增大而增大。與頻域分子帶C-FB、MVDR-FB方法進行比較,該方法具有較高的定位精度、較強的抗噪性和抗多徑干擾特性。

[1] Xu J F, Han S P, Shu X L, et al. Research on systemic accuracy of localization algorithm for the type of T-R bistatic sonar[J].SystemsEngineeringandElectronics, 2015, 37(10):2234-2238. (徐景峰, 韓樹平, 舒象蘭, 等. T-R型雙基地聲吶定位算法及系統精度研究[J].系統工程與電子技術, 2015, 37(10):2234-2238.)

[2] Aughenbaugh J M, La Cour B R, Gelb J M. Interference mitigation for multistatic active sonar[J].IEEEJournalofOceanicEngineering, 2015, 40(3):570-582.

[3] Liang J L, Xu L Z, Li J, et al. On designing the transmission and reception of multistatic continuous active sonar systems[J].IEEETrans.onAerospaceandElectronicSystems,2014,50(1):285-299.

[4] Cirpan H A, Cekli E. Deterministic maximum likelihood approach for localization of near-field sources[J].InternationalJournalofElectronicsandCommunications, 2002, 56(1):1-10.

[5] Hui J, Hu D, Hui J Y, et al. Researches on the measurement of distribution image of radiated noise using focused beamforming[J].ActaAcustica,2007,32(4):356-361.(惠娟,胡丹,惠俊英,等.聚焦波束形成聲圖測量原理研究[J].聲學學報,2007,32(4):356-361.)

[6] Mei J D, Hui J Y, Hui J. Research on simulation of near field passive ranging with underwater acoustic image by focused beam-forming[J].JournalofSystemSimulation, 2008, 20(5):1328-1333. (梅繼丹,惠俊英,惠娟.聚焦波束形成聲圖近場被動定位技術仿真研究[J].系統仿真學報,2008,20(5):1328-1333.)

[7] Shi J, Yang D S, Liu B S, et al. Radiated noise sources near-field location based on MVDR focused beamforming[J].JournalofDalianMaritimeUniversity, 2008, 34(3):55-58. (時潔, 楊德森, 劉伯勝, 等. 基于MVDR聚焦波束形成的輻射噪聲源近場定位方法[J].大連海事大學學報, 2008, 34(3):55-58.)

[8] Shi J, Yang D S. Coherent broadband MVDR focused beamforming based on vector sensor array processing[J].JournalofSystemSimulation, 2010, 22(2):473-477. (時潔, 楊德森. 矢量陣相干寬帶MVDR聚焦波束形成[J].系統仿真學報, 2010, 22(2):473-477.)

[9] Jiang J J, Duan F J, Chen J, et al. Mixed near-field and far-field sources localization using the uniform linear sensor array[J].IEEESensorsJournal, 2013, 13(8):3136-3143.

[10] Liang J L, Liu D. Passive localization of mixed near-field and far-field sources using two-stage MUSIC algorithm[J].IEEETrans.onSignalProcessing, 2010, 58(1):108-120.

[11] Jafari A H, Liu W, Morgan D R. Study of sensor positions for broadband beamforming[J].IEEESignalProcessingLetters, 2013, 20(8):779-782.

[12] Jian Y, Xi X Q, Yu Y. A implementation method to constant beamwidth and adaptive beamforming of broadband signal[C]∥Proc.oftheIEEEInternationalConferenceonSignalProcessing,CommunicationsandComputing,2014:394-397.

[13] Sam K A, Jesper R J, Mads G C. Fast joint DOA and pitch estimation using a broadband MVDR beamformer[C]∥Proc.ofthe21stEuropeanSignalProcessingConference, 2013: 1-5.

[14] Tao R, Zhou Y S. A novel method for the direction of arrival estimation of wideband linear frequency modulated sources based on fractional fourier transform[J].TransactionsofBeijingInstituteofTechnology,2005,25(10):895-899.(陶然,周云松.基于分數階傅里葉變換的寬帶LFM信號波達方向估計新算法[J].北京理工大學學報, 2005, 25(10):895-899.)

[15] Yang X M, Tao R. 2-D DOA estimation of LFM signals based on fractional fourier transform[J].ActaElectronicaSinica, 2008, 36(9):1737-1740. (楊小明, 陶然. 基于分數階傅里葉變換的線性調頻信號二維波達方向估計[J].電子學報, 2008, 36(9):1737-1740.)

[16] Liu S H, Shan T, Zhang Y M D, et al. A fast algorithm for multi-component LFM signal analysis exploiting segmented DPT and SDFrFT[C]∥Proc.oftheIEEEInternationalRadarConference, 2015: 1139-1143.

[17] Wang R, Ma Y. DOA estimation of wideband linear frequency modulated pulse signals based on fractional fourier transform[J].ActaArmamentarii,2014,35(3):421-427.(王瑞,馬艷.基于分數階傅里葉變換的線性調頻脈沖信號波達方向估計[J].兵工學報,2014,35(3):421-427.)

[18] Yin J W, Hui J Y, Cai P, et al. Underwater acoustic channel parameter estimation based on fractional fourier transforms[J].SystemsEngineeringandElectronics, 2007, 29(10):1624-1627.(殷敬偉,惠俊英,蔡平,等.基于分數階Fourier變換的水聲信道參數估計[J].系統工程與電子技術,2007,29(10):1624-1627.)

[19] Ozaktas H M, Arikan O, Kutay M A, et al. Digital computation of the fractional Fourier transform[J].IEEETrans.onSignalProcessing, 1996, 44(9):2141-2150.

[20] Feder M, Catipovic J A. Algorithm for joint channel estimation and data recovery-application to equalization in underwater communications[J].IEEEJournalofOceanicEngineering, 1991, 16(1):42-55.

Passive DOA and range estimation method for near-field broadband LFM signals

LI Jun1,2, LIN Qiu-hua1, YANG Xiu-ting2, KANG Chun-yu2

(1. Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology,Dalian 116024, China; 2. Department of Information Operations, Dalian Naval Academy, Dalian 116018, China)

A passive underwater acoustic localization method using fractional Fourier transform (FRFT) and focused beamformer is proposed. The direction-of-arrival and range information of the near-field broadband linear frequency modulated (LFM) signals can be estimated. The near-field spherical wave model for broadband LFM signal is introduced and the time-varying array manifold matrix corresponding to the LFM signals is changed into a fixed array manifold matrix by using FRFT. Combined with the near-field focused beamforming technology and multiple signal characteristics algorithm, the azimuths and distances of multiple broadband LFM signals will be estimated. The validity and efficiency of the proposed method is verified by simulations and experiments, and the performance of the proposed method is discussed in terms of signal-to-noise ratio, azimuth and distance of the underwater acoustic signal.

linear frequency modulated (LFM) signal; passive ranging; fractional Fourier transform (FRFT); multiple signal classification (MUSIC)

圖1均勻線列陣近場信號接收模型Fig.1The model of signal by uniform linear array in near-field

2015-04-20;

2016-02-16;網絡優先出版日期:2016-03-22。

國家自然科學基金(61271443,61471378)資助課題

TN 911.5

A

10.3969/j.issn.1001-506X.2016.08.05

李軍(1981-),男,講師,博士研究生,主要研究方向為信號與信息處理、水聲信號處理。

E-mail:lijunwk@163.com

林秋華(1969-),女,教授,博士,主要研究方向為信號與信息處理。

E-mail:qhlin@dlut.edu.cn

楊秀庭(1973-),男,副教授,博士,主要研究方向為水聲信號處理、聲吶作戰使用。

E-mail:yangxiuting2003@163.com

康春玉(1975-),男,副教授,博士,主要研究方向為水聲信號處理。

E-mail:dlkangcy@126.com

網絡優先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20160322.1453.008.html

猜你喜歡
信號方法
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
學習方法
孩子停止長個的信號
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 无码人中文字幕| 欧美成人区| 无码国产偷倩在线播放老年人| 国产成人综合日韩精品无码不卡| 国产激情无码一区二区APP| 国产又大又粗又猛又爽的视频| 男女猛烈无遮挡午夜视频| 永久免费精品视频| 色屁屁一区二区三区视频国产| 丁香婷婷久久| 9999在线视频| 69免费在线视频| 国产免费久久精品99re不卡| 亚洲国产日韩一区| 尤物成AV人片在线观看| 国产美女丝袜高潮| 欧美精品成人| 国产一二视频| 国产AV无码专区亚洲A∨毛片| 最新日本中文字幕| 青青草原国产精品啪啪视频| 国产日韩欧美一区二区三区在线| 呦女亚洲一区精品| 亚洲色图综合在线| 精品福利网| 免费毛片在线| 国产亚洲精品91| 一级毛片高清| 东京热高清无码精品| 国产成人久久综合777777麻豆| 91免费国产高清观看| 亚洲成aⅴ人在线观看| 欧美精品一二三区| 四虎国产精品永久一区| 四虎永久免费网站| 免费jizz在线播放| av尤物免费在线观看| 日本亚洲成高清一区二区三区| 国产精品19p| 综合网天天| 日本国产一区在线观看| 国产精品久久自在自线观看| 伊在人亚洲香蕉精品播放| 欧类av怡春院| 欧美a在线| 日韩午夜伦| 丝袜无码一区二区三区| 国产精品毛片一区| 91无码视频在线观看| 国产福利不卡视频| 91精品久久久久久无码人妻| 国内99精品激情视频精品| 欧美一区二区三区不卡免费| 欧洲日本亚洲中文字幕| 国产大片喷水在线在线视频| 免费日韩在线视频| 中文字幕在线永久在线视频2020| a级毛片毛片免费观看久潮| 国产精品天干天干在线观看| 国产成人免费| 国产一区二区三区视频| 波多野结衣无码AV在线| 国产精品va免费视频| a级毛片网| 国产日韩精品一区在线不卡 | 中文字幕2区| 女人爽到高潮免费视频大全| 日韩久草视频| 午夜在线不卡| 中文字幕免费在线视频| 国产亚洲欧美日韩在线一区| 丝袜亚洲综合| 亚洲国产成人久久77| 国产视频一二三区| 亚洲婷婷六月| 亚洲 欧美 日韩综合一区| 色有码无码视频| 911亚洲精品| 国产精品自在在线午夜区app| 免费无码在线观看| 亚洲色成人www在线观看| 亚洲人成网线在线播放va|