生态系统服务连接了社会和生态系统, 是人类直接或间接由生态系统的特征、功能或过程中收获的效益^[1—3]。生态系统服务的可持续供给关系着人类福祉与景观可持续性^[4—6]。而生态系统服务的可持续利用不可避免要面临不同生态系统服务的权衡, 因此评估生态系统服务的权衡关系是协调生态系统服务权衡以提高区域整体生态效益的先决条件。

土地利用作为人地交互的核心环节, 其不仅是生态系统服务的载体, 也是生态系统服务变化的重要驱动力之一^[7]。土地利用的决策、管理和优化是生态系统服务权衡评估的主要落脚点之一^[8]。2015年联合国可持续发展峰会通过了《2030年可持续发展议程》, 可持续发展目标(Sustainable Development Goals, SDGs)的范围涵盖社会、经济和环境三个维度^[9]。其中SDG1、2、3、5、9、11都直接或间接与土地利用有关, 可持续的土地利用方式在目标上与SDGs吻合, 其核心即是满足人类经济生活和环境需求并不断提升资源品质^[10—11]。因此, 如何通过优化土地利用方式, 进而保障生态系统服务的可持续供给, 是实现可持续发展的关键问题之一^[12]。

情景分析是从土地利用变化的驱动角度评估生态系统服务权衡的主要手段^[13], 在此基础上集成多目标优化算法得到的优化结果可以为土地利用政策提供参考以推进可持续发展。如Ramel等^[14]通过设置目标权重, 在山区识别了兼顾生物多样性保护热点及十项生态系统服务价值的生态保护情景。Kennedy等^[15]基于贪婪算法, 绘制了农产品供给服务、生物多样性保护和水质净化服务的效率前沿曲线, 找到了减少三者权衡效应的优化情景。Verhagen等^[16]使用遗传算法进行的土地利用优化, 在既保障区域农产品供给服务、景观美学价值和物种栖息地质量的同时又最小化乳业损失。多目标优化算法与土地利用变化模拟的耦合模型则可以基于优化后的各地类面积来模拟土地利用的未来情景, 进而对各生态系统服务的价值进行权衡^[17]。但是目前的研究多是关注优化后的景观配置, 很少在空间上分析土地利用的优化位置和生态系统服务不同权衡强度区域是否会表现出一致性, 也即是高生态系统服务权衡强度区域是否也具有很高的优化潜力。

云南省怒江流域面对严峻的可持续发展问题, 局部地区生态退化问题严重, 经济发展滞后^[18—19]。坡耕地资源是当地农户最主要的生计资本, 粗放的利用方式在提供产品供给服务的同时, 也对土壤保持服务造成了巨大损失。这两项生态系统服务利用的矛盾是制约区域可持续发展的重大阻碍之一。因此, 本文选取云南省流域中游的施甸县中粮食生产和土壤保持两项区域典型生态系统服务, 使用均方根偏差(Root Mean Square Error, RMSE)方法评估了2000—2020年每5年两项生态系统服务的权衡动态特征及分布模式。2020年施甸县森林覆盖率为55.34%, 未达到同期云南省森林覆盖率水平65.04%, 尚有提升空间。同时研究表明, 在生态保护情景下, 耕地向生态用地(如林地)转出的优化可以有效缓解除粮食供给外其他多个生态系统服务间的权衡效应^[20]。Peng等^[21]在云南省西北部地区开展的情景分析也证明坡度在20°以上的坡耕地退耕还林可以有效协调粮食生产和其他类型生态系统服务的权衡关系。基于此本研究以2020年作为基准年, 将决策变量设置为坡耕地向林地转出, 通过多目标线性规划绘制两项服务的效率前沿曲线, 评估优化情景下两项服务的供给特征, 继而识别高优先优化区域中各权衡强度等级坡耕地的分布情况。本文揭示了不同权衡强度坡耕地的优化潜力, 以期为边远山区坡耕地管理提供参考依据。

1 研究区和数据准备 1.1 研究区概况

施甸县地处云南省怒江流域中游, 隶属保山市, 下辖5镇8乡, 地处98°54′—99°21′E和24°16′—25°00′N。气候属中亚热带为主的低纬山地季风气候, 年平均气温17.6℃。地形属怒山尾翼山地峡谷区, 整体地势较为复杂, 海拔区间为480—2861m。植被类型包括半湿润常绿阔叶林、湿性常绿阔叶林、针阔混交林、高山灌丛草甸等^[22]。施甸县总国土面积为2009km²其中坝区面积为186km², 占总面积9.26%；山区面积为1517km², 占总面积75.51%；半山区面积为306km², 占总面积15.23%。截止2021年底施甸县总年平均人口3.48×10⁶人, 其中年平均乡村人口为3.23×10⁶人, 占比92.91%；规模以上工业总产值4.16×10⁶万元, 而农业总产值达5.18×10⁶万元。施甸全县都囊括在怒江流域内, 开展生态系统服务研究具有一定代表性。

1.2 数据采集

本文所用数据包括2000、2005、2010、2015、2020年5期施甸县土地利用/覆被图、数字高程模型(Digital Elevation Model, DEM)、气象数据、土壤质地类型数据、归一化植被指数(Normalized Difference Vegetation Index, NDVI)、社会经济数据。其中, 气象数据包括月平均温度、月平均降水；土壤数据包括砂粒、粉粒、黏粒和有机碳含量；社会经济数据为云南省、施甸县粮食产量, 各项数据来源见表 1。

2 研究方法 2.1 生态系统服务评估 2.1.1 土壤保持服务

土壤保持量通常由土壤潜在侵蚀量和实际侵蚀量确定^[23]。修正后的通用土壤流失方程模型(Revised Universal Soil Loss Equation, RUSLE)因结构简单、参数易得、计算便捷等优点在土壤侵蚀评估中得到广泛应用^[24], 本文采用该模型评估土壤保持服务, 如公式(1)所示：

(1)

式中, S_i是栅格i土壤保持模数(t hm^-2 a^-1)；R_i是降雨侵蚀因子(MJ mm hm^-2 h^-1 a^-1)；K_i是土壤可蚀性因子(t hm² h MJ^-1 hm^-2 mm^-1)；LS_i是坡长坡度因子；C_i为作物经营因子；P_i为水土保持措施因子。

降雨侵蚀因子R根据Wischmeier经验方程^[25]计算得到, 如公式(2)所示：

(2)

式中, R为降雨侵蚀因子(MJ mm hm^-2 h^-1 a^-1)；p_i为第i月平均降水(mm/m)；p为年平均降水(mm/a)。

土壤可蚀性因子K由侵蚀-生产力评价模型(Erosion-Productivity Impact Calculator, EPIC)^[26]计算得到, 如公式(3)所示：

(3)

式中, K为土壤可蚀性因子, SAN, SIL, CLA和C分别表示砂粒、粉粒、黏粒和有机碳的含量百分数。

坡长坡度因子LS的计算采用刘宝元等^[27]第一次全国水资源普查的方法, 作物经营因子C与水土保持措施因子P参考杨子生等^[28—30]2001年在云南省金沙江流域设置侵蚀试验小区的实测值。

2.1.2 粮食生产服务

植被指数是遥感作物监测的重要指标, 其中NDVI与作物产量具有显著线性相关关系^[31]。因此本文使用NDVI分配方法评估区域的粮食生产服务。云南省2000、2005、2010、2015、2020年粮食产量由统计年鉴提取, 然后根据逐栅格年最大NDVI与总年最大NDVI的比例分配至耕地中, 最后裁剪至研究区。粮食生产服务的计算见式(4)：

(4)

式中, GP_i是分配到栅格i的粮食产量, GP_sum是云南省年粮食总产量, NDVI_i是栅格i的NDVI值, NDVI_sum是云南省耕地NDVI的总值。

2.2 生态系统服务权衡识别

RMSE可以在空间上量化并可视化生态系统服务间权衡强度^[32]。RMSE的取值范围为0—1, 表征标准化后单个生态系统服务与生态系统服务期望之间的离散程度, RMSE越大则权衡强度越大^[33]。本文使用RMSE评估研究区两项生态系统服务的权衡强度, 并按照RMSE值0—0.2, 0.2—0.3, 0.3—0.4, 0.4—0.5, >0.5对权衡强度进行分级。其中RMSE＜0.3视为低权衡强度, RMSE在0.3—0.5视为中权衡强度, RMSE>0.5视为高权衡强度^[34]。RMSE的计算见式(5), 式(6):

(5)

(6)

式中, n是生态系统服务数量；ES_i是第i个生态系统服务的标准化值；ES为平均生态系统服务期望；W_i是各项生态系统服务的权重, 本文设定两项生态系统服务同样重要, 故W_i为0.5。

冷热点分析可以在空间上识别具有显著性的高/低值聚集区域^[35], 因此在RMSE结果的基础上再进行冷热点分析, 识别两项生态系统服务不同权衡强度区域在空间上的分布模式。

2.3 生态系统服务优化

为了实现两种服务权衡的优化, 通过线性规划算法来处理土壤保持服务和粮食生产服务的供给。本研究使用基于python的线性规划工具The Restoration Opportunities Optimization Tool (ROOT)^[36]进行优化操作, 其线性规划使用cvxpy库, 可以有效处理整数型决策变量的凸优化问题^[37]。线性规划在处理多目标问题时通过对各项目标设置不同的权重将多目标问题转化为单目标问题, 不同目标权重组合下目标函数的优化值构成了效率前沿, 而每次迭代计算中都容易被选中进行生态恢复的决策单元则可以被视为高优化潜力区域。首先生成空间决策单元, 在ArcGIS中利用泰森多边形工具产生覆盖研究区的六边形要素, 受限于计算机算力, 每个决策单元面积设置为15hm²。然后计算两项生态系统服务优化前后的差值输入决策单元进而构建目标函数。

目标函数为：

(7)

服从：

(8)

式中, w_k为第k项生态系统服务的权重值；F_ki为第k项生态系统服务在第i决策单元的值；x_i为决策变量, 取值为0或1, 表示决策单元中坡耕地是否转出。

本研究目标是探究高优化潜力区域中生态系统服务权衡强度的空间差异, 因此设定研究区所有坡耕地都能作为优化的潜在区域, 其总面积为13861.14hm²。

约束函数为：

(9)

式中, A_i为第i个决策单元中坡耕地的面积；x_i为决策变量。

3 结果与分析 3.1 生态系统服务权衡强度时空变化

研究区2000—2020年粮食生产与土壤保持服务的权衡强度时空分异明显(图 2), 各个年份的平均RMSE分别为0.466、0.482、0.483、0.473和0.499, 表现出逐渐增加的趋势。两项服务权衡高强度区长期集中在低海拔的中北部坝区部分, 这是由于坝区地势平缓, 坡长和坡度因子基本为0, 可以视为潜在土壤流失极小, 近似无土壤保持服务, 但是坝区内的水田具有高值的粮食生产服务, 因此产生了高强度的权衡。同时该区域权衡强度的空间动态变化明显, 从2015年开始权衡强度自南端出现衰减, 至2020年权衡强度衰减的部分扩散到中央位置。图 3为生态系统服务权衡强度的冷热点分析结果, 坝区及山区权衡强度的高值与低值空间分布均表现出显著的聚集特征, 其中大面积的聚集现象发生在研究区的中北部与东南缘的坝区, 小面积的聚集则零散分布在山区的耕地内。

2000—2020年施甸县RMSE大于0.4区域的面积远超0—0.4区域, 表明研究区20年间这两项生态系统服务在空间上以中、高权衡为主(表 2), 中权衡强度(RMSE: 0.3—0.5)区域面积在研究期间分别为8857.08、7440.57、8467.47、7621.92、5622.75hm², 除了2000年外均小于高权衡强度(RMSE>0.5)区域, 因此高权衡是研究期间两项生态系统服务关系的最主要特征。

3.2 优化情景下生态系统服务的权衡特征

研究区坡耕地两项生态系统服务在同时最大化的目标下产生的效率前沿曲线及优化情景如图 4所示。前沿曲线上的点集代表最优土地利用情景下, 也即研究区坡耕地转为林地后对应的两项生态系统服务的潜在权衡关系。其中点A和点E分别对应最小化粮食生产服务损失情景和最大化土壤保持收益情景, 点B、C、D则分别对应土壤保持服务提升至理论最大值的40%、50%和75%时的优化情景。当粮食生产目标的权重相对土壤保持目标的权重极小时, 需要以超过97%的坡耕地恢复为林地作为代价, 使得土壤保持服务的收益达到理论最大值26.69×10⁶t hm^-2 a^-1；而若优先保障粮食生产服务, 则土壤保持服务仅可以提升0.01×10⁶t hm^-2 a^-1；当土壤保持服务目标与粮食生产服务的目标权重分别为0.9和0.4时, 土壤保持服务将提升至13.35×10⁶t hm^-2 a^-1, 为理论最大值的一半, 此时粮食生产服务损失达到9.59×10⁶kg。各优化情景下坡耕地的转出面积可以视为优化的经济成本指标(表 3)。在前沿曲线上土壤保持服务收益由最低提升至理论最大值的50%, 仅需要恢复3388.51hm²坡耕地(点C), 而继续提升至75%则需要额外恢复3140.28hm²坡耕地(点D), 表明通过土地利用转换方式提升单一生态系统服务具有显著的阈值效应, 当土壤保持服务收益达到理论最大值50百分位数后的继续提升会显著增加成本。

3.3 高优先恢复区域生态系统服务权衡强度分布特征

在前沿曲线上选取的代表点对应的坡耕地转出图中可以看出, 随着多目标线性规划每次迭代, 总有一些决策单元始终被选为优化的对象, 这些优先选择的决策单元即为高优先恢复区域(图 5), 本研究识别的高优先恢复区域内坡耕地面积为2612.07hm², 占总坡耕地(潜在优化区域)面积的18.84%。通过对比高优先区域内坡耕地在潜在优化区域中RMSE的占比, 可以识别不同权衡强度等级坡耕地的优化潜力(表 4), 各权衡强度坡耕地的优化潜力排序为：中权衡强度>低权衡强度>高权衡强度。其中RMSE在0.4—0.5之间的坡耕地不仅具有最高的优化潜力, 同时面积占据了显著的优势, 说明这类坡耕地最适合进行生态恢复。而低权衡强度的坡耕地两项生态系统服务趋向于协同关系, 同时面积很小, 虽然有一定的优化潜力, 但是生态恢复的整体收益不高。高权衡强度的坡耕地则优化潜力不占优势, 只有在土壤保持目标的权重远超过粮食生产服务目标的权重时, 才应当考虑进行生态恢复。

4 讨论与结论 4.1 讨论

本研究中施甸县的粮食生产和土壤保持服务之间表现为权衡效应, 这与孙瑞等^[38], 余玉洋等^[39]在淮河流域和秦巴山区的研究结果一致。研究区的粮食生产与土壤保持服务的权衡强度在2000—2020年间表现出不断上升的趋势。一方面, 区域粮食产量的递增与这两项生态系统服务权衡的加剧具有一定的空间一致性^[40]；另一方面, 生态恢复工程的实施也可以缓解粮食生产服务与其他服务的权衡^[21]。2000—2020年间云南省总粮食产量提升428.06万余吨, 同时国家于2014年启动新一轮退耕还林工程, 这与研究区两项生态系统服务权衡强度的整体趋势与局部变化特征相吻合。

研究区两项生态系统服务的权衡特征表明进行土地利用优化的必要性, 本文通过线性规划识别的坡耕地生态恢复的高优先区域主要分布在研究区的边缘山区, 海拔相对较低(图 1), 平均坡度超过25.07°。我国主要水蚀区的研究表明, 土壤侵蚀量与坡度成正比关系^[41], 而25°以上陡坡地更是水土流失的主要发源地, 属于生态恢复的关键区域^[42]。这在一定程度上验证了本文使用模型的正确性和准确性。研究同时识别了各权衡强度等级坡耕地的优化潜力, 其中权衡强度过高或过低的坡耕地均不适宜直接进行生态恢复, 而应当以耕作、工程等措施进行调控。耕作措施包括等高耕作、沟垄种植、覆盖耕作、少耕免耕等, 其中横坡垄作投入少且实施便捷, 保土、保水、保肥作用明显^[43—44]。工程措施则有水平梯田、坡式梯田、隔坡梯田等, 其中水平梯田的土壤侵蚀控制能力最佳^[45]。本研究中两项生态系统服务的权衡曲线的特征表明当追求某一项服务提升超过一定的阈值时, 需要花费的成本激增, 这与其他生态系统服务多目标优化的研究相似^{[15, 46]}。说明在生态系统服务阈值范围内, 以低优化成本换取高收益是切实可行的。

坡耕地的转出容易威胁区域粮食安全, 平衡粮食生产和土壤保持服务是坡耕地利用优化面临的主要问题^[45], 且进行两项生态系统服务的权衡评估可以将结果直观展示在平面坐标系中, 与整合更多生态系统服务的高维空间图相比, 二维空间可视化效果更好。因此本文没有进行更多项生态系统服务的优化, 仅聚焦于研究区典型的粮食生产和土壤保持服务, 未来的研究中应整合更多项生态系统服务进行综合性的优化评估。

4.2 结论

本文评估并可视化了2000—2020年施甸县粮食生产与土壤保持服务的权衡强度, 然后以坡耕地转换为林地作为决策变量, 通过多目标线性规划方法, 分析了粮食生产与土壤保持两项生态系统服务在效率前沿曲线上的权衡关系及相应的坡耕地转出面积, 并且与权衡强度图层叠加分析, 识别了生态恢复高优先区与两项生态系统服务权衡强度在空间上的差异性。主要结论如下：(1)2000—2020年施甸县粮食生产与土壤保持服务在空间上以中、高权衡强度为主导且权衡强度整体在不断上升；(2)仅通过生态恢复手段提升坡耕地土壤保持服务具有阈值效应, 当土壤保持服务提升超过理论最大值的50%后会显著增加优化成本；(3)各权衡强度等级的坡耕地优化潜力为：中权衡强度>低权衡强度>高权衡强度, 因此中权衡强度坡耕地适合直接进行生态恢复, 而低、高权衡强度坡耕地则适合保护性耕作措施进行调控。