土地荒漠化作为一个重要的全球性生态环境及社会经济问题, 不仅造成生态环境的恶化和全球粮食产量的降低, 甚至会危害到人类生存和社会稳定^[1-3]。2015年《联合国防止荒漠化公约》统计显示, 全球约有30亿人口居住在土地荒漠化区；而且到2045年, 大约1.35亿人将会因荒漠化问题而流离失所。由于荒漠化和干旱, 全球农用地正以每年12万km²的速度转为荒地, 而通过修复生态系统退化区的土壤, 有望每年储存30亿t碳^[4]。中蒙俄经济走廊作为“一带一路”倡议中的六大经济走廊之一, 而中蒙俄跨境铁路、公路和原油管道作为支撑中蒙俄经济走廊区域合作发展的关键枢纽和通道, 其沿线区域的荒漠化问题直接关系到跨境铁路、公路和原油管道建设及安全运营重大生态环境^[5]。因此, 对中蒙俄经济走廊开展土地荒漠化程度识别与驱动机理研究, 不仅可以作为有效防止和治理荒漠化的基本前提, 而且有助于开展生物多样性、自然保护区及其他生态环境保护领域的技术交流合作。

自20世纪70年代以来, 国内外学者^[6-8]利用遥感数据覆盖范围广, 获取速度快等优势, 将不同平台或传感器的影像数据反演的植被指数^[9-11]与土壤^[12]、气候^[13]、地形^[14]等多种数据源结合, 对土地荒漠化程度识别。譬如, Lamchin等基于NDVI、表土粒径指数(TGSI)和地表反照率(Albedo), 并结合决策树方法识别了蒙古国的霍格诺汗自然保护区荒漠化^[15]。在此基础上, 大量学者尝试将遥感影像反演参数作为输入数据, 运用分类回归树(CART)^[16]、支持向量机(SVM)^[17]和随机森林(RF)^[18]等机器学习方法对荒漠化状况进行识别。其中, CART方法是Breiman在1984年提出的基于决策树构建的既可用于分类, 也可用于连续变量预测的算法^[19]；SVM模型是一种基于统计学习理论的重要机器学习算法, 即将非线性可分数据映射到高维空间, 完成线性可分, 已被广泛应用于遥感图像的分类^[20]；RF方法是Breiman在2001年提出的一种基于分类回归树的分类器, 能够处理较为复杂的遥感图像, 适用于多类别与多特征数据^[21]。另外, 还有学者在对地表反照率(Albedo)和归一化植被指数(NDVI)进行集成并构建Albedo-NDVI特征空间方法, 有效实现了荒漠化程度的等级划分^[22]。

虽然利用遥感数据和荒漠化程度评估模型可以识别荒漠化状况的动态, 但区分并测度荒漠化进程中气候变化和人类活动相对作用强度仍然是一个具有挑战性的工作^[23]。现有的基于实际与潜在的植被状况^[24-25]、残差趋势法^[26]、降水利用率指数^[27]虽然能揭示土地荒漠化过程中人类活动的潜在驱动作用, 但是无法定量评估气候变化和人类活动在荒漠化过程中的相对作用强度^[28]。基于净初级生产力的人类占用(Human Appropriation of NPP, HANPP)部分^[29], 即气候和遥感数据模拟的潜在净初级生产力(Potential NPP, PNPP)和实际存留在生态系统的实际NPP(Actual NPP, ANPP)之间的差异, 可分别用PNPP和HANPP表征气候变化和人类活动的影响^[30], 并利用线性趋势分析方法定量评估荒漠化过程中气候变化和人类活动所占比例。

论文旨在针对中蒙俄经济走廊区域的荒漠化格局及程度的定量分析问题, 利用Landsat、Sentinel-2等多源遥感影像数据构建荒漠化程度判别标准和指标体系, 基于Google Earth Engine(GEE)云平台, 在CART、SVM、RF和Albedo-NDVI四个模型对荒漠化动态识别精度对比的基础上, 筛选最优模型方法, 实现中蒙俄经济走廊荒漠化格局及程度的定量识别。并在此基础上, 引入线性趋势法计算PNPP和HANPP在2000—2015年的斜率, 实现该时段内PNPP和HANPP二者斜率的正负组合的情景分析, 进而实现气候变化和人文因素对中蒙俄经济走廊荒漠化驱动系数的定量测度。

1 数据与方法 1.1 数据收集与处理

用于荒漠化识别的多源遥感数据均来源于GEE平台。其中, 2000年荒漠化分类识别数据源为Landsat 5和Landsat 7数据, 2015年荒漠化分类识别数据源为Landsat 8数据。在荒漠化识别分类过程中, 首先, 利用GEE云计算平台对2000年和2015年的原始影像进行质量筛选和云掩膜分别获得3700景和2813景数据；其次, 采用NDVI最大值合成法进行镶嵌, 获得当前年份植被生长状况最好的地表反射率影像；再次, 由于荒漠化判别指标多为植被和土壤指数^[31], 为更好地对遥感影像进行荒漠化识别, 将水域、永久积雪和城镇建设用地统一归为“非植被区域”。另外, 为更好地构建机器学习算法中荒漠化识别分类的样本数据库, 结合Google earth的历史影像、Google照片及现有研究成果^{[17, 32-33]}, 采用欧空局的全球陆地覆盖数据和Sentinel-2 Level-1C产品作为荒漠化分级机器算法训练的样本数据的分类标准(表 1), 并基于随机算子选取70%的样本数据用于训练分类器, 30%样本数据用于验证。用于表征气候变化和人类活动的驱动数据, 主要包括MODIS的500m分辨率NPP产品(MOD17A3H和MOD17A2H)、分月份的WorldClim降雨数据、GLDAS-2.1气温数据、每8天的蒸散发数据和逐年份的GDP栅格数据^[34]。

1.2 荒漠化程度的判别指标

鉴于荒漠化程度判别指标应具备可反演性、可靠性和实用性的原则^[35], 选取植被、土壤、地形和气候等参数作为荒漠化程度的判别指标。根据提取的2000年和2015年的地表反射率影像数据, 采用各个指标求算模型(表 2), 分别获得：表征地表植被状态生物物理特征的归一化植被指数(NDVI)、植被覆盖度(FVC)和修正的土壤调节植被指数(MSAVI)；描述土壤信息的表土粒径指数(TGSI)、地表反照率(Albedo)、归一化水体指数(NDWI)和裸土指数(BSI)。另外, 地形指标和气候指标包括：高程、坡度、坡向和陆表温度(LST), 坡度和坡向是基于SRTM 30m数据, 经GEE平台的SRTM数字高程模型函数计算得出。

表 2中, Red为传感器红色波段, NIR为近红外波段, MIR为中红外波段, Blue为蓝色波段, Green为绿色波段, SWIR1为短波红外1波段, SWIR2为短波红外2波段；FVC中NDVI_S和NDVI_v分别表示裸土和植被的NDVI值, 为排除研究区NDVI的噪音现象, 选取区域内NDVI频率分布直方图的5%和95%的NDVI值分别作为NDVI_S和NDVI_v值。

在确定以上荒漠化程度提取指标及参数计算方法的基础上, 基于Google Earth Engine云计算平台上, 分别实现分类回归树(CART)、支持向量机(SVM)、随机森林(RF)和Albedo-NDVI(表 3)的模型算法, 提取中蒙俄经济走廊荒漠化程度的空间分布。

Albedo-NDVI理论公式：

式中, DDI为荒漠化差值指数。

1.3 荒漠化信息提取精度验证与评价

在分别运用CART、SVM、RF和Albedo-NDVI四个模型获取中蒙俄经济走廊荒漠化程度空间分布的基础上, 选取120个样点对荒漠化程度的分类和精度进行验证。其中包括研究区内24个全球通量观测塔站点和GEE平台随机生成的96个样点。在四种分类模型中, 基于GEE平台的随机算子抽取30%的样本点用于建立混淆矩阵, 得出总体精度(Overall Accuracy, OA)、生产者精度(Producer′s Accuracy, PA)、用户精度(User′s Accuracy, UA)及Kappa系数评价指标, 从不同的侧面对比荒漠化识别的精度, 其中生产者精度是一个样本点被正确分类的概率, 它是对丢失误差的测量。而用户精度是每一行中正确分类的样本点的个数与该行中所有样本点数的比值, 它表示了包含误差^[47-42]。其理论公式可表达如下:

式中, N代表样本点数的总和, X_ii代表混淆矩阵中第i种荒漠化类型反演结果与实际结果相同的样本点数, X_i+代表混淆矩阵中第i种荒漠化类型的列总和, X_+i代表混淆矩阵中第i种荒漠化类型的行总和。

通过建立混淆矩阵, 得出各个模型的识别精度估量指标(表 4)。结果表明, 在荒漠化分类结果中, CART模型是研究区荒漠化识别精度最高的分类方法, Kappa系数为0.754, 其次是Albedo-NDVI模型, 而SVM模型的精度最低。通过对生产者精度对比分析, Albedo-NDVI模型对中度、重度和极重度荒漠化具有较高的生产者精度, 说明对这三种荒漠化级别的漏测现象较小, 但是对轻度荒漠化的生产精度极低。CART模型对各个荒漠化级别的生产者精度都较高, 表明总体漏测现象较小。另外, CART模型对各个荒漠化级别的用户精度同样保持较高值, 且相比其他分类模型对非荒漠化区域的用户精度最高, 表明总体的多测误差较小, 而Albedo-NDVI模型对轻度荒漠化的用户精度极低, 说明在识别轻度荒漠化区域时, 具有很高的多测误差。SVM和RF的用户精度和生产者精度都较低, 表明两者在研究区的荒漠化识别效果并不好。

1.4 荒漠化驱动机理综合分析模型

基于MODIS数据和利用BIOME-BGC模型与光能利用率模型(CASA)模拟获得的NPP(ANPP)数据产品, 提取中蒙俄经济走廊区域的500m分辨率的NPP(ANPP)空间数据。在GEE云平台上, 实现Thornthwaite memorial模型^[43]的模型算法设计, 进而求解PNPP和HANPP。并结合降水和气温数据的空间数据, 在GEE平台上构建并实现了中蒙俄经济走廊荒漠化驱动机理综合分析模型, 其理论公式可表达如下：

式中, r表示年降水量(mm), t为年平均气温(℃), x表示年份, y为PNPP时得到Slope_PNPP(S_p), 为HANPP时得到Slope_HANPP(S_H), n为时间间隔。

对于研究区域, 首先基于给定时间段内荒漠化发展状况, 即荒漠化程度加重的区域定义为“退化区”, 荒漠化程度减轻的区域定义为“恢复区”, 将S_p和S_H的正负组合分为8种情景(表 5), 以测度气候变化和人类活动所占比例。在驱动因子定性方面, 当S_p为正值、或S_H为负值时, 表明在2000—2015年, 气候变化或者人类活动有利于荒漠化区域的恢复；而当S_p为负值、或S_H为正值时, 表明2000—2015年, 气候变化或者人类活动会加重土地荒漠化；在驱动因子定量方面, 如在荒漠化恢复区域, 情景1表示恢复完全归因于气候变化, 情景2表示恢复完全归因于人类活动, 情景3表示荒漠化恢复由气候变化和人类活动共同引起, 情景4表示这种状况超出模型的合理范围, 气候变化和人类活动都无法解释荒漠化恢复现象, 此时的“Error”区域是指被标记为错误的区域, 定义为气候变化和人类活动对植被的影响, 与研究区植被的实际生长状况并不符。即研究区气候变化和人类活动的驱动作用将加重土地荒漠化, 但实际上研究区的荒漠化现象被逐渐改善。

2 结果与分析 2.1 荒漠化时空分布格局及其转换趋势

根据荒漠化程度指标体系和最佳荒漠化程度评估模型获得2000年和2015年研究区荒漠化现状(图 1)。通过对比两期荒漠化现状, 发现研究区荒漠化空间分布具有很强的地带性特征, 极重度荒漠化区域主要分布在内蒙古自治区西部, 蒙古国的南戈壁盟、东戈壁盟和中戈壁盟, 重度和中度荒漠化区域主要分布在蒙古国的东戈壁盟、及其与内蒙古自治区交界处, 轻度荒漠化主要分布在蒙古国的东北部, 内蒙古自治区的南部, 而非荒漠化区域主要分布在俄罗斯区域和中国东北区域。另外, 在2000年到2015年间, 中蒙俄经济走廊区域的轻度荒漠化面积从71.97万km²增加到94.8万km², 以1.52万km²/a的速度扩张, 极重度、重度和中度荒漠化面积均呈现减小趋势, 其中极重度荒漠化面积减少最快, 以4873km²/a速度递减(表 6)。其中, 中国区域的荒漠化呈恢复趋势, 尤其是极重度、重度和中度荒漠化区域呈缩小趋势；蒙古国区域的荒漠化状况则呈加剧趋势, 尤其是极重度、中度和轻度荒漠化区域呈逐渐扩张趋势；而俄罗斯区域荒漠化的加重和减轻两个方向的状况均最为缓慢。

为更好地理解研究区2000—2015年荒漠化空间动态变化, 将荒漠化空间格局变化分为五类(图 2, 表 7), 分别为强烈退化(荒漠化加剧两级及其以上)、退化(荒漠化加剧到相邻一级)、稳定(荒漠化程度不变)、恢复(荒漠化降低到相邻一级)和显著恢复(荒漠化降低两级及其以上)。空间统计结果显示, 在2000—2015年间, 中蒙俄经济走廊的荒漠化程度保持稳定的区域主要分布在中国东部区域和俄罗斯区域, 退化和严重退化区主要分布在蒙古国的东部、中国区域与蒙古国区域的西部交界处, 恢复和显著恢复区域主要分布在内蒙古自治区和蒙古国的中部地区。另外, 在2000—2015年间, 中蒙俄经济走廊大部分区域的荒漠化程度未发生改变(占研究区面积的80.91%), 而荒漠化程度加重区域和恢复区域分别占总面积的9.68%和9.41%。其中, 在中国区域, 荒漠化恢复区域的面积比荒漠化程度加重区域的面积多15.18万km², 而在蒙古国区域和俄罗斯区域, 恢复区域的面积则分别比加重区域的总面积少13.20万km²和3.76万km²。

同时, 基于荒漠化转移矩阵(表 8)可以发现, 2000—2015年研究区荒漠化恢复区域主要来源于极重度荒漠化转为重度荒漠化和轻度荒漠化转为非荒漠化, 而荒漠化加重区域主要来源于重度荒漠化转为极重度荒漠化。

2.2 荒漠化驱动力定量评估

结合荒漠恢复区和退化区荒漠化程度数据, 运行荒漠化驱动机理综合分析模型, 获得气候变化和人类活动在中蒙俄经济走廊荒漠化“恢复区”和“退化区”的影响作用系数(图 3, 图 4, 表 9)。结果显示, 在2000—2015年间, 气候变化和人类活动在荒漠化恢复区的驱动作用在空间上具有明显的空间分异性。其中, 气候变化对蒙古国东部和内蒙古自治区的大部分区域荒漠化恢复起主导作用, 而且气候变化驱动作用>50%的荒漠化恢复区域占总恢复面积的68.8%, 而由人类活动驱动作用>50%的恢复区域仅占整个恢复区面积的27.36%。气候变化对荒漠恢复的驱动作用, 在中国区域、蒙古国区域和俄罗斯区域分别占荒漠恢复区面积的77%、58%和38%。另外, 在2000—2015年间, 在中蒙俄经济走廊的荒漠化程度加重区域, 人类活动的驱动效应则占主导作用, 其中, 人类活动驱动作用>50%的荒漠化加重区占所有荒漠化加重区域面积的比例高达69.68%, 主要分布在蒙古国的东部和内蒙古的中部地区。

2000—2015年间, 研究区荒漠化驱动机理结果中含有部分“Error”区域(图 3, 图 4), 其中恢复区的“Error”区域定义为气候变化和人类活动不利于荒漠化土地恢复, 但实际上有利于植被生长的区域, 这些“Error”区域主要分布在俄罗斯境内, 占所有恢复区域的3.83%。另外, 退化区的“Error”区域定义为气候变化和人类活动有利于荒漠化土地恢复, 但实际上抑制了植被生长的区域, 这些“Error”区域占所有退化区域的28.87%。事实上, 产生这种现象主要是因为遥感反演荒漠化结果存在一定的不确定性, 并且基于NPP计算驱动机理, 往往也受气候因子的噪音影响, 因此这些超出模型解释范畴的“Error”区域难以避免。

3 讨论

自21世纪以来, 中蒙俄经济走廊区域的极重度、重度和中度荒漠化面积减少, 轻度荒漠化面积增加；中国区域荒漠化土地恢复效果显著, 蒙古国区域荒漠化状况在加剧等结果与现有研究成果^{[15, 18, 44-46]}一致。譬如, 中蒙俄经济走廊区域的49.12%为草地, 裸地仅占14.59%。在气候变化和人类活动干扰的共同驱动下, 草地退化现象严重, 其中新增裸地面积的89.94%均自于草地退化, 但是裸地面积净减少近800km², 表明极重度荒漠化状况呈恢复趋势。

为了进一步解释线性趋势分析方法对气候变化和人类活动的荒漠化驱动作用测度的不足, 分别选取降水量、平均气温和年蒸散发量3个植被生长的关键气候参数^{[38, 40, 47]}和能够综合反映人类活动状态的人均GDP^[48], 对模型中未能解释的“Error”区域进行分析。空间统计分析显示, 在2000—2015年间, 中蒙俄经济走廊区域的年降水量减少29.7mm, 年平均气温增加1.57℃, 年蒸散发量减少17.3mm, 人均GDP增加6369＄。年降水量减少区域主要分布“非荒漠化”区(图 1), 尤其是俄罗斯区域的中部地区和中国区域的长白山脉部分地区。在中国区域和蒙古国区域的年降水量增加区域主要分布在荒漠化区, 其中中国区域的荒漠化“显著恢复”和“恢复”区与年降水量增加量大于50mm区的空间分布一致, 尤其是中国区域的东北平原、大兴安岭南侧地区和鄂尔多斯高原(图 1, 图 5)。

蒙古国区域的“极重度荒漠化”和“重度荒漠化”区的年降水量虽呈增加(0—50mm)趋势, 但其升温幅度(>2.8℃)和年蒸散发量增加幅度(>25mm), 尤其是荒漠化“严重退化”和“退化”区的升温幅度在3.0℃以上, 再加上该地区的人均GDP增加幅度(>8500＄)也大于其他区域。因此, 由于该区域升温和蒸散发增加在空间上的叠加, 使得荒漠化程度(尤其在“Error”区内)呈加重趋势。同时, 位于“Error”区的蒙古国南戈壁盟, 由于对该区域丰富矿产资源(如世界最大的铜矿奥尤陶勒盖金铜矿)大量开采^[49], 在致使人均GDP快速增长的同时, 也加重了土地沙化的程度。另外, 虽然中国区域的人均GDP也呈快速增长趋势, 但由于京津风沙源治理工程、退牧还草工程等系列国家重大生态工程的实施, 使得中国区域的荒漠化呈恢复趋势。

4 结论

基于GEE的云计算平台和多源数据, 结合植被、土壤、气候和地形等荒漠化程度的判别指标模型, 实现CART、SVM、RF和Albedo-NDVI 4个模型对荒漠化动态识别精度对比分析的基础上, 构建荒漠化驱动机理综合分析模型, 实现中蒙俄经济走廊荒漠化程度的时空变化及驱动效应分析。

精度验证表明, CART模型对中蒙俄经济走廊荒漠化程度信息提取和分类的精度高于SVM、RF和Albedo-NDVI模型, 其Kappa系数为0.754, 且对各个荒漠化程度分类的用户精度和生产者精度均比较稳定。虽然Albedo-NDVI模型的总体精度也比较高, 但是对于轻度荒漠化识别效果较差。

研究结果显示, 在2000—2015年间, 中蒙俄经济走廊区域的荒漠化程度整体呈减弱趋势, 尤其是极重度荒漠化、重度和中度荒漠化的区域呈缩减趋势, 但轻度荒漠化区域却呈增加趋势。中国区域的荒漠化呈恢复趋势, 净恢复面积达15.18万km², 而蒙古国区域和俄罗斯区域的荒漠化呈加重趋势, 其新增的荒漠化面积分别为13.2万km²和3.76万km²。

荒漠化程度时空变化及驱动作用分析表明, 中蒙俄经济走廊荒漠化程度减弱主要表现为极重度荒漠化转为重度荒漠化、轻度荒漠化转为非荒漠化两种状态, 而荒漠化程度加重区域则主要为重度荒漠化转为极重度荒漠化。气候变化在中蒙俄经济走廊荒漠化程度减弱和恢复种起到主导作用, 而局部的荒漠化加重区域, 人类活动的驱动作用超过了气候变化驱动作用。