陕北长城沿线植被恢复与生态系统防风固沙服务模拟分析

米朝娟,周自翔,刘 婷,武 佳,孙彦旭

西安科技大学测绘科学与技术学院,西安 710054

生态系统防风固沙服务是干旱和半干旱地区生态系统最重要的生态服务功能之一,实施防风固沙工程是抑制和固定风沙区沙尘不被再次移动以实现荒漠化逆转,有利于区域经济和人类福祉的可持续发展[1]。生态系统防风固沙服务通常从风蚀角度为防风固沙服务提供研究。学术界对风蚀的评估始于1940年代,风蚀模型是评估土壤风蚀状况的主要技术手段,包括德克萨斯模型(Texas erosion analysis model,TEAM)[2—3]、风蚀预报系统(WEPS)[4—5]、风蚀方程(WEQ)[6]、修正风蚀方程模型(Revised Wind Erosion Equation,RWEQ)[7—9]等。由于一些模型主要集中在风蚀发生的可能性,无法评估潜在风蚀量。而大多数学者考虑到REWQ模型具有良好的适用性,涉及气候条件、植被覆盖、土壤易蚀性、土壤结皮、地表粗糙度等综合因素来计算植被造成的沙粒滞留量所发生的风力侵蚀[10—12]。部分学者不断探索和验证该模型并通过调整参数,利用RWEQ模型对宁夏[7—12]、内蒙古[13—14]、京津冀[15—17]等中国干旱半干旱风沙区的风蚀状况进行了评估,取得了良好的效果。

植被是影响防风固沙生态功能的关键指标,也是检验防风固沙区生态保护成效的重要依据[11]。恢复植被有利于保持土壤水分、改变地表状态,从而降低风蚀发生的潜力[12]。风蚀的动态变化是自然气候变化和人类活动共同作用的结果,退耕还林使地表覆盖发生显著变化,进而与生态系统服务能力的变化相关度最高[18]。应用地理探测器模型分析自然和人为因子对宁夏地区防风固沙服务功能空间格局形成及演变过程中的贡献与交互作用[1]发现土壤类型和植被类型在防风固沙服务功能空间分布格局中发挥最重要的作用。而以归一化植被指数(NDVI)计算植被覆盖度,结合气象要素从不同尺度来评估宁夏灵武白芨滩自然保护区防风固沙功能时[12],植被覆盖度的增加与研究区的防风固沙服务成正比。朱趁趁等[14]利用RWEQ模型实现土地利用变化对内蒙古荒漠草原的防风固沙能力进行评估。结果表明：土地利用方式以林地恢复、建设用地扩张、不同覆盖度草地转化变化为主对防风固沙服务有一定的增强作用。徐洁等[19]基于RWEQ和HYSPLIT模型模拟了生态功能区防风固沙服务和防风固沙服务空间流动路径,从生态系统服务流动的角度建立了生态功能区及其防风固沙服务受益区之间的时空联系。而以往的研究大多是大区域尺度,很少评估土壤风蚀对风沙区边缘地带的影响,而且采用的气象数据精度不够。由于目前风沙边缘地区的气象环境变化、植被恢复与防风固沙能力之间的相互作用机制与如何保证防风固沙服务功能的可持续性还需进一步研究。

基于生态工程的防风固沙、植被恢复、土壤保育等生态功能监测与评价,可为促进地区经济和社会的可持续发展提供实践指导[20]。在我国土地沙化依然在扩展的趋势下,陕北长城沿线处于毛乌素沙漠的边缘地区,长期受到自然因素和人类活动扰动影响,是最脆弱的生态系统之一[21]。目前,榆林大地基本上消灭了沙漠化土地现象,在全球沙漠治理的意义上率先实现了荒漠化的逆转。但该地在干旱少雨、生态脆弱的自然条件下,导致频发干旱和风蚀现象。因此,从生态系统服务的角度探讨植被恢复、气象环境变化与防风固沙服务能力之间相互作用机质以及如何保证防风固沙服务的可持续性。本文研究内容包括：(1)分析陕北长城沿线2000—2018年植被覆盖时空变化。(2)利用RWEQ模型模拟长城沿线生态系统防风固沙服务物理量的时空模式,实现防风固沙服务的量化。(3)利用情景分析评估风蚀因子改变的情况下,植被恢复对防风固沙服务物理量和服务能力的影响。拟解决的关键问题：(1)以定量化研究的结果为长城沿线2000—2018年的风蚀现象的变化情况提供科学依据。(2)探讨风蚀因子改变,植被恢复对风沙边缘地区风蚀量的影响及生态系统防风固沙物质量和能力的影响。明确气象环境改变时,植被恢复与防风固沙服务的作用机制,并推进土地资源的保护和建设,实现可持续高质量发展，也为我国陕北地区防风固沙治理工作进程提供一定的科学参考。

1 研究区及数据源

1.1 研究区概况

陕北长城沿线(图1)位于陕西省最北端,包括定边、靖边、横山、榆阳、神木、府谷六个市县(区)[21—22],北邻鄂尔多斯草原,西邻宁夏自治区。位于36°45′—39°58′N,107°35′—111°29′E之间。作为陕北黄土高原丘陵沟壑区,地处干旱与半干旱地带,地势东低西高、年降水量偏少、土壤类型主要以沙土为主,为风力侵蚀提供了“有利”的发生条件。作为干旱频发区,面临着植被退化、水土流失等生态问题,生态系统脆弱易受人类活动干扰和气候变化的影响。

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

1.2 数据来源

本文使用的数据包括NDVI数据、数字高程模型(Digital Elevation Model,DEM)、土壤数据、雪盖数据和气象数据。(1)NDVI数据来源于美国地质勘探局(http://glovis.usgs.gov/)数据产品,选用2000—2018年时间段,空间分辨率为1km,利用MRT(MODIS Reprojection Tools)最大合成年数据。(2)DEM数据来源于中国科学院资源与环境数据云平台(http://www.resdc.cn)空间分辨率为1km。(3)土壤数据和积雪覆盖数据来源于中国科学院旱区寒区科学数据中心(http://westdc.westgis.ac.cn),土壤数据为世界土壤数据库(Harmonized World Soil Database,HWSD)所提供的中国土壤数据为1∶100万土壤数据库。中国雪深长时间序列数据集提供1979年1月1日到2020年12月31日逐日的中国范围的积雪厚度分布数据,空间分辨率为0.25°。(4)气象数据(风速、气温、降水)来源于时空三级环境大数据平台(https://data.tpdc.ac.cn/zh-hans/)提供的中国区域地面要素数据集(China meteorological forcing dataset,CMFD)[23—25],该数据时间范围为1979—2018年,时间分辨率为3小时,空间分辨率为0.1°,数据格式为NETCDF格式,运用Python语言提取出研究区内349个规则点气象要素,大于实际研究区4个国家气象站点,足以从栅格尺度上真实的反映出研究区内风速大于5m/s天数；(5)土壤湿度数据来源于世界气候中心(https://www.worldweather.cn/),空间分辨率为0.1°。最终将所有数据统一为CGCS_2000_Albers坐标系,利用ArcGIS平台重采样为1km分辨率。

2 研究方法

2.1 Theil-Sen斜率估算与M-K显著性检验结合

为了进一步体现出研究区植被变化的显著性,采用斜率估算和显著性检验相结合的方法。而Theil-Sen斜率估算和Mann-Kendall显著性检验是两种非参数检验方法,目前已经得到广泛的应用[26]。因其计算效率高,能够有效地测定时间变化趋势的起始位置,并具有检测范围广、定量化程度高等优点,比一元线性回归趋势检验空间变化更为敏感。Theil-Sen斜率估算通过计算序列中连续两期时间序列数据之间的斜率,将所有像元的数值对斜率的中值做时间序列的总变化趋势。公式如下：

(1)

式中,β是所有数据对斜率的中值。其大小表示平均变化率,正负表示时间序列的变化趋势。Kendall()为取中值函数。xj和xi分别表示时间序列中的第j项和第i项的值。但由于该斜率估算方法不能实现时间序列的趋势显著性判断。通常使用Mann-Kendall显著性检验来判断时间序列的趋势显著性。

Mann-Kendall方法是由彭曼首次提出并使用,并由Kendall进行改进,对于序列X=(X1,X2,X3...Xn)先确定所有对偶值,之后比较与检验统计量S的大小关系。该方法的原始假设是：H0：时间序列中的数据随机排列,即不存在显著趋势；H1：时间序列数据存在上升或下降的单调趋势。公式如下：

(2)

(3)

由于时间序列维度大小不同,显著性检验统计量的选取有所不同。当时间序列小于10时,直接使用统计量S进行双边趋势检验。在给定显著性水平α下,如果|S|≥Sα/2则拒绝原假设,则存在显著趋势。反之则不存在显著性趋势。当时间序列≥10时,统计量S近似服从标准正态分布,可将S标准化后得到检验统计量Z,Z值计算公式如下：

(4)

(5)

式中,n是时间序列中数据个数；m是序列中重复出现的数据组个数；ti是第i组重复数据组中的重复数据个数。同样,采用双边趋势检验,在给定显著性水平α下,在正态分布表中查得临界值Z1-α/2。如果|Z|≤Z1-α/2时,则接受原假设即趋势不显著；反之,则相反。

2.2 生态系统防风固沙服务物质量评价

植被所引起的风蚀减少量可视为防风固沙服务的物理量,其表示在裸土条件下潜在风蚀量与植被覆盖条件下实际风蚀量之间的差[27]。在充分考虑气候条件、植被覆盖、土壤等要素情况下,可根据RWEQ模型定量评估潜在风蚀量和实际风蚀量。利用两者的差值进行生态系统防风固沙服务物理量的评价[28]。其计算公式如下：

(6)

Qrmax=109.8×(WF×EF×SCF×K′)

(7)

Sr=150.71×(WF×EF×SCF×K′)-0.3711

(8)

(9)

Qmax=109.8×(WF×EF×SCF×K′×C)

(10)

S=150.71×(WF×EF×SCF×K′×C)-0.3711

(11)

G=SLr-SL

(12)

式中,SLr表示单位面积年潜在风蚀量(kg/m2)；Qrmax为潜在风力最大输沙能力(kg/m)；Sr为潜在关键地块长度(m)；SL表示单位面积年实际风蚀量(kg/m2)；Qmax为风力最大输沙能力(kg/m);S为关键地块长度(m)；G表示单位面积防风固沙量(kg/m2);z表示下风向距离(m),本文计算取50m；WF表示风蚀因子(kg/m)；EF表示土壤可蚀性因子(无量纲)；SCF表示为土壤结皮因子(无量纲)；K′表示为土壤粗糙度因子(无量纲)；C表示植被因子(无量纲)。具体表达式如下：

风蚀因子WF在考虑风速作用的前提下,还将积雪覆盖对地表的保护和土壤湿度纳入考量,计算公式为：

WF=wf×(ρ/g)×SW×SD

(13)

wf=u2(u2-u1)2×Nd

(14)

式中,wf表示风力系数(m/s3)；g表示重力加速度(m/s2),本文取9.8 m/s2；ρ表示空气密度(kg/m3),本文中采用20℃时的空气密度值,即1.205kg/m3；SW表示土壤湿度因子(无量纲)；SD表示积雪覆盖因子(无量纲),即无积雪覆盖的天数与研究总天数之比,积雪深度小于25.4mm时,表示无积雪；u1为最小起沙风速,模型中默认为5m/s;u2表示2m高度处的风速(m/s)；Nd为每月风速超过阈值风速的天数。

土壤可蚀性因子EF指在一定土壤理化条件下土壤受风蚀影响大小。土壤结皮因子SCF指在一定理化条件下土壤结皮抵抗风蚀能力的大小。其表达式为：

(15)

(16)

式中,sa为土壤砂粒含量(%);si为土壤粉粒含量(%)；cl为土壤粘粒含量(%)；OM为土壤有机质含量(%)；CaCO3为碳酸钙含量(%)。假定土壤可蚀性因子和土壤结皮因子随时间保持不变,在应用世界土壤数据库的基础上,结合陕北长城沿线风沙边缘区土壤质地特征,参考RWEQ模型内嵌土壤资料表,计算出土壤可蚀性因子和土壤结皮因子。

地表粗糙度因子K′是由地形引起的土地表面粗糙程度对土壤风蚀的影响,公式如下：

K′=cosα

(17)

式中,α表示为坡度。

植被因子C表示一定的植被覆盖对土壤风蚀的抑制程度,公式如下：

C=e-0.0483(FVC)

(18)

(19)

式中,FVC代表植被覆盖度(%),NDVIveg和NDVIsoil分别表示植被覆盖的最大和最小值。本文利用植被指数法,基于像元二分模型计算得出95%和5%的累积频率对应下的植被覆盖最大和最小值[29—30]。

生态系统防风固沙服务功能可以表明植被的实际固沙的能力[31]。利用研究区面积,定量化出防风固沙物理量,以此来评估区域生态系统防风固沙服务能力的大小。计算公式为：

QSR=10×G×A

(20)

ASR=QSR/A

(21)

式中,QSR为研究区防风固沙量(t/a)；A为研究区面积(hm2)；ASR为生态系统防风固沙能力(t/hm2)。

3 结果与分析

3.1 Theil-Sen斜率估算与M-K检验

采用Theil-Sen斜率估算和M-K显著性检验相结合的方法,利用NDVI数据得到的陕北长城沿线风沙边缘区植被变化情况(图2)。结果表明：2000—2018年研究区内98.74%的区域植被变化表现为改善。其中,约有90.34%植被覆盖区域的植被状况呈现出显著改善。结合哨兵二号影像数据对比来看,榆阳区、神木市等因城市扩展,导致植被出现退化。不显著改善集中表现在定边县占比较大,说明风沙边缘地区的植被还需进一步保护。

虽然Theil-Sen斜率估算和M-K检验方法有计算效率高、检测范围广、定量化程度高等优点,比一元线性回归趋势检验空间变化更为敏感。但本文检测到的陕北长城沿线风沙边缘地区植被变化空间分布模式与基于一元线性回归分析方法检测到的植被空间变化模式大体上一致。

3.2 风蚀模型因子时空格局

3.2.1风蚀因子结果

从表1可以看出,2000—2018年风蚀因子呈现出减少-增多-减少的趋势,且风蚀因子均值从24.02m/s波动式降至12.26m/s,整体上相对减小。在2010年风蚀因子均值达到最大,为74.13m/s,2005年最低,为6.32m/s。从各点大于5m/s的数值和天数来看,2010年相对多于其他年份,致使2010年的风蚀因子出现最大值。

表1 2000—2018年风蚀因子变化

在空间分布来看(图3),风蚀因子分布差异明显,且风蚀因子高值分布在榆阳区和神木县周边。对比分析,2010年较其他年份差异性更为显著。在2005年定边县风蚀因子高于其他五个县区,统计分析这是由于当年定边县周边气象要素点大于起沙风速阈值的天数和个数偏多。

图2 2000—2018年植被变化趋势及显著性检验Fig.2 Vegetation change trend and significance test from 2000 to 2018

图3 2000—2018年风蚀因子空间分布Fig.3 Spatial distribution of wind erosion factors from 2000 to 2018

3.2.2土壤及地形因子结果

土壤可蚀性因子表征土壤抵抗侵蚀能力的大小,土壤对风蚀作用敏感程度越大,可蚀性越高。陕北长城沿线风沙边缘区EF取值(图4)在0.26—0.61,平均值为0.48。在空间分布上,榆阳地区土壤对风蚀作用较敏感,EF取值较高,而定边县和府谷县EF取值较低。从土壤理化性质出发,红砂土的EF值最高,这是因为其沙粒含量达88%以上,容易发生风蚀现象。潜育土EF值最低且为0.266,因其粘粒和有机质含量达到43.7%相比于其他土壤类型不易起沙。

土壤结皮因子与土壤易受风力侵蚀威胁成反比,即土壤表层结皮越坚固,土壤结皮因子SCF值越大,土壤越不易受到风蚀威胁。整个研究区SCF的取值(图4)在0.06—1.00,平均值为0.56。SCF越高抵抗风蚀能力越差。从空间分布来看,SCF高值分布在长城沿线分割线以北,与内蒙古毛乌素沙漠相邻且抵抗风蚀能力较差。从土壤理化性质来看,红沙土粘粒含量和有机质含量最低,导致SCF最高,为0.837,抗风蚀能力最差。始成土的粘粒和有机质含量相对较高,SCF最低,说明其抗风蚀能力最强,不易发生风蚀。

地表粗糙度因子K′ 反映了地形对风蚀强度大小的影响。长城沿线K′(图4)取值范围为0.004—1.00,均值为0.65。空间分布上,长城沿线分割线以北地表粗糙度占大部分,主要分布在定边县北部、靖边县、榆阳区、神木县。根据海拔高度和地形来看,定边县和靖边县海拔北低南高,榆阳区、神木县西北高东南低,且两两县区高值区和低值区相近,造成地表粗糙度高,更易于侵蚀的发生。且定边县和靖边县南部山地海拔较高和起伏较大,不利于侵蚀的发生。

图4 土壤因子分布Fig.4 Distribution of soil factors

3.2.3植被因子结果

植被因子表征植被对风蚀抑制程度的大小,植被生长越好对风蚀的抑制程度越大进而影响生态系统防风固沙能力。研究区C因子取值范围为0.007—1.00,植被整体呈现好转趋势但均值呈轻微下降趋势,局部存在波动。空间分布上(图5),C因子较高值主要分布在神木市和榆阳区等地,表现出易被风力侵蚀,C因子较低值则相反。

图5 2000—2018植被因子Fig.5 The vegetation factor from 2000 to 2018

3.3 防风固沙物质量及时空格局

3.3.1潜在与实际风蚀时空变化

潜在风蚀能够反映出土壤可能的最大风蚀量,即土壤风蚀危害的“危险程度”。从表2可看出,2000—2018年潜在风蚀量在4569.18—64164.44万t,2010年达到潜在风蚀量峰值,2005年最少。与2000年相比,2018年减少了13160.93万t。可见陕北长城沿线土壤风蚀危害的“危险程度”在呈波动式降低。与2000年相比,2018年平均单位面积潜在风蚀量下降了3.92kg/m2,总体呈“减少-增加-减少”的波动式变化。2000—2018年实际风蚀量在1056.31—24211.89万t,在2010年达到实际风蚀量峰值(24211.89万t),2005年最少。与2000年相比,2018年增加了894.94万t。与2010年相比,2018年减少了20161.58万t,平均单位面积风蚀量减少了6.0kg/m2,减少幅度较大,整体呈“减少-增加-减少”的变化特点。

2000—2018年潜在风蚀出现显著的空间差异性。根据各个因子空间分布,在风蚀因子分布较大、地表较粗糙、土壤结皮较大的地区潜在风蚀强度最为剧烈且变化幅度最为明显,如榆阳区和神木市。而定边县、靖边县和府谷县潜在风蚀强度较轻。实际风蚀年际空间变化上,也存在明显的物质量变化。依据水利部《土壤风蚀分类标准(SL190—2007)》将土壤风蚀分为微度(<200t/a)、轻度(200—2500t/a)、中度(2500—5000t/a)、强烈(5000—8000t/a)、极强烈(8000—15000t/a)、剧烈(>15000t/a)[32—33]。根据分类统计,微度土壤风蚀面积增加了16.16%,轻度风蚀面积减少了15.94%,2010年强烈、极强烈、剧烈风蚀面积占比最大,与此对应该年的风蚀因子最大,导致风蚀现象严重。2010到2018年,风蚀因子逐渐减小,土壤风蚀强度也逐渐有了明显的改善。定边县、靖边县和府谷县微度风蚀强度占比较多；榆阳区和神木市潜在风蚀最为剧烈且变化幅度最为明显,也是单位面积实际风蚀量较严重的区域(图6)。从时间统计和空间分布比较来看,每年的潜在风蚀量>防风固沙量>实际风蚀量,说明在最大危害程度内,防风固沙服务发挥明显作用。而从整体情况来看,在2000—2018年陕北长城沿线防风固沙治理过程中,土壤风蚀现象已得到显著遏制。相比于潜在风蚀,实际风蚀多考虑植被因子,由于其他因定不变,主要是风蚀和植被因子对风蚀的影响最密切。

表2 2000—2018年防风固沙物质量统计/(kg/m2)

图6 2000—2018年风蚀强度空间分布Fig.6 Spatial distribution of wind erosion intensity from 2000 to 2018

3.3.2防风固沙服务物质量时空变化

潜在与实际风蚀差值表征土壤风蚀模数(单位面积防风固沙物理量)。研究表明：防风固沙物质量从2000年的21944.56万t减少到2018年的7987.44万t。2010年变化波动增长较大且达到最大固沙量。依据张彪等[31]将防风固沙服务能力ASR分为：低值区(ASR≤5t/hm2)、较低区(5≤ASR≤20t/hm2)、一般区(20≤ASR≤60t/hm2)、较高区(60≤ASR≤120t/hm2)和高值区(ASR≥120t/hm2)。从防风固沙服务能力分区统计来看,如表3所示,2000—2018年生态系统防风固沙服务能力高值区和较高值区出现减小-增多-减小的变化特点,且整体上是减小的趋势,且对比观察不同年份的防风固沙强度的分布规律即存在相似又存在差异。长城沿线生态系统防风固沙服务能力表现出明显的空间异质性,防风固沙服务能力高值区主要集中在神木市和榆阳区,如图7所示。根据模型分析,如果风蚀因子偏大且植被因子降低,即风场强度增加导致防风固沙量显著减少。风速的不稳定性会导致风蚀因子结果出现较大偏差,进而对防风固沙服务的空间分布变化产生较大的影响。在抑制土壤风蚀现象发生的基础上,固沙能力也在减小,后期需要加以重点保护植被。

表3 2000—2018年防风固沙能力空间分布/km2

图7 2000—2018年防风固沙能力空间分布Fig.7 Spatial distribution of wind prevention and sand fixation capacity from 2000 to 2018

3.4 植被恢复对土壤风蚀及防风固沙服务影响

为了进一步探究风蚀因子变化,植被恢复对风蚀现象及生态系统防风固沙服务的影响,利用情景模拟法,以最小(2005年)、最大风蚀(2010年)因子为基准,分设两种情景(即情景1和2)。由于计算潜在风蚀时没有考虑到植被覆盖,其他因子固定不变,即：2005年的潜在风蚀为情景1条件下的潜在风蚀,2010年的潜在风蚀为情景2条件下的潜在风蚀。从物质定量化分析角度出发(表4),横向对比情景1得出：最小风蚀因子情况下,植被长时间序列恢复,植被因子降低,抵抗风蚀能力增强,实际风蚀量从1066.56万t减少到825.38万t,相比减少了241.18万t,防风固沙量在逐年提升的同时防风固沙服务能力也显著提升。与情景1相比,情景2中风蚀因子的改变使得土壤风蚀量从2000年的11764.29万t波浪式增加到2018年的22649.44万t,说明植被不再是影响风蚀的主要因素,计算风蚀因子的土壤湿度、风力系数等因子的改变使得地表更不易抵抗风蚀。但纵向对比情景1和2：随着风蚀因子的显著变化,土壤风蚀危害的“危险程度”达到最大(潜在风蚀),以2000年为例,实际风蚀量从1066.56万t显著增长到11764.29万t,防风固沙量从3449.6万t变到52175.12}。由此可见：其他因素不变,气象条件恶劣的情况下,极易发生土壤风蚀。但两种情境下的每年防风固沙量>实际风蚀量,利用植被固沙的效果也越来越明显,生态系统防风固沙服务功能也在增加。

表4 2000—2018年分情景下防风固沙物质量统计/(kg/m2)

空间分布上(图8),在风蚀因子条件下(情景1),植被恢复,2000—2018年风蚀现象逐渐退化。在最大风蚀因子下(情景2),造成风蚀现象严重,且主要集中在榆阳区和神木县,导致定边县和靖边县受到殃及。这四个县区的防风固沙量也出现明显差异。可见风蚀因子是主导土壤风蚀的主导因素。整体来看：采取“锁绿”来防沙治沙是最为有效的,且随着风蚀因子的改变水土流失治理成效越为明显,生态系统防风固沙服务能力也越强。植被恢复对长城沿线的土壤风蚀起到有效的抑制,对生态系统防风固沙服务起到了一定的增强作用。但风蚀因子中风速属于不可防控的气象因素,只能依靠恢复植被来增强生态系统防风固沙能力。

图8 分情境下风蚀空间分布Fig.8 The amount of wind erosion under different scenarios

4 结论

本文基于RWEQ模型分析了陕北长城沿线防风固沙服务物理量的时空模式,实现风沙边缘地区生态系统防风固沙服务的量化,并评估风蚀和植被因子变化下对生态系统防风固沙服务物理量和服务能力的影响。得到以下结论：(1)通过分析2000—2018年陕北长城沿线六个县区年际NDVI的时空变化,风沙区内植被呈明显增加趋势；(2)2000—2018年整个风沙区实现了从“沙进人退”到“绿进沙退”的转变,土壤风蚀危害的“危险程度”也在逐步降低,实际风蚀和生态系统防风固沙总量2005年和2010年分别为4569.18万t和64164.44万t,且潜在风蚀量>防风固沙量>实际风蚀量,生态系统防风固沙服务发挥显著作用。(3)利用情景变化定量分析,探索植被恢复可以一定程度上抑制风蚀现象。在最小风蚀因子情况下,不受外力条件变化,植被防风固沙服务的作用发挥不明显,且一定程度抑制风蚀现象的发生。风蚀因子改变加重研究区的土壤风蚀,植被所发挥的生态系统防风固沙服务越为明显。而两因子相结合(即实际情况下的风蚀量),风蚀因子占主导地位,对土壤风蚀和防风固沙量影响作用较大。不论是真实情况还是分情景下的风蚀,生态系统防风固沙服务空间差异显著,榆阳区和神木市的土壤风蚀现象最为严重,也是生态系统防风固沙服物质量和能力高值集中区。总之,恢复植被是实现区域性荒漠化逆转和改善生态系统防风固沙服务与能力的必要条件,恢复植被不仅可以减少土壤风蚀,也可以增强生态系统防风固沙服务能力。

5 讨论

(1)当前土壤风力侵蚀的定量化验证还存在一定的困难。实测的土壤风蚀量是开展研究结果验证的基础,由于没有实测的土壤风蚀量数据,无法对利用模型模拟的实际风蚀量结果进行验证,后续可进一步采用沙尘暴发生频率来验证本文的可靠性。鉴于此,文中利用情景模拟法可以间接进行验证,对比分析植被恢复对风蚀现象有一定的遏制作用,与刘硕等[16]的研究结果一致。在改变风力因子时,风蚀现象依旧明显,植被发挥的生态系统防风固沙服务也显著提高。

(2)提高模型的优化性。在实现陕北长城沿线2000—2018年生态系统防风固沙量以及服务能力,由于实测气象站点少,而采用中国区域地面要素数据集,但该数据集时间范围限制,未能及时做到最新年份；考虑到该数据集时间分辨率较大,且有349个规则分布的点,在选用空间分辨率时采用1km空间分辨率,但对于土壤风蚀研究而言,1km分辨率足以反映其空间差异。在实现雪盖因子计算的过程中,计算的是年积雪覆盖天数的概率,得到年尺度风蚀量。在后期研究细化时间尺度,从月尺度考虑生态系统防风固沙物质量和服务能力。

(3)本文的研究结论仅从生态系统防风固沙量和服务能力角度出发,在情景1条件下植被的固沙作用效果不明显,但在最大风蚀因子条件下,植被对抑制风沙吹动,防风固沙效果越来越明显。三种情景每年的潜在风蚀量>防风固沙量>实际风蚀量。风蚀因子增大,潜在风蚀增加、实际风蚀增加、防风固沙量增加,生态系统防风固沙服务功能也增强。这与治沙研究者提出防沙治沙、植树造林工程以来的预想一致。从“沙进人退”到“绿进沙退”的转变是数十年治沙的功劳,“榆林绿”锁住了毛乌素沙地。作为自然地貌,毛乌素沙地不是真的消失,只是流动沙丘全部得到固定。但这些治理成果还属人工干预,尚需加快沙区生态系统的正向演替,实现自我循环发展,从根本上改变沙区生态面貌。

本文研究特色在于：采用中国区域地面要素数据集,更好的实现长城沿线风蚀因子的空间分布。但还存在不足之处：在考虑土壤理化性质随时间不发生变化的情况下,但实际情况中,由于陕北长城沿线六个县区耕地面积占比较大,长期使用化肥致使土壤理化性质发生轻微改变,导致文章实验结果的可靠性存在一定的误差,土壤理化性质的改变对土壤风蚀具有一定的影响,对风沙区土壤理化性质的测定研究也是今后的研究的一个重要科学问题。在可控的人为范围内,做好治沙工程,坚持不懈地防沙治沙、植树造林使得生态环境持续改善,为榆林经济社会的快速发展提供了强有力的基础保障。