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

基于P H D濾波的相控陣雷達多目標跟蹤算法

2016-09-20 08:20:12袁常順孫進平畢嚴先北京航空航天大學電子與信息工程學院北京100191
系統(tǒng)工程與電子技術 2016年3期

袁常順,王 俊,雷 鵬,孫進平,畢嚴先(北京航空航天大學電子與信息工程學院,北京100191)

基于P H D濾波的相控陣雷達多目標跟蹤算法

袁常順,王 俊,雷 鵬,孫進平,畢嚴先
(北京航空航天大學電子與信息工程學院,北京100191)

對于相控陣雷達方向余弦量測,采用擴展卡爾曼概率假設密度(extended Kalman probability hypothesis density,E K-P H D)濾波進行多目標跟蹤時,存在目標數(shù)估計偏高和目標狀態(tài)估計準確度低的問題。針對上述問題,提出了一種新的多目標跟蹤算法——無偏轉換量測概率假設密度(unbiased converted measurements P H D,U BC M-P H D)濾波算法。該算法采用方向余弦量測下的量測轉換方法,保留了更多的量測信息;同時對轉換后的量測偏差進行補償,使量測轉換誤差的均值、方差準確近似原始量測高斯分布的一、二階矩。仿真實驗表明,所提算法可提高目標數(shù)和目標狀態(tài)估計準確性。

相控陣雷達;多目標跟蹤;無偏轉換量測;隨機有限集;概率假設密度濾波

網(wǎng)址:www.sys-ele.co m

0 引 言

在多目標跟蹤中,目標數(shù)的變化和量測信息的不確定給多目標跟蹤帶來巨大的困難。如何有效、實時地實現(xiàn)多目標跟蹤,一直是目標跟蹤領域研究的熱點和難點[1 3]。傳統(tǒng)的多目標跟蹤方法采用數(shù)據(jù)關聯(lián)技術,但隨著目標數(shù)或者雜波數(shù)的增加,其計算量指數(shù)增長,如:多假設跟蹤(multiple hypothesis tracking,M H T)[2,4 5]、聯(lián)合概率數(shù)據(jù)互聯(lián)(joint probabilistic data association,JP D A)[1,6]、概率多假設跟_蹤(probabilistic_ M H T,P M H T)[7]和多目標粒子濾波[8]。

近年來,Mahler以有限集統(tǒng)計學為基礎,提出了隨機有限集(random finite set,RFS)的方法。該方法提供了一個簡潔的公式描述多目標跟蹤,避免了數(shù)據(jù)關聯(lián),很快成為了多目標跟蹤的研究熱點之一[3,9 23]。通過隨機有限集建模多目標狀態(tài)和量測,多目標跟蹤轉換為對多目標后驗密度估計的貝葉斯濾波。但由于多目標密度的組合特性和狀態(tài)、量測空間的多維積分,導致多目標貝葉斯濾波在許多實際應用中無法實現(xiàn)[9 10]。為了解決這個困難,概率假設密度(probability hypothesis density,P H D)濾波作為多目標貝葉斯濾波的一階矩近似被提出[9]。P H D濾波的優(yōu)點是在單目標狀態(tài)空間計算,避免數(shù)據(jù)關聯(lián)。現(xiàn)有PHD濾波實現(xiàn)主要包括:高斯混合PHD濾波、擴展卡爾曼PHD(extended Kalman PHD,EKPHD)濾波、不敏卡爾曼PHD(unscented Kalman PHD,UKPHD)濾波[11]、粒子P H D濾波[10],以及改進算法[12 14]。到目前為止,P H D濾波已經(jīng)得到廣泛應用,例如:即時定位與地圖構建[15]、機動目標跟蹤[16]、雷達目標跟蹤[17]、圖像序列跟蹤[18]、擴展目標跟蹤[19]、主動聲學跟蹤[20]、醫(yī)學視頻跟蹤[21]等。

但是,對于P H D濾波在相控陣雷達應用的研究較少。相控陣雷達中,量測坐標通常為方向余弦坐標,因在方向余弦坐標下,天線方向圖不隨掃描角變化。采用E K-P H D濾波進行多目標跟蹤時,由于E K-P H D濾波采用一階泰勒級數(shù)近似,忽略高階項,在濾波過程中引入了近似誤差。針對這一問題,本文提出一種方向余弦坐標下無偏轉換量測概率假設密度(unbiased converted measurements,U BC M-P H D)濾波的相控陣雷達多目標跟蹤算法。該算法充分利用P H D濾波和無偏轉換的優(yōu)點,提高了濾波和目標數(shù)估計精度。算法首先根據(jù)非線性量測下的無偏轉換方法,給出了方向余弦坐標下的無偏轉換方法;然后將上述無偏轉換方法與P H D濾波結合實現(xiàn)U BC M-P H D濾波。實驗仿真驗證了該算法對相控陣方向余弦坐標下多目標跟蹤的準確性和有效性。本文的方法主要通過P H D濾波來實現(xiàn),可直接被推廣到勢概率假設密度濾波[22]和勢平衡多目標多伯努利濾波[23]。

1 RFS及P H D濾波原理

假設k時刻有M(k)個狀態(tài)為xk,1,…,xk,M(k)的目標。傳感器在k 時刻接收到N (k)個量測zk,1,…,zk,N(k)。在這些量測中,一部分是真實目標產生的量測,一部分是虛假量測,同時有些真實目標產生的量測無法被檢測。基于RFS的多目標模型對目標狀態(tài)和量測沒有明確的順序要求,所以其可被表示為有限集的形式,即

若k-1時刻目標狀態(tài)為Xk-1,則k時刻的目標狀態(tài)Xk可表示為

式中,Sk|k-1(ζ)為k-1時刻到下一時刻存活目標狀態(tài)R FS;Bk|k-1(ζ)為k時刻由上一時刻狀態(tài)ζ衍生目標狀態(tài)R FS;Γk為k時刻新生目標狀態(tài)R FS[9]。

k時刻接收量測Zk可表示為

式中,Κk為雜波量測RFS;Θk(x)為真實目標產生量測RFS[9]。

根據(jù)最優(yōu)多目標貝葉斯濾波理論,基于R FS的多目標后驗概率密度的遞推表達式表示為

式中,fk|k-1(·|·)為狀態(tài)轉移概率密度函數(shù);gk(·|·)為量測似然函數(shù);μs為狀態(tài)空間近似Lebesgue測度[9]。

由于式(3)和(4)多維積分的計算量龐大,實現(xiàn)復雜,文獻[9]提出了利用多目標概率密度函數(shù)的一階矩P H D近似多目標后驗概率,實現(xiàn)對目標數(shù)時變且未知的多目標跟蹤。令vk|k-1和vk分別為k時刻的預測P H D和后驗P H D,則其可分別表示為

式中,βk|k-1(x)和γk(x)分別為衍生和新生目標狀態(tài)R FS的預測P H D;pS,k(ζ)為k-1時刻狀態(tài)ζ的目標在下一時刻存活概率;pD,k(x)為k時刻狀態(tài)x的目標檢測概率;κk(·)表示k時刻服從泊松分布雜波R FS的P H D[9]。

2 U BC M-P H D濾波

2.1 方向余弦坐標量測無偏轉換

在多目標跟蹤中,目標的狀態(tài)通常采用笛卡爾直角坐標系描述,但傳感器接收量測通常位于極坐標系或球坐標系。為了實現(xiàn)極坐標系或球坐標系的量測在笛卡爾直角坐標系下跟蹤,量測轉換方法被廣泛的使用。近年來,出現(xiàn)了大量關于量測轉換的改進算法[24 25],其差異在于量測轉換偏差和方差的估計方法。相控陣雷達由于陣列天線是不動的,當掃描波束偏離法線方向時,波速將會展寬,同時波束形狀也有變化,若采用極坐標來分析會十分復雜,因此提出了方向余弦坐標。針對相控陣雷達方向余弦坐標系下的量測,本文給出了一種無偏量測轉換方法,準確近似了原始量測高斯分布的一、二階矩。

相控陣雷達量測為Zmcos=(Rmαm)T,假設真值為R,α,有

方向余弦量測噪聲Vk=()T為零均值高斯白噪聲,其方差Rcos各分量為,。方向余弦坐標轉換到直角坐標為

類似于式(7),得到直角坐標真值為

將式(8)帶入式(9),化簡后得到的轉換誤差為因此,直角坐標系下誤差不再獨立,且該誤差與方向余弦、距離量測及原始量測誤差有關。此外,的統(tǒng)計特性未知,且β與α有如下的非線性關系為了使方向余弦坐標系量測轉換誤差的均值、方差準確近似原始量測高斯分布的一、二階矩,本文利用β在αm的四階展開β4來近似,其誤差為,則量測轉換偏差和方差如下:

式(12)中

最終得到的無偏轉換量測為Zcar=(αmRmβmRm)T+μ,方差為R。

2.2 U BC M-P H D

針對相控陣雷達的特點和E K-P H D濾波不足,本文提出了采用方向余弦坐標下無偏轉換量測的方法和基于隨機有限集的概率假設密度濾波對多目標進行跟蹤,在保證跟蹤精度的同時,實現(xiàn)對目標數(shù)精確的估計。該算法主要包含預測和更新兩部分:預測部分與標準的高斯混合P H D濾波相同;更新部分采用無偏轉換量測的方法將方向余弦坐標量測轉換為直角坐標量測更新目標狀態(tài),具體如下。

預測 假設k-1時刻多目標后驗概率密度的P H D為高斯混合形式,且可表示為

假設新生目標和衍生目標的P H D也分別為高斯混合形式,且表示為

根據(jù)前面的假設,k時刻預測的P H D也為高斯混合的形式,可表示為

式中

pS,k為目標存活概率。

更新 假設k-1時刻預測P H D為高斯混合形式

k時刻,對于每一觀測z∈Zk,采用第2.1節(jié)給出的無偏轉換量測的方式,求出該量測對應的轉換偏差μk(z)和方差Rk(z)。則更新后的P H D可表示為

式中

該方法相比于文獻[11]中的E K-P H D濾波,通過將量測誤差轉移到與目標狀態(tài)相同的坐標中,在轉換中由于各個變量之間不再是相互獨立,因此轉換后的方差可包含更多的信息。但由于在轉換的過程中存在線性化誤差的問題,采用去偏方式對轉換后的量測進行補償,同時給出了對應的轉換方差,這樣保證在轉換過程中信息丟失最少,準確近似原始量測高斯分布的一、二階矩,提高了算法跟蹤性能。

3 仿真分析

采用一個目標數(shù)目時變且未知的二維相控陣雷達多目標跟蹤場景來驗證本文算法。目標的運動和量測模型分別表示為

仿真場景中共包含12個目標,起始時刻存在2個勻速運動的目標,接著分別有10個勻速和直角轉彎運動的新生目標出現(xiàn),真實運動參數(shù)如表1所示。新生目標R FS的P H D為

表1 目標真實運動參數(shù)

采樣間隔取T=1 s,共采樣100個時刻,對于E K-P H D和U BC M-P H D兩種濾波算法的最多高斯分量數(shù)為100。仿真場景如圖1所示,其中和□分別表示目標運動起始和結束位置。

圖1 包含雜波量測的目標真實軌跡

圖2給出了U BC M-P H D濾波算法單次仿真實驗對多目標位置的估計結果。從仿真結果可以看出給出的算法可在大量虛警雜波中正確跟蹤單獨目標運動和不同的目標新生和消失量測。

圖2 U B C M-P H D濾波位置估計

通過采用目標航跡固定但雜波量測隨機產生的100次蒙特卡羅實驗來驗證算法的性能。圖3和圖4分別給出了E K-P H D和U BC M-P H D濾波算法的目標數(shù)估計均值和標準差比較結果。從圖中可以看出,U BC M-P H D濾波算法的目標數(shù)估計均值更接近于真實目標數(shù),E K-P H D濾波算法的目標數(shù)估計均值偏高;U BC M-P H D濾波算法的目標數(shù)估計標準差小于E K-P H D濾波算法。因為U BC M-P H D濾波算法在進行量測更新時,采用轉換量測的方式將量測的誤差全部轉移到與運動狀態(tài)相同的坐標系下,保留了更多的量測信息;同時本算法對轉換后量測偏差進行補償,使量測轉換誤差其均值、方差準確近似原始量測高斯分布的一、二階矩,改善了濾波性能。但E K-P H D濾波算法采用一階泰勒級數(shù)近似,忽略高階項,在濾波過程中引入了近似誤差,同時丟失部分量測信息。

圖3 E K-P H D和U B C M-P H D目標數(shù)估計均值

圖4 E K-P H D和U B C M-P H D目標數(shù)估計標準差

利用最優(yōu)子模式分配(optimal subpattern assign ment O SP A)距離對以上算法性能進行綜合評價,其定義為[26]

式中,p=1,c=300。

圖5給出了兩種算法的O SP A距離比較結果。從圖中可知U BC M-P H D濾波算法的O SP A距離小于E K-P H D濾波算法,進一步表明本文提出的算法具有更高的跟蹤精度;同時發(fā)現(xiàn)在目標數(shù)發(fā)生變化的時刻出現(xiàn)峰值,這主要是在目標新生的時刻,由于有新生目標加入,目標數(shù)出現(xiàn)變化,濾波器在該時刻對于目標數(shù)的估計偏離真實目標數(shù)(如圖3所示)。

圖5 E K-P H D和U BC M-P H D的O SP A距離

圖6給出了在不同雜波密度ρc下兩種濾波跟蹤性能,其中檢測概率pD,k=0.98。采用O SP A距離來進行比較,從比較結果可顯著的發(fā)現(xiàn)U BC M-P H D濾波算法的性能優(yōu)于E K-P H D濾波算法;同時隨著雜波密度ρc的增加,兩種算法的性能都降低。

圖7給出了不同檢測概率pD,k下兩種濾波跟蹤性能,其中雜波密度固定為ρc=31.25×10-3個/k m2。仍然采用O SP A距離來進行比較,從比較結果可以顯著的發(fā)現(xiàn)兩種濾波算法的跟蹤性能隨著pD,k的減小而變差。這主要是由于P H D濾波必須在解決比較高的檢測不確定性的基礎上接著解決目標數(shù)不確定,因此隨著檢測不確定的增加,目標的不確定性也增加了。但由于本文算法考慮了轉換量測后偏差的補償,因此跟蹤性能會優(yōu)于E K-P H D濾波算法。

圖6 不同雜波率下的O SP A距離

圖7 不同檢測概率下的O SP A距離

4 結束語

在相控陣雷達中,量測坐標通常為方向余弦坐標,采用E K-P H D濾波算法進行多目標跟蹤時,由于E K-P H D濾波采用一階泰勒級數(shù)近似,忽略高階項,在濾波過程中引入了轉換誤差。針對這一問題,本文提出了一種基于隨機有限集的U BC M-P H D濾波算法,該算法采用轉換量測的方式將量測誤差轉移到與運動狀態(tài)相同的坐標系下,保留了更多的量測信息;同時對轉換后的量測偏差進行補償,使量測轉換誤差的均值、方差準確地近似原始量測高斯分布的一、二階矩。仿真結果表明,本文所提算法可提高目標狀態(tài)和目標數(shù)估計準確度。該算法主要通過P H D濾波來實現(xiàn),可以直接被推廣到勢概率假設密度濾波和勢平衡多目標多伯努利濾波。

[1]Bar-Shalom Y,F(xiàn)ortmann T E.Tracking and data association[M]. San Diego:Academic,1988.

[2]Black man S.Multipletargettracking with radar applications[M]. Norwood:Artech House,1986.

[3]M ahler R.Statistical multisource-m ultitarget inform ation fusion[M].Norwood:Artech H ouse,2007.

[4]Reid D.A n algorith m for tracking m ultiple targets[J].IE E E Trans.on Autom atic Control,2004,24(6):843-854.

[5]Black man S.M ultiple hypothesis tracking for m ultiple target tracking[J].IE E E Aerospace and Electronic Systems M agazine,2004,19(1):5-18.

[6]Fortmann T E,Bar-Shalo m Y,Scheffe M.Sonar tracking of m ultiple targets using joint probabilistic data ssociation[J]. IE E E Journal of Oceanic Engineering,1983,8(3):173-184.

[7]Streit R L,Luginbuhl T E.M axim u m likelihood method for probabilistic m ulti-hypothesis tracking[C]∥Proc.of the S PIE,1994:5-7.

[8]H ue C,Lecadre J P,Perez P.Sequential M onte Carlo methods for m ultiple target tracking and data fusion[J].IE E E Trans.on Signal Process,2002,50(2):309-325.

[9]M ahler R.M ultitarget Bayes filtering via first-order m ultitarget m o ments[J].IE E E Trans.on Aerospace and Electronic Systems,2003,39(4):1152-1178.

[10]V o B N,Singh S,Doucet A.Sequential M onte Carlo methods for m ultitarget filtering with rando m finite sets[J].IE E E Trans.on Aerospace and Electronic Systems,2005,41(4):1224-1245.

[11]Vo B N,M a W.T he Gaussian mixture probability hypothesis density filter[J].IE E ETrans.on Signal Process,2006,54 (11):4091-4104.

[12]Yang J L,Ji H B,Liu J M.Gauss-H ermite particle P H D filter for bearings-only m ulti-target tracking[J].Systems Engineering and Electronics,2013,35(3):457-462.(楊金龍,姬紅兵,劉進忙.高斯厄米特粒子P H D被動測角多目標跟蹤算法[J].系統(tǒng)工程與電子技術,2013,35(3):457-462.).

[13]Chen L M,Chen Z,Yin F L,et al.Central difference Kalman probability hypothesis density filter for multi-target tracking[J]. Controland Decision,2013,28(1):36-42.(陳里銘,陳喆,殷福亮,等.基于中心差分卡爾曼-概率假設密度濾波的多目標跟蹤方法[J].控制與決策,2013,28(1):36-42.)

[14]Yan D L,Song Y D,Song Y,et al.T he application of squareroot cubature Kalman filter and probability hypothesis density in sim ultaneous localization and mapping for m obile robots[J]. Control Theory and A pplications,2014,31(8):1009-1017.(閆德立,宋永端,宋宇,等.平方根容積卡爾曼濾波概率假設密度_省略_機器人同時定位與地圖構建中的應用[J].控制理論與應用,2014,31(8):1009-1017.)

[15]A dams M,Vo B N,M ahler R,et al.SL A M gets a P H D:new concepts in map estimation[J].IE E E Robotics&A utomation M agazine,2014,22(1):26-37.

[16]Luo S H,Xu H,Xu Y,et al.U T based M M P H D filter for tracking maneuvering targets[J].Systems Engineering and Electronics,2012,34(4):666-672.(羅少華,徐暉,徐洋,等.基于U T變換的M M P H D機動目標跟蹤[J].統(tǒng)工程與電子技術,2012,34(4):666-672.)

[17]Jan P,To mas S,Zdenek N,et al.A n optimization of a P H D function for association of targets on m ultistatic radar[C]∥Proc.of the IE E E Radar Conference,2014:1084-1089.

[18]V o B N,Vo B T,Pham N T,et al.Joint detection and estimation of m ultiple objects fro m image observations[J].IE E E Trans.on Signal Procesing,2010,58(10):5129-5241.

[19]Tian S P,Zhou B,Qi Q F.Gaussian mixture P H D filter based tracking m ultiple maneuvering extended targets[J].Journalof Central South University(Science and Technology),2013,44 (12):4923-4929.(田森平,周波,戚其豐.基于高斯混合P H D濾波的多機動擴展目標跟蹤[J].中南大學學報(自然科學版),2013,44(12):4923-4929.

[20]Clark D,Vo B N,Bel l J.G M-P H D filter multitarget tracking in sonarimages[C]∥Proc.of the Defense and Security Symposium International Society for Optics and Photonics,2006:62350R-1 -62350R-8.

[21]Wood T M,Yates C A,Wilkinson D A,et al.Simplified multitarget tracking using the P H D filter for microscopic video data[J]. IE E E Trans.on Circuits and Systems for Video Technology,2012,22(5):702-713.

[22]Vo B T,Vo B N,Cantoni A.A nalytic im plementations of the cardinalized probability hypothesis density filter[J].IE E E Trans.on Signal Processing,2007,55(7):3553-3567.

[23]Vo B T,V o B N,Cantoni A.The cardinality balanced m ultitarget m ulti-bernoulli filter and its im plementations[J].IE E E Trans.on Signal Processing,2009,57(2):409-423.

[24]Zhao Z,Li X R,Ji lkov V P.Best l inear unbiased fi ltering with nonl inear measurements for target tracking[J].IE E E Trans.on Aerospace and Electronic Systems,2004,40(4):1324-1336.

[25]M ei W,H e Z,Liang G.Iterated debiased Kalman filterfor target tracking with converted measurements[C]∥Proc.of the International Conference on Inform ation Science and Technology,2012:185-189.

[26]Schuh macher D,Vo B T,Vo B N.A consistent metric for performance evaluation of multi-object filters[J].IE E E Trans.on Signal Processing,2008,56(8):3447-3457.

Multi-target tracking based on P H D filter for phased array radar

Y U A N Chang-shun,W A N G Jun,L EI Peng,S U N Jin-ping,BI Yan-xian
(School of Electronics and Information Engineering,Beihang University,Beijing 100191,China)

The extended Kalman probability hypothesis density(E K-P H D)filter has a higher biasin the estimation of the number of targets and a lower estimation accuracy oftheir states by using the direction cosine coordinate measurements for the phased array radar.To solve this problem,a novel multi-target tracking algorith m called unbiased converted measurements probability hypothesis density(U BC M-P H D)filter algorith mis proposed.The proposed algorith mutilizes the unbiased converted method to remain more information about the direction cosine coordinate measurements.Meanw hile,it compensates the bias caused by the converting direction cosine coordinate to Cartesian coordinate measurements,and the means and variances of the converted errors could accurately approximate the first-order and second-order moments of Gaussian distribution for original measurements.The simulation results indicate that the proposed algorith mim proves the estimation accuracy of both the number of targets and their states.

phased array radar;multi-target tracking;unbiased converted measurements;random finite set;probability hypothesis density filter

T N 953

A

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

1001-506 X(2016)03-0539-06

2015-01-06;

2015-09-11;網(wǎng)絡優(yōu)先出版日期:2015-11-18。

網(wǎng)絡優(yōu)先出版地址:http://www.cnki.net/kcms/detail/11.2422.T N.20151118.1208.012.html

國家自然科學基金(61171122,61201318,61471019,61501011,61501012);中央高校基本科研業(yè)務費專項資金(Y W F-15-GJS Y S-068)資助課題

袁常順(1987-),男,博士研究生,主要研究方向為雷達信號處理、隨機有限集多目標跟蹤。

E-mail:yuanchang61@126.com

王 俊(1972-),男,教授,博士,主要研究方向為雷達信號處理、實時信號處理。

E-mail:wangj203@buaa.edu.cn

雷 鵬(1985-),男,講師,博士,主要研究方向為信號處理、頻譜分析、目標識別。

E-mail:peng.lei@buaa.edu.cn

孫進平(1975-),男,教授,博士,主要研究方向為高分辨雷達信號處理、壓縮感知等。

E-mail:sunjinping@buaa.edu.cn

畢嚴先(1988-),男,博士研究生,主要研究方向為雷達信號處理、目標識別。

E-mail:biyanxian@126.com

主站蜘蛛池模板: 爱色欧美亚洲综合图区| 综合社区亚洲熟妇p| 精品无码人妻一区二区| 国产成人综合日韩精品无码不卡| 成人午夜在线播放| 亚洲区欧美区| 99久久亚洲综合精品TS| 9999在线视频| 欧美国产综合视频| 成人福利在线视频| 日韩免费成人| 国产极品美女在线播放| 91蝌蚪视频在线观看| 国内嫩模私拍精品视频| 亚洲av无码人妻| 亚洲最大综合网| 亚洲乱伦视频| 精品视频免费在线| 欧美一级专区免费大片| 久久semm亚洲国产| 91在线无码精品秘九色APP | 久久免费视频播放| 97视频免费看| 99国产精品免费观看视频| 亚洲美女一区| 麻豆精品在线播放| 婷婷色狠狠干| 成人91在线| 国产一级片网址| 一区二区偷拍美女撒尿视频| 午夜视频免费一区二区在线看| 日本色综合网| 国产男女免费完整版视频| 国产精品无码一二三视频| 青草91视频免费观看| 欧美a在线| 色亚洲成人| 久久无码高潮喷水| 欧美三级视频网站| 2021天堂在线亚洲精品专区| 制服丝袜无码每日更新| 男女男精品视频| 国产一区二区三区在线精品专区| 国产精品久久久久鬼色| 免费久久一级欧美特大黄| 一本大道视频精品人妻 | 久久中文无码精品| 五月综合色婷婷| 亚洲天堂免费| 久久国产精品无码hdav| 亚洲色精品国产一区二区三区| 91青青草视频在线观看的| 日韩精品毛片| 国产一区二区三区日韩精品| 久久香蕉国产线| 污污网站在线观看| av手机版在线播放| 国产色婷婷视频在线观看| 一级全黄毛片| 色婷婷天天综合在线| …亚洲 欧洲 另类 春色| 国产在线小视频| 在线欧美日韩| 日本一本在线视频| 久久不卡国产精品无码| 欧美人与性动交a欧美精品| 99热国产这里只有精品无卡顿" | 经典三级久久| 高清免费毛片| 国产成人精品免费视频大全五级| 日韩中文字幕免费在线观看| 午夜啪啪网| 全色黄大色大片免费久久老太| 国产性生大片免费观看性欧美| 极品尤物av美乳在线观看| 亚洲无线视频| 91免费观看视频| 亚洲精品麻豆| av无码久久精品| 亚洲成A人V欧美综合| 久久综合伊人 六十路| 男女男免费视频网站国产|