一种基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法

技术领域

本发明属于灌溉水管理领域，特别涉及一种基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法。

背景技术

灌溉用水量是世界淡水资源最大的消耗去向，约占世界用水总量的75％，但其利用全球约17％的耕地生产出的粮食占全球粮食总产量的40％左右。而我国利用全球7％的耕地养活了22％的世界人口，其中主要的功劳在于灌溉。根据水利部20世纪80年代初对于全国灌溉耕地与雨养耕地的粮食产量调查分析，灌溉耕地的粮食产量一般比雨养耕地的粮食产量高出1-3倍，尤其是对于干旱以及半干旱地区效果更加明显。所以，灌溉对于提高粮食产量具有重要意义。而灌溉面积的监测是灌区管理的重要环节，是获取灌区灌溉信息的重要来源。准确掌握灌溉面积、空间分布及其变化信息对灌区作物结构调整、灌溉信息化和管理、水资源优化配置、干旱监测和应急减灾均具有重要意义。大范围灌溉面积遥感监测使我们能够及时了解灌区的灌溉情况和干旱情势，为管理部门决策提供数据支持。

目前，获取灌溉面积的传统方法主要是通过灌溉用户上报，灌区工作人员记录整理的方式或利用灌溉用水量除以每亩所需水量的经验值的方式获得。但这两种方法不仅需要耗费大量的人力、物力、时间，而且会因漏记错记等问题造成获得的灌溉面积相关信息不准确，同时也不能获得较为准确的灌溉面积的空间分布情况，所以不能为灌区在水资源管理方面提供更为科学、准确的数据。此外，灌溉面积最直观的监测指标就是土壤含水量。如果耕地的土壤含水量骤增，则说明土壤发生了水的补给，即灌溉或降水。传统获取土壤含水量的方法有很多，如烘干称重法、电阻法、中子仪法、TDR测量法等，但这些方法只能获取某一点的不同土层深度的土壤含水量，不具有代表性，不仅费时、费力，还不能获取较大范围的土壤含水量信息。

发明内容

为了克服现有技术的问题，本发明提出了一种基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法，以遥感方法为技术支撑，采用MODIS遥感数据源，结合研究区物候、气象、灌水量等资料，以地表温度方法为技术手段，获取灌区灌溉面积时空信息。可通过实时有效的动态监测，定量掌握灌溉用水的去向，对灌溉用水的合理分配、调度和管理，提供数据支撑。

本发明的目的是这样实现的：

一种基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法，包括以下步骤：

步骤1，研究区地表温度信息的获取：

根据研究区选定研究作物的生育期，选用对应研究作物灌溉期间的遥感卫星MODIS地表温度数据，利用MRT工具对MODIS地表温度数据进行第1通道与第5通道波段的提取，提取后的数据再经过包括重投影、镶嵌、重采样、裁剪在内的预处理；

所述MODIS地表温度数据进行第1通道、第5通道的数据分别对应一天中最高的地表温度(Tmax)和最低地表温度(Tmin)；

步骤2，降雨数据的剔除：

获取步骤1选定的研究作物灌溉期内的历史气象数据，利用所述历史气象数据对步骤1获取的数据进行剔除处理；

首先，剔除降水当天和降水后第1天的数据，降水后第2天至第5天的数据依据不同情况进行判定，是否剔除；

步骤3，地表温度日较差曲线图的绘制：

基于经上述步骤1、2处理后的数据，计算地表温度日较差，并绘制地表温度日较差曲线图；

所述地表温度日较差定义为一天中最高的地表温度(Tmax)和最低地表温度(Tmin)之间的差值，即MODIS第1通道数据与第5通道数据的差值；

步骤4，多角度结果分析：

分别从空间、时间不同角度，对步骤3绘制的地表温度日较差曲线图进行分析，具体的：

从空间角度，进行单景影像空间分布分析：从不同像元间的空间温度差异判断像元是否发生灌溉；

从时间角度，进行一段时间多景影像纵向分析，根据发生灌溉后地表温度日较差曲线是否出现“先降低、后升高，最后趋于稳定”的变化趋势，判断发生灌溉的时间段；

步骤5，利用阈值率定判断是否发生灌溉：

首先，根据研究区选定研究作物的灌溉期气温变化特征和研究作物的物候特征(植被生长阶段、覆盖度变化等)，将灌溉期划分为若干时间段；

其次，结合步骤4获得的空间、时间趋势分析结果，利用步骤2所述历史气象数据，共同率定阈值△Ts，从而判断是否发生灌溉，并通过地表温度日较差变化曲线图率定不同时段的阈值，当地表温度日较差小于△Ts时则认为发生了灌溉；

所述阈值为地表温度日较差曲线中经历先降低后升高的最低值。

进一步的，步骤1中所述研究作物为水稻，选用每年5月～9月水稻灌溉期间的遥感卫星MODIS地表温度数据。

进一步的，步骤2中所述历史气象数据包括：日最高气温、日最低气温、气温日较差，从中国气象局数据网获取。

进一步的，步骤2中采用决策树方法进行降雨数据剔除，对降水后第2天至第5天的数据处理方法具体为：

若降水后第1天没有显示灌溉，而降水后第2天或以后几天显示灌溉，则说明未受降水影响，不需要进行剔除；

若降水后第1天显示灌溉，降水后第2天至第5天内显示灌溉且地表温度日较差逐渐变大，则说明受降水影响，需要进行剔除。

本发明的优点和有益效果是：

1、本发明所述基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法，采用MODIS遥感数据源，结合研究区物候、气象、灌水量等资料，以地表温度方法为技术手段，获取灌区灌溉面积时空信息。可通过实时有效的动态监测，定量掌握灌溉用水的去向，对灌溉用水的合理分配、调度和管理，提供数据支撑；

2、本发明基于MODIS卫星遥感数据，利用地表温度日较差等技术方法，建立基于遥感技术的多尺度的实际灌溉面积获取技术方案，制作研究区不同空间尺度的实际灌溉面积数据产品。

附图说明

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

图1是典型地物表面温度日变化曲线图；

图2是实施例步骤2剔除受降水影响数据的决策树分类图；

图3是实施例监测期间某日的地表温度日较差曲线图；

图4～8分别是2016-2020年灌溉强度和灌溉面积空间分布图；

图9是实施例精度验证部分样本点分布图；

图10是实施例精度验证部分IWMI灌溉面积提取结果样本点分布图；

图11是实施例精度验证部分基于地表温度日较差方法的提取结果与国外数据对比结果图。

具体实施方式

本申请所述提取阈值确定方法的理论依据为：

受太阳辐射、地物组分、物理状态、热特性、几何结构、生态环境、土壤物理参数等因素影响，不同地物地表温度及其日变化存在显著差异。图1示意了典型地物(沙地、草地、林地、湖泊)的相对辐射温度变化曲线。水体比热容最大，即放出/吸收相同的热量，其温度降低/升高的少，因此水体白天升温慢，夜间降温慢，昼夜温度变化比周围地物慢而小，其日温差比常见地物都小。对于裸土，干燥裸土的比热容远小于水体，在白天吸收热量迅速升温，晚上释放热量迅速降温，因此昼夜温度变化较大，其日温差远大于水体。因此，裸土中的水分能够增大裸土比热容，日温差比干燥裸土小，且裸土含水量越大，裸土日温差越小。对于有植被覆盖度的区域如草地、林地、耕地等，卫星遥感反演的地表温度受植被与土壤综合影响。土壤含水量影响植被蒸腾量，而蒸腾量又影响植被表层温度。白天，作物缺水时土壤水不能满足作物蒸腾需要，蒸腾吸热减少，叶面温度升高，导致其表层温度高于不缺水作物；夜间，由于水分的保温作用使缺水作物的表层温度比不缺水作物降温快；这使得缺水作物的地表温度日较差大于不缺水作物。基于上述原理，灌溉后耕地的土壤含水量显著增加，其地表温度日较差明显小于未灌溉耕地，因此可以通过确定一个地表温度日较差的阈值(△Ts)判断该耕地是否发生灌溉，即当地表温度日较差小于△Ts时则认为发生了灌溉。

实施例：

本实施例以水稻作为研究作物，对船行灌区进行水稻灌溉面积监测，提供了一种基于MODIS地表温度日较差的灌区灌溉面积提取阈值确定方法，包括以下步骤：

步骤1，研究区地表温度信息的获取：

根据水稻的生育期，购买2016年5月～9月、2017年5月～9月、2018年5月～9月、2019年5月～9月和2020年5月～9月水稻灌溉期间的遥感卫星MODIS地表温度数据产品——MOD11A1，空间分辨率1km，时间分辨率1d，水平编号和垂直编号为h27v05。

利用MRT工具(MODIS Reprojection Tool)对获取的MODIS地表温度数据MOD11A1进行第1通道与第5通道波段的提取，提取后的数据再经过包括重投影、镶嵌、重采样、裁剪在内的预处理。

所述MODIS地表温度数据进行第1通道、第5通道的数据分别对应一天中最高的地表温度(Tmax)和最低地表温度(Tmin)。

步骤2，降雨数据的剔除：

从中国气象局数据网(http://data.cma.cn/)获取2016年5月～9月、2017年5月～9月、2018年5月～9月、2019年5月～9月和2020年5月～9月的历史气象数据，包括日最高气温、日最低气温、气温日较差，利用该历史气象数据对步骤1获取的数据进行剔除处理，从而剔除因降雨导致的伪灌溉影响。

如图2所示，采用决策树方法进行降雨数据剔除，具体方法是：

首先，剔除降水当天和降水后第1天的数据。

其次，降水后第2天至第5天的数据依据不同情况进行判定，是否剔除。若降水后第1天没有显示灌溉，而降水后第2天或以后几天显示灌溉，则说明未受降水影响，不需要进行剔除；若降水后第1天显示灌溉，降水后第2天至第5天内显示灌溉且地表温度日较差逐渐变大，则说明受降水影响，需要进行剔除。

步骤3，地表温度日较差曲线图的绘制：

基于经上述步骤1、2处理后的数据，计算地表温度日较差，并绘制地表温度日较差曲线图，如图3所示。

所述地表温度日较差定义为一天中最高的地表温度(Tmax)和最低地表温度(Tmin)之间的差值，即MODIS第1通道数据与第5通道数据的差值。

步骤4，多角度结果分析：

我国农业灌溉以户为单位进行，单户耕地面积较小，且分布较为分散，因此在影像中寻找“灌溉纯像元”以率定阈值较为困难。本实施例利用经上述方法处理后的MODIS地表温度数据，分别从空间、时间不同角度，对步骤3绘制的地表温度日较差曲线图进行分析。

具体的方法为：

从横向上空间角度，进行单景影像空间分布分析：从不同像元间的空间温度差异判断像元是否发生灌溉；每次发生灌溉后，地表温度会减小，因此灌溉处的地表温度小于未灌溉处。

从纵向时间角度，进行一段时间多景影像纵向分析，根据发生灌溉后地表温度日较差曲线是否出现“先降低、后升高，最后趋于稳定”的变化趋势，判断发生灌溉的时间段。

这里需要说明的是，发生灌溉后一段时间内会有一个地表温度日较差变化曲线(即先将低后升高最后趋于稳定)。

步骤5，利用阈值率定判断是否发生灌溉：

首先，根据研究区选定研究作物的灌溉期气温变化特征和研究作物的物候特征(植被生长阶段、覆盖度变化等)，将灌溉期划分为若干时间段。

其次，结合步骤4获得的空间、时间趋势分析结果，利用步骤2获得的研究区日最高气温、日最低气温、气温日较差等历史气象数据，共同率定阈值△Ts，从而判断是否发生灌溉，并通过地表温度日较差变化曲线图率定不同时段的阈值，当地表温度日较差小于△Ts时则认为发生了灌溉。

所述阈值为地表温度日较差曲线中经历先降低后升高的最低值。

通过对2015-2020年灌溉期气温及天气状况分析发现，每年6月-9月日最高气温、日最低气温呈先升高后降低趋势，6月以后最低气温一般在21℃以上，8月下旬以后气温开始下降，最低气温20℃以下。研究区主要作物为水稻，一般具有幼苗期、插秧期、分蘖期、拔节期、孕穗期、抽穗期、扬花授粉期、灌浆期八个生长发育阶段，自幼苗期开始植被覆盖率逐渐变大。分蘖期期前水稻植株较小，植被覆盖率低，裸土较多；而分蘖期及以后阶段，水稻植株迅速生长，叶片数量逐渐增多并孕育稻穗，植被覆盖率变大，裸土减少。

在本实施例中根据研究区灌期气温变化和水稻物候特征，将灌溉期间划分为四个时段，即5月20日-6月29日、6月30日-7月23日、7月24日-8月20日和8月21日-9月30日，确定灌溉期间4个时段识别灌溉的日较差阈值△Ts1、△Ts2、△Ts3、△Ts4分别为10℃，9℃，8℃，10℃。

结果分析：

根据率定的阈值，剔除降雨数据的影响，对2016-2020年船行灌区区域分布情况进行遥感监测判别。灌溉时间主要在6月-9月进行。由于每年降雨量等气象因素的影响，每年的灌溉次数会有一定差异，2016-2020年的数据情况相对较好，可用数据量符合灌溉面积监测要求，根据降水、天气和数据质量等影响，不同年份的灌溉情况有所差别，遥感识别的结果能够有效的反映出灌溉次数的时空分布特征。

如图4～8所示，2016-2020年灌溉强度和灌溉面积空间分布图，遥感识别结果能够有效地反映出灌溉次数的分布特征，五年的灌溉次数在4次和5次以上，2016-2018年平均灌溉次数在5-7次范围内，2019-2020年平均灌溉次数为6-8次，其中2016-2018年由于灌溉期间降雨量较多，所以2016-2018年总体灌溉次数少于2019和2020年的总体灌溉次数。2016和2017年以5次和6次灌溉为主，2018年主要以6次和7次灌溉为主，2019年和2020年主要以灌溉7次和8次为主，其中有些年份出现灌溉次数少的部分区域，可能该地的水稻种植时间早一些，遥感监测的部分区域要灌溉次数少一些。

在灌溉次数识别结果的基础上，进一步计算得到船行灌区2016-2020年累计灌溉面积，如表1所示，2016-2020年累计灌溉面积约增加了65.02万亩次，其中2016和2017年的累计灌溉面积要小于2018-2020年累计灌溉面积。这是由于2016和2017年降水要多于2018-2020年。

表1船行灌区灌溉面积和累计灌溉面积统计表

方法精度验证：

1、基于样本点的实际灌溉面积精度评估

2020年10月野外勘察获取的149个样本点，其中主要包括水稻、小麦等，对2020年基于地表温度日较差方法提取的灌溉面积进行精度验证，如图9所示。

选取研究区主要作物水稻进行验证，其中水稻包括单季水稻和稻麦轮作水稻，2020年水稻的样本点共49个，有48个样本点落在光谱匹配内部，有1个样本点落在光谱匹配识别的灌溉区域外部，识别总体精度为97.96％。

2、基于国内外其他相关产品的灌溉面积空间分布合理性与精度评估

2010年世界水资源管理研究所(IWMI)绘制了亚洲2010年220m空间分辨率的灌溉与非灌溉区域分布图，用来进行对比分析。如图10，此数据共分为背景值、灌溉区域、雨养区域等，船行灌区研究区主要包含灌溉区域与非灌溉区域。

野外勘察获取149个样本点(包括水稻样本点49个，小麦样本点47个，)中，其中水稻和小麦的样本点一共96个，其中一共64个样本点(包括水稻样本点32个，小麦样本点32个)落在灌溉区域内，识别精度为66.67％。由此看出，光谱匹配结果样本点落在灌溉区域的要大于IWMI产品，因此光谱匹配方法获取的船行灌区实际灌溉面积精度要高于IWMI绘制的产品精度。

将光谱匹配识别的船行灌区灌溉面积与IWMI-Irrigated Area Map Asia(2000-2010)and Africa(2010)数据产品和GRIPC(2018)数据产品进行对比分析。图11分别为光谱匹配与其他数据产品的灌溉面积图，由图11可知光谱匹配识别的实际灌溉面积空间分辨率明显高于IWMI与GRIPC数据产品。实际灌溉面积识别结果精度较高，其灌溉区域分布范围与IWMI数据产品的灌溉范围基本一致，非灌溉区域主要分布在灌区的北部和东部某区域，灌溉区域主要分布在灌区的西部、南部、中部和东部区域，而GRIPC数据产品的只在灌区的中部与东部小范围识别出非灌溉区域，其余均为灌溉区域，识别精度较低。另外由局部图可知，实际灌溉面积的边界精度明显优于其他数据产品。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。