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

虛擬波束四階累積量DOA估計方法

2022-06-25 13:04:08劉甲磊馬佳智施龍飛
系統工程與電子技術 2022年7期
關鍵詞:信號

劉甲磊, 馬佳智, 施龍飛

(國防科技大學電子科學學院, 湖南 長沙 410073)

0 引 言

信號波達方向估計(direction of arrival,DOA)是陣列信號處理的重要研究內容,廣泛應用于雷達、通信、聲吶、無源定位等各個領域[1-2]。最具代表性的子空間多重信號分類(multiple signal classification,MUSIC)及衍生的多種算法[3-4],構成了陣列信號DOA領域的廣泛應用基礎,但其本質為基于信號二階統計特性構建的譜估計算法。相較信號二階統計特性,信號四階統計特性因其對任何形式的高斯過程的不敏感性,可以最大程度地抑制高斯噪聲的突出優點,使其在盲源分離和DOA估計上獲得許多開拓性的研究[5-7]。并且由四階累積量(fourth-order cumulant,FOC)計算性質,Liu等[8]證明了其在擴展陣列孔徑和陣列校準的有效性。陳建等[9-10]提出的基于FOC的MUSIC-like陣列擴展方法,應用陣列接收信號的全FOC構建擴展陣列,進行多源信號DOA估計,可有效提高陣列分辨率。文獻[11-12]分別提出了期望最大化(expectation maximization,EM)算法迭代和相關源信號下,利用信號FOC進行DOA估計的方法。齊棟等[13]提出了基于FOC的高斯色噪聲下二維DOA估計方法。Zhang等[14]提出了高斯色噪聲下基于FOC的無網格相干DOA估計方法,表現了FOC算法在高斯色噪聲背景下的有效性。Xue等[15]將FOC用于均勻圓陣DOA估計,獲得較好的估計精度。文獻[16-17]提出了一種對MUSIC-like算法去冗余的FOC的DOA方法,只利用部分無冗余的FOC切片構建FOC矩陣,矩陣階數大大下降,利于工程化實現。但是其構建的去冗余FOC矩陣,破壞了原有矩陣的對稱性,DOA估計精度有所下降,Shi等[18]提出了一種構建Toeplitz去冗余FOC陣的方法。但由于FOC對野值非常敏感,真實情況陣列誤差未知且未被校準時,估計性能下降十分嚴重。文獻[19-21]分別建立陣列幅相和互耦誤差下基于FOC的DOA估計模型,但每種方法只能部分消除一種誤差來源,較實際多誤差因素耦合時并未提及。

針對陣列雷達實際面臨復雜電磁環境,在強干擾場景和陣列誤差并存時,低信噪比(signal to noise ratio, SNR)目標信號DOA難估計的問題,本文提出了一種基于虛擬波束FOC(virtual beam forming-FOC,VBF-FOC)的DOA估計算法。算法通過陣列接收信號主成分提取構建初始虛擬波束,再利用參數化的方向波束高斯窗對初始波束進行修正,可比較有效提取強干擾下弱目標信號分量,一定程度減小陣列誤差對接收信號失真影響,降低SNR目標的估計損失。利用以上兩步形成新的聯合導向矢量,與FOC分解的噪聲子空間進行空間譜計算,可得到精確DOA估計。本文建立了算法流程框架,并推導該算法虛擬波束形成及其譜估計的公式,在非等功率信源場景下驗證了算法有效性。仿真對比了基于陣列信號全FOC的DOA估計的MUSIC-like算法,對MUSIC-like中FOC去冗余的簡化算法(reduced FOC,ReFOC),以及利用本文方法進行虛擬波束處理后的VBF-FOC和VBF-ReFOC兩種方法。仿真結果表明,虛擬波束處理后的VBF-FOC和VBF-ReFOC兩種方法,信號DOA檢測率和估計精度都有明顯提升,對于強干擾背景下的低SNR目標,表現出更強的適應性和魯棒性。

1 算法原理

論文虛擬波束形成含義是,對陣列單脈沖接收信號進行陣元級加權,利用合成波束信號作為輸入進行空間譜估計,算法模擬了陣列波束形成空間掃描的過程,構建本不存在的波束進行DOA估計,故為虛擬波束。VBF-FOC算法首先對陣列接收信號進行歸一化的虛擬波束域提取,將N路輸入提取為K(K≤N)路初始虛擬波束并對其中主要信號分量的進行DOA粗估計,計算DOA粗估計下主波束的方向性高斯修正窗函數,利用修正高斯窗與初始虛擬波束相乘,得到修正后的虛擬波束。將主波束信號代替陣列接收信號計算FOC矩陣,進而分解出FOC下的信號子空間和噪聲子空間,聯合上步計算得到的高斯窗修正的波束方向矢量,進行FOC譜估計,搜索譜峰從而得到信號的精確DOA估計。算法流程框圖如圖1所示,各步原理及分析下文一一敘述。

1.1 信號模型

以N元均勻等距線陣模型,陣元間距為d,雷達發射信號波長為λ,陣列波束指向為θ0。遠場空間存在M(M≤N)個入射角度為θm(i=1,2,…,M)的窄帶信源sm(t),其中目標信號位于θ1,剩余M-1個為干擾信源。則理想情況下,陣元接收信號為

X(t)=AS(t)+Vn(t)

(1)

式中:Vn(t)∈N(0,σ2)為高斯白噪聲矩陣;S(t)=[s1(t),s2(t),…,sM(t)]為信源的時域信息矩陣;[A=a(θ1),A=a(θ2),…,A=a(θM)]為陣列流形矩陣,a(θM)=[1,e-jφm,…,e-j(N-1)φm,]表示其中各分量,φm=(2πdsinθm)/λ為陣元間相移因子。

但實際應用,陣列誤差對信號觀測的影響不可忽略。陣列誤差來源分別為陣列幅相誤差、互耦誤差、空間位置擾動誤差3部分,且通常同時存在且伴隨觀測始終。為更真實模擬多源非理想特性場景,下面對3種誤差因素與信號觀測影響進行分析。

陣元通道幅相特性不一致是造成陣列幅相誤差的主要原因,各陣元通道幅相特性不一致時,其幅相誤差因子矩陣為

Γ=diag[ρ1ejφ1,ρ2ejφ2,…,ρNejφN]

(2)

式中:ρn表示第n(1≤n≤N)個陣元的幅度增益;φn表示第n個陣元的相位誤差。

工程應用中各陣元間的耦合效應產生的誤差為互耦誤差,互耦矩陣在通常情況下為對稱Toeplitz矩陣。其形式為

(3)

式中:cm,n為第m個與第n個陣元間的幅度互耦因子;γm,n為第m個與第n個陣元之間的相位互耦因子。由互易定理,cm,n=cn,m,γm,n=γn,m,由于隨著陣元間距的增大,互耦效應快速減弱,通常可將C建模為部分為0的稀疏矩陣。

N元陣列M個窄帶信號下,陣元n的位置擾動為(Δxn,Δyn,Δzn),c為光速,則第m個信源在陣元n上的延時擾動為

(4)

此時空間擾動誤差矩陣為

W=[w(θ1),w(θ2),…,w(θM)],w(θm)=[e-jw0Δτ1m,e-jw0Δτ2m,…,e-jw0ΔτNm]T

(5)

相似的,當3種誤差因素并存時,陣列的導向矢量及陣列流形可以表示為

(6)

(7)

此時真實信號接收模型為

(8)

1.2 VBF-FOC算法

1.2.1 FOC子空間分解

假設獨立信源發出的信號為零均值、非高斯的信號,則復信號陣列輸出的FOC矩陣為

(9)

式中:?k1,k2,k3,k4∈{1,2,…,N}。

根據式(9),定義陣列誤差影響下的FOC切片的展開式為

(10)

式中:將不同來波角度方向源信號FOC矩陣的簡記為

(11)

(12)

(13)

實際應用中,由于不同源信號來自于不同系統,可認為各信源間相互獨立,則源信號FOC矩陣CS的一個切片為

(14)

進一步可得,CS中僅有M個元素不為零,且非零元素均位于CS的對角線上。將CS的對角線特征值按大小重排后,可得新形式為

(15)

(16)

式(15)為陣列信號FOC矩陣CX的特征分解形式,式(16)為分解特征矩陣。其中,ΛM為M個大特征值(信源)所對應的特征值矩陣,ON2-M為N2-M個小特征值(噪聲)所對應的特征值矩陣,B=[b(θ1),b(θ2),…,b(θN2)]為分解得到的特征向量構成的矩陣,可理解為FOC形成虛擬陣列的導向矢量矩陣,其中對應方向θ的導向矢量為

(17)

式中:?表示Kronecker積。相應可得到信號子空間ES和噪聲子空間EN為

ES=span{b(θ1),b(θ2),…,b(θM)}

(18)

EN=span{b(θM+1),b(θM+2),…,b(θN2)}

(19)

則基于FOC直接分解的空間譜為

(20)

1.2.2 虛擬波束FOC計算

FOC作為信號的四階統計特性,其高階的計算過程決定了其對于數據野值很敏感。當場景中存在多個非等功率源信源,信號強度差異大時,SNR大的信源,其四階統計特性會遠強于小SNR信源,這就造成了復雜場景下低SNR目標DOA估計的困難。

虛擬波束形成,是對陣列接收信號X中各陣元信號進行加權,通過主成分提取方法分離陣列信號X中的主要信息量。實際權系數矩陣求取,是對X的自協方差矩陣進行特征分解,從中選取對應主信號分量的大特征值信號,對各來波方向信號進行粗提取。對各源信號提取時,為提升低SNR目標的DOA估計性能,虛擬波束會對各信號分量進行特征值歸一化處理,形成波束的權矩陣。波束的構建方法為

RXX=XXH=UΛUH

(21)

U=[u1,u2,…,uN]

(22)

Λ=diag[λ1,λ2,…λN]

(23)

將UH前乘至陣列數據陣X上得P=UHX,計算P的協方差陣為

RPP=[UHX][UHX]H=[UHX][XHU]=UH[UΛUH]U=Λ

(24)

由P的協方差陣為對角陣,可證明P的各行互相正交,且P各行為X特征分解的信號分量,其在P中排序與特征值能量大小順序相同,其中大特征對應信號分量構成信號主分量空間,譜估計時又叫信號子空間。

為使信號各分量能量相當,通過信號分量特征值歸一化的方法,減小強干擾壓制下對低SNR目標的混疊影響。取虛擬波束變換陣W=Λ-1/2UH,將其前乘至陣列數據陣X上,得構建波束為

(25)

(26)

同式(24),可證明構建的虛擬波束Z的協方差陣為單位陣,即Z中各行互相正交,且每行數據代表的信號分量能量都為1,通過虛擬波束構建的Z中各行信號為去除了二階統計相關下的正交信號。

以虛擬波束信號為輸入,進行FOC的子空間分解,則由式(13)得Z的FOC矩陣為

(27)

(28)

式中:D=[d(θ1),d(θ2),…,d(θK2)]為虛擬波束FOC矩陣特征分解,得到的特征向量構成的矩陣,代表擴展FOC矩陣的方向導向陣;K為提取的虛擬波束個數。當K=N時,Z為同維波束陣;當K

1.2.3 高斯窗修正的虛擬波束FOC譜估計

根據FOC擴展的方向導向陣,可得對應方向θ的導向矢量為

d(θ)=w(θ)?w*(θ)

(29)

(30)

式中:w(θ)定義為虛擬波束方向圖矢量,可理解為構建虛擬波束時,各波束對來波方向為θ的信號加權的權值系數。

w(θ)中,當θ遍歷角度區間(-90°,90°)取值時,記W(θ)=[w1(θ),w2(θ),…,wK(θ)]T表示各波束的全向方向圖矢量,表明構建的虛擬波束對來波信號具有主分量提取的作用,但其并不能完全提取出各單獨源信號。殘存的其余信號分量,會在構建的主波瓣外,形成小的鼓包,波束畸變會對譜估計時DOA估計性能產生影響。

為消除波束畸變影響,本文在虛擬波束基礎上,進一步利用w(θ)中源信號角度的粗估計信息,建立高斯窗來修正虛擬波束的方向圖矢量。具體方法為:根據信號子空間中大特征值的前M路信號,求取方向圖矢量w1(θ)~wM(θ)的極大值得到DOA粗估計,表示為

(31)

(32)

式中:δ為角度高斯分布標準差,δ取值一般為陣列3 dB主瓣寬度的平方根。以陣元數N,陣元間距為半波長時,可取

(33)

相應地,記W中前M路表示大特征值波束權陣記為W1,小特征值波束權陣記為W2。則高斯窗修正后的方向圖和陣列導向矢量分別為

(34)

(35)

綜上,可得在虛擬波束提取和高斯窗修正下的聯合譜估計函數為

(36)

(37)

1.3 虛擬波束提取與高斯窗修正效果分析

為更好說明虛擬波束提取和高斯窗修正方法,對源信號DOA估計的影響,下面以信號處理效果進行分析。

實驗場景:均勻線陣陣元數N=10,陣元間距半波長,3 dB主瓣寬度10.2°。空間3個獨立點源信號S1~S3方位角分別為:0°,-30°,15°;SNR分別為:0 dB,35 dB,45 dB。其中,S1模擬低SNR的弱目標信號,S2和S3為強干擾信號。S1~S3單位SNR下的信號樣式見圖2所示。其中,S1為線性調頻(linear frequency modulation, LFM)信號,帶寬B=5 MHz,脈寬τ=20 μs,采樣頻率fs=10 MHz,S2為對S1的轉發多假目標干擾信號,S3為高斯白壓制干擾信號。

計算空間譜函數時,每個角度θ的譜值由陣列導向矢量與噪聲子空間的正交關系可得。本文提出虛擬波束形成方法,陣列導向矢量為波束方向圖矢量Kronecker積,虛擬波束中主波瓣以外的混疊信號分量,會對譜估計性能產生影響。為進一步提升DOA估計性能,考慮修正信號主分量的波束方向圖,使虛擬波束主瓣內權值基本不變、抑制主波瓣外的權值起伏。

為保持各波束的連續性,選用高斯分布窗對波束方向圖進行修正。圖4所示為利用式(34)高斯窗修正后的歸一化波束方向圖,可見每個信號波束外的起伏被很好抑制,波束更加純凈。但是,高斯窗修正本質是為避免主瓣外波束起伏對空間譜估計產生影響,對DOA的粗估計結果不會有提升。

定量分析算法高斯窗修正效果其對信號譜估計精度的提升,取式(38)建立的波束誤差統計量,表征高斯窗修正前后,所抑制的權值起伏差異量在整個角度區間統計平均的歸一化誤差比。

(38)

式中:N表示θ遍歷角度區間(-90°,90°)區間的離散取值點數;下標序號i表示第i個虛擬波束。

以圖2仿真場景為例,進行100次蒙特卡羅仿真,如圖5所示。將3個高斯窗修正后波束虛擬波束的誤差做歸一化,可見高斯窗修正可平均對波瓣起伏均方根誤差(root mean square error, RMSE)抑制10%以上。

圖6所示為是否使用高斯窗修正,在100次蒙特卡羅仿真中,信號的DOA估計精度誤差統計平均。可見,使用高斯窗修正后的信號譜估計精度,較修正前有明顯提升,高斯窗修正前算法DOA估計精度為0.27°,修正后為0.06°,統計結果充分表明了高斯窗修正使信號等效導向矢量更加純凈,較大程度避免了信號混疊下對譜估計峰值計算的誤差影響。

2 仿真實驗

實驗參數設置見表1所示,空間場內共存在4個源信號,其中S1~S3設置同第1.3節不變,新增S4模擬主瓣干擾點源場景。陣列誤差由因子β整體控制,因子β越大,陣列誤差越大。相鄰陣元間互耦系數最大為0.1β,遠離依次減弱,5個陣元外互耦效應忽略不計。空間擾動誤差為乘性誤差,其大小可由幅相誤差、互耦誤差乘積近似代替,不再單獨設置。

表1 參數設置

實驗方法為基于陣列信號全FOC的DOA估計的MUSIC-like算法,以及對MUSIC-like中FOC去冗余的簡化算法。方便起見,分別記MUSIC-like及其簡化算法為FOC和ReFOC算法。其中,ReFOC算法運算更快,但犧牲了部分DOA精度。用本文構建高斯窗修正的虛擬波束進行DOA估計的方法,分別記為VBF-FOC和VBF-ReFOC。

2.1 信號DOA檢測性能分析

圖7所示為誤差因子β=1,即幅度、相位和陣元互耦有10%誤差時,4種方法的角度估計譜。其中,ReFOC估計精度較差,但經虛擬波束處理后的VBF-ReFOC精度明顯提升;且對比FOC和VBF-FOC(紅色與藍色線),0°的S1目標譜峰高度明顯改善,表明在波束構建時,利用信號空間特征值歸一化,彌補了非等功率信號的影響,有助于提升低SNR目標的DOA估計性能,與理論分析吻合。

DOA譜峰估計時,可能會因角度分辨率不夠、對非等功率信號和陣列誤差適應性差等因素,造成信號DOA漏檢。為考察幾種方法極限條件下信號檢出性能,設置場內含S1和S4兩點源信號,模擬主瓣干擾場景下目標DOA檢測性能。若區間內可形成兩個譜峰,則目標及干擾信號DOA均可成功估計,DOA檢測成功率為100%;若只能檢測到一個譜峰,則認為主瓣干擾信號S4淹沒目標信號S1,角度分辨率不夠致使目標檢測失敗,DOA檢測成功率為50%。圖8為β=1,兩點源角度間隔在5°~10°變化時,100次蒙特卡羅仿真的信號DOA檢測成功率。可見,利用信號全FOC進行估計的FOC方法,不容易產生漏檢,且經信號的虛擬波束化后,VBF-FOC、VBF-ReFOC兩種方法的目標信號DOA檢出成功率,較處理前性能有顯著提升。

2.2 DOA估計性能分析

除DOA的檢出性能外,更重要的評估指標是DOA估計精度。圖9為表1源信號參數下,陣列誤差因子β在0~1.5變化時,100次蒙特卡羅仿真下幾種方法DOA估計精度,評價指標選取DOA估計與真值間RMSE。其中,圖9(a)為場內4個源信號DOA估計的平均RMSE性能,考察算法DOA估計整體精度;圖9(b)為針對低SNR目標S1的DOA估計RMSE,考察算法對強干擾下弱目標DOA估計精度。結果可見,本文提出VBF-FOC算法的整體及對低SNR目標估計精度,都有明顯優勢;且在虛擬波束處理后,估計RMSE較之前方法都有很大下降。在陣列誤差不大于10%的影響下,VBF-FOC算法的源信號估計RMSE可穩定在0.5°以內,表現出對陣列誤差和非等功率信號情形下更強的適應性。

圖10為表1源信號參數下,弱目標信號SNR在-10~10 dB變化情況下,100次蒙特卡羅仿真的算法估計精度對比。顯而易見,文章提出基于VBF處理的兩種算法,低SNR目標DOA估計收斂更快。在都利用FOC進行計算的情況下,VBF-FOC在目標信號-6.5 dB時可收斂到1°估計精度,而對應FOC方法SNR需求是3 dB。相較而言,達到相同DOA精度,VBF-FOC對信號SNR的需求可降低平均10 dB,實現低SNR下的目標高精度DOA估計。

表2為兩種變化場景下,4種方法計算DOA估計平均精度值。場景1:固定SNR=0 dB,陣列誤差因子β在0~1.5變化;場景2:固定陣列誤差因子β=1,SNR在-10~10 dB變化。可見,VBF-FOC算法估計精度更高、適應性更強,且對比直接利用陣列接收數據的FOC進行DOA估計的算法,本文提出的VBF處理后算法的DOA估計精度平均可提升150%以上。

表2 不同算法DOA估計平均精度

2.3 運算速度分析

事實上,為保證虛擬波束方法DOA運算精度,第2.2節討論時,均采用波束數目K=N(陣元數)的等維波束構建方法。若K

通常對于一個K×K階矩陣進行特征值分解,其運算復雜性約為O(K3)。而FOC矩陣進行DOA估計時,相當于對原有陣列維數進行擴展,FOC矩陣階數越高,運算復雜度越大。而本文提出虛擬波束構建,雖然也需特征值分解,但卻是基于陣列信號的,對比FOC矩陣可視為低維矩陣。表3所示為100次蒙特卡羅仿真下,4種方法平均運算時間對比。當矩陣階數都為10時,VBF處理后算法運算時間只同比增加了約0.04 s,基本可忽略不計。

表3 算法運算時間

結合圖11分析,在計算精度和運算速度需要做平衡,取VBF=8時,VBF-FOC算法較FOC方法可在波束降維處理下,運算時間減小一半以上的情況下,仍可獲得DOA估計精度一倍的提升。

3 結 論

本文提出了一種在構建虛擬波束框架下,利用FOC進行DOA估計的改進算法(VBF-FOC和VBF-ReFOC)。算法通過陣列接收信號主成分提取和方向波束矢量高斯窗修正兩步,形成新的聯合導向矢量進行DOA空間譜估計。仿真結果表明,本文提出的VBF-FOC算法DOA估計精度最高,VBF-ReFOC精度次之但運算較快,二者在利用虛擬波束形成的處理后,可對源信號DOA估計精度較現有FOC方法提升150%以上,尤其對于場內同時存在多個非等功率源信號時,方法對低SNR目標精度提升更加明顯,表現出算法對復雜干擾環境DOA估計更好的適應性和魯棒性。在運算時間方面,虛擬波束形成對運算復雜度增加可忽略不計,且可通過虛擬波束降維處理,實現較高精度下的快速DOA估計。需要注意的是,本文在理論建模時采用均勻線陣,但是算法主體思想信號主成分提取和高階累積量特征分解,對于陣列形式約束較小,因此算法對于均勻圓陣等情形也適用。算法驗證時,考慮了陣列幅相、互耦等誤差因素影響,對工程應用中陣列誤差具備一定適應性。但是,考慮FOC的計算復雜性,在雷達等計算資源比較充分的平臺,推廣應用價值更好。

猜你喜歡
信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
7個信號,警惕寶寶要感冒
媽媽寶寶(2019年10期)2019-10-26 02:45:34
孩子停止長個的信號
《鐵道通信信號》訂閱單
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
基于Arduino的聯鎖信號控制接口研究
《鐵道通信信號》訂閱單
基于LabVIEW的力加載信號采集與PID控制
Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
主站蜘蛛池模板: 欧美第二区| 国产极品粉嫩小泬免费看| 久久伊人久久亚洲综合| 国产在线98福利播放视频免费| 园内精品自拍视频在线播放| 这里只有精品在线| 免费一级毛片不卡在线播放| 色成人亚洲| 国产一区二区三区在线观看视频| 免费一级大毛片a一观看不卡| 久久99国产乱子伦精品免| 欧美午夜小视频| 五月激情婷婷综合| 国产亚洲视频中文字幕视频| 天天综合色网| 免费a级毛片18以上观看精品| 青青热久麻豆精品视频在线观看| 毛片在线看网站| 国产视频资源在线观看| 欧美日本在线观看| 国产成人AV男人的天堂| 人妻丰满熟妇av五码区| 制服丝袜 91视频| 亚洲欧洲日产无码AV| 国产男人天堂| 成年人国产视频| 最新国语自产精品视频在| 亚洲AV无码久久天堂| 欧美日韩v| 国产成人做受免费视频| 亚洲最大福利视频网| 亚洲水蜜桃久久综合网站| 国产自无码视频在线观看| 青青操国产视频| 欧美精品在线观看视频| 国产欧美在线观看精品一区污| 韩国v欧美v亚洲v日本v| 欧美v在线| 欧美性精品| 亚洲品质国产精品无码| 亚洲AV无码久久精品色欲| 视频一本大道香蕉久在线播放| 亚洲成人在线网| 人人看人人鲁狠狠高清| 亚洲一区色| 天天色天天操综合网| 欧美激情视频一区| 亚洲欧美成人综合| 制服丝袜无码每日更新| 青青久视频| 久久久久久午夜精品| 国产亚洲高清在线精品99| 999国产精品永久免费视频精品久久| 在线精品欧美日韩| 国产精品无码影视久久久久久久 | 91免费国产高清观看| 国产对白刺激真实精品91| 亚洲无码A视频在线| 国产毛片不卡| 婷婷午夜影院| 2020国产精品视频| 亚洲精品国产日韩无码AV永久免费网| 华人在线亚洲欧美精品| 亚洲综合色婷婷中文字幕| 伊人福利视频| 国产成人精品一区二区不卡| 精品国产一区91在线| 亚洲天堂网视频| av手机版在线播放| 激情六月丁香婷婷| 91精品国产丝袜| 亚洲国产精品无码久久一线| 婷婷色在线视频| 国产成人a在线观看视频| 国产一二视频| 网友自拍视频精品区| 久久国产V一级毛多内射| 激情无码字幕综合| 国产乱子伦精品视频| 欧美性猛交xxxx乱大交极品| 91视频区| 欧美第九页|