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

合成孔徑雷達影像地形校正半經(jīng)驗?zāi)P?/h1>
2011-01-31 08:22:26鮑艷松王軍戰(zhàn)張友靜屈建軍張偉民
測繪學報 2011年4期
關(guān)鍵詞:模型

鮑艷松,王軍戰(zhàn),張友靜,屈建軍,張偉民,陳 晨

1.南京信息工程大學 氣象災(zāi)害省部共建教育部重點實驗室,江蘇南京210044;2.中國科學院 寒區(qū)旱區(qū)環(huán)境與工程研究所,甘肅蘭州730000;3.河海大學地球科學與工程學院,江蘇南京210098;4.江蘇省地質(zhì)調(diào)查院,江蘇南京210018

1 引 言

隨著遙感技術(shù)的發(fā)展,SAR(合成孔徑雷達)應(yīng)用越來越廣泛,如農(nóng)作物估產(chǎn)、土壤濕度反演、積雪參數(shù)反演、森林蓄積量反演、地表覆蓋分類等[1-4],但由地形因素引起的SAR后向散射系數(shù)變化嚴重影響了它的這些應(yīng)用。而在地形起伏地區(qū)本地入射角被認為是影響雷達后向散射的主要因素[5-7],因此,文獻[8—10]都建立了SAR圖像灰度值變化與本地入射角的余弦函數(shù)之間的簡單經(jīng)驗?zāi)P停孕U匦螌AR影像灰度值的影響。文獻[11]在總結(jié)前人建立的地形校正經(jīng)驗?zāi)P偷幕A(chǔ)上,利用雷達傳感器參數(shù)和DEM提取多個因子,建立圖像灰度值和多個因子間的7個新的地形校正經(jīng)驗?zāi)P停容^它們對不同地物覆蓋地表的應(yīng)用效果。文獻[12]將地形對SAR后向散射系數(shù)的影響分為兩個方面,一方面是地形變化引起面積變化而帶來的影響,另一方面是地形引起的本地入射角的變化而帶來的影響,并基于機載SAR數(shù)據(jù)推導(dǎo)去除地形對后向散射系數(shù)影響的方法。文獻[13—14]也將地形校正分為面積校正、輻射校正兩個方面,經(jīng)過嚴格地推導(dǎo)建立地形校正方法,取得較好的結(jié)果。文獻[15]建立以本地入射角為函數(shù)的線性模型校正后向散射系數(shù),并用RadarSat SAR影像對校正結(jié)果進行試驗,但這種簡單線性關(guān)系并不能充分說明地形起伏區(qū)域后向散射系數(shù)與本地入射角的復(fù)雜關(guān)系。

微波后向散射模型AIEM(advanced integrated equation model)模型是近年來模擬微波后向散射和土壤濕度反演的主要模型之一[16-18],該模型是IEM模型的發(fā)展與進一步完善[19],能夠很好地再現(xiàn)地表散射。因此本文首先基于理論知識,根據(jù)雷達側(cè)視(左視和右視)成像特點,分別簡單推導(dǎo)給出左視和右視兩種情況下影響SAR后向散射的重要因子——本地入射角計算模型;其次基于微波后向散射物理模型AIEM,模擬分析不同地形(坡度、坡向)對不同雷達入射角度后向散射的影響,建立以本地入射角為因子的地形校正半經(jīng)驗?zāi)P停阅M數(shù)據(jù)和ASAR數(shù)據(jù)作為試驗樣本,檢驗?zāi)P偷倪m用性;最后以地形校正前后的ASAR影像分類結(jié)果來檢驗?zāi)P偷膽?yīng)用效果。

2 模型推導(dǎo)

2.1 本地入射角計算

本地入射角是建立SAR影像地形校正模型的重要因子之一,因此本節(jié)首先給出左視和右視兩種情況下本地入射角的計算公式。

圖1 地形起伏地表SAR成像幾何示意圖(左視)Fig.1 SAR imaging geometry under the terrain undulations earth assumption(left-looking).

圖1中,HO為入射波方向;O′O為水準面垂線;OK為O′O的延長線;∠HOK即為雷達入射角θ;O Q′為坡面上點O的法線。本地入射角為入射波與地表法線夾角,如圖1中(假設(shè)為左視成像),OQ′為地表O點法線,則本地入射角β為∠HOQ′。本文中坡度角定義為:坡面與水平面的夾角,取值范圍為0°~90°,如圖1中的α;坡向定義為:地表面上一點的切平面的法線矢量在水平面的投影,與過該點的方位向逆向的夾角,順時針遞增,取值范圍為0°~360°,如圖1中的γ,而基于DEM利用ENVI提取的坡向為切平面的法線矢量在水平面的投影與過該點的正北向的夾角,因此需根據(jù)衛(wèi)星軌道傾角,將ENVI提取坡向轉(zhuǎn)換為本文定義的坡向。

假設(shè)已知雷達入射角θ,通過DEM計算出坡度角α,基于DEM提取的坡向和衛(wèi)星軌道傾角求出坡向γ。若求本地入射角β,利用立體幾何知識,先做以下輔助線:過底面垂足點O′分別作地表法線OQ′和入射向線段HO的平行線,分別交ON、OM、MH于點F、E、Q。最終簡化為如圖2所示。

圖2 地形起伏地表SAR成像幾何簡化圖(左視)Fig.2 Simplified SAR imaging geometry under the terrain undulations earth assumption(leftlooking).

圖2中,∠FO′E為本地入射角β。由圖2,利用角度關(guān)系、三角函數(shù)關(guān)系及正弦函數(shù)定理,可得到距離向坡度角的正切值(式(1))、方位向坡度角的正切值(式(2))、本地入射角余弦值計算公式(式(3))和本地入射角計算公式(式(4))

由于距離向坡度角、方位向坡度角都在0°~90°之間,它們的正切值始終為正,而此利用式(1)、式(2)計算時,由于γ角變化,它們的正切值可能為負,因此計算時取絕對值即可。

當雷達成像為右視時,同理可得距離向坡度角的正切值(式(5))、方位向坡度角的正切值(式(6))、本地入射角余弦值計算公式(式(7))和本地入射角計算公式(式(8))

同樣,利用式(5)、式(6)計算距離向坡度角正切值、方位向坡度角正切值時應(yīng)取絕對值。

至此,根據(jù)左視和右視兩種雷達成像情況,已經(jīng)簡單介紹給出了距離向坡度角正切值、方位向坡度角正切值、本地入射角余弦值和本地入射角的計算模型。由式(7)、式(8)可知,本地入射角主要受雷達入射角θ、地形坡度α、坡向γ以及左視或右視成像等因素的影響。下面將通過模型模擬分析本地入射角對后向散射的影響,并建立地形校正模型。

2.2 地形對后向散射影響模擬分析

利用AIEM模型模擬后向散射,需要多個輸入?yún)?shù),如土壤濕度、表示地表粗糙度的均方根高度和表面相關(guān)長度等[19],本文采用實測值作為AIEM模型的輸入?yún)?shù),參數(shù)范圍如表1。

表1 地表測量數(shù)據(jù)范圍Tab.1 Main characteristics of the ground truth data used in this study

假設(shè)為左視成像,則可利用式(4)分別計算雷達入射角度為24°、32°和40°時的本地入射角,式中坡度角分別設(shè)置為4°、8°、12°、16°、20°、24°,坡向分別設(shè)置為0°、45°、90°、135°、180°、225°、270°、315°。求取的本地入射角作為AIEM模型中的入射角參數(shù)。

基于AIEM模型,分別模擬雷達入射角度為24°、32°和40°時,不同坡度角、坡向計算得到的本地入射角的后向散射系數(shù)。其中土壤含水量值設(shè)為25%,均方根高度設(shè)為1.2cm,相關(guān)長度設(shè)為8.5cm。圖3為本地入射角與雷達后向散射系數(shù)散點圖。

由圖3可知,其他參數(shù)已定時,后向散射系數(shù)與本地入射角之間用二次函數(shù)關(guān)系擬合最佳,擬合結(jié)果如圖中所示。雷達入射角為24°時,坡度角從0°~24°變化、坡向從0°~360°變化,計算得到的本地入射角范圍為0°~48°,引起的后向散射系數(shù)變化約為-32dB;雷達入射角為32°時,計算得到的本地入射角范圍為8°~56°,引起的后向散射系數(shù)變化約為-43dB;雷達入射角為40°時,計算得到的本地入射角范圍為16°~64°,引起的后向散射系數(shù)變化約為-55dB。因此可以看出對于相對較小的雷達入射角,相同的坡度角、坡向變化引起的后向散射系數(shù)變化較小,這也從側(cè)面說明在地形起伏地表,相對較小的雷達入射角SAR數(shù)據(jù)是定量遙感如土壤濕度反演等的最佳數(shù)據(jù)源。右視時有相似的結(jié)論。因此,后向散射系數(shù)與本地入射角具有很好的二次函數(shù)關(guān)系,這也是本文建立SAR影像地形校正模型的依據(jù)。

圖3 不同雷達入射角下本地入射角與后向散射系數(shù)散點圖Fig.3 AIEM simulation of the backscattering at different slop and aspect

2.3 模型的建立

由2.2節(jié)模擬試驗可知,SAR后向散射系數(shù)與本地入射角之間用二次函數(shù)擬合最佳,表達式如下

式中,σg為地形變化時的SAR后向散射系數(shù);β為本地入射角;c為常數(shù)項;a、b為系數(shù)。a、b、c可由本地入射角和地形變化時的SAR后向散射系數(shù)值最小二乘擬合得到。

當?shù)匦螣o變化時,坡度角相當于0°,同樣滿足式(9),因此有

式中,σc為地形無變化時的SAR后向散射系數(shù);θ為雷達入射角或像元入射角;a、b、c同前。

式(9)除以式(10),即可得到計算地形無起伏時的后向散射系數(shù)σc的半經(jīng)驗?zāi)P?/p>

式中,σc為地形校正后的后向散射系數(shù);σg為地形變化時的后向散射系數(shù),即是獲取的SAR影像后向散射系數(shù);θ角可由已知的雷達入射角擬合插值求取;β角可由式(4)或式(8)求得;a、b、c同上,可由本地入射角和SAR影像后向散射系數(shù)最小二乘擬合得到。

所謂的半經(jīng)驗?zāi)P停侵咐碚撃P褪窃撃P退玫幕A(chǔ),而模型的一些系數(shù)則是由試驗測量數(shù)據(jù)擬合確定。式(11)基于物理理論模型AIEM模擬和嚴密推導(dǎo)得到,其系數(shù)可由實測DEM數(shù)據(jù)和SAR數(shù)據(jù)擬合得到,因此可稱為半經(jīng)驗?zāi)P汀T撃P托问胶唵巍?shù)少,且系數(shù)擬合容易,所以模型易于操作,下一節(jié)將分別以模擬數(shù)據(jù)和實測數(shù)據(jù)檢驗?zāi)P汀?/p>

3 模型的驗證與應(yīng)用

3.1 基于模擬數(shù)據(jù)的模型驗證

首先基于模擬數(shù)據(jù),利用式(11)對后向散射系數(shù)進行校正(雷達入射角為24°、32°、40°時,半經(jīng)驗?zāi)P椭械南禂?shù)和常數(shù)為圖3中擬合方程的系數(shù)和常數(shù)),校正結(jié)果如表2所示。同時利用常用的余弦模型[8]對模擬數(shù)據(jù)進行校正,結(jié)果如表2所示。均方根誤差為三個不同雷達入射角度下校正后的后向散射系數(shù)與無地形起伏時后向散射系數(shù)之間的計算結(jié)果。由表2可知,對于三個不同的雷達入射角,本地入射角半經(jīng)驗?zāi)P途让黠@好于余弦模型,且相對較小的雷達入射角校正效果更好。因此下面將用以本地入射角為因子的地形校正半經(jīng)驗?zāi)P蛯Λ@取的ASAR影像進行校正和應(yīng)用。

表2 地形校正結(jié)果比較Tab.2 Comparison of the correction results for the three incidence angles

3.2 基于ASAR數(shù)據(jù)的模型應(yīng)用

3.2.1 數(shù)據(jù)獲取與處理

本文使用的SAR數(shù)據(jù)是Envisat-1衛(wèi)星上搭載的高級雷達傳感器ASAR數(shù)據(jù),雷達工作頻率為5.36GHz,模式為可選擇極化精準模式(APP模式),其空間分辨率為30m。ASAR可選擇極化精準模式擁有7個入射角模式[20],本文獲取的ASAR影像為第二入射角模式(IS2,2005-05-08),入射角度為19.2°~26.7°,極化方式為HH/VV極化組合。在獲取的ASAR數(shù)據(jù)中選取試驗區(qū)如圖4所示。

圖4 研究區(qū)示意圖Fig.4 Location of the study area

本文獲取的ASAR-APP模式影像為1B級數(shù)據(jù),圖像記錄的灰度值(DN值)是振幅值,首先使用影像包含的經(jīng)緯度信息元數(shù)據(jù)對影像進行初步幾何糾正,然后配準到經(jīng)過精糾正的ETM+上,最后將灰度值(振幅值)轉(zhuǎn)化成后向散射系數(shù),計算如下

本文使用的DEM數(shù)據(jù)為NASA提供的2008年全球DEM數(shù)據(jù)(GDEM,global digital elevation model),是基于ASTER立體像對數(shù)據(jù),利用攝影測量技術(shù)生產(chǎn)的DEM,其分辨率為30m。將DEM數(shù)據(jù)精配準到ETM+影像,幾何校正后利用ENVI提取坡度角和坡向參數(shù),參數(shù)提取窗口設(shè)置為7×7[21]。然后根據(jù)ASAR的軌道傾角,將ENVI提取的坡向參數(shù)轉(zhuǎn)換為本文定義的坡向。

3.2.2 模型應(yīng)用

利用ASAR影像已知的入射角,擬合二次曲線,計算出每個像元的入射角,即得到半經(jīng)驗?zāi)P椭械南裨肷浣铅龋焕孟裨肷浣呛虯SAR傳感器運行軌道方向,即可判斷出傳感器成像時為右視,因此基于DEM獲取的坡度和計算得到的坡向,利用式(8)計算得到本地入射角β。然后利用本地入射角和原始影像后向散射系數(shù)的散點圖擬合半經(jīng)驗?zāi)P椭械南禂?shù),最后利用半經(jīng)驗?zāi)P屯瓿稍囼瀰^(qū)的地形校正。結(jié)果如圖5(b)所示,總體上看半經(jīng)驗?zāi)P托U^的影像后向散射系數(shù)值出現(xiàn)明顯好轉(zhuǎn)現(xiàn)象,由于地形原因造成的明暗對比變小,從同一區(qū)域局部放大圖中更容易看出。

圖5 影像地形校正前后比較Fig.5 ASAR image

為定量評價半經(jīng)驗?zāi)P偷匦涡UY(jié)果,選取地形校正前后影像方差較少的百分比作為定量評價因子[12,22],計算方法如式(13)。分別在未校正和已校正過的影像上,統(tǒng)計整個山體區(qū)域和隨機選取的10個30×30像元大小區(qū)域的方差,利用式(13)計算結(jié)果如表3所示為未校正影像方差;為地形校正后影像的方差。

表3 本文半經(jīng)驗?zāi)P偷匦涡UY(jié)果Tab.3 Correction results for ASAR image by the semi-model

由表3可知,使用本文地形校正半經(jīng)驗?zāi)P蛯ρ芯繀^(qū)進行校正,影像總體方差減小百分比為19%。方差減小百分比最小為11%,方差減小百分比最大達29%。

3.2.3 ASAR影像分類的結(jié)果比較

試驗區(qū)內(nèi)主要地物有林地、水體和建筑,結(jié)合同期的TM影像數(shù)據(jù),在地形校正前后的ASAR影像上選取三類地物相同的訓(xùn)練樣本,利用監(jiān)督分類中的最小距離分類法對地形校正前后的ASAR影像進行分類,然后利用選取的訓(xùn)練樣本對分類結(jié)果進行評價,可知地形校正后ASAR影像的分類總體精度達0.58,較地形校正前分類總體精度提高12個百分點。

不難看出無論地形校正前后,ASAR影像分類的總體精度都不高,主要有以下幾個原因:一是由于SAR影像受相干斑點噪聲的存在,嚴重影響地物分類精度[23],若使用適當?shù)臑V波方法濾除噪聲,有望提高分類精度;二是分類時只考慮后向散射信息,未考慮紋理特征等信息,若引入紋理特征等信息有望進一步提高分類精度。

4 結(jié)論與討論

地形起伏對SAR影像的影響直接關(guān)系到SAR數(shù)據(jù)的應(yīng)用效果,特別是在地形起伏較大地區(qū)進行定量遙感。本文基于模型模擬提出一種SAR影像地形校正半經(jīng)驗?zāi)P停饕Y(jié)論如下:

(1)簡單推導(dǎo)給出左視和右視兩種成像模式下本地入射角等的計算模型,然后基于微波后向散射模型AIEM,定量模擬分析地形起伏對不同雷達入射角度SAR的后向散射系數(shù)影響,分析可知:雷達入射角相對較小時,受到地形起伏因素的影響越小,因此地形起伏區(qū)域,雷達入射角相對較小的SAR數(shù)據(jù)是土壤濕度反演、森林蓄積量反演等定量遙感應(yīng)用的最佳數(shù)據(jù)源。

(2)建立后向散射系數(shù)與本地入射角的二次函數(shù)關(guān)系,并依此提出地形校正半經(jīng)驗?zāi)P停M數(shù)據(jù)和ASAR數(shù)據(jù)的校正結(jié)果表明本文模型在去除地形影響時的作用。對比校正前后ASAR影像監(jiān)督分類的結(jié)果,校正后ASAR數(shù)據(jù)的分類總體精度更高。

(3)植被覆蓋地表,植被是影響雷達后向散射的因素之一,下一步將研究地形起伏地表植被因素對后向散射的影響,以完善建立的地形校正模型;地表粗糙度也是影像影響雷達后向散射的重要因子,這是今后建立地形起伏地表SAR影像地形校正模型需考慮的又一因素。

致謝:感謝NASA提供的ASTER global digital elevation model VOO1數(shù)據(jù)。

[1] LIU L Y,WANG J H,BAO Y S,et al.Predicting Winter Wheat Condition,Grain Yield and Protein Content Using Multi-temporal Envisat-ASAR and Landsat TM Satellite Images[J].International Journal of Remote Sensing,2006,27(4):737-753.

[2] DENTE L,SATALINO G,MATTIA F,et al.Assimilation of Leaf Area Index Derived from ASAR and MERIS Data into CERES-Wheat Model to Map Wheat Yield[J].Remote Sensing of Environment,2008,112(4):1395-1407.

[3] BAGHDADI N,HOLAH N,ZRIBI M.Soil Moisture Estimation Using Multi-incidence and Multi-polarization ASAR Data[J].International Journal of Remote Sensing,2006,27(10):1907-1920.

[4] ZHANG Hailong,JIANG Jianjun,WU Hongan,et al.The BP Neural Net Work Classification Based on the Fusion of SAR and TM Images[J].Acta Geodaetica et Cartographica Sinica,2006,35(3):229-239.(張海龍,蔣建軍,吳宏安,等.SAR與TM影像融合及在BP神經(jīng)網(wǎng)絡(luò)分類中的應(yīng)用[J].測繪學報,2006,35(3):229-239.)

[5] FRANKLIN S E,LAVIGNE M B,HUNT E R,et al.Topographic Dependence of Synthetic Aperture Radar Imagery[J].Computer &Geosciences,1995,21(4):521-532.

[6] HINSE M,GWYN Q H J,BONN F.Radiometric Correction of C-band Imagery for Topographic Effects in Regions of Moderate Relief[J].IEEE Transactions on Geosicence and Remote Sensing,1988,26:122-132.

[7] RAUSTE Y.Incidence Angle Dependence in Forested and Non-forested Areas in SEASET SAR data[J].International Journal of Remote Sensing,1990,11(7):1267-1276.

[8] TEILLET P M,GUINDON B,MEUNIER J F,et al.Slope-aspect Effects in Synthetic Aperture Radar Imagery[J].Canadian J.Remote Sensing,1985,11(1):39-50.

[9] LECLERC G,BEAULIEU N,BONN F.A Simple Method to Account for Topography in the Radiometric Correction of Radar Imagery[J].International Journal of Remote Sensing,2001,22(17):3553-3570.

[10] REES W G,STEEL A M.Simplified Radar Mapping Equations for Terrain Correction of Space-borne SAR Images[J].International Journal of Remote Sensing,2001,22(18):3643-3649.

[11] THOMAS B,RUDOLF W,GUNTER S.Terrain Influences in SAR Backscatter and Attempts to Their Correction[J].IEEE Transactions on Geoscience and Remote Sensing,1991,29(3):451-462.

[12] VAN ZYL J J,CHAPMAN B D,DUBOIS P,et al.The Effect of Topography on SAR Calibration[J].IEEE Transactions on Geoscience and Remote Sensing,1993,31(5):1036-1043.

[13] LOEW A,MAUSER W.Generation of Geometrically and Radiometrically Corrected ScanSAR Images[J].Proceedings IEEE International Geoscience RS Symposium,2003,6(10):3995-3997.

[14] LOEW A,MAUSER W.Generation of Geometrically and Radiometrically Terrain Corrected SAR Image Products[J].Remote Sensing of Environment,2007,106(3):337-349.

[15] ZHANG Yonghong,ZHANG Jixian,LIN Zongjian,et al.Terra in Induced SAR Radiometric Distortions and Their Correction[J].Science of Surveying and Mapping,2002,27(4):23-26.(張永紅,張繼賢,林宗堅,等.地形引起的雷達輻射畸變及其校正[J].測繪科學,2002,27(4):23-26.)

[16] WU T D,CHEN K S.A Reappraisal of the Validity of the IEM Model for Backscattering from Rough Surfaces[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(4):743-753.

[17] CHEN K S,WU T D,TSANG L,et al.The Emission of Rough Surfaces Calculated by the Integral Equation Method with a Comparison to a Three-dimensional Moment Method Simulations[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(1):90-101.

[18] ZRIBI M,DECHAMBRE M.A New Empirical Model to Retrieve Soil Moisture and Roughness From C-band Radar Data[J].Remote Sensing of Environment,2002,84(1):42-52.

[19] REN Xin.A Surface Moisture Inversion Technique Using Multi-polarization and Multi-angle Radar Images[D].Beijing:Laboratory of Remote Sensing Information Sciences,Institute of Remote Sensing Applications,Chinese Academy of Sciences,2003:58-63.(任鑫.多極化、多角度SAR土壤水分反演研究[D].北京:中國科學院遙感應(yīng)用研究所,2003:58-63.)

[20] KUANG Yan,LI An,LI Ziyang,et al.Envisat Satellite Overview[J].Remote Sensing Information,2007,83(01):90-92.(匡艷,李安,李子楊,等.Envisat衛(wèi)星綜述[J].遙感信息,2007,83(01):90-92.)

[21] LI Z L,ZHU Q,GOLD C.Digital Terrain Modeling:Principles and Methodology[M].Florida:CRC Press,2005.

[22] GOYAL S K,SEYFRIED M S,O’NEILL P E.Correction of Surface Roughness and Topographic Effects on Airborne SAR in Mountainous Rangeland Areas[J].Remote Sensing of Environment,1999,67(2):124-136.

[23] HUANG Shiqi,LIU Daizhi.An Enhanced SAR Image Speckle Filter[J].Acta Geodaetica et Cartographica Sinica,2006,35(3):245-250.(黃世奇,劉代志.SAR圖像斑點噪聲抑制方法與應(yīng)用研究[J].測繪學報,2006,35(3):245-250.)

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包

主站蜘蛛池模板: 麻豆国产在线不卡一区二区| 国产91特黄特色A级毛片| 久久久久88色偷偷| 国产麻豆永久视频| 久久99国产视频| 白浆视频在线观看| 日本午夜在线视频| 亚洲天堂2014| 亚洲—日韩aV在线| AV不卡无码免费一区二区三区| 蝌蚪国产精品视频第一页| 在线播放真实国产乱子伦| 色播五月婷婷| 欧美日韩高清| 亚洲精品国产首次亮相| 欧美亚洲国产精品第一页| 欧美日韩国产在线观看一区二区三区| 另类专区亚洲| 91久久偷偷做嫩草影院| 114级毛片免费观看| 久久6免费视频| 国产永久在线视频| 久久99这里精品8国产| 亚洲成人在线网| 亚洲国产欧美国产综合久久 | 亚洲综合经典在线一区二区| 色悠久久综合| 毛片免费网址| 内射人妻无套中出无码| 色噜噜久久| 中文字幕无码制服中字| 精品久久久久久成人AV| 亚洲国内精品自在自线官| 国产免费一级精品视频| 亚洲第一区在线| 日韩a在线观看免费观看| 久久精品娱乐亚洲领先| 亚洲精品你懂的| 国产午夜福利亚洲第一| 欧洲高清无码在线| 激情影院内射美女| 四虎在线观看视频高清无码 | 青青操国产视频| 亚洲一区二区三区在线视频| 国产日韩欧美精品区性色| 成人亚洲国产| 午夜日韩久久影院| 欧美日韩综合网| 国产在线91在线电影| 国产午夜小视频| 国产精品99久久久| 国产欧美日韩91| 中文字幕在线看| 日韩一级毛一欧美一国产| 黄色三级网站免费| 亚洲精品777| 成人毛片免费观看| 欧美成人精品一级在线观看| 国产成人精品一区二区三区| 免费无码网站| 国产女同自拍视频| 乱系列中文字幕在线视频| 永久成人无码激情视频免费| 制服丝袜无码每日更新| 精品视频91| 成年免费在线观看| 少妇精品在线| 亚洲视频欧美不卡| 日韩专区欧美| 国产福利小视频在线播放观看| 久久性视频| 四虎影视永久在线精品| 久久熟女AV| 91无码人妻精品一区二区蜜桃| 欧美综合在线观看| 国产精品亚洲精品爽爽| 国内精品自在自线视频香蕉| 免费激情网站| av大片在线无码免费| 国产成人精品在线1区| 91成人在线观看| 亚洲清纯自偷自拍另类专区|