“三北”防护林工程（简称“三北”工程）是中国为防治西北、华北和东北地区风沙危害和水土流失等灾害而启动的第一项国家重点生态工程，也是全球规模最大的生态修复工程^［1］。目前正处于六期工程实施阶段。“三北”工程六期工程区分为核心攻坚区、协同推进区和巩固拓展区^［2］。其中，核心攻坚区是“三北”工程区内风沙灾害最严重的区域，也是“三北”工程六期核心治理区和三大标志性战役所在区域（黄河“几字弯”攻坚战、科尔沁和浑善达克沙地歼灭战、河西走廊-塔克拉玛干沙漠边缘阻击战，分别简称攻坚战工程区、歼灭战工程区和阻击战工程区）。区域任务量约占“三北”工程六期总任务量的2/3^［2］。认清区域土壤风蚀状况是确定风蚀防治重点、科学安排风蚀防治工程的前提。

风蚀模型是评估区域尺度土壤风蚀状况的主要手段。常用的风蚀模型主要有土壤风蚀方程（WEQ）、修正土壤风蚀方程（RWEQ）、风蚀预报系统（WEPS）和第一次全国水利普查风蚀模型。WEQ是美国学者为解决美国农田风蚀问题，构建的一个涵盖土块、土垄粗糙度、田块尺度、气候和植被5大因子的经验模型^［3］。为解决WEQ理论基础不足和实际应用效果差的问题，美国农业部同时启动了RWEQ和WEPS模型研究。RWEQ模型继承了WEQ模型结构，仍是针对特定区域的经验性模型^［4］。WEPS将风蚀影响因子以子模型形式表达，实现了高度程序化运行，但模型中的参数和常数针对美国制定，目前仍处于试验和完善阶段^［5］。第一次全国水利普查风蚀模型是针对中国北方风蚀地表类型建立的经验模型，模型结构简单、容易操作，可直接应用于中国北方地区的风蚀模拟，得到行业部门的认可和采纳，也成为水利部相关标准规范中使用的风蚀模型^［6-7］。后期部分学者根据该模型在实际应用中存在的问题进行了修正，提升了模型精度^［6，8］。

目前，已有不少学者利用风蚀模型在“三北”工程核心攻坚区内开展土壤风蚀时空演化及其驱动机制的相关研究。Du等^［9］利用RWEQ模型和控制变量法评估了生态工程对2001—2019年三北工程区土壤风蚀变化的贡献。邢瑜等^［10］利用WEPS模型，分析了河西走廊-塔克拉玛干沙漠边缘阻击战核心区2000—2023年土壤风蚀时空变化及其影响因素。同时，部分学者利用RWEQ模型或改进后第一次全国水利普查风蚀模型分别评估2022年以前核心攻坚区内部分地区（黄河流域^［11］、毛乌素沙地^［12］、塔河流域^［13］和京津风沙源^［6］等）土壤风蚀的时空演化。然而，目前尚未见到对整个“三北”工程核心攻坚区最新土壤风蚀状况的报道。

本研究选取修订的第一次全国水利普查风蚀模型、非参数趋势分析、极端梯度提升（XGBoost）模型和控制变量法，估算2000—2024年核心攻坚区逐年的土壤风蚀模数，分析区域土壤风蚀模数的空间分布格局及其时间演变规律，并进一步探讨了土壤风蚀时空分布格局的驱动因素。

核心攻坚区（34.78°—48.86°N，73.51°—126.18°E），东西长近4 300 km，南北宽近1 750 km，横跨中国北方13省（自治区、直辖市），总面积达302.07万km²（图1）。区域自然本底异质性强。地貌上，跨越了三级阶梯，涵盖平原、高原、山地、盆地等多种地貌类型。气候上，跨越了季风性气候和大陆性气候，2000—2024年区域年降水量14~1 200 mm，大部分区域年蒸发量超过1 500 mm^［14-16］。受降水和蒸发影响，核心攻坚区自东向西干旱程度不断加重，依次为半湿润、半干旱、干旱和极端干旱区。2000—2024年区域年平均风速为2.07~4.42 m·s^-1，中国北方多个大风区分布于此。攻坚战工程区和歼灭战工程区大风多发于冬春季节，阻击战工程区大风集中在春夏季节。受区域自然本底影响，其间分布着中国的四大沙地和七大沙漠，是中国土地沙漠化程度最高、问题最突出的区域。受气候等自然环境影响，核心攻坚区植被类型多样，主要包括荒漠、草原、栽培植被和草甸等。同时，核心攻坚区土壤类型也比较丰富，主要的土壤类型包括风沙土、栗钙土、棕漠土、灰棕漠土、棕钙土、黄绵土和寒钙土等。受自然条件影响，核心攻坚区生态环境脆弱、生态本底较差，是中国土壤风蚀问题最严重的区域。

核心攻坚区东西跨度大、自然环境复杂，人类活动强度差距显著。以2024年为例，区域人口密度从阿拉善盟的不足1人·km^-2到渭南市的近351人·km^-2。全区土地利用类型以未利用土地、草地和耕地为主。2020年，未利用土地、草地、耕地面积分别为1 449 267、995 608、345 738 km²，占区域总面积的47.99%、32.97%、11.45%，畜牧业和农业是区域主要的人类活动。近年来区域人口不断增多、社会经济得到快速发展。同时，区域内土地利用也发生了较大变化，2000—2020年区内耕地面积增加了9 397 km²、建设用地面积增大了14 356 km²。人类活动强度的增大加重了原本脆弱生态系统的压力，使之成为中国北方土壤风蚀和土地沙漠化的活跃地区。

截至2020年“三北”防护林工程已历时43年，完成了五期工程，目前正处于六期工程实施阶段。“三北”防护林工程的实施有效遏制了工程区沙化土地发展和风沙危害，基本控制住了水土流失，工程区生态环境得到了较大改善。

本研究所用数据包括遥感、气象、土地利用、土壤等。其中，遥感数据包括归一化植被指数（NDVI）和表土湿度数据。NDVI产品数据为覆盖全区2000—2024年MODIS 1 km/16 d NDVI产品数据，由MODIS网站（https：//ladsweb.modaps.eosdis.nasa.gov/search/）下载获得；表土湿度数据为2000—2024年美国航空航天局戈达德空间飞行中心与美国海洋和大气局国家环境预报中心联合发布的表土湿度产品数据，由EARTHDATA网站（https：//search.earthdata.nasa.gov/search）下载获得。气象数据包括研究区及周边地区气象台站2000—2024年的逐月降水、逐月气温和逐时风速数据，由国家气象数据中心（http：//data.cma.cn/）获得，气象数据获取后使用ANUSPLIN插值得到覆盖全区的栅格数据。空间分辨率1 km的土地利用数据（2000、2005、2010、2015、2020年）由资源环境科学与数据中心（https：//www.resdc.cn/）下载获得。土壤数据来自世界土壤数据库（HWSD）。所有数据均统一处理为空间分辨率1 km、Albers Equal-Area Conic/WGS84投影坐标系统的栅格数据。

本文采用修订的第一次全国水利普查风蚀模型^［8］。该模型是在水利部第一次水利普查风蚀模型基础上，增加了戈壁风蚀子模型。

式中：Q₁、Q₂、Q₃和Q₄分别为耕地、林草地、戈壁和沙地（漠）的土壤风蚀模数（t·km^-2·a^-1）；z₀为地表空气动力学粗糙度（cm），对于中国北方传统耕作方式，翻耕后裸露耕地表面的空气动力学粗糙度约为0.55 cm；A为与下垫面有关的风速修订系数；U_j 为气象站整点风速统计中高于临界侵蚀风速的第j级风速（m·s^-1）；T_j 为风速为U_i 的累积时间（min）；VC是春季植被覆盖度（%），利用像元二分法和春季NDVI数据计算获得；W为表土湿度因子，由表土湿度产品数据换算获得；GC为砾石覆盖度；a₁、b₁、c₁、a₂、b₂、c₂、a₃、b₃和c₃为常数项，分别取值-9.208、0.018、1.955、2.4869、-0.0014、-54.9472、6.1689、-0.0743和-27.9613。土壤风蚀模数核算后，依据《土壤侵蚀分类分级标准》（SL 190-2007）^［7］，确定土壤风蚀强度等级（表1）。

本研究采用Theil-Sen Median趋势分析和Mann-Kendall显著性检验对核心攻坚区2000—2024年土壤风蚀时间演化规律进行分析。

式中：β为变量x的线性变化趋势，β>0和β<0分别表示变量x随时间呈增加和减小的趋势；x_f 和x_e 分别表示第f和第e时间点变量x的像元值；S为所有正差异与负差异数量的统计值；n为时间序列的总长度（取值25）；sgn为符号函数；Z为标准化统计量，在0.05的显著性水平下，当|Z|≤1.96表明土壤风蚀变化趋势不显著，当|Z|>1.96表明土壤风蚀变化趋势显著；Var（S）为统计值S的方差（取值1833.3）。

极端梯度提升（eXtreme Gradient Boosting，XGBoost）模型是一种基于梯度提升决策树框架的高效机器学习算法^［15］。该模型通过引入损失函数的二阶泰勒，迭代地训练决策树模型将精度较低的多个模型组合成一个精度较高的模型，增加了模型对最优解的识别能力。

式中：Obj为目标函数；L为损失函数；Ω为正则化项。本文选用XGBoost模型分析土壤风蚀模数与自然和人为因素之间的关系。使用2 457 920个样本进行模型构建（其中70%用于训练，剩余30%用于验证）。随着训练次数增加，训练数据集与验证数据集的均方根误差均逐渐减小，并在1 653次后趋于稳定终止训练，计算获得最终结果。基于验证集的测试显示模型拟合良好（R²=0.82）。

SHAP（Shapley Additive Explanations）是一种基于博弈论、解释机器学习模型输出结果、因变量与各因素间复杂关系的方法。本文引入SHAP值来解释XGBoost模型输出结果中隐含的非线性关系，用于量化每个因素对模型预测结果的贡献大小和方向。

式中：y_base为整个模型的平均值；x_ij 为第i个样本的第j个特征，i∈（1，2，3，…，n），j∈（WIND，PRE，SAND，TEM，VC，PD，NLI）；f（x_ik ）为x_ij 的SHAP值。其中，f（x_ij ）为第i个样本中第j个特征对最终预测值yi^的贡献值。当f（x_ij ）>0时，表示该特征提升了预测值，具有正向影响；当f（x_ij ）<0时，表示该特征降低了预测值，具有负向影响。

控制变量法是目前定量评估各因素对土壤风蚀变化贡献率的常用方法^［15-16］。本文通过分别将风速、土壤湿度、植被覆盖度和土地利用固定为2000年的恒定值，其他因素变化条件下计算不同情形下的风蚀模数。不同情形下风蚀模数与实际风蚀模数之差即为不同因素对土壤风蚀变化的贡献率。

式中：CR_i 为因素i的贡献率；Q_Ci 和Q_Pi 分别为因素i变化情况下的土壤风蚀量和土壤风蚀变化量；Q_A 为实际风蚀量。

2000—2024年核心攻坚区年均土壤风蚀量1.19×10⁸ t，风蚀模数39.47 t·hm^-2·a^-1（图2）。核心攻坚区以微度和轻度风蚀为主，面积2 398 346 km²，占区域总面积的79.40%，广泛分布区域内。强烈、极强烈、剧烈风蚀区面积分别为83 277、67 668、195 407 km²，占区域面积的2.76%、2.24%、6.47%，集中分布在塔克拉玛干沙漠、库姆塔格沙漠、柴达木沙漠、巴丹吉林沙漠、腾格里沙漠、乌兰布和沙漠、浑善达克沙地西部等区域。中度风蚀区面积276 049 km²，占区域总面积的9.14%，主要分布在强烈等级以上风蚀区周边环境相对较好的区域。不同分区土壤风蚀强度不同。阻击战工程区风蚀强度最大，年均土壤风蚀量9.79×10⁷ t，风蚀模数53.06 t·hm^-2·a^-1。其次为攻坚战工程区，风蚀模数19.30 t·hm^-2·a^-1。歼灭战工程区风蚀强度最小，年均土壤风蚀模数为16.87 t·hm^-2·a^-1。

本文选取了平均风速、平均气温、降水量、植被覆盖度、土壤沙粒含量5项自然因子和夜间灯光指数、人口密度2项人类活动因子，结合XGBoost模型计算各因子的SHAP值，分析自然和人类活动因子对核心攻坚区土壤风蚀空间分布的影响。总体来说，自然因子是控制土壤风蚀空间分布的主导因素，人类活动的影响十分有限（图3）。不同因子的贡献有所不同。其中，平均风速的SHAP值最高，且对土壤风蚀多呈正向贡献。风速较高时，其对土壤风蚀的促进作用更大，表明风速越高的区域土壤风蚀模数越大；风速较低时，风蚀模数也较小（SHAP值<0）。其次为沙粒含量与降水量，SHAP值分别为16.63和15.98。土壤沙粒含量越大，土壤可蚀性越强，土壤风蚀模数也随之增大。降水量对风蚀的影响呈现出明显的差异，降水量较少的区域，降水对土壤风蚀的影响更显著。植被覆盖度和平均气温的SHAP值较为接近，分别为7.60和7.52，二者对土壤风蚀的贡献方向相反。植被覆盖度高的区域，地表裸露程度低，风蚀模数小；植被覆盖度较低的区域，地表裸露程度大，风蚀严重。气温较高区域因蒸发强烈、植被稀疏、表土干燥，风蚀程度较为严重。人口密度和夜间灯光指数对土壤风蚀表现出较弱的负向贡献（SHAP值较低），这可能与人口和经济活动多集中在水热条件好、植被覆盖度大的区域有关。

与核心攻坚区相似，自然因子尤其是平均风速和沙粒含量也是控制3个分区土壤风蚀空间分布的主导因素，SHAP值最高；人类活动因子对3个分区土壤风蚀空间分布的影响有限，SHAP值较低（图3）。降水量、植被覆盖度和平均气温对各分区土壤风蚀空间分布的影响不同。其中，从阻击战工程区、攻坚战工程区到歼灭战工程区，水热条件和植被生长状况趋于好转，植被覆盖度对土壤风蚀空间分布的贡献也相应增大。

2000—2024年核心攻坚区土壤风蚀模数总体呈显著减小趋势（P<0.01），风蚀模数年减小速率为0.77 t·hm^-2·a^-2（图4）。不同阶段，土壤风蚀模数变化趋势有所不同。其中，2000—2017年土壤风蚀模数呈显著减小趋势，2017—2024年土壤风蚀模数呈显著增大趋势。2010年是区域土壤风蚀模数最大的年份，平均土壤风蚀模数为53.43 t·hm^-2·a^-1；2017年是风蚀模数最小的年份，平均土壤风蚀模数为21.98 t·hm^-2·a^-1。

不同分区土壤风蚀模数变化趋势不同。其中，阻击战工程区和攻坚战工程区土壤风蚀模数变化趋势与核心攻坚区基本一致，2017年为区域土壤风蚀模数最小年份和土壤风蚀模数变化趋势转折年。两分区土壤风蚀模数变化速率有所不同。阻击战工程区土壤风蚀模数减小趋势较大，为1.09 t·hm^-2·a^-2。攻坚战工程区风蚀模数年减小速率仅为0.34 t·hm^-2·a^-2。歼灭战工程区整体变化趋势虽与其他分区相似，但土壤风蚀模数变化趋势的转折年份不同。2000—2014年风蚀模数呈减小趋势，年减小速率为1.15 t·hm^-2·a^-2；2014—2024年风蚀模数呈增大趋势，年增大速率为0.62 t·hm^-2·a^-2。

空间上，土壤风蚀模数呈减小趋势的区域面积1 484 609 km²，占核心攻坚区总面积的49.15%，广泛分布在核心攻坚区内（图5）。其中，呈显著减小趋势的区域面积837 875 km²，集中分布在塔克拉玛干沙漠西部和东部、柴达木盆地、巴丹吉林沙漠、浑善达克沙地西部、科尔沁沙地西部和内蒙古自治区北部的部分地区。土壤风蚀模数呈增大趋势的区域面积712 138 km²，占核心攻坚区总面积的23.57%，主要分布在新疆维吾尔自治区中部、甘肃省中部、内蒙古自治区北部和南部的部分地区（图5）。其中，呈显著增大趋势的区域面积289 323 km²，集中分布在塔克拉玛干沙漠部分地区、河西走廊、毛乌素沙地、库布齐沙漠和锡林郭勒盟北部。

控制变量分析结果显示，风速是驱动核心攻坚区土壤风蚀变化的主导因子（图6），风速变化对土壤风蚀变化贡献率绝对值主要位于80%~100%。风速变化对土壤风蚀变化总体呈负贡献，全区贡献率平均值为-12.83%。风速减小有利于土壤风蚀减小的区域面积占全区面积的56.47%，风速增大有利于土壤风蚀增大的区域面积占全区面积的43.53%。其次为植被覆盖度和土壤湿度，对土壤风蚀变化均主要呈负贡献。其中，植被覆盖度变化有利于土壤风蚀减小的区域面积占区域总面积的89.53%，土壤湿度变化有利于土壤风蚀减小的区域面积占区域总面积的88.47%。植被覆盖度和土壤湿度对土壤风蚀变化贡献率集中在-10%~0%，贡献率平均值分别为-6.00%和-4.57%，说明二者变化均有利于土壤风蚀的减弱，但影响程度相对有限。与上述因子不同，土地利用变化对土壤风蚀变化总体呈正贡献，贡献率平均值为1.39%，说明土地利用变化有利于土壤风蚀的增强，且影响十分有限。

风速对土壤风蚀变化的贡献呈现出明显的空间异质性。其中，强负贡献区（贡献率<-80%）主要分布在浑善达克沙地、阴山北部、乌兰布和沙漠、腾格里沙漠、巴丹吉林沙漠、柴达木盆地沙漠、塔克拉玛干沙漠东部、西部和中部部分地区。强正贡献区（贡献率>80%）主要分布在锡林郭勒、毛乌素沙地、河西走廊、塔克拉玛干沙漠中部部分地区和柴达木盆地沙漠南部。植被覆盖度贡献率普遍较低。强负贡献区（贡献率<-50%）集中在水热条件相对较好的科尔沁沙地西部、毛乌素沙地、库布齐沙漠东部、黄土高原部分地区和塔克拉玛干沙漠与柴达木盆地的绿洲区，说明植被恢复有利于这些区域土壤风蚀强度的降低。植被覆盖度正贡献区集中分布在浑善达克沙地和科尔沁沙地的部分地区，说明植被退化造成这些区域土壤风蚀的加重。土壤湿度和土地利用的贡献率普遍较低，强贡献区域零散分布在核心攻坚区内。

本研究计算的2000—2024年核心攻坚区多年平均土壤风蚀模数为39.47 t·hm^-2·a^-1，这与利用WEPS模型计算的河西走廊-塔克拉玛干沙漠边缘阻击战核心区（35.33 t·hm^-2·a^-1），利用RWEQ模型计算的内蒙古自治区（28.57 t·hm^-2·a^-1）、毛乌素沙地（5.4~26.5 t·hm^-2·a^-1）、库布齐沙漠（44.3~102.1 t·hm^-2·a^-1），改进的第一次全国水利普查风蚀模型计算的京津风沙源（15.25~22.38 t·hm^-2·a^-1）和河西走廊（27.74 t·hm^-2·a^-1）风蚀模数接近（表2）。然而，本研究结果明显高于利用RWEQ模型计算的“三北”工程区（3.92 t·hm^-2·a^-1）和塔河流域（10.71 t·hm^-2·a^-1）的风蚀模数，这可能有两方面原因：①RWEQ模型是基于美国农田松散土壤开发的风蚀模型。该模型适用于类似农田松散土壤的地表，应用在草地、戈壁等紧实土壤地表时可能会产生低估的现象^［18］。本研究使用的改进的第一次全国水利普查风蚀模型是中国政府行业主管部门指定并推广应用的模型^［19］，同时补充考虑了表土湿度和戈壁地表的影响，适用于中国风蚀区土壤风蚀模数估算^［6，8］。②风速是风蚀模型输入参数中的最敏感数据，风速数据时空分辨率越高，风蚀模型估算的结果越准确^［20］。本研究使用的风速数据是研究区内部及其周边335个气象站点记录的逐小时风速数据。日值数据会造成风蚀模数估算的显著降低^［20］。综上，本研究估算的风蚀模数结果是可信的，满足研究需要。

人类活动是引起土壤风蚀变化的重要因子^{［9， 16， 19］}。然而本文驱动力分析结果显示，自然因子是影响核心攻坚区土壤风蚀时空分布格局的主导因子，人类活动对区域土壤风蚀的影响十分有限。这可能主要有两方面原因：①核心攻坚区人口集中区域内气候条件好、土壤风蚀强度低的地区。大部分自然环境差、土壤风蚀强度高的地区人口密度，人类活动强度小，对土壤风蚀的影响有限。②围栏封育、造林种草、减小放牧强度等生态保护修复措施主要通过影响地表植被覆盖度影响土壤风蚀。植被覆盖度的变化有利于部分地区土壤风蚀的减弱。本文并未将植被覆盖度中人类活动的贡献分离出来，这也是造成人类活动对土壤风蚀影响有限的重要原因。未来研究应拆分气候变化和人类活动对植被覆盖度的影响，进而更加准确评估人类活动对土壤风蚀的贡献。

2000—2024年核心攻坚区年均土壤风蚀量1.19×10⁸ t，平均土壤风蚀模数39.47 t·hm^-2·a^-1。核心攻坚区以微度和轻度风蚀为主，面积2 398 346 km²，占区域总面积的79.40%。区内强烈及其以上等级风蚀区面积346 352 km²，占区域总面积的11.47%，主要分布在塔克拉玛干沙漠、库姆塔格沙漠、柴达木沙漠、巴丹吉林沙漠、腾格里沙漠、乌兰布和沙漠、浑善达克沙地西部等区域。三大分区中，河西走廊-塔克拉玛干沙漠边缘阻击战工程区风蚀强度最大，属于强烈风蚀水平；科尔沁和浑善达克沙地歼灭战工程区风蚀强度最小，平均土壤风蚀模数为16.87 t·hm^-2·a^-1。风速、沙粒含量和降水量等自然因子是控制土壤风蚀空间分布的主要因子，人类活动因子的影响十分有限。

2000—2024年核心攻坚区土壤风蚀模数总体呈显著减小趋势（P<0.01），其中2000—2017年土壤风蚀模数呈显著减小趋势，2017—2024年土壤风蚀模数呈显著增大趋势。不同分区土壤风蚀模数变化趋势不同。空间上，土壤风蚀模数呈减小趋势的区域面积占核心攻坚区总面积的49.15%，土壤风蚀模数呈增大趋势的区域面积占核心攻坚区总面积的23.57%。风速是引起区域土壤风蚀变化的主导因子，植被覆盖度是引起部分水热条件较好地区土壤风蚀变化的重要因子。