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

基于Landsat影像的南京市城市擴張與熱島效應響應關系研究

2023-11-02 09:46:40李夢秋石含寧楊樹文
農業與技術 2023年20期
關鍵詞:大氣

李夢秋石含寧楊樹文

(1.蘭州交通大學測繪與地理信息學院,甘肅 蘭州 730070;2.地理國情監測技術應用國家地方聯合工程研究中心,甘肅 蘭州 730070;3.甘肅省地理國情監測工程實驗室,甘肅 蘭州 730070)

引言

城市熱島效應(Urban Heat Island),是在城市化的過程中產生的,是由人類生產生活和城市生態系統共同作用而造成的[1]。隨著城市的擴張和人口的增加,城市熱島效應也會逐漸加劇,最明顯的表征就是城市區域的表層平均氣溫相對于周邊的溫度明顯升高,并且隨著時間的推進,溫度呈現出逐漸升高的趨勢。中國自改革開放以來,各個地區的城市建設進入了高峰期,尤其是沿海地區,如珠江三角洲、長江三角洲、環渤海經濟圈率先進入了城市化的進程當中。

熱島效應表現為城市上空的溫度比城郊的溫度要高,溫度曲線表現為城市高、兩端低的峰谷特征,溫度曲線此時就像是隆起的一座小島,因此被形象地稱為熱島效應。城市熱島效應一年四季都會存在,但是在夏季尤為明顯[2]。夏季氣溫相對較高,太陽輻射較高,由于城市的覆被類型大多為不透水面而且植被的覆蓋相對于城郊較少,所以城市對于太陽輻射的削弱作用就比較差。因此,在相同的太陽輻射條件下,城郊的溫度要比城市建成區低。另外,由于夏季天然的高溫,加上人類活動消耗了大量的能源,產生大量的熱能,這些熱能匯聚到城市近地面上空,導致了城市溫度的升高。

在城市熱島研究方面,早期研究的方式主要是利用氣象站獲取的氣象數據進行溫度的統計和計算。由于氣象站分布數量較少且分布的不均勻,加上氣象數據的處理工作較為繁瑣,所以對于區域性的城市熱島研究較為困難。1972年,Rao等[3]發現,可以利用熱紅外遙感進行城市熱島效應的研究,從此,航空航天器搭載傳感器進行城市溫度檢測稱為主流方法。隨著遙感技術和GPS、GIS技術的發展和成熟,遙感監測熱島效應更具有科學性。由于遙感具有宏觀、連續、分辨率高的特點,相比較傳統的方法能夠更加快捷精確地獲取溫度信息,從多角度進行城市熱島的分析。

早期遙感監測城市熱島主要是利用NOAA和AVHRR數據,但是這2種數據的空間分辨率較低,不能精確區分析中小城市的城市熱島分布情況。而Landsat系列數據由于具有較高的空間分辨率,彌補了前2種數據低空間分辨率的缺陷,能檢測中小尺度的城市熱島。

現階段對于城市熱島效應的研究主要集中在城市熱島的影響因素和熱島強度以及分布情況。Wen等[4]使用歸一化插值植被指數來研究城市熱島,發現熱島溫度與NDVI的大小呈現負相關,Chen等[5]通過歸一化處理裸露指數(NDBaI)和歸一化植被指數、歸一化水體指數、歸一化建筑指數研究了廣東地區的城市化和城市熱島的關系,發現除了NDVI在負值范圍內與溫度呈正相關,正值范圍內呈負相關,其他指數與溫度之間存在一定的相關性。Daniela Arnds等[6]通過研究分析了漢堡市夜間平均熱島的時空變化,觀測了40個不同的氣象站點后發現一年中4—9月為溫度的最高時間,并且隨著氣候的變化,城市的高溫情況會加劇。I N Kuznetsova等[7]發現,莫斯科的熱狀況主要受到人類活動的影響,并且溫度比寒冷的郊區要高的原因是供暖系統熱量的釋放和損失。在沒有積雪覆蓋的情況下,城市表面的熱輻射增強了這些熱源的人為影響。

我國利用遙感技術進行熱島效應的監測開始于20世紀80年代,但是由于數據來源的問題,發展速度較慢,最近十幾年才進入了快速發展的階段。發展初期,大量學者對于成熟時熱島的研究主要集中在北京、上海、廣州、深圳這些城市化程度比較高的區域。隨著城市化的進程不斷加快,一些二三線城市的城市熱島研究也逐漸增多。雷金睿等[8]利用1989年、1999年、2007年和2015年的4期Landat5 TM/8OLI-TIRS數據,通過景觀生態學和統計分析的方法研究了海口市的城市熱島和景觀格局的時空變化關系。孫永等[9]利用基于多層城市冠層方案BEP的建筑物能量模式BEM方案的WRF模式進行模擬研究重慶市的城市熱島的特征、形成原因和局部地區的熱力環流對于城市熱島形成的影響。樊志宇等[10]利用輻射傳輸方程對武漢市2001年、2007年和2016年的Landsat系列影像數據進行LST反演,并利用MODIS數據進行精度驗證,發現輻射傳輸方程目前依然適用于溫度反演,并且武漢市的城市熱島有從單中心向多中心的變化趨勢。魏冠軍等[11]同樣使用多期Landsat數據研究了蘭州市1994—2015年的城市熱島的空間分布,通過熱島強度分級發現熱島的擴展方向和城鎮用地的延伸方向一致,且強熱島對于城市下墊面的變化敏感程度最高。尹柯柯等[12]通過研究蘭州市的城市熱島與土壤調節植被指數(MSAVI)、增強型水體指數(EWI)、土壤亮度指數(NDSI)、土壤濕度指數(NDMI)之間的關系,發現熱島強度與這些指數呈現正向的線性關系,與修改后的土壤調節植被指數呈負相關。

基于以上的研究內容,本文選取了南京市作為研究區,采用單窗算法對南京市的地表溫度進行反演,并結合遙感影像計算出相關的指數,通過數學分析的方法來研究南京市的地表溫度與各種指數之間存在的數學關系。

1 研究方法

地表溫度的區域分布在研究城市氣候變化、植被的生長于分布、環境監測方面具有十分重要的作用。在城市應用上,由于城市建成區的面積較大,通過地面實時的監測困難較大且工作強度高,因此利用遙感數據進行地表溫度反演,將大大減少溫度監測的工作量。遙感數據的多種尺度,可以適用于不同尺度的溫度研究[13]。

熱輻射在傳遞的過程中,需要穿過大氣層才能被傳感器接收,在穿過大氣層的過程中,部分會被大氣吸收而發生衰減;此外,大氣自身也會向外發射一定強度的輻射,這部分熱輻射會向上層移動而被傳感器接收。向下的熱輻射在到達地表后會被地表反射回大氣層,經過大氣層后衰減以后到達傳感器。此外,地表本身也不是一個完全黑體,會向外產生熱輻射,所以目前復雜的熱場是進行熱遙感的困難之處。

根據地表熱輻射傳到方程,覃志豪等[14]推導出簡單且精度高的溫度推演方法——單窗算法。該方法把大氣和地表對熱輻射的影響全部包括在了推演公式之中,所以不需要進行大氣校正就能直接推算出地表的溫度。由于這種方法在計算時只需要利用遙感數據的熱紅外波段,所以這種方法被稱為單窗算法。相比于之前的溫度反演方法,單窗算法僅僅需要地表比輻射率、大氣透過率和大氣等效溫度這3個參數。

地表溫度反演的技術流程如圖1所示。

圖1 技術流程圖

圖2 研究區縮略圖

圖4 2013—2016年城市溫度分級圖

圖5 城市熱島面積散點及線性回歸圖

圖6 NDVI與地表溫度波譜對比圖

圖7 NDVI與溫度線性回歸圖

圖8 NDBI與溫度線性回歸圖

圖9 NDWI與溫度線性回歸圖

1.1 地表亮溫的計算

在進行地表溫度反演之前,需要對數據進行預處理,包括輻射定標、大氣校正、幾何校正、圖像增強、影像的鑲嵌、裁剪等。而本文選取的是Googleearth的數據源,影像的空間分辨率為0.5m。

亮度溫度是傳感器在衛星運行的軌道高度上觀測到的熱輻射強度所對應的溫度,這個溫度包括了大氣和地球表面對熱輻射傳導的影響。在推算影像的地表亮度溫度時,需要以下2步:將原始影像的DN值轉換為對應的熱輻射亮度值;根據熱輻射強度推算出對應的亮度溫度。

求算熱輻射強度的公式:

L(λ)=Lmin(λ)+[Lmax(λ)-Lmin(λ)]Qdn/Qmax

(1)

式中,L(λ)為輻射強度,W·m-2·sr-1·μm-1;Qmax為最大DN值(255);Qdn為像元灰度值;Lmax(λ)和Lmin(λ)為最大和最小輻射強度值。在計算出L(λ)以后,利用L(λ)可以計算出像元的亮度溫度:

T6=K2/Ln(1+K1/L(λ))-273

(2)

式中,T6為熱紅外影像的像元亮度溫度,℃;K1和K2均為常量。Landsat8 TIRS(band10)的K1為774.89,K2為1321.08,band11的K1為480.89,K2為1201.14。

1.2 植被覆蓋度的計算

在計算植被覆蓋度之前,需要計算出歸一化植被指數(NDVI)。使用經過大氣校正后的Landsat8全色波段數據進行NDVI的計算:

NDVI=(B5-B4)/(B5+B4)

(3)

式中,B5和B4分別代表了全色波段的近紅外和紅波段。

在求得了NDVI的結果以后,對NDVI進行重分類。當NDVI≤0.05時,可以將這些像元視為完全裸土覆蓋表面;當NDVI≥0.7時,可以將這些像元視為完全植被覆蓋表面;當0.05≤NDVI≤0.7時,像元視為植被與裸土按一定比例存在,可利用植被覆蓋度公式來計算這些像元的植被覆蓋度:

Pv=(NDVI-NDVIS)/(NDVIV-NDVIS)2

(4)

式中,NDVIS和NDVIV分別為完全裸土和完全植被時像元的NDVI值;Pv為影像中植被占據像元的比例。在計算時,根據重分類時的閾值,將NDVIS取值為0.05,NDVI取值為0.7。因此在使用ENVI計算Pv時,按照0.05和0.7對植被覆蓋度進行計算:

(b1gt0.7)×1+(b1lt0.05)×0+(b1ge0.05andb1le0.7)×[(b1-0.05)/(0.7-0.05)2]

(5)

由于Landsat8的熱紅外波段空間分辨率為400m,所以需要對計算出來的NDVI結果進行重采樣到400m。在得到NDVI計算結果以后,對其進行重分類。

1.3 地表比輻射率的計算

在計算地表比輻射率之前,需要對遙感影像進行地表覆被類型分類。本文按照計算所需利用監督分類將地表類型分成了3大類5個小類,分別是自然表面(林地、耕地、裸地)、水體和城鎮。

對于不同的地表覆被類型,求算比輻射率的方式也不同。

自然表面:

ε=-0.0462Pv2+0.0613Pv+0.9626

(6)

城鎮:

ε=-0.06712Pv2+0.086Pv+0.9589

(7)

水體:

ε=0.995

(8)

1.4 大氣等效溫度估計

影響大氣等效溫度的因素主要有大氣剖面的氣溫分布和大氣的狀態,在不同的溫度或區域,由于氣候條件的不同,在估算時也采用不同的估算模型。

夏季平均大氣:

Ta=16.0110+0.92621T0(2.4a)

(9)

冬季平均大氣:

Ta=19.2704+0.91118T0(2.4b)

(10)

熱帶平均大氣:

Ta=17.9769+0.91715T0(2.4c)

(11)

式中,T0是衛星過境時近地面附近的溫度,式中所有的參數的單位均為開爾文(K)。

1.5 大氣透射率估計

地表熱輻射在大氣中傳輸的過程中,大氣透射率會對其產生重要的影響,因此在地表溫度反演的過程中,大氣透射率是不可或缺的參數之一。研究發現,大氣透射率的大小收到了大氣中水分含量的影響,而其他大氣條件對大氣透射率的影響微乎其微。所以大氣水分含量是計算大氣透射率的關鍵參數。

τ=0.974290-0.08007w(2.5a)

(12)

τ=1.031412-0.11536w(2.5b)

(13)

τ=0.982007-0.09611w(2.5c)

(14)

τ=1.053710-0.14142w(2.5d)

(15)

其中,式(12)和式(13)式分別計算大氣水分含量在0.4~1.6時夏季和冬季的大氣透射率;式(14)和式(15)則是分別計算當大氣水分含量在1.6~3時夏季和冬季的大氣透射率。式中,w值衛星過境時近地面的大氣水分含量,g·cm-2。

2 研究區概況

本文選擇以南京市主城區作為研究對象。南京位于中國東部,長江下游地區,是江蘇省的省會,位于N31°14′~32°37′,E118°22′~119°14′。南京市總面積為6587km2,其中建成區面積為1398.69km2。截至2018年12月,南京市下轄11個區,其中南京主城區為鼓樓區、秦淮區、玄武區、建鄴區、雨花臺區、浦口區和棲霞區。目前,南京市常住人口為843.62萬人,戶籍人口為696.94萬人,城鎮常住人口為695.99萬人,城鎮化率為82.5%,南京主城區已經達到了完全城鎮化的水平。因此,選擇南京進行城市熱島效應的研究具有代表意義。

3 結果與分析

根據2013—2016年的地表覆被類型圖,在南京主城區,建設用地占據了大部分區域,且呈現出逐年增加的趨勢。水體以長江、秦淮河、玄武湖等大型水體和其他中小型湖泊為主。裸地隨著城市的建設逐漸減少,根據裸地的變化顯示,2013年開始,南京市江北地區開始開發建設,原先只在長江大橋兩側的建筑開始沿江向西南和東北地區延伸,并逐漸向北部發展。到2015年,江北地區的城市格局基本形成。而江南區域,由于南京青奧會的舉辦,奧體周圍大量的裸地開始出現,隨著建設進程的推進,奧體周邊的建筑開始出現并逐漸加密,到2016年,基本形成以奧體為中心的建筑群落。

主城轄區的林地比較固定,主要分為3塊,分別為浦口區的老山國家森林公園、棲霞區的鐘山風景區和雨花臺的將軍山風景區、牛首山文化旅游區。另外,轄區內還有很多規模較小的林地,如仙林地區和長江大橋北側的長江觀音景區。總體來看,南京市的城市綠地覆蓋面積較大,分布均勻,建筑密度較高的區域均有綠地分布,是比較宜居的城市環境。

耕地主要分布在八卦洲、江心洲和棲霞區與鎮江相鄰的沿江地帶,這些地區的耕地面積相對固定沒有出現明顯的下降和升高。而江北地區由于城市的建設導致一些耕地被侵占破壞,耕地面積下降,并存在不斷減少的趨勢。可以預見的是,隨著城市的不斷擴張,位于城市邊緣的耕地將會進一步縮減,被建設用地替代。

3.1 地表溫度分析

從南京市2013—2016年的溫度反演圖中看出,無論在哪一年,水體的反演溫度始終最低。在南京主城區的水體中,玄武湖的溫度最低,貫穿城區的長江溫度與玄武湖相比,相近或者更高一些。由于成像時間均為冬季,農田尚處于休耕時期,地表植被覆蓋情況不如耕種期理想,為裸地和耕地的混合覆被類型,所以溫度會略高。而南京市建成區域,大多為不透水面,所以溫度最高。

3.2 覆被類型溫度

通過溫度與用地類型的疊加分析,按照耕地、林地、裸地、水體、建設用地這5種類型對4期的溫度數據進行提取,結果發現,檢核用地的溫度是5種類型中最高的,裸地次之,林地和耕地的溫度相近,水體最低。

3.3 熱島面積分析

根據統計出的各覆被類型的溫度,提出相交的溫度區域,通過設置閾值提取出建設用地中溫度較高的區域,將這些區域定義為城市高溫熱島區域。利用Arcgis空間分析計算異常值的面積,從而獲取不同年份的城市熱島面積大小。根據統計信息發現,城市高溫熱島面積呈現出逐年增長的趨勢,通過線性擬合發現,時間與熱島面積存在正相關的關系。

表1 5種用地類型溫度表

表2 2013—2016年覆被類型轉移矩陣

綜合幾年的土地利用動態變化影像和轉移矩陣分析可以得出,2013—2016年,由于老舊城區的拆遷和重建,使得裸地的面積沒有發生太大變化,耕地從230.9508km2減少到97.0407km2,林地從107.2692km2增長到297.8667km2,建設用地從517.6808km2減少到了464.7726km2,水體面積基本保持不變。結合2013年和2016年的土地利用分類影像,發現城市建設用地減少,可能的因素是2013年待開發的區域較多,且有許多地區已經開始了城市的基礎建設,而到了2016年,這些區域大多已經建設完成,且城市綠地規劃也已完成,所以從分類結果中可以看到,城市的面積相比2013年明顯減少,城市中的綠地面積明顯增加。

3.4 溫度與相關指數關系分析

3.4.1 溫度與NDVI關系

通過地表溫度與NDVI的結果對比發現,NDVI值高的地方,地表溫度低,而NDVI值低的地方,地表溫度較高。同時,從空間分布來看,NDVI值高的區域主要集中在鐘山風景區和城市建成區以外的地方,城市區域內的NDVI值相對較低。閾值相反的是,對于地表反演溫度而言,城市建成區的值要高,而城市邊緣區域的值較低。此外比較特殊的是水體,無論是NDVI值還是地表溫度值,都較低。

為了更加直觀的展示這一特點,將地表溫度和NDVI的波譜曲線進行對比,通過對比發現,兩者之間的特性跟圖上反映規律一致。為了進一步研究各個像元的NDVI值和地表溫度之間的數學關系,以2013年的NDVI和溫度反演結果為例,在研究區內選擇了較小的一塊區域,區域內包括了自然水體、建設用地、利林地、裸地、耕地等地表類型,將區域內的像元對應的NDVI和溫度值進行提取,進行進一步的回歸分析。

經過線性回歸分析發現,地表溫度與NDVI值之間存在明顯的線性關系,且R2為0.410,即當NDVI值越大時,地表溫度越低。這與先前的影像和波譜曲線的對比結果一致。但由于進行擬合的值中沒有提出水體這部分的異常結果,導致R2的值并不理想。于是將提取到的NDVI值與地表溫度數據進行二次篩選,提出其中的水體對應值,再次進行回歸分析。

經過二次回歸分析,R2有了明顯的提升,表明拋開水體的特殊情況后,NDVI值與地表溫度之間的確存在明顯的線性負相關關系。進一步驗證了對比結果。

3.4.2 地表溫度與歸一化建筑指數的關系

從上文的NDVI與溫度之間的波譜關系,粗略的發現建設用地與溫度之間存在正相關的關系,為了進一步證明兩者之間存在的具體相關性,在此引入了歸一化建筑指數(NDBI)進行建筑物的指數計算。

歸一化建筑指數的原理與NDVI相似,但使用的波段是中紅外波波段和近紅外波段,具體的求算公式:

NDBI=(MIR-NIR)/(MIR+NIR)

(16)

在獲取了NDBI的結果之后,同樣選取相同的區域,對溫度和NDBI結果進行分析。

分析發現,NDBI與溫度之間的確存在正相關的線性關系,擬合出的線性表達式:

y=2.92+19.18×x

(17)

式中,y表示地表溫度;x表示NDBI。由此表明,當建筑物的面積不斷增加時,城市的溫度將會出現明顯的上升,城市熱島效應的影響也會增大。

3.4.3 地表溫度與歸一化水體指數的關系

上文水體與反演處的溫度之間沒有明顯的正相關的關系,也沒有明顯的負相關的關系,為了從理論方面揭示水體與溫度之間的關系,引入改進型的歸一化水體指數:

NDWI=(b2-b4)/(b2+b4)

(18)

式中,b2為綠波段;b4為紅波段。

通過對相同區域的NDWI和溫度之間進行散點的擬合,發現兩者之間的確不存在明確的線性關系。

4 緩解熱島效應的措施

在城市熱島效應越發嚴重的今天,如何有效的預防和治理是必須認真考慮的問題。結合研究成果和分析,本文認為可以從以下方面入手。

城市的擴張和人口的增長是導致熱島效應加劇的重要因素之一,在城市的建設過程中,前期應當進行合理的規劃,可以將城市進行功能區的劃定。如,在規劃住宅區與附屬的功能區時,應將人口密度與功能區的完備情況、便捷程度進行綜合考慮。

在城市建設的進程中,自然表面的破壞是不可避免的,但是自然表面對于還曾會熱島效應又具有抑制作用,所以在城市建設中,應當有所取舍,有些會不可避免的破壞原有的自然表面,而有些區域可以避免對這些地表覆蓋的破壞。在人口密集度比較高的區域,應當建設城市人工綠地,綠地的光合作用將有效改善城市的空氣質量,也能一定程度降低城市近地面的溫度,對改善城市熱島效應效果明顯。

另外,人類活動是影響城市熱島效應的最為重要的因素。以本文研究區南京為例,在南京城區,尤其在長江沿岸和棲霞山地區,分布著許多造船廠和工礦碼頭,這些地區大多是24h滿負荷運轉,這些企業在給城市帶來了GDP的同時,也在消耗著大量的熱能。由于這些企業大多是老舊企業,生產方式相對比較落后,大多位高耗能產業,在消耗了大量的能源的同時,卻不能提高資源利用率,所以在生產的過程中向外排出了大量的未利用熱能,導致區域性的地表航空的溫度激增。目前,改善這一問題的方法可以是進行產能升級,關停高耗能企業,改善生產工藝,有效提高資源利用率,在生產的過程中嚴格落實環保措施,杜絕廢氣廢水的直接排放。

5 結論與討論

通過對南京市進行覆被類型提取和地表溫度反演,得到了南京市2013—2016年的覆被類型圖和土地利用轉移矩陣和地表溫度分級圖,通過對這些數據的研究得出以下結論。

南京主城區的地表覆被類型從2013—2016年變化明顯,其中建設用地的面積出現先增后減的情況,究其原因是老城區的拆遷改造。水體面積基本保持不變,以長江、秦淮河、玄武湖等水系為主。耕地主要集中于江心洲和棲霞區東北部,并且隨著時間的推進逐漸縮減。綠地由于城市建設的完善,面積不斷增加。

在提取的5種覆被類型中,建設用地和裸地的溫度始終較高。對于林地和耕地,由于兩者具有相似的波譜信息,所以反演出的溫度相近,而水體的溫度始終最低。異常高溫的面積隨著時間的推進而不斷增加,通過時間和面積的線性回歸分析發現,兩者之間存在顯著的正向線性關系。

通過溫度與NDVI、NDWI、NDBI的線性回歸分析,得出地表溫度與NDVI值呈現較為顯著的線性負相關關系,與NDBI值呈現顯著的線性正相關關系。而地表溫度和NDWI之間則沒有明顯的線性關系,存在隨機性。

經過地表溫度和各種指數之間的回歸分析,可以明確在城市的擴張過程中,城市熱島效應是不斷加劇的。建筑物的不斷增加,導致城市不透水面積的增加,無法削弱太陽輻射,進而導致了城市近地表溫度的上升。緩解城市熱島問題,從抑制城市的建設來說是不科學的,如何做到城市的擴張和熱島效應的不加劇是一個比較矛盾的問題。在有限的空間里盡可能充分利用去進行綠化建設,是比較切實可行的方法。

目前的城市熱島研究已經到了比較成熟的地步,結合地表溫度進行各種指數的計算也有很多學者進行研究。本文選取了較小的樣區進行溫度和各類指數的線性回歸分析,相較于權研究區的指數回歸分析,其說服力還不夠強。此外,通過反演出地表溫度,與景觀生態學的各種指標結合進行研究,給城市的發展建設作出指導性的意見也是后面需要進行研究的方向。

猜你喜歡
大氣
大氣的呵護
軍事文摘(2023年10期)2023-06-09 09:15:06
首次發現系外行星大氣中存在CO2
科學(2022年5期)2022-12-29 09:48:56
宏偉大氣,氣勢與細膩兼備 Vivid Audio Giya G3 S2
太赫茲大氣臨邊探測儀遙感中高層大氣風仿真
有“心氣”才大氣
如何“看清”大氣中的二氧化碳
學生天地(2020年18期)2020-08-25 09:29:24
大氣穩健的美式之風Polk Audio Signature系列
稚拙率真 圓融大氣
中國篆刻(2017年3期)2017-05-17 06:20:46
大氣古樸揮灑自如
大氣、水之后,土十條來了
新農業(2016年18期)2016-08-16 03:28:27
主站蜘蛛池模板: 亚洲国产精品国自产拍A| 国产不卡网| 又大又硬又爽免费视频| 99偷拍视频精品一区二区| 国产jizzjizz视频| 国产亚洲视频免费播放| 亚洲天堂视频网| 77777亚洲午夜久久多人| 五月婷婷中文字幕| 亚洲精品视频免费看| 久久6免费视频| 亚洲视频黄| 日韩亚洲高清一区二区| 久久综合国产乱子免费| 亚洲中文字幕久久精品无码一区| 国产污视频在线观看| 制服丝袜在线视频香蕉| 久久人人爽人人爽人人片aV东京热| 国产91蝌蚪窝| 久久性妇女精品免费| 欧美黄网在线| 国产欧美日韩综合在线第一| 国产乱子伦精品视频| 国产精品v欧美| 青青草原国产精品啪啪视频| 91福利国产成人精品导航| 免费毛片在线| 一本大道香蕉中文日本不卡高清二区 | 精品一区二区三区波多野结衣| 亚洲国产清纯| 超碰色了色| 最新国产网站| 国产三区二区| 国产丰满大乳无码免费播放| 欧美一级在线看| 中文字幕色在线| 免费可以看的无遮挡av无码| 国产激爽爽爽大片在线观看| 亚洲精品图区| 狼友av永久网站免费观看| 国产女人18毛片水真多1| av午夜福利一片免费看| 国产精品手机视频| 国产69精品久久久久妇女| 日韩毛片免费视频| 国产成人精品视频一区二区电影| 国产免费久久精品99re不卡| 97久久精品人人| 亚洲愉拍一区二区精品| 日本伊人色综合网| 国产青青操| 亚洲日本在线免费观看| 特级做a爰片毛片免费69| 国产一区二区三区在线观看免费| 99久久免费精品特色大片| 欧美激情福利| 不卡国产视频第一页| 亚洲天堂在线免费| 国产乱论视频| 国内精品免费| 欧美精品另类| 日本免费福利视频| 精品国产亚洲人成在线| 在线播放91| 亚洲中文字幕在线观看| 中文字幕在线观| a天堂视频| 国产一级二级三级毛片| 国产91av在线| 51国产偷自视频区视频手机观看| 亚洲中文字幕手机在线第一页| av免费在线观看美女叉开腿| 精品国产自在现线看久久| 毛片国产精品完整版| 青青极品在线| 国产精品无码AV中文| 国产精品丝袜视频| 香蕉网久久| 91最新精品视频发布页| 国产视频大全| 亚洲,国产,日韩,综合一区| 88av在线看|