Farmland-scale evapotranspiration estimation method

Abstract

This invention belongs to the field of agricultural remote sensing technology and relates to a method for estimating farmland-scale evapotranspiration, comprising: 1) acquiring basic data of the study area based on a UAV; the basic data includes multispectral image data, thermal infrared image data, and meteorological data; 2) preprocessing the multispectral image data and thermal infrared image data obtained in step 1); 3) obtaining instantaneous farmland evapotranspiration based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2); 4) extending the instantaneous farmland evapotranspiration obtained in step 3) to a daily scale to complete the estimation of farmland-scale evapotranspiration. This invention provides a method for estimating farmland-scale evapotranspiration that can effectively estimate farmland-scale evapotranspiration while ensuring high estimation accuracy.

Description

Technical Field

[0001] This invention belongs to the field of agricultural remote sensing technology and relates to a method for estimating farmland-scale evapotranspiration, and more particularly to a method for estimating farmland-scale evapotranspiration based on multi-source remote sensing from unmanned aerial vehicles. Background Technology

[0002] Farmland evapotranspiration is the total flux of water vapor transported from crops and their surrounding soil into the atmosphere. It is an essential piece of information for agricultural water management, and determining its value can reduce ineffective water consumption during crop growth, providing a basis for developing scientific irrigation systems, guiding precise irrigation decisions, and improving water resource utilization efficiency. Therefore, it is crucial to establish monitoring methods suitable for farmland evapotranspiration. Traditional ground-based methods for farmland evapotranspiration observation mainly include eddy covariance, ripple ratio, and lysimeter methods. While these methods provide relatively accurate evapotranspiration observations, they are costly to implement, difficult to maintain, and have low universality. Currently, the highest spatiotemporal resolution of satellite imagery available for estimating ground evapotranspiration is 8 days and 30 meters, which is insufficient for estimating evapotranspiration at the farmland scale. Furthermore, the estimation process requires data from ground meteorological stations, making it difficult to apply. In contrast, unmanned aerial vehicles (UAVs) offer flexible operating times and spatial resolution down to the decimeter level. By using sensors mounted on UAVs to acquire the information needed for evapotranspiration estimation and combining it with the Earth's surface energy balance equation, the spatial distribution information of farmland evapotranspiration can be quickly obtained. Summary of the Invention

[0003] Purpose of the invention: In order to solve the above-mentioned technical problems in the background art, the present invention provides a method for estimating farmland-scale evapotranspiration that can effectively estimate farmland-scale evapotranspiration while ensuring high estimation accuracy.

[0004] Technical Solution: To achieve the above objectives, the present invention adopts the following technical solution: a method for estimating farmland-scale evapotranspiration, comprising the following steps:

[0005] 1) Acquire basic data of the study area using UAVs; the basic data includes multispectral image data, thermal infrared image data, and meteorological data;

[0006] 2) Preprocess the multispectral image data and thermal infrared image data obtained in step 1);

[0007] 3) Based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2), obtain the instantaneous farmland evapotranspiration.

[0008] 4) Extend the instantaneous farmland evapotranspiration obtained in step 3) to a daily scale to complete the farmland-scale evapotranspiration estimation.

[0009] Preferably, in step 1), the multispectral image data is acquired by a multispectral imager mounted on a UAV, and the multispectral image data includes blue light band image data, green light band image data, red light band image data, red edge band image data, and near-infrared band image data; the thermal infrared image data is acquired by a thermal infrared imager mounted on a UAV, and the thermal infrared image data includes surface temperature; the meteorological data is acquired by a meteorological sensor mounted on a UAV, and the meteorological data includes wind speed, air temperature, and air pressure.

[0010] Preferably, the preprocessing in step 2) is image registration, image stitching, and radiometric calibration, while ensuring that the sampling standards of the multispectral imager and the thermal infrared imager are consistent, for example, resampling to 0.2m resolution.

[0011] Preferably, step 3) is implemented as follows:

[0012] 3.1) Calculate the net surface radiation Rn based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2);

[0013] 3.2) Calculate the soil heat flux G based on the net surface radiation Rn obtained in step 3.1);

[0014] 3.3) Calculate the atmospheric sensible heat flux H;

[0015] 3.4) Based on the net surface radiation Rn obtained in step 3.1), the soil heat flux G obtained in step 3.2), and the atmospheric sensible heat flux H obtained in step 3.3), the surface energy balance equation is constructed, and the instantaneous farmland evapotranspiration is obtained according to the surface energy balance equation.

[0016] Preferably, the expression for the net surface radiation Rn in step 3.1) is:

[0017] R n =(1-α)R s +Lin-Lout-(1-ε)Lin;

[0018] in:

[0019] R s It is instantaneous incident shortwave radiation, and the unit is W / m 2 ;

[0020] α is the surface albedo, which is dimensionless;

[0021] Lin is the instantaneous incident longwave radiation, and its unit is W / m. 2 ;

[0022] Lout is the instantaneous emitted long-wave radiation, measured in W / m.2 ;

[0023] ε is the surface emissivity, which is dimensionless.

[0024] Preferably, the instantaneous incident shortwave radiation R s The expression is:

[0025] R s =G sc cosθd r 2 τ sw

[0026] G sc It is the solar constant;

[0027] θ is the solar zenith angle. It represents the local geographical latitude; δ represents the solar declination.

[0028] d r The Earth-Sun distance correction factor is dimensionless.

[0029] τ sw It is atmospheric transmittance, dimensionless, τ sw =0.75 + 2 × 10 -5 Z, where Z is the elevation of farmland in the study area, in meters;

[0030] The expression for the instantaneous incident long-wave radiation Lin is:

[0031] Lin = ε atm σT a 4

[0032] ε atm It is the atmospheric emissivity, dimensionless, ε. atm =1.08(-lnτ) sw ) 0.265 σ is the Stefan-Boltzmann constant; T a The air temperature in step 1) is measured in Kelvin (K).

[0033] The expression for the instantaneous emitted long-wave radiation Lout is:

[0034] Lout=εσT s 4

[0035] ε is the surface emissivity, which is dimensionless, and ε = 1.009 + 0.047lnNDVI, where NDVI is the normalized vegetation index. NIR is the reflectance in the near-infrared band, and RED is the reflectance in the red band.

[0036] T s It is the surface temperature in thermal infrared image data, and the unit is K.

[0037] Preferably, the expression for the soil heat flux G in step 3.2) is:

[0038]

[0039] Preferably, the expression for the atmospheric sensible heat flux H in step 3.3) is:

[0040]

[0041] ρ a It is air density, and the unit is kg / m³. 3 ;

[0042] C p It is the specific heat capacity of dry air;

[0043] r ah It is the aerodynamic impedance for heat transfer, measured in s / m. z1 represents the vegetation canopy height, z2 represents a rough reference height; k is the von Kármán constant; u * It is the frictional wind speed, and the unit is m / s. u z This is the wind speed at altitude z, measured in m / s; z om NDVI is the pixel surface roughness, measured in meters (m). When NDVI < 0, the value is 0.001; when NDVI > 0, the value is...

[0044] ΔT is the near-surface air temperature gradient. It is assumed that ΔT has a linear relationship with the surface temperature Ts, i.e., ΔT = aT s +b, where a and b are parameters. b = ΔT hot -aT cold ΔT hot and ΔT cold These represent the near-surface temperature gradients at the hottest and coldest points, respectively, T hot and T cold These are the surface temperatures at the hottest and coldest points.

[0045] Preferably, step 3.4) is implemented as follows:

[0046] 3.4.1) Based on the principle of surface energy balance, a surface energy balance equation is constructed using the net surface radiation Rn obtained in step 3.1), the soil heat flux G obtained in step 3.2), and the atmospheric sensible heat flux H obtained in step 3.3). The expression of the surface energy balance equation is:

[0047] LE=R n -HG;

[0048] in:

[0049] LE is latent heat flux, measured in W / m³. 2 ;

[0050] 3.4.2) Calculate the instantaneous evapotranspiration based on the surface energy balance equation obtained in step 3.4.1). The calculation method is as follows:

[0051] LE = ET a λρ w

[0052] in:

[0053] ET a It is the instantaneous evaporation rate, and the unit is mm / s;

[0054] Λ is the latent heat of vaporization, measured in J / kg, and λ = [2.501 - 0.00236 × (Ts - 273.15)] × 10 6 ;

[0055] ρ w It is the density of water.

[0056] Preferably, step 4) is specifically implemented by using the evaporation ratio method or the sine function method to perform a daily-scale expansion of the instantaneous farmland evapotranspiration obtained in step 3) to complete the farmland-scale evapotranspiration estimation.

[0057] Beneficial Effects: The advantages of this invention are: This invention provides a method for estimating farmland-scale evapotranspiration, comprising the following steps: 1) acquiring basic data of the study area based on UAVs; the basic data includes multispectral image data, thermal infrared image data, and meteorological data; 2) preprocessing the multispectral image data and thermal infrared image data obtained in step 1); 3) obtaining instantaneous farmland evapotranspiration based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2); 4) extending the instantaneous farmland evapotranspiration obtained in step 3) to a daily scale to complete the estimation of farmland-scale evapotranspiration. This invention utilizes the surface energy balance equation to calculate evapotranspiration, which has strong mechanistic validity. However, the operation of this equation requires more than ten types of ground monitoring information, such as solar radiation, soil heat flux, and air temperature and humidity, which are difficult to meet with existing agricultural production monitoring systems. This invention uses drones equipped with sensors to acquire the information needed to estimate farmland evapotranspiration, overcoming the problems of insufficient distribution and coverage of ground monitoring stations. The flexible operating time of drones can respond promptly to changes in farmland moisture, and the decimeter-level spatial resolution can effectively adapt to the characteristics of farmland of various sizes, providing real-time and effective data support for precise irrigation decisions, improving the practicality of farmland evapotranspiration estimation, and ensuring high estimation accuracy. Attached Figure Description

[0058] Figure 1 This is a flowchart of the method for estimating farmland-scale evapotranspiration provided by the present invention;

[0059] Figure 2 This is a farmland evapotranspiration estimation diagram based on the farmland-scale evapotranspiration estimation method provided by the present invention. Detailed Implementation

[0060] The implementation example of this invention selects the Qianke Village irrigation area in Zhuqiao Irrigation District, Yizheng City as the research area. The main agricultural planting method in this area is rice-wheat rotation. Through the implementation of this technical solution, the evapotranspiration of farmland can be effectively estimated. The specific implementation steps are as follows:

[0061] Step 1: A drone equipped with a multispectral imager, a thermal infrared imager, and meteorological sensors was used to acquire basic data of the study area. The multispectral imager acquired multispectral images in the bands of 450nm (blue light band), 550nm (green light band), 660nm (red light band), 680nm (red edge band), and 840nm (near-infrared band). The thermal infrared imager mainly acquired ground temperature information, and the meteorological sensors acquired data including wind speed, temperature, and air pressure.

[0062] Step 2: Design the flight route for the study area and preprocess the captured multispectral and thermal infrared image data, specifically including image registration, image stitching, radiometric calibration, and resampling to 0.2m resolution in sequence.

[0063] Step 3: Calculate the net surface radiation Rn using the surface radiation balance equation R. n =(1-α)R s +Lin-Lout-(1-ε)Lin,R s Instantaneous incident shortwave radiation (W / m 2 ); α is the surface albedo (dimensionless); Lin is the instantaneous incident longwave radiation (W / m²). 2 Lout represents the instantaneous emitted long-wave radiation (W / m²). 2 ); ε is the surface emissivity (dimensionless).

[0064] The instantaneous incident shortwave radiation R s =G sc cosθd r 2 τ sw G sc The solar constant has a value of 1367 W / m 2 θ is the solar zenith angle. δ represents the local geographical latitude; in this example, the latitude of the center point of the study area is 32.5°. δ represents the solar declination, which is the angle between the Earth's equatorial plane and the line connecting the Sun and the Earth's center, and can be calculated using an approximate formula: n is the day number, which is 108 on the day of the experiment; ω is the hour angle, ω = 15° (solar time - 12). The standard time in this example is 12:00 noon, the time of the UAV flight center. The longitude of the test area center is 119.2°, and the longitude of the center of the East 8th time zone is 120°. The calculated hour angle is approximately 0. r τ is the Earth-Sun distance correction factor (dimensionless), with a value of 1.015; sw It is atmospheric transmittance (dimensionless), τ sw =0.75 + 2 × 10 -5 Z represents the elevation (m) of farmland in the study area. The digital elevation model (DEM) of the study area was downloaded from 91 Satellite Map and its value ranges from 17.6 to 31.7 m.

[0065] Where the instantaneous incident long-wave radiation Lin = ε atm σT a 4 ;ε atm ε is the atmospheric emissivity (dimensionless). atm =1.08(-lnτ) sw ) 0.265 σ is the Stefan-Boltzmann constant, with a value of 5.67 × 10⁻⁶. -8 W / m 2 ·K 4Ta represents the air temperature (K), obtained from a meteorological sensor, with a value ranging from 302.7 to 304.9 K.

[0066] Where the instantaneously emitted long-wave radiation Lout=εσT s 4 ε is the surface emissivity (dimensionless), ε = 1.009 + 0.047lnNDVI, where NDVI is the normalized difference vegetation index. NIR represents the near-infrared reflectance, and RED represents the red reflectance, corresponding to the reflectance at 840nm and 660nm respectively for a multispectral imager. NDVI ranges from -0.95 to 0.89, where ε is assumed to be 1 when NDVI ≤ 0. s The surface temperature (K) was obtained from a thermal infrared imager and its value ranges from 301.26 to 307.11 K.

[0067] The calculated distribution of Rn is 667.93–897.55 W / m. 2 .

[0068] Step 4: Calculate soil heat flux G. Introducing multispectral data to estimate soil heat flux can reflect the differences in underlying surface for different crops.

[0069] The calculated distribution of G is 13.57–94.03 W / m. 2 .

[0070] Step 5: Calculate the atmospheric sensible heat flux H. ρ a air density (kg / m³) 3 ), calculate based on the ideal gas equation P is the pressure (Pa), obtained from a meteorological sensor, with a value of 100315.6 Pa; M is the average molar mass of air, taken as 0.02896 kg / mol; and R is the ideal gas constant, taken as 8.314 J·mol⁻¹. -1 ·K -1 T is taken as the average air temperature of 303.85K; C p It is the specific heat capacity of dry air, taken as 1004 J·kg -1 ·K -1 ;r ah For heat transfer aerodynamic impedance (s / m), z1 is taken as the vegetation canopy height of 0.9m, z2 is taken as a rough reference height, and in this example, it is taken as the flight altitude of the UAV when acquiring meteorological data, which is 30m. k is the von Kármán constant with a value of 0.41; u * Frictional wind speed (m / s) u zLet z be the wind speed (m / s). In this example, the wind speed at the drone's flight altitude of 30m is 4.85m / s. om Let represent the pixel surface roughness (m). When NDVI ≤ 0, the value is 0.001; when NDVI > 0, ... ΔT is the near-surface air temperature gradient. It is assumed that ΔT has a linear relationship with the surface temperature Ts, i.e., ΔT = aT s +b, where a and b are parameters. b = ΔT hot -aT cold The extreme points are determined by selecting two extreme pixels on the thermal infrared image, namely the hot spot and the cold spot. At the hot spot, it is assumed that all available energy can be used for surface heating, i.e., crop evaporation and transpiration are approximately zero. Therefore, H = R. n -G, assuming that at extremely cold points all available energy can be used for crop evaporation and transpiration, then H = 0, and the value of ΔT is obtained based on this assumption. ΔT hot and ΔT cold These represent the near-surface temperature gradients at the hottest and coldest points, respectively, T hot and T cold The surface temperatures at extreme hot and cold points were obtained using an unmanned thermal infrared imager and meteorological sensors. (T) hot The value is 301.26K, T cold The value is 307.11K, ΔT hot The value is 4.11K, ΔT cold The value is 1.44K.

[0071] The calculated distribution of H is 116.00–445.62 W / m. 2 .

[0072] Step Six: Based on the principle of surface energy balance, calculate the instantaneous data of each flux through the above steps, and then obtain the instantaneous value of evapotranspiration based on the relationship between flux and evapotranspiration. LE = R n -HG=ET a λρ w LE is the latent heat flux (W / m³) 2 ); ET a λ is the instantaneous evaporation rate (mm / s); λ is the latent heat of vaporization (J / kg), λ=[2.501-0.00236×(Ts-273.15)]×10 6 ;ρ w The density of water (1000 kg / m³) 3 ).

[0073] The calculated LE distribution is 122.96–410.46 W / m². 2 .

[0074] Step Seven: Farmland evapotranspiration is influenced by various factors throughout the day, including solar radiation, temperature, humidity, and wind speed, exhibiting a clear diurnal variation. Extending the instantaneous farmland evapotranspiration estimation results to a daily scale yields spatial distribution information of daily farmland evapotranspiration, facilitating comparison with measured evapotranspiration data and providing practical significance for water resource management and agricultural irrigation guidance. The instantaneous evapotranspiration calculated in Step Six is ​​extended to a daily scale using the sine function method. N ε The daily evapotranspiration time is calculated as the day's sunshine duration minus 2 (h), with a value of 11h. t is the time interval from sunrise to the drone's passage (h), with a value of 7h.

[0075] Step 8: Using the daily-scale expanded farmland evapotranspiration estimation results and geographic information data, generate a regional evapotranspiration map to obtain spatial distribution information of daily farmland evapotranspiration. Combined with irrigation and drainage canal distribution and control zone vector data, precise irrigation can be guided by zone, thereby improving water resource utilization efficiency. The specific steps of the farmland-scale evapotranspiration estimation method based on UAV multi-source remote sensing proposed in this invention are detailed below. Figure 1 .

[0076] Based on the method provided by this invention, the estimated daily evapotranspiration of farmland in the Qianke Village irrigation area of ​​Zhuqiao Irrigation District, Yizheng City, on April 18th is shown in the figure. Figure 2 The evaporation rate is 1.4–4.7 mm / d.

Claims

1. A method for estimating farmland-scale evapotranspiration, characterized in that: Includes the following steps: 1) Acquire basic data of the research area using UAVs; the basic data includes multispectral image data, thermal infrared image data, and meteorological data; 2) Preprocess the multispectral image data and thermal infrared image data obtained in step 1); 3) Based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2), obtain the instantaneous farmland evapotranspiration. 4) Extend the instantaneous farmland evapotranspiration obtained in step 3) to a daily scale to complete the farmland-scale evapotranspiration estimation; The specific implementation method in step 3) is as follows: 3.1) Calculate the net surface radiation Rn based on the meteorological data in step 1) and the multispectral image data and thermal infrared image data preprocessed in step 2); 3.2) Calculate the soil heat flux G based on the net surface radiation Rn obtained in step 3.1); 3.3) Calculate the atmospheric sensible heat flux H; 3.4) Based on the net surface radiation Rn obtained in step 3.1), the soil heat flux G obtained in step 3.2), and the atmospheric sensible heat flux H obtained in step 3.3), the surface energy balance equation is constructed, and the instantaneous farmland evapotranspiration is obtained according to the surface energy balance equation. The expression for the net surface radiation Rn in step 3.1) is: in: It is instantaneous incident shortwave radiation, and the unit is 1000 kilometres per minute (kJ / m ; Surface albedo, dimensionless; It is instantaneous incident longwave radiation, and the unit is . ; It is the instantaneous emission of long-wave radiation, and the unit is _____. ; It is the surface emissivity, which is dimensionless. The instantaneous incident shortwave radiation The expression is: , It is the solar constant; It is the zenith angle of the sun. , It is the local geographical latitude; Solar declination; The Earth-Sun distance correction factor is dimensionless. It is atmospheric transmittance, dimensionless. , where Z is the elevation of farmland in the study area, in meters; The expression for the instantaneous incident long-wave radiation Lin is: , It is atmospheric emissivity, dimensionless. σ is the Stefan-Boltzmann constant; T a The air temperature in step 1) is measured in Kelvin (K). The expression for the instantaneous emitted long-wave radiation Lout is: , It is the surface emissivity, dimensionless. NDVI is the Normalized Difference Vegetation Index; NIR is the reflectance in the near-infrared band, and RED is the reflectance in the red band. It is the surface temperature in thermal infrared image data, and the unit is K.

2. The method for estimating farmland-scale evapotranspiration according to claim 1, characterized in that: The multispectral image data mentioned in step 1) is acquired by a multispectral imager mounted on a drone. The multispectral image data includes blue light band image data, green light band image data, red light band image data, red edge band image data, and near-infrared band image data. The thermal infrared image data is acquired by a thermal infrared imager mounted on a drone. The thermal infrared image data includes the ground surface temperature. The meteorological data is acquired by a meteorological sensor mounted on a drone. The meteorological data includes wind speed, air temperature, and air pressure.

3. The method for estimating farmland-scale evapotranspiration according to claim 2, characterized in that: The preprocessing method in step 2) is image registration, image stitching, and radiometric calibration.

4. The method for estimating farmland-scale evapotranspiration according to claim 1, characterized in that: The expression for the soil heat flux G in step 3.2) is: .

5. The method for estimating farmland-scale evapotranspiration according to claim 4, characterized in that: The expression for the atmospheric sensible heat flux H in step 3.3) is: , It is air density, and the unit is kg / m³. 3 ; It is the specific heat capacity of dry air; It is the aerodynamic impedance for heat transfer, measured in s / m. z1 represents the height of the vegetation canopy, z2 represents a rough reference height; k is the von Kármán constant. It is the frictional wind speed, and the unit is m / s. , This is the wind speed at altitude z, measured in m / s. NDVI is the pixel surface roughness, measured in meters (m). When NDVI < 0, the value is 0.001; when NDVI > 0, the value is... , It is the near-surface air temperature gradient, assuming It has a linear relationship with surface temperature Ts, that is a and b are parameters. , , and These represent the near-surface temperature gradients at the hottest and coldest points, respectively. and These are the surface temperatures at the hottest and coldest points.

6. The method for estimating farmland-scale evapotranspiration according to claim 5, characterized in that: The specific implementation method of step 3.4) is as follows: 3.4.1) Based on the principle of surface energy balance, a surface energy balance equation is constructed using the net surface radiation Rn obtained in step 3.1), the soil heat flux G obtained in step 3.2), and the atmospheric sensible heat flux H obtained in step 3.3). The expression of the surface energy balance equation is: in: LE stands for latent heat flux, and its unit is... ; 3.4.2) Calculate the instantaneous evapotranspiration based on the surface energy balance equation obtained in step 3.4.1). The calculation method is as follows: ,in: It is the instantaneous evaporation rate, and the unit is mm / s; Λ is the latent heat of vaporization, measured in J / kg. ; It is the density of water.

7. The method for estimating farmland-scale evapotranspiration according to any one of claims 1-6, characterized in that: The specific implementation of step 4) is to use the evaporation ratio method or the sine function method to expand the instantaneous farmland evapotranspiration obtained in step 3) on a daily scale, and complete the farmland-scale evapotranspiration estimation.

Images