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

穿墻雷達墻體參數估計以及補償方法綜述

2022-11-24 01:53:46梁步閣楊德貴朱政亮
無線電工程 2022年11期
關鍵詞:方法

肖 駿,梁步閣,楊德貴,朱政亮

(1.中南大學 自動化學院,湖南 長沙 410083;2.廈門大學 水聲通信與海洋信息技術教育部重點實驗室,福建 廈門 361005)

0 引言

穿墻雷達(Through-Wall-Radar,TWR)能夠利用電磁波穿透磚墻、混凝土、木板等介質對遮擋目標進行探測和定位,因此在反恐安全、災后救援、醫療監護等領域得到廣泛應用[1-3]。在美國聯邦通信委員會向民用領域開放超寬帶(Ultra Wide Band,UWB)技術后,UWB TWR成為TWR應用的主流,UWB TWR工作時發射UWB信號,即發射帶寬比超過中心頻率25%的信號[4-6]。

與窄帶TWR相比,利用大帶寬信號UWB TWR的優勢[7-8]在于:超高距離分辨率、抗多徑衰落、抗干擾能力強、穿透性強。

根據信號體制不同,TWR可分為沖激脈沖(Impulse Radar,IR)、步進頻連續波(Stepped Frequency Continuous Waveform,SFCW)、調頻連續波(Frequency Modulated Continuous Wave,FMCW);按技術復雜度,TWR天線體制設計逐步由一發一收(Single Input Single Output,SISO)到一發多收(Single Input Multiple Output,SIMO),現以多發多收(Multiple Input Multiple Output,MIMO)為主;根據天線陣列設置不同,可分為固定陣列和可移動陣列。

TWR主要是探測墻后人體目標和分析墻體結構,二者都需要研究墻體的介質特性以及雷達信號穿透墻體后的衰減與延時補償。隨著收發通道的增多,雷達分辨率更高、探測距離更遠,能大幅提升TWR的性能,且探測模式更多樣化,但墻體參數估計、墻體補償等問題也更復雜。

本文主要對墻體參數估計方法和墻體補償方法的發展研究進行回顧。首先,簡要介紹TWR的發展以及探測墻后雷達目標的工作原理;然后,針對目前墻體參數方法和墻體補償方法的研究現狀進行分析與比較;最后,對TWR技術的發展趨勢進行展望。

1 TWR發展

TWR憑借在無接觸式探測領域的性能優勢,具有良好的應用前景,國內外對此進行了大量的研究,并取得了一定的成果。目前,由于國外公司和科研機構起步較早,產品較為成熟,在商業化上具有一定的優勢。但是,隨著國內各個領域對TWR的需求急劇增大,我國眾多科研機構與企業展開合作,也取得了較多的研究成果,并積極推動商業化進程。

在國外,美國Time Domain公司于2003年研發的SV2000A1[9-10]用于軍隊在人口密集地區的作戰需求,如圖1所示。其可通過檢測相鄰回波的能量變化來檢測目標的運動。以色列Camero-Tech公司于2004年研發的XaverTM系列[11],可穿透主要非金屬障礙物對人體目標進行探測。XaverTM-800是該系列中最先進的一款,可以實現三維成像和目標檢測、跟蹤,設備與運行結果如圖2所示。

圖1 Time Domain的穿墻雷達——SV2000A1

(a)XaverTM-800穿墻雷達

在國內,某公司研發的SJ6000+雷達[12],可以穿透42 cm的墻體,探測18 m內靜止人體和27 m內移動人體目標的呼吸信號。IR-UWB多通道雷達中心頻率為400 MHz[13],具有多個自由度,可根據探測場景調整收發單元的位置和探測方向。Radar Eye TWR系統[14]使用無載頻雙極短脈沖,可以探測3~5 m的目標,且對周圍環境電磁干擾小,適合部署在敏感地區。CE/CEM系列產品[15]中的CE200 雷達可以對墻體后方25 m內的靜止目標、30 m內的運動目標進行多目標探測和二維定位;CEM400雷達能夠實現墻后目標三維成像。TWR R300[16]屬于多通道UWB脈沖雷達,工作頻率為100~2 500 MHz,可穿透主要非金屬以及低含水量物體進行單/多目標探測。

2 TWR工作原理

TWR主要任務是對非金屬介質障礙物后方的人體目標進行探測,雷達通過發射天線向探測區域輻射電磁信號,穿透障礙物之后經目標反射并二次穿透障礙物并被接收天線接收進行處理。傳播過程中,探測區域的環境反射波、系統噪聲以及其他噪聲隨著目標回波一起被接收。

圖3 TWR回波模型

(1)

式中,a(t)為直耦波;b(t)為墻體表面的反射波;xt(t)為墻后目標回波;n(t)為外部噪聲。在進行目標檢測前,需要對回波信號進行預處理,去除直耦波、墻體反射波以及外部噪聲,提高信噪比(Signal to Noise Ratio,SNR)。

圖4 脈絡圖

3 墻體參數估計方法

目前,墻體參數估計方法主要分為基于時延估計、基于墻體特征值匹配和基于成像質量評估等3類。

3.1 基于時延估計的墻體參數估計方法

基于時延估計的墻體參數估計方法目前研究較多,按照理論模型不同分為2類:雙基地法[18]和單基地法[19],如圖5所示。在模型中,假設墻體為由單層均勻介質組成,厚度為d,相對介電常數為εr,天線與墻體前側的距離為r。

(a)雙基地

雙基地法:通過分置收發天線來獲取墻體的前、后側回波,天線水平放置,距離為2L。在進行墻體參數估計時主要使用圖5(a)中表示的后墻回波,根據該模型可將后墻回波的延時表示為:

(2)

式中,x為折射點P的水平位置。根據折射定理和回波模型,εr可表示為:

(3)

也可近似[20]為:

(4)

(5)

式中,Lm為第m組數據對應的發射天線和接收天線水平距離的一半。

單基地法:通過收發共用天線獲取多次反射回波,且通常考慮電導率σ的影響。在進行參數估計時,通常提取前3次回波進行處理,計算前3次回波的幅值相繼比來估計墻體參數d,εr和σ。在理論模型中,前3次回波的幅值可表示為:

(6)

(7)

(8)

式中,α為墻體衰減系數;Γ為墻體反射系數。所以,3次回波的幅值相繼比ρ1和ρ2可分別表示為:

(9)

(10)

通過求解式(9)和式(10),可得到d,εr和α的關系式:

(11)

(12)

(13)

式中,t1,t2分別為前2次回波的延時。σ可由α求出,而對于低損耗墻體,σ可表示為:

(14)

除此之外,有學者提出使用圖5(b)中的前2次墻體反射回波進行參數估計[21]。通過分析前2次反射回波與d,εr,σ的關系,得到代價函數:

(15)

再根據測得的回波最小化式(15),得到εr和σ的估計值,d可由式(12)求得。

單-雙基地混合法:相較于雙基地法,單基地法過程簡單,計算量較小,且能對電導率σ進行估計,但需要對后墻的2次回波進行提取,而因墻體導致的回波衰減,會使回波的提取難度遞增。所以,基于時延估計的參數估計方法通常采用單-雙基地混合的方式[22-24]:先利用雙基地法對墻體的d和εr進行估計,再利用單基地模型中前2個回波的比值ρ1由式(13)和式(14)計算出墻體電導率。這樣,既降低了對雷達技術的要求,又能對墻體電導率進行估計,且減少計算量[25-27]。

基于時延估計的參數估計方法除同側架設天線外,還可以在墻體兩側分別架設發射、接收天線且與墻面保持一定距離,利用傳播延時和振幅衰減系數與d,εr的非線性關系進行參數估計[28]。但該方法需要在兩側架設天線的方式與實際環境差別過大,僅用于實驗測試。

基于時延估計的參數估計方法準確度與時延估計結果相關,時延估計方法不同,最終參數估計值準確度也不同。常用的時延估計方法有基于相關矩陣特征分解的估計算法[22-23]和基于稀疏重建的估計算法[24-26],還可以直接使用矢量網絡分析儀(Vector Network Analyzer,VNA)進行測量。其中,基于相關矩陣特征分解的估計方法有基于旋轉不變技術的信號參數估計(Estimating Signal Parameter via Rotational Invariance Techniques,ESPRIT)方法、多信號分類(Multiple Signal Classification,MUSIC)方法等;基于稀疏重建的估計方法有稀疏盲反卷積、正交匹配跟蹤(Orthogonal Matching Pursuit,OMP)等算法。前者需要對回波數據進行解相干處理,且在低SNR時性能降低,導致參數估計準確度降低;后者在低SNR時仍有較好性能。在d=15 cm,εr=6和σ=0.012 S/m仿真數據下進行不同SNR環境的驗證,結果如表1所示。d和εr相對均方根誤差(Relative Root Mean-Squared Error,RRMSE)如圖6所示[26]。

表1 ESPRIT和OMP方法在不同SNR下估計結果

(a)厚度

隨著智能算法的發展與廣泛應用,支持向量機(Support Vector Machine,SVM)也被用于墻體參數估計,可分為2種訓練方式:第1種方式[29]的理論模型如圖5(a)所示,固定收發天線且與墻體保持一定距離,利用不同d和εr下前側與后側回波的最大幅值和延時分別訓練d,εr,σ;第2種方式[30]的理論模型如圖7所示。

圖7 基于SVM的墻體參數估計模型

天線緊貼墻面基于SAR的原理形成合成孔徑,利用墻后目標在不同d和εr回波的最大幅值和延時分別訓練d和εr。

總而言之,基于時延估計的墻體參數估計方法,主要適用于UWB脈沖信號體制TWR,多通道體制或可移動天線。其中,單基地法主要用于收發共用或天線間距較小的TWR;雙基地法和單-雙基地法主要用于多通道體制或可移動天線的TWR,且雙基地法對設備要求相對較低。而收發天線分置于墻體兩側的方法受應用場景影響較大,實用性較差,基于SVM的方法由于智能算法對芯片的性能要求高,主要用于離線處理或處理、探測分離的TWR。

此類方法雖然簡單直接,理論上能夠分辨多層介質,但因為實際操作環境和雷達硬件的限制,當前研究都是使用仿真、高性能試驗設備實測或VNA直接測量的數據,現有雷達產品難以準確地提取墻體回波,無法達到預期效果,所以現有雷達均不采用此類方法。基于時延估計的墻體參數估計方法也展現了獨特的優點及發展前景,隨著雷達技術的發展,該類方法將成為TWR的常用方法之一。

3.2 基于墻體特征值匹配的墻體參數估計方法

基于墻體特征值匹配的墻體參數估計方法通過不同的信號處理方法在雷達回波中提取參數與已知墻體參數進行匹配,得到墻體參數的估計值。此類方法大都通過改變天線陣列的布局來獲取多組數據進行特征值匹配,以其目標位置進行墻體參數的估計。

文獻[31-32]通過設計多種陣列結構進行墻體參數的估計,天線陣列緊貼墻面,在不同天線陣列下對墻后目標進行定位、成像,并假定多組參數對目標位置進行補償,得到每種天線陣列下的目標位置軌跡,再比較任意2種天線陣列形成的軌跡,其交點就是目標真實位置,對應的墻體參數即為估計值。同樣,也有固定天線陣列結構,改變天線陣列與墻面的間隔,然后假定不同參數對墻后目標進行定位、成像,任意2種距離形成的軌跡交點即為目標真實位置,其對應的墻體參數為估計值[33]。

這2種方法都是通過改變陣列或陣列與墻面位置進行墻體特征值匹配的,操作相對復雜,受場地限制。為了減少操作,文獻[34]進行了改進,將天線陣列與位置固定,假定幾個不同的εr,再對每個εr下不同的d對墻后目標進行目標定位、成像,形成目標位置軌跡,軌跡交點即為目標真實位置。

假設d和εr均未知,預測值分別表示為de=dT+Δd和εe=εT+Δε,其中dT,εT為真實值,Δd,Δε為誤差。厚度誤差Δd對目標定位的影響如圖8所示,目標位置p對應墻體參數(dT,εe),因為Δd,目標位置偏移到q點,對應墻體參數(dT+Δd,εe)。根據文獻[32-33]的結論,Δxpq和Δypq可表示為:

圖8 厚度誤差Δd對目標定位的影響

(16)

Δypq=yq-yp=

(17)

式中,θto,p,φto,p和θro,p,φro,p為目標位置p到發射、接收天線路徑在墻體-空氣分界處的入射角和折射角。

因此,在假定一個εr時,以不同d對目標進行定位,將目標位置擬合為線性軌跡,斜率k可表示為:

(18)

由式(18)可見,軌跡斜率k僅與εr,θto,p,φto,p和θro,p,φro,p有關;當天線固定時,θto,p,φto,p和θro,p,φro,p只與εr相關,所以εr是影響斜率k的唯一因素。

通過對固定陣列在至少2個假定的εr進行目標定位或成像,在每個假設的εr下,目標位置通過不同的d擬合成線性軌跡,軌跡斜率僅與εr有關,并且Δε引起的位置偏移和Δd引起的位置偏移會抵消,軌跡都會經過目標真實位置,所以這些軌跡交點為目標真實位置預測值。而對于每個假定的εr,d與目標位置存在一種線性關系,可表示為:

(19)

該改進方法減少了參數估計的操作,但仍存在該類方法的一些缺點:一是需多次實驗,計算復雜,實用性較差;二是估計的參數值雖能對墻體干擾進行補償,得到高質量的目標定位和成像效果,但大都是匹配值而非真實值,且存在多組。這是因為墻體參數之間在對雷達回波的干擾上具有一定相關性,它們的誤差在一些情況下會抵消,糾正了相關組合成像場景中每個像素波形返回的聚焦延遲集。所以,該類方法在應用時通常先通過肉眼等方式根據墻體類型縮小參數范圍來得到相對真實的參數估計值,極容易產生偏差,且僅適用于單層介質。因此,該類方法亟需對墻體參數間的耦合機理進行研究,以解決多組解和非真實解的問題。

3.3 基于成像質量評估的墻體參數估計方法

相較于前2類方法的高要求,基于成像質量評估的墻體參數估計方法不受外界環境影響,適用性較好,在現有雷達技術下能夠很好地應用。此類方法是通過對不同墻體參數組合計算出補償因子進行成像,再對圖像進行質量評估,找出成像效果最優的參數,其流程如圖9所示。此類方法的準確度與成像算法、圖像質量評估函數以及預估參數范圍相關。在此類方法中,成像算法主要使用后向投影(Back Projection,BP)算法、三角定位成像算法等具有代表性的方法。

圖9 基于成像質量評估的墻體參數估計方法

質量評估函數的選取上,文獻[35]通過分析理想墻體模型中目標和墻體的電磁波散射場,與真實墻體模型中目標和墻體的電磁波散射場進行對比,構建出真實墻體模型與理想墻體模型的目標函數,再針對墻體參數進行迭代優化,使目標函數最小化,從而得到墻體參數的估計值。

此類方法是對成像后圖像的質量進行評估,故常用評估方法大都是基于圖像自聚焦技術的,主要分為2類:熵和對比度。其中,熵類方法以圖像信息熵為評估指標[36-38],熵值越小成像效果越好,圖像信息熵表示為:

(20)

式中,I(rk,kk)為圖像在(k,kk)處的相對強度,可表示為:

(21)

對比度方法常以銳利度為評估指標[39],銳利度越高成像效果越好,銳利度為:

(22)

式中,n為第n銳利度。

同樣,也有利用圖像峰度作為評估指標[40],將Q個像素點圖像的峰度定義為:

(23)

(24)

(25)

為了減少計算量,即縮小參數估計范圍,文獻[37]通過分析SFCW雷達接收的墻體回波特性和電磁波在墻體中的多徑效應,得到d和εr與墻體回波的距離向壓縮關系:

(26)

式中,Δd為壓縮后2個峰值的延時差;Δf為信號的頻率階;θ1pq為折射角;L為IFFT點數。墻體的材質大多可以直接確定,得到εr的估計范圍,再根據式(26)計算d的估計范圍。然后,利用范圍內參數計算補償因子對墻后目標進行成像,再進行圖像質量評估,得到較精準的估計值。

除此之外,也可通過計算出d取值范圍內對應的折射率,使成像效果最佳,再擬合出參數對的最大銳利度曲線[39]。在這條曲線上取幾組參數成像,其成像效果相近,驗證了在成像方面d與εr之間存在一定相關性,即當滿足參數耦合條件時,不同參數組合能實現相似的成像效果。

總而言之,基于成像質量匹配的墻體參數估計方法實際上是通過多組墻體參數預測值的成像效果去找出最優的幾組預測值作為估計值,在此類方法中雖可用一些方法縮小預測值范圍,但由于成像的計算量較大,此類方法仍耗時較長,難以滿足當前需求。此類方法最終的參數估計值也不是唯一值,如第二類方法只是參數的匹配值而非真實值,需對在成像質量上墻體參數間的耦合機理進行研究,以確定墻體參數唯一、真實估計值。

4 墻體補償方法

墻體補償方法大都是通過尋找雷達信號在墻體與空氣交界處的折射點,再計算出每個位置目標回波時延,對雷達回波進行補償,實現位置修正和高聚焦成像。目前,國防科技大學、斯洛伐克技術大學等的專家學者對墻體補償方法進行研究,也取得了一定的成果。其中,傳統墻體補償方法都是基于牛頓-霍納迭代法[40]和牛頓-拉夫遜法[42]對四階多項式進行求解,需耗費巨大的計算量。為了降低計算量,專家們研究出多種方法進行替代。

最短路徑的墻體補償方法[44-45]是基于電磁波總是沿著最短延時路徑進行傳播的理論,即使在穿透介質發生折射時仍會尋找一條特定的路徑使傳輸時延最小,理論模型如圖10所示。

圖10 最短路徑法模型

通過對成像平面內所有像素點的最短延時進行計算:

(27)

再根據回波延時進行遍歷,確定目標位置,以實現墻體補償。該方法補償效果較好,減少了一部分計算量,但探測區域偏大時計算量仍較大,限制了在實際探測中的應用。

折射角近似法[46]的理論模型如圖11所示。

圖11 折射角近似法模型

將天線緊貼墻面置于A2(xA2,yA2),目標P(xP,yP)位于另一側。設天線在墻體另一側投影點為Q1(x1,y1),電磁波路徑折射點設為Q2(x2,y2)。

根據圖11,電磁波從天線A2傳輸到目標P的時延τA2P可表示為:

τA2P=τA2Q2+τQ2P=

(28)

式中,唯一不確定的參數是折射點的水平位置x2,可表示為:

x2=d/cotφ2+xA2,

(29)

式中,cotφ2可以通過折射定理的變式求得:

cot2φ2+1=(cot2θ2+1)·εr。

(30)

為了減少計算量,可假設存在一個虛擬天線A3(xA3,yA3),其折射點為A2與P連線和墻體后側的交點Q3(x3,y3),該點對應入射角如下:

(31)

圖12 不同εr下折射角及余弦值隨入射角變化

基于折射角近似的墻體補償方法能夠快速對墻后目標位置進行修正,但近距離目標難以聚焦成像,且效果比最短路徑法差,在實際中也存在一定限制。為了進一步提高目標位置校正效果,在最短路徑法和折射角近似法的基礎上,有學者提出一種基于折射角近似的最短路徑墻體補償方法[47]。該方法基于折射角近似的原理,確定折射點區間,在精度不變的情況下大幅減少了最短路徑法的計算量。

(32)

再求出該點的折射角φ1,找出對應的M1(xM1,yM1),使得A2M1‖A1Q1。這樣,M1可作為最短路徑法折射點搜索區間終點。在確定最短路徑法的折射點區間后,取N個折射點位置形成集合,再對成像區域所有像素點的最短時延進行計算。

當前,墻體補償方法雖對目標位置修正、成像聚焦具有一定的效果,但仍有著局限性,無法適用于所有場景。傳統方法能夠精確計算墻體的延時,但計算量過大,無法滿足應用需求;固定延時法能夠迅速對墻體延時進行粗略補償,但無法完全校正目標位置;最短路徑法能夠得到與傳統方法相似的效果,且計算量遠小于前者,但在探測范圍較大時仍需要較大的計算量;折射角近似法能夠快速校正目標位置,但近距離目標難以精準定位或聚焦成像;基于折射角近似的最短路徑延時補償方法在最短路徑法的基礎上減少一定計算量,但仍難以達到應用需求。不同于基于尋找折射點位置的墻體補償方法,也有學者提出一種基于圖像域的補償方法[41],首先,在忽略墻體影響下進行成像;然后,根據墻體參數計算出補償因子直接對圖像進行補償,只需要一次成像操作就能對多個目標進行補償,很大程度減少了計算量。

因此,目前對墻體補償方法的研究都旨在保證補償精度的情況下減少計算量,取得了一定的成果,展現了在TWR定位、成像領域的應用潛力,但計算量和精度仍難以同時達到需求,無法滿足所有場景,需進一步提升。對于天線體制,多發多收相較于一發多收的區別主要在于,不同發射天線的發射信號不同,需要在接收信號中區分出不同的發射信號,在補償時每個接收通道需要進行發射信號區分,并考慮到發射天線位置的不同。

5 結束語

本文圍繞TWR目標位置修正和成像聚焦,對墻體參數估計和補償方法近年的發展進行論述,通過對TWR現有墻體參數估計和補償方法的分析,介紹了雷達信號在穿透介質時的影響,講述了多種墻體參數估計和補償方法的優勢和適用領域,以及在某些方面的不足。

通過對國內外已發表的文獻進行歸納總結,可以發現大量學者和專家對墻體參數估計方法進行了多方面的研究,并在一些方面都取得了不錯的成果,能夠為目標位置修正、成像聚焦和介質檢測等方面提供較準確的數據支撐,具有一定的應用價值;而墻體補償方法能夠在保證一定精度的情況下大幅度減少計算量,展現了在TWR對墻后目標精準定位和成像領域的應用潛力。

對于未來的反恐安全、災后救援和醫療監護等領域,TWR對墻體特征參數估計和墻后目標定位、成像的速度、精準度有著很高的要求,結合目前的發展現狀以及TWR的應用需求,墻體參數估計和墻體補償領域的發展趨勢可能有:

① 復雜墻體的參數估計。目前,墻體參數估計方法大都基于單層均勻介質進行研究,在多層等復雜墻體的效果較差,而TWR應用場景中需要穿透的介質往往組成復雜,既有單層的也有多層的,既有均勻的也有非均勻的,需要加強對多層、非均勻介質參數估計的研究。

② 墻體參數估計方法的實用性和適用性。目前,由于TWR設備硬件和應用場景的限制,基于時延估計的墻體參數估計方法的實用性相對較差,需要對雷達硬件設備進行改進、升級,在保證設備便攜性的前提上提高雷達各模塊的性能,并嘗試使用其他信號體制或多模融合,與其他方式相結合,實現互補;基于成像質量評估的方法雖被認為是如今最有可能用于實際的方法,但此類方法得到的估計值存在多組且不一定為真實值,在多層、非均勻介質上更甚,無法有效解決多層、非均勻介質的參數估計問題,其根本原因是缺少對墻體參數之間耦合機理的研究;而基于特征匹配的墻體參數估計方法兩方面都存在一定問題,上述兩方面都需要進行改進。

③ 墻體補償方法計算速度、精度的提升與平衡。在墻體參數已知的情況下電磁波的傳播機理研究較為成熟,且相較于傳統的墻體補償方法,當前對墻體補償方法的研究旨在減少計算量,保持一定的結果精度。雖然目前的研究已大幅度減少了墻體補償的計算量,但仍存在一些缺點,難以適用于所有場景,且計算量和精度仍難以同時達到應用標準。所以,在當前研究的基礎上,可以結合多種方法進行改進或是與其他領域結合。

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲一级毛片免费观看| 无码福利日韩神码福利片| 国产肉感大码AV无码| 22sihu国产精品视频影视资讯| 久久中文电影| 999国产精品永久免费视频精品久久| 欧美色综合网站| 国产特级毛片| 噜噜噜久久| 91在线精品麻豆欧美在线| 国产精品久久久精品三级| 国产一二三区在线| 国产资源免费观看| 中文字幕欧美日韩| 国产成人综合亚洲欧美在| 五月婷婷综合网| 亚洲精品第一在线观看视频| 高潮爽到爆的喷水女主播视频| 91欧美在线| 爱色欧美亚洲综合图区| 在线人成精品免费视频| 亚洲日韩国产精品综合在线观看| 99久久国产综合精品2020| 日本91在线| a毛片基地免费大全| 免费av一区二区三区在线| 欧美午夜理伦三级在线观看| 亚洲国产综合精品中文第一| 91网址在线播放| 欧美日韩国产综合视频在线观看| 国产福利大秀91| 日韩欧美中文字幕一本| 亚洲成人高清在线观看| 精品少妇人妻一区二区| 亚洲日本精品一区二区| 国产综合亚洲欧洲区精品无码| 国产网站在线看| 91在线精品麻豆欧美在线| 成人毛片在线播放| 国产高清自拍视频| 91精品免费高清在线| 国产一级精品毛片基地| 在线亚洲精品自拍| 免费A级毛片无码无遮挡| 爱色欧美亚洲综合图区| 国产一线在线| 2021国产精品自产拍在线| 99九九成人免费视频精品| 亚洲精品久综合蜜| 欲色天天综合网| 色亚洲成人| 黄色网页在线观看| 日韩毛片免费观看| 午夜精品福利影院| 91精品专区国产盗摄| 人妻精品全国免费视频| 亚洲国产中文欧美在线人成大黄瓜| 久久国产精品77777| 99热免费在线| 久久九九热视频| 亚洲色无码专线精品观看| 国产精品永久久久久| 久久婷婷色综合老司机| 国产成人高清精品免费软件| 午夜老司机永久免费看片| 国产激爽大片高清在线观看| 精品国产美女福到在线直播| 国产h视频在线观看视频| 青青草久久伊人| 日本一区中文字幕最新在线| 国产美女在线免费观看| 欧美中出一区二区| 国产一在线| 久久77777| 91亚瑟视频| 乱色熟女综合一区二区| 国产成人精品男人的天堂| 久久夜夜视频| 亚洲国产日韩欧美在线| 亚洲永久色| 成人一级黄色毛片| 亚洲国产看片基地久久1024|