2001—2020年中国北方干旱半干旱区土壤风蚀的时空变化及驱动因素

doi:10.7522/j.issn.1000-694X.2025.00149

2001—2020年中国北方干旱半干旱区土壤风蚀的时空变化及驱动因素

王世宇^,, 王雪松^,

北京师范大学地表过程与水土风沙灾害风险防控全国重点实验室/防沙治沙教育部工程研究中心/地理科学学部，北京 100875

Spatiotemporal variation and driving factors of soil wind erosion in the arid and semi-arid areas of northern China in 2001-2020

Wang Shiyu^,, Wang Xuesong^,

State Key Laboratory of Earth Surface Processes and Disaster Risk Reduction / MOE Engineering Research Center of Desertification and Blown-Sand Control / Faculty of Geographical Science，Beijing Normal University，Beijing 100875，China

通讯作者: 王雪松（E-mail： xswang@bnu.edu.cn）

收稿日期: 2025-07-24 修回日期: 2025-09-08

基金资助:

国家自然科学基金青年科学基金项目. 42107343
地表过程与水土风沙灾害风险防控全国重点实验室重点项目. 2024-ZD-02

Received: 2025-07-24 Revised: 2025-09-08

作者简介 About authors

王世宇（2002—），男，新疆兵团人，研究生，主要从事土壤风蚀方面的研究E-mail：sywang1210@mail.bnu.edu.cn , E-mail：sywang1210@mail.bnu.edu.cn

摘要

中国北方干旱半干旱区地形复杂多样、生态环境较为脆弱，土壤风蚀是影响该地区生态系统稳定和社会经济发展的关键因素。本研究采用修正风蚀方程（RWEQ），定量分析了该地区2001—2020年土壤风蚀的空间分布格局和时间动态变化，并结合地理探测器进一步探究了土壤风蚀的驱动因素。结果表明：（1）研究时段内，研究区的多年平均实际土壤风蚀模数为13 449.59 t·km^-2·a^-1，以剧烈侵蚀为首，主要分布在塔里木盆地中心地带和内蒙古高原中西部；实际土壤风蚀模数总体呈现波动下降的态势，大致以每年102.78 t·km^-2的速率下降。（2）研究时段内，研究区不同土地利用类型的实际土壤风蚀模数存在显著差异。无论是发生转换的土地利用类型，还是未发生转换的土地利用类型，2020年的实际土壤风蚀模数：沙地>草地>耕地>建设用地>林地。此外，中国北方干旱半干旱区新增草地和耕地的实际土壤风蚀模数相对较高，分别为13 695.88 t·km^-2·a^-1和7 933.25 t·km^-2·a^-1，比研究时段内未发生转换的草地和耕地的实际土壤风蚀模数分别高19.90%和58.53%。（3）在影响研究区土壤风蚀的众多因素中，气候因子是主导因素（q=0.312，P<0.01）；其中，风因子与实际土壤风蚀模数呈现出较为显著的中等强度的正相关（r=0.43，P<0.05）。

关键词： 中国北方干旱半干旱区 ; 修正风蚀方程 ; 地理探测器

Abstract

The arid and semi-arid areas of northern China have complex and diverse topography and fragile ecological environment. Soil wind erosion is a key factor affecting the stability of the ecosystem and social and economic development in the areas. This study uses the Revised Wind Erosion Equation （RWEQ） to quantitatively analyze the spatial distribution and temporal dynamics of soil wind erosion in the study area from 2001 to 2020， and combined with geographic detector to further explore the driving factors of soil wind erosion. The results indicate that：（1） During the study period， the multi-year average actual soil wind erosion modulus of the study areas is 13 449.59 t·km^-2·a^-1， which was dominated by severe erosion and mainly distributed in the center of the Tarim Basin and the central and western parts of the Inner Mongolia Plateau. The actual soil wind erosion modulus shows a fluctuating downward trend， which decreases at a rate of 102.78 t·km^-2 per year. （2） During the study period， there are significant differences in the actual soil wind erosion modulus of different land use types in the study area. Whether it was the land use type that has been converted or the land use type that has not been converted， the order of the actual soil wind erosion modulus in 2020 was： sandy land > grassland > cultivated land > construction land > forestland. In addition， the actual soil wind erosion modulus of the newly-added grassland and newly-added cultivated land in the arid and semi-arid areas of northern China is relatively high， which were 13 695.88 t·km^-2·a^-1 and 7 933.25 t·km^-2·a^-1 respectively， and were higher than the actual soil wind erosion modulus of the grassland and cultivated land that did not change during the study period， which were 19.90% and 58.53% respectively. （3） Among the many factors affecting the soil wind erosion in the study areas， the weather factor was the dominant factor （q=0.312， P<0.01）. Among them， the wind factor and the actual soil wind erosion modulus show a relatively significant positive correlation of moderate intensity （r=0.43， P<0.05）.

Keywords： the arid and semi-arid areas of northern China ; RWEQ ; geographic detector

PDF (5073KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

王世宇, 王雪松. 2001—2020年中国北方干旱半干旱区土壤风蚀的时空变化及驱动因素. 中国沙漠[J], 2026, 46(2): 164-176 doi:10.7522/j.issn.1000-694X.2025.00149

Wang Shiyu, Wang Xuesong. Spatiotemporal variation and driving factors of soil wind erosion in the arid and semi-arid areas of northern China in 2001-2020. Journal of Desert Research[J], 2026, 46(2): 164-176 doi:10.7522/j.issn.1000-694X.2025.00149

0 引言

土壤不仅是重要的自然资源，还是维系生态环境和保障人类生存的重要基础^［1］，在维持生态平衡，提供植物生长所需的养分、水分和空气等方面起着不可替代的作用^［2］。土壤的健康状况对水文循环、气候调节及其他环境要素也具有深远的影响。然而，土壤会受到多种外营力的影响，从而出现退化。土壤风蚀作为土壤退化的主要驱动因素，指在风力作用下，表层土壤中的细颗粒和营养物质经历吹蚀、搬运与沉积的过程^［3］。在干旱半干旱及部分半湿润地区，土壤风蚀过程尤为显著，是土地沙漠化的关键环节^［4］。土壤风蚀引起的土地退化、土壤肥力下降^［5］、植物生长发育不良^［6-7］、泥沙淤积和生态环境恶化等问题^［8］也日益严重。此外，在土壤风蚀过程中，大量粉尘悬浮于大气中，对大气质量产生负面影响，造成所在地区及下风向地区大范围的沙尘暴^［9］，危害人体健康^［10］，还对各类基础设施造成危害^［11］。目前，全球约28%的土地正在遭受土壤风蚀的危害^［12］。根据《2024年中国水土保持公报》的数据，中国受土壤风蚀影响的土地面积约为155.09万km²，约占国土面积的16.20%，主要分布在新疆、内蒙古、青海、甘肃和西藏等15个省（自治区）^［13］。土壤风蚀已成为全球性的环境问题，并且一直都受到中国政府的高度关注。

基于土壤风蚀的发生过程，其研究方法主要有野外调查观测、风洞模拟实验、同位素示踪技术等^［14-16］。这些方法虽已被广泛应用并取得了一定的成果，但其开展往往涉及庞大的资源消耗，包括人力、物力、财力资源，并受到空间尺度的限制，难以将研究成果进一步推广。为了在大尺度上估算土壤风蚀量以及预测未来风蚀状况，研究者们基于计算机模拟技术和地理信息系统开发了一系列土壤风蚀模型^［17］。目前，已有的模型主要包括美国的风蚀方程（Wind Erosion Equation， WEQ）^［18］、捷克的Pasak风蚀模型^［19］、苏联的Bocharov模型^［20］、美国的侵蚀-生产力影响计算公式（Erosion-Productivity Impact Calculator， EPIC）^［21］、美国的风蚀预报系统（Wind Erosion Prediction System，WEPS）^［22］、澳大利亚的风蚀评价模型（Wind Erosion Assessment Model， WEAM）^［23］、美国的修正风蚀方程（Revised Wind Erosion Equation， RWEQ）^［24］、美国的德克萨斯侵蚀分析模型（Texas Erosion Analysis Model， TEAM）^［25］、欧洲的轻质土壤风蚀模型（Wind Erosion on European Light Soils， WEELS）^［26］、美国的风蚀评估模型（Single-event Wind Erosion Evaluation Program， SWEEP）^［27］等。在众多风蚀评估模型中，RWEQ模型综合考量了气候条件、土壤质地、地表状况等关键因素，估算精度较高，尤其在大尺度范围和长期序列的土壤风蚀评估方面表现出色。因此，该模型在土壤风蚀模拟及防风固沙研究领域得到了广泛的应用。中国学者对RWEQ模型进行了应用及优化，巩国丽等^［28］对RWEQ模型的相关参数进行率定，分析了内蒙古草原的风蚀状况及影响因素^［29］以及研究了中国北方风蚀区植被防风固沙的功能^［30］。江凌等^［31］采用RWEQ模型，完成了青海省土壤风蚀模数的估算，以及研究了内蒙古2001—2010年的土壤风蚀状况^［32］。申陆等^［33-34］则将RWEQ模型应用到浑善达克沙地，研究浑善达克沙地土壤风蚀的时空变化及主要影响因子。Zhang等^［35］将RWEQ运用在中国北方农牧交错带，估算其潜在风蚀量等。综上，将RWEQ模型运用于中国北方地区进行风蚀量的模拟具有较为显著的实用性。同时，应用该模型能够相对精确地评估中国北方地区植被在防风固沙方面的效用。因此，可以优先考虑将该模型应用于相关领域的研究和实践。目前，针对中国北方干旱半干旱地区土壤风蚀的系统研究仍显不足，特别是土壤风蚀的驱动因素亟需深入探究。

中国北方干旱半干旱区拥有丰富的生物多样性及能源矿藏等资源，是中国关键性的资源保障基地和生态安全屏障，同时也是“一带一路”倡议的核心区域，其生态地位及保护价值极为重要。然而，由于该区域独特的地理位置和气候因素，冬春季受到蒙古-西伯利亚高压的影响，该地区盛行从大陆中心吹来的干燥季风，从而使得该区域土壤风蚀和沙尘暴频繁发生，严重影响了该地区的环境质量，成为影响区域生态安全和社会发展的重要制约因素^［36］。因此，定量研究中国北方干旱半干旱区土壤风蚀的时空变化特征，提高在气候变化背景下土壤对气候反馈作用的认识，从而减少土壤风蚀以及沙尘暴的发生，对于保护中国北方干旱半干旱区的生态安全、维护生物多样性、促进可持续发展以及保障国家资源安全等具有重要意义。

1 研究区概况

中国北方干旱半干旱区位于中国的西北部（30°07′—50°44′N、73°26′—127°03′E）^［37］。本研究中，干旱半干旱区的范围依据多年平均干旱指数（年降水量与年潜在蒸散发量的比值，Aridity Index，AI）进行界定^［1］，干旱区根据干旱指数分为干旱半湿润区（0.5<AI≤0.65）、半干旱区（0.2<AI≤0.5）、干旱区（0.05<AI≤0.2）和极端干旱区（AI≤0.05）^［38］（图1）。该区域的主要气候类型为温带大陆性气候，年降水量较少，大部分地区为50~400 mm，年降水量呈现出由东南向西北递减的特征，且降水季节分配不均，集中在夏季；年平均气温相对较低，且气温年较差和气温日较差均较大。该区域地形复杂多样，地势西高东低，横跨三级地势阶梯。该区域的土地利用类型以沙地、裸地和草地为主，植被类型多为旱生草本和灌木。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 中国北方干旱半干旱区的地理位置

注：基于自然资源部标准地图服务网站标准地图（审图号：GS（2024）0650号）制作，底图边界无修改

Fig.1 Geographical location in the arid and semi-arid areas of northern China （ASANC）

2 数据与方法

2.1　数据来源

本研究使用的数据主要包括气象数据（风速、降水、雪深和潜在蒸散发）、土壤数据、地形数据、植被数据和土地利用数据（表1）。在气象数据中，风速通过MATLAB合成自美国国家航空航天局的戈达德地球科学数据和信息服务中心（https：//disc.gsfc.nasa.gov）的第二代全球再分析数据（Modern-Era Retrospective Analysis for Research and Applications， Version 2，MERRA-2）的M2T1NXSLV数据集中2 m处的东风（U2M）和2 m处的北风（V2M）；降水量（PRECTOTCORR）来自MERRA-2的M2T1NXFLX数据集；雪深（SNODP）来自MERRA-2的M2T1NXLND数据集；潜在蒸散发来自国家青藏高原科学数据中心（https：//www.tpdc.ac.cn）《中国1 km逐月潜在蒸散发数据集》^［39-40］。地形数据和植被数据分别来自资源环境科学数据平台（https：//www.resdc.cn）的全国DEM 500 m数据（SRTM 90 m）和《中国年度NDVI、EVI 1 km数据集》。土壤数据和土地利用数据分别来自国家冰川冻土沙漠科学数据中心（http：//www.ncdc.ac.cn）的《基于世界土壤数据库（HWSD）的中国土壤数据集（v1.1）》和《中国30 m年度土地覆盖数据集及其动态变化（1985—2022年）》^［41］。

表1 模型输入参数及数据说明

Table 1 Model input parameters and data description

模型名称	输入数据	参数	时间尺度	时间分辨率	空间分辨率	数据网址
RWEQ	气象数据	风速	2001—2020年	逐时	0.5°×0.625°	https://disc.gsfc.nasa.gov
		降水	2001—2020年	逐时	0.5°×0.625°	https://disc.gsfc.nasa.gov
		雪深	2001—2020年	逐时	0.5°×0.625°	https://disc.gsfc.nasa.gov
		潜在蒸散发	2001—2020年	逐月	1 km×1 km	https://www.tpdc.ac.cn
	土壤数据	土壤质地	不变量	不变量	1 km×1 km	https://www.ncdc.ac.cn
	地形数据	DEM	不变量	不变量	250 m×250 m	https://www.resdc.cn
	植被数据	NDVI	2001—2020年	逐月	1 km×1 km	https://www.resdc.cn
	土地利用	CLCD	2001—2020年	逐年	30 m×30 m	https://www.ncdc.ac.cn

新窗口打开| 下载CSV

2.2　研究方法

2.2.1　修正风蚀方程（RWEQ）

RWEQ模型是计算土壤风蚀的经典模型。在充分考虑气侯条件、土壤可蚀性、土壤结皮、地表粗糙度、植被状况等因子的情况下，RWEQ模型采用下式估算土壤风蚀量^［42］：

S_{L} = \frac{2 x}{S^{2}} \times Q_{m a x} \times e^{- {(\frac{x}{S})}^{2}}

（1）

Q_{m a x} = 109.8 \times (W F \times E F \times S C F \times K' \times C)

（2）

Q_{x} = Q_{m a x} (1 - e^{- {(\frac{x}{S})}^{2}})

（3）

S = 150.71 \times {(W F \times E F \times S C F \times K' \times C)}^{- 0.3711}

（4）

式中：S_L为实际土壤风蚀模数，kg·m^-2；x为距上风口距离，本研究取50 m；Q_x 为距上风口x处的土壤转移量，kg·m^-1；Q_max表示风力的最大输沙能力，kg·m^-1；当转运量Q_x 达到最大转运量Q_max的63.21%时的地块长度称为关键地块长度S，m。WF表示气候因子，kg·m^-1；EF表示土壤可蚀性因子，无量纲；SCF表示土壤结皮因子，无量纲； $K'$ 表示地表糙度因子，无量纲；C表示综合植被因子，无量纲。

气候因子WF。气候因子WF反映了各类气象要素对土壤风蚀的综合影响，主要包括风因子、土壤湿度因子、雪盖因子等指标。气候因子WF的基本表达式为：

W F = W f \times \frac{ρ}{g} \times S W \times S D

（5）

式中：WF为气候因子，kg·m^-1； $W f$ 为风因子，m³·s^-3；ρ表示空气密度，本研究取1.29 kg·m^-3；g表示重力加速度，本研究取9.80 m·s^-2；SW表示土壤湿度因子，无量纲；SD表示雪盖因子，无量纲。

W f = \frac{N_{d}}{N} \sum_{i = 1}^{N} U_{i, 2} {(U_{i, 2} - U_{t, 2})}^{2}

（6）

式中： $W f$ 为风因子，m³·s^-3； $U_{t, 2}$ 为2 m高处临界起动风速，m·s^-1，设定为5 m·s^-1， $U_{i, 2}$ 为2 m高处风速，m·s^-1；当2 m高处风速低于临界起动风速时，风因子应为0；N_d 为计算时段天数，d，本研究取各月日数；N为N_d 天内的风速观测次数（i），本研究取各月时数。

S W = \frac{E T_{p} - (R + I) \frac{R_{d}}{N_{d}}}{E T_{p}}

（7）

式中：SW为土壤湿度因子，无量纲；潜在蒸散发量ET_p采用Hargreaves公式^［43］，mm；R为各月降水量，mm；I为各月灌溉量，本研究取0 mm；R_d 为降水天数和灌溉天数，d；N_d 为观测天数，d。

S D = 1 - P

（8）

式中：SD为雪盖因子，无量纲； $P$ 表示计算时段内积雪覆盖深度大于25.4 mm的概率。

土壤可蚀性因子EF。对于不同粒级的土壤颗粒来说，颗粒越粗所需的剪切力越大，土壤可蚀性也越小。而土壤中有机质、黏粒、碳酸钙等物质的存在会使土壤颗粒形成微团聚体可降低土壤可蚀性。

土壤可蚀性因子EF通过下列方程来计算：

E F = \frac{29.09 + 0.31 S a + 0.17 S i + 0.33 \frac{S a}{C l} - 2.59 O M - 0.95 C a C O_{3}}{100}

（9）

式中：EF为土壤可蚀性因子，无量纲；Sa为砂粒含量，%；Si为粉粒含量，%；Cl为黏粒含量，%；OM为有机质含量，%；CaCO₃为碳酸钙含量，%。

土壤结皮因子SCF。土壤结皮为土壤颗粒物（特别是黏粒、粉粒与有机质颗粒）的胶结作用而在土壤表面生成一层物理、化学和生物性状均较特殊的土壤微层。在风蚀过程中，结皮的存在可以减少可侵蚀颗粒的含量，降低土壤颗粒的磨蚀作用，有利于沙丘的固定，防止和减弱土壤风力侵蚀。

通过对不同黏粒和有机质含量土壤的风洞试验，Hagen^［22］建立了土壤结皮因子SCF的定量方程。

S C F = \frac{1}{1 + 0.0066 {(C l)}^{2} + 0.021 {(O M)}^{2}}

（10）

式中：SCF为土壤结皮因子，无量纲；Cl为黏粒含量，%；OM为有机质含量，%。

地表粗糙度因子K'。在RWEQ中，地表粗糙度主要是指农田因耕作导致块状土形成及土垄存在，从而使地表条件发生改变，并对土壤风蚀产生影响。

地表粗糙度因子K'采用滚轴式链条法^［44］测定，土垄糙度K_r 和土地随机糙度C_rr 采用Smith-Carson方程^［45］计算。

K' = e^{(1.86 K_{r} - 2.41 K r^{0.934} - 0.127 C r r)}

（11）

K_{r} = 0.2 \times \frac{{(∆ H)}^{2}}{L}

（12）

式中：K'为地表粗糙度因子，无量纲；K_r为土垄糙度，cm；土地随机糙度C_rr难以获取，故其值取0，cm；ΔH表示距离L范围内的高程差，cm；L为地势起伏参数，cm。

植被覆盖因子C。植被覆盖度是影响风蚀发生的关键因子，其覆盖程度直接影响近地表风速以及地表粗糙度。植被能够提高地表粗糙度，从而增大临界起沙风速，同时对移动颗粒产生一定程度的阻碍效果^［46］。植被覆盖因子C通过下列方程来计算：

C = e^{a i (S C)}

（13）

式中：C为植被覆盖因子，无量纲；SC为植被覆盖度，%；ai为不同植被类型的系数，无量纲，林地取-0.1535，草地取-0.1151，灌丛取-0.0921，裸地取-0.0768，沙地取-0.0658，耕地取-0.0438。

SC由像元二分法计算：

S C = \frac{N D V I - N D V I_{m i n}}{N D V I_{m a x} - N D V I_{m i n}}

（14）

式中： $N D V I_{m a x}$ 和 $N D V I_{m i n}$ 分别为纯植被像元和纯土壤像元的NDVI值。

2.2.2　Mann-Kendall趋势检验

Mann-Kendall非参数检验不需要样本遵循特定分布，且不受少数异常值的影响，相较于其他方法，更易计算，能够体现时间序列趋势的上升和下降，并表现出趋势变化程度^［47］，其详细计算过程为：对时间序列变量x₁，x₂，…，x_n，其检验统计量S为：

S = \sum_{i = 1}^{n - 1} \sum_{j = i + 1}^{n} s g n (x_{j} - x_{i})

（15）

s g n (x_{j} - x_{i}) = \{\begin{matrix} 1 & x_{j} - x_{i} > 0 \\ 0 & x_{j} - x_{i} = 0 \\ - 1 & x_{j} - x_{i} < 0 \end{matrix}

（16）

式中：S呈现为正态分布，其均值为0，其方差公式为：

v a r (S) = \frac{n (n - 1) (2 n + 5)}{18}

（17）

当n>10时，标准正态统计变量Z为：

Z = \{\begin{matrix} \frac{S - 1}{\sqrt[]{v a r (S)}} & S > 0 \\ 0 & S = 0 \\ \frac{S + 1}{\sqrt[]{v a r (S)}} & S < 0 \end{matrix}

（18）

当Z值大于0时，表示时间序列呈增加趋势；当Z值小于0时，表示时间序列呈减少趋势。Z值的绝对值越大，表示该时间序列的变化趋势越显著，当|Z|≥1.96、|Z|≥2.58则说明该时间序列分别通过了置信水平为95%、99%的显著性检验。|Z|≥2.58代表时间序列为极显著变化；1.96≤|Z|≤2.58代表时间序列为显著变化；|Z|≤1.96代表时间序列为不显著变化。

2.2.3　线性相关系数

皮尔逊相关系数（r）量化变量之间的线性关系^［48］。接近1的相关系数表示强正相关关系，而接近-1的系数表示强负相关关系。系数接近0表示没有线性相关性。当|r|为0.0~0.1时，表明变量间几乎不相关；当|r|为0.1~0.3时，表明变量间为弱相关；当|r|为0.3~0.5时，表明变量间为中等强度相关；当|r|为0.5~0.7时，表明变量间为中等偏强相关；当|r|为0.7~1.0时，表明变量间为强相关。

r = \frac{\sum_{i = 1}^{n} (X_{i} - \bar{X}) (Y_{i} - \bar{Y})}{\sqrt[]{\sum_{i = 1}^{n} (X_{i} {- \bar{X})}^{2}} \sqrt[]{\sum_{i = 1}^{n} (Y_{i} {- \bar{Y})}^{2}}}

（19）

式中： $X_{i}$ 是变量X的第i个样本值； $\bar{X}$ 是变量X的样本均值； $Y_{i}$ 是变量Y的第i个样本值； $\bar{Y}$ 是变量Y的样本均值。

2.2.4　土地利用转移矩阵

土地利用转移矩阵反映了某一区域在某一时段内，期初和期末各地类之间的面积转化情况，提供了土地利用变化的动态过程信息^［49］，不仅包含了特定时间点的各地类面积数据，还包括了各地类之间的面积转移情况，即期初某地类面积如何转移到其他地类，以及期末各地类面积如何从其他地类转移而来的信息。一般来说，土地利用转移矩阵的通用形式可以表示为：

S_{i j} = [\begin{matrix} S_{11} & S_{12} & \cdot \cdot \cdot & S_{1 n} \\ S_{21} & S_{22} & \cdot \cdot \cdot & S_{2 n} \\ \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot & \cdot \cdot \cdot \\ S_{n 1} & S_{n 2} & \cdot \cdot \cdot & S_{n n} \end{matrix}]

（20）

式中： $S_{i j}$ 表示从转移前的i类土地转换成转移后的j类土地的面积；n代表土地利用类型的数量；i和j（其中i和j取自1，2，…，n）分别代表转移前和转移后的土地利用类型。

2.2.5　地理探测器

空间分异性是地理现象的基本特点，地理探测器是基于风险、因子、生态和交互等探测空间分异性，从而揭示驱动因子的新的空间统计方法^［50］。本研究基于R4.4.2平台，使用Song等^［51］改进后的R语言GD包，利用其中的因子探测和因子交互作用探测发掘中国北方干旱半干旱地区土壤风蚀的驱动因素。具体而言，本研究以中国北方干旱半干旱区实际土壤风蚀模数（Y）为因变量，以气候因子（X₁）、砂粒质量分数（X₂）、粉粒质量分数（X₃）、黏粒质量分数（X₄）、有机质质量分数（X₅）、碳酸钙质量分数（X₆）、海拔（X₇）、坡度（X₈）、土地利用类型（X₉）、归一化植被指数（X₁₀）为自变量。在ArcGIS10.8中，首先使用“创建渔网”工具生成20 km×20 km的渔网点；其次使用“按位置选择”工具以研究区边界为掩膜，裁剪出相应的渔网点；然后使用“多值提取到点”工具，将各因子进行采样；最后将采样点的属性表进行导出，每期生成约1.3万个采样点。在SPSS27.0.1.0中，使用“条件筛选”，剔除各因子的空值，每期剩余约1.1万个采样点。在R4.4.2中，导入筛选后的数据，选择等间距、自然间断、分位数、几何间隔、标准差等5种数据离散化方法，再选择4~6的分类数量，从而运行参数最优地理探测器。

空间分异与因子探测

q = 1 - \frac{\sum_{h = 1}^{L} N_{h} σ_{h}^{2}}{N σ^{2}} = 1 - \frac{S S W}{S S T}

（21）

S S W = \sum_{h = 1}^{L} N_{h} σ_{h}^{2}, S S T = N σ^{2}

（22）

式中：h=1，2，…，L为变量Y或因子X的分类； $N_{h}$ 和N分别为类h和全区的单元数； $σ_{h}^{2}$ 和 $σ^{2}$ 分别为类h和全区的Y值的方差；SSW和SST分别为类内方差和以及全区方差和。q的取值为［0，1］，值越大表明变量Y的空间分异性越强；如果分类是由解释变量X生成的，则q值越大表明解释变量X对响应变量Y的解释力越强，反之则越弱。

因子交互作用探测。交互作用探测用于识别因子之间的交互作用，即评估各解释变量共同作用和相互独立作用对响应变量的解释力（表2）。

表2 因子交互作用的类型^［50］

Table 2 Types of interaction between two covariates^［50］

交互类型	判定依据
非线性减弱	$q (X_{1} ⋂ X_{2}) < M i n (q (X_{1}), q (X_{2}))$
单因子非线性减弱	$M i n (q (X_{1}), q (X_{2})) < q (X_{1} ⋂ X_{2}) < M a x (q (X_{1}), q (X_{2}))$
双因子增强	$q (X_{1} ⋂ X_{2}) > M a x (q (X_{1}), q (X_{2}))$
独立	$q (X_{1} ⋂ X_{2}) = q (X_{1}) + q (X_{2})$
非线性增强	$q (X_{1} ⋂ X_{2}) > q (X_{1}) + q (X_{2})$

新窗口打开| 下载CSV

3 结果与分析

3.1　土壤风蚀的空间分布格局

2001—2020年，中国北方干旱半干旱区受土壤风蚀影响的面积约为4.47×10⁶ km²，约占总面积的79.90%；多年平均实际土壤风蚀模数为13 449.59 t·km^-2·a^-1。研究区土壤风蚀的等级以剧烈侵蚀（实际土壤风蚀模数>15 000 t·km^-2·a^-1）为首，约占研究区面积的24.59%；其次为轻度侵蚀（实际土壤风蚀模数200~2 500 t·km^-2·a^-1）和微度侵蚀（实际土壤风蚀模数<200 t·km^-2·a^-1），分别约占研究区面积的19.60%和10.83%。从空间分布来看，剧烈侵蚀横贯中国北方干旱半干旱地区，主要集中在塔里木盆地中心地带、青藏高原西部、内蒙古高原中西部和黄土高原中部；微度侵蚀和轻度侵蚀主要分布在中国北方干旱半干旱地区南部、新疆维吾尔自治区的天山北坡地区和塔里木盆地的边缘地区以及黑龙江省北部；极强烈侵蚀（土壤风蚀模数8 000~15 000 t·km^-2·a^-1）、强烈侵蚀（土壤风蚀模数5 000~8 000 t·km^-2·a^-1）、中度侵蚀（土壤风蚀模数2 500~5 000 t·km^-2·a^-1）大致分布在剧烈侵蚀向轻度侵蚀的过渡地区（图2）。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 2001—2020年中国北方干旱半干旱区多年平均实际风蚀等级空间分布

注：基于自然资源部标准地图服务网站标准地图（审图号：GS（2024）0650号）制作，底图边界无修改

Fig.2 Spatial distribution of multi-year average actual soil wind erosion modulus in the ASANC from 2001 to 2020

3.2　土壤风蚀的时空动态变化

2001—2020年，中国北方干旱半干旱区实际土壤风蚀模数总体呈现出波动下降态势，大致以每年102.78 t·km^-2的速率下降。其中，2001年实际土壤风蚀模数最大（16 281.30 t·km^-2·a^-1）；2019年的实际土壤风蚀模数最小（10 774.35 t·km^-2·a^-1，图3A）。2001—2020年，中国北方干旱半干旱区不同侵蚀强度的土地面积占比呈现波动变化的状况。具体而言，微度侵蚀、轻度侵蚀和中度侵蚀面积占比呈现波动上升的状况，强烈侵蚀、极强烈侵蚀和剧烈侵蚀面积占比呈波动下降趋势。其中，轻度侵蚀面积占比由2001年的19.32%波动上升至2020年的25.65%，占比增加了6.33%；剧烈侵蚀面积占比由2001年的36.47%波动下降至2020年的28.45%，占比减少了8.02%（图3B）。2001—2020年，中国北方干旱半干旱区实际土壤风蚀模数的变化趋势以减少为主，不显著减少、显著减少和极显著减少的累积面积约占研究区总面积的63.72%。其中，不显著减少区域面积占比为56.37%；显著减少区域面积占比为6.30%；极显著减少区域面积占比为1.05%，主要分布在黄土高原地区。此外，极显著增加、显著增加和不显著增加累积面积约占研究区总面积的17.22%。其中，极显著增加区域面积占比为0.01%；显著增加区域面积占比为0.28%；不显著增加区域面积占比为16.93%，主要分布在塔里木盆地的边缘地区、青藏高原西北部和南部地区以及大兴安岭西部和小兴安岭中部（图3C）。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 2001—2020年中国北方干旱半干旱区实际土壤风蚀模数的时空动态变化

注：基于自然资源部标准地图服务网站标准地图（审图号：GS（2024）0650号）制作，底图边界无修改

Fig.3 Spatiotemporal dynamics of actual soil wind erosion modulus in the ASANC from 2001 to 2020

3.3　土壤风蚀的驱动因素

3.3.1　土地利用变化对土壤风蚀的影响

2001—2020年，土地利用转移面积约为79 636 km²，约占研究区面积的1.42%。耕地和沙地的转出面积分别占总转出面积的20.92%、46.42%；建设用地和林地的转入面积分别占总转入面积的47.73%、36.95%。此外，28 868 km²的耕地被转化为建设用地，成为建设用地的主要来源，分布在京津冀地区中部和南部、黑龙江省西部、吉林省西北部和辽宁省北部。同时，林地增加了29 424 km²，其中最主要的来源是草地，其次为耕地（图4）。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 2001—2020年中国北方干旱半干旱区土地利用转移弦图

Fig.4 Landuse transfer chord diagram in the ASANC from 2001 to 2020

2001—2020年，发生转换的土地利用类型在2020年的实际土壤风蚀模数的排序：沙地（15 556.17 t·km^-2·a^-1）>草地（13 695.88 t·km^-2·a^-1）>耕地（7 933.25 t·km^-2·a^-1）>建设用地（5 024.63 t·km^-2·a^-1）>林地（3 896.55 t·km^-2·a^-1）；未发生转换的土地利用类型在2020年的实际土壤风蚀模数的排序：沙地（19 150.57 t·km^-2·a^-1）>草地（11 422.58 t·km^-2·a^-1）>耕地（5 004.18 t·km^-2·a^-1）>建设用地（4 271.00 t·km^-2·a^-1）>林地（2 460.38 t·km^-2·a^-1，图5）。研究表明，中国北方干旱半干旱区新增草地和新增耕地的实际土壤风蚀模数相对较高，分别为13 695.88 t·km^-2·a^-1和7 933.25 t·km^-2·a^-1，比研究时段内未发生转换的草地和耕地的实际土壤风蚀模数分别高19.90%和58.53%。通过进一步研究发现，新增耕地主要是由沙地转换而来，故该种土地利用类型实际土壤风蚀模数的增幅相对较大；而新增草地的主要来源为沙地，其次为耕地，故该种土地利用类型实际土壤风蚀模数的增幅相对较小（图4C）。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 2001年和2020年中国北方干旱半干旱区不同土地利用类型的实际土壤风蚀模数

Fig.5 The actual soil wind erosion modulus of different landuse in the ASANC in 2001 and 2020

3.3.2　土壤风蚀各影响因子的交互作用

本研究以实际土壤风蚀模数为响应变量，以气候因子、海拔、砂粒质量分数、土地利用类型和归一化植被指数等10个地理因子为解释变量，对中国北方干旱半干旱区2001—2020年的节点年份（2001、2005、2010、2015、2020年）进行单因子探测和因子交互探测。研究表明，各因子对实际土壤风蚀模数的解释力均有显著影响（P<0.01）。总体而言，2001—2020年实际土壤风蚀模数影响显著的因子按q值的均值排序为：气候因子（X₁）>黏粒质量分数（X₄）>砂粒质量分数（X₂）>粉粒质量分数（X₃）>归一化植被指数（X₁₀）>有机质质量分数（X₅）>土地利用类型（X₉）>海拔（X₇）>坡度（X₈）>碳酸钙质量分数（X₆）。在全部因素中，气候因子（X₁）对实际土壤风蚀模数的空间差异解释力最强（q=0.312），起主导作用；在土壤因素中，黏粒质量分数（X₄）对实际土壤风蚀模数的空间差异解释力最强（q=0.240），碳酸钙质量分数（X₆）对实际土壤风蚀模数的空间差异解释力最弱（q=0.040）；在其他下垫面因素中，NDVI（X₁₀）的解释力最强（q=0.127），土地利用类型（X₉）的解释力次之（q=0.086，图6A）。具体而言，2001—2020年各因子对实际土壤风蚀模数空间差异解释力（q值）的变化趋势各异，均呈现出波动变化的态势。对研究区实际土壤风蚀模数解释力下降的因子有粉粒质量分数（X₃）和坡度（X₈），其余因子均对研究区实际土壤风蚀模数解释力呈现上升趋势。其中，解释力上升最为明显的是海拔（X₇），从2001年的5.30%上升至2020年的11.80%；解释力上升最为微弱的是土地利用类型（X₉），从2001年的8.64%上升至2020年的8.76%。此外，从不同研究时段主导因子的位次来看，2001—2020年气候因子（X₁）始终占据第一主导因子的位置，碳酸钙质量分数（X₆）的解释力始终处于低位。砂粒质量分数（X₂）、粉粒质量分数（X₃）、黏粒质量分数（X₄）和有机质质量分数（X₅）对实际土壤风蚀模数的解释力在研究时段内的变化趋势一致，在2001—2005年均呈现上升趋势，在2005—2010年均呈现下降趋势，在2010—2015年均呈现上升趋势，在2015—2020年均呈现下降趋势，且与气候因子（X₁）对实际土壤风蚀模数解释力的变化呈现相反的状况（图6B）。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 2001—2020年中国北方干旱半干旱区实际土壤风蚀模数的单因子探测

注：X₁代表气候因子，X₂代表砂粒质量分数，X₃代表粉粒质量分数，X₄代表黏粒质量分数，X₅代表有机质质量分数，X₆代表碳酸钙质量分数，X₇代表海拔，X₈代表坡度，X₉代表土地利用类型，X₁₀代表归一化植被指数

Fig.6 Factor detection of the actual soil wind erosion modulus in the ASANC from 2001 to 2020

交互探测的结果表明（图7），研究区内的交互类型为非线性增强和双因子增强两种交互类型，即影响因子两两交互均会增强对实际土壤风蚀模数的解释力。交互作用解释力排在首位的是气候因子（X₁）和黏粒质量分数（X₄）的非线性增强作用，是土壤风蚀的显著控制因子，解释力达到62%。同时，气候因子（X₁）与所有因子交互作用后对实际土壤风蚀模数的影响均大于气候因子（X₁）的单独作用（q=0.31）。此外，相比于其他影响因子，碳酸钙质量分数（X₆）、海拔（X₇）、坡度（X₈）以及土地利用类型（X₉）单独对实际土壤风蚀模数的影响均相对较小，但与气候因子（X₁）进行交互后，解释力均上升至33%以上。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 2001—2020年中国北方干旱半干旱区实际土壤风蚀模数的交互探测

Fig.7 Interaction detection of the actual soil wind erosion modulus in the ASANC from 2001 to 2020

4 讨论

4.1　研究的可靠性

本研究采用中国生态环境部有关防风固沙评估模型，利用修正风蚀方程模拟了研究时段内中国北方干旱半干旱区的土壤风蚀模数。本研究结果表明2001—2020年中国北方干旱半干旱区的多年平均实际风蚀模数为13 449.59 t·km^-2·a^-1，该数值与前人研究^［52-54］相比略有差异，但鉴于输入数据源和尺度效应的影响，依然具备一定的参考价值。为了进一步验证模型结果的准确性，本文收集了新疆准东地区利用¹³⁷Cs示踪法^［55］和粒度对比法^［56］推算的土壤风蚀模数的研究成果，与本研究模拟的结果进行对比（表3）。结果显示，本研究模型反演的实际土壤风蚀模数略微偏大。通过进一步分析，通过MERRA-2合成的风速数据较之站点实测数据明显偏大，这直接影响了RWEQ的模拟效果。此外，在相关区域不同地类实际土壤风蚀模数的排序与前人^［57-59］的研究一致，均为沙地>草地>耕地>建设用地>林地。

表3 模型模拟验证结果对比^［55-56］

Table 3 Comparison results of RWEQ with measurement^［55-56］

研究区域	研究时段	研究方法	纬度	经度	原文结果 /(t·km^-2·a^-1)	本文结果 /(t·km^-2·a^-1)
新疆准东地区	2015年	¹³⁷Cs示踪法	44°16′14″N	90°21′00″E	2 644.08	3 678.03
			45°06′37″N	88°35′20″E	3 855.57	5 266.36
			44°29′58″N	89°54′54″E	1 871.13	4 340.87
			44°47′12″N	89°21′06″E	1 840.09	1 197.04
	2016年	粒度对比法	44°06′34″N	90°00′55″E	2 585.58	2 842.68
	2016年	粒度对比法	45°00′05″N	89°14′04″E	2 857.28	2 199.54

新窗口打开| 下载CSV

4.2　研究的局限性

本研究基于RWEQ模型对中国北方干旱半干旱区的土壤风蚀状况进行了定量模拟，同时与模型数据和¹³⁷Cs实测数据进行了对比验证，取得了一定研究成果，但仍存在一些不足需要进一步改进。

风速数据精度的提升。风力作为土壤风蚀发生与演进的核心驱动力^［15］，无论是在基于经验统计抑或基于过程机理的土壤风蚀模型中，风力侵蚀因子始终占据着模型构建的关键地位。因此，具备高时间分辨率与高空间分辨率的风力数据对于土壤风蚀模型中土壤风蚀模数的精确计算至关重要。然而，在实际风速观测过程中，受限于观测设备性能及观测条件等客观因素，基于大尺度空间范围获取高时间分辨率的风速数据颇具挑战。今后研究可以利用多模态信息融合来改进风速反演模型^［60］，结合机器学习或降尺度方法获得高分辨率风场^［61］等，进而为土壤风蚀模型提供更高精度的风速数据。

风蚀模型参数的修订。土壤风蚀模型大多源自国外，由于国内外土壤特性、气候条件和农田管理措施的差异，现有模型的参数需要修订以提高适用性^［62］。未来研究应重点关注根据中国本土土壤特性、气候条件和农业管理方式修订模型参数^［63］，特别在土壤粒径分类、作物覆盖度、耕作方式等方面^［64］。同时，需要验证模型参数在局地与大尺度区域中的适用性，并通过风洞试验探索基于实时数据动态调整模型参数的机制。

5 结论

本研究的主要结论如下：在研究时段内，研究区的多年平均实际土壤风蚀模数为13 449.59 t·km^-2·a^-1，以剧烈侵蚀为首，主要分布在塔里木盆地中心地带和内蒙古高原中西部。实际土壤风蚀模数总体呈现出波动下降的态势，大致以每年102.78 t·km^-2的速率下降。研究区不同土地利用类型的实际土壤风蚀模数存在显著差异。无论是发生转换的土地利用类型，还是未发生转换的土地利用类型，2020年的实际土壤风蚀模数：沙地>草地>耕地>建设用地>林地。此外，中国北方干旱半干旱区新增草地和新增耕地的实际土壤风蚀模数相对较高，分别为13 695.88 t·km^-2·a^-1和7 933.25 t·km^-2·a^-1，比研究时段内未发生转换的草地和耕地的实际土壤风蚀模数分别高19.90%和58.53%。研究区土壤风蚀的主导因子为气候因子（q=0.312， P<0.01）。其中，风因子与实际土壤风蚀模数呈现出较为显著的中等强度的正相关（r=0.43， P<0.05）。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

C J

， Fu

B J

， Wang

，et al.

Drivers and impacts of changes in China's drylands

［J］.Nature Reviews Earth & Environment，2021，2（12）：858-873.

[本文引用: 2]

[2]

Borrelli

， Panagos

， Alewell

，et al.

Policy implications of multiple concurrent soil erosion processes in European farmland

［J］.Nature Sustainability，2023，6（1）：103-112.