2000-2017年中巴经济走廊逐年荒漠化分布数据集

敏玉芳，冯克庭，康建芳，艾鸣浩

1. 中国科学院大学，北京 100049

2. 中国科学院西北生态环境资源研究院, 科学大数据中心，兰州 730000

3. 国家特殊环境、特殊功能观测研究台站共享服务平台, 平台服务中心，兰州730000

数据库（集）基本信息简介

引 言

荒漠化包括气候变化和人类活动在内的种种因素造成的干旱、半干旱和干旱亚湿润地区的土地退化，属于破坏性的生物地理过程，通常会造成生物多样性降低，土壤肥力下降乃至生态承载力丧失[1-2]。土地荒漠化是目前世界上最为严重的生态环境问题之一，据统计已有100多个国家正受到土地荒漠化的影响，荒漠化不仅破坏生态环境，也会削弱社会和经济的发展。监测土地荒漠化动态变化，掌握其变化规律对防治荒漠化有很重要的意义。

中巴经济走廊起点在中国新疆喀什，终点在巴基斯坦瓜达尔港，北接“丝绸之路经济带”、南连“21世纪海上丝绸之路”，是贯通南北丝路的关键枢纽，是一条包括公路、铁路、油气和光缆通道在内的贸易走廊，也是“一带一路”的重要组成部分。中巴经济走廊全长约3000公里，其中1300余公里为严寒高原、干旱和荒漠区段。尤其是在南段，干旱和大面积荒漠是其主要的生态环境约束因素。

遥感技术为人们提供了一种全新的荒漠化监测手段，它具有观测范围广、信息量大、数据更新快和精度高的特点，通过遥感图像解译或者定量反演都可以准确及时地获取荒漠化土地变化的信息。荒漠化在遥感影像上表现为裸地地表信息的增强和植被信息的减弱，可以采用地表反照率、地表温度、地表湿度、植被指数、植被覆盖度等指标因子表征。归一化植被指数（NDVI）是反映地表植被状态的重要生物物理参数，而地表反照率（Albedo）则是反映地表对太阳短波辐射反射特性的物理参量。随着荒漠化程度的加重，地表植被遭受严重破坏，地表植被盖度降低，生物量减少，地表粗糙度增加，在遥感图像上表现为NDVI值相应减少，地表反照率得到相应的增加。Li[3]等经过定位观测研究表明，当地表反照率达到一定数值时，会发生草地荒漠化，荒漠化发生的地表反照率阈值为0.3。Albedo 与 NDVI 这两个指标因子分别通过表征裸露地表信息、植被覆盖度信息，与荒漠化过程建立联系，具有较好的相关性。曾永年等[4]提出不单独依靠某一个参数，而是通过NDVI和Albedo的负相关关系构造“NDVI-Albedo”特征空间，并计算出荒漠化差值指数（DDI），利用DDI来进行荒漠化分级。该方法使用简单，易于获取指标，在荒漠化的分类及分级中，较仅使用遥感光谱信息进行分类的方法精度更高，所以在近年来不同地区的荒漠化定量评价中得到了一定的应用[4-7]。本文采用DDI评价中巴经济走廊荒漠化程度，研究区如图1所示。以NDVI与Albedo为监测指标，通过构造Albedo-NDVI特征空间，并利用Albedo和NDVI之间负相关的关系，构建中巴经济走廊DDI公式，并完成2000-2017年共18期中巴经济走廊荒漠化分类专题数据集，为相关部门的治理和决策提供技术支持。

1 数据采集和处理方法

1.1 数据来源与预处理

中巴经济走廊荒漠化分布数据的制备基于2种数据：MODIS植被指数产品MOD13A3数据和地表反照率产品MCD43A3数据，数据验证采用Landsat 7数据，详细列表如表1所示。

表1 研究采用数据列表

Albedo是遥感反演中的重要参数，指单位时间、单位面积内地表全部反射辐射通量与入射太阳总辐射通量之比。随着荒漠化程度增加，地表形态发生显著变化，集中表现在地表植被覆盖率下降，土壤有机质含量降低，土壤水分减少，地表粗糙度增加，反照率上升。MCD43A3是综合了Terra和Aqua卫星的数据产品，为天尺度L3级产品，时空分辨率是day/500 m。MCD43A3产品数据包含白空反照率和局地正午太阳角的黑空反照率，且白空反照率和黑空反照率都包含1-7波段共7个窄波段，以及3个宽波段可见光（0.3～0.7 μm）、近红外（0.7～5.0 μm）、短波（0.3～5.0 μm）。数据来源于美国国家航空航天局（NASA）陆地过程分布式数据档案中心（LPDAAC）。数据以HDF-EOS格式存储，采用正弦曲线投影。覆盖整个中巴走廊需要6景6个波段的数据。本研究获取2000-2017年18年MCD43A3数据共39 090景。首先使用MRT工具对数据进行镶嵌和重采样，将500 m分辨率的数据重采样为 1 km，并将所有数据正弦曲线投影转为 WGS84。然后利用中巴走廊矢量边界和Gdal的warp工具进行数据裁剪。

图1 研究区范围地理位置示意图

NDVI是反应植物量和植被检测的重要指标，主要应用于监测植被生长状态、植被覆盖度和消除部分辐射误差等。本文中使用1km的MODIS月合成植被指数产品MOD13A3。MOD13A3数据包括1 km月合成的NDVI、增强植被指数EVI、NDVI质量文件、EVI质量文件、可见光红光波段反射率、近红外波段反射率、可见光蓝光波段反射率、中红外波段反射率、平均观测天顶角、平均太阳天顶角、平均方位角。数据来源于LPDAAC。数据以HDF-EOS格式存储，采用Sinusoidal投影。本研究获取2000-2017年18年MOD13A3L3数据共1296景。首先使用MRT工具对数据进行镶嵌和重采样，将所有数据Sinusoidal坐标转为WGS84。然后利用中巴走廊矢量边界和Gdal的warp工具进行数据裁剪，最后得到216景研究区NDVI数据。

本文中还使用Landsat7 ETM+数据用于数据质量评估。Landsat 7 ETM+影像包括8个波段，band 1-band 5和band 7的空间分辨率为30 m，band 6的空间分辨率为60 m，band 8的空间分辨率为15 m，时间分辨率为16天。本研究中获取2010年7-8月，云量小于5%的Landsat 7影像8景。利用ENVI软件对数据进行几何校正和大气校正，NDVI计算与统计，最后利用像元二分模型进行植被覆盖度的遥感估算，得到研究区的植被覆盖度数据。

NDVI和Albedo原始数据有异常值和像元缺失的情况，质量不是非常可靠。需要进行无效值剔除和空间插补。NDVI数据的有效范围为-2000～10000，无效值用-3000填充。MCD43A3产品中所使用的1-6白空反照率波段的有效值范围为0～1000，无效值用32766填充。本文采用反距离加权（IDW）方法实现数据的空间差值，IDW方法参数设置如下：

（1）距离指数：通常取值范围在0.5～3之间，本文选择值为2；

（2）搜索半径：定义对缺失像元值进行插值的输入点，选择可变搜索半径方式，插值的最邻近输入采样点数量为12。

1.2 数据处理步骤

1.2.1 NDVI和Albedo数据年尺度重建

植被覆盖度的变化是荒漠化最直观的表现，需要在一年中植被最茂盛的情况下评价这一年的荒漠化程度。所以，基于年度NDVI和年度Albedo构建的荒漠化差值指数，需要以NDVI的年度最大值，Albedo的年度最小值作为基础数据。计算过程包括以下几个方面：

（1）Albedo数据计算

MODIS窄波段反照率向宽波段反照率的转化使用的是Liang[8-9]研究的算法，选用短波反照率计算方法，公式如下：

其中αshort是短波反照率，α1-α6分别代表MCD43A3中1-6波段，本研究中选用了白空反照率来计算Albedo，因为白空反照率是各个入射角的积分，更接近一般意义上的地表反照率。

（2）NDVI年度数据计算

采用最大值合成法（Maximum Value Compositing，MVC）。该方法在遥感数据处理方面，主要用于在一定时间段内像元的数据分析重构。计算公式：

其中，i为像元名，j为[1,n]时间区间的时间点，NDVIij指第i个像元在第j时间的NDVI值。

（3）Albedo年度数据计算

Albedo数据的年度重建采用最小值合成法。该方法是在某一指定的时间段内，选取像元的最小值作为新的像元值。年度Albedo的计算公式为：

其中，i为像元名，j为[1,n]时间区间的时间点，Albedoij指第i个像元在第j时间的Albedo值。

1.2.2 基于Albedo与NDVI构建荒漠化指数

NDVI值与植被覆盖度成正比，而Albedo与植被覆盖度成反比。根据NDVI与植被覆盖度的正相关关系，以及NDVI与Albedo的负相关关系可以得到Albedo-NDVI二维空间特征图。本文以2010年数据为例，先对NDVI和Albedo数据进行归一化处理，归一化处理公式如下：

然后在研究区域均匀选择1500个样本点，以NDVI作为X轴，Albedo作为Y轴，构建NDVIAlbedo特征空间。特征空间散点图如图2所示，根据线性统计回归分析结果，Albedo和NDVI间的负相关关系可以用以下线性关系式表示：

图2 中巴经济走廊2010年Albedo-NDVI二维特征空间图

如图2所示，由Albedo和NDVI构造的特征空间中，Albedo与NDVI的负相关关系式的不同位置，代表荒漠化不同阶段的状态和程度，荒漠化程度随着NDVI值的减少而增加，随着Albedo值的增加而增加，即Albedo与NDVI的负相关线性表达式可以反映荒漠化的变化趋势。根据Verstraete和 Pinty[10]的研究结论，如果在代表荒漠化变化趋势的垂直方向上划分 Albedo-NDVI特征空间，可以将不同的荒漠化土地有效区分开来，如图3所示，用荒漠化差值指数模型DDI来表示：

图3 NDVI-Albedo空间中荒漠化差值指数（DDI）图形表达

在具体应用中，常数α可根据Albedo-NDVI特征空间中斜率来确定，k为特征方程的斜率。

根据上面的计算结果，Albedo-NDVI的负相关关系斜率k=-0.2303，则α=4.3422，DDI的表达式如下：

1.2.3 荒漠化等级划分

为荒漠化评价与制图的需要，1984年联合国粮农组织（FAO）和联合国黄金规划署（UNEP）在《荒漠化评价与制图方案》中，从植被退化、风蚀、水蚀、盐碱化等4个方面，提出了荒漠化现状、发展速率、内在危险性评价的具体定量植被，并把荒漠化按其发展程度的不同，分为轻度、中度、重度、极重度。

通过上一节构建的 DDI模型，得到 2000-2017年的中巴经济走廊荒漠化指数专题数据。基于DDI进行荒漠化等级划分时，大部分研究人员都采用 ArcGIS重分类中的自然断裂法[5,7,11]。自然断裂法是基于统计学的Jank最优法得出的分级点[12]，这种分类法使各类别中差异最小，类之间差异最大。本文也采用自然断裂法和荒漠化程度四分法，将DDI值划分为5个区间，确定5级分级指标（表2），用户也可以结合实地调查数据，进一步精细化调整DDI分级表。

表2 DDI分级表

1.2.4 荒漠化指数专题数据模拟流程

首先重建NDVI和Albedo年数据，用最大值合成法计算年度NDVI值，用最小值合成法计算年度Albedo值；然后通过拟合Albedo和NDVI数据，构建Albedo-NDVI特征空间，得到拟合斜率k；最后构建荒漠化指数，进行荒漠化等级划分，完成基于荒漠化指数的荒漠化分级图。数据处理流程如图4所示。

图4 荒漠化指数专题数据模拟流程

2 数据样本描述

中巴经济走廊逐年荒漠化分布数据集有2000-2017年18个文件，每个文件夹包括年度NDVI、Albedo、DDI栅格数据和2010年荒漠化分级图示例。全部数据的空间分辨率均是1 km，年度NDVI、年度Albedo、DDI栅格数据保存格式为tif，荒漠化分级图采用jpg文件格式保存，用户可以根据DDI数据制作荒漠化分级数据。所有数据地理坐标系为WGS1984。数据结果展示如图5。

图5 2010年中巴经济走廊荒漠化分级图

3 数据质量控制和评估

本文采用高分辨率数据来评价荒漠化分级数据的质量。以2010年数据为例，选取8景Landsat 7 影像数据在小范围上进行验证。8景影像空间分辨率为30 m，平均云量小于5%，数据质量良好。DDI是通过NDVI年度最大值和Albedo年度最小值计算而来，所以本文用处于植被生长旺季的8月份Landsat数据作为验证数据，反演年度植被覆盖度最大值。运用ENVI软件对数据进行预处理，主要包括配准、几何校正、图像增强和FLAASH大气校正。首先计算8景影像的NDVI值，然后进行植被盖度的估算。植被覆盖度计算公式如下：

式中：vfc 为研究区的植被覆盖度，NDVIsoil为本研究区域纯沙漠的NDVI值，NDVIveg为完全植被覆盖区域的NDVI值。一般以5%置信度截取NDVI的上下阈值分别近似代表NDVIveg和NDVIsoil进行计算。

本文参考其他研究工作，将中巴经济走廊地区的荒漠化程度按植被覆盖度大概划分为：非荒漠化（高植被覆盖，vfc>60%）、轻度荒漠化（中等植被盖度，30%～60%）、中度荒漠化（中低植被覆盖，15%～30%）、重度荒漠化（低植被覆盖，<15%）。验证时，在每个影像上选取各级别荒漠化验证点30个，由于Landsat跟Modis数据分辨率不一致，所以在选取验证点时，在Landsat上选取小范围且植被覆盖度差别不大的区域作为一个验证点，与DDI分级结果进行比较（表3），并计算Kappa系数。表4是基于此方法验证的数据评价精度表，总体评价精度达到81.67%，Kappa系数为75.42%，两种分辨率的数据评价结果基本一致，反映出本文的数据在区域尺度上评价中巴经济走廊荒漠化具有可行性。

表4 单一类别精度及总体精度分析

4 数据价值

虽然荒漠化的影响因素复杂，使用目的多样化，造成指标因子繁杂，但所有的指标最终都可以归结到最基本的影响因子——植被。植被既是土地覆盖质量的综合体现，也是影响生态环境中水土保持等方面的主体因素。本文基于易获取的遥感影像数据，使用反映植被信息的NDVI与反映土壤、水分等信息的Albedo，构建中巴经济走廊的DDI模型，直观反映中巴经济走廊荒漠化程度。为定量评价荒漠化严重程度提供参考，为中巴经济走廊沿线国家进行荒漠化防治工作以及宏观决策的制定提供基础资料。

5 数据使用方法和建议

2000-2017年中巴经济走廊荒漠化分布相关数据保存为栅格tif格式。ArcGIS、QGIS、ENVI、ERDAS等常用的GIS与遥感软件均支持该数据的读取和操作。

致 谢

感谢LPDAAC提供的MODIS数据，感谢USGS提供的Landsat 7数据。感谢国家特殊环境、特殊功能观测研究台站共享服务平台给予的观测和数据方面指导与支持。