基于PLUS模型的山东省土地利用变化模拟分析

王洁宁，王文超，海蒙蒙

（山东建筑大学，山东 济南 250101）

面对当前日益加剧的人口-资源-环境问题，全球变化研究成为近年来国际上最为活跃的研究领域之一。在众多的全球变化问题中，土地利用/土地覆盖（The land-use and land-cover,LULC）变化研究显得尤为重要[1]。LULC变化是一个复杂的现象与过程，是人地耦合系统的核心载体，土地以不同功能的空间状态为人类社会提供各类产品与服务[2]。如何科学地引导城市发展，协调保护生态环境与保障城市建设，构建平衡关系已经成为当前城市规划与建设中急需解决的问题[3]。

建立LULC模型是深入了解LULC及全球变化研究的主要方法[1]。在土地利用空间动态模拟中常用的方法有基于经验统计的模型方法、基于多智能主体分析方法和以元胞自动机模型（Cellular automata,CA）为代表的基于栅格邻域关系分析的方法等[4]。CA模型作为时空动态模型，拥有复杂系统空间变化能力，逐渐成功应用于生物演化、环境变化、景观更替、林火扩散和城市系统等的模拟研究并取得许多有意义的研究成果[5]。CA模型作为一种“自下而上”的地理模拟系统，以其强大的空间运算能力和内在的实用性、合理性，动态地反映出城市系统的复杂结构和土地利用空间格局，但同时存在主要着眼于元胞局部相互作用的局限性，进行预测区域土地利用变化总需求量需要结合其他模型等问题[6]。

因此，CA与多种模型结合，形成现实复杂动态系统模拟的综合CA模型，如：SLEUTH-CA[7]、多准则评价CA[8]、GeoSOS(Logistic-CA、ANN-CA、DT-CA)[9-11]等，广泛应用于自然和社会经济影响下的时空土地利用动态建模。为了更好地理解土地利用潜在变化之间的关系并改善斑块增长模拟，PLUS模型是基于CA发展起来的新型土地利用变化模拟模型，能够有效处理人与自然共同因素影响下的土地利用变化过程中的不确定性，其将“自上而下”的用地需求空间分配和“自下而上”的总适宜性和转换设置相结合，适用于大尺度过程与局部动力学相互作用的其他区域和过程[12]。该模型应用随机森林（Random forest,RF）分类算法，可以更好地挖掘各类土地利用变化的诱因，包含一种新的多类随机斑块种子机制，可以更好地模拟多类土地利用斑块级的变化。

本文基于梁迅等人提出的PLUS v1.3.5，以山东省为研究范围，对各类土地利用扩张背后的驱动力因素进行挖掘；
验证PLUS模型在模拟研究区未来土地利用模拟精度的同时，进行山东省未来土地利用多情景预测，为山东省国土空间规划和生态安全保护提供决策参考。

1.1 研究区概况

山东省位于中国东部沿海，黄河下游，陆域位于34°25′～38°23′N，114°48′～122°42′E，东西长721.03千米，南北长437.28千米。境内中部山地突起，西南西北低洼平坦，东部缓丘起伏，形成以山地丘陵为骨架、平原盆地交错环列其间的地形大势。泰山雄踞中部，主峰海拔1 532.7米，为全省最高点，黄河三角洲一般海拔2～10米，为全省陆地最低处。境内地貌复杂，大体可分为平原、台地、丘陵、山地等基本地貌类型（图1）[13]。

图1 研究区位置示意图

近年来山东省农业发展稳固增益、林牧渔业发展平稳，农业高质量发展基础增强，工业发展向好，建筑业发展壮大。山东省在快速发展中产生了一系列生态问题，例如森林、湿地生态系统和物种、遗传多样性受到威胁，生态功能下降，水源涵养能力下降以及水土流失的问题依然严重等[14]。情景模拟预测对于城市增长边界划定、生态安全格局、生态系统功能和经济可持续发展具有重要意义。

1.2 数据来源

所用的数据主要包括：2010年和2020年30米LULC数据，来源于GlobeLand30（全球地理信息公共产品，http://www.globallandcover.com/），将数据镶嵌、裁剪到研究区范围后根据研究需要，对原有9类数据重分类为7类（表1）。为减少数据冗余，提高运算速度，将30米单元栅格数据重采样为100米单元栅格数据。根据CA邻域因子转换、人文因素与自然因素共同影响下的土地利用变化规则，确定人文因素影响下的交通、经济及人口分布等和自然因素影响下的河流水网，土壤类型、高程等14类驱动因子（表2）。

表1 基于GlobeLand30类型赋值的地表重分类

表2 驱动因子数据来源信息表

2.1 土地利用转移矩阵

在LULC变化研究中，掌握变化的数量和幅度，同时对土地利用类型间转移情况的了解十分必要[15]。转移矩阵能够较好地展示不同土地利用类型之间的转变情况[16]，将土地利用变化的类型转移面积按矩阵的形式加以列出，作为用地结构与变化方向分析的基础，可细致反应各类用地之间的转化关系，进而了解转移前后各地类的结构特征[17]。数学形式如公式（1）所示：

式中，S(t+1)、St分别为t+1、t时刻的土地利用状态；
Pij为状态转移矩阵；
n则为土地利用类型。

2.2 基于斑块生成的土地利用模拟模型（PLUS）

PLUS模型包含两个模块：（1）用地扩张策略分析（LEAS）；
（2）基 于 多 类 随 机 斑 块 种 子 的CA模 型（CARS）。框架一般结构如图2所示。

图2 斑块生成土地利用变化模拟模型

2.2.1 用地扩张策略分析（LEAS）。LEAS提取两个时期土地利用变化间各类用地扩张的部分并叠加多重驱动因素进行采样从而获得每种土地利用类型和驱动因素的增长样本，选取采样点后经过RF进行各类型土地用途扩张因素贡献和每类土地生长概率计算。LEAS在进行土地利用分析时，分析单类土地利用而忽略其他，进而得到所有土地利用类型的过渡规则，避免了对随着类别数量指数增长的转化类型进行分析，简化了土地利用变化的分析过程，而且保留了模型在一定时间段分析土地利用变化机理的能力，具有更好的解释性。

LEAS策略通过提取两期土地利用变化间各类用地扩张的部分，并从增加部分中采样。采用随机森林算法逐一对各类土地利用的驱动因素进行挖掘，获取各类用地的发展概率及该类驱动因素对该时段各类用地扩张的贡献。随机森林算法从原始数据集当中抽取随机样本，并最终确定k类土地利用类型在单元格i上出现的概率，其表达公式如下：

式中，d的取值为0或1，若d=1，表示有其他土地利用类型转变为k类土地利用类型；
当d=0，表示土地利用类型转变成了除k以外的其他土地利用类型。x是由若干驱动力因子组成的向量；
函数I是决策树集的指示函数；
hn(x)是向量x的第n个决策树的预测类型；
M为决策树的总量。

2.2.2 基于多类随机斑块种子的CA模型（CARS）。CARS模块是一个CA模型，包括基于多种土地利用随机种子的补丁生成机制。CARS模型是一种情景驱动的土地利用模拟模型，它综合了“自上而下”（土地利用需求）和“自下而上”（土地利用竞争）效应；
在“自下而上”效应中综合了邻域权重、转换成本矩阵和递减阈值等，在模拟过程中，土地利用需求通过一个自适应系数影响当地土地利用竞争，驱动土地利用量达到未来需求。

2.3 利用Markov模型生成土地需求预测

Markov模型[18]是一种基于马尔科夫过程理论预测事件发生概率的方法。它具有无后效的特点，即当前状态只与前时刻的状态有关，与其他因素无关。土地利用状况可以通过公式（2）[19]进行预测。

在Markov模型预测过程中，转移概率矩阵p必须逐期保持不变，由于Markov模型缺乏弹性空间，不适合对未来土地利用数量进行长期预测，为了避免较大误差的出现，所以本研究只对未来10年土地利用结构数量进行测算[20]。

2.4 PLUS模型与Markov模型耦合机理

Markov模型虽然可以有效地反应LULC时间变化，但不能有效表征数据的空间属性，也不能显示这些要素的空间变化分布。在进行土地利用变化的空间模拟时，往往需要先预测研究区的土地利用变化总量和各类土地利用结构从而进行“自上而下”地控制，因此，Markov模型和PLUS模型耦合可以更加精准地模拟和预测研究区域未来LULC的空间分布。

3.1 土地利用变化分析

3.1.1 各类型土地数量变化。通过对2010年与2020年山东省土地利用类型统计分析（图3），发现十年间耕地面积减少，减少量为10 283.01平方千米，人造地表面积增加，增加量为9 958.07平方千米。森林、草地灌木略有增加，增加量分别为412.83平方千米和349.56平方千米。湿地水体、裸地与海域略有减少，减少量分别为228.47平方千米、70.31平方千米、138.61平方千米。

图3 2010—2020年各类土地变化量（平方千米）

3.1.2 土地利用转移矩阵。利用ArcGIS对2010年与2020年2期土地利用数据进行处理分析并制作土地利用转移矩阵分析（见表3、图4）。土地利用转移矩阵中，耕地转人造地表总转移量为11 595.37平方千米，占转移总量的7.45%，为7类土地利用类型中转移变化量最大的类型，产生的原因是城市人口的增长、经济发展、城市环境的改善，城市建设用地不断对外扩张，侵占耕地[21-22]。其次为人造地表转为耕地，转移量为2 028.9平方千米，占转移总量的1.30%，产生的原因可能是山东省合村并居政策[23]，随着农民向城市地区单向流动引发农村人口逐渐减少，大量的“空心村”出现，政府将空心村合并，节约土地资源，让更多人造地表转变为耕地成为可能。再次，转移量较大的是湿地水体与耕地的相互转换，转移量分别为1 529.89平方千米（0.98%）、1 102.26平方千米（0.71%）。

图4 2010—2020年土地利用转移矩阵

表3 2000—2020年土地利用转移矩阵（平方千米）

从空间结构来看，城市发展的总体增长格局仍然以原有城区为核心，加之交通干线的驱动作用，延续当前城市发展趋势，城市周边非建设用地不断被城市外延建设用地侵占，快速扩张（图5）。受DEM高程和《山东省生态保护红线规划（2016—2020）》保护的共同作用，鲁中山区和胶东半岛山区草地灌木和森林转换并未受到过大影响，转移量均保持在0.27%之内。

图5 现状土地利用类型

3.2 斑块生成土地利用变化模拟（PLUS）

首先，将2010年和2020年重分类后的土地利用数据和14类驱动因子进行LEAS用地扩张策略分析。采样方式为随机采样，决策树的数目20，采样率0.1，训练RF的特征个数为14，并行线程数量设置为8，生成驱动因子贡献和各类用地发展概率，根据驱动因子贡献值进行用地扩张因子贡献值分析（图6）。

图6 各类用地贡献值分析

DEM和土壤类型对耕地的发展具有最强的影响作用，耕地受其他影响因素较为均衡显著。森林、草地灌木裸地受到DEM和坡度强烈影响。水域对湿地水体和人造地表呈现最强烈贡献，距高速距离对海域的发展具有重要关系，距铁路距离等其他驱动因子对各类土地利用类型具有较为均衡的作用。

其次，进行基于多类随机斑块种子的CA模型（CARS）运行。以2010年土地利用数据作为初始数据，邻域范围，默认值为3，并行线程数量设置为8，递减阈值的衰减系数、随机斑块种子的概率和随机种子最大比例在多次测试的条件下依次设置为0.9、0.1和0.000 1；
在土地利用需求界面输入2020年的真实土地利用数据需求量。邻域权重参数范围为0～1，越接近1代表该用地类型的扩张能力越强，依据研究区实际土地利用类型的转换量，将用地类型的扩张能力依次定义为人造地表＞森林＞草地灌木＞裸地＞海域＞湿地水体＞耕地。人造地表扩张能力最强，耕地扩张能力最弱，分别设置为1和0.01，按照2020年实际各类型用地扩张面积与减少面积来定义其他5类用地（表4）。转移矩阵定义：裸地可以转换为其他任意用地类型的条件下，其他用地类型转换方式按照3.1.2中土地利用转移矩阵进行简化设置。其中，0表示不允许发生转换，1表示可以发生转化（表5）。

表4 依据各类型用地扩张量确定邻域权重

表5 各类用地转移矩阵

3.3 模拟结果评价

对比两期的土地利用可知（表6），2020年的模拟结果与2020年实际现状相比总体上没有出现大幅度波动，土地利用结构较为稳定。具体到地类，相较2020年现状土地利用，湿地水体、裸地和海域像元有所增加，森林和草地灌木像元数量有所减少，变化幅度保持在15%以内。将模拟结果（图7）进行模拟精度验证，采用Kappa系数和FoM（Figure of Merit）指数进行分别验证评价。Kappa系数是一种较好地用来描述栅格数据集随时间变化以及变化方向的方法，常用于检查影像分类对于真实地物判断的正确性程度[24]。在分类精度评价中，不同的精度评价方法有不同的精度划分标准和含义，本文借用Cohen提出的Kappa系数分类评价标准，见表7[25]。FoM定量地检验元胞尺度上的模拟精度，值越大，用于城市土地利用发展模型的尺度模拟结果精度越高，其中值一般在0.01～0.25之间[26]，FoM的计算由公式（5）表示为：

图7 2020年模拟结果

表6 研究区2020年现状及2020年模拟结果的土地利用结构像元数

式中，A=变但变错；
B=变且变对；
C=预测变但实际没变；
D=预测不变但实际变。

2020年模拟土地变化的总体精度为0.861，Kappa=0.687，在评价中具有显著的一致性，各类用地类型“生产者准确度”分别为：耕地=0.911、森林=0.786、草地灌木=0.634、湿地水体=0.900、人造地表=0.717、裸地=0.425、海域=0.931；
在FoM中，A=1 529 275，B=536 193，C=5 301，D=634 099，FoM=0.198。表明在该段时间内利用PLUS模型模拟的各类用地发展概率达到比较理想的发展结果，各类用地发展概率具有较为显著的说明性。

依据2010年与2020年研究区土地利用数据，通过Markov模型进行2030年各类土地需求量预测，结果见表8。结合2010—2020年各类用地发展概率，通过CARS对2030年研究区土地利用空间布局进行自由发展情景预测。其次，在预测时加入相关规划约束因子——生态红线来控制研究区未来发展趋势和方向，得到研究区2030年生态红线保护规划政策导向下的土地利用空间分布。

表8 基于Markov模型土地利用需求预测（像元个数）

4.1 自由发展情景

如图8（2）所示，在自由发展情境下，到2030年前后，随着人造地表剧烈扩张，耕地面积将大幅度减少，减少量将达到2020的10%左右。森林、草地灌木、裸地和海域小幅度减少，湿地水体小幅度增加。在空间发展上，各市以建成区为中心向周围继续扩张与发展，城市周边非人造地表插针式建设，小范围人造地表聚集扩大，耕地面积持续减少，森林、草地灌木组成的绿地将进一步缩小。绿地受到DEM强烈影响，故2030年预测发展中，鲁中山区以及鲁东半岛山区绿地变化并未受到强烈影响。城市人口的增加是城市土地扩张的直接动因，城市人口的增加又表现为城市居民对交通等方面的需求加强，进而影响各组团人造地表之间的连接关系不断清晰[21]。根据2010年、2020年现状发展以及2030年预测结果发现，受南部山区、北部黄河等生态敏感区的制约，济南城市空间长期以东西条带状发展为主[27]，通过章丘与淄博相连接并成为一个东西带状组团整体。鲁东南地区地势平坦，各市县将继续呈现“摊大饼”式发展，交通可达性对城市发展具有重要作用，调节城市发展速度[28]。城市间交通可达性提升会促进城市全要素生产效率提升，道路等交通基础设施改善带来的交通可达性提升会促进要素在区域间的流动，优化资源配置促进经济高质量发展[29]。

图8 多情景预测下土地利用发展方式

4.2 生态保护规划情景

生态保护红线是指依法在重点生态功能区、生态环境敏感区和脆弱区等区域划定的严格管控边界，是国家和区域生态安全的底线，对于维护生态安全格局、保障生态系统功能、支撑经济社会可持续发展具有重要作用[14]。

在生态保护的目标下，研究中引入了山东省生态保护红线（图9）作为限制性条件，得到生态保护红线控制情景下的2030年土地利用预测模式，见图8（3）。在该政策影响下，胶东半岛、鲁中南山地、黄河三角洲等地的森林、草地灌木等土地利用得到了较好的保护。非生态保护区除受到强烈驱动因子影响禁止发展位置之外，各城市以该城市原有核心区域向周边辐射发展，各城市之间联系不断增强，尤其表现在鲁中南山地北侧济南与淄博之间。该情景将生态敏感/脆弱、禁止开发区的人为干扰降到最低，为自然情境演替创造了条件，对保护生物多样性、水源涵养以及土壤保持、防风固沙具有重要意义。

图9 山东省2016—2020年生态保护红线

本研究基于对山东省省域历史时期（2010—2020年）土地利用模拟及2030年预测，结论如下：

（1）山东省2010—2020年土地利用转移量最大的是耕地转人造地表与人造地表转耕地，分别占转移总量的7.45%与1.30%。

（2）土壤类型对耕地的发展具有最强的影响作用，受其他因素影响较为均衡显著。森林、草地灌木和裸地受到DEM强烈影响，水域对湿地水体和人造地表呈现最强烈贡献。GDP与海域的土地利用发展具有重要关系。铁路对各类土地利用类型具有较为均衡的作用。

（3）在2020年PLUS模型模拟与验证中，Kappa=0.687，FoM=0.198，具有显著一致性。Markov模型结合PLUS模型进行2030年自由发展与生态保护情景模拟中，说明规律挖掘和模拟结果显著性强、精度高。

（4）在自由发展情景下，耕地减少，城市呈现“摊大饼”式发展，小面积人造地表聚集扩大，人造地表之间联系增强；
在生态保护情景下，红线范围内绿地得到较好的保护，生态敏感区/脆弱区与禁止开发区的人为干扰影响降到最低。

在智慧城市建设和科学智能决策背景下，PLUS模型对山东省土地利用变化模拟的适用性较好，多情景模拟结果可以为山东省国土空间规划和生态安全保护提供决策参考。今后可尝试在PLUS模型中引入“三生”空间约束条件、生态正义约束条件等进一步探讨国土空间发展的多种可能性。