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

基于MODIS影像及不透水面積的珠江三角洲熱島效應時空分析

2017-12-20 03:21:00賀麗琴蘇琳琳
自然資源遙感 2017年4期
關鍵詞:效應區域

賀麗琴, 楊 鵬, 景 欣, 晏 磊, 蘇琳琳,2

(1.北京大學空間信息集成與3S工程應用北京市重點實驗室,北京 100087;2.首都師范大學信息工程學院,北京 100037)

基于MODIS影像及不透水面積的珠江三角洲熱島效應時空分析

賀麗琴1, 楊 鵬1, 景 欣1, 晏 磊1, 蘇琳琳1,2

(1.北京大學空間信息集成與3S工程應用北京市重點實驗室,北京 100087;2.首都師范大學信息工程學院,北京 100037)

珠江三角洲城市群是我國經濟快速發展的地區,在經濟發展的同時,其形成的熱島效應也日益明顯。采用熱紅外遙感方法研究熱島效應,可很直觀地了解城市熱島的空間分布狀況。基于MODIS影像,運用劈窗算法,對珠江三角洲地區1 a中4季的地表溫度(land surface temperature,LST)進行反演; 并用不透水面積(impermeable surface area,ISA)對城市中心區和郊區進行劃分,最終定量獲得城市熱島效應的大小。研究結果表明,珠江三角洲地區存在較嚴重的熱島現象,夏季最為嚴重,冬季最輕,且城市之間有相互連接形成大片城市熱島效應的趨勢; 熱島效應的大小與歸一化植被指數值成負相關,與城市的經濟發展程度成正相關。研究成果可為珠江三角洲地區的城市發展規劃提供一定的生態指導。

熱島效應; 珠江三角洲; MODIS; 地表溫度(LST); 不透水面積(ISA); 影響因素

0 引言

城市熱島效應是城市經濟快速發展帶來的一種特殊的城市生態現象。城市人口增加、下墊面改變和人工廢棄物排放是造成城市熱島效應的主要原因。改革開放以來,珠江三角洲城市經濟的快速發展加速了珠江三角洲的城鎮化進程,形成了以廣州、佛山、深圳以及周圍城市為中心的珠三角都市群; 而這種趨勢對環境最直接的影響就是形成以城市為中心的“熱島”現象。近十幾a來,很多學者們通過定點溫度監測研究珠三角城市群的熱島現象。張晶晶等[1]利用自動氣象站得到了以佛山為中心的周邊地區隨月份變化的熱島變化規律; 竇浩洋等[2]利用自動氣象站對珠三角城市群進行了熱島分析,得到白天和黑夜的熱島變化規律; 曾俠等[3]利用氣象站對珠江三角洲多年的地表溫度(land surface temperature,LST)數據進行了分析,得出該地區的熱島效應隨年份變化的規律。

傳統的通過定點溫度監測來研究熱島效應的方法難以全面反映城市熱島分布情況,而通過衛星遙感數據反演的LST分布可較直觀地了解城市熱島的空間分布狀況。MODIS是在NASA EOS的Terra/Aqua衛星上搭載的新一代地球觀測傳感器,其數據在全世界范圍內可免費接收,包含36個波段,其中第31波段和第32波段可用于反演LST,并用于對熱島效應進行分析。目前,國內已有很多學者利用MODIS數據對城市熱島效應進行了研究,如歷華等[4]對長珠潭地區熱島效應的研究; 閆峰等[5]對上海市熱島效應的研究; 楊鵬等[6]對石家莊市熱島效應的研究,其研究結果均很好地反映了當地的熱島現象,并對城市環境改造有一定的指導作用。

在以往的研究中,計算城市熱島效應時,對城市區域和郊區的位置劃分不是很明顯[4-6],因而會造成計算熱島效應大小的不確定性。本文采用建筑的不透水面積(impermeable surface area,ISA)正是為解決這個問題而引入的。建筑的不透水面(impermeable surface)指的是道路、停車場、建筑物和人行道等一系列現代化的瀝青或混凝土設施[7]。ISA在某一地區的比例大小往往反映該地區的城鎮化程度,ISA比例越大,則城鎮化程度越高,往往越接近城市中心區域[8]。ISA比例的大小對生態環境的影響也不同,當ISA比例在[1%,10%)時,會對生態造成輕度影響; 當ISA比例在[10%,25%)時,對生態的影響明顯; 當ISA比例達25%以上時,生態環境開始退化[8-9]。

本文利用MODIS影像數據確定當地的LST,并采用ISA比例確定城市中心區和郊區的位置,旨在更加準確地計算出當地的熱島效應強度隨時空的變化特征。利用MODIS數據的劈窗反演算法[10],對珠江三角洲地區1 a中秋、冬、春、夏4季的LST進行反演,得出該地區城市1 a中4季所形成的熱島效應,并對其表現出的熱島現象進行分析,以期對當地的城市發展提供一定的指導。

1 數據源及LST反演

1.1 數據源

本文使用的MODIS數據為2014年10月15日上午11: 20,2015年1月2日上午11: 15,2015年4月15上午11: 20和2015年7月11日上午11: 20獲取的MODIS1B影像,空間分辨率為1 km。LST算法中用到的波段包括第1,2,19,31和32共5個波段。反演得到的LST可大致反映珠江三角洲地區1 a中秋、冬、春、夏4季白天的LST。反演地區的面積大小為200 km×200 km,其在廣東省的分布示意圖見圖1紅框內。

圖1 珠江三角洲地區位置示意圖Fig.1 Location of Pearl River Delta region

從圖1可知,研究區域包括廣州、佛山、東莞和深圳等經濟發達地區,不僅包括陸地,還包括一部分海域。

目前得到ISA比例有2種常用的方法: ①利用30 m高空間分辨率的Landsat影像得到ISA比例產品,但這種產品只覆蓋美國; ②利用夜間燈光數據(nighttime light data)和全球人口空間分類數據(LandScan)得到空間分辨率為1 km的ISA比例產品[7]。本文采用第2種方法得到ISA比例產品,且該數據的空間分辨率正好與MODIS反演出的LST產品吻合,并可在美國NOAA國家地理數據中心網站上得到。本文使用研究區域最新(2010年)的ISA比例產品(圖2)。

圖2 研究區域的ISA比例Fig.2 ISA proportion of study area

圖2中把ISA比例分成了5個層次,其中,25%以上表示會造成生態環境的退化,同時可以把這個比例當作城市區域與郊區的分界線(即圖中棕色區域與黃色區域的分界線),可以看出城市區域與郊區的分界線明顯; 而把ISA在[50%,100%]的范圍(即圖中紅色區域)看成是城市最中心的區域[11],本文利用此區域的LST減去ISA為25%分界線以外20 km郊區的LST,得到熱島強度大小。

表1示出依據以上ISA比例數據得到的研究區內各個城市的城市區域面積及城市中心區域面積(廣州和佛山因相距很近,且熱島效應相連成片,故放在一起)。

表1 各城市的城市區域及城市中心區域面積Tab.1 Area of urban region and urban core region of each city km2

在利用MODIS數據進行LST反演之前,首先對MODIS數據進行預處理,以保證反演結果滿足精度和地理要求。本文中所做的數據預處理包括: ①投影變換。MODIS產品數據的投影一般是Sinusoidal,為便于與地形數據進行空間分析,需進行投影變換。所用的MODIS1B數據是hdf格式,自帶有經緯度坐標信息,可自動進行投影變換。本文利用MODIS Swath Tool工具,使用MODIS03經緯度坐標數據,分別對熱紅外數據集和反射率數據集進行自動投影變換; ②去云處理。在本文所選的4景圖像中,個別地區存在少量的云覆蓋,因被云覆蓋的區域不能反映地表的反射率與輻射亮度,故將被云覆蓋的區域去掉,以免得到錯誤的LST。本文所用的去云算法是多光譜綜合去云算法[12-13],可很好地去除大片的云層以及不容易被去除的薄云。所用到的數據有MODIS1B第1,6,8,26,29和31波段共6個波段。

1.2 LST反演

目前,國內外發展了很多利用熱紅外波段反演LST的算法,其中劈窗算法(split-window algorithm)因所需參數少、運算速度快、精度高而被廣泛運用。一些學者針對MODIS數據提出了相應的劈窗算法[14-15],并取得了良好的效果。本文采用Mao等[16]提出的劈窗算法計算LST,具有較高的精度,其計算公式為

LST=[C32(B31+D31)-C31(D32+B32)]/(C32A31-C31A32),

(1)

其中,

A31=0.137 87ε31τ31,

(2)

B31=0.137 87T31+31.656 77ε31τ31-31. 656 77 ,

(3)

C31=0.137 87(1-τ31)[1+(1-ε31)τ31] ,

(4)

D31=31.656 77(1-τ31)[1+(1-ε31)τ31] ,

(5)

A32=0.118 49ε32τ32,

(6)

B32=0.118 49T32+26.500 36ε32τ32-26.500 36 ,

(7)

C32=0.118 49(1-τ32)[1+(1-ε32)τ32] ,

(8)

D32=26.500 36(1-τ32)[1+(1-ε32)τ32] ,

(9)

式中:T31和T32分別為MODIS第31和32波段的亮度溫度,可通過輻射傳輸方程計算得出;τ31和τ32分別為第31和32波段的大氣透過率;ε31和ε32分別為第31和32波段的地表比輻射率。

2 熱島效應及分析

2.1 熱島效應結果

運用上述數據進行預處理和LST反演之后,得到的珠江三角洲LST結果如圖3所示(圖 3(d)中右上角的黑色陰影為云)。

(a) 2014年10月15日(秋季)(b) 2015年1月2日(冬季)

(c) 2015年4月15日(春季) (d) 2015年7月11日(夏季)

圖3反演得到的4季LST結果

Fig.3ResultsofretrievedLSTof4seasons

在得到的LST結果中,秋季的最高和最低溫度分別為40.3℃和15.5℃,冬季的最高和最低溫度分別為25.6℃和11.7℃,春季的最高和最低溫度分別為40.8℃和20.6℃,夏季的最高和最低溫度分別為41.1℃和24.6℃。因珠江三角洲地區是我國快速城鎮化地區之一,城市熱島現象比較突出; 而由圖3中的LST反演結果可以看出,城市地區的溫度明顯比周圍地區高。圖3中許多高溫地區都連接成片,形成了大范圍的高溫區,呈現明顯以大城市為中心的“區域熱島”現象,尤其是以廣州和佛山為中心的區域熱島。

為了更好地研究城市的熱島現象,對反演出的LST結果做1條剖面線,剖面線經過中山、江門、佛山、廣州、東莞和深圳等主要城市(圖4)。

圖4 剖面線位置Fig.4 Location of section lines

得到的剖面線經過地區LST變化見圖5。

(a) 2014年10月15日(秋季)(b) 2015年1月2日(冬季)

(c) 2015年4月15日(春季) (d) 2015年7月11日(夏季)

圖54季的LST剖面

Fig.5LSTinsectionlinesof4seasons

從圖5可以看出,不論秋、冬、春、夏,各主要城市都形成了比較明顯的熱島,且城市之間(中山和江門,佛山和廣州,東莞和深圳)有連片形成區域熱島的趨勢。其中,佛山和廣州形成的熱島現象與其他幾個城市相比,不論是在溫度上、還是范圍上都要嚴重。

為了更準確地計算各個城市熱島效應的大小,將得到的4季LST圖與ISA比例圖進行疊置,用ISA在[50%,100%]范圍城市中心區域的平均LST減去ISA為25%分界線以外20 km郊區的LST,得到熱島強度的準確大小。采用上述方法進行計算,得到具體熱島效應大小(表2)。從表2中可以看出,秋季(10月)熱島效應最強的廣州/佛山地區比郊區要高出5.1 ℃,東莞地區熱島效應也較明顯,比郊區高出4.2 ℃,中山、江門地區的熱島效應相對要小很多; 冬季(1月)熱島效應最強的廣州/佛山比郊區高出2.4 ℃,中山、江門、東莞和深圳地區表現出相當的熱島強度,但總體上都比較小; 春季(4月)熱島效應最強的廣州/佛山地區比郊區高出4.9 ℃,中山、江門和東莞地區也都比郊區高出3.5 ℃左右; 夏季(7月)熱島效應最強的廣州/佛山地區比郊區高出5.8 ℃,中山、江門和東莞地區也表現出較高的熱島效應。由此可見,在珠江三角洲地區,熱島效應的強度與季節有關,夏季的熱島效應強度最大,春季和秋季次之,冬季最弱。

表2 利用ISA計算得到的各城市熱島效應大小Tab.2 Intensity of urban heat island of each city calculated by using ISA (℃)

2.2 熱島效應的影響因素

眾所周知,植被對一個地區的LST影響很大。為了研究珠江三角洲地區植被對LST的影響作用,本文計算出MODIS影像的歸一化植被指數(normalized difference vegetation index,NDVI),NDVI能較好地反映出當地的植被特征、 覆蓋度和地表綠化程度。NDVI的計算公式為

NDVI=(ρ2-ρ1)/(ρ2+ρ1),

(10)

式中ρ1和ρ2分別為MODIS第1波段和第2波段的反射率值。珠江三角洲地區的NDVI計算結果如圖6所示。

(a) 2014年10月15日(秋季)(b) 2015年1月2日(冬季)

(c) 2015年4月15日(春季)(d) 2015年7月11日(夏季)

圖64季的NDVI計算結果

Fig.6ResultsofNDVIof4seasons

在得到的珠江三角洲地區NDVI結果中,水域部分的NDVI值都小于0,陸地部分的NDVI值在0以上。陸地部分的NDVI值中秋、夏2季最大(均為0.77),冬、春2季相對較小(但相差不大),說明這一地區的植被1 a中4季變化較小; 而城區和郊區的NDVI值相差較大,特別是在夏季(圖6(d)),城區普遍在0.17~0.30之間,郊區普遍在0.53~0.65之間。對比圖3和圖6可看出,LST與NDVI呈明顯的負相關關系,這證實了NDVI明顯影響著LST的高低,即NDVI值越高,LST值越低; NDVI值越低,LST值越高。這也很好地說明了城區NDVI值很低而LST值卻很高的事實。

熱島效應的形成不僅與植被覆蓋相關,更取決于一個地區的經濟發展程度。珠江三角洲地區之所以形成如此明顯的熱島效應,主要是因為城市經濟的快速發展,使得城市聚集了大量人口進行工業生產,從而消耗了大量能源,排放出較多的人為熱量。本文利用研究區各個城市的國內生產總值(gross domestic product,GDP)來反映其經濟發展程度。表3為珠江三角洲地區2014年城市GDP與城市LST值之間的對應關系(其中廣州和佛山市因相距太近,熱島已經連接成片,所以當成一個地區來研究; 表中的LST是各城區中心的LST均值)。

表3 各城市國內生產總值與LST的對應關系Tab.3 Relationship between GDP and LST of each city

從表3可以看出,1 a中4季廣州/佛山的LST值均明顯比其他城市高,而廣州/佛山的GDP也是最高的,由此可見經濟發展的程度會影響到LST值的高低。東莞的GDP比中山和江門高,4季中秋季、春季和夏季3個季節的LST值都比后2個城市高,而在冬季3個城市LST值的差別并不明顯,這也說明了冬季的熱島效應最小。除廣州以外,深圳的GDP明顯比其他城市高出很多,但其LST值并不比其他城市高,甚至春、夏2季的LST值還更低,這可能是因深圳距離海面太近,由海風及海水的原因導致LST值下降。

3 結論

本文基于MODIS數據,利用劈窗算法對經濟發達的珠江三角洲地區4個季節的地表溫度(LST)進行了反演,得到的LST結果反映出當地1 a中4個季節的城市LST分布情況。為了更準確地計算熱島效應的大小,引入不透水面積(ISA),區分出城市區域和郊區。得出以下結論:

1)秋、冬、春、夏4個季節中,夏季的熱島強度最大,春季和秋季次之,冬季最弱。

2)NDVI值直接影響LST的分布,NDVI值的大小與LST值成負相關,NDVI值低使得城區形成明顯的熱島效應。

3)經濟發展程度最高的廣州熱島現象最為嚴重,這證實了城市發展與熱島效應的正相關性。

4)經濟發展程度相對較高的深圳表現出的熱島效應并不強,其原因可能是深圳距離海面太近而由海風及海水導致LST值下降,可見影響熱島效應的因素還有除植被覆蓋和經濟發展程度之外的其他原因。城市熱島效應大小與植被指數和國內生產總值的關系可為珠江三角洲地區的城市發展規劃提供一定的生態指導。

本文對ISA的引入可為國內學者分析熱島效應提供一種更為精確的計算方法; 但也存在一些缺陷,由于珠江三角洲地區常年陰雨天氣,每個季節只能得到1 d的MODIS晴天影像,這對分析熱島效應隨季節的變化難免會存在一定的片面性。

[1] 張晶晶,竇浩洋,張恩潔,等.珠江三角洲城市群熱島特征研究[C]//中國氣象學會.中國氣象學會2008年年會城市氣象與城市可持續發展分會場論文集.北京:中國氣象學會,2008.

Zhang J J,Dou H Y,Zhang E J,et al.The features study of the Pearl River Delta urban heat island[C]//Chinese Meteorological Society.The collected papers of urban meteorological and urban sustainable development of Chinese meteorological society annual conference 2008.Beijing:Chinese Meteorological Society,2008.

[2] 竇浩洋,張晶晶,趙昕奕.珠江三角洲城市熱島空間分布及熱島強度研究[J].地域研究與開發,2010,29(4):72-77.

Dou H Y,Zhao J J,Zhao X Y.Study on spatial distribution and intensity of urban heat island in Pearl River Delta[J].Areal Research and Development,2010,29(4):72-77.

[3] 曾 俠,錢光明,潘蔚娟.珠江三角洲都市群城市熱島效應初步研究[J].氣象,2004,30(10):12-16.

Zeng X,Qian G M,Pan W J.Study on urban heat island effect in Pearl River Delta urban group[J].Meteorological Monthly,2004,30(10):12-16.

[4] 歷 華,曾永年,贠培東,等.基于MODIS數據的長株潭地區城市熱島時空分析[J].測繪科學,2007,32(5):108-110,116.

Li H,Zeng Y N,Yun P D,et al.Temporal and spatial characteristics of urban heat island in Changsha-Zhuzhou-Xiangtan area based on MODIS data[J].Science of Surveying and Mapping,2007,32(5):108-110,116.

[5] 閆 峰,覃志豪,李茂松,等.基于MODIS數據的上海市熱島效應研究[J].武漢大學學報(信息科學版),2007,32(7):576-580.

Yan F,Qin Z H,Li M S,et al.On urban heat island of Shanghai City from MODIS data[J].Geomatics and Information Science of Wuhan University,2007,32(7):576-580.

[6] 楊 鵬,陳 靜,高 祺,等.基于MODIS數據的石家莊城市熱島效應研究[J].河北遙感,2012(2):18-21.

Yang P,Chen J,Gao Q,et al.The research of Shijiazhuang urban heat island effect based on the MODIS data[J].Journal of Hebei Remote Sensing,2012(2):18-21.

[7] Elvidge C D,Tuttle B T,Sutton P C,et al.Global distribution and density of constructed impervious surfaces[J].Sensors,2007,7(9):1962-1979.

[8] Schueler T R.The importance of imperviousness[J].Watershed Protection Techniques,1994,1(3):100-111.

[9] Arnold Jr C L,Gibbons C J.Impervious surface coverage:The emergence of a key environmental indicator[J].Journal of the American Planning Association,1996,62(2):243-258.

[10] Qin Z H,Dall’Olmo G,Karnieli A,et al.Derivation of split window algorithm and its sensitivity analysis for retrieving land surface temperature from NOAA-advanced very high resolution radiometer data[J].Journal of Geophysical Research:Atmospheres(1984-2012),2001,106(D19):22655-22670.

[11] Zhang P,Imhoff M L,Wolfe R E,et al.Characterizing urban heat islands of global settlements using MODIS and nighttime lights products[J].Canadian Journal of Remote Sensing,2010,36(3):185-196.

[12] 何全軍,曹 靜,黃 江,等.基于多光譜綜合的MODIS數據云檢測研究[J].國土資源遙感,2006,18(3):19-22.doi:10.6046/gtzyyg.2006.03.05.

He Q J,Cao J,Huang J,et al.Cloud detection in MODIS data based on multi-spectrum synthesis[J].Remote Sensing for Land and Resources,2006,18(3):19-22.doi:10.6046/gtzyyg.2006.03.05.

[13] 楊鐵利,何全軍.MODIS數據的云檢測處理[J].鞍山科技大學學報,2006,29(2):162-166.

Yang T L,He Q J.Cloud detection in MODIS data[J].Journal of Anshan University of Science and Technology,2006,29(2):162-166.

[14] Wan Z M,Dozier J.A generalized split-window algorithm for retrieving land-surface temperature from space[J].IEEE Transactions on Geoscience and Remote Sensing,1996,34(4):892-905.

[15] Becker F.The impact of spectral emissivity on the measurement of land surface temperature from a satellite[J].International Journal of Remote Sensing,1987,8(10):1509-1522.

[16] Mao K B,Qin Z,Shi J,et al.A practical split-window algorithm for retrieving land-surface temperature from MODIS data[J].International Journal of Remote Sensing,2005,26(15):3181-3204.

Analysisoftemporal-spatialvariationofheatislandeffectinPearlRiverDeltausingMODISimagesandimpermeablesurfacearea

HE Liqin1, YANG Peng1, JING Xin1, YAN Lei1, SU Linlin1,2

(1.BeijingKeyLabofSpatialInformationIntegrationand3SApplication,PekingUniversity,Beijing100087,China;2.InformationEngineeringCollege,CapitalNormalUniversity,Beijing100037,China)

The Pearl River Delta urban group is a region with rapid economic development; nevertheless, with the economic development, the heat island effect becomes increasingly obvious. Contrast with traditional point surveillance, the thermal infrared remote sensing method can make us understand the spatial distribution of urban heat island more intuitively. In this paper, the authors retrieved the land surface temperature(LST)of the four seasons of the Pearl River Delta region using the split window algorithm based on MODIS images, and divided the urban core and suburban region based on impermeable surface area(ISA), and finally obtained the surface urban heat island intensity. The results show that the Pearl River Delta region has a serious heat island phenomenon, with the most severe season being summer and the weakest season being winter. There is a tendency that the connection of cities has led to the formation of large urban heat island, especially in the two most serious cities, Foshan and Guangzhou. The heat island intensity is negatively correlated with NDVI value and positively correlated with the degree of the city’s economic development. The research results would provide some ecological instructions for urban development planning of the Pearl River Delta region.

heat island intensity; Pearl River Delta; MODIS; land surface temperature(LST); impermeable surface area(ISA); influence factor

10.6046/gtzyyg.2017.04.21

賀麗琴,楊鵬,景欣,等.基于MODIS影像及不透水面積的珠江三角洲熱島效應時空分析[J].國土資源遙感,2017,29(4):140-146.(He L Q,Yang P,Jing X,et al.Analysis of temporal-spatial variation of heat island effect in Pearl River Delta using MODIS images and impermeable surface area[J].Remote Sensing for Land and Resources,2017,29(4):140-146.)

TP 79

A

1001-070X(2017)04-0140-07

2016-03-23;

2016-04-25

國家自然科學基金項目“遙感云圖-電磁波-熱紅外對地震等重大地質災害的預測機理”(編號: 41371492)資助。

賀麗琴(1990-),女,碩士研究生,主要從事熱紅外和高光譜遙感方面的研究。Email: heliqin_de@126.com。

(責任編輯:李瑜)

猜你喜歡
效應區域
鈾對大型溞的急性毒性效應
永久基本農田集中區域“禁廢”
今日農業(2021年9期)2021-11-26 07:41:24
分割區域
懶馬效應
今日農業(2020年19期)2020-12-14 14:16:52
場景效應
應變效應及其應用
關于四色猜想
分區域
基于嚴重區域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
偶像效應
主站蜘蛛池模板: 青青草原偷拍视频| 久久久久亚洲精品成人网| 亚欧美国产综合| jijzzizz老师出水喷水喷出| 欧美黄网在线| 日本免费精品| 91精品情国产情侣高潮对白蜜| 久久精品无码一区二区日韩免费| 欧美成人免费一区在线播放| 免费在线色| 国产产在线精品亚洲aavv| 重口调教一区二区视频| 四虎免费视频网站| 欧美有码在线观看| 亚洲视频二| 五月天丁香婷婷综合久久| 亚洲美女视频一区| 免费观看男人免费桶女人视频| 久久永久免费人妻精品| 国产乱子伦一区二区=| 中文字幕在线观| 欧美色视频网站| 精品久久人人爽人人玩人人妻| 国产微拍一区| 欧美翘臀一区二区三区| 福利片91| 97久久人人超碰国产精品| 国产午夜一级淫片| 欧美a在线| 国产日本视频91| 丰满少妇αⅴ无码区| 久久综合九色综合97网| 日韩A∨精品日韩精品无码| 国产精品女在线观看| 国产区在线观看视频| 国产毛片久久国产| 亚洲va精品中文字幕| 美女免费黄网站| 国产成人喷潮在线观看| 日韩午夜片| 亚洲福利网址| 精品一區二區久久久久久久網站| 欧美在线视频a| 国产大片喷水在线在线视频| 精品国产免费观看| 国产第一色| 亚洲人成影院在线观看| 国产99精品久久| 午夜天堂视频| 91色在线观看| 91国内外精品自在线播放| 亚洲人成在线精品| 国内毛片视频| 欧美激情视频二区| 久久久亚洲国产美女国产盗摄| 亚洲欧美天堂网| 亚洲美女一区| a天堂视频| 草逼视频国产| 国产专区综合另类日韩一区| 国产91色| 国产精品99久久久| 精品国产成人国产在线| 亚洲人成网站18禁动漫无码| 中文字幕欧美日韩高清| 青青青草国产| 丰满少妇αⅴ无码区| 中文字幕乱码中文乱码51精品| 国产精品一老牛影视频| 爽爽影院十八禁在线观看| 国产亚卅精品无码| 亚洲日韩精品无码专区| 在线综合亚洲欧美网站| 91探花国产综合在线精品| 成人中文在线| 美女一级免费毛片| 国产一级妓女av网站| 精品偷拍一区二区| 亚洲成人黄色在线| 免费国产不卡午夜福在线观看| 99精品国产自在现线观看| 黄色国产在线|