中巴经济走廊2000-2017年逐月温度植被干旱指数数据集

冯克庭，张耀南，田德宇，康建芳

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

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

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

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

国家科技基础条件平台“特殊环境特殊功能观测研究台站共享服务平台”（Y719H71006）、中国科学院信息化专项“寒旱区环境演变研究‘科技领域云’的建设与应用”（XXH13506）。数据集组成 数据集由216个数据文件组成，数据文件为月温度植被干旱指数数据，数据格式为tif，文件名为TVDI.AYYYYDDD.1_km_month.tif。

引 言

中国-巴基斯坦经济走廊（简称“中巴经济走廊”），范围包括巴基斯坦全国和中国新疆喀什地区及周边（图1），起点在中国新疆喀什，终点在巴基斯坦瓜达尔港，全长3000 km，北接“丝绸之路经济带”、南连“21世纪海上丝绸之路”，是贯通南北丝路关键枢纽，是一条包括公路、铁路、油气和光缆通道在内的贸易走廊，也是“一带一路”倡议的重要组成部分[1-2]，将对中巴两国的经济社会发展产生重大的影响，为“一带一路”建设实施发挥示范和推动作用[3]。中巴经济走廊沿线干旱灾害发生频繁，严重影响着沿线国家的安全和社会发展，制约着“一带一路”重大战略的实施[4]。因此，有必要利用各类数据对中巴经济走廊沿途干旱灾害开展监测研究，这将对抗旱减灾及风险评价提供有力的理论依据，为进一步掌握有效的综合干旱指标提供科技支撑，为抗旱生产实践提供决策参考，促进中国与“一带一路”沿线各国的灾害监测、预警、救灾、减灾的科技合作。

图1 研究区示意图

干旱的监测和分析，长期以来都是政府和学术界高度关注的热点问题[5]。传统的干旱监测方法是基于地面台站观测或实验观测，利用气象和水文观测站获得的降水、气温、蒸发、径流等气象和水文数据，以及农业气象观测的墒情等数据，依据各种干旱指标对观测数据进行统计分析来对干旱情况进行量化分析。由于观测站点空间密度有限，仅靠地面观测点的资料很难对干旱进行大范围、快速、连续的监测。随着遥感技术的发展和应用，遥感干旱监测已经成为全球抗旱减灾中不可或缺的手段，它与传统学科相结合，优势互补，可以提供区域、大陆乃至全球的旱情信息[6]，是一种宏观、快速、客观、经济的有效手段[7]。在遥感干旱监测中，将植被指数和地表温度相结合进行干旱监测的方法使用广泛[8-13]。其中，温度植被干旱指数（Temperature Vegetation Dryness Index，TVDI）方法应用最为广泛。Sandholt等[14]基于植被指数和地表温度的关系，提出了TVDI指数用以估测土壤表层水分状况。国内学者利用TVDI指数进行了全国级、区域级、省级等不同空间尺度的干旱监测。齐述华等[15]应用TVDI方法对全国进行了干旱监测，结果表明，该方法用于大范围评价旱情是合理的。杨秀海等[16]研究表明TVDI基本上能够反映表层土壤湿度状况，利用TVDI对西北地区进行夏季干旱动态监测是可行的。姚春生等[17]利用 TVDI方法反演了新疆地区的土壤湿度。张顺谦等[18]依据TVDI对2006年四川伏旱进行监测与评估，其结果与气候监测结果基本一致。杜灵通等[19]利用TVDI监测宁夏地区干旱，并分析了10年间的旱情变化特征。沙莎等[20]利用历史遥感数据构建了3种植被指数与地表温度特征空间，讨论了TVDI方法在甘肃省陇东地区的适宜性。但是，TVDI用于干旱监测时，由于研究区内地形起伏、南北纬度跨距的差异对地表温度数据的影响会带来TVDI的计算误差，从而降低TVDI反演的精度。冉琼等[21]考虑高程变化对地表温度的影响，对地表温度进行了高程校正，使得经过高程校正获取的TVDI能更好地反映土壤湿度。赵杰鹏等[22]利用地理纬度和地面高程，对地表温度进行校正，从而达到对大区域太阳辐射和大气背景差异校正的目的，使得TVDI监测土壤湿度的精度明显提高。

本数据集采用TVDI方法[14]，并利用高程和纬度对地表温度校正[21-22]，计算温度植被干旱指数，以MODIS植被指数和地表温度数据为基础，结合数字高程模型（Digital Elevation Model，DEM）数据，获取了 2000-2017年“中巴经济走廊”区域逐月温度植被干旱指数数据集，为区域灾害研究与决策提供基础数据。

1 数据采集和处理方法

本数据集采用的数据为MODIS植被指数产品MOD13A3、地表温度产品MOD11A2，SRTM DEM产品以及气象站观测降水数据和土壤湿度数据。其中MODIS数据来源于美国航空航天局（National Aeronautics and Space Administration，NASA）陆地产品处理分发数据中心 （Land Processes Distribution Active Archive Center，LPDAAC，https://lpdaac.usgs.gov）；DEM 数据来源于地理空间数据云（http://www.gscloud.cn）提供的SRTM数据集；降水和土壤湿度数据来源于国家气象数据共享服务平台（https://data.cma.cn）。数据生产流程包括数据预处理、数据重建、计算与评估，如图2所示。

其中数据预处理是进行数据拼接、投影转换、波段提取、重采样等，提取的波段包括归一化植被指数（Normalized difference vegetation index，NDVI）、地表温度（Land Surface Temperature，LST）及相应的质量控制文件（Quality Assurance，QA）；数据重建包括空间插值、时间序列滤波、地形校正、时间序列补齐、时间序列重建；计算与评估对 TVDI进行计算，并利用标准化降水指数（Standardized Precipitation Index，SPI）和土壤湿度对指数进行评估。

1.1 数据预处理

MODIS植被指数产品MOD13A3时间分辨率为月，空间分辨率为1 km，地表温度产品MOD11A2时间分辨率为8天，空间分辨率为1km，数据格式为hdf，投影方式为正弦曲线地图投影。利用MODIS再投影工具（MODIS Reprojection Tool，MRT）对数据进行拼接、投影转换、波段提取和重采样。MRT参数设置：输出格式为Geotiff，输出投影采用地理投影，水准面选择WGS84，采样方法为最近邻采样，像元大小为0.0083333333度（1 km）。获取NDVI、地表温度及相关质量控制数据。

图2 数据生产流程图

DEM数据格式为Geotiff，分辨率为90 m。利用ArcGIS工具，对DEM数据进行拼接、重采样、裁剪，生成分辨率为1 km的中巴经济走廊DEM数据。

采用Python语言编程，结合研究区矢量边界和MODIS质量控制数据，实现数据批量裁剪和掩膜，剔除质量不可靠像元，生成中巴经济走廊质量可靠的NDVI和地表温度数据集。

1.2 数据重建处理

遥感数据重建旨在利用多种统计和数值分析方法，模拟缺失数据或提高反演模型精度，从而实现插补缺失观测值，优化时间序列数据，为相关研究提供更加完备的基础数据。遥感数据重建方法分为空间重建和时间重建两类，本数据集生产过程中对数据进行空间和时间重建处理。

1.2.1 空间插值处理

数据预处理提取的NDVI和温度数据，如有质量不可靠像元缺失，需要进行空间插补，采用反距离加权（IDW）方法实现数据的空间插值。利用Python首先调用ArcGIS arcpy包的RasterToPoint模块将栅格数据集转换为点要素，其次调用ArcGIS arcpy包的IDW模块，将点插值成栅格表面，实现数据批量空间插值处理，IDW方法参数设置如下：

（1）距离指数：用于控制内插值周围点的显著性，值越高，距离较远的数据点影像越小，通常取值范围在0.5～3之间可以获得最合理结果，本次处理选择值为2；

（2）搜索半径：定义对缺失像元值进行插值的输入点，包括可变搜索半径和固定距离两种方式指定输入采样点，选择可变搜索半径方式，插值的最邻近输入采样点数量为12；

（3）像元尺寸：设置为与输入影像数据相同。

1.2.2 时间序列滤波处理

在空间重建的基础上，采用S-G滤波，在时间序列上对多期影像进行拟合重建。S-G滤波拟合方法是由Savitzky[23]等在1964年提出的一种基于平滑时间序列数据和最小二乘原理的卷积算法，它是一种移动窗口的加权平均算法，但其加权系数不是简单的常数窗口，而是通过在滑动窗口内对给定高阶多项式的最小二乘拟合得出，其表达式为：

在S-G滤波中，拟合效果指数取最小值时的迭代结果为最佳滤波效果，其计算公式为：

式（2）中，Fk为第k次迭代后的序列拟合结果指数，Yi0和Yik分别为未经迭代和第k次迭代后序列中第i个值，Wi为序列中第i个值的权重，N为滤波器长度。

处理流程如下：

（1）确定滤波窗口的一半宽度m和多项式拟合的阶数d，通常，m值取4～7，d值取2～4，本数据集处理中选择m=4，d=2；

（2）对空间插值后的时序数据进行S-G滤波处理，生成时序数据；

（3）计算序列中每个点的权重；

（4）根据权重、空间插值的时序数据和S-G滤波后的序列，生成新的时序数据，并利用S-G滤波拟合新的时序数据；

（5）计算拟合效果指数，若拟合效果指数达到最小，退出迭代；否则，返回步骤（4）继续迭代处理。

1.2.3 温度数据地形校正处理

研究表明[14]，TVDI模型反演的精度受地表温度、植被覆盖状况、地表参数、大气条件及太阳辐射等因子影响，而地理纬度和地面高程是影响大气背景差异和太阳辐射的两个重要因素。研究区内地形起伏、南北纬度跨距的差异对MODIS地表温度数据的影响会带来TVDI的计算误差，因此，需要对地表温度进行校正[20-22]。地表温度校正的公式如式（3）。

式（3）中Tc为校正后地表温度，Ts为校正前MODIS地表温度，H为高程，L为纬度，a为高程校正系数，b、c分别为纬度校正系数，其中a常取0.006（0.6℃/100 m）[21-22]，也有研究表明新疆地区a的最优取值为 0.003～0.005（0.3～0.5℃/100 m），纬度校正系数b、c分别取 0.3～0.5、-20～-12[22]，本数据集生产中a取值为0.003，b取值为0.4，c取值为-16。

1.2.4 数据时间序列处理

由于MODIS数据序列中部分影像期数缺失，需要对缺失期影像进行补全。处理过程中，采取多年同期数据平均值的方法补全缺失期影像。由于计算中需要的NDVI和地表温度数据时间分辨率为月，而MODIS 1 km地表温度无月分辨率数据，因此将MOD11A2数据转换为月尺度，采用方法为月内8天数据取平均。

1.3 温度植被指数计算

TVDI[14]是一种基于NDVI-LST特征空间的土壤水分监测方法，具有一定的物理意义，较单独使用NDVI或LST能够提供更加准确、丰富的干旱信息。LST与NDVI的斜率与土壤水分的负相关关系是特征空间中的重要统计特征。随着地表植被覆盖度的增加，地表温度开始下降。当地表干旱缺水时，地表温度会迅速升高；反之，土壤湿度较大，地表温度增加较少。由特征空间原理可知，计算TVDI 的关键是干湿边的拟合。现有的研究结果表明，特征空间中干边上NDVI与LST都呈现显著的负相关关系，这说明当植被受到水分胁迫时，地表温度随着植被覆盖度的增加而降低。大多数研究中湿边上NDVI与LST呈现正相关或者不相关，以两者的正相关关系居多[24]。TVDI的计算公式为：

式（4）、（5）、（6）中Tc为任意像元的地表温度，Tcmax和Tcmin分别为一定NDVI对应的最低和最高地表温度，a1、b1、a2、b2分别为待定系数，TVDI取值范围为[0, 1]。对于一组NDVI和LST遥感影像，对NDVI取步长0.01，求取每一NDVI对应的LST最高、最低值，用最小二乘法拟合得到干、湿边方程，进而将式（5）、（6）代入式（4），求得TVDI。图3为研究区2009年逐月干湿边拟合效果，表1为干湿边的拟合结果。

图3 2009年逐月NDVI-LST特征空间及干湿边拟合图

从图3的NDVI-LST特征空间可以看出，随着植被指数的增加，地表温度的最大值减小，而地表温度的最小值增加。从拟合效果来看（表1），干边斜率都小于0，湿边斜率都大于0，干边拟合相关系数R2最小值为0.8889，其他11个月份都大于0.90，NDVI与LST呈显著的负相关关系（p<0.001）；湿边拟合相关系数R2最小值为0.5113，最大值为0.9283，NDVI与LST呈显著的正相关关系（p<0.001），干湿边拟合的整体效果很好。此外，干湿边拟合方程的常数项，在特征空间中是干湿边在纵轴（地表温度）的截距，代表了裸土像元在水分充足和水分匮乏时的地表温度，从表1中可以看出，干湿边截距随着年内温度的变化发生相应变化，即冬季截距较小，夏季截距较大。

表1 2009年逐月干湿边拟合方程和相关系数

2 数据样本描述

2.1 命名规则

中巴经济走廊 2000-2017年逐月温度植被干旱指数数据集命名规则如下：TVDI.AYYYYDDD.1_km_month.tif，具体意义为：

（1）TVDI：表示温度植被干旱指数产品；

（2）AYYYYDDD：表示产品时间为YYYY年第DDD天（以每年1月1日为第一天）；

（3）1_km：表示产品空间分辨率为1 km；

（4）month：表示产品为月数据。

如TVDI.A2017032.1_km_month.tif，表示2017年2月，空间分辨率为1 km的温度植被干旱指数产品。

2.2 数据描述

中巴经济走廊2000-2017年逐月温度植被干旱指数产品信息如表2。其中，TVDI的范围为0～1，在存储时放大10000倍，像元数值大小范围为0～10000，填充值为-3000，数据类型为signed int16，使用时需乘以比例因子0.0001。

表2 温度植被干旱指数产品信息

序号 内容 数值5 填充值 -3 0 0 0 6 行数 2 1 2 0 7 列数 2 2 7 7 8 像元大小 0.0 0 8 3 3 3 3 3 3 3, 0.0 0 8 3 3 3 3 3 3 3 9 坐标系 W G S 8 4

以TVDI作为干旱分级指标，将干旱等级划分为5级[15]：湿润（0≤TVDI≤0.2）、正常（0.2＜TVDI≤0.4）、轻旱（0.4＜TVDI≤0.6）、中旱（0.6＜TVDI≤0.8）和重旱（0.8＜TVDI≤1）。图 4 所示，以TVDI分级指标划分的干旱等级图。

图4 中巴经济走廊干旱等级图

3 数据质量控制和评估

3.1 数据质量控制

TVDI数据产品的质量与NDVI和LST产品的质量有关，因此，通过对NDVI和LST数据产品的质量控制来保证TVDI产品的质量。按照TVDI数据生产流程（图2），质量控制包括提取质量可信数据处理和对提取的质量可信数据产品进行时空重建等处理。产品中的不可信像元一般都为云覆盖或者大气效应影响较强的区域，提取质量可信数据处理依据MODIS产品质量控制文件，剔除了产品的质量不可信像元。对于不可信像元，采用其周边可信像元值进行空间插值，以补全缺失像元，进而对插值产品进行时序滤波重建，从时空两方面提高数据质量。

其中，提取质量可信数据处理，对NDVI数据，按照NDVI产品数据说明，依据QA值将NDVI像元分为可信和不可信两种，提取质量可信像元。处理流程如下：

（1）像元可靠性波段（1 km monthly pixel reliability）中像元值为0的像元标定为可信像元；

（2）研究区内有大量冰川分布，因此将1 km monthly pixel reliability中像元值为2的像元标定为可信像元；

（3）1 km monthly pixel reliability中像元值为3的像元标定为不可信像元；

（4）1 km monthly pixel reliability中像元值为1的像元标定为待定像元；

（5）对于待定像元，比较相应的1km monthly VI Quality像元值0-5位数值确定该像元的可信度：当0-1位值为0时，确定为可信像元；当0-1位值为1时，依据2-5位值确定可信度；

（6）掩膜处理，剔除不可信像元，提取可信像元，生成新的NDVI数据。

对LST数据，依据LST产品质量控制（QC）文件说明，提取质量可信像元，处理流程如下：

（1）质量控制波段QC_Day中0-1位值为0的像元标定为可信像元；

（2）QC_Day中0-1位值为2和3的像元标定为不可信像元；

（3）QC_Day中0-1位值为1的像元标定为待定像元；

（4）对于待定像元，比较QC_Day中相应像元2-7位的值确定该像元的可信度：当2-3位的值为0时，标定为可信像元；当2-3位的值为1时，对应的4-5位值为0，6-7位值为0的像元标定为可信像元；其他像元标定为不可信像元；

（5）掩膜处理，剔除不可信像元，提取可信像元，生成新的地表温度数据。

3.2 TVDI与站点SPI和土壤湿度的相关性分析

SPI已经被证实能够很好地反映气象干旱，TVDI也是一种土壤水分监测方法，因此，用SPI和实测土壤水分对TVDI进行评价。由于未获取到巴基斯坦境内的气象站观测资料，这里仅利用中国境内7个气象站降水数据计算出的月尺度SPI和3个气象站的土壤湿度数据，对TVDI进行有效性验证。计算1960-2016年7个气象站的SPI，提取2000-2016年7个气象站对应像元的TVDI值，并对各气象站TVDI与同期逐月SPI和土壤湿度进行相关性分析。结果表明（图5），TVDI与SPI、土壤湿度呈负相关关系，TVDI与月尺度SPI相关系数最大为-0.389，最小为-0.158，除乌恰站（p=0.024）TVDI与SPI的相关性通过了p＜0.05的显著性检验外, 其他站TVDI与SPI的相关性都通过了p＜0.01的显著性检验；TVDI与实测土壤湿度相关系数最大为-0.656，最小为-0.217，TVDI与土壤湿度的相关性都通过了p＜0.01的显著性检验。

图5 TVDI与各气象站SPI、土壤湿度相关图

4 数据使用方法和建议

中巴经济走廊2000-2017年逐月温度植被干旱指数数据集在月尺度上反映区域干旱的变化特征，数据按月存储，格式为GeoTiff格式。数据读取和操作可以用ArcGIS、ENVI等常用的GIS与遥感软件。TVDI以影像的像元值表示，在使用过程中，用户可以根据自己的分级标准对数据进行分级，生成干旱等级图。

致 谢

感谢USGS提供MODIS数据，感谢地理空间云提供DEM数据，感谢国家气象数据共享服务平台提供降水和土壤湿度数据。