薛曾辉,高驭洋,卢枰达,周静,王可琳,张斌,马玉,刘梦云
(1.西北农林科技大学资源环境学院,陕西杨凌 712100;2.农业农村部西北旱地农业绿色低碳重点实验室,陕西杨凌 712100;3.农业农村部合阳农业环境与耕地保育科学观测实验站,陕西合阳 715300)
随着自然资源的开采,人类对生态系统产生了极大的影响,生态系统退化等问题日益突出[1]。生态系统服务(Ecosystem Services,ES)是指人类通过生态系统的结构、过程和功能直接或间接得到的生命支持产品和服务[2],将人类福祉与自然环境联系在一起,支持人类生存和经济发展[3]。各类ES随着地形和土地利用的变化出现明显空间分异,ES的空间异质性现象是生态评估的热点方向。不同ES间产生不同相互作用与联系[4],一个ES增强而其他ES随之增强、减弱或无变化,即产生协同、权衡或独立作用[5]。ES间通常是复杂且相关的[6],国内外学者使用遥感和地理信息科学相关手段,在不同区域不同尺度下取得了众多成果[5,7-11]。Hu等[5]通过斯皮尔曼(Spearman)秩相关性发现北欧国家ES间相关关系因空间而异;陈相标等[7]通过空间自相关等方法发现石林喀斯特岩溶区ES 间权衡与协同关系空间异质性显著;张静静等[8]采用空间叠置分析法确定了伏牛山地区森林生态系统的权衡协同区域;刘华超[9]、汪仕美[10]和武彦珍[11]等分别探究了大兴安岭、子午岭区和慈利县的ES 间相关关系,不过ES间权衡-协同-独立关系随地形和土地利用的变化机制有待进一步探究。
安塞区生态环境面临威胁,探究土地利用和地形下ES变化规律是精准开展生态治理的关键,亟需探究县域尺度下土地利用和地形对ES的影响。采用In VEST 和CASA 模型,定量评估安塞区2020年土壤保持(Soil Conservation,SC)、固碳释氧(Carbon Fixation and Oxygen Release,CO)、水源涵养(Water Harvesting,WH)和生境质量(Habitat Quality,HQ)4类ES,分析地形和土地利用下的各ES变化趋势,通过斯皮尔曼(Spearman)秩相关性系数计算ES间权衡-协同-独立关系,根据研究结果对安塞区土地综合利用规划及生态保护政策制定提供科学依据,并为县区尺度下生态文明建设提供参考。
安塞区位于陕西省延安市北部,东经108°5′44′′—109°26′18′′,北纬36°30′45′′—37°19′3′′,西毗志丹县,北靠靖边县,东接子长市,南与甘泉县、宝塔区相连;属于典型的中国西北黄土高原河谷丘陵区和河谷平原沟壑区,地形地貌丰富,面积2 950 km2,平均海拔1 372 m,境内有延河、大理河和清涧河3条水系;气候特征为中温带大陆性半干旱季风气候,四季长短不等,干湿分明,年平均气温8.8℃,年平均降水量505.3 mm,年日照时数2 390 h;土地利用类型包括草地、林地、耕地和未利用地等,草地面积占比超过50%,是安塞区最主要的土地利用方式之一,对生态环境质量有着极其重要的影响。研究安塞区土地利用和地形下各ES 及其相关关系,对其他受土壤侵蚀和水分流失影响的县区尺度ES间权衡-协同-独立研究具有典型示范性。安塞区地形图如图1所示。
图1 安塞区地形图Fig.1 Topographic map of Ansai District
数据主要包括:(1)土地利用数据,来源于欧洲航天局(ESA)World Cover 10 m 空间分辨率2020年土地覆盖数据(https:∥esa-worldcover.org/en),根据一级分类划分为耕地、林地、草地、建设用地、水域和未利用地。(2)归一化植被指数(NDVI)数据,来源于欧洲航天局(ESA)哨兵2 号(Sentinel-2)影像(https:∥scihub.copernicus.eu/dhus/#/home),在Google earth engine 中下载和处理得到2019 年5月、6月、10月、2020年2月、3月、4月、7月、8月、11月、12月、2021年1月、9月的10 m 空间分辨率NDVI数据。(3)土壤数据,来源于世界土壤数据库(HWSD)(http:∥data.tpdc.ac.cn/zh-hans/),数据空间分辨率为1 km,提取至点后进行插值。(4)气象数据,来源于国家气候中心(http:∥data.cma.cn/)提供的中国地面气候资料日值数据集(V3.0),选取研究区周边7个气象站点插值获得栅格数据。(5)数字高程模型(DEM),来源于美国国家航空航天局(NASA)网站的ALOS 12.5 m DEM 数据(https:∥search.asf.alaska.edu/#/)。(6)行政边界数据,来源于国家基础信息中心(http:∥www.ngcc.cn/ngcc/)。以上数据使用Arc Map软件统一坐标投影,栅格数据重采样至10 m 空间分辨率。
1.3.1 生态系统服务计算 SC,WH 和HQ 采用In-VEST 模型泥沙输移比(Sediment Delivery Ratio)、产水量(Annuai Water Yield)、生境质量(Habitat Quality)3个模块进行估算,具体计算公式参照刘华超[9]和吴志俊[12]等相关研究,降雨侵蚀力参照Wischmeier根据逐月平均降雨量提出的经验公式进行计算[13],土壤可蚀性参照EPIC(Erosion Productivity Impact Calculator)模型计算[14],植被可利用含水率参照周文佐总结的经验公式[15],In VEST 模型运行参数参考相关文献[16-20]和模型手册,并根据研究区实际情况进行取值。
植被净初级生产力(NPP)是植物单位时间单位面积上光合作用产生的净有机物量,表征植被对大气中二氧化碳的吸收能力[21],在ES 测算中用CASA模型估算NPP以表征CO,具体公式参照薛晓玉[22]和Jiang[23]等相关研究中NPP的计算。
1.3.2 生态系统服务权衡与协同关系分析 以各类ES评估结果为基础,使用SPSS进行Spearman秩相关系数分析,计算公式如下:
式中:a和b代表不同类别ES;R(a,b)代表a和b的相关系数;dab表示a和b的等级差值;n为样本个数;t为检验值。当R(a,b)>0时,代表a和b具有正相关,即ES间的协同,且数值越大代表协同越强;当R(a,b)<0时,代表a和b具有负相关,即ES间的权衡,且数值越小代表权衡越强。t检验方法可判断a和b相关的显著性p。
根据相关系数R及其显著性p进行权衡协同等级划分(表1),分析ES间的权衡-协同-独立关系,正值越大代表协同关系越强,负值越小代表权衡关系越强,0值则代表独立关系。
表1 权衡协同关系赋值表Table 1 Tradeoff and synergy assignment table
2020年安塞区ES空间分布如图2所示。SC分布整体均匀,低海拔区域数值较低,平均值根据土地利用排序为林地>草地>水域>建设用地>耕地>未利用地;CO 分布为南高北低,低海拔区域数值较低,高值集中在南部林地区域,平均值排序为林地>草地>耕地>未利用地>建设用地>水域;WH 分布为南高北低,水域明显低于其他区域,高值集中在低海拔区域,平均值排序为林地>建设用地>耕地>草地>未利用地>水域;HQ 分布中建设用地和耕地周围明显低于其他区域,平均值排序为林地>草地>水域>未利用地>耕地>建设用地。安塞区ES受地形和土地利用影响较大,呈现出明显的空间异质性。
图2 安塞区生态系统服务功能分布Fig.2 Distribution of ecosystem services in Ansai District
分析安塞区各ES随海拔、坡度和坡向等地形变化趋势(图3),结合自然间断法将海拔分级为低海拔区958~1 172 m、中低海拔区1 172~1 274 m、中高海拔区1 274~1 369 m 和高海拔区1 369~1 637 m,将坡度分级为平缓坡0°~12°、倾斜坡12°~24°、陡峭坡24°~36°和急悬坡36°~72°,将坡向分级为半阴坡50°~140°、阳坡140°~230°、半阳坡230°~320°和阴坡320°~50°。同时分析在安塞区不同土地利用下ES随海拔、坡度和坡向等地形变化趋势(图4)。
图3 安塞区生态系统服务随地形变化趋势Fig.3 The change trend of ecosystem services in Ansai District with topography
图4 安塞区不同土地利用生态系统服务随地形变化趋势Fig.4 The change trend of ecosystem services for different land use in Ansai District with topography
分析图3和图4可知,随着海拔升高,SC先升高后降低,在水域呈上升趋势,未利用地呈下降趋势;CO 先升高后降低再升高,在中低海拔区达到最高值,林地、建设用地和未利用地先上升再下降,水域先上升后下降;WH 呈下降趋势,耕地在高海拔区域略有上升,林地呈波浪形下降趋势,水域变化幅度较大;HQ 呈上升趋势,在林地和草地的数值明显高于其他地类,耕地和未利用地低值平缓波动,水域先下降后上升。随着坡度增加,SC先升高后降低,在急悬坡上升趋势放缓,在未利用地呈下降趋势,水域呈递上升趋势;CO 先上升后降低,在未利用地呈波浪上升趋势,林地呈高位下降趋势,耕地先降低后升高趋势,建设用地和水域呈上升趋势;WH 平缓上升后下降,在水域和建设用地先下降后上升,林、耕地和未利用地呈上升趋势;HQ 上升后平缓下降,在各土地利用上变化平缓但数值差异明显。随着坡向变化,SC在半阴坡达到最大值,半阳坡达到最小值,在建设用地呈相反趋势;CO 在阳坡达到最小值,阴坡达到最大值,在林地呈相反趋势,其他地类低水平波动;WH 在阳坡达到最大值,阴坡达到最小值,在耕地呈相反趋势;HQ 随坡度变化同CO基本一致,变化幅度较小。
各土地利用上ES空间分布规律差异明显,林、草和耕地ES值较高,且与整体一致性较高,对安塞区ES有着较强贡献,占据主导地位;建设用地、未利用地和水域受到自然条件约束,在高海拔和急悬坡分布较少,且各类ES值较低,与其他地类表现较大的差异。ES在低海拔和高海拔、平缓坡和急悬坡、阴坡和阳坡均有明显波动,以上区域存在较多的人类活动和自然条件限制,进而影响了ES的空间分布,通过政策法规和生态治理,受影响区域有较大的ES提升潜力。ES随海拔和坡度的变化幅度明显大于坡向,SC对于地形变化的敏感程度最强,HQ 对地形变化的敏感程度最弱。
通过表1赋值,计算各分类下SC-CO等6种关系的权衡-协同-独立值,正值相加为ES间协同度,负值相加取绝对值为ES间权衡度,为便于数据对比,将0值个数乘以4为ES间独立度,得到不同土地利用(图5)和地形(图6)ES间权衡-协同-独立关系和地形作用下ES间权衡-协同-独立度(表2)。
表2 安塞区地形作用下各生态系统服务间权衡-协同-独立度Table 2 Tradeoff-synergy-independence of ecosystem services under the influence of topography in Ansai District
图5 安塞区不同土地利用生态系统服务权衡-协同-独立关系雷达图Fig.5 Tradeoff-synergy-independence relationship of ecosystem services for different land use in Ansai District
图6 安塞区不同地形生态系统服务权衡-协同-独立关系雷达图Fig.6 Tradeoff-synergy-independence relationship of ecosystem services at different terrain in Ansai District
2.3.1 土地利用作用下的生态系统服务权衡-协同-独立分析 从土地利用来看(图5),协同度排序为不区分地类(258)>耕地(214)>草地(175)>林地(169)>水域(60)>未利用地(59)>建设用地(52);权衡度排序为未利用地(135)>不区分地类(91)>草地(67)>林地(24)>耕地(12)>建设用地(10)>水域(9);独立度排序为水域(236)>林地(144)>耕地(128)>草地(124)>未利用地(120)>建设用地(92)>不区分地类(48);累计值排序为耕地(202)>不区分地类(167)>林地(145)>草地(108)>水域(51)>建设用地(42)>未利用地(-76)。耕、草地和林地协同度明显高于其他地类,未利用地权衡度最高,水域独立度最高,不区分地类分析的相关关系结果与各土地利用有明显差异;耕地协同度和相关关系累计值最高,紧接着是林地和草地,表明了《耕地保护法》和退耕还林还草政策的实施成效;水域的权衡度和协同度均较低,但独立度最高,展示了水域生态的相对独立性;草地的协同和权衡能力均较强,展现了相关关系的复杂性,同时展现了草地对于安塞区ES的重要意义。
2.3.2 海拔作用下的生态系统服务权衡—协同—独立分析 在土地利用和海拔视角下(表2海拔部分),安塞区ES间以协同为主。协同度最高值为CO-WH(90),权衡度最高值为WH-HQ(60),独立度最高值为SC-WH(48),累计值排序为CO-WH(68)>SC-CO(54)>SCWH(52)>CO-HQ(10)>WH-HQ(-7)>SC-HQ(-27)。CO,WH 和SC 两两有较强的协同,HQ 与其他ES的相关关系以权衡为主。
从海拔区间来看(图6),协同度排序为海拔总区间(79)>中高海拔区(79)>高海拔区(78)>低海拔区(62)>中低海拔区(52);权衡度排序为海拔总区间(53)>低海拔区(50)>中低海拔区(42)>中高海拔区(41)>高海拔区(16);独立度排序为中低海拔区(64)>中高海拔区(44)>低海拔区(32)>高海拔区(28)>海拔总区间(24);累计值排序为高海拔区(62)>中高海拔区(38)>海拔总区间(27)>低海拔区(12)>中低海拔区(10)。海拔总区间与海拔分级区间下的ES 相关关系表现出差异性,且随着海拔升高,ES间权衡作用减弱,协同作用加强,相关生态治理应侧重于人类活动强烈的低海拔区域。
2.3.3 坡度作用下的生态系统服务权衡-协同-独立分析 在土地利用和坡度视角下(表2坡度部分),安塞区ES 间以协同为主。协同度最高值为CO-WH(102),权衡度最高值为WH-HQ(29),独立度最高值为SC-WH(72),累计值排序为CO-WH(93)>SC-HQ(57)>SC-WH(29)>CO-HQ(27)>SC-CO(23)>WH-HQ(14)。CO和WH 在ES间相关关系上起主导作用,HQ与其他ES间相关关系较为薄弱。
从坡度区间来看(图6),协同度排序为倾斜坡(86)>陡峭坡(79)>坡度总区间(73)>平缓坡(53)>急悬坡(43);权衡度排序为平缓坡(36)>急悬坡(17)=坡度总区间(17)>倾斜坡(11)>陡峭坡(10);独立度排序为急悬坡(88)>坡度总区间(60)=倾斜坡(60)=平缓坡(60)>陡峭坡(56);累计值排序为倾斜坡(75)>陡峭坡(69)>坡度总区间(56)>急悬坡(26)>平缓坡(17)。坡度总区间与坡度分级区间下的ES 相关关系表现出差异性,且随着坡度升高,ES间协同作用呈现先上升后下降变化趋势,权衡作用反之,相关生态治理应侧重于人类活动强烈的平缓坡和自然环境恶劣的急悬坡区域。
2.3.4 坡向作用下的生态系统服务间权衡-协同-独立分析 在土地利用和坡向视角下(表2 坡向部分),安塞区ES 间以独立为主。协同度最高值为CO-WH(85),权衡度最高值为WH-HQ(25),独立度最高值为SC-CO(96),累计值排序为CO-WH(80)>SC-HQ(56)>SC-WH(36)>CO-HQ(50)>WH-HQ(16)>SC-CO(14)。CO 和WH 主导ES间相关关系。
从坡向区间来看(图6),协同度排序为阳坡(72)>半阳坡(65)>坡向总区间(60)>半阴坡(53)>阴坡(52);权衡度排序为半阴坡(16)>阴坡(13)=坡向总区间(13)>半阳坡(8)>阳坡(5);独立度排序为阴坡(84)>半阴坡(76)=半阳坡(76)>坡向总区间(72)>阳坡(76);累计值排序为阳坡(67)>半阳坡(57)>坡向总区间(47)>阴坡(39)>半阴坡(37),坡向总区间与坡向分级区间下的ES相关关系表现出了差异性,坡向的变换对ES间相关关系产生了极大的影响,随着坡向向阳程度的加强,ES间协同能力随之提高,权衡能力随之降低,相关生态治理应侧重于阴坡区域。
土地利用和地形是影响ES分布、结构及功能的重要因素[24]。WH 随地形因子和土地利用变化情况与其他ES有着较明显的差异,参考胡砚霞[25]和高江波等[26]的相关研究可知,降水量为WH 最优解释力,通过降水量在地形和土地利用的分布间接影响了WH,人类活动和政策制定等对提升WH 值只能起到辅助作用;SC,CO 和HQ 在林地和草地较高,在耕地和未利用地较低,这与退耕还林还草的工作理念一致[27],表明了相关政策的合理性;SC,CO 和HQ 随海拔和坡度增加呈先增后减趋势,这与徐彩仙[24]和郜红娟等[28]研究结论一致,低海拔和平缓坡有强烈的人类活动,高海拔和急悬坡未利用地较多,导致ES平均值较低,应加快未利用地向其他土地利用方式转化;坡向变化影响热量差异[29],CO 和HQ 最高值在阴坡处,最低值在阳坡处,与陈奕竹等[30]研究结论一致,SC最高值在半阴坡,阳坡和半阳坡各类ES值较低,说明土壤水分和温度等对SC 驱动作用超过了太阳辐射。
国内外学者ES间权衡—协同—独立研究集中于各类ES间的相互作用、自然人文系统的分析以及尺度异质性上,对于在地形和土地利用变化下的研究分析较少[31]。相关学者在全流域和二级流域[32]和县域尺度与区域尺度[33]ES相关关系等均发现空间异质性在ES间相关性分析中表现明显,段宝玲等[34]在ES间相关关系研究中涉及地形部分也有类似结果。同Zhang[35]和Wu[36]等结论一致,ES间的协同能力远大于权衡,各类ES整体上呈协同作用,WH-HQ 在各地形因素下权衡程度均强于其他ES间权衡度。ES间协同作用随海拔升高而增强,随坡度增加呈先升高后降低变化趋势,与孙艺杰等[31]在延安市ES间相关关系研究结论基本一致,不同之处在于急悬坡处,原因可能在于ES选取和研究尺度不同。ES间协同作用随坡向向阳程度加强而增强,与张静静等[8]研究得到的南坡优于北坡结论一致,太阳辐射对人类活动与自然条件均产生了较大的影响,对于ES 间协同作用有促进效应。地形和土地利用变化下,气象条件、人类活动和植被类型等均发生了相应变化,对ES的空间分布和ES间相关关系产生了较大影响。
经过1999年和2014年两轮退耕还林还草政策的实施,安塞区ES持续向好发展,研究发现安塞区2020年生态系统服务整体表现为协同度>独立度>权衡度,各ES间随地形因素变化下协同—独立—权衡关系并存,协同作用常发生在SC,CO 和WH 之间,HQ 与其他ES易发生权衡作用。退耕还林还草方案要求,将陡坡耕地和梯田、重要水源地15°~25°坡耕地、严重沙化耕地进行退耕还林还草,实现了土地利用由耕地向草地和林地的转化,在改善植被的同时减少了土壤流失。林地和草地ES值高于耕地,随坡度等地形因子的升高,耕地SC 值降低,故在地形与土地利用角度下,退耕还林还草政策实现了SC 与CO 协同发展,王芸等[27]在陕北地区ES研究亦有类似结论。林地和草地的蒸散发量高于耕地,进而影响WH,WH 另一主要驱动因素为降水量[25-26],退耕还林还草背景下为提升WH 值及WH 与其他ES间的协同作用,应结合地形梯度变化,从减少蒸散发量和增强降水量等方面开展WH 生态治理工作。HQ 评估中生境威胁因子为建设用地、未利用地和耕地,随着经济社会发展,建设用地呈扩张趋势;随着退耕还林还草政策的继续实行,非平缓坡的未利用地和耕地呈缩减趋势,安塞区2020年HQ 与其他ES权衡作用明显,说明了建设用地扩张趋势对HQ 驱动作用大于退耕还林还草的生境影响,尤其在地形平缓区域。后续应平衡城市扩张与生态治理的平衡关系,在坚守耕地红线的基础上,继续大力开展退耕还林还草政策,力争SC,CO,WH 和HQ 等ES均朝协同方向发展,在地形和土地利用视角下实现县域尺度ES良性向好发展。
分析县域尺度下ES空间分布有助于提供科学依据支持政策制定,管理和优化人类活动。为提升安塞区ES水平及ES间协同能力,应根据区域生态提升潜力程度,因地制宜地实施不同强度的生态治理,加快未利用地向林地和草地等ES水平较高的土地利用转化,重视人类活动、生态治理和资源环境的和谐发展,在低海拔、高海拔、低坡度、阴坡、未利用地、建设用地和耕地等ES较弱或ES间权衡较强的区域开展相关生态治理,向提升ES值及减弱ES间权衡的方向进行政策制定和方案实施。ES低于区间中值和ES间为权衡作用的区域被定义为生态提升潜力区域,如表3所示。海拔视角下的生态提升潜力区域远多于坡度和坡向视角,安塞区地形起伏较大,生态治理中应着重注意地形对ES的相关影响。
表3 安塞区土地利用和地形视角下的生态提升潜力区域Table 3 Ecological improvement potential area from the perspective of land use and topography in Ansai District
安塞区ES研究亦存在些许不足,研究区附近气象站点数量较少,无法使用Anusplin等专业气象插值软件,克里金插值结果相对来说精度有限,影响ES结果精度;同时因分析侧重点等原因,本研究未对SC-CO,SC-WH,CO-WH,SC-HQ,CO-HQ,WH-HQ 等ES间权衡—协同—独立关系分别随地形和土地利用变化情况进行深入研究;ES间权衡—协同—独立关系中土地利用和地形因子的驱动机制研究不足。下一步可从SC-CO 等各类ES间相关关系随地形和土地利用变化、驱动力因素等进行深入探讨,为促进县域尺度下生态环境协调发展提出更加精准的土地管理建议及生态环境方案。
(1)安塞区林地和草地生态系统服务及其协同能力较强,未利用地和建设用地则较弱。地形和土地利用变化下,对生态系统服务空间分布及其相关关系产生了较大影响。
(2)安塞区生态系统服务以协同为主,土壤保持、固碳释氧和涵养水源生态系统服务均以协同为主,生境质量与其他生态系统服务权衡作用明显,退耕还林还草政策一定意义上促进了生态系统服务的协同作用。
(3)安塞区土地利用和地形因子总区间与分级区间下生态系统服务有较大差异,且随着海拔升高和坡向向阳程度的加强,生态系统服务协同作用加强,权衡作用减弱;随着坡度增加,协同作用呈先上升后下降趋势,权衡作用反之。生态系统服务权衡—协同—独立关系对土地利用与地形变化响应敏感,应根据生态提升潜力因地制宜开展生态治理。