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

基于Radon-Ambiguity變換的LFM信號時/頻差快速聯(lián)合估計

2016-03-07 08:52:23楊林森張子敬郭付陽
電波科學學報 2016年6期
關(guān)鍵詞:信號

楊林森 張子敬 郭付陽

(西安電子科技大學 雷達信號處理國家重點實驗室,西安 710071)

基于Radon-Ambiguity變換的LFM信號時/頻差快速聯(lián)合估計

楊林森 張子敬 郭付陽

(西安電子科技大學 雷達信號處理國家重點實驗室,西安 710071)

提出了一種基于Radon-Ambiguity變換(Radon-Ambiguity Transform, RAT)的線性調(diào)頻(Linear Frequency Modulated, LFM)信號時/頻差快速聯(lián)合估計的算法。根據(jù)LFM信號在多個不同角度上的RAT峰值位置建立一組以信號間時差和頻差為未知量的方程組,求解方程組即可得到時/頻差的估計值。對于存在噪聲的信號,RAT誤差會導致方程組不能直接求解,為了抑制噪聲干擾,采用最小二乘法估計時/頻差。本文算法無需計算二維平面上各點的模糊函數(shù)值,并且由于離散RAT可以通過快速傅里葉變換快速實現(xiàn),具有所需運算量低的優(yōu)點。仿真實驗表明,相比于常見的基于模糊函數(shù)峰值搜索的時/頻差估計算法,本文算法在保證時/頻差估計精度的同時能夠顯著地提高運算效率。

Radon-Ambiguity變換;最小二乘法;時/頻差;線性調(diào)頻信號

引 言

對于兩路信號間時差和頻差的估計問題廣泛地存在于雷達,聲納和導航的很多應用之中[1-4].在無源定位技術(shù)中,常常根據(jù)兩路接收信號間的時差和頻差對目標輻射源進行精確定位.雙星時/頻差定位是一種最常見的無源定位技術(shù)的應用,它利用兩顆衛(wèi)星接收同一地面輻射源的輻射信號,由于兩顆衛(wèi)星相對于該地面輻射源具有不同的距離和徑向速度,因此衛(wèi)星接收到的兩路信號間存在時差和多普勒頻差.由兩路接收信號間的時差和頻差可以確定出目標輻射源在空間中可能出現(xiàn)的區(qū)域范圍,再結(jié)合地球表面方程,即可最終實現(xiàn)地面輻射源的精確定位,顯然時/頻差估計技術(shù)直接影響到最終的定位性能.模糊函數(shù)是估計信號間時/頻差最常使用的工具[5],模糊函數(shù)是一種在時間和頻率二維平面上定義的廣義相關(guān).若兩路信號間僅存在時差,利用時域相關(guān)即可有效地估計該時差,然而當兩路信號間同時存在時差和頻差時,時域相關(guān)顯然無法滿足時差和頻差的同時估計,模糊函數(shù)作為時頻二維平面上的廣義相關(guān),可以實現(xiàn)兩路信號間時差和頻差的同時估計.類似于由時域相關(guān)峰值估計時差的過程,利用兩路信號間互模糊函數(shù)峰值位置也可以實現(xiàn)信號間時差和頻差的估計.然而,計算二維模糊平面上各點的模糊函數(shù)值并搜索峰值意味著巨大的運算量,信號間的時/頻差僅僅與模糊函數(shù)峰值點位置有關(guān),為了尋找峰值而計算整個模糊平面上各點的函數(shù)值顯然存在巨大的計算冗余.為此,研究人員提出了很多提高模糊函數(shù)峰值搜索效率的改進算法.文獻[6]提出了先粗估計模糊函數(shù)峰值位置再精估計的思路,減少了很多不必要的計算.文獻[7]根據(jù)多普勒頻移一般遠小于采樣頻率的特點,提出對兩路信號的混合積信號先進行低通濾波再抽取處理過程,顯著地減小了所處理信號的長度,從而降低了運算量.對于線性調(diào)頻(Linear Frequency Modulated, LFM)信號,文獻[8]根據(jù)其模糊函數(shù)圖存在一條過峰值點、斜率為調(diào)頻率的脊線的特點,將模糊函數(shù)二維峰值搜索轉(zhuǎn)化為沿脊線的一維峰值搜索,從而顯著地降低了時/頻差估計所需的運算量.

分數(shù)階傅里葉變換(fractional Fourier transform, FrFT)作為傅里葉變換的推廣,廣泛地應用于信號的檢測以及參數(shù)估計等問題中[9-12].文獻[13]中給出了離散FrFT基于快速傅里葉變換(Fast Fourier Transform, FFT)的快速實現(xiàn)方法,運算復雜度與FFT算法相當.Radon變換是一種沿直線積分的投影變換[14],Radon-Ambiguity變換(Radon-Ambiguity Transform, RAT)將信號的模糊函數(shù)在模糊平面上進行Radon變換.文獻[15]結(jié)合FrFT和模糊函數(shù)的性質(zhì),將RAT的定義式化簡為基于若干次FrFT和傅里葉變換的形式,因此離散RAT可以通過若干次FFT快速實現(xiàn),顯著地降低了RAT所需的運算量.

本文提出了一種基于RAT的LFM信號時/頻差快速聯(lián)合估計算法.首先根據(jù)兩路LFM信號在不同角度上RAT的峰值位置建立一組以時差和頻差為未知量的方程組,求解方程組即可獲得時差和頻差估計值.對于受到噪聲干擾的信號,RAT誤差會導致方程組不能直接求解.為了抑制噪聲的影響,采用最小二乘法求解時/頻差.最后,通過理論分析和仿真實驗,分別對比了本文算法和文獻[8]提出的Radon-Ambiguity單切片(single slice of Radon-Ambiguity transform,SSRAT)時/頻差估計算法所需的運算時間和估計精度.

1 模糊函數(shù)及Radon-Ambiguity變換

假設(shè)在噪聲環(huán)境中存在時差和頻差的兩路LFM信號為

(1)

對于兩路信號r1(t)和r2(t),其互模糊函數(shù)定義為[5]

(2)

若不考慮噪聲影響,將式(1)代入式(2),可以得到LFM信號的互模糊函數(shù)

(3)

由式(3)可以看出,LFM信號互模糊函數(shù)存在一個峰值點(Δτ,Δf),并且在二維模糊平面上還存在一條過峰值點的脊線,脊線所對應的方程為

f=mτ+Δf-mΔτ.

(4)

RAT變換是對信號的模糊函數(shù)χ(τ,f)在二維模糊平面上做Radon投影變換[14]的結(jié)果.圖1為RAT變換的示意圖,圖中τ軸和f軸分別為模糊平面的時延軸和多普勒頻率軸,將τ軸逆時針旋轉(zhuǎn)角度φ得到u軸,v軸垂直于u軸.信號在角度φ上的RAT定義為模糊函數(shù)χ(τ,f)沿著平行于v軸并且與u軸交于不同位置的直線PQ的積分變換,即

RAT(u,φ)=∫PQχ(τ,f)dv.

(5)

圖1 RAT變換示意圖

若根據(jù)式(5)計算RAT,則需先計算模糊函數(shù)在二維平面上各點的函數(shù)值,然后再對模糊函數(shù)沿不同的積分路徑進行積分,需要巨大的運算量,不利于快速實現(xiàn).文章[15]結(jié)合FrFT和模糊函數(shù)的性質(zhì),將信號的RAT化簡為基于FrFT和傅里葉變換的等價形式:

(6)

式中:Rφ(r)為信號r(t)在角度φ上的FrFT.由于離散FrFT能夠通過FFT快速實現(xiàn)[13],因此,根據(jù)式(6),離散RAT可以通過若干次FFT快速實現(xiàn),能夠顯著地提高RAT的計算效率.

2 利用RAT估計LFM信號時/頻差

(7)

根據(jù)式(4),LFM信號的模糊函數(shù)存在一條過峰值點且斜率為調(diào)頻率m的脊線.模糊函數(shù)上的主要能量都集中在脊線上,當φ=arctanm時,RAT的積分路徑與LFM信號模糊函數(shù)的脊線平行,脊線上的所有點將投影聚焦于一點,產(chǎn)生一個極大的峰值,因此信號在角度φ=arctanm的RAT峰值具有很強的抗噪聲性能,建立的方程式(7)具有較高的精度.然而,當φ≠arctanm時,RAT峰值點的產(chǎn)生過程沒有積累的作用,RAT峰值點的位置估計易受到噪聲的干擾,建立的方程式(7)易受噪聲干擾而產(chǎn)生誤差.因此,由兩個不同角度的RAT峰值位置建立的方程組中至少存在一個易受噪聲干擾方程,最終求解的時/頻差估計值無法保證精度.為了抑制噪聲對時/頻差估計精度的影響,本文將選取多個不同角度進行RAT,利用多次觀測來降低單次觀測的誤差對時/頻差估計值的影響.

(8)

Ax=b.

(9)

(10)

將LFM信號在各角度的RAT峰值位置觀測向量b代入式(10)即可得到時/頻差的最小二乘估計.

本文提出的基于RAT的LFM信號時/頻差快速聯(lián)合估計算法的步驟可以總結(jié)如下:

1) 在區(qū)間[0,π)均勻選取M個角度,記為φi(i=1,…,M).提前計算式(10)中的矩陣(ATA)-1AT.

4) 對于步驟3)中得到的乘積信號,計算其逆FFT,結(jié)果即為兩路信號在角度φi上的RAT.

6) 利用式(10)直接計算兩路LFM信號間時/頻差的最小二乘解.

3 算法運算量分析

模糊函數(shù)是估計信號間時/頻差最常用到的工具,為了提高模糊函數(shù)峰值搜索的運算效率,文獻[8]提出了SSRAT算法.SSRAT算法利用LFM信號模糊函數(shù)在二維平面上存在一條過峰值脊線的特點,首先估計出LFM模糊函數(shù)脊線的位置,然后沿著脊線搜索模糊函數(shù)峰值,由文獻[8]的結(jié)論可知SSRAT算法能夠顯著降低時/頻差估計所需的運算量.下面我們將對比本文算法和SSRAT算法所需的運算量.

假設(shè)離散采樣信號的長度為N,本文算法所需要的RAT角度個數(shù)為M.對于本文算法,由式(6)計算任意角度的RAT僅需計算兩次FrFT和三次傅里葉變換,計算離散采樣信號的FrFT和傅里葉變換都能夠利用FFT快速實現(xiàn),分別需要復乘次數(shù)約為3Nlog2N+N[13]和Nlog2N,則計算信號在M個角度的RAT所需的總復乘次數(shù)約為9MNlog2N+6MN.根據(jù)信號在M個角度上RAT的峰值位置,利用式(10)可直接求解時/頻差的最小二乘解.式(10)中矩陣A的取值只與RAT角度有關(guān),矩陣A可以認為是已知的,因此矩陣(ATA)-1AT也可以預先計算,在此不考慮其所需的運算量.將向量b代入式(10)計算信號時/頻差的最小二乘解所需的復乘次數(shù)為2M.綜上分析,本文算法的運算量約為9MNlog2N+6MN+2M.

考慮到實際問題中的離散采樣信號長度N一般都會較大,而本文算法所需的RAT角度個數(shù)M一般小于102,因此M往往遠遠小于信號的長度N,即M?N,則本文算法的運算量可以認為是O(Nlog2N),而SSRAT算法的運算量為O(N2)[8].因此相比于SSRAT算法,本文算法的運算量更小,運算效率更高.

4 仿真實驗

下面我們將通過仿真驗證本文算法的性能,假設(shè)兩路LFM信號帶寬為50 kHz,調(diào)頻率31.25 MHz/s,初始頻率κ=0,采樣頻率為200 kHz,信號持續(xù)時間為1.6 ms,觀測時間為2 ms,兩路信號間的時差和頻差分別為-10 μs和-3.2 kHz.

4.1 RAT角度搜索數(shù)量對本文算法性能的影響

圖2(a)仿真了信號在不同角度上RAT峰值位置與該信號模糊函數(shù)峰值點在對應角度的u軸上投影位置的關(guān)系.為了避免不同的量綱對實驗結(jié)果的影響,在仿真中對離散采樣信號的時間和頻率信息進行了去量綱處理[13].圖2(a)中的實線為該信號在[0,π)上各角度RAT峰值位置的仿真曲線,虛線為該信號模糊函數(shù)峰值點在各對應角度u軸上投影位置理論曲線.對于離散采樣信號,其RAT也是離散的,若選取的兩個RAT角度足夠接近,會使得不同角度的RAT峰值位置落入同一個離散分辨單元而無法區(qū)分,因此圖2(a)中的RAT峰值位置變化仿真曲線是一條存在階躍的曲線.可以看出,圖2(a)中的實線和虛線位置相互吻合,驗證了LFM信號在任意角度上RAT的峰值位置等于其互模糊函數(shù)峰值點在該角度所對應u軸上的投影位置.

圖2(b)和(c)分別仿真了在輸入信噪比為10 dB時,本文算法時差和頻差估計均方根誤差隨RAT角度數(shù)量的變化曲線.由圖2(b)和(c)中的曲線可以看出:當RAT角度數(shù)量小于30時,本文算法時/頻差估計均方根誤差隨RAT角度數(shù)量增加而迅速減小;當RAT角度數(shù)量大于60時,本文算法的時/頻差估計均方根誤差趨于穩(wěn)定,基本不再隨RAT角度數(shù)量的增加而降低.這是因為當RAT角度數(shù)量較小時,RAT角度數(shù)量的增加使得觀測樣本增加,有利于抑制單次觀測誤差對時/頻差估計精度的影響,然而當RAT角度數(shù)量已經(jīng)較大時,再增加RAT角度數(shù)量會導致相鄰角度間的差距過小,由前面的分析可知,對于差距較小的兩個角度,根據(jù)式(7)得到的兩個方程會接近線性相關(guān),進而求解的時/頻差估計值易受RAT峰值位置估計誤差的影響,因此當RAT角度數(shù)量增加到一定程度時,本文算法時/頻差估計精度不會再隨RAT角度數(shù)量增加而提高.考慮到RAT角度數(shù)量會影響本文算法所需運算量,為了平衡估計精度和所需的運算量,本文算法RAT角度數(shù)應在30~60選擇.

4.2 本文算法和幾種常見算法的性能比較

圖3仿真了本文算法和SSRAT算法時/頻差估計均方根誤差隨輸入信噪比的變化曲線.為了得到不同算法時差和頻差估計的均方根誤差,對于圖3中的每一個信噪比,我們分別進行了500次蒙特卡洛仿真,并且和克拉美-羅界[6]進行了比較.由圖3(a)和(b)可以看出,兩種不同算法的時/頻差估計均方根誤差曲線基本重合,并且隨信噪比的提高兩種不同算法的時/頻差估計均方根誤差曲線都不斷減小且接近克拉美-羅界,兩種算法都能夠得到較高的時/頻差估計精度.

圖3 不同算法的時/頻差估計均方誤根差隨信噪比的變化曲線

圖4仿真了本文算法和SSRAT算法所需的運算時間隨信號長度變化的仿真曲線.在該仿真中,信號長度在變化,其余參數(shù)與之前相同.根據(jù)4.1節(jié)中的仿真結(jié)果可知,為了同時滿足較高的時/頻差估計精度和較小的運算量,本文算法所需的RAT角度數(shù)量應在區(qū)間選擇,因此圖4針對本文算法的仿真僅選取了RAT角度數(shù)量為30和60的兩種情況.對比圖4中的3條不同曲線可以看出:在相同的信號參數(shù)條件下,本文算法在RAT角度數(shù)量為30和60時所需的運算時間都明顯小于SSRAT算法所需要的運算時間;另外,本文算法在RAT角度數(shù)為30時所需的運算時間略小于RAT角度數(shù)為60時所需的運算時間.

結(jié)合圖3和圖4的仿真結(jié)果可知,本文算法和SSRAT算法的時/頻差估計精度基本相同,都能實現(xiàn)時/頻差的精確估計,但是本文算法能夠顯著降低運算時間,運算效率具有明顯的優(yōu)勢.

圖4 不同算法時/頻差估計運算時間比較

5 結(jié) 論

對于兩路存在時差和頻差的LFM信號,本文提出了一種基于RAT的時/頻差快速聯(lián)合估計算法.根據(jù)LFM信號在任意一個角度上RAT的峰值位置能夠得到一個以信號間時差和頻差為未知量的二元方程,選取多個不同的角度能夠建立一組二元方程組,求解方程組即可得到時/頻差的估計值.對于存在噪聲干擾的信號,其RAT誤差會導致得到的方程組不能直接求解.為了抑制噪聲的影響,本文算法采用最小二乘法實現(xiàn)時/頻差的估計.本文算法無需計算二維平面上各點的模糊函數(shù)值,并且存在基于FFT的快速實現(xiàn)方法,由理論分析可知,本文算法所需的運算量遠小于SSRAT算法,同時實驗仿真驗證了,在相同條件下,本文算法和SSRAT算法時/頻差估計精度相同,都能夠獲得較高的時/頻差估計精度,然而本文算法所需的運算時間明顯小于SSRAT算法.由于時/頻差估計技術(shù)是無源定位的關(guān)鍵技術(shù),本文提出的時/頻差估計算法同時具有估計精度高和運算復雜度低的特點,具有一定的實際應用價值,特別適合于對時/頻差估計的速度和精度都有要求的實時處理問題,例如雙星利用時/頻差對雷達或通信基站的定位問題.

[1] KIM Y H, KIM D G, KIM H N. Two-step estimator for moving-emitter geolocation using time difference of arrival/frequency-difference of arrival measurements[J]. IET Radar, Sonar & Navigation, 2015, 9(7): 881-887.[2] 李昕. 基于雙通道DFRFT互譜法的Chirp信號時延估計[J]. 電子學報, 2014, 42(6): 1068-1073.

LI X. Time delay estimation of Chirp signals based on doubled-channel DFRFT cross-spectrum[J]. Acta Electronica Sinica, 2014, 42(6): 1068-1073. (in Chinese)

[3] MUSICKI D, KAUNE R, KOCH W. Mobile emitter geolocation and tracking using TDOA and FDOA measurements[J]. IEEE Transactions on Signal Processing, 2010, 58(3): 1863-1874.

[4] 朱文濤, 蘇濤, 楊濤, 等. 低信噪比下線性調(diào)頻連續(xù)波信號的參數(shù)估計[J]. 電波科學學報, 2013, 28(6): 1158-1164.

ZHU W T, SU T and YANG T, et al. Parameter estimation of linear frequency modulated continuous wave signal in low SNR[J]. Chinese Journal of Radio Science, 2013, 28(6): 1158-1164.

[5] STEIN S. Differential delay/Doppler ML estimation with unknown signals[J]. IEEE Transactions on Acoustic, Speech, and Signal Processing, 1993, 41(8): 2717-2719.

[6] STEIN S. Algorithms for Ambiguity function processing[J]. IEEE Transactions on Acoustic, Speech, and Signal Processing, 1981, 29(3): 588-599.

[7] TOLIMIERI R, WINOGRAD S. Computing the Ambiguity surface[J]. IEEE Transactions on Acoustic, Speech, and Signal Processing, 1985, 33(4): 1239-1245.

[8] SHARIF M R, ABEYSEKERA S S. Efficient wideband signal parameter estimation using a Radon-Ambiguity transform slice[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(2): 673-688.

[9] 劉鋒, 徐會法, 陶然, 等. 分數(shù)階Fourier域多分量LFM信號間的分辨研究[J]. 中國科學(信息科學), 2012, 42(2): 136-148.

LIU F, XU H F, TAO R, et al. Research on resolution among multi-component LFM signals in the fractional Fourier domain[J]. Science China (Information Science), 2012, 42(2): 136-148. (in Chinese)

[10] 戰(zhàn)立曉, 湯子躍, 朱振波. 基于分數(shù)階傅里葉變換的加速微弱目標檢測與估計[J]. 電波科學學報, 2013, 28(2): 296-304.

ZHAN L X, Tang Z Y, ZHU Z B. FRFT based detection and estimation of accelerating weak target[J]. Chinese Journal of Radio Science, 2013, 28(2): 296-304. (in Chinese)

[11] SEJDIC E, DJUROVIC I, STANKOVIC L. Fractional Fourier transform as a signal processing tool: An overview of recent developments[J]. Signal Processing, 2010, 91(6): 1351-1369.

[12] OONINCX P J. Joint time-frequency offset detection using the fractional Fourier transform[J]. Signal Processing, 2008, 88: 2936-2942.

[13] OZAKTAS H M, ARIKAN O, KUTAY M A, et al. Digital computation of the fractional Fourier transform[J]. IEEE Transactions on Signal Processing, 1996, 44(9): 2141-2150.

[14] WOOD J C, BARRY D T. Tomographic time-frequency analysis and its application toward time-varying filtering and adaptive kernel design for multicomponent linear-FM signals[J]. IEEE Transactions on Signal Processing, 1994, 42(8), pp: 2094-2104

[15] YANG L S, ZHANG Z J, GUO F Y. Fast algorithm for Radon-Ambiguity transform[J]. IET Radar, Sonar & Navigation, 2016, 10(3), pp: 553-559

[16] 張賢達. 矩陣分析與應用[M]. 北京:清華大學出版社, 2004: 71-96.

ZHANG X D. Matrix analysis and applications[M]. Beijing: Tsinghua University Press, 2004: 71-96 (in Chinese)

Joint time-frequency offset estimation for LFM signals based on the Radon-Ambiguity transform

YANG Linsen, ZHANG Zijing, GUO Fuyang

(NationalLaboratoryofRadarSignalProcessing,XidianUniversity,Xi’an710071,China)

A fast method for joint estimation of the time-frequency offset for linear frequency modulated signals based on the Radon-Ambiguity transform (RAT) is proposed in this paper. According to peak positions of the RAT of a LFM signal on different angles, a set of equations with the time-frequency offset as the unknowns can be established and the time-frequency offset can be estimated by solving the equations. For signals disturbed by noise, errors of the RAT will cause the equations having no solutions. In order to eliminate the noise, the least square method is used to estimate the time-frequency offset. Because the proposed algorithm does not need to calculate the value of each point on the 2-D Ambiguity plane and the RAT for sampled signals can be realized rapidly by several processes of the fast Fourier transform, the proposed method has advantage of low computational cost. Simulation results show that the proposed algorithm can ensure the accuracy of the estimation of the time-frequency offset as well as is computationally efficient compared with common methods based on peak searching of the Ambiguity function.

Radon-Ambiguity transform; Least square method; Time-frequency offset; LFM

10.13443/j.cjors.20160916

2016-09-16

國家自然科學基金(No.61571349, No.61201325);國家留學基金資助項目

TN911.72

A

1005-0388(2016)06-0-0

楊林森 (1988-),男,陜西人,博士研究生,研究方向為時頻分析、雷達信號處理.

張子敬 (1967-),男,陜西人,教授,博士生導師,研究方向為無源定位、雷達信號處理.

郭付陽 (1991-),男,江西人,博士研究生,研究方向為時頻分析、雷達信號處理.

楊林森, 張子敬, 郭付陽. 基于Radon-Ambiguity變換的LFM信號時/頻差快速聯(lián)合估計[J]. 電波科學學報,2016,31(6):頁碼.

YANG L S, ZHANG Z J, GUO F Y. Joint time-frequency offset estimation for LFM signals based on the Radon-Ambiguity transform[J]. Chinese journal of radio science,2016,31(6):頁碼.(in Chinese). DOI: 10.13443/j.cjors.20160916

聯(lián)系人: 楊林森 E-mail:yls198802@163.com

DOI 10.13443/j.cjors.20160916

猜你喜歡
信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
7個信號,警惕寶寶要感冒
媽媽寶寶(2019年10期)2019-10-26 02:45:34
孩子停止長個的信號
《鐵道通信信號》訂閱單
基于FPGA的多功能信號發(fā)生器的設(shè)計
電子制作(2018年11期)2018-08-04 03:25:42
基于Arduino的聯(lián)鎖信號控制接口研究
《鐵道通信信號》訂閱單
基于LabVIEW的力加載信號采集與PID控制
Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: 好吊日免费视频| 一级毛片免费不卡在线| 欧美高清三区| 国产精品免费p区| 亚洲床戏一区| 久久久久免费精品国产| 国产极品粉嫩小泬免费看| 国产成人免费高清AⅤ| 国产精品免费p区| 婷婷色丁香综合激情| 亚洲国产天堂久久九九九| 欧美伊人色综合久久天天| 中文字幕免费在线视频| 日韩专区第一页| a级毛片免费播放| 久久国产亚洲偷自| 色视频国产| 亚洲一级无毛片无码在线免费视频| 91精品久久久久久无码人妻| 99伊人精品| 在线观看免费黄色网址| 呦视频在线一区二区三区| 青青草原国产免费av观看| 国产一区在线视频观看| 亚洲人成网站色7777| 五月综合色婷婷| 中文字幕无码电影| 福利一区在线| 一级成人a毛片免费播放| 亚洲va在线观看| 精品国产乱码久久久久久一区二区| 免费99精品国产自在现线| 亚洲无码A视频在线| 亚洲国产欧美自拍| 国产一区在线观看无码| 国产主播在线一区| 国产成人综合亚洲网址| 婷婷综合缴情亚洲五月伊| 国产区在线观看视频| 国产视频久久久久| 亚洲va欧美ⅴa国产va影院| 亚洲免费人成影院| 午夜免费视频网站| 亚洲91在线精品| 亚洲大尺码专区影院| 国产麻豆精品久久一二三| 国产精品福利尤物youwu| 综合色88| 国产国产人成免费视频77777 | 一级毛片高清| 国产浮力第一页永久地址| 日韩在线1| 亚洲最大福利视频网| 国产欧美日韩免费| a级高清毛片| 超碰色了色| 午夜丁香婷婷| 麻豆国产精品一二三在线观看| 四虎国产精品永久一区| 成人免费一级片| 精品丝袜美腿国产一区| 久久综合九九亚洲一区| 久久久久亚洲AV成人人电影软件| 亚洲天堂成人| 国产丝袜91| 狠狠五月天中文字幕| 粉嫩国产白浆在线观看| 国产在线一区二区视频| 国产精品欧美亚洲韩国日本不卡| 欧美国产在线看| 亚洲精品在线91| 人人妻人人澡人人爽欧美一区| 亚洲第一黄色网址| 色婷婷狠狠干| 99久久精品久久久久久婷婷| 亚洲福利网址| 草草影院国产第一页| 综合色88| 色婷婷色丁香| 国产视频大全| 免费国产高清视频| 美女高潮全身流白浆福利区|