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

摘要

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

主权项

1.一种基于温度植被干旱指数TVDI的农业旱灾等级监测方法，其特征在于，包括如下步骤：/n一、数据的准备：/n利用遥感技术获取的数据进行农业旱灾监测研究，其中所用遥感数据为NASA，http：//ladsweb.nascom.nasa.gov/提供的MODIS LST产品、植被指数产品、土地利用覆盖产品和DEM数据，以及耕地复种指数数据；统计数据为中国统计年鉴中的耕地面积数据、受灾面积数据和中华人民共和国农业农村部种植业管理司，http://www.zzys.moa.gov.cn/提供的农作物播种面积数据和农作物物候历；/n二、基于背景值与波动值的陆地表面温度数据重构：/n基于获取的MODIS LST产品依据相应的方法对陆地表面温度数据LST进行重构；/n所述步骤二的基于背景值与波动值的陆地表面温度数据重构包括重构方法和重构操作步骤，具体方法和操作流程如下：/n(1)LST重构方法/n基于Cressman客观分析法中利用初始场值与订正值共同逼近观测值的思想，将像元LST的多年平均水平，即背景值视作像元的初始场值；将像元LST的波动值，即周边影响半径内像元的观测值与背景值之差作为订正值，以此对插值像元LST进行一次订正插值；像元a(i，j)插值的具体算法可表示如下：/nLST_insert＝LST_background+LST_variance 式1-1/n /n /n式中，LST_insert为像元a LST插值结果，LST_background为像元a多年背景值；/n由于LST数据中存在缺失值和低质量数据，在计算背景值时应将这些点去除，因此需要对LST时间序列数据进行重构；选择所需设定参数较少的非对称高斯函数拟合法对LST时间序列数据进行重构，并对拟合后的时间序列数据LST_Gaussian求取各期多年背景值，n为数据集中某一期数据所含年数，如上式1-2所示；若某些像元在时间序列数据中缺失值或低质量数据过多，无法进行拟合，该像元值则由周边影响半径内同类地表覆盖类型像元的平均值替代，在历史年份中至少超过半数地表覆盖类型相同，其中对于缺失像元与周边像元存在高程差异的问题，利用海拔每升高1000m温度下降约6K的关系，去除参与计算像元的高程对LST的影响；/nLST_variance为像元a处的波动值，由地表覆盖类型和影响半径决定，ΔLST为影响半径内高质量像元的观测值与背景值之差，K为影响半径内相同地表覆盖类型的高质量像元个数，W_K为其对应权重，由下式计算得到：/n /n式中，d_ijk为插值像元到相同地表类型的高质量像元的距离；R为影响半径，由于温度在水平距离6000m范围内的变化一般小于0.6K，因此，将影响半径R设为3，即在7x7窗口内的高质量像元参与插值计算；/n(2)LST重构操作步骤/n1)数据预处理；选取MODIS产品用于干旱监测区LST数据重构；重构中所用到的MOD11A2、MOD12Q1产品首先通过MODIS重投影工具MODISReprojectionTool，MRT重投影为WGS84坐标系统，并提取产品中的昼夜LST数据、昼夜LST质量控制文件和LAI/fPAR体系土地利用/覆盖分类结果；研究中还用到由NASA提供的数字高程数据Digital ElevationModel，DEM，先将DEM数据同MODIS数据进行配准，再将DEM数据的空间分辨率从30m重采样为1000m；/n由于所用数据的种类及数量较多，在处理前对上述研究数据进行剪裁，考虑重构算法中涉及空间窗口卷积运算，因此选择矩形掩膜文件对数据进行裁切，空间分辨率为1000m，尽可能确保监测区域处于矩形区域中央；经上述预处理后所用研究数据为：昼夜LST数据、昼夜LST质量控制文件、DEM数据和MOD12Q1数据，对以上数据进行波段融合，构成昼夜LST数据集、昼夜LST质量控制文件数据集、地表覆盖类型数据集和DEM数据；/n2)计算LST多年背景值；将昼夜LST数据集分别进行非对称高斯函数拟合；之后，检验拟合后的时间序列数据集：由于，数据集中LST的像元值为开式温度，所以将未拟合成功的像元所在的时间序列所有值设为0；其后，对所有像元计算该像元同期多年平均值，作为该期背景值，获得LST背景值数据集，昼夜各46个波段；对LST背景值数据集中像元值为0的像元，在该像元周围7x7窗口内进行插值：首先，参照地表覆盖类型数据集，选取7x7窗口内与待插值像元地表覆盖类型相同的像元，在历史年份中与待插值像元类型相同次数超过半数即认为相同地表覆盖类型，若没有则选取全部像元；随后，利用Cressman客观分析法中的权重计算方法根据选取像元与待插值像元的距离计算各像元权重；之后，利用DEM数据，计算选取像元与待插值像元的高程差，利用海拔每升高1000m，温度下降6K的关系，将所有选取像元的LST统一到待插值像元所在高程温度；最后，对窗口内选取像元“找平”后的LST与其权重进行加权求和，获取时空连续的LST背景值数据集；/n3)利用昼夜LST质量控制文件数据集，分别对昼夜LST数据进行逐项元筛选：依据昼夜LST质量控制文件，仅保留各波段中高质量像元值，将其余像元值设为0，获得昼夜待插值LST数据集；/n4)利用计算所得LST背景值数据集、昼夜待插值LST数据集以及地表覆盖类型数据集进行插值，具体步骤为a.利用地表覆盖类型数据集，选取像元值为0的像元周边7x7窗口内，与0值像元地表覆盖类型相同像元，若没有则选取全部像元；b.利用Cressman客观分析法中的权重计算方法计算选取像元的权重；c.计算选取像元的实际值，即高质量像元观测值与背景值之差，作为波动值；d.对窗口内选取像元的波动值与其对应权重进行加权求和，获取最终LST插值结果；/n三、农作物种植区归一化植被指数-陆地表面温度特征空间的构建：/n提取干旱监测区域耕地，利用多年同期数据共同构建农作物各生长期农作物种植区归一化植被指数-陆地表面温度NDVI-LST散点图，并拟合各期干湿边方程，构建NDVI-LST特征空间；/n四、农作物温度植被干旱指数的计算：/n基于步骤三的结果，利用Price于1990年提出了温度植被干旱指数计算监测区域农作物种植区农作物温度植被干旱指数C-TVDI；/n五、基于农作物温度植被干旱指数的旱灾等级监测：/n基于有监督分类器的设计，通过历史数据确定监测区域基于农作物温度植被干旱指数的旱灾等级监测模型的参数，进行旱灾的监测；/n所述步骤五基于农作物温度植被干旱指数的旱灾等级监测，包括基于有监督分类器的设计，通过历史数据确定监测区域基于农作物温度植被干旱指数的旱灾等级监测模型的参数，进行旱灾的监测；具体如下：/n基于有监督分类器的设计，以往年作物生长期的遥感数据与各年实测旱灾等级数据作为训练样本，通过对训练样本的学习，确定旱灾等级监测模型中的参数；之后再利用最新一年的作物生长期遥感数据，对该年年终旱灾等级进行监测；该模型主要由受灾面积估算函数、参数寻优和旱灾等级监测三部分组成：/n(1)受灾面积估算函数：受灾面积定义为遭受各种自然灾害的农作物播种面积，且在同年遭受几种或几次自然灾害的，以其中危害最大的一次计算受灾面积，不重复计灾，旱灾受灾面积则为该年一次或多次遭受旱灾影响的最大农作物播种面积；考虑受灾面积值为年内旱灾对农作物的最大影响面积，考虑监测区如果为多次复种的种植模式，且农作物经收割之后一次旱灾的影响结束，因此，将受灾面积估算函数分为不同种植期，设定如下：/nS_i＝Max[Mean(s_ik，...，s_il)，...，Mean(s_im，...，s_in)]式1-8/n式中，S_i为第i年受灾面积估算值；k～l为每年第一期收割作物关键生长期，m～n为每年最后一期收割作物关键生长期，s_ij(j＝k～l)为第i年第一期收割作物各关键生长期的受旱面积，s_ij(j＝m～n)为最后一期收割作物各生长期的受旱面积；由于从受旱到旱灾的发生需要一定时间的积累，单期数据的旱情未必引发旱灾，而短期内干旱的空间变异性不大，因此将各期受旱面积的平均值作为受灾面积；全年的受灾面积值则为各期受灾面积的最大值；/n受旱面积s_ij(j＝k～l或m～n)由第i年农作物关键生长期农业遥感指数数据和该生长期阈值参数t_j(j＝k～l或m～n)决定；以农作物温度植被干旱指数为例，指数值越接近1则表示越干旱，则第i年第j个关键生长期的受旱面积S_ij为该年该期数据中C-TVDI值大于该期阈值t_j的像元所占面积，阈值t_j为模型参数，参数与所用生长期一一对应；/n(2)参数寻优，即利用寻优算法对设定的目标函数求取最值的过程；旱灾等级监测模型以多年遥感数据与农业受灾面积数据作为训练样本进行模型的参数寻优；模型参数的最优解可认为是在该组参数下，计算得到多年受灾面积估算值与多年训练样本中实际受灾面积的最小平均偏差，因此将目标函数设定为：/n /n式中，S_i为第i年受灾面积估算值，S_i0为第i年实际受灾面积值，n为训练样本个数，对目标函数求得最小值时对应的模型参数t_j(j＝k～l或m～n)即为最优参数；/n选用网格搜索算法进行参数寻优，网格搜索算法的寻优速度较快，获得全局最优解，不会陷入局部最优解；但在搜索前要设定参数范围，在参数范围较大，搜索步长较小的搜索中，训练耗时较长；/n参数寻优过程表述如下：(i)将训练样本中的遥感数据进行农作物遥感指数提取，并选择参与模型计算的作物生长期；(ii)设定最小目标函数值、各生长期参数t_i0的搜索范围、搜索步长和参数t_i0的起始值；(iii)将提取的农业遥感指数数据输入到受灾面积估算函数中，获得受灾面积估算值；(iv)将受灾面积估算值与训练样本中受灾面积数据带入目标函数中，计算目标函数值，若该值小于当前最小目标函数值，则保留当前参数，并将最小值更新；(v)通过网格搜索算法获得新的参数值，并重复上述(iii)、(iv)步骤，若全部参数均已遍历则完成寻优，记录最优参数；/n(3)旱灾等级监测，即通过待监测年农作物生长期遥感数据，对该年旱灾等级进行监测；由于参与训练的样本数较少，训练样本对不同程度灾情的代表性有限，因此利用训练样本受灾面积估算值与实际受灾面积间的线性关系，对待监测年的受灾面积估算值进行调整；旱灾等级监测的具体步骤如下：(i)计算训练样本中实际受灾面积与最优参数下受灾面积估算值的线性回归方程；(ii)提取待监测年农作物关键生长期遥感指数；(iii)利用受灾面积估算函数求取最优参数下的受灾面积估算值；(iv)将上一步中估算值作为线性方程自变量，获得调整后受灾面积估算值；(v)利用调整后受灾面积估算值计算受灾率I，并参照旱灾等级标准划分表，获得旱灾等级监测结果；旱灾等级标准划分如下：/n无旱：I≤5％；轻度干旱：5％＜I≤10％；中度干旱：10％＜I≤20％；严重干旱：20％＜I≤30％；特大干旱：I≥30％。/n