Farmland scale evapotranspiration estimation method

Abstract

The invention belongs to the technical field of agricultural remote sensing, and relates to a farmland scale evapotranspiration estimation method, which comprises the following steps: 1) acquiring basic data of a research area based on an unmanned aerial vehicle; the basic data comprises multispectral image data, thermal infrared image data and meteorological data; 2) preprocessing the multispectral image data and the thermal infrared image data obtained in the step 1); (3) acquiring instantaneous farmland evapotranspiration based on the meteorological data in the step (1) and the multispectral image data and the thermal infrared image data preprocessed in the step (2); 4, the instantaneous farmland evapotranspiration obtained in the step 3 is expanded to the daily scale, and farmland scale evapotranspiration estimation is accomplished. The farmland scale evapotranspiration estimation method can effectively estimate the farmland scale evapotranspiration, and meanwhile high estimation precision can be guaranteed.

Description

Technical Field

[0001] The present invention belongs to the field of agricultural remote sensing technology and relates to a method for estimating farmland-scale evapotranspiration, in particular to a method for estimating farmland-scale evapotranspiration based on multi-source remote sensing using unmanned aerial vehicles (UAVs). Background Art

[0002] Farmland evapotranspiration (ET) is the total water vapor flux from crops and the soil they adhere to to the atmosphere. It is essential fundamental information for agricultural water management. Determining its value can reduce ineffective water consumption during crop growth, provide data support for formulating scientific irrigation schedules, guiding precise irrigation decisions, and improving water resource utilization efficiency. Therefore, it is imperative to develop monitoring methods suitable for farmland ET. Traditional ground-based methods for observing farmland ET mainly rely on eddy covariance, ripple ratio, and lysimeters. While these methods provide relatively accurate ET measurements, they are subject to high equipment costs, difficulty in ongoing maintenance, and limited universal applicability. Currently available satellite imagery for estimating ground-based ET has a maximum temporal and spatial resolution of 8 days and 30 meters, making it difficult to apply to farmland-scale ET. Furthermore, the estimation process requires data from ground-based meteorological stations, making it difficult to apply. In contrast, drones offer flexible operation times and can achieve decimeter-level spatial resolution. Using drone-mounted sensors to obtain the information necessary for ET estimation, combined with the surface energy balance equation, the spatial distribution of farmland ET can be rapidly derived. Summary of the Invention

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

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

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

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

[0007] 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);

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

[0009] Preferably, the multispectral image data in step 1) is collected by a multispectral imager carried on an unmanned aerial vehicle, 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 collected by a thermal infrared imager carried on an unmanned aerial vehicle, and the thermal infrared image data includes surface temperature; the meteorological data is collected by a meteorological sensor carried on an unmanned aerial vehicle, and the meteorological data includes wind speed, air temperature and air pressure.

[0010] Preferably, the preprocessing method 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, and exemplarily resampling to 0.2m resolution.

[0011] Preferably, the specific implementation of step 3) is:

[0012] 3.1) Calculating the surface net 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 surface net radiation Rn obtained in step 3.1);

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

[0015] 3.4) Based on the surface net 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), a surface energy balance equation is constructed. The instantaneous farmland evapotranspiration is obtained according to the surface energy balance equation.

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

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

[0018] in:

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

[0020] α is the surface albedo, dimensionless;

[0021] Lin is the instantaneous incident longwave radiation, in W / m 2 ;

[0022] Lout is the instantaneous outgoing longwave radiation, in W / m2 ;

[0023] ε is the surface emissivity, 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 is the solar constant;

[0027] θ is the solar zenith angle, is the local geographic latitude; δ is the solar declination;

[0028] d r is the correction factor for the distance between the Sun and the Earth, dimensionless;

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

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

[0031] Lin=ε atm σT a 4

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

[0033] The expression of the instantaneous outgoing long-wave radiation Lout is:

[0034] Lout=εσT s 4

[0035] ε is the surface emissivity, dimensionless, ε = 1.009 + 0.047 ln NDVI, NDVI is the normalized difference vegetation index; NIR is the reflectance of the near infrared band, and RED is the reflectance of the red band;

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

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

[0038]

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

[0040]

[0041] ρ a is the air density in kg / m 3 ;

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

[0043] r ah is the aerodynamic resistance to heat transfer, in s / m, The value of z1 is the height of the vegetation canopy, and the value of z2 is the reference height; k is the von Karman constant; u * is the friction wind speed, in m / s, u z is the wind speed at height z, in m / s; om It is the surface roughness of the pixel, the unit is m, when NDVI < 0, the value is 0.001, when NDVI > 0,

[0044] ΔT is the near-surface air temperature gradient. Assume that ΔT is linearly related to the surface temperature Ts, that is, ΔT = aT s +b, a, b are parameters, b=ΔT hot -aT cold , ΔT hot and ΔT cold are the near-surface temperature gradients at the hottest and coldest points, respectively, T hot and T cold The surface temperatures of the hottest and coldest points.

[0045] Preferably, the specific implementation of step 3.4) is:

[0046] 3.4.1) Based on the surface energy balance principle, the surface energy balance equation is constructed using the surface net 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 the latent heat flux, in W / m 2 ;

[0050] 3.4.2) Calculate instantaneous evapotranspiration based on the surface energy balance equation obtained in step 3.4.1), the calculation method being:

[0051] LE=ET a λρ w

[0052] in:

[0053] ET a is the instantaneous evapotranspiration, in mm / s;

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

[0055] ρ w is the density of water.

[0056] Preferably, the specific implementation method of step 4) is to use the evaporation ratio method or the sine function method to expand the instantaneous farmland evaporation obtained in step 3) to a daily scale to complete the farmland scale evaporation estimation.

[0057] Beneficial effects: The advantages of the present invention are: the present invention provides a method for estimating evapotranspiration at the farmland scale, comprising the following steps: 1) obtaining basic data of the study area based on an unmanned aerial vehicle; 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) expanding the instantaneous farmland evapotranspiration obtained in step 3) to a daily scale to complete the estimation of farmland-scale evapotranspiration. The present invention utilizes the surface energy balance equation to calculate evapotranspiration, which has a strong mechanistic rationality. However, the operation of this equation requires more than ten types of ground monitoring information, such as solar radiation, soil heat flux, air temperature and humidity, which are difficult to meet in the existing agricultural production monitoring system. The present invention uses drones equipped with sensors to obtain the information needed to estimate farmland evapotranspiration, overcoming the problems of insufficient distribution and coverage of ground monitoring stations. The flexible operation time of drones can promptly respond to changes in farmland moisture, and the decimeter-level spatial resolution can effectively adapt to the characteristics of farmland of various specifications, providing real-time and effective data support for precise irrigation decision-making, improving the practicality of farmland evapotranspiration estimation, and ensuring a high estimation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 It is a flow chart of the method for estimating evapotranspiration at the farmland scale provided by the present invention;

[0059] Figure 2 This is a farmland evapotranspiration estimation diagram of an embodiment based on the farmland scale evapotranspiration estimation method provided by the present invention. DETAILED DESCRIPTION

[0060] The implementation case of the present invention selects the Qianke Village irrigation area of ​​​​Yizheng City's Zhuqiao Irrigation District 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: Use a drone equipped with a multispectral imager, a thermal infrared imager, and a meteorological sensor to obtain basic data for the study area. The multispectral imager acquires multispectral images at wavelengths including 450nm (blue light), 550nm (green light), 660nm (red light), 680nm (red edge), and 840nm (near infrared). The thermal infrared imager primarily acquires ground temperature information, while the meteorological sensor acquires information including wind speed, temperature, and air pressure.

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

[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 is the instantaneous incident shortwave radiation (W / m 2 ); α is the surface albedo (dimensionless); Lin is the instantaneous incident longwave radiation (W / m 2 ); Lout is the instantaneous outgoing 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 is the solar constant, which is 1367W / m 2 ; θ is the solar zenith angle, is the local geographic latitude. In this example, the latitude of the center of the study area is 32.5°. δ is the solar declination, which is the angle between the Earth's equatorial plane and the line connecting the Sun and the Earth's center. It can be calculated using the 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 noon at the drone flight center time. The central longitude of the test area is 119.2°, and the central longitude of the East 8th District is 120°. The calculated hour angle is approximately 0; d r is the correction factor for the distance between the Sun and the Earth (dimensionless), which is 1.015; τ sw is the atmospheric transmittance (dimensionless), τ sw =0.75+2×10 -5 Z is the altitude of farmland in the study area (m). The digital elevation DEM of the study area was downloaded from 91 satellite images, and its value ranged from 17.6 to 31.7 m.

[0065] 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, which is 5.67×10 -8 W / m 2 ·K 4; Ta is the air temperature (K), obtained by the meteorological sensor, and its value is 302.7~304.9K.

[0066] The instantaneous outgoing long-wave radiation Lout=εσT s 4 ; ε is the surface emissivity (dimensionless), ε=1.009+0.047lnNDVI, NDVI is the normalized difference vegetation index, NIR is the near-infrared band reflectance, RED is the red band reflectance, corresponding to the reflectance of 840nm and 660nm of the multispectral imager, respectively. The distribution of NDVI is -0.95 to 0.89. When NDVI ≤ 0, ε = 1 is assumed; T s is the surface temperature (K), obtained by a thermal infrared imager, and its value ranges from 301.26 to 307.11K.

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

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

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

[0070] Step 5: Calculate the atmospheric sensible heat flux H, ρ a is the air density (kg / m 3 ), calculated according to the ideal gas equation P is the pressure (Pa), obtained from the meteorological sensor, and its value is 100315.6 Pa. M is the average molar mass of air, which is 0.02896 kg / mol. R is the ideal gas constant, which is 8.314 J·mol. -1 ·K -1 , T takes the average air temperature as 303.85K; C p is the specific heat capacity of dry air, which is 1004 J·kg -1 ·K -1 ; r ah is the aerodynamic resistance of heat transfer (s / m), The value of z1 is 0.9m, the value of z2 is the reference height. In this example, the value is 30m, the flight height of the drone when acquiring meteorological data. k is the von Karman constant, which is 0.41. u * is the friction wind speed (m / s), u zis the wind speed at the z height (m / s). In this example, the wind speed at the drone's flight altitude of 30m is 4.85m / s; om is the surface roughness of the pixel (m), when NDVI≤0, the value is 0.001, when NDVI>0, ΔT is the near-surface air temperature gradient. Assume that ΔT is linearly related to the surface temperature Ts, that is, ΔT = aT s +b, a, b are parameters, b=ΔT hot -aT cold , is determined by selecting two extreme pixels on the thermal infrared image, namely the hottest point and the coldest point. At the hottest point, it is assumed that all available energy is used for surface heating, that is, the crop evapotranspiration is approximately 0, then H = R n -G, assuming that the available energy is completely used for crop evaporation and transpiration at the extreme cold point, then H = 0, and the ΔT value is obtained based on this assumption. ΔT hot and ΔT cold are the near-surface temperature gradients at the hottest and coldest points, T hot and T cold is the surface temperature at the hottest and coldest points, obtained by unmanned thermal infrared imagers 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.62W / m 2 .

[0072] Step 6: Based on the surface energy balance principle, 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. 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 evaporation (J / kg), λ = [2.501-0.00236×(Ts-273.15)]×10 6 ρ w is the density of water (1000kg / m 3 ).

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

[0074] Step 7: The evapotranspiration process of farmland will be affected by multiple factors such as solar radiation, temperature, humidity, wind speed, etc. in a day, showing obvious diurnal variation characteristics. The instantaneous evapotranspiration estimation results of farmland are expanded to the daily scale to obtain the spatial distribution information of daily evapotranspiration of farmland, which is convenient for comparison with the measured evapotranspiration data and has more practical significance for water resources management and agricultural irrigation guidance. The sine function method is used to expand the instantaneous evapotranspiration calculated in step 6 to the daily scale, and the daily evapotranspiration is N ε is the daily evapotranspiration hours, specifically the sunshine duration of the day minus 2 (h), which is 11h; t is the time interval from sunrise to the UAV passing (h), which is 7h.

[0075] Step 8: Combine the daily scaled farmland evapotranspiration estimation results with geographic information data to achieve regional evapotranspiration mapping output, and obtain the spatial distribution information of the farmland daily evapotranspiration. Combined with the irrigation and drainage canal distribution and control zone vector data, it can guide the zone-by-zone precision irrigation, thereby improving the efficiency of water resource utilization. The specific steps of the farmland scale evapotranspiration estimation method based on multi-source remote sensing of unmanned aerial vehicles proposed in this invention are shown in Figure 1 .

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

Claims

1. A method for estimating evapotranspiration at the farmland scale, characterized by: The following steps are involved: 1) Obtain basic data of the study area using drones; 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) Expand the instantaneous farmland evapotranspiration obtained in step 3) to the daily scale to complete the farmland-scale evapotranspiration estimation.

2. The method for estimating farmland-scale evapotranspiration according to claim 1, wherein: The multispectral image data in step 1) is collected by a multispectral imager mounted on an unmanned aerial vehicle, 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 collected by a thermal infrared imager mounted on an unmanned aerial vehicle, and the thermal infrared image data includes surface temperature; the meteorological data is collected by a meteorological sensor mounted on an unmanned aerial vehicle, and the meteorological data includes wind speed, air temperature, and air pressure.

3. The method for estimating farmland-scale evapotranspiration according to claim 2, wherein: The pre-processing methods in step 2) are image registration, image stitching and radiometric calibration.

4. The method for estimating farmland-scale evapotranspiration according to claim 3, wherein: The specific implementation method in step 3) is: 3.1) Calculating the surface net 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 surface net radiation Rn obtained in step 3.1); 3.3) Calculate the atmospheric sensible heat flux H; 3.4) Based on the surface net 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), a surface energy balance equation is constructed. The instantaneous farmland evapotranspiration is obtained according to the surface energy balance equation.

5. The method for estimating farmland-scale evapotranspiration according to claim 4, wherein: The expression of surface net radiation Rn in step 3.1) is: R n =(1-α)R s +Lin-Lout-(1-e)Lin; in: R s is the instantaneous incident shortwave radiation, in W / m 2 ; α is the surface albedo, dimensionless; Lin is the instantaneous incident longwave radiation, in W / m 2 ; Lout is the instantaneous outgoing longwave radiation, in W / m 2 ; ε is the surface emissivity, dimensionless.

6. The method for estimating farmland-scale evapotranspiration according to claim 5, wherein: The instantaneous incident shortwave radiation R s The expression is: R s =G sc cosθd r 2 t sw G sc is the solar constant; θ is the solar zenith angle, is the local geographical latitude; δ is the solar declination; d r is the correction factor for the distance between the Sun and the Earth, dimensionless; τ sw is the atmospheric transmittance, dimensionless, τ sw =0.75+2×10 -5 Z, where Z is the altitude of the farmland in the study area, in m; The expression of the instantaneous incident long-wave radiation Lin is: Lin=e atm σT a 4 ε atm is the atmospheric emissivity, dimensionless, ε atm =1.08(-lnτ sw ) 0.265 ; σ is the Stefan-Boltzmann constant; T a The air temperature of the meteorological data in step 1) is in K; The expression of the instantaneous outgoing long-wave radiation Lout is: Lout=εσT s 4 ε is the surface emissivity, dimensionless, ε = 1.009 + 0.047 ln NDVI, NDVI is the normalized difference vegetation index; NIR is the reflectance of the near infrared band, and RED is the reflectance of the red band; T s It is the surface temperature in thermal infrared image data, in K.

7. The method for estimating farmland-scale evapotranspiration according to claim 6, wherein: The expression of soil heat flux G in step 3.2) is:

8. The method for estimating farmland-scale evapotranspiration according to claim 7, wherein: The expression of atmospheric sensible heat flux H in step 3.3) is: ρ a is the air density in kg / m 3 ; C p is the specific heat capacity of dry air; r ah is the aerodynamic resistance to heat transfer, in s / m, The value of z1 is the height of the vegetation canopy, and the value of z2 is the reference height; k is the von Karman constant; u * is the friction wind speed, in m / s, u z is the wind speed at height z, in m / s; om It is the surface roughness of the pixel, the unit is m, when NDVI < 0, the value is 0.001, when NDVI > 0, ΔT is the near-surface air temperature gradient. Assume that ΔT is linearly related to the surface temperature Ts, that is, ΔT = aT s +b, a, b are parameters, b=ΔT hot -aT cold , ΔT hot and ΔT cold are the near-surface temperature gradients at the hottest and coldest points, respectively, T hot and T cold The surface temperatures of the hottest and coldest points.

9. The method for estimating farmland-scale evapotranspiration according to claim 8, wherein: The specific implementation of step 3.4) is: 3.4.1) Based on the surface energy balance principle, the surface energy balance equation is constructed using the surface net 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: THE R n -HG; in: LE is the latent heat flux, in W / m 2 ; 3.4.2) Calculate instantaneous evapotranspiration based on the surface energy balance equation obtained in step 3.4.1), the calculation method being: THE AND a λρ w in: ET a is the instantaneous evapotranspiration, in mm / s; Λ is the latent heat of vaporization, in J / kg, λ = [2.501-0.00236×(Ts-273.15)]×10 6 ; ρ w is the density of water.

10. The method for estimating farmland-scale evapotranspiration according to any one of claims 1 to 9, characterized in that: The specific implementation method 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) to a daily scale to complete the farmland scale evapotranspiration estimation.

Images