近20 a 黄河源园区土地退化时空演化分析

秦 彤， 李功权， 范佳晨

（1.长江大学地球科学学院，湖北 武汉 430000；2.西北师范大学地理与环境科学学院，甘肃 兰州 730070）

土地退化被定义为土地生产能力的暂时或永久衰退［1］，主要分为土壤物质的迁移和土壤内部化学性质的变化［2］。由于其对环境、生物与农业等方面带来的严重破坏作用，得到国内外相关部门广泛重视。然而，对土地退化的评估与监测尚未达到一个共识［3］，特别是面向区域及以上空间尺度的研究［4］。

土地退化问题研究主要以遥感手段为主，通过遥感反演建立土地退化指数［3-7］的方法进行土地退化评价。土壤［4-7］、土地利用/土地覆被［5-6］、植被［3-7］、地表反照率［3,7］、气候［4］等因子均被前人应用于建立土地退化指数。近年来，为了得到更准确的评估、监测与预测结果，土地退化研究更加趋向于采用多源时空数据，应用机器/深度学习算法充分挖掘多源时空数据中隐含的信息。Yue 等［3］采用反向传播神经网络（Back propagation-artificial neural network，BPANN）得到了鄂尔多斯高原的土地退化监测结果；许明杰等［8］以一种加权自组织映射神经网络（Self or⁃ganizing map，SOM）对黄冈市土地生态进行了评估；Grinand等［9］采用随机森林算法进行了土地退化建模与预测。研究结果表明，机器/深度学习算法在土地退化分类与预测方面均取得了较好的效果［3,8-9］。但是，综合时空数据进行土地退化时空演化分析的研究仍有不足，结合地理空间人工智能（Geospatial ar⁃tificial intelligence，GeoAI）进行后续研究则更为匮乏。

三江源黄河源园区作为黄河的源头，对整个黄河流域至关重要，但由于水土流失与土地沙化等问题的存在导致生态发展受限［10-11］，国家对此也越来越重视。特别是近20 a以来，为了保护这片敏感的地区，《青海三江源自然保护区生态保护和建设总体规划》等政策相继出台。为了响应国家的号召，本文旨在对当地的土地退化治理提供理论依据。根据研究区现存的水土流失与土地沙漠化问题［10-11］，以2000、2005、2010、2015年和2020年5期的土壤侵蚀模数、沙化特征指数与土壤湿度指数为评价指标，应用空间托普利茨逆协方差聚类［12］（Spatial To⁃eplitz inverse covariance-based clustering，STICC）进行土地退化空间分布模式分析，并结合土地退化指数分析其时空演化过程。在此基础上，采用极致梯度提升（Extreme gradient boosting，XGBoost）算法［13］进行土地退化时空演化驱动机制探索。

1 研究区概况

三江源黄河源园区（97°01′20″～99°14′57″E，33°55′05″～35°28′15″N）位于青海省果洛藏族自治州玛多县境内［14］（图1），地处青藏高原，海拔4125～5257 m，中部偏南与东北地区以山地地貌为主，西部与北部地区以湖泊湿地为主，拥有扎陵湖、鄂陵湖2个最大的天然湖泊。植被主要分布在研究区东南、西南地区，北部相对稀疏。气温、降雨年度差异不大，大致呈现由南向北逐渐降低的趋势。年降雨量100～200 mm，集中在6—10 月，年平均气温维持在1°左右。但是随着全球气候变化、高原生态环境暖干化趋势加剧，再加上人类不合理活动的双重作用，三江源黄河源园区生态环境日益恶化、水土流失加剧、土地沙漠化面积扩大［10-11］。

图1 三江源黄河源园区地形图Fig.1 Topographic map of Sanjiangyuan Yellow River Source Park

2 数据与方法

2.1 数据来源

数字高程模型（Digital elevation model，DEM）数据：该数据来源于美国国家航天局地球科学数据系统（http://Earthdata.nasa.gov/）2019年发布的DEM V3数据，其空间分辨率为30 m，垂直精度为20 m。

降雨数据：降雨数据来源于ERA5数据集［15］，由于2020年的降雨数据不全，本文基于该数据集选取了2000、2005、2010、2015 年和2019 年共5 a 的逐日降雨量数据。

土壤数据：土壤数据来自于中国科学院资源环境科学数据中心（http://www.resdc.cn/）的1 km 分辨率的土壤质地类型数据，该数据集采集于全国第二次土壤普查期间；以及国家青藏高原科学数据中心（http://data.tpdc.ac.cn/）的中国土壤有机质含量数据［16］，其元数据更新于2021年。

植被覆盖度数据：植被覆盖度数据来源于MOD13A1的NDVI数据集，由于7—9月中旬植被旺盛，因此分别统计2000、2005、2010、2015 年和2020年5 a 的NDVI 均值。并在Google Earth Engine 平台上，通过像元二分法进行处理后得植被覆盖度数据。

遥感数据：利用Landsat系列卫星的蓝（BLUE）、短波红外1（SWIR1）以及短波红外2（SWIR2）提取沙化特征指数。为了避免植被与雨雪的干扰，在影像选取时以10—11月为时间区间，以云量低于15%为宜。由于2015年和2017年的Landsat 7影像存在损坏，2015 年的Landsat 5 和2020 年的Landsat 8 影像在该时间段内云量较高，最终选择了2000 年的Landsat 7、2004 年和2010 年的Landsat 5、2015 年和2019年的Landsat 8为基础数据源。参考郑子豪等［17］提出的土壤湿度指数提取方法，选取MOD09A1 卫星中的第1～7波段提取土壤湿度指数。由于北方每年7—9月降雨较充分、土壤湿度较高，且2000年的卫星影像存在异常，分别提取2001、2005、2010、2015年和2020年7—9月中旬土壤湿度指数均值作为研究数据。

2.2 研究方法

2.2.1 土地退化评价指标参考前人研究成果［5-6,10-11］，选用土壤侵蚀模数表征水土流失问题，以沙化特征指数与土壤湿度指数探究土地沙化问题。共线性检验得任何2 个因子之间的方差膨胀系数（Vari⁃ance inflation factor，VIF）均处于10 以下，所选因子之间不存在多重共线性问题。

（1）土壤侵蚀模数

自20世纪起，国内外研究广泛采用土壤侵蚀模数来反映水土流失问题。在区域尺度下，多采用通用土壤流失方程（Universal soil loss equation，USLE）进行水土流失监测［18］，具体计算公式如下：

A=R×K×LS×C×P（1）

式中：A为土壤侵蚀模数（t·km-2·a-1）；R为降雨侵蚀力因子（MJ·mm·hm-2·h-1·a-1）；K为土壤可蚀性因子（t·hm2·h·hm-2·MJ-1·mm-1）；LS 为坡长坡度因子；C为植被覆盖因子；P为水土保持工程措施因子。

在实际计算过程中，可根据不同地区的地理条件进行相应修正［19-26］。根据降雨数据，R值计算采用章文波等［20］提出的计算方法。根据张科利等［21］针对中国地区进行了修正后的K值计算方法，基于土壤数据进行土壤可蚀性计算。以DEM数据为准，在LS因子计算过程中，缓坡坡度因子采用McCool［22］的方法，陡坡坡度因子采用Liu等［22］提出的方法，结合Foster 等［23］提出的非流量累积坡长计算方法与Wischmier 等［24］提出的坡长因子计算公式提取坡长因子。基于植被覆盖度数据，C值采用蔡崇法等［25］针对中国地区提出的方法。P值采用章影等［26］使用的方法。

（2）沙化特征指数

由于已沙化地物的光谱特征与其他地物存在明显差别，可根据Landsat系列卫星的BLUE、SWIR1及SWIR2 波段反演得到土地沙化程度［5］。在Google Earth Engine 平台上对Landsat 系列数据集进行去云、拼接、裁剪以及波段运算后提取沙化特征指数，计算公式如下：

SFI=(SWIR1-BLUE)/(200-SWIR1-SWIR2)（2）式中：SFI 为沙化特征指数；BLUE、SWIR1 与SWIR2分别对应着Landsat的蓝、热红外1及热红外2波段。

（3）土壤湿度指数

土壤湿度指数既反映了土地干旱程度，又反映了土地沙化的风险［5-6］。MODIS 作为具有高时间分辨率的高光谱数据，在反演土壤湿度时起到了较好的 效 果［17］。在Google Earth Engine 平 台 上 对MOD09A1数据集采用缨帽变换计算土壤湿度指数，计算公式如下：

式中：WET 为土壤湿度指数；b1～b7 分别为MOD09A1卫星的波段1到波段7。

（4）评价指标归一化

由于所选土地退化影响因子存在量纲不一致性，因此需对其进行归一化。土壤侵蚀模数与沙化特征指数是正向指标，需进行正向归一化；土壤湿度指数是逆向指标，需进行逆向归一化。具体归一化公式如下：

式中：Xscale、Xscale-分别为单因子正向、逆向归一化后的结果；X为原始数据；Xmax为5期数据中的最大值；Xmin为5期数据中的最小值。

2.2.2 土地退化空间分布模式

（1）土地退化评价单元提取

流域作为地貌的最小单元，最合适作为土地退化的评价单元。为了提取到合适的流域，需要找到最佳汇水阈值［27］。并基于该阈值，以DEM 为数据源，通过水文分析进行流域提取。前人研究表明［27］，通过对不同汇水阈值下的河网密度进行均值变点分析，求取变点后，可得到最佳汇水阈值（图2）。分析可得，汇水阈值为300的点为变点，以此将整个黄河源园区划分为44816 个流域，并将此作为土地退化空间分布模式分析的基础。

图2 最佳汇水阈值提取Fig.2 Extraction of optimal catchment threshold

（2）土地退化空间分布模式

为了更优地显示土地退化空间分异现象，开展了SOM 神经网络、基于密度的聚类方法（Densitybased spatial clustering of applications with noise，DB⁃SCAN）和STICC 聚类算法对比分析。STICC 聚类结果的空间连续性强于基于属性的SOM神经网络，相似属性对象的聚集优于基于密度的DBSCAN。可见，STICC 聚类方法更适合进行地理现象空间分布问题的研究。该方法采用Delaunay 三角剖分法保持了空间上的连续性，使用多源高斯分布的逆协方差矩阵定义子区域属性组成的马尔科夫随机场，同时考虑空间连续性与重复地理模式发现问题［12］。

考虑到样本数据之间的不确定性与随机性存在差异，采用熵值法［28］对正向、逆向归一化后的土壤侵蚀模数、沙化特征指数与土壤湿度指数分别确权后（表1），以44816个流域作为聚类评价单元进行样本值的提取，最终得到一个3×44816的矩阵，并以此进行STICC聚类。为了避免异常值对聚类结果的干扰，首先采用局部Moran’sI指数［29］进行空间异常值的去除。聚类时，将分类个数统一设置为3，将逆协方差矩阵稀疏性水平调整指数（lambada）统一设置为10-3，将空间一致性调整指数（beta）统一设置为10-10，迭代次数分别设置为5、10、…、30。为了得到最佳迭代次数，采用分类适确性指数（Davies-Bould⁃in index，DBI）指数［8］对聚类结果进行评价，DBI指数越小表明聚类效果越好。结合DBI指数得到，2000、2005、2010、2015 年和2020 年的迭代次数分别为25、30、10、15和30时聚类结果较为适宜。

表1 土地退化指数权重确定Tab.1 Weight determination of land degradation indices

2.2.3 土地退化时空演化探索由于聚类结果仅仅能保证将相似的样本分为一类，无法得到该类所代表的土地退化程度高低，也无法在时空尺度上对土地退化程度进行比较。为了探索土地退化时空演化过程，还需要结合土地退化指数。该指数分布在0～1 之间，值越高，土地退化程度越高。整合5期土地退化指数数据，采用自然间断点分级法可将土地退化程度分为5级，土地退化指数处于0～0.098为微度土地退化、0.098～0.162 为轻度土地退化、0.162～0.217 为中度土地退化、0.217～0.290 为微重土地退化、＞0.290为重度土地退化。

LDI=a×USLEscale+b×SFIscale+c×WETscale-（6）式中：LDI 为土地退化指数；USLEscale、SFIscale为正向归一化土壤侵蚀模数、沙化特征指数；WETscale-为逆向归一化土壤湿度指数；a、b及c分别代表对应权重。

综合分析5 期STICC 聚类结果中3 类的类别内土地退化指数均值，可将3类进一步定性评价为高、中值及低值聚类区。在此基础上分别统计5 期STICC聚类结果中每一类内的土地退化指数分级结果面积占比，根据面积占比随时间演变的结果可进行土地退化时空演化过程探索。

2.2.4 土地退化驱动机制分析为了深入探究土地退化时空演化驱动机制，本文采用XGBoost算法［13］，以每一类中的正向、逆向归一化土壤侵蚀模数、沙化特征指数与土壤湿度指数与LDI分级结果为数据源，对不同类进行特征重要性评估。

3 结果与分析

3.1 土地退化空间分布模式分析及结果验证

STICC 聚类将整个研究区划分为3 类，第Ⅰ类分布面积最少，仅占10%左右，处于中部偏南与小部分东北地区；除2005 年外，第Ⅱ类占整个研究区面积大约50%，大体分布在偏北地区；第Ⅲ类除2005年面积占比较高外，其余年份占整个研究区面积比约为40%，以西南、东南地区分布最广（图3）。

第Ⅰ类所处区域土地退化程度相对较高，属于高值聚类区；第Ⅱ类为中值聚类区；第Ⅲ类表现为低值聚类区（图3）。结合2020 年Landsat 8 真彩色合成影像［3,5,7］与2020年的聚类结果进行对比后得到实际的土地退化情况与真实遥感反演的结果基本一致。示例3和5表现为土地退化程度轻或者不存在土地退化的一类，植被覆盖度高，整体呈现为绿色或者深绿色；示例1、2 和6 基本处于土地退化中值聚类区，植被覆盖度低、以裸土分布为主，地表呈现为较亮的颜色，地表分布较为平整，整体色差不大；示例4和7表现为土地退化程度最重的一类，地表破碎、沟壑纵横，土壤颜色偏红、偏黑（图4）。所以利用STICC聚类进行土地退化评价是可靠的。

图3 土地退化空间分布Fig.3 Spatial distribution of land degradation

图4 土地退化空间分布结果对比Fig.4 Comparison of spatial distribution results of land degradation

3.2 土地退化时空演化探索分析

基于空间分布模式分析结果，分别统计5 期STICC 聚类结果中第Ⅰ、Ⅱ类及第Ⅲ类内的土地退化指数分级结果（图3、表2）。2000年土地退化程度以轻度为主，占整个研究区面积的92.68%。第Ⅲ类占整个研究区面积的52.66%，主要分布在西南地区，在东南地区也有少量分布；第Ⅱ类分布在偏北地区，占据整个研究区面积的36.56%；第Ⅰ类主要分布在中部偏南地区与东北偏北地区，但分布面积仅占整个研究区面积的10.78%，且以轻、中度为主。2000—2005年，土地退化程度由以轻度为主转化为以微重为主，土地退化程度加深，第Ⅲ类大幅增长，第Ⅰ类与第Ⅱ类则存在一定程度的减少。这表明土地退化程度在偏西南地区增长趋势高于其他地区，第Ⅱ类中的中部偏北地区与西南地区之间的差异减小，中部偏北地区转化为第Ⅲ类，第Ⅲ类土地退化面积增长了33.29%。到2010 年时，土地退化程度由以微重为主好转至以中度为主，89.28%的区域处于中度土地退化。由于偏南地区的土地退化好转趋势相对较快，第Ⅲ类中的偏北与偏南地区之间的差异逐渐增大，偏北地区的第Ⅲ类向第Ⅰ、Ⅱ类转化，第Ⅱ类面积大幅增长，占据整个研究区面积的56.93%。同时，偏南地区的第Ⅰ类由于土地退化指数的降低，逐渐与偏北、偏东地区的一部分第Ⅱ、Ⅲ类融合，面积增长了13.22%，但第Ⅰ类土地退化程度好转至以中度为主。2010—2015年，土地退化指数继续下降，并转化为以微度为主，处于微度的区域占据整个研究区的94.12%。土地退化好转趋势由南向北扩展，由于中部略偏北土地退化现象好转，有18.79%的中部偏北地区的第Ⅰ、Ⅱ类转化为第Ⅲ类。2015—2020年，土地退化程度基本稳定，继续以微度为主，但第Ⅱ类面积增加，5.71%的第Ⅰ类与17.54%的第Ⅲ类转化为第Ⅱ类，这主要是由于中部偏西的第Ⅲ类土地退化指数存在小幅度的提高、中部偏东的第Ⅰ类土地退化指数存在小幅度的下降。同时，由于东北地区的土地退化指数存在小幅上涨，3.42%的区域土地退化程度加重至微度以上。

表2 土地退化程度时空演化统计Tab.2 Spatial and temporal evolution statistics of land degradation degree

3.3 土地退化驱动机制分析

本文采用XGBoost算法［13］，对5期数据的3个类别分别进行土地退化特征因子重要性分析（图5）发现，土壤侵蚀模数是造成土地退化程度空间分异的主要原因，但2005年的第Ⅲ类、2015年和2020年的第Ⅱ类土壤湿度指数对土地退化空间分异的影响力有所提高。结合单因子分布情况（图6）可知，2000—2005年土地退化程度加重，水土流失风险提高、土地沙化风险小幅降低。由于降雨侵蚀力由19.92～276.79 MJ·mm·hm-2·h-1·a-1增长至76.39～695.06 MJ·mm·hm-2·h-1·a-1，且由中部偏南地区向偏西南地区偏移，植被覆盖度仅在中部偏北地区有小幅增长。因此，土壤侵蚀模数在中部偏南、特别是偏西南地区的土壤侵蚀模数大幅增长，逐渐趋近于土地退化程度较高的偏北地区，第Ⅲ类分布面积大幅增加。同时，土壤湿度在中部偏东南地区有所增长，但是在中部偏北地区变化不大，土壤湿度指数对第Ⅲ类的土地退化程度空间分异的影响力大幅提高。偏东南地区的土地沙化风险降低，土地退化程度加重趋势相对较缓。2010年，降雨侵蚀力在中部偏西南地区的高值由695.06 MJ·mm·hm-2·h-1·a-1下降至268.59 MJ·mm·hm-2·h-1·a-1，偏西北地区的降雨侵蚀力相对提高，植被覆盖度在中部偏北地区轻微下降。因此，偏南地区的土地退化好转趋势高于北方，中部偏西南地区与偏北地区的空间分异拉大，分异成2 类，这也表明了2005 年中部偏西南地区的土壤侵蚀模数提高并不具备长期发展趋势。同时，偏北地区沙化特征指数下降、土壤湿度提高，土地沙化风险降低。2010—2015年，中部地区的植被覆盖度显著提高，土壤侵蚀模数继续保持下降趋势。但偏北地区的土壤湿度下降、沙化特征指数提高，土地沙化风险加大，因此第Ⅱ类分布向北偏移。第Ⅱ类中的偏东地区土地沙化现象更严重，造成了土壤湿度指数与沙化特征指数对第Ⅱ类的空间分异影响力提高。到2020年，降雨量在偏东地区有所增长，植被覆盖度继续提高，但土壤侵蚀模数整体变化不大。随着气温升高，地表蒸散量在偏北地区显著提高，偏北地区的土地沙化风险增加。东北偏北地区的土壤侵蚀模数还存在小幅上涨，造成了土地退化程度的轻微恶化。

图5 土地退化影响因子重要性Fig.5 Importance of land degradation impact factors

图6 土地退化影响因子分布Fig.6 Distributions of land degradation impact factors

4 讨论

黄河流域生态保护与高质量发展近年来逐渐被国家重视，黄河源头更应得到更多的关注与重视。研究结果指出，降雨侵蚀力因子分布在200～300 MJ·mm·hm-2·h-1·a-1，与章文波等［20］针对全国所作的降雨侵蚀力分布一致；土壤可蚀性因子分布在0.01～0.04 t·hm2·h·hm-2·MJ-1·mm-1之间，与张科利等［21］的研究结果一致；其余因子均来自官方发布的DEM 数据与遥感数据直接计算，并未进行二次处理。同时，胡光印等［30］对黄河流域沙漠化问题分析发现，黄河源园区沙漠化问题从偏东南地区向偏西北地区逐渐加重，与本文研究结果相符。

相比仅依靠土地退化指数［5-7］进行土地退化时空演化分析，本文从GeoAI 角度入手，引入了STICC聚类算法［12］。聚类后的空间分异性更显著，与许明杰等［8］研究一致。同时，本文分析得土地退化程度整体呈现波动性下降，与Ye等［7］针对黄河源园区的分析相符。但局部地区仍然存在问题，中部偏南地区与东北偏北地区的高水土流失问题，是黄河源园区土地退化空间分异的主要原因。在国家政策扶持下，水土流失问题逐渐缓解，土壤侵蚀模数对空间分异的影响力也不断降低。不过，中部偏北地区的土地沙化风险近10 a来有所增加，主要是由于地表蒸散量的增加，这与气候变暖息息相关［31］。特别是东北地区，水土流失风险与土地沙化风险均存在小幅提高，需要引起重视。

5 结论

本文根据三江源黄河源园区2000—2020 年存在的土地退化问题，综合应用土壤侵蚀模数、沙化特征指数与土壤湿度指数探究了该区域土地退化的时空演化过程，主要结论如下：

（1）第Ⅰ类主要分布在中部偏南与东北地区，属于高值聚类区，分布面积较少，以水土流失问题为主；第Ⅱ类主要分布在偏北地区，为中值聚类区，多年分布面积占半数左右，主要表现为土地沙化问题；西南与东南地区主要分布为土地退化程度较轻的第Ⅲ类。

（2）土地退化指数整体处于波动性下降，近20 a基本呈现轻度-微重-中度-微度-微度分布。土地退化程度在黄河源园区整体好转，在偏南地区的好转趋势更明显。

（3）黄河源园区土地退化空间分异现象主要受土壤侵蚀模数影响，但随着水土流失问题逐渐被治理，其对土地退化空间分异的影响力不断降低。同时，由于偏北地区近些年来渐渐显示出土地沙化的趋势，该区域内的沙化特征指数与土壤湿度指数对土地退化空间分异的影响力有轻微提高。