土地利用规划与管理的目标是保障社会、经济和生态的协同可持续发展, 即在保障社会和经济发展的同时, 不影响生态系统服务的可持续供给^[1]。越来越多的研究证明, 土地利用对生态系统作用的影响不容忽视, 在进行土地利用规划时应当考虑到生态系统的影响。然而, 生态系统服务价值具有测量的困难性、驱动因素的复合性和评价标准的片面性, 导致生态系统服务的相关研究成果难以应用到土地利用规划管理当中^[2], 如何将生态系统服务充分融入到土地利用规划、管理和决策, 还存在着许多尚未解决的问题。在规划方面, 最主要的问题是如何决定众多不同土地利用类型的最佳分配和管理^[3]。De Groot等^[4]指出, 生态系统服务与土地利用规划、管理和设计中的关键问题有：(1) 土地与生态系统特征及其相关功能和服务之间关系的量化；(2) 土地功能和生态系统服务的空间定义(映射)和可视化；(3) 土地和生态系统服务的经济和社会价值保持一致和可比性的方法；(4) 在体现成本效益问题上适当考虑生态系统服务和利益相关者(在时间和空间上)价值变化的所有成本和收益(包括生态、社会和经济方面)。

虽然诸多的土地利用模拟方法为土地利用决策提供了有力而高效的工具, 但自上而下和自下而上的土地利用模型都无法完全解释跨多个层次的人类社会与生态环境相互作用的复杂性, 包括空间、时间和社会维度上的相互权衡^[5]。例如, 系统动力学模型在表示社会生态动力学和时间权衡方面很强, 但在解决社会生态异质性和空间权衡方面却相当弱。空间优化模型能够检测空间权衡, 但预测时间权衡的潜力较小。空间贝叶斯网络模型非常适合应对社会生态相互关系的不确定性, 但理解或表示复杂权衡的能力很差, 单独使用时不包含反馈回路^[6]。这些传统的空间建模方法没有明确代表单个参与者, 因此原则上无法捕捉异质人类社区和景观中的权衡。在土地利用管理中, 不同生态系统服务之间往往存在权衡或协同关系。例如, 增加林地有助于提升碳固存和水源涵养功能, 而扩展农田则可能导致水源涵养功能的下降。有助于在土地利用规划中进行更明智的取舍和平衡, 并且能够同时评估多个生态系统服务的价值及其之间的相互关系, 提供对生态系统功能的全面理解。识别和增强土地利用中的协同效应, 通过优化配置不同的土地利用类型, 提高整体生态服务效益^[7]。

社会与生态系统之间的平衡一直是可持续发展中面临的关键问题^[8]。黄土高原的可持续发展对中国整个生态安全和自然资源供应具有重要意义^[9—10]。研究结果表明, 黄土高原的生态修复和环境保护工作取得了显著成果^[11], 其中, 草地向林地的转化是改善地区生态环境质量的主要因素之一^[12]。《黄河流域生态保护和高质量发展规划纲要》指出, 黄土高原地区生态治理任务的重点是水土保持。然而, 土壤保持服务依赖植被恢复, 自然植被的蒸腾作用和水源涵养作用又会导致区域供水量的减少^[13], 对用水安全造成威胁, 影响经济和社会的发展。因此, 明晰黄土高原地区供水服务和土壤保持服务的权衡和协同作用, 以及这两种生态系统服务受土地利用的影响, 对该地区生态与社会的协调发展至关重要, 未来黄土高原地区的土地利用规划, 需要将这两个方面作为重点考虑。

针对上述问题, 本研究旨在建立一种基于土地利用结构的生态系统服务权衡协同的量化方法, 研究土地利用结构对水土生态系统服务的非线性影响。将生态系统服务和土地利用之间的非线性关系应用于土地利用规划模拟。根据土地利用结构对水土生态系统服务的协同区间确定合适的土地利用占比, 将约束效应的研究成果直接体现在约束条件的设定上, 进行多情景下的土地利用模拟。

1 材料与方法 1.1 研究区概况

黄土高原位于32°N—41°N, 107°E—114°E之间, 夏秋季炎热多暴雨, 全年≥10℃的积温2300—4500 ℃, 无霜期120—250 d, 日照时数1900—3200 h, 均比同纬度的华北平原高, 是我国辐射能高值区之一。按县级行政区界划分, 黄土高原地区总面积为64.87万km², 占全国土地总面积的6.8% (图 1)。

阳泉市全市面积4559 km², 地理坐标为37°40′′N—38°31′N, 112°54′E—114°04′E, 位于中国山西省东部, 黄土高原东部边缘和太行山脉西侧的交汇地带, 黄土的堆积从东到西逐渐增加, 海拔西高东低、北高南低, 高差接近1700 m。阳泉市地形以山地为主, 占总面积的80%左右(图 1)。

1.2 数据获取与处理 1.2.1 自然地理数据

(1) 1990—2020年降水1 km栅格数据、气温1 km栅格数据, 来自国家青藏高原科学数据中心(https://data.tpdc.ac.cn/)；(2) 1990—2020年潜在蒸散发1 km栅格数据, 来自国家青藏高原科学数据中心(https://data.tpdc.ac.cn/)；(3) 植被根系限制层栅格250 m数据, 植物有效含水量栅格250 m数据, 来自世界土壤数据库(HWSD, https://data.isric.org/)；(4) 数字高程模型(DEM), 来自地理数据云(https://essd.copernicus.org/)；(5) 土壤可蚀性因子1 km栅格数据, 来自国家青藏高原科学数据中心(https://data.tpdc.ac.cn/)；(6) 土壤类型1 km栅格数据, 来自资源环境科学与数据中心(https://www.resdc.cn/)；(7) 1985—2015年全国初级净生产力1 km数据, 来自全球变化科学研究数据出版系统(https://www.geodoi.ac.cn/)；(8) 1990—2022年黄土高原水文站点年径流和年输沙量数据, 来自国家地球系统科学数据中心(http://www.geodata.cn/)；(9) 1990—2020年的降雨侵蚀力采用李静等^[14]开发的黄土高原年降雨侵蚀力的二次多项式拟合公式。

1.2.2 社会经济数据

(1) 1990—2020年土地利用30 m栅格数据, 来自Earth System Science Data (https://essd.copernicus.org/)；(2) 2010年人口和GDP1 km栅格数据, 来自资源环境科学与数据中心(https://www.resdc.cn/)；(3) 一级道路、二级道路、高铁站位置和政府位置矢量数据, 来自Open Street Map (https://www.openstreetmap.org/)。

1.3 研究方法 1.3.1 K-means聚类

K-means聚类算法是一种无监督学习方法, 基本思想是以迭代的方式, 将给定的样本集划分为K个簇, 令每个簇中的数据对象与其簇中心尽可能相似^[15]。对于给定的一个包含n个d维数据点的数据集X={x₁, x₂, …, x_i, …, x_n}, 其中x_i∈R^d, 以及要生成的数据子集的数目K, 将数据组织为K个划分：C={c_k, i=1, 2, …, K}。每个划分代表一个类c_k, 每个类有一个中心μ_i。一般选取欧氏距离作为衡量数据对象间相似度的指标, 计算该类内各点到聚类中心μ_i的距离平方和, 以及各类的距离平方总和。

1.3.2 InVEST模型

InVEST模型基于生产函数来量化和评估生态系统服务, 包括陆地、淡水和海洋三类生态系统服务评估模型。该模型在空间上使用地图作为信息源, 并生成地图输出, 主要包括碳储存和固定、蓝碳储存和固定、水产量、养分持留、沉积物持留、异花传粉等服务的评估。本研究将InVEST模型用于黄土高原产水量和土壤保持的模拟, 主要使用以下两个模块。

(1) 年产水量模块

InVEST模型产水量模块(Water Yield, WY), 水循环被简化为降水、蒸散发和地下水入渗, 忽略了侧向入流。该模块基于Budyko水热耦合平衡假设和年平均降水量数据。确定研究区每个栅格单元x的年产水量Y(x), 公式如下：

(1)

式中, AET(x)表示栅格单元x的年实际蒸散量、P(x)表示栅格单元x的年降水量。水量平衡公式中, 土地利用/覆被类型的植被蒸散发计算, 采用Budyko水热耦合平衡假设公式计算。

(2) 沉积物输送模块

对于每个栅格, 使用模型先计算该栅格的年度土壤流失量, 然后计算沉积物输送量, 即实际到达溪流的土壤流失比例。一旦沉积物到达溪流, 模型将假设它们被输送到集水口。因此, 该模型无法模拟出可能影响沉积物负荷量的溪流过程。

年土壤流失量和土壤保持量由修订的通用土壤流失方程给出：

(2)

(3)

式中, SC是土壤保持量, USLE是实际土壤损失(t hm² a^-1), R是降雨侵蚀力(MJ mm hm^-2 h^-1 a^-1), K是土壤可蚀性因子(t h MJ^-1 mm^-1), L为斜率长度因子, S是斜坡陡度系数, C为植被覆盖管理因子, 表示土壤和作物管理系统在控制土壤流失方面的相对有效性, 是在特定作物和管理制度下土地的土壤流失量与连续休耕和耕作土地相应流失量的比值。P是保护支持实践因素, 反映了减少径流的数量和速率从而减少侵蚀量, 表示为支护实践条件下与斜坡上执行耕作条件下的土壤流失之比。

1.3.3 改进的约束线法

相比于传统的线性回归和相关分析方法, 约束线方法能够更好地刻画多因素影响的复杂生态系统中限制变量对响应变量的限制作用^[16], 绘制方法主要有参数法、散点云网格法、分位数回归法和分位数分割法^[17]。约束线上的点代表两个生态变量受到其他变量的影响最小^[18], 约束线的单调变化是阈值两侧主导因素不同的结果^[19], 线型上的波动则是多个主导因素交替的作用。

DBSCAN (Density-Based Spatial Clustering of Applications with Noise)是一种基于密度的空间聚类算法, 它在数据挖掘和机器学习等领域中得到广泛应用^[20]。该算法的核心思想是通过检查数据点在给定半径ε内包含点的数量来形成簇, 如果一个点的邻域内包含足够多的点(大于等于所设置的MinPts), 则这些点被视为一个簇的一部分, 并继续扩展该簇, 直到没有新的点满足条件。DBSCAN聚类可以在带有“噪声”的空间数据库中发现任意形状的聚类^[21], 能够识别并排除被视为噪声的低密度点。

在数据量较少的情况下, 本研究使用DBSCAN聚类对数据归一化后的散点进行异常值的判定。然后使用散点云网格法(dividing scatter cloud into size classes), 将散点云的限制变量与响应变量划分为若干区间, 根据区间最大值来确定边界点的位置, 最后用多项式函数对边界点拟合得到约束线(图 2)^[22]。根据约束线的变化趋势判断线型, 在表述时选择最为贴近的形状, 变量之间的约束线形状涉及到了线性线、凸曲线、驼峰曲线等。除了Hao等^[23]定义的12种约束曲线的线型外, 本研究补充了新的约束线线型：半凹波曲线, 即不完整的凹波曲线(图 3)。

1.3.4 多目标规划

多目标规划(MOP)的核心概念是在给定区域上需要同时考虑多个目标函数的最优化, 目标函数之间通常存在矛盾关系, 需要进行权衡和折中, 使各个子目标都尽可能达到最优^[24]。因此, MOP的解通常是一组妥协解, 没有任何解能在不损害其他目标的情况下改进任何一个目标, 即Pareto最优解集。

1.3.5 PLUS模型

PLUS模型(Patch-level Land Use Simulation Model)有两个组成部分：基于土地扩张分析策略(Land Expansion Analysis Strategy, LEAS)的规则挖掘框架和基于多类随机斑块种子的元胞自动机模型(CA based on multi-type random patch seeds, CARS)。

LEAS提取两期土地利用变化间各类用地扩张的部分, 并从中进行一定比例的随机采样, 根据土地利用类型划分子集。采用随机森林分类(random forest classification, RFC)算法逐一对各类土地利用扩张和驱动力因素进行挖掘, 通过对历史土地利用数据和相关驱动因素的分析, 获取各类用地的发展概率和驱动因素对该时段各类用地扩张的贡献。RFC算法能够处理高维数据, 可以根据随机噪声引起的袋外误差(out-of-bag error)测量出自变量对因变量变化的重要性, 并处理变量之间的多重共线性, 最终输出土地利用类型k在像元i上的增长概率P_{i, k}^d。

(4)

式中, d值为0或1；1表示有其他土地利用类型转变为土地利用类型k, 0表示其他的转变；x是由多个驱动因素组成的向量；I是决策树集的指示函数；h_n(x)是第n个决策树对向量x的预测类型；M是决策树的总数。

CARS模型能够基于预期的土地利用需求、土地利用转移矩阵、土地利用类型扩张权重, 结合随机斑块种子生成和阈值递减机制, 通过自适应系数使土地利用需求影响局部土地利用竞争, 在发展概率的约束下模拟土地利用斑块的动态演化过程, 从而预测未来土地利用的变化情况。

1.3.6 多情景土地利用模拟

对阳泉市的土地利用进行提取, 该地区的土地利用类型有耕地、林地、灌木、草地、水域和城建用地, 依次对应决策变量x₁-x₆ (表 1)。这些变量将作为优化对象, 通过设置目标函数和约束条件, 使用遗传算法NSGA-Ⅱ求解多目标规划模型, 以获得水土生态系统协同管理和可持续发展的土地利用占比, 为PLUS模型模拟土地利用分布提供初始数据。

一、目标函数

基准情景：基于马尔科夫链模型预测2050年的土地利用变化趋势。

对照情景：未充分考虑地区土地利用结构对水土生态系统服务特点影响的情景。追求经济和生态的双重目标, 但不考虑水土生态系统服务之间的协同作用。目标函数为F₁和F₂(表 2)。

协同情景：水土生态系统服务协同的土地利用规划, 目标函数为F₁, F₂, F₃, F₄(表 2)。

经济效益通过GDP栅格数据依照土地利用分类进行赋值, 反映不同土地利用类型的经济价值(表 3)。生态系统生态服务价值系数依据谢高地等^[25]的生态服务价值当量因子表进行计算。为了与阳泉市的土地利用类型相匹配, 本研究参考了张瑜等^[26]的研究成果, 将其折合为一级分类(表 4), 再求和得到生态系统生态服务价值系数。由于协同情景已经在目标函数层面考虑了水供给能力函数F₃(x)和土壤保持能力函数F₄(x), 因此在协同情景的生态服务价值系数里, 剔除了土壤保持服务系数和水资源供给服务系数, 以保证协同情景和对照情景之间在方法理论上的对比更加确切。

二、约束条件

基于阳泉市2050年人口和粮食单产预测, 设置约束方程(表 5)。求解时决策变量满足下列要求：(1) 所有土地利用类型之和为阳泉市的总面积；(2) 阳泉市土地利用类型能满足未来人口发展的需要；(3) 阳泉市耕地面积能实现90%的粮食自给率, 以保障粮食安全。复种指数MCI参照全国复种指数分布图, 采用阳泉市内的最高值, 取值为2。

2 结果与讨论 2.1 黄土高原土地利用结构聚类及分析

黄土高原包含334个县级行政区域, 面积20.6 km²—33893.4 km²不等, 平均面积为2957 km²。本研究采取和县域类似大小的空间尺度, 使用ArcGIS10.7对研究区的DEM进行水文分析, 划分了226个小流域, 平均面积为4078.4 km²。基于K-means聚类划分了六种土地利用结构类型(图 4), 分别农地结构、林地结构、农地-草地结构(以下简称“农草地结构”)、草地结构、以耕地和自然植被(包括林地和草地)为主导的混合结构(农地-林地-草地结构, 以下简称“农林草地结构”)、裸地结构。

土地利用结构的转换以农草地结构向草地结构为主, 其次是农草地结构向农林草地结构。农地结构的转化方向为草地结构和农林草地结构；农草地结构大部分转化为农林草地结构, 其次是草地结构, 少部分转变为农地结构；草地结构的转变源自农地比例的增加, 少数的草地结构转化为了农地结构和农草地结构, 但绝大多数的草地结构草地面积呈增加趋势(图 4)；农林草地结构少部分转化为了林地结构, 没有转化的部分也能观察到林地比重的增长。数量上, 黄土高原土地利用结构的转换以农草地结构向草地结构为主, 其次是农草地结构向农林草地结构, 分别代表黄土高原西部的退耕还草和东部的退耕还林。其次是草地结构和农地结构向农草地结构转化, 分别代表黄土高原西部对草地的开垦, 与东部耕地的自然植被恢复(图 4)。

在对黄土高原的土地利用结构进行划分时, 应用K-means聚类算法得到了明确的土地利用类型。尽管阳泉市存在建设用地, 但由于建设用地在整体区域中的占比相对较小, 因此在划分中未予单独列出。这种处理方式在一定程度上简化了分析, 但可能忽略了建设用地对生态系统服务的潜在影响。与已有研究相比, 本研究的划分方法较为简洁, 但缺乏对建设用地的详细分析, 可能导致某些区域的土地利用类型未能全面反映。例如, 李松林和王萍^[27]在对类似区域的研究中, 强调了城市化进程对生态环境的影响, 并建议在土地利用分类中纳入建设用地以提升评估的精度。未来研究应考虑将建设用地纳入分类体系, 以提供更为全面和准确的土地利用结构分析。此外, 阳泉市的建设用地虽占比不大, 但其对当地生态系统服务的影响不容忽视。建设用地通常伴随人口密集和经济活动, 对水源供给、土壤保持等生态系统服务有显著影响。例如, 城市扩张可能导致水源涵养能力下降和土壤侵蚀增加。因此, 未来研究应细化土地利用类型, 特别是对建设用地进行详细分类和评估, 以更准确地反映其对生态系统服务的影响^[28]。

2.2 土地利用结构对水土生态系统服务的影响 2.2.1 土地利用结构对水土生态系统服务的约束效应

从约束线的形状和边界值判断, 土地利用结构对碳固存的约束并不如对水供给、土壤保持的明显。在6种土地利用结构中, 林地结构的3种生态系统服务(以下简称“服务”)都处于较高水平, 而裸地结构生态系统服务极低(图 5)。

在农地结构中, 耕地的占比对水供给和土壤保持约束线为驼峰曲线, 对碳固存为负凸曲线, 当耕地占比为0.6—0.7时, 能较好地均衡三种服务。林地结构对水供给、土壤保持和碳固存的约束线分别为U型、对数曲线和正线性线, 当林地占比约为0.5时, 能较好地均衡三种服务。森林占比较低时(0.5—0.6), 就已经表现出很好的土壤保持作用, 随着森林比例的增加, 碳固存的边界值增加, 土壤保持反而下降。林地占比小于0.8时, 对水供给的约束随森林比例的增加而增加。农草地结构以耕地和草地为主导, 耕地占比对水供给、土壤保持和碳固存的约束线型分别为近似驼峰型的凸波、下降的凸波、凸波型。草地占比对3种服务的约束线分别为近似负凸型、驼峰趋势的凸波和驼峰型。随着草地占比的增加, 水供给和碳固存的变化不明显, 而土壤保持呈现双峰的波动趋势, 总体上略高于单一的耕地结构。当耕地占比约为0.45、草地占比约为0.4时, 能较好地均衡三种服务。草地结构中草地的占比对水供给、土壤保持和碳固存的约束线分别为负凸型、驼峰型和波动不大的凸波型。在草地占比接近于1时观察到碳固存的下降, 当草地占比在0.75—0.85时, 能较好地均衡三种服务。随着耕地占比的增加, 农林草地结构对3种服务的约束线分别呈凸波型、下降的凸波、平坦的半凹波, 林地占比对水供给的约束线呈平坦的半凹波, 对土壤保持和碳固存则呈现正凸曲线, 草地占比对3种服务的约束线分别呈负线型、半凹波曲线和负凸曲线。农林草地结构的成分最为复杂, 约束因素较多, 当耕地、林地和草地占比在0.15、0.4、0.25—0.35时, 能较好地均衡三种服务。裸地结构中的水供给和土壤保持的约束呈驼峰型, 而碳固存接近于负线型。

2.2.2 土地利用结构对水土生态系统服务的权衡协同作用

各土地利用结构中生态系统服务之间的约束曲线基本上是驼峰型、凹波、半凹波和凸波(图 6)。水供给-土壤保持的约束关系在农地结构、农林草地结构中呈现明显的驼峰或凸波, 即超过某个阈值之后, 土壤保持随水供给的增加而增加转变为随着水供给的增加而减少。而裸地结构的土壤保持-水供给约束线的趋势接近于对数曲线, 即随着水供给的增加土壤保持减少, 并且超过某个阈值后, 趋于一个较为稳定的值。土壤保持-碳固存的关系在农地结构和草地结构中的约束线形状接近, 呈现不明显的负凸型。林地结构和农林草地结构的形状较为接近, 都为半凹波型, 即在中间呈现一个凹部。农草地结构中的土壤保持则随碳固存的增加呈现波动的下降趋势(凹波型)。碳固存和水供给的关系, 除农草地结构和草地结构外, 其他土地利用结构约束线都呈现半凹波型, 即先下降, 而后在某个阈值上水供给能达到最大值。而农草地结构和草地结构这两种土地利用结构在黄土高原地区占比较大, 影响约束关系的因素较多, 呈现下降的凹波。

本研究使用InVEST模型评估了黄土高原的生态系统服务功能, 重点关注了水供给和土壤保持两项主要服务。结果显示, 这两项服务在不同土地利用类型之间存在显著差异, 农田和林地在水供给和土壤保持方面表现突出。林子雁等^[29]也发现农田和林地在生态系统服务提供方面具有重要作用。然而, 本研究的评估中未包括碳固存能力, 这一方面的缺失可能导致对生态系统服务综合价值的低估。碳固存是生态系统服务的重要组成部分, 尤其在应对气候变化方面具有重要意义。未来研究应考虑将碳固存能力纳入评估, 以提供更为全面的生态系统服务分析。通过多目标规划法, 本研究分析了三类生态系统服务(水供给、土壤保持和生态系统服务价值)之间的权衡与协同关系。结果表明, 这些服务之间存在显著的协同关系, 特别是在土地利用类型和空间布局方面。优化结果显示, 多目标规划能够有效提升生态系统服务的整体价值。本研究中的多目标规划方法具有显著的优势, 但未能单独量化碳固存能力, 可能影响优化结果的全面性。孟欣宇和吴远翔^[30]提出综合指标的构建能够更好地反映不同生态系统服务之间的复杂关系, 从而为决策提供更为科学的依据。因此, 未来研究应考虑构建更为综合的生态系统服务指标, 包括碳固存能力在内, 以提升优化结果的准确性和实用性。

2.3 典型区多情景土地利用模拟

基于本研究对黄土高原地区土地利用结构的划分, 阳泉市的土地利用结构为农地-林地-草地混合结构, 当耕地、林地、草地的占比区间分别为0.1—0.2、0.32—0.42、0.25—0.35时, 水供给和土壤保持服务可以同时处于较高水平(图 7), 因此依据生态系统服务的约束, 设置协同情景的耕地x₁、林地x₂、草地x₄决策变量边界。对于非主导土地利用类型, 设置上下限为(1±0.2)为基准情景上下限(表 6)。

为了在方法上与协同情景形成对照, 尽可能控制变量, 确保对照情景与协同情景在决策变量的区间长度上保持一致, 区间中点参照2020年的土地利用和基准情景的土地利用取中间值(表 6)。

在确定不同土地利用类型的决策变量边界后, 采用遗传算法NSGA-Ⅱ来求解构建的多目标规划模型。NSGA-Ⅱ是一种高效的非支配排序遗传算法, 适用于处理多目标优化问题, 能够在搜索过程中找到一组非支配解, 即Pareto最优解集。由NSGA-Ⅱ求解了协同情景和对照情景的土地利用占比, 由马尔科夫链对基准情景的土地利用占比进行了预测(表 7)。

从未来的土地利用数量上看, 三种情景下的水域和灌木面积均处于较低水平。其中, 对照情景中的城建用地面积最大, 基准情景次之, 协同情景最小, 但三者均超过了2020年的水平(图 8)。这一趋势表明, 无论是基于需求的多目标优化, 还是依赖历史趋势的马尔科夫链, 阳泉市在未来30年内都倾向于增加城市建设用地的面积。在人口容纳方面, 基准情景、对照情景和协同情景分别可容纳167万人、196万人和145万人, 均能满足阳泉市的人口增长需求。从未来耕地需求来看, 三种情景的需求依次递减, 但均低于2020年阳泉市的耕地面积, 预示着阳泉市在未来将面临耕地数量持续减少的趋势。三种情景的粮食自给率均超过1.5, 在未来都能充分确保粮食安全。在基准情景中草地面积减少, 与另外两种情景相反, 这一差异的主要原因在于基准情景在规划时沿用了过去的土地利用变化趋势, 在未来计划中布置了更多的林地, 因此使得草地的面积相对减少。

从未来的土地分布格局来看, 阳泉市的土地利用呈现出以城市建设为中心, 向外依次为耕地、草地和林地的层次分布特征。该分布模式在三种不同情景下得到了明显体现, 三者最显著的差异在于城市建设的规模和布局。在对照情景中, 城市建设的分布相对更为分散, 从空间发展趋势来看, 城市建设倾向于扩展至原有的耕地和草地区域(例如, 图 9中红圈内部), 表明对照情景下城市建设与耕地、草地之间的竞争关系。而在基准情景中草地占据了最大的面积, 其扩张方向主要是向林地和耕地延伸(例如, 图 9中紫圈内部)。因此, 可以粗略地观察到土地利用的发展过程中各种土地利用类型之间的替代规律, 大致上遵循它们分布的层次顺序, 即城市建设倾向于取代耕地和草地, 草地则可能取代林地和耕地, 而林地则可能与草地发生替代关系(例如, 图 9中蓝圈内部), 这种替代关系体现了阳泉市土地利用变化的内在逻辑和动态平衡。

经济发展情景下, 建设用地显著增加, 而生态保护情景下, 林地和草地面积有所扩大, 这与经济发展和生态保护之间的权衡关系相一致。然而, 模拟结果中某些土地利用类型的变化趋势可能受到模型参数设定和输入数据的限制, 未来研究应进一步验证和调整模型以提高模拟精度。不同情景的模拟结果为政策制定提供了有价值的参考。例如, 在生态保护优先情景下, 林地和草地的增加有助于提升区域的生态系统服务功能, 如水源涵养和土壤保持。然而, 这一情景下的经济发展可能受到一定限制。因此, 政策制定者需要在生态保护和经济发展之间寻找平衡, 通过合理的土地利用规划, 实现生态与经济的协调发展。尽管多情景模拟提供了不同政策和管理措施下土地利用变化的可能性, 但仍存在一些局限性。首先, 模拟结果依赖输入数据的准确性和完整性, 特别是经济、社会和环境数据的精度对模拟结果有重要影响^[31]。其次, 模型参数的设定和情景设计需要更加精细和符合实际, 特别是对于政策调控的情景, 应考虑更加复杂的社会经济因素。此外, 模拟结果的空间分辨率和时间尺度也需要进一步优化, 以提高模拟结果的实用性和可靠性。本研究的多情景模拟结果与其他研究在总体趋势上具有一致性。例如, Ha等^[32]发现生态保护措施对提升区域生态系统服务具有显著效果, 而经济发展则会加剧土地资源的紧张。本研究在情景设计和模型应用上具有一定创新性, 如引入PLUS模型进行精细化模拟, 并结合多目标规划法进行优化分析。这为未来的研究提供了新的思路和方法。

3 结论

本研究以土地利用和水土生态系统服务为研究对象, 在小流域尺度上划分土地利用结构。通过InVEST模型计算了黄土高原的水供给和土壤保持服务。采用K-means聚类、改进的约束线法获得了土地利用结构对水土生态系统服务的约束效应以及它们之间的协同-权衡关系, 使用多目标优化结合PLUS模型将阳泉市作为典型区域规划2050年的土地利用, 验证了土地利用和水土生态系统服务的非线性关系在土地利用规划方面的应用成效。在小流域尺度上汇总土地利用结构和水土生态系统服务, 基于改进的约束线法评估了每种土地利用结构能兼顾水土生态系统服务的主导土地利用占比：农地结构、林地结构、农地-草地结构、草地结构、农地-林地-草地结构、裸地结构分别对应0.6—0.7, 0.5, 0.45/0.4 (农地/草地), 0.75—0.85, 0.15/0.4/0.25—0.35 (农地/林地/草地)。阳泉市未来的土地利用发展倾向于增加城市建设用地面积、减少耕地面积, 呈现以城市建设为中心, 向外依次为耕地、草地和林地的层次分布特征。