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

空時頻域中欠定混合條件下的波達方向估計*

2015-06-21 12:39:37朱立為汪亞王翔黃知濤國防科技大學(xué)電子科學(xué)與工程學(xué)院湖南長沙40073電子信息系統(tǒng)復(fù)雜電磁環(huán)境效應(yīng)國家重點實驗室河南洛陽47003
國防科技大學(xué)學(xué)報 2015年5期
關(guān)鍵詞:信號

朱立為,汪亞,王翔,黃知濤(.國防科技大學(xué)電子科學(xué)與工程學(xué)院,湖南長沙40073;.電子信息系統(tǒng)復(fù)雜電磁環(huán)境效應(yīng)國家重點實驗室,河南洛陽47003)

空時頻域中欠定混合條件下的波達方向估計*

朱立為1,2,汪亞1,王翔2,黃知濤2
(1.國防科技大學(xué)電子科學(xué)與工程學(xué)院,湖南長沙410073;
2.電子信息系統(tǒng)復(fù)雜電磁環(huán)境效應(yīng)國家重點實驗室,河南洛陽471003)

波達方向估計是陣列信號處理領(lǐng)域的熱點問題,但經(jīng)典的波達方向估計方法通常要求陣元數(shù)大于源信號個數(shù),即滿足超定條件,而在實際中往往面臨的是源信號個數(shù)大于陣元數(shù)的欠定條件。基于此,提出了一種基于空間時頻分布的多重信號分類擴展算法,通過將空間時頻分布矩陣進行擴展,實現(xiàn)了欠定條件下的波達方向估計。相比時頻多重信號分類算法,所提算法能同時適應(yīng)超定和欠定條件;相比已有的欠定波達方向估計方法,其不但保證了波達方向估計的精度,而且放寬了對源信號稀疏性的要求,同時還降低了對快拍數(shù)的要求。仿真實驗結(jié)果證明了該方法的有效性。

波達方向;時頻分布;欠定

在陣列信號處理領(lǐng)域,波達方向(Direction Of Arrival,DOA)估計是聲吶、雷達、地震學(xué)的一個重要研究課題。傳統(tǒng)的DOA估計方法,如多重信號分類(MUltiple SIgnal Classification,MUSIC)方法[1]、改進MUSIC方法[2-3]、時頻MUSIC方法[4]、極大似然方法[5]、空間譜估計方法[6]、波束成形方法[7]等,一般要求信號環(huán)境滿足超定條件,即接收天線的陣元個數(shù)要大于潛在的源信號個數(shù)。但在實際的信號環(huán)境下,這一條件并不總能滿足。例如,在機載或星載的非合作通信應(yīng)用中,由于地面雷達、通信等各種輻射源的大量使用,并且使用的頻段不斷擴展、相互重疊,加上各種自然輻射產(chǎn)生的無意干擾,使得機載或星載接收設(shè)備面臨時域高度密集、頻域嚴(yán)重混疊、空間相互交錯的復(fù)雜電磁環(huán)境,而實際上機載或星載設(shè)備本身受體積限制,陣元數(shù)不能隨意增加,且陣元數(shù)目越多,接收設(shè)備結(jié)構(gòu)越復(fù)雜,價格也越昂貴。因此需要利用有限的接收陣元對盡可能多的源信號的DOA進行估計,即需要研究可以適用于欠定混合條件下的DOA估計方法。

近年來,相繼出現(xiàn)了一些適應(yīng)欠定混合條件的DOA估計算法。主要可分為兩類:一是基于高階累積量的欠定DOA估計算法[8-9]。文獻[9]利用四階累積量代替協(xié)方差矩陣,將MUSIC算法擴展成4-MUSIC算法(4-MUSIC)。為了進一步提高陣列對多個潛在源信號的DOA估計能力,文獻[10]將四階累積量推廣到任意偶數(shù)階累積量,提出了2q-MUSIC(q>1)算法,并詳細(xì)分析了采用不同陣列結(jié)構(gòu)等價的虛擬陣元數(shù)目。雖然這些高階累積量算法可以通過擴展產(chǎn)生虛擬的陣元以適應(yīng)欠定混合條件,但是高階累積量的計算過程復(fù)雜且需要較多的樣本數(shù)目才能保證算法的估計精度。二是其他一些基于時頻域稀疏性的DOA估計算法,如退化分解估計技術(shù)(Degenerate Unmixing Estimation Technique,DUET)等也能適應(yīng)欠定混合條件[11-15],但這些方法對信號的時頻域稀疏性要求過于嚴(yán)格,即任意時頻點只有一個源信號起主導(dǎo)作用,其余源信號取值為0[16],這在實際應(yīng)用中并不能完全滿足。

本文提出了基于空間時頻分布(Spatial Time-Frequency Distributions,STFD)的MUSIC算法(STFD-MUSIC),通過擴展空間時頻混合矩陣[17]來實現(xiàn)欠定條件下的DOA估計。該方法并不需要假設(shè)信號是稀疏的,且能夠降低算法對快拍數(shù)的要求。

1 算法分析

假設(shè)N個窄帶遠(yuǎn)場信號S(t)=[s1(t),…,sN(t)]T入射到由M個陣元組成的天線陣上,觀測信號可以寫為:

其中,X(t)=[x1(t),…,xM(t)]T和V(t)=[v1(t),…,vM(t)]T分別表示M個陣元的輸出即觀測信號和噪聲,θ=[θ1,…,θN]代表各源信號的入射角度,混合矩陣A(θ)=[a1(θ1),…,aN(θN)]表示傳輸矩陣,作為線性矩陣,其K列可以表示為:

本文研究的欠定條件指的是陣元數(shù)小于源信號數(shù)目,即N>M;研究目的是估計θ。

1.1 空間時頻分布基本理論

本文算法的核心思路是利用觀測信號的空間時頻分布矩陣實現(xiàn)陣列的虛擬擴展,實現(xiàn)欠定混合條件下的DOA估計。

空間時頻分布矩陣定義如式(3)、式(4):

聯(lián)立式(4)和式(5),可得接收信號的空間時頻分布矩陣與源信號的空間時頻分布矩陣之間的關(guān)系,其表示如式(6):

由文獻[18]易知,接收信號的空間時頻分布矩陣的對角元素稱為自源點(auto-terms),也稱自項,其對應(yīng)著X(t)中代表各個信號的向量平方;而反對角元素稱為互源點(cross-terms),對應(yīng)著多個信號向量的線性組合。

1.2 算法假設(shè)條件

為了能夠在時頻域上估計出源信號的DOA,采用了Wigner-Vill分布。同時,假設(shè)源信號滿足以下條件:

假設(shè)1:在時頻平面上信號的自源時頻點與互源時頻點幾乎是不混疊的;

根據(jù)假設(shè)2可知本文算法在天線陣元數(shù)一定時可接收處理的最大源信號個數(shù)。

1.3 算法描述

1.3.1 自源點選擇

對于任意自源點(t,f)∈Ωs(Ω表示整個時頻平面),接收信號X(t)的空間時頻分布矩陣可以表示為:

式中,DˉSS(t,f)表示對角矩陣,其表示為DˉSS(t,f)=diag[Ds1s1(t,f),…,DsNsN(t,f)]。

先找出自源點,同時抑制互源點,則可以實現(xiàn)對各個源信號的DOA估計。利用的準(zhǔn)則如式(8)所示:

在超定混合條件下,直接對DXX(t,f)進行奇異值分解可以很容易地估計出源信號的DOA。但在欠定混合條件下,不能直接分解來估計源信號的DOA。下面將介紹一種能適用于欠定混合條件下基于空間時頻分布的DOA估計算法。1.3.2基于空間時頻分布的擴展矩陣構(gòu)造

根據(jù)式(8)所述的準(zhǔn)則,求出自源點的個數(shù),設(shè)在式(8)準(zhǔn)則下找出K個自源點,則K個自源點對應(yīng)的時頻分布矩陣表示如式(9)所示:

定義矩陣D∈K×M,其中(D)km=(Dk)mm(k =1,…,K,m=1,…,M),則C可以表示成如式(11)所示的形式:

式中,A⊙A*=[a1θ1),…,aN(θN)?θN)],其中*,⊙和?分別表示復(fù)數(shù)的轉(zhuǎn)置、ri-Rao乘積和Kronecker乘積。因此,要解決的問題就是當(dāng)N>M時,如何利用矩陣C來估計角度θ。

一般來說,假設(shè)N≤min(M2,K)。根據(jù)文獻[19]中的證明可知,當(dāng)潛在源信號個數(shù)與陣元數(shù)滿足假設(shè)2的關(guān)系式時,則A⊙A*和D是列滿秩的。而A⊙A*和D是列滿秩的,意味著C也是列滿秩的,并且C的秩就等于源信號個數(shù)N。

1.3.3 子空間求解與角度估計

與傳統(tǒng)的MUSIC方法類似,可以直接對自相關(guān)陣RC=CCH進行奇異值分解來獲得噪聲子空間。在本文算法中,雖然信號子空間和噪聲子空間可以通過任意一個RCi進行奇異值分解求得,但是當(dāng)RCi奇異性較強時,子空間的估計精度會很差,進而影響DOA的估計結(jié)果。因此,提出對擴展自相關(guān)矩陣集合{RCi|i=1,…,K}進行聯(lián)合對角化的方法求解信號子空間和噪聲子空間。相比單一矩陣的處理方法,利用多個矩陣的估計方法可以同時保證算法的估計精度和信噪比適應(yīng)能力。采用聯(lián)合對角化的方法分解自相關(guān)陣RC= CCH,其分解后的主要形式如式(12):

式中,RDi=(,Λ為×對角陣。Λ的N個較大對角對應(yīng)著信號子空間,M2-N個較小對角線元素代表著噪聲子空間。因此,可以利用V矩陣相應(yīng)的M2-N列向量構(gòu)成的噪聲子空間來估計源信號的DOA。根據(jù)文獻[20]中的研究,可以使用式(13)所示的代價函數(shù)進行聯(lián)合對角化。

當(dāng)式(13)的值最小時,可以估計出信號子空間與噪聲子空間。因為RCi是Hermitian陣,所以RCi中組成信號子空間的列向量與組成噪聲子空間的列向量是正交的。詳細(xì)來說,Span{VS}= Span{Α⊙Α*},其中VS是N維信號子空間的基,則Α⊙Α*中所有列向量都與VN中的列向量正交,而且是M2-N維噪聲子空間的基。因此,所有{ai(θi)(θi),1≤i≤N}向量都與VN中的列向量正交,因此有式(14)所示的關(guān)系式:

在實際運用過程中,常常定義空間偽譜的概念來運用上述正交性進行DOA估計,如式(15):

從式(15)可知,只要找出最大“譜線”P(θ)的位置處對應(yīng)的角度,則這個角度就等于要估計的某一源信號的DOA。

綜上所述,本文所提出的基于空間時頻分布的欠定混合條件下DOA估計的核心算法步驟如下:

(i-1)M+j,k)ij計算拓維后的新矩陣,同時求出C^的協(xié)方差矩陣=,然后對協(xié)方差矩陣進行聯(lián)合對角化,估計出VN。

2 仿真分析

2.1 評價準(zhǔn)則

為了驗證算法的性能,仿真時采用平均測向均方根誤差作為指標(biāo)來衡量DOA的估計性能,定義平均測向均方根誤差(Root-Mean-Square Error,RMSE):

2.2 仿真實驗

仿真分析中,用4個線性調(diào)頻信號作為源信號,其頻率分別為:1000MHz,1000.1MHz,1000.4MHz,1000.5MHz,碼速率都為200kbit/s,角度分別為:-43°,-13°,26°,48°。仿真實驗1到4中接收天線是陣元數(shù)目為3的均勻線陣,仿真實驗5中接收天線是陣元數(shù)目分別為4和5的均勻線陣,相鄰陣元間距為半波長。為降低采樣率和計算量,先將射頻信號的頻率變到中頻,中頻頻率分別為350kHz,450kHz,750kHz和850kHz,采樣率為2MHz。仿真將[-90°,90°]的空域以0.1°間隔均勻采樣,得到離散的假設(shè)角度集θ~,進行500次蒙特卡洛試驗。

在仿真實驗中,對本文算法與4-MUSIC算法進行了性能比較,而沒有將本文算法與基于時頻稀疏性的DOA估計算法進行比較,這是因為實驗設(shè)置的源信號在時頻域是重疊的,不滿足時頻域稀疏性,基于稀疏性的DOA估計算法無法完成源信號DOA的估計。

在圖1中,仿真實驗給出了信噪比(Signal to Noise Ratio,SNR)為10dB,采樣點數(shù)分別為1024和512時,本文算法與4-MUSIC算法的任意10次仿真實驗的空域偽譜估計結(jié)果。從圖1(a)可以看出本文算法的空間譜中出現(xiàn)了4個顯著的譜峰,能夠清晰地分辨出四個源信號的DOA;圖1(b)是4-MUSIC算法的空間譜,圖中出現(xiàn)了4個大致的譜峰,勉強能分辨出四個源信號的DOA。當(dāng)采樣點數(shù)降低時,從圖1(c)中可以看出,本文算法仍然可以清晰地分辨出源信號的DOA;而圖1(d)中4-MUSIC算法的空間譜上的譜峰已很不明顯,難以分辨出源信號的DOA。

圖1本文算法與4-MUSIC算法空間偽譜Fig.1 Spatial pseudo-spectrums of the proposed and 4-MUSIC algorithm

圖2 為信噪比為10dB時,本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨采樣點數(shù)變化的曲線(圖中CRLB為克拉美羅下限)。從圖中可看出,隨著采樣點數(shù)的增加,平均測向均方根誤差隨之減少;在相同的采樣點數(shù)條件下,本文算法的平均測向均方根誤差更小,DOA估計精度更高;而在相同的DOA估計精度下,本文算法需要采樣點數(shù)要少,因為在相同的估計性能條件下,相比4-MUSIC的高階累積量的估計,本文算法的空間時頻分布矩陣估計需要的采樣數(shù)據(jù)更少。

圖2不同采樣點數(shù)的RMSE估計結(jié)果Fig.2 Estimation results of RMSE for different sample points

圖3 為采樣點數(shù)等于1024、信噪比從-5dB到30dB變化時,本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化曲線。從圖中可以看出,隨著信噪比的變化,本文算法的平均測向均方根誤差更接近CRLB,DOA估計精度更高。當(dāng)DOA估計精度相同時,本文算法比4-MUSIC算法適應(yīng)的信噪比更低。

圖3 不同的SNR的估計結(jié)果(欠定條件)Fig.3 Estimation results of RMSE for different SNR(underdetermined)

假設(shè)陣列誤差矢量ei(1≤i≤2)是相互獨立的高斯隨機變量,即滿足E{eieHj}=σ2δIN,令σ2= 0.000 3,圖4給出了本文算法和4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化的曲線。從圖中可以看出,相比4-MUSIC算法,本文算法的陣列誤差適應(yīng)能力更強。

圖4不同的SNR的估計結(jié)果(考慮陣列誤差)Fig.4 Estimation results of RMSE for different SNR(model error)

圖5 為采樣點數(shù)設(shè)為1024,信噪比從-5dB到30dB變化,陣元數(shù)分別為5和6時,本文算法與4-MUSIC算法的平均測向均方根誤差RMSE隨信噪比變化曲線。結(jié)合圖3的結(jié)果可以看出,本文的DOA估計算法可以有效應(yīng)用于適定和超定的接收條件,極大地擴展了算法的實際應(yīng)用范圍,并且隨著陣元數(shù)的增加,DOA估計性能也隨之提高。

3 結(jié)論

針對欠定混合條件下DOA估計問題,提出了一種基于空間時頻域的估計算法。該方法通過Khatri-Rao乘積和Kronecker乘積對空間時頻混合矩陣進行擴展,以適應(yīng)欠定混合條件。提出的欠定DOA估計算法中,采用聯(lián)合對角化的方法來處理空間時頻矩陣,增強了算法的穩(wěn)健性,提高了算法在低信噪比條件下的DOA估計能力。仿真結(jié)果表明,所提算法的DOA估計性能優(yōu)于4-MUSIC算法。同時,隨著天線陣元數(shù)的增加,所提算法對DOA的估計性能還會大幅度提升。

圖5 不同的SNR的估計結(jié)果(接收陣元為5,6)Fig.5 Estimation results of RMSE for different SNR(5 and 6 array elements)

References)

[1]Schmidt R O.Multiple emitter location and signal parameter estimation[J].IEEE Transactions on Antennas and Propagation,1986,34(3):276-280.

[2]Zhou Q C,Gao H T,Wang F,etal.Modified DOA estimation methods with unknown source number based on projection pretransformation[J].Progress in Electromagnetics Research B,2012,38:387-403.

[3]Zhou Q C,Gao H T,Wang F.A high resolution DOA estimatingmethod without estimating the number of sources[J].Progress in Electromagnetics Research C,2012,25:233-247.

[4]Belouchrani A,Amin M G.Time-frequency MUSIC[J].IEEE Signal Processing,Letters,1999,6(5):109-110.

[5]Ziskind I,Wax M.Maximum likelihood localization of multiple sources by alternating projection[J].IEEE Transactions on Acoustics Speech&Signal Processing,1988,36(10):1553-1560.

[6]Malioutov D,Cetin M,Willsky A S.A sparse signal reconstruction perspective for source localization with sensor arrays[J].IEEE Transactions on Signal Processing,2005,53(8):3010-3022.

[7]Li J,Stoica P.Robust adaptive beamforming[J].Series in Telecommunications&Signal Processing,2006,15(7): 2345-2348.

[8]Chevalier P,Albera L,F(xiàn)erréol A,et al.On the virtual array concept for higher order array processing[J].IEEE Transactions on Signal Processing,2005,53(4):1254-1271.

[9]Porat B,F(xiàn)riedlander B.Direction finding algorithms based on high-order statistics[J].IEEE Transactions on Signal Processing,1991,39(9):2016-2024.

[10]Chevalier P,F(xiàn)erreol A,Albera L.High-resolution direction finding from higher order statistics:the 2-MUSIC algorithm[J].IEEE Transactions on Signal Processing,2006,54(8):2986 -2997.

[11]Rickard S,Dietrich F.DOA estimation of many W-disjoint orthogonal sources from two mixtures using DUET[C]// Proceedings of the Tenth IEEEWorkshop on Statistical Signal and Array Processing,IEEE,2000:311-314.

[12]Araki S,Sawada H,Mukai R,et al.DOA estimation for multiple sparse sources with normalized observation vector clustering[C]//Proceedings of 2006 IEEE International Conference on Acoustics,Speech and Signal Processing,ICASSP 2006,IEEE,2006.

[13]Zhang W Y,Rao B D.A twomicrophone-based approach for source localization of multiple speech sources[J].IEEE Transactions on Audio Speech&Language Processing,2010,18(8):1913-1928.

[14]Zhou Z,Lu S J,Zhang E Y,et al.Underdetermined DOA estimation of LFM signals[C]//Proceedings of 2012 Second International Conference on Instrumentation,Measurement,Computer,Communication and Control(IMCCC),IEEE,2012:869-872.

[15]Lie JP,Ng B P,See C M S.Multiple UWB emitters DOA estimation employing time hopping spread spectrum[J].Progress in Electromagnetics Research,2008,78:83-101.

[16]王翔,黃知濤,任嘯天,等.基于時頻單源點檢測和聚類驗證技術(shù)的欠定混合盲辨識算法[J].國防科技大學(xué)學(xué)報,2013,35(2):69-74.WANG Xiang,HUANG Zhitao,REN Xiaotian,et al.Blind identification of underdetermined mixtures based on detection of time frequency single source point and cluster validation technique[J].Journal of National University of Defense Technology,2013,35(2):69-74.(in Chinese)

[17]Lu F B,Huang Z T,Jiang W L.Underdetermined blind separation of non-disjoint signals in time-frequency domain based onmatrix diagonalization[J].Signal Processing,2011,91(7):1568-1577.

[18]Belouchrani A,Amin M G,Thirion-Moreau N,et al.Source separation and localization using time-frequency distributions: an overview[J].IEEE Signal Processing Magazine,2013,30(6):97-107.

[19]De Lathauwer L.A link between the canonical decomposition in multilinear algebra and simultaneousmatrix diagonalization[J].Siam Journal on Matrix Analysis and Applications,2006,28 (3):642-666.

[20]Golub G H,Van Loan C F.Matrix computation[M].2nd ed.USA:the John HopKins University Press,1989.

Underdeterm ined direction of arrival estimation based on spatial time-frequency distributions

ZHU Liwei1,2,WANG Ya1,WANG Xiang2,HUANG Zhitao2
(1.College of Electronic Science and Engineering,National University of Defense Technology,Changsha 410073,China;2.The State Key Laboratory of Complex Electromagnetic Environment Effects on Electronics and Information System,Luoyang 471003,China)

In the field of array signal processing,direction of arrival(DOA)estimation is a hotspot problem.Classical DOA estimation methods usually require the number of sensor should be larger than the source signals’(which the so-called over-determined case is).However,whatwe encounter in practice is always the underdetermined case in which the number of source signal is larger than the sensors’.To solve the problem,amultiple signal classification(MUSIC)extension algorithm based on spatial time-frequency distribution was proposed to achieve the underdetermined DOA estimation by expanding the dimension of the spatial time-frequency distributionsmatrices.Compared with the existing timefrequency MUSIC,the proposed algorithm can be applied to both the over-determined and the underdetermined cases.The proposed algorithm also has advantages over the existing underdetermined DOA estimation methods for it guarantees the estimation precision,relaxes the requirements for source signal sparseness and lowers standards of the number of snapshots.Simulation results confirm the validity and high performance of the proposed algorithm.

direction of arrival;time-frequency distributions;underdetermined

TN911

A

1001-2486(2015)05-149-06

10.11887/j.cn.201505023

http://journal.nudt.edu.cn

2014-09-28

CEMEE國家實驗室開放課題基金資助項目(2014K104B);國家自然科學(xué)基金資助項目(61401490)

朱立為(1985—),男,湖南郴州人,博士研究生,E-mail:iendwin@163.com;黃知濤(通信作者),男,教授,博士,博士生導(dǎo)師,E-mail:taldcn@yahoo.com.cn

猜你喜歡
信號
信號
鴨綠江(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信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: www成人国产在线观看网站| 亚洲欧洲综合| 青青青视频免费一区二区| 国产91精品调教在线播放| 国产精欧美一区二区三区| 欧美区国产区| 国产成人精品午夜视频'| 漂亮人妻被中出中文字幕久久| 欧美狠狠干| 欧美精品在线观看视频| 高清大学生毛片一级| 2048国产精品原创综合在线| 久久人体视频| 九九九久久国产精品| 亚洲综合色区在线播放2019| 亚洲色中色| 亚洲天堂视频网| 狠狠色婷婷丁香综合久久韩国| 99偷拍视频精品一区二区| 欧美a在线| 亚洲成年网站在线观看| 蜜臀av性久久久久蜜臀aⅴ麻豆| 婷婷亚洲综合五月天在线| 国产一级在线播放| 国产高潮流白浆视频| 久久综合九色综合97网| 91探花国产综合在线精品| 91蝌蚪视频在线观看| 97狠狠操| 久久亚洲日本不卡一区二区| P尤物久久99国产综合精品| 日韩欧美国产精品| 久久这里只有精品23| 午夜免费小视频| 毛片卡一卡二| 99久久性生片| 色婷婷丁香| 精品无码一区二区三区电影| 中国毛片网| 国产人成在线观看| 99久久亚洲精品影院| AV无码国产在线看岛国岛| AⅤ色综合久久天堂AV色综合| 日韩麻豆小视频| 婷婷六月综合网| 又大又硬又爽免费视频| 国产精品.com| 亚洲aⅴ天堂| 中文字幕 91| 欧美日韩在线第一页| 国产va免费精品观看| 亚洲一级毛片免费看| 99久久精品免费看国产免费软件| 91精品最新国内在线播放| 国产jizz| 国产一区二区网站| 无码 在线 在线| 国产亚洲视频在线观看| 激情综合激情| 天天色综网| 青青草原国产免费av观看| 亚洲一区二区三区香蕉| 国产欧美成人不卡视频| 欧美成人在线免费| 一区二区三区在线不卡免费| 国产精品白浆在线播放| 成人噜噜噜视频在线观看| 国产成人亚洲无吗淙合青草| 国产日韩精品欧美一区灰| 亚洲第一网站男人都懂| www.精品国产| 亚洲乱伦视频| 亚洲色图综合在线| 国内精品自在自线视频香蕉 | 国产精品综合久久久| 国产男人的天堂| 亚洲无码在线午夜电影| 成人福利在线免费观看| 欧美亚洲一二三区| AV无码国产在线看岛国岛| 国产人免费人成免费视频| 69综合网|