防风固沙服务指植被生态系统对风沙起到的抑制和固定作用^[1], 是干旱地区荒漠生态系统提供的最为重要的服务^[2]。中国是世界上荒漠化最严重的国家之一, 约有660万km²旱地面临荒漠化的风险^[3]。第五次全国荒漠化和沙化状况公报显示, 截至2014年, 我国荒漠化土地面积261.16万km², 其中沙化土地面积达172.12万km²。风蚀是造成干旱、半干旱及半湿润地区土地荒漠化的主要因素^[4], 严重的风蚀荒漠化不仅造成流沙堆积、地表风蚀粗化、土壤水分散失、土地生产力下降和沙尘暴灾害频发等生态环境问题^[5], 同时对人民的生命健康和财产安全产生了严重威胁^[6—7]。荒漠化地区生态脆弱, 为治理土壤侵蚀、遏制土地退化, 中国实施了多样化的治理措施和管理政策^[8], 先后开展了一系列生态恢复工程措施, 如三北防护林、退耕还林还草、天然林保护工程、京津冀风沙源治理工程等, 有效改善区域生态环境, 提高了生态系统服务功能^[9—11]。然而全球气候变化和人口的快速增长将会增加荒漠化的风险^[12], 科学定量评估区域防风固沙功能、开展长时间序列防风固沙服务研究可为应对荒漠化、促进生态修复和可持续发展提供理论依据和参考。

目前, 涉及植被生态系统防风固沙服务的研究主要包括物质量与价值量的评估和变化^[2]、空间流动^{[1, 13]}以及防风固沙服务与其他服务间权衡协同关系^[14—15]等三个方面, 其中防风固沙服务空间流动、服务间权衡协同关系的研究较多, 主要侧重于分析区域生态系统所提供的服务功能效益, 而物质量评估和变化方面的研究作为基础性工作较多集中于区域土壤风蚀和防风固沙服务时空演变, 对于驱动机制的研究尚不完善, 近年来已受到诸多学者广泛关注。Teng等^[16]从气候变化和人类活动两方面研究得出气候变化对青藏高原风蚀的影响具有大尺度和空间连续性, 而人类活动对风蚀的影响只发生在局部地区并呈斑块分布; Li等^[17]基于气候和社会经济因子采用多元线性回归方法分析得出温度和造林工程是内蒙古防风固沙服务功能提高的最重要因素; 彭婉月等^[18]定量评估了2000—2017年黑河中下游防风固沙服务变化并基于灰色关联法分析得出风力因子是影响防风固沙功能变化的最主要因子; 张玥等^[19]应用地理探测器模型分析了自然和人为因素对锡林郭勒盟防风固沙功能的贡献及交互作用。以上研究集中于各类因子对防风固沙服务的影响程度和交互作用, 虽进一步加深了对于防风固沙服务的理解, 但这些研究方法多从全局性角度出发, 侧重各因子的整体性影响, 掩盖了变量间关系的局部特性, 忽略了各驱动因子作用的区域差异和空间异质性^[20—21]。地理加权回归模型(Geographically weighted regression, GWR)^[22]是一种局域空间分析方法, 可用于检验回归结果的空间变异性, 揭示自变量与因变量之间的空间非平稳性, 有助于更加全面了解防风固沙服务的驱动机制。

黄河流域是我国重要的生态屏障和经济区, 其生态环境脆弱, 沙漠化土地广泛分布于该区域^[23]。黄河流域生态保护和高质量发展已被确定为国家战略, 然而生态保护与发展矛盾突出, 荒漠化问题一直是黄河流域生态保护修复的突出问题^[24]。因此, 本文以黄河流域为研究区, 基于修正土壤风蚀模型(Revised wind erosion equation, RWEQ)计算2000—2018年黄河流域防风固沙服务, 探究其时空变化特征, 结合相关分析和GWR模型揭示研究区防风固沙服务驱动因素影响的区域差异和空间异质性, 以期为黄河流域荒漠化防治、生态保护修复和可持续发展提供一定的科学依据和理论参考。

1 研究区概况与数据来源 1.1 研究区概况

黄河是我国的第二长河, 发源于青藏高原巴彦喀拉山北麓, 自西向东流经青海、四川、甘肃、宁夏、内蒙古、陕西、山西、河南、山东9省区, 最后注入渤海。黄河流域位于我国中北部(32°10′—41°50′N、95°53′—119°05′E), 总面积约79.5万km², 地势西高东低, 地形差异较大, 自西向东横跨青藏高原、内蒙古高原、黄土高原以及黄淮海平原4个地貌单元(图 1)。流域气候复杂多样, 多年平均气温约9.5℃, 多年平均降雨量约466.5mm^[25], 自北向南横跨干旱、半干旱、半湿润3个气候分区。黄河流域生态环境脆弱, 沙漠化土地广泛分布于该区域, 流域内存在库布奇沙漠和毛乌素沙漠。

1.2 数据来源及处理

本研究采用多源数据集, 包括气象数据集、土壤数据集、社会数据集以及土地利用、DEM、NDVI等相关辅助数据集, 具体数据来源见表 1。

其中温度、降水、风速等气象数据首先利用Matlab检验去除缺失严重的站点, 对部分站点进行插补后以DEM为协变量, 利用ANUSPLIN专业气象插值软件进行插值, 空间分辨率为1km; 粗砂、粉砂、黏粒和有机质等土壤属性数据来自面向陆面模拟的中国土壤数据集, 在Arcmap10.2中栅格计算相应土壤因子, 空间分辨率为1km; 2000—2018年雪深数据采用Arcpy进行批量格式转换后, 用Matlab计算雪盖因子, 空间分辨率为1km; 归一化植被指数数据集已经过最大值合成处理, 有效去除大气云层、建筑阴影和太阳高度角的影响, 本研究以像元二分法计算植被覆盖度应用于植被因子计算, 空间分辨率为1km; SRTM V4.1 DEM数据经过镶嵌裁剪后用来提取地表粗糙度因子, 空间分辨率为1km; 土地利用数据经裁剪后将二级分类合并为一级类, 空间分辨率为1km; 各因子空间参考设置为WGS_1984_Albers, 以此进行研究区防风固沙服务研究。

2 研究方法 2.1 RWEQ模型

参考《生态红线划定指南》和相关学者研究^{[16—19, 32]}, 本文采用修正土壤风蚀模型(RWEQ)计算黄河流域土壤风蚀模数和防风固沙服务物质量。RWEQ模型在考虑气候、植被、地形和土壤等因素的基础上, 以裸土条件下的潜在土壤风蚀量与植被覆盖条件下实际土壤风蚀量的差值定量评估黄河流域防风固沙服务功能量的变化, 计算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中, SP为单位面积固沙量(kg/m²); Z表示下风向距离(m), 取50m^{[1, 32—33]}; S_LP为潜在风蚀量(kg/m²); S_L为实际风蚀量(kg/m²); Q_maxP为潜在风沙最大转运量; Q_max为实际风沙最大转运量; S为关键地块长度(m); S_P为潜在关键地块长度(m); WF为气候因子(kg/m); EF为土壤可蚀性因子(无量纲); SCF为土壤结皮因子(无量纲); C为植被覆盖因子(无量纲); K′为地表粗糙度因子(无量纲)。

RWEQ模型中各因子计算公式如下：

(1) 气候因子WF综合考虑风力、降雪、温度、降雨、太阳辐射等多种气候条件对土壤风蚀产生的影响, 计算公式如下：

(8)

(9)

(10)

式中, Wf为风力因子(m³/s³); g为重力加速度(取9.8m/s²); ρ为空气密度(kg/m³), 由海拔高度H(km)和绝对温度T(开尔文)计算得到; SW为土壤湿度因子; SD为雪盖因子; u₁为距地面2m处的临界起沙风速(m/s); u₂为距地面2m处风速(m/s); N_d为测定风速的时间周期; N为测定周期内风速次数。

(2) 土壤可蚀性因子EF表示一定理化条件下土壤受风蚀影响的程度^[1], 计算公式如下：

(11)

式中, sa为土壤粗砂含量(%); si为土壤粉砂含量(%); cl为土壤黏粒含量(%); OM为土壤有机质含量(%); CaCO₃为碳酸钙含量(%), 本文不考虑。

(3) 土壤结皮因子SCF指一定理化条件下土壤结皮抵抗风蚀能力的大小^[1], 计算公式如下：

(12)

式中, cl为土壤黏粒含量(%); OM为土壤有机质含量(%)。

(4) 在RWEQ模型中, 地表粗糙度因子K′主要指农田耕作使地表条件发生变化, 从而对土壤风蚀产生影响, 分为随机糙度C_rr和土垄糙度K_r。由于研究区空间跨度较大, 土地覆被类型多样, 计算地表粗糙度的随机糙度因子和土垄糙度因子获取困难, 因此本文参考江凌等^[32]的研究, 采用土垄造成的地表糙度因子K_r以地形起伏产生的地表粗糙度予以替代, 计算公式如下：

(13)

(14)

式中, K_r为地形糙度因子(cm); C_rr为随机糙度因子(cm), 本次研究取0;L为地势起伏参数(m); ΔH为距离L范围内的海拔高程差(m)。

(5) 植被覆盖因子C是影响土壤风蚀过程的重要因子, 一定程度的植被覆盖可以降低有效风速, 削弱风力搬运土壤颗粒的能力^{[1, 18, 32]}, 计算公式如下：

(15)

式中, SC为植被覆盖度(%); a_i为不同植被类型的系数, 分别为林地0.1535、草地0.1151、灌丛0.0921、裸地0.0768、沙地0.0658、农田0.0438。

2.2 防风固沙功能保有率

单位面积防风固沙量可以表示生态系统中植被作用引起的风蚀减少量, 但由于实际风蚀量与潜在风蚀量均与风蚀力有关, 并不能有效表征生态系统本身对固沙的贡献率, 因此, 将单位面积防风固沙量SP与潜在土壤风蚀量S_Lp的比值作为防风固沙服务保有率F^[34]来表征黄河流域植被防风固沙功能。公式如下：

(16)

2.3 地理加权回归模型

地理加权回归(GWR)是一种局域空间分析方法, 是对普通线性回归的拓展, 通过获取局部参数评估自变量与因变量关系在空间尺度上的变异^[35]。本文引入GWR模型探究黄河流域防风固沙生态系统服务的驱动机制, 并分析其空间异质性分布成因。公式如下：

(17)

式中, y_i为样点i的因变量, (u_i, v_i)为样点的中心位置坐标, α₀(u_i, _i)为截距, α_k(u_i, v_i)为第i个点的回归系数, x_ik为第i个点上第k个变量的观测值, ε_i为误差项。

3 结果分析 3.1 土壤风蚀及防风固沙服务时空变化特征

2000—2018年黄河流域土壤风蚀模数以及防风固沙量呈明显下降趋势, 而植被覆盖度和防风固沙保有率呈现上升趋势(图 2)。土壤风蚀模数呈先上升后下降的趋势, 由2000年的4.93 t hm^-2 a^-1下降到2012年的2.55 t hm^-2 a^-1, 后上升至2018年的4.19 t hm^-2 a^-1, 19年来平均以0.04 t hm^-2 a^-1的速率下降; 而黄河流域防风固沙物质量从2000年的28.05 t hm^-2 a^-1下降到2018年的25.37 t hm^-2 a^-1, 最低值为2011年的18.41 t hm^-2 a^-1, 多年来平均以0.14 t hm^-2 a^-1的速率下降。土壤风蚀模数与单位面积防风固沙物质量的数值差距较大, 说明黄河流域荒漠化治理成效显著。防风固沙物质量的变化受到风速、植被覆盖等气候和人为因素的双重影响, 物质量的下降并不能说明黄河流域防风固沙能力减弱。防风固沙保有率与植被覆盖度的时间变化均处于上升趋势, 两者分别以0.05%·a^-1和0. 26%·a^-1的速率增加, 说明2000—2018年黄河流域植被防风固沙功能增强。

通过对2000—2018年19年数据计算平均值得到黄河流域土壤风蚀模数和防风固沙服务空间分布(图 3), 土壤风蚀模数多年平均值为3.56 t hm^-2 a^-1, 而防风固沙物质量19年平均值达23.08 t hm^-2 a^-1, 黄河流域土壤风蚀和防风固沙服务分布表现出明显的空间异质性。整体上看, 黄河流域土壤风蚀以微度(< 2 t hm^-2 a^-1)和轻度(2—25 t hm^-2 a^-1)为主占研究区面积的96.28%, 侵蚀强度为中度以上(>25 t hm^-2 a^-1)的土壤风蚀区域占比较少, 主要集中于黄河流域西北部, 侵蚀面积随侵蚀强度的增加而减小。而防风固沙服务高值区则主要分布于黄河流域西北部的内蒙和宁夏地区、以及黄河下游的河南地区。

从2000—2018年黄河流域土壤风蚀模数和防风固沙变化趋势空间分布来看, 黄河流域土壤风蚀模数增加区域多于减少区域, 增加区域占研究区面积的22.22%, 零星分布于流域各处, 主要包括甘肃、青海、宁夏西部、阴山山脉、太原盆地以及河南部分区域; 而风蚀模数减少区域占研究区面积的12.91%, 主要分布于流域西北部以及东部部分区域, 包括内蒙古中南部、宁夏东部、陕西北部、山西北部以及山东济南、东营等地。研究区单位面积防风固沙量增加区域多于减少区域, 增加区域占黄河流域面积的17.12%, 分布于流域西部、南部以及东部部分区域, 包括甘肃、宁夏南部、陕西榆林南部、关中平原、山西南部以及河南部分区域; 单位面积防风固沙量减少区域占研究区面积的15.06%, 主要分布于流域西北部和黄河下游, 包括内蒙古、山西中北部以及山东济南、东营等地。

3.2 黄河流域防风固沙服务驱动机制

土壤风蚀受到人类活动、气候、植被、土壤、地形等多种因素综合影响^{[17, 36]}。参考相关学者的研究^{[16—19, 32—33, 37—38]}, 基于防风固沙服务相关影响因子以及数据的可获取性, 本文从社会、气候、土壤、植被、地形等角度选取人口密度(Pop)、夜间灯光(NTL)、降水(Prec)、温度(Temp)、风速(Wind)、气压(Pres)、实际蒸散量(ET)、太阳辐射总量(SR)、归一化植被指数(NDVI)、土壤粗砂含量(SA)、土壤粉砂含量(SI)、土壤黏粒含量(CL)、土壤有机质含量(OM)、坡度(Slope)、高程(DEM)等15种因素作为自变量, 以防风固沙物质量为因变量来探究黄河流域防风固沙生态系统服务的驱动机制。

3.2.1 相关性分析

相关性是筛选模型解释变量的初步条件, 在进行地理加权回归分析之前先对自变量和因变量进行相关分析, 初步筛选出降水、温度、风速、气压、实际蒸散量、太阳辐射总量、归一化植被指数、土壤粗砂含量、土壤粉砂含量、土壤黏粒含量、土壤有机质含量、坡度、高程等13个与防风固沙物质量在0.001级别置信水平显著相关的因子作为地理加权回归模型的解释变量(图 4)。因子间多重共线性会导致地理加权回归模型结果出现偏差, 应选用共线性较低, 即方差膨胀因子(Variance inflation factor, VIF) < 7.5的因子进行地理加权回归分析, 因此本文基于R语言car包进行冗余因子的剔除最终选择降水、温度、风速、气压、实际蒸散量、归一化植被指数、土壤粗砂含量、土壤黏粒含量、土壤有机质含量、坡度等10个因子作为驱动因素进一步分析影响黄河流域防风固沙物质量空间异质性分布的因素。

根据皮尔逊相关分析结果可知, 与防风固沙物质量正相关的驱动因素为：土壤粗砂含量(SA)、温度(Temp)、风速(Wind)、气压(Pres), 与防风固沙物质量呈负相关的驱动因素为：土壤黏粒含量(CL)、归一化植被指数(NDVI)、实际蒸散量(ET)、坡度(Slope)、土壤有机质含量(OM)、降水(Prec)。

3.2.2 地理加权回归分析

地理加权回归分析之前需要诊断模型拟合效果, 对筛选出的驱动因子分别利用全局最小二乘法(OLS)模型与GWR模型构建回归模型。由表 2可知, GWR模型拟合结果中AICc值低于OLS模型拟合结果, 且差值大于3, 这表明GWR模型模拟结果优于OLS模型。GWR模型校正后的局部决定系数R²达0.6070, 相较传统的OLS模型提高了48.45%, 并且GWR模型残差平方和低于OLS, 这说明GWR模型相较于OLS模型能更好地解释各驱动因素对黄河流域防风固沙服务空间分布异质性的影响。

图 5为各驱动因子地理加权回归系数拟合结果, 回归系数可以反映各因子驱动作用强度的空间差异性。局部决定系数R²在-0.14—0.72间, 高值主要集中于黄河流域中北部, 说明GWR模型在黄河流域中部和北部的拟合精度更高。各因子回归系数空间分布具有非均匀性和局部性, 表明各因子对黄河流域不同区域防风固沙服务作用强度呈明显异质性。

降水量对防风固沙服务的影响以负面效应为主, 回归系数小于0的区域占研究区面积的67.22%, 主要分布于黄河流域西部和中游区域, 而在黄河中游流域西北部和下游流域则主要表现为正向促进。温度对防风固沙服务的影响以负面效应为主, 研究区77.59%区域的防风固沙服务受到温度的负面抑制作用, 集中分布于黄河流域上游和中游。风速对防风固沙服务的影响以正面效应为主, 占研究区总面积的60.53%, 风速越高, 防风固沙量也越高。气压对防风固沙服务的影响与风速类似, 以正向促进作用为主。实际蒸散量对防风固沙服务的变化产生正负两面效应, 呈多中心分散分布, 回归系数范围为-2.27—1.51。归一化植被指数对研究区防风固沙服务的分布产生正负两面效应, 在黄河流域南部NDVI对防风固沙服务的影响主要表现为微弱的促进作用, 而在黄河流域北部和东部则表现为抑制作用。土壤黏粒含量对防风固沙服务的影响以负面效应为主, 约占研究区总面积的82.50%。土壤粗砂含量对防风固沙服务的影响以正面效应为主, 约占研究区总面积的81.57%, 仅在黄河流域中下游部分区域表现为微弱的抑制作用。85.68%区域土壤有机质含量对防风固沙服务为负面影响, 仅在黄河流域中北部零散表现为促进作用。坡度对防风固沙服务的影响表现出明显的负面效应, 回归系数范围为-2.66—0.12, 小于0的区域占研究区面积的94.18%, 仅在黄河流域中游部分区域零星表现为正面效应。

4 讨论

本文采用RWEQ模型定量评估了2000—2018年黄河流域土壤风蚀模数以及防风固沙服务的时空变化并基于相关分析和GWR模型探究了防风固沙物质量空间分布异质性驱动因素。黄河流域土壤风蚀模数与单位面积防风固沙量均呈下降趋势, 主要与风蚀力的减弱有关, 在研究时段内研究区风蚀力从19.36 m³/s³下降至18.89 m³/s³(图 6), 风蚀力减弱导致潜在土壤风蚀量与实际土壤风蚀量均下降, 防风固沙量也呈现下降趋势, 但这并不能说明生态系统防风固沙功能在降低, 从防风固沙保有率以及植被覆盖度的的变化趋势来看, 生态系统植被对固沙的贡献率整体是上升, 这与前人的研究相一致^{[9, 39]}。相关分析结果表明研究区防风固沙服务主要受到气候、植被、地形、土壤等自然因素的作用。

驱动因子间存在相互作用, 并不是单独对防风固沙服务产生影响, 并且各驱动因素对防风固沙服务影响具有明显的空间差异, 驱动因素的作用导致防风固沙服务呈现出空间异质性。从回归系数整体上看, 降水、温度、风速、实际蒸散量、气压等气候因素对研究区防风固沙服务的影响普遍较强, 气候变化是影响区域土壤风蚀的主导因素^[40—41]。降水量与温度不仅直接影响土壤湿度、土壤颗粒粘力, 同时影响植被生长状况, 达到调节风蚀的作用^[42]。降水具有抑制风蚀的能力, 雨水增加了土壤湿度并允许土壤颗粒间形成毛细作用力增强土壤凝聚力^[43], 更多的降水可以促进植被生长, 有利于土壤表面结皮从而增强地表对风蚀的抵抗力^[44], 黄河流域东南部降水量高于西北部, 具有更好的植被覆盖与土壤湿度, 可以有效抵抗风蚀, 促使东南部的实际风蚀与潜在风蚀明显低于西北部, 因而降水对防风固沙量影响在西北部和东南部呈现明显差异。温度对防风固沙服务的影响与降水相似, 体现在温度升高一定程度上促进植被生长抑制风蚀发生, 从而提高防风固沙服务功能, 但同样会导致蒸散量升高土壤湿度降低, 使得地表更加干燥, 容易遭受风蚀^[17]。风是产生风蚀的重要气候驱动力^[45], 而风蚀的产生需要达到临界风速, 在风速较大的一些区域如内蒙古南部和榆林市北部, 风蚀模数与防风固沙服务均较高, 而风速大小与水平气压梯度力成正比, 气压场的变化会引起地面风速的改变^[46], 对防风固沙量的影响呈现出相似性。

除气候因素外, 植被覆盖和土壤性质对防风固沙服务也具有较大影响。植被覆盖可增强土壤抗侵蚀能力同时会影响区域环境如温度、土壤湿度、蒸散发等^[47—48], 在研究区南部较高的植被覆盖减少蒸发增加了含水量, 同时增加了土壤有机质含量, 这可能是植被覆盖对防风固沙量的影响在南北区域空间分异的原因。土壤性质对风蚀有显著影响, 相关研究表明有机质含量较低的土壤在自然条件下更容易遭受侵蚀^[49], 黏粒含量和有机质含量较高时土壤颗粒间的黏着力较强^[50]。研究区土壤性质存在明显的空间差异, 中南部土壤黏粒和有机质含量明显高于北部, 因而在中南部的促进作用强于北部。有学者认为地形起伏度与土壤风蚀存在一定关系, 地势平坦地区土壤风蚀较为强烈, 而在坡度较大的区域, 地表粗糙度增大, 增强了土壤抗风蚀能力^[51], 也有研究表明随着坡度上升沙化土地面积呈快速下降趋势^[52], 坡度较低区域防风固沙能力要弱于较高区域, 坡度的负面效应在西北区域表现最为明显。

本文针对黄河流域土壤风蚀及防风固沙服务的研究进行了有益探索和尝试。从模型角度看, 进行研究时多参考相关学者^{[1, 18, 32—33]}的成果对其中一些参数和因子进行处理, 如下风向距离Z因缺乏实测数据而参考相关研究^{[1, 32—33]}取值50m, 但相较于1km×1km的栅格研究单元而言这样的取值是否合适仍然需要进一步研究; 植被因子C处理时受限于数据, 以土地利用为基础考虑了不同植被类型防风固沙能力不同进行相应计算^[18], 但这与RWEQ模型中植被因子的原始计算方式有所差别。RWEQ模型最初主要是为了估算美国农田的土壤流失状况而创建的^[53], 综合目前研究来看, 该模型在其他土地类型的计算结果与实际情况间存在的不确定性即可能高估或者低估土壤风蚀的问题还需进一步根据实测数据量化研究。从数据角度来讲, 本文所使用的多源数据集中各种数据格式分辨率均不一致, 考虑到研究区范围较大, 最终将数据统一为1km分辨率进行计算, 与高分辨率数据相比存在一定的误差。防风固沙服务受多种因素共同作用, 本研究从多角度选择了10种因子探究黄河流域防风固沙服务驱动机制, 但由于GWR模型不适用于分析土地利用、土壤类型、植被类型等类型量数据, 未将这些数据纳入驱动因素中进行分析。受限于数据的可获取性, 在人类活动方面本次研究仅选取了人口密度和夜间灯光进行相应分析并没有全面考虑例如经济、造林面积、牲畜等等, 这对研究结果会造成一定影响, 同时植被对防风固沙服务的影响并不能简单地归结为自然因素, 大规模生态工程的实施是植被恢复的有力措施, 应当结合自然以及人为影响两方面考虑, 对于自然因素和人为因素耦合影响下防风固沙服务变化将是下一步研究的重点, 如预测未来气候情景下土壤风蚀和防风固沙服务的变化、未来土地利用情景下各土地利用类型防风固沙能力变化。

5 结论

(1) 2000—2018年黄河流域土壤风蚀模数与单位面积防风固沙量均呈下降趋势, 防风固沙保有率和植被覆盖度则呈上升趋势。土壤风蚀模数和单位面积防风固沙量分别以0.04 t hm^-2 a^-1和0.14 t hm^-2 a^-1的速率下降, 而防风固沙保有率与植被覆盖度分别以0.05%·a^-1和0. 26%·a^-1的速率上升, 说明黄河流域风蚀状况减轻, 植被防风固沙功能增强。

(2) 黄河流域多年土壤风蚀模数与单位面积防风固沙量平均值分别为3.56 t hm^-2 a^-1和23.08 t hm^-2 a^-1。黄河流域土壤风蚀以微度(< 2 t·hm^-2·a^-1)和轻度(2—25 t hm^-2 a^-1)为主占研究区面积的96.28%, 侵蚀强度为中度以上(>25 t hm^-2 a^-1)的土壤风蚀区域占比较少, 空间分布呈现西北高东南低的特点, 而防风固沙服务高值主要分布于西北部和下游流域。

(3) 各驱动因素对黄河流域防风固沙服务的影响具有明显的空间异质性。降水、温度、实际蒸散量、土壤黏粒含量、土壤有机质含量和坡度对防风固沙服务的影响均呈现负面效应强于正面效应, 而风速、气压、归一化植被指数、土壤粗砂含量对防风固沙服务的影响则表现出正面效应强于负面效应。