河西走廊-塔克拉玛干沙漠边缘阻击战核心区位于中国西北干旱地区，是“三北”工程攻坚战和北方生态安全屏障建设的主战场。作为中国最强烈的风沙活动线和最严重的风沙灾害区，该地风蚀起尘、起沙问题对中国西北乃至全国的生态安全具有深远影响^［1-2］。同时，这一区域不仅见证了古代“丝绸之路”的辉煌历史，还是现代“欧亚大陆桥”和现代“一带一路”倡议中的生态安全关键枢纽，承载着促进区域乃至国际合作与交流的重要功能^［3］。然而，针对该区域风蚀和起尘量的研究相对匮乏。

风蚀是漫长的地质过程^［4］，深刻影响着生物圈的多个层面。风蚀带走大量富含土壤养分的地表颗粒物，导致土壤退化、生产力下降^［5-6］。此外，风蚀形成的细颗粒物是大气气溶胶的主要来源，为沙尘暴事件提供了丰富的物质来源^［7］，造成大气环境污染，危害交通、通信和水利设施等，严重影响人类健康和经济发展。尤其是包括河西走廊-塔克拉玛干沙漠边缘阻击战核心区的东亚沙尘源区，其风蚀影响范围之广、程度之深，不容忽视^［8］。因此，准确估算地表风蚀量和起尘量是防治土壤风蚀和区域荒漠化的重要需求。为此，国内外学者基于大量野外观测^［9］、室内风洞及可移动风洞试验工作^［10］，先后提出了多种土壤风蚀模型，可开展区域尺度的评估^［11］。常用的风蚀模型主要包括风蚀方程（Wind Erosion Equation， WEQ）^［12］、修正风蚀方程（Revised Wind Erosion Equation， RWEQ）^［13］和风蚀预报系统（Wind Erosion Prediction System， WEPS）^［4］等，其中WEPS能够以每小时步长有效模拟不同气候、地形和植被覆盖条件下的土壤风蚀动态过程，尤其是可以给出PM₁₀结果，是目前公认的较为先进和全面的风蚀模型。

全球气候变化加剧为研究区的生态环境带来了新的挑战^［14］。土壤风蚀对气候变化高度敏感^［15-16］，例如气温上升、降水模式变化及极端干旱事件的增多^［17］，不仅会增加土壤水分蒸发，缩短冻结期，加剧土壤风蚀，还影响植被的生长分布与防风固沙能力。在此背景下，“三北”防护林等重点生态建设工程在应对气候变化挑战、维持研究区生态环境中的作用更为重要。

本研究采用土壤风蚀预报系统WEPS的核心算法，使用气象、土壤和植被等多源地理数据分析了2000—2023年河西走廊-塔克拉玛干沙漠边缘阻击战核心区的风蚀量及PM₁₀变化趋势，探讨其主要影响因素，以期深化对该研究区风沙活动及其影响的认识，为“三北”重点生态建设工程布局与优化、提升防风固沙生态系统服务功能提供科学依据和决策参考。

河西走廊-塔克拉玛干沙漠边缘阻击战核心区主要包括河西走廊地区（甘肃省酒泉市、嘉峪关市、武威市、金昌市、张掖市、白银市西部）和新疆五地州（阿克苏地区、喀什地区、克孜勒苏柯尔克孜自治州、和田地区、巴音郭楞蒙古自治州），总面积为80.67 万km^2［18］（图1）。核心区内沙漠、戈壁广布，是中国北方强度最大的沙尘暴源地。其中，塔克拉玛干沙漠是中国第一大沙漠、世界第二大流动沙漠。河西走廊地区东起乌鞘岭，西接吐哈盆地，南依祁连山，北偎腾格里、巴丹吉林两大沙漠，不仅是中国北方地区主要沙尘源区，也是风沙东移南下的主要通道，北方冷空气南下在这一区域易形成“狭管效应”。

在本研究中，我们使用了2000—2023年的气象数据、土壤数据和植被数据（表1）估算风蚀量及PM₁₀量。风速数据来自第五代欧洲中期天气预报中心（ECMWF，https：//www.ecmwf.int/en/era5-land67）。土壤数据来自世界土壤数据库（HWSD），由于土壤的物理和化学性质在短时间内不会发生变化，因此我们假定土壤质地和化学性质在研究期间保持不变。归一化植被指数（NDVI）和增强植被指数（EVI）数据来自MOD13C1产品（https：//search.earthdata.nasa.gov/search）。叶面积指数（LAI）来自GEE平台的MODIS/061/MOD15A2H数据集，植被冠层高度来自全球冠层高度2020数据集（https：//earthengine.google.com/）。降水量数据来自全球降水气候项目（GPCP）数据集（https：//www.ncei.noaa.gov/access/metadata/landing-page/bin/iso？id=gov.noaa.ncdc：C00979）。气温数据来自GLDAS Noah Land Surface Model L4月度数据集。土壤水分数据来自GLDAS Noah Land Surface Model L4日数据集（https：//search.earthdata.nasa.gov/search）。

本研究采用WEPS中的风蚀测算方案，结合气象、土壤质地、植被特征等数据，使用ArcGIS10.2在研究区内建立5 km×5 km格网，在小时级的时间步长内分别计算每个栅格的风蚀和PM₁₀总量，并叠加各栅格结果，根据区域面积给出整体风蚀模数和PM₁₀释放速率。基本方程如下：

式中：q为单宽输沙量，kg·m^-1；u_* 为摩阻风速，m·s^-1；u_*t为临界摩阻风速，m·s^-1。可见输沙的驱动力在于摩阻风速u_* 大于临界摩阻风速u_*t^［19］。

当地表有直立植被时，需计算直立植被对摩阻风速的降低作用，因此观测地点摩阻风速由气象站摩阻风速u_*f和有无植被摩阻风速u_* 共同计算得到：

式中：u_*f为气象站摩阻风速，m·s^-1；u为高度为z时风速；z为观测风速的高度（通常情况获取的标准观测风速的高度是10 m）；z_0f为空气动力学粗糙度，取25 mm；z₀为当地空气动力学粗糙度，取地表垄作粗糙度与随机粗糙度的最大值；u_*b为观测地点无植被时的摩阻风速，m·s^-1。

观测地点有植被覆盖的摩阻风速u_* 依据有效生物量拖曳系数大小，采用下列公式来进行计算。

式中：u_*v为植被冠层上部摩阻风速，m·s^-1；WZ_ov为直立植被冠层空气动力学粗糙度，m；BR_cd 为植被有效拖曳力系数；BR_lai为直立植被叶面积指数，m²·m^-1；BR_sai为植被茎面积指数，是茎面积与水平地表面积的比值，m²·m^-1；BZ为直立植被高度，m。

临界摩阻风速u_*t由光滑平坦地表临界起动摩阻风速u_*tb，地表有倒放植物引起的临界起动摩阻风速增加量u_*tve和含水率引起的临界摩阻风速增加量u_*twc计算得到：

式中：SF_cv为地表覆盖土块、结皮或石块的面积百分比（不起沙面积百分比）；SFC_cv为免受排放影响的土壤表面积的变化率；BFF_cv为倒放植被覆盖面积；M_wc为地表含水率；M_1.5wc为1.5 MPa地表含水率^［20］。

沙尘物质中，能够长距离传输的PM₁₀含量可作为起尘量的有效估计。WEPS中PM₁₀计算方法如下：

该式为二维模拟区域的PM₁₀的质量守恒方程。式中：C10¯为PM₁₀颗粒物平均浓度，mg·m^-3；H10¯为格网内PM₁₀能达到的高度，m；x，y分别为格网的长度和宽度，即平行于边界垂直方向的水平距离，m； t为时间，s；q₁₀为PM₁₀单位时间和宽度内的排放量，kg·m^-1·s^-1；G_10en为松散土壤释放的PM₁₀垂直方向净通量，kg·m^-2·s^-1；G_10an为土块和地壳磨损而产生的PM₁₀垂直方向净通量，kg·m^-2·s^-1；G_10ss为蠕移和跃移破裂而产生的PM₁₀垂直方向净通量，kg·m^-2·s^-1；SF_10en为从地表排放的悬移物质中PM₁₀的占比；SF_SSen为可流失松散土壤中悬移颗粒比例；C_en为排放系数；SF₁为直径小于0.01 mm的土壤占比；SF₁₀则为直径小于0.10 mm的土壤占比；SF_10an为地表释放的磨损悬移物质中PM₁₀的占比；SF_cla为土壤质地中黏土含量；SF_10bk为蠕移和跃移撞击破碎所产生的悬移物质中PM₁₀的占比；F_ani 表示第i个产生跃移磨损的土壤占比；C_ani 为第i个产生跃移磨损相应系数；SF_si为土壤质地中粉砂土含量；其他参数意义同上。

本研究采用Theil-Sen Median趋势分析和Mann-Kendall显著性检验相结合的方法，对研究区2000—2023年土壤风蚀和各气象因素的年际变化趋势进行分析。该方法能够有效应对数据收集过程中的测量误差和不规则性，在地理研究中得到了广泛应用^［21］。Theil-Sen Median 斜率β的计算公式为：

式中：β为相关变量线性变化趋势的幅度，β0表示时间序列内呈趋势上升，反之下降；x_j 和x_i 分别为第j次和第i次（ji）的相应数据值。

Mann-Kendall显著性检验过程为：

式中： S表示所有正差异与负差异数量的差值；n为时间序列的总长度；x_i 和x_j 分别是时间序列x的第i和j个数值；sgn（）为统计符号；m为序列中结（重复出现的数据组）的个数；t_i 为第i组重复数据组中的重复数据个数。根据正态分布统计表，在给定显著性水平（α=0.05）下，当|Z|≤Z_1-α/2时，接受原假设，即趋势不显著；若|Z|Z_1-α/2，则拒绝原假设，即认为趋势显著。

之后，通过相关性分析方法，计算了土壤湿度（SM）、风速（WS）和归一化植被指数（NDVI）这3个主要影响因子与风蚀之间的相关系数^［22］。线性相关的计算公式如下：

式中：R_xy 为两变量x，y的相关系数；x¯，y¯分别为变量x，y的均值；n为样本量。

为了剥离其他变量的影响，仅关注风蚀与单个影响因子之间的关系，进一步计算了偏相关系数。偏相关系数的计算公式如下：

式中：r_12⋅34表示在控制变量x₃和x₄的影响后，变量x₁与x₂之间的偏相关系数；r_12·3，r_14·3，r_24·3分别为排除所计算的影响因素x₃后，x₁与x₂、x₁与x₄、x₂ 与x₄之间的偏相关系数。

最后，将上述3个主要影响因子，在每个格点上建立多元线性回归模型，并对比得到标准化系数绝对值，回归系数绝对值最大的因子即为格点中的主导影响因子，模型可以表示为：

式中：b₀是回归常数项；b_i （i=1，2，3）是使用最小二乘法计算的标准化回归系数。分析通过95%的显著性检验，每个格点上的拟合方程对应的判定系数计算公式为：

式中：WE_i 为第i个风蚀量值；WE_ave是风蚀量平均值；WE_MLR是根据回归模型计算得到的风蚀量结果。

2000—2023年，研究区风蚀总量为6.88×10¹⁰ t，年风蚀模数为2 232~5 116 t·km^-2，多年平均风蚀模数为3 553 t·km^-2（图2）。按照水利部《土壤侵蚀分类分级标准》（SL190-2007）^［23］，将研究区风蚀强度分为微度、轻度、中度、强烈、极强烈和剧烈6个等级，年风蚀模数大于200 t·km^-2、强度达到轻度以上的区域面积为6.00×10⁵ km²，约占研究区总面积的69%，强烈及以上风蚀强度面积约占研究区总面积的21%。风蚀最强烈的区域主要分布在塔克拉玛干沙漠东南缘和罗布泊边缘地带的且末县、若羌县以及河西地区的敦煌市、瓜州县、肃北蒙古族自治县、玉门市、金塔县和民勤县（图3）。在这些地区，多年平均风蚀模数超过5 000 t·km^-2。

研究区各季节风蚀模数空间分布与年平均值相似，但差异明显。春季（3—5月）风蚀模数最大（1 652 t·km^-2），且影响范围最广，占年内风蚀总量的47%。夏季（6—8月）风蚀模数较春季有所降低（1 020 t·km^-2），对年内风蚀量贡献为29%，强风蚀区域主要在塔克拉玛干沙漠东缘若羌县和河西走廊地区瓜州县和玉门市。秋季（9—11月）和冬季（12月至次年2月）风蚀模数较低（分别为507 t·km^-2和358 t·km^-2），且范围较小，对年内风蚀量贡献分别为14%和10%，且集中在罗布泊边缘地带和河西地区肃北蒙古族自治县与玉门市（图4）。

2000—2023年研究区平均风蚀模数变化趋势为减弱—增强—减弱（图5）。21世纪初期，尤其是2001年风蚀模数较大，达到4 424 t·km^-2，之后呈下降趋势，2005年风蚀模数降到最低值，为2 232 t·km^-2。随后研究区风蚀模数呈上升趋势，2010年风蚀模数为最大值，达到5 116 t·km^-2。此后研究区风蚀模数呈下降趋势，并且均未超过初期2001年以及2010年的水平。总体来看，2000—2023年，研究区风蚀程度减轻，强烈及以上程度风蚀区域面积减少了2.15×10⁴ km²。空间上，风蚀减少区域占主导地位，约占研究区总面积的48%，年代下降速率达到41 t·km^-2，主要分布在塔克拉玛干沙漠中部和东南缘，以及河西走廊中部、北部地区，仅在塔克拉玛干沙漠西北部阿克苏市、喀什市，东北部尉犁县和河西地区肃北蒙古族自治县北部等局部地区风蚀模数呈现增加趋势（图6）。

从研究区风蚀模数季节变化来看，2000—2023年，风蚀季节变化特征明显（图7A），风蚀模数减少最快的是秋季（12 t·km^-2），其次是冬季（7 t·km^-2）和夏季（3 t·km^-2），而春季风蚀模数呈现增加趋势（8 t·km^-2）。从月平均风蚀模数变化来看（图7B），研究区4月平均风蚀模数最大（686 t·km^-2），1月平均风蚀模数最小（94 t·km^-2）。其中4月平均风蚀模数增加趋势最大（38 t·km^-2），而6月平均风蚀模数下降趋势明显（48 t·km^-2）。

在所有被吹蚀物质中，尘的传输距离要长得多，特别是直径小于10 μm的颗粒（PM₁₀），它们能在大气中停留足够长的时间，在顺风方向上传输数千千米^［24］。研究区2000—2023年平均PM₁₀释放量为3.11×10⁷ t，平均释放速率为38.53 t·km^-2·a^-1。其中，春季PM₁₀释放量是其他季节的1.6~4.3倍。研究区PM₁₀与风蚀量高值区空间分布基本一致，主要位于塔克拉玛干沙漠东缘和河西走廊北部戈壁地区，而罗布泊周边盐碱地地区较低（图8）。

2000—2023年，研究区风蚀模数总体上呈现出下降趋势，风速、植被覆盖度和土壤湿度是其主要影响因素。风蚀量在研究区与归一化植被指数和土壤湿度负相关（负相关区域面积占比分别为69%、73%），与风速正相关（正相关区域面积占比为91%）。这种相关关系是不同因子对风蚀起尘影响的反映，大风是土壤风蚀最直接的动力来源^［25］，密集的植被覆盖和较高的土壤湿度则会增加临界摩擦速度，防止风蚀。

图2和图9显示，风蚀高值区集中在风速高、植被稀少、土壤干燥的区域。为量化各因素的独立贡献，我们进行了偏相关分析以剔除其他因素的干扰。结果表明，风速是影响风蚀量的关键因素，尤其是在阿克苏地区和罗布泊、若羌地区以及河西走廊肃北、瓜州、敦煌和民勤等地，风速与风蚀量的偏相关系数超过0.9。而归一化植被指数和土壤湿度与风蚀量的显著负相关仅出现在塔克拉玛干沙漠南部和河西走廊北部的部分区域（图10）。

多元线性回归分析显示（图11A），风速是风蚀高值区的主导因素，特别是在塔克拉玛干沙漠东南缘和河西走廊北部，而土壤湿度和归一化植被指数主要影响研究区西部的风蚀低值区。在上述风速主导的区域，回归方程拟合效果较好（R²0.8，图11B）。

从2000—2023年的影响因素年际变化趋势来看，研究区风速总体呈下降趋势（减少区域面积占比为63%），其中中西部区域下降趋势明显，仅南疆和田地区、阿克苏地区和河西走廊肃北县局部区域风速上升（图12A）。与此同时，归一化植被指数、降水量和温度均呈现增加趋势（增加区域面积占比分别为94%、98%和99%，图12B~D）。这些变化会改善区域生态状况并减弱风蚀^［26］，这与上述影响因子与风蚀量的相关性分析结果一致。但我们需要注意的是，气温增加会缩短土壤的冻结期，导致其抗侵蚀性降低，从而增加春季风蚀风险。此外，温度升高加剧土壤蒸发，导致部分区域土壤湿度下降（减少区域面积占比54.38%），主要在塔克拉玛干沙漠北缘和河西走廊中部地区，干燥松散地表可能会加剧风蚀风险（图12E）。

本文计算的2000—2023年研究区多年平均风蚀模数（3 553 t·km^-2）接近于利用RWEQ风蚀模型计算的1990—2015年中国干旱区多年平均风蚀模数（3 222 t·km^-2）^［27］，并与邻近内蒙古地区基于风蚀经验公式计算的2000—2017年平均风蚀模数（3 192 t·km^-2）相符，且变化均呈下降趋势^［28］。此外，本研究区2001—2014年PM₁₀多年平均释放量（3.29 ×10⁷ t）在Wang等^［29］使用沙尘排放模型所计算的中国北方同一时期起尘量（1.10 ×10⁸ t）范围内。Liu等^［30］计算的2016年南疆地区PM₁₀释放总量为1.2×10⁷ t，这与本文计算的塔克拉玛干沙漠范围内多年PM₁₀释放均值（1.47×10⁷ t）结果基本一致。综上，本文得到的研究区风蚀量和PM₁₀动态变化是合理的。同时，由于本研究基于小时时间尺度，能够更精确地捕捉风蚀和PM₁₀释放过程及其总量。未来，随着观测数据的不断积累，我们将继续优化模型方案，并提出更适用于中国北方荒漠区的计算方法。

研究区内广袤的沙漠、戈壁及沙漠化土地，为风蚀和起尘活动提供了物质基础。此外，其气候四季分明，进一步影响了风蚀活动的强度变化。春季是风蚀和起尘的高发期，尤其是4月，因大气不稳定、气旋活动频繁，风速与大风日数显著增加，成为风蚀与起尘的高发期。夏秋季风蚀活动虽有所减弱，但仍不可忽视。而冬季土壤冻结，土壤易蚀性降低^［31］，但寒冷干燥的气流仍会对区域环境产生影响，增加风蚀与起尘威胁。

从风蚀变化的驱动因素分析可见，风速是影响研究区风蚀变化的主要气候因素，其动态变化直接影响着风蚀强度的改变。同时，区域风蚀及起尘活动受到特定风沙路径的显著影响。研究区大风入侵的主要路径表现为盛行风在塔克拉玛干沙漠的“西进”和“东灌”两条气流路径。“西进”沿昆仑山北麓，影响南疆盆地中西部地区；而“东灌”则是冷空气由罗布泊和库鲁克塔格沙漠至若羌和喀什地区沿线，直接影响塔克拉玛干沙漠腹地^［32-33］。河西走廊区域则受“西路”和“西北路”冷空气的影响更为突出。河西走廊为东西走向，南北介于祁连山和马鬃山之间，其特殊地理位置使“西路”冷空气由南疆东移，穿过塔里木盆地并沿河西走廊东行，显著影响瓜州、玉门和敦煌等沿线地区。“西北路”冷空气则从西西伯利亚和蒙古国一带吹向巴丹吉林沙漠和腾格里沙漠，进而影响河西走廊的金昌、武威和民勤等地^［34-36］。这些路径不仅决定了风沙的传输方向，还加剧了特定区域的风蚀与起尘强度，如塔克拉玛干沙漠、腾格里沙漠、巴丹吉林沙漠边缘及其绿洲地带，成为风蚀与起尘活动的重点区域。

地表状况方面，研究区的植被覆盖度变化总体呈上升趋势，这一趋势主要受到自然降水增加和生态建设工程的双重驱动^［37］。近几十年来，生态建设工程在研究区取得显著成效^［38-39］，逆转了沙漠化趋势^［40］。然而，塔克拉玛干沙漠的西缘和北缘，归一化植被指数呈下降趋势，风蚀活动加剧。该区域的植被退化可能受人类活动等多种因素影响。在塔克拉玛干沙漠北缘，不合理的耕作和灌溉导致荒漠植被减少，耕地荒漠化问题突出^［2］。此外，由于上游水库建设和用水量增加，河流的洪泛作用减弱，植被可利用水资源减少，胡杨林逐渐衰退，加剧了区域风蚀与起尘现象。

研究区风蚀及起尘的分布和变化是气候背景与地表状况等多重因素共同作用的结果。在干旱荒漠地带，地表的稳定性通常仅在风速减弱的情况下得以维持。面对强风的侵蚀，往往需要依赖人工干预和植被种植等措施^［41］。因此，制定防治措施时，必须全面考虑这些复杂因素的综合影响。鉴于植被恢复受到研究区水资源有限的制约，科学的防治策略应以建立防护林体系为核心，重点关注风蚀和起尘严重的区域，如塔克拉玛干沙漠、河西走廊中部，以及腾格里沙漠与巴丹吉林沙漠边缘的绿洲地区。在保护荒漠地表自然状态的同时，合理调控人类活动的强度，减少人为干预可能带来的风蚀加剧^［42］。此外，防治措施应根据季节变化进行动态调整。春季应加强防风固沙的工作，通过植被和防护林减缓风速，降低土壤侵蚀风险。夏季则需注重土壤保护，以保持土壤湿度和稳定性。冬季应减少地表扰动，防止扬尘增加。通过实施这些针对性的防治措施，可以实现风蚀与起尘的精准控制，保障区域生态环境的持续稳定与健康发展。

2000—2023年，河西走廊-塔克拉玛干沙漠边缘阻击战核心区的土壤风蚀总风蚀量达6.88×10¹⁰ t，多年平均风蚀模数为3 553 t·km²。该区域的风蚀强度在不同区域和季节存在差异，强烈风蚀主要在塔克拉玛干沙漠东南缘和罗布泊边缘，以及河西走廊的敦煌、瓜州等地区。季节性风蚀与年风蚀空间分布相似，其中春季风蚀模数最高，秋冬季则较低。同期，研究区PM₁₀的平均释放量为3.11×10⁷ t，其中春季释放量为其他季节的1.6~4.3倍。

21世纪以来，研究区风蚀强度经历了减弱—增强—再减弱的变化，风蚀模数在2005年降至最低值，2010年达到峰值后逐渐下降。整体上，风蚀强度呈下降趋势，强烈及以上程度的风蚀区域面积减少了2.15×10⁴ km²，减少区域主要在塔克拉玛干沙漠中部和东南缘以及河西走廊中部地区。

研究区91%的范围内风速与风蚀模数呈正相关，风速的变化直接影响风蚀强度；而植被覆盖与土壤湿度则分别在69%和73%的区域内与风蚀量呈现负相关关系。研究期间，63%的区域风速呈下降趋势，94%的区域植被覆盖增加，98%的区域降水量上升，这些变化共同推动了风蚀量的减少。然而，全域温度升高可能导致土壤冻结期缩短，从而增加冬春干旱季节的潜在风蚀风险。因此，在风蚀防治过程中，应重视区域气候变化与生态恢复对风蚀动态的复杂影响。