基于作物水分指示线的农田灌溉概率指数计算方法

1.本发明涉及农业遥感监测领域，尤其涉及一种基于作物水分指示线的农田灌溉概率指数计算方法。


背景技术：

2.灌溉是人为利用自然降水以外的其他水源补给植被生长所需水分的措施，其能有效缓解气候变化和干旱气象对农业生产的不利影响；保证作物正常生长，实现农业的高产稳产(与雨养农业相比)。灌溉在改善降水量不足或分布不匀、保证作物正常生长、提高粮食产量和保障粮食安全中发挥了重要作用；同时灌溉也是淡水资源消耗的主要途径，全球灌溉农业用水量约占淡水抽取总量的70％(许多发展中国家的灌溉甚至消耗了超过90％的取水量)；更进一步，不合理的灌溉还会对地下水环境、地表水质、土壤、社会经济甚至人体健康造成不利的影响。在当前全球淡水资源紧缺背景下，科学掌握农业灌溉的分布、需水量、发生时段及其动态变化对于水资源管理、粮食安全保障和经济发展规划等具有重要意义。
3.灌溉面积及其空间分布是农田灌溉的基础信息。传统费时费力的野外调查方法已无法满足现代农业对大区域覆盖与及时更新的灌溉空间分布信息的需要。随着空间信息技术的飞速发展，具有覆盖范围广、信息时效高和作业成本低等优点的地理空间信息制图技术成为区域农业灌溉分布信息获取的主要手段。当前常用的耕地灌溉分布信息提取方法主要包括空间分配技术与遥感分类技术两大类。其中空间分配技术是借助一定的空间分配规则，将将行政单元尺度上统计灌溉面积空间分配至特定的网格单元，从而生成灌溉耕地分布图；核心在于找出与灌溉高度相关的某个或某些变量(如灌溉设施分布、作物复种指数等)，并建立合理的空间分配模型。这种方法是目前全球大尺度灌溉耕地制图的主要方法，代表产品有全球灌溉耕地数字地图gmia(global map of irrigation areas)和rica(monthly irrigated and rainfed crop areas)。基于遥感分类技术获取得到灌溉耕地的分布，往往选择多个能区分灌溉耕地与雨养耕地的特征变量，比如波段指标(光谱反射率、植被指数、水体指数等)、物候指标(生长季起止时间、峰值等)、气候和环境变量(地表温度。蒸散发等)等，采用监督或者非监督的分类方法来提取灌溉耕地。随着遥感成像与解译技术的发展，面向灌溉耕地提取的遥感数据源与分类特征从单一发展到多源、技术方法从人工解译到自动分类乃至多分类器集成、制图尺度也从局地典型农业区扩展至流域、平原、国家甚至全球。
4.在上述两种方法实施中，选择合适有效的特征变量是重要步骤。目前用来表征耕地受灌溉可能性大小的特征变量主要有两类：第一类参量描述灌溉引起的土壤水分或作物水分的变化(如作物水分指数)；第二类参量关注灌溉引起的植被长势的变化(如归一化植被指数峰值等)；也有研究结合了植被因素和非植被因素来表征耕地的灌溉潜力。然而大部分的特征参量对于灌溉可能性的表征都比较间接，物理机理不够明确、灌溉指示性较弱。
5.在灌溉耕地分布的遥感提取中，结合多源遥感数据与丰富光谱特征，借助机器学习方法(诸如支持向量机、随机森林、甚至深度学习等)进行灌溉制图是当前的主流方法。然
而这类方法需要大量的实地样本或真值，样本缺失给区域灌溉耕地制图带来较大挑战；且基于局部区域实验所构建的分类或提取规则普适性不高，限制了方法的区域推广。
6.最近有学者尝试从雨养耕地的气象干旱与农业干旱的高度相关出发，寻找雨养指示线、并构建耕地灌溉概率指数。这给农业灌溉遥感监测提供了新的思路，然而其存在的不足主要表现为：(1)对气象干旱的表达过于简化，即未考虑作物缺水的累积效应，也未考虑作物生长环境的蒸散发条件；(2)雨养指示线仅考虑了上边线所表达的气象干旱与农业干旱的相关性，而未能考虑下边线对气象干旱与农业干旱的边界极值约束。整体上使得对干旱的估算偏高、产生过饱和现象，进而影响对作物灌溉状况的估算精度和模型的适用性。
7.综上所述，在目前可见的专利和文献中，基于遥感的农田灌溉信息提取方法较少；在遥感灌溉制图相关的特征选择、计算模型构建、算法推广应用等方面尚存在的问题较多；缺乏物理机制明确、指标指示性强的灌溉指数计算方法。


技术实现要素：

8.发明目的：本发明的目的在于提供一种基于作物水分指示线的农田灌溉概率指数计算方法，准确、简便表达作物水分亏缺程度，提升农田灌溉概率指数的物理意义和指示意义。
9.技术方案：基于作物水分指示线的农田灌溉概率指数计算方法，包括以下步骤：
10.(1)确定灌溉监测的地理空间参考和时空分辨率，并预处理多源遥感数据集，包括时间序列降水、潜在蒸散发、实际蒸散发以及农作物分布图，生成具有相同地理空间参考和时空分辨率的、且像元对齐的栅格数据集；
11.(2)在作物分布区域内，基于像元的降水量、实际蒸散发量和潜在蒸散发量，计算其作物水分亏缺指数与作物水分胁迫指数，利用实际蒸散发量与潜在蒸散发量的加权和替代作物生长需水量，计算像元每时相的作物水分亏缺值，并以加权和的方式累计其前期多时相的作物水分亏缺，得到当期作物水分亏缺指数，利用实际蒸散发量与潜在蒸散发量的相对差异百分比，计算作物水分胁迫指数；
12.(3)构建以像元作物水分亏缺指数为横坐标、作物水分胁迫指数为纵坐标的作物水分散点图；
13.(4)在作物水分散点图中，假设横坐标分割区间内的纵坐标值符合正态分布，并基于其均值与标准差来确定纵坐标值的置信区间，同步提取散点图的上边界点集、下边界点集，直线拟合其上边界点集、生成作物水分的上指示线，直线拟合其下边界点集、生成作物水分的下指示线；
14.(5)在作物水分散点图中，分别计算农田区域像元到上、下指示线的竖直上距离、竖直下距离，构建以上距离与下距离为变量的像元灌溉概率指数，综合利用作物水分的上指示线、下指示线，构建农田像元灌溉概率指数计算模型；
15.(6)计算作物分布区域内的像元灌溉概率指数，生成区域农田灌溉概率指数分布图。
16.所述步骤(1)具体为，在确定农田灌溉遥感监测的地理空间参考r、空间分辨率s、时间分辨率t参数的基础上，首先对时间序列降水、潜在蒸散发量、实际蒸散发量进行t天像元值累计合成，生成时间分辨率为t的时间序列数据集；然后将合成时间序列数据集与农作
物分布图转换到地理空间参考r空间中；接着将合成时间序列数据重采样到空间分辨率s，并严格对齐像元空间位置；最后统计每个像元范围内的农作物面积占比，并将面积占比超过50％的像元判定为耕地像元，获取农田区域数据集。
17.所述步骤(2)中选取作物水分胁迫指数(crop water stress index,cwsi)来表征作物的农业干旱状态，选取作物水分亏缺指数(crop water deficit index,cwdi)来表征作物的气象干旱状态；作物水分亏缺指数cwdi计算方法如下式：cwdi＝a
×
cwdii+b
×
cwdi
i-1
+c
×
cwdi
i-2
+d
×
cwdi
i-3
+e
×
cwdi
i-4
18.其中cwdi为作物生育时段的累计水分亏缺指数；考虑作物干旱过程的累积效应，作物水分亏缺指数一般取连续5个时相：当前及其前4个时相，cwdii、cwdi
i-1
、cwdi
i-2
、cwdi
i-3
、cwdi
i-4
；作物亏缺的加权和；进一步考虑到累积效应随间隔时长的衰减，5个时相的累计权重系数分别为a、b、c、d、e，分别取值为0.3、0.25、0.2、0.15、0.1；
19.作物水分亏缺指数cwdi表征作物水分的亏缺程度，表达为作物需水量与实际供水量之差，以百分率表示；本发明简化作物水分亏缺的计算、并考虑作物生长环境的实际蒸散条件，用像元的实际蒸散发量与潜在蒸散发量的加权和代替参考作物蒸散量，其权重分别为0.5，简化的cwdi计算如下：
20.cwdii＝1-prei/(0.5
×
aeti+0.5
×
peti)
21.其中aet和pet分别表示实际蒸散发量与潜在蒸散发量，pre为同期降水量，作物水分胁迫指数cwsi的计算方法如下式：
22.cwsi＝1-aet/pet
23.作物水分亏缺指数和胁迫指数的数值越大，表示干旱越严重，数值越小，表示干旱情况越轻。
24.所述步骤(3)在作物水分散点图构建中，以作物水分亏缺指数cwdi为横坐标轴、作物水分胁迫指数cwsi为纵坐标轴，建立像元气象干旱与农业干旱的分布空间，并绘制农田像元散点图。
25.所述步骤(4)首先将横坐标气象干旱指数cwdi以0.01间隔进行区间划分；然后提取每个划分区间q中的散点，统计其农业干旱指数cwsi的均值、标准差、以及置信区间，为去除由噪声等因素导致的异常值，取95％的置信区间；接着查找置信区间内的cwsi最大值与最小值，并分别记其对应的散点为划分区间q的作物水分上边界点与下边界点；最后提取所有划分区间的上边界点构成散点图的上边界线，并以最小二乘的方式直线拟合生成作物水分上指示线；同样方法生成作物水分下指示线，作物水分的上指示线与下指示线分别记为fu和fb：
26.fu＝au×
cwdi+bu27.fb＝ab×
cwdi+bb28.其中au和bu为拟合得到的上指示线的斜率和截距，ab和bb为拟合得到的下指示线的斜率和截距。
29.所述步骤(5)在作物水分散点空间中，以作物水分上、下指示线为约束，计算每个农田像元与作物水分上、下指示线在纵轴方向上的距离irr、pir，像元坐标为(cwdi,cwsi)，其中irr表征灌溉对作物水分胁迫的缓解程度，irr与pir之和mi表征农业灌溉对作物水分胁迫的最大缓解程度，irr与mi的比值表征该农田像元受到灌溉的可能性，农田像元的灌溉
概率指数建模如下：
30.irr＝fu(cwdi)-cwsi
31.pir＝cwsi-fu(cwdi)或mi＝fu(cwdi)-fb(cwdi)
[0032][0033]
对于位于作物水分上指示线之上的像元，直接对其灌溉概率指数赋值为0；对于位于作物水分下指示线之下的像元，直接对其灌溉概率指数赋值为1.0。
[0034]
有益效果：本发明综合考虑作物生长环境的蒸散发条件与作物缺水的累积效应，利用实际蒸散发量与潜在蒸散发量的加权和来替代作物需水量，准确、简便表达作物水分亏缺程度；在基于作物水分散点图的作物水分指示线提取中，能够同步拟合作物水分的上指示线、下指示线，进一步约束作物水分(干旱指数)的散点空间，准确表达灌溉对作物水分的影响程度、提升农田灌溉概率指数的物理意义和指示意义。
附图说明
[0035]
图1是基于作物水分指示线的农田灌溉概率指数计算方法示意图
[0036]
图2是基于作物水分指示线的农田灌溉概率指数计算实现流程图
[0037]
图3是多源遥感数据预处理的流程图
[0038]
图4是像元作物水分指数计算与散点图建立的流程图
[0039]
图5是作物水分散点图中作物水分上下指示线提取的流程图
[0040]
图6是基于作物水分上下指示线的耕地灌溉概率指数建模流程图
[0041]
图7是某地农田灌溉概率指数空间分布示意图
具体实施方式
[0042]
下面结合附图对本发明的技术方案作进一步说明。
[0043]
如图1所示，在构建作物水分指数散点图步骤中，利用预处理后的多源遥感数据集计算农田像元的作物水分亏缺指数和作物水分胁迫指数，并构建农田区域的作物水分散点图；在像元灌溉概率指数建模步骤中，提取作物水分散点图的上边界点集与下边界点集，分别直线拟合出作物水分上指示线与下指示线，进而构建像元灌溉概率指数模型。
[0044]
如图2所示基于作物水分指示线的农田灌溉概率指数计算方法的实现流程图，包含有多源数据预处理、作物水分(亏缺与胁迫)指数计算和散点图建立、作物水分上下边界点提取、作物水分上下指示线拟合、像元灌溉概率指数建模、农田灌溉概率指数分布图生成等6个实现单元。
[0045]
如图1-2，本发明提供一种技术方案：基于作物水分指示线的农田灌溉概率指数计算方法，包括以下步骤：
[0046]
(1)确定灌溉监测的地理空间参考和时空分辨率，并预处理多源遥感数据集，包括时间序列降水、潜在蒸散发、实际蒸散发以及农作物分布图，生成具有相同地理空间参考和时空分辨率的、且像元对齐的栅格数据集；
[0047]
(2)在作物分布区域内，基于像元的降水量、实际蒸散发量和潜在蒸散发量，计算其作物水分亏缺指数与作物水分胁迫指数，利用实际蒸散发量与潜在蒸散发量的加权和替
代作物生长需水量，计算像元每时相的作物水分亏缺值，并以加权和的方式累计其前期多时相的作物水分亏缺，得到当期作物水分亏缺指数，利用实际蒸散发量与潜在蒸散发量的相对差异百分比，计算作物水分胁迫指数；
[0048]
(3)构建以像元作物水分亏缺指数为横坐标、作物水分胁迫指数为纵坐标的作物水分散点图；
[0049]
(4)在作物水分散点图中，假设横坐标分割区间内的纵坐标值符合正态分布，并基于其均值与标准差来确定纵坐标值的置信区间，同步提取散点图的上边界点集、下边界点集，直线拟合其上边界点集、生成作物水分的上指示线，直线拟合其下边界点集、生成作物水分的下指示线；
[0050]
(5)在作物水分散点图中，分别计算农田区域像元到上、下指示线的竖直上距离、竖直下距离，构建以上距离与下距离为变量的像元灌溉概率指数，综合利用作物水分的上指示线、下指示线，构建农田像元灌溉概率指数计算模型；
[0051]
(6)计算作物分布区域内的像元灌溉概率指数，生成区域农田灌溉概率指数分布图。
[0052]
如图3所示，在多源数据预处理中，本发明提出的农田灌溉概率指数模型涉及到时间序列降水、潜在蒸散发、实际蒸散发以及农作物分布图等多源栅格数据，其往往存在着地理空间参考和时空分辨率等不一致的问题，数据预处理将为后续建模提供具有相同地理空间参考和时空分辨率的、且像元对齐的栅格数据集。
[0053]
所述步骤(1)具体为，在确定农田灌溉遥感监测的地理空间参考r、空间分辨率s、时间分辨率t参数的基础上，首先对时间序列降水、潜在蒸散发量、实际蒸散发量进行t天像元值累计合成，生成时间分辨率为t的时间序列数据集；然后将合成时间序列数据集与农作物分布图转换到地理空间参考r空间中；接着将合成时间序列数据重采样到空间分辨率s，并严格对齐像元空间位置；最后统计每个像元范围内的农作物面积占比，并将面积占比超过50％的像元判定为耕地像元，获取农田区域数据集。
[0054]
如图4所示，在像元作物水分指数计算中，利用预处理后的时间序列降水量、实际蒸散量和潜在蒸散量数据来计算像元作物水分指数。
[0055]
所述步骤(2)中选取作物水分胁迫指数(crop water stress index,cwsi)来表征作物的农业干旱状态，选取作物水分亏缺指数(crop water deficit index,cwdi)来表征作物的气象干旱状态；作物水分亏缺指数cwdi计算方法如下式：
[0056]
cwdi＝a
×
cwdii+b
×
cwdi
i-1
+c
×
cwdi
i-2
+d
×
cwdi
i-3
+e
×
cwdi
i-4
[0057]
其中cwdi为作物生育时段的累计水分亏缺指数；考虑作物干旱过程的累积效应，作物水分亏缺指数一般取连续5个时相：当前及其前4个时相，cwdii、cwdi
i-1
、cwdi
i-2
、cwdi
i-3
、cwdi
i-4
；作物亏缺的加权和；进一步考虑到累积效应随间隔时长的衰减，5个时相的累计权重系数分别为a、b、c、d、e，分别取值为0.3、0.25、0.2、0.15、0.1；
[0058]
作物水分亏缺指数cwdi表征作物水分的亏缺程度，表达为作物需水量与实际供水量之差，以百分率表示；本发明简化作物水分亏缺的计算、并考虑作物生长环境的实际蒸散条件，用像元的实际蒸散发量与潜在蒸散发量的加权和代替参考作物蒸散量，其权重分别为0.5，简化的cwdi计算如下：
[0059]
cwdii＝1-prei/(0.5
×
aeti+0.5
×
peti)
[0060]
其中aet和pet分别表示实际蒸散发量与潜在蒸散发量，pre为同期降水量，作物水分胁迫指数cwsi的计算方法如下式：
[0061]
cwsi＝1-aet/pet
[0062]
作物水分亏缺指数和胁迫指数的数值越大，表示干旱越严重，数值越小，表示干旱情况越轻。
[0063]
所述步骤(3)在作物水分散点图构建中，以作物水分亏缺指数cwdi为横坐标轴、作物水分胁迫指数cwsi为纵坐标轴，建立像元气象干旱与农业干旱的分布空间，并绘制农田像元散点图。
[0064]
如图5所示，在作物水分散点图的作物水分上下指示线提取中，所述步骤(4)首先将横坐标气象干旱指数cwdi以0.01间隔进行区间划分；然后提取每个划分区间q中的散点，统计其农业干旱指数cwsi的均值、标准差、以及置信区间，为去除由噪声等因素导致的异常值，取95％的置信区间；接着查找置信区间内的cwsi最大值与最小值，并分别记其对应的散点为划分区间q的作物水分上边界点与下边界点；最后提取所有划分区间的上边界点构成散点图的上边界线，并以最小二乘的方式直线拟合生成作物水分上指示线；同样方法生成作物水分下指示线，作物水分的上指示线与下指示线分别记为fu和fb：
[0065]fu
＝au×
cwdi+bu[0066]
fb＝ab×
cwdi+bb[0067]
其中au和bu为拟合得到的上指示线的斜率和截距，ab和bb为拟合得到的下指示线的斜率和截距。
[0068]
如图6所示，在农田灌溉概率指数建模中，气象干旱是土壤干旱的主要诱因，当土壤含水率变少，作物根系难以从土壤中吸收到足够的水分来补偿蒸腾消耗的水分时，就进一步诱发农业干旱。在没有灌溉补给时，气象干旱和农业干旱是高度相关的；而当灌溉发生时，气象干旱向农业干旱的演变过程就会受阻，农业干旱程度降低。在作物气象干旱程度相同的条件下，假设不同耕地之间的农业干旱差异主要是由灌溉补给决定的(忽略局部区域内诸如作物品种、地形、土壤类型等因素的影响)，则农业干旱程度越轻的像元是灌溉作物(或发生灌溉)的可能性越高，农业干旱程度越严重的农田像元是雨养作物的可能性越高。所述步骤(5)在作物水分散点空间中，以作物水分上、下指示线为约束，计算每个农田像元与作物水分上、下指示线在纵轴方向上的距离irr、pir，像元坐标为(cwdi,cwsi)，其中irr表征灌溉对作物水分胁迫的缓解程度，irr与pir之和mi表征农业灌溉对作物水分胁迫的最大缓解程度，irr与mi的比值表征该农田像元受到灌溉的可能性，农田像元的灌溉概率指数建模如下：
[0069]
irr＝fu(cwdi)-cwsi
[0070]
pir＝cwsi-fu(cwdi)或mi＝fu(cwdi)-fb(cwdi)
[0071][0072]
对于位于作物水分上指示线之上的像元，直接对其灌溉概率指数赋值为0；对于位于作物水分下指示线之下的像元，直接对其灌溉概率指数赋值为1.0。
[0073]
图7示意了利用本发明中的方法实施的某地农田灌溉概率指数的计算效果。从整体灌溉概率指数分布来看，本发明方法计算的灌溉概率指数整体上呈现东部高、西部低的
趋势，与该地区的气温和降水的分布趋势以及农业活动的变化梯度类似；灌溉概率指数较高的区域也基本是沿着河网分布(靠近水源)的。这恰好与我们对农业灌溉空间分布的判断和认识相一致。
[0074]
本发明的实例是在pc平台上实现的。经实验证明，本发明能够较简洁、准确地表征农田像元的灌溉可能性，并生成区域的农业灌溉空间分布图；后续的相关性验证表明：本发明所给实例的灌溉概率指数与县域灌溉面积的相关性达到0.82。本发明所提及的方法可广泛应用于区域农业灌溉信息的提取与分析程中，以支撑水资源管理、农业生产活动等决策。