亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于互質(zhì)陣虛擬陣列空間平滑的相干信號DOA估計方法

        2022-04-07 12:10:16嚴世安寧方立
        關(guān)鍵詞:互質(zhì)協(xié)方差信噪比

        韋 娟, 嚴世安, 寧方立

        (1. 西安電子科技大學(xué)通信工程學(xué)院, 陜西 西安 710071; 2. 西北工業(yè)大學(xué)機電學(xué)院, 陜西 西安 710072)

        0 引 言

        波達方向(direction of arrival,DOA)估計是一種陣列接收信號獲取信號入射方向信息的技術(shù),被廣泛應(yīng)用在雷達、聲納、無線通信等領(lǐng)域。目前多以非相干信號為研究對象,但在實際信號傳播過程中,由于多徑傳播和同頻干擾造成大量相干信號源存在,因此相干信號的DOA估計亦是一研究熱點。

        均勻線陣相干信號DOA估計算法中,以空間平滑、Toeplitz矩陣重構(gòu)為基礎(chǔ)的解相干算法,均存在陣列孔徑損失問題??梢酝ㄟ^增加陣元數(shù)來提高DOA估計算法的角度自由度(degrees of freedom,DOF),但該方法會增加系統(tǒng)部署成本;亦可引入如互質(zhì)陣列、嵌套陣列等非均勻稀疏陣列構(gòu)建虛擬陣列以拓展陣列孔徑?;ベ|(zhì)陣列是由兩個陣元數(shù)分別為互質(zhì)數(shù)的均勻子陣相互穿插形成,結(jié)合了稀疏陣列的優(yōu)勢和質(zhì)數(shù)的性質(zhì),能夠突破奈奎斯特采樣定理的限制,陣列部署靈活,具有廣闊的應(yīng)用前景。

        前向空間平滑(forward spatial smoothing,FOSS)算法是最早被提出來估計相干信號DOA的,但該方法孔徑損失嚴重。前后向空間平滑(forward and backward spatial smooth,FBSS)算法利用了后向子陣列,增加了平滑次數(shù),增大了陣列孔徑,擁有比FOSS算法更高的DOF,但對信號沒有進行完全的解相干。Han等提出了一種類旋轉(zhuǎn)不變子空間信號參數(shù)估計(estimating signal parameters via rotational invariance techniques-like,ESPRIT-like)算法,將數(shù)據(jù)協(xié)方差矩陣的任意一行構(gòu)造成Toeplitz矩陣,能夠完全消除信號的相干性,利用ESPRIT-like算法進行DOA估計,但該算法估計精度低,亦存在陣列孔徑損失。Stoica等提出的最大似然(maximum likelihood,ML)估計算法以及Ottersten等提出的加權(quán)子空間擬合(weighted subspace fitting, WSF)算法均可以在不損失陣列孔徑的前提下對相干信號進行DOA估計,但這兩種算法包含多維非線性搜索,計算復(fù)雜度高,且初值設(shè)置對于DOA測向影響較大。Wang等提出的通過四階累積重構(gòu)Toeplitz矩陣進行相干信號DOA估計,DOA估計精度較高,但計算復(fù)雜度高。Zhang等提出了一種基于多重Toeplitz矩陣重構(gòu)的DOA估計方法,對相關(guān)矩陣的平方進行加權(quán)求和,以形成滿秩等效數(shù)據(jù)協(xié)方差矩陣。該算法無需提前消噪,但低信噪比、低快拍條件下估計性能不佳。唐曉杰等提出了一種前后向多重Toeplitz矩陣重構(gòu)(forward and backward multi-Toeplitz,FBMT)算法,構(gòu)造了 Toeplitz 矩陣,通過誤差最小準則構(gòu)造代價函數(shù),實現(xiàn)相干信號的 DOA 估計。該方法無需信源數(shù)先驗信息,但陣列孔徑損失嚴重,低信噪比條件下估計性能較差。Peng等提出一種將協(xié)方差矩陣子矩陣對角線元素之和的共軛作為加權(quán)因子的自加權(quán)空間平滑(self-weighted spatial smoothing,SWSS)算法,信號估計精度較高,由于采用空間平滑預(yù)處理,依舊存在陣列孔徑損失。

        基于此,對互質(zhì)陣虛擬陣列的連續(xù)部分進行空域平滑,提出一種基于互質(zhì)陣列虛擬陣列空間平滑(coprime array virtual array spatial smoothing,CASS)的相干信號DOA估計方法,該算法不需要信源數(shù)先驗信息,利用互質(zhì)陣虛擬陣列可以實現(xiàn)DOF拓展的特點,搭載空間平滑算法進行相干信號DOA估計。在低信噪比環(huán)境下對于信號的DOA估計精度以及魯棒性擁有較好表現(xiàn)。

        1 陣列信號模型

        1.1 互質(zhì)陣列接收信號模型

        互質(zhì)陣列利用虛擬域信號處理方法,在虛擬域上形成一個增廣的虛擬陣列,對該虛擬陣列信號進行與均勻線陣等價的統(tǒng)計信號處理,能夠有效提升DOA估計的DOF。

        互質(zhì)陣列結(jié)構(gòu)如圖1所示,該陣列由兩個均勻線陣嵌套構(gòu)成。

        圖1 互質(zhì)陣列Fig.1 Coprime array

        其中,,為互質(zhì)的兩個數(shù),虛擬陣元間距=2?;ベ|(zhì)子陣1為2元均勻線陣,陣元間距為,子陣1陣元位置為{|0≤≤2-1}?;ベ|(zhì)子陣2為元均勻線陣,陣元間距為,子陣2陣元位置為{|0≤≤-1}。

        互質(zhì)陣列物理陣元真實位置為=[,,…,2+-2]。

        互質(zhì)陣虛擬等效差聯(lián)合陣列如圖2所示。

        圖2 互質(zhì)陣列虛擬陣列Fig.2 Virtual array of coprime array

        虛擬陣元位置為={ ±(-)|0≤≤2-1, 0≤≤-1}。其中,虛擬陣列陣元位置-(+-1)與(+-1)之間無空洞,可以視為均勻線陣。

        個信號入射到互質(zhì)陣列接收信號為

        (1)

        對于第個相干信號可表示為

        ()=(),=2,3,…,

        式中:()為原始參照信號;為相干系數(shù)。

        1.2 空間前后向平滑處理

        互質(zhì)陣列協(xié)方差矩陣為

        (2)

        對接收協(xié)方差矩陣進行向量化處理:

        (3)

        式中:=vec (),為單位矩陣;=[()?(),…,()?()],?表示Kronecker內(nèi)積,(·)表示共軛。

        通過在Toeplitz矩陣結(jié)構(gòu)中重新排列的元素可獲得虛擬陣列信號的等價協(xié)方差。由于信噪比和快拍數(shù)對協(xié)方差矩陣的影響,同一波程差對應(yīng)的協(xié)方差矩陣元素不同,故對同一波程差對應(yīng)的元素取平均值:

        (4)

        式中:(diff)表示同一波程差diff對應(yīng)的第個協(xié)方差矩陣元素;(diff)表示波程差diff相同的協(xié)方差元素總和。

        在互質(zhì)虛擬陣列空洞位置插入天線,將聯(lián)合差陣列變?yōu)榫鶆蚓€陣后,假設(shè)空洞處天線無接收信號。根據(jù)波程差元素對協(xié)方差矩陣進行擴展,構(gòu)成Toeplitz矩陣(2-+1)×(2-+1)。

        (5)

        式中:1<,<2-+1。

        盡管矩陣包含所有互質(zhì)陣列信號信息,但由于稀疏陣列特性,是低秩矩陣,把看作部分數(shù)據(jù)信息缺失的協(xié)方差矩陣,利用低秩矩陣恢復(fù)理論,對矩陣中的零元素進行近似填充。

        創(chuàng)建優(yōu)化問題:

        (6)

        式中:()是以復(fù)向量為第一列的Hermitian Toeplitz矩陣,(2-+1)×1;為協(xié)方差矩陣擬合誤差。

        由于矩陣秩函數(shù)是非連續(xù)、非凸函數(shù),利用跡范數(shù)最小化對式(6)進行凸松弛表示:

        (7)

        構(gòu)建映射矩陣,(2-+1)×(2-+1)。

        (8)

        由于,向量的首元素無法進行共軛拓展,故將第一個元素初值先定義為實數(shù),經(jīng)CVX工具箱計算得到的整體優(yōu)化向量后,再對其初值進行優(yōu)化。

        經(jīng)凸優(yōu)化重構(gòu)后的矩陣()可看作由2-+1個陣元組成的均勻陣列的協(xié)方差矩陣。為提高算法估計精度,對重構(gòu)后的Toeplitz矩陣()進行前后向空間平滑分成個子陣列。為避免虛擬陣列內(nèi)插零陣元造成協(xié)方差秩虧損,滑窗長度須不小于連續(xù)差聯(lián)合子陣正向長度(+-1)?;伴L度設(shè)置為=+-1。

        FOSS處理如圖3所示。

        圖3 FOSSFig.3 FOSS

        前向平滑時,定義第個子陣列接收數(shù)據(jù)為

        (9)

        第個子陣接收數(shù)據(jù)協(xié)方差矩陣可表示為

        (10)

        FOSS后的協(xié)方差矩陣可通過各子陣的協(xié)方差矩陣求平均而得

        (11)

        (12)

        式中:為置換矩陣,其反對角線元素均為1,其余位置全為0。

        將FOSS矩陣及后向空間平滑矩陣取平均得到FBSS矩陣((+-1)×(+-1));

        (13)

        由于除2以外的質(zhì)數(shù)均為奇數(shù),當=2時,為奇數(shù),則+-1為奇數(shù);當>2時,,均為奇數(shù),則+-1為奇數(shù)。即對于互質(zhì)陣列,子陣陣元數(shù)<時,連續(xù)位置+-1均為奇數(shù)。則構(gòu)建的Toeplitz矩陣為奇數(shù)階。

        Toeplitz矩陣可以看作陣元位置為{-,…,0,…,}均勻線性陣列的輸出協(xié)方差矩陣:

        (14)

        式中:=(-1)2。

        對于相干信號,又有

        (15)

        式中:

        (16)

        2 基于后向平滑的相干信號估計

        2.1 重構(gòu)后向平滑Toeplitz矩陣

        選取的第行:

        (17)

        式中:=[,1,…,,],且信號協(xié)方差矩陣中不會存在全零行,因此中不存在零元素。利用長度為+1的滑窗截取的元素,得到行向量:

        ,=[(,-),(,1-),…,(,-)], 0≤≤

        (18)

        取第一個行向量,0:

        (19)

        式中:

        利用向量重構(gòu)Toeplitz矩陣:

        (20)

        (21)

        式中:+1,表示第條對角線元素為1,其余元素均為零的+1階方陣;=diag{}表示矩陣重構(gòu)后的虛擬信號協(xié)方差矩陣。

        忽略噪聲影響,可以寫為

        (22)

        2.2 凸優(yōu)化去噪處理

        對于第個信號,必存在1×(+1)和其余-1個導(dǎo)向矢量張成的信號子空間正交:

        (23)

        式中:range{·}表示張成的子空間。

        (24)

        將式(24)代入式(22)得

        (25)

        根據(jù)正交性質(zhì),有

        (26)

        考慮到噪聲影響,式(26)可表示為

        (27)

        為減小誤差,利用每次平滑的子陣列構(gòu)造代價函數(shù),得到凸優(yōu)化模型:

        (28)

        (29)

        構(gòu)建映射矩陣:

        (30)

        (31)

        (32)

        (33)

        式(33)的根為

        (34)

        式中:(·)表示偽逆。

        將式(34)代入式(28)得

        (35)

        ()()進行特征值分解:

        (36)

        (37)

        ()=+1-max eig{()()}

        (38)

        式中:max eig{·}表示特征值分解后最大特征值。

        構(gòu)建譜函數(shù)():

        (39)

        對()進行譜搜索,對每個()()求其最大特征值,()譜峰對應(yīng)角度即為DOA估計方向。

        算法步驟歸納如下:

        根據(jù)式(2)計算陣列數(shù)據(jù)協(xié)方差矩陣=E{()()};

        根據(jù)式(4)對虛擬陣列元素進行整合,根據(jù)式(5)構(gòu)建矩陣;

        利用CVX工具箱,根據(jù)式(7)重構(gòu);

        矩陣進行后向平滑矩陣。用+1長度的滑窗截取,根據(jù)式(17)得到;

        根據(jù)式(30)、式(31)計算映射矩陣();

        根據(jù)式(39)計算空間譜函數(shù)(),通過譜峰搜索得到信號DOA。

        2.3 算法DOF分析

        由2+-1個陣元組成的互質(zhì)陣列經(jīng)Toeplitz重構(gòu)產(chǎn)生(),因平滑子陣陣元數(shù)目大于等于相干信號數(shù)時,可以有效解相干。則對()進行空間平滑時,>。且由于對互質(zhì)陣虛擬陣列進行了空洞填充處理,為保證信號信息的完整,滑窗長度不小于虛擬連續(xù)陣列正向長度(+-1),即≥+-1。則取滑窗長度=+-1,經(jīng)前后向空間平滑預(yù)處理后,生成矩陣(2+1)×(2+1),其中=(+-2)2。

        矩陣不存在秩虧損,對矩陣第行進行空間平滑處理生成,此時滑窗長度只需不小于信源數(shù)即可完成對的解相干?;伴L度取最小值,即假設(shè)噪聲子空間僅由一條導(dǎo)向矢量構(gòu)成,滑窗長度為+1。

        綜上,2+-1個陣元組成的互質(zhì)陣列,其虛擬陣經(jīng)空間平滑后最多可以完成(+)2個相干角度估計。

        3 實驗分析

        本節(jié)通過仿真實驗將CASS算法與文獻[15]提出的FBSS算法、文獻[18-20]提出的ML算法、文獻[28]提出的FBMT算法、文獻[29]提出的SWSS算法進行比較,對不同條件下的均方根誤差(root mean square error, RMSE)進行分析。RMSE可以通過500次蒙特卡羅實驗得到:

        (40)

        式中:為信源數(shù)。

        設(shè)信號為零均值高斯信號,噪聲為零均值高斯白噪聲。由于FBSS算法、ML算法、SWSS算法需要信號源數(shù)先驗信息,而FBMT算法和CASS算法無需信號源數(shù)。故仿真實驗中,FBSS算法、ML算法、SWSS算法的信號源數(shù)是已知的,FBMT算法和CASS算法信源數(shù)先驗信息均未知。

        為方便DOF的比較,CASS算法構(gòu)建互質(zhì)陣列子陣陣元數(shù)分別為2=6,=5,總陣元數(shù)為10。作為不損失陣列孔徑算法對照,ML算法采取10元均勻線陣。由于CASS算法、ML算法可以實現(xiàn)相干信號的最大測向數(shù)為9。FBSS算法、FBMT算法、SWSS算法構(gòu)建的均勻陣列陣元設(shè)置為17,即完成9個相干信號DOA估計的FMBT算法所需要的最小陣元數(shù)。

        3.1 全相干信號的估計性能分析

        3.1.1 不同信噪比條件下的算法性能

        假設(shè)9個入射信號均相干,入射角度分別為[-66°,-42°,-30°,-15°, 0°, 10°, 25°, 42°, 56°]。相干系數(shù)分別為=[1, 02ejπ3, ejπ4, 04ejπ6, ejπ5, 06ej2π5, ejπ7, 08ej2π7, ej3π7],正則化參數(shù)=1.5×10,快拍數(shù)為500,掃描間隔為0.1°。全相干信號中信噪比對算法估計精度的影響如圖4所示。分析圖4可知, FMBT算法由于已經(jīng)達到其最大信號測量角度極限,高信噪比條件下測向性能與FBSS算法接近。ML算法由于測向信號較多,難以保證經(jīng)多次角度代入處理收斂到的局部最小值是真實信號值。SWSS算法由于利用信號協(xié)方差信息,在低質(zhì)量信號中提取出有效信息能力較強,估計精度較高。在信噪比為-12~-4 dB條件下,CASS算法由于對原互質(zhì)陣列進行同波程差數(shù)據(jù)整合以及前后向空間平滑,增強了信號強度,估計精度高。CASS算法由于填補虛擬陣列空洞,存在信息缺失,故在8 dB信噪比條件下估計性能略差。

        圖4 全相干信號信噪比對不同算法估計精度的影響Fig.4 Influence of signal to noise ratio of fully coherent signal on the estimation accuracy of different algorithms

        全相干信號-12 dB條件下各算法的歸一化功率譜如圖5所示。分析圖5可知,在低信噪比條件下,FBSS算法由于沒有進行完全解相干,估計精度較低;ML算法在對多個相干信號進行DOA估計時,需要選定合適的初值,DOA估計精度受初值設(shè)置影響。FBMT算法沒有對信號信息進行預(yù)處理,噪聲影響較大,不能準確估計DOA;SWSS算法在較大噪聲干擾下,無法通過加權(quán)子空間算法估計準確的DOA。CASS算法經(jīng)過預(yù)處理,去噪相對徹底,算法精度較高。

        圖5 全相干信號-12 dB條件下的歸一化功率譜Fig.5 Normalized power spectrum of fully coherent signal at -12 dB

        3.1.2 不同快拍數(shù)條件下的算法性能

        設(shè)置信噪比為0 dB,其余仿真環(huán)境同第3.1.1節(jié)。全相干信號中快拍數(shù)對算法估計精度的影響如圖6所示。

        圖6 全相干信號快拍數(shù)對不同算法估計精度的影響Fig.6 Influence of snapshots number of fully coherent signal on the estimation accuracy of different algorithms

        分析圖6可知,FBSS算法估計性能受快拍數(shù)影響較大,且隨著快拍數(shù)增加,估計精度提升明顯。在快拍數(shù)為20時,ML算法、FBMT算法、SWSS算法以及CASS算法估計精度遠高于FBSS算法。且CASS算法在快拍數(shù)為20時依舊保持良好的算法精度,估計性能最佳。

        3.2 部分相干信號的估計性能分析

        3.2.1 不同信噪比條件下的算法性能

        假設(shè)入射角度為[-66°, -42°, -30°, -15°, 0°, 10°, 25°, 42°, 56°]的9個信號,前6個信號為相干信號,后3個信號是與其他信號完全不相干的信號。正則化參數(shù)τ=1.5×10,快拍數(shù)為500,掃描間隔為0.1°。部分相干信號中信噪比對算法估計精度的影響如圖7所示。

        圖7 部分相干信號信噪比對不同算法估計精度的影響Fig.7 Influence of signal to noise ratio of partially coherent signal on the estimation accuracy of different algorithms

        分析圖7可知,對于部分相干信號,FBSS算法在低信噪比條件下估計精度依舊較低。ML算法、FBMT算法與SWSS算法估計精度均低于CASS算法。CASS算法在低信噪比條件下估計角度均方誤差較小。

        部分相干信號-12 dB條件下各算法的歸一化功率譜如圖8所示。

        圖8 部分相干信號-12 dB條件下的歸一化功率譜Fig.8 Normalized power spectrum of partially coherent signal at -12 dB

        分析圖8可知,在部分相干信號中,FBSS波峰尖銳,完成所有DOA估計,但與真實角度偏差較大。 ML算法經(jīng)多維非線性搜索,收斂值不能保證為全局最佳。FBMT算法在已達到其相干信源數(shù)估計極限的情況下,功率譜函數(shù)波峰較平;SWSS算法經(jīng)自加權(quán)處理,在低信噪比條件下能夠有效去噪,波峰較尖銳。CASS算法由于凸優(yōu)化去噪相對徹底,其混合信號估計性能與全相干信號估計性能相對穩(wěn)定。

        3.2.2 不同快拍數(shù)條件下的算法性能

        設(shè)置信噪比為0 dB,其余仿真環(huán)境同第3.2.1節(jié)。部分相干信號中快拍數(shù)對算法估計精度的影響如圖9所示。

        圖9 部分相干信號快拍數(shù)對算法估計精度的影響Fig.9 Influence of snapshots number of partially coherent signal on the estimation accuracy of algorithm

        分析圖9可知,FBSS算法估計精度隨快拍數(shù)增加而提高。ML算法、FBMT算法、SWSS算法、CASS算法隨快拍數(shù)增加估計精度變化曲線較為平穩(wěn)。CASS算法估計性能最佳,且在低快拍條件下性能穩(wěn)定。

        3.3 運算時間比較

        實驗仿真環(huán)境為Matlab R2020a平臺,Inter i7-9570H處理器,16 G內(nèi)存。仿真條件同第3.1.2節(jié)。各算法運算時間比較如圖10所示。

        圖10 不同算法運算時間比較Fig.10 Comparison of operation time of different algorithms

        分析圖10可知, FBSS算法只經(jīng)過空間平滑處理,運算速度最快;ML算法由于需進行多維搜索,且DOA估計數(shù)目愈多,運算時間愈長;FBMT算法由于進行多次前后向空間平滑處理及特征分解處理,運算時間也相對較長。SWSS算法由于需要對信號協(xié)方差進行加權(quán)累加,隨著精度要求的提高,迭代次數(shù)變多,計算復(fù)雜度增高;CASS算法采用互質(zhì)陣列作為接收陣列,陣元數(shù)較少,且相較FBMT算法省去多次前向平滑處理以及映射陣列的特征值分解,計算復(fù)雜度較低,運算時間相對較短。

        4 結(jié) 論

        針對現(xiàn)有基于完全相干信號的DOA估計算法陣列孔徑損失嚴重,且在低信噪比環(huán)境下DOA計算精度較低的問題,提出了一種基于互質(zhì)陣虛擬陣列空間平滑的DOA估計算法。該算法將互質(zhì)陣列進行空洞填充,協(xié)方差矩陣按照波程差大小重構(gòu)Toeplitz矩陣,然后再通過后向空間平滑算法二次構(gòu)建Toeplitz矩陣,利用構(gòu)建的與信號子空間最大信號特征值相關(guān)的代價函數(shù)來搜索譜函數(shù),得到DOA。該算法優(yōu)點如下:

        (1) CASS算法在比FBSS、FBMT、SWSS算法少陣元的條件下完成了比以上算法更準確的相干信號的DOA估計。同時,CASS算法估計精度又高于同所需陣元數(shù)相同的ML算法。CASS算法系統(tǒng)架置成本低、部署靈活,且計算復(fù)雜度低、性價比高;

        (2) CASS算法在低快拍數(shù)、低信噪比環(huán)境下對完全相干、部分相干信號DOA估計精度更高,具有更好的魯棒性;

        (3) CASS算法估計相干信號DOA時,不需要信源數(shù)先驗信息,更符合實際應(yīng)用環(huán)境。

        由于互質(zhì)陣列虛擬陣列法在信號分離時利用空洞填充,存在一定信號損失,且互質(zhì)陣列虛擬陣列算法是在虛擬域內(nèi)進行的二階等價虛擬陣列信號處理,故該算法在較高信噪比環(huán)境下DOA估計精度遜色于FBMT算法、SWSS算法。高信噪比條件下CASS算法的DOA估計精度仍存在提升空間。

        猜你喜歡
        互質(zhì)協(xié)方差信噪比
        基于互質(zhì)陣列的信號波達方向估計算法
        航空兵器(2023年2期)2023-06-25 03:04:39
        基于深度學(xué)習(xí)的無人機數(shù)據(jù)鏈信噪比估計算法
        低信噪比下LFMCW信號調(diào)頻參數(shù)估計
        電子測試(2018年11期)2018-06-26 05:56:02
        低信噪比下基于Hough變換的前視陣列SAR稀疏三維成像
        Short-range Radar Detection with(M,N)-Coprime Array Configurations
        不確定系統(tǒng)改進的魯棒協(xié)方差交叉融合穩(wěn)態(tài)Kalman預(yù)報器
        一種基于廣義協(xié)方差矩陣的欠定盲辨識方法
        保持信噪比的相位分解反褶積方法研究
        縱向數(shù)據(jù)分析中使用滑動平均Cholesky分解對回歸均值和協(xié)方差矩陣進行同時半?yún)?shù)建模
        關(guān)于協(xié)方差的U統(tǒng)計量檢驗法
        亚洲第一页综合av免费在线观看| 少妇特黄a一区二区三区| 巨大欧美黑人xxxxbbbb| 少妇高潮无码自拍| 国产一区二区三区 在线观看| 欧美丰满老熟妇aaaa片| 亚洲一区精品无码色成人| 午夜福利影院不卡影院| 亚洲av本道一本二本三区| 无码人妻一区二区三区免费看| 狠狠色综合网站久久久久久久| 久久久久久无中无码| 黄色潮片三级三级三级免费| 国产精品白浆在线观看免费| 精品熟女日韩中文十区| 亚洲欧美久久婷婷爱综合一区天堂| 国产夫妻自偷自拍第一页| 亚洲精品久久激情国产片| 台湾佬娱乐中文22vvvv| 爆乳无码AV国内| 日本高清视频在线观看一区二区| 欧美成人猛交69| 日本VA欧美VA精品发布| 国产日韩精品视频一区二区三区| 人妻少妇哀求别拔出来| 欧美操逼视频| 国产黑色丝袜在线观看视频| 国产一区二区三区不卡视频| 免费无码高潮流白浆视频| 欧美中文字幕在线| 中文字幕人妻少妇精品| 日韩人妻不卡一区二区三区| 日韩亚洲av无码一区二区不卡| 小13箩利洗澡无码免费视频| 加勒比东京热一区二区| 国语对白做受xxxxx在线| 中文字幕av一区二区三区| 伊人婷婷综合缴情亚洲五月| 无码人妻精品一区二区三区9厂 | 亚洲国产精品国自拍av| 特级a欧美做爰片第一次|