金海珍,于德永,*,郝蕊芳,黄 婷

1 北京师范大学地理科学学部,北京 100875 2 北京林业大学水土保持学院,北京 100083

生态系统服务是指人类从自然生态系统中获得的直接和间接利益[1],是生态系统所形成并维持的人类赖以生存的自然环境条件与效用[2],包括供给服务(如维持人类生存的粮食、淡水等),调节服务(如洪水调节、气候调节等),支持服务(如维持地球生命生存环境的养分循环)和文化服务(如精神、娱乐和文化收益)[3]。各项生态系统服务间是相互影响的,它们的可持续供给对维持整个生态系统稳定和提高人类福祉具有重要意义。由于生态系统服务的类型丰富,空间分布的异质性和人类对生态系统服务的需求偏好[4],生态系统服务之间形成了复杂的关系。科学认知和理解生态系统服务之间的关系是区域生态系统服务管理和优化的关键[5]。目前,生态系统服务间关系的研究多为二元化方法,即此消彼长的权衡关系与相互促进的协同关系[6-7]。研究内容多包括权衡与协同特征的识别、时空变化的分析和驱动因子的探究等[8-11]。例如,Shen等采用空间叠加和统计方法定性分析京津冀城市群生态系统服务在不同服务簇中权衡与协同的空间分布格局[12]。Raudsepp-Hearne等利用相关系数反映加拿大魁北克城市景观供给服务和调节服务之间权衡关系的强度[13]。Kubiszewski等通过模拟未来情景下土地利用情况,估算生态系统服务的未来价值及服务间可能存在的关系[14]。然而国内外关于权衡与协调的研究存在一些局限。权衡与协调关系的相关分析是假设生态系统服务之间的关系是单调的,这种二元化的分析忽略了生态系统受多种因素共同影响的复杂性,样本点往往呈现散点云的分布,而不全是围绕着某一条线分布；其次,二元化的权衡与协同分析表征生态系统服务间的定量关系,却难以科学地解释背后的生态学机理。约束线分割法,被证明是分析生态系统服务非线性特征的适宜方法[5,15-17],但现阶段关于约束关系的量化及驱动机制的案例研究还有待于进一步深化。科尔沁沙地作为我国四大沙地(毛乌素沙地,浑善达克沙地,科尔沁沙地和呼伦贝尔沙地)之一,人地矛盾突出。近二十年来,政府致力于实施一系列生态工程调整土地利用类型,提高生态环境质量[18-19]。这些工程取得良好的效果,但2000—2015年间科尔沁沙地地区生态赤字仍在不断增加,大面积的农业灌溉对水资源造成压力,导致水资源可利用性减小[20-22]。李佳鸣等的研究指出,科尔沁地区生态恢复效率较低,其生态系统服务价值的提升滞后于内蒙古整体水平[23]。而有关科尔沁沙地的研究中,多以价值量评估生态系统服务,不仅较少利用模型评估生态系统服务物质的量,还缺乏服务间的非线性关系的讨论。

鉴于此,本研究基于遥感数据评估了2000—2018年科尔沁沙地长时间序列下的四种关键生态系统服务(净初级生产力、产水服务、土壤保持服务和风力侵蚀负服务)；揭示沙地关键生态系统服务的非线性约束关系,阐明约束作用关系的生态学意义,以期为实现科尔沁沙地的土地利用优化提供决策支持。

1 研究区概况

科尔沁沙地(117°49′—123°42′E,41°41′—46°05′N)位于内蒙古自治区东南部,地处东北平原向内蒙古高原的过渡地带,是北方农牧交错带的典型地区之一。总面积为12.29 × 104km2(沙化土地面积为1.79×104km2)[21],包括内蒙古自治区的15个县级行政县(图1)。沙地地势西高东低,海拔高度范围为86 m 至 1935 m,北部是大兴安岭及其丘陵地区,中部是冲积平原,南部是黄土丘陵区。土地利用主要以耕地、草地和林地为主。科尔沁沙地主要分布的河流有西辽河、老哈河和乌尔吉木伦河等；主要土壤类型有栗钙土、沙土、草甸土和盐碱土等；原生植被景观为典型草原到森林草原的过渡类型。该区域位于干旱与半干旱区,对气候变化敏感,年降雨量350 mm左右,主要集中在夏季(6、7、8月)；年平均风速3.5—4.5 m/s。在干旱多风等气候因素和过度的草地开垦等人类活动影响下,科尔沁沙地生态环境脆弱。除了国家层面的生态工程,在科尔沁沙地还特别地实施了2004年“双百万亩”综合治理工程和2014年“双千万亩”综合治理工程[24]。总的来说,沙地得到一定改善,但是仍有部分“远沙大沙”没有得到有效治理,沙漠化形势依旧严峻[25]。

图1 研究区地理位置Fig.1 Location of the study area

2 数据与方法

2.1 数据来源与处理

本文主要的数据说明与来源如表1所示。利用MRT(Modis Reprojection Tool)投影工具软件对归一化植被指数(normalized difference vegetation,NDVI)遥感数据进行镶嵌和定义投影等批处理,再通过最大值合成得到每月合成影像；选取研究区以及周边16个气象台站的降水、气温、风速和日照时数等地面气候资料,采用克里金插值法进行插值得到栅格数据；土地利用数据按照一级分类系统划分,包含林地、草地、水域、农田、建设用地和未利用地六大类。所有数据经重采样,空间分辨率为250 m,坐标与投影均保持一致(UTM Zone 50N,WGS 1984)。

表1 数据说明及来源Table 1 Description and resource of the study data

2.2 研究方法

2.2.1生态系统服务评估

本文定量估算了2000—2018年科尔沁沙地4种重要的生态系统服务功能,分别是净初级生产力(NPP),产水服务(WY),土壤保持服务(SC)和风力侵蚀负服务(WE)。

NPP是地表碳循环的重要组成,表征生态系统的质量状况。基于改进的CASA(Carnegie-Ames-Stanford-Approach)模型[26],NPP由植物吸收的光合有效辐射和实际光能利用率进行估算。CASA模型的输入数据包括月均温度、月降雨量、月总日照辐射、月NDVI、土地利用分类数据和相关生物参数[27]。

产水服务用产水量来表征,指降水量除去蒸腾作用部分的水源供给部分[28]。采用InVEST模型产水模块[29],基于budyko水量平衡原理方程得到产水量。产水模块输入数据包括月总降水量、月潜在蒸散发(PET)、土壤深度、土壤有效水分(PAWC)、土地利用/覆盖数据、流域矢量数据、生物参数以及Zhang系数。

土壤保持服务用土壤保持量来表征[30],即潜在土壤流失与实际土壤流失之间的差值。采用修订土壤通用流失方程(RULSE)[30],模型输入数据包括降水侵蚀力因子、土壤侵蚀性因子、地形因子、植被覆盖因子和水土保持措施因子。

风力侵蚀负服务采用修订的土壤风蚀方程(RWEQ)计算土壤风蚀量[31]。风力侵蚀是导致科尔沁沙地生态环境恶化的主要原因。模型输入包括气候因子、土壤可蚀性因子、土壤结皮因子、土壤糙度因子和结合残差因子。

2.2.2生态系统服务约束线提取方法

约束线能够刻画在多因素影响的复杂生态系统中限制变量对响应变量的限制作用,从而使响应变量不超过某一个范围或能够达到的最大值[32-33]。目前主要有4种绘制方法：参数法、散点云网格法、分位数回归法和分位数分割法[34]。本研究采取分位数分割法(图2),首先绘制成对生态系统服务的散点图,将x轴按照其值域平均分为100组,选取每一组中99.9%的分位数作为边界点,最后在Origin 2019中用提取到的边界点拟合得到约束线[16]。约束线上的点表示变量y受变量x影响时,变量y的理论最大值。

图2 约束线提取示意图Fig.2 Diagram of constraint line extraction

2.2.3约束线特征值分析方法

本文通过箱线图和离散度分析约束线特征值的变化情况。箱线图以四分位数和四分位距为基础,能够直观明了的识别异常值。在R语言软件中,使用plotbox函数确定超过上限和下限的约束线阈值。离散度即变异系数,用来分析线性约束线的特征值斜率和常数项。变异系数是标准偏差的绝对值除以平均值,变异系数越小其偏离程度越小。对多年线性约束线特征值变化趋势分析采用Sen 趋势分析[35]并结合Mann-Kendall 检验对变化趋势进行显著性分析[36-37],判断其是否在长时间序列上发生显著变化。

3 结果与分析

3.1 生态系统服务时空特征

2000—2018年科尔沁沙地4种生态系统服务均值空间分布及变化量如图3所示。植被净初级生产力多年平均值为256.5 gC/m2,空间分布具有明显的异质性。植被净初级生产力的高值区主要集中在研究区西北部的草地和林地,低值区主要分布在中部和西南部。从变化量来看,植被净初级生产力在东南部和中部增量高,最高可达25.3 gC/m2；在西部有所减少。产水服务的多年供给平均值为211.9 mm,空间格局差异明显。产水量高值区主要集中在科尔沁沙地东部,最高值可达522.2 mm；低值区分布在西北部。从变化量来看,中部偏东南地区及北部部分地区产水服务呈增加趋势；相反地,东部边缘地区和西部为产水量减少区域。土壤保持服务多年平均值达到196.84 t/km2,空间分布差异较小,高值区主要分布在坡度较高的西北和西南地区。从变化量来看,土壤保持服务整体呈增加趋势,增量最高值分布在西北大片区域和西南小片区域；局地地区土壤保持服务减少量稍小,约26.3 t/km2。风力侵蚀量的多年供给平均值为3.8 t/m2,空间分布明显差异。风力侵蚀高值区呈条带状分布在南部的沙地,土壤类型多为风沙土,涉及翁牛特旗、奈曼旗和科尔沁左翼后旗等旗县。从变化量来看,大部分地区风力侵蚀量有所减少,其中中部地区明显改善,风力侵蚀强度较大减少；南部沙地是风力侵蚀量增加的高值区,该地区的风力侵蚀情况仍在加剧。

图3 科尔沁沙地四种生态系统服务多年平均值及变化量空间分布Fig.3 Spatial distribution of the mean value and change slope of four ecosystem services in Horqin Sandy Land NPP:net primary productivity;WY:water yield;SC:soil conservation;WE:wind erosion

3.2 生态系统服务约束关系

不同水平下的生态系统服务两两之间的约束线提取结果如图4所示。植被净初级生产力与产水量、土壤保持量和风力侵蚀量,产水量与土壤保持量和风力侵蚀量,在景观水平和类水平上的约束线类型一致,两两之间呈现为抛物线类型。土壤保持量与风力侵蚀量,在景观水平和类水平上均呈现线性约束线。

图4 景观水平和类水平(森林、草地和耕地)生态系统服务约束关系Fig.4 Constraint line of the paired ESs at landscape level and class level (forest,grassland and farmland)

3.2.1植被净初级生产力和产水量的约束关系

在景观水平和类水平上,2000—2018年植被净初级生产力与产水量的约束线均呈抛物线类型(图4),约束线的形状在各个水平上几乎没有变化。随着植被净初级生产力的增加,其对产水量的约束作用先减少后增加,即当植被净初级生产力未超过阈值时,植被净初级生产力逐渐增加,其对产水量的约束作用逐渐减小；当植被净初级生产力超过阈值时,植被净初级生产力逐渐增加,植被生长消耗较多水分,在降水量总量有限的情况下,产水量减少,因此植被净初级生产力对产水量的约束作用逐渐增强。

3.2.2植被净初级生产力和土壤保持量的约束关系

植被净初级生产力与土壤保持量的约束线在景观水平和类水平上均呈抛物线类型(图4),但大多数年份的约束线峰值不够明显,且约束线上的值主要处于阈值左侧,植被净初级生产力对土壤保持量的约束作用随着植被净初级生产力增加而逐渐减小。在一定程度上,植被净初级生产力增加,能够提升水土保持能力,有效减少水分对土壤的侵蚀作用。但是植被净初级生产力超过阈值时,降水一方面促进植被生长,维持较高的植被净初级生产力,但另一方面较多降水也会增加土壤侵蚀程度,使得土壤保持能力降低。植被净初级生产力与土壤保持量约束线形状在耕地和林地、草地两个水平上有较明显差异；耕地水平的约束线阈值低于林地和草地水平的阈值,表明在相同植被净初级生产力的情况下,耕地提供的土壤保持能力较低。

3.2.3植被净初级生产力和风力侵蚀量的约束关系

植被净初级生产力与风力侵蚀量的约束线与产水量相似,在景观水平和类水平上均呈开口向下的抛物线,驼峰阈值显著(图4)。一般而言,高植被生产力代表较好的植被覆盖状况,一定程度上能够减少风速,同时发达的植被根系能够增强土壤抵抗风力侵蚀的能力[38]。但是在植被净初级生产力阈值的左侧,植被净初级生产力较低,土壤植被覆盖状况较差,因此受风场强度等气候因子驱动,尽管植被净初级生产力增加,风力侵蚀量也仍然增加。

3.2.4产水量和土壤保持量的约束关系

各个水平上产水量和土壤保持量的约束线均呈抛物线类型,但是形状差异明显,其中景观水平上和草地水平上的约束线相似度高,在耕地水平上个别年份约束线形状变化剧烈,这与耕地的耕作方式和种植类型有关(图4)。大部分年份,随着产水量的增加,土壤保持量呈现下降趋势,产水量对土壤保持量的约束作用逐渐减弱。降水是影响产水量和土壤保持量的共同关键因素,当降水达到产生径流量时,会降低土壤保持能力,尤其在砂质土壤分布区域。在同等产水量的条件下,林地和草地的土壤保持量较耕地的土壤保持量大。

3.2.5产水量和风力侵蚀量的约束关系

产水量和风力侵蚀量在各个水平上的约束线均为抛物线类型,大多数年份的驼峰扁平,少部分年份的驼峰明显(图4)。除了耕地,其他类水平的产水量和风力侵蚀量之间约束作用较弱,随着产水量的增加,风力侵蚀量变化不明显。一般情况下,在土地退化严重的沙地,风力侵蚀作用较强。风驱动土壤发生侵蚀,同时风促进地表径流蒸发,且砂质土壤有较强的渗透作用,使得产水量减少。研究区降水量少,常年干旱缺水,产水量的波动较小。风力侵蚀是该区域存在的突出问题,产水量的变化不足以影响风力侵蚀的趋势性变化。耕地水平上产水量与风力侵蚀量约束关系的变化与当地部分农田采用灌溉补水、追求高作物产量的耕作方式有关[39]。在该区域应当采取保护性耕作措施,增加土壤表面粗糙度,减少风力侵蚀量；另一方面,使用节水灌溉、垄作耕作方式可以提高水资源的利用率,减少风力侵蚀量[20]。

3.2.6土壤保持量和风力侵蚀量的约束关系

在各个水平上,土壤保持量与风力侵蚀量的约束线均呈现负向线性关系,随着土壤保持量的增加,风力侵蚀量减少 (图4)。科尔沁沙地土壤主要为砂质土壤,流动性强,保水性差,很容易被水力和风力侵蚀破坏。在降水少的地区,植被覆盖度低,土壤保持服务很低,风力侵蚀很容易发生,因此土壤的水蚀和风蚀这两者很难同时发生。在降水多的地方,土壤湿度大,地表植被覆盖度高且土壤保持量高,会遏制风蚀的发生,因此在各个水平上这两种服务之间呈现负向的线性约束关系。

3.3 约束线特征值

3.3.1阈值特征

在景观水平上,呈抛物线类型的非线性关系的约束线阈值统计如表2所示。由表2可知,抛物线类型的约束线阈值表现出一定的时间差异。

表2 景观水平上抛物线型约束线阈值Table 2 The thresholds on the hump-shaped constraint lines at landscape level

由于影响产水量的土壤类型、土壤深度和地形等基本保持不变,植被净初级生产力与产水量阈值产生时间差异的主要原因与降水量、潜在蒸散量的变化有关。植被净初级生产力和土壤保持量受到多种因素的影响,在计算土壤保持量的RUSLE模型中,土壤可蚀性因子和地形因子在一定时间内保持稳定,降水量、土地利用变化和植被覆盖度的变化均可能引起净初级生产力和产水量阈值在时间上的差异。植被净初级生产力和风力侵蚀量约束关系的影响因素主要为土壤类型、风速和植被覆盖度。

在植被净初级生产力有关的3对生态系统服务约束线阈值的箱线图中(图5),每对生态系统服务约束关系中的植被净初级生产力的多年阈值变化范围较为集中,相对稳定,尤其在与产水量、风力侵蚀量的约束关系上表现更为明显。而植被净初级生产力阈值对应的产水量、土壤保持量和风力侵蚀量阈值变化幅度大,其中土壤保持量阈值的偏离程度最大。本研究发现,植被净初级生产力作为重要的支持服务,对其他生态系统服务起着非线性放大或缩小的重要调节作用。植被净初级生产力微小的变化会引起产水量、土壤保持量、风力侵蚀量较大幅度的变化。

图5 景观水平的生态系统服务约束线(NPP-WY,NPP-SC,and NPP-WE)阈值的箱线图Fig.5 The box plots of thresholds on the constraint lines of NPP-WY,NPP-SC and NPP-WE at landscape levelNPP1、NPP2和NPP3分别对应在NPP-WY,NPP-SC,NPP-WE中NPP的阈值;NPP:(gC/m2),WY:(mm),SC:(t/hm2),WE:(t/m2)

3.2.2线性约束关系

土壤保持量在各个水平上对风力侵蚀量有负向线性约束关系,其特征值斜率(k)和截距(b)如表3所示,其中斜率代表着约束关系的强弱；截距表示不存在限制因子时,响应因子能够达到的最大值。2000—2018年,各个水平上k值具有显著的增加趋势,土壤保持量和风力侵蚀量之间的约束关系逐渐减弱；除了耕地水平上的截距较小,其他水平上的截距平均值几乎一致。结合耕地在其他服务约束关系中提供服务的能力,研究区的耕地应当注重产水量与土壤保持量之间的关系,选择需水量较少、耐旱的作物能够有效改善耕地的整体生态系统服务。

表3 生态系统服务线性约束线的斜率和截距Table 3 Slopes (k) and constant terms (b) of the linear constraint lines between paired ecosystem services

4 讨论

目前关于沙地生态系统服务的研究结果表明生态恢复政策对环境质量改善有积极作用[40],但具体的政策实施方案仍需不断优化。本文通过探究沙地生态系统服务的非线性关系,即各项服务间的约束线及其特征值,为沙地土地利用优化提供了定量依据,也为预测服务的潜在最大值提供了有效手段[41]。植被净初级生产力与产水量的抛物线类型约束线一方面很好地解释了科尔沁沙地在荒漠化逆转过程中出现水资源短缺的问题,当植被净初级生产力增加到超过阈值时,产水量反而会减小,约束作用变强。另一方面,其阈值对于优化区域整体生态系统服务的供给,合理分配资源至关重要。Hao等[42]在锡林郭勒盟的生态系统服务约束关系研究中同样指出,植被净初级生产力的阈值是实现植被净初级生产力、产水、土壤水蚀控制和土壤风蚀控制协同的重要参考指标。官冬杰等[43]以重庆市湿地为研究区的结果表明,在景观水平上水资源供给服务和水土保持服务间呈抛物线型约束关系,降水和气温作为主要的自然因素直接影响湿地的生态系统服务。相比于水资源高供给的湿地,科尔沁沙地产水量低导致本研究中的产水量与土壤保持量的抛物线型约束关系较弱。戴路炜等[11]关于北方农牧交错带多伦县的研究结果表明,土壤保持与防风固沙服务之间为权衡关系,但无法解释生态系统服务间相互作用的阈值特征；本研究很好地揭示了土壤保持量和风力侵蚀量线性约束关系的斜率与截距,可为政策制定者提供定量依据。

由于数据可得性限制,本文重点探讨了支持服务(净初级生产力)和调节服务(土壤保持、产水和风力侵蚀),但对供给服务和文化服务没有涉及,政策制定者应充分注重景观的多功能性,以整体生态系统服务最优为目标,而不是某项生态系统服务达到最大化。

5 结论

通过对科尔沁沙地生态系统服务时空变化及约束关系的研究,本文得到以下结论：

(1)2000—2018年科尔沁沙地植被净初级生产力、产水量、土壤保持量和风力侵蚀量的多年均值分别为256.5 gC/m2,211.9 mm,26.3 t/km2和3.8 t/m2。翁牛特旗、奈曼旗和科尔沁左翼后旗等旗县是风力侵蚀的重度区,应当重点治理。

(2)植被净初级生产力与产水量、土壤保持量、风力侵蚀量的约束线均呈开口向下抛物线类型,随着植被净初级生产力的增加,三种服务均呈先增加后减小的变化趋势,表明植被净初级生产力对三种服务的约束效应均呈先减小后增加的变化特征。产水量与土壤保持量、风力侵蚀量的约束线也呈抛物线类型,但抛物线形状扁平趋于直线。土壤保持量与风力侵蚀量的约束关系呈线性类型。

(3)植被净初级生产力的变化对其他生态系统服务的动态起到放大器的作用,它的微小变化可能会引起产水量、土壤保持服务、风力侵蚀量较大幅度的变化。以植被净初级生产力与其他类型生态系统服务的约束线阈值为基础,合理利用或提高区域植被净初级生产力,可以减少其动态变化对其他服务的负面约束作用,从而促进整体生态系统服务达到较优状态。

(4)林地和草地提供土壤保持服务的能力较耕地高,应重点发展草地、林地生态系统。同时,应在耕作条件适宜的地区选择需水量较少的农作物并实施节水灌溉,有效改善和提高耕地的生态系统服务,以提高科尔沁沙地景观总体生态系统服务的供给能力。