一种基于遥感数据的参考作物腾发量计算方法

摘要

本发明提出一种基于遥感数据的参考作物腾发量计算方法，属于定量灌溉设计与农业水资源管理领域。本发明首先获取研究区域训练气象数据并计算对应的训练参考作物腾发量ET0作为训练集的输出数据组，获取研究区域的训练遥感数据作为训练集的输入数据组，建立研究区域的参考作物腾发量计算模型；然后通过验证ET0和验证遥感数据检验模型精度；使用时，获取所需时间段的遥感数据并输入验证完毕的模型，模型输出该研究区域选定时间段内参考作物腾发量。本发明结合少量的气象站资料与易获取的遥感资料，可实际应用于气象站数据缺失地区和气象站稀疏地区的区域参考作物腾发量计算，有效的协助当地灌溉规划和农业用水管理。

权利要求

1.一种基于遥感数据的参考作物腾发量计算方法，其特征在于，包括以下步骤：

1)建立研究区域的参考作物腾发量计算模型；具体步骤如下：

1-1)确定研究区域，在该区域内及周边选取获取训练气象数据的气象站点；

选定任意研究区域，在研究区域范围内或在距离区域边界100公里内的研究区域周边具有相近气候、地理状况的区域选取气象站点；统计气象站点的总数量，随机选取数量为总数量80％的气象站点作为获取训练气象数据的气象站点；

1-2)获取训练气象数据；

从步骤1-1)选取的获取训练气象数据的气象站点获取气象数据与气象站点地理信息，包括：日最高气温Tmax、日最低气温Tmin、日平均气温Tmean、平均水汽压ea、平均风速u、实际日照时数N、相对湿度RH、风速测点高度H，气象站点高程Z、气象站点经度β和气象站点纬度每个气象站点选取的气象数据时间长度大于等于1年；将获取的气象数据记为训练气象数据；

1-3)利用步骤1-2)获取的训练气象数据，计算训练ET0；

参考作物腾发量ET0通过彭曼一蒙特斯公式计算，表达式如下：式中，Rn为净辐射，G为土壤热通量，u2为2米高度处的风速，es为饱和水汽压，ea为平均水汽压，Δ为饱和水汽压—温度曲线斜率，γ为湿度计常数；

1-4)获取研究区域的训练遥感数据；

确定遥感数据类型，根据步骤1-2)中每个获取训练气象数据的气象站点的经度和纬度，选取对所有获取训练气象数据的气象站点进行覆盖的遥感数据，所述遥感数据采用栅格形式的数据，遥感数据观测的时间与气象站点观测的时间对应；将获取的遥感数据作为训练遥感数据；

1-5)建立研究区域的参考作物腾发量计算模型；

将步骤1-4)获取的训练遥感数据作为训练集的输入数据组，将步骤1-3)得到的训练ET0作为训练集的输出数据组，利用机器学习模型进行训练，建立研究区域的参考作物腾发量计算模型；

2)对步骤1)建立的模型进行验证；具体步骤如下：

2-1)选取获取验证气象数据的气象站点；

获取验证气象数据的气象站点为步骤1-1)中未选中作为获取训练气象数据气象站点的剩余气象站点；

2-2)获取验证气象数据；

重复步骤1-2)，从步骤2-1)选取的获取验证气象数据的气象站点获取气象数据与气象站点地理信息；每个气象站点选取的气象数据时间长度大于等于1年；将获取的气象数据记为验证气象数据；

2-3)重复步骤1-3)，利用步骤2-2)获取的验证气象数据计算验证ET0；

2-4)重复步骤1-4)，获取验证遥感数据；验证遥感数据选取与步骤1-4)中训练遥感数据相同类型的遥感数据；

2-5)对步骤1)建立的模型进行验证；

将验证遥感数据作为验证集的输入数据组输入步骤1)建立的模型，该模型输出对应的模型估算ET0；使用精度评价指标，将模型估算ET0与验证ET0进行对比，评价模型精度；所述精度评价指标包括：均方根误差RMSE，相关系数CC和相对偏差BIAS，计算公式分别如下：式中，n为验证ET0的个数，与模型估算ET0的个数相同；Pi为第i个模型估算ET0，Ai为第i个验证ET0，Pmean为模型估算ET0的均值，Amean为验证ET0的均值；

将每个精度评价指标的计算结果与研究区域对于该精度评价指标设定的阈值进行比较并判定：若三个精度评价指标均符合对应的阈值要求，则模型验证完毕，步骤1)建立的模型即为最终得到的该研究区域参考作物腾发量计算模型，进入步骤3)；若任一精度评价指标不符合其对应的阈值要求，则重新返回步骤1-4)，重新选取遥感数据类型，或重新返回步骤1-

5)更换机器学习模型，直到三个精度评价指标均符合对应的阈值要求，模型验证完毕，进入步骤3)；

3)利用步骤2)验证完毕的模型，计算研究区域选定时间段内的参考作物腾发量；具体步骤如下：

3-1)选定实际使用需求的时间段，根据步骤1-4)中确定的遥感数据类型，获取该时间段内对应的遥感数据；

3-2)将步骤3-1)获取的遥感数据输入步骤2)验证完毕的模型，模型输出对应时间段内的ET0，即为该研究区域选定时间段内遥感数据覆盖范围内的参考作物腾发量。

说明书

一种基于遥感数据的参考作物腾发量计算方法

技术领域

[0001] 本发明属于定量灌溉设计与农业水资源管理领域，具体涉及一种基于遥感数据的参考作物腾发量计算方法。

背景技术

[0002] 参考作物腾发量(ET0)是水文循环和水资源管理的重要组成部分。准确估计ET0对于提高灌溉效率，干旱监测，水资源重复利用，农业水资源调度具有重要意义。ET0被定义为具有特定特征的假设作物参考表面的蒸散量，在农业应用中，ET0常用于计算实际条件下的作物蒸散量以指导灌溉设计。许多研究表明，详细而合理的灌溉系统可以提高灌溉水利用效率。因此，对区域ET0进行可靠的估算是灌溉规划和农业用水管理的重要基础。

[0003] 当前的水资源问题尤其农业用水管理问题中需要进行ET0空间尺度的估算，然而许多发展中和贫困地区由于缺乏必要的气象数据，难以良好估算ET0的时空变化。粮农组织提出彭曼一蒙特斯公式(FAO-PM)方法作为估算ET0的标准方法，FAO-PM方法需要大气温度，风速，相对湿度和太阳辐射数据来估计ET0。然而能够观察完整参数的气象站数量常常是有限的，且在发展中国家，风速，湿度，辐射数据的收集常存在缺失。而当气象站不够密集时，气象站对区域的代表性较弱，也难以获得可代表整个区域的参考作物腾发量值。因此，这些地区的灌溉管理难以达到最佳状态，增加了水资源短缺和水资源冲突的风险。

[0004] 由于传统ET0估算方法的局限性，利用较少的气象数据估算区域ET0的方法成为当前研究热点。近年来一些研究使用地表温度、地表温度与风速、地表温度与日照时数等数据结合人工神经网络方法估算ET0，以期使用更少的气象数据计算ET0；但这些方法局限于气象站点，无法在气象站点以外的区域进行合理应用，难以获取区域性ET0数据。近年来遥感技术的发展也使得遥感数据在水文及水资源领域的运用范围越来越广，Maeda等(2011)应用MODIS地表温度替代Hargreaves模型，Thornthwaite模型和Blaney-Criddle方程的大气温度进行ET0计算。结果表明Hargreaves模型最为合适，均方根误差平均为0.47mm/d，相关系数为0.67。但该方法简单的利用MODIS地表温度数据替代气温，相关系数未达到较好的水平，精度难以满足实际运用要求。

发明内容

[0005] 本发明的目的是为克服已有技术中气象站数据缺失地区和气象站稀疏地区的区域参考作物腾发量难以获得的问题，提出一种基于遥感数据的参考作物腾发量计算方法。本发明结合少量的气象站资料与易获取的遥感资料，利用机器学习方法，建立计算区域性的参考作物腾发量的模型；可实际应用于气象站数据缺失地区和气象站稀疏地区的区域参考作物腾发量计算，有效的协助当地灌溉规划和农业用水管理。

[0006] 本发明提出一种基于遥感数据的参考作物腾发量计算方法，其特征在于，包括以下步骤：

[0007] 1)建立研究区域的参考作物腾发量计算模型；具体步骤如下：

[0008] 1-1)确定研究区域，在该区域内及周边选取获取训练气象数据的气象站点；

[0009] 选定任意研究区域，在研究区域范围内或在距离区域边界100公里内的研究区域周边具有相近气候、地理状况的区域选取气象站点；统计气象站点的总数量，随机选取数量为总数量80％的气象站点作为获取训练气象数据的气象站点；

[0010] 1-2)获取训练气象数据；

[0011] 从步骤1-1)选取的获取训练气象数据的气象站点获取气象数据与气象站点地理信息，包括：日最高气温Tmax、日最低气温Tmin、日平均气温Tmean、平均水汽压ea、平均风速u、实际日照时数N、相对湿度RH、风速测点高度H，气象站点高程Z、气象站点经度β和气象站点纬度 每个气象站点选取的气象数据时间长度大于等于1年；将获取的气象数据记为训练气象数据；

[0012] 1-3)利用步骤1-2)获取的训练气象数据，计算训练ET0；

[0013] 参考作物腾发量ET0通过彭曼一蒙特斯公式计算，表达式如下：

[0014]

[0015] 式中，Rn为净辐射，G为土壤热通量，，u2为2米高度处的风速，es为饱和水汽压，ea为平均水汽压，Δ为饱和水汽压—温度曲线斜率，γ为湿度计常数；

[0016] 1-4)获取研究区域的训练遥感数据；

[0017] 确定遥感数据类型，根据步骤1-2)中每个获取训练气象数据的气象站点的经度和纬度，选取对所有获取训练气象数据的气象站点进行覆盖的遥感数据，所述遥感数据采用栅格形式的数据，遥感数据观测的时间与气象站点观测的时间对应；将获取的遥感数据作为训练遥感数据；

[0018] 1-5)建立研究区域的参考作物腾发量计算模型；

[0019] 将步骤1-4)获取的训练遥感数据作为训练集的输入数据组，将步骤1-3)得到的训练ET0作为训练集的输出数据组，利用机器学习模型进行训练，建立研究区域的参考作物腾发量计算模型；

[0020] 2)对步骤1)建立的模型进行验证；具体步骤如下：

[0021] 2-1)选取获取验证气象数据的气象站点；

[0022] 获取验证气象数据的气象站点为步骤1-1)中未选中作为获取训练气象数据气象站点的剩余气象站点；

[0023] 2-2)获取验证气象数据；

[0024] 重复步骤1-2)，从步骤2-1)选取的获取验证气象数据的气象站点获取气象数据与气象站点地理信息；每个气象站点选取的气象数据时间长度大于等于1年；将获取的气象数据记为验证气象数据；

[0025] 2-3)重复步骤1-3)，利用步骤2-2)获取的验证气象数据计算验证ET0；

[0026] 2-4)重复步骤1-4)，获取验证遥感数据；验证遥感数据选取与步骤1-4)中训练遥感数据相同类型的遥感数据；

[0027] 2-5)对步骤1)建立的模型进行验证；

[0028] 将验证遥感数据作为验证集的输入数据组输入步骤1)建立的模型，该模型输出对应的模型估算ET0；使用精度评价指标，将模型估算ET0与验证ET0进行对比，评价模型精度；所述精度评价指标包括：均方根误差RMSE，相关系数CC和相对偏差BIAS，计算公式分别如下：

[0029]

[0030]

[0031]

[0032] 式中，n为验证ET0的个数，与模型估算ET0的个数相同；Pi为第i个模型估算ET0，Ai为第i个验证ET0，Pmean为模型估算ET0的均值，Amean为验证ET0的均值；

[0033] 将每个精度评价指标的计算结果与研究区域对于该精度评价指标设定的阈值进行比较并判定：

[0034] 若三个精度评价指标均符合对应的阈值要求，则模型验证完毕，步骤1)建立的模型即为最终得到的该研究区域参考作物腾发量计算模型，进入步骤3)；若任一精度评价指标不符合其对应的阈值要求，则重新返回步骤1-4)，重新选取遥感数据类型，或重新返回步骤1-5)更换机器学习模型，直到三个精度评价指标均符合对应的阈值要求，模型验证完毕，进入步骤3)；

[0035] 3)利用步骤2)验证完毕的模型，计算研究区域选定时间段内的参考作物腾发量；具体步骤如下：

[0036] 3-1)选定实际使用需求的时间段，根据步骤1-4)中确定的遥感数据类型，获取该时间段内对应的的遥感数据；

[0037] 3-2)将步骤3-1)获取的遥感数据输入步骤2)验证完毕的模型，模型输出对应时间段内的ET0，即为该研究区域选定时间段内遥感数据覆盖范围内的参考作物腾发量。

[0038] 本发明的特点及有益效果在于：

[0039] 本发明提出的一种基于遥感数据的参考作物腾发量计算方法，利用易获取的遥感数据，减少已有的ET0计算方法对气象数据的依赖。使用机器学习方法，可以有效地找到遥感数据与ET0之间的关系，弥补了气象站数据缺失地区和气象站稀疏地区的区域参考作物腾发量难以估算的问题，为灌溉规划和农业用水管理提供有效的参考作物腾发量资料，具有使用简便、可操作性强、容易得到实际应用和推广的优点。

具体实施方式

[0040] 本发明提出一种基于遥感数据的参考作物腾发量计算方法，下面结合具体实施例进一步详细说明如下。

[0041] 本发明提出的一种基于遥感数据的参考作物腾发量计算方法，包括以下步骤：

[0042] 1)建立研究区域的参考作物腾发量计算模型；具体步骤如下：

[0043] 1-1)确定研究区域，在该区域内及周边选取获取训练气象数据的气象站点；

[0044] 本发明中研究区域可为任意区域，气象站点的选取要求为气象站点的位置在研究区域范围内，或在研究区域周边(距离区域边界100公里内)具有相近气候、地理状况的区域。统计可供选择的气象站点的总数量，而后在可供选择的气象站点中，随机选取数量为总数量80％的气象站点(若有小数则四舍五入取整)，作为获取训练气象数据的气象站点。本实施例中，选取的研究区域为中国甘肃省的疏勒河流域，研究区域范围内气象站点6个，研究区域周边(距离区域边界100公里内)具有相近气候、地理状况的气象站点5个，总数量为11个；随机选取9个(数量为总数量的80％，取整后为9个)气象站点，作为训练气象数据的气象站点来源。

[0045] 1-2)获取训练气象数据；

[0046] 从步骤1-1)选取的获取训练气象数据的气象站点获取气象数据与气象站点地理信息，包括：日最高气温(Tmax)、日最低气温(Tmin)、日平均气温(Tmean)、平均水汽压(ea)、平均风速(u)、实际日照时数(N)、相对湿度(RH)、风速测点高度(H)，气象站点高程(Z)、气象站点经度(β)和气象站点纬度 气象数据简称、名称与单位如表1所示。

[0047] 表1气象数据简称、名称与单位列表

[0048]

[0049] 将获取的气象数据记为训练气象数据，训练气象数据需满足时间长度要求，每个气象站点选取的气象数据时间长度需大于等于1年。本实施例中的9个气象站点均选取2010年1月1日至2013年12月31日的气象数据，符合时间长度要求。

[0050] 1-3)利用步骤1-2)获取的训练气象数据，计算训练ET0；

[0051] 本发明中，ET0通过FAO-PM公式计算，表达式如下：

[0052]

[0053] 式中，ET0为参考作物腾发量，单位：mm/day，Rn为净辐射，单位：MJ/(m2day)，G为土壤热通量，单位：MJ/(m2day)，由于日尺度的土壤热通量相对很小，一般可以忽略。Tmean为日平均气温，单位：℃，u2为2米高度处的风速，单位：m/s，es为饱和水汽压，单位：kPa，ea为平均水汽压，单位：kPa，Δ为饱和水汽压—温度曲线斜率，单位：kPa/℃，γ为湿度计常数，单位：kPa/℃。

[0054] 1-3-1)饱和水汽压es和饱和水汽压—温度曲线斜率Δ的计算方法如下：

[0055] 饱和水汽压es为日最低气温Tmin时的饱和水汽压和日最高气温Tmax时的饱和水汽压的平均值，表达式如下：

[0056]

[0057] 其中，气温为T时的饱和水汽压e0(T)可以根据下式计算：

[0058]

[0059] 式中，e0(T)为气温为T时的饱和水汽压，单位：kPa。

[0060] 饱和水汽压—温度曲线斜率Δ可通过下式计算：

[0061]

[0062] 1-3-2)湿度计常数γ计算表达式如下：

[0063]

[0064] 式中，P为气压，单位：kPa；λ为水的汽化潜热，单位：J/g；

[0065] 其中，气压P计算表达式如下：

[0066]

[0067] 式中，Z为气象站点高程；

[0068] 1-3-3)2m处风速u2计算方法如下：

[0069] 通常气象站点的风速测点高度不一定为2m，为获得2m处风速，需进行换算。通过气象站点获得的平均风速u与风速测点高度H换算2m处风速，风速换算公式：

[0070]

[0071] 式中，u为气象站点平均风速；

[0072] 1-3-4)净辐射Rn计算方法如下：

[0073] 日地相对距离的倒数由下式得到：

[0074]

[0075] 式中，dr为日地相对距离的倒数；J为日序数(令1月1日为1，逐日增加)。

[0076] 日倾角计算公式：

[0077]

[0078] 式中，δ为日倾角。

[0079] 日照时数角计算公式：

[0080]

[0081] 式中，ωs为日照时数角；为地理纬度。

[0082] 最大可能日照时数计算公式：

[0083]

[0084] 式中，Ns为最大可能日照时数。

[0085] 大气顶太阳辐射计算公式：

[0086]

[0087] 式中，Ra为大气顶太阳辐射，单位：MJ/(m2·d)。

[0088] 到达地球表面的实际太阳辐射计算公式：

[0089]

[0090] 式中，Rs为到达地球表面的实际太阳辐射，单位：MJ/(m2·d)；N为实际日照时数。

[0091] 到达地球表面的晴空太阳辐射计算公式：

[0092]

[0093] 式中，Rso为到达地球表面的晴空太阳辐射，单位：MJ/(m2·d)；Z为气象站点高程高程。

[0094] 净短波辐射计算公式：

[0095] Rns＝0.77Rs(15)

[0096] 式中，Rns为净短波辐射，单位：MJ/(m2day)。

[0097] 净长波辐射计算公式：

[0098]

[0099] 式中，Rnl为净长波辐射，单位：MJ/(m2day)；σ为斯蒂芬-玻尔兹曼(Stefan-Boltzmann)常数，σ＝4.903×10-9MJ/(K4·m2·d)；Tmax,K为最高绝对温度，Tmax,K＝Tmax+273.16；Tmin,K为最低绝对温度，Tmin,K＝Tmin+273.16。

[0100] 净辐射计算公式：

[0101] Rn＝Rns-Rnl   (17)

[0102] 1-4)获取研究区域的训练遥感数据；

[0103] 本发明中，机器学习的输入数据组由遥感数据组成，遥感数据采用栅格形式的数据。

[0104] 确定遥感数据类型，选取对步骤1-1)中所有获取训练气象数据的气象站点进行覆盖的遥感数据，根据步骤1-2)中每个获取训练气象数据的气象站点的经纬度，使用Matlab或其它编程方法提取每个气象站点对应的遥感栅格值作为训练遥感数据。遥感数据观测的时间，与气象站点观测的时间对应，气象站点观测的时间为一段时间内的日数据，则遥感数据也为这一段时间的日数据。遥感数据可使用的类型包含但不限于Modis，Landsat，SPOT，高分一号等。本实施例中采用的训练遥感数据为Modis遥感白天地表温度数据、遥感夜晚地表温度数据、遥感地表反射率三种数据产品。

[0105] 1-5)建立研究区域的参考作物腾发量计算模型；

[0106] 将步骤1-4)获取的训练遥感数据作为训练集的输入数据组，将步骤1-3)得到训练ET0作为训练集的输出数据组，利用机器学习模型进行训练，建立研究区域的参考作物腾发量计算模型。其中，通过从每个气象站点获取的每日气象数据计算得到的该气象站点每日ET0与该日覆盖该气象站点的遥感数据是一一对应的。可使用的机器学习模型包含但不限于支持向量机(SVM)、自适应神经模糊推理系(ANFIS)、BP神经网络(BP)等。将训练集的输入数据组和输出数据组输入到机器训练模型后，机器学习模型就能够建立起输出数据和输入数据之间的关系，生成参考作物腾发量计算模型。

[0107] 2)对步骤1)建立的模型进行验证；具体步骤如下：

[0108] 2-1)选取获取验证气象数据的气象站点；

[0109] 获取验证气象数据的气象站点为步骤1-1)中未选中作为获取训练气象数据气象站点的剩余气象站点。本实施例中研究区域的气象站点总数量为11个，已随机选取9个作为训练气象数据的来源，未选取的气象站点为2个，则以这2个气象站点作为验证气象数据的来源。

[0110] 2-2)获取验证气象数据；

[0111] 重复步骤1-2)，从步骤2-1)选取的获取验证气象数据的气象站点获取气象数据与气象站点地理信息；每个气象站点选取的气象数据时间长度大于等于1年，时间的起始时间与结束时间可随机选取，将获取的气象数据记为验证气象数据。本实施例中一个获取验证气象数据的气象站点选取的气象数据时间长度为2014年9月5日至2016年12月31日，另一个气象站点选取的气象数据时间长度为2014年1月1日至2014年9月29日和2015年7月14日至2016年12月31日。

[0112] 2-3)计算验证ET0；

[0113] 重复步骤1-3)，利用步骤2-2)获取的验证气象数据通过FAO-PM公式计算验证ET0。

[0114] 2-4)获取验证遥感数据；

[0115] 重复步骤1-4)，验证遥感数据选取与步骤1-4)中训练遥感数据相同类型的遥感数据，验证遥感数据对步骤2-1)中所有获取验证气象数据的气象站点进行覆盖，根据气象站点的经纬度，使用Matlab或其它编程方法提取气象站点对应的遥感栅格值。验证遥感数据观测的时间也与获取验证气象数据的气象站点观测的时间对应，气象站点观测的时间为一段时间内的日数据，则遥感数据也为这一段时间的日数据。本实施例中采用的验证遥感数据为Modis遥感白天地表温度数据、遥感夜晚地表温度数据、遥感地表反射率三种数据产品。

[0116] 2-5)对步骤1)建立的模型进行验证；

[0117] 将验证遥感数据作为验证集的输入数据组输入步骤1)建立的模型，该模型输出对应的模型估算ET0；将模型估算ET0与步骤2-3)获取的验证ET0进行比较；

[0118] 本发明使用精度评价指标，将模型估算ET0与验证ET0进行对比，评价模型精度。可使用的精度评价指标包括：均方根误差(RMSE)，相关系数(CC)，相对偏差(BIAS)，计算公式分别如下：

[0119]

[0120]

[0121]

[0122] 式中，n为验证ET0的个数，与验证ET0与所选取的天数对应，与模型估算ET0的个数相同。Pi为第i个模型估算ET0，Ai为第i个验证ET0，Pmean为模型估算ET0的均值，Amean为验证ET0的均值。

[0123] 依据具体区域使用ET0的精度要求，根据精度评价指标均方根误差(RMSE)，相关系数(CC)，相对偏差(BIAS)判断模型精度是否符合对应的阈值要求。每个精度评价指标的对应阈值通常使用实际区域的规范与历史经验。本实施例中要求相关系数大于等于0.9，均方根误差小于等于1.5mm，相对偏差小于等于10％。

[0124] 若三个精度评价指标均符合对应的阈值要求，则模型验证完毕，步骤1)建立的模型即为最终得到的该研究区域参考作物腾发量计算模型，进入步骤3)；若任一精度评价指标不符合其对应的阈值要求，则重新返回步骤1-4)，重新选取遥感数据类型，或重新返回步骤1-5)更换机器学习模型，直到三个精度评价指标均符合对应的阈值要求，模型验证完毕，进入步骤3)。

[0125] 3)利用步骤2)验证完毕的模型，计算研究区域选定时间段内的参考作物腾发量；具体步骤如下：

[0126] 3-1)选定实际使用需求的时间段(可为任意时间长度)，根据步骤1-4)中确定的遥感数据类型，获取该时间段内对应的遥感数据(遥感数据可覆盖研究区的全部范围或部分范围)。本实施例中使用研究区域全部范围的2017年1月1日至2017年6月31日遥感数据。

[0127] 3-2)将步骤3-1)获取的遥感数据输入步骤2)验证完毕的模型，模型输出对应时间段内的ET0，即为该研究区域选定时间段内遥感数据覆盖范围内的参考作物腾发量。本实施例中输出研究区域全部范围的2017年1月1日至2017年6月31日的每日ET0。