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

基于6S模型的MODIS影像逐像元大氣校正及其應用

2015-03-07 03:42:56姜琦剛
吉林大學學報(地球科學版) 2015年5期
關鍵詞:大氣模型

徐 言,姜琦剛

吉林大學地球探測科學與技術學院,長春 130026

?

基于6S模型的MODIS影像逐像元大氣校正及其應用

徐 言,姜琦剛

吉林大學地球探測科學與技術學院,長春 130026

以獲取地物真實反射率為目的,介紹了逐像元大氣校正的方法。應用6S模型逐像元綜合考慮太陽天頂角、傳感器天頂角、相對方位角、大氣氣溶膠厚度、觀測波段以及地表海拔這6個參數,生成查找表,統計分析各參數關于反射率的敏感度;并與以往單一參數校正的方法比較了校正的效果。結果表明,逐像元大氣校正算法更加接近地物的真實反射率。利用MODIS地表反射率產品對大氣校正的結果進行驗證,單一參數、逐像元大氣校正的相對誤差分別控制在26.9%和12.7%以內;在以植被指數(NDVI)為例的后續遙感定量化反演過程中,逐像元比單一參數大氣校正方法計算的NDVI平均高出14.4%。

逐像元;6S;查找表;大氣校正

0 引言

遙感技術的應用越來越廣泛,而隨著遙感應用的定量化、精細化發展,越來越多的定量化遙感涉及到光學分析[1-5]。因傳感器獲得的信號是地面信號和大氣信號的疊加,大氣影響不僅造成成像傳感器圖像模糊、對比度下降、信噪比降低以及細節損失[6],還使得遙感數據在后續的植被指數、地表溫度、葉綠素濃度等定量化反演中所獲得的結果嚴重偏離實際值[7];因此需要重視大氣對遙感數據的影響。

大氣校正是遙感輻射校正的重要步驟,是由遙感影像獲取反映地表真實反射率的一個必不可少的環節,對于定量遙感而言尤為重要[8-9]。這是由于傳感器在獲取信息過程中受到大氣分子、氣溶膠和云粒子等大氣成分吸收與散射的影響,使其獲取的遙感信息中帶有一定非目標地物的成像信息,數據預處理的精度達不到定量分析的高度。一方面,大氣的吸收、散射作用削弱了地表反射能量,對圖像信號有衰減作用;另一方面,大氣程輻射對圖像信號有增加作用[10]。消除這些大氣影響的處理,稱為大氣校正[11]。

以往進行大氣校正過程中,假設大氣條件均一,整幅影像僅使用一套參數進行校正。但在實際情況下,衛星觀測地面生成的影像每一個像元的大氣參數、觀測角度、海拔高度等參數均不相同,特別是在霧霾天氣的情況下,氣溶膠分布的不均勻更較為直觀地影響著衛星影像大氣校正的結果。

因此綜合上述情況,筆者嘗試逐像元考慮各方面參數進行大氣校正,選用MODIS(moderate-resolution imaging spectroradiometer)數據作為研究對象,基于6S輻射校正模型統計分析各個參數對反射率的敏感性,生成查找表,逐像元考慮大氣、角度、高程等參數進行大氣校正,以提高大氣校正的精度。

1 6S模型原理

6S模型由Eric Vermote等人在5S(the simulation of the satellite signal in the solar speetrum radiative code)模型的基礎上進行改進,采用SOS(successive order of scattering)的散射計算方法,建立在輻射傳輸理論基礎之上,準確模擬太陽--目標物--傳感器路徑上的大氣影響,應用范圍廣,受研究區特點及目標類型等的影響較小,精度高,與LOWTRAN、MODTRAN等模型同為目前發展比較成熟的大氣訂正模型[12]。6S模型適用于0.25~4.00 μm波長范圍內電磁波的大氣輻射傳輸的模擬。利用該模型進行大氣校正的工作流程是[13]:按順序將幾何參數、大氣模式、氣溶膠模式、可見度、海拔高度、觀測波段等參數輸入6S模型;通過計算模擬太陽輻射在大氣中的輻射狀況;給出大氣校正系數Xa、Xb和Xc[14]。

假設目標表面為均勻朗伯體表面,不考慮氣體吸收,在6S模型中衛星觀測的表觀反射率可用式(1)[15]表達:

(1)

式中:ρTOA為大氣表觀反射率;ρ0為大氣的路徑輻射項等效反射率;ρS為地表反射率;θs為太陽天頂角;θv為衛星天頂角;φs為太陽方位角;φv為衛星方位角;φs--φv代表相對方位角;S代表大氣球面反照率;T(θs)代表大氣下行透過率;T(θv) 代表上行輻射總透過率。

運行6S輻射傳輸模型后,可得到大氣校正系數Xa、Xb、Xc,然后按式(2)計算得到地表反射率[16]:

(2)

其中:

(3)

式中:Y為中間變量;L為定標后的輻射亮度值。

2 研究區概況

本文選用的MODIS遙感影像以安徽省為例,區域范圍為東經114.9°--119.6°,北緯29.4°--34.6°,成像日期為2013年7月3日。

根據研究區遙感影像圖,6S模型的初始輸入參數見表1,同時假定地表具有均一的朗伯反射特性。其中,氣溶膠厚度利用MOD04氣溶膠產品得出,地面高程利用SRTM(shuttle radar topography mission)地形產品數據得出,二者均經過投影、剪裁、重采樣、輻射定標等預處理。

3 敏感性分析

對于處理大范圍遙感影像圖,參數一致的假設顯然是導致地表反射率校正誤差的原因[17]。因此通過對敏感性的分析,我們能更深層次地了解各種參數對大氣校正結果的影響程度,為大氣校正參數選取提供參考。根據6S模型模擬,影響大氣校正獲得地表反射率的參數包括幾何參數、目標和傳感器的高程參數、大氣模式、氣溶膠模式以及光譜條件。因此,筆者選擇了太陽天頂角、傳感器天頂角、相對方位角、大氣氣溶膠厚度、觀測波段以及地表海拔這6個參數的變化來分析不同參數對大氣校正的影響。利用IDL語言線下調用6S輻射模型生成查找表,統計分析各參數關于地表反射率敏感度的優劣(圖1),各參數的取值范圍如下:太陽天頂角0°~70°,衛星天頂角0°~80°,氣溶膠光學厚度0.00~1.95,海拔高程0.0~1.9 km。

表1 6S模型的初始輸入參數

據圖1分析,氣溶膠光學厚度對反射率的影響穩定且明顯,氣溶膠厚度越大對反射率的敏感度也就越好(圖1a)。太陽天頂角在角度小于25°時,對地表反射率的影響不明顯,而太陽天頂角大于25°時,太陽天頂角對反射率的敏感度明顯增高,甚至高于氣溶膠光學厚度對于反射率的敏感度(圖1b)。衛星天頂角在MODIS的波段3和波段4的敏感度明顯高于其他波段,隨著角度的增大反射率略有下降,衛星天頂角整體敏感度不高(圖1c)。海拔高度對反射率的敏感度同樣較低,筆者僅計算了海拔高度在1.9 km以下的地表高程,海拔高度僅在MODIS的波段3和波段7對反射率有一定的影響,其他波段的敏感度幾乎可忽略不計(圖1d)。另外,筆者也注意到當太陽天頂角變化較大時與變化較小時觀測的海拔高程對反射率的敏感度幾乎一致。而相對方位角對反射率的敏感度最差,相對方位角對應的大氣校正系數僅Xb發生細微變化(圖1e);根據公式(2)、(3)可得出相對方位角對反射率的影響幾乎可忽略。

對地表反射率影響敏感度較高的有氣溶膠光學厚度及太陽天頂角,而相對方位角敏感度最差。經過針對上述各個參數的敏感性分析,我們在對某一景影像做大氣校正時,可以根據這景衛星圖像的各參數情況,選擇性地忽視一些敏感度較低的參數,這樣在提高大氣校正效率的同時也一定程度地縮小了大氣校正的工作量,同時也擴展了利用輻射傳輸模型進行大氣校正的應用。

4 查找表的建立和逐像元校正方法

筆者在對研究區做逐像元大氣校正時,考慮到因參數較多導致數據量偏大、耗時較多,因此針對研究區利用IDL語言調用6S輻射模型做查找表,利用6S模型離線計算不同太陽天頂角、衛星天頂角、氣溶膠光學厚度以及地表海拔高程情況下的校正系數Xa、Xb、Xc,生成研究區上述各個參數范圍內的大氣校正系數。由于大氣校正是逐像元分析而且又加入考慮了多個影響大氣反射率的參數,所以針對工作區做查找表的好處在于,輸入大氣校正參數的范圍被縮小,根據各參數范圍分別調整步長,這樣在縮小工作量的同時可以更加精細化地計算大氣校正的參數。

圖1 不同觀測波段上的氣溶膠光學厚度(a)、太陽天頂角(b)、衛星天頂角(c)、海拔高程(d)及相對方位角(e)對反射率的敏感度Fig.1 AOD(a), solar zenith(b), sensor zenith(c), ground altitude(d) and relative azimuth(e) of the sensitivity to the reflectivity on different observation bands

本文研究區面積較小,研究區內參數的范圍比較狹小,所以本文參數的取值范圍如下:太陽天頂角為0°~17°時步長設為3°,太陽天頂角為24°~60°時步長設為12°,共11個值;衛星天頂角0°~75°以15°為步長,共6個值;氣溶膠光學厚度0~1.95,步長不規則,共11個值;地表海拔0~1.9 km,步長不規則,共9個值;以及MODIS的前7個波段。將上述參數組合代入6S輻射模型,得到不同參數下大氣校正系數的查找表。

逐像元大氣校正步驟:

1)根據遙感圖像的成像時間,獲取研究區那一時刻的大氣參數、角度信息以及高程數據,分別進行投影、剪裁、重采樣等預處理,使得獲取的各類參數對應于遙感影像信息。

2)利用獲取的研究區內的各類參數信息,分別取均值作為輸入6S輻射模型的初始參數,得到6S輻射校正輸出表。

3)在輻射校正輸出表中確定需要參數的位置,根據研究區各參數范圍分別設計步長,并利用IDL編程循環調用6S輻射校正模型,以獲取不同參數不同組合的大氣校正系數查找表[14]。

4)對大氣校正系數查找表進行線性插值,并根據研究區各像元的大氣、角度、高程等信息查詢查找表中對應的大氣校正系數,根據像元對應的大氣校正系數計算該像元處的地表反射率。

5 結果與驗證

筆者將查找表逐像元大氣校正所得出的結果與統一參數大氣校正以及校正前的表觀反射率結果進行了對比,采用人工目視解譯的方法隨機選取了陸地和水體做為典型地類分別獲取像元處的地表反射率,結果如圖2所示。由于氣溶膠及觀測角度等參數的不同,統一參數大氣校正與逐像元大氣校正的精度差值比較明顯。

a.陸地;b.水體;c.頻率。圖2 三種不同方法的校正結果在各觀測波段上的比較Fig.2 Comparison of the correction result of three method for each observation band

校正前衛星定標后得到的表觀反射率雖然消除了太陽天頂角、日地距離和太陽輻射量的差異,但是不能消除由大氣程輻射、瑞利散射和吸收帶來的影響[18]。這些影響造成的大氣衰減使得衛星傳感器接收到的信號很大部分不是地面反射率的光譜信息,校正前后差值很大。地表反射率圖像中地物在可見光波段(波段1--4)反射率減小,其中藍光波段(波段3)的瑞利散射最強,反射率減小最為明顯,綠光(波段4)次之(圖2a、b)。

圖3 掃描線上各像元的氣溶膠光學厚度(a)、太陽天頂角(b)、衛星天頂角(c)、海拔高程(d)、相對方位角(e)情況Fig.3 AOD(a), solar zenith(b), sensor zenith(c), ground altitude(d) and relative azimuth(e) of the pixels on the scan line

根據研究區校正前后反射率頻率變化圖(圖2c)可以看出,校正后的反射率峰值向低反射率方向移動,整體上校正前的表觀反射率值大于校正后的反射率值。

圖4 第497行反射率(a)、植被指數(b)在三種大氣校正方法下的結果對比Fig.4 The 497 line reflectivity and vegetation index results were compared in three atmospheric correction methods

由于6S模型計算耗時太長,取研究區像元第497行來進行驗證。圖3給出了研究區第497行掃描線上自西向東各個像元的氣溶膠光學厚度、太陽天頂角、衛星天頂角、海拔高程和相對方位角的情況,結合各參數對反射率的敏感度,對比校正前后反射率的變化情況。本文以紅光波段為例,根據掃描線上的像元參數進行逐像元大氣校正,將校正結果與統一參數大氣校正以及表觀反射率的結果進行對比(圖4a)。大氣校正削弱了大氣瑞利散射和氣溶膠散射所導致的紅光波段信號增強,而地形因素又使得校正后的反射率大于表觀反射率,校正后的反射率與實際情況更為接近。利用MODIS地表反射率產品對大氣校正的結果進行驗證,統一參數、逐像元大氣校正的相對誤差分別控制在26.9%和12.7%以內。因此,對于影像范圍較大,或氣溶膠和地面高程分布不均勻的山區或者農田等地逐像元大氣校正方法更具優勢。

以植被指數為例證明大氣校正精度對后續遙感定量計算的影響程度[19](圖4b)。從圖4b可以看出,校正前后同一剖面(第497行)兩條NDVI曲線的分布和變化趨勢吻合,峰谷和峰值基本一致,逐像元比單一參數大氣校正方法計算的NDVI平均高出14.4%。究其原因,經過大氣校正后,可見光的反射率降低,近紅外波段的反射率升高,兩波段的對比度增高,光在大氣傳播過程中所受到的衰減得到彌補,則影像的植被指數呈現增加的趨勢[18]。

6 結論

分析了6S模型中各種參數對反射率的敏感度,并詳細敘述了逐像元大氣校正方法。結論如下:

1)利用查找表逐像元大氣糾正的算法,根據各個像元的大氣狀況、觀測角度、地面高程等條件差異,綜合考慮了各種敏感度不同的參數,與以往常用的基于整景遙感影像獲取平均參數并進行大氣校正的算法相比,逐像元大氣校正算法與實際情況較為接近,能夠更加精確地對遙感影像進行大氣校正并獲取地物的真實反射率。

2)利用MODIS地表反射率產品對大氣校正的結果進行驗證,單一參數、逐像元大氣校正的相對誤差分別控制在26.9%和12.7%以內;在以植被指數(NDVI)為例的后續遙感定量化反演過程中,逐像元比單一參數大氣校正方法計算的NDVI平均高出14.4%。

3)查找表的建立是大氣校正精度提高的關鍵,在一定范圍內,參數的步長過長影響校正的精度,而步長過小則使得數據量過大,冗余數據較多。因此可以根據研究區情況設定參數的范圍,并分別設計各參數的步長,根據研究區來建立查找表,而不是固定用一套查找表。

4)經過大氣校正后的植被指數整體上明顯提高,本文的逐像元大氣校正方法可以用于分析影像范圍較大,或氣溶膠和地面高程分布不均勻的山區或者農田等地。

[1] Cairns B, Carlson B E, Ying R X, et al. Atmospheric Correction and Its Application to an Analysis of Hyperion Data[J]. IEEE Transactions on Geoscience & Remote Sense, 2003, 41(6):1232-1245.

[2] 杜鑫,陳雪洋,蒙繼華,等. 基于6S模型的環境星CCD數據大氣校正[J]. 國土資源遙感, 2010(2):22-25. Du Xin, Chen Xueyang, Meng Jihua, et al. Atmospheric Correction of HJ-1 CCD Data Based on 6S Model[J]. Remote Sensing for Land and Resources, 2010(2):22-25.

[3] Liang S L, Fang H L, Chen M Z, et al. Validating MODIS Land Surface Reflectance and Albedo Products: Methods and Preliminary Results[J]. Remote Sensing of Environment, 2002, 83(1): 149-162.

[4] 程偉,王黎明,田慶久.一種基于陰影像元的光學遙感大氣校正方法[J]. 測繪學報,2008,37(4):469-475. Cheng Wei, Wang Liming, Tian Qingjiu. A Method of Atmospheric Correction Based on Shadow-Pixel for Optical Satellite Data[J]. Acta Godaetica et Cartographica Sinica, 2008, 37(4):469-475.

[5] 申茜,張兵,李俊生,等.航天高光譜遙感器CHRIS的水體圖像大氣校正[J]. 測繪學報,2008,37(4):476-481. Shen Qian, Zhang Bing, Li Junsheng, et al. Atmospheric Correction for Waterbody Images Acquired by Spaceborne Hyperspectral Sensor CHRIS[J]. Acta Geodatetica et Cartographica Sinica, 2008, 37(4): 476-481.

[6] 王釗.6S輻射模型算法解析及在MODIS大氣校正中的應用[J]. 陜西氣象,2006(5):34-37. Wang Zhao. 6S Radiation Model Parsing Algorithm and Its Application in the MODIS Atmospheric Correction[J]. Journal of Shaanxi Meteorology, 2006(5):34-37.

[7] 李衛國,蔣楠,王紀華.基于薄云霧去除的ETM+影像大氣校正[J]. 農業工程學報,2013, 29(增刊1):82-88. Li Weiguo, Jiang Nan, Wang Jihua. Atmospheric Correction for ETM+ Image Based on Thin Cloud Removal[J]. Transactions of the Chinese Society of Agricultural Engineering, 2013, 29(Sup.1): 82-88.

[8] Liang S L. Advanced Remote Sensing[M]. Oxford: Academic Press, 2012: 111-126.

[9] 王明常,王亞楠,陳圣波,等. 多角度高光譜CHRIS/Proba植被模式數據大氣校正[J]. 吉林大學學報:地球科學版,2011,41(2):609-614. Wang Mingchang, Wang Yanan ,Chen Shengbo, et al. Atmospheric Correction Research on Vegetation Patterns of Multi-Angle Hyperspectral CHRIS/Proba Data[J]. Journal of Jilin University: Earth Science Edition, 2011,41(2):609-614.

[10] 劉其悅,余濤,方莉,等.基于6S模型的HJ-1/CCD影像逐像元大氣校正[J]. 黑龍江科技信息,2010(31):11. Liu Qiyue, Yu Tao, Fang Li, et al. Based on the 6S Model HJ-1/CCD Pixel by Pixel Atmospheric Correction[J]. Heilongjiang Science and Technology Information, 2010(31):11.

[11] 亓雪勇,田慶久.光學遙感大氣校正研究進展[J]. 國土資源遙感,2005(4):1-6. Qi Xueyong, Tian Qingjiu. The Advances in the Study of Atmospheric Correction for Optical Remote Sensing[J]. Remote Sensing for Land & Resources, 2005(4):1-6.

[12] 阿布都瓦斯提·吾拉木,秦其明,朱黎江.基于6S模型的可見光、近紅外遙感數據的大氣校正[J]. 北京大學學報:自然科學版,2004,40(4):611-618. Ghulam A, Qin Qiming, Zhu Lijiang, et al. 6S Model Based Atmospheric Correction of Visible and Near-Infrared Data and Sensitivity Analysis[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2004, 40(4):611-618.

[13] 徐永明,覃志豪,陳愛軍.基于查找表的MODIS逐像元大氣校正方法研究[J]. 武漢大學學報:信息科學版,2010,35(8):959-962. Xu Yongming, Qin Zhihao, Chen Aijun, et al. A Pixel-by-Pixel Atmospheric Correction Algorithm for MODIS Data Based on Look-up Table[J]. Geomatics and Information Science of Wuhan University, 2010,35(8):959-962.

[14] 甘文霞,沈煥鋒,張良培,等.采用6S模型的多時相MODIS植被指數NDVI歸一化方法[J]. 武漢大學學報:信息科學版,2014,39(3):300-304. Gan Wenxia, Shen Huanfeng, Zhang Liangpei, et al. Normalization of Multi-Temporal MODIS NDVI Based on 6S Radiative Transfer Model[J]. Geomatics and Information Science of Wuhan University, 2014,39(3):300-304.

[15] Vermote E F, Tanre D, Deuze J L, et al.The Se-cond Simulation of the Satellite Signal in the Solar Spectrum (6S) User’s Guide: Version 3[Z]. [S. l.]: NASA Godard Space Flight Center, 2006:5-16.

[16] 孫源,顧行發,余濤,等.環境星CCD數據大氣校正研究[J]. 國土資源遙感,2010(4):6-9. Sun Yuan, Gu Xingfa, Yu Tao, et al. The Environmental Satellite Data Atmospheric Correction Research[J]. Remote Sensing for Land & Resources, 2010(4):6-9.

[17] 鄭文武,曾永年.利用MODIS數據和查找表的城區TM影像大氣校正[J]. 測繪科學,2012,37(3):23-25. Zheng Wenwu, Zeng Yongnian. An Atmospheric Correction Algorithm for Urban Area TM Data Base on MODIS Data and Look-up Table[J]. Science of Surveying and Mapping, 2012, 37(3):23-25.

[18] 孫長奎,孫林,麻盛芳,等. HJ-1 CCD數據大氣校正方法研究[J]. 遙感學報,2012,16(4):826-836. Sun Changkui, Sun Lin, Ma Shengfang, et al. Atmospheric Correction Method Based on HJ-1 CCD Data[J]. Journal of Remote Sensing, 2012, 16(4):826-836.

[19] Franch B, Vermote E F, Sobrino J A, et al. Analysis of Directional Effects on Atmospheric Correction[J]. Remote Sensing of Environment, 2013, 128(1):276-288.

A Pixel by Pixel Atmospheric Correction Algorithm and Its Application for MODIS Data Based on 6S Model

Xu Yan,Jiang Qigang

CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China

In order to obtain the real reflectance of earth objects, we introduced a new method to correct atmospheric parameters pixel by pixel. Considering with the six parameters, including solar zenith angle, sensor zenith angle, azimuth angle, aerosol optical depth, observation band, and surface elevation, we applied the 6S model pixel by pixel to generate a lookup table for analyzing the sensitivity of each parameter to reflectance; and compared the correction results with the ones by the old method using single parameter. The results show that the reflectance by 6S model pixel by pixel is much closer to the true reflectance of objects. By using the MODIS surface reflectance products, we verified the accuracy both of the 6S model pixel by pixel and the single parameter method, and found that their relative errors are 12.7% and 26.9% respectively. Taking vegetation index (NDVI) as an example for the subsequent quantitative remote sensing inversion, theNDVIby 6S pixel by pixel method is 14.4% higher than the single parameter atmospheric correction method.

pixel by pixel; 6S; lookup table; atmospheric correction

10.13278/j.cnki.jjuese.201505304.

2015-01-17

中國地質調查局項目(1212010510613)

徐言(1987--),女,博士研究生,主要從事GIS與遙感研究,E-mail:xuyanby@163.com

姜琦剛(1964--),男,教授,博士生導師,主要從事GIS與遙感地學環境的研究,E-mail:jiangqigang_rs@163.com。

10.13278/j.cnki.jjuese.201505304

P407.8;TP751

A

徐言,姜琦剛. 基于6S模型的MODIS影像逐像元大氣校正及其應用.吉林大學學報:地球科學版,2015,45(5):1547-1553.

Xu Yan,Jiang Qigang. A Pixel by Pixel Atmospheric Correction Algorithm and Its Application for MODIS Data Based on 6S Model.Journal of Jilin University:Earth Science Edition,2015,45(5):1547-1553.doi:10.13278/j.cnki.jjuese.201505304.

猜你喜歡
大氣模型
一半模型
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
太赫茲大氣臨邊探測儀遙感中高層大氣風仿真
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
莊嚴大氣 雄媚兼備
主站蜘蛛池模板: 久久不卡国产精品无码| 九色91在线视频| 久久综合伊人77777| 欧美日韩国产在线人| 色欲不卡无码一区二区| 日本在线视频免费| 久久精品无码国产一区二区三区 | 亚洲精品男人天堂| 综合久久五月天| 精品国产黑色丝袜高跟鞋| 久久99国产精品成人欧美| 91亚洲视频下载| 三上悠亚在线精品二区| 嫩草在线视频| 国产微拍一区| 亚洲性视频网站| 国产一级视频久久| 国产亚洲精久久久久久无码AV| 久久人人97超碰人人澡爱香蕉| 亚洲一区二区视频在线观看| 波多野结衣第一页| 国产精品久久精品| 精品国产一二三区| 女人一级毛片| 97精品久久久大香线焦| 黄色网在线| 秋霞午夜国产精品成人片| 久久免费精品琪琪| 成人一区在线| 97se亚洲综合在线天天| 69精品在线观看| 日韩黄色大片免费看| 91国内在线观看| 91精品啪在线观看国产| 免费看美女自慰的网站| 亚洲va在线∨a天堂va欧美va| 亚洲色图综合在线| 国产成人精品无码一区二| 国产网站免费观看| 欧美精品黑人粗大| 18禁影院亚洲专区| 一级成人a做片免费| 99尹人香蕉国产免费天天拍| 国产精品九九视频| 亚洲第一区欧美国产综合| 亚洲精品无码久久久久苍井空| 国产精品亚欧美一区二区| 欧美亚洲国产精品第一页| 51国产偷自视频区视频手机观看| 97视频免费在线观看| 天天色综合4| 男人天堂伊人网| AⅤ色综合久久天堂AV色综合 | 国产在线视频导航| 亚洲高清日韩heyzo| 亚洲欧洲日产国产无码AV| 99视频精品全国免费品| 成年人免费国产视频| 国产成人超碰无码| 国产成人精品优优av| 波多野结衣国产精品| 波多野结衣一区二区三区88| 蜜臀av性久久久久蜜臀aⅴ麻豆| 国产福利微拍精品一区二区| 亚洲天堂精品视频| 久久免费观看视频| 日韩高清一区 | 色天堂无毒不卡| 思思热精品在线8| 91九色视频网| 国产无码网站在线观看| 欧美激情第一欧美在线| 伊伊人成亚洲综合人网7777| 国产精品无码久久久久久| 免费无码AV片在线观看国产| 国产极品美女在线| 女人天堂av免费| 国产精品片在线观看手机版| 国产毛片不卡| 九色在线视频导航91| 国产农村妇女精品一二区| 蜜臀AV在线播放|