生态系统服务使自然环境与人类福祉相连，同时也为相关政策制定提供基本框架^[1]。生态系统服务类型繁多，由于气候和土地利用/覆盖的变化，以及生态系统服务间直接和间接的相互作用，导致当前全球约60%的生态系统服务正在退化，对全球生态安全以及人类福祉产生重要威胁。不同生态系统服务间存在权衡或协同关系，权衡通常会增强1个生态系统服务，导致其他生态系统服务损耗^[2]。与权衡相反，协同关系是2个生态系统服务同时增强或减弱的情况，深入了解生态系统服务与驱动因素间的复杂关系，为实施可持续管理战略奠定基础^[3]。

当多个生态系统服务受相同驱动因素时，生态系统服务之间的关系就会发生变化；同时，生态系统服务之间权衡和协同关系在很大程度上受政策干预和环境变异性影响。例如，旱地过度的造林计划威胁到生态系统的可持续性^[4]；联合国发布的评估报告中指出24个生态系统中有15个正在持续退化，且退化趋势加剧^[5]。因此，探索生态系统服务的权衡与协同关系及其驱动因素成为研究热点。通常识别生态系统服务驱动因素的方法包括相关性分析^[6]、主成分分析和回归分析^[7-8]。王凯平等^[9]认为人口、国内生产总值等社会经济因素对生态系统服务有重要影响。针对生态脆弱的北方农牧交错带的研究大多从大尺度进行分析，而生态系统服务的尺度关联性使不同区域的同种生态服务功能存在较大差异，大尺度上的评估结果可能在小空间尺度地域并不适用。

小滦河流域位于阴山北麓—浑善达克沙地国家重点生态功能区，是京津冀水源重要供给区，承载着京津冀地区的涵养水源、水土保持、防风固沙等重要生态功能。近年来，伴随着气候变化，小滦河流域水资源供需矛盾日益突出，严重制约着京津冀地区的用水安全。因此，笔者以小滦河流域为研究对象，采用InVEST及RUSLE模型对其2000—2020年的土壤保持、水源涵养、生境质量、水质净化4项生态系统服务功能进行评估，以期厘清各服务间权衡与协同效应的表现类型及时空差异，并探讨生态系统服务的驱动因素。研究结果可为京津冀生态系统服务评估和管理提供理论参考。

1 研究区概况

小滦河流域(E 116°32′~118°14′，N 41°35′~42°40′)位于河北省承德市围场满族自治县境内(图 1)，流域总面积约1 646.42 km²，海拔909~1 871 m。小滦河是滦河上游的一级支流，北起围场县塞罕坝机械林场，南至围场县石桌子乡，东北部为以山地森林为主，向西南高原地势平坦过渡。研究区为半干旱半湿润的大陆性季风气候，年均温5.2 ℃，冬季寒冷干燥，夏季多雨。降水时空分布不均，年均降水量439 mm，多集中在6—9月，蒸发量991.8 mm。土壤类型主要有灰色森林土、风沙土、沼泽土、草甸土等。河北省围场县地处内蒙古高原东南缘与冀北山地的过渡带，是典型的农林牧交错区，流域地表分布着较多固定及移动沙丘，生态环境极为脆弱。

2 材料与方法 2.1 数据来源与处理

数据来源详见表 1，关键参数及数值见表 2。所有数据均通过ArcGIS10.3进行定义投影、掩膜提取、裁剪、重采样、重分类等，栅格数据分辨率统一为30 m×30 m，投影坐标系统一为WGS_1984_Albers。气象数据基于Anusplin4.2进行插值。

1) 蒸散系数K_c范围在0~1.5之间，参考指导手册3.9.1中自带数据库进行赋值。

2) 植被最大根系深度，mm：无植被覆盖区域取值为1，参考相似区域文献^{[17, 19]}及指导手册3.9.1自带数据库进行赋值。

3) 水土保持措施因子P，参考刘灿均等^[10]研究，将耕地、林地、灌丛地、草地、水域、建设用地和未利用土地分别赋值为0.3、0.9、1、1、1、0、0、1。

4) 氮过滤系数、氮输出负荷，kg/hm² ·a、氮养分滞留临界距离，m；参考InVEST模型手册3.9.1及吴瑞等^[11]、杨荣金等^[12]的研究，并进行本土化后赋值。

2.2 生态系统服务评估 2.2.1 水源涵养

水源涵养量(water retention，W_R)是生态系统服务的重要环节，也是生态用水量的主要标志。基于InVEST模型计算出产水量后，需再用地形指数、土壤饱和导水率和流速系数对计算结果进行修正，以获得最终水源涵养量。

$ \begin{gathered} W_{\mathrm{R}}=\min \left(1, \frac{249}{V}\right) \min \left(1, \frac{0.9 \times T_{\mathrm{I}}}{3}\right) \\ \min \left(1, \frac{K}{300}\right) \times Y；\end{gathered} $

(1)

$ Y_{x j}=\left(1-\frac{\mathrm{AET}_{x j}}{\operatorname{PRE}_x}\right) \operatorname{PRE}_x 。$

(2)

式中: W_R为单位栅格水源涵养量，mm；V为流速系数；T_I为地形数，由DEM处理得到，量纲均为1；K为土壤饱和导水率，mm/d，由spaw软件结合土地利用得出；Y为产水量，mm；通过式(2)得到：Y_xj为地类j栅格单位x的产水量，mm；PRE_x为栅格x的年均降水量，mm；AET_xj为j地类栅格x的年均蒸散发量，mm。

2.2.2 土壤保持

土壤保持(soil retention，SR)是生态系统防止自然或人为因素引起土壤流失而侵蚀调控能力及对泥沙储积的保持能力。笔者利用修正的通用土壤流失方程(RULSE)进行估算土壤保持作用。计算式为：

$ \begin{gathered} \operatorname{RKLS}_x=R_x K_x L S_x ; \\ \operatorname{USLE}_x=R_x K_x L S_x C_x P_x ; \\ \operatorname{SEDR}_x=R_{\mathrm{KLS} x}-U_{\mathrm{SLE} x} \end{gathered} $

(3)

式中：RKLS_x为基于地貌气候下栅格x的潜在土壤侵蚀量，t/km²；USLE_x为栅格x实际土壤侵蚀量，t/km²；SEDR_x为栅格x的泥沙持流量，t/km²；R_x为栅格x降雨侵蚀力，MJ ·mm/(km²h ·a)；K_x为栅格x土壤可蚀因子，t ·km² ·h/(km² ·MJ ·mm)；L为坡长系数；S为坡度系数，量纲均为1；C_x为栅格x的覆盖和管理因子；P_t为栅格x的措施因子；二者量纲均为1。

2.2.3 生境质量

生境质量(habitat quality，HQ) 是指为个体或群体提供在生态系统中生存的能力，是支持其他生态系统服务的基础。笔者选择耕地、城市用地、铁路和不同道路作为威胁因素。公式为：

$ Q_{x j}=H_j\left(1-\left(\frac{D_{x j}^z}{D_{x j}^z+K^z}\right)\right) 。$

(4)

式中：Q_xj为生境类型j中栅格x的生境质量指数；H_j为土地利用类型j的生境适宜性；z为归一化常数，取值2.5；K为半饱和常数，取值0.5；D_xj^z为土地利用类型j中网格x的生境退化程度。

2.2.4 水质净化

氮输出(nitrogen export，NE)主要基于InVEST模型的水质净化模块(NDR，Nutrient delivery ratio)中得到，通过植被和土壤可以通过储存和转化拦截径流中的污染物概念来计算每个网格中的氮输出含量。氮产量可表征水净化服务，输出量越高，则证明其水质净化能力越弱，反之亦然。所有像素上氮输出量之和则为整个流域的氮输出量。计算式如下：

$ \begin{gathered} X_{\text {export }}=\text { Load_surf }_i \times \mathrm{NDR}_{\text {surf }, i}+\operatorname{load}_{\text {subs }, i} \times \mathrm{NDR}_{\text {subs }, i} \\ X_{\text {export tot }}=\sum\limits_i X_{\text {export }, i} ; \end{gathered} $

(5)

$ \operatorname{ALV}_x=\operatorname{HSS}_x \times \operatorname{pol}_x \text { 。} $

(6)

式中：X_export为流域中氮输出养分总量，t/(hm² ·a)；load_surf，i为地表氮输出的负荷量；NDR_surf，i为地表氮传输率; Load_surf_i表示地下氮总负荷量，t/(hm² ·a)；NDR_subs，i为栅格单元上地表氮传输率，量纲为1；X_{export tot}为研究区划分子流域的氮养分总输出量，t/(hm² ·a)；X_{export, i}为单位像素i上的输出量，t/(hm² ·a)；ALV_x为像素点x的调节负荷值；HSS_x为像素点x计算方法的水文敏感性得分，pol_x为像素点x的输出系数，量纲为1。

2.3 生态系统服务的权衡/协同关系

通过SPSS和ArcGIS 10.3软件进行Pearson积矩相关系数检验判定权衡协同关系。正系数表示协同关系，负系数表示权衡关系，绝对值越大代表协调或权衡关系越显著，相关系数的显著性水平设置为0.01。

2.4 生态系统服务驱动因素分析

为探讨自然和社会经济因素对生态系统服务的影响，笔者选取13个解释变量进行冗余分析(RDA)(表 3)。在进行RDA分析之前，引入人类活动强度模型^[13]对生态系统受人类活动影响程度进行量化。

$ \mathrm{HAI}=\sum\limits_{k=1} \frac{S_k C_k}{S} $

(7)

式中：HAI为人类活动强度；S为研究区的总面积，km²；S_k为用地类型k的面积，km²；C_k为用地类型k的人类活动强度系数，量纲为1。

3 结果分析 3.1 不同土地利用类型生态系统服务变化

小滦河流域土地利用结构变化差异明显(表 4)。2000—2020年间，研究区耕地转入和转出面积基本持平；草地面积共增加85.37 km²；林地主要向草地转移；灌丛因为向林地和草地转换，其面积减少27.9 km²；水域和居民、城乡建设用地面积在不同地类的转化中均增长，面积增幅分别为31.89%和69.34%；未利用土地减少52.59 km²，降幅达到33.54%。

2000—2020年，林地、水域和草地的生境质量指数较高，林地的生境质量指数>0.9。林地土壤保持、水源涵养服务功能在2000—2020年间波动上升，草地、耕地水源涵养功能骤增。未利用地的水源涵养能力下降明显，土壤保持能力则以2005年为节点，呈现先增强后减弱变化。整体来看(图 2)，各生态系统服务功能的排序为：土壤保持服务>生境质量>水源涵养>水质净化。20年间，小滦河流域城镇化扩张和社会经济活动的加强使居民、建设用地迅速扩张，同时周边林、草、灌丛地面积有不同程度减少，改变了研究区土地利用的空间结构，从而影响不同生态系统服务间供给能力的强弱和相关转化。

3.2 不同生态系统服务时空演变特征

2000—2020年小滦河流域水源涵养、土壤保持、生境质量与水质净化4项服务总体呈上升趋势(图 3)。生境质量指数由0.77增至0.79，空间上增加区域主要位于流域中部及河流沿岸人类活动干扰强烈的区域。城镇化的快速发展使小滦河流域部分林、草地向建设用地转化，水土保持能力有所降低，土壤保持服务以2010年为时间段先降后升，减弱的区域主要分布在流域中、南部。水源涵养服务减弱区域主要集中在流域中部草牧场区域。水质净化服务在20年间增长3.58%，2010年流域水质净化服务能力最强，下降区域主要集中河流下游沿岸，主要由于耕地转化为林地、灌丛地，这两种地类都具有极强的固氮能力，流域中北部区域水质净化显著增强，是由于农村周边居民点的垃圾堆放、牲畜粪便堆积以及耕地施肥的非点源污染所致。

3.3 生态系统服务间权衡与协同关系

不同生态系统服务间的权衡与协同关系用Pearson相关系数来衡量。在2000—2010年间，水源涵养服务与其他3项服务相关系数为负值(表 5)，证明各服务之间存在明显权衡关系。而土壤保持与生境质量、水质净化服务间相关系数为正值，证明这些生态系统服务之间存在明显协同关系。2010—2020年间水源涵养服务与其他3项服务相关系数转变为负值，表明不同生态服务间的权衡与协同关系具备转换性。生境质量与土壤保持服务间的正相关系数最小，权衡关系稳定中呈现出小幅增强的趋势。以2005年为拐点，水源涵养与土壤保持、生境质量服务间的权衡关系从减弱转为增强的协同关系；而土壤保持与水质净化服务间的协同关系转为权衡关系。水源涵养与水质净化服务间关系不稳定，2010年为协同关系，其余年份则呈现权衡关系。2000—2020年，生境质量和水质净化的相关系数均最高，证明其相关性最紧密。此外，随着时间推移，土壤保持和生境质量服务间呈现递增的协同关系。

3.4 生态系统服务驱动因素分析

由于研究时间段内各因素变化波动甚微，因此笔者重点选取2000年和2020年进行冗余分析。

不同驱动因子对不同生态系统服务影响程度不同，锐角表明变量之间呈正相关，夹角越大，协同影响越弱；钝角表明变量之间呈负相关，夹角越大，权衡影响越弱；直角时，无相关关系；箭头长短代表该因子对响应变量的贡献率。由图 4可见，在2000年和2020年，小滦河流域生态系统服务与自然、社会因子间存在极显著相关性(P < 0.001)，社会和自然因子对生态系统服务的累积解释比例为56.2%和63.9%。2000年，未利用地比例、土壤有机质和坡度对4种生态系统服务的相关性较大，而林地面积比例、耕地面积比例及人口密度则显著相关性较小。2020年时，除土壤有机质、坡度、降水、草地面积比例外，其余影响因子对4种生态系统服务的贡献率有所降低，降水、坡度、NDVI指数、林地比例对水源涵养、土壤保持、生境质量服务表现为协同影响，对水质净化服务表现为权衡影响。土壤有机质含量与水源涵养呈权衡相关，与土壤保持、生境质量、水质净化服务呈协同相关；坡度和降水与水质净化呈权衡相关，与土壤保持、生境质量、水源涵养服务呈协同相关。同时林地、草地、耕地、未利用地对各生态系统服务的解释率较高，证明城镇化的发展引起建设用地扩张，进而破坏生态系统服务之间关系的平衡，但林草地和耕地面积的增加可以有效维持其关系稳定。

4 讨论

笔者在评估流域生态系统服务时空演变时，采用Pearson系数分析研究区生态系统服务的权衡与协同关系。结果表明，小滦河流域4种生态系统服务的权衡与协同关系共存，权衡协同关系之间相互转化。其中，水质净化与水源涵养和生境质量呈增强的权衡关系，与土壤保持呈现协同关系向权衡关系转变，可能是由于建设用地挤压生态用地，两者矛盾不断增多。此外，笔者发现小滦河流域生态系统服务间的权衡/协同关系可能在不同时期发生变化^[14]。2000—2020年间，小滦河流域水质净化与土壤保持、水源涵养协同效应减弱，权衡效应增强。以往研究表明，生态系统服务间权衡/协同关系较弱时，权衡/协同关系很容易发生相互转换^[15]。生态系统服务间关系的演变往往与区域自然条件和社会经济发展有关。在小滦河流域，水质净化和土壤保持间的关系从协同效应转变为权衡关系，主要是由于植被覆盖度显著增加。小滦河流域地处农牧交错带，耕地承担着水源涵养和粮食产量的双边责任，因此不适宜大面积规划退耕还林政策，还需要根据流域经济发展实际不断调整。

通过RDA分析发现，地形因子一定程度上决定着土地利用的分布和功能，大量研究证实，土地利用结构的变化和植被覆盖度的增加可引起生态系统的结构和功能发生变化^[16]。RDA结果还表明，土壤有机质对小滦河流域生态系统服务有重要影响，气候、土壤、人为等多因素共同作用造成生态系统服务的空间异质性。植被覆盖、耕地比例和未利用地比例对小滦河流域生态系统服务产生重要影响，其中NDVI的重要性不断增强，而耕地与未利用地比例在下降，主要原因在于随着植树造林及生态工程建设的推进使流域植被覆盖度增加，对应耕地及未利用地则不断减少。建议小滦河流域应进行长期的土地规划和配置优化，加强对流域内林地和草地资源的保护工作，实行弹性放牧积极恢复草地生态系统的生产力，有效提升小滦河流域水源涵养、土壤保持力及生境质量，减少氮输出量，促进流域内生态用地的良性转换和生态系统的稳定性。本研究还发现，降水的影响在不断下降，降水作为干旱半干旱区植被需水量的主要来源，是影响生态系统服务的关键因素^[17]，但在适当的工程措施的配合下，可提高区域的生态系统服务功能^[18]。

本研究在模型参数本土化的基础上成功运用在小滦河流域，对其4项生态系统服务表述较为清晰。但是仍存在部分重要参数获取方法单一，模型计算过程未考虑更多实际复杂因素等不足，结果精度仍有提升空间。

5 结论

1) 研究区不同土地利用结构的生态系统服务功能差异明显，林地、草地的生境质量指数最高，耕地的水质净化服务供给能力最强，未利用土地的土壤保持能力呈先增后降的变化趋势。城镇化和社会经济活动的增加改变了土地利用结构的变化，影响不同生态系统服务间供给能力强弱的转化。

2) 2000—2020年小滦河流域4项服务均呈增加趋势，水源涵养服务与其他3项服务存在明显权衡关系，土壤保持与生境质量、水质净化服务存在明显的协同关系。生态系统服务间权衡与协同关系的演变往往与区域自然条件和社会经济发展有关。

3) 不同驱动因子对各生态系统服务的影响程度不同，小滦河流域生态系统服务与自然、社会因子间呈极显著相关。降水、NDVI指数、林地比例对水源涵养、土壤保持、生境质量服务表现为协同影响，对水质净化服务表现为权衡影响。城镇化的发展破坏生态系统服务之间关系的平衡，但林草和耕地面积的增加可以有效维持其关系稳定。