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

基于張量分解的互質陣MIMO雷達目標多參數估計方法

2015-07-12 13:55:59樊勁宇顧紅蘇衛民陳金立
電子與信息學報 2015年4期

樊勁宇顧 紅蘇衛民陳金立

①(南京理工大學電子工程與光電技術學院 南京 210094)

②(南京信息工程大學電子與信息工程學院 南京 210044)

基于張量分解的互質陣MIMO雷達目標多參數估計方法

樊勁宇*①顧 紅①蘇衛民①陳金立②

①(南京理工大學電子工程與光電技術學院 南京 210094)

②(南京信息工程大學電子與信息工程學院 南京 210044)

該文提出了一種基于雙基地互質陣列(CPA)多輸入多輸出(MIMO)雷達的多目標波離方向角(DOD)、波達方向角(DOA)和多普勒頻率估計算法。收發陣列各由兩個滿足互質結構的稀疏均勻子線陣組成。時域的快拍序列同樣由兩個互質的稀疏均勻采樣構成。算法利用張量因子分解得到分別包含DOD, DOA和多普勒頻率信息的3個流形矩陣,再從中構造出具有范德蒙德矩陣結構的虛擬流形矩陣。為了提高估計精度,還提出了一種基于特征值分解的誤差抑制算法,并通過旋轉不變子空間算法(ESPRIT)求取各目標的3個待估參數。與傳統算法相比,該算法通過構造均勻虛擬陣列和虛擬快拍提高參數估計性能,且不會產生模糊,避免了譜峰搜索和額外的配對過程。仿真實驗驗證了該算法有效性。

雙基地MIMO雷達;互質數對;張量分解;Swerling-I模型

1 引言

由無線通信中多輸入多輸出(MIMO)技術發展而來的的MIMO雷達近些年來引起了人們的廣泛關注,文獻[1]首次提出這一概念。與相控陣雷達相比,該雷達具有更大的陣列孔徑、更高的陣列自由度和更好的低截獲性能[2?4]。得益于發射端的波形分集特性,雙基地MIMO雷達可以同時估計多個目標的波離方向角(Direction Of Departure, DOD)和波達方向角(Direction Of Arrival, DOA),從而實現交叉定位。文獻[5]將2維多重信號分類(2 Dimension-MUltiple SIgnal Classification, 2D-MUSIC)技術引入雙基地MIMO雷達,能夠同時估計多個目標的DOD和DOA,并提出了一種降維MUSIC算法降低2維譜峰搜索的運算量。文獻[6]通過2次獨立的旋轉不變子空間參數估計(Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT)算法分別估計目標DOD和DOA,算法避免了譜峰搜索,但需要進行額外的參數配對。文獻[7]對其做出了改進,提出一種自動配對的ESPRIT算法。文獻[8]將回波數據構造成測量張量,利用張量因子分解得到目標的各個流形矩陣,再通過傳統子空間算法(如ESPRIT等)求得目標的DOD和DOA。在目標雷達截面積(Radar Cross Section, RCS)滿足Swerling-I模型的條件下還可求得各目標的多普勒頻率,且這些參數均自動配對。上述算法要求陣元間距滿足半波長限制,且脈沖重復頻率(Pulse Recurrence Frequency, PRF)大于2倍多普勒頻率,否則估計的參數會存在模糊。為了在相同陣元個數下得到更大的陣列孔徑,非均勻陣列被引入MIMO雷達信號處理領域。

非均勻線陣的各陣元間距組成一個特定序列,且均為半波長的整數倍。常見的非均勻線陣有:最小冗余陣列[9](Minimum Redundancy Array, MRA)、嵌套陣列[10](Nested Array, NA)和互質陣列[11](Co-Prime Array, CPA)等。文獻[12]將非均勻線陣應用于雙基地MIMO雷達,從匹配濾波輸出數據的自相關矩陣中構造出虛擬MIMO陣列的導向矢量。再通過2維空間平滑算法估計出各目標的DOD和DOA。該算法提高了MIMO雷達的角度估計精度且可以實現估計參數的自動配對。當物理陣元個數較多時形成的虛擬陣元數將大大增加,此時2維空間平滑算法得到的2維虛擬陣列流形矩陣過于龐大不利于后續運算。此外,自相關運算會導致運算過程中丟失時域信息,因而無法估計目標的多普勒頻率。

本文提出了一種基于非均勻線陣的目標參數估計算法。選用了互質陣作為收發陣列,并對文獻[8]中均勻采樣的快拍序列做出改進,僅對其中構成互質序列的若干非均勻快拍進行采樣。從而將非均勻陣列角度估計算法拓展到了非均勻快拍下的多普勒頻率估計中。算法首先將匹配濾波后數據構造成三階張量,并利用張量分解從中求取互質結構下的收發陣列流形矩陣和多普勒信息矩陣。然后通過矩陣的Khatri-Rao乘積運算估計出均勻虛擬收發線陣各自的陣列流形矩陣和均勻虛擬快拍下的多普勒信息矩陣。針對運算過程中誤差放大問題提出了一種利用特征值分解進行誤差抑制的方法。最后利用ESPRIT算法從誤差抑制后的3組導向矢量中逐一估計出各目標的DOD, DOA和多普勒頻率。算法利用非均勻線陣形成的虛擬陣列獲得更高的角度估計性能,通過構造數量遠高于實際采樣快拍數的均勻虛擬快拍提高了多普勒頻率估計精度,并降低了接收端信號處理的存儲成本。該算法能夠實現參數間自動配對,且適用于MRA等其他非均勻陣列,文中僅以互質陣為例是為了便于非均勻陣列和非均勻快拍序列的構造。

2 系統模型

雙基地MIMO雷達的收發陣列分別為N和M個陣元的非均勻線陣。M路正交窄帶發射信號S=[s1,s2,…,sm]T∈?M×P滿足SSH=IM,其中(·)T表示轉置,(·)H表示共軛轉置,P為每個脈沖的采樣點數,IM為M維單位矩陣。假設有K個不相關的遠場目標,收發陣列流形矩陣分別為Ar=[ar1, ar2,…,arK]和At=[at1,at2,…,atK],其中ark和atk分別為第k個目標的接收和發射導向矢量。則接收端第q個快拍的輸出信號可以表示為其中Λq是以λkq=ckqej2πfdk(q?1)T(k=1,2,…,K)為對角元素的K維對角矩陣,ckq和fdk分別為第k個目標在第q個快拍的RCS和多普勒頻率,T為脈沖重復周期。Nq為零均值加性高斯白噪聲,Q為每個相干處理間隔(Coherent Processing Interval, CPI)內的快拍數。假設發射陣列由兩個具有互質結構的子線陣構成,具有圖1所示結構。陣元間距較大的子陣用●表示,陣元間距較小的子陣用○表示,兩個子陣共用的最左端參考陣元用?表示,陣元個數分別為2Mt和Nt(Mt<Nt),陣元間距分別為Mtλ/2和 Mtλ/2。其中λ為信號波長,Mt和Nt為兩個互質的正整數,即Mt和Nt不含1以外的公約數,且滿足總陣元數M=2Mt+Nt?1。

圖1 互質發射陣列結構

設第k個目標的DOD為φtk∈[0,2π),其發射導向矢量可以表示為集合IM,N包含兩個子陣所有陣元位置。該導向矢量不滿足范德蒙德結構,無法使用常規的角度估計算法。因此構造矩陣Rtk=atka,其元素可以表示為。根據互質數性質,由p=im?in組成的集合?包含了區間[?MtNtMtNt]上的所有整數[13]。去除Rtk中冗余和非均勻的元素即可構造向量

其等效于位于[?MtNtMtNt]處陣元間距為λ/2的均勻虛擬線陣的導向矢量,再利用ESPRIT算法即可估計出目標發射角。綜上所述,虛擬陣列技術可以利用2Mt+Nt?1個陣元的稀疏互質陣列構造出2MtNt+1個陣元的虛擬滿陣。同理,在接收端可以得到2MrNr+1個陣元的虛擬接收陣列。

將上述互質陣列的思想拓展到目標多普勒頻率估計中,假設目標滿足Swerling-I模型,其RCS在一個CPI內保持不變,即ckq≡ck(q=1,2,…,Q)。通過將式(1)中Yq列向量化得到yq=vec(Yq),令Y = [y1,…,yQ],即可以構造出整個CPI內Q個快拍的數據矩陣。

其中⊕表示Khatri-Rao乘積。矩陣ΛT的第q列為矩陣Λq的對角線元素構造成的列向量λq=[λ1q,…,λkq,…,λKq]T。N=[vec(N1),…,vec(Nq),…,vec(NQ)], vec(·)表示矩陣的列向量化。當q取[1 Q]上的均勻序列,且相鄰兩快拍間隔滿足奈奎斯特采樣條件T≤1/(2fd)時,即可利用ESPRIT算法從矩陣Λ的第k個列向量λk中得到該目標多普勒頻率的不模糊估計。為了降低數據存儲成本,僅采樣L個快拍的數據,快拍序列具有類似于收發陣列的互質結構,滿足lm,n∈LM,N={nlMl,0≤nl≤Nl?1}∪{mlNl, 0≤ml≤2Ml?1},Ml和Nl為一對互質數。此時,式(4)中矩陣Y的維數降低為(MN)×L,但仍可從2MlNl+1個連續的虛擬快拍中估計目標的多普勒頻率。

基于上述回波模型,3.1節提出了一種在非均勻陣列以及非均勻快拍條件下估計目標DOD, DOA和多普勒頻率的算法,并在3.2節中進行優化以降低估計誤差。

3 互質陣MIMO雷達參數估計算法

3.1 目標多參數估計

為了從式(4)中得到tA,Ar和Λ的估計,在此引入張量分解算法。文獻[14]給出了張量的矩陣表示法。

定義1[14]N階張量X∈?I1×I2×…×IN的P模式展開矩陣X(P)可以表示為

其中A(1)~A(N)稱為張量X的N個因子矩陣。由上述定義可得,式(4)中矩陣Y為三階張量的二模式展開矩陣,該張量包含3個因子矩陣At,Ar和Λ。因此可將矩陣Y重構成三階張量Y∈?M×N×L,其元素滿足ymnl=[Y](m?1)·N+n,l, [·]a,b表示矩陣第a行第b列的元素。可見ymnl表示第m路發射信號經過多目標散射在第n個接收陣元處得到的第l個快拍的輸出數據,因此可以改寫為

其中[·]a表示向量的第a個元素,emnl為噪聲張量中的相應元素。

定義2[14]若張量X∈?I1×I2×…×IN為N個向量的外積,即張量的每個元素等于這N個向量相應元素的乘積,則X為秩1張量。

根據上述定義和式(6),張量Y可以表示為K個秩1張量的和:

“°”表示向量外積。atk,ark和λk分別為矩陣At,Ar和Λ的第k列。張量分解的目的就是從具有式(7)結構的張量Y中得到3個因子矩陣的估計,和Λ?。交替最小二乘算法(Alternating Least Squares, ALS)可以求解這一問題[15]。求得的矩陣與At的列向量存在幅度和順序上的模糊。記為=AtΓ1Γ2, Γ1為列置換矩陣,Γ2=diag(c1,…,ck,…,cK)為列加權矩陣,ck為隨機復常數。為了構造形如式(3)的虛擬導向矢量,令

其中Πt為(2MtNt+1)×(2Mt+Nt?1)2維選擇矩陣,使得選取后的矩陣 等效于[?MtNtMtNt]處虛擬發射陣列的陣列流形。令Πt1,Πt2分別為矩陣Πt的前2MtNt行和后2MtNt行,[·]fr和[·]lr分別表示矩陣的第一行和最后一行,由式(3)可得

滿足Πt2=Πt1·Φt,Φt是以(e?jπcosφt1,…, e?jπcosφtk,…,e?jπcosφtK)為對角元素的對角矩陣。則第k個目標的DOD可以利用ESPRIT算法求得

其中arg[·]表示取復數幅角,(·)+表示矩陣偽逆。以此類推,可以利用估計得到的因子矩陣和構造出形如式(3)的向量和,求得第k個目標的DOA和多普勒頻率

討論1 均勻快拍下的目標多普勒頻率估計算法要求發射信號PRF滿足奈奎斯特采樣條件fPRF≥2fdk。本文采用的非均勻快拍序列由2個均勻稀疏快拍序列構成,具有類似圖的互質結構。其PRF分別為fPRF1=fPRF/Ml和fPRF2=fPRF/Nl,均不滿足上述采樣條件,常規算法下得到的速度估計會存在周期性模糊。但得益于Ml和Nl的互質關系,構造出的2MlNl+1個虛擬快拍均勻分布,其PRF滿足=fPRF。因此,僅要求≥2fdk即可解決由非均勻稀疏快拍造成的速度模糊問題。

3.2 改進的參數估計算法

由式(7)可知,在噪聲環境下通過ALS算法得到3個因子矩陣的估計t,r和存在一定誤差。令tk=ck(atk+σntk),其中σntk為該導向矢量的估計誤差,且σ?1。虛擬均勻發射線陣的導向矢量的估計值為

?表示kronecker積,該運算使得張量分解產生的估計誤差被放大,從而影響式(10)中DOD估計精度。針對這一問題,提出了一種基于特征值分解的誤差抑制算法。

其中0a×b表示a×b維全零矩陣,0≤i≤MtNt。由式(3)和式(13)可以構造:

討論2 張量分解參數估計算法的最大可估目標數由張量因子矩陣的秩決定。當目標個數K滿足2K+2≤kAt+kAr+kΛ時,式(7)中張量Y分解出的K個秩1張量是唯一的。其中ktA,kAr和kΛ分別表示式(4)中矩陣tA,rA和Λ的秩。每個秩1張量僅包含單一目標的待估參數信息,因此式(10)式(11)估得的參數不需要額外的配對算法。3.2節中的優化算法同樣也是基于單個秩1張量的運算,所以式(16)中的各估計值也是一一對應的。

討論3 由式(9)和式(10)可知旋轉不變因子ηtk=ej?=ej2πdtcosφtk /λ的估計誤差影響了方向余弦utk=cosφtk的估計精度,從而決定了算法的角度估計性能。而ηtk的估計值又由相鄰兩虛擬陣元的信號相位誤差Δ?和陣元間距誤差Δd決定,即= ej2π(dt+Δd)utk /λ+Δ?。不難得到方向余弦utk的估計誤差為

可見參數估計誤差將不再僅僅受到回波噪聲的影響。而3.2節所提改進算法主要作用是抑制加性噪聲,因此本文需要建立在兩項測量誤差Δ?和Δd均為0,僅存在加性高斯白噪聲時的理想情況下。

表1列出了算法的基本步驟。

4 仿真實驗

本節通過蒙特卡羅實驗驗證算法有效性,并對比了文獻[8],文獻[12],本文3.1節和本文3.2節所提算法的估計性能。文獻[8]采用均勻線型收發陣列,文獻[12]和本文算法采用相同結構的非均勻互質線陣。文獻[8]和文獻[12]采用均勻采樣,快拍數為Q,本文算法采樣的快拍序號服從第2節所述的互質結構,共計L個快拍。雷達工作載頻為1 GHz,發射信號為正交的Hadamard編碼信號,脈沖內編碼數為128,脈沖重復周期均為100 μs。3個遠場不相關目標的DOD, DOA和多普勒頻率分別為:φt= {70.2738°,65.4826°,55.5839°},φr={65.2485°,55.3847°, 44.3847°}和fd={2 kHz, 3.5 kHz, 1.5 kHz}。噪聲為獨立零均值加性高斯白噪聲。假設收發陣列結構相同,定義角度估計均方根誤差(Angle Root Mean Square Error, ARMSE)和多普勒頻率估計均方根誤差(Doppler Root Mean Square Error, DRMSE)為

表1 算法基本步驟

實驗次數Nm=100。各圖中短虛線、長虛線、實線和點畫線分別表示本文3.1節、本文3.2節、文獻[12]和文獻[8]算法的仿真數據。

圖2為4種算法估計誤差隨信噪比變化曲線。文獻[8]和文獻[12]算法的快拍數為20,本文所提2種算法的快拍序列具有互質結構,互質數對為Ml=5和Nl=11,實際采樣的快拍數也等于20。文獻[8]采用6陣元均勻收發陣列,陣元間距為λ/2。文獻[12]和本文2種算法均采用圖1結構的收發陣列,互質數對為Mt=2和Nt=3。信噪比變化范圍為[0 dB, 30 dB]。由圖2(a)可知,本文及文獻[12]由于收發陣列的互質結構形成了更多的虛擬陣元,獲得了比文獻[8]中均勻線陣更高的角度估計精度。本文3.1節的算法由于運算過程中誤差被放大,角度估計精度低于文獻[12]中的2維空間平滑算法。經過本文3.2節算法的優化,最終的估計精度略高于文獻[12]所述算法。由圖2(b)可知,文獻[12]的算法無法估計目標的多普勒頻率。得益于快拍序列的互質結構,本文3.1節算法具有比文獻[8]更高的多普勒估計精度。并且在3.2節算法優化后估計誤差進一步減小。

圖3為4種算法估計誤差隨實際采樣的快拍數變化曲線。收發陣列結構同上一實驗相同,信噪比為0 dB。實際采樣快拍數的變化范圍為[20, 100]。本文2種算法快拍序號具有互質結構,為了保證同文獻[8]和文獻[12]的快拍數相同,互質數對分別為Ml={5,10,15,20,25}和Nl={11,21,31,41,51}。由圖3(a)可知,本文3.2節的算法可以獲得遠高于文獻[8]且略高于文獻[12]角度估計精度,這一優勢在低快拍條件下仍然保持。由圖3(b)可知,在低快拍下本文3.2節的算法可以利用虛擬快拍提高多普勒頻率估計精度,且提升效果隨著快拍數的增加而增大。

5 結束語

本文提出了一種基于非均勻線型收發陣列的雙基地MIMO雷達目標DOD, DOA和多普勒頻率估計算法。收發陣列均采用非均勻互質線陣,且一個CPI內采樣的快拍序號也具有互質結構。算法首先利用張量分解從測量張量中求取3個因子矩陣,從中構造出虛擬均勻線陣的陣列流形矩陣和虛擬均勻采樣的多普勒信息矩陣。然后構造虛擬導向矢量的自相關矩陣,通過特征值分解得到誤差抑制后的虛擬導向矢量。最后通過ESPRIT算法求得目標的DOD, DOA和多普勒頻率。仿真結果表明算法在相同陣元數和快拍數條件下能得到比普通MIMO雷達更高的參數估計精度,且各參數間無需額外的配對算法。

圖2 估計誤差與信噪比關系曲線

圖3 估計誤差與快拍數關系曲線

[1] Rabideau D J and Parker P. Ubiquitous MIMO multifunction digital array radar[C]. Proceedings of the Thirty-seventh Asilomar Conference on Signals, Systems and Computers, California, USA, 2003: 1057-1064.

[2] Kai L and Manikas A. Superresolution multitarget parameter estimation in MIMO radar[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(6): 3683-3693.

[3] Gong Peng-cheng, Shao Zhen-hai, Tu Guang-peng, et al.. Transmit beampattern design based on convex optimization for MIMO radar systems[J]. Signal Processing, 2014, 94(3): 195-201.

[4] Li L. Cramer-Rao bound for parameter estimation in narrowband bistatic MIMO radar[J]. Applied Mechanics and Materials, 2014, 56(2): 5034-5037.

[5] Zhang Xiao-fei, Xu Ling-yun, Xu Lei, et al.. Direction of departure (DOD) and direction of arrival (DOA) estimation in MIMO radar with reduced-dimension MUSIC[J]. IEEE Communications Letters, 2010, 14(12): 1161-1163.

[6] Chen Duo-fang, Chen Bai-xiao, and Qin Guo-dong. Angle estimation using ESPRIT in MIMO radar[J]. Electronics Letters, 2008, 44(12): 770-771.

[7] Chen Jin-li, Gu Hong, and Su Wei-min. Angle estimation using ESPRIT without pairing in MIMO radar[J]. Electronics Letters, 2008, 44(24): 1422-1423.

[8] Chen Yuan-bing, Gu Hong, and Su Wei-min. Joint 4-D angle and doppler shift estimation via tensor decomposition for MIMO array[J]. IEEE Communications Letters, 2012, 16(6): 917-920.

[9] Chen Chun-yang and Vaidyanathan P P. Minimum redundancy MIMO radars[C]. Proceedings of the IEEE International Symposium on Circuits and Systems, California, USA, 2008: 45-48.

[10] Pal P and Vaidyanathan P P. Nested arrays: a novel approach to array processing with enhanced degrees of freedom[J]. IEEE Transactions on Signal Processing, 2010, 58(8): 4167-4181.

[11] Pal P and Vaidyanathan P P. Coprime sampling and the music algorithm[C]. Proceedings of the Digital Signal Processing Workshop and IEEE Signal Processing Education Workshop (DSP/SPE), California, USA, 2011: 289-294.

[12] Yao Bo-bin, Wang Wen-jie, and Yin Qin-ye. DOD and DOA estimation in bistatic non-uniform multiple-input multipleoutput radar systems[J]. IEEE Communications Letters, 2012, 16(11): 1796-1799.

[13] Vaidyanathan P P and Pal P. Sparse sensing with co-prime samplers and arrays[J]. IEEE Transactions on Signal Processing, 2011, 59(2): 573-586.

[14] S?rensen M and Comon P. Tensor decompositions with banded matrix factors[J]. Linear Algebra and Its Applications, 2013, 38(2): 919-941.

[15] Lü Hui and Zhang Meng. Angle and Doppler estimation using alternating least squares method in bistatic MIMO radar[C]. Proceedings of the IET International Radar Conference 2013, Xi,an, China, 2013: 987-990.

樊勁宇: 男,1985年生,博士生,研究方向為MIMO雷達陣列信號處理、目標定位等.

顧 紅: 男,1967年生,博士,教授,博士生導師,研究方向為噪聲雷達、MIMO雷達信號處理、雷達成像、目標識別等.

蘇衛民: 男,1959年生,博士,教授,博士生導師,研究方向為陣列信號處理、雷達成像等.

Co-prime MIMO Radar Multi-parameter Estimation Based on Tensor Decomposition

Fan Jin-yu①Gu Hong①Su Wei-min①Chen Jin-li②

①(School of Electronic Engineering and Optoelectronic Technology, Nanjing University of Science and Technology, Nanjing 210094, China)

②(Colllege of Electronic and Information Engineering, Nanjing University of Information Science and Technology, Nanjing 210044, China)

A novel algorithm for estimation of Direction Of Departure (DOD), Direction Of Arrival (DOA), and Doppler frequency based on bistatic MIMO radar with Co-Prime Array (CPA) is presented. The transmit and receive arrays are both composed of a pair of sparse uniform subarrays. Similarly, a pair of snapshot sequences with co-prime intervals constitutes the sampling of temporal. Three manifold matrices which contain multi-targets' DODs, DOAs and Doppler frequencies respectively are estimated through tensor decomposition. From which a group of Vandermonde matrices of virtual manifold are constructed. To improve the estimation accuracy, an error depressing algorithm based on eigenvalue decomposition is proposed. Finally, the above three parameters are estimated by an Estimation of Signal Parameters via Rotation Invariant Techniques (ESPRIT) algorithm. The proposed algorithm offers better performance through virtual array and virtual snapshot without parameter ambiguous. It requires neither peak searching nor pairing processes, and the simulation results are presented to verify the effectiveness of the proposed algorithm.

Bistatic MIMO radar; Co-prime pair; Tensor decomposition; Swerling-I model

TN958

: A

:1009-5896(2015)04-0933-06

10.11999/JEIT140826

2014-06-23收到,2014-11-06改回

國家部委基金,教育部博士點基金(20113219110018),中國航天科技集團公司航天科技創新基金(CASC04-02),國家自然科學基金(61302188)和江蘇省自然科學基金(BK20131005)資助課題

*通信作者:樊勁宇 chisame904@gmail.com

主站蜘蛛池模板: 热99精品视频| 91精品国产一区自在线拍| 国产无码性爱一区二区三区| 97成人在线观看| 久久激情影院| 九九热免费在线视频| 99久久精品久久久久久婷婷| 911亚洲精品| 国产真实乱人视频| 国产视频大全| 欧美成人午夜视频免看| 色综合激情网| 色国产视频| 日本国产精品| 黄色国产在线| 天堂成人av| 久久精品aⅴ无码中文字幕| 欧美精品1区2区| 九九精品在线观看| 91在线视频福利| 欧洲高清无码在线| 91精品国产情侣高潮露脸| 国产手机在线ΑⅤ片无码观看| 中国精品久久| 国产高清又黄又嫩的免费视频网站| 40岁成熟女人牲交片免费| 国产免费羞羞视频| 亚洲黄色片免费看| 亚洲AV无码久久精品色欲| 亚洲欧美极品| 成年A级毛片| 91国内视频在线观看| 91欧美在线| 国产午夜无码专区喷水| 这里只有精品在线| 精品小视频在线观看| 99热这里只有精品2| 国产91熟女高潮一区二区| 亚洲精品在线观看91| 国产精品免费福利久久播放| 亚洲一级毛片在线观| 精品一区二区三区视频免费观看| 亚洲第一视频网| 精品无码一区二区在线观看| 欧美精品一区二区三区中文字幕| 国产视频久久久久| P尤物久久99国产综合精品| 久久视精品| 欧美69视频在线| 精品人妻系列无码专区久久| 国产美女无遮挡免费视频网站| www中文字幕在线观看| 又污又黄又无遮挡网站| 伊人中文网| 国产不卡一级毛片视频| 国产午夜在线观看视频| 婷婷六月在线| 国产午夜精品一区二区三| 亚洲人成人伊人成综合网无码| 热热久久狠狠偷偷色男同| 日日拍夜夜操| 天天综合网站| www成人国产在线观看网站| 最近最新中文字幕在线第一页 | 国产在线麻豆波多野结衣| 日韩在线网址| 国产精品成人一区二区不卡 | 成人综合网址| 亚洲最猛黑人xxxx黑人猛交| 丁香六月综合网| 四虎亚洲国产成人久久精品| 国产三级毛片| 精品亚洲麻豆1区2区3区| 在线免费亚洲无码视频| 在线播放精品一区二区啪视频| 国产无码高清视频不卡| 一本大道在线一本久道| 四虎永久免费在线| 欧美成人精品在线| 国产精品妖精视频| 久久国产精品娇妻素人| 99视频精品全国免费品|