近20 a蒙古国土壤风蚀变化特征及主要影响因素分析

郭 茵， 雷加强,3， 范敬龙,3， 王海峰,4， 吕振涛

（1.中国科学院新疆生态与地理研究所国家荒漠-绿洲生态建设工程技术研究中心，新疆 乌鲁木齐 830011；2.中国科学院大学，北京 100049；3.中国科学院新疆生态与地理研究所塔克拉玛干沙漠研究站，新疆 库尔勒 841000；4.中国科学院新疆生态与地理研究所莫索湾沙漠研究站，新疆 石河子 832000）

土壤风蚀是土地退化最主要的原因之一。世界上约三分之一的土地受到风蚀的影响，其中严重风蚀土地占比达50%以上［1-3］。土壤风蚀会对陆地生态系统以及人类社会发展产生诸多不利影响［4］，如显著降低土壤肥力和植物生产力［5］，导致空气质量降低等［6］。此外，地表风蚀形成的细粒物质是大气气溶胶的主要来源［7］，为沙尘暴提供了物源。只要有发生风蚀的条件，就可能产生较大规模的风沙灾害［8］，如2021年3月14日蒙古国发生的特大沙尘暴事件就危及到东亚广大地区。

土壤风蚀的定量评估方法包括野外风蚀直接测量［9］、风洞模拟实验［10］、基于同位素示踪法的风蚀速率计算等［11］。利用风蚀监测设备或风洞模拟实验具有一定的局限性［12］，难以对较大空间范围进行风蚀评估。随着遥感技术的发展，越来越多的学者利用地理信息系统进行风蚀建模，发展出WEQ、RWEQ、WEPS、SWEEP 等众多模型，其中以RWEQ模型应用最为广泛。RWEQ 模型以WEQ 模型为基础，融合了自然过程和实际经验，能够将野外数据与计算机模型相结合，来描述物理的风蚀过程［13］，在实际中得到了广泛的应用。RWEQ 模型起初是美国农业部（USDA）为估算耕地土壤流失而开发的一种基于经验的模型［14］，之后经过数次改进和验证，在其他各地得到了广泛的应用。Zhang 等［15］使用RWEQ模型，对中国内蒙古地区1990—2015年的土壤风蚀进行了模拟，评估了气候变化对风蚀的影响；Lin 等［12］利用RWEQ 模型估算了1982—2015 年中国河西地区的潜在风蚀率，Chi 等［16］利用RWEQ模型评估了中国土壤风蚀的时空格局，Li等［17］利用RWEQ 模型模拟了1986—2099 年中亚土壤风蚀的时空动态。蒙古国是东亚地区沙尘暴事件发生频率最高的地区［18］，是东亚地区沙尘暴灾害的主要起源地［19］。前人使用基于专家的方法评估了蒙古国的风蚀速率和程度，结果表明蒙古国大约30%的领土受到风蚀的影响，其中戈壁沙漠被认为是受风蚀影响最严重的区域［20］。Natsagdorj 等［21］利用气象观测和激光雷达，发现蒙古国各地的风蚀总体上有所增加。师华定等［22］利用地理信息系统和遥感技术，通过模糊c 均值聚类（FCM）方法对植被覆盖率、地形起伏度、土壤干燥度和风场强度4 个主要蒙古高原风蚀环境因子数据进行模糊聚类，结合专家知识，将环境组合与特定的风蚀危险度相对应，得到风蚀危险度图，以此对蒙古高原土壤风蚀危险度进行评价研究，并分析了蒙古高原土壤风蚀危险度的空间格局特征。FCM方法可有效地获取风蚀-环境关系知识，为预测性风蚀危险度制图提供依据，此研究并没有对蒙古高原土壤风蚀危险度进行定量评价。Mandakh 等［23］在ArcGIS 环境中使用风蚀方程模型（WEQ）评估了蒙古国各地的风蚀情况，尽管此研究对蒙古国的风蚀进行了定量分析，但是缺乏对蒙古国长时间序列的风蚀趋势的分析。目前，关于蒙古国的风蚀研究较为有限，亟待开展蒙古国长期风蚀状况定量描述和评估研究。本文使用RWEQ 模型，计算蒙古国全域的风蚀模数，揭示风蚀的时空变化规律和趋势，分析识别主要影响因素，有助于深入了解蒙古国土壤风蚀整体特征，为重点地区风沙灾害和荒漠化防治工作提供重要参考。

1 研究区概况

蒙古国位于亚洲东中部（41°32′～52°15′N，87°44′～119°56′E），面积156.65×104km2［24］，是世界第二大内陆国家。蒙古国平均海拔约1580 m，最高点（西部的khuiten 峰）海拔4653 m，最低点（东部的Khookh Nuur 湖洼地）海拔518 m，首都乌兰巴托市海拔约1310 m。蒙古国远离海洋，具有显著的大陆性气候特征，四季分明，年均温和日均温变化较大；年平均降水量约120～250 mm，降水量从北到南、从东到西递减，降水多集中在7—8月［25］。根据蒙古国家气象和环境监测局（https://eic.mn/）提供的蒙古国自然地带分布数据显示，南部为沙漠和荒漠草原地区，东南部是戈壁沙漠区；典型草原主要分布在蒙古国中部东西延伸；森林集中分布在蒙古国北部山区（图1）。其中，荒漠草原、半荒漠和沙漠地区约占41.3%，这些地区大风频发，土壤风蚀发生强烈，生态环境脆弱，年平均风速为4～6 m·s-1［23］。蒙古国南部的荒漠地区由于土壤疏松，更易被风蚀并形成沙尘暴［23,26］，年沙尘暴日数可达30～120 d；蒙古国是一个沙漠化严重的国家，沙漠化的加剧直接增加了风蚀和沙尘排放源，威胁到东亚其他国家的空气质量［27］。

图1 研究区概况示意图Fig.1 Geographical distribution map of Mongolia

2 研究方法与数据来源

2.1 研究方法

2.1.1 风蚀模型 RWEQ模型充分考虑了包括气候条件、土壤特性、表面粗糙度、植被覆盖和其他环境要素，主要基于式（1）～（3）来估算风蚀量。

式中：SWEQ表示风蚀模数（t·hm-2·a-1）；Qmax表示最大输沙能力(kg·m-1)；s表示关键地块长度(m)；x是模拟区域不可侵蚀边界到下风口的距离(m)；WF表示气候因子(kg·m-1)；EF表示土壤可蚀性因子（无量纲）；SCF表示土壤结皮因子（无量纲）；K′表示土壤粗糙度影响因子（无量纲）；COG表示植被影响因子。具体计算过程如下式所示：

式中：Wf表示风影响因子(kg·m-1·s-1)；ρ表示空气密度(kg·m-3)；g表示重力加速度(m·s-2)；SW表示土壤湿度因子；SD表示积雪覆盖影响因子；W表示风力值（m3·s-3）；U2指2 m 高的风速（m·s-1）；Ut指风速阈值，一般取5 m·s-1；Nd是风速测量的时间间隔（d）；N代表风速观测频率（d）；EL表示从数字高程模型（DEM）数据获得的高程（km）；T表示绝对温度（K）；SD由积雪覆盖度计算获得；P表示测定时段内积雪覆盖深度大于25.4 mm的发生概率。

土壤可蚀性因子EF计算过程如下式所示：

式中：Sa为土壤砂粒含量(%)；Si为土壤粉粒含量(%)；Cl是黏土含量（%）；Sa/Cl为土壤砂粒与黏土之比；OM为土壤有机质含量(%)；CaCO3为碳酸钙含量（%）。

土壤结皮影响因子SCF计算过程如下式所示：

土壤质地和有机质含量的年际变化非常缓慢，可将土壤可蚀性因子和结皮因子视为常量值，假设土壤可蚀性因子和结皮因子不随时间改变。

地表粗糙度影响因子K′，根据Xu 等［28］的研究结果，计算如下式所示：

式中：α表示地形坡度，可从DEM数据中提取。

植被影响因子COG表示一定植被覆盖条件下对土壤风蚀的抑制程度，采用最大合成法（MVC）得到半月NDVI 数据，并用像元二分法求得2001—2020年的半月植被覆盖度，计算过程如下式所示：

式中：FVC 为植被覆盖率(%)；NDVIsoil为裸土像素的NDVI值；NDVIveg为植被的NDVI值。

2.1.2 变化趋势分析方法 利用Sen’s趋势分析法，计算了2001—2020 年蒙古国土壤风蚀以及气象因子（风速、气温、降水）的变化趋势。该方法是非参数检验，只需要保证样本的独立性，不需要拟合正态分布［29］，方法如下式所示：式中：xi和xj是第i年和第j年的风蚀模数值或风速、气温、降水的数值；n是研究期的总年数。Slope表示时间序列数据的变化率，正值表示风蚀或气象因子（风速、气温、降水）随时间变化总体呈增加趋势，负值表示总体呈减少趋势。

2.2 数据来源

（1）气象数据

气象数据包括风速、气温、降水、积雪深度、土壤水分，均来自全球陆地数据同化系统GLDAS2.1产品［30］。GLDAS数据是美国航空航天局（NASA）戈达德空间飞行中心（GSFC）和美国海洋和大气局（NOAA）国家环境预报中心（NCEP）联合发布的基于卫星、陆面模式和地面观测数据的同化产品，它能够提供多种驱动数据，这些数据来源于大气同化产品、再分析和实际观测。GLDAS2.1产品时间分辨率为3 h，空间分辨率为0.25°，与地面气象观测站实测数据具有较高的一致性［30］。

（2）土壤数据

土壤数据包括土壤有机质含量、黏粒含量、粉粒含量、沙粒含量、CaCO3含量，其中土壤有机质含量、黏粒含量、粉粒含量、沙粒含量来源于Open Land Map(OLM)数据集，CaCO3含量来源于世界土壤数据库（HWSD）。

（3）植被数据

选用由美国航空航天局（NASA）发布的NDVI中分辨率成像光谱仪获得的陆地2级标准数据产品MOD13Q1，时间分辨率为16 d，空间分辨率为250 m，可用于计算植被覆盖度。

（4）地形数据

选择由美国航空航天局（NASA）和美国测绘局（NIMA）联合测量的数字高程模型（SRTM）版本4，空间分辨率为90 m。这是目前质量最高的SRTM数据集，水平标准误差为1 m，垂直标准误差为0 m，可用于计算坡度和地表粗糙度。

（5）其他数据

牲畜数量数据来自于蒙古国国家统计局（www.1212.mn）。土地利用数据使用的是由中国国家基础地理信息中心（NGCC）提供的GlobeLand30产品［31］，数据集使用了包括20000 多张美国陆地资源卫星（Landsat）和中国环境减灾卫星（HJ-1）的多光谱影像，具有较高的空间分辨率和数据精度［32］。

将以上所有栅格数据的空间分辨率重新采样至500 m，时间分辨率为16 d，全年有23 个时间段。对计算结果进行汇总，得到每年的风蚀总量。

3 结果与分析

3.1 风蚀模数的空间变化特征

从蒙古国土壤风蚀模数空间分布来看（图2），土壤风蚀变化空间差异较大，风蚀模数的最大值与最小值的差值在0～60 t·hm-2·a-1之间。整体上，土壤风蚀模数呈现北部低、南部高的空间分布模式，由于南部的戈壁沙漠区土壤质地较粗，物源丰富，且可蚀性较高，加之降水少，地表植被稀疏，因此该区成为蒙古国土壤风蚀最为严重的地区。

图2 2001—2020年蒙古国土壤风蚀模数空间分布特征Fig.2 Spatial distribution characteristics map of wind erosion in Mongolia from 2001 to 2020

由图3、图4可知，对风蚀模数进行多年平均，从年内空间分布可以看出，风蚀整体上春季（3—5月）最强，约占全年风蚀模数的45%；其次为秋季（9—11月），约占全年风蚀模数的33%；再次为夏季（6—8月），约占全年风蚀模数的12%；冬季（12月—次年2月）最弱，约占全年风蚀模数的10%。蒙古高原春季降水稀少，土壤干燥，植物尚未进入快速生长期，植被覆盖度低，因此春季风蚀非常强烈。冬季土壤表层冻结，且多有积雪覆盖，因此风蚀最弱。夏季和秋季植被条件均较好，但夏季降水集中，土壤湿度明显大于秋季，因此夏季风蚀强度小于秋季。气象观测表明，蒙古国的沙尘事件主要发生在春季［18］，与本文风蚀主要发生在春季的结论一致。

图3 2001—2020年蒙古国季节平均土壤风蚀模数空间分布Fig.3 Spatial distribution map of seasonal average soil wind erosion in Mongolia from 2001 to 2020

图4 2001—2020年蒙古国季节平均土壤风蚀模数变化Fig.4 Variation of seasonal average soil wind erosion modulus in Mongolia from 2001 to 2020

3.2 风蚀模数的时间变化特征

由图5 可知，在2001—2020 年间，蒙古国土壤风蚀模数总体呈显著上升趋势（R2=0.2258，P＜0.05），其中2018年最强（4.4463 t·hm-2·a-1），2008年最弱（1.369 t·hm-2·a-1）。从时间序列的变化特征来看，2005—2009 年的风蚀区域明显小于其他年份，而2009 年后风蚀范围明显增大。根据风蚀模数的年际变化，2001—2020 年期间风蚀变化分为2 个明显的阶段，2001—2008 年为风蚀模数下降时期（R2=0.536，P＜0.05），2010—2020 年为风蚀模数增加时期。

从土壤风蚀变化速率空间差异来看（图6），蒙古国南部戈壁沙漠地区的风蚀变化速率约为2 t·hm-2·a-1，明显高于其他区域；蒙古国北部风蚀最低，变化最小；蒙古国西南部风蚀速率变化为负，表明此区域风蚀呈下降趋势，平均每年的风蚀变化速率约-0.01～-0.6 t·hm-2·a-1。总体来看，蒙古国土壤风蚀模数高值区变化速率快，而低值区变化速率慢。蒙古国北部以森林生态系统为主，植被类型以高山泰加林、高山森林和森林草原为主，植被覆盖度较高，土壤风蚀阈值不大，整体森林生态系统稳定，受气候变化及人为活动影响小，风蚀很小，且变化率很小；蒙古国南部荒漠生态系统脆弱，植被稀少，受气候变化影响大，因此风蚀年际变化率也较大；蒙古国西南部风蚀呈下降趋势的主要原因是此区域风速呈下降趋势，风蚀量常年较低。

图6 2001—2020年蒙古国土壤风蚀模数变化趋势的空间分布Fig.6 Spatial distribution of variation trend of soil wind erosion in Mongolia from 2001 to 2020

3.3 土壤风蚀的影响因子分析

3.3.1 气候因子变化及影响 影响土壤风蚀的主要气候因子有风速、气温、降水等。本文采用Sen’s趋势分析法分析了蒙古国主要气候因子变化趋势及空间分布。从图7 可以看出，近20 a 来蒙古国年最大风速呈明显波动变化，2008年最低（14.04 m·s-1），2001 年最高（16.26 m·s-1）；2001—2008 年，年最大风速呈下降趋势（R2=0.427，P＞0.05）；2009—2020年，年最大风速呈显著增加趋势（R2=0.3959，P＜0.05）。土壤风蚀模数随着风速变化趋势的波动而波动，二者变化规律基本呈现对应关系（图5，图7），表明风速是影响土壤风蚀最重要的自然动力条件。

图5 2001—2020年蒙古国多年平均土壤风蚀模数变化趋势Fig.5 Variation trend graph of average soil wind erosion in Mongolia from 2001 to 2020

图7 2001—2020年蒙古国年最大风速变化趋势Fig.7 Variation trend of annual average maximum wind speed in Mongolia from 2001 to 2020

结合年均风速变化趋势空间分布（图8）来看，蒙古国南部戈壁沙漠区风速显著上升，其风蚀量呈增加趋势；中部风速下降，风蚀量呈减少趋势。此外，将蒙古国2001—2020年逐年平均风速与逐年风蚀模数进行相关性分析（图9）发现，在蒙古国仅3.3%左右的区域平均风速与风蚀模数呈负相关关系，绝大数区域（96.6%）风蚀模数都与平均风速呈正相关。尤其是在蒙古国南部及西南部的戈壁区域，风速与风蚀量呈显著正相关关系，说明风速是影响土壤风蚀最重要的自然动力条件。这也解释了蒙古国西南部风蚀量变化速率下降是受风速下降的影响。

图8 2001—2020年蒙古国年均风速变化趋势空间分布Fig.8 Spatial distribution of annual average wind speed change trend in Mongolia from 2001 to 2020

图9 蒙古国年均风速与风蚀模数相关关系Fig.9 Correlation between annual average wind speed and wind erosion modulus in Mongolia

气温和降水会通过影响土壤湿度和植物生长而间接影响风蚀［33］。从图10 可以看出，近20 a 来，蒙古国大部分地区年均温呈上升趋势，仅有北部少部分区域年均温呈下降趋势。南部增温最高，增温速率达到0.06 ℃·a-1，升温促进了地表水的蒸发，导致土壤湿度降低，可蚀性增加，特别是在缺乏植被的地区，升温增加了风蚀的发生。降水与表层土壤含水量密切相关，降水量减少增加了风蚀的发生。整体来看，近20 a来蒙古国降水呈显著增加趋势（图11b），但降水变化存在较大空间差异（图11a），北部地区降水明显增加，中部和南部部分地区降水减少，尤其南部戈壁地区呈明显减少趋势（图11c）。南部地区干旱灾害发生的可能性增加，加上风力的作用，风蚀量也呈增大趋势。

图10 2001—2020年蒙古国年均温变化趋势时空分布Fig.10 Temporal and spatial distribution of annual mean temperature change trend in Mongolia from 2001 to 2020

图11 2001—2020年蒙古国年降水量变化趋势时空分布Fig.11 Temporal and spatial distribution of annual rainfall variation trend in Mongolia from 2001 to 2020

3.3.2 植被覆盖度变化及影响 植被覆盖可增加地

表粗糙度和削弱近地表风速，从而降低地表风蚀危险性［15,34］。图12 为近20 a 来蒙古国FVC 的变化，可以看出，蒙古国大部分区域FVC 整体呈增长趋势，东部部分地区植被覆盖度快速增加，因此土壤风蚀降低；但部分地区尤其是南部和西部的荒漠草原区多地植被覆盖度减少，导致风蚀发生风险增大。

图12 2001—2020年蒙古国植被覆盖度变化趋势Fig.12 Change trend of annual vegetative cover in Mongolia from 2001 to 2020

4 讨论

4.1 风蚀模数的结果验证

RWEQ 模型起初是美国农业部（USDA）为估算耕地土壤流失而开发的一种基于经验的模型［14］，之后经过数次改进和验证，不仅是将其用于评估农田发生的土壤风蚀，如Li等［35］将其用于评估内蒙古地区的固沙服务功能，Chi 等［36］评估了中国全域的土壤风蚀，也有一些学者将其用于中亚地区的土壤风蚀评估［17,37］。然而其结果往往需要其他方法估算的结果来对比验证。

综上所述，蒙古国2001—2020年多年平均土壤风蚀模数为2.736 t·hm-2·a-1。137Cs示踪法是估算土壤侵蚀速率的常用方法［38］。前人基于137Cs 监测方法计算了蒙古高原的风蚀模数，蒙古国境内的监测样点从北向南涵盖了森林草原区、典型草原区、荒漠化草原区、草原化荒漠区、戈壁区、农牧交错区［39］。将本研究中基于RWEQ 模型模拟的20 a 的土壤风蚀模数与刘纪远等［39］使用137Cs 实测的土壤风蚀模数基于相同的地理坐标点进行了相关性分析（图13），结果表明，RWEQ 模型模拟结果与137Cs 测量值显著相关（R=0.91，P＜0.01）。Mandakh 等［23］使用WEQ模型评估的蒙古国风蚀模数在2.7～27.5 t·hm-2·a-1之间，比本文用RWEQ 估算结果略偏低，但数量值大体上相当。土壤风蚀的发生和发展依赖于侵蚀动力（气候因子如风速）和可蚀性因子（地表状况如植被盖度、土壤特性等)之间的相互作用，土壤可蚀性因子可以衡量土壤表层被风力侵蚀的相对可能性。Jugder等［40］从土壤可蚀性的角度评价了蒙古国土壤被风力侵蚀的相对可能性，发现从蒙古国南部的沙漠和荒漠草原地区到西部的沙地，存在着广泛的高可蚀性区域，与本文中的严重风蚀区域也大致相同。中等可蚀性区域主要覆盖西部的干旱草原区和荒漠草原区，低可蚀性区域从蒙古国西部延伸到北部。Mandakh 等［23］使用WEQ 模型评估了蒙古国风蚀模数，分析了风蚀的空间分布规律，与本文结论一致。

图13137Cs测量和RWEQ模型模拟风蚀模数的相关系数Fig.13 Correlation coefficient of137Cs measurement and RWEQ model simulation wind erosion modulus

4.2 人为活动因素变化及影响

土壤风蚀的发生与气候条件和下垫面条件有密切关系，人为活动通过影响下垫面而影响风蚀的发生。过度放牧以及煤矿开采等人类活动加剧了蒙古国土地退化［41］。由图14 可知，在2001—2020年间，蒙古国的牲畜数量由2001年的26075千头增加到2020年的67068千头，每年大约以2400千头的速度显著增加。牲畜数量的增加使草原压力增大，而宽松的放牧政策造成多数草场过度放牧，牲畜踩踏和过牧行为破坏了土壤表层结构，降低了冬、春季地表的植被盖度，增加了土壤的风蚀风险。因此，过度放牧可能是导致蒙古国土壤风蚀增加的一个人为活动因素。

图14 2001—2020年蒙古国牲畜数量变化趋势Fig.14 Variation trend of livestock quantity in Mongolia from 2001 to 2020

近20 a来，蒙古国经济社会快速发展，为了适应生存和发展的需求，土地利用类型也发生了较大变化。从土地利用类型转化（图15）来看，蒙古国境内有1733 km2的草地转化为耕地，部分未利用地和林地也有转化为耕地的情况。与传统畜牧业生产方式相比，发展种植业，会导致更严重的土壤风蚀发生。其他土地类型转化为耕地后，人为翻动表土，土壤结构、水分和养分状况发生了很大变化，导致局部地区风蚀加剧。

图15 2001—2020年蒙古国主要土地利用类型转化Fig.15 Transformation of main land use types in Mongolia from 2001 to 2020

4.3 影响因素耦合的空间异质性作用

气温和降水通过影响土壤湿度和植被覆盖度而间接地影响风蚀。2001—2020 年蒙古国全域气温和降水均呈增加趋势，植被覆盖度也呈增加趋势，对地表风蚀应发挥一定的抑制作用，但全域风蚀仍整体呈增加趋势，似乎存在一定的矛盾。但分析发现，原因在于影响风蚀因素耦合存在较大的时空异质性。对于一些风蚀相对严重的地区，植被覆盖度的增加对减少风蚀没有明显影响，这是因为高风蚀区域的植被覆盖度较低，且增加不显著（图12），风速等其他因素的变化对风蚀有更显著的影响，而且在正值风季的春季作用较为强烈。学者们对蒙古国土壤风蚀的研究中发现，不同的气候因素影响不同的区域，干旱对戈壁和沙漠地区影响更大［23］。因此，从蒙古国全域的尺度来看，风速是土壤风蚀发生的主导因素，风速越大，风蚀越强。而且强风蚀的干旱荒漠土地面积占比极高，因此蒙古国整体风蚀仍呈增加趋势。

5 结论

本文通过RWEQ 模型对蒙古国2001—2020 年的土壤风蚀时空变化状况及对主要影响因子进行了分析，得到的结论如下：

（1）蒙古国土壤风蚀年内季节变化明显，全年45%的风蚀发生在春季，秋季其次，冬季最少。

（2）从2001—2020 年，蒙古国土壤风蚀呈年际波动变化，总体上呈显著增加趋势，风蚀模数以0.06 t·hm-2·a-1的速度增长。

（3）土壤风蚀具有明显的空间异质性，整体呈南部高、北部低的空间分布格局，风蚀严重区域位于蒙古国南部与中国接壤的戈壁和沙漠地区；风蚀变化速率空间变异较大，土壤风蚀高值区大于低值区。

（4）蒙古国土壤风蚀变化整体上与气候因子及引起的植被盖度变化有关，也与过牧、耕地开垦等人为活动增强有密切关系。为了有效防治沙漠化危害，蒙古国应在土地资源开发利用方面加强科学管理，对于严重沙漠化地区应积极采取可行的治理措施。