CN108760643A - 一种适宜于高海拔地区的干旱遥感监测方法 - Google Patents

Abstract

本发明公开一种适宜于高海拔地区的干旱遥感监测方法，包括以下步骤：步骤S10、收集目标区域的数据资料；步骤S20、将步骤S10中得到的数据进行预处理，得到增强型植被指数、地表温度指数、土地覆盖类型、降尺度降水数据；步骤S30、再通过步骤S20中的数据计算得到植被状态指数、温度条件指数、降水状态指数、重分类后的土地覆盖类型、高程；步骤S40、基于空间主成分分析构建干旱遥感监测模型。本发明综合考虑影响干旱的多方面因素，包括植被因子、地表温度因子、降水因子、土地覆盖类型因子和地形因子，利用空间主成分分析方法构建干旱监测模型可以有效消除所选变量之间相关性较大的变量，提取出较少的几项互不相关的综合性指标。

Description

一种适宜于高海拔地区的干旱遥感监测方法

技术领域

本发明属于环境遥感应用领域，尤其涉及一种适宜于高海拔地区的干旱遥感监测方法。

背景技术

干旱遥感监测是以干旱过程中的各致旱因子信息为基础，建立干旱遥感综合监测模型并进行干旱监测。干旱作为一种常见的，并且发生频率较高的灾害，一直以来给农业生产、粮食安全、生态环境及经济社会的发展造成了严重的威胁。中国是一个农业大国，农业的兴衰对国民经济有着非常严重的影响。而对于农业来说，旱灾比水患更为严重，每年遭受各类自然灾害的农田面积和粮食作物减产损失中，旱灾要占一半以上。因此，提高对干旱灾害的监测和评估，从而增强干旱灾害的防御和预警能力已经成为国民经济中亟待解决的重大科学问题。

现在的干旱监测多为气象干旱监测和遥感干旱监测，气象干旱监测方法运用较为成熟，传统的气象干旱监测是对气象站点获取的气象数据进行统计分析，但是每个区域的气象站点数量有限，空间分布不均且覆盖范围较小，气象数据不能及时获取，同时气象干旱监测只考虑了气温、降水等干旱致灾因素，很难全面获取干旱过程的所有信息，使干旱监测存在局限性。传统的以气象站点监测数据来监测干旱的方法已远不能满足现代生产和管理的需要，迫切要求采用遥感、GIS技术等现代技术手段对干旱灾害进行大面积、实时的动态监测，尤其是遥感技术可以快速高效的获取大范围连续空间面上的信息，弥补了传统气象站点干旱监测以点带面的不足，同时利用遥感技术获取的地表土壤和植被等水分载体对干旱响应的信息能更全面的反映干旱的真实情况。这些方法的应用有利于加强干旱灾害监测与评估的能力，可以切实提供能够支持抗灾决策和干旱灾害风险管理的理论和技术方法。

目前大多数的研究都集中在气象干旱、农业干旱或水文干旱某一单一的类型上，常规的干旱监测中，考虑最多的影响因素通常是降水量和温度，且往往采用气象站点观测数据，虽然较为准确，但是气象站点的分布不均匀，尤其是高海拔地区，难以获得较为准确的面状信息，而目前的大多数干旱遥感监测模型仅考虑了植被因素、地表温度某单一因素或温度-植被综合指数，但干旱的发生发展是复杂的过程，降水量的异常变化、植被的生长状况、蒸发量异常、土壤含水量都是衡量干旱的重要指标，单独依赖某一种指标很难衡量干旱这一自然过程。因此，在干旱监测时应综合考虑大气降水、蒸散、植被的生长状况以及土壤含水量等因素，只有研究这些因素的综合耦合过程，才能更准确地对干旱进行监测，本文结合气候环境、植被生长等多种因子利用多源遥感数据构建综合干旱监测模型有利于增强干旱监测模型的适用性。

目前利用多因素变量构建干旱监测模型的方法主要分为三种：(1)基于半经验的模型变量权重赋值法，该方法是主要对模型所选择的各变量进行3～4种经验性的权重赋值，通过构建的不同权重的干旱监测模型与土壤相对湿度和气象干旱指数的对比验证结果比较不同权重监测模型的干旱监测精度，从而选择适宜于研究区的最佳权重组合，构建最终的干旱监测模型；(2)基于数据挖掘的机器算法，当前应用较多的算法有分类回归树(CART)、人工神经网络(ANN)、支持向量机(SVM)、随机森林(RF)算法等，该方法主要是利用气象站点计算的干旱指数与站点对应的遥感数据进行模型的综合构建，以气象站点干旱指数作为因变量，各类对应的遥感变量作为自变量，通过机器算法进行训练构建模型，并预测其余区域的干旱指数，得出最终的干旱监测结果；(3)利用现有的商业模型构建软件建立干旱监测模型，该方法对前两种方法进行综合运用，也是通过气象站点监测数据计算的干旱指数作为因变量，各模型构建变量作为自变量，通过商业软件做趋势拟合构建多个不同监测模型，然后通过比较不同模型的监测精度选择适合研究区的干旱监测模型。但研究区位于高海拔地区，气象站点分布较少且极不均匀，无法表现研究区域内气象的整体情况，以极少的样本训练数据预测整个研究区范围内的整体情况可能会存在较大的差异且实现难度较大。

发明内容

本发明主要是克服现有技术中的不足之处，提出一种适宜于高海拔地区的干旱遥感监测方法，本方法采用空间主成分分析对研究区进行干旱综合监测模型的构建，该方法可以有效消除所选变量之间相关性较大的变量，提取出较少的几项互不相关的综合性指标，并且这些综合性指标能最大限度的反映出原来较多指标所反映的信息，同时利用空间主成分分析方法可以直接利用遥感数据进行计算，降低了以点带面所存在的点数据分布不均匀所带来的影响，同时该方法的权重为各主成分对应的贡献率，它反映了该主成分包含原始数据的信息量占全部信息量的比重，这样确定权重是客观的、合理的，克服了某些方法以人为经验主观赋权的缺陷。

本发明解决上述技术问题所提供的技术方案是：一种适宜于高海拔地区的干旱遥感监测方法，包括以下步骤：

步骤S10、收集目标区域的数据资料，所述数据资料包括MODIS遥感数据、 TRMM遥感数据和SRTM-DEM数据；

步骤S20、将步骤S10中得到的数据进行预处理，得到增强型植被指数、地表温度指数、土地覆盖类型、降尺度降水数据；

步骤S30、再通过步骤S20中的数据计算得到植被状态指数、温度条件指数、降水状态指数、土地利用类型、高程；

步骤S40、基于空间主成分分析构建干旱遥感监测模型。

进一步的是，所述MODIS遥感数据包括地表温度产品数据、植被指数产品数据和地表覆盖产品数据。

进一步的是，所述步骤S20中的具体处理过程为：

步骤S201、MODIS遥感数据预处理：

①利用MODIS处理工具对地表温度产品数据、植被指数产品数据和地表覆盖产品数据进行拼接、投影转换，统一为Albers投影坐标系，并通过MODIS 处理工具提取各产品单波段数据；

②然后对影像数据进行研究区的裁剪，同时使用质量控制文件剔除图像中的无效值，并使用多年均值充填算法对无效值区域进行填充修补；

③再对地表产品数据进行重采样统一空间分辨率为1km，利用最大值合成法得到地表温度月值；

④利用栅格计算器将像元值转化为摄氏度，转换关系如下：

New_Value＝Old_Value×0.02-273.15

式中：New_Value为转换后温度值，℃；Old_Value为像元原始值；

步骤S202、TRMM遥感数据预处理：

1)TRMM遥感数据投影变换及月值计算

对TRMM遥感数据进行投影变换，统一为Albers投影坐标系，0.25°分辨率，对影像数据进行研究区的裁剪，先对TRMM遥感数据进行月降水量的计算，计算公式如下：

P_month＝Value×24×Days_month

式中：P_month为该月总降水量，mm；Value为初始像元值，mm/h；Days_month为该月对应的天数；

3)再通过上述数据计算降水数据降尺度。

进一步的是，所述降水数据降尺度的具体计算步骤为：

步骤1、将1km空间分辨率的增强型植被指数重采样为与降水数据一致的空间分辨率，然后提取增强型植被指数与降水数据相同位置的点，建立数据库用于低分辨率GWR模型的训练数据；

步骤2、将低分辨率降水数据作为因变量，增强型植被指数作为自变量建立低分辨率GWR模型。在本研究中，选择了ADAPTIVE作为核函数和CV作为选择标准，因为具有这些参数建立的GWR模型可以模拟具有最高确定系数，此步骤中，得到每个位置下由增强型植被指数的值预测低分辨率的遥感预测降水，系数和残差；

步骤3、将步骤2得到的系数和残差按照spline方法插值为1km分辨率的空间数据；

步骤4、将1km分辨率的系数和残差代回步骤2的GWR模型中，其中以 1km空间分辨率的增强型植被指数为自变量得到1km的遥感预测降水；

步骤5、计算各气象台站观测降水与1km遥感预测降水的误差；

步骤6、运用克里金插值方法将步骤5的误差插值变为1km分辨率；

步骤7、将步骤4的1km遥感预测降水加上步骤6的误差表面，得到最终的降尺度降水数据。

进一步的是，所述步骤S30的具体计算步骤为：

步骤S301、利用增强型植被指数计算植被状态指数，其计算公式如下：

式中：VCI_i为某年i月的植被状态指数；EVI_i为i月的增强型植被指数值； EVI_max为多年相应月份增强型植被指数的最大值；EVI_min为多年相应月份增强型植被指数的最小值；

步骤S302、利用地表温度指数计算温度条件指数，其计算公式如下：

式中：TCI_i为某年i月的温度条件指数；LST_i为i月的温度条件指数值；LST_max为多年相应月份温度条件指数的最大值；LST_min为多年相应月份温度条件指数的最小值；

步骤S303、利用降尺度降水数据计算降水状态指数，其计算公式如下：

式中：PCI_i为某年i月的降水状态指数；TRMM_i为i月的降尺度降水数据值； TRMM_max为多年相应月份降尺度降水数据的最大值；TRMM_min为多年相应月份降尺度降水数据最小值；

步骤S304、根据不同土地覆盖类型对干旱的响应程度对土地覆盖类型按分等级赋值法进行重分类，得到重分类后的地表覆盖类型；

步骤S305、对目标区域SRTM-DEM数据和重分类后的地表覆盖类型数据进行极差标准化处理，统一量纲得到DEM数据，其计算公式如下：

式中：i为经度；j为纬度；x_ij′为(i,j)为处像元标准化后的数据；x_ij为(i,j) 为处像元原始值；x_max为同期影像像元最大值；x_min为同期影像像元最小值。

进一步的是，所述步骤S40中具体过程为：

步骤S401、利用空间主成分分析工具对植被状态指数、温度条件指数、降水状态指数、地表覆盖类型和DEM数据五个因子数据进行主成分分析，得出单月计算的主成分波段数据和对应的计算结果文本数据：

步骤S402、根据上述数据计算干旱综合监测指数，其计算公式如下：

式中：PC_i为第i个主成分；α_i为第i个主成分对应的贡献率的值；

步骤S403、最后根据干旱综合监测指数划分干旱等级，其中：

当干旱综合监测指数≤0.2时，干旱等级为极端干旱；

当0.2＜干旱综合监测指数≤0.3时，干旱等级为重度干旱；

当0.3＜干旱综合监测指数≤0.4时，干旱等级为中度干旱；

当0.4＜干旱综合监测指数≤0.5时，干旱等级为轻度干旱；

当0.5＜干旱综合监测指数≤0.6时，干旱等级为正常；

当干旱综合监测指数＞0.6时，干旱等级为湿润。

本发明的有益效果是：本发明综合考虑影响干旱的多方面因素，包括植被因子、地表温度因子、降水因子、土地覆盖类型因子和地形因子，同时利用空间主成分分析方法构建干旱监测模型可以有效消除所选变量之间相关性较大的变量，提取出较少的几项互不相关的综合性指标，且该方法直接利用遥感数据进行计算，降低了以点带面所存在的点数据分布不均匀所带来的影响。

附图说明

图1是本发明的线路图；

图2是干旱遥感监测模型构建流程图；

图3是研究区干旱监测结果图；

图4是研究区整体干旱等级划分结果图；

图5是研究区部分干旱等级划分结果图。

具体实施方式

下面结合实施例和附图对本发明做更进一步的说明。

如图1所示，本发明的一种适宜于高海拔地区的干旱遥感监测方法，包括以下步骤：

步骤S10、收集目标区域的数据资料，所述数据资料包括MODIS遥感数据、 TRMM遥感数据和SRTM-DEM数据，所述MODIS遥感数据包括地表温度产品数据、植被指数产品数据和地表覆盖产品数据；

步骤S20、将步骤S10中得到的数据进行预处理，对上述数据进行基本的投影转换统一坐标系，同时统一时间分辨率和空间分辨率并进行数据质量控制，其中对降水数据进行降尺度计算统一空间分辨率，得到增强型植被指数、地表温度指数、土地覆盖类型、降尺度降水数据；

步骤S30、再通过步骤S20中的数据计算得到植被状态指数、温度条件指数、降水状态指数、土地利用类型、高程；

步骤S40、基于空间主成分分析构建干旱遥感监测模型。

其中，所述步骤S20中的具体处理过程为：

步骤S201、MODIS遥感数据预处理：

①利用NASA官网提供的MODIS处理工具MRT(MODIS Reprojection Tool) 软件对地表温度产品数据、植被指数产品数据和地表覆盖产品数据进行拼接、投影转换，统一为Albers投影坐标系，并通过MRT工具提取各产品单波段数据；

②然后利用ArcGIS10.2对影像数据进行研究区的裁剪，同时使用质量控制文件剔除图像中的无效值，并使用多年均值充填算法对无效值区域进行填充修补；

③地表温度产品和植被指数产品空间分辨率为1km，地表覆盖产品的空间分辨率为500m，所以利用ArcGIS10.2对地表产品数据进行重采样统一空间分辨率为1km，因为白天的LST产品在生产时对大气，气溶胶散射和薄卷云进行大气校正，能够有效消除云噪声，所以本文选取白天LST产品作为数据源， MOD11A2为8天地表温度数据，因此利用最大值合成法得到地表温度月值；

④利用ArcGIS10.2栅格计算器将像元值转化为摄氏度，转换关系如下：

New_Value＝Old_Value×0.02-273.15

式中：New_Value为转换后温度值，℃；Old_Value为像元原始值；

步骤S202、TRMM遥感数据预处理：

1)TRMM遥感数据投影变换及月值计算

利用ENVI5.1对TRMM3B43数据进行投影变换，统一为Albers投影坐标系，0.25°分辨率，利用ArcGIS10.2对影像数据进行研究区的裁剪。同时， TRMM3B43数据中，数据初始像元值的大小对应的是每小时的降水量，单位为 mm/h，所以先对TRMM数据进行月降水量的计算，计算公式如下：

P_month＝Value×24×Days_month

式中：P_month为该月总降水量，mm；Value为初始像元值，mm/h；Days_month为该月对应的天数；

4)再通过上述数据计算降水数据降尺度，由于TRMM3B43空间分辨率为 0.25°，因此对该降水数据进行降尺度处理，统一为1km空间分辨率，本研究采用地理加权回归模型GWR进行降水数据的降尺度运算，地理加权回归基本思想是变量间的关系随着空间位置的变化而变化，通过估算每一位置的因变量与相关变量的参数来建立回归模型，GWR是常用线性回归模型的拓展。GWR使用空间x，y坐标直接构建位置和参数之间的关系以及因变量和自变量之间的局部拟合关系。作为描述空间非平稳性的工具，GWR中的回归系数不是基于全局信息；而是随着位置而变化，通过使用来自最近的邻近观测值的次采样数据的局部回归估计产生的。GWR的原理如下：

式中：(μ_i,ν_i)为i点坐标；β_j(μ_i,ν_i)——i观测点的第j个回归系数；ε_i——第 i个样点的随机误差。

很多研究表明植被指数NDVI和降水之间有着密切的联系，而且增强型植被指数EVI与降水的响应关系比NDVI与降水的响应关系更为强烈，尤其是降水量较大时其优势更为明显，所以本方法选取EVI指数参与降水数据的降尺度计算，降水数据降尺度具体方法步骤如下：

步骤1、将1km空间分辨率的增强型植被指数重采样为与降水数据一致的空间分辨率，然后提取增强型植被指数与降水数据相同位置的点，建立数据库用于低分辨率GWR模型的训练数据；

步骤2、使用ArcGIS10.2，将低分辨率降水数据作为因变量，增强型植被指数作为自变量建立低分辨率GWR模型。在本研究中，选择了ADAPTIVE作为核函数和CV作为选择标准，因为具有这些参数建立的GWR模型可以模拟具有最高确定系数，此步骤中，得到每个位置下由增强型植被指数的值预测低分辨率的遥感预测降水，系数和残差；

步骤3、将步骤2得到的系数和残差按照spline方法插值为1km分辨率的空间数据；

步骤4、将1km分辨率的系数和残差代回步骤2的GWR模型中，其中以 1km空间分辨率的增强型植被指数为自变量得到1km的遥感预测降水；

步骤5、计算各气象台站观测降水与1km遥感预测降水的误差；

步骤6、运用克里金插值方法将步骤5的误差插值变为1km分辨率；

步骤7、将步骤4的1km遥感预测降水加上步骤6的误差表面，得到最终的降尺度降水数据。

本发明优选的方案是，所述步骤S30的具体计算步骤为：

步骤S301、利用增强型植被指数EVI计算植被状态指数VCI，VCI是EVI 的基于像素的归一化，其中EVI的短期气候信号通过将其与长期生态信号分离而被滤除，VCI对干旱监测的假设是，植被活跃受气候变化的影响，在一定的气候地区和季节，极端干旱将削弱植被生长，并导致多年观测中EVI最低，相反，最高的EVI代表了最佳的气候条件，在植被干旱胁迫指标应用中VCI比EVI 更好，其计算公式如下：

式中：VCI_i为某年i月的植被状态指数；EVI_i为i月的增强型植被指数值； EVI_max为多年相应月份增强型植被指数的最大值；EVI_min为多年相应月份增强型植被指数的最小值；

步骤S302、用地表温度指数LST不能完全有效的监测干旱，所以利用LST 计算基于地表温度时间序列的温度条件指数TCI，该指数强调温度条件与植被长势的关系，可以弥补VCI受植被长势及生长季节的限制，在植被稀疏及作物播收期均可监测，TCI主要用于评估植被在较高温度条件下的生长状况以及干旱对水分的胁迫作用，被视为旱情监测的初始指标，其数值在0～1之间，是一个无量纲的量，其计算公式如下：

式中：TCI_i为某年i月的温度条件指数；LST_i为i月的温度条件指数值；LST_max为多年相应月份温度条件指数的最大值；LST_min为多年相应月份温度条件指数的最小值；

步骤S303、降水的多少代表干湿的程度，干旱年份的降水量与正常年份同一年的降水量会有较大差异，弱降水量表示不利或干旱条件，而较强的降水量则表现出湿润的状态，为有利条件；因此，提出的植被状态指数和温度条件指数，对降水数据进行相应的变换处理，得到降水状态指数，其计算公式如下：

式中：PCI_i为某年i月的降水状态指数；TRMM_i为i月的降尺度降水数据值； TRMM_max为多年相应月份降尺度降水数据的最大值；TRMM_min为多年相应月份降尺度降水数据最小值；

步骤S304、根据不同土地利用类型对干旱的响应程度对土地利用类型按分等级赋值法进行重分类，得到地表覆盖类型；

较荒地而言，植被能够改善土壤结构，使土壤的持水能力和水分调节能力增强，森林和林地通常对土壤水分具有较高的保持能力，而草地和农田由于只能利用表层土壤的水分，如果出现长时间的降水短缺或高温蒸散，便会导致土壤表层含水量下降、植被供水短缺，草地和农田会更容易发生干旱灾害，而农田受人为因素影响较大其水分的保持能力较弱于草地。因而本发明根据不同土地利用类型对于水分的保持能力从小到大和对干旱的响应程度对MCD12Q1数据的IGBP土地利用分类方案按照分等级赋值法将土地重分类为荒地、农耕地、草地、林地、水体和城镇建设用地6大类。具体分类规则见表2。

表2土地利用重分类

步骤S305、对目标区域SRTM-DEM数据和重分类后的地表覆盖类型数据进行极差标准化处理，统一量纲得到DEM数据，其计算公式如下：

式中：i为经度；j为纬度；x_ij′为(i,j)为处像元标准化后的数据；x_ij为(i,j) 为处像元原始值；x_max为同期影像像元最大值；x_min为同期影像像元最小值。

其中本发明优选的实施方式是，所述步骤S40中具体过程为：

步骤S401、利用ArcGIS10.2进行空间主成分分析工具对逐月植被状态指数 VCI、温度条件指数TCI、降水状态指数PCI、标准化后的当月对应的年土地利用类型数据LC和DEM数据五个因子数据进行主成分分析，得出单月计算的主成分波段数据和对应的计算结果文本数据：

其中空间主成分分析(Spatial Principal Component Analysis，SPCA)是指在地理信息系统软件平台支持下，将输入多元属性空间中的输入波段内的数据变换到相对于原始空间对轴进行旋转的新的多元属性空间，其轴相对于原始空间旋转，新空间中的轴(属性)互不相关，通过对原始空间轴的旋转，将相关的多变量空间数据转化为少数几个不相关的综合指标，将空间变量对因变量的影响程度分配到相应的主成分因子上，实现利用尽量少的综合指标最大限度地保留原来较多空间变量所反映的信息，完成主成分分析；

步骤S402、本发明提取累计贡献率大于85％的前n个主成分，对主成分进行加权求和计算研究区干旱综合监测指数，其中权重用每个主成分对应的贡献率来表示，即：

式中：PC_i为第i个主成分；α_i为第i个主成分对应的贡献率的值；

步骤S403、最后根据干旱综合监测指数划分干旱等级，其中：

当干旱综合监测指数≤0.2时，干旱等级为极端干旱；

当0.2＜干旱综合监测指数≤0.3时，干旱等级为重度干旱；

当0.3＜干旱综合监测指数≤0.4时，干旱等级为中度干旱；

当0.4＜干旱综合监测指数≤0.5时，干旱等级为轻度干旱；

当0.5＜干旱综合监测指数≤0.6时，干旱等级为正常；

当干旱综合监测指数＞0.6时，干旱等级为湿润。

实施例1

以处于高海拔地区的西藏东南部拉萨市、日喀则市、昌都市、山南市和林芝市为研究区为研究区，用本发明方法对该研究区进行干旱遥感综合监测。

研究区位于位于西藏自治区的藏南谷地和藏东峡谷地带两大区域，包括拉萨市、昌都市、林芝市、山南市和日喀则市五个地区56个县(区)，研究区主要分布于西藏自治区东部和南部，面积约51.69万km²，约占全区土地面积的 42.37％。该区域是西藏主要的农业区，拉萨市农业主要分布在拉萨河谷的羊八井、堆龙德庆等地，日喀则市耕种土壤面积全区最大，而山南市素有“西藏粮仓”之称，主要种植小麦、青棵、油菜等。根据1983～2013年西藏自治区气象灾害事件普查数据统计，研究区五个市在1983～2013年共计发生干旱灾害308次，约占全区总数的86.52％。因此对研究区域进行干旱监测的研究对于农业生产具有重要意义。

为了对研究区进行干旱遥感监测，采用本发明方法对研究区进行方法验证，步骤如下：

第一步，研究区遥感数据和基础地理数据收集及预处理。

(1)对区域进行干旱综合监测模型的构建，数据的获取及处理的准确性直接影响模型的监测精度以及干旱分布规律的统计分析。收集整理干旱遥感监测模型所需研究区2001～2015年MODIS和TRMM遥感数据、SRTM-DEM数据和基础地理信息数据。本方法所用到的基础数据如表1所示。

表1干旱综合监测模型构建所用基础数据表

考虑本研究研究区域较大且研究时间为2001～2015年共计180个月，时序较长，所以选用时间分辨率能满足要求的MODIS产品数据，本文根据研究区域大小将所有遥感数据空间分辨率都统一为1km。

(2)MODIS遥感数据预处理：

(3)TRMM遥感数据预处理：

1)TRMM遥感数据投影变换及月值计算

2)TRMM月降水数据降尺度计算

第二步，干旱遥感监测模型变量的计算。

(1)利用EVI指数、LST指数和TRMM降尺度结果计算研究区域2001～2015 年逐月植被状态指数VCI、温度条件指数TCI和降水状态指数PCI。

(2)对土地利用数据进行重分类。

(3)统一量纲。研究区地处青藏高原南部，青藏高原由于其特殊的地形，高海拔雪线以上地区常年冰雪覆盖，并不会发生干旱灾害，而海拔较低的河谷地带则大多处于干旱半干旱地区。对于地形因素而言，研究区域内海拔较高的区域发生干旱灾害几率较小，而高程相对较低的河谷地区则是干旱灾害的高发区域，所以将高程也纳入干旱影响因子。土地利用类型和高程数据与植被状态指数、温度条件指数和降水状态指数具有不同的量纲，进行综合计算时对计算结果会有较大的影响，因此对两者采用极差标准化处理转化为和前三个变量统一量纲的数据。

第三步，干旱遥感监测模型的构建。

本方法利用ArcGIS10.2进行空间主成分分析。工具生成的是波段数与指定的成分数相同的多波段数据。第一个主成分将具有最大的方差，第二个主成分将具有未通过第一个主成分描述的第二大方差，以此类推。工具输出的文本文件包含了输入的各空间数据间的协方差矩阵、相关系数矩阵、特征值、特征向量以及特征值对应的贡献率和累计贡献率。

本方法将2001～2015年每个月共计180期对应的植被状态指数VCI、温度条件指数TCI、降水状态指数PCI、标准化后的当月对应的年土地利用类型数据 LC和DEM数据五个因子数据进行主成分分析，得出单月计算的主成分波段数据和对应的计算结果文本数据，本研究提取累计贡献率大于85％的前n个主成分，对主成分进行加权求和计算研究区干旱综合监测指数，其中权重用每个主成分对应的贡献率来表示，即：

式中：PC_i为第i个主成分；α_i为第i个主成分对应的贡献率的值。具体流程图如图2。干旱遥感监测指数值越小，表示越干旱。部分(2015年1月～12月) 计算结果特征值、贡献率和累计贡献率如表3。其余时段干旱遥感监测指数按上述流程进行计算。部分计算结果遥感图见图3。

表32015年1月～12月模型计算特征值、贡献率和累计贡献率

第四步，干旱遥感监测计算结果干旱等级分类。

通过参考植被健康指数VHI干旱等级划分以及前人研究得出的干旱分类方案，对本研究所计算的干旱综合监测指数进行干旱等级划分，分类结果见图4 和图5。

第五步，结果验证。

(1)据中国天气网：2009年4月至6月中旬，西藏东部降水量为近10年来同期最少，特别是6月上中旬降水量为1956年以来历史同期最少值，出现了严重初夏旱。如果仅针对6月上中旬，也是比较罕见的极端天气，但是，如果比较4月至6月，极端性就不是那么明显。截至6月18日，西藏全区有五个地区27个县(市)部分地区不同程度受旱，受旱总面积近37万亩。其中，拉萨、日喀则地区中部、山南市、林芝市、昌都市大部有中等强度以上的气象干旱，在历史上属10～20年一遇的干旱。拉萨市和山南市的部分地区有重度气象干旱，在历史上属30年一遇的干旱。历来多雨丰水的“藏上江南”林芝地区入夏来降雨减半，边境察隅17日最高气温近32℃，超过历史极值，旱情有进一步扩大趋势。据统计，改则、狮泉河、尼木、南木林、泽当、浪卡子、日喀则和拉萨9个气象站已连续226天以上无有效降水，其中改则站监测结果显示，这一地区已有 299天无有效降水。

由图4可以看出2009年5月和6月西藏自治区东南五市都发生大面积干旱，图4a中在2009年5月昌都市大面积发生干旱现象，且大部分区域干旱等级在中等强度及以上；拉萨市大部分区域也发生中等轻度及以上的干旱现象；日喀则市西部边缘区域和东北部区域发生干旱显现；在山南市中部和北部与拉萨市相邻的区县都发生较为严重的干旱现象；林芝市发生干旱的区域主要集中于东部与昌都市相邻的区域，察隅县部分地区确有中等强度现象出现。图4b中显示 2009年6月西藏东南部区域发生的干旱强度有所减弱，大部分区域都发生轻度干旱，只有日喀则市中部和西部、拉萨市中南部和昌都市小部分区域有中等强度干旱，东南部察隅县干旱也大部分为轻度干旱。所以2009年5月和6月干旱综合监测指数计算结果与实际干旱事件发生区域范围基本吻合，且计算结果所显示的干旱发生区域也较为广泛。

(2)据西藏自治区气象灾害事件普查数据，2005年自入夏以来，拉萨市各县(区)持续出现高温少雨天气，致使各县(区)出现不同的旱情，同期5月下旬至6月中旬，日喀则地区出现高温少雨天气；5月下旬至6月份山南市贡嘎县降水异常偏少，尤其是6月份降水量只有18.0mm，出现持续性高温少雨天气，导致全县受到不同程度的旱灾；2005年入春以来山南市浪卡子县、昌都市左贡县出现高温少雨天气，降水量明显偏少，6月1日～7月6日浪卡子县降水量不足20mm，6月份的降水与历史同期值相比偏少67％，该时间段内左贡县降水总量不足40mm，干旱使得全县农牧业受到不同程度的旱灾。

通过对2005年6月干旱综合监测指数进行等级划分，并对上述干旱灾害事件的发生区域进行分析，结果如图5。由图5可以看出2005年6月，西藏自治区昌都市八宿县、察雅县、贡觉县、左贡县和芒康县部分地区发生干旱灾害，大部分区域为轻旱，只有极少数区域发生中等以上干旱灾害；图5中拉萨市尼木县、曲水县、堆龙德庆区、城关区和墨竹工卡县均有轻度干旱出现；而图5 中日喀则市则大部分区域发生干旱现象，其中仲巴县、萨嘎县、吉隆县、聂拉木县、定日县和康马县等地区还出现中等强度及以上的干旱情况；图5中山南市干旱主要发生在贡嘎县、浪卡子县、措美县、洛扎县、隆子县西部和错那县西北部，其中措美县和浪卡子县干旱发生区域较大，部分地区出现中等强度及以上干旱。通过对比2005年6月的实际干旱灾害事件与干旱综合监测结果可以看出干旱综合监测模型计算结果显示的干旱发生区域基本能够与实际干旱事件相对应，同时本研究所构建的干旱综合监测模型计算结果所划分的干旱发生区域较大于实际资料所统计的结果。

以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (6)

1.一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，包括以下步骤：

步骤S10、收集目标区域的数据资料，所述数据资料包括MODIS遥感数据、TRMM遥感数据和SRTM-DEM数据；

步骤S20、将步骤S10中得到的数据进行预处理，得到增强型植被指数、地表温度指数、土地覆盖类型、降尺度降水数据；

步骤S30、再通过步骤S20中的数据计算得到植被状态指数、温度条件指数、降水状态指数、土地覆盖类型、高程；

步骤S40、基于空间主成分分析构建干旱遥感监测模型。

2.根据权利要求1所述的一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，所述MODIS遥感数据包括地表温度产品数据、植被指数产品数据和地表覆盖产品数据。

3.根据权利要求2所述的一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，所述步骤S20中的具体处理过程为：

步骤S201、MODIS遥感数据预处理：

①利用MODIS处理工具对地表温度产品数据、植被指数产品数据和地表覆盖产品数据进行拼接、投影转换，统一为Albers投影坐标系，并通过MODIS处理工具提取各产品单波段数据；

②然后对影像数据进行研究区的裁剪，同时使用质量控制文件剔除图像中的无效值，并使用多年均值充填算法对无效值区域进行填充修补；

③再对地表产品数据进行重采样统一空间分辨率为1km，利用最大值合成法得到地表温度月值；

④利用ArcGIS栅格计算器将地表温度数据的像元值转化为摄氏度，转换关系如下：

New_Value＝Old_Value×0.02-273.15

式中：New_Value为转换后温度值，℃；Old_Value为像元原始值；

步骤S202、TRMM遥感数据预处理：

1)TRMM遥感数据投影变换及月值计算

对TRMM遥感数据进行投影变换，统一为Albers投影坐标系，0.25°分辨率，对影像数据进行研究区的裁剪，先对TRMM遥感数据进行月降水量的计算，计算公式如下：

P_month＝Value×24×Days_month

式中：P_month为该月总降水量，mm；Value为初始像元值，mm/h；Days_month为该月对应的天数；

2)再通过上述数据计算降水数据降尺度。

4.根据权利要求3所述的一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，所述降水数据降尺度的具体计算步骤为：

步骤1、将1km空间分辨率的增强型植被指数重采样为与降水数据一致的空间分辨率，然后提取增强型植被指数与降水数据相同位置的点，建立数据库用于低分辨率GWR模型的训练数据；

步骤2、将低分辨率降水数据作为因变量，增强型植被指数作为自变量建立低分辨率GWR模型。在本研究中，选择了ADAPTIVE作为核函数和CV作为选择标准，因为具有这些参数建立的GWR模型可以模拟具有最高确定系数，此步骤中，得到每个位置下由增强型植被指数的值预测低分辨率的遥感预测降水，系数和残差；

步骤3、将步骤2得到的系数和残差按照spline方法插值为1km分辨率的空间数据；

步骤4、将1km分辨率的系数和残差代回步骤2的GWR模型中，其中以1km空间分辨率的增强型植被指数为自变量得到1km的遥感预测降水；

步骤5、计算各气象台站观测降水与1km遥感预测降水的误差；

步骤6、运用克里金插值方法将步骤5的误差插值变为1km分辨率；

步骤7、将步骤4的1km遥感预测降水加上步骤6的误差表面，得到最终的降尺度降水数据。

5.根据权利要求1所述的一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，所述步骤S30的具体计算步骤为：

步骤S301、利用增强型植被指数计算植被状态指数，其计算公式如下：

式中：VCI_i为某年i月的植被状态指数；EVI_i为i月的增强型植被指数值；EVI_max为多年相应月份增强型植被指数的最大值；EVI_min为多年相应月份增强型植被指数的最小值；

步骤S302、利用地表温度指数计算温度条件指数，其计算公式如下：

式中：TCI_i为某年i月的温度条件指数；LST_i为i月的温度条件指数值；LST_max为多年相应月份温度条件指数的最大值；LST_min为多年相应月份温度条件指数的最小值；

步骤S303、利用降尺度降水数据计算降水状态指数，其计算公式如下：

式中：PCI_i为某年i月的降水状态指数；TRMM_i为i月的降尺度降水数据值；TRMM_max为多年相应月份降尺度降水数据的最大值；TRMM_min为多年相应月份降尺度降水数据最小值；

步骤S304、根据不同土地覆盖类型对干旱的响应程度对土地覆盖类型按分等级赋值法进行重分类，得到重分类后的地表覆盖类型；

步骤S305、对目标区域SRTM-DEM数据和重分类后的地表覆盖类型数据进行极差标准化处理，统一量纲得到DEM数据，其计算公式如下：

式中：i为经度；j为纬度；x_ij′为(i,j)为处像元标准化后的数据；x_ij为(i,j)为处像元原始值；x_max为同期影像像元最大值；x_min为同期影像像元最小值。

6.根据权利要求5所述的一种适宜于高海拔地区的干旱遥感监测方法，其特征在于，所述步骤S40中具体过程为：

步骤S401、利用空间主成分分析工具对植被状态指数、温度条件指数、降水状态指数、地表覆盖类型和DEM数据五个因子数据进行主成分分析，得出单月计算的主成分波段数据和对应的计算结果文本数据：

步骤S402、根据上述数据计算干旱综合监测指数，其计算公式如下：

式中：PC_i为第i个主成分；α_i为第i个主成分对应的贡献率的值；

步骤S403、最后根据干旱综合监测指数划分干旱等级，其中：

当干旱综合监测指数≤0.2时，干旱等级为极端干旱；

当0.2＜干旱综合监测指数≤0.3时，干旱等级为重度干旱；

当0.3＜干旱综合监测指数≤0.4时，干旱等级为中度干旱；

当0.4＜干旱综合监测指数≤0.5时，干旱等级为轻度干旱；

当0.5＜干旱综合监测指数≤0.6时，干旱等级为正常；

当干旱综合监测指数＞0.6时，干旱等级为湿润。