基于Fragstats4的景观格局指数与地表温度的相关性——以石家庄市为例

杨鹏 高祺* 张艳品 刘思廷 齐晓华 曹春莉 程巳阳

(1 石家庄市气象局，石家庄 050081; 2 中国气象科学研究院灾害天气国家重点实验室，北京 100081)

引言

随着城市化的发展，土地利用/覆盖和景观格局都发生变化，水泥、沥青等不透水地表组成的建成区持续扩张，水体、天然植被等构建的生态空间不断萎缩，交通及路网用地密度提升，景观格局完整性遭到破坏[1-5]。近年来，从景观生态学的角度开展土地利用尺度效应的研究是一个逐渐兴起的研究方向。土地利用景观格局变化研究对优化城市布局和加快美丽城市建设至关重要，景观的组成和结构特征是影响热场分布的主要因素，组成特征对热场的影响主要是由于各地类热效应不同，结构特征主要影响斑块间能量交换[6-8]。

国内外许多学者从不同方面、利用不同方法开展景观格局指数与地表温度的相关性研究。Callo等[9]最先利用NOAA/AVHRR数据证明NDVI与LST有负相关；Weng Qihao等[10]利用景观格局指数研究不同土地利用景观格局对热环境的影响；Coseo等[11]通过对芝加哥城市热岛影响因素分析得出，不透水地表比率和数目可以解释68%的气温变化；江学顶[12]采用遥感和MM5模式研究珠江三角洲城市群的城市热力景观及其日变化规律；马勇刚等[13]以乌鲁木齐为研究区，将景观分类和景观格局变化图与地表温度图进行对比，研究景观格局变化对城市热岛的影响；徐双等[14]以Landsat7 ETM+数据和土地利用现状数据为数据源，分析了2004年和2010年长沙市热力景观空间格局演变；吴未等[15]以无锡市2010年TM数据研究了快速城市化平原地区不同景观指数随粒度的变化特征及适宜尺度的选择。

本文以石家庄市区作为研究区，选择1987年、2004年和2019年Landsat遥感影像数据，以城市景观格局为视角，深入研究景观格局指数与地表温度的相关性，可以掌握城市化发展对热环境的影响，可以为城市的合理规划、生态环境改善以及可持续发展提供理论依据和技术支撑。

1 研究区概况

本文以石家庄市区作为研究区(图1)。石家庄市地处河北省中南部，太行山东麓，位于37°27′～38°47′N，113°30′～115°20′E之间，属于温带半湿润半干旱大陆性季风气候， 太阳辐射的季节性变化显著，四季分明。近年来，由于城市工业化进程加快，使得城市绿地越来越少，城市景观结构发生改变，城市工业、交通运输和居民向大气中排放的人为热和废热日益增加，城郊之间某些气象要素已出现明显差异， 城市热岛效应愈发显著。通过研究景观格局指数与地表温度的相关关系，可以掌握城市景观对城市环境的影响，为石家庄城市建设、城市规划和生态平衡的改善提供一定的理论依据和技术支撑。

图1 石家庄市区地理位置

2 资料和方法

2.1 Landsat遥感数据及处理

本文所需数据资料主要涉及了Landsat卫星遥感资料、气象资料以及相关的地理信息资料等。遥感影像数据来源于国际科学数据服务平台[16]，条带号/行编号为124/34，数据格式为 .TIF，数据选取时间为1987年7月31日、2004年8月30日和2019年7月23日，其中1987年和2004年为Landsat5 TM数据，2019年为Landsat8 OLI_TIRS数据，具体遥感数据信息见表1所示；气象资料来源于河北省气象数据共享服务平台；地理信息资料由河北省地理信息局提供，数据格式为.shp格式。

表1 Landsat遥感数据信息

由于遥感影像在采集过程中不可避免地会受到传感器自身光学性质和结构特点等内部因素以及太阳辐射、地形条件等外部因素的影响，为保证研究结果的可靠性，需要对遥感数据进行预处理，主要包括波段合成、辐射定标、大气校正、几何校正和边界图裁剪等[17]。遥感数据预处理操作过程均在ENVI 5.2下进行，经过处理后，1987年、2004年和2019年遥感影像图见图2所示。

图2 1987年(a)、2004年(b)和2019年(c)石家庄卫星遥感影像

2.2 监督分类方法

监督分类方法又称训练场地法，是以建立统计识别函数为理论基础，依据典型样本训练方法进行分类的技术[18]。即根据已知训练区提供的样本，通过选择特征参数，求出特征参数作为决策规则，建立判别函数以对各待分类影像进行的图像分类，是模式识别的一种方法。判别准则若满足分类精度要求，则此准则成立；反之，需重新建立分类的决策规则，直至满足分类精度要求为止。目前比较常见的监督分类方法主要有平行六面体法、最大似然法、最小距离法、马氏距离法、神经网络法、支持向量机法等。本文利用最大似然法的监督分类方法对遥感影像数据进行分类处理，根据研究需要和4种类型的影像图的表征特点(表2)，将景观类型分为绿地、水体、不透水地表和未利用地等四种类型，分类结果见图3所示。通过精度计算得到1987年、2004年和2019年下垫面分类的总体精度分别为98.23%、98.79%和96.48%，Kappa系数分别为95.92%、97.43%和95.12%，总体精度和Kappa系数均在95%以上，表明分类精度很高，能够满足研究要求。

表2 4种下垫面类型的样本影像

图3 1987年(a)、2004年(b)和2019年(c)石家庄4类下垫面分类结果

2.3 地表温度反演算法

本文针对Landsat5 TM数据，采用单窗算法反演地表温度[19]；针对Landsat8 OLI数据，采用劈窗算法反演地表温度[20]。具体方法及参数请参考文献。其中， Landsat5 TM数据的热红外通道为第6波段，Landsat8 OLI数据是第10和第11波段。本文通过相同时间点的地面气象观测站点的地温数据作为验证数据源对反演精度进行检验， 发现单窗算法和劈窗算法的误差绝对值分别在1.2 K和0.8 K之内，反演结果能够满足研究需要，数据可用。

2.4 景观格局指数

景观格局指数能够比较精确地反映景观的结构组成和空间配置信息，是一个简单的定量指标。根据描述对象的等级，可分为斑块水平指数、斑块类型水平指数和景观水平指数[21-22]。本文利用Fragstats 4.2软件计算景观格局指数，在类型水平上选取8种景观格局指数，在景观水平上选取9种景观格局指数，具体描述见表3、表4所示。

表3 类型水平上选取的景观指数定义

表4 景观水平上选取的景观指数定义

2.5 移动窗口

移动窗口(Moving window)是通过对窗口内所选中景观指标进行统计计算，输出对应于所选景观指数栅格图，移动窗口算法[23]可以从空间角度展示出景观格局的动态变化过程。移动窗口算法的基本原理是将窗口(用户指定形状和大小)放置在每个栅格像元上，计算选定的度量并将值返回给栅格像元。因此，栅格像元周围的每个窗口都被视为子景观，并对每个单元重复此过程，移动窗口算法原理的示意,如图4所示，窗口a、b、c计算结果分别写入栅格1、2、3中，以次类推。

图4 移动窗口算法原理(窗口a、b、c计算结果分别写入栅格1、2、3中)

3 结果与分析

3.1 地表温度反演

图5给出了1987年、2004年和2019年地表温度反演结果。可以看出，1987年、2004年和2019年研究区地表温度空间分布均表现出中心高、周边低的现象，城市热岛效应比较明显。随着石家庄城市建设的快速发展，市区热岛发展较快，热岛强度和面积不断增加，东-西向发展明显；1987年和2004年滹沱河河道及附近存在一条比较显著的“高温带”，随着滹沱河的综合治理及蓄水工程的实施，“高温带”逐渐消失，现已成为石家庄市区的“生态冷源”。

图5 反演地表温度:(a)1987年，(b)2004年，(c)2019年

3.2 景观粒度的选择

景观粒度是景观中可识别的最小斑块单元，有关研究发现，景观格局指数会随粒度的变化而发生变化[24-25]。本文将研究区重采样为空间分辨率为30 m、60 m、90 m、120 m、150 m、180 m、210 m、240 m、

270 m和300 m的粒度大小，选取蔓延度指数(CONTAG)、景观导度指数(COHESION)、香农多样性指数(SHDI)和香农均度指数(SHEI)等景观格局指数，研究不同粒度下景观格局指数的变化。

图6给出了1987年、2004年和2019年不同粒度下景观格局指数的变化。从图6可以看出，1987年、2004年和2019年不同年份的景观格局指数大小不同，但变化趋势基本一致；四种景观格局指数均随粒度的不同而发生变化，蔓延度指数和景观导度指数变化较为明显，香农多样性指数和香农均度指数变化相对较小。研究表明，景观指数随粒度变化的第1尺度域是确定粒度大小，进行景观格局分析的适宜取值范围，根据研究区各景观指数变化曲线的拐点和跃变区间确定最佳适宜粒度[26]，从图6可以看出，在120 m、180 m、240 m和270 m处出现拐点，从图6b景观导度指数(COHESION)看，90 m～150 m范围为第1尺度域，210 m～270 m范围为第2尺度域，由此可以判定，90 m～150 m范围内为最佳适宜粒度，由于热红外波段的空间分辨率为120 m，经综合考虑，最终确定该研究区的适宜粒度为120 m。

图6 不同空间粒度下景观格局指数的变化:(a)蔓延度，(b)景观导度，(c)香农多样性，(d)香农均度

3.3 景观格局指数计算

景观格局指数能够比较精确地反映景观的结构组成和空间配置信息，是一个简单的定量指标。根据描述对象的等级，可分为斑块水平指数、斑块类型水平指数和景观水平指数。斑块水平指数往往被作为计算其他景观指数的基础，而其本身对了解整个景观的结构并没有很大的解释价值，通常采用的是斑块类型水平指数和景观水平指数。在获得了各种类型水平指数和景观水平指数的基础上，对单个时相的景观格局和不同时相的景观格局变化进行分析，说明景观的变化情况及变化趋势。利用Fragstats 4.2软件分别计算斑块类型水平指数和景观水平指数，结果见表5、表6所示。从表5可以看出，1987—2019年绿地的斑块类型面积(CA)整体呈逐年减少趋势，不透水地表和水体的斑块类型面积(CA)在逐年增加。1987—2019年绿地的斑块数量(NP)在逐年增加，绿地的破碎化指数(SPLIT)基本呈逐年增加趋势，绿地的最大斑块面积指数(LPI)在1987年、2004年和2019年分别为68.27%、12.40%和9.54%，呈逐渐减少趋势，这都说明研究区范围内绿地斑块在不断地裂化和减少。绿地的斑块形状指数(LSI)在1987年、2004年和2019年分别为19.69、22.15和36.37，说明绿地斑块形状不规则化程度加剧。绿地的聚集度指数(AI)在1987年、2004年和2019年分别为89.44%、86.51%和71.66%，说明研究区内绿地景观斑块的离散程度增大，绿地斑块的空间优势性大大降低。不透水地表的最大斑块面积指数(LPI)在1987年、2004年和2019年分别为12.09%、29.58%和54.18%，呈逐年增加趋势，不透水地表的破碎化指数(SPLIT)在1987年、2004年和2019年分别为68.15、11.41和3.40，呈逐年减少趋势，说明研究区范围内不透水地表斑块在不断的增加和聚合，形成优势斑块，说明城市化进程对城市景观格局影响较大。

表5 1987、2004、2019年类型水平上的景观指数

景观水平指数可以从全局的角度对研究区的景观水平进行分析，通过表6可以看出，1987—2019年最大斑块面积指数(LPI)呈现先减少、后增加的趋势，破碎化指数(SPLIT)与之相反，呈现先增加、后减少的趋势，说明随着城市化进程的加快，绿地斑块面积在逐渐减少和裂化，绿地的景观优势在不断下降，不透水地表的面积在逐渐增加和聚合，不透水地表的景观优势在不断加强，逐渐形成优势景观，这一结论合理解释了研究区范围内优势景观由绿地景观逐渐向不透水地表景观演变的过程。另外，1987—2019年聚集度(AI)逐渐减少，表明景观中不同类型的斑块离散程度在逐渐增大。

表6 1987、2004、2019年景观水平上的景观指数

3.4 移动窗口尺寸选取

景观格局指数的分析仅以数值描述景观格局，难以直接描述景观格局的空间形态，利用移动窗口可以实现景观格局指数的空间可视化，这就更加有利于统计分析景观格局指数与对应地表温度的关系。本文利用Fragstats 4.2的移动窗口功能得到景观格局指数空间分布图，为了更加科学地确定研究区移动窗口的网格尺度大小，本研究以2019年为例，统计分析了不同网格尺度下不同景观类型占比(PLAND)与地表温度的关系(图7)。从图7可以看到，绿地和不透水地表与地表温度的相关性较高，在360×360～1080×1080(m)范围内，基本呈稳定增加趋势，到1320×1320(m)则开始下降；水体则呈逐渐下降趋势，未利用地呈波动趋势，在600×600(m)处相关性最高，但未利用地与地表温度的相关性最差，这是与未利用地的不规则分布有关，在本文分析中将不再考虑未利用地的分析部分。经反复测试，600×600和360×360(m)相差不大，由于水体的跨度在1400 m左右，相关研究[23-25]表明一般选取其尺度的1/2作为窗口尺寸为宜，经综合考虑，最终确定以600×600(m)的正方形移动窗口进行景观格局分析。

图7 不同网格尺度下4种下垫面类型占比与地表温度的相关系数(绝对值)

3.5 不同景观格局指数之间的相关关系

以2019年为例，利用移动窗口法对2019年的遥感影像图进行计算，得到9种景观格局指数栅格图，并对其进行相关性分析，结果见表7。从表7可以看出，很多景观格局指数之间呈现出较高的相关性，如香农多样性指数(SHDI)和最大斑块面积指数(LPI)的相关性为-0.910，SHDI值越高，表明景观种类越多，斑块破碎化程度越高，其最大斑块面积指数LPI越小；斑块破碎化指数(SPLIT)和最大斑块面积指数(LPI)的相关性为-0.946，说明破碎化程度与最大斑块呈显著负相关，斑块越破碎导致LPI指数越小；斑块形状指数(LSI)和聚集度(AI)的相关性为-0.897，表明LSI与AI呈显著负相关，斑块形状指数增加使其不规则度增加，进而导致聚集度减少；斑块形状指数(LSI)和香农多样性指数(SHDI)的相关性为0.891，LSI增加使得其不规则度增加，最终导致斑块破碎，香农多样性指数SHDI增加。根据表7并考虑不同景观格局指数的生态学意义，本文选取景观占比(PLAND)、最大斑块指数(LPI)、斑块破碎化指数(SPLIT)和聚集度指数(AI)进行分析。

表7 2019年不同景观格局指数之间的相关系数

3.6 景观格局指数与地表温度的相关性分析

本文选取景观占比(PLAND)、最大斑块指数(LPI)、斑块破碎化指数(SPLIT)和聚集度指数(AI)，对2019年的遥感影像图进行计算和分析，从景观优势度、破碎化以及聚集程度等角度，分析了研究区范围内绿地、水体和不透水地表的景观格局指数与平均地表温度之间的相关性，结果见表8所示。

表8 4种类型景观格局指数与平均地表温度的相关系数

从表8可以看出，绿地、水体和不透水地表的景观格局指数与平均地表温度均表现出较好的相关关系，均为极显著性相关，其中绿地和不透水地表的相关性最好，水体次之。斑块百分比指数(PLAND)、最大斑块指数(LPI)和聚集度(AI)与平均地表温度表现出一致的极显著相关关系(绿地和水体为负相关，不透水地表为正相关)，破碎化指数(SPLIT)则恰恰相反(绿地和水体为正相关，不透水地表为负相关)。绿地和水体占比与平均地表温度相关系数分别为-0.844和-0.598，不透水地表占比与平均地表温度相关系数为0.911，绿地和不透水占比与平均地表温度的相关性高于水体，表明绿地为城市地表热环境效应抑制的主要贡献源，不透水地表是造成城市热岛的主要因素。平均地表温度与斑块百分比指数(PLAND)和最大斑块指数(LPI)的相关系数明显高于平均地表温度与聚集度(AI)和破碎化指数(SPLIT)的相关系数，这一研究结果可以说明，一个优势斑块(景观)对平均地表温度的影响效果明显大于几个比较分散或破碎的斑块(景观)。

由此可以得出，一个较大规模的绿地(或水体)产生的冷岛效应明显强于几个较小的分散绿地(或水体)；同理，规模较大或连片的不透水地表会比几个较为分散的不透水地表产生更强的热岛效应。由此可见，合理优化城市景观的空间布局是改善城市热环境、减轻城市热岛效应最有效的途径和手段。

4 结论

本文以石家庄市为例，基于1987年、2004年和2019年Landsat遥感影像数据，运用遥感、景观生态学、地理信息系统和统计分析等理论方法，分析景观格局指数与地表温度的相关性，得出结论如下：

(1)1987—2019年绿地CA逐年减少，NP、LSI和SPLIT逐年增加，LPI和AI逐年减少，说明研究区范围内绿地斑块在不断地裂化和减少，绿地斑块形状不规则化程度加剧，绿地景观斑块的离散程度增大，空间优势性不断降低；不透水地表和水体CA逐年增加，LPI逐年增加，SPLIT逐年减少，说明研究区范围内不透水地表斑块在不断地增加和聚合，形成优势斑块，城市化进程对城市景观格局影响较大。

(2)1987—2019年LPI呈现先减少、后增加的趋势，SPLIT先增加、后减少，说明绿地斑块面积在逐渐减少和裂化，绿地的景观优势在不断下降，不透水地表的面积在逐渐增加和聚合，不透水地表的景观优势在不断加强，解释了优势景观由绿地景观逐渐向不透水地表景观演变的过程。

(3)LST与PLAND和LPI的相关系数明显高于LST与AI和SPLIT的相关系数，说明一个优势斑块(景观)对LST的影响效果明显大于几个比较分散或破碎的斑块(景观)，一个较大规模的绿地(或水体)产生的冷岛效应明显强于几个较小的分散绿地(或水体)。