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

城市地表溫度反演及其與下墊面定量關系分析—以西安市為例

2015-02-21 03:47:54毛文婷王旭紅祝明英蔡靜程德強

毛文婷,王旭紅,祝明英,蔡靜,程德強

西北大學城市與環境學院,陜西西安710127

城市地表溫度反演及其與下墊面定量關系分析—以西安市為例

毛文婷,王旭紅*,祝明英,蔡靜,程德強

西北大學城市與環境學院,陜西西安710127

城市地理與氣象條件、人為因素共同決定了當今城市的熱島現象,其中后者的影響更為重要,而下墊面又是主要的人為因素之一。本文以西安市為研究區域,利用2014年的landsat8數據,反演出地表溫度、植被和不透水面等信息,對城市熱島效應與下墊面關系進行定量分析,結果表明:作為下墊面重要表征組分的植被、不透水面與地表溫度有著密切的關系,其中植被指數與地表溫度之間的相關系數為-0.912,呈負相關關系;而不透水面指數與地表溫度之間的相關系數是0.747,呈正相關關系;此外,還對地表溫度、植被指數和不透水面指數進行了多元回歸分析。由此說明,植被對城市熱島效應的改善起到積極的作用,而不透水面會增加地表溫度,從而加劇熱島效應,因此可根據植被、不透水面這兩大因子對城市生態環境進行評估,從而制定相應的生態決策。

城市熱島效應;下墊面;地表溫度;不透水面;西安市

城市熱島(Urban Heat Island,UHI)反映的是城市與郊區之間的溫度差異,城市化是促進當今城市熱島形成的最主要原因。中國經濟網2015年1月20日報道,國家統計局發布報告顯示,2014年我國城鎮化率達到54.77%。隨著我國城市化水平的不斷提高,城市熱島現象日益加劇,城市熱島效應對城市生態系統產生各種影響,如城市高溫、暖冬、暴雨等異常城市氣象,與其他城市氣候現象相互作用會產生諸如“城市雨島”、“城市霧島”、“城市污島”等現象,城市持續高溫對大氣臭氧層也會造成相應破壞。城市熱島效應不僅影響人們正常的生產生活和工作學習,還阻礙人類生存質量提高和城市發展。因此,城市熱島現象已成為專家學者研究的熱點之一。

城市熱島現象反映的是城市的局地氣候特征,主要自然和人為兩大因素決定。自然因素主要指當地地理、氣象等,人為因素主要包括下墊面、人為熱源及大氣污染、城市形態及景觀格局。其中下墊面是人類活動最為直接的對象,因而是影響城市地表溫度(Land surface temperature)進而成為影響城市熱島效應最為重要的因子之一。不少學者對城市熱島效應及其影響因子之間的關系進行了研究,例如:唐菲等以福州市蒼霞片區為例,分析了在舊城改造下不透水面地表溫度呈指數相關關系,高不透水面比例地區溫度上升比低不透水面比率地區更快[1];錢樂祥等使用“單窗算法”利用3個不同時段的TM/ETM+數據獲取了珠江三角洲核心區域的地表溫度分布圖,分別分析了歸一化水汽指數、歸一化植被指數與地表溫度之間的關系并對兩指數進行了對比[2];Kok Chooi Tan利用Landsat數據使用經過兩種大氣校正法之后影像分別對歸一化植被指數(NDVI)和地表溫度(Land surface temperature,LST)之間的關系進行了評價[3];Zhang等利用ETM+(Enhanced thematic mapper plus)數據反演南京地表溫度,研究在不同尺度下提取的植被和不透水面與地表溫度的線性關系[4];賈寶全等[5]利用TM影像研究了西安市2006年和2010年的城市熱島的變化并分析了城市綠地對城市熱島效應的影響。以上研究較注重于單一因子對城市地表溫度的影響,而通過分析多因子對地表溫度的綜合影響的研究相對較少。本文在上述研究的基礎上,以西安市作為主要研究區,探討城市下墊面組分包括植被和不透水面兩大因子對城市地表溫度的影響,從而為降低城市地表溫度,減緩城市熱島效應,改善城市生態環境提供相應的參考對策。

1 研究區與數據源

西安市地處黃河流域中部的關中平原地區,107.40°E~109.49°E和33.42°N~34.45°N之間,是陜西省省會城市,現轄新城、碑林、蓮湖、雁塔、灞橋、未央、閻良、臨潼、長安、高陵10個區,藍田、周至、戶縣3個縣,氣候屬暖溫帶半濕潤大陸性季風氣候,是國家重要的科研、教育和工業基地,我國西部地區重要的中心城市,世界歷史文化名城。本文主要通過城市地表溫度的反演,研究地溫與下墊面組分之間的關系。因此,選取西安市繞城高速公路以內的中心城區作為主要研究區,如圖1。該區是西安市城市發展的核心區域,道路網為環狀與放射狀相結合,人口密度大。

研究數據為2014年1月3日、5月11日和8月15日的三期Landsat8影像,PATH/ROW為127/36,數據來源于NASA官方網站,主要采用陸地成像儀(OLI)的4個波段,分別是綠波段(Band3)、紅波段(Band4)、近紅外波段(Band5)和中紅外波段(Band6),空間分辨率為30 m;熱紅外傳感器(TIRS)的熱紅外波段(Band10),空間分辨率是120 m。本文對Landsat8數據進行了預處理,包括使用FLAASH大氣校正法做輻射校正,各項定標數據可從影像頭文件中獲取,以1:5000地形圖為參考做幾何校正。

圖1 研究區示意圖Fig.1 The schematic study area

2 研究方法

2.1 植被指數

自勞斯(Rouse)等人[6]引入歸一化差值植被指數(Normalized Difference Vegetation Index,NDVI)在遙感影像上提取植被信息以來,該方法得到了廣泛應用。該指數的原理是植被在多光譜遙感數據的近紅外波段具有高反射性,在紅光波段具有強吸收性。NDVI計算公式如下:

式(1)中,NIR、Red分別對應于landsat8數據的第五波段、第四波段,結果見圖2((a)是2014年1月3日的NDVI,(b)是2014年8月15日的NDVI)。從圖2(a)和(b)都可以看出植被指數高的區域集中于該研究區的外部邊沿,而低的區域主要集中于城市中心地區;對比兩圖也可以看出,圖2(b)代表的夏季NDVI最大值較圖2(a)代表的冬季最大值大很多;圖2(a)和(b)兩期影像反演出的NDVI平均值分別為0.085、0.299,這是因為夏季氣溫高,植被茂密,而冬季氣溫低,植被稀疏。

2.2 不透水面指數

不透水面[7](Impervious surfaces)是指水不能直接通過且不能下滲到土壤中的物質,包括瀝青、水泥、瓦片等材料構成的建筑物、路面和停車場等。利用遙感手段提取不透水面信息的方法有很多,本文采用徐涵秋[8,9]提出的歸一化差值不透水面指數(normalized difference impervious surfure index,NDISI),該指數能夠較好的提取不透水面信息。

式中,NIR、MIR1、TIR分別為遙感影像的近紅外波段、中紅外1波段和熱紅外波段,分別對應landsat8的第5波段、第6波段、第10波段;這里還引入了另一個歸一化差值指數即改進型歸一化水體指數[10,11]MNDWI(modified normalized difference water index)。

其中,Green為綠光波段,即為Landsat8的第3波段,結果見圖2((c)是2014年1月3日的MNDISI,(d)是2014年8月15日的MNDISI)。對比圖2(c)和(d)可以看出,圖2(c)代表的冬季MNDISI較圖2(d)代表的夏季MNDISI值高的范圍分布廣,這可從兩個季節的MNDISI的均值可以看出,1月和8月MNDISI均值分別為0.5696、0.470。

圖2 NDVI和MNDISIFig.2 NDVI and MNDISI

2.3 地表溫度的反演

地表溫度的反演方法有很多種,常見的有大氣校正法[12](也稱“輻射傳輸方程法”;Radiative Transfer Equation,RTM)、單通道算法[13,4](Single-channel method)、單窗算法[15](Mono-window Algorithm)、劈窗算法[16](也稱“分裂窗算法”,Split-window Algorithm),其中前三種方法都是針對單一熱紅外通道進行地溫反演的方法,后一種針對多通道的熱紅外地溫反演。大氣校正法由于需要詳細大氣剖面參數,而這些參數不容易獲取,因而這種方法較少使用;單窗算法和單通道算法估測大氣參數的情況下也可獲取較高的地溫反演精度,但是在有實時氣象資料提供的前提下,單窗算法的反演精度優于單通道算法[17]。對于Landsat8TIRS有兩個熱紅外波段,是和MODIS、NOAA/AVHRR等具有相類似的多通道熱紅外波段的傳感器,故而既可以采用劈窗算法;也可以采用已被證實具有較高反演精度的單窗算法。本文采用了單窗算法[15],對熱紅外10波段進行了地表溫度反演,公式如下:

式中,Ts是地表溫度(K);i為第i波段;ai和bi是常量,在這里選擇0~40℃時的情況,其中a10=-60.919,b10=0.428;C和D是中間變量,C=ετ,,其中ε是地表比輻射率,τ是大氣透射率;Ti是衛星傳感器上傳感器所探測到的像元亮度溫度(K);Ta為大氣平均作用溫度(K)。2.3.1大氣平均作用溫度Ta大氣平均作用溫度[18]與地面附近(一般為2 m處)氣溫T0(K)存在以下線性關系(如表1)。根據在中國氣象科學數據共享服務網(http://cdc.cma.gov.cn/home.do)提供的西安市實測溫度數據計算得出2014年1月3日和2014年8月15日氣溫分別為5℃、25℃,即T0為278.15 K、298.15 K,由此計算Ta為272.888 K、292.161 K。

表1 大氣平均作用溫度計算Table1 Calculation of the average temperature of the atmosphere

2.3.2 大氣透射率τ估計大氣透射率[18]是指通過大氣(或某氣層)后的輻射強度與入射前輻射強度之比,主要取決于大氣中的水分含量。根據西安市的地理位置,通過水分含量估算2014年1月3日和2014年8月15日大氣透射率τ分別約為0.909、0.883。

2.3.3 亮度溫度Ti亮度溫度反演采用普朗克函數[19]計算:

式中,Li為輻射亮度值,landsat8 TIRS熱紅外波段的Ki,1、Ki,2值可通過查詢影像頭文件查詢到,其中K10,1=774.89 W/m2·μm·sr,K10,2=1321.08 K。

2.3.4 地表比輻射率ε地表比輻射率的估算比較常用的一種方法是先對遙感圖像進行分類,將地表分為不同的覆蓋類型,在根據實測或者經驗值的地物比輻射率給各個地表覆蓋類型賦予不同的值,從而生成地表比輻射率圖像。目前,已有一些比輻射率數據庫,如MODIS UCSB比輻射率庫(http://www.icess.ucsb.edu/modis/EMIS/html/em.html)等。此外,還可利用NDVI計算地表比輻射率,這是由于NDVI的對數與地表比輻射率存在線性相關性,利用NDVI的閾值進行地表分類,然后給各個地表覆蓋類型賦予不同的值,該方法得到受到了很多專家學者的關注并加以運用[20][21]。本文采用覃志豪等[20]提出的現將地表分成水體、自然表面和城鎮區,分別針對3種地表類型計算地表比輻射率:

式中,NDVIs為完全被裸土或無植被覆蓋區域的NDVI值,NDVIv則代表完全被植被所覆蓋的像元的NDVI值,即純植被像元的NDVI值。取經驗值NDVIv=0.70和NDVIs=0.05,即當某個像元的NDVI值大于0.70時,Pv取值為1,;當NDVI小于0.05,Pv取值為0。

根據以上對單窗算法計算參數的計算,對Landsat8 THRIS Band10進行反演,反演出亮度溫度、地表溫度,并對二者進行差值處理后的分級結果如圖3所示。2014年1月3日和2014年8月15日反演的地表溫度與亮度溫度差值的平均值分別為2.32℃、2.689℃,可以看出利用單窗算法反演出的地表溫度比亮度溫度高,夏季反演的溫度差值比冬季反演的溫度差值高;由圖3(e)(f)可以看到兩期影響反演的地表溫度與亮度溫度差值分別集中在2.4~3.92℃、1.4~3.2℃。

3 結果與分析

3.1 城市地表溫度的空間分布特征

地表溫度的反演結果:研究區兩期影像地表溫度范圍分別在-1.46~13.84℃、24.67~46.78℃,平均溫度分別為6.55℃、31.26℃。地表溫度的空間分布特征如圖3所示:按均值±整數倍標準差將溫度分割為5個等級。

圖3(b)、(d)中,研究區地表溫度分別集中于3.5~9.5℃之間、24.6~34℃之間,兩景不同時間的影像反演出來的地表溫度高溫分布區域和低溫分布區域大致相同;其中高溫區主要集中于該研究區的中部、中北部以及西部,而沿繞城高速的外圍區則主要屬于該研究區低溫區,這與城市本身中心氣溫高于郊區的特征相符。此外,該研究區的低溫主要集中于灞河流域附近,高溫區主要集中在新城區東部、蓮湖區中部,這是因為這些區域建有一些工廠和加油站。

圖3 亮度溫度、地表溫度反演結果與二者差值Fig.3 Inversion effect of brightness temperature and surface temperature and their difference

3.2 地表溫度與NDVI、MNDISI的定量關系

目前,已經有很多學者對地表溫度與植被、不透水面這兩大因子之間的定量關系進行了研究。張小飛等[22]以深圳市為例,從土地利用角度將下墊面進行詳細劃分、使用線性光譜分離模型提取植被蓋度、反演亮度溫度,確立了植被覆蓋度與地表溫度呈負相關關系且呈現一定的線性關系,并分析了不同分辨率和下墊面類型對二者關系的影響;Xu[23]提出一種新的可以有效提取不透水面的指數,并將其與反演出的地表溫度進行相關性分析,比較好的揭示了二者之間的關系;歷華等[24]以MODIS為數據源,提取長株潭地區的NDBI、NDVI和地表溫度,分析各自的時空分布特征并對地表溫度與NDBI、NDVI的定量關系分別進行了研究。本文借鑒以上學者的研究方法,在上文中已經采用指數法提取了植被指數NDVI、不透水面指數MNDISI,通過單窗算法反演了地表溫度,并分析了三者之間的空間分布特征。在此基礎上,本文采用如下方式對地表溫度與NDVI、MNDISI之間的定量關系進行研究。

本文利用2014年5月11日Landsat 8影像反演出地表溫度、NDVI和MNDISI進行隨機采樣,樣區尺度為60*60個柵格(每個柵格尺度為30 m*30 m),樣本個數為30個,統計每一個樣區的平均值。對NDVI、MNDISI與地表溫度Ts分別計算其相關系數,做統計回歸分析;然后利用回歸方法來綜合分析NDVI、MNDISI對地表溫度的影響,從而實現從二維和三維的角度分別研究下墊面組分對地表溫度的影響(表2)。

NDVI與地表溫度Ts之間相關系數R為-0.912,二者進行一元線性回歸分析(圖4),發現用指數模型擬合二者的關系比較好,如圖4所示。從相關系數及圖4可以看出地表溫度與NDVI高度相關并呈負相關關系,即植被覆蓋率高的地區地表溫度Ts低,反之,地表溫度Ts高,具體統計NDVI增加與地表溫度的遞減幅度如表3。由表3可以看出,隨著NDVI值增加,地表溫度降低的幅度先大后小,說明小幅度的增加植被覆蓋率就可以達到降低地表溫度的效果。地表溫度Ts與MNDISI之間的相關系數為0.747,二者之間用對數模型可以達到較好的擬合效果,如圖4。從相關系數和對數模型可以看出,地表溫度與MNDISI呈正相關關系,即隨著不透水面指數值的增大,地表溫度也會隨之上升。由表3可以看出,地表溫度Ts的上升幅度是隨著不透水面指數值的增加而減小的。因此,應當重點關注不透水面這一下墊面組分對地表溫度的影響。同時在表3中對比植被和不透水面對地表溫度的影響,在植被指數或不透水面指數增加相同比例時,后者對地表溫度的影響更為顯著。

表2 植被指數(NDVI)、不透水面指數(MNDISI)和地表溫度(Ts)的關系方程Table 2 Regression models between NDVI,MNDISI and land surface temperature(Ts)

圖4 植被指數(NDVI)、不透水面指數(MNDISI)與地表溫度(Ts)的回歸關系Fig.4 Regression analysis of the relationship of Ts with NDVI and MNDISI

表3 植被指數(NDVI)值、不透水面指數(MNDISI)值及其所對應的地表溫度(Ts)Table 3 The value of NDVI,MNDISI and their corresponding Ts

多元回歸分析進一步揭示了作為下墊面重要組分的植被和不透水面都是影響地表溫度的重要因素,相關系數達到0.917。植被指數NDVI每增加0.1,地表溫度Ts下降0.67℃;不透水面指數MNDISI每增加0.1,溫度將上升1.97℃;如果在植被指數增加0.1,同時不透水面指數減少0.1,則地表溫度可降低2.64℃,這說明,植被和不透水面共同影響地表溫度,前者增加同時后者減少才能有效降低城市地表溫度。

4 結論

以上研究表明:(1)西安市繞城高速以內的主城區植被與地表溫度相關系數為-0.912,呈指數型負相關關系,且植被指數的值在增加的初期對地表溫度的降低效果更顯著;(2)不透水面與地表溫度相關系數為0.747,呈對數型負相關關系,同時在不透水面指數的值在增加的初期對地表溫度的影響更為顯著;(3)植被指數和不透水面指數在增加或減少相同比例時對地表溫度的影響更顯著;(4)三者之間進行多元回歸,相關系數為0.917,說明作為下墊面組分的植被和不透水面都是影響地表溫度的重要因素。因此,在增加城市植被覆蓋率的基礎上同時減少不透水面可以綜合降低城市地表溫度,并且采用這種方式達到的降溫效果優于改善植被或不透水面單一要素分布狀況時所達到的降溫效果。此外,通過此研究發現,由于影響城市地表溫度的因子是復雜多變的,對多因子的綜合分析對降低地表溫度,改善城市熱島效應的效果可能更加顯著。

[1]唐菲,徐涵秋.舊城改造與城市熱島效應關系的遙感研究——以福州市蒼霞片區為例[J].地理科學,2011,31(10):1228-1234

[2]錢樂祥,崔海山.歸一化水汽指數與地表溫度的關系[J].地理研究,2008,27(6):1358-1366

[3]Kok Chooi Tan,Hwee San Lim,Mohd Zubir MatJafri,et al.A comparison of radiometric correction techniques in the evaluation of the relationship between LST and NDVI in Landsat imagery[J].Environmental Monitoring and Assessment, 2012,184(6):3813-3829

[4]Zhang X,Zhong T,Wang K,et al.Scaling of impervious surface area and vegetation as indicators to urban land surface temperature using satellite data[J].International Journal of Remote Sensing,2009,30(4):841-859

[5]賈寶全,邱爾發.基于TM衛星遙感影像的西安市城市熱島效應變化分析[J].干旱區研究,2013,30(2):347-355

[6]Rouse JW,Haas RH,Schell JA,et al.Monitoring Vegetation Systems in the Great Plain with ERTS[C]//Proceedings of the Third Earth Resource Technology Satellite—1 Symposium.Greenbelt:NASA SP-351,1974:3010-317

[7]Arnold CL,Gibbons CJ.Impervious Surface Coverage:The Emergence of a Key Environmental Indicator[J]. Journal of the American Planning Association,1996,62(2):243-258

[8]徐涵秋.一種快速提取不透水面的新型遙感指數[J].武漢大學學報:信息科學版,2008,33(11):1150-1153

[9]Xu Hanqiu.Rule-based impervious surface mapping using high spatial resolution imagery[J].International Journal of Remote Sensing,2013,34(1):27-44

[10]徐涵秋.利用改進的歸一化差異水體指數(MNDWI)提取水體信息的研究[J].遙感學報,2005,9(5):589-595

[11]Xu Hanqiu.Modification of normalised difference water index(NDWI)to enhance open water features in remotely sensed imagery[J].International Journal of Remote Sensing,2006,27(14):3025-3033

[12]覃志豪,Arnon Karnieli.用NOAA-AVHRR熱通道數據演算地表溫度的劈窗算法[J].國土資源遙感,2001(2):33-42

[13]Jiménez-Mu noz JC,Sobrino JA.A generalized single-channel method for retrieving land surface temperature from remote sensing data[J].Journal of Geophysical Research,2003,108(D22):4688

[14]Jiménez-Mu noz JC,Cristobal J,Sobrino JA,et al.Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval From Landsat Thermal-Infrared Data[J].IEEE Transactions on Geoscience and Remote Sensing,2009,47(1):339-349

[15]覃志豪,Zhang Minghua,Arnon Karnieli,等.用陸地衛星TM6數據演算地表溫度的單窗算法[J].地理學報,2001,56(4):456-466

[16]Sobrino JA,Li Zhaoliang,Stoll MP,et al.Improvements in the split-window technique for land surface temperature determination[J].IEEE Transactions on Geoscience and Remote Sensing,1994,32(2):243-253

[17]丁鳳,徐涵秋.TM熱波段圖像的地表溫度反演算法與實驗分析[J].地球信息科學,2006,8(3):125-130

[18]覃志豪,Li Wenjuan,Zhang Minghua,等.單窗算法的大氣參數估計方法[J].國土資源遙感,2003(2):37-43

[19]Sospedra F,Caselles V,Valor E.Effective wave number for thermal infrared bands application to Landsat T M[J].International Journal of Remote Sensing,1998,19(11):2105-2117

[20]覃志豪,李文娟,徐斌,等.陸地衛星TM6波段范圍內地表比輻射率的估計[J].國土資源遙感,2004,(3):28-32

[21]Sobrino JA,Raissouni N,Li Zhaoliang.A Comparative Study of Land Surface Emissivity Retrieval from NO AA Data[J].Remote Sensing of Environment,2001,75(2):256-266

[22]張小飛,王仰麟,吳健生,等.城市地域地表溫度-植被覆蓋定量關系分析——以深圳市為例[J].地理研究,2006,25(3):369-377

[23]Xu Hanqiu.Analysis of Impervious Surface and its Impact on Urban Heat Environment using the Normalized Difference Impervious Surface Index(NDISI)[J].Photogrammetric Engineering and Remote Sensing:Journal of the American Society of Photogrammetry,2010,76(5):557-565

[24]歷華,柳欽火,鄒杰.基于MODIS數據的長株潭地區NDBI和NDVI與地表溫度的關系研究[J].地理科學,2009,29(2):262-267

Analysis on the Urban Land Surface Temperature Inversion and Its Quantitative Relationship with Underlying Surface-Taking Xi'an City as a case

MAOWen-ting,WANGXu-hong*,ZHUMing-ying,CAIJing,CHENGDe-qiang
College of Urban and Environmental Science/Northwest University,Xi'an 710127,China

Urban geographical and meteorological conditions and human factors jointly determines the current urban heat island phenomenon.The more important effect comes from human factors,however the underlying surface is one of the main human factors.This research took Xi'an City,provincial capital of Shanxi China as research area and applied the 2014 landsat8 data to extract land surface temperature,vegetation and impervious surface information and analyzed the quantitative relationship between urban heat island effect and vegetation,underlying surface.The results showed that vegetation and impervious surface which were important components of underlying surface had close relations with land surface temperature.The correlation coefficient between vegetation index and land surface temperature was negative 0.912, and the correlation coefficient between the land surface temperature and impervious surface index was 0.747.In addition,the land surface temperature,vegetation index and impervious surface index had been carried on a multiple regression analysis. Above research indicated that vegetation played a positive role in the improvement of the urban heat island effect,while impervious surface increased the land surface temperature to aggravate the heat island effect.So on the basis of vegetation and impervious surface,the urban ecological environment could be evaluated,and the corresponding ecological decisions could be made.

Urban heat island effect;underlying surface;land surface temperature;impervious surface;Xi'an City

P423.2

A

1000-2324(2015)05-0708-07

2015-04-14

2015-05-14

國家自然基金:不同地貌類型區的遙感圖像信息容量空間差異性研究(41071271);陜西省自然科學基礎研究計劃:基于遙感圖像信息容量的城市熱島效應研究(2015JM4132)

毛文婷(1988-),女,山東泰安人,碩士研究生,研究方向為主要從事地理信息系統與遙感信息分析研究.E-mail:wentingmao@163.com

*通訊作者:Author for correspondence.E-mail:jqy_wxh@163.com

主站蜘蛛池模板: 亚洲综合激情另类专区| 久久成人免费| 婷婷中文在线| 好久久免费视频高清| 91国内外精品自在线播放| 成人年鲁鲁在线观看视频| 亚洲色欲色欲www在线观看| 国产专区综合另类日韩一区| 青青青草国产| 成AV人片一区二区三区久久| 免费可以看的无遮挡av无码| 精品国产一区91在线| 国产jizzjizz视频| 热99精品视频| 国产丝袜91| 国产经典免费播放视频| 亚洲无码37.| 中文字幕1区2区| 1769国产精品视频免费观看| 中文字幕自拍偷拍| 久久亚洲国产一区二区| 午夜无码一区二区三区| 黄色国产在线| 久久性妇女精品免费| 黄色国产在线| 国产主播喷水| 国产日韩欧美一区二区三区在线| 国产精品亚洲一区二区在线观看| 激情乱人伦| 国产欧美日韩专区发布| 亚洲三级影院| 日韩av手机在线| 欧美一级色视频| 亚洲天堂啪啪| 婷婷色狠狠干| 亚洲国产成人精品无码区性色| 干中文字幕| 这里只有精品在线播放| 欧美日韩在线成人| 国产91熟女高潮一区二区| 欧美成人怡春院在线激情| 免费国产一级 片内射老| 亚洲成人免费在线| 性做久久久久久久免费看| 国产免费黄| 激情在线网| 五月婷婷丁香综合| 黄色网址免费在线| 国产精品无码AV片在线观看播放| 久久无码av三级| 国产高颜值露脸在线观看| 99热国产这里只有精品9九| 日韩专区欧美| 日本高清在线看免费观看| 国产情侣一区二区三区| 3D动漫精品啪啪一区二区下载| 波多野结衣在线一区二区| 久久一本日韩精品中文字幕屁孩| 欧美97欧美综合色伦图 | 精品成人一区二区| 中文精品久久久久国产网址| 成人小视频在线观看免费| 天天综合网站| 日韩在线播放欧美字幕| 日韩精品无码不卡无码| 全免费a级毛片免费看不卡| 1769国产精品免费视频| 亚洲视频在线网| 深夜福利视频一区二区| 中文字幕在线一区二区在线| 99精品免费欧美成人小视频| 天天婬欲婬香婬色婬视频播放| 欧美激情视频二区| 欧美成人怡春院在线激情| 一本一道波多野结衣av黑人在线| 波多野结衣中文字幕久久| 美女裸体18禁网站| 日韩在线永久免费播放| 999在线免费视频| 老司国产精品视频91| 一级全免费视频播放| 久久午夜夜伦鲁鲁片无码免费|