生态系统服务是指人类从生态系统中所获得的收益, 包括供给、调节、支持和文化等服务。由于服务种类多样、空间分布不均衡以及人类使用的选择性, 生态系统服务之间表现为此消彼长的权衡和相互增益的协同关系^[1]。生态系统管理需要兼顾多种服务, 增强协同效应并减弱权衡效应, 最大化综合效益而非只追求某一种服务的效益^[1—4], 因此开展权衡与协同作用机制、类型及其特征的分析研究很有必要^[5]。近年来, 已有大量研究关注区域尺度的生态系统服务权衡与协同关系, 研究方法主要包括统计学方法、空间分析、情景模拟与权衡模型等^[3], 然而目前关于生态系统服务相互作用机理研究及尺度效应尚不明确。

东北地区黑土肥沃, 粮食产量和商品粮分别占全国的1/4和1/3, 是我国重要的商品粮基地^[6—8], 每年能满足超过2.16亿居民的粮食需求^[9—10], 是国家粮食安全的“稳定器”和“压舱石”^[11], 此外东北是我国湿地资源的主要分布区之一, 被列入国际湿地公约保护区名录的重要湿地数量占全国的1/3, 同时东北还拥有我国面积最大的天然林(面积占全国的37%), 不仅为居民提供农业生产等供给服务, 是我国的粮仓, 还提供了水源涵养、土壤保持和防风固沙等调节服务, 是我国重要的生态屏障。东北生态系统服务是近年来研究的热点, 其时空演变及权衡与协同关系是气候变化和人类活动综合作用的结果。干旱、热浪、强降水等极端气候事件的持续增加, 导致东北粮食生产潜力波动^[12], 粮食生产风险区面积增大, 部分地区粮食减产量可达10%^[13]。此外, 随着人口快速增长和经济社会快速发展, 东北城市化进程加快, 城乡建设用地不断扩展, 极大地改变了陆地生态系统的格局、过程、功能, 导致生态系统服务及其相互关系发生极大变化。气候变化和人类活动导致了土壤侵蚀^[14—21]、土壤有机质下降^{[8, 22—26]}、土壤酸化^[7]等土地退化问题^{[6, 27—28]}, 造成了耕地生产潜力下降^[29]、增产质量不高、种植结构单一化和粗放管理、农业生产不确定性增加等供给和产能方面的问题^[30], 同时影响了森林、草地和湿地的水源涵养、土壤保持和防风固沙等生态功能, 生态系统服务之间矛盾凸显, 生产和生态功能协同发展受到限制。当各种服务之间存在冲突时, 应进行权衡分析以减少冲突而增加协同^[31], 否则可能会导致某些服务功能供给能力下降, 降低生态系统的稳定性^[32]。

东北生态系统服务的权衡与协同研究涉及粮食供给、水源涵养、土壤保持、防风固沙、固碳释氧、生物多样性等服务。已有研究主要聚焦于局部地区, 分析的基本单元多为县域, 且多为起始独立年份的简单变化分析, 然而东北黑土区大尺度、长时间序列的研究相对较少, 在宏观层面上对东北黑土区生态系统服务的权衡与协同关系还缺乏统一的认识。此外, 研究方法主要是统计学方法, 例如相关分析^[33—36]、空间自相关分析^[37—38]、聚类分析^[39—40]、均方根误差^[41]、权衡协同度^[42—43]等。这些方法虽然能有效揭示服务间权衡与协同关系的科学规律, 但对其空间分异特征的表达有所欠缺。空间制图方法可以弥补这一不足, 但通常只考虑前后两个年份而非连续的长时间序列的生态系统服务, 容易产生较大误差。逐像元相关分析方法是评价区域生态系统服务权衡与协同关系的可行方法, 该方法基于长时间序列的生态系统服务的变化趋势, 能够准确揭示服务间权衡与协同关系的空间分异特征, 又能够清晰地反映服务功能连续的时间变化规律^[44]。综上, 本文基于CASA(Carnegie-Ames-Stanford Approach)模型、降水贮存量法、修正的通用土壤流失方程(Revised Universal Soil Loss Equation, RUSLE)和修正的风蚀方程(Revised Wind Erosion Equation, RWEQ)量化了东北黑土区供给、水源涵养、土壤保持和防风固沙4种关键生态系统服务并分析其时空演化特征, 采用多元回归残差分析法厘定气候变化与人类活动因素对供给服务的贡献率, 采用逐像元相关分析法探讨了供给服务与调节服务之间的权衡与协同关系, 明确了气候变化与人类活动因素对权衡与协同的影响。研究创新之处在于建立了一套东北黑土区长时间序列的1 km分辨率的栅格数据集, 在宏观层面上揭示了不同尺度下东北黑土区生态系统服务的时空分布、驱动因素以及权衡与协同关系, 可为东北黑土区生态保护修复相关政策的制定提供重要的科学依据, 对保障粮食和生态安全、实现黑土地生态高效可持续利用具有重要意义。

1 数据与方法 1.1 研究区概况

研究区的范围包括黑龙江、吉林、辽宁3省以及内蒙古自治区东部, 面积约为125万km², 占全国总面积的13%(图 1)。地形以平原、丘陵和山地为主, 东、西、北部三面环山, 山地之间分布着东北平原。气候以温带大陆性气候和温带季风气候为主, 不同地区气温和降水差异明显。东北黑土地土壤性状好、肥力高, 是适宜农耕的优质土地, 但由于长期超负荷利用, 耕地质量逐年下降, 黑土地“变薄、变瘦、变硬”的问题日趋严重, 黑土地肥力下降、水土流失、土壤酸化等土地退化问题未得到充分的重视。此外由于林地开垦等人类活动的影响, 切沟侵蚀呈加剧趋势, 且在未来气候变化情景下, 东北黑土区土壤侵蚀具有增大风险^[20]。为了缓解东北生态系统服务间的矛盾与冲突, 促进生产与生态功能的协同发展, 有必要深入开展生态系统服务权衡与协同研究。

1.2 数据来源与预处理

本文所用到的数据包括土地利用(LUCC)数据、气象数据、数字高程模型(DEM)数据, 归一化植被指数(NDVI)数据。

土地利用数据：来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn)土地利用遥感监测数据集^[45], 该数据集是以美国陆地卫星Landsat遥感图像为主要信息源, 通过人工目视解译获得, 数据空间分辨率为30 m, 时间跨度为2000—2020每5年一期。提取研究区范围内的耕地、林地、草地以及水体和湿地, 进行生态系统服务的量化、驱动因素分析以及权衡与协同分析。

气象数据：来源于国家气象科学数据共享平台(http://data.cma.cn)2000—2020年的日值观测数据, 包括日均风速、风向、降水、温度、日照时数等。采用ANUSPLIN进行插值得到空间分辨率为1 km的栅格数据。

DEM数据：来源于CIAT-CSI SRTM网站(https://srtm.csi.cgiar.org)SRTM数据, 空间分辨率为90 m, 用于计算坡度和坡长。

NDVI数据：获取MODIS 2000—2020年的NDVI时间序列数据产品(MOD13Q1), 空间和时间分辨率分别为250 m和16 d。经过格式转换、重投影、拼接、重采样和S-G滤波处理, 采用最大合成法得到月度NDVI数据, 然后计算植被覆盖度：

式中, f为植被覆盖度；NDVI_veg和NDVI_soil分别是纯植被和裸土像元的NDVI值, 依据土地利用类型确定纯植被和裸土的NDVI。

数据预处理：在研究区范围建立1 km×1 km网格, 基于该网格对已有栅格数据分辨率和坐标系进行统一, 以便后续的数据处理和分析。

1.3 研究方法 1.3.1 生态系统服务量化

(1) 供给服务

生态系统净初级生产力(Net Primary Productivity, NPP)是绿色植被在单位面积、单位时间内累积的有机物数量, 由光合作用获得的有机质总量扣除植被自身生长和代谢消耗后的剩余部分, 它直接反映了自然环境条件下生态系统的供给能力, 可作为供给服务量化的指标^[46]。以NPP为指标来表征东北黑土区供给服务, 使用基于光能利用率原理的CASA模型对NPP进行估算^[47], 具体模型如下：

式中, x为像元；t为时间；APAR为植被所吸收的光合有效辐射(MJ/m²)；ε为光能转化率(gC/MJ), NPP(x, t)为像元x在t时间内的净初级生产力(gC m^-2 a^-1), APAR(x, t)为空间位置x在时间t内吸收的光合有效辐射(MJ/m²)；ε(x, t)表示像元x在时间t内的实际光能利用率(gC/MJ)^[48]。根据上述公式最终得到2000—2020年逐年1 km的NPP栅格数据。

(2) 水源涵养服务

采用降水贮存量法^[49]估算东北森林、湿地、草地生态系统通过调节水量和净化水质所产生的水源涵养服务量, 公式表示为：

式中, W_r为与裸地相比森林、湿地和草地生态系统涵养水分的增加量(m³), A为生态系统面积(km²), J为产流降水量(mm), J₀为年均降水量(mm), K₀为产流降水量占降水总量的比例；R₀为与裸地相比生态系统减少径流的效益系数。通过已发表文献搜集获得实测降雨产流临界值, 并以临近气象台站实测的日降水量修正同时期得到的TRMM逐日3 h降水量, 累计单次降雨量大于降雨产流临界值的数值, 得到单点产流降雨量占降雨总量的比例, 进而与多年平均河川径流系数建立线性关系, 得到K₀值的空间分布。森林R₀值主要通过搜集已有文献资料, 草地R₀值根据植被覆盖度^[50—51]计算。根据上述公式最终得到2000—2020年逐年的1 km水源涵养服务量栅格数据。

(3) 土壤保持服务

利用生态系统土壤保持量(SC)衡量生态系统保育土壤的能力, 即生态系统在极度退化状况下或没有土地覆盖管理措施情景下土壤侵蚀量(SE_d)与现实状况下土壤侵蚀量(SE_a)的差值。基于RUSLE模型计算土壤水蚀模数：

式中, SE为水蚀模数(t/hm²)；R为降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1), 利用基于日降雨量的半月降雨侵蚀力模型计算得出^[52]；K为土壤可蚀性因子(t h MJ^-1 mm^-1), 由Nomo模型^[53]根据1 ∶ 100万土壤数据估计得出；L(m)为坡长因子, S(%)为坡度因子, 由McCool等^[54]和Liu等^[55]的核心算法计算；C为覆盖和管理因子, 利用植被覆盖度f计算^[56]；P是侵蚀控制实践因子, 数值范围在0—1之间, 主要根据相关文献和专家咨询确定, 草地、森林和裸地的P值定为1, 水田为0.15, 旱地为0.352, 定居点为0.01, 湿地/水体为0^[57—58]。根据上述公式最终得到2000—2020年逐年的1 km土壤保持服务量栅格数据。

(4) 防风固沙服务

防风固沙量(SF)为裸土条件和地表覆盖植被条件下的土壤风蚀量的差值：

式中, SF表示防风固沙服务量, SL_d表示裸土条件下的潜在土壤风蚀量, SL_a表示植被覆盖条件下的现实土壤风蚀量。

充分考虑气候条件、植被状况、地表土壤的粗糙度、土壤可蚀性和土壤结皮的情况下, 利用RWEQ模型量化土壤风蚀：

式中, SL为风蚀模数(t/hm²)；x为地块长度(m)；Q_x为x处的沙通量(kg/m)；Q_max为风力的最大输沙能力(kg/m)；s为关键地块长度(m)；WF为气象因子, 反映了风因子、空气密度、重力加速度、土壤湿润度因子和雪盖因子对风蚀的综合影响；EF是土壤可蚀性因子；SCF是土壤结皮因子；K′是土壤糙度因子；COG是综合植被因子^[57—60]。WS₂为2 m高度风速(m/s), WS_t为2 m高度临界风速, N为风速的观测次数, N_d为试验的天数, ρ为空气密度(kg/m³), g为重力加速度(m/s²), SW为土壤湿度因子, SD为雪覆盖因子；Sa为土壤砂粒含量, Si为土壤粉砂含量, Cl为土壤粘土含量, OM为有机质含量, CaCO₃为碳酸钙含量。根据上述公式最终得到2000—2020年逐年的1 km防风固沙服务量栅格数据。

1.3.2 驱动因素分析

本文利用多元回归残差分析法^[61—63]厘定气候变化与人类活动对供给服务的贡献率, 识别各地区供给服务的主导驱动因素。以气温(T)和降水量(P)为自变量, 供给服务(NPP_obs)为因变量建立二元线性回归方程；假设供给服务仅受气候变化影响, 根据回归方程以及气温和降水量数据, 计算得到气候变化驱动下供给服务(NPP_cc)；通过NPP_obs与NPP_cc的差值NPP_ha表示人类活动对生态系统服务的影响。

式中, T和P分别为气温与降水量；a、b和c为回归模型参数。通过一元线性回归法计算NPP_cc与NPP_ha的趋势Slope_nppcc和Slope_nppha, 趋势为正表示驱动因素对供给服务起到促进作用, 反之则表示抑制作用；两者与供给服务实际变化趋势(Slope_nppobs)的比值, 即气候变化或人类活动对供给服务变化的相对贡献率(表 1), 然后通过比较两者的大小来确定主导驱动因素。

1.3.3 生态系统服务权衡与协同分析

以东北黑土区供给服务量与水源涵养量、土壤保持量和防风固沙量的多年变化趋势为基础, 基于逐像元空间相关分析法, 通过判断两两服务间相关系数的绝对值大小及正负方向, 来确定栅格和县域尺度的权衡与协同关系, 并进行空间明晰的定量制图^[44]。相关系数r计算公式如下：

式中, x_i和y_i分别是某像元第i年的生态系统服务量, n为总年数, 共21年。相关系数为正, 则两种服务之间表现为协同关系；反之则为权衡关系。根据相关系数的零假设T检验方法判断显著性：

式中, r是相关系数, n是样本观测数, k是可控制变量的数目, n-k-2是自由度, 当|T|>T_0.05n-k-2时即P＜0.05, 拒绝原假设, 相关性结果显著, 当| T |>T_0.01n-k-2时即P＜0.01, 拒绝原假设, 相关性结果极显著。T_0.05n-k-2 和T_0.01n-k-2 查阅T检验表确定。查表时, 根据经验两种服务间可能存在正相关和负相关, 所以用双侧界线。最终将权衡与协同关系按照不同程度分为强协同(r>0, P≤0.05)、中协同(r>0, 0.05＜P≤0.1)、弱协同(r>0, P>0.1)、强权衡(r＜0, P≤0.05)、中权衡(r＜0, 0.05＜P≤0.1)和弱权衡(r＜0, P>0.1)。

2 结果与分析 2.1 生态系统服务时空演变特征 2.1.1 供给服务

东北黑土区NPP多年平均值约为406.79 gC m^-2, 多年平均总量为473.98 Tg, 其中吉林多年平均值最高(454.58 gC m^-2), 内蒙古东部最低(347.10 gC m^-2)。2000—2020年, 东北黑土区NPP总体上呈现小幅提升趋势, 平均每年提升4.24 gC m^-2, 其中内蒙古东部和黑龙江增势最明显, 年际变化量达到5.41 gC m^-2 a^-1和4.11 gC m^-2 a^-1, 吉林和辽宁增势较小, 年际变化量仅为2.92 gC m^-2 a^-1和2.46 gC m^-2 a^-1。从县域来看, 96.6%的县域NPP表现为上升趋势, 仅有3.4%的县域为下降趋势(图 2, 表 2)。

2.1.2 水源涵养服务

东北黑土区森林、草地和湿地生态系统多年平均水源涵养总量约为11.53×10¹⁰ m³, 单位面积水源涵养量为11.67×10⁴ m³/km², 其中黑龙江和吉林较高, 分别约为15.25×10⁴和14.40×10⁴ m³/km², 而辽宁和内蒙古东部较低。2000—2020年, 东北黑土区水源涵养量总体呈现小幅提升趋势, 平均每年上升0.08×10⁴ m³/km², 其中黑龙江和吉林增势较明显, 年际变化量分别为0.11×10⁴ m³/km²和0.10×10⁴ m³/km², 而内蒙古东部和辽宁上升趋势较小。县域尺度来看, 76.9%的县域水源涵养量表现为上升趋势, 而23.1%的县域为下降趋势(图 2, 表 2)。

2.1.3 土壤保持服务

东北黑土区多年平均土壤保持总量约为10.54×10⁸ t, 单位面积土壤保持量为9.08 t/hm²。其中辽宁单位面积土壤保持量最高(28.72 t/hm²), 其次是内蒙古东部(9.70 t/hm²)和吉林(6.37 t/hm²)。2000—2020年, 东北黑土区单位面积土壤保持量平均每年下降0.01 t/hm², 从各行政区来看, 黑龙江、吉林下降, 而辽宁、内蒙古东部上升。从县域来看, 60.7%的县域土壤保持量下降, 而39.3%的县域上升(图 2, 表 2)。

2.1.4 防风固沙服务

东北黑土区多年平均防风固沙总量约为161.56×10⁸ t, 单位面积防风固沙量为139.18 t/hm², 其中吉林和内蒙古东部单位面积防风固沙量较高, 分别为194.56和172.32 t/hm², 而黑龙江和辽宁则相对较低。2000—2020年, 东北黑土区防风固沙量总体上表现为轻微上升趋势, 每年增量为0.02 t/hm², 其中辽宁增势最明显(0.30 t hm^-2 a^-1), 其次是吉林(0.15 t hm^-2 a^-1), 而黑龙江和内蒙古东部则表现为下降趋势。从县域来看, 24.6%的县域防风固沙量表现为下降趋势, 而75.4%的县域为上升趋势(图 2, 表 2)。

2.2 气候变化与人类活动驱动因素分析 2.2.1 驱动类型分析

东北黑土区供给服务变化的驱动类型主要包括6类：CH-D、C-D、H-D、H-I、C-I以及CH-I(图 3)。表现为CH-I的面积最大, 约为65.95 km², 约占总面积的53.13%, 其中吉林和辽宁表现为CH-I的面积占比最大, 分别为71.94%和66.76%。其次是H-D, 面积为23.17 km², 占18.67%, 主要分布于内蒙古东部和黑龙江。表现为CH-D和C-I面积分别占11.07%和13.27%, 表现为C-D和H-I的面积较少, 面积仅占0.50%和3.36%(表 3)。从县域尺度来看, 约有75.1%的县域生态系统服务变化驱动类型表现为CH-I；10.3%的县域表现为H-D, 主要分布于内蒙古呼伦贝尔；而表现为CH-D、H-I和C-I的县域数量占比分别为5.0%、3.1%和6.5%(图 3)。因此, 东北黑土区生态系统服务变化的驱动类型以气候变化与人类活动的共同促进作用为主。

2.2.2 主导驱动因素

从供给服务变化的主导驱动因素及作用率来看(图 4, 表 4), 总体上主导驱动因素为人类活动, 贡献率为62.65%, 气候变化的贡献率为37.35%, 各省人类活动作用程度均超过60%。从面积统计来看, 人类活动主导驱动的面积约为85.81 km², 占总面积的69.6%；气候变化为主导驱动因素的面积约为37.49 km², 占总面积的30.4%。各省供给服务主导驱动因素及作用率差异不大, 人类活动主导面积占比均在70%左右。从县域尺度来看, 约有23.4%的县域主导驱动因素为气候变化, 其中黑龙江鹤岗市、伊春市的部分县域气候变化贡献率超过了75%；76.6%的县域主导因素为人类活动, 其中人类活动贡献率超过75%的县域数量约占38.9%。由此可见, 东北黑土区供给服务变化的主导驱动因素为人类活动。

2.3 生态系统服务权衡协同特征 2.3.1 供给与水源涵养

东北黑土区供给与水源涵养服务总体上以协同关系为主, 面积约占79.17%, 表现为强权衡和强协同的区域面积分别约占总面积的6.32%和47.54%(图 5, 表 5)。从各省表现为强权衡的面积占比来看, 吉林为各省最大(9.16%), 其次是内蒙古东部(7.08%)和黑龙江(6.27%), 辽宁占比最小(2.69%)。从各省强协同面积占比来看, 内蒙古东部为各省最大(69.41%), 吉林最小, 仅为19.42%。

从县域尺度来看, 分别有88.67%和11.33%的县域表现为协同与权衡, 表现为强协同和强权衡的县域数量分别占县域总数的53.72%和1.29%。其中, 强协同面积占比过半的县域数量约占42.72%, 强权衡面积高占比县域主要分布于吉林白山、延边以及黑龙江绥化、大庆等地(图 6), 这些县域多表现为供给服务量升高、水源涵养量降低, 因此需要保护森林、草地、湿地生态系统, 实施退耕还林、封山育林等工程, 增加地表蓄水, 减少径流和蒸发, 提高水源涵养能力以增加协同效应。

从不同驱动因素影响下的权衡与协同特征来看, 驱动类型为CH-D的区域主要表现为权衡, 而其他驱动类型表现为协同。人类活动主导区表现为权衡的面积占22.63%, 表现为强权衡面积占8.18%, 均高于气候变化主导驱动的区域, 因此人类活动因素加剧了供给与水源涵养服务之间的权衡关系(图 7)。

2.3.2 供给与土壤保持

东北黑土区供给与土壤保持服务总体上以权衡关系为主, 面积约占51.89%, 表现为强权衡和强协同的区域面积分别约占总面积的3.25%和4.62%(图 5, 表 5)。从各省表现为强权衡的面积占比来看, 辽宁为各省最大, 达到了7.27%, 其他各省相差不大。而从各省表现为强协同的面积占比来看, 辽宁和内蒙古东部面积占比较大, 分别占该省的8.45%和6.02%, 而黑龙江的强协同面积占比最小(2.65%)。

从县域尺度来看, 分别有31.75%和68.25%的县域表现为协同与权衡, 其中大部分县域表现为弱协同和弱权衡, 表现为强协同和强权衡的县域数量分别仅占县域总数的3.17%和0.01%。强权衡面积高占比县域主要分布于辽宁大连、盘锦, 黑龙江齐齐哈尔以及内蒙古通辽等地；强协同面积高占比县域分布于辽宁朝阳、内蒙古赤峰市等地(图 6), 这些县域多表现为供给服务量升高、土壤保持量降低, 因此可以通过实施水土保持工程, 增加植被覆盖, 减少水土流失。

从不同驱动因素影响下的权衡与协同特征来看, 驱动类型为H-I、C-D和CH-I的区域主要表现为权衡, 而其他驱动类型表现为协同。人类活动主导区和气候变化主导区表现为权衡的面积分别占54.41%和46.41, 表现为强权衡面积分别占3.61%和2.45%, 对比发现人类活动因素加剧了供给与土壤保持服务之间的权衡关系(图 7)。

2.3.3 供给与防风固沙

东北黑土区供给与防风固沙服务总体上以协同关系为主, 面积约占61.77%, 表现为强权衡和强协同的区域面积分别约占总面积的4.25%和25.23%(图 5, 表 5)。从各省表现为强权衡的面积占比来看, 黑龙江为各省占比最大(6.82%), 其次是内蒙古东部(3.92%), 而辽宁最小(0.21%)。从各省表现为强协同的面积占比来看, 辽宁和吉林的强协同面积占比最大(46.90%和34.67%), 而黑龙江的强协同面积占比最小, 仅为14.53%。

从县域尺度来看, 分别有77.14%和22.86%的县域表现为协同与权衡, 表现为强协同和强权衡的县域分别占总数的39.68%和4.13%, 其中表现为强权衡的县域主要分布在黑龙江鸡西、牡丹江以及内蒙古呼伦贝尔等地(图 6), 这些县域表现为权衡的原因主要是防风固沙量的下降, 因此可以通过建设防护林, 利用林木的阻挡和滞留作用减少风速, 降低风沙运动能力, 提高风沙沉降率。此外草原地区可以通过控制放牧强度, 防止草地退化, 提高生态系统防风固沙的能力。

从不同驱动因素影响下的权衡与协同特征来看, 驱动类型为CH-D、C-D和H-D的区域主要表现为权衡, 而其他驱动类型表现为协同。人类活动和气候变化主导驱动的区域均表现为协同关系, 人类活动主导区域表现为协同关系的面积占比略低于气候变化主导区域, 因此人类活动因素减弱了供给与防风固沙服务之间的协同关系(图 7)。

3 讨论与结论 3.1 讨论

本文对东北黑土区生态系统服务的时空演变进行了深入分析。从生态系统服务量化来看, 已有研究结果显示, 黑龙江近20年水源涵养总量提升46.92%^[64], 三江平原和大兴安岭寒温带地区水源涵养量显著提升^[65], 大兴安岭南部、松嫩平原、三江平原土壤保持量下降^[66], 沈阳市各生态系统服务提升^[40], 这些结果都与本文生态系统服务的变化趋势结果一致。此外, 土壤保持和防风固沙量的空间分布结果与Wang等^[67]关于水蚀和风蚀的研究结果是契合的。从各服务间权衡与协同特征结果来看, 已有研究表明, 大兴安岭北部、长白山南部供给与水源涵养服务之间表现为协同关系, 松嫩平原供给与土壤保持服务之间表现为权衡关系^[37], 宾县生态系统服务之间以协同关系为主且空间分异性显著^[36], 内蒙古东部各服务之间以协同关系为主, 而东北平原粮食生产区供给与土壤保持服务之间存在显著的权衡^[39], 哈长城市群北部供给与土壤保持服务之间多为权衡关系^[42], 这些研究都验证了本文结果的准确性。

通过分析东北黑土区生态系统服务的时空演化特征及权衡与协同关系, 可以看出部分地区生态系统服务下降, 服务之间的权衡关系显著, 研究区仍有1/5的县域服务之间主要表现为权衡关系, 这些都将是东北黑土区生产—生态协同发展的负面因素。权衡程度较高的县域, 政府管理部门应在保证黑土地产能的同时, 尽可能多关注该地区的生态功能。不同主导驱动因素的县域需要采取相应的措施, 例如气候变化驱动下的县域应加强农业水利基础设施建设, 减少强降雨、旱涝等极端气候事件的不利影响；人类活动驱动下的县域应加强农田周边森林、草地的保护, 开展植树造林、退化草地修复、土地退化防治工作, 减少土地盐碱化、沙化、土壤侵蚀的不利影响。

近年来, 政府部门高度重视东北黑土地保护和生态修复工作。2021年颁布文件指出“实施国家黑土地保护工程, 推广保护性耕作模式”, 通过实施免耕、少耕及秸秆还田等, 通过提高土壤结构的稳定性, 提升土壤蓄水保水能力, 有效减少水土流失^[68], 从而有利于土壤养分的积累、显著提高土壤肥力和粮食产量。政府部门发布了文件《东北森林带生态保护和修复重大工程建设规划(2021—2035年)》, 提出到2035年, 要全面提升东北森林带生态系统质量和稳定性, 提升优质生态产品供给能力, 建设维护东北地区生态安全、粮食安全、气候安全的重要屏障。此后国家颁布实施了《黑土地保护法》和《关于加强新时代水土保持工作的意见》, 对推进水土流失治理、农田建设、保护性耕种、生态保护等提出了新的要求。

本文的不确定性和不足主要表现在几个方面：(1)本文在东北黑土区生态系统服务的量化工作中, 仅对20年来各服务时空演化规律进行了分析, 并对供给服务受气候变化与人类活动的影响进行了初步探讨, 但对更具体的驱动因子定量分析有待系统和深入的研究。(2)为了进一步提高生态系统服务的估算精度与可信度, 需要进一步结合地面调查对方程中的管理因子进行参数本地化。(3)本文分别单独讨论了各生态系统服务的状况, 但对它们之间的相互影响并未深入研究。(4)本文采用残差分析法对东北黑土区供给服务驱动因素及贡献程度进行了研究, 由于各调节服务驱动因素复杂且服务类型较多, 因此并未考虑全部服务。本文以气温和降水量为自变量与供给服务拟合, 没有考虑高程、坡度和距离等其他因子, 因此存在一定的不确定性。

3.2 结论

本文量化了东北黑土区NPP、水源涵养、土壤保持和防风固沙4种关键生态系统服务并分析了其时空演化特征, 厘定了气候变化与人类活动因素对供给服务的贡献率, 刻画了不同地区、不同主导驱动因素下供给服务与调节服务间的权衡与协同特征, 结果显示：

(1) 2000—2020年, 东北黑土区NPP、水源涵养和防风固沙服务量总体上表现为小幅提升趋势, 分别有96.6%、76.9%和75.4%的县域生态系统服务量表现为上升趋势, 内蒙古东部和黑龙江的供给服务量、黑龙江和吉林的水源涵养量以及辽宁的防风固沙量上升趋势较为明显。土壤保持量总体表现为小幅下降趋势, 60.7%的县域土壤保持量下降, 吉林下降趋势最为明显。

(2) 东北黑土区供给服务变化的驱动类型以气候变化与人类活动的共同促进作用为主。人类活动是东北黑土区供给服务变化的主导驱动因素, 总体贡献率达到了62.65%, 县域数量占76.6%。

(3) 东北黑土区供给与水源涵养、供给与防风固沙服务之间总体上表现为协同关系, 约有88.67%和77.14%的县域表现为协同, 而供给与土壤保持服务之间表现为权衡关系, 县域数量占比为68.25%。人类活动因素对服务间关系的作用表现为加重了权衡、减弱协同作用。