CN113470175B - 基于光学梯形模型的灌溉区域制图方法 - Google Patents

Abstract

本发明公开了一种基于光学梯形模型的灌溉区域制图方法，该方法中将待绘制区域划分为多个网格，解算每个网格中的参数，即认为网格中的每个像素点的参数一致，在保证准确的情况下提高了计算效率；利用所述参数，解算获得每个像素点在一段时间内的土壤含水量变化情况，进而与该像素点处的降水量相比较，设定合理的阈值，即可判断出该像素点是否为灌溉区域。

Description

基于光学梯形模型的灌溉区域制图方法

技术领域

本发明涉及一种灌溉区域测算技术领域，具体涉及一种基于光学梯形模型的灌溉区域制图方法。

背景技术

农耕文明是人类历史中不可或缺的一部分。伴随着耕作技术的诞生和发展，人类从条件适宜的大河流域走向世界的每个角落，这个过程离不开灌溉农业；特别是在降水量远小于蒸发量的干旱-半干旱地区，灌溉给农作物的发芽、生长和成熟提供了必要的湿度条件，很大程度上使得自然条件并不适宜耕种的土地能够产出足够多的粮食。

为了实现粮食安全和土地利用可持续发展的目标，需要使用大范围高分辨率的实际灌溉数据作为决策资料。在过去的半个世纪，全球灌溉农田面积翻了一倍，占全部耕地面积18％的灌溉农业区生产了全球粮食总产量40%，每年消耗的淡水量达到了人类淡水使用总量的70％。一方面，在人口增长和全球变化的双重影响下，部分地区的灌溉情况发生了变化，过去生产的一些数据在时间和空间尺度上的应用价值有所降低。另一方面，相比于传统的“灌区”概念，实际发生灌溉的区域更加符合水资源管理和农业管理的需要。

现有技术中，灌溉区域遥感制图方法主要分为三类：一类是完全使用统计数据的GIS方法，数据缺失多、更新慢、共享差；第二类使用以MODIS为代表的中等分辨率遥感数据，通过监测地表光谱特征和时间变化，最终获得空间分辨率较低的全球灌区产品；另一类使用Landsat等高分辨率遥感数据在流域范围内以监督分类的方式进行灌区提取，严重依赖于当地的统计资料，因而不能推广到统计资料缺乏的地区。目前来说全球灌区数据多为低分辨率，尚没有一套完整的全球30m分辨率的真实农田灌溉数据。

GIAM低分辨率全球灌区数据的分辨率低，且较为陈旧。M18中分辨率灌区数据中的灌区提取的方法利用了NDVI表征的作物生长周期，并和理论周期数进行比较推断灌溉区域，但是这种方法并没有考虑到不同作物本身需水量的差异，仍旧是间接的手段获取灌溉区域，无法直接用于分析区域灌溉强度，对于深入的灌溉研究难以提供有效信息。

LANID美国灌溉农业数据集获得的数据集为区域数据，使用了大量的地区统计数据；在制作过程中需要一定程度的人工干预和样本选择，在全球范围推广方面存在较大不确定性；

多源数据联合制图中，不同年份的数据训练的分类器在其他年的数据上的迁移效果差；获取的训练数据质量限制了方法的进一步应用。

利用遥感技术确定灌溉区域是大范围、高效率灌区制图的主要方法。为了得到多时相、大范围的实际灌溉区域分布图，实现完全基于高分辨率遥感影像的动态化、自动化制图，本发明人对遥感图像分析方法做了深入研究，以期待提出一种能够解决上述问题的基于高分辨率光学遥感影像的实际灌溉区域制图方法，通过采用光学梯形模型与时间序列分析的变化检测，来实现实际灌溉区域的自动制图。

发明内容

为了克服上述问题，本发明人进行了锐意研究，设计出一种基于光学梯形模型的灌溉区域制图方法，该方法中将待绘制区域划分为多个网格，解算每个网格中的参数，即认为网格中的每个像素点的参数一致，在保证准确的情况下提高了计算效率；利用所述参数，解算获得每个像素点在一段时间内的土壤含水量变化情况，进而与该像素点处的降水量相比较，设定合理的阈值，即可判断出该像素点是否为灌溉区域，从而完成本发明。

具体来说，本发明的目的在于提供一种基于光学梯形模型的灌溉区域制图方法，该方法包括如下步骤：

步骤1，将待绘制区域划分为多个网格，分别获得每个网格中的参数；

步骤2，针对遥感图像中的每个像素点，获得其在一段时间内的土壤含水量W；

步骤3，针对遥感图像中的每个像素点，比较其相邻两个时刻的土壤含水量变化量，并与对应的降水量做比较，从而判断该像素点是否为灌溉区域。

其中，在步骤1中，所述参数包括该网格中的干边截距i_d、干边斜率s_d、湿边截距i_w、湿边斜率s_w、局地最大湿土壤含水量W_w和局地最小干土壤含水量W_d。

其中，在步骤1中，获得干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w的过程包括如下子步骤：

子步骤1，选取网格对应的剔除不含云、云影、水体、冰雪像素点之后的卫星图像，据此利用红光波段、近红外波段、短波红外波段的反射率观测值获得该网格的归一化植被指数NDVI和短波红外转换反射率STR，，针对每个网格都获得对应的NDVI和STR时间序列影像；

子步骤2，在0.3到0.7范围内将特征空间的NDVI轴均匀分为宽度为0.01的40个区间，选取取值点，得到特征空间的干边和湿边；

子步骤3，对子步骤2得到的干边和湿边，分别用最小二乘法做线性回归拟合，得到干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w。

其中，所述STR通过下式(一)获得：

所述NDVI通过下式(二)获得：

其中，R_SWIR表示短波红外波段的地表反射率，R_NIR表示近红外波段地表反射率，R_red表示红光波段地表反射率。

其中，在步骤1中，获得最大湿土壤含水量W_w和局地最小干土壤含水量W_d的过程为：

调取一段时间内主被动土壤湿度任务SMAP提供的表层土壤含水量信息，对于每个网格，挑选出其在该时间段内所观测到的最大值和最小值，其中最大值即为所述局地最大湿土壤含水量W_w，所述最小值即为局地最小干土壤含水量W_d。

其中，在步骤2中，通过下式(三)获得土壤含水量W：

W＝W_d+θ(W_w-W_d) (三)

其中，θ表示归一化土壤湿度；

优选地，所述归一化土壤湿度θ通过下式(四)获得：

其中，在步骤3中，通过下式(五)获得土壤湿度变化量：

Δ＝W_n-W_n-1-P_(n-1，n) (五)

其中，Δ表示扣除降水的土壤含水量变化量，

W_n表示n时刻的像素点土壤含水量，W_n-1表示n-1时刻的像素点土壤含水量，P_(n-1，n)表示n-1时刻到n时刻的降水量。

其中，当Δ＞5mm时，该像素点为灌溉区域。

本发明所具有的有益效果包括：

(1)根据本发明提供的基于光学梯形模型的灌溉区域制图方法，从灌溉最直接的效果出发，即引起局部土壤含水量的剧烈变化，将高分辨率光学遥感影像反演土壤湿度的手段应用于灌溉区域的土壤含水量变化检测，结合降水数据对灌溉区域的灌溉事件的空间分布进行定量化分析，在大数据分析平台的支撑下，实现对实际灌溉区域的提取分析；

(2)根据本发明提供的基于光学梯形模型的灌溉区域制图方法，是一种完全基于遥感数据的方法，所有参数都是通过遥感观测来获得；相比于大量使用地面观测资料的现有方法，本方法的效率更高，可靠性更强；

(3)根据本发明提供的基于光学梯形模型的灌溉区域制图方法，该方法中选用光学遥感的方式获得观测资料，相对于使用热红外观测资料的方法，其分辨率更高，最终绘制的灌溉区域图的精度更高；

(4)根据本发明提供的基于光学梯形模型的灌溉区域制图方法，该方法中检测土壤含水量变化，因此降低了土壤湿度反演模型精度的影响，并且对于不同分辨率和不同区域的遥感图像有很好的适应性：只要卫星能在红光波段、近红外波段、短波红外波段对地表成像，即可使用改方法在卫星工作期间提取实际灌溉区域。

附图说明

图1示出根据本发明一种优选实施方式的基于光学梯形模型的灌溉区域制图方法的整体逻辑过程示意图；

图2示出根据本发明一种优选实施方式的基于光学梯形模型的灌溉区域制图方法中光学梯形模型原理示意图；

图3示出本发明实施例中得到的黑河流域的灌溉区域图。

具体实施方式

下面通过优选实施方式和实施例对本发明进一步详细说明。通过这些说明，本发明的特点和优点将变得更为清楚明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。

根据本发明提供基于光学梯形模型的灌溉区域制图方法，如图1中所示，该方法包括如下步骤：

步骤1，将待绘制区域划分为多个网格，分别获得每个网格中的参数；

优选地，将待绘制区域划分为边长为0.1°的正方形网格，即该网格的长为0.1纬度，宽为0.1经度，都约为11公里。卫星遥感图像的精度为30米，即其像素点一般对应于边长为30米的正方形区域，一个网格中最多能够包含十多万个像素点。

优选地，在步骤1中，所述参数包括该网格中的干边截距i_d、干边斜率s_d、湿边截距i_w、湿边斜率s_w、局地最大湿土壤含水量W_w和局地最小干土壤含水量W_d。

本申请中，所述干边和湿边，都是光学梯形模型中的参数，所述光学梯形模型的原理如图2中所示，其特征空间一般采用：横轴为归一化植被指数NDVI，纵轴为短波红外转换反射率STR的形式。同一地面目标在不同土壤湿度条件下，其(STR，NDVI)数据对的分布会构成一个梯形；土壤含水量与植被含水量线性相关，即土壤含水量越大，同等植被覆盖条件下植被含水量越高，STR越大。NDVI和STR线性相关，NDVI越大，同等土壤含水量条件下的STR越大。将最小干土壤含水量对应的STR-NDVI分布的下包络称为干边，将最大湿土壤含水量对应的STR-NDVI分布的上包络称为湿边。

优选地，在步骤1中，获得干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w的过程包括如下子步骤：

子步骤1，选取网格对应的剔除不含云、云影、水体、冰雪像素点之后的卫星图像，据此利用红光波段、近红外波段、短波红外波段的反射率观测值获得该网格的归一化植被指数NDVI和短波红外转换反射率STR，，针对每个网格都获得对应的NDVI和STR时间序列影像；

子步骤2，在0.3到0.7范围内将特征空间的NDVI轴均匀分为宽度为0.01的40个区间，选取取值点，得到特征空间的干边和湿边；具体来说，对于任意一个小区域的任意一张影像，裁切以后将所有像素点按照这些NDVI区间分为40组，仅保留每组中具有最大STR值和最小STR值的点；对于任意一个小区域，合并从多张影像得到的点，仅保留每个NDVI区间范围内中具有最大STR值和最小STR值的点；每个NDVI区间上具有最小STR值的点构成了特征空间的干边，每个NDVI区间上具有最大STR值的点构成了特征空间的湿边。

子步骤3，对子步骤2得到的干边和湿边，分别用最小二乘法做线性回归拟合，得到干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w。

在一个优选的实施方式中，所述STR通过下式(一)获得：

所述NDVI通过下式(二)获得：

其中，R_SWIR表示短波红外波段的地表反射率，R_NIR表示近红外波段地表反射率，R_red表示红光波段地表反射率，都能够通过卫星遥感图像中直接获得。

在一个优选的实施方式中，在步骤1中，获得最大湿土壤含水量W_w和局地最小干土壤含水量W_d的过程为：

调取一段时间内主被动土壤湿度任务SMAP提供的表层土壤含水量信息，对于每个网格，挑选出其在该时间段内所观测到的最大值和最小值，其中最大值即为所述局地最大湿土壤含水量W_w，所述最小值即为局地最小干土壤含水量W_d。W_w和W_d具有时间不变性。所述一段时间优选为5～6年，可以任选主被动土壤湿度任务SMAP运行期间连续的5～6年，如2015年～2020年。

本申请中，所述W_w和W_d可以将光学梯形模型估算的归一化土壤湿度转换为土壤含水量，进而满足后续步骤的要求。在常规做法中，一般将W_d定义为0、将W_w定义为土壤饱和含水量，但实际上灌溉区域土壤很少达到这些极端情况。本申请中采用的上述方法建立了一种从光学卫星观测值到土壤湿度卫星观测值的归一化土壤湿度转换方式，弥补了常规做法的不足，提高了土壤含水量估算的精度。

步骤2，针对遥感图像中的每个像素点，获得其在一段时间内的土壤含水量W；在解算任意一个像素点土壤含水量的过程中涉及到的参数，都选用该像素点所在网格的参数，即同一个网格中的任意像素点的参数都是彼此相同的。

在步骤2中，通过下式(三)获得土壤含水量W：

W＝W_d+θ(W_w-W_d) (三)

其中，θ表示归一化土壤湿度；

优选地，所述归一化土壤湿度θ通过下式(四)获得：

由于针对同一个像素点位置，有多帧不同时刻拍摄获得的卫星图像，在去除云、云影、水体、冰雪遮盖的情况后，还会有至少两帧卫星图像，此处需要分别针对每一帧卫星图像中的像素点分别予以计算，获得不同时刻的同一像素点处的土壤含水量。

步骤3，针对遥感图像中的每个像素点，比较其相邻两个时刻的土壤含水量差，并与对应的降水量做比较，从而判断该像素点是否为灌溉区域。

具体来说，在步骤3中，通过下式(五)获得一个像素点中扣除降水的土壤含水量变化量：

Δ＝W_n-W_n-1-P_(n-1，n) (五)

其中，Δ表示扣除降水的土壤含水量变化量，

W_n表示n时刻的像素点土壤含水量，W_n-1表示n-1时刻的像素点土壤含水量，P_(n-1，n)表示n-1时刻到n时刻的降水量，该值通过CHIRPS逐日降水数据从n-1时刻到n时刻逐日累积得到。所述具体时刻n的取值需要根据该像素点经筛选剩余的图像的卫星过境时间来确定，Δ＞Tmm能够满足n的一个取值就能认为该像素点发生一次灌溉事件，该像素点区域属于灌溉区域；当Δ＞Tmm能满足多个n的取值，就能认为该像素点发生多次灌溉事件，该像素点区域属于高频灌溉区域。所述阈值T表示扣除降水的土壤含水量变化量阈值，是一个经验的非负实数。

优选地，阈值T取5，即当Δ＞5mm时，该像素点在n-1到n时刻之间发生一次灌溉事件，该像素点为灌溉区域；

待绘制区域中的每一个像素点都经过上述方式，即都得到对应的Δ，统计每个像素点满足Δ＞5mm的时刻n的取值的个数，将个数大于0的像素点导出储存，即为灌溉区域图。

本申请中，将该临界值设置为5mm，既能够排除偶然因素导致的数据波动等意外因素，由于灌溉时的水量都远高于5mm，设置该临界值也不会将真实的灌溉区域排除在外。

实施例

针对黑河流域的灌溉区域，采用基于光学梯形模型的灌溉区域制图方法进行灌溉区域图绘制；

步骤1，调取黑河流域从2014年至2019年间来自Landsat的卫星图片，在其上按照横竖各0.1°划分出多个网格，解算每个网格中对应的参数，

所述解算参数的过程包括：子步骤1，选取网格对应的剔除不含云、云影、水体、冰雪像素点之后的卫星图像，据此利用红光波段、近红外波段、短波红外波段的反射率观测值获得该网格的归一化植被指数NDVI和短波红外转换反射率STR，针对每个网格都获得对应的NDVI和STR时间序列影像；

所述STR通过下式(一)获得：

所述NDVI通过下式(二)获得：

子步骤2，在0.3到0.7范围内将特征空间的NDVI轴均匀分为宽度为0.01的40个区间，选取取值点，得到特征空间的干边和湿边；

子步骤3，分别用最小二乘法做线性回归拟合，得到干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w。

再调取2015年至2020年内SMAP表层土壤含水量信息，对于每个网格，挑选出其在该时间段内所观测到的最大值和最小值，其中最大值即为局地最大湿土壤含水量W_w，所述最小值即为局地最小干土壤含水量W_d。

步骤2，通过下式(三)获得土壤含水量W：

W＝W_d+θ(W_w-W_d) (三)

其中，θ表示归一化土壤湿度，通过式(四)获得；

步骤3，通过下式(五)获得扣除降水的土壤含水量变化量，统计每个像素点扣除降水的土壤含水量变化量大于5mm的次数，将次数大于0的像素点作为灌溉区域并存储输出：

Δ＝W_n-W_n-1-P_(n-1，n) (五)

其中，Δ表示扣除降水的土壤含水量变化量，W_n表示n时刻的像素点土壤含水量，W_n-1表示n-1时刻的像素点土壤含水量，P_(n-1，n)表示n-1时刻到n时刻的降水量。

得到的黑河流域的灌溉区域图如图3中所示，该图与其他来源即低分辨率灌溉区域分布图的黑河流域灌溉区域面积统计结果一致，但是该图提供了更多的空间细节，并反映了不同地区灌溉频率的差异。在地面土壤湿度观测站点处使用步骤3所述的过程，得到真实灌溉事件的时间，与卫星得到的灌溉事件，即扣除降水的土壤含水量变化量大于5mm的时间段相比较，在黑河流域地表过程综合观测网的大满超级站在2014年至2020年的长期观测和SoilNET土壤湿度观测网络在2012年夏季的大范围观测中，发现所有卫星得到的灌溉事件的时间段内都发生了地面站点提供的真实灌溉事件，说明本申请方案得到的灌溉事件在验证站点上准确度为100％，基于此得到的灌溉区域地图具有很高的可靠性。

以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。

Claims (2)

1.一种基于光学梯形模型的灌溉区域制图方法，其特征在于，该方法包括如下步骤：

步骤1，将待绘制区域划分为多个网格，分别获得每个网格中的参数，将待绘制区域划分为边长为0.1°的正方形网格，即该网格的长为0.1纬度，宽为0.1经度，所述参数包括该网格中的干边截距i_d、干边斜率s_d、湿边截距i_w、湿边斜率s_w、局地最大湿土壤含水量W_w和局地最小干土壤含水量W_d；

步骤2，针对遥感图像中的每个像素点，获得其在一段时间内的土壤含水量W；

步骤3，针对遥感图像中的每个像素点，比较其相邻两个时刻的土壤含水量差，并与对应的降水量做比较，从而判断该像素点是否为灌溉区域；

在步骤1中，获得最大湿土壤含水量W_w和局地最小干土壤含水量W_d的过程为：

调取一段时间内主被动土壤湿度任务SMAP提供的表层土壤湿度信息，对于每个网格，挑选出其在该时间段内所观测到的最大值和最小值，其中最大值即为所述局地最大湿土壤含水量W_w，所述最小值即为局地最小干土壤含水量W_d；所述一段时间为主被动土壤湿度任务SMAP运行期间连续的5～6年；

在步骤2中，通过下式(三)获得土壤含水量W：

W＝W_d+θ(W_w-W_d) (三)

其中，θ表示归一化土壤湿度；

所述归一化土壤湿度θ通过下式(四)获得：

R_SWIR表示短波红外波段的地表反射率，R_NIR表示近红外波段地表反射率，R_red表示红光波段地表反射率；

在步骤3中，通过下式(五)获得扣除降水的土壤含水量变化量：

Δ＝W_n-W_n-1-P_(n-1，n) (五)

其中，Δ表示扣除降水的土壤含水量变化量，

W_n表示n时刻的像素点土壤含水量，W_n-1表示n-1时刻的像素点土壤含水量，P_(n-1，n)表示n-1时刻到n时刻的降水量；

Δ＞Tmm能够满足n的一个取值就能认为该像素点发生一次灌溉事件，该像素点区域属于灌溉区域；当Δ＞Tmm能满足多个n的取值，就能认为该像素点发生多次灌溉事件；阈值T表示扣除降水的土壤含水量变化量阈值；

当Δ＞5mm时，该像素点为灌溉区域；

统计每个像素点满足Δ＞5mm的时刻n的取值的个数，将个数大于0的像素点导出储存，即为灌溉区域图；

在步骤1中，获得干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w的过程包括如下子步骤：

子步骤1，选取网格对应的剔除不含云、云影、水体、冰雪像素点之后的卫星图像，据此利用红光波段、近红外波段、短波红外波段的反射率观测值获得该网格的归一化植被指数NDVI和短波红外转换反射率STR，针对每个网格都获得对应的NDVI和STR时间序列影像；

子步骤2，在0.3到0.7范围内将特征空间的NDVI轴均匀分为宽度为0.01的40个区间，选取取值点，得到特征空间的干边和湿边；其中，对于任意一个小区域的任意一张影像，裁切以后将所有像素点按照这些NDVI区间分为40组，仅保留每组中具有最大STR值和最小STR值的点；对于任意一个小区域，合并从多张影像得到的点，仅保留每个NDVI区间范围内中具有最大STR值和最小STR值的点；每个NDVI区间上具有最小STR值的点构成了特征空间的干边，每个NDVI区间上具有最大STR值的点构成了特征空间的湿边；

子步骤3，对子步骤2得到的干边和湿边，分别用最小二乘法做线性回归拟合，得到干边截距i_d、干边斜率s_d、湿边截距i_w和湿边斜率s_w。

2.根据权利要求1所述的基于光学梯形模型的灌溉区域制图方法，其特征在于，

所述STR通过下式(一)获得：

所述NDVI通过下式(二)获得：

其中，R_SWIR表示短波红外波段的地表反射率，R_NIR表示近红外波段地表反射率，R_red表示红光波段地表反射率。

Images