河西走廊是中国“两屏三带”中的北方防沙带的重要组成部分，对北方生态安全屏障建设至关重要。土壤风蚀作为干旱、半干旱区最广泛的地表过程，以及土地沙漠化和沙尘暴的首要环节，在河西走廊地区尤为典型，使得河西走廊和相邻的内蒙古阿拉善盟地区成为中国沙尘暴策源区。风蚀模型是定量评价土壤风蚀强度，进而为区域风蚀防治提供决策依据的重要手段。常用的风蚀模型主要有美国农业部组织开发的修正土壤风蚀方程（Revised wind erosion equation，RWEQ）^［1］和风蚀预报系统（Wind erosion prediction system，WEPS）^［2］，以及Shao等^［3］提出的风蚀评价模型（Wind erosion assessment model，WEAM）。其中，RWEQ和WEPS是以农田风蚀为研究对象的经验模型，WEAM采用的粒度分布函数以疏松表土为前提，因此这些模型主要应用于农田等疏松表土的风蚀预报，针对草地、戈壁等非耕作地表的适用性和结果的可靠性尚待检验。尽管如此，这些模型仍被直接应用于估算其他类型地表风蚀量^［4-6］。最近开始有学者尝试对RWEQ进行扩展，以更好地应用于草地风蚀量估算^［7］。

第一次全国水利普查风蚀模型由农田、草灌地、沙地3个风蚀方程组成^［7］，该模型作为纯经验模型，虽然缺乏力学基础，难以获得国际同行的认可，但其形式简洁，而且可直接应用于中国北方各类主要风蚀地表，因此得到中国行业部门认可和采纳，成为全国风蚀监测中被水利部规定为唯一使用的模型。然而，该模型计算结果显示，西部广大戈壁区风蚀强度普遍达到极强烈甚至剧烈水平^［8］。按照水利部《土壤侵蚀分类分级标准》（SLl90—2007），极强烈风蚀对应的风蚀模数为8 000~15 000 t·km^-2·a^-1（对应风蚀厚度为50~100 mm·a^-1），剧烈风蚀对应的风蚀模数大于15 000 t·km^-2·a^-1（风蚀厚度大于100 mm·a^-1）^［9］。对于砾石覆盖度普遍较高的戈壁，尽管其地质历史时期曾经历了强烈的风蚀，但由于长期的地表风蚀粗化，早已形成了砾石保护层，因此极强烈和剧烈风蚀对现代戈壁来说很难达到。邢瑜等^［10］采用WEPS计算了河西走廊-塔克拉玛干沙漠土壤风蚀时空变化，结果同样显示戈壁区土壤风蚀普遍达到极强烈和剧烈水平。这意味着第一次全国水利普查风蚀模型和WEPS均可能显著高估了戈壁土壤风蚀量。鉴于此，本文拟通过改进和完善第一次全国水利普查风蚀模型，使其能够应用于戈壁风蚀量计算，提高该模型的适用性和推广应用价值。同时，利用改进和完善的第一次全国水利普查风蚀模型，评估河西走廊土壤风蚀空间分布和年际变化，为客观认识河西走廊土壤风蚀现状和土壤风蚀防治提供依据。

河西走廊位于中国西部地区，地处黄河以西、祁连山和巴丹吉林沙漠中间的甘肃省西北部，是一个呈北西-南东走向的狭长地带，行政区划包括嘉峪关、酒泉、张掖、金昌和武威5个地级市（图1）。该区属于温带大陆性气候，由于昆仑山、祁连山、阿尔金山等高山的阻挡，来自印度洋和太平洋的湿润空气难以进入，造成降水极度匮乏，年降水量只有 40~410 mm，年蒸发量1 500~3 311 mm，气候干燥多风。区内沙漠、戈壁广布，风蚀荒漠化是最严重的生态环境问题。

第一次全国水利普查风蚀模型是在高尚玉等^［11］2008年提出的风蚀模型的基础上，补充表土湿度因子后的改进模型，包括以下3个子模型。

耕地风蚀子模型：

式中：Q₁为耕地风蚀模数（t·km^-2·a^-1）；W为表土湿度（%）；T_j 为一年内风蚀期10 m高度处，超过临界起沙风速的第j个等级风速的累计时间（min）；z₀为地表空气动力学粗糙度；A为风速修订系数（取值为0.897）；U_j 为第j个等级风速（m·s^-1）；a₁、b₁、c₁分别为-9.208、0.018、1.955。

草灌地风蚀子模型：

式中：Q₂为草灌地风蚀模数（t·km^-2·a^-1）；V为植被盖度（%）；a₂、b₂、c₂分别为2.4869、-0.0014、-54.9472；其他变量含义同式（1）。不同植被盖度条件下草灌地临界起沙风速按照文献［11］提供的方程计算。

沙地风蚀子模型：

式中：Q₃为沙地风蚀模数（t·km^-2·a^-1）；a₃、b₃、c₃分别为6.1689、-0.0743、-27.9631；其他变量含义同式（1）。不同植被盖度条件下沙地临界起沙风速按照文献［11］提供的方程计算。

上述模型组没有包含戈壁地表风蚀量的计算方法。本文基于以下事实，在上述模型的基础上，补充戈壁风蚀子模型：①风蚀模数与砾石盖度负相关；②戈壁土壤具有与草灌地土壤相似的紧实度特征，其结构稳定性明显不同于沙漠和耕作土壤；③戈壁区尽管植被稀疏，但稀疏植被仍对地表风蚀具有削弱作用。据此，在草灌地风蚀子模型的基础上，补充砾石覆盖因子，作为戈壁风蚀子模型。即：

式中：Q₄为戈壁风蚀模数（t·km^-2·a^-1）；G为砾石盖度；其他变量含义同式（2）。

按以下步骤，利用上述补充完善后的风蚀模型，按年度计算河西走廊土壤风蚀模数。

确定风蚀地表类型。利用河西走廊2020年土地利用现状图（数据来源：资源环境科学与数据平台https：//www.resdc.cn/Default.aspx），将研究区地表类型重新划分为沙地（沙漠）、耕地、草灌地、戈壁和不可蚀地表5类。其中，草灌地为扣除沙地（沙漠）范围内林地和草地后的其余各类林地和草地，不可蚀地表包括水域、建设用地和裸岩（图1）。

计算植被盖度。河西走廊土壤风蚀主要发生在每年7、8、9月以外的其他月份，其中冬春季节风蚀最为强烈，而这一阶段植被大多处于非生长期，植被NDVI难以反映地表植被覆盖的真实情况。为准确反映风蚀期植被覆盖状况，借鉴高尚玉等^［11］的方法，选取治理区大多植被开始返青的4—5月植被NDVI（数据来源：NASA数据中心https：//ladsweb.modaps.eosdis.nasa.gov/search/），利用像元二分法计算风蚀期植被盖度^［12］。本文以4月中旬至5月中旬的植被盖度作为风蚀模型中植被盖度因子输入参数。图2为河西走廊2000—2023年平均植被盖度分布图。

统计风蚀期各等级风速累计时间。统计河西走廊各气象站点2000—2023年逐年1—6月、10—12月5 m·s^-1以上各等级风速的累计时间，对统计结果进行空间插值，获得该区历年各等级风速累计时间空间分布数据。风速数据来源于国家气象科学数据中心。以2023年为例，河西走廊各等级风速累计时间空间分布见图3。

确定戈壁地表砾石盖度：中国西北地区戈壁表面砾石盖度普遍为0.32~0.85^［13］，结合作者等的野外调查结果，河西走廊戈壁表面砾石盖度取值为0.75。

计算表土湿度。基于1—6月和10—12月土壤湿度数据（数据来源：NASA数据中心https：//ldas.gsfc.nasa.gov/gldas/gldas-get-data），计算获得历年各月平均土壤湿度。

利用上述数据和沙地（沙漠）、耕地、草灌地、戈壁风蚀子模型，利用GIS方法，按像元计算历年风蚀模数。依据《土壤侵蚀分类分级标准》（SL 190—2007）^［9］，确定土壤风蚀强度。

2000—2023年河西走廊平均风蚀模数为27.74 t·hm^-2·a^-1，整个研究区总体属于中度风蚀水平。土壤风蚀空间差异显著：东南部、中部风蚀较弱，以微度和轻度风蚀为主；东北部民勤绿洲及沙漠地区风蚀普遍强烈，以中度-极强烈风蚀为主；西部沙漠与戈壁镶嵌分布，土壤风蚀极为不均，轻度至剧烈风蚀均有分布，但以轻度风蚀为主，其他等级的风蚀斑块镶嵌分布于轻度风蚀区（图4）。除微度风蚀（占总面积的29.9%）外，全区轻度风蚀面积最大（134 395 km²），其次为剧烈风蚀（14 840 km²），强烈风蚀面积最小（7 354 km²）。轻度、中度、强烈、极强烈、剧烈风蚀面积占研究区总面积的比例依次为54.65%、3.11%、2.99%、3.27%和6.03%。

土壤风蚀强度的上述空间分布，与区域气候、植被、土壤等自然地理要素的空间异质性和土地利用方式有关。河西走廊土地利用类型多为沙漠、戈壁和草地，植被覆盖度整体呈现出从西北向东南逐渐增加的趋势，高植被盖度区域主要在肃南、民乐、山丹、天祝；中部和东部主要为草地，植被覆盖率较高；整个西部及东北部民勤为低植被盖度区域^［14-15］。高植被盖度区域风蚀模数小于2 t·hm^-2·a^-1，为微度风蚀。中、中高盖度区域（嘉峪关、高台、临泽、张掖、武威、古浪）多为轻度风蚀。区内绿洲多位于疏勒河、黑河、石羊河的冲积平原，土地利用类型主要为草地、耕地、灌木林地和人类生活用地，植被盖度较高，土壤风蚀主要为微度、轻度风蚀。整个西部为极干旱戈壁地区，地表砾石覆盖度大，虽然风力极强，但砾石覆盖显著抑制了土壤风蚀，风蚀强度主要为轻度，剧烈、极强烈、强烈风蚀主要发生在戈壁区的流沙斑块。其中，疏勒河流域绿洲区（敦煌、瓜州、玉门、肃北部分区域）以轻度、中度风蚀为主，敦煌以西百余千米外为哈顺戈壁西部、库姆塔格沙漠东端，地貌从风蚀地貌（砾漠）过渡为风积地貌（沙漠），冬春季西路冷空气经敦煌并沿“狭管”东行^［16］，在高风速条件下，极易发生风蚀，风蚀模数大于150 t·hm^-2·a^-1，剧烈风蚀面积较大且连续分布。酒泉市肃州区、张掖市甘州区、临泽、高台位于黑河流域中下游平原洪积粗砾质戈壁亚地区（植被盖度低于5%），戈壁、绿洲相间分布，绿洲边缘植被稀疏，地表类型由耕地、草地逐渐转为戈壁和沙地，轻度风蚀和中度风蚀相间分布。东部石羊河流域水源丰富处形成连片绿洲，如古浪绿洲、武威绿洲、永昌绿洲、民勤绿洲，风蚀模数小于25 t·hm^-2·a^-1，绿洲沙漠过渡带山前洪积砂砾戈壁与高原剥蚀砾质戈壁广泛发育^［17］，冬春季河西走廊东部西北路冷空气吹向腾格里沙漠和巴丹吉林沙漠^［16］，对民勤造成极大影响，区域内风速高，降水相对较少^［6］，极强烈、强烈、中度风蚀遍布整个区域。表1为河西走廊各土地利用类型的平均风蚀模数统计。

河西走廊各市各等级风蚀强度面积统计见表2。①酒泉市位于河西走廊西部，是河西走廊风蚀面积最大、强度最大的区域。全市多年平均风蚀模数为35.65 t·hm^-2·a^-1。其中轻度风蚀土地面积最大，占比60.3%；中度、强烈风蚀土地面积零星分布；存在2.9%的极强烈风蚀土地和8.9%的剧烈风蚀土地。②嘉峪关市区域面积较小，土地利用类型多为耕地、建设用地和未利用土地，中西部地区地表多为砾石覆盖，因此土壤风蚀相对较弱。全市多年平均风蚀模数为3 t·hm^-2·a^-1，主要为微度风蚀土地，占嘉峪关区域面积的63.0%，其次为轻度风蚀土地，占比35.3%；其他强度风蚀土地面积很小，中度、强烈、极强烈风蚀土地均有分布，但面积合计只有20 km²，只占全市总面积的1.6%。③张掖市土地利用类型多为森林和耕地，相对于河西走廊西部地区，年降水量相对较高，一定程度上抑制了土壤风蚀的发生。全市多年平均风蚀模数为3.75 t·hm^-2·a^-1，微度风蚀和轻度风蚀面积占比高达96.5%，存在3.3%的中度风蚀土地和0.2%的强烈风蚀土地，无极强烈风蚀土地和剧烈风蚀土地。④金昌市面积较小，北部与巴丹吉林沙漠、腾格里沙漠相接，风蚀相对强烈，南部绿洲地区植被覆盖良好，风蚀较弱。全市多年平均风蚀模数为9.15 t·hm^-2·a^-1，土壤风蚀多为轻度风蚀，占比52.2%；其次为轻度风蚀，占40.7%；中度风蚀、强烈风蚀、极强烈风蚀也有一定比例，无剧烈风蚀。⑤武威市北部属于腾格里沙漠和巴丹吉林沙漠的边缘地带，东北部的民勤是河西走廊区域沙尘暴暴发频次最多的地区^［18］，土壤风蚀模数长时间处于高水平。全市多年平均风蚀模数为21.95 t·hm^-2·a^-1，轻度风蚀面积最大，占比45.1%；微度风蚀土地面积次之，占比27.9%；中度风蚀和极强烈风蚀也有一定占比，但都不超过10%；无剧烈风蚀。

2000—2023年河西走廊风蚀模数总体呈波动下降趋势（图5A）。2010年是风蚀最强烈的年份，2001年次之，这两个年份全区风蚀模数分别高达46.52 t·hm^-2·a^-1和38.52 t·hm^-2·a^-1，比这一时期的平均风蚀模数分别高出67.7%和38.9%。与之相对应，2001年和2010年也是该区2000年以来沙尘暴最强烈的两个年份。2001年河西走廊地区共发生沙尘暴20~30次，其中强沙尘暴5~10次^［19］；2010年河西走廊地区共发生沙尘暴10~20次，其中强沙尘暴3~6次^［20］，分别以2001年4月6—8日的特强沙尘暴和2010年4月24日特强沙尘暴最具代表性，被认为是21世纪初中国北方最严重的两次沙尘天气事件。2019年是2000年以来河西走廊风蚀最弱的年份，全区平均风蚀模数只有15.57 t·hm^-2·a^-1，仅为风蚀最强年份（2010年）风蚀模数的三分之一。2020年以来风蚀强度有所回升，2023年全区风蚀模数接近多年平均值。

河西走廊风蚀模数年际变化与2000年以来的气候和植被覆盖变化密切相关。研究显示，2000年以来研究区年平均气温呈显著上升趋势，年降水量总体呈增加趋势^［21］。植被覆盖年际变化相对复杂，其中，敦煌、肃北、山丹、民勤、武威植被覆盖无显著变化，肃北、民勤、武威、张掖、敦煌、古浪少部分区域植被发生显著、极显著退化，但显著改善和及显著改善区域遍及整个河西走廊^［15］，总体表现为裸露土壤面积减少、植被覆盖逐渐变好的趋势^［22-23］。本文计算结果同样显示研究区植被总盖度呈波动增大的趋势（图5A），其中，2001年是这一时期植被覆盖最差的年份。这是风蚀模数减小总趋势的主要原因。风速方面，研究区平均风速总体呈增大趋势（图5A），而2010年风速异常增大，之后风速波动性降低，2020年以来风速有增大趋势。显然，2001年风蚀强烈的主要原因是植被覆盖极差，2010年风蚀显著加剧的主要原因是风速的异常增大，2020年后风蚀模数回升的主要原因也在于风速增大。可见，河西走廊土壤风蚀仍在很大程度上受控于气候条件的变化，这是干旱、极端干旱区生态环境脆弱敏感、植被恢复困难的必然结果。

研究区各市土壤风蚀年际变化趋势总体相近，细节有所不同（图5B）。其中，酒泉市土壤风蚀年际变化最为剧烈，2010年的峰值与2020年的最小值相差近2.6倍。嘉峪关市、张掖市风蚀模数均处于0~10 t·hm^-2·a^-1，年际波动幅度最小，风蚀模数在2016年出现高峰。金昌市2005年风蚀最弱，2005—2015年除2010年风蚀明显加剧外其余年份变化不大，2016年之后风蚀模数逐渐回升，2021年达到2000年以来的峰值。武威市2015年之前土壤风蚀年际变化与金昌市类似，2019年风蚀模数减小到最低值之后迅速回升，同样在2021年达到峰值。该市2000—2023年风蚀模数变化曲线与沙尘暴日数年际变化曲线十分相似^［24］。

第一次全国水利普查风蚀模型是中国第一个由政府行业主管部门指定并推广应用于全国风蚀监测的风蚀模型^［8］。然而，该模型所包含的方程只适用于草灌地、耕地和沙地（沙漠），不包含针对戈壁地表的风蚀方程。上述问题导致两方面的后果，一是国家决策部门与学术界计算的土壤风蚀数据脱节，二是戈壁风蚀量计算结果与实际不符。本文以提高第一次全国水利普查风蚀模型适用性和推广应用价值为基本出发点，在此基础上，研究河西走廊土壤风蚀时空变化规律。从补充的戈壁风蚀方程计算结果看，酒泉市北部戈壁区多年平均风蚀模数约为33.6 t·hm^-2·a^-1，按照戈壁地区表层土壤平均容重约为2.0 t·m^-3计，折合风蚀厚度约为1.68 mm·a^-1。考虑到该区域包含流沙斑块，我们认为这样的平均风蚀强度是合理的。若扣除其中的流沙斑块，戈壁地表平均风蚀模数按5.54 t·hm^-2·a^-1，年风蚀厚度约0.27 mm。而利用原模型计算得到的第一次全国水利普查风蚀强度图^［8］，以及邢瑜等^［10］采用WEPS计算得到的河西走廊-塔克拉玛干沙漠边缘阻击战核心区风蚀强度图则显示，河西走廊境内酒泉市北部戈壁区全部由强烈、极强烈、剧烈风蚀斑块所组成。由于无法从这两个文献获取该区域范围内的风蚀模数原始数据，若取强烈风蚀下限即50 t·hm^-2·a^-1，则意味着采用第一次全国水利普查原始风蚀模型和WEPS计算的酒泉北部戈壁区年风蚀厚度最小可达2.5 mm，其余极强烈、剧烈风蚀斑块风蚀厚度更大，平均风蚀厚度将远远高于本文计算结果。对于砾石覆盖度普遍较高的戈壁，这样的风蚀厚度显然有所高估。

目前尚未发现利用RWEQ、WEAM等模型计算的河西走廊土壤风蚀模数。由于这些模型从机理上讲只能适用于疏松表土的风蚀预报，对于戈壁区砾石覆盖下的紧实结构土壤，模型模拟结果也将不可避免地高估风蚀量。

鉴于此，我们建议将本文提出的戈壁风蚀方程补充到第一次全国水利普查风蚀模型中，形成包含耕地、草灌地、沙地（沙漠）、戈壁等4个风蚀方程的模型组，作为改进的第一次全国水利普查风蚀模型，在中国进行推广应用，以提高风蚀评估精度并消除政府与学术界计算的土壤风蚀数据差。

2000—2023年，河西走廊平均风蚀模数为27.74 t·hm^-2·a^-1，土壤风蚀空间差异显著。东南部和中部以微度和轻度风蚀为主；东北部民勤绿洲及沙漠地区风蚀以中度-极强烈风蚀为主；西部沙漠与戈壁镶嵌分布，土壤风蚀极度不均，轻度至剧烈风蚀均有分布，但以轻度风蚀为主。2000—2023年风蚀模数总体呈波动下降趋势，2020年以来风蚀强度有所回升。2010年和2001年是这一时期内风蚀最强烈同时也是沙尘暴最强烈的两个年份，体现了两者之间的因果关系。