清水河流域生态系统服务评估及预测

裴宏伟，张红娟，张良，李雅丽，张高伟，刘孟竹*

1.河北建筑工程学院市政与环境工程系

2.河北省水质工程与水资源综合利用重点实验室

3.河北省张家口水文勘测研究中心

生态系统是人类赖以生存的环境与基础，在生态过程中产生的物质与服务有着无法替代的价值[1]。生态系统服务是人类从生态系统中获得的直接或间接的所有惠益，包括有形的物质供给和无形的服务提供，主要分为供给服务、调节服务、文化服务和支持服务[2]。随着人类活动对生态系统的干扰日益加强，加上不稳定的气候变化，全球超过2/3的生态系统服务出现下降，生态系统退化的程度已经超过了预期水平[3-4]。世界范围内，生态系统服务已经成为地理学、生态学、测绘学、社会科学等相关学科的研究前沿与热点。目前，生态系统服务评估已经上升到国家尺度[5-7]，生态系统服务价值的研究也已在全球尺度上开展[8-9]。生态系统服务与价值的评估对自然资源的合理配置与利用，生态补偿、生态恢复和生态红线政策制定等具有重大意义[10-11]。

由人类活动引起的土地利用变化是影响生态系统服务与价值的主要原因，情景分析法能够准确量化土地利用变化对生态系统服务的具体影响。例如，Bai等[12]基于情景比较法区分了土地利用对生态系统服务的相对贡献；王培俊等[13]通过土地预测模型预测了未来的生态系统服务价值，为未来的生态资产评估提供了可行的方案。然而，基于预测模型或假设情景建立的情景分析法往往难以量化某一种或某几种特定的土地利用变化对生态系统服务的影响。这是因为预测模型模拟出的土地情景产生了全局的变化，无法准确解释最终的生态系统服务变化主要受哪种土地利用转变的影响。为了解决这一问题，部分研究通过设计未来的土地情景来判断某一种人为的土地管理方式是否能够达到预期的生态系统服务供给水平[14-15]，研究的方法是将对应某一土地利用类型的栅格改变为另一种特定的土地利用栅格，从而实现准确地量化某一种土地转变对生态系统服务的影响，但这种设计情景法缺少了对未来的预测。综上，通过模型预测未来的情景，并以此对未来情景进行土地栅格的特定转变的设计（预测+设计），能够达到准确地模拟未来不同发展方向变化的土地情景，进而实现对未来不同情景下生态系统服务和价值的评估。

清水河流域地处河北省西北部，与张家口市崇礼区行政边界较为一致。该流域为半湿润向半干旱过渡的生态脆弱区，也是土地利用受人类活动影响较大的典型区域。为揭示清水河流域未来的生态系统服务和价值的变化，本研究以该流域2000—2018年的产水服务、固碳服务、土壤保持服务以及生态系统服务价值为研究内容，采用预测+设计的情景分析法，在预测未来土地格局自然发展的基础上加入了人为管理的土地利用设计方案，以预测未来自然、人为影响下的生态系统服务与价值变化，促使实现生态系统服务与价值的最大化供给，以期为生态系统服务和价值的预测提供评估框架，同时也为生态可持续发展、土地规划决策提供科学依据。

1 研究区与研究方法

1.1 研究区概况

清水河流域（114°47′E～115°31′E，40°49′N～41°16′N）位于海河上游，处于内蒙古高原和华北平原过渡地区。流域内400 mm等降水量线穿过，多年平均降水量为442 mm，降水主要集中在6—9月，多年平均气温为8.6 ℃，属于大陆性季风气候。流域总面积约 2 200 km2，高程为 780～2 165 m（图 1）。2018年，清水河流域总人口约13万人，人口密度约50人/km2，农、林、牧产值分别占经济总产值的1/2、1/5、1/4，其中仅蔬菜的经济产值就达到了2/5。高强度的农业活动和城镇化快速扩张给当地的水资源、生态环境带来了较大影响。此外，流域气候较为恶劣，暴雨、冰雹灾害时有发生，水土流失较为严重、生态环境极为脆弱[16]。

图 1 清水河流域高程Fig.1 Elevation of Qingshui River Basin

1.2 数据来源与处理

土地利用数据来源于遥感数据解译，在前期研究[16]中已经完成，分辨率为30 m，选取了2000年、2009年、2018年3个年份数据，根据清水河流域实际情况分为耕地、林地、草地、水域、建设用地5个一级土地利用类型。DEM数据（30 m）来源于地理空间数据云（http://www.gscloud.cn）。气象数据来源于中国地面气候资料日值数据集（V3.0）（http://data.cma.cn），包括清水河流域及周边地区4个气象站2000—2018年降水量、相对湿度、温度等数据。由Penman-Monteith公式计算得到每个站点作物参考蒸散量，最后由克里金插值法得到30 m的蒸散量栅格数据（与降水量数据分辨率一致）。年降水量栅格数据（1 km）来源于国家地球系统科学数据中心（http://www.geodata.cn/），相比于插值法得到的降水量数据，栅格数据具有更好的空间分布，空间分辨率经过重采样达到了30 m。流域及子流域矢量图由ArcMap软件中水文分析工具提取得到。粮食产量数据来源于《张家口经济年鉴》（2000—2019年），粮食价格来源于《中国农产品价格调查年鉴2019》，均下载于中国知网（https://www.cnki.net/）。

1.3 研究方法

1.3.1 技术路线

元胞自助机（cellular automata，CA）模型是在空间上相互作用、时间上又具有因果关系的一种网格动力学模型，具有处理复杂空间系统的能力。Markov模型基于Markov链，因其较好的稳定性和无后效性，能够预测土地利用变化中各时刻的变动过程。CA-Markov模型综合了Markov模型的时序预测和CA模型的空间分布模拟，已经被广泛应用于土地利用数量变化、空间格局变化的模拟中[16-17]。不同于CLUE-S模型以及FLUS模型，CA-Markov模型可以将上一阶段的土地利用转移矩阵作为Markov链，进而根据以往时期的土地流转趋势作为预测未来土地格局的方向，因此该模拟方法符合本研究中的自然发展情景，即土地利用变化向着原本的趋势发展。

分析2000—2018年清水河流域的土地利用时空变化，利用CA-Markov模型预测未来自然发展情形下的土地利用情景；在此基础上，通过ArcMap软件中重分类工具设计未来的退耕还林、还草发展情景，得到不同时空尺度下的土地利用情景；根据InVEST（integrated valuation of ecosystem services and trade-offs）模型、当量因子法和当量系数表来模拟不同时空尺度下的产水量、土壤保持量、碳储量以及生态系统服务价值，最后实现清水河流域近20年来的生态系统服务及价值的评估。研究技术路线如图2所示。

图 2 研究技术路线Fig.2 Research technology roadmap

1.3.2 产水服务

使用InVEST模型中的产水量模块计算产水量，以表示产水服务，计算公式如下：

式中：Yij为栅格i中土地利用类型j的年产水量，mm/a；AETij为栅格i中土地利用类型j的年实际蒸散量，mm/a；Pi为栅格i的年降水量，mm/a；PETi为栅格i的年潜在蒸散量，mm/a；ET0,i为栅格i的年参考作物蒸散量，mm/a；Kc，j为土地利用类型j的蒸散系数，无量纲；w为气候-土壤的非物理参数；Z为季节性常数；max(layer_depthi)为栅格i的最大土壤深度，mm；root_depthi为栅格i的最大根系限制深度，mm；AWCi为栅格i的土壤有效含水量，mm；PAWCi为栅格i的植物可利用含水量，其值处于[0，1]。

1.3.3 固碳服务

InVEST碳储量模型中需要的输入数据为土地利用数据以及基本的四大碳库（地上碳库、地下碳库、土壤碳库和死亡有机质碳库）数据，输出结果为总碳储量以及不同土地利用类型的碳储量。本研究只考虑基本的四大碳库，碳密度数据参考同地区研究结果（表1）[18]，用总碳储量表示固碳服务，计算公式如下：

表 1 清水河流域不同土地利用类型碳密度参数Table 1 Carbon density parameters of land use types in Qingshui River Basin t/hm2

式中：C为总碳储量，t；Ca为地上碳密度，t/hm2；Cb为地下碳密度，t/hm2；Cc为土壤碳密度，t/hm2；Cd为死亡有机质碳密度，t/hm2；S为各土地利用类型面积，hm2。

1.3.4 土壤保持服务

InVEST泥沙运移模型在土壤流失流失方程（USLE）的理论基础上进行了改进，在沉积物持留方面，该模型考虑了植物对土壤侵蚀的减缓以及对上坡沉积物具有拦截作用这一重要水文过程，使得模拟结果更科学合理。本研究采用改进的USLE方程，具体计算方法如下：

式中：RKLSi为栅格i基于地貌和气候计算的潜在土壤侵蚀量，t；Ri、Ki、LSi、Ci、Zi为单元栅格i的降水侵蚀力因子、土壤可蚀性因子、坡长坡度因子、植被覆盖与管理因子、土壤保持措施因子；USLEi为栅格i的土壤实际侵蚀量，t；SEDRi为栅格i的土壤持留量，t；SEi为栅格i的泥沙持留率；SEDRETi为栅格i的土壤保持量，t；Ci对应于栅格i中耕地、林地、草地、水域、建设用地的取值分别为0.05、0.03、0.04、0、0， 类似地，Zi对应栅格i中的耕地、林地、草地、水域、建设用地分别取值1、1、1、0、0。其他因子根据模型推荐参数输入。

1.3.5 生态系统服务价值计算

根据谢高地等[19-20]发展的当量因子法来计算生态系统服务价值，计算公式如下：

式中：Vj为土地利用类型j的生态系统服务价值，元；Cjn为土地利用类型j对应的二级生态系统服务价值n（食物生产价值、气候调节价值等），元；Sj为土地利用类型j的总面积，hm2；ESV为研究区总的生态系统服务价值，元。

以单位面积农田生态系统粮食生产的净利润当作1个标准当量因子的生态系统服务价值量，标准当量因子的值参考以往研究定义为实际粮食产值的1/7[21]，当量系数表参照文献[19]。根据《全国农产品成本收益资料汇编》，2018年全国3种稻谷粮食收购价格平均值为2.51元/kg，清水河流域2000—2018年粮食产量均值为19 314 t/a，粮食种植面积均值为8 374 hm2/a，其中2001年、2002年粮食种植面积数据缺失，采用2000年、2003年的数据利用线性插值法得到。参考谢高地等[19]发展的植被净初级生产力（net primary productivity，NPP）修正方法，崇礼区2018年平均NPP为422.38 g/(m2·a)，全国水平为338.16 g/(m2·a)，修正系数为1.25。最终计算得到标准当量因子为1 033.82元/hm2，进而得到2000—2018年清水河流域不同土地利用类型单位面积生态系统服务价值（表2）。

1.3.6 未来情景构建

结合CA-Markov模型以及ArcMap软件中的重分类工具来实现土地情景预测及设计。按照清水河流域自然发展情景，通过CA-Markov模型模拟出2027年流域的土地利用分布（参考前期研究[16]），在ArcMap软件中对该栅格实施不同规则的栅格属性重分类，得到不同力度生态修复措施下的土地利用情景。重分类规则参考清水河流域当地的退耕还林、还草政策而制定，不同情景设计具体见表3。

2 结果和分析

2.1 清水河流域2000—2018年生态系统服务时空变化

2000—2018年清水河流域碳储量、土壤保持量、产水量及其空间变化分别如表4和图3所示。碳储量由2000年的2.23×107t减至2009年的2.21×107t，2018年进一步减至2.20×107t。整个时期清水河流域单位面积碳储量减少约1.10 t/hm2（表4），固碳服务呈下降趋势。研究区土壤保持总量在2000—2009 年由 175.92×106t增至 176.03×106t，2018 年增至176.04×106t，单位面积土壤保持量在近20年增加约0.54 t/hm2，土壤保持服务呈轻微增强趋势。研究区产水量由2000年的123.77×106m3增至2009年的 127.18×106m3，2018 年又增至 130.78×106m3，单位面积产水量增加了3.19 mm，产水服务相较于固碳和土壤保持服务明显增强。

表 2 清水河流域不同土地利用类型单位面积的生态系统服务价值Table 2 Ecosystem service values per unit area in different land use types in Qingshui River Basin 元/hm2

表 3 清水河流域不同力度生态修复措施下的土地利用情景Table 3 Land use scenarios under different ecological restoration measures in Qingshui River Basin

表 4 清水河流域单位面积生态系统服务变化Table 4 Changes of ecosystem services per unit in Qingshui River Basin

近20年来，清水河流域碳储量变化空间差异较为显著，固碳功能较强的区域分布在流域西北部和南部等林地覆盖度较高的区域，而碳储量减少的区域主要分布在城镇扩张以及草地转变为耕地的区域（图3）。土壤保持量空间格局与高程、坡度以及土地利用分布格局相关，在海拔低、坡度平缓的耕地区域，土壤保持量较低，在海拔高、坡度陡的林草地区域土壤保持量较高。

产水量空间分布在全流域尺度上与降水量分布有关，在区域尺度上则与土地利用类型分布相关。研究区东北部多年平均降水量较高，因此该区域产水量高于其他地区。从不同土地利用类型的产水量而言，建设用地产水量较高，通常建设用地下垫面缺少土壤无法存续水，导致降水直接形成了径流产生更多的水量，产水量显著增加的区域也主要分布于建设用地扩张区域；林地区域的产水量整体偏低，这是由于林地区域的年降水量远远小于潜在蒸散量（潜在蒸散量约是降水量的2倍以上），林地的高蒸散发增加了更多的水分消耗。尽管林地的增加减少了生态系统产水服务，但林地的降水截留以及固土、持水能力也能够调节径流、涵养水源，提供了潜在的额外生态系统服务与价值。

2.2 不同情景下生态系统服务预测

通过土地利用预测及设计的方式得到清水河流域未来不同发展情景下的生态系统服务及其变化（图4）。以2018年碳储量为基准，清水河流域在自然发展情景下2027年单位面积碳储量减少了0.56 t/hm2，退耕还林情景下单位面积碳储量增长了0.59 t/hm2，而加大还林力度后则增长了4.32 t/hm2，在退耕还草和深度退耕还草情景下分别减少了2.63、3.70 t/hm2。在退耕还草后，碳储量有所降低，这是由于在该区域碳库中草地的碳密度小于农田的碳密度，因此在耕地转变为草地后总碳储量减少。

图 3 清水河流域2000—2018年生态系统服务变化Fig.3 Changes of ecosystem services in Qingshui River Basin from 2000 to 2018

图 4 不同预测情景下2027年相较于2018年生态系统服务变化Fig.4 Changes of ecosystem services in 2027 compared with 2018 under different forecast scenarios

在自然发展情景下，土壤保持量呈减少趋势，2027年单位面积土壤保持量相较2018年减少了0.08 t/hm2，这主要是由于自然情景下植被覆被土地面积减少，耕地和建设用地面积增多，导致土地的固土能力减弱。情景分析结果显示，退耕还林措施能够逆转未来土壤保持服务减弱的趋势，将坡度大于25°的耕地改造成林地后，2027年清水河流域单位面积土壤保持量相较于2018年增加了0.19 t/hm2，而将坡度大于15°的耕地改造成林地后，单位面积土壤保持量增加了0.71 t/hm2；退耕还草情景下单位面积土壤保持量增加了0.10 t/hm2，深度退耕还草情景下单位面积土壤保持量增长了0.37 t/hm2。

在气候不变的情况下，2018—2027年清水河流域植被面积将会减少，建设用地和耕地面积增加；相应地，流域的实际蒸散发量将减小，径流、土壤中水以及下渗水等地表水将增多，最终表现流域生态系统的产水服务增强。在退耕还林情景下，2027年产水量相较于还林前减少了0.67 mm，而深度退耕还林情景下产水量相较还林前减少了2.84 mm；在轻度、重度退耕还草情景下产水量相较于还草前分别减小了0.41、1.77 mm。总体而言，力度较小的还草措施会增大流域生态系统产水服务，而力度较大的还草措施则会减少产水服务。

2.3 不同情景下生态系统服务价值变化

2000—2018年以及2027年不同情景下清水河流域生态系统服务价值变化如图5所示。2000年、2009年和2018年清水河流域生态系统服务价值分别为24.97亿、25.38亿、24.42亿元，单位面积生态系统服务价值分别为11 360、11 550、11 110元/hm2。2018年，研究区生态系统服务价值分别是该区当年农林渔牧业产值（10.59亿元）和GDP（31.37亿元）的2.31、0.78倍。流域生态系统服务价值在 2000—2009年增长了0.41亿元（+1.66%），在2009—2018年减少了0.96亿元（-3.77%）。经分析，2000—2009年研究区生态系统服务价值增加是由于水域、林地等具有高生态价值的土地利用类型扩张，2009—2018年减小则是因为植被、水域等土地面积减少以及城镇化扩张。总体来讲，2000—2018年清水河流域生态系统服务价值呈轻微减小趋势，减小幅度约为0.55亿元，主要是由于城镇化的扩张导致原本具有多种生态功能和高生态价值的土地萎缩进而减小了生态系统服务价值。

图 5 2000—2018年以及不同预测情景下单位面积生态系统服务价值的变化Fig.5 Changes of ecosystem service values per unit area in 2000-2018 and under different forecast scenarios

与基础情景(2018年)相比较，在自然发展情景下2027年清水河流域单位面积生态系统服务价值降低了590元/hm2，退耕还林情景下单位面积生态系统服务价值降低了150元/hm2；深度退耕还林情景下单位面积生态系统服务价值出现显著增高趋势，相较于2018年水平增长了1 850元/hm2；退耕还草、深度退耕还草情景下，流域单位面积生态系统服务价值相较于2018年分别减少约30、110元/hm2。

3 结论与建议

（1）2000年、2009年和2018年，清水河流域单位面积碳储量平均水平为100.68 t/hm2，流域生态系统固碳服务近20年呈下降趋势，单位面积碳储量下降了1.10 t/hm2；单位面积土壤保持量平均水平为800.62 t/hm2，土壤保持服务近20年呈上升趋势，单位面积土壤保持量增加了0.54 t/hm2；3个年份的产水量平均值为57.884 mm，流域生态系统产水服务增强，近20年产水量增加了3.19 mm。

（2）2000年、2009年、2018年清水河流域生态系统服务价值分别为24.97亿、25.38亿、24.42亿元，呈先增加后减少的趋势。2018年，流域生态系统服务价值是同时期GDP的0.78倍，这体现了清水河流域生态系统服务价值的高度稀缺性。

（3）在大于15°的陡坡区域实施退耕还林将增加175.98 km2的林地，同时也会显著提升清水河流域固碳、土壤保持服务和生态系统服务价值，但扩张的林地同时会产生更多蒸散量从而削弱生态系统的产水服务。

（4）针对清水河流域近20年来生态系统服务评估结果，提出如下建议：研究区整体坡度较高，应加大15°以上陡坡的植树造林力度，以防止水土流失，增加生态系统固碳功能；增加对草地修复的人力物力投资，严禁过度草地开垦、放牧等；进一步推进生态退耕、生态补偿，减少耗水蔬菜种植规模，以保障流域生态可持续发展。