基于遥感的迁西县土地利用景观格局变化及脆弱性分析

曹应举 邢英梅 曹 琳 王 娜 郝雪景

(河北水文工程地质勘察院 石家庄 050000)

景观是由多个生态系统构成的异质性地域或不同土地覆被方式的镶嵌体[1]；景观格局是大小和形状各异的不同景观要素在空间上的分布与组合，是自然、生态过程和人为因素相互作用的结果[2-3]。景观格局变化能够体现空间上的斑块和景观变化，与区域生态环境变化密切相关[4]。因此，对景观格局动态变化的研究成为当今学者研究区域生态环境变化的必要前提。多年来，不少学者通过景观格局的数量和空间变化分析景观的演变特征，通过结合多种生态敏感性因子构建景观脆弱性指数来评价景观生态的好坏[5]。例如，曾光[6]等借助RS和GIS技术，对山西省1975～2014年间四个时间段的多源遥感数据采用景观指数的方法，进行湿地景观格局变化特征分析。结果表明，40年来山西省湿地面积减少、退化速率加快，湿地景观的破碎度指数先增后减，山西省湿地景观格局的变化归根于人类活动的影响。张荣天[7]等借助Landsat TM/ETM+数据，以镇江市为研究案例地，运用CA-Markov模型、动态度模型、景观指数以及PCA方法，对1994～2012年研究区景观格局演变特征、驱动机制以及变化趋势进行研究，通过定量化分析了解了区域土地利用景观格局演变的基本特征。董雅雯[8]等基于克里金插值法，通过构建的景观脆弱度指数对三亚市景观脆弱度进行评价，结果以不同等级景观脆弱区的分布特征和数量大小反映出三亚市整体受外界干扰的适应程度大小。李慧颖[9]等以吉林省长吉示范区为例，利用1976～2015年多源遥感数据的景观分类结果，对研究区森林景观的变化特征以及驱动因素进行分析，结果表明人类活动是导致研究区森林景观破碎化程度加剧的主要原因。

总结前人的研究成果，发现景观格局变化的研究大多是对人口、经济和气候等外在变化驱动力的分析，通过景观格局变化进行景观脆弱性评价的研究还略显不足。基于此，本研究以迁西县为例，基于土地利用动态度模型、景观格局指数分析方法，对2009年、2013年以及2017年遥感影像景观分类结果进行空间叠置对比分析，构建景观脆弱性指数对三期数据的景观脆弱性进行综合评价，旨在进一步分析研究区景观变化的动态特征以及发展方向。

1 材料与研究方法

1.1 研究区概况

迁西县位于河北省东北部，地处东经118°6′～118°37′，北纬39°57′～40°27′，东西横距39 km，南北纵距51 km，四季分明，属于暖温带大陆性半湿润的季风气候。迁西县是唐山市山区较多的县，地势南北高中间低，海拔70～831.3 m，境内有滦河、洒河、长河、清河等7条河流分布其中，土地利用类型多样，林地覆盖率达68 %，蕴含丰富的水资源、林地资源和矿产资源，其地貌呈典型的低山丘陵景观，因其独特的地形地貌特征被评选为国家级生态示范县。该县属环渤海经济圈，毗邻北京、天津、秦皇岛、承德等市，优越的地理位置和便捷的交通使得迁西县经济得以快速发展(图1)。

图1 研究区位置及高程

1.2 研究数据

研究所采用的数据来源于USGS美国地质调查局网站的Landsat TM/OLI影像，时段分别为2009年、2013年和2017年，获取月份分别为8月、7月和6月，云量覆盖率均小于1%且研究范围内无云量覆盖，影像投影分度带为UTM Zone 50N，基准面为WGS-84。其他辅助数据包括迁西县矢量边界确定研究区精确范围，Google Earth高分辨率历史影像用于目视解译选取高精度验证样本，ASTER GDEM 30 m分辨率高程数据，统计数据主要来源于2010～2016年《唐山统计年鉴》，用于对研究区土地利用分类结果进行对照分析。

为了精确反演地表辐射亮度值和反射率值，在ENVI5.1软件中对遥感影像进行辐射定标和大气校正，Landsat8-OLI影像的大气校正结果与全色波段利用Gram-Schmidt变换方法进行图像融合，再利用迁西县矢量文件进行边界裁剪，最终得到研究区范围。后期为达到增强和提取有用信息的目的，对研究区影像进行了分段线性拉伸和主成分分析，主成分分析保留了最大主成分波段。

1.3 景观分类及精度评价

为了系统分析迁西县景观格局动态特征，参照《土地利用现状分类》(GB/T 21010-2017)，结合迁西县景观格局的实际情况，将研究区景观类型分为：居民地、工矿用地、林地、水体、耕地、园地和其他用地七大景观类型。研究以eCognition软件为平台，运用面向对象和决策树相结合的分类方法，通过建立分类规则，区分开特征参数差异较大的地类，之后从易到难依次进行分类，根据制定的分类标准和分类规则，综合考虑影像的光谱、纹理、形状等特征，对迁西县的3期遥感影像进行多尺度分割与分类，最后对分类结果进行后处理(图2)。

图2 2009～2017年迁西县景观分类结果

借助谷歌地球高分辨率历史影像，随机选取足够数量且分布均匀的样本点对2009～2017年分类结果逐一进行精度评价。精度验证结果表明，迁西县土地利用总体分类精度分别为91.3 %(2009年)、92.0 %(2013年)和90.6 %(2017年)，均达到最低允许判读的精度，满足研究需要。

1.4 研究方法

(1)土地利用动态度模型

土地利用动态度反映的是研究区在一定时间范围内，土地利用类型变化的数量和速度[10]。应用土地利用动态度分析土地利用类型的动态变化，可以真实反映区域土地利用类型的变化快慢程度[11]。根据研究对象，可将土地利用动态度分为单一土地利用动态度和综合土地利用动态度。

单一土地利用动态度表达的是在特定时期内某一土地利用类型变化的程度和速率[12]，其计算公式如式1所示。

(1)

式中：K为研究时段内某一土地利用类型的动态度，ua和ub分别为该研究时段初期和末期的土地利用类型面积，t为研究时间间隔，土地利用动态度K反映了研究区在t时间间隔内的土地利用变化速率。

综合土地利用动态度表达的是土地利用变化的剧烈程度[13]，其计算公式如式2所示。

(2)

式中：LC为研究区土地利用类型的综合土地利用动态度，Si为研究初期第i类土地利用类型的面积，ΔSi-j为研究时段内第i类土地利用类型转为非i类土地利用类型面积的绝对值，t为研究时间间隔。

(2)景观格局指数的选择

景观格局指数能够高度浓缩景观空间格局信息，并且能够反映景观结构组成和空间配置等方面的特征，是研究景观格局时空变化的关键数据[14]。根据研究需要同时结合研究区特点，筛选出适合本研究的多种景观指数。在斑块类型水平上选取斑块数量、最大斑块指数、边缘密度、斑块密度、景观形状指数以及分散和并列指数；在景观水平上选取蔓延度指数、景观连接度指数、聚合指数和香农多样性指数。以期对研究区的景观破碎程度进行简单度量，同时计算出丰度密度指数、景观均匀度指数、面积加权平均分维数、斑块类型面积4个重要的景观指数，用以对研究区脆弱性进行时空变化分析。研究所用的指数通过景观结构分析软件Fragstats 4.2计算得到。

(3)景观干扰度指数

景观干扰度指数由景观类型破碎度、分维数倒数和优势度共同决定，其计算公式如式3所示。

Si=aFNi+bFDi+cDOi

(3)

式中：景观类型破碎度FNi，分维数倒数FDi，优势度DOi，a,b,c分别为各个指数对景观干扰度指数的贡献量大小，本研究三个参数分别赋值为0.5、0.3和0.2。

(4)景观敏感性指数

景观敏感性指数是通过景观干扰度指数和景观类型易损度指数的乘积相叠加，用来表征景观受外界干扰敏感程度的指标[16]，景观敏感性指数的计算公式如式4所示。

(4)

式中：LSI为景观敏感性指数，i为景观类型，n为景观类型数目，Si为景观干扰度指数，Vi为景观类型易损度指数。参考他人研究成果[17]，将各景观类型按易损程度由低到高分别赋值：居民地为1，工矿用地为2，林地为3，园地为4，耕地为5，水体为6，其他用地为7，最后进行标准化处理得到每个景观类型的易损度指数。

(5)景观适应性指数

景观脆弱性不仅与景观敏感性有关，还与景观的适应和恢复能力有关。一般情况下，景观多样性越大，结构分布越均匀，其景观适应性越强。本文借鉴孙才志[18]等的研究，引入景观适应性指数用来描述景观的抗外界干扰及适应能力，景观适应性指数LAI计算公式如式5所示。

LAI=PRD×SHDI×SHEI

(5)

式中：PRD为斑块丰度密度指数，SHDI为香农多样性指数，SHEI为香农均匀性指数。

(6)景观格局脆弱性评价指数

为了直观反映研究区2009～2017年三个时期的景观格局空间变化对生态环境产生的负面影响，利用Fragstats 4.2软件计算出的景观指数构建景观脆弱性指数。

系统的脆弱性通过系统面对外界干扰的敏感程度与适应能力来体现，根据引入景观敏感性指数和景观适应性指数构建可以评价景观生态安全的景观脆弱性指数LVI[15]，其计算公式如式6所示。

LVI=LSI×(1-LAI)

(6)

式中：LSI表示景观敏感性指数，LAI表示景观适应性指数。

2 结果与分析

2.1 景观格局变化特征分析

基于获得的土地利用分类结果，采用GIS空间分析功能统计各时期土地利用类型的变化面积，根据公式(1)和公式(2)，可以得出2009～2013年、2013～2017年迁西县单一土地利用动态度和综合土地利用动态度(表1)。

表1 2009～2017年迁西县土地利用结构及其变化

从迁西县土地利用动态变化表可以看出：2009～2013年间，迁西县土地利用空间变化面积最大的地类是工矿用地，变化总量为2489.37 hm2，其次是林地，变化总量为-2375.29 hm2，其他用地面积变化最小，变化总量为-44.96 hm2。在单一动态度方面，工矿用地动态度最大，其数值为14.44%，其次分别是园地、水体、居民地、耕地和其他用地，其数值分别为10.50%、3.92%、3.91%、-2.36%和-1.35%，由于林地面积比重最大，使得其面积变化量较高、单一土地利用动态度最小，林地的单一土地利用动态度为-0.59%。2013～2017年间，面积变化最大的地类是水体，变化总量为-2033.35 hm2，其次是其他用地，变化总量为1916.06 hm2，居民地面积变化最小，变化总量为-254.52 hm2。在单一动态度方面，其他用地动态度最大，其数值为61.02%，其次分别是水体、园地、工矿用地、耕地和居民地，其数值分别为-8.77%、-3.91%、-3.16%、1.05%和-0.4%，林地的单一土地利用动态度最小为0.21%。以上数量变化表明在2009～2013年期间工矿用地增加速度最快，园地减少速度最快，工矿用地和居民地面积的大幅增加，主要由林地、耕地和园地转化而来，其转化原因与经济增速发展和人口数量的增加有关。在2013年到2017年期间其他用地增加速度最快，水体减少速度最快。

从总体来看：2009～2017年间，迁西县各土地利用类型之间变化较为频繁，综合土地利用动态度为1.64%。其中，耕地、水体、林地、园地的面积呈减少趋势，且减少量依次上升，工矿用地、其他用地、居民地的面积呈增加趋势，且增加量依次上升。从变化数量看，居民地的大幅增加和园地的大幅减少最为显著，8年来居民地面积净增1911.68 hm2，年变化量为238.96 hm2，园地净减少1611.94 hm2，年变化量为201.49 hm2；从变化速度来看，其他用地的变化速度最快，平均每年增加28.18%，其次是园地，平均每年减少6.39%，工矿用地平均每年增加4.73%，水体平均每年减少3.11%，居民地平均每年增加1.72%，耕地平均每年减少0.7%，林地由于占总体面积比重较大，是变化速度最小的用地类型，平均每年减少0.19%。就整个研究区而言，2009～2013年土地利用行为相对稳定有序，2013～2017 年土地利用类型间流转程度较为剧烈。

2.2 景观指数分析

利用ArcGIS将土地利用分类结果转化为栅格图，为了保留原有数据的细节信息，栅格数据的像元大小设置为原始影像数据的分辨率大小，利用Fragstats 4.2软件进行景观指数的计算，最终绘制斑块类型水平下各景观指数的雷达图(图3)。

图3 斑块类型指数雷达变化图

(1)林地

在研究时段内，林地的最大斑块指数和边缘密度最大，并且在2017年达到峰值，而林地的分散和并列指数变化不大，说明林地景观一直是迁西县的优势景观。同时，林地与其他斑块类型的邻接程度较高，斑块数量持续增加说明林地景观的破碎化程度也在增加。

(2)耕地

耕地的斑块数量、边缘密度、斑块密度和景观形状指数一直上升，而分散和并列指数一直下降，表明耕地图斑受地形起伏的影响形状更加复杂，与其它景观类型的嵌合度大，分散性强，团聚性差，破碎化程度逐年增加，相比较耕地与林地，虽然耕地的占有面积较少，但是由于整体破碎度高，空间结构更复杂，所以景观形状指数高于林地。

(3)水体

水体的分散和并列指数最大并且保持上升趋势，表明水体景观的破碎化程度较低，人类活动对水体景观的影响不明显，水体景观维持在一个较为稳定的状态，其他类型的斑块受水体分布广的影响，与水体斑块邻接程度很高，其分散和并列指数值也普遍较高。另外，2009年到2013年研究区内斑块密度和边缘密度都有所增加，主要是因为研究所用的2013年遥感影像获取时间正值雨量充沛的7月份，沟渠、水库蓄水量增大，淹没了周围区域，使得该景观斑块面积增大。从2013年分类图上可以看出水体形状变得狭长曲折，且比较复杂，这使得水体的景观形状指数增大；而2013年到2017年研究区内斑块密度和边缘密度都有所降低，形状指数转而减小，说明水体自身形成了较大的斑块，水体的形状趋于规则和简单。

(4)居民地和工矿用地

居民地和工矿用地的斑块数量、边缘密度、斑块密度和景观形状指数先上升后下降总体趋于上升，但分散和并列指数先下降后上升总体趋于下降，说明人类活动的干扰造成不同斑块间的相互镶嵌。居民地和工矿用地的破碎化程度一直在增加、斑块形状越来越复杂。但从2013年开始，居民地和工矿用地的破碎化程度有降低的趋势，斑块形状趋向规则、简单，说明城镇化发展中经过合理的规划引导和空间资源的有序配置，居民地和工矿用地景观格局朝着向好的方向发展。

(5)指数分析

景观水平指数主要识别整个景观镶嵌体的空间特征，为了了解迁西县的整体景观空间构型，计算2009、2013和2017年的景观水平指数如表2所示。

表2 不同时期的景观水平指数

从景观水平上来看，研究时段内蔓延度指数呈递减趋势，在2017年达到最小，其值为65.083，说明研究区中林地景观斑块类型连接性减弱，研究区内不同斑块类型的团聚程度在下降，空间分布逐渐趋向分散，研究区景观的异质性增强；景观连接度指数均在99.4%以上，并呈现出逐年递增的趋势，说明整个研究区是以少数大斑块为主体构建起来的，其他不同斑块之间的连接度较为紧密，研究区内景观破碎度有所改善；2009年到2013年聚合指数和Shannon多样性指数呈增长趋势，2013年以后呈降低趋势，说明2013年前迁西县景观类型趋于丰富多元化，各种景观类型的发展变化，使得景观中存在较多的小斑块，景观破碎化程度增加；2013年后景观格局连通性向好的趋势发展，景观中不同斑块类型聚集程度有所减小，各种景观类型之间逐渐形成较好的连接性，且景观格局破碎度缓和，朝着均匀化和多样化方向发展。

2.3 景观脆弱性分析

利用网格分析法对研究区景观脆弱性进行评价，综合考虑研究区的大小和采样区的工作量，将研究区分割为4 km*4 km的单元格网，按照从上到下从左至右的顺序依次编号，生成编号为1-119的119个景观类型图文件。同时，每个方格网中的景观类型图在ArcGIS中批量转换为栅格文件，在Fragstats 4.2中导入前期在Excel中制作好的batch文件批量计算景观指数，计算每个格网的脆弱性指数并将其作为每个格网中心点的脆弱性值，以每个格网中心点的ID标识码和格网编号为连接字段建立属性连接关系，利用统计学中的克里金空间插值法绘制景观格局脆弱性分布图，最后按照相等间隔将研究区划分为5个等级分区：低脆弱区、较低脆弱区、中等脆弱区、较高脆弱区和高脆弱区。该尺度下的景观脆弱性分布如图4所示，统计研究区5个等级分区所占面积比例见表3。

图4 4 km*4 km采样区景观格局脆弱性分布图

表3 景观脆弱分区面积比例

结合图4和表3可以看出：2009～2017年研究区景观脆弱性整体呈现出以较低脆弱区、中等脆弱区和较高脆弱区分布较多，而低脆弱区和高脆弱区分布较少的空间布局，且整体脆弱性程度从西北到东南大致呈现降低趋势。

从景观脆弱性不同等级来看，研究区2009年低脆弱区和较低脆弱区都不具优势，占整个研究区的面积比例分别为5.72%、17.86%，主要分布在研究区东部及南部边缘区域，该区域零散分布有林地、耕地、居民地和工矿用地，景观适应性提高，敏感性和脆弱性较低；中等脆弱区占研究区面积比例为21.71%，主要呈半包围状分布于研究区中部地区，该区域林地、居民地和工矿用地相互嵌套，耕地景观较为聚集，景观适应性处于一般水平，多样性不明显，导致该区域有一定的脆弱性；较高脆弱区和高脆弱区占研究区面积比例分别为37.34%、17.37%，主要分布在研究区北部、西北部及中部地区，以林地为主，耕地为辅，其景观类型相对单一，该区域也是滦河流经的区域，滩涂等其他用地面积变化较大，导致景观脆弱性值较大。

与2009年相比，研究区2013年低脆弱区和较低脆弱区面积比例都有所增加，空间分布从研究区东部及东南部边缘区域向中部及南部扩展，占据2009年部分中等脆弱区和较高脆弱区；中等脆弱区面积比例较2009年有所增加，其变化走向朝西挤压并向东北方向延伸，占据了一些景观较高脆弱区和高脆弱区；处于较高脆弱区和高脆弱区的面积比例大幅减少，较高脆弱区空间位置整体向西部靠拢，高脆弱区变动剧烈，由北部及西北部区域迁移至西南区域，在东南部区域也出现小范围的较高脆弱区和高脆弱区。

2017年，研究区内处于低脆弱区面积比例较2013年大幅下降，由14.67%下降至3.10%，较低脆弱区面积比例有所增加，由23.34%上升到25.19%，两者的空间位置从中部及南部区域逐步向东部边缘地区扩展，出现在2013年东南部区域的较高脆弱区和高脆弱区完全消失，被低脆弱区和较低脆弱区占据。结合景观分类图可以发现，造成该片区域景观脆弱性发生剧烈变化的原因是水体景观和耕地景观变化导致的生态环境失衡、脆弱性升高，经过综合整治后的水体景观和耕地景观朝向好的方向发展，景观适应性提高，从而脆弱性降低；研究区中等脆弱区的面积比例较2013年持续增加，由27.48%上升到29.69%，造成该变化的原因是，自2013年以来耕地景观的合理规划加上对土地资源的高效整合，造成中等脆弱区面积增大，空间分布呈半包围趋势朝西南方向扩展；相较2013年、2017年研究区处于较高脆弱区的面积比例减小，由30.22%下降至21.16%，高脆弱区面积比例大幅增加，由4.29%上升到20.86%，迁西县中西部山峦环抱，滦河干流纵贯迁西县境，两山夹一河的地形条件在加快当地旅游业发展的同时也在影响着各类景观要素的稳定性，该区域受人类活动的影响较大。因此，对外界的干扰较为敏感，较高脆弱区大多分布于该区域；同时结合景观分类图可以看出，研究区高脆弱区域分布有大量的其他用地景观，景观类型的频繁更替导致该区域的景观脆弱性达到最高。

3 结语

对迁西县景观变化及脆弱性分析应考虑以下方面：虽然描述景观格局的指数有很多，但许多景观指数之间存在很大的相关关系，因而在景观格局研究中，选取信息不重复的景观指数是值得探讨的问题。自然因素和人为因素是迁西县景观格局脆弱性发生变化的主要原因，随着一系列生态保护政策的实施，迁西县景观脆弱性会逐渐稳定。此外，本文所构建的景观脆弱性评价指标中，低脆弱性并不表示生态环境的质量优良，仅代表了区域生态处于相对稳定的状态，可恢复性较好。

迁西县土地利用从变化数量来看，以居民地的大幅增加和园地的大幅减少最为显著；从变化速度来看，其他用地的变化速度最快，然后依次为园地、工矿用地、水体、居民地、耕地、林地，林地是变化速度最小的用地类型。但大体上迁西县耕地、水体、林地和园地的面积仍呈逐渐减少趋势，工矿用地、其他用地和居民地的面积逐步增加。

通过斑块水平和景观水平对迁西县景观特征进行分析，迁西县林地景观和耕地景观的斑块数量一直处于上升趋势，景观边界越来越复杂，破碎化程度逐渐增长；水体景观受降雨量等的影响，形状趋于简单和规则，破碎化程度逐渐降低；居民地景观和工矿用地景观在人类活动的干预下，形状呈现出由复杂向简单过渡的趋势，破碎化程度逐渐得到缓和；其它景观破碎化程度变化不大。在景观水平上，连接度指数和Shannon多样性指数增加，蔓延度指数和聚合指数减少，表明迁西县景观趋于丰富多元化，景观异质性减少，破碎化程度朝着均匀化和多样化方向发展。

迁西县景观脆弱性呈现出由西北到东南逐步降低的趋势，且脆弱性西部高于东部，北部高于南部，中西部分布有大量的其他用地景观，景观类型的频繁更替导致高脆弱区面积有所扩张。但从整体来看，较低脆弱区和中等脆弱区面积比例上升，较高脆弱区所占面积比例下降，景观脆弱性较严重的地区趋于分散，逐渐被脆弱性较低的区域所占据，迁西县景观脆弱性整体有所缓和，但局部区域仍需加强治理，保证该区域的可持续发展。