科尔沁沙地防风固沙服务时空变化

doi:10.7522/j.issn.1000-694X.2025.00083

科尔沁沙地防风固沙服务时空变化

张慧琳^,¹^,², 王卫国^,³, 薄一览¹, 晋子振¹

1.忻州师范学院地理系，山西忻州 034000

2.山西大学黄土高原研究所，山西太原 030006

3.山西师范大学地理科学学院，山西太原 030031

Spatiotemporal dynamics of wind erosion prevention and sand fixation service in the Horqin Sandy Land

Zhang Huilin^,¹^,², Wang Weiguo^,³, Bo Yilan¹, Jin Zizhen¹

1.Department of Geography，Xinzhou Teachers University，Xinzhou 034000，Shanxi，China

2.Institute of Loess Plateau，Shanxi University，Taiyuan 030006，China

3.College of Geographical Sciences，Shanxi Normal University，Taiyuan 030031，China

通讯作者: 王卫国（E-mail： wangwg@sxnu.edu.cn）

收稿日期: 2025-05-12 修回日期: 2025-06-16

基金资助:

山西省基础研究计划（自由探索类）青年项目.  202303021222182
山西省基础研究计划（自由探索类）青年项目.  202403021212316
山西省哲学社会科学规划课题.  2023YJ055

Received: 2025-05-12 Revised: 2025-06-16

作者简介 About authors

张慧琳（1993—），女，山西大同人，在读博士研究生，讲师，主要从事生态系统服务流变化E-mail：zhanghlmm@163.com , E-mail：zhanghlmm@163.com

摘要

科尔沁沙地土地沙化严重，研究其防风固沙服务的时空演变规律对区域生态安全与可持续发展具有重要意义。本研究利用修正风蚀方程（RWEQ）模型、XGBoost-SHAP、HYSPLIT（Hybrid Single-Particle Lagrangian Integrated Trajectory Model）等方法，系统分析了2000—2020年科尔沁沙地防风固沙服务时空演变、驱动机制。结果表明：（1）科尔沁沙地防风固沙服务呈现显著的时空异质性特征，单位面积固沙量呈先降后升态势，空间格局表现为中南高、北部低的分异特征。（2）XGBoost-SHAP模型显示，植被覆盖度是实际风蚀作用的核心驱动因子，而风速、降水与气温气象要素通过非线性作用影响风蚀过程。（3）HYSPLIT模拟揭示2000—2020年共产生3 559条风沙传输路径，轨迹数量整体下降，表明固沙能力提升。聚类分析识别出4条主导路径：31.67%向东入太平洋，29.92%向东北影响吉林省和黑龙江省，20.46%跨境扩散至朝鲜、日本等区域，17.95%南下抵达东南沿海。本文研究结果可为优化沙地治理策略、协调区域生态补偿机制提供科学依据，同时为干旱半干旱区生态系统服务流动研究提供方法参考。

关键词： 防风固沙服务 ; HYSPLIT ; XGBoost-SHAP ; 科尔沁沙地

Abstract

The Horqin Sandy Land suffers from severe land desertification. Research into the spatiotemporal evolution patterns of the wind erosion prevention and sand fixation service holds significant importance for regional ecological security and sustainable development. This study systematically analyzed the spatiotemporal dynamics， driving mechanisms of wind erosion prevention and sand fixation service in the Horqin Sandy Land from 2000 to 2020 using the Revised Wind Erosion Equation （RWEQ） model， the XGBoost-SHAP algorithm， and HYSPLIT trajectory simulations. The results indicate：（1） The wind erosion prevention and sand fixation service exhibited significant spatiotemporal heterogeneity. The sand fixation per unit area showed a trend of initial decrease followed by increase， with a spatial pattern characterized by a higher in the central-southern regions and lower in the northern areas. （2） The XGBoost-SHAP model quantified vegetation coverage as the core driving factor， while interactions among wind days， precipitation， and temperature nonlinearly influenced wind erosion processes. （3） HYSPLIT simulations identified 3 559 sand transport trajectories over the 20-year period， with an overall decline in trajectory frequency， indicating enhanced sand fixation capacity. Cluster analysis revealed four dominant pathways： 31.67% of trajectories moved eastward into the Pacific Ocean， 29.92% influenced northeastern regions （Jilin and Heilongjiang）， 20.46% crossed borders to areas such as North Korea and Japan， and 17.95% migrated southward to China’s southeastern coast. These findings provide a scientific basis for optimizing desertification control strategies and coordinating regional ecological compensation mechanisms， while offering methodological references for studying ecosystem service flows in arid and semi-arid regions.

Keywords： wind erosion prevention and sand fixation service ; HYSPLIT model ; XGBoost-SHAP model ; Horqin Sandy Land

PDF (9702KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

张慧琳, 王卫国, 薄一览, 晋子振. 科尔沁沙地防风固沙服务时空变化. 中国沙漠[J], 2025, 45(4): 227-240 doi:10.7522/j.issn.1000-694X.2025.00083

Zhang Huilin, Wang Weiguo, Bo Yilan, Jin Zizhen. Spatiotemporal dynamics of wind erosion prevention and sand fixation service in the Horqin Sandy Land. Journal of Desert Research[J], 2025, 45(4): 227-240 doi:10.7522/j.issn.1000-694X.2025.00083

0 引言

风蚀是全球陆地生态系统重要的环境胁迫因子，其发生与临界风导致的土壤颗粒运移密切相关^［1］，主要发生于干旱半干旱地区^［2］，影响范围广阔^［3］，覆盖全球陆地面积的三分之一^［4］。风蚀通过剥离表土、释放粉尘等，引发土地退化、碳库受损和生物多样性衰减等多重生态危机^［5-7］。特别是其衍生的沙尘暴现象，已成为威胁3.5亿人口健康和区域可持续发展的重大环境挑战^［8］。因此，及时开展风蚀评估和监测工作十分必要。在风蚀模拟研究领域，已有的研究模型包括Bangnold等^［9］提出的输沙方程、德克萨斯模型（Texas Erosion Analysis Model，TEAM）^［10］、风蚀方程（Wind Erosion Equation，WEQ）^［11］、修正风蚀方程（Revised Wind Erosion Equation，RWEQ）^［12］、风蚀预报系统（Wind Erosion Prediction System）^［13］等。RWEQ模型充分考虑了气候、植被、土壤等因素，被越来越多地应用到中国的风蚀和防风固沙研究中^［14-15］。李琪等^［16］利用RWEQ模型评估分析京津风沙源区防风固沙功能时空变化及其区域内部差异。柴茵超等^［17］利用RWEQ模型模拟毛乌素沙地实际风蚀量、潜在风蚀量、防风固沙量的时空变化特征，并借助地理因子探测器识别驱动因素。黄孟冬等^［18］利用RWEQ模型对浑善达克地区防风固沙时空变化进行分析。

防风固沙服务（Wind Erosion Prevention and Sand fixation Service， WEPS）指生态系统中各类自然植被对风沙抑制和固定作用，是干旱风蚀区生态系统的关键性防护型服务^［19］，有助于水土保持和降低风速，更能保护人体健康、农业生产和交通设施免受风沙危害^［20］，对提升环境质量和保障人类福祉具有关键作用。防风固沙服务供给区和受益区的空间解耦，决定了防风固沙服务在空间上的流动性，可以映射生态系统服务流传输的动态过程，识别生态系统服务流的流动路径、方向和距离，确定受益区和供给区^［21］。防风固沙服务的空间流动实际上是沙尘在大气与地表覆被界面之间的输移量（起沙量）以及沙尘在大气中的扩散过程，大气环流形势决定了沙尘输送距离及其影响范围^［22］。HYSPLIT（Hybrid Single-Particle Lagrangian Integrated Trajectory）模型由美国国家海洋与大气管理局（National Oceanic and Atmosphere Administration，NOAA）的空气资源实验室（Air Resource Laboratory，ARL）研发，用于计算和分析大气污染物输送、扩散轨迹的专业模式^［23］，可以模拟风沙轨迹与粉尘运输途径^［24］、模拟PM₁₀和PM_2.5扩散的时空变化和传播情况^［25-26］等。王国宇等^［19］利用RWEQ模型与HYSPLIT模型模拟无定河流域土壤风蚀量以及防风固沙服务流的时空变化。Yu等^［27］利用HYSPLIT模拟塔里木河流域防风固沙服务的空间流动。Gao等^［21］采用RWEQ风蚀模拟、HYSPLIT预测等多项模型，量化了锡林郭勒草原国家级自然保护区防风固沙服务存量并识别流向，提出了生态补偿方案。Xu等^［28］和肖玉等^［29］对中国防风固沙型重点生态功能区和宁夏回族自治区等区域的防风固沙量与防风固沙服务流的时空分布进行了深入的研究，对受益区进行识别。

科尔沁沙地是中国北方半干旱农牧交错带的典型区域、国家防风固沙型重点生态功能区^［30］，承载着遏制沙尘东侵、维护东北-华北生态屏障的重要使命。然而，持续的人口压力与高强度的人类活动（过度放牧、开垦和樵采等），导致该区顶级植被群落类型——疏林草原植被遭到严重破坏，诱发了不同类型和程度的土地沙漠化^［31-32］。但现有关于科尔沁沙地的研究中，定量解析该区WEPS服务时空格局、服务传输路径等方面仍存在明显不足，这严重制约着区域生态补偿机制构建与沙尘暴联防联控策略制定。因此本文基于RWEQ模型计算防风固沙服务供给的时空格局，并且运用机器学习解析风蚀的主要驱动因素，通过HYSPLIT模型模拟沙尘传输路径，确定受益区范围，研究成果将为优化沙地治理模式、建立跨区域生态补偿制度提供科学依据，对保障北方生态安全具有重要实践价值。

1 材料与方法

1.1　研究区概况

科尔沁沙地位于内蒙古高原向东北平原的过渡地带（图1）。行政区域主要包括内蒙古自治区通辽市、赤峰市和兴安盟科尔沁右翼中旗，此外，还涉及辽宁省西部彰武县和吉林省通榆县部分，本文的研究范围位于内蒙古自治区的15个旗（县、区、市）^［33］，总面积125.51万km²。年降水量为300~500 mm，集中于7—9月，年蒸发量1 800~2 000 mm。地势自西向东逐渐下降。冬、春季受蒙古高压影响，大风天气多，地表干燥，常出现强沙尘暴^［34］。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 研究区位置示意图

注：基于自然资源部标准地图服务网站标准地图（审图号为GS（2016）0650号）制作，底图边界无修改

Fig.1 Location of the study area

1.2　数据来源及预处理

本文使用的数据包括气象数据、土壤数据、雪深数据、植被数据等，具体见表1。气象数据中温度、降水、风速，均来源于国家气象信息中心（https：//data.cma.cn），利用ANUSPLIN软件对研究区内及邻近的45个气象站数据进行插值。地表潜在蒸散发数据来自国家青藏高原科学数据中心发布的1901—2022年高精度气象指标栅格数据，分辨率为1 km。数字高程模型（DEM）数据来自地理空间数据云平台（https：//www.gscloud.cn/），分辨率为30 m。土地利用（LUCC）数据来自中国科学院资源与环境数据中心（https：//www.resdc.cn），分辨率为30 m。土壤数据，包括土壤有机碳、沙粒、粉砂、黏土以及碳酸钙含量数据，来自世界土壤数据库（http：//www.fao.org），分辨率为1 km。雪深数据来自国家青藏高原科学数据中心（https：//data.tpdc.ac.cn）的中国雪深长时间序列数据集（1979—2023年），空间分辨率为0.25°。植被覆盖度（FVC）数据来源于MODIS-NDVI（MOD13Q1）数据集，该数据集提供了16天合成的无云产品，并通过Google Earth Engine（GEE）进行计算。在模拟轨迹时，采用ARL的再分析资料（https：//ready.arl.noaa.gov/），该数据包括U风速（m·s^-1）、V风速（m·s^-1）、相对湿度（%），不同高度的位势高度（gpm）。所有栅格数据集均被重采样至1 km，重投影为WGS-1984阿尔伯斯等面积投影。

表1 数据来源

Table 1 Data source

数据分类	数据名称	时、空分辨率	数据来源
RWEQ模型驱动	气象数据	逐日	国家气象信息中心
	土壤数据	1 km	世界土壤数据库
	DEM	30 m	美国国家航空航天局
	地表潜在蒸散发	逐日/1 km	国家青藏高原科学数据中心
	雪深数据	逐日/0.25°	国家青藏高原科学数据中心
HYSPLIT模型驱动	再分析资料	6 h/2.5°	美国国家海洋和大气管理局（https://ready.arl.noaa.gov/）
XGBoost-SHAP 模型驱动	气温	逐日	国家气象信息中心（https://data.cma.cn）
	降水
	大风天数
	植被覆盖度	16天/500 m	GEE

新窗口打开| 下载CSV

1.3 研究方法

1.3.1 防风固沙服务功能评估

本文通过修正风蚀方程（RWEQ）模型评估科尔沁沙地2000—2020年防风固沙量^［35-37］。

G = S L_{r} - S L

（1）

式中：G表示防风固沙量（kg·m^-2）；SL_r和SL分别表示潜在风蚀量和实际风蚀量（kg·m^-2）；SL_r为无植被覆盖情况下的最大可能风蚀量；SL则是在当前植被覆盖条件下的实际风蚀量。

S L_{r} = \frac{2 z}{S_{r}^{2}} Q_{r m a x} \cdot e^{- (z / S_{r})^{2}}

（2）

S L = \frac{2 z}{S^{2}} Q_{m a x} \cdot e^{- {(z / S)}^{2}}

（3）

式中：Q_rmax和Q_max分别表示潜在条件下风力下最大输沙量和当前风力条件下最大输沙量（kg·m^-1）；S_r和S分别表示潜在关键地块长度和实际临界地块长度（m）；z是计算得到的下风向距离（m），通常取50 m。

Q_{r m a x} = 109.8 \times (W F \times E F \times S C F \times K')

（4）

S_{r} = 150.71 \times {(W F \times E F \times S C F \times K')}^{- 0.3711}

（5）

Q_{m a x} = 109.8 \times (W F \times E F \times S C F \times K' \times C)

（6）

S = 150.71 \times {(W F \times E F \times S C F \times K' \times C)}^{- 0.3711}

（7）

式中：WF是气象因子（kg·m^-1）；EF是土壤可蚀性因子；SCF是土壤结皮因子；K'是土壤粗糙度因子；C表示植被因子。

W F = w f \times (ρ / g) \times S W \times S D

（8）

w f = u_{2} \times (u_{2} - u_{1}) \times N_{d}

（9）

式中：wf为风力因子（m·s^-³）；g为重力加速度（m·s^-2）；ρ为空气密度（kg·m^-³）；SW为土壤湿度因子；SD为雪盖因子，由雪深数据计算得到；u₁表示沙粒开始被风吹起的临界风速，取值为5 m·s^-1；u₂为气象站测得的月平均风速（m·s^-1）；N_d为每月风速超过5 m·s^-1的天数，该数值通过气象数据的统计分析得出。

E F = \frac{29.09 + 0.31 s a n d + 0.17 s i l t + 0.33 (s a n d / c l a y) - 2.59 O M - 0.95 C a C O_{3}}{100}

（10）

S C F = 1 / [1 + 0.0066 \times {(c l a y)}^{2} + 0.021 \times {(O M)}^{2}]

（11）

式中：EF是土壤可蚀性因子；sand表示土壤中砂粒含量（%）；silt表示土壤中粉粒含量（%）；clay表示土壤中黏粒含量（%）；OM表示土壤中的有机质含量（%）；由于碳酸钙（CaCO₃）数据缺失，因此未将其纳入计算。

C = e^{- 0.0483 \times F V C}

（12）

F V C = (N D V I - N D V I_{s o i l}) / (N D V I_{v e g} - N D V I_{s o i l})

（13）

式中：C是植被因子；FVC为植被覆盖度；NDVI_soil为完全是裸土或无植被覆盖区域的NDVI值；NDVI_veg为完全被植被所覆盖像元的NDVI值。

K' = c o s α

（14）

式中：K'是土壤粗糙度因子；α表示坡度，通过ArcGIS中的数字高程模型（DEM）计算得出。

为避免气候因素对防风固沙服务的影响，本文使用防风固沙量（G）与潜在土壤风蚀量（SLr)的比值作为防风固沙保有率^[15]。

1.3.2　气团轨迹模拟

本文利用HYSPLIT模型的前向轨迹模块模拟科尔沁沙地沙尘的传输。HYSPLIT模型是大气传输和扩散模型，已被广泛地应用于多种污染物在各地区传输和扩散的研究中^［38-39］，Xu等^［28］将其用于模拟防风固沙服务研究。HYSPLIT模型计算气团轨迹使用拉格朗日法：根据空气中随风飘动粒子的运动参数随时间变化规律推算粒子在一定时间内的移动轨迹。移动轨迹由初始位置 $P (t)$ 和第1个推测位置P′（t+Δt）三维速度矢量的平均值计算^［40］。第1个推测位置的计算方法为：

P' (t + Δ t) = P (t) + V (p, t) + Δ t

（15）

最终位置的计算公式如下：

P (t + Δ t) = P (t) + 0.5 \times [V (p, t) + V (p', t + Δ t)] + Δ t

（16）

式中： $P (t)$ 表示初始位置； $V (p, t)$ 表示初始位置处的速度向量；Δt表示积分时间步长； $V (p', t + Δ t)$ 表示第一次预测位置处的速度向量。

采用MeteoInfo软件中的TrajStat模块，依据HYSPLIT模拟轨迹之间总空间变异度显著增加这一准则，对其进行聚类分析^［41］，识别主要的沙尘输送路径。

1.3.3　影响因素分析

XGBoost（eXtreme Gradient Boosting）是基于梯度提升框架的集成学习算法^［42］，通过迭代构建决策树组合优化预测性能，在回归、分类及排序任务中具有广泛应用。在传统梯度提升算法基础上改进，包括加权处理数据分布、正则化防止过拟合、支持计算加速训练，并能自动处理缺失值。凭借强大的学习能力和优良的泛化性能，XGBoost在大规模高维数据集处理方面具有显著优势^［43］。本研究采用XGBoost模型评估影响因素对风蚀的作用，并在模型训练过程中进行了超参数调优。

为解析XGBoost模型的黑箱特性，本研究采用基于合作博弈论Shapley值的SHAP（SHapley Additive exPlanations）解释框架^［44］。该方法通过计算特征边际贡献的加权平均值。

Φ_{i} = \sum_{S \subseteq F \ \{j\}} \frac{|S|! (|F| - |S| - 1)!}{|F|!} [f (S ⋃ \{j\}) - f (S)]

（17）

式中：Φ_i 表示SHAP值，反映特征在模型预测中的重要性，值越大，说明该特征对模型预测结果的贡献越大，反之则越小；F为特征全集；S为特征子集； $| S |$ 是集合S中的特征数； $f (S)$ 是特征集S对模型预测输出的贡献； $f (S \cup {j}$ ）是包含特征j的特征集 $S \cup {j}$ 对模型预测输出的贡献。此外，SHAP细化到每个样本的特征贡献，直观展示各特征在风蚀变化中的作用。

2 结果分析

2.1　防风固沙服务时空变化特征

科尔沁沙地2000—2020年防风固沙服务呈显著的时间异质性。从单位面积服务能力来看，年固沙量为5 232.45~8 194.69 t·km^-2·a^-1，年平均值6 528.9 t·km^-2·a^-1（图2A）。从变化趋势来看，单位面积固沙量呈现阶段性演变特征，2000—2010年单位面积固沙量呈现递减趋势，平均变化率为146.51 t·km^-2·a^-1，而2010年后趋势逆转，以65.97 t·km^-2·a^-1的速率回升，反映近10年生态治理成效逐步显现（图2A）。2000—2020年防风固沙保有率的年际变化呈上升趋势，以每年0.17%的速率增加（图2B），表明2000—2020年研究区防风固沙服务能力增强。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 2000—2020年单位面积固沙量和防风固沙保有率

Fig.2 Unit area and preservation rate of WEPS from 2000 to 2020

从空间分布来看（图3），科尔沁沙地单位面积防风固沙量呈中南高北部低的特征。高值区（>8 000 t·km^-2）主要分布于研究区中部和南部，集中在科尔沁左翼后旗、科尔沁区、库伦旗北部、奈曼旗北部、敖汉旗北部等区域。低值区主要分布于研究区北部，分布在突泉县、科尔沁右翼中旗、霍林郭勒市、扎鲁特旗等区域。2000—2020年防风固沙总量呈现先降后升的V形波动曲线，研究区防风固沙总量在10.6亿~16.6亿t，均值为13.1亿t，最大值和最小值分别出现在2000年和2012年。具体而言，2000年固沙总量达到研究期峰值16.6亿t，此后10年间急剧下降至2010年的12.4亿t（降幅达25.3%），反映出此阶段人类活动干扰与气候变化叠加导致的生态退化风险，2010年后生态治理成效凸显，固沙总量于2015年降低至11.8亿t，2020年小幅增加至13.7亿t。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 2000—2020年防风固沙服务空间分布

注：基于自然资源部标准地图服务网站标准地图（审图号为GS（2016）0650号）制作，底图边界无修改

Fig.3 The spatial distribution of WEPS from 2000 to 2020

2000—2020年科尔沁沙地微度和轻度单位面积固沙量分别为52.61、1 699.70 t·km^-2·a^-1，占比面积合计最大，达到53.87%。中度和强烈单位面积固沙量分别为3 300.62、6 455.74 t·km^-2·a^-1，占比面积合计16.16%。极强烈和剧烈单位面积固沙量分别为11 653.15、18 989.65 t·km^-2·a^-1，占比面积合计29.96%。2000—2020年微度等级面积占比稳定在8%~9%，轻度等级占比从2000年27%跃升至2020年42%，说明生态修复工程显著提升了初级固沙效能（表2）。与之形成对比，中度等级占比由31%降至17%，揭示中度侵蚀区正向低等级转化的积极态势；强烈、极强烈和剧烈等级占比在降低，显示高强度侵蚀得到有效遏制。

表2 侵蚀强度分级及占比 (%)

Table 2 Classification and proportion of WEPS

等级	侵蚀强度 /（t·km^-²）	年份
等级	侵蚀强度 /（t·km^-²）	2000	2005	2010	2015	2020
微度	<500	8	8	8	9	8
轻度	1 500~2 500	27	54	51	52	42
中度	2 500~5 000	31	7	10	8	17
强烈	5 000~8 000	4	2	2	3	3
极强烈	8 000~15 000	6	13	9	7	5
剧烈	>15 000	25	17	21	21	24

新窗口打开| 下载CSV

从时间演变规律分析，2000—2020年研究区平均防风固沙保有率呈现波动上升特征，由2000年的81.25%提升至2020年的86.19%，增幅达4.94%。具体表现为三阶段演化：2000—2005年快速提升期（年均增长率1.10%），2005—2010年短暂回调期（保有率下降2.56%），2010年后进入稳态增长期（年均增速0.15%）。这种波动轨迹与生态工程的阶段性推进及降水年际变率密切相关，特别是2010年后的平缓增长，反映生态修复进入提质增效阶段。空间分异格局揭示出与固沙量分布相反的特征（图4），防风固沙保有率较高的区域主要位于北部地区，中部沙地的防风固沙保有率较低。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 2000—2020年防风固沙保有率空间分布

注：基于自然资源部标准地图服务网站标准地图（审图号为GS（2016）0650）制作，底图边界无修改

Fig.4 Preservation rate of WEPS from 2000 to 2020

2.2　防风固沙服务驱动因素

2.2.1　基于 XGBoost 回归模型分析

风蚀受多种因素的影响，包括风速、地表覆盖度和粗糙度、气候条件、植被条件、土地利用/覆盖变化等^［45］，这些因素共同影响风蚀程度。本文选取风速（WIN）、气温（TEM）、降水（PRE）、植被覆盖度（FVC）4个自然因子，探讨其与实际风蚀的关系。多重共线性检验结果发现WIN、TEM、PRE、FVC的膨胀方差因子（VIF）均小于10，表明各因子之间不存在共线性问题。

XGBoost回归模型经过超参数调优，得到最优参数组合，这些模型表现出较高的解释能力。其中2000、2010、2020年的R²分别为0.8911、0.8853、0.8423，均方根值MSE分别为4.3555、1.9361、3.9810，均方根误差RMSE分别是2.0870、1.3914、1.9952，高R²值和较低的误差指标表明模型能够有效捕捉数据的模式和关系，所选的因子能够解释实际风蚀的空间变异。

2.2.2　基于SHAP模型的解释性分析

基于XGBoost机器学习算法，使用SHAP算法影响因子进行重要性评价，较高的SHAP值表示该因素对风蚀度的影响较大。由图5可知，2000年主要的驱动因素是FVC和TEM，2010年主要的驱动因素是FVC和WIN，2020年主要的驱动因素是FVC和PRE。总之，FVC是抑制风蚀最重要的因素，因此在所有年份都有最高的SHAP值。其他3个因素是风蚀发生的环境条件，对应的SHAP值在不同年份波动，反映当年特定的气候条件放大了该因子对风蚀过程的影响强度。各影响因素对风蚀的作用方向在不同年份大致相同，FVC在各年份均对实际风蚀的影响最大，是核心驱动因子，其SHAP值为5~20，跨度较大。第二是WIN，其SHAP值为-5~5，对风蚀主要起正向影响。其次是TEM，对风蚀具有正向影响，其SHAP值为-5~5。最后是PRE，其SHAP值为-5~5，对风蚀具有负向影响。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 SHAP重要性蜂群图和条形图

注：WIN：风速；TEM：气温；PRE：降水量；FVC：植被覆盖度

Fig.5 SHAP importance swarm plot and bar chart in 2000， 2010 and 2020

为了更好地明确影响因子对风蚀的影响，本研究分析了影响因子与风蚀之间的关系（图6）。具体来说，FVC值与风蚀呈现指数型关系，当FVC值增至0.4左右时，SHAP急剧下降，表明植被覆盖度增加，可使风蚀风险降低，这印证了植物根系固结土体、冠层削减近地表风速的双重防护机理，随着FVC的增加，SHAP值逐渐趋于平稳，意味着较高的植被覆盖率对减少风蚀有积极作用。气象因子与风蚀的关系更为复杂，通过非线性作用影响风蚀过程。PRE的SHAP值显示降水量对风蚀的影响不是简单的线性关系，当PRE较小时，风蚀的SHAP值较大，随着PRE的增加，风蚀趋于稳定，反映充沛降水通过促进植被生长间接抑制风蚀。TEM的SHAP值的变化显示，当TEM较小时，SHAP值变化较为稳定，但是当TEM升高时，SHAP值急剧上升，波动较大。WIN和风蚀的关系无显著线性关系，表明风速对风蚀的影响存在非线性或阈值效应。风速较小时，风蚀的SHAP值较稳定，但是当风蚀增加，SHAP值急剧增大。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 影响因素与风蚀的SHAP散点图

Fig.6 SHAP dependence plots of influencing factors and wind erosion

2.3　沙尘传输路径模拟

本研究以科尔沁沙地内的站点为HYSPLIT模型模拟的起始点，利用ARL气象再分析资料，通过设定风速阈值（≥5 m·s^-1）筛选符合起沙条件的天气过程，模拟气团在起沙事件发生后120 h（5 d）的前向扩散路径。2000—2020年共有3 559条风沙轨迹路径（图7），年际分布特征显著，轨迹数量最大值出现在2001年（282条，占比7.92%），最小值出现在2014年（77条，占比2.16%）。2000—2020年线性回归分析表明，轨迹数量年际变化呈现显著负相关趋势（图7），风沙轨迹数量呈现下降的趋势，表明防风固沙能力的增强。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 2000—2020年风沙轨迹数量图

Fig.7 Number of wind sand trajectories from 2000 to 2020

本研究通过MeteoInfo软件的TrajStat插件对2000—2020年风沙前向轨迹进行聚类分析，揭示沙尘输运路径的时空演变规律（图8A~F）。2000—2020年研究区存在4条主导路径（图8A），31.67%的轨迹向东输送至太平洋，29.92%向东北方向影响吉林和黑龙江，20.46%沿朝鲜半岛扩散至日本等境外区域，17.95%经黄海、渤海抵达上海、浙江、福建等东南沿海。图8B~F分别表示2000、2005、2010、2015年和2020年的轨迹聚类。2000年轨迹以东南方向为主（90.94%），其中，52.83%的轨迹沿着辽宁到达朝鲜半岛，26.79%的轨迹继续延伸到达日本，11.32%的轨迹吹向中国东南沿海地区，而东北路径仅占9.06%，吹向吉林和黑龙江方向。2005年34.52%的轨迹影响黄海、东海区域，30.46%延伸至朝鲜半岛及日本，24.87%的轨迹吹向了东北部的吉林和黑龙江，而其余10.15%的轨迹吹向了西北部的内蒙古。2010年73.19%的轨迹向东扩散，26.81%扩散至黄海、东海；2015年，87.64%的轨迹沿辽宁—朝鲜半岛—日本轴线输送，仅12.36%进入东海。2020年路径复杂性增强，48.83%的轨迹覆盖东北亚跨境区域，35.21%影响内蒙古-吉林，15.96%沿中国海岸线延伸至南海。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 2000—2020年沙尘轨迹聚类

注：基于自然资源部标准地图服务站标准地图（审图号为（2016）1600号）制作，底图边界无修改

Fig.8 Trajectory clustering from 2000 to 2020

3 讨论

3.1　防风固沙量时空变化

本文的研究结果表明科尔沁沙地单位面积年均防风固沙量为6 528.9 t·km^-2·a^-1，年均总量为13.1亿t，研究结果高于邵国媚等^［46］、朱梦媛等^［47］的研究结果，前者的年均总量为3.71亿t，后者的单位面积年均防风固沙量为2 160 t·km^-2·a^-1，总量为1.40亿t，主要是由于研究区范围的不一致。邵国媚等^［46］以科尔沁生态功能区为研究区域，主要包括内蒙古自治区的11个旗县，朱梦媛等^［47］以科尔沁沙地南缘为研究区，主要是辽西北的狭长地带，而本文的研究区包括的是内蒙古自治区的15个旗（县、区、市），使得下垫面状况、土壤可蚀性、风力因子等发生深刻变化^［48］，导致防风固沙量不同。

本研究揭示，研究区防风固沙能力呈现先降后升的动态演变特征。这一趋势与黄麟等^［49］对北方风蚀区退耕还林工程、迟文峰等^［50］对内蒙古高原的评估结果一致，这是气候暖湿化和生态工程治理双重驱动的结果。为改善区域生态状况，各区域实施“三北”防护林体系建设、退耕还林还草等重大生态修复项目，工程实施初期，受干旱气候和水资源短缺制约^［51］，新造林地普遍存在幼苗成活率偏低、植被覆盖度不足等问题，导致地表植被防护效能阶段性弱化，防风固沙功能出现短期回落。但随着京津风沙源治理工程的持续推进，区域植被盖度显著提升，生物多样性逐步恢复，叠加土地利用结构的优化调整，有效抑制了土壤风蚀过程^［45］。

3.2　防风固沙的驱动因素

本文基于XGBoost-SHAP分析了风蚀的驱动因素，结果表明植被覆盖度是21世纪以来科尔沁沙地风蚀演变的主导调控因子。这一结论与研究区实施的系列生态恢复工程密切相关，自2000年“三北”防护林等生态恢复工程实施以来，加之对于天然林保护、生态公益林保护和退耕还林还草等政策协同效应^［52］，区域植被覆盖度大幅度提升，数据显示20年来植被明显恢复地区面积超过64.91%，主要分布在科尔沁沙地北部、中部以及东南缘^［53］。SHAP值的蜂群图显示植被覆盖度和风蚀呈正相关关系，主要是植被通过直接覆盖地表、吸收气流动量、拦截阻滞跃移颗粒、吸附细颗粒，减弱土壤风蚀^［54］，植被覆盖度越大，对土壤保护作用越大^［55］。在气候要素交互作用方面，研究揭示了非线性的响应特征。风速与风蚀强度呈正相关关联，但二者未呈现简单线性关系，可能依赖于植被覆盖度等其他变量，植被覆盖度发挥着关键缓冲调节作用，当植被覆盖度达阈值时，其防风固沙效应可部分抵消风速增强带来的侵蚀。温度对风蚀有着正向驱动作用，这与Ma等^［56］和Gao等^［57］的研究结果一致，全球范围内以及中国北方，气温上升趋势导致风蚀增加，主要是气温升高加速土壤水分蒸散，降低黏粒胶结作用，增强大气湍流，提升风动力条件。多因素耦合作用凸显了气候变暖对风蚀过程的复杂影响机制。降水要素呈现双刃剑效应，主要是正向驱动，研究区大部分地区属于干旱半干旱地区，水是制约植被长势的主要因子^［48］。降水量增加，有利于植被恢复和土壤湿度增大，抑制风蚀的发展^［58］，但本文的研究表明，降水和风蚀之间并不是线性关系，部分区域也呈负相关关系，可能是由于高降水地表径流冲刷破坏地表结构，反而加剧风蚀，这种非线性响应印证了干旱区生态-水文过程的脆弱复杂特性。本研究结果对生态工程效益评估具有重要启示，在生态工程实施20年后，植被恢复已显著改变区域风蚀动力系统，但气候变化引发的极端天气频率增加可能削弱生态防护效能。未来需建立动态适应性管理机制，在植被重建中优先选择抗旱耐风蚀物种，同时构建水土保持工程与植被措施协同的复合防护体系，以应对气候变化带来的新型风蚀风险。

3.3　沙尘的传输路径

沙尘以常呈湍流状态的风为载体，其运动方向在各时刻均有很大不同，因此，沙尘传输属于全向型。在各方向气流环境差异背景下，大气环境存在一定的不稳定性，但尚未有更好的解决方案，本文将沙尘传输进行简化，假设风速大于临界起沙风速时，沙源地的沙尘即随气流运动。科尔沁沙地位于东亚季风活跃区，距离蒙古高压常年平均中心位置不足1 500 km^［59］，沙尘传输以西北-东南向为主，加之流域的盛行风向为西北风，东南方向的传输路径比西北方更多。研究结论符合这一客观事实，研究时段东南方向的沙尘传输频率较西北方更为频繁。

3.4　不确定性

本文的不确定性主要体现在这3个方面：①基础数据层面，风速数据采用的方法是结合气象站点风速数据与插值软件ANUSPLIN得到风速的栅格数据，因此，在风速插值的结果中存在不确定性。②动力机制层面，将沙尘的流动路径简化，假定只要在风速大于等于临界起沙风速的条件下，有充足的沙源可供流动，没有考虑大气层结的不稳定性。③驱动因素方面，鉴于数据的可获得性，当前模型主要考虑了风速、温度、植被覆盖度和降水4个核心自然参数，但对放牧强度、牲畜数量、土地利用动态、农业生产方式调整、生态保护政策变化等人类活动因子缺乏定量表征。未来可通过构建高密度气象观测网络、发展多源遥感数据同化技术、引入大气边界层数值模式耦合等方法，系统提升模型对复杂地表风蚀过程的多尺度模拟能力，为干旱区风蚀防治提供更精准的决策支持。

4 结论

本文利用RWEQ模型计算了2000—2020年科尔沁沙地的防风固沙服务，并分析了其时空格局变化特征，利用XGBoost-SHAP算法分析了风蚀的驱动影响因素，最后利用HYSPLIT模型解释了风沙的流动轨迹。结论如下：防风固沙服务呈现出显著的时空异质性特征，年固沙量为5 232.45~8 194.69 t·km^-2，年平均值6 528.9 t·km^-2，单位面积固沙量呈先降后升趋势，空间格局表现为中南高-北部低的分异特征。

使用SHAP算法对影响因子进行重要性评价，植被覆盖度是实际风蚀的核心驱动因子，植被覆盖度值与风蚀呈现指数型关系，随着植被覆盖度的增加，SHAP值逐渐降低。其次的影响因素是风速，主要是正向影响，气象因子主要是通过非线性交互作用对风蚀产生影响。

HYSPLIT模拟揭示2000—2020年共产生3 559条风沙传输路径，轨迹数量整体下降，表明固沙能力提升，聚类分析识别出4条主导路径：31.67%向东入太平洋，29.92%向东北影响吉林和黑龙江，20.46%跨境扩散至朝鲜、日本等区域，17.95%南下抵达东南沿海。本文的研究结果可为优化沙地治理策略、协调区域生态补偿机制提供科学依据，同时为干旱半干旱区生态系统服务流动研究提供方法参考。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Ravi

， D.Odorico

， Over

T M

，et al.

On the effect of air humidity on soil susceptibility to wind erosion：the case of air-dry soils

［J］.Geophysical Research Letters，2004，31（9）：L09501.