一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法

摘要

本发明公开了一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法，包括如下步骤：(1)数据准备；(2)陆地表面温度(LST)数据重构；(3)农作物种植区归一化植被指数-陆地表面温度(NDVI-LST)特征空间构建；(4)TVDI计算；(5)基于TVDI的旱灾等级监测。本发明提出了基于多年背景值与区域波动值的LST数据重构方法，并针对农作物构建耕地区域多年NDVI-LST特征空间，计算了温度植被干旱指数，设计了基于监督分类思想的农作物旱灾等级监测模型进行旱灾等级遥感监测，该模型能够实时、较为准确地反映不同条件下作物受到的干旱胁迫程度，在农业旱灾的监测、预警及防范中具有重要意义。

说明书

技术领域

本发明涉及一种以MODIS遥感数据为主要数据，通过重构陆地表面温度数据并结合归一化植被指数建立农作物温度干旱植被指数，基于农作物温度植被干旱指数的农业旱灾等级监测模型进行农业干旱等级监测方法，相似度评估模型实现自然灾害灾情快速评估的方法，具体为一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法。

背景技术

旱灾是人类生产生活中最主要的自然灾害之一，与其他自然灾害相比，旱灾具有发生频率高、持续时间长、波及范围广的特点。干旱由水分亏缺所引起，在承灾体存在的条件下，较大范围的长时间干旱就可引发旱灾。旱灾是一个复杂、累积的过程，其表现主要为地表水干涸、地下水位下降、地表植被生长受到影响等，其影响主要有农作物减产或绝收，人畜饮水困难，以及持续高温缺水所导致的工业生产受阻。我国是人口和农业大国，粮食安全问题是关乎社会稳定及经济发展的重要问题。在我国旱灾对农业生产带来的危害最为严重。自1949年以来，由干旱引起的粮食受灾面积，平均每年可达3.2亿亩，其中成灾率近30％，不仅给我国农业带来损失，也严重制约了社会和经济发展。

旱灾成因多，机理复杂，具有波及范围大，持续时间长，空间及时间上变异性强等特点，目前在农业旱灾监测上没有成熟有效的方法；传统的气象干旱监测：气象站点稀疏，且分布不均，数据空间精度不高，不能有效的将旱情与农业灾害建立对应关系；农业旱灾灾情统计：虽然可以较精确的得到农业干旱的受灾程度，但统计数据获取周期较长，耗费大量人力物力，且数据多以行政区划为单位，空间精度有限，此外，还存在统计口径不统一等问题。针对上述问题，遥感技术具有多平台、多传感器、高分辨率等特点，可以快速、高精度地获取旱情信息，在短时间内可对受灾区域进行重复性观察；根据不同波段的地物光谱特征，遥感技术可以直接或间接地得到地表水分变化、植被生长、土壤墒情等信息。相比气象干旱监测与农业旱灾灾情统计，农业旱灾遥感监测可以高精度、大范围、快速地获取土壤水分变化及作物生长状态，能够实时、准确地反映不同条件下作物受到的干旱胁迫程度，在农业旱灾的监测、预警及防范中具有重要意义。

发明内容

遥感技术可以快速、多次、大范围地获取地表信息，在旱情监测中具有重要意义。其中利用光学与热红外遥感数据的温度植被干旱指数(TVDI)，有效地结合了植被状态的指示因子归一化植被指数(NDVI)，与地表水热平衡的重要参数陆地表面温度(LST)，在遥感干旱监测研究中应用广泛。但TVDI在实际农业旱情监测中存在如下问题：(1)缺失值较多使LST数据在时空上不连续，无法进一步提取TVDI；(2)在NDVI-LST特征空间构建中存在非耕地像元的影响，且传统的TVDI在多年数据间缺乏可比性。针对如上两个问题，本发明提出了基于多年背景值与区域波动值的LST数据重构方法，相比现有重构方法，可重构单景影像中较大区域的空间缺值，和单像元长时间序列的数据缺失，并能在一定程度上保留LST的变化细节；以农作物为研究对象，构建耕地区域多年NDVI-LST特征空间，并提出农作物温度植被干旱指数(C-TVDI)，在此基础之上，计算监测区农作物关键生长期的C-TVDI，并设计基于监督分类思想的农作物旱灾等级监测模型进行河南省旱灾等级遥感监测，为抗灾救灾决策提供参考。

为实现上述目的本发明采用如下技术方案：一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法，包括如下步骤：

(1)数据的准备：

本发明利用遥感技术获取的数据进行农业旱灾监测研究，其中所用遥感数据为NASA(http://ladsweb.nascom.nasa.gov/)提供的MODISLST产品、植被指数产品、土地利用覆盖产品和DEM数据，以及耕地复种指数数据；统计数据为中国统计年鉴中的耕地面积数据、受灾面积数据和中国农业部种植业管理司(http://www.zzys.moa.gov.cn/)提供的农作物播种面积数据和农作物物候历。

(2)基于背景值与波动值的陆地表面温度数据重构：

基于获取的MODISLST数据产品依据相应的方法对陆地表面温度数据(下面简称LST)进行重构。

(3)农作物种植区归一化植被指数-陆地表面温度特征空间的构建：

提取干旱监测区域耕地，利用多年同期数据共同构建农作物各生长期农作物种植区归一化植被指数-陆地表面温度(下面简称NDVI-LST)散点图，并拟合各期于湿边方程，构建构建NDVI-LST特征空间。

(4)农作物温度植被干旱指数的计算：

基于步骤三的结果，利用Price于1990年提出了温度植被干旱指数计算监测区域农作物种植区农作物温度植被干旱指数(下面简称C-TVDI)。

(5)基于农作物温度植被干旱指数的旱灾等级监测：

基于有监督分类器的设计思想，通过历史数据确定监测区域基于农作物温度植被干旱指数的旱灾等级监测模型的参数，进行旱灾的监测。

作为本发明进一步的方案，所述步骤(2)的基于背景值与波动值的陆地表面温度数据重构包括重构方法和重构操作步骤，具体方法和操作流程如下：

(1)LST重构方法

本发明基于Cressman客观分析法中利用初始场值与订正值共同逼近观测值的思想，将像元LST的多年平均水平，即背景值视作像元的初始场值；将像元LST的波动值，即周边影响半径内像元的观测值与背景值之差作为订正值，以此对插值像元LST进行一次订正插值。像元a(i，j)插值的具体算法可表示如下：

LST_insert＝LST_background+LST_variance式1-1

$>{\mathrm{LST}}_{\mathrm{background}}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{\mathrm{LST}}_{\mathrm{Gaussian}}$>式1-2

$>{\mathrm{LST}}_{\mathrm{variance}}=\frac{\Sigma W_K\Delta {\mathrm{LST}}_K}{\Sigma W_K}$>式1-3

式中，LST_insert为a点LST插值结果，LST_background为a点多年背景值。

由于LST数据中存在缺失值和低质量数据，在计算背景值时应将这些点去除，因此需要对LST时间序列数据进行重构。本发明选择所需设定参数较少的非对称高斯函数拟合法对LST时间序列数据进行重构，并对拟合后的时间序列数据LST_Gaussian求取各期多年背景值，n为数据集中某一期数据所含年数，如上式1-2所示；若某些像元在时间序列数据中缺失值或低质量数据过多，无法进行拟合，该像元值则由周边影响半径内同类地表覆盖类型(在历史年份中至少超过半数地表覆盖类型相同)像元的平均值替代，其中对于缺失像元与周边像元存在高程差异的问题，本发明中利用海拔每升高1000m温度下降约6K的关系，去除参与计算像元的高程对LST的影响。

LST_variance为a点处的波动值，由地表覆盖类型和影响半径决定，ΔLST为影响半径内高质量像元的观测值与背景值之差，K为影响半径内相同地表覆盖类型的高质量像元个数，W_K为其对应权重，由下式计算得到：

$>W_{\mathrm{ijK}}=\left(\begin{array}{cc}\frac{R^2-d_{\mathrm{ijK}}^2}{R^2+d_{\mathrm{ijK}}^2}& (d_{\mathrm{ijK}}<R)\\ 0& (d_{\mathrm{ijK}}\ge R)\end{array}\right)$>式1-4

式中，d_ijk为插值像元到相同地表类型的高质量像元的距离；R为影响半径，由于温度在水平距离6000m范围内的变化一般小于0.6K^[57]，因此，本文将影响半径R设为3，即在7x7窗口内的高质量像元参与插值计算。

(2)LST重构操作步骤

1)数据预处理。本发明选取MODIS产品用于干旱监测区LST数据重构。重构中所用到的MOD11A2、MOD12Q1产品首先通过MODIS重投影工具(MODISReprojectionTool，MRT)重投影为WGS84坐标系统，并提取产品中的昼夜LST数据、昼夜LST质量控制文件和LAI/fPAR体系土地利用/覆盖分类结果；研究中还用到由NASA提供的数字高程数据(DigitalElevationModel，DEM)，先将DEM数据同MODIS数据进行配准，再将DEM数据的空间分辨率从30m重采样为1000m。

由于研究所用数据的种类及数量较多，在处理前对上述研究数据进行剪裁，考虑重构算法中涉及空间窗口卷积运算，因此选择矩形掩膜文件对数据进行裁切，空间分辨率为1000m，尽可能确保监测区域处于矩形区域中央。经上述预处理后所用研究数据为：昼夜LST数据、昼夜LST质量控制文件、DEM数据和MOD12Q1数据，对以上数据进行波段融合，构成昼/夜LST数据集、昼/夜LST质量控制文件数据集、地表覆盖类型数据集和DEM数据。

2)计算LST多年背景值。将昼夜LST数据集分别进行非对称高斯函数拟合。之后，检验拟合后的时间序列数据集：由于，数据集中LST的像元值为开式温度，所以将未拟合成功的像元所在的时间序列所有值设为0；其后，对所有像元计算该像元同期多年平均值，作为该期背景值，获得背景值数据集，昼夜各46个波段。对背景值数据集中像元值为0的像元，在该像元周围7x7窗口内进行插值：首先，参照地表覆盖类型数据集，选取7x7窗口内与待插值像元地表覆盖类型相同的像元(在历史年份中与待插值像元类型相同次数超过半数即认为相同地表覆盖类型)，若没有则选取全部像元；随后，利用Cressman客观分析法中的权重计算方法根据选取像元与待插值像元的距离计算各像元权重；之后，利用DEM数据，计算选取像元与待插值像元的高程差，利用海拔每升高1000m，温度下降6k的关系，将所有选取像元的LST统一到待插值像元所在高程温度；最后，对窗口内选取像元“找平”后的LST与其权重进行加权求和，获取时空连续的LST多年背景值数据集。

3)利用昼夜LST质量控制文件数据集，分别对昼夜LST数据进行逐项元筛选：依据昼夜LST质量控制文件，仅保留各波段中高质量像元值，将其余像元值设为0，获得昼夜待插值LST数据集。

4)利用计算所得LST背景值数据集、待插值LST数据集以及地表覆盖类型数据集进行插值，具体步骤为a.利用地表覆盖类型数据集，选取像元值为0的像元周边7x7窗口内，与0值像元地表覆盖类型相同像元，若没有则选取全部像元；b.利用Cressman客观分析法中的权重计算方法计算选取像元的权重；c.计算选取像元的实际值(即高质量像元观测值)与背景值之差，作为波动值；d.对窗口内选取像元的波动值与其对应权重进行加权求和，获取最终LST插值结果。

作为本发明进一步的方案，所述步骤(3)农作物种植区归一化植被指数-陆地表面温度特征空间的构建，包括利用提取干旱监测区域耕地，利用多年同期数据共同构建农作物各生长期农作物种植区归一化植被指数-陆地表面温度(下面简称NDVI-LST)散点图，并拟合各期干湿边方程，构建构建NDVI-LST特征空间。具体如下：

(1)本文发明利用MODIS植被指数16天合成产品MOD13A2，经本发明方法重构后的白天LST8天数据以及监测区域农作物复种指数数据，各数据空间分辨率均为1000m。

1)将MODIS植被指数产品MOD13A2用MODIS重投影工具进行投影转换，并提取其中的NDVI数据及质量评估数据(QualityAssurance，QA)；

2)对监测区逐年数据进行耕地像元提取，其中由于在耕地复种指数数据中，选择像元面积与统计年鉴中耕地面积数据最为吻合的像元值，各年面积偏差控制在10％以下，以该复种指数像元值的像元作为监测区域耕地像元，对监测区域耕地进行提取；

3)参照QA数据中NDVI质量的16级分类标准，如表1所示，筛选NDVI数据质量为中等以上的像元，即在QA数据2～5位的值小于4的像元，将NDVI数据中其余像元赋为0；

4)利用Timesat3.1.1对经QA数据筛选后的NDVI数据，进行时间序列重构，重构方法选择Savitzky-Golay滤波法，滤波窗口设为5；

5)将白天LST8天数据每2景合并为1景16天LST数据，合并中原始LST数据的高质量像元替代插值像元，当两像元同为高质量像元或插值像元时，合并算法参考MOD11A2用户手册，对两像元值求取平均值，作为合并数据像元值。

(2)首先，对多年数据按期分组，将NDVI、LST数据分成若干组，各组中包含相同数量同期NDVI、LST数据，对各组数据的所有农作物像元构建NDVI-LST特征空间；然后，以各组数据中的NDVI值最大、最小值之差的百分之一为步长，逐步长获取不同NDVI值所对应的组内LST最大、最小值；最后利用NDVI与LST最大、最小值拟合各组数据中的干湿边方程。

作为本发明进一步的方案，所述步骤(4)农作物温度植被干旱指数的计算，包括基于Price于1990年提出了温度植被干旱指数计算监测区域C-TVDI值。具体的温度植被干旱指数的计算方法表示如下：

$>\mathrm{TVDI}=\frac{T_s-T_{s\_\min }}{T_{s\_\max }-T_{s\_\min }}$>式1-5

其中，

T_{s_max}＝a₁+b₁NDVI式1-6

T_{s_min}＝a₂+b₂NDVI式1-7

式中，T_s为像元(i，j)的LST值，T_{s_max}与T_{s_min}分别为像元(i，j)的NDVI值对应的干边和湿边上的LST值，干边由对特征空间中各类NDVI值对应的最大值构成的散点图，进行线性回归求得，湿边则由对应的最小值散点图的线性回归求得；a₁，a₂，b₁，b₂分别为干湿边方程中的待定系数。

作为本发明进一步的方案，所述步骤(5)基于农作物温度植被干旱指数的旱灾等级监测，包括基于有监督分类器的设计思想，通过历史数据确定监测区域基于农作物温度植被干旱指数的旱灾等级监测模型的参数，进行旱灾的监测。具体如下：

基于有监督分类器的设计思想，以往年作物生长期的遥感数据与各年实测旱灾等级数据作为训练样本，通过对训练样本的学习，确定旱灾等级监测模型中的参数；之后再利用最新一年的作物生长期遥感数据，对该年年终旱灾等级进行监测。该模型主要由受灾面积估算函数、参数寻优和旱灾等级监测三部分组成：

(1)受灾面积估算函数：受灾面积定义为遭受各种自然灾害的农作物播种面积，且在同年遭受几种或几次自然灾害的，以其中危害最大的一次计算受灾面积，不重复计灾，旱灾受灾面积则为该年一次或多次遭受旱灾影响的最大农作物播种面积。考虑受灾面积值为年内旱灾对农作物的最大影响面积，考虑监测区如果为多次复种的种植模式，且农作物经收割之后一次旱灾的影响结束，因此，将受灾面积估算函数分为不同种植期，设定如下：

S_i＝Max[Mean(s_ik，...，s_il)，...，Mean(s_im，...，s_in)]式1-8

式中，S_i为第i年受灾面积估算值；k～1为每年第一期收割作物关键生长期，m～n为每年最后一期收割作物关键生长期，s_ij(j＝k～1)为第i年第一期收割作物各关键生长期的受旱面积，s_ij(j＝m～n)为最后一期收割作物各生长期的受旱面积；由于从受旱到旱灾的发生需要一定时间的积累，单期数据的旱情未必引发旱灾，而短期内干旱的空间变异性不大，因此将各期受旱面积的平均值作为受灾面积；全年的受灾面积值则为各期受灾面积的最大值。

受旱面积s_ij(j＝k～1或m～n)由第i年农作物关键生长期农业遥感指数数据和该生长期阈值参数t_j(j＝k～1或m～n)决定。以农作物温度植被干旱指数为例，指数值越接近1则表示越干旱，则第i年第j个关键生长期的受旱面积S_ij为该年该期数据中C-TVDI值大于该期阈值t_j的像元所占面积，阈值t_j为模型参数，参数与所用生长期一一对应。

(2)参数寻优，即利用寻优算法对设定的目标函数求取最值的过程。旱灾等级监测模型以多年遥感数据与农业受灾面积数据作为训练样本进行模型的参数寻优。模型参数的最优解可认为是在该组参数下，计算得到多年受灾面积估算值与多年训练样本中实际受灾面积的最小平均偏差，因此将目标函数设定为：

$>\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}|S_i-S_{i0}|$>式1-9

式中，S_i为第i年受灾面积估算值，S_i0为第i年实际受灾面积值，n为训练样本个数，对目标函数求得最小值时对应的模型参数t_j(j＝k～1或m～n)即为最优参数。

本研究选用网格搜索算法进行参数寻优，网格搜索算法的寻优速度较快，可以获得全局最优解，不会陷入局部最优解。但在搜索前要设定参数范围，在参数范围较大，搜索步长较小的搜索中，训练耗时较长。

参数寻优过程表述如下：(i)将训练样本中的遥感数据进行农作物遥感指数提取，并选择参与模型计算的作物生长期；(ii)设定最小目标函数值、各生长期参数t_i0的搜索范围、搜索步长和参数t_i0的起始值；(iii)将提取的农业遥感指数数据输入到受灾面积估算函数中，获得受灾面积估算值；(iv)将受灾面积估算值与训练样本中受灾面积数据带入目标函数中，计算目标函数值，若该值小于当前最小目标函数值，则保留当前参数，并将最小值更新；(v)通过网格搜索算法获得新的参数值，并重复上述(iii)、(iv)步骤，若全部参数均已遍历则完成寻优，记录最优参数。

(3)旱灾等级监测，即通过待监测年农作物生长期遥感数据，对该年旱灾等级进行监测。由于本研究中参与训练的样本数较少，训练样本对不同程度灾情的代表性有限，因此利用训练样本受灾面积估算值与实际受灾面积间的线性关系，对待监测年的受灾面积估算值进行调整。旱灾等级监测的具体步骤如下：(i)计算训练样本中实际受灾面积与最优参数下受灾面积估算值的线性回归方程；(ii)提取待监测年农作物关键生长期遥感指数；(iii)利用受灾面积估算函数求取最优参数下的受灾面积估算值；(iv)将上一步中估算值作为线性方程自变量，获得调整后受灾面积估算值；(v)利用调整后受灾面积估算值计算受灾率(I)，并参照旱灾等级标准划分表，获得旱灾等级监测结果。旱灾等级标准划分表如下：

表2旱灾等级标准划分表

与现有技术相比，本发明进行了如下改进：(1)提出了基于多年背景值与区域波动值的LST数据重构方法，相比现有重构方法，可重构单景影像中较大区域的空间缺值，和单像元长时间序列的数据缺失，并能在一定程度上保留LST的变化细节；(2)以农作物为研究对象，构建耕地区域多年NDVI-LST特征空间，并提出农作物温度植被干旱指数(C-TVDI)，在此基础之上，计算监测区农作物关键生长期的C-TVDI，并设计基于监督分类思想的农作物旱灾等级监测模型进行河南省旱灾等级遥感监测，为抗灾救灾决策提供参考。

附图说明

图1是根据本发明说明书和具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的陆地表面温度数据重构技术路线图；

图2是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的2005年1月7日白天河南陆地表面温度数据重构结果对比图；

图3是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的2009年1月1日夜间河南陆地表面温度数据重构结果对比图；

图4是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的陆地表面温度数据重构原始值重构结果散点图；

图5是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的多年农作物NDVI-LST部分散点图；

图6是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的C-TVDI与农业受灾面积行政区划分布对比图；

图7是根据本发明具体实施例的一种基于温度植被干旱指数(TVDI)的农业旱灾等级监测方法的旱灾等级监测流程图；

具体实施方式

下面结合附图和具体实施例对本发明作进一步阐述。

本案例以河南为研究区，主要用到植被指数数据产品、LST数据产品、土地覆盖/土地覆盖变化数据产品。本文所用的各类MODIS产品如下表：

表1MODIS陆地产品

耕地复种指数数据主要用于获得河南耕地像元的空间分布，该数据由刘建红博士提供，数据空间分辨率为500m，本案例所用为河南地区2001-2011年共11景数据，像元值为每年该像元的农作物种植次数，非耕地像元标记为0，耕地复种指数数据具体提取算法如下：

(1)根据多年平均气象资料将全国划分一熟区、两熟区和三熟区；

(2)在每个熟制区内选择适当的训练样本，根据训练样本确定每个熟制区内农作物的最短生长季长度、最长生长季长度、最小生长幅度；

(3)利用MODIS数据获得每个像元增强型植被指数(EnhancedVegetationIndex，EVI)时间序列曲线，并提取植被生长季个数、生长季长度、生长幅度等物侯参数；

(4)将像元提取的物侯参数与像元所在熟制区内的农作物的最短生长季长度、最长生长季长度、最小生长幅度进行比较，判别一个植被生长季是否属于农作物生长季，如果是则保留，如果不是则删除；

(5)最终得到的农作物生长季个数即为该像元的复种指数。

在本案例中，还用到了统计年鉴数据、农事历数据和各类农作物播种面积数据3类统计数据；其中农事历数据和各类农作物播种面积数据，用于确定主要农作物的关键生长期，并筛选与农作物关键生长期一致的遥感数据用于农业旱灾监测研究；统计年鉴数据构建遥感旱灾等级监测模型与对模型监测结果进行检验。研究所用统计年鉴数据包括河南省各年旱灾受灾面积数据和耕地面积数据，该数据全部来自2001-2011年中国统计年鉴；农事历数据和各类农作物播种面积数据来自于中国农业部种植业管理司(http://www.zzys.moa.gov.cn/)，根据各类农作物的播种面积，河南省主要农作物为冬小麦、夏玉米、中稻、大豆、棉花、油菜和花生，各类主要农作物在研究年限内的农事历如下表：

本案例主要流程包括：1)利用多年背景值与区域波动值的LST数据重构方法，对河南2000-2011年LST数据进行重构，随机抽选了5000个高质量像元，对重构数据进行精度验证；2)利用重构后的LST数据与MODISNDVI数据，构建河南地区农作物NDVI-LST特征空间；3)以2001-2007年农作物温度植被干旱指数数据与实际受灾面积数据作为训练样本，以2008-2011年农作物温度植被干旱指数数据、实际受灾面积数据以及旱灾等级数据作为检验样本，对农作物旱灾监测模型进行评价。

本案例选取MODIS产品用于河南地区LST数据重构研究。重构中所用到的MOD11A2、MOD12Q1产品首先通过MODIS重投影工具(MODISReprojectionTool，MRT)重投影为WGS84坐标系统，并提取产品中的昼夜LST数据、昼夜LST质量控制文件和LAI/fPAR体系土地利用/覆盖分类结果；研究中还用到由NASA提供的数字高程数据(DigitalElevationModel，DEM)，先将DEM数据同MODIS数据进行配准，再将DEM数据的空间分辨率从30m重采样为1000m。

由于研究所用数据的种类及数量较多，在处理前对上述研究数据进行剪裁，考虑重构算法中涉及空间窗口卷积运算，因此选择矩形掩膜文件对数据进行裁切，大小设为641x622，空间分辨率为1000m，研究区河南处于矩形区域中央。经上述预处理后所用研究数据为：昼夜LST数据各506景，昼夜LST质量控制文件各506景，其中每年46景，共11年；DEM数据1景；MOD12Q1共11景，每年1景，对以上数据进行波段融合，构成昼/夜LST数据集、昼/夜LST质量控制文件数据集、地表覆盖类型数据集和DEM数据，便于进一步重构算法的处理。

本案例陆地表面温度数据重构(LST)包括三个步骤(流程图见说明书附图1)：

(1)计算LST多年背景值。将昼夜LST数据集分别进行非对称高斯函数拟合，该步骤由Timesat3.1.1实现。之后，检验拟合后的时间序列数据集：由于，数据集中LST的像元值为开式温度，所以将未拟合成功的像元所在的时间序列所有值设为0；其后，对所有像元计算该像元同期多年平均值，作为该期背景值，获得背景值数据集，昼夜各46个波段。对背景值数据集中像元值为0的像元，在该像元周围7x7窗口内进行插值：首先，参照地表覆盖类型数据集，选取7x7窗口内与待插值像元地表覆盖类型相同的像元(在11年中与待插值像元类型相同次数超过6次即认为相同地表覆盖类型)，若没有则选取全部像元；随后，利用Cressman客观分析法中的权重计算方法根据选取像元与待插值像元的距离计算各像元权重；之后，利用DEM数据，计算选取像元与待插值像元的高程差，利用海拔每升高1000m，温度下降6k的关系，将所有选取像元的LST统一到待插值像元所在高程温度；最后，对窗口内选取像元“找平”后的LST与其权重进行加权求和，获取时空连续的LST多年背景值数据集。

(2)利用昼夜LST质量控制文件数据集，分别对昼夜LST数据进行逐项元筛选：按照3.1节中像元质量等级的划分，仅保留各波段中高质量像元值，将其余像元值设为0，获得昼夜待插值LST数据集，各506景。

(3)利用计算所得LST背景值数据集、待插值LST数据集以及地表覆盖类型数据集进行插值，具体步骤为1.利用地表覆盖类型数据集，选取像元值为0的像元周边7x7窗口内，与0值像元地表覆盖类型相同像元，若没有则选取全部像元；2.利用Cressman客观分析法中的权重计算方法计算选取像元的权重；3.计算选取像元的实际值(即高质量像元观测值)与背景值之差，作为波动值；4.对窗口内选取像元的波动值与其对应权重进行加权求和，获取最终LST插值结果(说明书附图2)。

本案例为验证基于背景值与波动值LST数据重构算法的插值精度，设计精度验证实验，在原始昼夜LST数据集中随机选取了250景LST数据，在选取的250景数据中再随机选取5000个高质量像元，将像元值赋为0，生成昼夜LST验证数据集，验证数据集与原始LST数据集除像元值被修改外，其余均相同。将验证数据集进行基于背景值与波动值LST数据重构，并将重构结果中的插值结果与原始像元值进行对比。对比结果中，LST插值与原始像元值的最大偏差为15.44K，平均偏差为0.81K，偏差超过2K的像元占像元总数的5.2％；可见，基于背景值与波动值LST数据重构算法的插值精度较高，可对MODISLST数据进行较好的修复。说明书附图4为插值结果与原始数据的散点图，原始像元值与插值结果有较强的线性相关性，其中，线性斜率为0.986，截距为1.21，R2为0.960。(说明书附图2、图3、图4)

本案例利用MODIS植被指数16天合成产品MOD13A2，经上述重构后的白天LST8天数据以及河南地区2001-2011年农作物复种指数数据，各数据空间分辨率均为1000m。数据处理过程如下：

(1)将MODIS植被指数产品MOD13A2用MODIS重投影工具进行投影转换，并提取其中的NDVI数据及质量评估数据(QualityAssurance，QA)；

(2)对研究区逐年数据进行耕地像元提取，其中由于在耕地复种指数数据中，像元值为2的像元面积与河南统计年鉴中耕地面积数据最为吻合，各年面积偏差在5％～10％不等，因此以复种指数值为2的像元作为研究区耕地像元，对研究区耕地进行提取；

(3)参照QA数据中NDVI质量的16级分类标准，如表6所示，筛选NDVI数据质量为中等以上的像元，即在QA数据2～5位的值小于4的像元，将NDVI数据中其余像元赋为0；

(4)利用Timesat3.1.1对经QA数据筛选后的NDVI数据，进行时间序列重构，重构方法选择Savitzky-Golay滤波法，滤波窗口设为5；

(5)将白天LST8天数据每2景合并为1景16天LST数据，合并中原始LST数据的高质量像元替代插值像元，当两像元同为高质量像元或插值像元时，合并算法参考MOD11A2用户手册，对两像元值求取平均值，作为合并数据像元值。

本案例利用预处理后的NDVI、LST数据，对河南农作物种植区构建NDVI-LST特征空间。首先，对多年数据按期分组，将NDVI、LST数据分成23组，各组中包含11景同期NDVI、LST数据，对各组数据的所有农作物像元构建NDVI-LST特征空间；然后，以各组数据中的NDVI值最大、最小值之差的百分之一为步长，逐步长获取不同NDVI值所对应的组内LST最大、最小值；最后利用NDVI与LST最大、最小值拟合各组数据中的干湿边方程，部分结果如说明书附图5所示。

在多期结果中，多年农作物NDVI-LST散点分布大致可分为三类：其中，大部分时期的多年农作物NDVI-LST散点分布符合理论上的三角分布关系，如说明书附图5(a)第11期(6月下旬)，5(b)第18期(10月中旬)，5(c)第21期(11月上旬)，其中干边斜率均小于0，而湿边并不全为平行于NDVI轴的直线，如在说明书附图5第11期、18期中，湿边斜率值均大于0，在说明书附图5第21期中斜率小于0，从整体分布看，湿边斜率值均大于干边的斜率值，干湿边随NDVI增长逐渐靠拢；部分时期的多年农作物NDVI-LST散点图分布与理论三角形分布相反，如说明书附图5(d)第8期(5月上旬)、5(e)第13期(7月下旬)所示，其中干边斜率均大于0，湿边斜率均小于0，干湿边随NDVI的增长逐渐背离；另有个别时期的多年农作物NDVI-LST散点图分布为一矩形，干湿边相互平行，且斜率均接近于0，如说明书附图5(f)第3期(2月中旬)。

在上述结果中，多年农作物NDVI-LST散点分布的区别主要由散点图各期所对应的农作物生长期的不同决定：呈理论三角形分布的散点图，如11期、18期和21期，时间分别为6月下旬、10月中旬和11月上旬，分别对应为秋粮播种/出苗期、秋粮收获期和夏粮播种期，在上述三个时段，农作物种植区域NDVI值分布较为均匀，裸地、低覆盖作物和高覆盖作物共存，因此NDVI-LST散点图与理论三角形相近；而呈反三角形的散点图，如8期、13期，时间分别为5月上旬和7月下旬，分别对应为夏粮生长期和秋粮生长期，这两个时期农作物长势达到最大，在农作物种植区的绝大部分像元为高覆盖作物，仅有少数像元为裸地或为裸地与作物的混合像元，NDVI-LST散点图中，点多集中于NDVI的中高值区域，NDVI低值区域只有少数点，对应的LST值不具备代表性且变化差异有限，因此干湿边在NDVI低值区域更为接近；呈矩形分布的散点图，如3期，时间为2月中旬，该时期为夏粮冬小麦的越冬期，此时农作物种植区域像元有一定的植被覆盖度，散点图中点集中于NDVI的中低值区，但处于越冬期作物蒸腾作用很弱，不具备调节冠层温度的能力，因此干湿边均呈平行于NDVI轴的直线。

本案例利用上述过程求得的干湿边方程进一步计算河南多年农作物种植区域的C-TVDI值，计算方法如下：

$>\mathrm{TVDI}=\frac{T_s-T_{s\_\min }}{T_{s\_\max }-T_{s\_\min }}$>式1

其中，

T_{s_max}＝a₁+b₁NDVI式2

T_{s_min}＝a₂+b₂NDVI式3

式中，T_s为像元(i，j)的LST值，T_{s_max}与T_{s_min}分别为像元(i，j)的NDVI值对应的干边和湿边上的LST值，干边由对特征空间中各类NDVI值对应的最大值构成的散点图，进行线性回归求得，湿边则由对应的最小值散点图的线性回归求得；a₁，a₂，b₁，b₂分别为干湿边方程中的待定系数。

本案例计算河南多年农作物种植区域的C-TVDI值部分结果与同年农业受灾面积行政区划分布图见说明书附图6.

本案例河南省各市受灾面积空间分布与C-TVDI指数空间分布对比结果(说明书附图6)可见，C-TVDI在整体上能够正确地反映出农业旱灾的空间分布，但在地市级区域上的准确性不足，其原因可能为旱灾灾情受多种因素共同决定；为进一步将C-TVDI用于农业旱灾遥感监测，并将遥感数据与农业旱灾灾情建立关系，本案例提出一种基于多期遥感数据(C-TVDI)的旱灾等级监测模型，并结合河南省旱灾受灾面积数据对模型监测结果进行验证。

本案例旱灾等级监测模型是基于有监督分类器的设计思想，以往年作物生长期的遥感数据与各年实测旱灾等级数据作为训练样本，通过对训练样本的学习，确定旱灾等级监测模型中的参数；之后再利用最新一年的作物生长期遥感数据，对该年年终旱灾等级进行监测。该模型主要由受灾面积估算函数、参数寻优和旱灾等级监测三部分组成：

(1)受灾面积估算函数：受灾面积定义为遭受各种自然灾害的农作物播种面积，且在同年遭受几种或几次自然灾害的，以其中危害最大的一次计算受灾面积，不重复计灾^[65]，旱灾受灾面积则为该年一次或多次遭受旱灾影响的最大农作物播种面积。考虑受灾面积值为年内旱灾对农作物的最大影响面积，研究区河南主要为两次复种的种植模式，且农作物经收割之后一次旱灾的影响结束，因此，将受灾面积估算函数分为夏粮受灾面积和秋粮受灾面积两部分，设定如下：

S_i＝Max[Mean(s_ik，...，s_il)，Mean(s_im，...，s_in)]式4

式中，S_i为第i年受灾面积估算值；k～1为夏粮作物关键生长期，m～n为秋粮作物关键生长期，s_ti(j＝k～1)为第i年夏粮作物各关键生长期的受旱面积，s_ij(j＝m～n)为秋粮作物各生长期的受旱面积；由于从受旱到旱灾的发生需要一定时间的积累，单期数据的旱情未必引发旱灾，而短期内干旱的空间变异性不大，因此将夏/秋粮各期受旱面积的平均值作为夏/秋粮的受灾面积；全年的受灾面积值则为夏粮与秋粮受灾面积的最大值。

受旱面积s_ij(j＝k～1或m～n)由第i年农作物关键生长期农业遥感指数数据和该生长期阈值参数t_j(j＝k～1或m～n)决定。以农作物温度植被干旱指数为例，指数值越接近1则表示越干旱，则第i年第j个关键生长期的受旱面积S_ij为该年该期数据中C-TVDI值大于该期阈值t_j的像元所占面积，阈值t_j为模型参数，参数与所用生长期一一对应。

(2)参数寻优，即利用寻优算法对设定的目标函数求取最值的过程。旱灾等级监测模型以多年遥感数据与农业受灾面积数据作为训练样本进行模型的参数寻优。模型参数的最优解可认为是在该组参数下，计算得到多年受灾面积估算值与多年训练样本中实际受灾面积的最小平均偏差，因此将目标函数设定为：

$>\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}|S_i-S_{i0}|$>式5

式中，S_i为第i年受灾面积估算值，S_i0为第i年实际受灾面积值，n为训练样本个数，对目标函数求得最小值时对应的模型参数t_j(j＝k～1或m～n)即为最优参数。

常用的寻优算法有遗传算法、模拟退火算法、粒子群算法以及网格搜索算法。由于本文模型涉及参数较少，且训练样本量较小，因此选用网格搜索算法进行参数寻优，网格搜索算法的寻优速度较快，可以获得全局最优解，不会陷入局部最优解。但在搜索前要设定参数范围，在参数范围较大，搜索步长较小的搜索中，训练耗时较长。

参数寻优过程表述如下：(i)将训练样本中的遥感数据进行农作物遥感指数提取，并选择参与模型计算的作物生长期；(ii)设定最小目标函数值、各生长期参数t_i0的搜索范围、搜索步长和参数t_i0的起始值；(iii)将提取的农业遥感指数数据输入到受灾面积估算函数中，获得受灾面积估算值；(iv)将受灾面积估算值与训练样本中受灾面积数据带入目标函数中，计算目标函数值，若该值小于当前最小目标函数值，则保留当前参数，并将最小值更新；(v)通过网格搜索算法获得新的参数值，并重复上述(iii)、(iv)步骤，若全部参数均已遍历则完成寻优，记录最优参数。

(3)旱灾等级监测，即通过待监测年农作物生长期遥感数据，对该年旱灾等级进行监测。由于本研究中参与训练的样本数较少，训练样本对不同程度灾情的代表性有限，因此利用训练样本受灾面积估算值与实际受灾面积间的线性关系，对待监测年的受灾面积估算值进行调整。旱灾等级监测的具体步骤如下：(i)计算训练样本中实际受灾面积与最优参数下受灾面积估算值的线性回归方程；(ii)提取待监测年农作物关键生长期遥感指数；(iii)利用受灾面积估算函数求取最优参数下的受灾面积估算值；(iv)将上一步中估算值作为线性方程自变量，获得调整后受灾面积估算值；(v)利用调整后受灾面积估算值计算受灾率(I)，并参照旱灾等级标准划分表，获得旱灾等级监测结果^[64]。旱灾等级标准划分表如下：

表4旱灾等级标准划分表

本案例旱灾登记监测利用2001-2011年河南省农作物温度植被干旱指数数据与中国统计年鉴河南省受灾面积、耕地面积数据，对旱灾等级监测模型进行评价；根据河南夏秋粮主要作物，冬小麦和夏玉米在各时段的对水分的不同需求，将冬小麦的拔节、抽穗和灌浆期，以及夏玉米的抽雄和乳熟期作为旱灾等级监测模型中的农作物关键生长期，该生长期对应每年中的6期遥感数据，每期数据对应一个阈值参数tj(j＝1～6)。为确定模型的6个阈值参数，将2001-2007年数据作为模型训练样本，进行参数寻优并确定受灾面积初始估算值与实际受灾面积的回归方程；为评价模型的监测结果，将2008-2011年数据作为检验样本，利用模型的最优参数与线性回归方程对各年受灾面积进行估算，并对旱灾等级进行监测，并对比监测结果与实际旱灾等级。具体流程见说明书附图7。

本案例利用训练样本2001-2007年河南省农作物温度植被干旱指数数据与河南受灾面积数据，对旱灾等级监测模型各生长关键期阈值参数寻优结果如下：

表5模型各期阈值寻优结果

利用最优参数求得训练样本受灾面积估算值，并与实际受灾面积值进行线性回归，回归方程及决定系数如下：

S_实测＝0.8448×S_估算+2693.4式6

R²＝0.9370

根据上述参数及回归方程，将检验样本2008-2011年农作物温度植被干旱指数数据分别输入旱灾等级监测模型中，求得检验样本各年受灾面积估算值，并根据旱灾等级标准划分表获得各年旱灾等级监测结果，如下表所示：

表6旱灾等级监测结果

表6中，受灾率为受灾面积与耕地面积之比，受灾面积偏差率为受灾面积估算值与实际值绝对偏差占实际受灾面积的比值。通过表6可见，作为检验样本的2008-2011年，涵盖了从无旱到重旱四个等级，其中受灾面积偏差率随实际旱灾等级的下降而增高，如发生重旱的2009年，实际受灾面积超过1500khm²，而受灾面积偏差率仅为6％，而在无旱灾发生的2010年，实际受灾面积不足100khm²，受灾面积偏差为实际受灾面积4倍多，从受灾面积偏差值可以看出，上述结果是因为各年受灾面积的估算值与实际值均存在200khm²左右的偏差，偏差率则随着实际受灾面积的减少而升高；从旱灾等级的监测结果上看，检验样本的旱灾监测等级与旱灾实际等级均一致，在灾害等级的定性监测上具有一定的可行性。

本案例详细介绍了基于农作物温度植被干旱指数的农业旱灾等级监测技术流程，并以河南省为例，收集了2001年到2011年的该省农业旱灾相关数据，将2008-2011年数据作为检验样本，利用模型的最优参数与线性回归方程对各年受灾面积进行估算，并对旱灾等级进行监测，并对比监测结果与实际旱灾等级。从旱灾等级的监测结果上看，检验样本的旱灾监测等级与旱灾实际等级均一致，在灾害等级的定性监测上具有一定的可行性，能够满足干旱灾害等级监测的基本要求。基于农作物温度植被干旱指数的农业旱灾等级监测方法数据可获取性强，方法简单，操作容易，具有一定业务化应用前景。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。