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

        ?

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

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

        劉甲磊, 馬佳智, 施龍飛

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

        0 引 言

        信號波達方向估計(direction of arrival,DOA)是陣列信號處理的重要研究內容,廣泛應用于雷達、通信、聲吶、無源定位等各個領域[1-2]。最具代表性的子空間多重信號分類(multiple signal classification,MUSIC)及衍生的多種算法[3-4],構成了陣列信號DOA領域的廣泛應用基礎,但其本質為基于信號二階統(tǒng)計特性構建的譜估計算法。相較信號二階統(tǒng)計特性,信號四階統(tǒng)計特性因其對任何形式的高斯過程的不敏感性,可以最大程度地抑制高斯噪聲的突出優(yōu)點,使其在盲源分離和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估計方法,表現(xiàn)了FOC算法在高斯色噪聲背景下的有效性。Xue等[15]將FOC用于均勻圓陣DOA估計,獲得較好的估計精度。文獻[16-17]提出了一種對MUSIC-like算法去冗余的FOC的DOA方法,只利用部分無冗余的FOC切片構建FOC矩陣,矩陣階數大大下降,利于工程化實現(xiàn)。但是其構建的去冗余FOC矩陣,破壞了原有矩陣的對稱性,DOA估計精度有所下降,Shi等[18]提出了一種構建Toeplitz去冗余FOC陣的方法。但由于FOC對野值非常敏感,真實情況陣列誤差未知且未被校準時,估計性能下降十分嚴重。文獻[19-21]分別建立陣列幅相和互耦誤差下基于FOC的DOA估計模型,但每種方法只能部分消除一種誤差來源,較實際多誤差因素耦合時并未提及。

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

        1 算法原理

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

        1.1 信號模型

        以N元均勻等距線陣模型,陣元間距為d,雷達發(fā)射信號波長為λ,陣列波束指向為θ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,由于隨著陣元間距的增大,互耦效應快速減弱,通??蓪建模為部分為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子空間分解

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

        (9)

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

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

        (10)

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

        (11)

        (12)

        (13)

        實際應用中,由于不同源信號來自于不同系統(tǒng),可認為各信源間相互獨立,則源信號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作為信號的四階統(tǒng)計特性,其高階的計算過程決定了其對于數據野值很敏感。當場景中存在多個非等功率源信源,信號強度差異大時,SNR大的信源,其四階統(tǒng)計特性會遠強于小SNR信源,這就造成了復雜場景下低SNR目標DOA估計的困難。

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

        RXX=XXH=UΛUH

        (21)

        U=[u1,u2,…,uN]

        (22)

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

        (23)

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

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

        (24)

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

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

        (25)

        (26)

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

        以虛擬波束信號為輸入,進行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(θ)中,當θ遍歷角度區(qū)間(-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)

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

        (36)

        (37)

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

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

        實驗場景:均勻線陣陣元數N=10,陣元間距半波長,3 dB主瓣寬度10.2°??臻g3個獨立點源信號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的轉發(fā)多假目標干擾信號,S3為高斯白壓制干擾信號。

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

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

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

        (38)

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

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

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

        2 仿真實驗

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

        表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檢測性能。若區(qū)間內可形成兩個譜峰,則目標及干擾信號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目標估計精度,都有明顯優(yōu)勢;且在虛擬波束處理后,估計RMSE較之前方法都有很大下降。在陣列誤差不大于10%的影響下,VBF-FOC算法的源信號估計RMSE可穩(wěn)定在0.5°以內,表現(xiàn)出對陣列誤差和非等功率信號情形下更強的適應性。

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

        猜你喜歡
        信號
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        7個信號,警惕寶寶要感冒
        媽媽寶寶(2019年10期)2019-10-26 02:45:34
        孩子停止長個的信號
        《鐵道通信信號》訂閱單
        基于FPGA的多功能信號發(fā)生器的設計
        電子制作(2018年11期)2018-08-04 03:25:42
        基于Arduino的聯(lián)鎖信號控制接口研究
        《鐵道通信信號》訂閱單
        基于LabVIEW的力加載信號采集與PID控制
        Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
        久久天天躁夜夜躁狠狠躁2022| av色一区二区三区精品| 免费人成小说在线观看网站| 国产在线一区二区三区av| 国产午夜福利精品| 日本精品熟妇一区二区三区| 精品久久久少妇一区二区| 人妻少妇偷人精品无码| 国产精品麻豆综合在线| 亚洲一码二码在线观看| 国产在线视频91九色| 一本色道久久88综合日韩精品| 国产女人成人精品视频| 国产三级在线观看性色av| 最新日本人妻中文字幕| 免费无码又黄又爽又刺激| 国产成人影院一区二区| 日韩国产有码精品一区二在线| 日本高清一道本一区二区| 东北女人毛多水多牲交视频| 日韩h网站| 国产一级自拍av播放| 极品少妇xxxx精品少妇偷拍| 欧洲成人午夜精品无码区久久| 色播中文字幕在线视频| 国产美女高潮流白浆视频| 亚洲va欧美va日韩va成人网| 国产精品 视频一区 二区三区| 日本少妇爽的大叫高潮了| 国产精品国产三级国产专播下| 久久久久人妻一区精品色欧美| 麻豆久久五月国产综合| 中文字幕视频一区二区| 高h小月被几个老头调教| 免费黄色电影在线观看| 人妻av一区二区三区高| 国产av剧情刺激对白| 欧美黑人xxxx又粗又长| 日韩欧美国产丝袜视频| 81久久免费精品国产色夜| 亚洲av无码乱码国产精品久久|