基于USLE模型的祁连山国家公园土壤水力侵蚀评价
陈豪
,
丁文广
,
Tanjia Binte Zafar
收稿日期：2020-02-29; 修回日期：2020-03-17
项目名称：中国科学院战略性先导科技专项"祁连山生态系统变化归因与善治对策"(XDA20100102)
摘要：基于通用土壤流失方程(USLE)模型对2005、2010和2015年祁连山国家公园土壤水力侵蚀情况进行计算和分析, 为本地区的土壤侵蚀防治提供决策建议。经计算方法选择、插值方法对比、同类研究类比、统计分析, 获得每年及平均土壤水力侵蚀空间强度与等级分类。结果表明: 2005—2015年, 祁连山国家公园的土壤水力侵蚀水平有所下降, 土壤侵蚀量由2005年的6 314.26万t/a下降至2015年的5 755.59万t/a; 就土壤侵蚀等级来说, 研究区范围内70%的区域属于微度侵蚀, 整体上土壤侵蚀强度水平较低; 土壤水力侵蚀主要发生在高程为≥3 600~4 600 m的区域和≥8°~15°的缓坡, 而高程≥4 600 m的区域贡献25.67%的侵蚀量, 如果冰川积雪退缩, 裸土面临暴露, 可能加重土壤侵蚀。因此在国家公园未来发展与土壤侵蚀防治过程中, 需要严格保护3 600~4 600 m的高程带, 以及5°~15°缓坡地带的自然环境, 同时保护冰川与永久积雪。
关键词：USLE GIS 土壤水力侵蚀 祁连山国家公园
Evaluation of soil water erosion in Qilian Mountain National Park based on USLE model
College of Earth and Environmental Sciences, Lanzhou University, 730000, Lanzhou, China
土壤流失通用方程USLE是在多个国家和地区已经验证的土壤侵蚀模型，具有较为便捷，评估量化等优点。在我国，对于西北地区的土壤侵蚀研究在20世纪70年代已经开始，最典型的是围绕天水等地进行土壤侵蚀径流小区的研究，形成了对USLE的修正模型[1]。
祁连山地区是西北重要的生态安全屏障，1988年在国务院的批准下，于该地区建立了祁连山国家级自然保护区，但在随后的30年中保护区生态安全状况总体恶化[2-3]，直到2018年祁连山国家公园管理机构正式成立，地区环境保护进入新的阶段。对国家公园过往的土壤水力侵蚀状况进行分析有助于为地区土壤侵蚀防治提供依据，并有补充国家公园相关研究空白的意义，可对本地区的土壤侵蚀防治提供决策建议。
1 区域背景
国家公园最初的建设目的是保护珍稀的生态系统、自然景观以及文化遗产[4]，随后国家公园的理念在世界范围产生了广泛影响。2017年9月，中共中央办公厅、国务院印发了《建立国家公园体制总体方案》，2018年，祁连山国家公园管理机构成立[5]。
祁连山国家公园位于青海省与甘肃省交界，青藏高原北缘，海拔1684~5604m。其南部为青藏高原，北部为河西走廊。国家公园总体为西北至东南走向，区域狭长，总面积5.02万km2，其中甘肃片区3.44万km2，青海片区1.58万km2。从土壤侵蚀类型分区的角度，祁连山国家公园处于风力侵蚀类型区、冻融侵蚀类型区和水力侵蚀类型区的交汇处。
2 数据与方法
2.1 数据来源
本研究中所采用的数据包括：由资源环境数据云平台(http://www.resdc.cn)提供的1km空间分辨率的数字高程模型(DEM)、2005、2010、2015年1km空间分辨率的土地利用数据及1:100万全国土壤类型属性数据，土地利用数据解译精度大于90%[6]，具有较高的精度；由中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)提供的由TREEA卫星获得的MODIS植被指数旬合成产品，空间分辨率为500m。气象数据由中国气象数据共享网(https://data.cma.cn/)以及国家气象信息中心(http://www.nmic.cn/)获得的研究区周围现有28个站点的日值降水数据。
2.2 模型概述
本研究应用由Wischmeier[7]提出的通用土壤流失方程(USLE)为基础，并逐个选择合适的因子计算方法，以修正后的结果计算土壤水力侵蚀模数。USLE模型公式为
式中：A为土壤侵蚀模数，t/(hm2·a)，R为降水侵蚀力因子，MJ·mm/(hm2·h·a)；K为土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm)；LS为地形因子(L为坡长因子，S为坡度因子)，量纲为1；C为植被覆盖与作物管理因子，量纲为1；P为水土保持措施因子，量纲为1。
2.3 土壤侵蚀分类标准
根据SL190—2007《土壤侵蚀分类分级标准》[8]以及其他相似研究区做法，如祁连山南坡、三江源地区、大兴安岭冻融侵蚀区、大别山区等[9-12]，将土壤侵蚀强度分级标准确定为微度侵蚀≤500t/(km2·a)、轻度侵蚀(500, 2500]t/(km2·a)、中度侵蚀(2500, 5000]t/(km2·a)、强烈侵蚀(5000, 8000]t/(km2·a)、极强烈侵蚀(8000, 15000]t/(km2·a)、剧烈侵蚀＞1万5000t/(km2·a)。不同坡度带分级标准确定为[0~5°)、[5°~8°)、[8°~15°)、[15°~25°)、[25°~35°)、[35°, 90°)6个等级。高程带划分考虑到祁连山国家公园所涉及高程在1600~5600m之间，故将其分均分为[1600, 2600)m、[2600, 3600)m、[3600, 4600)m、≥4600m 5个范围。
3 土壤水力侵蚀模数计算
3.1 降水侵蚀力因子R
本研究采用Wischmeier在1978年提出的简化算法[13]，该方法需要年降水量和月降水量，均衡了需要的数据量和对全年降水分布的反映。其具体公式为
|
$
R=\sum\limits_{i=1}^{12} 1.735 \times 10^{\left(1.5 \lg \left(P_{i} / P\right)-0.8188\right)}。$
|
(2) |
式中：Pi为月降水量，mm；P为年降水量，mm。抽取研究区内不相邻的3个站点，选用不同插值方法模拟降水量，与实测值相较选择相对平均误差MRE(mean relative error)最小的反距离权重法对式(2)结果进行插值，各插值方法对比结果见表 1。
表 1(Tab. 1)
表 1 各插值方法对比结果
Tab. 1 Compared results of interpolation methods
| 方法Methods |
MRE |
| 普通克里金Ordinary kriging |
4.26 |
| 泛克里金Universal kriging |
11.79 |
| 反距离权重Inverse distance weighting |
3.19 |
| 样条函数Splines |
6.08 |
| 自然邻域Natural neighbor |
6.22 |
| Notes: MRE stands for mean relative error. |
|
表 1 各插值方法对比结果
Tab. 1 Compared results of interpolation methods
|
3.2 土壤土壤可蚀性因子K
本研究区的土壤可蚀性因子K采用Williams等在EPIC模型中提出的计算公式[14]，此公式需要提供土壤颗粒组成以及土壤有机碳质量分数。通过1:100万全国土壤类型数据，与土壤属性数据库说明，可以获得各类土壤上述组分资料。
|
$
\begin{array}{c}
K=(0.2+0.3 \exp (-0.0256 \mathrm{SAN}(1-\mathrm{SIL} / 100))) \cdot \\
\left(\frac{\mathrm{SIL}}{\mathrm{CAL}+\mathrm{SIL}}\right)^{0.3}\left(1.0-\frac{0.25 C}{C+\exp (3.72-2.95 C)}\right) \cdot \\
\left(1.0-\frac{0.7 \mathrm{SN}_{1}}{\mathrm{SN}_{1}+\exp \left(-5.51+22.9 \mathrm{SN}_{1}\right)}\right)
\end{array}
$
|
(3) |
式中：SAN为砂粒质量分数，%；SIL为粉粒质量分数，%；CAL为黏粒质量分数，%；C为土壤有机碳质量分数，%；SN1=1-SAN/100。通过1:100万全国土壤类型数据，确定研究区内的土壤类型分布，计算不同类型土壤的砂粒、粉粒、黏粒质量分数。
3.3 地形因子LS
本研究中应用基于黄土高原侵蚀情况，经过改进的坡长因子和坡度因子[15]公式计算研究区每个位置的地形因子。此算法为《第4次全国土壤侵蚀普查技术规程》所采用的计算方法，由于其坡度适用区间较大，而祁连山国家公园普遍存在较大坡度的地形，故适合研究区使用。具体公式[16]为：
|
$
S=\left\{\begin{array}{ll}
10.8 \sin \theta+0.03 & \theta <5^{\circ} \\
16.8 \sin \theta-0.50 & 5^{\circ} \leqslant \theta<10^{\circ} \\
21.9 \sin \theta-0.96 & \theta \geqslant 10^{\circ}
\end{array}\right.;
$
|
(4) |
|
$
L=\left(\frac{\lambda}{22.1}\right)^{n}, n=\left\{\begin{array}{ll}
0.2 & \theta \leqslant 1^{\circ} \\
0.3 & 1^{\circ} <\theta \leqslant 3^{\circ} \\
0.4 & 3^{\circ} <\theta \leqslant 5^{\circ^{\circ}} \\
0.5 & \theta>5^{\circ}
\end{array}\right.
$
|
(5) |
式中：θ为坡度，(°)；n为坡长指数；量纲为1；λ为坡长，m。采用分辨率为1km的DEM数据，提取计算坡长因子和坡度因子：先对影像进行填洼，再得到径流方向(角度)，根据流向计算单元栅格坡长λ；然后使用ENVI软件中的band math工具，利用分段函数对式(4)、(5)分别进行计算，得到坡长因子L和坡度因子S，相乘得到地形因子LS。
3.4 植被覆盖与作物管理因子C
本研究应用线性混合像元分解法计算植被覆盖因子[17]，C值的估算公式为：
|
$
\left\{\begin{array}{ll}
C=1 & G=0 \\
C=0.6508-0.3436 \lg G & 0 <G <78.3 \% \\
C=0 & G \geqslant 78.3 \%
\end{array}\right.。$
|
(6) |
式中G为植被覆盖度，%，植被覆盖率与归一化植被指数(normalized difference vegetation index, NDVI)的关系可以用下式表示：
式中：N为像元NDVI值；Nsoil为裸地像元NDVI值；Nveg为高纯度植被像元NDVI值。C值计算过程为，先对所获取的MODIS植被指数影响进行提取，生成全年各旬NDVI影响，再对其进行最大化合成，获得作物生长期NDVI。在ArcGIS中计算植被覆盖度G，而后根据式(6)进行分段函数计算，获得植被覆盖与作物管理因子C分布。
3.5 水土保持措施因子P
水土保持因子P赋值方法有2种，一种为根据土地利用类型直接赋值[18]，一种为根据经验公式进行计算赋值[19]，经验公式为
式中$\phi $为坡度，%。本研究采用赋值和公式计算结合的方法，对水体、岩石、冰川、建设用地等地类的P值赋值为0，其他地类根据式(8)进行计算。
4 结果与分析
4.1 土壤侵蚀量时空分布统计分析
通过式(2)~(8)计算式(1)ULSE中6个因子得到2005、2010、2015年祁连山国家公园区域土壤侵蚀模数空间分布，2005年，研究区平均侵蚀模数为1257.82t/(km2·a)，2010年平均侵蚀模数为1094.20t/(km2·a)，2015年平均侵蚀模数为1146.53t/(km2·a)。研究区各年年均侵蚀模数为1166.90t/(km2·a)，年均侵蚀量为5857.86万t/a。根据土壤侵蚀强度分类标准划分对应等级(图 1)，可以看出侵蚀重点出现在国家公园中部地区，各年均保持这一态势。
从各侵蚀强度等级面积比例(图 2)来看，研究区各年份70%面积处于微度侵蚀等级，10%~20%属于轻度侵蚀，中度侵蚀面积比例在7%左右，其他等级均在5%以下。其中2010年微度等级面积比例最大，其他更高等级面积比例均为最小，反映出土壤侵蚀总体强度的降低，2015年轻度侵蚀面积比例出现回升，但总体而言仍比2005年有改善。
从总侵蚀量变化(图 3)来看，除微度侵蚀外，各年各等级侵蚀量呈现U形变化，2005年最高，2010年最低，2015年有少量回升，但仍比2005年下降了约10%。2010年微度侵蚀强度对应侵蚀量最高，是因为更高侵蚀强度等级向微度等级转换后，此等级面积明显增加，累计侵蚀量也有所增加。总体上看，轻度、中度、强烈、极强烈4个等级构成了侵蚀量的主体，其中极强烈等级侵蚀量最大(1500万t左右)，而其仅有3%左右的等级对应面积。
4.2 不同地形因素下土壤水力侵蚀时空分布
在不同高程带的侵蚀量比例中，比例最高的高程带是≥3600~4600m，这是由于祁连山区域降水受坡度海拔影响较大。由于祁连山地区雪线位于4673m(2012年)[20]，因此在高程≥4600m，仅占整个研究区8.10%面积的区域发生总侵蚀量的25.67%，达到1539.59万t(图 4)，在这个高程带植被覆盖度较低或完全没有植被覆盖，常存在较大坡度，导致了较高的土壤侵蚀模数。因此当冰川和积雪层出现萎缩，则会有更多裸地暴露出来，加剧土壤侵蚀的情况。
从坡度等级来说，不同坡度等级所占面积随坡度等级上升而下降，侵蚀量主要集中在0~25°的地带，而平均侵蚀模数(等级对应侵蚀量与面积之商)则随着坡度上升而上升(图 5)。在1km尺度上没有出现≥35°的地带。在分级侵蚀量方面，[8°, 15°)的地带所占侵蚀量比例达到44.36%，为2598.55万t，是土壤侵蚀的最大构成部分。[15°, 25°)区间土壤侵蚀量为938.27万t，[25°, 35°)的陡坡区域侵蚀量仅为18.31万t，这2个坡度区间侵蚀量虽然很少，但侵蚀模数为研究区最高水平，大于2900t/(km2·a)。
5 结论与讨论
从2005年至2015年，祁连山国家公园的土壤水力侵蚀水平有所下降，总侵蚀量下降约10%，各侵蚀强度等级有向低等级转移的趋势，反映出地区土壤侵蚀情况有所缓解。研究区近85%的面积均处于微度和轻度侵蚀强度，近年来比例有所上升，总体上国家公园土壤侵蚀处于较低水平。侵蚀主要集中在面积比例小于10%的强烈与极强烈等级，因此在未来土壤侵蚀防治工作中，应有针对性的对重点侵蚀区进行治理。从地形因素上看，≥3600~4600m高程带以及≥8°~15°的缓坡地区是侵蚀发生的重点地区。国家公园成立之后，应当依法依规停止公园范围内的生产生活活动，以保持并恢复这类地区的植被覆盖，从而降低侵蚀模数。
特别需要注意的是，≥4600m高程带仅占总面积的8%，但此高程带的平均侵蚀模数接近3700t/(km2·a)，在各高程带中最高。祁连山地区的雪线位于4673m(2012年，上升速率4.23m/a)，如果冰川萎缩、雪被消融，则可能暴露出更多裸土，造成土壤侵蚀的加重。
基于ULSE的土壤水力侵蚀研究的先例中，有许多与本研究类似的将研究区选择在典型水力侵蚀区之外。祁连山南坡、三江源地区、大兴安岭冻融侵蚀区、大别山区等[8-11]，以上研究区属于高海拔地区、高寒地区、或者山区，与祁连山国家公园有许多可以类比的相似之处，在上述地理位置、地形、气候、研究年份等条件相似的地区，利用USLE模型计算的侵蚀模数处于1000~3000t/(km2·a)范围之内，而本文所获得的结果为研究区平均侵蚀模数1166.90t/(km2·a)。故而说明本研究结果具有一定可信度。
2020, Vol. 18 
