京津冀城市群生態(tài)空間格局變化與地表溫度關(guān)系研究

王 戈 于 強 YANG Di 趙曉婷 趙桂芳 岳德鵬

(1.北京林業(yè)大學精準林業(yè)北京市重點實驗室， 北京 100083；2.蒙大拿大學空間分析實驗室， 米蘇拉 MT 59808)

0 引言

城市群是城市發(fā)展到成熟階段的最高空間組織形式[1]。隨著城市群內(nèi)城市范圍的擴張與城市數(shù)量的增加，建筑用地比例持續(xù)增加、生態(tài)用地不斷減少，景觀趨于破碎，產(chǎn)生了大量生態(tài)環(huán)境問題[2]，如熱環(huán)境、水環(huán)境、城市生態(tài)服務功能下降等問題[3]。其中，城市群熱環(huán)境問題尤為突出。熱環(huán)境問題是絕大多數(shù)城市群所共有的問題[4]，由熱環(huán)境效應影響區(qū)域所構(gòu)成的“熱空間”面積不斷增加，破壞了大氣環(huán)境，導致空氣污染物不易擴散、城郊對流增強、單次降雨增加，易引發(fā)城市洪澇災害[5]。

城市群生態(tài)空間可分為綠色生態(tài)空間與藍色生態(tài)空間(簡稱綠色空間與藍色空間)。綠色空間是指可以維持城市群生態(tài)環(huán)境穩(wěn)定的所有植被出現(xiàn)的地方，可以為自然界提供生境，保護生物多樣性[6]；城市群藍色空間是指城市群中水體區(qū)域[7]。由藍色空間與綠色空間構(gòu)成的生態(tài)空間對城市群生態(tài)穩(wěn)定意義重大，建設(shè)城市群藍綠空間是緩解熱環(huán)境問題的重要手段[8]。然而，我國人多地少，城市群綠地與水體的規(guī)劃建設(shè)受到諸多方面制約[9]，生態(tài)空間的降溫效果不僅與植被、水體的面積和空間分布有關(guān)，更受生態(tài)空間本身格局的影響[10]。地表溫度作為有效表征熱環(huán)境影響程度的指標而被廣泛應用于城市群熱環(huán)境問題研究中[11]。景觀指數(shù)可以高度濃縮景觀信息、反映其結(jié)構(gòu)組成和空間配置。目前，生態(tài)空間與地表溫度關(guān)系研究大多集中在綠地景觀與地表溫度的關(guān)系方面，有關(guān)學者基于遙感數(shù)據(jù)與氣象觀測數(shù)據(jù)的研究表明，城市植被綠地可通過蒸騰作用和陰影來降低地表溫度，植被覆蓋組成、NDVI與地表溫度間存在顯著相關(guān)性，并且植被比例越高、地表溫度越低[12]。綠地景觀空間結(jié)構(gòu)對地表溫度也具有顯著影響，地表溫度與區(qū)域內(nèi)綠地景觀斑塊大小、形狀指數(shù)、聚集度指數(shù)等景觀指數(shù)間存在正相關(guān)關(guān)系，綠地格局參數(shù)對地表溫度影響差異較大[13]。陳愛蓮等[14]探究了北京市城市綠地格局與地表溫度及冷島強度的相關(guān)關(guān)系，程好好等[15]研究了深圳市綠地類型、屬性特征與地表溫度分異之間的關(guān)系，袁振等[16]將地溫空間格局與城市綠地斑塊進行了空間耦合分析。以上研究僅局限于單一研究區(qū)，缺乏對不同景觀組成地區(qū)的比較，采用相關(guān)性分析與線性回歸忽略了景觀的空間異質(zhì)性[17]。本研究采用空間自相關(guān)與空間自回歸模型探究京津冀城市群內(nèi)不同地區(qū)生態(tài)空間格局與地表溫度關(guān)系[18]，為京津冀城市群熱環(huán)境問題緩解提供科學依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

京津冀城市群位于華北平原北部(北緯36°5′～42°40′，東經(jīng)113°27′～119°50′)，北靠燕山山脈，南接華北平原，西倚太行山，東臨渤海灣，西北和北部地形較高，南面和東面地形較為平坦，氣候?qū)儆谂瘻貛Т箨懶约撅L型氣候，降水分配不均勻，春秋短促且干旱多風，夏季高溫多雨，冬季寒冷干燥，分布較廣的植被類型為草地、常綠針葉林和落葉闊葉林。京津冀城市群是我國三大城市群之一，主要包括北京市、天津市和河北省在內(nèi)共13個城市(圖1)，區(qū)域面積為 18.85萬km2，人口數(shù)量1.18億。2018年京津冀城市群GDP為8.5×104億元,占全國GDP總量的9.44%，作為國家“十三五”規(guī)劃重點發(fā)展城市群，對我國華北地區(qū)具有重要意義。近些年來，北京市、天津市和河北省發(fā)展迅速，京津冀城市群的下墊面改變使得地表溫度發(fā)生變化，加上人口聚集密度高等因素，導致地表溫度升高，城市群熱環(huán)境狀況進一步加劇。

1.2 數(shù)據(jù)來源

本文選用的地表溫度數(shù)據(jù)是搭載在Terra衛(wèi)星上MODIS陸地系列產(chǎn)品中的地表溫度8 d合成產(chǎn)品(MOD11)。研究選用的2003—2018年樹木覆蓋率、短植被覆蓋率與裸地覆比率由MOD44B植被覆蓋轉(zhuǎn)換產(chǎn)品(VCF)提供，VCF是指植被冠層或葉面積在地表的垂直投影面積占植被區(qū)總面積的比例，其中包括3種類型: 樹木冠層(Tree canopy, TC)覆蓋率、短植被(Short vegetation, SV)覆蓋率與裸地(Bare ground, BG)覆比率，數(shù)值范圍為0～200，0～100表示該區(qū)域的樹木、短植被或裸地覆蓋率，100～200表示該區(qū)域為水體。植被指數(shù)(NDVI)由MODIS陸地標準數(shù)據(jù)產(chǎn)品MOD13提供，MNDWI由MOD09數(shù)據(jù)中第4、6波段計算得到。MODIS數(shù)據(jù)的儲存格式為HDF，地圖投影為正弦投影，利用MODIS Reprojection Tool(MRT)工具進行批量幾何校正與重采樣處理，批處理后數(shù)據(jù)分辨率為1 000 m×1 000 m，可以滿足京津冀城市群研究的需要。

1.3 研究方法

1.3.1MODIS地表溫度處理

本文研究重點為探究生態(tài)空間格局與地表溫度的耦合關(guān)系，夏季白天是高溫對城市群造成影響最嚴重的時間段，故選取2018年夏季典型高溫月(7月)的白天地表溫度數(shù)據(jù)。選取8 d合成全球地表溫度和發(fā)射率產(chǎn)品(MOD11)，對MOD11中Day_view_time(白天數(shù)據(jù)觀測時間)子數(shù)據(jù)集進行預處理，然后對其進行輻射定標并轉(zhuǎn)換為地表溫度。地表溫度計算公式為

Ts=0.02DN-273.15

(1)

式中Ts——地表溫度,℃

DN——像元灰度

受云層影響，地表溫度產(chǎn)品無法記錄云層覆蓋區(qū)域的地表溫度數(shù)據(jù)，MOD11A2缺失值為0或低于正常值。運用月最大值合成法記錄影像中同一坐標像元在合成范圍內(nèi)最大值來有效填補缺失數(shù)據(jù)，將所有數(shù)據(jù)進行重采樣至分辨率為1 000 m×1 000 m。

1.3.2趨勢分析法

為了能夠掌握研究區(qū)2003—2018年的樹木覆蓋率(TC)、耕地覆蓋率(SV)和裸地覆比率(BG)在時間與空間的變化特征，采用一元線性回歸法模擬每個柵格的年際變化趨勢，運用ArcGIS軟件中空間分析模塊來計算變化趨勢，計算公式為

(2)

式中n——某一特征參數(shù)在研究時間段內(nèi)累計時間,a

Xi——某一種特征參數(shù)第i年的平均值

Slope——趨勢線斜率，為負表示該參數(shù)在整個時間段呈下降趨勢，為正呈上升趨勢

對每一個柵格進行趨勢斜率分析，來研究整個區(qū)域各個特征參數(shù)的變化情況[19]。

1.3.3景觀格局指數(shù)

景觀格局可以反映景觀的結(jié)構(gòu)組成與空間配置，景觀格局指數(shù)可以將景觀格局定量化并且高度濃縮景觀信息。根據(jù)文獻[20]關(guān)于運用景觀格局指數(shù)評價景觀格局與熱環(huán)境關(guān)系的研究，從類型水平與景觀水平上共選取7個常用景觀格局指數(shù)(表1)，利用Fragstats 4.2軟件計算2018年各樣區(qū)景觀格局指數(shù)。

1.3.4移動窗口分析

景觀格局指數(shù)之間存在一定關(guān)聯(lián)性[21]，為全面反映生態(tài)空間的景觀格局[18]，選取景觀類型比例(PLAND)、斑塊密度(PD)、最大斑塊指數(shù)(LPI)、聚合度指數(shù)(AI)、平均斑塊面積(AREA_MN)、面積加權(quán)平均斑塊分維數(shù)(FRAC_AM)、形狀指數(shù)(LSI)來研究2018年的生態(tài)空間與地表溫度間的關(guān)系。結(jié)合實際情況選取Fragstats 4.2軟件，采用邊長為10 km的正方形窗口自研究區(qū)左上角移動，提取窗口內(nèi)生態(tài)空間格局指數(shù)，最終得到生態(tài)空間景觀指數(shù)的柵格圖像[22]。通過ArcGIS 10.4軟件，計算各個窗口內(nèi)平均地表溫度與景觀指數(shù)并進行歸一化處理。

1.3.5空間自相關(guān)分析

空間自相關(guān)分析可以解釋景觀的聚集特性[23]，分為全局空間自相關(guān)與局部空間自相關(guān)，兩者通過Moran’sI來衡量空間自相關(guān)程度。ANSELIN提出雙變量空間自相關(guān)分析來解釋空間變量與鄰近區(qū)域其他變量之間的相關(guān)性[24]。

運用GeoDa1.6軟件進行空間自相關(guān)分析，全局自相關(guān)分析用于衡量樣區(qū)內(nèi)各個空間單元與臨近單元之間的聚集性，計算公式為

表1 景觀格局指數(shù)Tab.1 Landscape pattern index

(3)

式中n——空間單元總數(shù)

I——全局空間自相關(guān)值

wij——位置i、j觀測值的權(quán)重

xi、xj——位置i、j的觀測值

Wij——通過K鄰接關(guān)系建立的空間權(quán)重矩陣

S2——所有空間單元中變量觀測值x的方差

局部自相關(guān)分析用于分析變量在局部區(qū)域的空間相關(guān)性，基于z檢驗的LISA分布圖能直觀呈現(xiàn)局部區(qū)域變量的聚集性及分異特征，計算公式為

(4)

式中I——局部空間自相關(guān)值

1.3.6空間自回歸模型

本研究通過GeoDa1.6軟件采用空間自回歸模型對2018年7月地表溫度與生態(tài)空間景觀指數(shù)進行回歸分析，公式為

y=ρW1y+βx+μ+a

(5)

其中

μ=λW2μ+ε

式中y——因變量，即LST

x——自變量，即各類生態(tài)空間景觀指數(shù)

β——自變量回歸系數(shù)

a——截距μ——隨機誤差項

ε——服從均值為0、方差為δ2的隨機誤差[25]

W1、W2μ——因變量與殘差的空間鄰接權(quán)重矩陣

ρ——空間滯后項W1y的回歸系數(shù)

λ——空間殘差項回歸系數(shù)

通過最大似然對數(shù)(Maximum likelihood logarithm, LIK)、赤池信息量準則(Akaike information criterion, AIC)、Schwartz指標(Schwartz criterion, SC)與回歸模型誤差項的Moran’sI值檢驗不同空間自回歸模型的擬合效果[26]。模型LIK越大， AIC、SC越小，模型殘差的Moran’sI值越接近于0，說明模型擬合效果越好[27]。

2 結(jié)果與分析

2.1 樹木覆蓋率、短植被覆蓋率和裸地覆比率的變化趨勢分析

年際變化空間趨勢分布可反映參數(shù)趨勢分布情況和研究區(qū)不同位置趨勢分布情況。基于2003—2018年研究區(qū)VCF變化情況，用斜率表示近15年的變化幅度，正或負值表示VCF的增加或減少(圖2)。從圖2a可看出，在2003—2018年，京津冀地區(qū)的中部、東北以及西南邊界的TC變化趨勢值為正，其中大部分在0.5以上，表明該地區(qū)的TC呈明顯的上升趨勢，東南部有零散的正值區(qū)域，但數(shù)值較小，上升趨勢不明顯；京津冀地區(qū)的西北及東南大片區(qū)域TC的變化趨勢值為負，且在-2～0之間，表明該部分地區(qū)TC呈下降趨勢，但幅度較小，在承德市西北、東南以及天津東部沿海區(qū)域變化趨勢值達-2以下，TC呈明顯下降趨勢。2003—2018年間，京津冀地區(qū)的東北、西南邊界和東部沿海地區(qū)SV變化趨勢值為正(圖2b)，表明該地區(qū)的SV呈上升趨勢，其中張家口市東部、承德市中部、唐山市南部地區(qū)變化趨勢值達2以上，SV上升趨勢明顯，京津冀地區(qū)西北、中部和東部區(qū)域的SV變化趨勢值為負，表明該地區(qū)SV呈下降趨勢，其中張家口市、廊坊市、石家莊市中部區(qū)域變化趨勢值達1.7以上，SV上升趨勢明顯。從圖2c可以看出，2003—2018年間，承德市西北和南部、保定市、石家莊市、邢臺市和邯鄲市的大部分區(qū)域變化趨勢值為正，表明BG呈上升趨勢，其中邯鄲市東南部上升趨勢明顯；張家口市、承德中部、北京市、廊坊市、秦皇島市、唐山市、天津市、滄州市和衡水市的大部分區(qū)域變化趨勢值為負，表明BG呈下降趨勢。

年際變化的空間趨勢可以在空間上清楚反映樹木、短植被和裸地的變化趨勢分布情況和不同位置的趨勢及差異情況。基于研究區(qū)樹木、短植被和裸地趨勢分析，選取2018年的7個樣區(qū)(圖3)進行研究，樣區(qū)1為張北地區(qū)，全域內(nèi)最大面積景觀為草地，建成區(qū)面積較小。樣區(qū)2為京津地區(qū)，全域內(nèi)建成區(qū)分布較為集中，草地分布在樣區(qū)2西北部，水體主要分布在西南與東南部。樣區(qū)3位于承德以南唐山以北，樣區(qū)內(nèi)分布大量林地景觀，建成區(qū)主要分布在南部。樣區(qū)4分布在邯鄲市附近，地勢平緩海拔較低，樣區(qū)內(nèi)分布著大量耕地。樣區(qū)5位于邢臺以南衡水以北，樣區(qū)內(nèi)耕地呈現(xiàn)增長趨勢，并伴隨土地裸露風險。樣區(qū)6主要位于石家莊市，草地分布在樣區(qū)西部，中部與東部分布大量建筑用地與耕地，樣區(qū)7位于承德市北部，海拔較高，地勢起伏較大，樣區(qū)北部與南部分布大量林地景觀，少量耕地分布在樣區(qū)東南方向。

2.2 生態(tài)空間格局分析

基于2018年7月京津冀月平均NDVI與MNDWI數(shù)據(jù)，提取NDVI或MNDVI值大于0.85的區(qū)域作為2018年7月綠色空間或藍色空間(圖4)，綠色與藍色空間共同組成了京津冀地區(qū)的生態(tài)空間。比較7個樣區(qū)的生態(tài)空間指數(shù)，從表2可以看出，樣區(qū)4生態(tài)空間最少，樣區(qū)3最多，達到30.58%。樣區(qū)7 LPI最大，為55.600%，樣區(qū)6最小，為20.007%。樣區(qū)2斑塊密度最大，PD為0.169，樣區(qū)5斑塊密度最小，PD為0.068。樣區(qū)7 AREA_MN最大，為4 090.6 m2，樣區(qū)4 AREA_MN最小，為650 m2。樣區(qū)2 LSI最大，景觀復雜程度高，樣區(qū)7 LSI最小。樣區(qū)1～7的FRAC_AM差別不大，各樣區(qū)內(nèi)斑塊形狀對生態(tài)過程影像差別不大。樣區(qū)7景觀聚合度指數(shù)最大，樣區(qū)內(nèi)景觀斑塊團聚程度較高。

PLAND作為評價景觀豐度的指標，可以反映生態(tài)空間分布情況。以生態(tài)空間的PLAND指數(shù)作為變量進行全局空間自相關(guān)分析，樣區(qū)5生態(tài)空間景觀的Moran’sI值為0.743(表2)，表明具有明顯空間聚集性。

表2 生態(tài)空間格局指數(shù)全局自相關(guān)分析Tab.2 Global autocorrelation analysis of ecological spatial pattern index

2.3 地表溫度格局分析

對京津冀城市群的地表溫度數(shù)據(jù)進行空間自相關(guān)分析，結(jié)果如圖5所示，樣區(qū)1～7的全局空間自相關(guān)指數(shù)分別為0.475、0.737、0.607、0.376、0.476、0.735和0.638，且Moran’sI顯著性P值均小于0.001，樣區(qū)1～7的地表溫度具有較強的空間正相關(guān)性。由圖5可知，在樣區(qū)1中西北與東南方向出現(xiàn)LST的高-高聚集區(qū)，樣區(qū)2的西南部出現(xiàn)LST高溫聚集區(qū)，樣區(qū)3的北部與東部、樣區(qū)4西北部、樣區(qū)5中部、樣區(qū)6西南部、樣區(qū)7西北與東北部均出現(xiàn)LST的高-高聚集區(qū)。LST的低-低聚集區(qū)分布在樣區(qū)1中部、樣區(qū)2北部與東部、樣區(qū)3的中部與西部、樣區(qū)4的東部、樣區(qū)5與樣區(qū)6的北部、樣區(qū)7的北部與南部。圖5中移動窗口尺寸為10 km×10 km。

2.4 生態(tài)空間格局與地表溫度耦合關(guān)系分析

2.4.1相關(guān)性與雙變量空間自相關(guān)分析

Pearson 相關(guān)性及雙變量空間自相關(guān)分析如表3所示，樣區(qū)1 PLAND、PD、LSI、LPI、 FRAC_AM AREA_MN 與LST均呈負相關(guān)，AI與LST呈正相關(guān)；樣區(qū)2 PD、LSI、LPI、FRAC_AM、AREA_MN、AI與LST均呈正相關(guān)；樣區(qū)3 PLAND、LPI、AREA_MN、AI在空間上與LST呈正相關(guān)；樣區(qū)4 PLAND、PD、LSI、LPI、FRAC_AM、AREA_MN、AI與LST均呈負相關(guān)；樣區(qū)5中PLAND、LSI、LPI、FRAC_AM、AREA_MN、AI與LST呈現(xiàn)正相關(guān)，PD與LST呈負相關(guān)；樣區(qū)6中PLAND、LPI、AREA_MN、AI與LST呈現(xiàn)正相關(guān)，PD、LSI、FRAC_AM與LST呈負相關(guān)；樣區(qū)7中PLAND、LSI、LPI、FRAC_AM、AREA_MN、AI與LST均呈負相關(guān)，PD與LST呈正相關(guān)。7個樣區(qū)的生態(tài)空間景觀格局與LST相關(guān)系數(shù)以及雙變量空間自相關(guān)Moran’sI差異明顯。樣區(qū)5、7位于河北省北部，林地景觀比例較高，相關(guān)性以及雙變量空間自相關(guān)性高于其他樣區(qū)，與景觀優(yōu)勢度、斑塊破碎度有關(guān)。由于樣區(qū)1、4生態(tài)空間比例較低，對LST影響有限。PLAND對LST影響較大，樣區(qū)7生態(tài)空間比例高，生態(tài)空間斑塊化且集中連片對地表溫度影響明顯。

表3 生態(tài)空間格局指數(shù)與地表溫度相關(guān)性及雙變量空間自相關(guān)分析Tab.3 Correlation between ecological spatial pattern index and surface temperature and bivariate spatial autocorrelation analysis

2.4.2空間自回歸分析

生態(tài)空間景觀指數(shù)與地表溫度具有空間相關(guān)性，在此基礎(chǔ)上以景觀指數(shù)為自變量，地表溫度為因變量分別采用3種模型進行空間回歸分析，結(jié)果見表4，由于生態(tài)空間LST的影響受景觀優(yōu)勢度、破碎度以及聚集程度等因素共同作用，通過對比可得多元回歸擬合效果。對比樣區(qū)1、7 3種回歸模型的殘差Moran’sI值如表5所示，OLS模型普遍高于SLM與SEM模型，樣區(qū)6 OLS模型最高，數(shù)值為0.729，不可以有效解釋變量間的空間關(guān)系，樣區(qū)1～7的SLM與SEM模型擬合效果遠優(yōu)于OLS。大部分樣區(qū)SEM模型R2大于SLM模型，SEM模型解釋變量的能力更強，樣區(qū)7 SEM模型R2最大，為0.834，樣區(qū)5 SEM模型R2最小，為0.389。各個樣區(qū)SEM模型的LIK值較大，AIC、SC數(shù)值和Moran’sI值較小，整體上模型SEM的擬合效果優(yōu)于SLM。

3 結(jié)論

(1)基于2003—2018年研究區(qū)的MODIS遙感數(shù)據(jù)，京津冀地區(qū)的中部、東北以及西南邊界地區(qū)的森林覆蓋度呈現(xiàn)增長趨勢，該地區(qū)TC變化趨勢值為正，耕地覆蓋率呈上升趨勢；承德市西北和南部、保定市、石家莊市、邢臺市和邯鄲市的大部分區(qū)域BG變化趨勢值為正，該地區(qū)具有土地裸露風險。

(2)根據(jù)研究區(qū)VCF變化情況，分析提取了7個樣區(qū)，并提取了2018年7月的綠色空間或藍色空間。各樣區(qū)生態(tài)空間與地表溫度的空間分布具有顯著的空間自相關(guān)性。

表4 3種空間回歸模型參數(shù)對比Tab.4 Comparison of parameters of three spatial regression models

表5 3種空間回歸模型誤差對比Tab.5 Error comparison of three spatial regression models

(3)7個樣區(qū)的生態(tài)空間景觀格局與LST相關(guān)系數(shù)以及雙變量空間自相關(guān)Moran’sI差異明顯。樣區(qū)5、7位于河北省北部，林地景觀比例較高，相關(guān)性以及雙變量空間自相關(guān)性高于其他樣區(qū)，這與景觀優(yōu)勢度、斑塊破碎度有關(guān)。樣區(qū)1、4生態(tài)空間比例較低，對LST影響有限。PLAND對LST影響較大，樣區(qū)7生態(tài)空間比例高，生態(tài)空間斑塊集中連片對地表溫度影響明顯。

(4)地表溫度受生態(tài)空間的景觀優(yōu)勢度、破碎度以及聚集程度等因素共同影響，通過對比得出，多元回歸擬合效果較好。對比樣區(qū)1、7的殘差Moran’sI值，樣區(qū)6的普通線性回歸模型的殘差Moran’sI值為0.729，無法有效解釋變量間的空間關(guān)系，樣區(qū)1、7空間滯后模型與空間誤差模型擬合效果遠優(yōu)于OLS模型。各個樣區(qū)空間誤差模型的R2大于空間滯后模型，空間誤差模型解釋變量的能力更強。各個樣區(qū)空間誤差模型的LIK值較大，AIC、SC以及模型殘差Moran’sI值較小，空間誤差模型擬合效果優(yōu)于空間滯后模型。