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

基于Landsat8遙感影像和SEBS模型的呼圖壁縣蒸散量時空格局分析

2016-06-05 14:15:31張圓鄭江華劉志輝姚俊強
生態科學 2016年2期
關鍵詞:區域模型

張圓, 鄭江華,*, 劉志輝,3, 姚俊強

1. 新疆大學資源與環境科學學院, 新疆 烏魯木齊 830046 2. 綠洲生態教育部重點實驗室, 新疆 烏魯木齊 830046 3. 新疆大學干旱與半干旱生態研究所, 新疆 烏魯木齊 830046

基于Landsat8遙感影像和SEBS模型的呼圖壁縣蒸散量時空格局分析

張圓1,2, 鄭江華1,2,*, 劉志輝1,2,3, 姚俊強1,2

1. 新疆大學資源與環境科學學院, 新疆 烏魯木齊 830046 2. 綠洲生態教育部重點實驗室, 新疆 烏魯木齊 830046 3. 新疆大學干旱與半干旱生態研究所, 新疆 烏魯木齊 830046

張圓, 鄭江華, 劉志輝, 等. 基于Landsat8遙感影像和SEBS模型的呼圖壁縣蒸散量時空格局分析[J]. 生態科學, 2016, 35(2): 26-32.

ZHANG Yuan, ZHENG Jianghua, LIU Zhihui, et al. Spatial and temporal distribution of evapotranspiration in the hutubi County based on Landsat8 data and SEBS model[J]. Ecological Science, 2016, 35(2): 26-32.

利用Landsat8影像, 采用SEBS模型, 結合呼圖壁縣氣象站觀測的溫度風速、日照時數等氣象數據, 對新疆昌吉回族自治州呼圖壁縣2013年4月22日、6月9日、8月28日 10月15日的蒸散發量進行了估算。從時間上看, 估算結果存在明顯的季節變化規律: 夏季最大, 春季次之, 秋季最小, 以耕地為例四天蒸散發量分別為: 1.938, 3.136, 2.641和1.314 mm·day–1。從空間上看, 縣域蒸散發量整體變化趨勢為: 北部荒漠區<中部平原區<南部山區。四天當中最大值出現在南部山區6月9日達到了4.128mm·day–1。對SEBS的估算結果與呼圖壁縣氣象站的觀測結果和利用彭曼公式計算的結果進行比較,表明SBES 模型的結果是合理的, 可以在實踐中用來反映天山北坡典型縣域蒸散量的時空變化特征。

SEBS; 蒸散發; Landsat8; 天山北坡; 縣域

1 前言

蒸散發是地表水量和熱量平衡的重要參量, 包含土壤蒸發、水面蒸發和植被的蒸騰, 也是衡量植被生長狀況和作物產量的重要指標[1]。對于區域尺度上蒸發的估算,遙感技術是最為經濟和最為準確的手段[2], 不同學者從不同角度對蒸散發的研究進展進行了系統總結[3–4], 經過20多年的發展, 出現了許多的遙感反演蒸散發的模型, 荷蘭Wageningen大學為主發展起來的地表能量平衡算法SEBAL (surface energy balance algorithm for land)和地表能量平衡系統SEBS[5](surface energy balance system)是兩種目前使用最為廣泛的基于地表能量平衡原理的單層模型。其中, SEBS模型是由荷蘭籍華人蘇中波提出的[3], 它是在SEBAL基礎上發展起來的遙感蒸散單層模型,相對SEBAL而言, SEBS具備更明確的物理概念和較高的通量估算精度,具有更好的實用價值,因而近年來在國內外獲得了較廣泛的應用[6],例如李發鵬等[7]基于MODIS 數據, 應用表面能量平衡(SEBS)模型, 對黃河三角洲的區域陸面蒸散發量進行了估算, 分析區域蒸散發量的時空分布特征;張雨航等[8]以海流兔河流域為例, 利用MODIS數據,結合氣象資料, 運用遙感模型中的表面能量平衡法對該流域蒸散量進行估算; 金曉媚[9]等基于MODIS遙感數據, 應用表面能量平衡系統(SEBS), 對柴達木盆地及8個水資源區2001—2011年的區域蒸散量進行了計算, 并分析其影響因素。田國珍[10]等利用風云三號、風云二號氣象衛星數據, 結合自動站氣象數據, 基于SEBS模型對山西省進行干旱監測研究;唐婷[11]等以京津唐地區為例, 基于SEBS模型, 利用MODIS遙感數據和氣象數據, 計算了2000、2005和2010年四季代表月份的平均日蒸散發量, 并結合3期土地利用圖, 定量評估了由城市擴張引起的日蒸散發量的變化。Ma W, Hafeez M等[12], 使用ASTER衛星數據, 利用SEBS模型估算澳大利亞新南威爾士西南的Coleamball灌區的蒸散量, 結合實地測量數據證明了SEBS模型在該研究估算的合理性。這些研究為區域水資源管理科學可靠的耗水和需水量估算提供了較為可行的理論和方法基礎, 但他們的研究區域均為較大空間范圍區域, 應用的多為中低空間分辨率的遙感影像, 對于較小區域、較小流域的遙感蒸散量估算卻少有涉及, 而在最嚴格水資源管理背景下, 較小區域(縣級行政區劃)或較小流域的更精確蒸散量估算成為水資源合理有效分配和減少資源紛爭的必然要求。

呼圖壁縣位于新疆天山北坡, 屬于干旱半干旱氣候區, 呼圖壁縣擁有軍塘湖河和呼圖壁河兩個小流域, 具有從冰川~高山草甸~林帶~高山草原~平原灌區~荒漠區的地貌特征。該縣是論文依托的最嚴格水資源管理示范項目的研究區。研究采用2014年3月18日服務用戶的30米空間分辨率的Landsat8影像, 結合傳統的SEBS模型, 以期計算出更為合理精確的每年四季典型時期縣域蒸散量, 服務研究區的最嚴格水資源管理示范。

對四個時期的蒸散量結果進行驗證, 表明該模型在估算呼圖壁縣蒸散發上具有一定的精度, 可滿足區域日蒸散發估算的需要, 能為生態需水的計算提供一定基礎。

2 研究方案與研究方法

2.1 研究方案

本文利用Landsat8衛星影像反演一系列地表參數, 包括寬波段反射率、比輻射率、NDVI、利用單窗算法反演出的地表溫度, 結合地面觀測的氣象資料(中國氣象科學數據共享服務網http://cdc.cma. gov.cn/home.do), 主要有氣溫, 日照時間, 風速, 氣壓等, 以及ASTER DEM 高程, 將以上參數在ENVI當中轉換成TIFF格式, 導入ILWIS軟件中計算得到呼圖壁縣蒸散發的分布結果, 技術路線圖詳見圖 1。

圖1 蒸散量計算流程簡圖Fig. 1 Computational flow diagram of evapotranspiration

2.2 SEBS模型原理及主要參數

SEBS模型應用對遙感數據處理獲得的一系列地表參數(如反照率、比輻射率、地表溫度、NDVI 等), 結合地面觀測的氣象數據(如氣溫, 降雨, 相對濕度, 風速, 氣壓等), 對區域尺度空氣動力學參數(d0,z0m)、地表凈輻射、感熱進行估算, 進而根據實際蒸發比(潛熱/(凈輻射-地熱))得到區域尺度的蒸散發[13]。

任意時刻的地表能量平衡方程為

式中:Rn為凈輻射(w·rn–2);G0為土壤熱通量(w·m–2);λE為潛熱通量(w·rn–2) (其中λ=2.49×106為水的汽化熱(J·kg–1);E為蒸散率(kg·m–2·s–1);H為感熱通量(w·m–2)。

(1) 凈輻射(Rn)公式如下:

式中:a為地表反照率; ε為地表發射率;T0為地表溫度(℃); δ為斯蒂芬—波爾茲曼常數(5.67e10-8)。

(2) 土壤熱通量(G0)

式中:Γc為植被覆蓋較好區域的參數為0.05;Γs為裸土區域的參數為0.315;fc為植被覆蓋度。

(3) 感熱通量(H)。

在大氣近地層中, 根據大氣邊界層相似理論,有以下的關系

式中:z為參考高度(m);u為風速(m·s–1);u﹡為摩擦風速(m·s–1);d0為零平面位移高度(m);z0m和z0h分別為動力學粗糙長度(m)和地表熱傳輸粗糙長度(m);Ψm和Ψk分別為動力學和熱力學傳輸的穩定度訂正函數;θ0和θα分別為觀測面和參考面高度的虛溫(℃);L為莫寧霍夫長度(m);H為感熱通量(w·m–2);k為卡爾曼常數;ρ為空氣密度(g·m–3);Cp為空氣的熱容(J·g–1·℃);θv為近地表的位溫(℃);g為重力加速度(m·s–2)。

摩擦風速(u)、感熱通量(H)、莫寧霍夫長度(L)可以通過迭代求解方程(4)、(5)、(6)得到。其他的變量可通過氣象觀測信息結合遙感觀測信息求得。

(4) 相對蒸發比的提出。

根據地表能量平衡方程, 在土壤水分虧缺的干燥地表環境下, 由于沒有土壤水分供給蒸發, 潛熱通量約為零, 此時感熱通量達到最大值:

式中:Hdry為干燥地表環境下的感熱通量(w·m–2)。

土壤充分供水情況下的感熱通量是通過結合Penman-Monteith由公式得出:

2.3 Landsat8影像特點

Landsat 數據的波譜信息豐富、空間分辨率高,數據源穩定, 是通過遙感技術反演蒸散發的理想數據源[14]。2013年2月11號, NASA 成功發射了Landsat 8衛星, 為四十年的Landsat計劃注入了新鮮的血液。Landsat8上攜帶有兩個主要載荷: OLI和TIRS。Landsat8衛星在設計的特征上與之前Landsat系列基本相同, 因此Landsat8影像與前期的Landsat數據保持很高的一致性和可比性。基于Landsat8自身的特點, 研究者們已經開始挖掘Landsat8影像對于生態環境和地表監測的不同意義[15–17]。

3 研究區域與數據

3.1 研究區概況

呼圖壁縣屬于新疆昌吉回族自治州, 境內主要有呼圖壁河、軍塘湖河兩大水系。地勢南高北低, 由南向北的地形分布依次為:南部多為高山和丘陵, 平均海拔在2400米左右, 占全縣總面積的31.6%; 中部是整個呼圖壁縣主要的農作物種植區, 是多年形成的沖積平原, 平均海拔在580米之間, 占總面積43.2%; 北部為沙漠和戈壁, 海拔在360-460米之間,占總面積的25.2%。呼圖壁縣受到中緯度西風帶的控制下, 是溫帶大陸性氣候, 因此晝夜溫差大, 冬長夏短。氣候隨著南北部的地形差異而變化明顯, 其中低山、平原和沙漠地區屬中溫帶, 南部中山和高山地區屬于寒溫帶。溫度由北向南逐漸降低, 年最高氣溫為 36.0 ℃-43.1 ℃, 日照時數2900小時, 年降水量為110-400毫米, 蒸發量2300毫米, 年均風速3.1 m·s-1[18]。研究區如圖所示:

3.2 遙感數據與氣象數據

下載http://glovis.usgs.gov/ USGS網站中4月22日、6月9日、8月28日, 10月15日Landsat8遙感影像, 并且進行拼接、裁剪、大氣校正等預處理。

根據ILWIS軟件中關于SEBS模塊對參數的需要,結合實際獲取的數據, 所輸入的數據分為遙感和氣象兩類, 遙感數據為: 由Landsat8數據反演得到的地表溫度(K), 比輻射率, 反射率, NDVI和ASTERM DEM高程, 以及自中國氣象科學共享網的氣象數據,包括衛星過境當天的溫度, 氣壓, 風速, 日照時數等。SEBS模型所輸入數據詳見表 1。

圖2 研究區示意圖Fig. 2 Study area

4 結果分析

利用SEBS模型估算呼圖壁縣4月22日、 6 月13日、8月28日 、10月15日4天的蒸散量, 對應植被生長隨季節變化的規律, 并對蒸散發結果與氣象站的蒸發皿測量值和彭曼公式計算的參考蒸散量進行對比。

4.1 呼圖壁縣蒸散發量的時空規律分析

模型計算的呼圖壁縣單日蒸散量分布圖, 如圖2所示, 圖3(a)、圖3(b)、圖3(c), 圖3(d)分別代表4月22日、6月13日、8月28日, 10月15日呼圖壁縣蒸散量,圖中可以看出估算結果存在明顯的季節變化規律:夏季最大, 春季次之, 秋季最小。空間上整體變化趨勢為: 北部荒漠區<中部平原區<南部山區。

該結果反映出的呼圖壁河流域蒸散發空間變化特征是: 1)南部高原及丘陵地區多被林地覆蓋,蒸散發量較大2)中部平原地區蒸散發量能夠和農作物的生長具有較好對應, 受到人類活動影響最為明顯 3)呼圖壁縣南部為戈壁和沙漠覆蓋, 供給蒸發的水分不足, 蒸散量較小。

從時間上來看: 4月22日蒸散結果圖3(a)對應著初春, 是農作物生長的初期, 大部分地區還未開始播種, 此時氣溫較低限制蒸發能力, 中部的平原區以耕地為主其中夾雜部分城鎮, 南部高山和丘陵大多被森林覆蓋是呼圖壁河的產流區, 因此蒸散量相對較大。

表1 輸入信息Tab. 1 The input parameters for SEBS

圖3 SEBS模型計算結果Fig. 3 Evapotranspiration results

6月9日蒸散結果對應初夏, 由于大量引水灌溉,水分供給充足, 作物生長比較旺盛, 較 4月份相比南部山區與中部平原區蒸散量顯著升, 大水面如水庫等地的蒸散量值也大幅增加, 呼圖壁縣蒸散量的差異由北向南也達到最大, 位于南部山區的呼圖壁河產流區蒸散量極高, 北部荒漠區蒸此時散量仍為三個部分的最低。

8月28日的蒸散結果對應于夏末, 整個呼圖壁縣蒸散量仍處于較高水平, 中部平原區與南部山區蒸散量值略低于6月。

10月15日的蒸散結果對應于秋季, 隨著溫度下降, 植被生長季進入尾聲, 部分農作物已經收割完成, 蒸散量較8月28日有較大回落, 中部平原區蒸散量較8月28日大大減少。

為近一步分析不同地物類型下的蒸散發量特征,利用GEOEYE-1 影像將地表分為耕地、林地、水體、裸地、沙地、其他共計6類, 使用ArcGIS軟件中的區域分析功能, 得到各類地物蒸散發量統計直方圖如下圖4所示:

我們可以得出以下結論: 林地和水體的蒸散量在所計算的 4天當中的值都較高, 主要原因是北部的林地屬于高海拔山區, 受到的太陽輻射多, 植被覆蓋度較高, 蒸散量較大。中部平原區受到人類活動影響, 蒸散量的季節變化很明顯。

4.2 模型驗證

4.2.1 蒸發皿數據驗證

SEBS模型蒸散發結果的精度驗證數據為中國氣象科學數據共享服務網(http://cdc.cma.gov.cn/ home.do)呼圖壁氣象站(44.08N,86.49E)4月22日、6 月9日、8月26日、10月15日的蒸發皿水面蒸散量值與呼圖壁遙感蒸散量分布圖上統計出的最大值進行比較。水面蒸發反映了一定區域特定時段內蒸發潛能, 可視為實際蒸散發的上限, 如果估算的陸面蒸散發高于蒸發皿蒸發, 則結果不合理[19], 由下表可知, 遙感估算的蒸散發量大都低于蒸發皿觀測值,但在10月15日的模型計算結果當中, 蒸散量值出現了高于水面蒸發的情況, 有可能是受到了遙感影像中的云的影響。

4.2.2 空間位置上的可靠性驗證

在蒸散發分布圖上提取氣象站點所對應經緯度的像元值和利用該氣站數據以彭曼公式計算的蒸散量值進行對比。

一般, 蒸發皿觀測的蒸散量總是保持最大, 其次是彭曼方法估算的蒸散量。在半干旱區, 陸面水分條件大多時候不能滿足蒸散需要, 實際蒸發要明顯比彭曼法估算的蒸散量低。

圖4 蒸散發量直方圖Fig. 4 Histogram of evapotranspiration

表2 呼圖壁縣水面蒸散量值對比Tab. 2 Comparision of Hutubi water surface on evapotranspiration map

表3 三種蒸散發計算結果對比Tab. 3 Comparision of three kinds of evapotranspiration calculated results

由以上兩種方式可以看出, SEBS模型所計算的蒸散發量與氣象站所測水面蒸發量和彭曼公式的計算結果三者變化趨勢相吻合, 盡管在驗證方面實測數據較為缺乏, 但表中仍然可以認為 SEBS模型在呼圖壁縣的蒸散量計算結果基本符合實際。

5 結論

(1) 使用 Landsat8影像, 利用SEBS 模型計算出2013年4月22日、6月9日、8月28日、10月15日的呼圖壁縣蒸散量, 并對SEBS的模擬結果與呼圖壁縣氣象站的觀察結果和利用彭曼公式計算的結果進行了比較, 表明SBES模型的結果是合理的,可以適用于較小區域蒸散發估算有一定適用性。

(2) 呼圖壁縣蒸散發量時空變化顯著, 春季的蒸散量較低, 夏季最大。南部山區具有較高的蒸發量, 中部農業灌溉地區包含城鎮部分, 夾雜著植被和水體等具有較高蒸散發的地類物型, 也呈現出較高的蒸發量, 北部地區多為戈壁和沙漠, 供給蒸發的水分不足, 該區域蒸散量在各個時期都較小, 估算結果合理地反映了蒸發量的時空差異。

(3) 利用 Landsat8影像周期性獲取可視化蒸散發量的空間分布, 對于縣域水資源管理和分配, 保護水資源平衡, 促進農業發展都具有重要意義。

[1] 張榮華, 杜君平, 孫睿. 區域蒸散發遙感估算方法及驗證綜述[J]. 地球科學進展, 2012, 27(12): 1295–1307.

[2] 傅國斌, 劉昌明. 遙感技術在水文學中的應用與研究進展[J]. 水科學進展, 2001, 12(4): 547–559.

[3] 高彥春, 龍笛. 遙感蒸散發模型研究進展[J]. 遙感學報, 2008, 12(3): 515–528.

[4] 郭曉寅, 程國棟. 遙感技術應用于地表面蒸散發的研究進展[J]. 地球科學進展, 2004, 19(1): 107–114.

[5] SU Z. The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes[J]. Hydrology and Earth System Sciences, 1999, 6(1): 85–100.

[6] 何延波, 王石立. 遙感數據支持下不同地表覆蓋的區域蒸散[J]. 應用生態學報, 2007, 18(2): 288–296.

[7] 李發鵬, 徐宗學, 李景玉. 基于MODIS數據的黃河三角洲區域蒸散發量時空分布特征[J]. 農業工程學報, 2009, 25(2): 113–120.

[8] 張雨航, 王曉林, 胡光成. 基于MODIS數據的海流兔河流域蒸散量的計算[J]. 地球科學--中國地質大學學報, 2012, 37(2): 376–380.

[9] 金曉媚, 郭任宏, 夏薇. 基于MODIS數據的柴達木盆地區域蒸散量的變化特征[J]. 水文地質工程地質, 2013, 40(6): 8–13.

[10] 田國珍, 武永利. 基于遙感蒸散發模型的山西省干旱監測研究[J]. 中國農學通報, 2013, 29(36): 160–166.

[11] 唐婷, 冉圣宏, 談明洪. 京津唐地區城市擴張對地表蒸散發的影響[J]. 地球信息科學學報, 2013, 15(2): 233–240.

[12] MA W, HAFEEZ M, ISHIKAWA H, et al. Evaluation of SEBS for estimation of actual evapotranspiration using ASTER satellite data for irrigation areas of Australia [J]. Theoretical and Applied Climatology, 2013, 112(3/4): 609–616.

[13] 周劍, 程國棟, 李新, 等. 應用遙感技術反演流域尺度的蒸散發[J]. 水利學報, 2009, 40(6): 679–687.

[14] 曾麗紅, 宋開山, 張柏, 等. 應用Landsat數據和SEBAL模型反演區域蒸散發及其參數估算[J]. 遙感技術與應用, 2008, 23(3): 255–263.

[15] 徐涵秋, 唐菲. 新一代Landsat系列衛星: Landsat8遙感影像新增特征及其生態環境意義[J]. 生態學報. 2013, 33(11): 3249–3257.

[16] 汲旭生, 袁昆, 筱佳, 等. 利用Landsat 8 影像制作解譯樣本點輔助高分影像分類可行性的研究[J]. 測繪與空間地理信息, 2014, 37(6): 51–55.

[17] 焦志敏, 張曉麗, 李法玲, 等.Landsat8多光譜波段紋理對葉面積指數的影響分析[J]. 地理與地理信息科學, 2014, 30(3): 42–45.

[18] 孫挺. 基于農業生物質利用下的土地潛力分析——以呼圖壁縣為例[D]. 新疆農業大學碩士畢業論文, 2013, 6.

[19] 周峰, 王文, 王曉剛. 基于植被指數-地表溫度特征空間的伊河流域蒸散發量估算[J]. 地理與地理信息科學, 2013, 29(002): 116–120.

[20] 張強, 張之賢, 問曉梅, 等. 陸面蒸散量觀測方法比較分析及其影響因素研究[J]. 地球科學進展, 2011, 26(5): 538–547.

Spatial and temporal distribution of evapotranspiration in the hutubi County based on Landsat8 data and SEBS model

ZHANG Yuan1,2, ZHENG Jianghua1,2,*, LIU Zhihui1,2,3, YAO Junqiang1,2
1.College of Resources and Environment Sciences,Xinjiang University,Urumqi,Xinjiang830046,China2.Key Lab for Oasis Ecosystem of MOE,Xinjiang University,Urumqi,Xinjiang830046,China3.Institute of Arid Ecology and Environment,Xinjiang University,Urumqi Xinjiang830046,China

The evapotranspiration of Hutubi County was estimated by Landsat8 and meteorological observations based on SEBS model. Estimated evapotranspiration of four days on April 22th, June 9th, August 28thand October 15th, and the distribution maps of were retrieved. The results could reflect season change; for instance, the results of farmland were 1.938, 3.136, 2.641, 1.314 mm·day-1on four days respectively. The overall trend was the northern desert plain <the central region < the southern mountains. The maximum value in four days appeared on June 9th in southern mountain area which was 4.128 mm·day-1. The results showed the SEBS model has accuracy in estimating evapotranspiration of arid and semi-arid climate and it can capture temporal and spatial change on typical northern slope of Tianshan Mountain.

SEBS model; evapotranspiration; Landsat8; northern slope of Tianshan Mountain; County

10.14108/j.cnki.1008-8873.2016.02.005

K90

A

1008-8873(2016)02-026-07

2014-11-24;

2015-09-22

水利部公益性行業科研專項(201301103)資助; 教育部促進與美大地區科研合作與高層次人才培養項目(117-40101)

張圓(1990—)女, 漢, 籍貫: 山西, 碩士研究生, 主要從事水資源管理及GIS應用等方面研究, E-mail: zy13659958670@aliyun.com

*通信作者:鄭江華(1973—), 男, 教授, 碩士生導師, 主要從事資源監測和3S技術應用研究, E-mail: zheng_jianghua@126.com

猜你喜歡
區域模型
一半模型
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
關于四色猜想
分區域
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 在线国产资源| 国产系列在线| 国产精品人莉莉成在线播放| 亚洲无码高清免费视频亚洲| 中文字幕有乳无码| 久久影院一区二区h| 全裸无码专区| 麻豆AV网站免费进入| 高清免费毛片| 亚洲成人在线播放 | 老色鬼久久亚洲AV综合| 国产欧美性爱网| 日本国产在线| 蜜臀AVWWW国产天堂| 欧美精品在线免费| 国产专区综合另类日韩一区 | 国产成人午夜福利免费无码r| 2021最新国产精品网站| 国产伦片中文免费观看| 国产成人亚洲综合a∨婷婷| 亚洲精品777| 极品性荡少妇一区二区色欲| 思思热在线视频精品| 欧美亚洲一区二区三区导航| 久久国产精品电影| 伊人欧美在线| 国产精品所毛片视频| 亚洲天堂视频在线免费观看| 久久精品中文字幕免费| 久久国产高潮流白浆免费观看| 国产制服丝袜无码视频| 狠狠色综合久久狠狠色综合| 中文字幕中文字字幕码一二区| 精品久久蜜桃| 亚洲一区免费看| 国产第一页屁屁影院| 黄色网站不卡无码| 久久久久久国产精品mv| 91娇喘视频| 国产色婷婷| 免费观看亚洲人成网站| 亚洲欧美成aⅴ人在线观看| 国产h视频在线观看视频| 欧美一级夜夜爽www| 国产靠逼视频| 婷婷午夜影院| 国产精品成人一区二区不卡 | 伊人久久婷婷五月综合97色| 午夜视频免费试看| 免费可以看的无遮挡av无码| 成人免费网站在线观看| 欧美翘臀一区二区三区| 精品国产网站| 99偷拍视频精品一区二区| 婷婷丁香色| 九九久久精品免费观看| 91久久偷偷做嫩草影院电| 欧美一级高清视频在线播放| 在线a视频免费观看| 久久精品视频亚洲| 国产欧美日韩视频怡春院| 亚洲日本一本dvd高清| 欧美国产精品不卡在线观看| 97在线观看视频免费| 无码电影在线观看| 五月天福利视频| 91区国产福利在线观看午夜| 69av在线| 亚洲黄色激情网站| 国产精品尤物在线| 久久婷婷六月| 国产产在线精品亚洲aavv| 91无码人妻精品一区二区蜜桃| 欧美黄网站免费观看| 99精品一区二区免费视频| 美女高潮全身流白浆福利区| 国产91特黄特色A级毛片| 国产女人在线| 日韩小视频在线播放| 最新亚洲人成无码网站欣赏网| 成人在线综合| 色噜噜在线观看|