基于RWEQ的毛乌素沙地防风固沙功能及其影响机理研究

Windbreak and sand fixation function and its driving mechanisms in Mu Us Sandy Land based on revised wind erosion equation

摘要

目的 毛乌素沙地作为我国北方生态安全屏障建设的核心区域，研究其防风固沙功能的时空变异特征及驱动机制，为荒漠化防治工程效益评估提供量化依据，为生态恢复区划策略的制定提供理论支撑。 方法 基于修正风蚀方程（Revised Wind Erosion Equation，RWEQ），分析了2000—2023年毛乌素沙地防风固沙功能的时空变化特征，并运用地理探测器定量分析了11个因子对防风固沙功能的影响机理。 结果 （1） 2000—2023年毛乌素沙地年均防风固沙量为3.32×10⁷ t，单位面积年均防风固沙量约3.82 t/hm²，多年来以0.045 t/（hm² · a）的速率上升。（2） 空间上呈现东南高西北低的分布特征，单位面积防风固沙量高值区〔>10 t/（hm² · a）〕主要分布于横山区、靖边县、神木市中部、榆阳区东南部、定边县东部和伊金霍洛旗的中部地区；低值区〔0~2 t/（hm² · a）〕主要分布于乌审旗的北部、鄂托克旗的东南部。（3） 2000—2023年林地、耕地和草地单位面积固沙量变化量分别上升了138.12%，98.56%，95.40%。（4） 单因子分析结果表明：归一化植被指数（NDVI）、降水、土壤类型是对防风固沙功能变化解释力最强的因子。双因子交互作用分析结果表明：NDVI与土壤类型、土地利用类型、风速和高程的交互作用具有更强的解释力。 结论 毛乌素沙地防风固沙功能的显著提升是生态工程和气候变化多因子协同作用的结果，且生态工程和气候变化对区域防风固沙功能的作用机制存在显著差异，未来需要持续关注并优化多因子交互关系以进一步提升防风固沙能力。

Abstract

Objective The Mu Us Sandy Land serves as a core region for ecological security barrier construction in northern China. This study examines the spatiotemporal variation characteristics and driving mechanisms of its windbreak and sand fixation function， providing a quantitative basis for evaluating the effectiveness of desertification control projects and offering theoretical support for the formulation of ecological restoration zoning strategies. Methods Based on the revised wind erosion equation （RWEQ）， the spatiotemporal variation characteristics of windbreak and sand fixation function in the Mu Us Sandy Land from 2000 to 2023 were analyzed， and the driving mechanisms of 11 factors on windbreak and sand fixation function were quantitatively analyzed using the geographical detector method. Results （1） From 2000 to 2023， the annual average windbreak and sand fixation amount in the Mu Us Sandy Land was 3.32×10⁷ t， and the annual average amount per unit area was about 3.82 t/hm²， with an increasing rate of 0.045 t/（hm² · a） over the years. （2） Spatially， the function exhibited higher values in the southeast and lower values in the northwest. High-value areas 〔>10 t/（hm² · a）〕were mainly distributed in Hengshan District， Jingbian County， central Shenmu City， southeastern Yuyang District， eastern Dingbian County， and the central region of Ejin Horo Banner. Low-value areas 〔0~2 t/（hm² · a）〕 were mainly distributed in northern Uxin Banner and southeastern Otog Banner. （3） From 2000 to 2023， windbreak and sand fixation amount per unit area in forest land， cultivated land， and grassland increased by 138.12%， 98.56%， and 95.40%， respectively. （4） Single-factor analysis results showed that the normalized difference vegetation index （NDVI）， precipitation， and soil type were the most influential factors in explaining changes in windbreak and sand fixation function. Two-factor interaction analysis results showed that the interactions between NDVI and soil type， land use type， wind speed， and elevation had stronger explanatory power. Conclusion The significant improvement of the windbreak and sand fixation function in the Mu Us Sandy Land results from the synergistic effects of ecological engineering and climate change， and their respective influence mechanisms show significant differences. In the future， the ongoing monitoring and optimization of multi-factor interactions are required to further enhance windbreak and sand fixation capacity.

Graphical abstract

关键词

防风固沙 / RWEQ模型 / 地理探测器 / 毛乌素沙地

Key words

windbreak and sand fixation / revised wind erosion equation model / geographical detector / Mu Us Sandy Land

沙漠化是全球重大生态问题之一^［1］。毛乌素沙地地处我国北方农牧交错带中部、北方防沙带的西部，它既是京津冀地区的主要沙尘源、也是黄河几字湾风沙防治“攻坚战”主战场、国家重要能源战略基地、北方生态安全屏障重点建设区域^［2］。历史上，毛乌素沙地由于过度放牧、气候变化、战争等多种原因，土地沙化严重，造成流沙堆积、地表风蚀粗化、土壤水分散失、土地生产力下降和沙尘暴灾害频发等生态问题^［3］。因此，其防风固沙功能提升对筑牢北方生态屏障和黄河流域生态保护和高质量发展具有重要意义。为遏制土地沙化，我国政府自1978年以来陆续启动了一系列生态工程，包括“三北防护林工程”“退耕还林”“退牧还草”等，取得了显著的生态功效，具体表现为林草面积占比从1980年的10.2%增长到2020年的73.7%，绿化土地面积从1980年的3.25万 km²增长到 2020年的3.39万 km²，极重度、重度和中度荒漠化面积大幅下降^［4］。防风固沙是生态系统通过其结构与过程减少由于风蚀所导致的土壤侵蚀作用，这种保持土壤，抑制风蚀过程的功能即为防风固沙功能^［5］。人工植被是土地沙化防治最有效的方法。植被能够通过拦截风沙流、改变土壤物质组成等作用，阻滞沙尘和减缓风沙对土壤的侵蚀，起到防风固沙功能^［6］。防风固沙功能受气象、土壤、植被、土地利用等多种因素的影响，土地利用类型的改变是人类活动对自然改造过程的具象化体现，能够直接反映人类活动的方向和强度。生态工程区土地利用变化必然引起区域防风固沙功能的相应变化。生态工程的防风固沙功能与植被覆盖度密切相关，植被的面积和类型变化势必会造成防风固沙生态功能的增减变化^［7］。毛乌素沙地作为国家防风固沙的重点生态功能区，监测其生态工程的防风固沙功能动态变化及其发展趋势，不仅可以准确、及时、全面地掌握荒漠化土地的消长变化信息^［8］，还可为后续生态工程分区施策和设计提供依据。

目前有多种风蚀模型可用于防风固沙功能的评估。包括德克萨斯模型、风蚀方程、Bocharov模型、WEPS （Wind Erosion Prediction System）、 RWEQ （Revised Wind Erosion Equation）等^［9］。由于受研究对象、尺度、适用范围及考虑因子的多寡等，不能兼顾科学性和实用性，使各风蚀模型的理论局限性和适用性差异较大。如德克萨斯模型考虑因子较少，不能全面反映风蚀过程；WEQ （Wind Erosion Equation）模型主要针对农田开发，适用范围有限；Bocharov模型科学性较强，但缺乏各影响因子的定量关系表达^［5］。基于过程模拟的经验性模型RWEQ则充分考虑了土壤、气候、植被等因素，且影响因子数据易获取，与GIS系统具有较强的亲和力，因此被广泛应用到我国防风固沙功能的时空变化评估中^［10］。在驱动因素方面，地理探测器可利用空间数据的分层异质性揭示其背后的驱动机制，并被成功应用于各类生态系统功能的驱动因素分析^［11］。

尽管已有诸多研究对毛乌素沙地的生态环境变化及防风固沙功能进行了监测与评估，如李艳娇^［12］应用RWEQ模型分析了毛乌素沙地东南部的防风固沙功能对土地利用变化的响应。张秀娟等^［13］对包含毛乌素沙地荒漠草原在内的荒漠草原在禁牧前后草地生态系统的防风固沙功能进行了估算。马苏等^［14］对位于毛乌素沙地东北部的鄂尔多斯市防风固沙功能时空变化及驱动因素进行了分析。上述研究或集中于人为管理措施对防风固沙服务的影响，或集中于某一生态系统类型或行政区域的防风固沙功能变化及其影响因素分析，缺少因生态工程实施导致的土地利用变化对毛乌素沙地防风固沙功能空间分异和驱动因素的量化分析。因此，本文基于毛乌素沙地遥感影像数据，以毛乌素沙地为研究区，采用RWEQ与GIS空间统计技术相结合的研究方法，评估分析2000—2023年毛乌素沙地防风固沙功能的时空变化及其区域差异，并基于地理探测器模型，定量研究11个自然因素对毛乌素沙地防风固沙功能的单独影响及交互作用，旨在探究生态工程对防风固沙功能的影响，并为后期工程分区施策治理提供参考依据。

1 研究区概况

毛乌素沙地是中国四大沙地之一，面积约4.22万km²。它既是北方沙尘暴的主要沙源地，也是国家确定的以防风固沙为核心的重点生态功能区。行政区划上涵盖了陕西、宁夏和内蒙古自治区。地理坐标为37°27.5′—39°22.5′N、107°20′—111°30′E，位于我国北方农牧交错带的中部，地处鄂尔多斯高原与黄土高原的过渡地带，属于荒漠草原—干草原—森林草原的过渡带，具有明显的生态脆弱性和空间异质性。研究区气候特征为典型的大陆季风气候，兼具干旱和半干旱气候特点^［15］。年均气温6.5~10.3 ℃，年均降水量为150~500 mm，干燥度在空间分布上呈自东南向西北递减趋势^［16］。冬春季节盛行西北风，夏秋季节则以东南风为主，起沙风主要出现在冬末、春季，年大风时间20~40 d。毛乌素沙地海拔高度950~1 600 m，整体呈自西北向东南倾斜趋势^［17］，根据营力作用可划分为风成地貌、湖成地貌、流水地貌、干燥地貌、黄土地貌和其他地貌^［18］。沙地内部呈现流动沙丘、半固定沙丘和固定沙丘镶嵌分布。毛乌素地带性土壤以棕钙土、栗钙土和黑垆土为主，非地带性土壤则以盐碱土、草甸土及风沙土为主。毛乌素沙地从西北向东南依次由荒漠向草原再到森林过渡，植被覆盖度呈现增加的趋势。地带性草原植被以大针茅（Stipa grandis）、戈壁针茅（Stipa tianschanica）、百里香（Thymus mongolicus）等为主；地带性荒漠植被以藏锦鸡儿（Caragana spinifera）、红砂（Reaumuria Pall. Maxim.）、柠条（Caragana korshinskii Kom.）等为主；参照1∶100万中国植被图集可将当地植被类型分为9类（图1），分别为草丛、栽培植被、沼泽、灌丛、阔叶林、荒漠草原、草甸、沙漠和其他。其中，栽培植被多以植被恢复和防风固沙建植的人工林，包括长柄扁桃 （Amygdalus pedunculata），紫穗槐 （Amorpha fruticosa L.），樟子松 （Pinus sylvestris var. mongholica Litv.），油蒿 （Artemisia ordosica），沙柳 （Salix psammophila）和柠条等^［15］。

2 研究方法与数据来源

2.1 数据来源及处理

本研究的时间范围为2000—2023年共计24 a。鉴于本研究使用RWEQ模型估算研究区防风固沙功能，该模型涉及主要参数包括气候、土壤、植被和地形等因素。所需基础数据主要包括2000—2023年的土地利用/覆盖、DEM、NDVI、气象、土壤类型数据。上述数据来源及内容如表1所示。其中，气象数据首先利用Matlab检验去除缺失严重的站点，对部分站点进行插补后以DEM为协变量，利用ANUSPLIN软件进行插值；大风天数利用MATLAB软件统计研究时段内毛乌素沙地及周围的10个气象站≥临界起沙风速（5 m/s）^［19］的各等级风速累计时间，10个气象站分别为鄂托克旗、鄂托克前旗、伊金霍洛旗、乌审旗、榆林市、神木市、盐池县、定边县、靖边县和横山区。粗砂、粉砂、黏粒、碳酸钙含量和有机质等土壤属性数据在Arcmap10.2中使用栅格计算器进行相应计算；雪深数据采用Arcpy进行批量格式转换后，用Matlab计算雪盖因子；NDVI数据集已经过最大值合成处理，有效去除大气云层、建筑阴影和太阳高度角的影响，本研究以像元二分法计算植被覆盖度应用于植被因子计算；ASTER GDEM数据经过镶嵌裁剪后用来提取地表粗糙度因子；土地利用数据经裁剪后将二级分类合并为一级类；在模型运算过程中，本研究所用栅格数据均采用WGS_1984_UTM_Zone_48 N投影，并统一使用NEAREST重采样为1 km空间分辨率。

2.2 研究方法

2.2.1 RWEQ模型

风蚀模型是评估土壤风蚀状况的主要技术手段。防风固沙功能一般用防风固沙量来体现，可通过风蚀模型估算地表无植被状况下的土壤风蚀量（潜在风蚀量，用SL_r表示）与植被覆盖条件下土壤风蚀量（实际风蚀量，用SL表示）的差值来实现。参照我国生态红线划定指南及目前风蚀模型应用的广泛程度^［20］，本研究采用RWEQ模型进行土壤风蚀量的估算。计算公式如下：

\mathrm{S}{\mathrm{L}}_r=\frac{\mathrm{2}z}{s_r^{\mathrm{2}}}{Q}_{r\mathrm{m}\mathrm{a}\mathrm{x}}\cdot {\mathrm{e}}^{-{\left(z/s_r\right)}^{\mathrm{2}}}（1）

Q_{r\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{109.8}\times 〔\mathrm{W}\mathrm{F}\times \mathrm{E}\mathrm{F}\times \mathrm{S}\mathrm{C}\mathrm{F}\times K\mathrm{\text{'}}〕（2）

S_r=\mathrm{150.71}\times {\left(\mathrm{W}\mathrm{F}\times \mathrm{E}\mathrm{F}\times \mathrm{S}\mathrm{C}\mathrm{F}\times K\mathrm{\text{'}}\right)}^{-\mathrm{0.3711}}（3）

\mathrm{S}\mathrm{L}=\frac{\mathrm{2}z}{S^{\mathrm{2}}}{Q}_{\mathrm{m}\mathrm{a}\mathrm{x}}\cdot {\mathrm{e}}^{-{\left(z/s\right)}^{\mathrm{2}}}（4）

Q_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{109.8}\times 〔\mathrm{W}\mathrm{F}\times \mathrm{E}\mathrm{F}\times \mathrm{S}\mathrm{C}\mathrm{F}\times K\mathrm{\text{'}}\times C〕（5）

S=\mathrm{150.71}\times {\left(\mathrm{W}\mathrm{F}\times \mathrm{E}\mathrm{F}\times \mathrm{S}\mathrm{C}\mathrm{F}\times K\mathrm{\text{'}}\times C\right)}^{-\mathrm{0.3711}}（6）

式中：\mathrm{S}{\mathrm{L}}_r为潜在风蚀量（kg/m²）； Q_rmax为潜在风力的最大输沙能力（kg/m）； S_r为潜在关键地块长度（m）；SL为土壤侵蚀模数（kg/m²）； Q_max为风力的最大输沙能力（kg/m）；S为关键地块长度（m）；SR为防风固沙量（kg/m²）；Z为下风向距离（m），取50 m； WF、EF、SCF、K′、C分别为影响土壤风蚀状况的气象因子、土壤可蚀性因子、土壤结皮因子、地表粗糙度因子和植被因子。

（1） WF：气象因子，它表征了综合气象因素下风力对土壤颗粒的搬运能力；

\mathrm{W}\mathrm{F}=\mathrm{W}\mathrm{E}\times \frac{\rho }{g}\times \mathrm{S}\mathrm{W}\times \mathrm{S}\mathrm{D}（7）

\mathrm{W}\mathrm{E}=u_{\mathrm{2}}\times {\left(u_{\mathrm{2}}-u_{\mathrm{1}}\right)}^{\mathrm{2}}\times N_{\mathrm{d}}（8）

\mathrm{S}\mathrm{W}=\frac{\mathrm{E}{\mathrm{T}}_{\mathrm{p}}-\left(R+I\right)\left(R_{\mathrm{d}}/N\right)}{\mathrm{E}{\mathrm{T}}_{\mathrm{p}}}（9）

\mathrm{E}\mathrm{F}=\frac{\mathrm{29.09}+\mathrm{0.31}\mathrm{s}\mathrm{a}+\mathrm{0.17}\mathrm{s}\mathrm{i}+\mathrm{0.33}\left(\frac{\mathrm{s}\mathrm{a}}{\mathrm{c}\mathrm{l}}\right)-\mathrm{2.59}\mathrm{O}\mathrm{M}-\mathrm{0.95}\mathrm{C}\mathrm{C}}{\mathrm{100}}

\rho =\mathrm{348.01}\left(\frac{\mathrm{1.013}-\mathrm{0.1183}H+\mathrm{0.0048}{H}^{\mathrm{2}}}{T}\right)（10）

式中：WF表示气象因子（kg/m）；WE表示风场强度因子；ρ表示空气密度（kg/m³）；g表示重力加速度，取值9.8 m/s²；H为海拔高度（km）；T为绝对温度（开氏度）；SW表示土壤湿度因子；SD表示雪盖因子（无积雪覆盖天数/研究总天数），定义雪盖深度<25.4 mm为无积雪覆盖；u₂表示近地面2 m风速（m/s）；u₁表示临界起沙风速（m/s）；ET_p为潜在蒸散量（mm）；R为观测时段内的降水量（mm）；I为观测时段内的灌溉量（mm），（本文取0）；R_d为观测时段内的降水次数和灌溉天数；N为观测时段天数。

本研究整理下载的气象站点风速数据是在距地面高度为10 m处观测记录所得，但本研究RWEQ模型输入的风速数据需距地面高度为 2 m^［21］，故采用以下公式进行转换：

u_{\mathrm{2}}=u_{\mathrm{10}}{\left(\frac{z_{\mathrm{2}}}{z_{\mathrm{10}}}\right)}^{\frac{\mathrm{1}}{\mathrm{7}}}（11）

式中：u_{\mathrm{2}}为2 m高处风速（m/s）； u_{\mathrm{10}}为10 m高度风速（m/s）；z_{\mathrm{2}}和z_{\mathrm{10}}分别为距离地面2 m和10 m高处的粗糙度（无量纲）。

（2）EF：土壤可蚀性因子，它表征了一定土壤理化条件下土壤受风蚀的程度；

（12）

式中：EF表示土壤可蚀性因子；sa，si 和cl分别表示土壤粗砂、粉砂和黏粒含量（%）；OM表示土壤有机质含量（g/kg）；CC表示土壤中碳酸钙含量（%）。其中，土壤质地数据需采用对数正态分布模型将土壤粒径转化为RWEQ模型输入需要的美国标准。

（3） SCF：土壤结皮因子，它表征了一定土壤结皮盖度下土壤受风蚀的程度；

\mathrm{S}\mathrm{C}\mathrm{F}=\frac{\mathrm{1}}{\mathrm{1}+\mathrm{0.006}{\left(\mathrm{c}\mathrm{l}\right)}^{\mathrm{2}}+\mathrm{0.021}{\left(\mathrm{O}\mathrm{M}\right)}^{\mathrm{2}}}（13）

式中：SCF表示土壤结皮因子；cl表示土壤黏粒含量（%）；OM表示土壤有机质含量（g/kg）。

（4） K′：地表粗糙度因子，它用以模拟地表粗糙度对风蚀的影响；

K\mathrm{\text{'}}=\mathrm{c}\mathrm{o}\mathrm{s}\alpha（14）

式中：K′表示地表粗糙度因子（无量纲）；α为地形坡度（°），利用 ArcGIS软件中坡度（slope）模块计算获得。

（5） C：植被因子，它表征一定植被盖度下对土壤风蚀的抑制程度；

C={\mathrm{e}}^{-\mathrm{0}.{\mathrm{0483}}_{\mathrm{S}\mathrm{C}}}（15）

\mathrm{S}\mathrm{C}=\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}（16）

式中：\mathrm{S}\mathrm{C}为植被覆盖度（%）；\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}、NDVI_max、NDVI_min分别为NDVI实际值、最大值和最小值。

为分析毛乌素沙地实施生态工程前后防风固沙功能的时空分异，本研究参考中国北方固沙量划分等级，将研究区单位面积固沙量划分为5个等级（表2）。

2.2.2 趋势分析

基于ArcGIS软件利用空间栅格数据，以像元为单元，采用一元线性回归法计算防风固沙功能的变化趋势和累积变化幅度。计算公式为：

S=\frac{n\times {\displaystyle \sum _{i=\mathrm{1}}^n}i\times v_i-{\displaystyle \sum _{i=\mathrm{1}}^n}i{\displaystyle \sum _{i=\mathrm{1}}^n}{v}_i}{n\times {{\displaystyle \sum _{i=\mathrm{1}}^n}{i}^{\mathrm{2}}-〔{\displaystyle \sum _{i=\mathrm{1}}^n}i〕}^{\mathrm{2}}}（17）

式中：n代表年份，取值范围为［1，24］。n为研究时段长度，n=24；v_i为第i年的年均单位面积防风固沙量；S为各像元年均单位面积防风固沙量的变化趋势斜率。该斜率值可以反映出防风固沙功能在不同时期和不同区域的变化速度和趋势。S>0，表明防风固沙能力呈增加趋势；S<0，表明防风固沙能力呈下降趋势。使用F检验进行显著性检验。将线性回归与F检验的结果相结合，分为极显著减少（S<0，p≤0.01）、显著减少（S<0，0.01<p≤0.05）、不显著降低（S<0，p>0.05）、极显著增加（S≥0，p≤0.01）、显著增加（S≥0，0.01<p≤0.05）、不显著增加（S≥0，p>0.05）。

2.2.3 防风固沙功能驱动因素分析

地理探测器可基于统计学原理，对研究对象的空间分异进行探测和驱动力分析^［22］。本研究主要通过地理探测器对影响防风固沙功能变化的多个解释变量进行解释力分析，厘定解释各变量在交互作用下的解释力，并揭示防风固沙功能变化对驱动因子响应的区域差异性。本文所用到的地理探测器为分异及因子探测和交互探测器。

分异及因子探测：可以解释某X因子在多大程度上可以解释因变量Y的空间分异，结果用q值表示，q值越大表示X因子对因变量Y的解释力越强，反之就越弱。

q=\mathrm{1}-{\displaystyle \sum _H^L}{N}_h{\sigma }_h^{\mathrm{2}}/N{\sigma }^{\mathrm{2}}（18）

式中：ℎ=1，…，L为变量Y或因子X的分层，即分类或分区；N_h和N分别为层ℎ和全区的单元数；{\sigma }_h^{\mathrm{2}}和{\sigma }^{\mathrm{2}}分别是层ℎ和全区Y值的方差。

交互探测器：用以计算任意两个变量的交互效应。交互作用如表3所示。

本研究选取影响防风固沙功能的气候（年均降水量X₁，年均温度X₂，年均风速X₃，年均蒸散发X₄）、地形（高程X₅，坡度X₆，地貌类型X₇）、土壤（土壤类型X₈）、植被（NDVI X₉和植被类型X₁₀）和土地利用类型X₁₁共5类自然因素11个因子作为探测因子，分析毛乌素沙地2000年、2005年、2010年、2015年、2020年和2023年年均防风固沙功能物质量变化的驱动因素。

3 结果与分析

3.1 防风固沙功能时空变化特征

2000—2023年毛乌素沙地土壤风蚀模数呈下降趋势，防风固沙量则整体呈上升趋势但存在明显的年际波动（图2）。

从时间尺度上看，2000—2023年毛乌素沙地土壤风蚀模数呈降低趋势，变幅为241.64~409.26 t/（hm² · a），平均值为347.12 t/（hm² · a），其年际变化呈减少趋势，平均变化量为-2.09 t/（hm² · a）。说明研究区整体土壤风蚀活动减弱，生态工程建设成效显著。单位面积防风固沙量整体呈波动上升趋势，平均增速为0.045 t/（hm² · a）。具体表现为从2000年的2.18 t/hm²，先波动上升至2019年的4.68 t/hm²，后于2021年急速下降到3.35 t/hm²。多年平均单位面积防风固沙量约为3.82 t/（hm² · a）。2000—2023年毛乌素沙地年均防风固沙量为3.32×10⁷ t/a，2019年毛乌素沙地防风固沙量最高，为4.63×10⁷ t，与2000年相比，2023年防风固沙总量增加了2.06×10⁷ t，并以0.86×10⁶ t/a的多年平均增速增加。

通过与以往的研究结果进行对比验证以判断RWEQ模型计算结果的可靠性。本文计算的多年平均土壤风蚀模数为3.82 t/（hm² · a），高于京津风沙源治理工程区的2.83 t/（hm² · a），但与在中国北方^［23-25］、黄土高原北部风蚀区^［10］的计算结果范围基本保持一致。

参考中国北方固沙量划分等级，近24年毛乌素沙地单位面积防风固沙量表现出明显的等级差异性和空间异质性（图3）。空间上，毛乌素沙地单位面积防风固沙量整体上呈东南高西北低的趋势，且其高值区范围随时间推移有明显扩大趋势（图4），即 2023年明显高于2010年和2000年，且逐渐从东南向西北扩张。具体表现为单位面积防风固沙量高值区〔>10 t/（hm² · a）〕主要分布于横山区、靖边县、神木市中部、榆阳区东南部、定边县东部和伊金霍洛旗的中部地区；低值区〔0~2 t/（hm² · a）〕主要分布于乌审旗的北部、鄂托克旗的东南部与西北角。

3.2 土地利用类型与防风固沙功能的关系

基于土地利用类型和不同植被类型的防风固沙能力差异进行分时段计算，结果表明（表4），相较于2000年，2023年毛乌素沙地建设用地、林地面积分别增加326.62%，23.51%，而草地、耕地、水域和未利用地面积分别减少1.28%，4.68%，39.47%和2.12%，林地的增加主要来源于生态工程的实施。不同土地利用类型在控制土壤风蚀、发挥防风固沙作用方面存在显著差异，从表4可以看出，2000年各土地利用类型单位面积固沙量序列依次为耕地〔3.26 t/（hm² · a）〕>草地〔2.33 t/（hm² · a）〕>林地〔2.09 t/（hm² · a）〕>建设用地〔1.66 t/（hm² · a）〕>水域〔1.58 t/（hm² · a）〕>未利用地〔1.11 t/（hm² · a）〕；2023年依次为耕地〔6.48 t/（hm² · a）〕>林地〔4.99 t/（hm² · a）〕>草地〔4.55 t/（hm² · a）〕>水域〔1.95 t/（hm² · a）〕>建设用地〔1.72 t/（hm² · a）〕>未利用地〔1.33 t/（hm² · a）〕。相较于2000年，2023年林地、耕地和草地单位面积固沙量变化较大，分别上升138.12%，98.56%，95.40%。草地在2000年和2023年分别贡献了61.82%和64.21%的固沙量，说明草地是起主要防风固沙作用的土地利用类型。

3.3 防风固沙功能驱动因素分析

3.3.1 单因子分析

在单因子探测计算结果当中，按q值大小排序，对影响防风固沙功能的气候、地形、土壤、植被和土地利用类型5类自然因素11个因子进行排序。如表5所示，年均防风固沙量空间变化的影响程度由大到小依次为NDVI（30.30%）>降水（17.90%）>土壤类型（15.90%）>地貌类型（15.90%）>植被类型（10.90%）>土地利用（10.10%）>高程（8.40%）>温度（7.60%）>风速（6.60%）>潜在蒸散发（6.10%）>坡度（4.90%）。选取2000年、2005年、2010年、2015年、2020年、2023年6个典型年份，基于地理探测器的探测因子q值变化探究毛乌素沙地防风固沙功能的驱动因素，综合2000年、2005年、2010年、2015年、2020年和2023年的探测结果，NDVI、降水、土壤类型、地貌类型、植被类型和土地利用是影响防风固沙功能空间变化的主要驱动因素，且其解释力均稳定在10%以上。整体看来，不同驱动因子在不同时期对于毛乌素沙地荒漠化的影响存在较大差异，且由于自然环境与社会经济条件不同，各影响因子的解释力存在较大的年际变化。

3.3.2 多因素交互作用分析

鉴于单一因子无法完全解释防风固沙功能的空间分异性，需考虑多因子之间的协同和交互作用。本研究采用地理探测器中的交互探测器可用以确定各因子之间的交互作用，评价两两因子之间的交互作用。2000年、2005年、2010年、2015年、2020年和2023年交互探测器分析结果显示，各因子之间对毛乌素沙地防风固沙功能的影响均存在增强型交互作用（图5）。除蒸散发、高程、风速、温度、坡度之间大部分为非线性增强关系外，其他因子之间普遍存在双因子的增强关系。相较于单因子分析，各因子间的交互作用均强于单因子q值的影响。2000年NDVI与其他各因子的交互作用均明显高于其他双因子交互作用对防风固沙功能空间分异的解释能力，其中解释能力最高的为NDVI∩土壤类型（45.10%）和NDVI∩风速（44.30%）；2005年解释能力最高的为NDVI∩土壤类型（35.60%）和NDVI∩地貌类型（35.40%）；2010年解释能力最高的为NDVI∩土地利用（30.70%）和NDVI∩土壤类型（30.20%）；2015年解释能力最高的为NDVI∩风速（41.20%）和NDVI∩高程（40.10%）；2020年解释能力最高的为NDVI∩高程（47.00%）和NDVI∩地貌类型（46.80%）；2023年解释能力最高的为风速∩地貌类型（24.30%）和风速∩降水（22.50%）。整体上来看，NDVI、降水、风速与其他因子组合产生的交互作用均明显强于其他双因子交互作用对防风固沙功能空间分异的解释能力。

4 讨 论

4.1 毛乌素沙地防风固沙功能的时空分异

2000—2023年毛乌素沙地防风固沙功能呈波动增长趋势，具体表现为土壤风蚀模数呈明显下降趋势，而防风固沙量呈增加趋势。这与王旭洋等^［26］对北方防沙带防风固沙能力明显增强的结论一致。本研究参照同研究区已有成果和同类型区的侵蚀模数来验证土壤风蚀模型（RWEQ）计算结果。本研究结果数值上与中国北方、黄土高原北部风蚀区^{［10，23-25］}的估算结果接近，说明RWEQ模型在毛乌素沙地防风固沙功能评估中具有较好的适用性。从时间尺度来看，侵蚀模数的变化趋势也与同区域的研究结果一致，均呈现出整体波动中增加的变化趋势。在空间分布上，表现为单位面积防风固沙量呈减少区域少于增加区域，且其东南部高于西北部，表现出明显的空间分异性。低值区的产生可能与人工造林密度、物种选择超出了当地土壤水分阈值，导致出现土壤干层、植物衰败死亡，从而使防风固沙功能降低有关。而高值区的产生与发展可能与近年来北方风蚀区平均风速下降有关，也可能与研究区近年来因实施一系列生态工程使整体NDVI有所提升有关。近40年来，毛乌素沙地经过地方和国家工程和项目的实施，使NDVI从1985年的45.1%上升到2020年的73.7%，沙漠化土地面积减少1 963.75 km^2［4］。尤其是2000年后的封山禁牧^［27］、退耕还林还草^［28］、绿洲化建设和自然保护区建设^［4］等措施均能在一定程度上逆转沙漠化进程^［29］。此外，由于退耕还林引入的灌丛化不仅影响草原生态系统植被组成，还通过灌丛斑块的“沃岛效应”加剧了草地退化和土壤风蚀^［30］，从而使草地防风固沙功能低于林地（表4）。

4.2 毛乌素沙地防风固沙功能驱动因素

基于防风固沙功能相关影响因子以及数据的可获取性，本文从驱动因子的单独作用和交互作用入手，分别分析了自然因素对防风固沙功能变化和空间分异的解释力。结果表明，不同自然因素对防风固沙功能的解释力存在差异性。单因子分析结果表明，气候因子如降水、风速等对研究区防风固沙功能的影响普遍较强，这与王晓峰等对黄河流域防风固沙功能驱动因子的分析结论一致。交互作用分析结果表明，除气候因素外，NDVI、土壤类型和DEM对防风固沙功能具有较大影响，可通过各类影响因素之间的交互作用使防风固沙功能增强。如植被在可增强土壤抗蚀性的同时还会影响区域环境如温度、土壤湿度、蒸散发等。土地利用作为具象化的人类活动，其类型转化会使土壤结构、功能、粗糙度、NDVI等发生显著变化，导致固沙功能存在差异。本研究中，城市化、农业活动与生态工程作用的土地利用变化是影响植被空间格局的关键因素。毛乌素沙地近年来实施的大型生态工程如退耕还林还草、天然林保护工程、三北防护林建设等会通过土地利用类型转化影响研究区的生态系统结构、过程，使林草面积不断扩大，NDVI呈显著上升趋势，防风固沙功能也趋于增强^［25］。增加植被覆盖是应用最广泛的生物治沙措施，它不仅可以提高土壤抗蚀性，还可通过植被覆盖地表减小风力对地表的直接吹蚀作用^［31］。NDVI和土地利用类型都会影响土壤风蚀，这与王旭洋等^［26］对中国北方农牧交错带土壤风蚀影响因子的分析结论一致。本研究中，草地是发挥防风固沙作用的主要土地利用类型，这与王金凤等^［25］对黄土高原北部风蚀区防风固沙功能驱动因素分析结论一致。但毛乌素沙地部分区域因人工灌丛化的加剧，使部分草地防风固沙功能弱于林地。气候变化和人类活动共同解释了毛乌素沙地植被变化，但由于毛乌素沙地存在相对显著的气候、植被、地貌和土壤属性的空间分异，自然要素组合和配置使防风固沙功能空间分异更为显著，如位于毛乌素沙地东南缘的NDVI与降水的交互作用就显著强于毛乌素沙地西北部。在降水相对充足的东南部，高温会促进NDVI增大，进而使防风固沙功能增强。而毛乌素沙地西北部则因降水稀少且植被稀疏，高温会增加干旱发生强度和频率，限制植被生长，使得地表更加干燥，容易遭受风蚀^［32］。

4.3 不确定性

本文的不确定性主要体现在4个方面：（1） 基础数据层面，风速是RWEQ模型中最重要的气象驱动力数据。实际的风速时刻在变化，模型中代入不同时段的平均风速会产生不同的模拟结果。本研究为避免单个年份的气象数据带来的随机误差，以多年研究为主。且本研究中风速数据采用的方法是结合气象站点风速数据与插值软件ANUSPLIN得到风速的栅格数据，因此，在风速插值的结果中存在不确定性。（2） 参数选择方面，本研究在筛选影响防风固沙功能变化的驱动因素时，考虑了气候变化和人类活动两大方面的因素，但仍存在一些不足之处。如本研究仅聚焦于年均气温、年降水量和年均风速三大关键因子，没有考虑更广泛的社会经济因素（包括GDP、人口变化），而这些因素对防风固沙功能的影响同样不容忽视。（3） 驱动因素方面，鉴于数据的可获得性，当前模型主要考虑了气象、地形、土壤和植被等自然因素，但对除土地利用动态之外的人类活动（如放牧、农业生产方式调整、生态保护活动等）影响强度缺乏定量表征。（4） 模型适用性层面，RWEQ模型是基于平原的统计值建立的经验模型，具体应用时需要修正。一些重要参数的计算如结皮因子、粗糙度因子均为基于大量野外试验的一种经验模型，这可能导致模型预测结果与研究区实际情况之间存在一定的偏差。为了进一步验证和提升RWEQ模型的准确性，未来需结合野外观测数据，通过提高气候、土壤、植被覆盖等参数的分辨率，提升模拟精度，并结合实地风蚀监测数据或其他经过验证的风蚀模型，以期为RWEQ模型的持续优化与广泛应用提供理论支持。

5 结 论

本文基于遥感数据，以毛乌素沙地为研究区，采用RWEQ与GIS空间统计技术相结合的研究方法，评估分析2000—2023年毛乌素沙地防风固沙功能的时空变化及其区域差异，并基于地理探测器模型，定量研究了11个自然因素对毛乌素沙地防风固沙功能的单独影响及交互作用，主要结论如下：

（1） 从时间尺度上看，2000—2023年毛乌素沙地土壤风蚀模数呈明显下降趋势，单位面积防风固沙量和防风固沙总量则呈波动上升趋势。从空间尺度上看，近24年毛乌素沙地单位面积防风固沙量表现出明显的等级差异性和空间异质性。毛乌素沙地防风固沙量整体上呈现东南高西北低的趋势。

（2） 不同土地利用类型在控制土壤风蚀、发挥防风固沙作用方面存在显著差异，其中，草地是毛乌素沙地起主要防风固沙作用的土地利用类型。

（3） 驱动因素方面，毛乌素沙地防风固沙功能在2000—2023年呈现显著提升趋势，主要受NDVI和降水的影响。单因子分析结果表明，NDVI对防风固沙功能的空间变化影响最大，其次为降水。双因子交互作用分析结果表明，NDVI与其他因子（如土壤类型、风速等）的交互效应尤为突出，并在多个年份中表现出最高解释能力。

参考文献

[1] 王涛，朱震达.我国沙漠化研究的若干问题：1.沙漠化的概念及其内涵［J］.中国沙漠，2003，23（3）：209-214.

[2] Wang T， Zhu Z D. Study on sandy desertification in China：1. definition of sandy desertification and its connotation［J］. Journal of Desert Research， 2003，23（3）：209-214.

[3] 卢琦，肖春蕾，包英爽，等.打赢“三北”攻坚战，再造一个“新三北”：实现路径与战略规划［J］.中国科学院院刊，2023，38（7）：956-965.

[4] Lu Q， Xiao C L， Bao Y S， et al. Implementation path and strategic planning of winning the battle of “Three-North” and reconstructing “New Three-North”［J］. Bulletin of Chinese Academy of Sciences， 2023，38（7）：956-965.

[5] 原媛，母艳梅，邓钰洁，等.植被覆盖度和物候变化对典型黑沙蒿灌丛生态系统总初级生产力的影响［J］.植物生态学报，2022，46（2）：162-175.

[6] Yuan Y， Mu Y M， Deng Y J， et al. Effects of land cover and phenology changes on the gross primary productivity in an Artemisia ordosica shrubland［J］. Chinese Journal of Plant Ecology， 2022，46（2）：162-175.

[7] 杨皓钦，王海兵，左合君，等.毛乌素沙地绿洲化土地变化模式及稳定性［J］.应用生态学报，2024，35（3）：687-694.

[8] Yang H Q， Wang H B， Zuo H J， et al. Change pattern and stability of oasisization land in Mu Us Sandy Land［J］. Chinese Journal of Applied Ecology， 2024，35（3）：687-694.

[9] 邹学勇，张春来，程宏，等.土壤风蚀模型中的影响因子分类与表达［J］.地球科学进展，2014，29（8）： 875-889.

[10] Zou X Y， Zhang C L， Cheng H， et al. Classification and representation of factors affecting soil wind erosion in a model［J］. Advances in Earth Science， 2014，29（8）： 875-889.

[11] 巩国丽，刘纪远，邵全琴.草地覆盖度变化对生态系统防风固沙服务的影响分析：以内蒙古典型草原区为例［J］.地球信息科学学报，2014，16（3）：426-434.

[12] Gong G L， Liu J Y， Shao Q Q. Effects of vegetation coverage change on soil conservation service of typical steppe in Inner Mongolia［J］. Journal of Geo-Information Science， 2014，16（3）：426-434.

[13] 薛永盛，杨雪梅.北方防沙带典型县域主要生态系统服务变化及权衡协同关系［J］.草业科学，2022，39（7）：1363-1374.

[14] Xue Y S， Yang X M. Changes in major ecosystem services and their trade-offs and synergy in typical counties of the northern desert belt［J］. Pratacultural Science， 2022，39（7）：1363-1374.

[15] 常文静，丛士翔，王融融，等.气候变化和人类活动对毛乌素沙地NDVI变化的量化分析［J］.干旱区地理，2025，48（1）：63-74.

[16] Chang W J， Cong S X， Wang R R， et al. Quantitative analysis of NDVI changes in Mu Us Sandy Land by climate change and human activities［J］. Arid Land Geography， 2025，48（1）：63-74.

[17] 王彬，韩庆杰，刘冰，等.沙面温度对风蚀动力过程的影响［J］.中国沙漠，2023，43（6）：10-19.

[18] Wang B， Han Q J， Liu B， et al. The influence of sand surface temperature on wind erosion dynamic process［J］. Journal of Desert Research， 2023，43（6）：10-19.

[19] Gong G L， Liu J Y， Shao Q Q， et al. Sand-fixing function under the change of vegetation coverage in a wind erosion area in northern China［J］. Journal of Resources and Ecology， 2014，5（2）：105-114.

[20] Gao J B， Zuo L Y， Liu W L. Environmental determinants impacting the spatial heterogeneity of karst ecosystem services in Southwest China［J］. Land Degradation & Development， 2021，32（4）：1718-1731.

[21] 李艳娇.陕北地区土地利用变化对防风固沙效益的影响［J］.西部大开发（土地开发工程研究），2020，5（3）：56-63，70.

[22] Li Y J. The effect of land use change on windbreak and sand fixation in northern Shaanxi［J］. Land Development and Engineering Research， 2020，5（3）：56-63，70.

[23] 张秀娟，马珂昕，刘晓斌，等.实施禁牧前后宁夏草地防风固沙的时空变化特征［J］.草业科学，2023，40（3）：638-653.

[24] Zhang X J， Ma K X， Liu X B， et al. Spatiotemporal patterns of windbreak and sand fixation before and after grazing prohibition in Ningxia［J］. Pratacultural Science， 2023，40（3）：638-653.

[25] 马苏，刘军会，康玉麟，等.鄂尔多斯市防风固沙功能时空变化及驱动因素分析［J］.环境科学研究，2022，35（11）：2477-2485.

[26] Ma S， Liu J H， Kang Y L， et al. Spatio-temporal changes of sand-fixing function and its driving factors in the Ordos［J］. Research of Environmental Sciences， 2022，35（11）：2477-2485.

[27] Zhu Z C， Shao M A， Jia X X， et al. Rainfall partitioning characteristics and simulation of typical shelter forest in Chinese Mu Us Sandy Land［J］. Science of the Total Environment， 2024，945：174091.

[28] Feng K， Wang T， Liu S L， et al. Path analysis model to identify and analyse the causes of aeolian desertification in Mu Us Sandy Land， China［J］. Ecological Indicators， 2021，124：107386.

[29] 杨梅焕，靳小燕，王涛.毛乌素沙地植被物候变化及其对气候变化的响应［J］.水土保持通报，2022，42（2）：242-249.

[30] Yang M H， Jin X Y， Wang T. Vegetation phenology change of Mu Us sandy land and its response to climate change［J］. Bulletin of Soil and Water Conservation， 2022，42（2）：242-249.

[31] 温仰磊，蒿承智，谭利华，等.1∶25万《毛乌素风沙地貌图》的编制［J］.中国沙漠，2018，38（3）：508-515.

[32] Wen Y L， Hao C Z， Tan L H， et al. Compilation of geomorphic map of the Mu Us desert［J］. Journal of Desert Research， 2018，38（3）：508-515.

[33] Zhang H B， Peng J， Zhao C N， et al. Wind speed in spring dominated the decrease in wind erosion across the Horqin Sandy Land in northern China［J］. Ecological Indicators， 2021，127：107599.

[34] 王洋洋，肖玉，谢高地，等.基于RWEQ的宁夏草地防风固沙服务评估［J］.资源科学，2019，41（5）：980-991.

[35] Wang Y Y， Xiao Y， Xie G D， et al. Sand-fixing function of the grassland ecosystem in Ningxia based on the revised wind erosion model［J］. Resources Science， 2019，41（5）：980-991.

[36] 闫峰，吴波，王艳姣.2000—2011年毛乌素沙地植被生长状况时空变化特征［J］.地理科学，2013，33（5）：602-608.

[37] Yan F， Wu B， Wang Y J. Spatial and temporal variations of vegetation growth status in mu us sandy land in 2000-2011［J］. Scientia Geographica Sinica， 2013，33（5）：602-608.

[38] 王劲峰，徐成东.地理探测器：原理与展望［J］.地理学报，2017，72（1）：116-134.

[39] Wang J F， Xu C D. Geodetector： Principle and prospective［J］. Acta Geographica Sinica， 2017，72（1）：116-134.

[40] 黄麟，祝萍，肖桐，等.近35年三北防护林体系建设工程的防风固沙效应［J］.地理科学，2018，38（4）：600-609.

[41] Huang L， Zhu P， Xiao T， et al. The sand fixation effects of three-north shelter forest program in recent 35 years［J］. Scientia Geographica Sinica， 2018，38（4）：600-609.

[42] 刘珺，郭中领，常春平，等.基于RWEQ和WEPS模型的中国北方农牧交错带潜在风蚀模拟［J］.中国沙漠，2021，41（2）：27-37.

[43] Liu J， Guo Z L， Chang C P， et al. Potential wind erosion simulation in the agro-pastoral ecotone of northern China using RWEQ and WEPS models［J］. Journal of Desert Research， 2021，41（2）：27-37.

[44] 王金凤，刘小玲，李庆，等.黄土高原北部风蚀区防风固沙服务时空分异及驱动因素［J］.中国沙漠，2023，43（4）：220-230.

[45] Wang J F， Liu X L， Li Q， et al. Spatio-temporal differentiation and driving factors of windbreak and sand fixation services in wind erosion area of the northern Loess Plateau［J］. Journal of Desert Research， 2023，43（4）：220-230.

[46] 王旭洋，郭中领，常春平，等.中国北方农牧交错带土壤风蚀时空分布［J］.中国沙漠，2020，40（1）：12-22.

[47] Wang X Y， Guo Z L， Chang C P， et al. Spatiotemporal distribution of soil wind erosion modulus in the agro-pastoral ecotone of North China［J］. Journal of Desert Research， 2020，40（1）：12-22.

[48] 程积民，井赵斌，金晶炜，等.黄土高原半干旱区退化草地恢复与利用过程研究［J］.中国科学：生命科学，2014，44（3）：267-279.

[49] Cheng J M， Jing Z B， Jin J W， et al. Restoration and utilization mechanism of degraded grassland in the semi-arid region of Loess Plateau［J］. Scientia Sinica （Vitae）， 2014，44（3）：267-279.

[50] Gao W D， Zheng C， Liu X H， et al. NDVI-based vegetation dynamics and their responses to climate change and human activities from 1982 to 2020： a case study in the Mu Us Sandy Land， China［J］. Ecological Indicators， 2022，137：108745.

[51] Liu M， Dries L， Heijman W， et al. The impact of ecological construction programs on grassland conservation in Inner Mongolia， China［J］. Land Degradation & Development， 2018，29（2）：326-336.

[52] Shen X J， Liu Y W， Liu B H， et al. Effect of shrub encroachment on land surface temperature in semi-arid areas of temperate regions of the Northern Hemisphere［J］. Agricultural and Forest Meteorology， 2022，320：108943.

[53] 高函，张宇清，吴斌，等.柠条林防风阻沙效益模拟研究初报［J］.北京林业大学学报，2010，32（4）：175-180.

[54] Gao H， Zhang Y Q， Wu B， et al. Simulating efficiency of wind-speed reduction and sand-break of Caragana korshinskii coppice［J］. Journal of Beijing Forestry University， 2010，32（4）：175-180.

[55] Li D J， Xu D Y， Wang Z Y， et al. The dynamics of sand-stabilization services in Inner Mongolia， China from 1981 to 2010 and its relationship with climate change and human activities［J］. Ecological Indicators， 2018，88：351-360.

基金资助

国家自然科学基金区域创新发展联合基金资助项目 “宁夏中部干旱带典型固沙植物高效抚育利用的水力调控机理”(U23A20223)

国家自然科学基金“白刺灌丛沙堆演化过程中冠层结构和根系构型对土壤水分变化的响应与适应”(42561011)