基于PLUS模型的金佛山自然保护地毗邻区土地利用变化模拟与多情景预测

杜文武 胡 瑶 刘渝杰 眭 淼

（西南大学园艺园林学院，重庆 400715）

位于自然保护地边界以外的毗邻区，土地超载、破坏性利用及管控乏力现象普遍，对自然保护地本体造成严重威胁。科学预测土地利用变化趋势，是理解毗邻区土地利用变化对自然保护地本体生态安全影响，进而制定针对性管控措施的重要前提。运用PLUS模型，基于2000年、2010年、2020年三期土地利用数据，分析2000 - 2020年重庆金佛山自然保护地毗邻区土地利用变化时空格局及其驱动力，并设置自然发展、城镇发展、生态保护三种情景，模拟预测2040年保护地毗邻区土地利用情况。结果表明：（1）金佛山自然保护地毗邻区土地利用类型以林地、耕地为主；2000 - 2020年间，耕地、林地面积减少，建设用地面积大幅增加，转入来源主要为林地和耕地；（2）PLUS模型对金佛山自然保护地毗邻区土地利用驱动因子的解释能力较好，其中，对耕地、林地、建设用地面积变化解释力度最大的驱动因子均为NDVI指数和人口密度；（3）PLUS模型在金佛山自然保护地毗邻区的适用性较强，总体精度为0.95，Kappa系数为0.92，远高于0.75；（4）三种情景下，地类变化趋势均为林地、耕地面积减少，建设用地面积由保护地毗邻区逐渐向本体边界处蔓延。研究表明，PLUS模型适用于同类型自然保护地毗邻区土地利用模拟研究；仅对毗邻区生态类用地进行保护限制不足以遏制建设用地向本体扩张，亟需对毗邻区人类活动进行更为精细化的管控。

自然保护地；毗邻区；PLUS模型；土地利用变化；模拟预测

自然保护地是国土生态安全的基石，在维护国家生态安全中居于首要地位[1]。位于自然保护地边界外、作为自然保护地本体生态屏障的毗邻区，因土地超载、破坏性利用及管控乏力现象对自然保护地本体的生态安全造成严重威胁，已成为国内自然保护地体系构建过程中的现实难题[2-3]。在自然保护地建设发展过程中，由于城镇化扩张、旅游开发、乡村居民点建设等情况，自然保护地毗邻区人类活动加剧，对自然保护地内部的生态过程产生干扰。

土地利用/土地覆被变化能够直接表现人类与自然生态环境的相互作用关系，是衡量和判断毗邻区受人类活动扰动状况的关键[4]。目前，土地利用变化相关的研究内容主要集中于土地利用格局动态变化[5-6]、土地利用时空格局变化[7-8]、土地利用驱动机制[9-10]、土地利用多情景模拟预测[11-13]等。随着地理信息系统和遥感技术的发展，与土地利用研究相关的模型不断发展，代表性的数量预测模型有马尔可夫模型（Markov）、人工神经网络模型（ANN）等，代表性的空间预测模型有元胞自动机（CA）、CLUE-S模型、FLUS模型等[14-16]。

文章选用在数量预测与空间预测上具有双重优势的斑块生成土地利用变化模拟模型（Patch-generating Land Use Simulation Model，PLUS）[17]，对2000、2010、2020年金佛山自然保护地毗邻区土地利用时空变化特征及驱动力进行分析，并设置自然发展、城镇发展、生态保护三种情景，预测2040年金佛山自然保护地毗邻区土地利用情况，以探究毗邻区土地利用变化对自然保护地本体的影响，为制定针对性管控措施提供参考。

1 研究区概况

金佛山位于重庆市南川区境内28°50′-29°20′N，107°00′-107°20′E，地处渝、黔两省交界处，是四川盆地东南缘与云贵高原的过渡地带，由金佛、柏枝、箐坝、三元（庙坝）4片108座山峰组成，总面积达1 300 km2，海拔在700～2 251 m之间。气候属亚热带温润季风气候，年平均气温8.3℃，年平均降水量约1 400 mm，相对湿度可达90%，云雾较多。金佛山是中国早期建立的自然保护区之一，生物多样性丰富，是中国珍贵的生物物种基因库和天然植物园[18]。金佛山地质古老，属于典型的喀斯特台原地貌，2014年作为典型代表被列为世界自然遗产。

金佛山存在数种自然保护地类型[19]，根据各类型边界重叠情况及具体生态保护目标。综合考虑人类活动干扰情况、各保护地功能分区情况、管控要求等，最终选择金佛山国家级自然保护区的核心区和缓冲区、金佛山风景名胜区的核心景区及金佛山世界自然遗产的核心区作为金佛山自然保护地本体，并以此范围边界缓冲15 km的区域作为金佛山自然保护地毗邻区[20]（图1）。

图1 研究区域Fig.1 Study area

2 研究方法

2.1 数据收集和处理

本研究的数据主要包括土地利用数据、自然条件数据及社会经济数据，数据具体情况见表1。土地利用/土地覆被基础数据参照《土地利用现状分类标准》（GBT21010-2017），运用三级分类系统，根据研究区实际情况及研究的具体内容，土地利用类型分为耕地、林地、草地、水域及建设用地；DEM（Digital Elevation Model，数字高程模型）数据为ASTER GDEMV2数据，坡度数据在其基础上获得；归一化植被指数（NDVI）是基于Landsat遥感影遥感影像数据，经过辐射定标、大气校正、条带修复、拼接裁剪等处理后，根据公式计算得来；人口、GDP数据经过重采样至30 m分辨率；建筑数据根据不透水面数据、土地利用数据及卫星影像数据进行校正获得。

表1 数据信息Tab.1 Data information

2.2 土地利用/覆被变化时空格局分析

2.2.1 土地利用转移矩阵

土地利用转移矩阵可用于分析研究区域在一定时间段内，一种土地类型向其他用地类型转换的情况，可以直观且定量地分析出各地类转化数量和转换方向[22]。计算见公式（1），式中S为面积；i，j分别为研究初期与研究末期的土地利用类型；n为研究区土地利用类型的数量。

2.2.2 土地利用动态度

土地利用动态度可用于揭示研究区域内各类土地面积的综合变化速率[23]，体现各类用地之间变化强度情况及之间存在的差异，其中包含了单一土地利用动态度和综合土地利用动态度。

单一土地利用动态度计算见公式（2），式中K为所求时间段内其中一种用地类型的土地利用单一动态度；Ua、Ub分别为该地类研究时段初期和末期土地利用类型面积；T为研究时段的间隔时长，单位为年。

综合土地利用动态度计算见公式（3），式中LU为所求时段内综合土地利用动态度；LUi为研究初期第i类土地利用类型面积；LUi-j为第i类土地利用类型转为第j类土地利用类型面积的绝对值；T为研究时段长度，单位为年。

2.3 PLUS模型

PLUS模型是利用斑块模拟土地利用状况的模型，其耦合了一种新的用地扩张策略分析（Land Expansion Analysis Strategy，LEAS）和基于多类随机斑块种子的CA模块（CA Model Based on Multi-type Random Patch Seeds，CARS），通过选择土地利用变化的驱动因子并计算其解释力，进一步测度土地利用的变化规律，并以此模拟预测未来土地利用格局[24]。PLUS模型作为土地利用模拟模型具有精度高、速度快的特点[17]。

2.3.1 LEAS

LEAS能提取分析各类土地利用类型中扩张部分及其发展概率情况。其通过随机森林模型，基于驱动因子栅格数据得出各个因子对研究区土地利用格局演变的贡献度。LEAS模块融合了多种算法优势，可以更好地处理多重共线性问题，对土地利用扩张具有更好的解释作用。

土地利用变化的驱动机制复杂多样。本文综合考虑金佛山自然保护地毗邻区土地利用变化的自然环境影响因素和社会经济影响因素[25]，结合前人研究，选取高程、坡度、NDVI指数、年均降水、年均气温、GDP、人口密度、建筑面积指数、距道路距离共9个驱动因子。

2.3.2 CARS

CARS结合CA模型、随机种子生成机制和阈值递减机制，使PLUS模型在模拟运行过程中受到土地利用扩张发展概率约束，从而模拟出研究区未来土地利用数据。CARS模型的参数设置主要涉及邻域权重因子、随机斑块种子模型、自适应系数等，本文基于前人经验[26]和研究区实际情况，将CARS各部分参数设置如下：（1）邻域效应，以30 m×30 m栅格数据为研究区域范围，将领域范围设置为3，递减阈值设置为0.5，膨胀系数设置为0.1；（2）转化规则，转化成本表格见表2；（3）邻域权重，经过多次修正测试和精度检测评估，最终邻域因子参数设置见表3。

表2 转化成本规则Tab.2 Conversion cost rules

表3 邻域权重设置Tab.3 Neighborhood weight setting

2.3.3 模拟结果与精度检验

为综合评估PLUS模型模拟预测结果的可信度，本研究利用2020年模拟预测数据与实际现状数据进行以下两种方式的验证：（1）混淆矩阵分析。将模拟结果与实际现状相对应栅格进行单元一一对应的混淆矩阵模型，通过混淆矩阵了解栅格数据中模拟预测产生的差异；（2）Kappa系数检验。通过计算模拟预测之后对比实际现状的总体精度以及Kappa系数大小，测算出模拟预测及实际现状在空间栅格单元中的一致性大小。Kappa系数具体计算公式见公式（4），式中P1为模型模拟预测得到正确结果的像元栅格比例，P2表示栅格数据在随机条件状态下期望正确的像元栅格比例。在Kappa系数检验中其数值取值范围为0 ～1，Kappa系数越大表示模型检验结果与实际现状的一致性越高，也反映出模拟预测结果精度越高。其判断标准通常为：Kappa系数大于0.75时一致性较高，模拟预测精度较高；0.4 ～0.75范围内一致性一般；小于0.4时一致性较差，模拟精度较差。

2.4 情景设置

为探究不同发展目标下金佛山自然保护地及其毗邻区土地利用变化情况，根据前人经验[27-28]，本文设置自然发展情景、城镇发展情景和生态保护情景，从这三种未来发展情景角度模拟预测2040年金佛山自然保护地毗邻区土地利用格局。各情景成本设置如表4所示。

表4 各情景土地利用转化成本规则矩阵Tab.4 Cost rule matrix of land use transformation in different scenarios

3 结果分析

3.1 2000-2020年金佛山自然保护地毗邻区土地利用变化时空格局分析

3.1.1 土地利用空间分布特征

经过修正的2000 - 2020年金佛山自然保护地及其毗邻区土地利用分类数据如图2及表5所示。其中林地面积占比最高，其次为耕地，草地、水域、建设用地面积占比较小。20年间自然保护地本体区域土地利用类型以林地为主，占比持续在90%以上；耕地面积持续减少，占比在7%左右；建设用地面积少量增加，但总占比较少。自然保护地毗邻区土地利用类型仍以林地为主，占比约61%；耕地面积减少较多，总体占比约31%；建设用地面积占比增加至2.23%，空间上向东北方向逐步迁移至南城街道及三泉镇等靠近自然保护地边界的区域，不利于自然保护地本体自然保护，亟需重视毗邻区土地利用变化动态。

表5 土地利用地类面积变化分析表Tab.5 Analysis table of land use land type area change

图2 2000 - 2020年土地利用类型图Fig.2 Land use type map from 2000 to 2020

3.1.2 土地利用转移矩阵

根据表6，2000 - 2020年间，建设用地、耕地、林地的地类转化关系较为显著。其中，建设用地明显增加，2020年的面积较2000年增长两倍多，主要流入来源为林地，转移面积为2 486.25 hm2，占建设用地转入面积的75.37%；耕地面积减少，主要转化为建设用地和林地，转移面分别为2 486.25 hm2、1 378.35 hm2，分别占耕地面积转出面积的58.54%、32.46%。根据2000 - 2020年间自然保护地毗邻区空间分布地类转换图像（图3），2000 - 2010年期间转为建设用地的部分主要位于毗邻区北部的水江镇和西部的万东镇，转向耕地的部分主要位于东北部三桥镇下辖的乡村，可能为乡村开垦土地，复耕等；2010 - 2020年期间建设用地转入十分明显，不仅包括原有的水江镇和万东镇的进一步村镇建设，还有南城街道的进一步扩张，在靠近自然保护地本体边界的建设用地也在进一步增加，包括黑山镇和金山镇的建设扩张。

表6 自然保护地毗邻区土地利用转移矩阵（单位：hm2）Tab.6 Land use transfer matrix for adjacent areas of protected areas

图3 自然保护地毗邻区空间分布地类转换图像Fig.3 Spatial distribution of adjacent areas of protected areas and conversion images of land types

图4 驱动因子贡献度Fig.4 Driving factors contribution

3.1.3 土地利用动态度

根据表7，金佛山自然保护地毗邻区综合土地利用动态度整体较低，2010 - 2020年间土地利用类型变化频繁，综合土地利用动态度为0.079%。2000 - 2010年间，各类用地变化幅度为建设用地＞草地＞耕地＞水域＞林地，建设用地变化最为显著，单一土地利用动态度为3.32%，草地、耕地、水域面积减少，变化速度较缓，林地面积增加，但增速仅为0.05%，由于毗邻区林地占较多，因此变化速率不明显；2010 - 2020年间，各类用地变化幅度为建设用地＞耕地＞林地，其中耕地和林地面积减少，单一土地利用动态度分别为-0.18%和-0.03%，建设用地增长速率进一步加快，达到5.47%，几乎为前10年的近一倍。

表7 自然保护地毗邻区土地利用动态度（单位：%）Tab.7 Land use dynamic degree in the adjacent areas of protected areas

3.2 基于PLUS模型的金佛山自然保护地及其毗邻区土地利用变化驱动力分析

PLUS模型中LEAS模块计算各个驱动因子对每种土地利用类型具体解释力程度，如图5所示。对于耕地来说各类驱动因子相对比较平均，其中影响贡献度最大的是NDVI指数，其次为人口密度，影响力最小的是年均气温；对于林地来说，影响力较大的驱动因子依次为NDVI指数、人口密度以及GDP，贡献度相对较小的驱动因子为距道路距离以及高程；对于建设用地来说，各驱动因子之间最值差距最大，其中NDVI指数贡献度最大，可以达到0.25以上，人口密度与高程、建筑面积指数同样对建设用地贡献度较大。

图5 2020年土地利用模拟对比图Fig.5 Comparison chart of land use simulation in 2020

3.3 PLUS模型精度验证

根据PLUS模型得出2020年金佛山自然保护地毗邻区土地利用模拟数据，并与2020年土地利用实际数据对比（图5，表8）。然后选取其中三个区域，将模拟与实际的土地利用数据扩大分析进行比对，可以看出模拟得到的土地利用类型数据在空间分布中较为一致，可信度较高。

表8 2020年实际土地利用模式与预测模式的混淆矩阵（单位：hm2）Tab.8 The confusion matrix between the actual land use model and the forecast model in 2020

依据土地利用模拟结果数据与实际现状数据进行精度检验计算，得出Kappa系数和总体精度分别为0.91和0.95，其中Kappa系数超过0.9，大于0.75，可以反映出模型模拟预测精度较高，使用模型预测可信度高。由此可见，在PLUS模型中进行的相关参数设定具有合理性，可以在后续的土地利用模拟预测中使用。

3.4 多情境模拟下金佛山自然保护地及其毗邻区土地利用分析

3.4.1 土地利用多情景模拟

本文结合筛选出的驱动因子及设定参照规则，运用PLUS模型得到2040年土地利用模拟预测的结果（表9）。自然发展情景下（图6-a），对比2020年，2040年耕地、林地面积在20年间将分别减少2 871.99 hm2、1 114.92 hm2，草地及水域面积变化相对较小，建设用地面积扩张较为明显，由6 365.79 hm2增至10 517.40 hm2，共增加4 151.61 hm2，表明在不考虑政策变迁、重大生态变化等干扰因素的情况下，建设用地将在原有村镇的基础上不断地扩散式增长，增加自然保护地本体的保护压力；城镇发展情景下（图6-b），对比2020年，2040年耕地、林地、草地面积下降，其中耕地及林地面积变化较为显著，分别降低4 975.47 hm2、1 114.92hm2，建设用地面积增加十分明显，增至2020年的两倍，共6 255.09 hm2，表明在对建设用地减少转出、增加转入，同时不考虑调整生态保护政策的情况下，建设用地将会侵占耕地、林地、草地，威胁自然保护地本体的生态安全和毗邻区的粮食安全；生态保护情景下（图6-c），对比2020年，2040年林地、草地、水域面积分别增加了7 613.46 hm2、591.39 hm2、5.31 hm2，建设用地面积增速较缓，共增加4 151.61 hm2，耕地面积明显减少，共减少12 361.77 hm2，表明建设用地的蔓延趋势得到缓解，林地等生态类用地得到了保护。

表9 2040年三种情景下土地利用需求预测表（单位：hm2）Tab.9 Land use demand forecast table under three scenarios in 2040

图6 2020 - 2040年三种情景下土地利用类型模拟图Fig.6 Simulation map of land use types in three scenarios from 2020 to 2040

3.4.2 多情景下建设用地变化趋势

对比2020年，本文进一步分析了三种情景下2040年金佛山自然保护地毗邻区的建设用地对于自然保护地本体趋近压迫的情况，分别得出建设用地重心迁移和标准化椭圆（图7）。分析可知，三种情景下，建设用地面积由自然保护地毗邻区周边逐渐向本体边界处蔓延，其中主要以原有的村镇如南城街道、水江镇、黑山镇等逐步向内连片发展，金山镇、大有镇、头渡镇、德隆镇等靠近自然保护地本体的村镇也随着经济社会发展逐步向边界蔓延逼近。这意味着自然保护地边界将会有部分区域生态质量受到人类活动的影响，若不加以重视，将会严重威胁自然保护地本体生的生态安全。

4 讨论与结论

本文基于2000 - 2020年金佛山自然保护地毗邻区土地利用现状，分析金佛山自然保护地毗邻区土地利用的时空变化，使用PLUS模型进行各类用地扩张的驱动力分析、土地利用模拟验证及2040年土地利用多情景预测。金佛山自然保护地毗邻区作为保护地本体的生态屏障和战略缓冲，土地利用类型以林地、耕地为主，占99%以上，草地、水域、建设用地较少。2000 - 2020年间，耕地、林地、建设用地之间的转换最为强烈，在自然因素和社会经济因素的多重驱动下建设用地占用大量林地和耕地，2020年面积增长至2000年的两倍。在自然发展、城镇发展、生态保护三种情境下，2040年金佛山自然保护地毗邻区土地利用变化均为林地、耕地面积减少，建设用地增加。其中，生态保护情景限制林地、草地、水域转化为其他用地，降低了建设用地扩张概率，但仍由保护地毗邻区逐渐向本体边界处蔓延，呈破坏性土地利用趋势，说明此情景下生态保护力度不足以遏制建设用地向本体扩张，不利于当前及长远自然保护地的发展。

PLUS模型在金佛山自然保护地毗邻区模拟的总体精度为0.95，Kappa系数为0.92，远高于0.75，对同类型自然保护地毗邻区土地利用的模拟研究提供了有效案例。但仍有不足之处，PLUS模型在参数设置时根据已有研究经验经过多次调试后设定，带有一定主观性；驱动因子的选择根据代表性和数据的可获得性而定，难以解释复杂的土地利用变化原因，在未来的研究中应探寻更加科学客观的方法。

自然保护地毗邻区的破坏性土地利用将切断自然保护地本体与毗邻区之间生态过程的连续性，加剧自然保护地本体孤岛化和毗邻区破碎化的风险。仅对毗邻区生态类用地进行保护限制不足以遏制建设用地向本体扩张，亟需对毗邻区人类活动进行更为精细化的管控。PLUS模型在金佛山自然保护地毗邻区的适用性较强，可以运用于同类型自然保护地毗邻区土地利用的模拟研究。

注：文中图表均由作者自绘。