京津冀城市群生态空间格局变化与地表温度关系研究

王 戈 于 强 YANG Di 赵晓婷 赵桂芳 岳德鹏

(1.北京林业大学精准林业北京市重点实验室， 北京 100083；2.蒙大拿大学空间分析实验室， 米苏拉 MT 59808)

0 引言

城市群是城市发展到成熟阶段的最高空间组织形式[1]。随着城市群内城市范围的扩张与城市数量的增加，建筑用地比例持续增加、生态用地不断减少，景观趋于破碎，产生了大量生态环境问题[2]，如热环境、水环境、城市生态服务功能下降等问题[3]。其中，城市群热环境问题尤为突出。热环境问题是绝大多数城市群所共有的问题[4]，由热环境效应影响区域所构成的“热空间”面积不断增加，破坏了大气环境，导致空气污染物不易扩散、城郊对流增强、单次降雨增加，易引发城市洪涝灾害[5]。

城市群生态空间可分为绿色生态空间与蓝色生态空间(简称绿色空间与蓝色空间)。绿色空间是指可以维持城市群生态环境稳定的所有植被出现的地方，可以为自然界提供生境，保护生物多样性[6]；城市群蓝色空间是指城市群中水体区域[7]。由蓝色空间与绿色空间构成的生态空间对城市群生态稳定意义重大，建设城市群蓝绿空间是缓解热环境问题的重要手段[8]。然而，我国人多地少，城市群绿地与水体的规划建设受到诸多方面制约[9]，生态空间的降温效果不仅与植被、水体的面积和空间分布有关，更受生态空间本身格局的影响[10]。地表温度作为有效表征热环境影响程度的指标而被广泛应用于城市群热环境问题研究中[11]。景观指数可以高度浓缩景观信息、反映其结构组成和空间配置。目前，生态空间与地表温度关系研究大多集中在绿地景观与地表温度的关系方面，有关学者基于遥感数据与气象观测数据的研究表明，城市植被绿地可通过蒸腾作用和阴影来降低地表温度，植被覆盖组成、NDVI与地表温度间存在显著相关性，并且植被比例越高、地表温度越低[12]。绿地景观空间结构对地表温度也具有显著影响，地表温度与区域内绿地景观斑块大小、形状指数、聚集度指数等景观指数间存在正相关关系，绿地格局参数对地表温度影响差异较大[13]。陈爱莲等[14]探究了北京市城市绿地格局与地表温度及冷岛强度的相关关系，程好好等[15]研究了深圳市绿地类型、属性特征与地表温度分异之间的关系，袁振等[16]将地温空间格局与城市绿地斑块进行了空间耦合分析。以上研究仅局限于单一研究区，缺乏对不同景观组成地区的比较，采用相关性分析与线性回归忽略了景观的空间异质性[17]。本研究采用空间自相关与空间自回归模型探究京津冀城市群内不同地区生态空间格局与地表温度关系[18]，为京津冀城市群热环境问题缓解提供科学依据。

1 材料与方法

1.1 研究区概况

京津冀城市群位于华北平原北部(北纬36°5′～42°40′，东经113°27′～119°50′)，北靠燕山山脉，南接华北平原，西倚太行山，东临渤海湾，西北和北部地形较高，南面和东面地形较为平坦，气候属于暖温带大陆性季风型气候，降水分配不均匀，春秋短促且干旱多风，夏季高温多雨，冬季寒冷干燥，分布较广的植被类型为草地、常绿针叶林和落叶阔叶林。京津冀城市群是我国三大城市群之一，主要包括北京市、天津市和河北省在内共13个城市(图1)，区域面积为 18.85万km2，人口数量1.18亿。2018年京津冀城市群GDP为8.5×104亿元,占全国GDP总量的9.44%，作为国家“十三五”规划重点发展城市群，对我国华北地区具有重要意义。近些年来，北京市、天津市和河北省发展迅速，京津冀城市群的下垫面改变使得地表温度发生变化，加上人口聚集密度高等因素，导致地表温度升高，城市群热环境状况进一步加剧。

1.2 数据来源

本文选用的地表温度数据是搭载在Terra卫星上MODIS陆地系列产品中的地表温度8 d合成产品(MOD11)。研究选用的2003—2018年树木覆盖率、短植被覆盖率与裸地覆比率由MOD44B植被覆盖转换产品(VCF)提供，VCF是指植被冠层或叶面积在地表的垂直投影面积占植被区总面积的比例，其中包括3种类型: 树木冠层(Tree canopy, TC)覆盖率、短植被(Short vegetation, SV)覆盖率与裸地(Bare ground, BG)覆比率，数值范围为0～200，0～100表示该区域的树木、短植被或裸地覆盖率，100～200表示该区域为水体。植被指数(NDVI)由MODIS陆地标准数据产品MOD13提供，MNDWI由MOD09数据中第4、6波段计算得到。MODIS数据的储存格式为HDF，地图投影为正弦投影，利用MODIS Reprojection Tool(MRT)工具进行批量几何校正与重采样处理，批处理后数据分辨率为1 000 m×1 000 m，可以满足京津冀城市群研究的需要。

1.3 研究方法

1.3.1MODIS地表温度处理

本文研究重点为探究生态空间格局与地表温度的耦合关系，夏季白天是高温对城市群造成影响最严重的时间段，故选取2018年夏季典型高温月(7月)的白天地表温度数据。选取8 d合成全球地表温度和发射率产品(MOD11)，对MOD11中Day_view_time(白天数据观测时间)子数据集进行预处理，然后对其进行辐射定标并转换为地表温度。地表温度计算公式为

Ts=0.02DN-273.15

(1)

式中Ts——地表温度,℃

DN——像元灰度

受云层影响，地表温度产品无法记录云层覆盖区域的地表温度数据，MOD11A2缺失值为0或低于正常值。运用月最大值合成法记录影像中同一坐标像元在合成范围内最大值来有效填补缺失数据，将所有数据进行重采样至分辨率为1 000 m×1 000 m。

1.3.2趋势分析法

为了能够掌握研究区2003—2018年的树木覆盖率(TC)、耕地覆盖率(SV)和裸地覆比率(BG)在时间与空间的变化特征，采用一元线性回归法模拟每个栅格的年际变化趋势，运用ArcGIS软件中空间分析模块来计算变化趋势，计算公式为

(2)

式中n——某一特征参数在研究时间段内累计时间,a

Xi——某一种特征参数第i年的平均值

Slope——趋势线斜率，为负表示该参数在整个时间段呈下降趋势，为正呈上升趋势

对每一个栅格进行趋势斜率分析，来研究整个区域各个特征参数的变化情况[19]。

1.3.3景观格局指数

景观格局可以反映景观的结构组成与空间配置，景观格局指数可以将景观格局定量化并且高度浓缩景观信息。根据文献[20]关于运用景观格局指数评价景观格局与热环境关系的研究，从类型水平与景观水平上共选取7个常用景观格局指数(表1)，利用Fragstats 4.2软件计算2018年各样区景观格局指数。

1.3.4移动窗口分析

景观格局指数之间存在一定关联性[21]，为全面反映生态空间的景观格局[18]，选取景观类型比例(PLAND)、斑块密度(PD)、最大斑块指数(LPI)、聚合度指数(AI)、平均斑块面积(AREA_MN)、面积加权平均斑块分维数(FRAC_AM)、形状指数(LSI)来研究2018年的生态空间与地表温度间的关系。结合实际情况选取Fragstats 4.2软件，采用边长为10 km的正方形窗口自研究区左上角移动，提取窗口内生态空间格局指数，最终得到生态空间景观指数的栅格图像[22]。通过ArcGIS 10.4软件，计算各个窗口内平均地表温度与景观指数并进行归一化处理。

1.3.5空间自相关分析

空间自相关分析可以解释景观的聚集特性[23]，分为全局空间自相关与局部空间自相关，两者通过Moran’sI来衡量空间自相关程度。ANSELIN提出双变量空间自相关分析来解释空间变量与邻近区域其他变量之间的相关性[24]。

运用GeoDa1.6软件进行空间自相关分析，全局自相关分析用于衡量样区内各个空间单元与临近单元之间的聚集性，计算公式为

表1 景观格局指数Tab.1 Landscape pattern index

(3)

式中n——空间单元总数

I——全局空间自相关值

wij——位置i、j观测值的权重

xi、xj——位置i、j的观测值

Wij——通过K邻接关系建立的空间权重矩阵

S2——所有空间单元中变量观测值x的方差

局部自相关分析用于分析变量在局部区域的空间相关性，基于z检验的LISA分布图能直观呈现局部区域变量的聚集性及分异特征，计算公式为

(4)

式中I——局部空间自相关值

1.3.6空间自回归模型

本研究通过GeoDa1.6软件采用空间自回归模型对2018年7月地表温度与生态空间景观指数进行回归分析，公式为

y=ρW1y+βx+μ+a

(5)

其中

μ=λW2μ+ε

式中y——因变量，即LST

x——自变量，即各类生态空间景观指数

β——自变量回归系数

a——截距μ——随机误差项

ε——服从均值为0、方差为δ2的随机误差[25]

W1、W2μ——因变量与残差的空间邻接权重矩阵

ρ——空间滞后项W1y的回归系数

λ——空间残差项回归系数

通过最大似然对数(Maximum likelihood logarithm, LIK)、赤池信息量准则(Akaike information criterion, AIC)、Schwartz指标(Schwartz criterion, SC)与回归模型误差项的Moran’sI值检验不同空间自回归模型的拟合效果[26]。模型LIK越大， AIC、SC越小，模型残差的Moran’sI值越接近于0，说明模型拟合效果越好[27]。

2 结果与分析

2.1 树木覆盖率、短植被覆盖率和裸地覆比率的变化趋势分析

年际变化空间趋势分布可反映参数趋势分布情况和研究区不同位置趋势分布情况。基于2003—2018年研究区VCF变化情况，用斜率表示近15年的变化幅度，正或负值表示VCF的增加或减少(图2)。从图2a可看出，在2003—2018年，京津冀地区的中部、东北以及西南边界的TC变化趋势值为正，其中大部分在0.5以上，表明该地区的TC呈明显的上升趋势，东南部有零散的正值区域，但数值较小，上升趋势不明显；京津冀地区的西北及东南大片区域TC的变化趋势值为负，且在-2～0之间，表明该部分地区TC呈下降趋势，但幅度较小，在承德市西北、东南以及天津东部沿海区域变化趋势值达-2以下，TC呈明显下降趋势。2003—2018年间，京津冀地区的东北、西南边界和东部沿海地区SV变化趋势值为正(图2b)，表明该地区的SV呈上升趋势，其中张家口市东部、承德市中部、唐山市南部地区变化趋势值达2以上，SV上升趋势明显，京津冀地区西北、中部和东部区域的SV变化趋势值为负，表明该地区SV呈下降趋势，其中张家口市、廊坊市、石家庄市中部区域变化趋势值达1.7以上，SV上升趋势明显。从图2c可以看出，2003—2018年间，承德市西北和南部、保定市、石家庄市、邢台市和邯郸市的大部分区域变化趋势值为正，表明BG呈上升趋势，其中邯郸市东南部上升趋势明显；张家口市、承德中部、北京市、廊坊市、秦皇岛市、唐山市、天津市、沧州市和衡水市的大部分区域变化趋势值为负，表明BG呈下降趋势。

年际变化的空间趋势可以在空间上清楚反映树木、短植被和裸地的变化趋势分布情况和不同位置的趋势及差异情况。基于研究区树木、短植被和裸地趋势分析，选取2018年的7个样区(图3)进行研究，样区1为张北地区，全域内最大面积景观为草地，建成区面积较小。样区2为京津地区，全域内建成区分布较为集中，草地分布在样区2西北部，水体主要分布在西南与东南部。样区3位于承德以南唐山以北，样区内分布大量林地景观，建成区主要分布在南部。样区4分布在邯郸市附近，地势平缓海拔较低，样区内分布着大量耕地。样区5位于邢台以南衡水以北，样区内耕地呈现增长趋势，并伴随土地裸露风险。样区6主要位于石家庄市，草地分布在样区西部，中部与东部分布大量建筑用地与耕地，样区7位于承德市北部，海拔较高，地势起伏较大，样区北部与南部分布大量林地景观，少量耕地分布在样区东南方向。

2.2 生态空间格局分析

基于2018年7月京津冀月平均NDVI与MNDWI数据，提取NDVI或MNDVI值大于0.85的区域作为2018年7月绿色空间或蓝色空间(图4)，绿色与蓝色空间共同组成了京津冀地区的生态空间。比较7个样区的生态空间指数，从表2可以看出，样区4生态空间最少，样区3最多，达到30.58%。样区7 LPI最大，为55.600%，样区6最小，为20.007%。样区2斑块密度最大，PD为0.169，样区5斑块密度最小，PD为0.068。样区7 AREA_MN最大，为4 090.6 m2，样区4 AREA_MN最小，为650 m2。样区2 LSI最大，景观复杂程度高，样区7 LSI最小。样区1～7的FRAC_AM差别不大，各样区内斑块形状对生态过程影像差别不大。样区7景观聚合度指数最大，样区内景观斑块团聚程度较高。

PLAND作为评价景观丰度的指标，可以反映生态空间分布情况。以生态空间的PLAND指数作为变量进行全局空间自相关分析，样区5生态空间景观的Moran’sI值为0.743(表2)，表明具有明显空间聚集性。

表2 生态空间格局指数全局自相关分析Tab.2 Global autocorrelation analysis of ecological spatial pattern index

2.3 地表温度格局分析

对京津冀城市群的地表温度数据进行空间自相关分析，结果如图5所示，样区1～7的全局空间自相关指数分别为0.475、0.737、0.607、0.376、0.476、0.735和0.638，且Moran’sI显著性P值均小于0.001，样区1～7的地表温度具有较强的空间正相关性。由图5可知，在样区1中西北与东南方向出现LST的高-高聚集区，样区2的西南部出现LST高温聚集区，样区3的北部与东部、样区4西北部、样区5中部、样区6西南部、样区7西北与东北部均出现LST的高-高聚集区。LST的低-低聚集区分布在样区1中部、样区2北部与东部、样区3的中部与西部、样区4的东部、样区5与样区6的北部、样区7的北部与南部。图5中移动窗口尺寸为10 km×10 km。

2.4 生态空间格局与地表温度耦合关系分析

2.4.1相关性与双变量空间自相关分析

Pearson 相关性及双变量空间自相关分析如表3所示，样区1 PLAND、PD、LSI、LPI、 FRAC_AM AREA_MN 与LST均呈负相关，AI与LST呈正相关；样区2 PD、LSI、LPI、FRAC_AM、AREA_MN、AI与LST均呈正相关；样区3 PLAND、LPI、AREA_MN、AI在空间上与LST呈正相关；样区4 PLAND、PD、LSI、LPI、FRAC_AM、AREA_MN、AI与LST均呈负相关；样区5中PLAND、LSI、LPI、FRAC_AM、AREA_MN、AI与LST呈现正相关，PD与LST呈负相关；样区6中PLAND、LPI、AREA_MN、AI与LST呈现正相关，PD、LSI、FRAC_AM与LST呈负相关；样区7中PLAND、LSI、LPI、FRAC_AM、AREA_MN、AI与LST均呈负相关，PD与LST呈正相关。7个样区的生态空间景观格局与LST相关系数以及双变量空间自相关Moran’sI差异明显。样区5、7位于河北省北部，林地景观比例较高，相关性以及双变量空间自相关性高于其他样区，与景观优势度、斑块破碎度有关。由于样区1、4生态空间比例较低，对LST影响有限。PLAND对LST影响较大，样区7生态空间比例高，生态空间斑块化且集中连片对地表温度影响明显。

表3 生态空间格局指数与地表温度相关性及双变量空间自相关分析Tab.3 Correlation between ecological spatial pattern index and surface temperature and bivariate spatial autocorrelation analysis

2.4.2空间自回归分析

生态空间景观指数与地表温度具有空间相关性，在此基础上以景观指数为自变量，地表温度为因变量分别采用3种模型进行空间回归分析，结果见表4，由于生态空间LST的影响受景观优势度、破碎度以及聚集程度等因素共同作用，通过对比可得多元回归拟合效果。对比样区1、7 3种回归模型的残差Moran’sI值如表5所示，OLS模型普遍高于SLM与SEM模型，样区6 OLS模型最高，数值为0.729，不可以有效解释变量间的空间关系，样区1～7的SLM与SEM模型拟合效果远优于OLS。大部分样区SEM模型R2大于SLM模型，SEM模型解释变量的能力更强，样区7 SEM模型R2最大，为0.834，样区5 SEM模型R2最小，为0.389。各个样区SEM模型的LIK值较大，AIC、SC数值和Moran’sI值较小，整体上模型SEM的拟合效果优于SLM。

3 结论

(1)基于2003—2018年研究区的MODIS遥感数据，京津冀地区的中部、东北以及西南边界地区的森林覆盖度呈现增长趋势，该地区TC变化趋势值为正，耕地覆盖率呈上升趋势；承德市西北和南部、保定市、石家庄市、邢台市和邯郸市的大部分区域BG变化趋势值为正，该地区具有土地裸露风险。

(2)根据研究区VCF变化情况，分析提取了7个样区，并提取了2018年7月的绿色空间或蓝色空间。各样区生态空间与地表温度的空间分布具有显著的空间自相关性。

表4 3种空间回归模型参数对比Tab.4 Comparison of parameters of three spatial regression models

表5 3种空间回归模型误差对比Tab.5 Error comparison of three spatial regression models

(3)7个样区的生态空间景观格局与LST相关系数以及双变量空间自相关Moran’sI差异明显。样区5、7位于河北省北部，林地景观比例较高，相关性以及双变量空间自相关性高于其他样区，这与景观优势度、斑块破碎度有关。样区1、4生态空间比例较低，对LST影响有限。PLAND对LST影响较大，样区7生态空间比例高，生态空间斑块集中连片对地表温度影响明显。

(4)地表温度受生态空间的景观优势度、破碎度以及聚集程度等因素共同影响，通过对比得出，多元回归拟合效果较好。对比样区1、7的残差Moran’sI值，样区6的普通线性回归模型的残差Moran’sI值为0.729，无法有效解释变量间的空间关系，样区1、7空间滞后模型与空间误差模型拟合效果远优于OLS模型。各个样区空间误差模型的R2大于空间滞后模型，空间误差模型解释变量的能力更强。各个样区空间误差模型的LIK值较大，AIC、SC以及模型残差Moran’sI值较小，空间误差模型拟合效果优于空间滞后模型。