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

含曲面介質結構復雜目標電磁散射計算的射線追蹤方法

2016-11-29 03:44:34張磊侯兆國董純柱王超殷紅成
電波科學學報 2016年3期
關鍵詞:方法模型

張磊 侯兆國 董純柱 王超 殷紅成

(1.電磁散射重點實驗室,北京 100854;2.中國航天科工二院,北京 100854)

?

含曲面介質結構復雜目標電磁散射計算的射線追蹤方法

張磊1,2侯兆國1董純柱1王超1殷紅成1,2

(1.電磁散射重點實驗室,北京 100854;2.中國航天科工二院,北京 100854)

針對含曲面介質結構的電大復雜目標電磁散射計算問題,提出一種基于平面元網格模型曲率重構與射線密度歸一化概念相結合的快速射線追蹤方法.該方法通過曲率重構計算復雜目標表面的主曲率半徑,考慮從光疏介質到光密介質和從光密介質到光疏介質時電磁波照射凹凸曲面所具有的不同擴散或聚焦效應,并利用射線密度歸一化計算射線追蹤過程中每一根射線對總散射場的貢獻.當射線與介質表面的碰撞點位于焦散處時,通過引入功率追蹤成功克服了傳統幾何射線管在焦散處的奇異性.仿真結果驗證了該方法的正確性和高效性.

電磁散射;射線追蹤;曲面介質結構;射線密度歸一化;雷達散射截面

DOI 10.13443/j.cjors.2015072702

引 言

曲面介質結構廣泛存在于飛機、導彈等復雜目標,研究其電磁散射建模方法對目標低散射設計、目標特征提取與識別等應用具有重要意義.

傳統的數值方法,如矩量法(Method of Moment, MoM)、多層快速多極子算法(Multilevel Fast Multi-pole Algorithm,MLFMA)、時域有限差分(Finite Difference Time Domain,FDTD)方法等,通過求解積分方程或微分方程,能夠獲得含曲面介質結構復雜目標散射問題的高精度計算結果[1-10].然而,曲面介質結構剖分尺寸通常為介質波長的1/20,且需進行體剖分,隨著目標電尺寸的增大,未知量數目將變得異常龐大,使得計算資源和效率成為限制數值方法實際應用的瓶頸問題.

高頻漸近方法從電磁散射機理出發,利用射線追蹤技術快速獲取電磁波的傳播路徑和場分布,已成功用于高效分析平面結構或簡單曲面結構介質目標的電磁散射特性,且具有較高的計算精度[11-17].文獻[10]利用改進的幾何光學(Geometrical Optics, GO)法計算電大均勻介質目標的電磁散射特性;文獻[11]提出了一種射線追蹤方法計算電大非均勻介質目標的電磁散射特性;文獻[13]提出了一種射線密度歸一化(Ray-Density Normalization,RDN)[14]與物理光學(Physical Optics,PO)法相結合的混合算法計算均勻介質目標的電磁散射特性;文獻[15]提出了一種改進的等效流近似算法求解介質表面的等效電磁流,結合積分方程計算目標的電磁散射特性;文獻[16]利用PO與射線追蹤方法分析了頻率選擇表面(Frequency Selective Surface,FSS)的電磁特性.上述研究很少涉及含曲面介質結構的電大復雜目標電磁散射計算.

為準確高效地計算含曲面介質結構的電大復雜目標高頻電磁散射問題,需要解決曲面介質結構上射線反射與折射時引入的擴散或聚焦效應,以及射線與介質分界面的碰撞點位于焦散處時的奇異性等問題.為此,本文提出一種基于平面元網格模型曲率重構與RDN概念相結合的快速射線追蹤方法.通過將目標平面元網格模型擬合為局部二次曲面,再根據微分幾何公式和線性插值算法計算目標表面上任一點的主曲率半徑來完成曲率重構.與此同時,對比分析了從光疏介質到光密介質和從光密介質到光疏介質時電磁波分別照射凹凸曲面所具有的擴散或聚焦效應.區別于傳統的幾何射線管方法,本文方法利用RDN計算每一根射線對總散射場的貢獻,不僅降低了射線總數,而且當射線與介質分界面的碰撞點位于焦散處使得傳統的幾何射線管方法失效時,通過引入功率追蹤的思想成功解決了焦散處的奇異性問題.

1 平面元曲率重構

在復雜目標電磁散射建模過程中,為提高擬合精度并實現高效計算,曲面結構通常被描述為三角平面元網格模型,這種離散形式缺少連續的法向矢量和曲率.當利用射線追蹤技術計算介質目標多次反射或折射貢獻時,需要利用介質表面的曲率信息計算電磁波的擴散系數.因此,在對含曲面介質結構復雜目標的電磁建模過程中,需要快速、精確地求解平面元網格模型各點的曲率信息.

1.1 頂點曲率重構

圖1所示為頂點vi的鄰域,圖中與vi共面的頂點集合記為Vi,Vi中頂點個數記為|Vi|[18].頂點曲率重構的具體步驟如下:

圖1 頂點vi的鄰域

首先,計算頂點vi的法向矢量,方法是將與其相鄰的三角面元的法向矢量進行面積和夾角的加權平均[18].

然后,利用頂點vi及其鄰域中其他頂點坐標信息,求解式(1)中二次曲面的多項式系數.分析可知|Vi|≥3,當|Vi|>3時,通過最小二乘擬合得到最優的二次曲面:

(1)

最后,根據微分幾何公式計算頂點vi的主曲率1/ρ1和1/ρ2為[19]

(2)

為驗證頂點曲率重構算法的計算精度,表1給出了在不同剖分尺寸下計算得到的直徑1 m球的頂點主曲率均值,以及與理論真值(2 rad/m)的相對誤差.計算結果表明,頂點曲率重構算法具有較高的計算精度,且隨著剖分尺寸的減小,重構精度不斷提高.

表1 直徑1 m的球經重構得到的主曲率均值與相對誤差

1.2 擴散和聚焦效應

電磁波照射實際復雜介質目標,可分為從光疏介質入射到光密介質和從光密介質入射到光疏介質兩種狀態,而每種狀態下介質分界面又包含凹凸曲面兩種情況,對應的擴散和聚焦效應各不相同.圖2所示為凹凸曲面的反射與折射示意圖,圖中ε1、μ1為光疏介質的電磁參數,ε2、μ2為光密介質的電磁參數,即ε2μ2>ε1μ1,曲面法向方向由光密介質指向光疏介質.

(a) 從光疏介質入射到光密介質

(b) 從光密介質入射到光疏介質圖2 凹凸曲面的反射與折射示意圖

表2給出了圖2中四種狀態下反射與折射的擴散和聚焦效應,可知擴散和聚焦效應與曲面凹凸性、介質材料分布密切相關.

表2 四種狀態下反射與折射的擴散和聚焦效應

圖3 凸曲面邊界上的射線追蹤示意圖

P0點處的反射波波前和折射波波前的曲率矩陣分別為[20]

(3)

(4)

式中,h=k1cos θi-k2cos θt,k1、k2分別為入射波和折射波所在空間的波數.

根據曲率矩陣計算反射波波前與折射波波前的主曲率半徑,可表示為

(5)

距離P0點s處反射波或折射波的擴散系數為

(6)

式中,ρ1、ρ2分別為反射波波前或折射波波前的主曲率半徑.

2 RDN概念

Didascalou[14]假定每一條實際的射線傳播路徑上存在著許多相關聯的射線,通過確定該路徑上被接收面接收的關聯射線的總數,可對該實際射線對總散射場的貢獻進行歸一化.定義射線密度為單位面積上的射線數,則被接收的射線總數等于射線密度與接收面在傳播方向上的投影面積的乘積.因此,每根射線不僅具有幅度、相位、極化等信息,還攜帶射線密度信息.

2.1 射線密度

平面波入射條件下,射線密度的初始值為

(7)

針對曲面介質表面,反射和折射射線發生擴散和聚焦,根據幾何光學理論,反射與折射射線密度可分別表示為:

(8)

(9)

2.2 場追蹤

在射線追蹤過程中,賦予每根射線一個初始的電場強度,計算射線在傳輸路徑上的場分布,直至射線離開目標或電場能量衰減至足夠小后停止,該過程即為場追蹤.

P0點的反射場和折射場分別表示為:

exp(-jk1sr);

(10)

exp(-jk2st) .

(11)

式中: Ei(P0)為P0點處入射場; Er(sr)為距離P0點sr處的反射場; Et(st)為距離P0點st處的折射場; R(P0)為并矢反射系數; T(P0)為并矢折射系數[20].當反射或折射射線再次與介質分界面發生碰撞,反射場或折射場變為新的入射場.

場追蹤完成后,總散射場求解需要對各面元的貢獻進行疊加.基于上述RDN概念[14],利用各自面元上接收到的射線總數歸一化每根射線對總散射場的貢獻,具體為:

ERDN(r′)=XFE(r′) ;

(12)

HRDN(r′)=XFH(r′) ;

(13)

(14)

式中: E(r′)和H(r′)為一條實際傳播路徑在面元上的總場; ERDN(r′)和HRDN(r′)為歸一化后的總場; XF為歸一化加權因子; M為面元上接收到的射線總數; A為面元的面積.

當射線追蹤過程中碰撞點位于焦散處時,擴散系數和射線密度均趨于無窮大,式(10)、式(11)計算得到的電場為無窮大,與實際情況不符,場追蹤失效.對此,引入功率追蹤的思想,將焦散處的場追蹤轉變為功率追蹤.

2.3 功率追蹤

功率追蹤的思路是賦予每根射線一個初始的功率取代電場強度,其余與場追蹤類似[14].兩者區別為:功率追蹤在傳輸過程中不考慮擴散系數,因為除了有耗介質的傳輸衰減外,功率保持不變.

功率追蹤過程中反射場和折射場可表示為:

Er(sr)=R(P0)Ei(P0)exp(-jk1sr) ;

(15)

Et(st)=T(P0)Ei(P0)exp(-jk2st).

(16)

基于上述RDN概念,功率追蹤的歸一化加權系數為[14]

(17)

式中, TD為場追蹤過程中的擴散系數.

對比式(14)和式(17)可知,射線追蹤過程中無碰撞點位于焦散處時,場追蹤與功率追蹤計算結果完全相同.反之,場追蹤失效.基于射線總數量有限的前提,功率追蹤的歸一化加權系數可進一步表示為

(18)

對此,功率追蹤為幾何光學場追蹤在焦散處失效的問題提供了一種有效的近似解決方案,避免了焦散處的奇異性.需要強調的是,功率追蹤僅作為場追蹤在焦散處失效時的補充,當碰撞點不處于焦散處時,以場追蹤為主.

3 仿真示例與分析

驗證本文提出的射線追蹤方法的正確性和高效性,以均勻介質球、介質球頭柱和某導彈模型為例開展仿真計算,并分別與Mie級數解、傳統的幾何射線管方法和MLFMA計算結果進行對比.

3.1 均勻介質球

均勻介質球幾何模型如圖4所示,半徑為0.5 m,電磁參數為εr=3.7-j0.1,μr=1.0.計算條件為:頻率f=2~10GHz,步長Δf=10MHz,垂直極化.圖5給出了射線追蹤方法計算的介質球后向散射RCS與Mie級數解對比.

圖4 均勻介質球幾何模型

圖5 射線追蹤方法與Mie級數解計算的后向散射RCS

由圖5可知,文中提出的快速射線追蹤方法計算結果與Mie級數解吻合良好,尤其隨著頻率的不斷增大,兩種方法計算得到的RCS趨于一致.因此,針對電大曲面介質體目標,本文方法具有較高的計算精度.

3.2 球頭柱模型

球頭柱模型及其射線追蹤過程如圖6所示,其中半球部分為均勻損耗介質體,半徑為0.2m,電磁參數為εr=3.8-j0.24,μr=1.0;圓柱部分為理想導體,高為1.0m,半徑為0.2m.計算條件為:俯仰角為0°(與Z軸夾角為90°),方位角為0°~180°,頻率f=2.5GHz,平行極化.圖7給出了利用本文射線追蹤方法、傳統幾何射線管方法、以及MLFMA計算的球頭柱模型后向散射RCS對比.

圖6 球頭柱模型及其射線追蹤示意圖

圖7 射線追蹤方法與MLFMA計算的后向散射RCS對比

由圖6可知,介質半球的射線追蹤過程遠比導體圓柱復雜,入射射線在P1點處發生反射和折射,折射射線在半球與圓柱的分界面上P2點處僅發生反射,反射射線再依次在P3、P4、P5點處發生反射與折射,而照射到導體圓柱上的電磁波,在圓柱表面Q1、Q2處僅發生反射.圖中射線的粗細變化表示射線攜帶能量的衰減情況.

由圖7可知,本文射線追蹤方法計算的球頭柱模型后向散射RCS曲線與傳統幾何射線管方法和MLFMA計算結果在0°~180°方位角范圍內均吻合較好.同時,圖7還給出了不考慮曲面介質結構擴散系數時射線追蹤方法的計算結果,與MLFMA計算結果對比可知,在0°~90°方位角范圍內存在較大偏差,其原因在于忽略了三角面元網格模型的局部幾何特征,即擴散或聚焦效應.

3.3 導彈仿真模型

為進一步驗證本文射線追蹤方法計算含曲面介質結構復雜目標時的精度和效率,以導彈仿真模型為例開展仿真計算.幾何模型如圖8所示,彈頭為均勻損耗介質,電磁參數為εr=3.8-j0.24,μr=1.0,彈體與尾翼為理想導體,整體長為51.2λ0,平面波水平方向入射,方位角為0°~180°,平行極化.圖9給出了射線追蹤方法與MLFMA計算的后向散射RCS對比,其中縱坐標為后向散射RCS與波長平方的比值.

圖8 導彈仿真模型

圖9 射線追蹤方法與MLFMA計算的后向散射RCS

由圖9可知,射線追蹤方法與MLFMA計算得到的導彈仿真模型后向散射RCS在20°~90°方位角范圍內吻合較好,而在0°~20°方位角范圍內存在一定誤差,主要是因為在該方位角范圍內為掠入射,且未考慮行波等貢獻,高頻漸近方法難以獲得準確結果.類似的問題出現在100°~120°方位角范圍內.

最后,為驗證本文方法的高效性,定量對比了射線追蹤方法與MLFMA的計算效率和資源需求.表3給出了針對該模型兩種方法所需的計算時間和峰值內存使用情況,其中MLFMA采用12核并行計算,射線追蹤方法僅采用1核計算.分析可知,射線追蹤方法的計算效率為MLFMA的93.6倍,使用的峰值內存僅為MLFMA的1.25%,可直接在普通計算機上開展仿真計算,無需配置高性能計算機工作站.

表3 兩種方法的計算時間和峰值內存使用

4 結 論

本文提出了基于三角平面元網格模型曲率重構與RDN概念相結合的快速射線追蹤方法,實現了對含曲面介質結構電大復雜目標電磁散射特性的高效、準確計算;同時,通過引入功率追蹤思想,避免了碰撞點位于焦散處時的奇異性,為幾何光學場追蹤在焦散處失效的問題提供了一種有效的近似解決方案.仿真結果驗證了本文方法的正確性和高效性.

[1]UMASHANKARKR,TAFLOVEA,RAOSM.Electro-magneticscatteringbyarbitraryshapedthree-dimensionalhomogeneouslossydielectricobject[J].IEEEtransactionsonantennasandpropagation, 1986, 34: 758-766.

[2]RAOSM,CHACC,CRAVEYRL,etal.Electro-magneticscatteringfromarbitraryshapedconductingbodiescoatedwithlossymaterialsofarbitrarythickness[J].IEEEtransactionsonantennasandpropagation, 1991, 39(5): 627-631.

[3]PETERSONAF.EfficientsolenoidaldiscretizationofthevolumeEFIEforelectromagneticscatteringfromdielectricobjects[J].IEEEtransactionsonantennasandpropagation, 2014, 62(3): 1475-1478.

[4]UMASHANKARK,TAFLOVEA.Anovelmethodtoanalyzeelectromagneticscatteringofcomplexobjects[J].IEEEtransactionsonelectromagneticcompatibility, 1982, 24(4): 397-405.

[5]DEMARESTK,PLUMBR,HUANGZhubo.FDTDmodelingofscatterersinstratifiedmedia[J].IEEEtransactionsonantennasandpropagation, 43(10), 1995: 1164-1168.

[6]LUCC.Afastalgorithmbasedonvolumeintegralequationforanalysisofarbitrarilyshapeddielectricradomes[J].IEEEtransactionsonantennasandpropagation, 2003, 51(3): 606-612.

[7]YANGPG,NIEZP,QUEXF,etal.ApplicationofCFIEwithhierarchicalvectorbasisfunctionsforEMscatteringfromhybridmetal-dielectricobjectswithMLFMA[C]//2012 10thInternationalSymposiumonAntennas,Propagation&EMTheory, 2012: 967-970.

[8]ERGULO.Fastandaccuratesolutionsofelectro-magneticsproblemsinvolvinglossydielectricobjectswiththemultilevelfastmultipolealgorithm[J].Engineeringanalysiswithboundaryelements, 2012, 36: 423-432.

[9] 聶在平, 胡俊, 闕肖峰, 等. 面向工程應用能力提升的電磁散射高效數值分析: 進展與挑戰[J]. 電波科學學報, 2014, 29(1): 1-11.

NIEZP,HUJ,QUEXF,etal.ApplicationcapabilitypromotionorientedefficientnumericalanalysisofEMscattering:progressesandchallenges[J].Chinesejournalofradioscience, 2014, 29(1): 1-11. (inChinese)

[10] 闕肖峰, 聶在平, 胡俊. 導體介質組合體電磁分析的建模與計算[J]. 電波科學學報, 2008, 23(3): 396-401.

QUEXF,NIEZP,HUJ.Electromagneticmodelingandcalculationforcompositeconductinganddielectricobjects[J].Chinesejournalofradioscience, 2008, 23(3): 396-401. (inChinese)

[11]KOUYOUMJIANRG,PETERSL,THOMASD.Amodifiedgeometricalopticmethodforscatteringbydielectricbodies[J].IEEEtransactionsonantennasandpropagation, 1963, 11(6): 690-703.

[12]KIMH,LINGH.Electromagneticscatteringfromaninhomogeneousobjectbyraytracing[J].IEEEtransactionsonantennasandpropagation, 1992, 40(5): 517-525.

[13]WEINMANNF.PO/PTDraytracingforarbitrarymetallicanddielectricobjects[C]//Europeanconferenceonantennas&propagation, 2006: 6-10.

[14]DIDASCALOUD,SCHAFERTM,WEINMANNF,etal.Ray-densitynormalizationforray-opticalwavepropagationmodelinginarbitrarilyshapedtunnels[J].IEEEtransactionsonantennasandpropagation, 2000, 48(9): 1316-1325.

[15]MEANAJG,MARTINEZ-LORENZOJA,LAS-HERASF,etal.Wavescatteringbydielectricandlossymaterialsusingthemodifiedequivalentcurrentapproximation(MECA)[J].IEEEtransactionsonantennasandpropagation, 2010, 58(11): 3757-3761.

[16]KIMJH,CHUNHJ,HONGIP,etal.AnalysisofFSSradomesbasedonphysicalopticsmethodandraytracingtechnique[J].IEEEantennasandwirelesspropagationletters, 2014, 13: 868-871.

[17] 張磊, 王超, 董純柱, 等. 分層介質目標電磁散射計算的快速射線追蹤方法[J]. 電波科學學報, 2015, 30(2): 208-216.

ZHANGL,WANGC,DONGCZ,etal.Fastray-tracingmethodforelectromagneticscatteringcomputationofmulti-layereddielectricstructure[J].Chinesejournalofradioscience, 2015, 30(2): 208-216. (inChinese)

[18] 神會存, 李建華, 周來水. 三角網格模型頂點法矢與離散曲率計算[J]. 計算機工程與應用, 2005, 26: 12-15.

SHENHC,LIJH,ZHOULS.Estimationoftriangularmeshvertexnormalvectoranddiscretecurvature[J].Computerengineeringandapplications, 2005, 26: 12-15.(inChinese)

[19] 汪茂光. 幾何繞射理論[M]. 西安電子科技大學出版社, 1994: 95-104.

[20]JAMESGL.Geometricaltheoryofdiffractionforelectromagneticwaves[M].TheInstitutionofEngineeringandTechnology, 1979.

張磊 (1989-),男,安徽人,現為中國航天科工二院研究生院在讀博士生,主要研究方向為目標與環境電磁散射計算方法及計算機加速技術.

侯兆國 (1983-),男,甘肅人,高級工程師,博士,主要研究方向為雷達目標特性建模與模型評估.

董純柱 (1981-),男,河南人,高級工程師,博士,主要研究方向為雷達目標特性建模與特征提取.王超 (1979-),男,陜西人,高級工程師,博士,主要研究方向為雷達目標特性建模與特征分析.

殷紅成 (1967-),男,江西人,研究員,中國航天科工二院研究生院博士生導師,主要研究方向為雷達目標特性研究與應用.

Ray tracing method for electromagnetic scattering computation from complex target with curved surface dielectric structure

ZHANG Lei1,2HOU Zhaoguo1DONG Chunzhu1WANG Chao1YIN Hongcheng1,2

(1.Science and Technology on Electromagnetic Scattering Laboratory, Beijing 100854, China;2.TheSecondAcademyofChinaAerospaceScienceandIndustryCorporation,Beijing100854,China)

A ray tracing method combining the curvature reconstruction of faceted model with the ray-density normalization (RDN) is presented for electromagnetic(EM) scattering computation from electrically large complex target with curved surface dielectric structure. This method calculates the principal radii of curvature by curvature reconstruction, considers the different divergence or convergence effect for convex or concave surface illuminated by EM waves from optically thin medium to optically dense medium and from optically dense medium to optically thin medium, and calculates the contribution of each ray to the total scattered field by using RDN. The method also overcomes the singularity of traditional ray tubes at caustics by introducing the power tracing when the impact point between the ray and the interface is located at caustics. Simulation results verify the accuracy and efficiency of the proposed method.

EM scattering; ray tracing; curved surface dielectric structure; ray-density normalization (RDN); radar cross section (RCS)

10.13443/j.cjors.2015072702

2015-07-27

國家自然科學基金重大項目(No.61490690)

TN011

A

1005-0388(2016)03-0546-08

張磊, 侯兆國, 董純柱, 等. 含曲面介質結構復雜目標電磁散射計算的射線追蹤方法[J]. 電波科學學報,2016,31(3):546-552+596.

ZHANG L, HOU Z G, DONG C Z, et al. Ray tracing method for electromagnetic scattering computation from complex target with curved surface dielectric structure[J]. Chinese journal of radio science, 2016,31(3):546-552+596.(in Chinese). DOI:10.13443/j.cjors.2015072702

聯系人: 張磊 E-mail: zhangleigcss@126.com

猜你喜歡
方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
學習方法
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲大尺码专区影院| 国产中文一区二区苍井空| V一区无码内射国产| 国产日韩欧美一区二区三区在线 | 亚洲欧美人成人让影院| 国产日韩久久久久无码精品| 国产国语一级毛片在线视频| 精品国产自| 亚洲人成网站观看在线观看| 亚洲午夜天堂| 午夜不卡福利| 亚洲国产精品一区二区第一页免| 久久久亚洲国产美女国产盗摄| 国产91无码福利在线| 欧美成人二区| 99热这里只有精品2| 一级爱做片免费观看久久| 无码人中文字幕| 男人的天堂久久精品激情| 欧美成人综合视频| 手机精品视频在线观看免费| 114级毛片免费观看| 成人欧美日韩| 亚洲国产看片基地久久1024| 99九九成人免费视频精品 | 在线观看无码av五月花| 精品人妻无码中字系列| 亚洲欧美成人影院| 又黄又湿又爽的视频| 人妻精品全国免费视频| 亚洲国产精品无码AV| 国产电话自拍伊人| 91在线一9|永久视频在线| 3D动漫精品啪啪一区二区下载| 女人18毛片一级毛片在线 | 99精品影院| 丁香婷婷激情网| 欧美成人一级| 日韩精品亚洲人旧成在线| 999国内精品视频免费| 伊人久久影视| 国产精品熟女亚洲AV麻豆| 国产精品9| 中文纯内无码H| 在线另类稀缺国产呦| 婷婷伊人久久| 欧美不卡视频一区发布| 中国国产一级毛片| 久久婷婷综合色一区二区| a天堂视频| 无码'专区第一页| 国产成人你懂的在线观看| 中文一区二区视频| 久热精品免费| 国产一二三区视频| 国产va在线| 91综合色区亚洲熟妇p| 国产区91| 2019年国产精品自拍不卡| 精品欧美一区二区三区久久久| 一级黄色欧美| 毛片三级在线观看| 国产区精品高清在线观看| 亚洲国产成人久久精品软件| 国产成人区在线观看视频| 色婷婷成人| 国内精品小视频福利网址| 日韩精品无码免费一区二区三区| 国产熟女一级毛片| 亚洲美女久久| 2020精品极品国产色在线观看 | 日韩a在线观看免费观看| 成人福利一区二区视频在线| 女人18毛片一级毛片在线 | 成人国内精品久久久久影院| 福利一区三区| 在线观看网站国产| 中文字幕欧美日韩高清| 久无码久无码av无码| 欧美日本激情| 亚洲欧美成人| 欧美在线国产|