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

基于雙向拋物方程逆算法的障礙物定位技術研究

2015-02-17 02:55:56王昆龍云亮劉震宇
電波科學學報 2015年6期

王昆 龍云亮 劉震宇

(1.廣東工業大學信息工程學院,廣州 510006;2.中山大學電子與通信工程系,廣州 510006)

?

基于雙向拋物方程逆算法的障礙物定位技術研究

王昆1龍云亮2劉震宇1

(1.廣東工業大學信息工程學院,廣州 510006;2.中山大學電子與通信工程系,廣州 510006)

摘要在單向拋物方程法的逆算法即逆繞射拋物方程法的基礎上,研究了雙向拋物方程法的逆算法,對障礙物引起的后向場進行了逆繞射運算,并根據該逆繞射場分布來確定障礙物的位置和高度.基于此逆算法對刃峰形障礙物的定位問題進行了數值仿真求解,分析了天線高度和刃峰位置的變化對定位精度的影響.數值結果表明雙向拋物方程逆算法能用于對障礙物進行定位,且對單刃峰的定位精度很高.

關鍵詞雙向拋物方程;逆算法;定位;刃峰

資助項目: 廣東省自然科學基金(S2013040013643); 國家自然科學基金(61401106,41376041); 教育部博士點基金(20130171110024); 廣州市科技計劃項目科學研究專項(2014J4100206)

聯系人: 王昆 E-mail:wkkelly@gdut.edu.cn

引言

最早由Leontovich和Fock提出的拋物方程法(Parabolic Equation, PE)是波動方程的一種前向近軸傳播近似[1], 許多研究者對其進行了深入研究,并用其求解了各種復雜環境的電波傳播問題[2-10].傳統的PE法是一種單向算法,只考慮了前向傳播波,忽略了后向傳播波,而當有障礙物存在時,障礙物的反射可能會引起較強的后向傳播波.Oraizi最先研究了電波傳播問題中的雙向拋物方程法(Two-Way Parabolic Equation,2WPE)的理論和應用,用其計算了短距離路徑的單刃峰和多刃峰環境的前向和后向傳播的電磁場,并用一致性繞射理論(Uniform Theory of Diffraction, UTD)和時域有限差分(Finite-Differenie Time-Domain, FDTD)法驗證了2WPE的有效性和準確性[11]. Ozgun研究了雙向分步步進拋物方程法(Two-Way Split-Step Parabolic Equation,2W-SSPE),計算了對流層遠距離多刃峰環境的雙向電波傳播問題,考慮了大氣波導的影響[12]. Apaydin提出了基于有限元法求解的雙向拋物方程模型(Two-Way Finite Element Based Parabolic Equation, 2W-FEMPE),在計算步長、仿真時間、源設置等方面對2W-SSPE和2W-FEMPE兩種模型做了比較[13]. 我們曾提出了改進的雙向拋物方程法,推導了修正的后向拋物方程,采用雙向分步步進的離散混合傅里葉變換處理有限導電地面,對混合的任意不規則地形環境的雙向電波傳播問題進行了分析[14-16].

單向拋物方程法(One-Way Parabolic Equation, 1WPE)是從天線源的初始場出發,沿前向求解電磁波的空間分布,屬于正向問題;而逆算法或者稱為逆繞射拋物方程法則與之相反,是由某一接收位置的場分布出發,對電磁波的空間分布進行逆繞射反演求解,從反演場分布中判斷和估計源的位置信息,屬于逆問題. Spencer和Walker最先研究了單向拋物方程法的逆算法,提出了逆繞射拋物方程定位系統(Inverse Diffraction Parabolic Equation Localisation System,IDPELS),給出了逆算法仿真結果,并在澳大利亞進行了實測[17-18]. IDPELS采用大天線陣或可移動裝置測量某一位置處的接收信號場剖面,將此接收場剖面進行逆繞射傳播,由反演的逆繞射電磁場分布中的收斂點給出發射源的位置. IDPELS系統對于復雜環境下的輻射源的無源定位有著很好的適用性和穩定性,可以廣泛應用在很多領域,如神經網絡定位、自動移動機器人導航、全球定位系統、航空電子戰、雷達定向武器系統和安全通信系統等. Hawkes等人研究了三維拋物方程法的逆繞射問題,模型中考慮了大氣折射率與不規則地形的影響[19]. 國內的郭建炎也研究了拋物方程的逆算法問題,并用逆算法模型對單個或者多個發射源進行定位,分析了森林對發射源定位的影響,改進和完善了誤差橢圓的表示法[20]. 李德鑫等人在理想大氣輻射源定位的基礎上,分析了標準大氣、地形遮蔽條件下的一個和兩個輻射源的定位問題,對定位誤差、數據點取值、誤差橢圓的繪制等問題進行研究,給出了參數取值的方法及建議,提高了算法的定位精度[21-22].

上述均為針對單向拋物方程法的逆算法的研究,并將其應用于發射源定位問題;而本文研究雙向拋物方程法的逆算法,或稱為逆繞射雙向拋物方程法,并將其用于障礙物探測和定位問題.以刃峰形障礙物為例,當發射天線輻射的電磁波在前向傳播過程中遇到刃峰時會產生明顯的后向反射,形成后向傳播波,用2WPE法可以計算前向和后向傳播波疊加的總電磁場. 如果在天線源和刃峰間選擇某一位置作為接收點,用2WPE法計算得到的總電磁場減去1WPE法計算得到的前向場,就得到刃峰反射的后向場剖面,對此后向場進行逆繞射傳播計算,得到反演的逆繞射場分布,根據反演的逆繞射場分布情況可以分析和確定刃峰的位置和高度. 這種方法也可以推廣到任意不規則形狀障礙物的檢測與定位.

1雙向拋物方程法的逆算法模型

假設電磁場的時諧因子為e-iωt,且所有場分量與方位角無關,x為水平距離,z為垂直高度,則麥克斯韋方程組可簡化為二維標量亥姆霍茲(Helmholtz)方程,在直角坐標系中,方程為[8]

(1)

式中: U在水平極化和垂直極化中分別代表電場或者磁場的橫向分量; k是波數; n=n(x,z)是隨距離和高度而變化的折射率.

根據微分算子理論對式(1)進行因式分解,得到前向和后向兩個拋物方程,分別引入前向和后向簡化函數uf(x,z)=e-ikxUf(x,z)和ub(x,z)=eikxUb(x,z),然后采用泰勒近似法對方程中的偽微分算子進行近似,就得到標準雙向拋物方程,或稱為窄角雙向拋物方程(Two-WayNarrowAnglePE, 2W-NAPE)[14-15]:

(2)

(3)

此處引入了地球曲率的影響,ae為地球半徑,f表示前向,b表示后向. 前向和后向傳播波疊加的總場為U=uf·eikx+ube-ikx.

采用雙向分步步進混合傅里葉變換法求解窄角雙向拋物方程,在x+Δx和x-Δx處的場解分別為:

(4)

(5)

式(4)和式(5)表明,x+Δx處的前向場和x-Δx處的后向場可以分別由x處的前向和后向場步進遞推得到,因此,也可以通過反演的方法由x+Δx和x-Δx處的場逆推得到x處的前向和后向場分布,即雙向拋物方程的逆算法可表示為:

(6)

(7)

由上可見,如果測量得到某點的前向場分布剖面,則可以根據式(6)的逆繞射前向拋物方程法反演得到空間中的前向場分布,且逆繞射前向場分布中的收斂點就代表發射源或其鏡像的位置,這就是單向拋物方程法的逆算法. 如果測量并計算得到某點的障礙物引起的后向傳播場剖面,則通過式(7),由后向場在x-Δx處的值可以反演得到后向場在x處的值,這樣通過多次步進計算就能得到障礙物引起的整個后向場空間分布,在反演的后向場分布中包含了產生后向場的障礙物的位置和高度信息,據此可以進行障礙物的探測和定位,此即雙向拋物方程法的逆算法.

綜上所述,總結雙向拋物方程法的逆算法的步驟如下:① 選擇和設置天線源,向傳播路徑上的障礙物輻射電磁波,并分別用單向拋物方程法和雙向拋物方程法計算空間電磁場分布;② 在天線源和障礙物之間選擇某一接收位置(在實測中為便于測量,選擇的接收位置盡量靠近發射天線),在接收位置處,用2WPE法計算得到的總場減去1WPE法計算得到的前向場,得到的差值就是障礙物引起的在接收位置處的后向場剖面. 需要注意和說明的是,在實驗和實際運用中,需要用大天線陣或者可移動設備來測量接收信號剖面,由于很難在較大高度范圍上連續地測量,所以只能得到部分測量值,這將造成逆算法的定位精確度降低;③ 從后向場剖面出發,通過逆繞射后向拋物方程(7)步進逆推后向電磁場的空間分布;④ 分析用逆算法反演得到的后向電磁場空間分布,根據其中的場轉折點的信息來識別和判斷障礙物的位置和高度,從而實現對障礙物的探測和定位.

2數值算例

在以下所有算例中,如無特殊說明,均假設發射天線為水平極化的高斯天線,其3 dB垂直波束寬度為3°,天線仰角為0°. 假設標準大氣環境,有限導電地表為中等干燥地面,工作頻率為900 MHz,電磁波傳播的最大距離為10 km,距離步長為10 m.

考慮障礙物是單刃峰的情況,假設刃峰是理想導體,位于x=xe處,高度為he,如圖1所示.

圖1 單刃峰地形

在數值計算中,采用1WPE算法時,從x=0處的天線初始場開始,用式(4)將場在前向步進,將單刃峰當做理想導體處理,直至最大距離處;采用2 WPE算法時,當步進的前向場在xe處遇到刃峰時被分解為向前向和后向傳播的兩個分量(分別是+x和-x方向). 在刃峰處應用適當的邊界條件可以得到后向反射初始場,用式(5)將后向反射初始場在后向步進傳播,直至到達天線處,而x>xe區域的前向波仍以原方式在+x方向傳播,將前向場和后向場疊加起來得到總場. 最后用逆算法確定刃峰的位置和高度.

2.1 發射天線高度變化對定位精度的影響

首先假設在距離發射天線5 km的位置上存在一個高度為100 m的單刃峰,即xe=5 km,he=100 m. 分別考查發射天線高度為50 m,100 m,150 m時的單刃峰定位的效果和精確度.

當發射天線高度為50 m時,在圖2(a)中畫出了用2WPE算法計算得到的雙向傳播因子覆蓋圖. 在距離發射天線100 m,即與單刃峰相距4.9 km的位置上,將2WPE算法在此處計算得到的總場值和1WPE算法在此處計算得到的前向場值相減,得到單刃峰在此處產生的后向場剖面,如圖2(b)所示.

(a) 單刃峰的雙向電波傳播因子分布圖

(b) 在x=100 m處的后向場剖面圖2 單刃峰環境的雙向電波傳播(天線高50 m,f=900 MHz)

從圖2(b)的后向場剖面出發進行逆繞射遞推計算,在圖3(a)中給出了單刃峰后向場的逆繞射傳播因子分布,圖中實際的單刃峰頂點用圓與十字線標出. 可見,逆繞射傳播因子在單刃峰頂點附近有明顯的轉折點,且轉折點與實際單刃峰頂點位置十分接近,由此可以估計單刃峰的位置和高度.

為了更清楚地顯示場的變化,圖3(b)中畫出了單刃峰頂點附近的局部逆繞射傳播因子分布,并用“×”標示出轉折點,由此可判定單刃峰位于距離接收點4 915.7 m處,即距離天線5 015.7 m處,相對誤差為0.32%;刃峰高度約為100.446 m,相對誤差為0.446%.

(a) 逆繞射傳播因子分布圖

(b) 刃峰頂點附近的局部逆繞射傳播因子分布圖圖3 單刃峰環境2WPE逆算法的數值結果(天線高50 m,f=900 MHz)

當天線高度增加到100 m時,其他條件不變,在圖4(a)畫出了x=100 m處的后向場剖面,圖4(b)給出了單刃峰頂點附近的局部逆繞射傳播因子,并用“×”標示出轉折點,可見在單刃峰頂點附近仍有明顯轉折點,確定單刃峰位于距離接收點4 919.3 m處,即距離天線5 019.3 m處,相對誤差為0.395%;刃峰高度約為100.556 m,相對誤差為0.556%.

(a) 在x=100 m處的后向場剖面

(b) 刃峰頂點附近的局部逆繞射傳播因子分布圖圖4 單刃峰環境2WPE及逆算法的數值結果(天線高100 m,f=900 MHz)

當發射天線高度增加到150 m時,其他條件不變,在圖5(a)中畫出了x=100 m處的后向場剖面,圖5(b)給出了單刃峰頂點附近的局部逆繞射傳播因子分布,并用“×”標示出轉折點,可確定單刃峰位于距離接收點4 923.1 m處,即距離天線5 023.1 m處,相對誤差為0.47%;刃峰高度約為100.998 m,相對誤差為0.998%.

由圖3(b)、4(b)、5(b)可見,當發射天線高度從50 m增加到100 m和150 m時,由于接收到的單刃峰反射的后向場信息減少,所以單刃峰的定位誤差隨天線高度增加而略有增大,如表1所示. 因此,在仿真和實測中,天線高度不宜取高,應低于障礙物高度,此時在天線前的接收點處能獲得豐富的后向場信息,從而可以對障礙物進行精確定位,特別是在實驗測量中,較低的天線高度更有利于實測系統的建立和操作.

(a) 在x=100 m處的后向場剖面

(b) 刃峰頂點附近的局部逆繞射傳播因子分布圖圖5 單刃峰環境2WPE及逆算法的數值結果(天線高150 m,f=900 MHz)

表1 不同高度天線對應的單刃峰定位結果

2.2 刃峰位置變化對定位精度的影響

下面考察刃峰位置的變化對逆算法定位精度的影響.假設天線高度為50 m,單刃峰的位置分別變化到距離發射天線xe=3 km和xe=8 km處,高度仍為100 m,其它條件不變.

當刃峰位于xe=3 km時,圖6(a)和(b)中分別給出了全空間的逆繞射傳播因子分布和單刃峰頂點附近的局部逆繞射傳播因子分布. 可見,逆繞射傳播因子在刃峰頂點附近有明顯的轉折點. 從圖6(b)中的轉折點可以確定刃峰位于距離接收點2 901.4 m處,即距離發射天線3 001.4 m處,相對誤差為0.048%;刃峰高度約為100.4 m,相對誤差為0.4%.

(a) 逆繞射傳播因子分布圖

(b) 刃峰頂點附近的局部逆繞射傳播因子分布圖圖6 2WPE逆算法的數值結果(天線高50 m,xe=3 km,f=900 MHz)

當刃峰位于xe=8 km時,在圖7 (a)和(b)中分別給出了全空間的逆繞射傳播因子分布和單刃峰頂點附近的局部逆繞射傳播因子分布. 從圖7(b)中的轉折點可以確定單刃峰位于距離接收點7 967.8 m處,即距離發射天線8 067.8 m處,相對誤差為0.859%;刃峰高度約為101.879 m,相對誤差為1.879%.

由圖6、3、7可見,單刃峰的位置從xe=3 km變化到5 km,8 km時,刃峰高度和位置的定位誤差都逐漸增大,如表2所示. 可見障礙物離天線越近,定位越準確,障礙物離天線越遠,定位精度越低. 其主要原因是隨著距離的增加,后向場擴散了,有更多的能量輻射損耗掉了,接收點接收到的信息減少,因此定位誤差增加.

(a) 逆繞射傳播因子分布圖

(b) 刃峰頂點附近的局部逆繞射傳播因子分布圖圖7 2WPE逆算法的數值結果(天線高50 m,xe=8 km,f=900 MHz)

表2 不同刃峰位置對應的單刃峰定位結果

上述討論的是單刃峰定位問題,對多刃峰環境,定位誤差會增大,定位精度會下降,這是因為接收測量點處會接收到多個刃峰的后向傳播場疊加的總后向場,在應用逆繞射傳播運算對每個刃峰進行定位的過程中,多個后向場會彼此互為干擾,從而導致各個刃峰的定位誤差增大,此時的定位精度不但受刃峰之間的距離位置關系影響,還和刃峰之間的高度關系密切相關.作者將在后續工作中對此進行研究和討論.

3結論

在逆繞射單向拋物方程法的基礎上,提出并研究了雙向拋物方程法的逆算法,以窄角雙向拋物方程為例推導了逆算法的公式,給出了雙向拋物方程法的逆算法的實現過程與步驟,并用該算法分析了單刃峰的探測與定位問題. 由數值結果可知,采用雙向拋物方程法的逆算法可以對電波傳播路徑上的障礙物的位置和高度進行探測和確定. 分析了天線高度變化和刃峰位置變化對定位精度的影響,數值結果表明,隨著天線高度增加,定位誤差略有增大;而隨著刃峰位置變化,當接收測量點與刃峰位置之間距離增加時,定位誤差增大. 綜合可知,雙向拋物方程法的逆算法可用于障礙物的探測和定位,精度較高,所以該算法在地理環境探測、各種目標定位等方法中將會有較好的應用前景.

參考文獻

[1] LEONTOVICH M A, FOCK V A. Solution of the problem of propagation of electromagnetic waves along the Earth’s surface by method of parabolic equations [J]. Journal Phys, USSR, 1946,10(1):13-23.

[2] DOCKERY G D. Modeling electromagnetic wave propagation in the troposphere using the parabolic equation [J]. IEEE Trans Antennas Propagation, 1988, 10: 1464-1470.

[3] KUTTLER J R, DOCKERY G D. Theoretical description of the parabolic approximation/Fourier split-step method of representing electromagnetic propagation in the troposphere [J]. Radio Science,1991,26(2):381-393.

[4] BARRIOS A E. A terrain parabolic equation model for propagation in the troposphere[J]. IEEE Trans Antennas Propagation, 1994, 42: 90-98.

[5] DOCKERY G D, KUTTLER J R. An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation [J]. IEEE Trans Antennas Propagation, 1996,44(12): 1592-1599.

[6] DONOHUE D, KUTTLER J. Propagation modeling over terrain using the parabolic wave equation [J]. IEEE Trans Antennas Propag, 2000, 48(2): 260-277.

[7] KUTTLER J, JANASWAMY R. Improved Fourier transform methods for solving the parabolic wave equation [J]. Radio Science, 2002, 37(2): 5.1-5.11.

[8] LEVY M. Parabolic Equation Methods for Electromagnetic Wave Propagation [M]. London: IEE, 2000.

[9]胡繪斌,柴舜連,毛鈞杰.寬角拋物方程在阻抗邊界條件下的應用[J]. 電波科學學報, 2005, 20(4):500-504.

HU Huibin, CHAI Shunlian, MAO Junjie. Application of the wideangle parabolic equation under impedance boundary condition[J]. Chinese Journal of Radio Science, 2005, 20(4): 500-504. (in Chinese)

[10]郭建炎, 王劍瑩, 龍云亮, 等. 基于拋物方程法的部分森林覆蓋山區電波傳播分析[J].電波科學學報, 2008, 23(6): 1045-1050.

GUO Jianyan, WANG Jianying, LONG Yunliang, et al. Analysis of radio propagation in partly forested terrain environment using parabolic equation approach [J]. Chinese Journal of Radio Science, 2008, 23(6):1045-1050. (in Chinese)

[11]ORAIZI H, HOSSEINZADEH S. Radio-wave-propagation modeling in the presence of multiple knife edges by the bidirectional parabolic-equation method [J]. IEEE Trans Vehicular Technology, 2007, 56 (3): 1033-1040.

[12]OZGUN O. Recursive two-way parabolic equation approach for modeling terrain effects in tropospheric propagation[J]. IEEE Trans Antennas Propagation, 2009, 57 (9): 2706-2713.

[13]APAYDIN G, SEVGI L. The split step Fourier and finite element based parabolic equation propagation prediction tools: canonical tests, systematic comparisons, and calibration [J]. IEEE Antennas and Propagation Magazine, 2010, 52(3): 66-79.

[14]王昆, 楊永欽, 龍云亮, 等. 多刃峰環境無線電波傳播預測的雙向拋物方程法[J]. 電波科學學報, 2011, 26(6): 1058-1064.

WANG Kun, YANG Yongqin, LONG Yunliang, et al. Twoway parabolic equation approach for modeling radio wave propagation in the presence of multiple knife edge[J]. Chinese Journal of Radio Science, 2011, 26(6): 1058-1064. (in Chinese)

[15]WANG Kun, LONG Yunliang. Propagation modeling over irregular terrain by the improved two-way parabolic equation method [J]. IEEE Trans Antennas Propagation, 2012, 60(9): 4467-4471.

[16]WANG Kun, GUO Jianyan, LONG Yunliang. Two-way wide-angle parabolic equation method for modeling electromagnetic waves propagation in urban areas [C]//The International Conference on Automatic Control and Artificial Intelligence. London: IET Press, 4010-4013.

[17]SPENCER T A, WALKER R A, HAWKES R M. Inverse diffraction parabolic wave equation localisation system (IDPELS) [C]//The 2004 International Symposium on GNSS/GPS. Sydney, 2004.

[18]SPENCER T A, WALKER R A, HAWKES R M. Inverse diffraction parabolic wave equation localisation system (IDPELS) [J]. Journal of Global Positioning Systems, 2005, 4(1/2): 245-257.

[19]HAWKES R M, SPENCER T A, WALKER R A. Three dimensional models for propagation in the troposphere and inverse diffraction [C]//IPS Radio and Space Services Proceedings, 2006: 6.1-6.6.

[20]GUO Jianyan, JIANG Hongyan, LONG Yunliang. Localization of transmitter in forest environments using inverse diffraction parabolic equation [C]//Asia-Pacific Microwave Conference, 2008.

[21]LI Dexin, LI Bo, WANG Xingbo. Passive localization of emitter source using inverse diffraction parabolic equation [C]//IEEE 10th International Conference on Signal Processing (ICSP), 2010: 111-114

[22]李德鑫, 楊日杰, 王鴻吉, 等. 基于逆繞射拋物方程法的輻射源定位技術研究[J]. 電波科學學報, 2011, 26(4): 683-687.

LI Dexin, YANG Rijie, WANG Hongji, et al. Passive location based on inverse diffraction parabolic equation[J]. Chinese Journal of Radio Science, 2011, 26(4): 683-687. (in Chinese)

王昆(1977-),女,遼寧人,博士,廣東工業大學信息工程學院講師,研究方向為無線通信、電波傳播理論和電磁數值計算等.

龍云亮(1963-),男,重慶人,中山大學教授,博士生導師,研究方向為天線理論與設計、電波傳播理論和電磁數值計算等.

劉震宇(1976-),男,湖南人,廣東工業大學信息工程學院副研究員,研究方向為物聯網技術、通信網信息安全處理和數字信號處理.

樊振宏, 譚延君, 陳如山. 一種旋轉對稱涂覆導體電磁散射高效分析方法[J]. 電波科學學報,2015,30(6):1116-1122. doi: 10.13443/j.cjors. 2014121701

FAN Zhenhong, TAN Yanjun, CHEN Rushan. Efficient approach for electromagnetic scattering by coated conducting bodies of revolution [J]. Chinese Journal of Radio Science,2015,30(6):1116-1122. (in Chinese). doi: 10.13443/j.cjors. 2014121701

Localization of obstacles based on the inverse algorithm of

two-way parabolic equation approach

WANG Kun1LONG Yunliang2LIU Zhenyu1

(1.SchoolofInformationEngineering,GuangdongUniversityofTechnology,

Guangzhou510006,China;2.DeptofElectronicsandCommunication

Engineering,SunYat-SenUniversity,Guangzhou510006,China)

AbstractIn this paper, an inverse algorithm of two-way parabolic equation (2WPE) method is presented and investigated on the basis of the inverse algorithm of one-way parabolic equation (1WPE) method, i.e. inverse diffraction parabolic equation method. At a certain receiving point, the backward-fields reflected by obstacles are obtained and inversely propagated, and then the location and height of the obstacles can be determined by the inverse diffraction fields. Simulation results of the inverse algorithm of 2WPE method for the scenario of single knife-edge is provided. For single knife-edge, the effects of the antenna height and the knife-edge location on the positioning accuracy is analyzed, and the accuracy is very high.

Key wordstwo-way parabolic equation; inverse algorithm; localization; knife-edges

作者簡介

收稿日期:2014-12-22

中圖分類號TN011

文獻標志碼A

文章編號1005-0388(2015)06-1108-08

主站蜘蛛池模板: 国产精品不卡片视频免费观看| 亚洲国产欧美自拍| 极品国产一区二区三区| 国产网站一区二区三区| 中文字幕伦视频| 亚洲午夜福利在线| 亚洲欧美另类中文字幕| 国产成人综合久久| 亚洲精品中文字幕午夜 | 91成人在线观看视频| 亚洲日本中文字幕天堂网| 性69交片免费看| 最新国产高清在线| 亚洲国产成人久久精品软件| www.av男人.com| 国产福利大秀91| 91破解版在线亚洲| 中日无码在线观看| 亚洲国产成人久久精品软件| 亚洲,国产,日韩,综合一区 | 亚洲欧美日韩另类在线一| 国产波多野结衣中文在线播放| 刘亦菲一区二区在线观看| 亚洲综合婷婷激情| 99久久性生片| 亚洲第一色网站| 99re这里只有国产中文精品国产精品 | 在线观看亚洲人成网站| 全免费a级毛片免费看不卡| 国产高清在线观看91精品| 蜜桃视频一区| 欧美国产日本高清不卡| 国产专区综合另类日韩一区 | 亚洲AⅤ综合在线欧美一区| 一个色综合久久| 在线播放精品一区二区啪视频 | 亚洲精品日产AⅤ| 91无码人妻精品一区| 最新国产精品第1页| 热99精品视频| 黄色一及毛片| 日韩欧美中文| 精品午夜国产福利观看| 97超碰精品成人国产| 99久久免费精品特色大片| 小说区 亚洲 自拍 另类| 18禁色诱爆乳网站| 亚洲成人播放| 欧美性精品不卡在线观看| 色综合激情网| 亚洲国产亚综合在线区| 亚洲国产精品日韩欧美一区| 国产99视频在线| 精品视频一区在线观看| 日韩天堂网| 亚洲欧洲国产成人综合不卡| 免费毛片a| 最新亚洲人成网站在线观看| 国产在线第二页| 久久成人免费| 久久精品人妻中文系列| 亚洲一区二区在线无码| 99精品在线看| 国产69囗曝护士吞精在线视频| 亚洲中文在线视频| 亚洲欧美精品一中文字幕| 国产精品视频白浆免费视频| 国产成人免费手机在线观看视频 | h网站在线播放| 在线色国产| 亚洲无线一二三四区男男| 日韩高清成人| 精品无码一区二区三区在线视频| 日韩在线观看网站| 香蕉综合在线视频91| 亚洲精品天堂自在久久77| 久久免费视频播放| 国产免费观看av大片的网站| a在线观看免费| 无遮挡国产高潮视频免费观看 | 国产精品无码翘臀在线看纯欲 | 深夜福利视频一区二区|