夏兴生，田绍梅，吕圣慧，潘耀忠

(1.青海省人民政府—北京师范大学 高原科学与可持续发展研究院，青海 西宁 810016；2.北京师范大学 地理科学学部 遥感科学国家重点实验室，北京 100875；3.青海师范大学 地理科学学院，青海 西宁 810016)

1 引言

土壤是支持生态系统产品供给和功能运转的基础，因此，土壤保持能力的大小对人类从生态系统中获得惠益至关重要[1].青藏高原作为世界海拔最高的独特地理单元，也作为东亚的生态屏障，其土壤保持能力所创造的经济价值占其农林牧业收入总量的10%左右[2]，但受全球气候变化和过渡的人类活动影响，青藏高原部分地区的土壤保持能力处于持续下降的状态，威胁着生态系统服务功能的正常运转和社会经济的稳定发展[3].因此，关注青藏高原生态系统土壤保持能力的变化，对当地生态环境的修复和养护策略的制定具有重要的参考意义.

目前，针对高原土壤保持研究评估大致可分为基于模型的定量评估和基于土壤保持因子，如气候变化、植被覆盖等的定性评估.具体用于土壤保持研究的定量模型主要有USLE模型[4]、In VEST模型[5]等[6-8].在基于土地利用/植被覆盖度、气候变化等定性评价土壤保持服务方面，Jackson等[9]基于地形因素对土壤侵蚀影响的考虑，提出复合地形评价指数来识别易受土壤侵蚀的脆弱区；朱青等[10]采用趋势分析法，重点分析了黄土高原植被覆盖状况对土壤保持的作用，这些研究为后续的深入探讨奠定了方法论基础.

青藏高原因为兼具生态屏障和生态服务的双重功能，针对其土壤保持研究也有丰富的成果.例如，于格等[11]基于风洞实验方法，定量研究了青藏高原高寒草甸土壤保持能力的经济效益；朱永刚等[12]针对青藏高原生态系统脆弱、环境易破坏难恢复的问题，总结了青藏高原水土保持养护规划的重要技术和管理措施；闻亮等[13]基于In VEST模型采用叠加分析法客观评价了三江源土壤保持能力的时空变化.但是，目前对于青藏高原土壤保持仍然集中在特定时空下的静态研究方面[14]，基于不同因子变化下的土壤保持能力动态变化关注较少，导致驱动土壤保持能力时空变化的机理和核心要素不明晰.因此，有必要在三江源自然保护区和国家公园建设的背景下，进一步结合多源空间数据，探讨多因子共同作用下的青藏高原典型区土壤保持能力的时空变化，为当地的土壤侵蚀防治工作和生态环境的修复工作提供参考.

本研究选择三江源区域的玛多县为研究区，结合2000年、2010年、2020年三期土地覆盖、降雨及植被指数等影响土壤保持能力的数据产品，基于In VEST模型的土壤保持模块进行潜在和实际土壤流失量的估算，分析玛多县土壤保持能力的时空变化特征及植被和降水的驱动作用，以期为研究区生态保护与修复措施的制定和实施提供参考.

2 研究区与数据

2.1 研究区概况

玛多县位于青海省果洛藏族自治州西北部，地处青藏高原腹地，全县总面积约2.53万km2，平均海拔在4200m，年平均气温-4℃，属高寒草原气候(图1)[15].玛多县也占据了三江源自然保护区、国家公园黄河源园区的大部分地区，是高寒生态脆弱地带和国家重点生态功能区的核心区域[16].高山草甸土、高山草原土组成了该县的主要土壤类型，二者皆为天然的草甸、草原生态系统发育提供了良好的支撑[17，18].但是，由于高原地理环境的影响，该区的土壤成土过程缓慢，土层薄，保水性能差，发育的草种一般都长势矮小、鲜草产量比较低，土壤-草地耦合的脆弱生态系统对高原严峻的气候条件和极端天气/气候事件的响应敏感[19，20]，因此，该区域的土壤保持能力变化对生态系统变化、气候变化也具有一定的指示作用，关注这一变化对推进三江源国家公园的生态建设，筑牢我国的生态屏障和核心水源涵养地的保护具有重要的作用.

图1 研究区概况

2.2 数据及预处理

本研究所用的数据主要包括土地利用数据、数字高程模型(DEM)数据、土壤数据、气象数据、遥感反演的植被指数数据、基础地理数据，具体及预处理如下：

(1)土地利用数据：来源于中国自然资源部(https：//www.webmap.cn/commres.do?method=globeIndex)发布的Globeland30数据产品，时间范围为2000年、2010年、2020年，主要用于提取研究区的水土保持措施因子.预处理主要是将前期下载好的数据进行拼接并按照地类代码转为栅格图层.

(2)DEM数据：来源于地理空间数据云(http：//www.gscloud.cn)，空间分辨率为30m，用于In VEST模型参数的提取.通过对原始DEM的填洼处理得到研究区的坡度坡向图层，用于模型坡度坡长因子(LS)的计算.

(3)土壤数据：来源于寒区旱区科学数据中心(http：//westdc.westgis.ac.cn/)世界土壤数据库(HWSD)中的中国土壤数据集(v1.1)，空间分辨率为 1km，采集时间为2009年，主要用于In VEST模型参数的计算及土壤保持服务的评价.预处理主要是基于土壤质地属性代码转换为栅格图层.

(4)气象数据：来源于温室数据共享平台(http：//data.sheshiyuanyi.com/ Weather Data/)包括玛多县及其周边的9个站点，主要为降水数据，用于In VEST模型参数的提取.参考中国气象数据网(http：//data.cma.gov.cn)中降水产品的空间插值方法，通过简单样条函数插值得到研究区2000年、2010年、2020年的降水量分布.

(5)遥感反演的植被指数：主要是MODIS 归一化植被指数(NDVI)，来源于地球数据中心(https：//earthdata.nasa.gov/)16天合成的250m空间分辨率数据产品，时间为2000年、2010年、2020年，用于研究区数据预处理及In VEST模型参数的计算.对原始的NDVI数据利用MODIS Reprojection Tool(MRT)软件进行预处理(包括拼接、投影转换、重采样等)，之后按照元数据剔除无效值并进行0～1的取值范围换算.

(6)基础地理数据：包括玛多县行政区划数据，来源于全国地理信息资源目录服务系统1∶100万全国基础地理数据库(http：//www.webmap.cn/main.do?method=index)，用于研究区范围的提取.

3 In VEST土壤保持模型

In VEST模型是由美国斯坦福大学、大自然保护协会(TNC)以及世界自然基金会(WWF)于2007年为生态系统服务功能的空间表达、动态分析和定量评估而共同研发的量化模型[20].唐尧[21]等通过对相关文献的整理，认为该模型在生物多样性、土壤保持、水环境、碳储量等方面的应用具有较好的适用性[22-26].因此，本研究直接选择In VEST模型的土壤保持模块[20]，具体公式如下：

USLE=R×K×LS×C×P;

(1)

RKLS=R×K×LS;

(2)

SD=RKLS-USLE,

(3)

式中：USLE是实际土壤流失量[t/(km2·a)]，RKLS是潜在土壤流失量[t/(km2·a)]，R是降水侵蚀力因子[(MJ·mm)/(km2·h·a)]，K是土壤可蚀性因子(t·km2·MJ·mm)，LS是坡度坡长因子(无量纲)，C是植被覆盖因子(无量纲)，P是水土保持措施因子(无量纲)，SD是土壤保持量[t/(km2·a)].各个因子计算过程如下：

(1)降水侵蚀因子R

R是横量降水影响土壤侵蚀能力大小的一个重要指标[27]，降水侵蚀因子的计算受研究区地理环境的影响较大，玛多县河流湖泊密布，所以采用适用于河海流域的降水侵蚀因子计算方法[28]，即由Wischmerier[29]提出的基于月尺度的降雨侵蚀因子的计算方法(式4).

,

(4)

式中：R降雨侵蚀力因子[(MJ·mm)/(hm2·h·a)]；Pi代表第i月降雨量(mm)；P代表年降雨量(mm)，基于上述方法，得到玛多县2000年、2010年、2020年降水侵蚀因子空间分布如图2所示.

图2 玛多县降水侵蚀因子空间分布

(2)土壤可蚀性因子K

K是土壤被降雨侵蚀，迁移运输敏感程度的重要指标，是影响土壤侵蚀的内在因素[28].其计算方法可分为直接测定法、诺谟图法以及公式法[30]，其中公式法较前面两者应用比较普遍，而Sharpley和Williams等[31]在EPIC模型中发展的K值估算方法经过多次修正应用比较广泛[4，6，32]，所以本文选择Sharpley等[31]提出的K值计算方法(式5、式6)，获得的土壤可蚀性因子的空间栅格分布如图3a所示.

(5)

,

(6)

式中：K是土壤可蚀性因子(t·km2·MJ·mm)，SAN是土壤砂粒含量(%)，SIL是土壤粉粒含量(%)，CLA是土壤粘粒含量(%)，C是土壤有机质含量(%)

(3)坡度坡长因子LS

LS主要表示地形因素对降雨径流侵蚀强度的影响.已有研究中多采用RKLS模型中的坡度因子算法[4，33]，但是，张芷温等[34]研究表明RKLS模型仅在坡度<10°的地区具有较好的适用性[35]，而对坡度≥20°的情形，刘斌涛等[36]提出的修正坡度因子算法适用性较好.玛多县大多属于坡度范围大于25°的山区，因此选择修正的坡度坡长因子算法(式7、式8)，具体得出的LS空间分布如图3b、图3c.

;

(7)

,

(8)

式中：L表示坡长因子(无量纲)，S表示坡度因子(无量纲)，λ表示水平坡长(m)，m代表坡长指数(无量纲)，θ表示地面坡度.

图3 a.土壤可蚀性因子；b.坡长因子；c.坡度因子

(4)植被覆盖与管理因子C

C被定义为在相同生态系统生物组分条件下，一定植被覆盖和管理措施下土地的土壤流失量和农耕空闲时的土地土壤流失量之比[37].魏健美等[38]、楚冰洋等[39]、王敏等[40]根据蔡崇法等[41]对C因子的计算方法(式9、式10)开展了相应研究，并在自己的研究领域进行应用实践，并确定了C因子的值介于[0，1]，其中，0表示不产生水土流失.鉴于此，本文也采用蔡崇法等[41]的研究方法，基于NDVI数据对C因子进行计算，结果如图4所示.

,

(9)

式中：C是植被覆盖与管理因子，f是植被覆盖度(%)，其计算公式为：

,

(10)

式中：NDVI是归一化植被覆盖指数，NDVImax和NDVImin是研究区NDVI的最大值与最小值.

图4 植被覆盖与管理因子C空间分布图(注：小于等于0.09的取值为水体)

(5)水土保持因子P

P指采取一定水土保持措施的土壤流失量与未采取水土保持措施的耕地土壤流失量之比，其值介于0～1之间.不同地域单元因为自然地理环境和人类活动强度的不同，P的取值不同.例如，邱春霞等[42]根据经验方法估算了黄土高原不同土地利用类型的P值，而陈琼等[43]同样在前人研究的经验基础上，给出了三江源地区不同土地利用类型的水土保持与措施因子P值，本研究区恰好在三江源范围，所以直接取三江源的P值(表1、图5).

表1 玛多县不同土地利用类型的P值

图5 土地利用及水土保持因子P空间分布图

4 结果与分析

4.1 降水侵蚀因子和植被覆盖与管理因子时空变化分析

图2是基于式4计算得到的玛多县三期降水侵蚀因子空间分布图.从图可以看出，研究区2000年降水侵蚀因子(图2a)大体上呈从西北到东南逐渐增大的趋势.2010年降水侵蚀因子(图2b)与2000年(图2a)相比，分布趋势具有明显差异，其在研究区北部形成了一个低值中心，由低值中心向东南部和西南部降水侵蚀因子逐步增大，且在西南部形成一个最高值区域.2020年降水侵蚀因子(图2c)的分布规律与2000年(图2a)较为接近，但西北的低值区域进一步缩小，中值区域覆盖研究区的东北部、东部、东南部和南部部分区域.图2的结果说明，玛多县2000—2010年降水侵蚀因子可能呈现增加的趋势，2010—2020年的降水侵蚀因子可能呈现减小的趋势，整体上2000—2020年的降水侵蚀因子表现为先增加后减小的态势.

图4是基于式9、式10计算得到的玛多县三期植被覆盖与管理因子空间分布图.从图可以看出，研究区2000年植被覆盖与管理因子(图4a)按纬度分布较为明显，即由北到南大致呈低值-中值-高值分布的情况，其中值大都集中在35°线附近，是低值和高值的过渡地带.低值地带覆盖面积较小，中值和高值地带覆盖面积较大，且在某一地带中夹杂着呈块状或带状聚集的其他值的地块.2010年植被覆盖与管理因子(图4b)与2000年(图4a)相比，中值地带和高值地带略微北移，且在35°线以北的东北部地区形成明显的带状高值区.2020年植被覆盖与管理因子(图4c)相比2010年(图4b)，其低值地带未发生明显变化，但位于35°线以南的中部中值地带部分转变为高值地带.图4的结果说明，玛多县2000—2010年、2010—2020年两个阶段的植被覆盖与管理因子均可能呈现增加的趋势.伴随植被的变化，土地覆盖也发了明显的变化，主要表现为湿地大面积的消失而裸地大面积增加(图5)，由此使得土壤保持因子的空间格局发生了一定的变化.综上，2000-2020年，玛多县植被/土地覆盖所呈现的土壤保持能力在减弱.

4.2 实际土壤流失模数时空变化分析

图6是基于式1计算得到的玛多县三期实际土壤流失量空间分布图.从图可以看出，研究区2000年实际土壤流失量(图6a)的高值大致集中在东北部的东西两侧区域和南部偏东的大部分区域，且块状聚集比较明显，中部、西部地区实际土壤流失量较小，且大部分区域在0～6t/(km2·a).2010年实际土壤流失量(图6b)与2000年(图6a)相比，明显存在整体的增加趋势，无论是2000年的高值分布区还是低值分布区，实际土壤流失量均有显著的增加.2020年实际土壤流失量(图6c)相比2010年(图6b)明显得到改善，但是相比2000年(图6a)还是有微弱的增加.图6的结果说明，基于降水和植被/土地覆盖的变化，玛多县2000—2010年实际土壤流失量可能呈现增加的趋势，2010—2020年的实际土壤流失量可能呈现减小的趋势.

图6 玛多县实际土壤流失空间分布

通过统计分析显示(表2)玛多县2000年实际土壤流失量为1.38×107(t/a)，均值为610.92t/(km2·a)，2010年实际土壤流失量为2.01×107(t/a)，均值为889.81t/(km2·a)，2020年实际土壤流失量为1.59×107(t/a)，均值为703.88t/(km2·a).2010年实际土壤流失量比2000年显著增加，2020年实际土壤流失量相比2010年则有所降低，但相比于2000年还是有所增加，与其空间变化的趋势(图6)一致，说明玛多县的实际土壤保持能力近20年可能处于波动上升的趋势.

表2 玛多县多年土壤侵蚀量和土壤保持量

4.3 潜在土壤流失模数时空变化分析

图7是基于式2计算得到的玛多县三期潜在土壤流失量空间分布图.从图可以看出，与实际土壤流失空间分布(图6)基本一致，2000年研究区潜在土壤流失量(图7a)高值主要集中在东北部的东西两侧区域和南部偏东的大部分区域，且呈块状聚集，中部以西的大部分地区潜在土壤流失在0.1～8t/(km2·a).2010年研究区潜在土壤流失量(图7b)高值分部区域大致与2000年接近，但2010年潜在土壤流失量呈现由中部向西部及西北端逐渐增加空间变化态势，其中中部及西部地区变化明显，西南端潜在土壤流失量相比2000年也明显增加.2020年研究区潜在土壤流失(图7c)空间分布与2010年(图7b)十分接近，但相比于2010年有略微的改善，但和2000年(图7a)相比仍然表现为明显增加的状态.图8的结果说明，玛多县2000—2010年潜在土壤流失量也可能呈现增加的趋势，2010—2020年的潜在土壤流失量同样可能呈现较弱的减小趋势，但整体上2000—2020年的潜在土壤流失量表现为增加的态势.

图7 玛多县潜在土壤流失空间分布

通过统计分析显示(表2)玛多县2000潜在土壤流失量为1.6×107(t/a)，均值为730.1t/(km2·a)；2010年潜在土壤流失量为2.51×107(t/a)，均值为1110.7t/(km2·a)；2020年潜在土壤流失量为2.01×107(t/a)，均值为889.4t/(km2·a).2010年潜在土壤流失量比2000年也显著增加，2020年潜在土壤流失量比2010年也同样有所降低，但相比于2000年还是明显增加，与其空间变化的趋势(图7)一致，说明玛多县的潜在土壤流失形势比较严峻.

4.4 土壤保持模数时空变化分析

图8是基于式3计算得到的玛多县三期土壤保持模数空间分布图.从图可以看出，玛多县2000年的土壤保持量(图8a)高值集中在东部及西南端的小区域内，其他区域的土壤保持量比较小，且大部分在0～2t/(km2·a).相比于2000年(图8a)，2010年土壤保持量(图8b)明显提高，中部偏北的倾斜区域土壤保持量没有变化，其他区域均明显增加.2020年土壤保持量(图8c)相比2010年土壤保持量(图8b)则又明显降低，但相比2000年(图9a)的空间分布还是呈现明显的增加.图8的结果表明，在玛多县实际土壤流失量(图6)和潜在土壤流失量(图7)均呈现先增后减的趋势下，土壤保持量也呈现先增后减的态势.

图8 玛多县土壤保持空间分布

通过统计分析显示(表2)玛多县地区2000年土壤保持量为2.63×106(t/a)，均值为116.43t/(km2·a)；2010年土壤保持量为4.996×106(t/a)，均值为221.17t/(km2·a)；2020年土壤保持量为4.2×106(t/a)，均值为185.93t/(km2·a).这种变化趋势与其空间变化态势(图8)也是一致的，且仅从土壤保持量的角度考虑，玛多县2010之前的土壤保持能力可能处于不断提升的态势，但近10年则可能呈现弱的波动变化或减弱的态势.

5 讨论

本研究以玛多县为研究区，基于In VEST模型的土壤保持模块，利用遥感产品数据、降水数据、土壤数据等，计算了2000年、2010年和2020年三期的潜在土壤流失量、实际土壤流失量以及土壤保持量.综合本研究的结果，无论是潜在土壤流失量、实际土壤流失量还是土壤保持量，均表现为先增加后减少的变化，考查In VEST模型土壤保持模块的因子，本研究中主要的变化因子为降水和植被/土地覆盖，降水侵蚀因子反映了降水冲刷对土壤流失的侵蚀效应，理论上与土壤保持能力成反比，本研究中降水侵蚀因子在2000—2020年表现为先增加后减小的态势变化(图2)，与潜在土壤流失量、实际土壤流失量的变化态势相符合；而植被/土地覆盖的变化明显是持续减弱了土壤保持能力，这一结果与土壤保持能力的变化态势并不一致，但是也解释了在降水“先增后减”的条件下，土壤保持能力“后减”的幅度相较“先增”的幅度弱的原因.因此，仅考虑降水和植被/土地覆盖变化的条件下，主导玛多县土壤保持能力的因素是降水，但植被/土地覆盖变化则提供了环境条件，而考虑植被/土地覆盖变化明显的是湿地大面积消失而裸地大面积的增加，推测可能的原因是全球变暖的背景导致青藏高原冻土消融，地表水分随之流失致使植被退化，从而导致土壤保持能力下降.但是降水和植被/土地覆盖变化具体如何耦合土壤保持能力及冻土消融导致的湿地大面积消失和裸地增加如何影响青藏高原土壤保持能力的机制还有待进一步探讨.

另外，本研究的结果也可能存在不确定性，主要是In VEST模型源于Kavocs[44]和Shaw等[45]对美国明尼苏达洲和内华达山脉区域的研究，其应用必须要求有完整的水文观测数据，且需要一些经验数据进行模型参数的调整，直至模型输出的绝对值在正确范围内[21].本研究降水的数据首先存在一定的缺失，且站点密度较低，空间分布仅参考中国气象数据网(http：//data.cma.cn)网格化降水数据产品生产的方法，进行简单样条插值获得，所以降水侵蚀因子的计算结果可能并不精准，获取精准的空间化降水数据可能是进一步研究的关键.更重要的是，In VEST模型虽然在中国的研究应用已经取得了丰富的成果[22-26]，但由于中国地理环境的复杂性，可能仍处于初级的探索应用阶段，尚未形成权威性覆盖全国的经验参数.所以，还需要进一步研究分析和优化此模型在中国不同地理环境应用下的参数以进行细致的模型校验.即使如此，对比刘敏超等[16]对三江源地区的土壤保持能力的评估结果，本文研究的结果仍然具有一定的合理性，可为相关研究和应用提供参考.

6 结论

(1)玛多县土壤保持能力的三个指标(潜在土壤流失量、实际土壤流失量、土壤保持量)2000—2020年均呈现先增加后减少的态势.其中，2000年、2010年、2020年的实际土壤流失量均值分别为610.92t/(km2·a)、889.81t/(km2·a)、703.88t/(km2·a)；潜在土壤流失量均值分别为730.1t/(km2·a)、1110.7t/(km2·a)、889.4t/(km2·a)；土壤保持量均值分别为116.43t/(km2·a)、221.17t/(km2·a)、185.93t/(km2·a).

(2)空间分布上，潜在土壤流失量、实际土壤流失量在三期均表现为呈点状或块状稀疏分布，高值主要集中在东北部的东西两侧区域和南部偏东的大部分区域，中部以西地区潜在土壤流失量、实际土壤流失量较小.土壤保持量则呈块状分布，高值集中在东部及西南端小部分区域，其他区域土壤保持量较小.

(3)玛多县2000—2020年的降水侵蚀因子表现为先增加后减小的态势，植被/土地覆盖所呈现的土壤保持能力则表现为减弱趋势.在仅考虑降水和植被/土地覆盖变化的条件下，玛多县土壤保持能力变化以降水为主导，但植被/土地覆盖的变化则为其提供环境条件，而植被环境变化背后的机制性因素可能来自于全球变暖导致的青藏高原冻土消融，但需要进一步探讨验证.