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

自旋目標(biāo)平動(dòng)補(bǔ)償與微多普勒參數(shù)估計(jì)研究

2016-07-05 09:28:33楊文革屈文星張若禹
裝備學(xué)院學(xué)報(bào) 2016年3期

楊文革, 屈文星, 張若禹

(1. 裝備學(xué)院 光電裝備系, 北京 101416; 2. 裝備學(xué)院 研究生管理大隊(duì), 北京 101416;3. 北京跟蹤與通信技術(shù)研究所, 北京 100094)

自旋目標(biāo)平動(dòng)補(bǔ)償與微多普勒參數(shù)估計(jì)研究

楊文革1,屈文星2,張若禹3

(1. 裝備學(xué)院 光電裝備系, 北京 101416;2. 裝備學(xué)院 研究生管理大隊(duì), 北京 101416;3. 北京跟蹤與通信技術(shù)研究所, 北京 100094)

摘要導(dǎo)彈目標(biāo)在飛行過程中的自旋運(yùn)動(dòng)會(huì)對(duì)連續(xù)波雷達(dá)信號(hào)產(chǎn)生微多普勒調(diào)制,對(duì)調(diào)制后信號(hào)進(jìn)行參數(shù)估計(jì)可獲得目標(biāo)的自旋頻率,而目標(biāo)因平動(dòng)產(chǎn)生的多普勒頻率會(huì)對(duì)微動(dòng)參數(shù)估計(jì)產(chǎn)生影響。針對(duì)平動(dòng)項(xiàng)引起的微多普勒頻率估計(jì)參數(shù)偏移問題,通過建立同時(shí)具有平動(dòng)與微動(dòng)項(xiàng)的載波信號(hào)模型,采用延遲共軛相乘處理補(bǔ)償平動(dòng)項(xiàng)獲得微多普勒信息,并通過匹配相關(guān)算法對(duì)微動(dòng)參數(shù)進(jìn)行估計(jì),最后給出了該算法的具體試驗(yàn)步驟。仿真結(jié)果驗(yàn)證了該方法的有效性。

關(guān)鍵詞微多普勒;連續(xù)波雷達(dá);平動(dòng)補(bǔ)償;延遲共軛相乘;參數(shù)估計(jì)

導(dǎo)彈目標(biāo)在飛行過程中為維持較好的定向性,會(huì)在平動(dòng)的同時(shí)伴有自旋,自旋運(yùn)動(dòng)會(huì)在雷達(dá)回波中產(chǎn)生邊帶調(diào)制,這種現(xiàn)象稱為微多普勒效應(yīng)[1]。在連續(xù)波雷達(dá)體制下此類邊帶調(diào)制表現(xiàn)為多普勒頻率上有規(guī)律的波動(dòng),其時(shí)頻變化規(guī)律同目標(biāo)的自旋頻率相同[2],可用于目標(biāo)旋轉(zhuǎn)頻率的估計(jì)。

關(guān)于微多普勒頻率的參數(shù)估計(jì)及提取,早期研究為便于信號(hào)處理算法的應(yīng)用,均設(shè)置目標(biāo)位于一固定點(diǎn),平動(dòng)速度為零[1-3];而導(dǎo)彈目標(biāo)飛行過程中會(huì)產(chǎn)生多普勒頻偏,多普勒頻偏在微多普勒參數(shù)估計(jì)中會(huì)引起微多普勒頻率的折疊與偏移[4]48,從而無法獲得準(zhǔn)確的自旋頻率估計(jì)值,因此在微動(dòng)參數(shù)估計(jì)前應(yīng)先完成平動(dòng)多普勒的補(bǔ)償。

文獻(xiàn)[4]49較早地針對(duì)速度對(duì)微多普勒的影響進(jìn)行了分析,該文獻(xiàn)采用中心法對(duì)速度進(jìn)行等效補(bǔ)償,但所采用的頻譜搬移方法并未考慮速度的多階變化率,也未在信號(hào)層做出多普勒頻率補(bǔ)償。文獻(xiàn)[5]1173采用延遲共軛相乘的方法,對(duì)信號(hào)中的平動(dòng)項(xiàng)有較好的補(bǔ)償效果,但在參數(shù)搜索中采用的Hough變換適合多分量信號(hào)處理,搜索參數(shù)計(jì)算量巨大,算法也較為復(fù)雜。文獻(xiàn)[6]采用了基函數(shù)匹配相關(guān)搜索譜峰的方法來估計(jì)微多普勒參數(shù),但未對(duì)調(diào)頻指數(shù)進(jìn)行估計(jì),且未考慮目標(biāo)平動(dòng)對(duì)微動(dòng)參數(shù)估計(jì)帶來的影響。

本文受已有研究結(jié)果啟發(fā),提出了一種基于延遲共軛相乘和匹配相關(guān)思想的平動(dòng)補(bǔ)償與微多普勒參數(shù)搜索方法,并利用卡森準(zhǔn)則與最優(yōu)門限法估算信號(hào)帶寬[7]355,降低了參數(shù)搜索的維度,減少了搜索算法的復(fù)雜度。仿真結(jié)果表明:該方法可實(shí)現(xiàn)對(duì)平動(dòng)多普勒的補(bǔ)償,參數(shù)估計(jì)算法較為靈活,適用于連續(xù)波雷達(dá)信號(hào)中的微多普勒參數(shù)估計(jì)。

1平動(dòng)-微動(dòng)復(fù)合信號(hào)模型

在目標(biāo)的速度測(cè)量中常采用精度較高的非調(diào)制單頻連續(xù)波,通過鎖相接收回波信號(hào)獲得多普勒信息。對(duì)于平動(dòng)多普勒頻率fT(t)可通過多項(xiàng)式近似表達(dá)[8]28

(1)

式中,ai為各次項(xiàng)系數(shù),也為第i階多普勒變化率;P為多普勒階數(shù),其各次項(xiàng)系數(shù)均與目標(biāo)運(yùn)動(dòng)狀態(tài)有關(guān)。自旋運(yùn)動(dòng)引起的微多普勒頻率fmD(t)可表示為三參數(shù)的余弦函數(shù)[8]19

(2)

式中,A為微多普勒頻偏;fm為目標(biāo)自旋頻率;φm為微多普勒初相。平動(dòng)項(xiàng)與微多普勒頻率對(duì)信號(hào)的調(diào)制作用均可看做是信號(hào)的相位調(diào)制,頻率與調(diào)制相位函數(shù)之間的關(guān)系可表示為[9]

(3)

因此通過對(duì)平動(dòng)多普勒頻率域微多普勒頻率進(jìn)行積分,可得平動(dòng)分量與微多普勒分量的復(fù)合信號(hào),表示為

(4)

式中,fc為接收載頻;系數(shù)bi=ai/(i+1);調(diào)頻指數(shù)mf=A/fm;sc(t),sT(t),sm(t)分別為載波分量、平動(dòng)多普勒分量和微多普勒調(diào)制分量。通過基帶變換可得基帶信號(hào)sb(t):

(5)

2平動(dòng)補(bǔ)償

平動(dòng)補(bǔ)償可將接收信號(hào)中的平動(dòng)分量進(jìn)行補(bǔ)償,從而提取微多普勒信號(hào)分量,在實(shí)際應(yīng)用中平動(dòng)補(bǔ)償方法較多,根據(jù)前文分析可采用延遲共軛相乘[5]1173的方法對(duì)平動(dòng)項(xiàng)進(jìn)行補(bǔ)償。補(bǔ)償原理為:延遲共軛相乘處理可使信號(hào)中的平動(dòng)項(xiàng)降階一次,通過多階延遲共軛相乘處理即可完成對(duì)高階平動(dòng)項(xiàng)的補(bǔ)償。延遲共軛相乘算法可表示為

(6)

式中,τ為時(shí)間延遲量;s*(t)為原信號(hào)的共軛形式;sD(t)為延遲共軛相乘處理后的信號(hào)。因此基帶信號(hào)的延遲共軛相乘函數(shù)sDb(t)可表示為

(7)

式中,sDT(t)與sDm(t)分別為平動(dòng)多普勒項(xiàng)與微多普勒項(xiàng)的延遲共軛相乘函數(shù),對(duì)兩項(xiàng)分別展開求解。首先平動(dòng)多普勒項(xiàng)的延遲共軛相乘函數(shù)可表示為

(8)

式中,b1,i,c1,i均為化簡后的常數(shù)系數(shù);τ為時(shí)間延遲。由于在延遲共軛相乘處理中多普勒頻率最高階ti+1得以消除,最終其sDT(t)中僅剩t的i次方項(xiàng)與直流分量,即一次延遲共軛相乘處理可使平動(dòng)多普勒分量降低一階。

因此平動(dòng)項(xiàng)的k階延遲共軛相乘函數(shù)可表示為

(9)

式中,bk,i,ck,i均為化簡后的常數(shù)系數(shù)。通過P次延遲共軛相乘即可使平動(dòng)多普勒多項(xiàng)式降低P階,因此sDT(t)可化為

(10)

式中,t的一階項(xiàng)表現(xiàn)為多普勒頻率上的直流分量2πbP,i,在頻域已不影響微多普勒頻率的幅值與變化規(guī)律,只需消除此直流分量即可。根據(jù)延遲共軛相乘的性質(zhì),一階處理可使多普勒分量降階一次,因此需要消除一階多普勒直流分量,要對(duì)信號(hào)再進(jìn)行一次延遲共軛相乘處理。可得最終補(bǔ)償后的信號(hào)為

(11)

關(guān)于微多普勒調(diào)制項(xiàng)的延遲共軛相乘函數(shù)可表示為

(12)

式中,l1為一階幅值增益因子;φ1為一階相位延遲因子,可分別表示為

(13)

可看出一階延遲共軛相乘只改變了微多普勒調(diào)頻項(xiàng)中微多普勒頻率的幅值與相位,其自旋速度并無改變。因此微多普勒項(xiàng)的P+1階延遲共軛相乘函數(shù)可表示為

(14)

(15)

(16)

根據(jù)文獻(xiàn)[10]的研究表明,彈道目標(biāo)在飛行中段只受地球引力的作用,彈道較為平緩,平動(dòng)多普勒近似為三階多項(xiàng)式足夠滿足精度要求,即四階延遲共軛相乘處理足以滿足對(duì)平動(dòng)多普勒補(bǔ)償?shù)男枨蟆?/p>

3微多普勒信號(hào)參數(shù)估計(jì)

微多普勒信號(hào)參數(shù)估計(jì)可采用與基函數(shù)匹配相關(guān)處理,搜索譜峰的方法來完成[11]44。對(duì)平動(dòng)補(bǔ)償后的信號(hào)采用門限檢測(cè)的方法估計(jì)信號(hào)帶寬,利用卡森準(zhǔn)則換算得調(diào)頻指數(shù),設(shè)置2項(xiàng)相關(guān)搜索參數(shù)分別為:旋轉(zhuǎn)頻率與微多普勒相位,取信號(hào)實(shí)部與基函數(shù)進(jìn)行匹配相關(guān)處理并搜索譜峰,通過換算可得所估參數(shù)。

根據(jù)上一節(jié)所得平動(dòng)補(bǔ)償后的信號(hào)可知,現(xiàn)需要估計(jì)的參數(shù)有3項(xiàng),分別為調(diào)頻指數(shù)mfs、自旋頻率fm與微多普勒初始相位φms。由于對(duì)3項(xiàng)參數(shù)同時(shí)搜索,迭代次數(shù)多,搜索譜峰算法復(fù)雜。因此需要對(duì)調(diào)頻指數(shù)進(jìn)行換算,由卡森(Carson)公式[12]可換算得調(diào)頻指數(shù)為

(17)

式中,B為調(diào)頻信號(hào)帶寬。采用帶寬采樣門限檢測(cè)的方法可較為方便準(zhǔn)確地估計(jì)出信號(hào)帶寬,并將搜索維度降低為二維,很大程度上減少了相關(guān)運(yùn)算的次數(shù),降低了譜峰搜索的難度。

關(guān)于門限檢測(cè)法估計(jì)信號(hào)帶寬的算法流程如下[7]354:

1) 對(duì)待估信號(hào)進(jìn)行FFT運(yùn)算計(jì)算其頻譜,并對(duì)頻譜做歸一化處理。

2) 尋找最大譜峰與次大譜峰,記其譜峰歸一化幅值分別為E1與E2,計(jì)算能量比η=E2/E1。根據(jù)能量比η判斷此調(diào)頻信號(hào)為弱調(diào)制還是深度調(diào)制(根據(jù)經(jīng)驗(yàn)值可知當(dāng)η<1時(shí),信號(hào)為弱正弦調(diào)制,否則為深度正弦調(diào)制)。

3) 當(dāng)判斷信號(hào)為深度正弦調(diào)制后,設(shè)置檢測(cè)門限μ,統(tǒng)計(jì)頻譜歸一化幅度高于檢測(cè)門限的頻率,取極值頻率之差即為估計(jì)所得帶寬B,通過多次測(cè)量剔除大于平均值的野值點(diǎn)。由于白噪聲對(duì)信號(hào)頻譜的影響常表現(xiàn)為邊帶頻率,因此檢測(cè)門限μ應(yīng)隨著信號(hào)載噪比環(huán)境的優(yōu)劣性而隨之改變。

結(jié)合式(17)構(gòu)建搜索基函數(shù)

(18)

(19)

(20)

(21)

其在有限時(shí)間T內(nèi)的數(shù)字信號(hào)形式可表達(dá)為

(22)

由于在做延遲共軛相乘處理時(shí),改變了微多普勒信號(hào)的調(diào)頻指數(shù)mf與初相φm,因此在構(gòu)建微多普勒補(bǔ)償函數(shù)前,需對(duì)搜索到的參數(shù)進(jìn)行縮放與移相處理。

(23)

由于在P階延遲共軛相乘處理中對(duì)信號(hào)進(jìn)行過無用信息的截除,截除處理改變了理論相位的初始位置,因此在計(jì)算實(shí)際微多普勒初始相位時(shí)應(yīng)做相位補(bǔ)償。截除均為從0時(shí)刻開始,分P次共截取Pτ(單位為s)的無用信息,此段時(shí)長內(nèi)的信號(hào)相位變化量為

(24)

(25)

4仿真試驗(yàn)與分析

4.1平動(dòng)補(bǔ)償

設(shè)置信號(hào)時(shí)長為4 s,目標(biāo)自旋頻率0.5 Hz,微多普勒頻偏40 Hz,微多普勒初始相位為π/3,根據(jù)標(biāo)準(zhǔn)測(cè)控信號(hào)取信噪比為20dB[13],一階多普勒變化率500Hz,通過式(4)可得多普勒理論時(shí)頻特征。

對(duì)信號(hào)S(t)做短時(shí)傅里葉變換(Short-TimeFourierTransform,STFT),可得基頻信號(hào)時(shí)頻特征,如圖1所示。

可看出STFT變換結(jié)果與理論時(shí)頻特征是一致的,微多普勒頻率表現(xiàn)為多普勒頻率上有規(guī)律的波動(dòng)。取延遲量τ=0.5s,進(jìn)行二階延遲共軛相乘處理,并對(duì)處理后的無用信息做截除,其補(bǔ)償結(jié)果如圖2所示。

從圖2a)中可以看出微多普勒頻率表現(xiàn)為以固定頻率為中心的周期性波動(dòng),說明一階延遲共軛相乘處理已補(bǔ)償多普勒頻率的一階分量。

從圖2b)中可看出二階延遲共軛相乘處理已將信號(hào)中的固定頻率予以補(bǔ)償,微多普勒頻率已被搬移至零頻。仿真結(jié)果與前文分析結(jié)果是一致的,滿足微動(dòng)參數(shù)估計(jì)的條件。

a) 理論時(shí)頻特征

b) STFT變換結(jié)果圖1 基頻信號(hào)多普勒時(shí)頻特征

a) 一階延遲共軛相乘處理

b) 二階延遲共軛相乘處理圖2 平動(dòng)補(bǔ)償信號(hào)時(shí)頻特征

4.2微多普勒參數(shù)搜索

根據(jù)上文分析通過對(duì)調(diào)頻指數(shù)的換算可減少參數(shù)搜索的維度以及運(yùn)算復(fù)雜度,因此首先對(duì)平動(dòng)補(bǔ)償后的信號(hào)帶寬進(jìn)行估計(jì),其帶寬如圖3所示。之后設(shè)置自旋頻率搜索范圍為0~2Hz,微多普勒相位搜索范圍為0~2π,整體搜索步長為1/100,結(jié)合式(22)進(jìn)行參數(shù)搜索,結(jié)果如圖 4所示。

圖3 平動(dòng)補(bǔ)償后的信號(hào)頻譜

圖4 微多普勒信號(hào)參數(shù)搜索結(jié)果

通過頻譜門限過濾與譜峰搜索,結(jié)合式(20)、式(23)與式(25)可計(jì)算對(duì)應(yīng)的微多普勒頻偏、自旋頻率與微多普勒相位。其結(jié)果與真實(shí)參數(shù)與未進(jìn)行平動(dòng)補(bǔ)償?shù)墓烙?jì)結(jié)果進(jìn)行比較,如表 1所示。

表1 微多普勒參數(shù)估計(jì)結(jié)果比較

可看出對(duì)未進(jìn)行平動(dòng)補(bǔ)償?shù)男盘?hào)進(jìn)行參數(shù)估計(jì),其估計(jì)值與真值出現(xiàn)了較大的偏差,而本方法通過延遲共軛相乘處理,完成了平動(dòng)多普勒的補(bǔ)償,所估參數(shù)精度得到了較大改善。

本算法結(jié)合卡森準(zhǔn)則進(jìn)行微多普勒信號(hào)帶寬估計(jì),根據(jù)式(20)可得本文方法的乘法運(yùn)算復(fù)雜度Ca為

(26)

式中,M為單參數(shù)循環(huán)次數(shù);K為復(fù)雜度系數(shù);N為信號(hào)采樣點(diǎn)數(shù);KN為單次匹配相關(guān)乘法運(yùn)算次數(shù);NF為FFT點(diǎn)數(shù);相比文獻(xiàn)[11]51中對(duì)微多普勒參數(shù)的三維搜索,本算法運(yùn)算復(fù)雜度減少M(fèi)KN-NFlbNF,其中NF<

為驗(yàn)證不同搜索步長下的參數(shù)估計(jì)效果,分別設(shè)置搜索步長為1/25,1/50與1/400進(jìn)行參數(shù)搜索,可得不同搜索步長下的估計(jì)精度,如表 2所示。

表2 不同搜索步長下的參數(shù)搜索時(shí)長與

可看出當(dāng)搜索步長較大時(shí),估計(jì)誤差較大;反之步長較小時(shí),估計(jì)誤差隨之減小。雖然短步長條件下估計(jì)誤差小,所估參數(shù)可取較高的精度,但搜索時(shí)間相比其他搜索補(bǔ)償較大,因此在實(shí)際應(yīng)用中通過結(jié)合測(cè)量精度與運(yùn)算時(shí)間的需求來權(quán)衡步長的選擇。

5結(jié) 束 語

針對(duì)平動(dòng)對(duì)旋轉(zhuǎn)目標(biāo)微多普勒參數(shù)估計(jì)帶來的頻率偏移問題,本文基于連續(xù)波雷達(dá)體制,通過延遲共軛相乘的方法完成了對(duì)信號(hào)平動(dòng)分量的補(bǔ)償,對(duì)微多普勒信號(hào)采用相關(guān)匹配的方法實(shí)現(xiàn)了微多普勒參數(shù)的估計(jì),并給出了仿真試驗(yàn)的具體步驟。仿真結(jié)果表明:該方法搜索維度低,算法原理簡單且易于實(shí)現(xiàn),在搜索步長為1/400時(shí),各項(xiàng)參數(shù)估計(jì)值的相對(duì)誤差均能控制在1%以下,適用于動(dòng)目標(biāo)的微動(dòng)參數(shù)估計(jì)。但在較低信噪比環(huán)境或其他干擾環(huán)境下(如常見的多徑干擾),如何完成平動(dòng)多普勒的補(bǔ)償與微動(dòng)參數(shù)的估計(jì),仍需結(jié)合實(shí)際情況對(duì)搜索步長的選擇進(jìn)行更好的權(quán)衡,進(jìn)一步完善估計(jì)算法。

參考文獻(xiàn)(References)

[1]CHENVC,LIF,HOSS,etal.Analysisofmicro-Dopplersignatures[J].IEEEProc.RadarSonarNavig,2003,150(4):271-276.

[2]屈文星,楊文革,張若禹.應(yīng)答式旋轉(zhuǎn)目標(biāo)的微多普勒研究[J].電子測(cè)量技術(shù),2015,38(11):32-36.

[3]CHENVC,LIF,HOSS,etal.Micro-Dopplereffectinradar:phenomenon,model,andsimulationstudy[J].IEEETransactionsonAerospaceandElectronicSystems,2006,42(1):2-20.

[4]高紅衛(wèi),謝良貴,文樹梁,等.速度對(duì)微多普勒的影響及其補(bǔ)償研究[J].航天電子對(duì)抗,2008(4):46-50.

[5]楊有春,童寧寧,馮存前,等.彈道中段目標(biāo)回波平動(dòng)補(bǔ)償與微多普勒提取[J].中國科學(xué):信息科學(xué),2013(9):1172-1182.

[6]李寶柱,袁起,何佩琨,等.目標(biāo)自旋引起微多普勒的補(bǔ)償新方法[J].現(xiàn)代雷達(dá),2008(10): 49-51.

[7]熊輝,呂遠(yuǎn),曾德國,等.利用卡森準(zhǔn)則的正弦調(diào)頻信號(hào)參數(shù)估計(jì)方法[J].電子測(cè)量與儀器學(xué)報(bào),2010,24(4):353-358.

[8]劉進(jìn).微動(dòng)目標(biāo)雷達(dá)信號(hào)參數(shù)估計(jì)與物理特征提取[D].長沙:國防科學(xué)技術(shù)大學(xué),2010:19-28.

[9]ZHAOY,WANGG,BAIY.AnewmethodoftranslationalcompensationforspacespinningtargetbasedonEMD[C]// 2013IEEEInternationalConferenceonSignalProcessingCommunicationandComputing.Kunming:IEEE,2013:1-5.

[10]陳行勇.微動(dòng)目標(biāo)雷達(dá)特征提取技術(shù)研究[D].長沙:國防科學(xué)技術(shù)大學(xué),2006:131-132.

[11]魏萌.微多普勒效應(yīng)分析及參數(shù)提取[D].成都:電子科技大學(xué),2010:44-51.

[12]樊昌信,曹麗娜.通信原理[M].北京:國防工業(yè)出版社,2008:109-110.

[13]李偉,楊文革,趙江.基于VFD濾波器的測(cè)控信號(hào)動(dòng)態(tài)信息加載方法[J].現(xiàn)代防御技術(shù),2015,43(3):146 -150.

(編輯:李江濤)

Research on Translational Motion Compensation and Micro-Doppler Parameter Estimation of Spinning Target

YANG Wenge1,QU Wenxing2,ZHANG Ruoyu3

(1. Department of Optical and Electronic Equipment, Equipment Academy, Beijing 101416, China;2. Department of Graduate Management, Equipment Academy, Beijing 101416, China;3. Beijing Institute of Tracking and Telecommunication Technology, Beijing 100094, China)

AbstractThe spinning motion of flying missile will produce micro-Doppler modulation on continuous wave radar echo signal. Parameter estimation for this modulated signals will obtain the spinning frequency of the target and the Doppler frequency of translational motion will have influence on the estimation of the micro parameters. Aiming at the problem of deviation of estimation on the micro-Doppler frequency parameter, through building up a carrier signal model with translational motion and micro-motion items , the paper obtains the information of the micro-Doppler by compensating translational motion with delayed conjugated multiplication, then estimates the micro motion parameters by matching relevant algorithm, and finally concludes the specific steps of this algorithm. The simulation result verifies the effectiveness of this method.

Keywordsmicro-Doppler; continuous wave radar; translational motion compensation; delayed conjugated multiplication; parameters estimation

收稿日期2015-12-01

作者簡介楊文革(1966-),男,教授,博士生導(dǎo)師,主要研究方向?yàn)楹教鞙y(cè)控。

中圖分類號(hào)TN959.6

文章編號(hào)2095-3828(2016)03-0107-06

文獻(xiàn)標(biāo)志碼A

DOI10.3783/j.issn.2095-3828.2016.03.021

主站蜘蛛池模板: 国产欧美日韩视频一区二区三区| 无码日韩精品91超碰| 国产精品永久在线| 久久午夜夜伦鲁鲁片不卡| 四虎成人精品| 精品在线免费播放| 天堂成人在线| 久久婷婷色综合老司机| 日本草草视频在线观看| 无码一区二区波多野结衣播放搜索| 久久精品66| 色综合成人| 中文字幕在线日本| 在线观看免费人成视频色快速| 天天视频在线91频| 亚洲人成日本在线观看| 午夜色综合| 毛片久久网站小视频| 一区二区午夜| 国产白浆一区二区三区视频在线| 日韩麻豆小视频| 欧美在线综合视频| 久久综合五月婷婷| 亚洲性影院| 久久天天躁狠狠躁夜夜2020一| 亚洲国产精品成人久久综合影院| 欧洲日本亚洲中文字幕| 2022国产无码在线| 日本一本正道综合久久dvd | 婷婷色一二三区波多野衣| 福利在线不卡| 性视频久久| 亚洲 成人国产| 亚洲水蜜桃久久综合网站| 免费在线看黄网址| 国产激情无码一区二区免费| 一区二区偷拍美女撒尿视频| 毛片免费观看视频| 欧美性猛交一区二区三区| 美女被狂躁www在线观看| 精品无码人妻一区二区| 色噜噜综合网| 精品国产香蕉在线播出| 日韩精品无码免费专网站| 成人免费黄色小视频| 黄色污网站在线观看| 青青草欧美| 无码 在线 在线| 国产尤物在线播放| 91视频青青草| 欧美国产日韩在线| 在线综合亚洲欧美网站| 国产精品美女自慰喷水| 国产福利在线观看精品| 久久这里只有精品66| 黄色国产在线| 91精品国产自产在线老师啪l| 2021国产在线视频| 亚洲综合一区国产精品| 亚洲午夜综合网| 亚洲第一极品精品无码| 久久动漫精品| 91热爆在线| 全部免费毛片免费播放| 久久亚洲高清国产| 亚洲欧洲日韩综合| 日韩免费无码人妻系列| 欧美成人怡春院在线激情| 国产精品女人呻吟在线观看| 高清免费毛片| 亚洲国产看片基地久久1024| 国产毛片片精品天天看视频| 亚洲国产看片基地久久1024| 欧美日韩久久综合| 亚洲VA中文字幕| 国产精品大尺度尺度视频| 好紧好深好大乳无码中文字幕| 中文字幕1区2区| 国产成人精品优优av| 亚洲中文在线看视频一区| 欧美三级不卡在线观看视频| 日韩欧美中文|