CN111079256B - 一种基于遥感的灌溉水有效利用系数测算方法 - Google Patents

Abstract

本发明公开了一种基于遥感的灌溉水有效利用系数测算方法，包括步骤：(1)获取含有可见光、热红外波段的遥感影像，并进行预处理；(2)获取气象数据、数字高程模型数据；(3)获取灌区每一次灌溉用水量数据；(4)灌区总蒸散发量计算；(5)获取土壤类型数据、土地利用数据、植被数据、降水量数据，通过水文模型计算灌区内的总地表径流量；(6)测算灌溉水有效利用系数。本发明弥补了目前首尾测算法对于农业计量率要求高，复杂灌区实际操作过程中典型田块选取困难，人工监测范围有限，以点带面现象突出，并且测量结果受监测人员水平影响较大等不足，具有大尺度、高效、精确的优点。

Description

一种基于遥感的灌溉水有效利用系数测算方法

技术领域

本发明涉及农田水利研究领域，特别涉及一种基于遥感的灌溉水有效利用系数测算方法。

背景技术

灌区在灌水时，灌溉水除了一部分为被灌溉的农作物耗用外，还有一部分在输水、配水和灌水过程中损失掉，这部分水量成为非生产性水量损失。这些损失包括渗漏损失、蒸发损失、田面损失、泄水损失和跑水损失等。灌溉水有效利用系数是指灌溉时能够被农作物利用的净水量除以灌溉系统取用的灌溉总水量。灌溉水有效利用系数作为农业用水效率红线的主要指标，是最严格水资源管理考核的重要指标之一，对于节约用水具有重要意义。

目前灌溉水有效利用系数的测算主要采用的是首尾测算法(如图1)，该方法是指直接测量统计灌区从水源引入的毛灌溉用水总量，通过分析测算得到田间实际净灌溉用水总量，田间实际净灌溉用水总量与毛灌溉用水总量的比值即为灌溉水有效利用系数。但是这种方法对于灌区的计量要求较高。以中国为例，其农业用水计量率较低，整体不足55％，并且在其南方灌区具有地形复杂、灌溉图斑破碎、种植结构复杂、退排水很难监测等特征，导致该方法在实际操作过程中，存在典型田块选取困难、监测范围小、以点带面精度低、监测结果很难验证等不足，在考核期内较难及时对灌溉水有效利用系数进行准确确定，导致考核缺乏行之有效的抓手。

因此，研究一种大尺度、高效、精确的灌溉水有效利用系数测算方法具有很重要的研究意义和实用价值。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种基于遥感的灌溉水有效利用系数测算方法，该方法弥补了目前首尾测算法对于农业计量率要求高，对复杂灌区实际操作过程中典型田块选取困难，人工监测范围有限，以点带面现象突出，并且测量结果受监测人员水平影响较大等不足。

本发明的目的通过以下的技术方案实现：一种基于遥感的灌溉水有效利用系数测算方法，包括步骤：

(1)获取含有可见光、热红外波段的遥感影像，并进行预处理；

(2)获取气象数据、数字高程模型数据；

(3)获取灌区每一次灌溉用水量数据；

(4)灌区总蒸散发量(ET，Evapo-Transpiration)计算；

(5)获取土壤类型数据、土地利用数据、植被数据、降水量数据，通过水文模型计算灌区内的总地表径流量(R，Runoff)；

(6)测算灌溉水有效利用系数，公式是：

其中，IW为灌溉水有效利用系数；ET为灌区内的总蒸散发量；R为灌区内的总地表径流量；w为灌区的总引水量。

优选地，所述步骤(1)中含有可见光、热红外波段的遥感影像可以采用MODIS卫星影像，或者Landsat系列卫星影像。

优选地，所述步骤(1)中对遥感影像进行预处理，步骤是：根据卫星影像公布的定标参数，对遥感影像进行辐射定标、大气校正。

优选地，所述步骤(2)中获取的气象数据包括降水、大气气温、风速和大气湿度。

优选地，所述步骤(4)灌区总蒸散发量的计算方法是基于SEBAL模型的能量平衡算法。

更进一步的，基于SEBAL模型的能量平衡算法，步骤是：

λET＝R_n-G-H

其中，λ为蒸发潜热；ET为蒸散发量；R_n为净辐射量；G为土壤热通量；H为显热通量。

更进一步的，所述的净辐射量R_n的计算公式如下：

R_n＝(R_S↓-R_S↑)+(R_L↓-R_L↑)

其中，R_S↑和R_S↓为地面向上和向下的短波辐射项；R_L↑和R_L↓为地面向上和向下的长波辐射项。

优选的，在计算净辐射量R_n过程中，净短波辐射R_Sn为：

R_Sn＝R_S↓-R_S↑

其中，R_S↓的计算公式为：

G_SC为太阳常数，为1367w m^-2；d_e-s为日地相对距离，为天文单位；θ为太阳倾斜角，用弧度表示；τ_SW代表单向大气传输，采用高程函数进行计算，公式为：

τ_SW＝0.75+0.00002z

z代表海拔高度，可以通过步骤(2)中的数字高程模型获取；

其中，R_S↑的计算公式为：

R_S↑＝αR_S↓

α表示给定的比照率(albedo)，α_toa表示大气顶部反照率，α_ρ表示路径辐射反照率，其取值范围在0.025-0.04之间，这里取固定值0.03。

最终，净短波辐射可以表示为：

优选的，在计算净辐射量R_n过程中，净长波辐射R_Ln为：

R_Ln＝R_L↓-R_L↑

其中，R_L↓＝σε_αΤ_α ⁴，R_L↑＝σε_SΤ_S ⁴，σ是斯蒂芬玻尔兹曼常数(＝5.67×10-8Wm^-2K^-4)；Τ_α是大气温度；ε_α为大气比辐射率，计算公式为ε_α＝-0.85(lnτ_SW)^0.09；Τ_S是地表温度；ε_S为地表比辐射率，计算公式为ε_S＝1.009+0.047(lnNDVI)，NDVI取值大于0，当NDVI小于0时，取值为1。

更进一步的，地表温度Τ_S的计算公式为：

其中，R代表了热红外波段数据，K₁和K₂为定标常量。

更进一步的，所述的土壤热通量G的计算公式如下：

其中，R_n为净辐射量。

更进一步的，所述的显热通量H的计算公式如下：

H＝ρC_p(T₀-T_a)/r_a'

其中，ρ是空气密度(kg·m^-3)；C_p为标准大气压下的空气比热(1004J·K^-1·kg^-1)T₀为空气动力学温度，Τ_α是大气温度，SEBAL模型引入T₀与大气高度Τ_α的温差dT与Ts的线性关系假设，即：dT＝a+bT_s；空气动力学阻抗r_a′可表示为：

其中，z1是空气动力学阻抗r_a′对应的上边界，一般取值0.01m；z2为r_a′对应的上边界和下边界的高度；φ(z2,L)和φ(z1,L)为高度为z2和z1的热量校正长度；k为常量，u_*(200)是在100-200m高的掺混层的摩阻风速，可由该层高度的风速u200计算得到：

其中，下标200代表200m高度对应的参数值，k为常量，φm(z200)为迭代参数。

优选地，所述步骤(5)中植被数据采用获取的Landsat影像的近红外波段和红波段数据计算得到，计算公式：

其中，NIR为近红外波段，R为红波段。

本发明与现有技术相比，具有如下优点和有益效果：

针对目前首尾测算法对于农业计量率要求高，复杂灌区实际操作过程中典型田块选取困难，人工监测范围有限，以点带面现象突出，并且测量结果受监测人员水平影响较大等不足，本发明提出了一种基于遥感的灌溉水有效利用系数测算方法，相对于传统的首尾测算法，克服了对于灌区计量设施过分依赖、监测周期长、监测范围小、监测频次低等不足，通过现有网络数据以及遥感数据即可实现较为准确的灌溉水有效利用系数测算，具有效率高的优点。

附图说明

图1是传统首尾测算法模型示意图。

图2是本发明方法的流程图。

图3是本发明基于遥感的灌区灌溉水有效利用系数测算方法模型示意图。

图4(a)-(f)是本实施例计算得到的植被指数NDVI、地表温度、净辐射通量、土壤热通量、瞬时蒸散发量、日蒸散发量示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

如图1所示，本实施例给出一种基于遥感的灌溉水有效利用系数测算方法，该方法根据遥感影像、灌区灌溉用水量记录资料、土壤类型数据、土地利用数据、植被数据、降水量等数据，即可实现灌溉水有效利用系数测算。具体的，本实施例以中国广州市示例，详细说明利用本发明方法进行灌溉水有效利用系数测算的过程。根据图1所示的方法流程图，方法包括步骤：

S1：获取含有可见光、热红外波段的遥感影像，并进行预处理。

所述遥感影像主要可采用MODIS卫星影像和Landsat系列卫星影像，数据可以从各自官方网站下载：MODIS卫星(http://e4ftl01.cr.usgs.gov/)、Landsat系列卫星(http://data.cma.cn/)。本实施例中从USGS官网(http://data.cma.cn/)下载Landsat8卫星遥感影像，并根据卫星影像公布的定标参数，对遥感影像进行辐射定标、大气校正。处理完成的数据作为步骤S4中的输入数据。

S2：获取气象数据、数字高程模型数据。

本案例的研究区内安装有小型气象监测站，因此没有采用中国气象网提供的数据，主要监测的参数包括降水(mm)、大气气温(℃)、风速(m/s)和大气湿度(％)等，可通过中国气象数据网(http://data.cma.cn/)下载。监测的时间跨度为2017年1月1日至2017年12月31日，记录频率为分钟。

数字高程模型数据从地理空间数据云(http://www.giscloud.cn/)下载，数据图像空间分辨率30米。然后通过ArcGIS对数据进行拼接，得到研究区的DEM。

S3：获取灌区每一次灌溉用水量数据。

通过流溪河灌区管理处获取灌溉用水记录，累计灌溉用水总量为3863.55万方。

S4：灌区总蒸散发量计算。

蒸散发量的计算模型很多，本发明中选择目前全世界范围内使用的最广泛的能量平衡算法SEBAL模型。

SEBAL模型的具体计算公式为：

λET＝R_n-G-H (1)

其中，λ为蒸发潜热，J/m³；ET为蒸散发量，m/s；R_n为净辐射量，W/m²；G为土壤热通量，W/m²；H为显热通量，W/m²。

R_n可以由以下公式计算得到：

R_n＝(R_S↓-R_S↑)+(R_L↓-R_L↑) (2)

其中，R_S↑和R_S↓为地面向上和向下的短波辐射项，单位W·m^-2；R_L↑和R_L↓为地面向上和向下的长波辐射项，单位W·m^-2。

净短波辐射(R_Sn)为：

R_Sn＝R_S↓-R_S↑ (3)

净长波辐射(R_Ln)为：

R_Ln＝R_L↓-R_L↑ (4)

其中R_S↓的计算公式为：

其中，G_SC为太阳常数，为1367w·m^-2；d_e-s为日地相对距离，为天文单位；θ为太阳倾斜角，用弧度表示；τ_SW代表单向大气传输，采用高程函数进行计算，公式为：

τ_sW＝0.75+0.00002z (6)

其中，z代表了海拔高度(m)，可以通过步骤S2中的数字高程模型获取。

对于给定的比照率(albedo，α)R_S↑可以表述为：

R_S↑＝αR_S↓ (7)

α表示给定的比照率(albedo)，α_toa表示大气顶部反照率，α_ρ表示路径辐射反照率，其取值范围在0.025-0.04之间，这里取固定值0.03。

综合公式(3)、(5)、(7)，短波净辐射可以表示为：

向下长波辐射R_L↓计算公式为：

R_L↓＝σε_αΤ_α ⁴ (10)

其中，σ是斯蒂芬玻尔兹曼常数(＝5.67×10-8W·m^-2·K^-4)；Ta是大气温度(K)；εa为大气比辐射率(无量纲)，计算公式为：

ε_α＝-0.85(lnτ_SW)^0.09 (11)

向上长波辐射R_L↑计算公式为：

R_L↑＝σε_SΤ_S ⁴ (12)

其中，σ是斯蒂芬玻尔兹曼常数(＝5.67×10-8W·m^-2·K^-4)；Ts是地表温度(K)；εs为地表比辐射率(无量纲)，计算公式为：

ε_S＝1.009+0.047(lnNDVI) (13)

在公式(13)中NDVI取值大于0，当NDVI小于0时，取值为1。

其中，R代表了热红外波段数据，K1和K2为定标常量。

H＝ρC_p(T₀-T_a)/r_a' (16)

其中，ρ是空气密度(kg·m^-3)；C_p为标准大气压下的空气比热(1004J·K^-1·kg^-1)，T₀为空气动力学温度，℃，T₀不能直接用遥感反演的Ts代替，因此SEBAL模型引入T₀与参考高度Ta的温差dT与Ts的线性关系假设，即：

dT＝a+bT_s (17)

空气动力学阻抗r_a′可表示为：

其中，z1是空气动力学阻抗r_a′对应的上边界，一般取值0.01m；z2为r_a′对应的上边界和下边界的高度，一般取值2m；φ(z2,L)和φ(z1,L)为高度为z2和z1的热量校正长度；k为常量0.41，u_*(200)是在100-200m高的掺混层的摩阻风速，可由该层高度的风速u200计算得到：

其中，下标200代表200m高度对应的参数值。k为常量0.41，Z200取值200，zom取值0.5，φm(z200)为迭代参数。本实施例基于上述SEBAL模型，计算得到植被指数NDVI、地表温度、净辐射通量、土壤热通量、日蒸散发量，并对计算的结果进行了校正，进一步采用时间尺度拓展的方式获得整个研究区内的蒸散发总量2074.67万方。

S5：获取土壤类型数据、土地利用数据、植被数据、降水量等数据，通过水文模型计算地表径流量。

本实施例中，包括步骤如下：

S5-1：获取土壤类型数据(来源：世界土壤数据库(Harmonized World SoilDatabase version 1.1))，

S5-2：获取土地利用数据(来源：MODIS卫星土地利用产品数据库(http://e4ftl01.cr.usgs.gov/))；

S5-3：植被数据采用获取的Landsat影像的近红外波段和红波段数据，计算公式：

其中，NIR为近红外波段，R为红波段；

S5-4：降水量数据通过中国气象数据网获取(http://data.cma.cn/)；

S5-5：最后输入进入常用的水文模型，如SWAT、SCS-CN等计算得到地表径流量。

本实施例中，基于上述数据采用SWAT水文模型，根据研究区内的土地利用数据、降水量、植被指数NDVI、土壤类型计算得到地表径流总量235.62万方。

S6：测算灌溉水有效利用系数。

融合可见光-热红外遥感的灌溉水有效利用系数计算公式：

其中，IW为灌溉水有效利用系数；ET为灌区内的总蒸散发量；R为灌区内的总地表径流量；w为灌区的总引水量。

根据本发明提出的灌溉水有效系数的计算公式得到研究区内灌溉水有效利用系数为0.476。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (8)

1.一种基于遥感的灌溉水有效利用系数测算方法，其特征在于，包括步骤：

(1)获取含有可见光、热红外波段的遥感影像，并进行预处理；

遥感影像采用MODIS卫星影像或Landsat系列卫星影像，根据卫星影像公布的定标参数，对遥感影像进行辐射定标、大气校正预处理；

(2)获取气象数据、数字高程模型数据；

气象数据通过灌区小型气象监测站监测或通过中国气象数据网下载，气象监测站监测的参数包括降水、大气气温、风速和大气湿度；

数字高程模型数据从地理空间数据云下载，通过ArcGIS对数据进行拼接，得到灌区的数字高程模型DEM；

(3)获取灌区每一次灌溉用水量数据；

(4)灌区总蒸散发量计算；

(5)获取土壤类型数据、土地利用数据、植被数据、降水量数据，通过水文模型计算灌区内的总地表径流量；

通过世界土壤数据库获取土壤类型数据；

通过MODIS卫星土地利用产品数据库获取土地利用数据；

植被数据采用获取的Landsat影像的近红外波段和红波段数据计算得到，计算公式：

其中，NIR为近红外波段，R为红波段；

降水量数据通过中国气象数据网获取；

将土壤类型数据、土地利用数据、植被数据、降水量数据输入SWAT或SCS-CN水文模型，计算得到地表径流量；

(6)测算灌溉水有效利用系数，公式是：

其中，IW为灌溉水有效利用系数；ET为灌区内的总蒸散发量；R为灌区内的总地表径流量；w为灌区的总引水量。

2.根据权利要求1所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，所述步骤(4)灌区总蒸散发量的计算方法是基于SEBAL模型的能量平衡算法。

3.根据权利要求2所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，基于SEBAL模型的能量平衡算法，步骤是：

λET＝R_n-G-H

其中，λ为蒸发潜热；ET为蒸散发量；R_n为净辐射量；G为土壤热通量；H为显热通量。

4.根据权利要求3所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，所述的净辐射量R_n的计算公式如下：

R_n＝(R_S↓-R_S↑)+(R_L↓-R_L↑)

其中，R_S↑和R_S↓为地面向上和向下的短波辐射项；R_L↑和R_L↓为地面向上和向下的长波辐射项。

5.根据权利要求4所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，在计算净辐射量R_n过程中，净短波辐射R_Sn为：

R_Sn＝R_S↓-R_S↑

其中，R_S↓的计算公式为：

G_SC为太阳常数；d_e-s为日地相对距离；θ为太阳倾斜角；τ_SW代表单向大气传输，计算公式为：

τ_SW＝0.75+0.00002z

z代表海拔高度，通过步骤(2)中的数字高程模型获取；

其中，R_S↑的计算公式为：

R_S↑＝αR_S↓

α表示给定的比照率，α_toa表示大气顶部反照率，α_ρ表示路径辐射反照率；

最终，净短波辐射可以表示为：

6.根据权利要求4所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，在计算净辐射量R_n过程中，净长波辐射R_Ln为：

R_Ln＝R_L↓-R_L↑

其中，R_L↓＝σε_αT_α ⁴，R_L↑＝σε_ST_S ⁴，σ是斯蒂芬玻尔兹曼常数；T_α是大气温度；ε_α为大气比辐射率，计算公式为ε_α＝-0.85(lnτ_SW)^0.09；T_S是地表温度；ε_S为地表比辐射率，计算公式为ε_S＝1.009+0.047(ln NDVI)，NDVI取值大于0，当NDVI小于0时，取值为1。

7.根据权利要求6所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，地表温度T_S的计算公式为：

其中，R代表了热红外波段数据，K₁和K₂为定标常量。

8.根据权利要求3所述的基于遥感的灌溉水有效利用系数测算方法，其特征在于，所述的土壤热通量G的计算公式如下：

其中，R_n为净辐射量；

所述的显热通量H的计算公式如下：

H＝ρC_p(T₀-T_a)/r_a'

其中，ρ是空气密度；C_p为标准大气压下的空气比热，T₀为空气动力学温度，T_α是大气温度，SEBAL模型引入T₀与大气高度T_α的温差dT与Ts的线性关系假设，即：dT＝a+bT_s；空气动力学阻抗r_a′可表示为：

其中，z1是空气动力学阻抗r_a′对应的上边界；z2为r_a′对应的上边界和下边界的高度；φ(z2,L)和φ(z1,L)为高度为z2和z1的热量校正长度；k为常量，u_*(200)是在100-200m高的掺混层的摩阻风速，可由该层高度的风速u200计算得到：

其中，下标200代表200m高度对应的参数值，φm(z200)为迭代参数。