李 婷,乔志宏,冯玮含,杨亚东 (.西安科技大学测绘科学与技术学院,陕西 西安 70054;.中国农业科学院农业资源与农业区划研究所,北京 0008)

旱地覆盖陆地表面41%的土地面积,支撑着全球超过20亿人口的生存与发展[1]。旱地生态系统不仅是草料生产和畜牧发展的基础,而且具有气候调节、防止土壤侵蚀和维持生物多样性等多项生态系统功能,由此成为开展生态恢复的关键区域[2]。随着全球人口持续增长及畜牧产品需求不断提升,预计到21世纪中叶,放牧活动频率及放牧强度增加的风险将持续提高。放牧干扰与水资源短缺、气候变化等因素交互叠加,造成旱地生态恢复成效的多重不确定性[3-4]。因此,有必要明确不同放牧强度对生态系统服务的影响,为旱地生态系统恢复与可持续牧业管理提供关键信息。

放牧牲畜通过采食、践踏和代谢物排放改变植物群落结构,进而影响生态系统生产力和其他功能过程[5]。过度放牧会对多项生态系统服务产生负面影响,造成植被和土地功能退化。PETZ等[6]证实高强度放牧将导致生物多样性和土壤防侵蚀能力降低,且碳排放量比低度放牧高4倍以上;FAN等[7]通过内蒙古草原的放牧试验发现,逐渐增大的放牧强度将减少草原土壤蓄水,增加土壤侵蚀风险。但合理的放牧活动在生物多样性和遗传资源维持、动植物产品供给等方面产生着积极作用[8]。例如,基于旱区地形梯度下放牧强度的情景模拟表明,适度放牧可促进水分胁迫下低山草地和平原草地植被净初级生产力显著提高[9]。针对黄土丘陵区典型坡面的研究结果也显示,适度放牧有助于增加地表枯落物盖度,促进坡面土壤水分入渗,并能够改善退耕草地植被群落发育状况[10-11]。综合来看,放牧活动对生态系统的影响随环境梯度和放牧压力而产生变化,因此放牧强度估算是开展生态系统功能响应研究的必要前提[12]。

以往关于放牧强度的量化研究以小尺度控制试验居多,即通过控制牧场面积、牲畜数量和牧草利用率等指标直接量化放牧强度[13],并开展放牧强度对植被群落特征、生产力、物种多样性等影响研究[10,14]。尽管控制试验可以相对准确地量化放牧强度,但难以充分反映大范围放牧引起的地面植被变化。近年来,遥感数据因其大范围、长时序、高时效等优势成为进行大尺度放牧强度估算的重要手段,利用归一化植被指数(normalized difference vegetation index,NDVI)、地上生物量(aboveground biomass,AGB)和地上净初级生产力(aboveground net primary productivity,ANPP)等遥感植被指数与实测数据建立函数模型是研究区域尺度放牧强度的常见方法。例如,邵明轩等[15]利用实测AGB数据与NDVI建立反演模型,并通过分析生长季NDVI变化特征开展高寒草原放牧强度估算;翟星等[16]则以NDVI数据反演了区域尺度ANPP空间分布,并结合气象数据与实测AGB估算了若尔盖高原放牧强度。然而,基于植被指数与实测数据建立的放牧强度经验模型在不同地区的适用性难以确定。对于生态恢复地区,植被变化直接受到放牧活动的干扰,并影响着生态系统服务的产生与传递[17]。如何快速、有效地确定生态恢复以来的放牧水平仍有待进一步研究;此外,不同放牧强度下生态系统服务的响应方式尚不清晰,这也是生态恢复区开展适度放牧面临的挑战之一。

2014年起,我国进入第二轮退耕还林还草阶段。由于农牧民所获直接经济补偿减少,伴随生活生产成本提高,造成偷牧活动禁而不止,形成禁牧—偷牧—土地退化的恶性循环,阻碍了生态恢复效益的实现和维持[18]。但长期封禁也带来草地生态系统多样性降低、生产力下降等新的生态问题;部分封禁地植被及其生态功能恢复已进入相对缓慢阶段,从长远来看将造成自然资源的浪费[19-20]。基于上述背景,以黄土高原为研究区,对2000—2020年研究区产水量、草料供给和土壤保持3项关键生态系统服务空间分布与变化特征进行量化;基于净初级生产力(NPP)和全球牲畜密度数据量化放牧强度,并通过提取约束线确定放牧强度对这3项服务的制约作用;在此基础上探究不同放牧强度下生态系统服务权衡/协同时空变化特征。研究结果可为旱区生态恢复与可持续牧业管理提供科学依据。

1 研究区与数据来源

1.1 研究区概况

黄土高原位于黄河中游地区〔图1(a)〕,包括宁夏和山西全境,陕西中北部,甘肃陇中和陇东地区,青海东北部以及内蒙古河套平原和鄂尔多斯高原,面积约64万km2;海拔分布于1 200～1 600 m。研究区属于半湿润向半干旱、干旱气候区的过渡地带,降水事件多集中于6—9月,年平均降水量分布自东南向西北内陆递减〔图1(b)〕,多年平均气温为9 ℃左右〔图1(c)〕。严重的水资源匮乏和水土流失是黄土高原可持续发展的两大障碍。行政单元上,黄土高原横跨7省(区),其中,青海省、甘肃省、宁夏回族自治区和内蒙古自治区是黄土高原畜牧产品供给的主要产区,山西省和陕西省是黄土高原主要粮食及其他农产品供给区。但在农业区划上,黄土丘陵区大面积地区包含于农牧交错带,畜牧发展在农民经济收入中占据重要地位,放牧干扰对黄土高原植被恢复影响实际仍然存在。退耕封禁20年来,黄土高原地区土地利用/覆被变化以林地、草地和人工用地面积增加为主要特征,植被盖度显著提升,生态环境明显改善。但如何有序、科学地开展适度放牧以维护生态恢复成效的持续性和稳定性已成为该地区重要的研究议题。

图1 研究区概况Fig.1 Study area overview

1.2 数据需求及来源

放牧强度及生态系统服务量化涉及牲畜密度、遥感植被指数以及地形因子、土壤属性、气候因子等数据。主要数据包括:(1)2000、2010和2020年30 m分辨率土地利用/覆被数据来源于中国科学院资源环境科学与数据中心(https:∥www.resdc.cn/);(2)NPP、植被覆盖度(FVC)为MODIS产品数据集,由中国科学院遥感与数字地球研究所提供,空间分辨率为250 m,时间分辨率为16 d;(3)黄土高原及其周边113个气象站台数据来源于国家气象信息中心(http:∥data.cma.cn/),用于各时期年均降水量空间插值和降雨侵蚀力计算;(4)1 km分辨率土壤属性数据来源于中国西部环境与生态科学数据中心(http:∥westdc.westgis.ac.cn)提供的世界土壤数据库(harmonized world soil database,HWSD)——中国土壤属性数据集,用于估算土壤可蚀性、土壤深度、植物可利用含水量等生态参数;(5)30 m分辨率DEM数据来源于地理空间数据云(http:∥www.gscloud.cn/),用于子流域提取、耕地坡度提取及坡度坡长因子等计算;(6)不同土地利用/覆被类型蒸散相关系数参考联合国粮食及农业组织(FAO)在线目录(http:∥www.fao.org/docrep/X0490E/x0490e0b.htm)取值;(7)1 km分辨率牲畜密度数据来源于FAO世界牲畜格网数据集(https:∥dataverse.harvard.edu/dataverse/glw),用于放牧强度量化。

2 研究方法

2.1 生态系统服务量化

2.1.1产水量

采用InVEST模型中Water Yield模块进行产水量计算。根据Budyko曲线和年平均降水量,Water Yield模块通过年降水量减去年实际蒸散发量计算每个像素〔Y(x)〕的年产水量,评估各子流域对区域产水量的相对贡献,计算原理如下:

(1)

(2)

PET(x)=Kc(lx)·ET0(x),

(3)

(4)

式(1)～(4),AET(x)为栅格单元x年实际蒸散量,mm;P(x)为栅格单元x年降水量,mm;PET(x)为栅格单元x年潜在蒸散量,mm;Kc(lx)为栅格单元x植物蒸散系数;ET0(x)为栅格单元x参考作物蒸散发量,mm;ω为表示自然气候-土壤性质的非物理参数,由植物可利用水含量和反映降水季节分布的经验常数确定;Z为反映降水季节分布的经验常数;AWC(x)为土壤有效含水量,由植物可利用水分含量、最大根系埋藏深度和植物根系深度的最小值确定。笔者研究根据黄土高原及周边113个气象站点数据,通过空间插值生成不同年份降水量栅格图层;ET0(x)采用修正后Hargreaves方程[21]进行量化得到;土壤最大深度由HWSD数据库获取,植物可利用水分含量根据HWSD中砂粒、粉粒、黏粒和土壤有机质进行估算得到[22]。Kc(lx)和植物根系深度参考FAO在线目录对土地利用/覆盖类型进行赋值,并参考黄土高原相关研究[23]进行修正,最终数据以生物物理表格形式输入。Z分布范围为1～30,由各年份降水事件次数经验公式[24]计算得到。

2.1.2草料供给

产草量是衡量草地生产力高低的重要指标,也是制定畜牧业生产和管理规划的基础。除草地外,牲畜对灌丛同样具有较高的采食积极性,因此灌丛也是牲畜重要的草料来源[25]。FLOMBAUM等[26]通过对阿根廷干旱生态系统开展实地研究,建立干旱环境中草地和灌丛绿色生物量(Gb,g·m-2)与FVC的线性回归方程(表1[26])。其中,Gb为草地所有绿叶以及灌丛所有绿叶和当年生枝条,是一种快速估算大尺度生态系统草料供给的有效指标,并在旱地生态系统相关研究[27-28]中得到借鉴。笔者亦采用该方法,根据2000、2010和2020年土地利用/覆被数据二级分类和FVC数据,提取草地与灌丛斑块FVC;利用ArcMap 10.3平台进行不同时期草料供给的线性回归计算,最终形成各时期研究区草料供给服务空间格局。

表1 灌丛及草地绿色生物量(Gb,g·m-2)与植被覆盖度(FVC,CFV,%)线性拟合关系[26]Table 1 Linear fitting relationship betweenthe green biomass (Gb) and FVC (CFV) in shrubs and grasses

采用SUN等[27]建立的全球草地生态系统NPP实测样点数据集对草料供给量化结果进行验证。该数据集包含1957—2018年全球草地生态系统地上净初级生产力(ANPP)、地下净初级生产力(BNPP)和总净初级生产力(TNPP)共2 985条实测数据。其中,ANPP代表生态系统地上生物量,也是反映生态系统产草能力的直接指标[28]。因此,采用2000和2010年实测ANPP数据对草料供给服务结果进行检验。由于草料供给服务仅涉及草地和灌丛斑块,排除空值样点后,最终筛选保留黄土高原地区28个样点进行结果验证〔图2(a)〕。结果显示,2000和2010年绿色生物量(Gb)与实测ANPP的决定系数(R2)均大于0.5〔图2(b)～(c)〕,表明该方法能够较好地反映研究区草料供给功能。

Gb为绿色生物量,ANPP为地上净初级生产力。图2 ANPP实测样点空间分布及数据验证Fig.2 Spatial distribution of measured sample points of ANPP and the data verification

2.1.3土壤保持

采用InVEST模型Sediment Delivery Ratio模块量化年均潜在土壤侵蚀,并以潜在土壤流失量乘以植被覆盖度因子(C)和控制实践因子(P)得到实际土壤流失量,最后通过计算两者之差得到各时期年均土壤保持量,计算公式为

S=R·K·LS·(1-C·P)。

(5)

式(5)中,S为年均土壤保持量,t·hm-2·a-1;R为降水侵蚀力因子,MJ·mm·hm-2·h-1·a-1,基于月平均降水量计算多年平均降雨侵蚀力,并乘以17.02进行国际制单位转换[29];K为土壤可蚀性因子,t·hm2·h·MJ-1·hm-2·mm-1,根据EPIC模型进行计算,并乘以0.131 7转换为国际制单位[22];LS为坡度坡长因子,根据在模型中输入无洼地DEM运算得到;C和P分别为植被覆盖因子和水土保持因子,分别基于植被覆盖度[30]和耕地坡度方法[31]计算得到;InVEST模型的汇流累积量是形成河网的阈值条件,汇流累积阈值根据黄土高原实验样区分析结果[32]设定。

2.2 放牧强度量化

考虑牲畜活动对灌草食用强度的影响,参考PETZ等[6]对大尺度放牧强度量化思路,以放牧生物量与牲畜可食用最大生物量之比量化研究区禁牧后放牧强度的空间分布:

I=D/A。

(6)

式(6)中,I为放牧强度指数;D为放牧生物量,取决于牲畜分布密度;A为可供放牧的最大生物量,定义为在最高放牧强度下地上NPP平均百分比,取值为64%[33],但不同生态系统之间及内部生物群落可食用性存在普遍差异,后续研究多以地上生态系统NPP的60%进行估算,并使用多年NPP平均值以避免年际波动引起的误差[33-34]。

FAO世界牲畜格网数据集(gridded livestock of the world,GLW)[35]提供了全球及国家尺度牛、绵羊、山羊和马等牲畜的空间分布密度数据,数据更新时间为2015年。GLW数据产品的推出为开展大范围土地利用强度格局[36]、放牧牲畜影响因素[37]以及草地生态系统承载力[33]等研究提供了重要支撑。已有研究利用中国省级和地市级农业主管部门山羊和绵羊统计资料对GLW数据集进行验证,证实该数据集对国家尺度放牧牲畜分布及预测研究具有良好的适用性[37]。笔者根据黄土高原地市级和县区级统计资料,选取牛、羊和马作为产生放牧干扰的主要牲畜类型。鉴于数据可获取性和县区边界完整性,收集研究区2015年194个县区牛存栏和羊存栏统计值用于GLW牲畜放牧密度分布验证。由GLW获取的放牧牛和放牧羊(山羊和绵羊)空间分布密度具有明显差异(图3),牛高密度放牧区集中于青海省州县,羊高密度放牧区则主要分布于内蒙古自治区各旗县、山西省和陕西省北部县区。GLW两种放牧牲畜密度数据与县区统计值的回归系数分别达到0.774和0.616,R2分别为0.718和0.638(图3),说明GLW数据集能够有效反映研究区放牧牲畜空间分布的梯度差异。在此基础上,综合GLW牛、马和羊(山羊和绵羊)空间放牧密度数据用于放牧强度的量化。

图3 GLW放牧牛、放牧羊(山羊和绵羊)分布密度及其与统计值对比Fig.3 Spatial distribution on the grazing density of cattle and the total of goat and sheep, and their comparisons with the statistical data at county level

为实现不同类型牲畜的单位换算,参考FAO《畜牧业发展与减贫》研究报告[38],采用国际通用热带牲畜单位(TLU)将各类牲畜折算为统一标准,1 TLU等于250 kg,每头牛或马换算为0.6 TLU,每只山羊或绵羊换算为0.1 TLU。最终,以牲畜密度预处理结果和2000—2020年平均NPP数据(取值为60%)量化截至2020年研究区放牧强度指数(I),并参照PETZ等[6]研究结果将I划分为低(0～0.4)、中(>0.4～0.6)和高(>0.6～1.0)3种类型,用于不同放牧强度下生态系统服务关系的进一步分析。其中,低度放牧表示生态系统受到放牧生产干扰极小,畜牧业生产消耗远低于自然生态系统生产力;中度放牧表示放牧活动遵循季节规律,植被结构保留原始生态系统结构和物种组成,且畜牧业生产消耗略低于自然生产能力;而高度放牧表示集约化的放牧活动,严重依赖原始生态系统的饲料投入,且畜牧业生产消耗接近或超过自然生产力。

2.3 放牧强度对生态系统服务的约束效应

基于分段分位数回归的约束线提取能够充分表征限制变量对响应变量的主要约束过程,已成为探索生态系统服务之间约束关系及地理要素对生态系统服务影响的成熟方法[39-40]。根据约束线信息属性,约束线上方的点表示限制变量与响应变量之间很少或几乎不受其他因素影响,约束线下方的点则表示限制变量与响应变量之间仍存在其他因素影响[40]。以I作为限制变量,将2000、2010和2020年产水量、草料供给和土壤保持3项生态系统服务作为响应变量,采用约束线提取方法对比退耕后不同时期放牧强度对生态系统服务的制约影响。具体方法为:(1)对限制变量即放牧强度进行标准化处理,将取值范围平均分成100份作为x轴100列数据;(2)为减少异常值的影响,采用MATLAB R2020a软件提取每列数据99%分位数值作为边界点;(3)对边界点进行趋势线拟合,并根据散点云图分布形状和分位数拟合优度(R2)确定放牧强度对生态系统服务的最优约束线类型。

2.4 生态系统服务之间权衡/协同关系分析

根据放牧强度和生态系统服务量化结果,以不同放牧强度区分别对3项生态系统服务进行空间采样取值,采用“3σ准则”排除异常值[41],并剔除草料供给服务空值点,作为各时期服务权衡/协同关系分析样本数据;其次,对各项生态系统服务进行无量纲化处理,并采用R v4.1.2软件包Performance Analytics工具进行生态系统服务之间的Pearson相关分析,开展不同时期、不同放牧强度下生态系统服务权衡/协同变化研究。

3 结果与分析

3.1 放牧强度空间格局及牲畜承载量

以县区放牧牲畜统计数据验证GLW数据在黄土高原的放牧牲畜空间分布密度,在此基础上量化研究区放牧强度指数。从放牧强度空间格局上来看,生态恢复20年来黄土高原地区禁牧效果显著。截至2020年,高度放牧区面积仅占研究区总面积的5.30%,主要分布于内蒙古自治区及青海省部分旗或县区,且大部分均为传统畜牧区(图4);中度放牧区面积占比为31.65%,主要分布于山西省和陕西省以北,以及甘肃省和宁夏回族自治区以西范围;研究区其他地区均属于低度放牧区,面积占比达63.05%。

图4 2020年黄土高原放牧强度空间格局Fig.4 Spatial pattern of grazing intensity of the Loess Plateau in 2020

根据牲畜密度及不同牲畜热量的单位换算,不同放牧强度下牲畜活动密度范围见表2。由于地表可食用最大生物量分布具有显著空间差异,研究区单位面积牲畜质量及相应牲畜数量并非连续分布,而受NPP影响呈断点式分布。以山羊或绵羊为例,低度放牧区每km2放羊数量小于100只;中度放牧区每km2放羊数量分布范围为150～200只;高度放牧区每km2放羊数量通常大于400只,如研究区西部西宁市部分州县区(图4)。

表2 不同放牧强度单位面积牲畜承载范围Table 2 Livestock carrying range per unit area under different grazing intensity

3.2 放牧强度对生态系统服务的约束效应

3.2.12000—2020年生态系统服务空间分布及变化

2000—2020年单位面积年平均产水量、草料供给和土壤保持均有不同程度增长,其中,产水量经历先增后减的变化;草料供给服务增长更为明显,至2020年增幅达37.12%;植被恢复措施成为降低区域土壤侵蚀的主要贡献因素,研究期间单位面积年均土壤保持量持续增长(表3)。

表3 2000—2020年黄土高原年平均生态系统服务量Table 3 Annual average ecosystem services of the Loess Plateau from 2000 to 2020

生态系统服务时空分布主要取决于降水、河网密度和生态系统类型等自然禀赋,不同时期各项服务空间分布具有一致性特征〔图5(a)～(c)〕。产水量自东南至西北梯度减少,低值区域广泛分布在黄土高原西北部地区。灌丛生物量对区域草料供给贡献更为显著,单位面积草料供给高值区域主要分布在陕西北部及黄土高原东部地区,而不是以草地生态系统为主的传统牧区,如内蒙古自治区或青海省部分旗县。土壤保持服务空间分布主要与植被覆盖度、水土保持措施和土地利用/覆被类型相关,高值区集中分布在陕西北部以及黄土高原东部和南部地区。

图5 2000—2020年黄土高原生态系统服务空间格局及变化Fig.5 Spatial pattern and changes of ecosystem services in the Loess Plateau from 2000 to 2020

为反映2000—2020年各项生态系统服务空间变化,以子流域为分析单元进一步识别各项服务呈现增加或减少趋势的空间分布,并统计其占黄土高原面积的百分比〔图5(d)〕。受生态恢复工程实施和降水分布的影响,区域尺度产水量增长面积占比为60.05%,主要分布于中部和西北部地区。从草料供给空间变化上来看,该项服务改善面积占比达到76.46%,但内蒙古自治区、甘肃省和陕西省南部部分旗县区草料供给功能有所下降。此外,显著的植被恢复使得各子流域土壤保持服务整体有所提高。

3.2.2放牧强度对生态系统服务的约束效应

分段分位数方法可以揭示不同时期生态系统服务受放牧强度的约束效应。图6显示,放牧强度对3项生态系统服务均呈负向约束关系,即各项服务随放牧强度的增加整体呈减少趋势,决定系数(R2)表明提取的约束线具有良好的拟合优度。

红色线为约束线,R2为决定系数。灰色点为数据散点,彩色点为边界点。图6 2000、2010和2020年放牧强度对生态系统服务的约束效应Fig.6 Scatter plots and constraint lines between grazing intensity and ecosystem services in 2000, 2010 and 2020

不同时期放牧强度对各项服务的约束线具有相似性。其中,产水量受放牧强度制约均呈反S型约束曲线。2000和2010年放牧强度-产水量关系的拐点值接近(表4),放牧强度分别小于14.30%和18.05%时产水量随之增加,且产水量第1拐点均大于500 mm,放牧强度超过上述阈值后产水量随之减少,产水量下降至约275 mm时出现第2拐点。2020年,放牧强度第1拐点值提前至8.86%,且产水量拐点值分别降低至440.55和250.32 mm。放牧强度增加对草料供给产生显著且单调的负线性制约效应。土壤保持服务受放牧强度影响则均呈指数递减约束曲线,超过放牧强度-土壤保持阈值后曲线负斜率增加,表明放牧强度高于约40%时,土壤保持量更易受放牧强度制约而减少;并且在放牧干扰下,2020年土壤保持服务阈值下降至150.81 t·hm-2,表明随着时间尺度延长,放牧强度对土壤保持服务的制约作用将更加明显。

表4 放牧强度对产水量及土壤保持服务约束阈值Table 4 Thresholds of the constraint lines of grazing intensity on water yield and soil conservation

3.3 不同放牧强度下生态系统服务权衡/协同变化

结合放牧强度空间格局与生态系统服务量化结果,以研究区3项生态系统服务权衡/协同变化为参照,进一步开展不同放牧强度下3项生态系统服权衡/协同变化特征分析。结果表明,2000—2020年黄土高原3项生态系统服务整体呈现协同关系,相关系数均通过0.01水平显著检验(表5),但20年间正相关系数均呈减少趋势。其中,产水量与草料供给之间正相关系数由0.60降低为0.47,协同关系削弱程度最大,产水量与土壤保持服务之间正相关系数由0.40减少至0.36;草料供给与土壤保持之间相关系数由0.31波动减少至0.29。

表5 2000、2010和2020年黄土高原及不同放牧强度下生态系统服务相关系数变化Table 5 Correlation coefficients of ecosystem services of the Loess Plateau in the different grazing areas in 2000, 2010 and 2020

高度放牧区3项服务之间相互关系与区域尺度研究结果具有一致性,即正相关系数呈减少趋势,但其协同关系削弱程度较区域尺度更为剧烈。2000—2020年,产水量与草料供给、土壤保持之间正相关系数由0.34和0.24分别减少至0.16和0.19;草料供给与土壤保持之间正相关系数亦呈波动减少趋势,由0.21减少至0.17。总体来看,高于阈值范围的放牧干扰仍是制约生态系统功能恢复的重要因素(表5)。诸多研究结果也证实,高强度牲畜活动对地表灌草等植被的频繁啃食和践踏将导致表层土壤紧实度增加和孔隙度减少,造成草原土壤蓄水减少和土壤侵蚀风险增加[34]。

中度放牧区3项服务之间相关系数持续增长。其中,2000年产水量与草料供给之间为不显著权衡关系,而2010—2020年间相关系数由0.04增加至0.13;且2000—2020年产水量与土壤保持之间相关系数由0.25增至0.36;草料供给与土壤保持之间正相关系数增长趋势亦较为明显。低度放牧区面积占研究区面积的50%以上,但3项生态系统服务之间权衡关系加剧。其中,2000—2020年产水量与草料供给关系变化最为突出,由协同关系(相关系数为0.29)转变为权衡关系(相关系数为-0.14)。在封禁最佳年限制约下,严格禁牧或低度放牧干扰对生态系统功能的改善已呈现边际效应递减趋势。

4 讨论

黄土高原地区全面推行退耕禁牧至今已20余年。生态恢复工程的实施深刻影响区域生态系统服务空间格局与关系变化。林草地面积增加有助于改善固碳服务,改良土壤物理性质,减少土壤侵蚀[42]。在此影响下,研究区产水量、草料供给和土壤保持整体呈现东南高西北低分布特征,这与前人研究结果[43]具有一致性。但研究区不同放牧强度下3项服务间关系强弱及方向存在明显差异,且服务权衡/协同关系随封禁时间产生变化。

尽管单位面积产水量差异并不显著(表3),但放牧强度对产水量的约束拐点值由2000年的518.42和275.92 mm分别降低至2020年的440.55和250.32 mm(表4),进一步证实区域产水功能下降的趋势,且下降区域主要分布于低度放牧区〔图4,图5(d)〕。作为典型农牧交错区,黄土高原低度放牧区生态系统服务关系的变化更多地是由大面积退耕造林引起的生态效应。人工林草生长阶段伴随的植被冠层蒸发和蒸腾效应的增强,使得区域产水量与碳固存、侵蚀控制间协同作用随时间变化呈现削弱趋势[44]。但典型退耕区产水量和其他生态系统服务之间的权衡关系更加凸显[45]。草料供给服务体现灌草生态系统固碳功能,其与产水量之间的协同关系明显减弱,且低度放牧区产水量与草料供给之间已由协同关系转为权衡关系(表5)。基于长期定位观测的试验[46]同样表明,退耕区群落盖度、物种丰富度和地上生物量等特征指标增长峰值均出现在封禁20年时,随后显现下降趋势。因此,忽略最佳封禁时限将对生态系统自然更新、生物量累积造成一定负面影响,并可能加剧植被恢复带来的水资源利用负担。

放牧强度-产水量的反S型约束关系显示一定阈值范围内产水量将随放牧强度增大而增强(图6)。对于中度放牧区,尽管可能存在局地产水量和草料供给同时减少的情况,但中度放牧区产水量和草料供给均以增长为主〔图4,图5(d)〕,因此生态系统服务空间采样分析结果整体呈现增强的协同关系(表5)。已有研究[47]证实,适度放牧能够使牧草凋落物和立枯物减少,激发牧草超补偿性生长机制。例如,适度放牧减少了过于繁密的植物叶片面积,有利于缓解植被蒸散发,从而降低植被生长所需水分[48]。因此,放牧活动在一定程度上削弱了禁牧最佳时限与可持续水资源利用短缺的消极影响[33]。尽管如此,笔者研究中反S型约束线放牧强度拐点值均小于20%,且2020年放牧强度对产水量约束效应拐点出现更早(表4),表明放牧强度对研究区水资源的限制已更加明显。因此,针对不同区域特征进行放牧活动的适应性管理时,必须考虑可供放牧的最大地上生物量及其空间异质性特征。

现阶段不同研究对于放牧强度的定义和划分标准存在较大差异,仅以强度等级进行结果对比可能会引起偏差。例如,黄土区退化草地试验结果表明封育15年后可开展每km2200只羊的低度放牧活动,以提高草地质量,促进物种多样性[46],与研究区中度放牧区放牧牲畜种类及数量类似,即每km225～33头牛或150～200只羊(表2);也有研究[10]通过羊群粪球密度划分放牧干扰强度,并发现粪球密度低于10个·m-2时可有效促进退耕草地植物群落发育。因此,除考虑放牧牲畜采食生物量的空间异质性,放牧强度的适应性管理还须以实际放牧的牲畜类型及数量来细化指标,以避免由定义或类别划分差异引起的管理偏差。划分标准的差异也使得区域放牧强度难以被直接验证[13]。因此,笔者参考相似研究采用实际统计数据对GLW放牧牲畜密度空间分布进行评估。尽管GLW数据集能够有效反映区域放牧牲畜空间梯度差异,但由于数据算法及时空分辨率缺陷,仍可能对局地放牧密度估算造成一定误差。后续将选择典型样区开展不同类型放牧牲畜的实测校验以降低GLW数据可能存在的误差及不确定性;对可供放牧的最大生物量(A)也将进一步综合地形、温度和植被高度等关键因素进行细化[28,33]。未来,在气候变化与畜牧产品需求的双重压力下,植被恢复与可利用水资源权衡范式仍将不可避免地对旱地生态系统功能关系造成长期影响。面对上述挑战,探索适度放牧模式并发挥放牧干扰的正向作用将是促进旱地生态恢复成效与畜牧业可持续管理的有效措施。

5 结论

该文对退耕后黄土高原地区放牧强度进行估算,通过分段约束线方法明确放牧强度对生态系统服务的制约作用,并探索不同强度放牧下黄土高原生态系统服务权衡/协同关系的时空异质性。结果表明:(1)生态恢复20年以来,黄土高原地区高度放牧区面积仅占研究区总面积的5.30%,低度放牧区面积占比达63.05%。(2)放牧强度对产水量、草料供给和土壤保持服务均呈显著负向约束关系,但在一定阈值范围内产水量随放牧强度增加呈现增加趋势。(3)放牧约束条件下,高度放牧区生态系统服务之间协同关系削弱程度更为剧烈,中度放牧区服务之间协同关系整体呈增长趋势,其中,产水量与草料供给之间由权衡关系转为协同关系,但低度放牧区服务之间权衡关系有所加剧。(4)严格禁牧模式下生态系统功能的改善已呈现边际效应递减趋势。因此,掌握放牧牲畜类型和恢复生态系统空间异质性是开展适度放牧的必要前提。