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

二維數字陣列三階泰勒展開單脈沖角度估計

2016-08-15 08:09:36馬曉峰馮丹萍盛衛星韓玉兵張仁李
系統工程與電子技術 2016年8期
關鍵詞:方向方法

馬曉峰, 馮丹萍,2, 盛衛星, 韓玉兵, 張仁李

(1. 南京理工大學電子工程與光電技術學院, 江蘇 南京 210094;2. 中航工業雷華電子技術研究所,江蘇 無錫 214063)

?

二維數字陣列三階泰勒展開單脈沖角度估計

馬曉峰1, 馮丹萍1,2, 盛衛星1, 韓玉兵1, 張仁李1

(1. 南京理工大學電子工程與光電技術學院, 江蘇 南京 210094;2. 中航工業雷華電子技術研究所,江蘇 無錫 214063)

全數字或相控陣子陣結構的數字單脈沖測角一般都采用基于單脈沖比一階泰勒展開的測角公式,主瓣干擾下的自適應數字單脈沖測角也采用鑒角曲線斜率約束的方法,這些都會導致和波束主瓣寬度內偏離波束指向較遠目標的測角偏差增大。為了解決上述問題,針對任意二維數字陣列結構,推導了基于單脈沖比三階泰勒展開的二維單脈沖測角公式,分析了公式求解的收斂問題。仿真結果表明所提方法對和波束主瓣范圍內的目標均有較小的測角偏差,可以有效地提高主瓣干擾下的測角性能。該方法在提高測角性能的同時運算量增加不大,適合工程實現。

二維數字陣列; 單脈沖角度估計; 泰勒展開; 測角偏差; 主瓣干擾

0 引 言

和差單脈沖技術能夠精確地測量目標的角度信息,因此被廣泛應用于預警和跟蹤雷達等系統。全數字或相控陣子陣結構數字陣列采用的數字單脈沖測角方法一般都是基于單脈沖比一階泰勒展開(下面簡稱一階展開)的[1-2],鑒角曲線用一定斜率的直線逼近,由于理想鑒角曲線在偏離波束指向角度較遠處存在較大的非線性,這些位置的測角偏差較大。為了擴大線性測角范圍,并減小存在主瓣干擾的條件下,自適應單脈沖鑒角曲線的非線性,有文獻提出在偏離波束指向處采用三點或者多點斜率約束的方法優化差波束方向圖[3-5],但是額外增加的約束會減少自適應算法的自由度,導致自適應差波束輸出的信干噪比下降,可能改變差波束方向圖的形狀,使得差波束的副瓣電平惡化,而且該方法擴大目標線性測角區間的范圍也比較有限。針對鑒角曲線不精確引起的測角偏差問題,傳統的單脈沖測角系統也有采用鑒角曲線直接存儲或者近似的解決方法[6],這種方法對于指向精度不高的一維單脈沖測角系統比較適用,而高指向精度的二維大陣列單脈沖測角系統,波束指向為二維角度,且每個角度均需要存儲理想的鑒角曲面,存儲數據量大,靈活性比較差[7]。

文獻[8]針對一維數字陣列結構,提出了采用單脈沖比三階泰勒展開(下面簡稱三階展開)的自適應單脈沖測角方法,使得一維鑒角曲線以三次多項式的形式逼近理想鑒角曲線,和波束主瓣寬度內,在目標偏離波束指向較遠處仍保持較小的測角偏差。還有文獻給出單脈沖比更高階泰勒展開的自適應單脈沖測角方法[9]。本文針對任意二維數字陣列結構,詳細推導了基于三階展開的二維單脈沖測角公式,考慮到陣列結構的對稱性,簡化測角公式,接著,采用牛頓迭代法計算角度估計值,并分析了采用牛頓迭代法僅需要幾次迭代就可以獲得精確角度估計值的原因。仿真結果表明,基于三階展開的二維單脈沖角度估計方法,在和波束主瓣寬度內具有均勻的測角偏差,測角偏差優于一階展開的方法;該方法還十分適用于主瓣干擾情況下,全數字或相控陣子陣結構數字陣列的自適應單脈沖角度估計。

1 問題模型

1.1陣列模型與測角公式

如圖1所示,假設平面陣列天線的N個各向同性的單元天線位于極坐標系的yoz平面內,第i個單元天線的位置為ri=(0,yi,zi)T。平面波來波的單位方向性矢量為U=(w,u,v), 其中,w=sinθcosφ,u=sinθsinφ,v=cosθ,θ為與z軸夾角,φ為xoy平面內與x軸夾角。x軸(θ,φ)=(90°,0°)為波束指向法線方向,目標角度估計范圍為θ:0°~180°,φ:-90°~90°。則第i個單元天線接收到的復基帶信號為

(1)

式中,b為回波復幅度;λ為工作波長;ni為第i個單元天線輸出的復高斯白噪聲信號。回波信號復向量可以表示為

(2)

式中,a(u,v)為(u,v)方向的導向性矢量,ai(u,v)=ej(2π/λ)(yiu+ziv)(i=1,…,N);復高斯白噪聲向量n的均值為0,方差為σ2I;b的最大似然估計可以通過最大化如下的高斯概率密度函數獲得

(3)

因此,b的最大似然解為

(4)

式中,(·)H為共軛轉置;U的最大似然估計值為式(5)給出的和波束掃描功率對數方向圖的最大值對應的角度[1-2]。

(5)

式中,w(u,v)為和波束的權重系數。當和波束為靜態波束時,w(u,v)=a(u,v);當和波束為自適應波束時,w(u,v)=μR-1a(u,v),μ=1/aH(u,v)R-1a(u,v)為一個標量系數,R-1為陣列接收復基帶信號的協方差矩陣估計的逆。

圖1 均勻平面陣列結構示意圖

記 F(u,v)的一階偏導數為Fα:

(6)

式中,Fu,Fv就是u,v方向的單脈沖比;wu,wv為和波束權重系數在u,v方向的偏導數,即u,v方向差波束的權重系數。

假設,目標方向(us, vs)在陣列天線和波束指向(u0, v0)附近,將單脈沖比Fα在(u0, v0)處進行三階泰勒展開,可得

(7)

式中,Fαβ,Fαβγ,Fαβγδ(α,β,γ,δ可取值為u,v)分別為和波束對數功率方向圖F(u,v)的二階、三階和四階偏導數。通過求解公式(7)就可以得到目標方向(us, vs)的精確估計,本論文第三部分將推導公式(7)所有系數的計算方法,給出采用牛頓迭代法求解目標角度估計值的過程,并分析牛頓迭代法的收斂性。

1.2二維角度估計的克拉美羅界

為了對本文提出方法的測角性能進行有效評估,下面給出二維角度估計的克拉美羅界。克拉美羅界確定了任何無偏估計量的最小方差。

以估計u值為例,克拉美羅不等式為

(8)

式中

(9)

將式(9)代入式(8),可以得到測量所得u偏離實際角度us的方差的下界。

同理可以得到v偏離實際角度vs的方差的下界。

2 單脈沖比三階泰勒展開方法

由式(6)給出的和波束對數功率方向圖F(u,v)的一階偏導數Fα的計算公式可以知道,當忽略接收機的噪聲分量,即z≈ba(u,v)時,若設波束指向角度的導向矢量為權矢量,則Fα在(u0,v0)處的值為0。文獻[10]指出,對于高增益大規模陣列天線(如902陣元的陣列),當陣元信噪比低于-10dB時,忽略噪聲分量會有較大的影響,但是在太低的信噪比條件下,任何單脈沖角度估計方法都是不可用的。因此,我們僅考慮實際單脈沖角度估計可以使用的場合,比如經過積累處理后的和波束輸出信噪比大于10dB的情況,此時可以忽略噪聲分量。在后續推導F(u,v)的二階、三階和四階偏導數Fαβ,Fαβγ,Fαβγδ(α,β,γ,δ可取值為u,v)時,都將忽略接收機噪聲的影響。推導得到的各項系數值都是在波束指向角度(u0,v0)處的值。

(10)

h=wHzzHw

(11)

式中,g[α](α=u,v)表示α方向的差波束功率。

則F(u,v)的二階偏導數為

(12)

F(u,v)的三階偏導數為

(13)

四階偏導數可以表示為

(14)

其中

(15)

(16)

(17)

(18)

式(15)~(18)中

(19)

(20)

(21)

(22)

(23)

(24)

為了后面計算方便,令

(25)

(26)

將計算所得系數代入式(7)的三階展開測角公式中,解二元三次方程就可以計算目標方向(us,vs)。

需要指出,式(25)和式(26)中的各個值,在波束指向角度(u0,v0)處的取值僅與陣元的坐標位置ri=(0,yi,zi)(i=1,2,…,N)有關,當二維陣列陣元位置中心對稱時(一般都能滿足),可以得到如下結論

(27)

(28)

(29)

(30)

式中,α,β,γ,δ=u,v,且u的個數為奇數;o,p,q,s=x,y。

(31)

(32)

(33)

(34)

式中,α,β,γ,δ=u,v,且u的個數為奇數;o,p,q,s=x,y。

因此,將式(27)~式(34)代入式(19)~式(24)后可以看到,F(u,v)的二階、三階和四階偏導數中,除以下系數外其余系數均為0。

(35)

(36)

(37)

(38)

(39)

(40)

所以,基于三階展開測角公式(7)可以簡化為

(41)

可以看到,式(41)中二次項系數全為0,因此,單脈沖比二階泰勒展開的測角結果與一階展開的結果是一樣的。

3 角度估計值求解與收斂性分析

三階展開式(41)是一個典型的非線性方程組,可以采用牛頓迭代法求解,定義函數向量

(42)

則式(41)的解即為l(us,vs)=0的解。

(43)

于是,得到牛頓迭代法的迭代公式為

(44)

迭代公式的初始值可以設定為和波束指向(u0,v0),由于目標方向(us,vs)靠近和波束指向(u0,v0),經過幾次迭代就可以收斂。

使用式(44)計算得到的第1次迭代結果為

(45)

式中

(46)

(47)

將式(46),式(47)代入式(45),得

(48)

從式(48)可以看到,第一次牛頓迭代得到的結果等價于一階展開測角公式求解得到的結果。

下面以函數l1(us,vs)為例,分析一下采用上述牛頓迭代法的收斂性。

定理 1[11]設函數f(x)滿足f(x*)=0,f′(x*)≠0且f(x)在x*的領域內有二階連續導數,則當初值x0足夠接近于x*時,由牛頓迭代公式得到的序列至少是二階收斂的。

l(us,vs)的二階導數仍為初等函數,故其在us的領域Bδ(us)=[us-δ,us+δ]是連續的,并且

(49)

(50)

則以任意初值us0∈Bδ(us),由牛頓迭代法得到的序列至少二階收斂于us。

由于牛頓迭代法的局部收斂性,故初值的選取非常重要。在解式(41)時,我們取和波束指向方向(u0,v0)為初值,第一次牛頓迭代的結果與一階展開的測角結果完全相同,已經很接近實際值,所以牛頓迭代法能快速收斂。

4 仿真結果分析

4.1無干擾條件下測角性能

仿真設置的陣元信噪比為20dB,噪聲是均值為0,方差為1的高斯白噪聲。

仿真 1假設和波束中心指向為(θ0,φ0)=(60°,0°),即(u0,v0)=(0,0.5),圖2(a)和圖2(b)分別給出了目標處于和波束3dB波束寬度內及其附近不同位置,三階展開和一階展開的測角性能。箭頭的起點為目標真實角度,每個角度進行500次蒙特卡羅仿真,每次仿真得到的測角結果在圖中用小圓圈標出,箭頭的終點為500次角度估計結果的均值。從小圓圈的分布范圍可以看出測角方法的測角精度,而從箭頭的長短可以看出測角方法的測角偏差。

圖2 和波束指向(u0, v0)=(0, 0.5)時的測角性能比較

從圖2可以看出,和波束波束寬度范圍內,當目標角度偏離和波束指向角度較遠時,三階展開具有更小的測角偏差,這是因為三階展開測角方法的鑒角曲面采用三次多項式的形式逼近理想的鑒角曲面,當目標偏離波束指向角度較大時仍具有較好的近似。圖3給出了u方向的鑒角曲面切圖,當目標角度偏離波束指向角度較遠(仍在主瓣3 dB寬度范圍內)時,一階展開鑒角曲線已經偏離了理想鑒角曲線,而三階展開仍能較好的逼近理想情況,因此測角性能更好。

另外,從每個角度500次測角結果的分布來看,測角精度差別不大,u方向或者v方向的測角精度可由如下計算公式估算得到[12]

(51)

式中,k為常數,只與陣元結構有關,u3 dB為和波束3 dB波束寬度,SNRout為和波束輸出信噪比(signal to noise ratio,SNR)。兩種測角方法在和波束指向附近的500次測角方差與式(51)吻合,差別不大,但一階展開方法在偏離和波束指向較遠處的測角方差逐漸增大,主要還是由于鑒角曲面誤差卷入引起的。

圖3 和波束指向(u0, v0)=(0, 0.5),固定 v=0.5時, u方向的鑒角曲線 

為了更清楚地比較一階展開和三階展開測角性能差異,表1給出了和波束3 dB波束寬度內及其附近的部分角度位置,500次統計平均得到的實測角度和實際角度夾角。

表1 無干擾情況下一階展開和三階展開測角結果

從表1中可以看到,目標離波束指向很近時,三階展開和一階展開兩種方法測角偏差均值相當,隨著目標偏離波束指向,三階展開的測角偏差均值明顯小于一階展開。

仿真 2為了比較一階展開和三階展開在整個測角區間的測角性能,我們設置波束指向方位角φ=0°,俯仰角θ在45°~135°范圍內以1°為步進掃描,計算每個波束指向角度下,目標處于和波束3 dB波束寬度內及其附近任意位置的平均測角偏差,如圖4所示。

圖4 φ=0°時θ方向的平均測角誤差比較

從圖4可以看到,兩種方法在θ變化時測角性能變化不大,三階展開的測角偏差在整個測角區間始終比一階展開好,更逼近克拉美羅界。

4.2干擾條件下測角性能

仍然以圖1所示二維陣列為例,給出三階展開測角方法應用于自適應單脈沖系統的測角性能。仿真條件為,陣元信噪比20 dB,陣元干噪比40 dB,干擾和信號互不相關,噪聲為均值為0,方差為1的高斯白噪聲。和波束與差波束的自適應權重系數均采用線性約束最小方差準則優化得到[10]。

仿真 3旁瓣干擾條件下測角性能仿真。和波束中心指向為(θ0,φ0)=(60°,0°),即(u0,v0)=(0,0.5);干擾方向角度(θJ,φJ)=(60°,35.26°),即(uJ,vJ)=(0.5,0.5),圖5給出了v=0.5時,u方向的鑒角曲面切圖。可以看到,存在旁瓣干擾情況下,三階展開鑒角曲線可以更好的逼近理想鑒角曲線,由于旁瓣干擾的抑制對主瓣區域的形狀影響較小,鑒角曲線與無干擾情況下的鑒角曲線接近,因此角度估計精度和偏差也與無干擾情況下的類似,此處不再給出仿真結果。

圖5 和波束指向(u0,v0)=(0,0.5),旁瓣干擾方向(uJ,vJ)=(0.5,0.5),固定v=0.5時,u方向的鑒角曲線 

仿真 4主瓣干擾條件下測角性能仿真。和波束中心指向為(θ0,φ0)=(60°,0°),即(u0,v0)=(0,0.5);干擾方向角度(θJ,φJ)=(60°,13.35°),即(uJ,vJ)=(0.2,0.5),圖6給出了v=0.5時,u方向的鑒角曲面切圖。可以看到存在主瓣干擾的情況下,由于和波束及差波束方向圖在主瓣區域,特別是在干擾角度附近嚴重變形,鑒角曲線發生較大的畸變,一階展開方法逼近理想鑒角曲線的誤差變大,而三階展開方法除了干擾角度附近外逼近仍然良好。

圖6 和波束指向(u0, v0)=(0, 0.5),主瓣干擾方向(uJ, vJ) (0.2, 0.5),固定v=0.5時,u方向的鑒角曲線 Fig.6 Monopulse ratio of u when sum beam points to (u0, v0)=(0, 0.5), v=0.5 and the mainlobe interface located at (uJ, vJ)=(0.2, 0.5)

圖7(a)和圖7(b)分別給出了目標在和波束3 dB波束寬度范圍內及其附近時,三階展開和一階展開的測角性能。每個角度同樣進行500次蒙特卡羅仿真。

圖7 和波束指向(u0, v0)=(0, 0.5),主瓣干擾方向(uJ, vJ) (0.2, 0.5)時的測角性能比較        Fig.7 Comparison of angle estimation performace when sum beam points to (u0, v0)=(0, 0.5) and the mainlobe interface located at (uJ, vJ)=(0.2, 0.5)

從圖7可以看出,當存在主瓣干擾時,三階展開方法測角性能保持較好,僅在干擾角度附近測角偏差和測角精度有所降低,而一階展開方法在主瓣區域內偏離和波束指向角度,特別是主瓣干擾附近角度,測角偏差和測角精度都明顯降低,測角穩健性較差。

同樣,表2給出了和波束3 dB波束寬度內及其附近的部分角度位置,500次統計平均得到的實測角度和實際角度夾角。

表2 主瓣干擾情況下一階展開和三階展開測角結果

從表2中可以看到,當存在主瓣干擾時,三階展開的測角偏差更小。另外從表1和表2的測角結果可以看到,當存在干擾時,兩種方法各自的測角誤差都有所增大。

5 結 論

本文針對任意二維數字陣列結構,推導了基于單脈沖比三階泰勒展開的二維單脈沖測角簡化公式及其快速迭代求解算法,并分析了迭代算法的收斂性。采用該方法的單脈沖鑒角曲面與理想鑒角曲面逼近程度優于單脈沖比一階泰勒展開的方法,因此針對和波束主波束寬度范圍內的目標均具有明顯優于單脈沖比一階泰勒展開的方法的測角偏差。同時,該方法還十分適用于提高主瓣干擾情況下自適應單脈沖算法的角度估計精度及其穩健性。

[1] Nickel U. Overview of generalized monopulse estimation[J].IEEEAerospaceandElectronicSystemsMagazine,2006,21(6):27-55.

[2] Nickel U. Monopulse estimation with adaptive arrays[J].IEEProceedings-Radar,SonarandNavigation, 1993, 140 (5): 303-308.

[3] Fante R L. Synthesis of adaptive monopulse patterns[J].IEEETrans.onAntennasandPropagation, 1999, 47(5): 773-774.

[4] Yu K B, David J M. Adaptive digital beamforming for angle estimation in jamming[J].IEEETrans.onAerospaceandElectronicSystems, 2001, 37(2): 508-523.

[5] Li R F, Rao C, Dai L Y, et al. Algorithm for constrained adaptive sum-difference monopulse among sub-arrays[J].JournalHuazhongUniversityofScienceandTechnolgy, 2013, 41(9): 6-10. (李榮鋒, 饒燦, 戴凌燕, 等. 子陣間約束自適應和差單脈沖測角算法[J].華中科技大學學報, 2013, 41(9): 6-10. )

[6] Niu B J, Li Y B. Study and application of angle measurement for 2-D phased array monopluse tracking[J].ModernRadar, 2003, 25(5): 16-18.(牛寶君, 李延波. 二維相控陣單脈沖跟蹤測角方法的研究與應用[J].現代雷達, 2003, 25(5): 16-18.)

[7] Zhu W, Chen B X, Zhou Q. Angle measurement method with digital monopulse for 2-dimensional digital array radar[J].SystemsEngineeringandElectronics, 2011, 33(7): 1503-1509.(朱偉, 陳伯孝, 周琦. 兩維數字陣列雷達的數字單脈沖測角方法[J].系統工程與電子技術, 2011, 33(7): 1503-1509.)

[8] Chen C Z, Wang T,Wu J X. Monopulse angle estimation with adaptive array based on the third-order taylor series[J].ModernRadar,2013,35(8):32-36.(陳成增,王彤,吳建新.基于三階泰勒展開的自適應單脈沖測角方法[J].現代雷達,2013,35(8):32-36.)

[9] Chen C Z, Wang T, Wu J X, et al. Monopulse angle estimation with adaptive array based on monopulse response curve fitting[J].SystemsEngineeringandElectronics, 2013, 35(7): 1404-1408.(陳成增,王彤,吳建新,等.基于鑒角曲線擬合的自適應單脈沖測角方法[J].系統工程與電子技術,2013,35(7):1404-1408.)

[10] Nickel U. Monopulse estimation with subarray-adaptive arrays and arbitrary sum and difference beams[J].IEEProceedings-Radar,SonarandNavigation, 1996, 143(4): 232-238.

[11] Guan Z, Lu J P.Basicofnumericalanalysis[M]. Beijing: Higher Education Press, 2002. (關治, 陸金甫. 數值分析基礎[M]. 北京:高等教育出版社, 2002.)

[12] Xu Z H, Huang Y G, Xiong Z Y, et al. On the consistency of mono-pulse and maximal likelihood estimation with array radar[J].ModernRadar, 2013, 35(10): 32-35.(徐振海, 黃艷剛, 熊子源,等. 陣列雷達單脈沖與極大似然估計的一致性[J].現代雷達, 2013, 35(10): 32-35.)

Third-order taylor expansion based monopulse angle estimation for 2-dimension digital array

MA Xiao-feng1, FENG Dan-ping1,2, SHENG Wei-xing1, HAN Yu-bin1, ZHANG Ren-li1

(1. School of Electronic and Optical Engineering, Nanjing University of Science and Technology, Nanjing 210094,China; 2. The Radar and Avionics Institute of Aviation Industry Corporation of China, Wuxi 214063, China)

Digital monopulse angle estimation for full digital or subarray phased array system is generally based on the monopulse ratio first-order taylor expansion, while the monopulse ratio is always linearly constrained in adaptive digital monopulse angle estimation for mainlobe jamming suppression. However, above methods will cause large angle estimation deviation when the target visibly departs from the antenna look direction. In view of the above problems, the 2-dimension monopulse angle estimation formula for arbitrary 2-dimension digital array structure is presented, which is based on the monopulse ratio third-order taylor expansion. The convergence analysis of formula solution is also given. Numerical simulations illustrate that less angle estimation deviation can be obtained for the target within the mainlobe region. It also demonstrates that the target angle estimation performance is improved in the presence of mainlobe jamming. The presented method is suitable for engineering implementation for its less calculation amount.

2-dimension digital array; monopulse angle estimation; Taylor expansion; angle estimation deviation; mainlobe jamming

2015-05-22;

2016-02-21;網絡優先出版日期:2016-06-07。

國家自然科學基金(61401207)資助課題

TN 953+.5

A

10.3969/j.issn.1001-506X.2016.08.04

馬曉峰(1981-),男,講師,博士研究生,主要研究方向為陣列信號處理、MIMO雷達信號處理。

E-mail:maxiaofeng@njust.edu.cn

馮丹萍(1991-),女,助理工程師,碩士,主要研究方向為數字陣列角度估計、波束賦形。

E-mail:enny28fdp@126.com

盛衛星(1966-),男,教授,博士,主要研究方向為陣列天線和陣列信號處理、雷達目標電磁散射特性建模及應用、微波成像。

E-mail:shengwx@njust.edu.cn韓玉兵(1972-),男,副教授,博士,主要研究方向為圖像處理、陣列信號處理。

E-mail:hanyb@njust.edu.cn

張仁李(1986-),男,講師,博士,主要研究方向為雷達信號處理。

E-mail:zhangrenli_nust@163.com

網絡優先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20160607.1138.004.html

猜你喜歡
方向方法
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 欧美 亚洲 日韩 国产| 国产菊爆视频在线观看| 亚洲人成网站在线观看播放不卡| 精品无码人妻一区二区| 国产一二三区视频| 国产成人精品免费视频大全五级| 91午夜福利在线观看精品| 特级毛片8级毛片免费观看| 国产成人精品免费av| 亚洲欧州色色免费AV| 亚洲日韩Av中文字幕无码| 亚洲一欧洲中文字幕在线| 怡红院美国分院一区二区| 亚洲无限乱码一二三四区| 国产精品亚洲五月天高清| 国产欧美中文字幕| 国产麻豆精品久久一二三| 国产亚洲精品va在线| 亚洲国内精品自在自线官| 午夜福利视频一区| 啪啪国产视频| 国产手机在线ΑⅤ片无码观看| 国产成人a毛片在线| 亚洲 欧美 中文 AⅤ在线视频| 国产不卡在线看| 在线无码九区| 免费欧美一级| 中文一级毛片| 亚洲伊人久久精品影院| 国产区91| 国产玖玖玖精品视频| 国产成人高清精品免费5388| 国产精品第5页| 无码网站免费观看| 欧美日韩国产在线播放| 国产白浆在线| 国产高颜值露脸在线观看| 久久毛片网| 国产精品所毛片视频| 亚洲人成网线在线播放va| 免费人成视频在线观看网站| 国产精品无码影视久久久久久久| 国产无码在线调教| 久久这里只有精品免费| 91青青草视频| AV不卡国产在线观看| 国产精品任我爽爆在线播放6080 | 美女国内精品自产拍在线播放| 国产久草视频| 欧洲在线免费视频| 在线观看免费AV网| 久久一级电影| 一级毛片在线播放免费观看| 国产精品视频白浆免费视频| 国产手机在线小视频免费观看| 欧美中文字幕一区二区三区| 中文字幕人妻av一区二区| 国产福利2021最新在线观看| 99久久精品免费看国产免费软件| 亚洲第一页在线观看| 亚洲人成亚洲精品| 亚洲网综合| 午夜欧美理论2019理论| 亚洲精品日产精品乱码不卡| 超清人妻系列无码专区| 久久亚洲国产一区二区| 亚洲精品你懂的| 欧美在线视频a| 在线免费看黄的网站| 国产日韩精品一区在线不卡| 久久久久中文字幕精品视频| 97国产在线播放| 欧美丝袜高跟鞋一区二区| 国产男人天堂| 国产69精品久久久久孕妇大杂乱| 99久久婷婷国产综合精| 国产成人综合亚洲欧美在| 亚洲AⅤ无码国产精品| 最近最新中文字幕免费的一页| 九九九国产| 亚洲男人的天堂视频| 国产一级毛片yw|