土壤风蚀是很多国家和地区农业土地利用中的重要问题^［1-2］，主要发生在干旱、半干旱以及部分半湿润地区^［3-5］。大约28%的全球陆地面积受土壤风蚀的影响^［6-8］。中国是世界上受风蚀危害严重的国家，沙化土地面积达172.12万km²，占国土面积的17.93%^［9］。风蚀带走大量富含土壤养分的地表颗粒物，导致土壤退化、生产力下降，区域生态环境服务功能降低，同时造成大气环境污染，危害交通、通讯和水利设施等，严重危害人类健康和经济发展^［10-15］。因此，研究土壤风蚀模数的时空变化规律对保护生态环境和生态文明建设至关重要。

地表风蚀量的准确估算是进行土壤风蚀研究的基础。为了准确估算地表风蚀量，国内外科学家做了大量研究并先后提出了多种土壤风蚀模型。目前相对成熟的风蚀模型主要包括两类：一是经验模型，主要有风蚀方程（WEQ）^［16］、修正风蚀方程（RWEQ）^［17］、风蚀预报系统（WEPS）^［18］、第一次全国水利普查模型（NWESMC）^［19］等。经验模型往往受建模区域的限制较大，在模型构建区域应用效果较好，推广应用则存在很大的局限性。二是物理模型，以风蚀评价模型（WEAM）和综合风蚀模拟系统（Integrated wind-erosion modeling system，IWEMS）为代表^［20-21］。WEAM模型以风沙流和大气尘输移的实验与理论成果为基础，构建用以估算农田风沙流及大气尘输移量的模型。该模型物理过程明确，但模型中的变量未能涵盖影响风蚀过程的各种主要因素及因素之间的相互作用^［22］。基于WEAM模型的弊端，Shao^［21］提出了综合风蚀模拟系统（IWEMS）。该模型包括大气模型、地表结构模型、风蚀过程模型、输送和沉积模型以及地表信息数据库。其理论基础是沙粒碰撞理论，认为土壤风蚀的机制主要是跃移颗粒冲击引起的土壤团粒磨蚀和溅射以及沙粒的分离和释放，而冲击效率与跃移颗粒的直径、沙粒直径以及起动风速有关。IWEMS模型将沙尘扩散和输移尺度扩大到了澳洲全境，同时基于多种不容易侵蚀物质和地表可蚀性颗粒的研究成果对临界摩阻风速的计算方式进行了修正，并且对于不同的侵蚀颗粒分别建立了独立的拖曳力模型，在粉尘排放量和风蚀量的计算上应用效果较好^［23-28］。该模型在中国应用较广，已经被应用于黄土高原、内蒙古等干旱半干旱地区的土地退化与草地生态系统管理研究中^［29-30］。该模型被用于计算中国北方地区的粉尘释放量、输沙通量和宁蒙河段的风沙入河量，并得到校准，应用比较成熟^［27-28］。

河北坝上地区地处京津风沙源上风向，是典型的风蚀区，特殊的地理位置使本区成为京津冀的重要生态屏障。由于自然环境影响和不合理的土地利用方式，生态环境日益恶化，对当地和整个京津冀地区人民的生产生活产生严重影响。许多学者在该区进行了风蚀研究，不论是在风蚀机理方面还是在风蚀模型方面，均取得了丰硕的研究成果^［31-33］。但是经济发展导致的土地利用方式变化以及全球变暖背景下出现的气候暖干趋势的影响，使这一地区的土壤风蚀危害仍不容忽视。因此，本文借鉴前人的研究结果，选择河北坝上地区为研究区，以遥感、气象、土壤等多源数据为基础，采用IWEMS模型估算该地区2000—2018年的逐月风蚀模数，分析风蚀模数的时空变化规律，为坝上地区风沙防护体系的调整与优化提供科学依据。

坝上地区位于河北省西北部，地处北方农牧交错带中段，包含尚义县、沽源县、张北县、康保县、围场满族蒙古族自治县（简称围场县）北部和丰宁满族自治县（简称丰宁县）北部共6个区域（图1），面积约1.9万km²，占河北省总面积的9.8%。该区域为东亚大陆性季风气候，是半干旱区向干旱区的过渡带，年平均气温1~2 ℃，年降水量400 mm左右，70%集中在6—9月，年蒸发量高达1 800 mm。特殊的地形和气候条件使该地区生态系统相对脆弱，自我修复能力差。主要土地利用类型为耕地和草地，土壤以栗钙土、沙质栗钙土为主，有机质含量少且土质松软，易受侵蚀风化。当地的耕作习惯为秋季或早春进行农田翻耕，从而形成大面积的裸露农田，冬春季又有强风天气，使该地区成为中国农田风蚀严重地区，造成土地沙化退化、土壤肥力损失和风沙灾害频发等一系列严重的生态问题。

基础数据主要包括土壤数据、气象数据、植被数据、DEM（数字高程模型）数据、遥感数据、土壤湿度数据和不同时期的土地利用数据。其中土壤数据来源于中国1∶100万土壤数据库，利用数据库中砂粒、粉粒和黏粒的含量，计算不同土壤粒径的粒级占比p（d_s）（粒径为d_s的沙粒在该区域自然地表中所占比例），并计算不同粒径在没有植被和土壤湿度影响下的临界摩阻起动风速（u_*t）：

式中：ρ_a为空气密度（1.29 kg·m^-3）；ρ_p为沙粒的密度（2 600 kg·m^-3）；g为重力加速度（9.8 m·s^-2）；a₁为无量纲参数；a₂为量纲参数。

气象数据来源于中国气象科学数据共享服务网（http：//data.cma.cn）提供的2000—2018年国家台站数据，均为每小时的平均风速。利用风速数据计算整个坝上地区每小时的摩阻风速u*分布情况：

式中：u_z表示气象站点风速；z为气象站高度；z₀为空气动力学粗糙度；k为卡曼系数。

遥感数据来源于USGS网站（https：//www.usgs.gov）提供的2000—2018年MOD13A2数据产品，采用像元二分法模型计算植被盖度（VC）和迎风面积指数（λ），进而计算植被对摩阻起动风速的影响：

式中：σ_r表示植被基部面积与迎风面积的比值；β_r表示作用于单株植被上的阻力与作用于地表阻力的比值；m_r为调整参数。

土壤湿度数据使用美国航空航天局戈达德空间飞行中心和美国海洋和大气局国家环境预报中心联合发布的2000—2018年土壤湿度数据（GLDAS）。土壤湿度对摩阻起动风速的影响函数为：

式中：A、b均为无量纲参数；θ_r表示风干土壤的含水率；θ表示土壤湿度，使用GLDAS数据。

土地利用数据采用的是根据2000、2005、2010、2015年Landsat遥感影像人机交互解译的数据（https：//www.resdc.cn/）。

IWEMS是由Shao^［21］提出的用于预测区域尺度的风蚀预报模型，该模型是基于过程的物理模型，采用独立起动的概念计算输沙通量，能够计算区域尺度的土壤风蚀模数。本文主要采用校准后的模型进行计算^［27］。

IWEMS有两个关键计算步骤：摩阻起动风速u_*t 的计算和跃移风蚀量的计算。摩阻起动风速主要考虑植被、土壤湿度和土壤粒径的影响：

式中：u_*t （d_s，λ，θ）表示粒径为d_s的土壤颗粒在植被和土壤湿度影响下的摩阻起动风速；f_v（λ）、 f_w（θ）分别表示植被和土壤湿度对摩阻起动风速的影响。均一地表条件下某一粒径的跃移通量计算公式为：

式中：A_c表示可蚀性土壤所占比例；c₀为Owen系数；u_*t 为摩阻风速，根据气象站风速数据计算获得。IWEMS模型采用沙粒独立起动的概念，认为每个粒级的颗粒在土壤风蚀过程中互不影响，因此自然地表的跃移通量可表示为：

式中：d₁和d₂分别表示土壤粒径的下限和上限；p（d_s）表示粒径为d_s的沙粒在该区域自然地表中所占比例。

2000—2018年坝上地区的平均年风蚀模数为4.89~39.89 t·hm^-2·a^-1，风蚀以轻度和中度侵蚀为主，大部分地区土壤风蚀程度较低（图2）。风蚀严重的地区主要分布在坝上中西部，尤其是康保县、张北县北部和丰宁县西部，其次为沽源县和尚义县的北部，其余地区风蚀模数较小，属于微度侵蚀，有些地区风蚀模数甚至为0。

2000—2018年，风蚀模数的空间分布情况有所变化，主要表现在坝上地区西部和中部。其中坝上地区西部，即张北县北部和康保县的风蚀模数，从2000年到2009年显著减小，2014年之后再次增加。而坝上地区中部，即丰宁县西部和沽源县东部，自2009年以来，风蚀模数逐渐增加，逐渐超过坝上地区西部，成为整个坝上地区风蚀最为强烈的区域。

2000—2018年平均风蚀模数的变化趋势为减弱-增强-减弱（图3）。2001年平均风蚀模数最高，为39.89 t·hm^-2·a^-1，之后连续4年降低，2005年降到最低，为4.91 t·hm^-2·a^-1，比2001年降低34.98 t·hm^-2·a^-1。2005年之后平均风蚀模数先增大后减小，在2012年和2015年出现了两个小高峰，分别为30.79 t·hm^-2·a^-1和26.68 t·hm^-2·a^-1，但均未超过2001年的风蚀水平。2015年之后，平均风蚀模数持续降低，到2018年降低到7.99 t·hm^-2·a^-1。

从2000—2018年的平均月风蚀模数可以看出，12个月内均有土壤风蚀发生，但是土壤风蚀程度差异很大，1月土壤风蚀模数最高，之后呈降低趋势，6—9月风蚀模数最低，尤其是7月，只有0.01 t·hm^-2·month^-1，10月再次升高。这可能与植被盖度的变化有关。因为植被可通过固定表层土壤、增加地表粗糙度来减弱风蚀。坝上地区植被的非生长季大致为11月至次年4月，期间草原地区植被大部分处于枯萎状态，耕地处于裸露状态，加之冬春季节强风盛行，沙尘天气最为集中，从而导致11月至次年4月风蚀最严重，土壤风蚀模数最高。而5月草原地区植被逐渐返青，耕地开始耕作，因此自5月始，植被盖度逐渐增加。另外，由于坝上地区为季风气候，降水集中在7—9月，这期间植被生长状况最好，植被盖度最高，可有效降低风蚀。因此，5—10月土壤风蚀模数最小。

根据中华人民共和国水利部制定的风蚀强度分级标准（表1）^［34］，统计了坝上地区2000—2018年不同侵蚀强度所占的面积比例（表2），并对坝上地区的土壤风蚀风险进行评价。

坝上地区土壤侵蚀强度在中度侵蚀以上（含中度侵蚀）所占的比例总体呈减少的趋势，轻度侵蚀以下（含轻度侵蚀）所占比例呈增加趋势。本文将侵蚀级别高于中度侵蚀（风蚀模数高于50 t·hm^-2·a^-1）的区域设置为风蚀风险区，自2000年至2018年，只有2001、2002、2012年和2014年风蚀风险区比例高于5%，其余年份均低于1%。另外，根据风蚀风险区面积所占比例的大小，可将风蚀过程大致分为3个阶段：2000—2002、2003—2008年和2009—2018年，其风蚀程度第一阶段>第三阶段>第二阶段。

下垫面条件是影响土壤风蚀的重要因素。河北坝上地区地表性质复杂，按土地利用类型划分，既有易风蚀的土地利用类型（林地、草地、耕地、沙地等），又广泛分布着无法风蚀的土地利用类型（水域、建设用地等）。为了便于计算不同土地利用类型对土壤风蚀的影响，本文利用2000、2005、2010、2015年4期土地利用类型数据，将土地利用类型重新划分为耕地、林地、草地、沙地和非风蚀地，并分别计算了2000—2004、2005—2009、2010—2014、2015—2018年4个时间段内（以下分别用“第1、2、3、4年段”来表示）不同土地利用类型所占比例及其平均风蚀模数。

图4显示，不同土地利用类型所占比例从大到小依次为耕地（55.57%±0.68%）、草地（25.25%±0.36%）、林地（9.82%±1.01%）、沙地（5.46%±0.10%）、非风蚀地（3.90%±0.17%）。从第1年段到第4年段，草地和沙地面积所占比例分别从25.78%和5.56%减小到25.17%和5.37%，非风蚀地面积所占比例从3.78%增加到4.14%，耕地和林地面积所占比例虽略微增加，但变化很小。4个年段内，5种土地利用类型中，除非风蚀地不产生风蚀之外，其余4种土地利用类型的平均风蚀模数总体呈下降趋势。从第1年段到第4年段，耕地、林地、草地和沙地的平均风蚀模数分别从23.21、14.85、22.02、24.36 t·hm^-2·a^-1下降到16.81、10.49、15.94、17.62 t·hm^-2·a^-1。4种可风蚀的土地利用类型风蚀程度均为轻度侵蚀，其风蚀模数从大到小依次为沙地（18.89±4.39 t·hm^-2·a^-1）、耕地（17.14±5.07 t·hm^-2·a^-1）、草地（16.54±4.72 t·hm^-2·a^-1）、林地（10.81±4.18 t·hm^-2·a^-1）。

从4个年段的平均风蚀模数可以看出，第2年段的平均风速模数最低，耕地、林地、草地和沙地的平均风蚀模数分别为10.82、5.14、10.59、13.83 t·hm^-2·a^-1，显著低于其余3个年段。这是因为2003年国家开始实施京津风沙源治理工程，在第2年段，工程效果开始凸显，土壤风蚀模数开始降低。这一结果与相关研究结果一致^［35-36］。从不同土地利用类型的平均风蚀模数可知，耕地的土壤风蚀程度仅次于沙地，与沙地相差只有1.75 t·hm^-2·a^-1。这是因为坝上地区的耕地均为旱作农田，主要农作物为莜麦、亚麻、马铃薯等，由于当地长期形成的耕作习惯是秋季或早春进行农田翻耕，从而形成大面积的裸露农田，与冬春季节强风叠加，使该地区成为中国农田土壤风蚀沙化最严重的地区。林地的平均风蚀模数最小，这是由于林地枝叶茂盛、防风效果明显，且地表覆盖一定的凋落物，可固定地表土壤并增加地表粗糙度，进而削弱土壤风蚀。草地的平均风蚀模数略小于耕地，显著大于林地，可能由于部分草地被破坏，使表层土壤疏松，加速水分蒸发，使地表土壤含水量和植被盖度均小于林地。

植被盖度是定量评价土壤风蚀的重要指标。植被通过直接覆盖地表、截获风沙颗粒物、吸收气流动量以及改变气流场分布等途径抑制土壤风蚀^［37-39］。有关植被对土壤风蚀的影响，从定性分析、实验模拟到模型评估，国内外学者做了大量研究^{［36，40-44］}。本文使用NDVI数据逐月计算坝上地区的植被盖度，得到年平均植被盖度和月平均植被盖度。从图5可以看出，植被覆盖的高值区出现在坝上地区东北部和东南部，即丰宁县东北部和围场县，这些地区风蚀程度以微度侵蚀为主，甚至有些地区风蚀模数为0。而其余地区植被盖度均低于50%，土壤风蚀相对严重。2000—2018年植被盖度总体呈上升趋势，从2000年的27.48%增加到2018年的36.20%，与土壤风蚀模数的变化总体上为负相关。比如2003—2008年，平均植被盖度为31.41%，较2000—2002年平均增加3.71%，这6年的土壤风蚀模数均在9 t·hm^-2·a^-1以下，平均为7.13 t·hm^-2·a^-1，相比2000—2002年平均降低19.07 t·hm^-2·a^-1；而从2008年到2009年，植被盖度从32.81%降到28.17%，土壤风蚀模数则从4.89 t·hm^-2·a^-1增加到16.64 t·hm^-2·a^-1。

植被盖度和风蚀模数存在很好的负相关关系，不仅体现在两者年平均值的变化上，也体现在月平均值变化上。12月至次年4月植被盖度最低，均在20%以下，其中，1—4月植被盖度从9.39%缓慢升高到18.44%。5月为坝上地区植被返青期，植被盖度迅速上升到30.32%。6—8月植被盖度持续增加，8月增加到最大值64.63%。9月之后气温降低，植被枯萎，植被盖度开始下降，9—12月，植被盖度从50.41%下降到15.48%。

植被盖度与土壤风蚀模数之间的关系说明植被在一定程度上能够抑制土壤风蚀，这主要是由于在干旱、半干旱地区，土壤易发生风蚀，植被可固定表层土壤、保持土壤水分、增加地表粗糙度，进而抑制土壤风蚀。

土壤湿度是影响土壤风蚀强弱的重要因素，土壤表层水分的增加能有效增加土壤颗粒间的黏聚性，降低风力对土壤的破坏和剥离，减弱直至抑制风蚀过程^［45］。根据土壤表层0~10 cm的土壤湿度数据，得到2000—2018年平均土壤湿度的空间分布和逐年变化（图6）。土壤湿度分布具有很好的空间变化规律，即从东南向西北土壤湿度逐渐减小。其中，尚义县和丰宁县南部以及围场县东北部土壤湿度最高，这些地区也是土壤风蚀的低风险区。与之相反，康保县北部、围场县西北部以及沽源县和丰宁县北部的小部分区域土壤湿度最低，这些地区也是风蚀最为严重的地区。2000—2018年超过半数的年份康保县北部土壤风蚀最为严重；其次为沽源县和丰宁县北部地区，在2003—2004年成为研究区内风蚀最严重的区域；围场县西北部只在少数几个年份土壤风蚀程度较大，但均低于前面两个地区。结合这3个区域的植被盖度，其大小为围场县西北部>沽源县和丰宁县北部>康保县北部，可以看出，土壤湿度基本相同的情况下，植被盖度越高，土壤风蚀模数越小。这说明土壤湿度和植被盖度均是影响土壤风蚀的主要因素。

土壤湿度的年际变化不显著，为14.85%~18.21%，总体上略呈上升趋势。但从月平均土壤湿度可以看出，土壤湿度随月份存在很好的变化规律。3—7月，月平均土壤湿度逐渐增大，从11.93%增加到19.17%；自8月至次年3月，月平均土壤湿度逐渐减小，从18.75%减小到12.29%。月平均土壤湿度和风蚀模数之间存在很好的负相关关系，R²为0.85，表明土壤湿度是影响土壤风蚀的主要因素。同时，土壤湿度的月变化和植被盖度的月变化具有一定的相似性，两者均是在冬春季节较低，夏秋季节较高。这是因为坝上地区为温带大陆性气候，雨热同期，夏季降水多，气温高，有利于植被生长，使植被盖度增加，而植被盖度的增大又能够提高土壤保水能力；秋季之后，降水稀少，加之气候极其寒冷，导致大部分植被枯萎，植被盖度逐渐降低。

在干旱半干旱地区，风是土壤风蚀形成与发展的自然动力，是造成风沙危害的直接驱动力^［46］。研究表明，风蚀速率与风速呈正相关关系，随着风速增大，风蚀速率相应增加^［47-48］。起沙风（>5 m·s^-1）的大小和累计时间对输沙量有很大影响，起沙风累计时间越长，越容易发生风蚀。图7显示，2000—2018年坝上地区大于5 m·s^-1的起沙风速累计时间总体呈下降趋势，与风蚀模数的变化正相关，R²为0.69。2000—2018年，起沙风累计时间由892.76 h下降到716.93 h，风蚀模数则从16.44 t·hm^-2·a^-1减小到7.99 t·hm^-2·a^-1。起沙风累计时间和风蚀模数的最大值和最小值均出现在2001年和2005年，最大值分别为1 171.94 h和39.89 t·hm^-2·a^-1，最小值分别为456.83 h和4.91 t·hm^-2·a^-1。2003—2007年，起沙风累计时间较短，为685.39 h，相应的风蚀模数也较低，为7.58 t·hm^-2·a^-1；而2009—2013年，起沙风持续时间增加至938.81 h，风蚀模数也增加到18.96 t·hm^-2·a^-1。风蚀模数随起沙风累计时间的这种变化说明起沙风累计时间是影响土壤风蚀的主要因素。

由此可见，土地利用方式、植被盖度、土壤湿度和起沙风累计时间均是影响土壤风蚀的重要因素，这些因素之间又相互影响。不同土地利用方式和植被盖度变化均会导致地表土壤理化性质和空气动力学粗糙度的变化，继而对土壤湿度产生影响；土壤湿度的变化又会反过来影响植被盖度和土地利用方式；而土壤湿度和植被盖度的变化又直接影响地表沙物质的临界起动风速，进而影响起沙风累计时间。比如退耕还林，土地利用方式由耕地转为林地，植被盖度增加，较高植被盖度不仅可以有效降低风速，使起沙风累计时间降低，还可降低土壤表层的水分蒸发，提高林地所覆盖区域的土壤湿度，增大地表沙物质的临界起动风速，最终降低土壤风蚀。因此，未来可通过改善土地利用方式或增加植被盖度等措施在一定程度上改善土壤风蚀。但是影响土壤风蚀的因素众多，如人类活动、气候变化以及各影响因素之间的相互作用等。虽然坝上地区实施了一系列生态保护工程，如三北防护林工程、京津风沙源治理工程等，但仍有部分区域受降水、风速等气候因素以及过度樵采、放牧等人类活动的影响，植被恢复效果不好，土壤风蚀依然严重。因此，在干旱半干旱地区，应遵循植被的地带性规律，对于天然草原应以自然修复为主，而耕地应采取保护性耕作，尽量免耕或少耕，减少土壤扰动，并在条件适宜的区域适当种植人工林。

本文利用IWEMS模型对河北坝上地区的土壤风蚀模数进行了评估，但是缺少野外实地验证，下一步计划对坝上地区不同地点进行¹³⁷Cs土壤样品采集，并根据所测数据进一步修正IWEMS模型，以期提高风蚀预报的准确性。

本文针对河北坝上地区2000—2018年土壤风蚀模数的时空变化及其影响因素，结合气象、土壤、遥感和土地利用等数据，利用IWEMS定量评估了河北坝上地区土壤风蚀模数的时空分布状况，分析了土地利用类型、植被盖度、土壤湿度和起沙风累计时间对土壤风蚀模数的影响，并探究了坝上地区的风蚀治理措施。

坝上地区平均年风蚀模数变化为减弱-增强-减弱，但总体呈减弱趋势，平均月风蚀模数变化为每年的11月至次年4月土壤风蚀模数最高，其余月份风蚀模数较低。土壤风蚀程度以轻度和中度侵蚀为主，风蚀灾害主要发生在坝上地区中西部的康保县、张北县北部和丰宁县西部。自2009年以来，康保县和张北县的风蚀模数逐渐降低，丰宁县西部的风蚀模数逐渐增加，成为坝上地区风蚀最严重的区域。

风蚀风险区所占面积比例总体呈减小趋势，只有少数年份的风蚀风险区比例高于5%，其余年份均低于1%。尤其是京津风沙源治理工程实施之后的6年，土壤风蚀一直处于较低状态，风蚀风险区面积所占比例均低于3%。

土地利用类型、植被盖度、土壤湿度和起沙风累计时间均是影响土壤风蚀的重要因素，不同土地利用类型的土壤风蚀模数存在很大差异（沙地>耕地>草地>林地），风蚀模数与植被盖度和土壤湿度呈负相关关系，与起沙风累计时间呈正相关关系。因此，可通过调整土地利用类型、增加植被覆盖度、保护性耕作等措施降低土壤风蚀。