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

對(duì)稱穩(wěn)定分布噪聲下基于廣義相關(guān)熵的DOA估計(jì)新方法

2016-08-30 11:57:29邱天爽任福全李景春譚海峰大連理工大學(xué)電子信息與電氣工程學(xué)部大連116024國(guó)家無(wú)線電監(jiān)測(cè)中心北京100037北京郵電大學(xué)信息與通信工程學(xué)院北京100876
電子與信息學(xué)報(bào) 2016年8期
關(guān)鍵詞:定義方法

王 鵬 邱天爽* 任福全 李景春 譚海峰(大連理工大學(xué)電子信息與電氣工程學(xué)部大連116024)(國(guó)家無(wú)線電監(jiān)測(cè)中心北京100037)(北京郵電大學(xué)信息與通信工程學(xué)院北京100876)

?

對(duì)稱穩(wěn)定分布噪聲下基于廣義相關(guān)熵的DOA估計(jì)新方法

王鵬①邱天爽*①任福全①李景春②譚海峰②③①
①(大連理工大學(xué)電子信息與電氣工程學(xué)部大連116024)
②(國(guó)家無(wú)線電監(jiān)測(cè)中心北京100037)
③(北京郵電大學(xué)信息與通信工程學(xué)院北京100876)

針對(duì)穩(wěn)定隨機(jī)變量有限二階矩不存在的特點(diǎn),該文定義了一種新的廣義相關(guān)熵,并從理論上證明了對(duì)稱穩(wěn)定分布隨機(jī)變量廣義相關(guān)熵的有界性。此外,提出了一種穩(wěn)定分布噪聲下基于最小廣義相關(guān)熵準(zhǔn)則的DOA估計(jì)新方法,給出了一種迭代優(yōu)化算法并通過仿真實(shí)驗(yàn)分析了算法的收斂性。仿真結(jié)果表明,與現(xiàn)有基于分?jǐn)?shù)低階矩的FLOM-MUSIC、基于類相關(guān)熵的CRCO-MUSIC以及基于lp范數(shù)的ACO-MUSIC算法相比,所提方法可以獲得更好的估計(jì)結(jié)果,尤其是在高脈沖性噪聲環(huán)境下具有更加明顯的優(yōu)勢(shì)。

波達(dá)方向估計(jì);相關(guān)熵;廣義相關(guān)熵;穩(wěn)定分布噪聲;MUSIC算法

1 引言

DOA(Direction O f Arrival)估計(jì)是陣列信號(hào)處理中的基本問題之一,廣泛應(yīng)用于雷達(dá)、聲吶以及無(wú)線電通信等領(lǐng)域[1]。多重信號(hào)分類[2](MU ltiple SIgnal Classification,MUSIC)算法能夠?qū)崿F(xiàn)DOA的超分辨率估計(jì),但是傳統(tǒng)算法多假設(shè)背景噪聲服從高斯分布。實(shí)際上,由于受到自然因素(如大氣噪聲、海雜波等)以及人為因素(如電動(dòng)機(jī)等電磁設(shè)備)的影響,噪聲可能呈現(xiàn)較強(qiáng)的脈沖性,此時(shí)利用A lpha穩(wěn)定分布[3]進(jìn)行描述更加合適。與高斯分布隨機(jī)變量不同,A lpha穩(wěn)定分布隨機(jī)變量不具有有限二階矩,傳統(tǒng)MUSIC方法不再適用。

為抑制A lpha穩(wěn)定分布噪聲的影響,文獻(xiàn)[4-8]提出了基于分?jǐn)?shù)低階統(tǒng)計(jì)量(Fractional Lower Order Statistics,F(xiàn)LOS)的DOA估計(jì)方法。該類方法雖取得了較好的估計(jì)效果,但存在一定的局限性:(1)階次p必須滿足1≤p<α[4-6]或0

為了更好地抑制脈沖噪聲,提高DOA估計(jì)算法的魯棒性,本文定義了一種新的廣義相關(guān)熵,提出了一種基于最小廣義相關(guān)熵準(zhǔn)則的DOA估計(jì)新方法。該方法無(wú)需構(gòu)造協(xié)方差矩陣,通過優(yōu)化算法直接實(shí)現(xiàn)信號(hào)子空間的估計(jì),為利用相關(guān)熵實(shí)現(xiàn)DOA估計(jì)提供了一種新思路。仿真結(jié)果表明,本文算法在高脈沖性噪聲環(huán)境下能夠獲得更好的估計(jì)結(jié)果。

2 穩(wěn)定分布與相關(guān)熵

2.1穩(wěn)定分布

對(duì)稱A lpha穩(wěn)定(Symmetric A lpha Stable,SαS)分布隨機(jī)變量常用其特征函數(shù)進(jìn)行描述[3]。

其中,2α<≤0稱為特征指數(shù),用來(lái)度量概率密度函數(shù)的拖尾厚度,0γ>是分散系數(shù),用以度量樣本的分散程度,μ是位置參數(shù),當(dāng)1 2α<≤時(shí),該參數(shù)表示該隨機(jī)變量的均值,當(dāng)0 1α<<時(shí),該參數(shù)則表示該隨機(jī)變量的中值。特別地,當(dāng)2α=時(shí),SαS分布退化為高斯分布。當(dāng)2α<時(shí),S Sα隨機(jī)變量不存在有限的二階矩,因此基于二階統(tǒng)計(jì)量的傳統(tǒng)DOA估計(jì)算法在穩(wěn)定分布噪聲環(huán)境下性能嚴(yán)重退化。

2.2相關(guān)熵

相關(guān)熵作為一種隨機(jī)變量局部相似性的度量[13],廣泛應(yīng)用于時(shí)延估計(jì)[14]、自適應(yīng)濾波[15,16]、圖像處理[17]以及醫(yī)學(xué)信號(hào)處理[18]等領(lǐng)域,其定義為

3 基于廣義相關(guān)熵的DOA估計(jì)

3.1問題描述

假設(shè)均勻等距線陣由M個(gè)陣元組成,陣元間距為d,存在P個(gè)頻率已知的窄帶信源,信源之間相互獨(dú)立,入射角度分別為θ1,θ2,…,θP,以第1個(gè)陣元為參考陣元,則第m個(gè)陣元t時(shí)刻的輸出可以表示為

其中,si(t)為第i個(gè)信源的復(fù)包絡(luò),λ表示信號(hào)波長(zhǎng)且滿足d≤λ/2,nm(t)表示服從SαS分布的加性噪聲,各陣元噪聲相互獨(dú)立,噪聲與信號(hào)之間相互獨(dú)立。進(jìn)一步將式(5)寫成矩陣形式:

根據(jù)低秩近似理論[19],X可以利用低秩矩陣X進(jìn)行近似:

其全局最優(yōu)解如下:

其中,SVD(·)表示奇異值分解,Ds是由P個(gè)最大奇異值組成的對(duì)角矩陣,Us和Vs分別為對(duì)應(yīng)的左特征向量和右特征向量。由矩陣相關(guān)理論可知,矩陣Us與A張成相同的子空間。將式(8)的代價(jià)函數(shù)進(jìn)一步寫為

由于服從S Sα分布的隨機(jī)變量不具有有限的二階矩,基于式(11)算法的性能必然下降。

3.2廣義相關(guān)熵

為抑制穩(wěn)定分布噪聲的影響,受相關(guān)熵啟發(fā),本文定義了一種新的廣義相關(guān)熵:

可以證明,對(duì)稱穩(wěn)定分布隨機(jī)變量的廣義相關(guān)熵是有界的。

定理1假設(shè)存在兩個(gè)獨(dú)立同分布的對(duì)稱穩(wěn)定分布隨機(jī)變量X和Y,定義r X Y=-,則廣義相關(guān)熵是有界的。

證明根據(jù)廣義相關(guān)熵的定義可知

由S Sα分布的性質(zhì)[3]可知,r X Y=-同樣服從S Sα分布。根據(jù)數(shù)學(xué)期望的定義,式(15)可進(jìn)一步寫為

其中,()f r表示隨機(jī)變量的概率密度函數(shù)。

根據(jù)概率密度函數(shù)與特征函數(shù)之間的關(guān)系,將式(1)所示的特征函數(shù)代入,可得

又因?yàn)?/p>

所以,式(17)可以進(jìn)一步簡(jiǎn)化為

3.3基于最小廣義相關(guān)熵準(zhǔn)則的DOA估計(jì)算法

利用廣義相關(guān)熵替換式(11)中的二階矩,本文提出如式(21)所示的代價(jià)函數(shù)。

為了更好說明廣義相關(guān)熵對(duì)穩(wěn)定分布噪聲的抑制作用,將式(21)進(jìn)一步寫為

其中,B C⊙表示矩陣對(duì)應(yīng)元素相乘,1/2B表示對(duì)矩陣中每個(gè)元素進(jìn)行開方運(yùn)算,R為殘差矩陣,W為加權(quán)矩陣,其第m行n列元素

圖1給出了核長(zhǎng)1σ=時(shí)的加權(quán)因子。從圖中可以看出,當(dāng)殘差較大時(shí),加權(quán)因子迅速趨于0,可以有效抑制穩(wěn)定分布噪聲中的野點(diǎn)。

利用最小廣義相關(guān)熵準(zhǔn)則(M inimum Generalized Correntropy Criterion,MGCC),最小化式(21)就可以實(shí)現(xiàn)信號(hào)子空間的估計(jì)進(jìn)而實(shí)現(xiàn)空間譜的估計(jì):

圖1 加權(quán)因子幅度圖

3.4算法實(shí)現(xiàn)步驟

對(duì)比式(22)和式(11)發(fā)現(xiàn)可利用奇異值分解實(shí)現(xiàn)最小化式(21)的求解,但是由于加權(quán)矩陣W與待求解矩陣Y和Z有關(guān),無(wú)法單次獲得最優(yōu)解,因此本文采用IR-SVD[12]方法,主要步驟如下:步驟1隨機(jī)生成列滿秩矩陣Y(0)與行滿秩矩陣設(shè)定迭代次數(shù)k=0,初始誤差為,迭代停止誤差sε以及最大迭代次數(shù)K;

步驟6按式(24)計(jì)算空間譜,通過搜索峰值實(shí)現(xiàn)P個(gè)信源DOA的估計(jì)。

4 仿真實(shí)驗(yàn)

其中,s(t)表示信號(hào),γ為噪聲的分散系數(shù)。

本文采用兩個(gè)指標(biāo)對(duì)算法性能進(jìn)行評(píng)價(jià):DOA估計(jì)成功率和均方根誤差(Root M ean Square E r r o r,RMSE)。當(dāng)2個(gè)信源的入射角度估計(jì)誤差都不超過3°時(shí),則認(rèn)為此次DOA估計(jì)是成功的。估計(jì)成功率是估計(jì)成功次數(shù)與隨機(jī)實(shí)驗(yàn)次數(shù)之比。DOA估計(jì)的RMS?E定義為

本文同時(shí)仿真了文獻(xiàn)[2]中的經(jīng)典MUSIC算法、文獻(xiàn)[5]中基于分?jǐn)?shù)低階矩的FLOM-MUSIC算法、文獻(xiàn)[10]中基于類相關(guān)熵的CRCO-MUSIC算法、文獻(xiàn)[12]中基于范數(shù)的ACO-MUSIC算法以及本文的MGCC-MUSIC算法。

實(shí)驗(yàn)1算法收斂性能設(shè)定迭代停止條件εs=,K=100。圖2給出了不同噪聲環(huán)境下本文算法的迭代誤差、子空間距離以及代價(jià)函數(shù)的變化曲線。其中,子空間距離定義為從圖中可以看出,當(dāng)?shù)`差滿足時(shí),子空間距離以及代價(jià)函數(shù)基本穩(wěn)定,即以此為迭代停止條件是合適的。

實(shí)驗(yàn)2核長(zhǎng)的選擇圖3給出了不同噪聲環(huán)境下算法性能隨核長(zhǎng)σ的變化曲線,其中快拍數(shù)N=100。可以發(fā)現(xiàn),當(dāng)核長(zhǎng)滿足σ∈[1.3,2.5]時(shí),本文算法可以達(dá)到較好的估計(jì)結(jié)果。因此,在后續(xù)仿真實(shí)驗(yàn)中,核長(zhǎng)參數(shù)均取σ=1.5。此外,還可以發(fā)現(xiàn),當(dāng)噪聲的特征指數(shù)對(duì)算法性能影響很小。

圖2 算法迭代收斂曲線

實(shí)驗(yàn)3 GSNR對(duì)算法性能的影響圖4給出了不同GSNR下的估計(jì)結(jié)果,其中噪聲的特征指數(shù)α=1.5,快拍數(shù)N=100。可以看出,隨著GSNR的提高,所有算法性能均有顯著提高。與其它4種算法相比,經(jīng)典MUSIC算法在穩(wěn)定分布噪聲環(huán)境下的性能較差。在低信噪比環(huán)境下,MGCC-MUSIC算法與CRCO-MUSIC算法可以獲得更高的估計(jì)成功率,而且MGCC-MUSIC算法具有更低的估計(jì)誤差;而在高信噪比環(huán)境下,MGCC-MUSIC算法則與ACO-MUSIC算法性能相當(dāng)。

實(shí)驗(yàn)4快拍數(shù)對(duì)算法性能的影響圖5給出了快拍數(shù)N對(duì)算法性能的影響。設(shè)定SαS噪聲特征指數(shù)α=1.5,GSNR為4 dB。可以發(fā)現(xiàn),除經(jīng)典MUSIC算法外,其它4種算法的RMSE均隨著快拍數(shù)的增加而減小,其中MGCC-MUSIC,CRCO-MUSIC以及ACO-MUSIC 3種算法的估計(jì)成功概率相近,但是MGCC-MUSIC算法可以獲得更低的估計(jì)誤差。

實(shí)驗(yàn)5穩(wěn)定分布噪聲特征指數(shù)的影響圖6給出了不同噪聲特征指數(shù)下算法的估計(jì)結(jié)果,其中快拍數(shù)N=100,GSNR為4 dB。可以看出,在高脈沖性噪聲環(huán)境下,本文算法具有非常明顯的優(yōu)勢(shì)。此外,與圖2和圖3現(xiàn)象一致,噪聲特征指數(shù)α對(duì)本文算法性能的影響較小。

實(shí)驗(yàn)6角度分辨率圖7給出了不同算法的角度分辨能力。信源1的入射角度1θ從0°變化至18°,信源2的入射角度固定為220θ=°,噪聲特征指數(shù)1.5α=,GSNR為4 dB。從圖中可以看出,本文算法能夠分辨的最小角度差為6°,優(yōu)于其它算法。此外,還可以發(fā)現(xiàn),當(dāng)兩個(gè)信源的角度差較小時(shí),本文方法的估計(jì)成功率更高,而當(dāng)兩個(gè)信源的角度差較大時(shí),本文方法則在RMSE方面具有微弱的優(yōu)勢(shì)。

圖3 算法性能隨核長(zhǎng)的變化曲線

圖4 算法性能隨GSNR的變化曲線

圖5 算法性能隨快拍數(shù)的變化曲線

圖6 算法性能隨特征指數(shù)的變化曲線

圖7 算法的角度分辨率

5 結(jié)論

本文定義了一種新的廣義相關(guān)熵,證明了對(duì)稱穩(wěn)定分布隨機(jī)變量廣義相關(guān)熵的有界性,并基于最小廣義相關(guān)熵準(zhǔn)則提出了一種穩(wěn)定分布噪聲下的DOA估計(jì)新方法。該方法無(wú)需構(gòu)造協(xié)方差矩陣,通過求解優(yōu)化問題直接實(shí)現(xiàn)信號(hào)子空間的估計(jì)。為實(shí)現(xiàn)優(yōu)化問題的求解,本文給出了一種迭代優(yōu)化算法,并通過仿真實(shí)驗(yàn)分析了算法的收斂性。最后,本文分析了核長(zhǎng)、廣義信噪比、快拍數(shù)、噪聲特征指數(shù)以及信源間入射角度差對(duì)算法性能的影響。仿真結(jié)果表明,本文算法可以獲得更好的估計(jì)結(jié)果,尤其是在高脈沖性噪聲環(huán)境下具有更加明顯的優(yōu)勢(shì)。

[1]KRIM H and VIBERGM.Twodecadesofarray signalprocessing research:the parametric approach[J].IEEE Signal Processing Magazine,1996,13(4):64-97.doi:10.1109/79.526899.

[2]SCHM IDE R O.Multiple emitter location and signal parameter estimation[J].IEEE Transactions on Antennas and Propagation,1986,34(3):276-280.doi:10.1109/tap.1986.1143830.

[3]SHAO M and NIKIAS C L.Signal processing with fractional lower ordermoments:stable processes and their applications[J]. Proceedings of the IEEE,1993,81(7):986-1010.doi: 10.1109/5.231338.

[4]TSAKALIDESP and NIKIASC L.The robust covariation-based MUSIC(ROC-MUSIC)algorithm for bearing estimation in impulsive noise environments[J].IEEE Transactions on Signal Processing,1996,44(7):1623-1633.doi:10.1109/78.510611.

[5]LIU T and MENDEL JM.A subspace-based direction finding algorithm using fractional lower order statistics[J].IEEE Transactions on Signal Processing,2001,49(8):1605-1613.doi: 10.1109/78.934131.

[6]ZHANG JF and QIU T S.A novel covariation based noncircular sources direction finding method under impulsive noise environments[J].Signal Processing,2014,98:252-262.doi: 10.1016/j.sigpro.2013.11.031.

[7]BELKACEM I H and MARCOS S.Robust subspace-based algorithms for joint angle/Doppler estimation in non-Gaussian clutter[J].Signal Processing,2007,87(7):1547-1558.doi: 10.1016/j.sigpro.2006.12.015.

[8]YOUG H,QIU T S,and SONG AM.Noveldirection findings for cyclostationary signals in im pulsive noise environments[J]. Circuits,System s,and Signal Processing,2013,32(6):2939-2956. doi:10.1007/s00034-013-9597-0.

[9]張金鳳,邱天爽,宋愛民,等.Alpha穩(wěn)定分布噪聲環(huán)境下類M估計(jì)相關(guān)的DOA估計(jì)新算法[J].通信學(xué)報(bào),2013,34(5):71-78.doi: 10.3969/j.issn.1000-436x.2013.05.008.

ZHANG Jinfeng,QIU Tianshuang,SONG Aim in,et al. M-estimate like correlation based algorithm for direction ofarrival estimation under alpha-stable environm ents[J].Journal on Communications,2013,34(5):71-78.doi: 10.3969/j.issn.1000-436x.2013.05.008.

[10]ZHANG J F,QIU T S,SONG A M,et al.A novel correntropy based DOA estimation algorithm in impulsive noise environments[J].Signal Processing,2014,104:346-357.doi: 10.1016/j.sigpro.2014.04.033.

[11]邱天爽,張金鳳,宋愛民,等.脈沖噪聲下基于廣義類相關(guān)熵的DOA估計(jì)新方法[J].信號(hào)處理,2012,28(4):463-466.

QIU Tianshuang,ZHANG Jinfeng,SONG Aimin,et al.The generalized correntropy-analogous statistics based direction of arrival estimation in impulsive noise environments[J].Signal Processing,2012,28(4):463-466.

[12]ZENG W J,SO H C,and HUANG L.lp-MUSIC:Robust direction-of-arrivalestimator for impulsive noise environments[J]. IEEETransactionson SignalProcessing,2013,61(17):4296-4308. doi:10.1109/tsp.2013.2263502.

[13]LIU W F,POKHAREL P P,and PRINCIPE J C.Correntropy: properties and app lications in non-Gaussian signal processing[J]. IEEETransactionson SignalProcessing,2007,55(11):5286-5298. doi:10.1109/tsp.2007.896065.

[14]宋愛民,邱天爽,佟祉諫.對(duì)稱穩(wěn)定分布的相關(guān)熵及其在時(shí)間延遲估計(jì)上的應(yīng)用[J].電子與信息學(xué)報(bào),2011,33(2):494-498.doi: 10.3724/SP.J.1146.2010.00309.

SONG Aim in,QIU Tianshuang,and TONG Zhijian.Correntropy of the symmetric stable distribution and its application to the time dealy estimation[J].Journal of Electronics&Information Technology,2011,33(2):494-498.doi: 10.3724/SP.J.1146.2010.00309.

[15]WU Z Z,SHIJ H,ZHANG X,et al.Kernel recursivemaximum correntropy[J].Signal Processing,2015,117:11-26.doi: 10.1016/j.sigpro.2015.04.024.

[16]CHEN B D,LEIX,LIANG J L,etal.Steady-statemean-square error analysis for adaptive filtering under the maximum correntropy criterion[J].IEEE Signal Processing Letters,2014,21(7):880-884.doi:10.1109/LSP.2014.2319308.

[17]WANG L F and PAN C H.Robust level set image segm entation via a local correntropy-based K-means clustering[J].Pattern Recognition,2014,47(5):1917-1925.doi: 10.1016/j.patcog.2013.11.014.

[18]MELIA U,GUAITA M,VALLVERDU M,et al.Correntropy measures to detect daytim e sleepiness from EEG signals[J]. Physiological Measurements,2014,35(10):2067-2083.doi: 10.1088/0967-3334/35/10/2067.

[19]MARKOVSKY I.Low Rank Approximation:Algorithms,Im plementation,Applications[M].London,UK,Springer-Verlag London,2011:60-64.

王鵬:男,1989年生,博士生,研究方向?yàn)殛嚵行盘?hào)處理、非高斯非平穩(wěn)信號(hào)處理.

邱天爽:男,1954年生,博士,教授,博士生導(dǎo)師,主要研究方向?yàn)榉瞧椒€(wěn)非高斯統(tǒng)計(jì)信號(hào)處理、數(shù)字信號(hào)處理等.

任福全:男,1984年生,博士生,研究方向?yàn)閳D像處理、成像與壓縮感知等.

李景春:男,1966年生,博士,高級(jí)工程師,主要研究方向?yàn)闊o(wú)線電監(jiān)測(cè)理論與應(yīng)用等.

譚海峰:男,1977年生,高級(jí)工程師,主要研究方向?yàn)闊o(wú)線電頻譜工程以及無(wú)線電監(jiān)測(cè)等.

A Novel Generalized Correntropy Based Method for Direction of Arrival Estimation in Symmetric Alpha Stable Noise Environments

WANG Peng①Q(mào)IU Tianshuang①REN Fuquan①LIJingchun②TAN Haifeng②③①
①(Faculty ofElectronic Information and Electrical Engineering,Dalian University of Technology,Dalian 116024,China)
②(State Radio Monitoring Center,Beijing 100037,China)
③(Schoo l of Inform ation and Communication Engineering,Beijing University of Posts and Telecommunications,Beijing 100876,China)

To overcome the lim itation that the alpha stab le distributed variable possesses infinite second-order moment,a novelgeneralized correntropy is defined and thebounded property of the generalized correntropy for the symmetric alpha stab le variable is proved.Furthermore,a novelm inimum generalized correntropy criterion based DOA estim ation method for im pu lsive noise is proposed,and an iterative optim ization algorithm is presented,the convergence ofwhich is analyzed by simulation experiments.The simulation resu lts demonstrate that the proposed method can get better estimation resu lts than the fractional lower order moments based FLOM-MUSIC,the correntropy-like based CRCO-MUSIC and the lp norm based ACO-MUSIC methods,especially in the highly im pulsive noise environments.

Direction O f Arrival(DOA)estimation;Correntropy;Generalized correntropy;A lpha stable distribu ted noise;MUSIC algorithm

s:The National Natural Science Foundation of China(61139001,61172108,81241059)

TN911.7

A

1009-5896(2016)08-2007-07

10.11999/JEIT 151217

2015-11-03;改回日期:2016-03-03;網(wǎng)絡(luò)出版:2016-05-09

邱天爽qiutsh@d lut.edu.cn

國(guó)家自然科學(xué)基金(61139001,61172108,81241059)

猜你喜歡
定義方法
永遠(yuǎn)不要用“起點(diǎn)”定義自己
海峽姐妹(2020年9期)2021-01-04 01:35:44
定義“風(fēng)格”
學(xué)習(xí)方法
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
成功的定義
山東青年(2016年1期)2016-02-28 14:25:25
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
修辭學(xué)的重大定義
山的定義
主站蜘蛛池模板: 精品国产免费观看| 久久国产精品影院| 国产精品无码作爱| 久久精品亚洲专区| 欧美一区二区丝袜高跟鞋| av尤物免费在线观看| 国产精品久久精品| 国产成人精品免费视频大全五级| 丁香五月激情图片| 中文字幕在线免费看| 欧美激情综合一区二区| 精品久久久久成人码免费动漫| 精品国产免费观看一区| 亚洲精品无码成人片在线观看| 91久久国产热精品免费| 亚洲成人一区二区三区| 国产特级毛片aaaaaaa高清| 美女啪啪无遮挡| 亚洲精品无码高潮喷水A| 欧美中文字幕无线码视频| 国产精品亚洲αv天堂无码| 国产免费网址| 久久这里只有精品免费| 日韩欧美视频第一区在线观看| V一区无码内射国产| 国产99在线| 日韩成人在线网站| 欧美伦理一区| 色香蕉影院| 97超爽成人免费视频在线播放| 国产91视频免费观看| 国产激情影院| 国产精品区视频中文字幕| 无码福利视频| 全免费a级毛片免费看不卡| 久久综合色视频| 狠狠做深爱婷婷久久一区| 亚洲AV无码久久精品色欲| 亚洲最猛黑人xxxx黑人猛交 | 男人天堂伊人网| 成年女人a毛片免费视频| 黄片在线永久| 91日本在线观看亚洲精品| 精品国产免费观看一区| 久久精品国产免费观看频道| 国产美女人喷水在线观看| 亚洲天堂免费在线视频| 国产成人精品优优av| 欧美α片免费观看| 成人国产精品网站在线看| 国产精品亚洲欧美日韩久久| 精品無碼一區在線觀看 | 丁香婷婷久久| www.av男人.com| 91视频首页| 香蕉伊思人视频| 久久这里只有精品23| 久久无码av三级| 国产欧美日韩综合在线第一| 日本国产精品| 欧美不卡视频在线| 亚洲色欲色欲www网| 久久免费视频播放| 四虎影视库国产精品一区| 欧美一区精品| 97亚洲色综久久精品| 青青青视频免费一区二区| 九九热在线视频| 日韩欧美中文在线| 国产十八禁在线观看免费| 亚洲欧美在线精品一区二区| 久久夜色撩人精品国产| 中文无码毛片又爽又刺激| 亚洲日韩AV无码一区二区三区人| 看av免费毛片手机播放| 欧美A级V片在线观看| 99久久精品国产精品亚洲| 国产免费黄| 国禁国产you女视频网站| 国产一级视频久久| 精品伊人久久久久7777人| 亚洲二区视频|