Irrigation prescription map inversion method based on unmanned aerial vehicle spectrum data

Abstract

The invention relates to the technical field of image processing, in particular to an irrigation prescription map inversion method based on unmanned aerial vehicle spectrum data. The method comprises the following steps: carrying out unmanned aerial vehicle data acquisition on a farmland area to obtain a multispectral image, DEM data and meteorological interpolation data; performing micro-landform parameter calculation according to the multispectral image, DEM data and unmanned aerial vehicle GPS data, and constructing a terrain-spectrum coupling model to obtain a terrain relief degree map and the terrain-spectrum coupling model; performing vegetation multi-dimensional state analysis by using a terrain-spectrum coupling model to obtain a vegetation index map, a canopy temperature map and a water stress preliminary judgment map; performing net radiation calculation by using the meteorological interpolation data to obtain a net radiation diagram; according to the net radiation diagram, the vegetation index diagram and the topographic relief degree diagram, soil heat flux estimation is conducted, and a soil heat flux diagram is obtained. According to the method, collaborative inversion is carried out by eliminating the terrain effect, and the precision and reliability of irrigation prescription map inversion are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and particularly to a method for inverting an irrigation prescription map based on unmanned aerial vehicle (UAV) spectral data. Background Art

[0002] The development of satellite remote sensing, aerial remote sensing, and UAV remote sensing technologies has provided new means for obtaining farmland information. The traditional extensive irrigation method causes waste of water resources, and precision agriculture requires fine irrigation according to the actual water demand of crops.

[0003] The undulating terrain leads to uneven distribution of solar radiation on different slopes and aspects, thus affecting the energy received by the ground surface. Vegetation on different slopes and aspects is also affected by the terrain in terms of its growth status and spectral characteristics. Traditional remote sensing inversion methods usually assume that the ground surface is horizontal, ignoring the influence of the terrain, resulting in reduced inversion accuracy in areas with complex terrain (such as mountains and hills). The terrain effect will affect the calculation of vegetation indices, the inversion of temperature, the estimation of parameters of the energy balance model, etc., ultimately affecting the accuracy of the irrigation prescription map.

[0004] Traditional remote sensing inversion methods usually rely only on a single type of remote sensing data, such as multispectral data or thermal infrared data. Multispectral data mainly reflects the growth status and coverage of vegetation, but is not sensitive to soil moisture information. Thermal infrared data mainly reflects the surface temperature, but is easily affected by factors such as the atmosphere and vegetation coverage. The information provided by a single data source is limited and it is difficult to comprehensively and accurately reflect the water shortage status of crops. Summary of the Invention

[0005] Based on this, it is necessary to provide a method for inverting an irrigation prescription map based on UAV spectral data to solve at least one of the above technical problems.

[0006] To achieve the above object, a method for inverting an irrigation prescription map based on UAV spectral data includes the following steps:

[0007] Step S1: Collect UAV data for the farmland area, and perform multi-source data registration and fusion to obtain a fusion data set, where the fusion data set includes a multispectral image, a thermal infrared image, DEM data, and meteorological interpolation data;

[0008] Step S2: Calculate microtopography parameters based on the multispectral image, DEM data, and UAV GPS data, and construct a terrain-spectral coupling model to obtain a terrain undulation map and a terrain-spectral coupling model; use the terrain-spectral coupling model to perform multi-dimensional state analysis of vegetation to obtain a vegetation index map, a canopy temperature map, and a preliminary water stress judgment map;

[0009] Step S3: Use meteorological interpolation data to calculate net radiation and obtain a net radiation map; estimate soil heat flux based on the net radiation map, vegetation index map, and topographic relief map to obtain a soil heat flux map; calculate apparent thermal inertia based on the thermal infrared image and canopy temperature map, and estimate sensible heat flux to obtain an apparent thermal inertia map and a sensible heat flux map; perform evapotranspiration inversion based on the soil heat flux map, apparent thermal inertia map, and sensible heat flux map to obtain an evapotranspiration distribution map;

[0010] Step S4: Extract sample characteristics based on the initial water stress judgment map and evapotranspiration distribution map, and perform water shortage judgment to obtain a phenological feature dataset and a water shortage diagnosis index;

[0011] Step S5: Divide the irrigation management area based on the phenological feature data and water shortage diagnosis index, and calculate the irrigation amount to obtain an irrigation prescription map.

[0012] The present invention obtains high-resolution and multi-source farmland data through UAV remote sensing and meteorological observations, and generates a fusion dataset containing multi-spectral, thermal infrared, topographic, and meteorological information through a series of preprocessing and registration fusion operations, providing a comprehensive, accurate, and consistent data basis for all subsequent analyses. Micro-topographic parameters are extracted using DEM data, a topographic-spectral coupling model is constructed, and based on this model, multi-spectral and thermal infrared data are corrected, effectively eliminating the influence of topography on key parameters such as vegetation indices and canopy temperature, improving the inversion accuracy of these parameters in areas with complex topography, and laying a foundation for subsequent energy balance analysis and water shortage diagnosis. Based on the energy balance principle, net radiation, soil heat flux, sensible heat flux, and latent heat flux are inverted using multi-source data, and apparent thermal inertia is innovatively introduced to assist in the inversion of evapotranspiration, finally obtaining a high-precision evapotranspiration distribution map. Evapotranspiration is the main way of crop water consumption, and the results of this step provide key information for accurately evaluating crop water requirements. By extracting multi-dimensional features, a water shortage diagnosis model is constructed using machine learning methods, and a water shortage diagnosis index (WDI) is generated. WDI can quantitatively reflect the severity of crop water shortage, realizing the transformation from qualitative evaluation to quantitative diagnosis, and providing a more refined and accurate basis for irrigation decision-making. According to WDI and meteorological data, irrigation zoning is carried out, and combined with potential evapotranspiration, crop coefficient, actual evapotranspiration, and soil water holding capacity correction coefficient, the irrigation amount of each zone is calculated in detail, and finally an irrigation prescription map is generated. This prescription map can guide the precise irrigation system to achieve on-demand and quantitative irrigation, which is the core result of the entire method and has significant water-saving potential. Therefore, the present invention provides a method for inverting an irrigation prescription map based on UAV spectral data, eliminates the topographic effect by constructing a fine topographic-spectral coupling model, and realizes the collaborative inversion of multi-spectral, thermal infrared, DEM, and meteorological data through multi-source data fusion, solving the disadvantages of existing methods in terms of landform and data sources, and improving the accuracy and reliability of the inversion of the irrigation prescription map.

[0013] Preferably, step S1 includes the following steps:

[0014] Step S11: Uniformly arrange at least 10 ground control points in the farmland, measure the coordinates of the control points using an RTK GNSS receiver, with horizontal and vertical accuracies both better than 2 cm. Plan the UAV flight path to fly in an east-west strip shape, set the flight altitude to 100 meters, the forward overlap to 80%, and the side overlap to 70%. Collect UAV data for the farmland area through the UAV and synchronously obtain meteorological data to obtain original spectral data, original digital elevation model data, UAV GPS data, and original weather station data;

[0015] Step S12: Perform spectral data preprocessing on the UAV spectral data to obtain multi-spectral images and thermal infrared images;

[0016] Step S13: Perform DEM data preprocessing on the digital elevation model data to obtain DEM data;

[0017] Step S14: Perform meteorological data interpolation processing on the original meteorological station data to obtain meteorological interpolation data.

[0018] The present invention obtains high-resolution, multi-spectral / thermal infrared farmland images, providing a data basis for refined agricultural monitoring; synchronizing GPS and ground control point data ensures the geometric accuracy of the images; synchronizing meteorological data provides necessary environmental inputs for subsequent energy balance models and crop growth models, improving the model accuracy. Through radiometric calibration, geometric correction, atmospheric correction (thermal infrared), and image registration, the effects of sensors, geometric distortion, and the atmosphere are eliminated, obtaining multi-spectral and thermal infrared images with physical significance, spatially consistent, and accurately registered, laying a foundation for subsequent quantitative analysis and multi-source data fusion. Through denoising, interpolation, smoothing, and coordinate transformation, the quality and accuracy of the DEM data are improved, and the spatial consistency between the DEM data and the remote sensing images is ensured, providing reliable data for subsequent terrain analysis. Converting discrete meteorological station data into continuous spatially distributed data to match the spatial resolution of the remote sensing images provides spatially distributed meteorological input parameters for energy balance calculations and crop growth models. Accurately registering and fusing multi-spectral, thermal infrared, DEM, and meteorological data into a multi-dimensional data cube realizes the effective integration of multi-source data, providing comprehensive and consistent data support for subsequent comprehensive analysis, which is the key to improving the accuracy of irrigation prescription maps.

[0019] Preferably, the calculation of the micro-topographic parameters in step S2 and the construction of the terrain-spectral coupling model are specifically as follows:

[0020] Perform micro-topographic parameter calculation on the DEM data. According to the elevation differences in the horizontal and vertical directions, calculate the angle between the slope normal vector and the due north direction, in degrees, with a range of 0 - 360 degrees, where 0 degrees represents due north, 90 degrees represents due east, 180 degrees represents due south, and 270 degrees represents due west. Define a 5x5 pixel window and calculate the difference between the maximum elevation value and the minimum elevation value within the window as the terrain undulation value of the central pixel of the window, in meters, to obtain the slope map, aspect map, and terrain undulation map;

[0021] Perform solar geometric parameter calculation based on the UAV GPS data, slope map, and aspect map to obtain the solar altitude angle map, solar azimuth angle map, and incidence angle map;

[0022] Perform direct radiation correction coefficient calculation based on the solar altitude angle map, solar azimuth angle map, and incidence angle map to obtain the total direct radiation correction coefficient map;

[0023] Calculate the scattering radiation correction coefficient according to the slope map and aspect map to obtain the scattering radiation correction coefficient map;

[0024] Calculate the adjacency effect correction coefficient for the multispectral image to obtain the adjacency effect correction coefficient map;

[0025] Construct a band - by - band radiative transfer model based on the total direct radiation correction coefficient map, the scattering radiation correction coefficient map, and the adjacency effect correction coefficient map to obtain a terrain - spectral coupling model.

[0026] In the present invention, by using the processed DEM data, three key micro - geomorphic parameters, namely slope, aspect, and terrain undulation degree, are calculated. These parameters quantitatively describe the micro - topographic characteristics of the farmland area, providing basic data for subsequent analysis of the impact of terrain on radiation and vegetation indices, and are the key inputs for terrain correction. Combining the UAV GPS data and the micro - geomorphic parameters, the solar altitude angle, solar azimuth angle, and incident angle of each pixel are accurately calculated. These parameters accurately describe the geometric position of the sun relative to the earth's surface, taking into account the influence of terrain, providing key geometric information for subsequent radiation correction and improving the accuracy of radiation correction. Based on the solar geometric parameters, the total direct radiation correction coefficient is calculated, which quantitatively describes the impact of terrain on direct solar radiation. Through this coefficient, the direct radiation difference caused by terrain can be effectively eliminated, providing a guarantee for the accurate calculation of subsequent vegetation indices. Based on the slope and aspect, the scattering radiation correction coefficient is calculated, which quantitatively describes the impact of terrain on sky - scattered radiation. Through this coefficient, the scattered radiation difference caused by terrain can be effectively eliminated, further improving the accuracy of radiation correction. Based on the corrected multispectral image, the adjacency effect correction coefficient is calculated, which quantitatively describes the radiation interaction between adjacent pixels. Through this coefficient, the radiation interference between pixels in the terrain - complex area can be effectively weakened, improving the accuracy of reflectance data. Integrating the direct radiation correction coefficient, the scattering radiation correction coefficient, and the adjacency effect correction coefficient, a band - by - band terrain - spectral coupling model is constructed. This model comprehensively considers various impacts of terrain on solar radiation, can correct the observed terrain - affected reflectance to the horizontal - plane reflectance, providing a key guarantee for the accurate calculation of subsequent vegetation indices and moisture indices, and is the core for improving the inversion accuracy.

[0027] Preferably, the multi - dimensional state analysis of vegetation in step S2 is specifically as follows:

[0028] Calculate the initial spectral index by using the multispectral image; correct the initial spectral index by using the terrain - spectral coupling model to obtain the vegetation index map and the moisture index map;

[0029] Calculate the canopy temperature by using the thermal infrared image and the vegetation index map to obtain the canopy temperature map;

[0030] Based on the vegetation index map, water index map and canopy temperature map, a preliminary judgment of water stress is made to obtain a preliminary water stress judgment map.

[0031] In the present invention, by using the corrected multi-spectral images, initial spectral indices such as NDVI and NDWI are calculated. These indices can reflect the growth status and water content of vegetation, providing preliminary quantitative indicators for subsequent vegetation status analysis. The initial spectral indices are corrected using a terrain-spectral coupling model to eliminate the influence of terrain on the spectral indices, obtaining more accurate vegetation index maps and water index maps. This improves the reliability of the vegetation index and water index in areas with complex terrain, laying a foundation for subsequent analysis. Using the corrected thermal infrared image and the corrected vegetation index map, a canopy temperature map is calculated. Canopy temperature is an important indicator of crop water status, and this step provides key thermal infrared information for water stress analysis. By comprehensively using the corrected vegetation index map, water index map and canopy temperature map, and adopting a decision tree method, a preliminary water stress judgment map is generated. This map preliminarily divides the water stress levels of farmland, providing a basis for more refined water shortage diagnosis in the future and realizing a preliminary assessment of crop water status.

[0032] Preferably, the estimation of soil heat flux in step S3 is specifically as follows:

[0033] Calculate the vegetation coverage according to the vegetation index map to obtain a vegetation coverage map;

[0034] Calculate the terrain correction coefficient according to the terrain undulation map to obtain a terrain correction coefficient map;

[0035] Perform a vegetation coverage correction calculation on the vegetation coverage map and the terrain correction coefficient map to obtain a corrected vegetation coverage map;

[0036] Use the corrected vegetation coverage map and the net radiation map to estimate the ratio of soil heat flux to net radiation to obtain a G / Rn ratio map, where G is the soil heat flux and Rn is the net radiation value;

[0037] Calculate the soil heat flux according to the G / Rn ratio map to obtain a soil heat flux map.

[0038] In the present invention, by using the corrected vegetation index map (NDVI), a vegetation coverage map is obtained through the pixel dichotomy model. Vegetation coverage is a key parameter affecting the surface energy balance, and this step provides important vegetation information for the subsequent estimation of soil heat flux. The terrain correction coefficient is calculated using the terrain undulation map. The terrain undulation reflects the roughness of the surface, and this coefficient takes into account the influence of terrain on the calculation of vegetation coverage, improving the accuracy of vegetation coverage estimation, especially in areas with complex terrain. By combining the vegetation coverage map with the terrain correction coefficient map, a corrected vegetation coverage map is obtained. This map more accurately reflects the actual vegetation coverage situation and further improves the accuracy of subsequent soil heat flux estimation. Using the corrected vegetation coverage map and the net radiation map, the ratio of soil heat flux to net radiation (G / Rn) is estimated. G / Rn is a key parameter for estimating soil heat flux, and this step establishes the relationship between G / Rn and vegetation coverage through an empirical formula, providing a basis for the calculation of soil heat flux. Using the G / Rn ratio map and the net radiation map, a soil heat flux map is calculated. Soil heat flux is an important part of the surface energy balance, and the calculation result of this step provides a key input for the subsequent inversion of latent heat flux and evapotranspiration, which is a key link in the energy balance closure.

[0039] Preferably, the calculation of the apparent thermal inertia and the estimation of the sensible heat flux in step S3 are specifically as follows:

[0040] Extract temperature data of different time phases from the thermal infrared image to obtain a temperature data set;

[0041] Calculate the apparent thermal inertia according to the meteorological interpolation data and the temperature data set to obtain an apparent thermal inertia map;

[0042] Estimate the sensible heat flux according to the canopy temperature map, the meteorological interpolation data and the vegetation index map to obtain a sensible heat flux map.

[0043] In the present invention, by extracting the brightness temperature data of at least two different time phases from the corrected thermal infrared image, the necessary time series temperature information is provided for calculating the apparent thermal inertia. Selecting time phases with a large temperature difference helps to improve the accuracy of apparent thermal inertia calculation. Using the temperature data of different time phases, meteorological data (solar radiation) and surface albedo, an apparent thermal inertia map is calculated. Apparent thermal inertia reflects the soil thermal properties and moisture conditions, providing important soil information for the subsequent inversion of evapotranspiration and calculation of irrigation amount. Using the canopy temperature map, meteorological interpolation data (wind speed, air temperature) and the corrected vegetation index map, the sensible heat flux map is estimated by the aerodynamic method. Sensible heat flux is another important part of the surface energy balance, and the calculation result of this step, together with the net radiation and soil heat flux, provides a key input for the inversion of latent heat flux and evapotranspiration, ensuring the closure of the energy balance.

[0044] Preferably, the evapotranspiration inversion described in step S3 is specifically as follows:

[0045] Calculate the latent heat flux based on the net radiation map, the apparent thermal inertia map, and the sensible heat flux map to obtain the latent heat flux map;

[0046] Calculate the latent heat of vaporization based on the air temperature in the meteorological interpolation data to obtain the latent heat of vaporization map;

[0047] Calculate the initial evapotranspiration based on the latent heat flux map and the latent heat of vaporization map to obtain the initial evapotranspiration map;

[0048] Use the apparent thermal inertia map to assist in selecting cold pixels and the initial evapotranspiration map to assist in selecting hot pixels to obtain a cold pixel set and a hot pixel set;

[0049] Calculate the average latent heat flux of cold pixels and the average sensible heat flux of hot pixels, and perform evapotranspiration ratio factor calculation to obtain the evapotranspiration ratio factor map;

[0050] Calculate the average evapotranspiration of cold pixels, and perform final evapotranspiration calculation based on the evapotranspiration ratio factor map to obtain the evapotranspiration distribution map.

[0051] The present invention calculates the latent heat flux map based on the net radiation map, the soil heat flux map, and the sensible heat flux map according to the energy balance principle. The latent heat flux is the energy consumed in the evapotranspiration process, and the calculation result of this step is directly related to the inversion accuracy of evapotranspiration and is the final link of energy balance closure. Using the air temperature in the meteorological interpolation data, the latent heat of vaporization map is calculated. The latent heat of vaporization is the energy required for water to change from liquid to gas, and this parameter links the latent heat flux with the evapotranspiration amount and is the key to the conversion between energy units and water quantity units. Using the latent heat flux map and the latent heat of vaporization map, the initial evapotranspiration map is calculated. This map provides a preliminary estimate of evapotranspiration and provides a basis for subsequent refined evapotranspiration calculation. Use the apparent thermal inertia map to assist in selecting cold pixels (high ATI, high NDVI), and use the initial evapotranspiration map to assist in selecting hot pixels (low ETi, low NDVI). This method is more objective and reliable than traditional methods, improves the accuracy of cold and hot pixel selection, and provides key end-member pixels for the calculation of the evapotranspiration ratio factor. Based on the average latent heat flux and sensible heat flux of cold and hot pixels, the evapotranspiration ratio factor map is calculated. This factor reflects the proportion of evapotranspiration at the current moment to potential evapotranspiration and is a key parameter for calculating actual evapotranspiration. Using the evapotranspiration ratio factor map and the average evapotranspiration of cold pixels (as potential evapotranspiration), the final evapotranspiration distribution map is calculated. This map provides high-spatial-resolution information on actual evapotranspiration of farmland and is a key basis for water shortage diagnosis and formulating irrigation prescription maps.

[0052] Preferably, step S4 includes the following steps:

[0053] Step S41: Obtain historical irrigation data; perform preprocessing on the historical irrigation data, the initial moisture stress judgment map, and the evapotranspiration distribution map to obtain a sample data set;

[0054] Step S42: Extract multi-dimensional features from the sample data set to obtain a phenological feature data set;

[0055] Step S43: Select and train a model based on the feature data set to obtain a water shortage diagnosis model;

[0056] Step S44: Input the feature data set into the water shortage diagnosis model to generate a water shortage diagnosis index, and obtain the water shortage diagnosis index.

[0057] In the present invention, by performing spatio-temporal matching and integration on historical irrigation data, soil moisture data, crop yield data, the initial moisture stress judgment map, and the evapotranspiration distribution map, and performing data cleaning and standardization processing, a high-quality sample data set is constructed. This data set provides reliable and representative samples for the training of the water shortage diagnosis model and is the basis for model construction. Multi-dimensional features such as the cumulative evapotranspiration, water deficit index, time features, historical irrigation features, and statistical values of meteorological features are extracted from the sample data, and topographic features are combined. These features comprehensively consider various factors affecting crop water shortage, increase the input information volume of the model, and provide more comprehensive information for the accurate prediction of the model. The random forest regression model is selected as the water shortage diagnosis model, and the model is trained, parameter-tuned, and evaluated using the training set, validation set, and test set. The random forest model has strong non-linear fitting ability and anti-noise ability, can effectively process multi-dimensional feature input, and improves the prediction accuracy and generalization ability of the model. The feature data set is input into the trained water shortage diagnosis model to generate a water shortage diagnosis index (WDI). WDI is a continuous numerical value that can quantitatively represent the severity of crop water shortage, provides a more refined and accurate basis for irrigation decision-making, and realizes the transformation from qualitative evaluation to quantitative diagnosis.

[0058] Preferably, step S5 includes the following steps:

[0059] Step S51: Set and classify the WDI threshold according to the water shortage diagnosis index and meteorological interpolation data to obtain a WDI classification map;

[0060] Step S52: Divide the irrigation management area according to the WDI classification map to obtain an irrigation management area division map;

[0061] Step S53: Calculate the irrigation amount according to the phenological feature data and the irrigation management area division map to obtain an irrigation amount allocation table;

[0062] Step S54: Assign the irrigation amount information in the irrigation amount distribution table to grid pixels according to the spatial positions in the irrigation management zoning map to obtain an irrigation prescription map.

[0063] In the present invention, by setting a WDI threshold according to the water shortage diagnosis index (WDI), meteorological data (especially future rainfall forecasts), in combination with crop water requirement characteristics and historical irrigation experience, farmland is divided into different water shortage levels. This method takes into account future weather conditions, realizes dynamic adjustment of irrigation decisions, and improves the pertinence and timeliness of irrigation. According to the WDI classification map, using the region growing algorithm, farmland is divided into different irrigation management areas. This zoning management method takes into account the spatial differences in crop water shortage conditions and combines the actual layout of the irrigation system, facilitating the implementation of precise irrigation and improving irrigation efficiency. According to the phenological characteristic data and the irrigation management zoning map, a method based on crop water requirements is used to calculate the irrigation amount for each management area. This method comprehensively considers crop type, growth stage, meteorological conditions, and soil water holding capacity, realizes refined calculation of the irrigation amount, avoids over-irrigation or under-irrigation, and improves water resource utilization efficiency. The information in the irrigation amount distribution table is assigned to grid pixels according to the spatial positions in the irrigation management zoning map to generate an irrigation prescription map. This map intuitively shows the irrigation requirements of different areas of farmland and can be directly used to guide the operation of the precise irrigation system, realizing the transformation from data analysis to practical application, which is the core result of the whole method.

[0064] Preferably, step S53 is specifically as follows:

[0065] Calculate potential evapotranspiration according to meteorological interpolation data to obtain a potential evapotranspiration map;

[0066] Calculate the crop coefficient according to the phenological characteristic data and the vegetation index map to obtain a crop coefficient map;

[0067] Estimate actual evapotranspiration according to the evapotranspiration distribution map, the crop coefficient map, and the potential evapotranspiration map to obtain an actual evapotranspiration map;

[0068] Calculate the soil water holding capacity correction coefficient according to the apparent thermal inertia map to obtain a soil water holding capacity correction coefficient map;

[0069] Preliminarily calculate the irrigation amount according to the actual evapotranspiration map to obtain a preliminary irrigation amount map;

[0070] Use the soil water holding capacity correction coefficient map to correct the preliminary irrigation amount map to obtain a corrected irrigation amount map;

[0071] Allocate the irrigation amount according to the corrected irrigation amount map and the irrigation management zoning map to obtain an irrigation amount distribution table.

[0072] In the present invention, by utilizing meteorological interpolation data and calculating with the Penman-Monteith formula, a potential evapotranspiration map is obtained. Potential evapotranspiration is the benchmark for calculating crop water requirements, and this step provides an important meteorological basis for the calculation of irrigation amounts. According to the phenological characteristic data and the corrected vegetation index map, a crop coefficient map is calculated using the piecewise linear function method. The crop coefficient reflects the water requirement differences of different crops at different growth stages. This step combines the crop water requirement characteristics with remote sensing data, improving the pertinence of the irrigation amount calculation. Using the evapotranspiration distribution map obtained in step S3 as the measured evapotranspiration, this map provides high-spatial-resolution information on the actual evapotranspiration of farmland, avoiding errors caused by model estimation. A soil water holding capacity correction coefficient map is calculated using the apparent thermal inertia map. This coefficient reflects the water retention capacity of the soil, taking into account the influence of soil factors on irrigation demand, making the irrigation amount calculation more in line with the actual situation. A preliminary irrigation amount map is calculated using the actual evapotranspiration map, which reflects the actual water requirements at the farmland scale. The preliminary irrigation amount map is corrected using the soil water holding capacity correction coefficient map to obtain the corrected irrigation amount map. This correction takes into account the influence of soil factors on irrigation demand, making the irrigation amount calculation more in line with the actual situation and avoiding over-irrigation of areas with strong water retention capacity. According to the corrected irrigation amount map and the irrigation management zoning map, an irrigation amount allocation table is calculated. This table corresponds the irrigation amount to the irrigation management area, providing the final quantitative data for the generation of the irrigation prescription map and realizing the conversion from the continuous irrigation amount map to the zoned management irrigation amount. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 FIG. is a schematic flow chart of the steps of a method for inverting an irrigation prescription map based on UAV spectral data;

[0074] Figure 2 FIG. is a detailed schematic flow chart of step S5 in the present invention.

[0075] The implementation, functional features, and advantages of the present invention will be further described with reference to the embodiments and the accompanying drawings. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0076] The technical method of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art within the scope of the present invention without creative efforts belong to the scope of protection of the present invention.

[0077] In addition, the accompanying drawings are only schematic illustrations of the present invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and thus repeated descriptions thereof will be omitted. Some of the block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. The functional entities may be implemented in software form, or in one or more hardware modules or integrated circuits, or in different networks and / or processor methods and / or microcontroller methods.

[0078] It should be understood that although the terms "first", "second", etc. may be used herein to describe various units, these units should not be limited by these terms. These terms are only used to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, the first unit may be referred to as the second unit, and similarly the second unit may be referred to as the first unit. The term "and / or" used herein includes any and all combinations of one or more of the listed associated items.

[0079] In the embodiments of the present invention, reference is made to Figure 1 As shown, it is a schematic diagram of the step flow of the method for inverting the irrigation prescription map based on the drone spectral data of the present invention. In this example, the method for inverting the irrigation prescription map based on the drone spectral data includes the following steps:

[0080] Step S1: Collect drone data for the farmland area, and perform multi-source data registration and fusion to obtain a fusion data set, where the fusion data set includes multi-spectral images, thermal infrared images, DEM data, and meteorological interpolation data;

[0081] In the embodiments of the present invention, a drone is used to carry a multispectral camera and a thermal infrared camera to conduct aerial photography of farmland areas, obtain original multispectral images and original thermal infrared images, and record the drone's GPS data at the same time. Ground control points are arranged in the farmland and their precise coordinates are measured. Meteorological data (temperature, humidity, wind speed, radiation, rainfall) synchronized with the drone flight are obtained from a nearby meteorological station. The original multispectral images are subjected to radiometric calibration, geometric correction, and atmospheric correction (optional) to obtain multispectral images. The original thermal infrared images are subjected to radiometric calibration and geometric correction, and atmospheric correction is performed using the ground-measured temperature to obtain thermal infrared images. The obtained DEM data are denoised, interpolated, and smoothed, and converted to the same coordinate system and resolution as the images to obtain DEM data. Spatial interpolation (such as inverse distance weighting interpolation or Kriging interpolation) is performed on the meteorological station data to obtain meteorological interpolation data with the same resolution as the images. Finally, based on the calibrated multispectral images, the thermal infrared images, DEM data, and meteorological interpolation data are registered to the same coordinate system using an image registration algorithm (such as registration based on mutual information), and all the data are stacked by band to generate a fusion dataset.

[0082] Step S2: Calculate microtopographic parameters based on the multispectral images, DEM data, and drone GPS data, and construct a terrain-spectral coupling model to obtain a terrain undulation map and a terrain-spectral coupling model; use the terrain-spectral coupling model to conduct a multi-dimensional state analysis of vegetation to obtain a vegetation index map, a canopy temperature map, and a preliminary water stress judgment map;

[0083] In the embodiments of the present invention, using the DEM data, the slope, aspect, and terrain undulation of each pixel are calculated. According to the drone GPS data (flight time, longitude and latitude) and the slope and aspect data, the solar altitude angle, solar azimuth angle, and incident angle of each pixel are calculated. A radiation transfer model based on a physical model is constructed to calculate the direct radiation correction coefficient, scattered radiation correction coefficient, and adjacent effect correction coefficient respectively, and the model parameters are adjusted according to different bands. Using the multispectral images, spectral indices such as NDVI and NDWI are calculated, and these indices are corrected using the constructed terrain-spectral coupling model to obtain a vegetation index map and a water index map. Using the thermal infrared images, according to Planck's law and the Stefan-Boltzmann law, combined with the surface emissivity estimated by NDVI, the canopy temperature is calculated to obtain a canopy temperature map. By comprehensively analyzing the vegetation index map, water index map, and canopy temperature map, a decision tree method is used to generate a preliminary water stress judgment map.

[0084] Step S3: Calculate the net radiation using the meteorological interpolation data to obtain a net radiation map; estimate the soil heat flux based on the net radiation map, vegetation index map, and topographic undulation map to obtain a soil heat flux map; calculate the apparent thermal inertia based on the thermal infrared image and canopy temperature map, and estimate the sensible heat flux to obtain an apparent thermal inertia map and a sensible heat flux map; perform evapotranspiration inversion based on the soil heat flux map, apparent thermal inertia map, and sensible heat flux map to obtain an evapotranspiration distribution map;

[0085] In the embodiment of the present invention, the net radiation is calculated using the albedo calculated from the meteorological interpolation data (solar radiation, air temperature, etc.) and the multispectral image. Based on the vegetation index map and the topographic undulation map, the ratio of the soil heat flux to the net radiation is estimated, and then the soil heat flux is calculated. The apparent thermal inertia is calculated using at least two-phase thermal infrared images and meteorological data, combined with the surface albedo and solar altitude angle. The sensible heat flux is estimated using the aerodynamic method based on the canopy temperature map, meteorological interpolation data (wind speed, air temperature), and vegetation index map. According to the energy balance equation (Rn = G + H + LE), the latent heat flux is calculated using the net radiation, soil heat flux, and sensible heat flux. The latent heat of vaporization of water is calculated using the air temperature in the meteorological interpolation data. The latent heat flux is divided by the latent heat of vaporization to obtain the initial evapotranspiration. The apparent thermal inertia map is used to assist in selecting cold pixels (high ATI, high NDVI), and the initial evapotranspiration map is used to assist in selecting hot pixels (low ETi, low NDVI). The average latent heat flux of the cold pixels and the average sensible heat flux of the hot pixels are calculated, and then the evapotranspiration scaling factor is calculated. The average evapotranspiration of the cold pixels is calculated as the potential evapotranspiration, and the evapotranspiration scaling factor is multiplied by the potential evapotranspiration to obtain the final evapotranspiration distribution map.

[0086] Step S4: Extract sample characteristics based on the initial water stress judgment map and the evapotranspiration distribution map, and perform water shortage judgment to obtain a phenological feature dataset and a water shortage diagnosis index;

[0087] In the embodiment of the present invention, historical irrigation data, soil moisture data, and crop yield data are collected, and spatiotemporally matched with the initial water stress judgment map and the evapotranspiration distribution map to construct a sample dataset. Features are extracted from the sample data, including cumulative evapotranspiration, water deficit index, time features, historical irrigation features, and statistical values of meteorological features. These features are combined with the original features (initial water stress judgment level, evapotranspiration) to form a feature dataset. A random forest regression model is selected as the water shortage diagnosis model, and the feature dataset is divided into a training set, a validation set, and a test set. The model is trained using the training set, the model parameters are adjusted using the validation set, and the model performance is evaluated using the test set. The feature dataset is input into the trained model to predict the water shortage diagnosis index (WDI) for each pixel, generating a water shortage diagnosis index map.

[0088] Step S5: Divide the irrigation management areas based on the phenological characteristic data and the water shortage diagnosis index, calculate the irrigation amount, and generate an irrigation prescription map;

[0089] In the embodiment of the present invention, according to the water shortage diagnosis index (WDI) and the meteorological interpolation data (future rainfall forecast), combined with the water demand characteristics of local crops and historical irrigation experience, a WDI threshold is set, and the WDI map is divided into four levels: severe water shortage, moderate water shortage, mild water shortage, and sufficient water, to obtain a WDI classification map. The region growing algorithm is used to merge the pixels with the same WDI level and spatially adjacent into the same irrigation management area, and considering the actual layout of the irrigation system, an irrigation management area division map is obtained. Calculate the potential evapotranspiration according to the meteorological interpolation data, calculate the crop coefficient according to the phenological characteristic data and the vegetation index map, and estimate the actual evapotranspiration in combination with the actual evapotranspiration map. Calculate the soil water holding capacity correction coefficient according to the apparent thermal inertia map. Initially calculate the irrigation amount according to the actual evapotranspiration amount, and correct it using the soil water holding capacity correction coefficient. Calculate the average corrected irrigation amount for each irrigation management area, store the management area ID and the corresponding irrigation amount in the irrigation amount distribution table. According to the irrigation management area division map, assign the irrigation amount value in the irrigation amount distribution table to the corresponding pixels to generate an irrigation prescription map.

[0090] Preferably, step S1 includes the following steps:

[0091] Step S11: Uniformly arrange at least 10 ground control points in the farmland, use an RTK GNSS receiver to measure the coordinates of the control points, with both horizontal and vertical accuracies better than 2 cm. Plan the UAV flight path to fly in an east-west strip shape, set the flight altitude to 100 meters, the heading overlap to 80%, and the side overlap to 70%. Collect UAV data for the farmland area through the UAV, and synchronously obtain meteorological data to obtain the original spectral data, the original digital elevation model data, the UAV GPS data, and the original weather station data;

[0092] Step S12: Perform spectral data preprocessing on the UAV spectral data to obtain a multispectral image and a thermal infrared image;

[0093] Step S13: Perform DEM data preprocessing on the digital elevation model data to obtain DEM data;

[0094] Step S14: Perform meteorological data interpolation processing on the original weather station data to obtain meteorological interpolation data.

[0095] In the embodiments of the present invention, in the target farmland area, a DJI Matrice 300 RTK drone is used, equipped with a Zenmuse P1 multispectral imaging system and a Zenmuse H20T thermal infrared camera. The P1 multispectral imaging system includes visible light bands (blue, green, red) and a near-infrared band, with central wavelengths of 450 nm, 560 nm, 650 nm, and 840 nm respectively, and the spectral resolution is better than 10 nm. The wavelength range of the H20T thermal infrared camera is 8 - 14 μm, and the thermal sensitivity (NETD) ≤ 50 mK @ f / 1.0. Before flight, at least 10 ground control points are evenly arranged in the farmland, and the coordinates of the control points are measured using an RTK GNSS receiver, with the horizontal and vertical accuracies both better than 2 cm. The drone flight path is planned to fly in an east-west strip shape, the flight altitude is set at 100 meters, the forward overlap is 80%, and the side overlap is 70% to ensure the quality of image stitching. The flight time is selected between 10:00 and 12:00 in the morning on a sunny and cloudless day with a wind speed less than 5 m / s and a solar altitude angle greater than 45°. The drone GPS data, including longitude, latitude, altitude, and time information, is synchronously recorded, and the sampling frequency is 1 Hz. At the same time, meteorological data synchronized with the drone flight period is obtained from the meteorological station closest to the farmland (not exceeding 5 kilometers), including temperature (accuracy ±0.5 °C), humidity (accuracy ±3%), wind speed (accuracy ±0.2 m / s), total radiation (accuracy ±5 W / m 2 ²) and rainfall (accuracy ±0.1 mm), and the data time resolution is 1 hour. Finally, the original multispectral images, original thermal infrared images, drone GPS data, and original meteorological station data are obtained.

[0096] Preprocess the original multispectral images and original thermal infrared images obtained in step S11. First, using the standard gray scale board reflectance data synchronously recorded during image shooting, perform radiometric calibration on each band of the original multispectral images to convert the DN values into surface reflectance. The radiometric calibration uses an empirical linear model: L = Gain * DN + Offset, where Gain and Offset are calculated from the gray scale board reflectance and the corresponding DN values. Then, using the drone GPS data and the coordinates of the ground control points, use a polynomial geometric correction model to perform geometric fine correction on the multispectral images and thermal infrared images. The polynomial order is selected as quadratic, and the resampling method uses bilinear interpolation. The geometric error of the corrected images is less than 1 pixel. For the thermal infrared images, use the Planck formula to convert the DN values into brightness temperature, and use the synchronously measured ground target (such as asphalt pavement, water body) temperature data for atmospheric correction, using a single-channel atmospheric correction algorithm to eliminate the influence of atmospheric absorption and scattering. Finally, register the corrected multispectral images and thermal infrared images. Taking the multispectral images as the reference, use a feature point-based registration method to extract SIFT features, and use the RANSAC algorithm to eliminate mis-matched points. The registration accuracy is controlled at the sub-pixel level. Output the multispectral images and thermal infrared images.

[0097] Preprocess the original digital elevation model (DEM) data obtained in step S11. First, check whether there are holes and outliers in the DEM data. For the hole areas, use the Kriging interpolation method for filling. Kriging interpolation uses a semivariogram model, and the model parameters are obtained by fitting the existing data. For the abnormal elevation values (elevation values that deviate significantly from the surrounding areas), use the median filtering method for smoothing, and set the filter window size to 3x3 pixels. Then, perform a projection transformation on the DEM data to convert its coordinate system to the same UTM projection coordinate system as the corrected multispectral image. Finally, resample the DEM data to make its spatial resolution consistent with that of the corrected multispectral image (for example, both are 1 meter), and use the bilinear interpolation method for resampling. Output the DEM data.

[0098] Perform spatial interpolation processing on the original weather station data obtained in step S11. Since the weather station data are discrete point data, and the irrigation prescription map requires continuous spatial distribution data, interpolation is needed. Taking the temperature data as an example, use the inverse distance weighted interpolation method (IDW) to interpolate the temperature data observed at the weather stations into raster data with the same resolution as the corrected multispectral image. The weight of the IDW method is inversely proportional to the square of the distance, and the search radius is set to 5 kilometers to ensure that each pixel is affected by at least 3 weather station data. The same interpolation method is also used for the humidity, wind speed, and total radiation data, but the search radius and weight index are adjusted according to the spatial variability characteristics of each element. Since the rainfall data has large spatial variability, use the Thiessen polygon method for interpolation, and the rainfall within each polygon is equal to the observed value of its corresponding weather station. Output the meteorological interpolation data, including raster data of temperature, humidity, wind speed, total radiation, and rainfall. Perform multi-source data registration and fusion on the multispectral image, thermal infrared image, DEM data, and meteorological interpolation data. First, ensure that all data have the same spatial reference (UTM projection coordinate system) and the same spatial resolution (for example, 1 meter). Then, taking the corrected multispectral image as the reference, use the image registration method to accurately register the thermal infrared image, DEM data, and meteorological interpolation data into the coordinate system of the multispectral image. The registration uses a registration algorithm based on mutual information, which has good robustness for images of different modalities. Finally, stack the registered multispectral image (multiple bands), thermal infrared image (single band), DEM data (single band), and meteorological interpolation data (multiple bands) in the band order to form a multi-dimensional data cube, that is, the fusion dataset. Each pixel of this dataset contains multispectral, thermal infrared, terrain, and meteorological information.

[0099] Preferably, the calculation of the microtopography parameters and the construction of the terrain-spectral coupling model described in step S2 are specifically as follows:

[0100] Calculate the microtopographic parameters of the DEM data. According to the elevation differences in the horizontal and vertical directions, calculate the angle between the slope normal vector and the due north direction, in degrees, with a range of 0 - 360 degrees, where 0 degrees represents due north, 90 degrees represents due east, 180 degrees represents due south, and 270 degrees represents due west. Define a 5x5 pixel window, and calculate the difference between the maximum elevation value and the minimum elevation value within the window as the topographic undulation value of the pixel at the window center, in meters, to obtain the slope map, aspect map, and topographic undulation map.

[0101] Calculate the solar geometric parameters based on the UAV GPS data, slope map, and aspect map to obtain the solar altitude angle map, solar azimuth angle map, and incident angle map.

[0102] Calculate the direct radiation correction coefficient based on the solar altitude angle map, solar azimuth angle map, and incident angle map to obtain the total direct radiation correction coefficient map.

[0103] Calculate the diffuse radiation correction coefficient based on the slope map and aspect map to obtain the diffuse radiation correction coefficient map.

[0104] Calculate the adjacency effect correction coefficient for the multispectral image to obtain the adjacency effect correction coefficient map.

[0105] Construct a band - by - band radiative transfer model based on the total direct radiation correction coefficient map, diffuse radiation correction coefficient map, and adjacency effect correction coefficient map to obtain the terrain - spectral coupling model.

[0106] In the embodiment of the present invention, using the DEM data obtained in step S13, calculate the slope, aspect, and topographic undulation. The slope calculation uses the Horn algorithm, which is based on the elevation differences between each pixel in the DEM data and its eight adjacent pixels around it, calculates the maximum rate of change of the pixel in the horizontal and vertical directions, and then takes the square root of the sum of their squares to obtain the slope value, in degrees. The aspect calculation also uses the Horn algorithm. According to the elevation differences in the horizontal and vertical directions, calculate the angle between the slope normal vector and the due north direction, in degrees, with a range of 0 - 360 degrees, where 0 degrees represents due north, 90 degrees represents due east, 180 degrees represents due south, and 270 degrees represents due west. The topographic undulation calculation uses the neighborhood analysis method. Define a 5x5 pixel window, and calculate the difference between the maximum elevation value and the minimum elevation value within the window as the topographic undulation value of the pixel at the window center, in meters. Traverse the entire DEM data to obtain the slope map, aspect map, and topographic undulation map.

[0107] Based on the UAV GPS data obtained in step S11 (recording the image shooting time, longitude and latitude information) and the slope map and aspect map obtained in the previous step, calculate the solar altitude angle, solar azimuth angle and incident angle for each pixel. First, based on the image shooting time (day of the year and hours, minutes and seconds) and geographical location (longitude and latitude), calculate the solar zenith angle and solar azimuth angle at that moment. The calculation of the solar zenith angle and azimuth angle uses the standard solar position calculation formula, considering the influence of the earth's rotation and revolution. Then, based on the slope (β) and aspect (α) of each pixel, as well as the solar zenith angle (θ) and solar azimuth angle (φ), calculate the incident angle (i). The incident angle calculation formula is: cos(i) = cos(β)cos(θ) + sin(β)sin(θ)cos(φ - α). Traverse the entire study area to obtain the solar altitude angle map (90° - solar zenith angle), solar azimuth angle map and incident angle map.

[0108] Based on the solar altitude angle map, solar azimuth angle map and incident angle map calculated in the previous step, calculate the total direct radiation correction coefficient for each pixel. The direct radiation correction coefficient takes into account the influence of the solar altitude angle and incident angle on the direct solar radiation received by the ground surface. Use the C-factor method to calculate the direct radiation correction coefficient (Cd): Cd = cos(i) / cos(θ), where i is the incident angle and θ is the solar zenith angle. For the shadow area (cos(i) < 0), Cd is set to 0. Traverse the entire study area to obtain the total direct radiation correction coefficient map.

[0109] Based on the slope map and aspect map obtained from the microtopography parameter calculation step, calculate the diffuse radiation correction coefficient for each pixel. The diffuse radiation correction coefficient takes into account the influence of the terrain on the sky diffuse radiation. Use the method proposed by Dozier and Frew to calculate the sky view factor (Vd): Vd = (1 + cos(β)) / 2, where β is the slope. Assuming isotropic atmosphere, the diffuse radiation correction coefficient (Cs) is equal to the sky view factor: Cs = Vd. Traverse the entire study area to obtain the diffuse radiation correction coefficient map.

[0110] For the multispectral image obtained in step S12, calculate the adjacent effect correction coefficient for each pixel. The adjacent effect refers to the fact that the radiation received by the target pixel not only comes from itself, but is also affected by the reflected radiation of the surrounding adjacent pixels. Use a simplified adjacent effect model, assuming that the adjacent effect is proportional to the weighted average of the reflectivities of the surrounding pixels. Define a 5x5 pixel window, calculate the difference between the reflectivity of each pixel in the window and the reflectivity of the target pixel, and then use the Gaussian function to perform weighted averaging on these differences, with the weights inversely proportional to the square of the distance. The adjacent effect correction coefficient (Cp) is equal to 1 plus this weighted average. Calculate the adjacent effect correction coefficient for each band of the multispectral image separately. Traverse the entire study area to obtain the adjacent effect correction coefficient map for each band.

[0111] Based on the total direct radiation correction coefficient map (Cd), diffuse radiation correction coefficient map (Cs), and proximity effect correction coefficient map (Cp) calculated in the previous steps, a sub-band topographic-spectral coupling model is constructed. This model is used to correct the observed terrain-affected reflectance (ρt) to the horizontal plane reflectance (ρ0). The model uses the following formula: ρ0 = (ρt - Cs * ρd) / (Cd * τ + Cp), where ρd is the ratio of downward diffuse radiation to total radiation in the atmosphere, and τ is the atmospheric transmittance. ρd and τ can be obtained by simulating with an atmospheric radiative transfer model (such as MODTRAN) or calculated using synchronous atmospheric parameter measurement data. For each band (blue, green, red, near-infrared) of the multispectral image, the above correction model is established separately to obtain the sub-band topographic-spectral coupling model. The model parameters (ρd and τ) are adjusted according to the central wavelength of different bands and atmospheric conditions.

[0112] Preferably, the multi-dimensional state analysis of vegetation described in step S2 is specifically as follows:

[0113] Calculate the spectral index using the multispectral image to obtain the initial spectral index; use the topographic-spectral coupling model to correct the initial spectral index to obtain the vegetation index map and the moisture index map;

[0114] Calculate the canopy temperature using the thermal infrared image and the vegetation index map to obtain the canopy temperature map;

[0115] Based on the vegetation index map, moisture index map, and canopy temperature map, make a preliminary judgment on water stress to obtain the preliminary water stress judgment map.

[0116] In the embodiments of the present invention, the normalized difference vegetation index (NDVI) and the normalized difference water index (NDWI) are calculated using the multispectral images obtained in step S12. The formula for calculating NDVI is: NDVI = (ρNIR - ρRed) / (ρNIR + ρRed), where ρNIR is the reflectance in the near-infrared band and ρRed is the reflectance in the red band. The formula for calculating NDWI is: NDWI = (ρGreen - ρNIR) / (ρGreen + ρNIR), where ρGreen is the reflectance in the green band. The entire multispectral image is traversed to obtain the initial NDVI map and the initial NDWI map. Then, the initial NDVI and the initial NDWI are corrected using the sub-band terrain-spectral coupling model constructed in the previous step. Since NDVI and NDWI are calculated based on the ratio of band reflectances, the correction method is to apply the model to each reflectance term in the numerator and denominator. Taking NDVI as an example, the corrected NDVI = (ρNIR_corrected - ρRed_corrected) / (ρNIR_corrected + ρRed_corrected), where ρNIR_corrected and ρRed_corrected are the reflectances in the near-infrared and red bands corrected using the terrain-spectral coupling model, respectively. Similarly, NDWI is corrected. The entire image is traversed to obtain the vegetation index map (NDVI) and the water index map (NDWI).

[0117] Using the thermal infrared image obtained in step S12 and the vegetation index map (NDVI) obtained in the previous step, the canopy temperature is calculated. First, according to Planck's law, the brightness temperature of the thermal infrared image is converted into radiance. Then, the land surface emissivity (ε) of each pixel is estimated using NDVI. The empirical formula proposed by Vande Griend and Owe is used: ε = 1.0094 + 0.047 * ln(NDVI). For pixels with NDVI less than 0.2 (water bodies or bare soil), ε is set to 0.98; for pixels with NDVI greater than 0.7 (fully vegetated), ε is set to 0.99. Finally, according to the Stefan-Boltzmann law, the radiance is converted into the canopy temperature (Tc): Tc = (L / (ε * σ))^(1 / 4) - 273.15, where L is the radiance and σ is the Stefan-Boltzmann constant (5.67x10^-8 W / m 2 / K 4 ^4), and the calculation result is in degrees Celsius. The entire image is traversed to obtain the canopy temperature map.

[0118] Based on the vegetation index map (NDVI), water index map (NDWI), and canopy temperature map (Tc) obtained in the previous two steps, a preliminary judgment of water stress is made. The decision tree method is adopted, considering these three indicators comprehensively. First, set the threshold of NDVI to 0.4, the threshold of NDWI to 0.1, and the threshold of Tc to 35 °C. For each pixel, if NDVI < 0.4 and NDWI < 0.1 and Tc > 35 °C, it is determined as severe water stress; if NDVI < 0.6 and NDWI < 0.2 and Tc > 30 °C, it is determined as moderate water stress; if NDVI < 0.8 and NDWI < 0.3 and Tc > 25 °C, it is determined as mild water stress; otherwise, it is determined as sufficient water. Traverse the entire image to obtain the preliminary judgment map of water stress, which divides the farmland into four levels: severe stress, moderate stress, mild stress, and sufficient water.

[0119] Preferably, the estimation of soil heat flux described in step S3 is specifically as follows:

[0120] Calculate the vegetation coverage according to the vegetation index map to obtain the vegetation coverage map;

[0121] Calculate the terrain correction coefficient according to the terrain undulation map to obtain the terrain correction coefficient map;

[0122] Perform vegetation coverage correction calculation on the vegetation coverage map and the terrain correction coefficient map to obtain the corrected vegetation coverage map;

[0123] Use the corrected vegetation coverage map and the net radiation map to estimate the ratio of soil heat flux to net radiation to obtain the G / Rn ratio map, where G is the soil heat flux and Rn is the net radiation value;

[0124] Calculate the soil heat flux according to the G / Rn ratio map to obtain the soil heat flux map.

[0125] In the embodiments of the present invention, the fractional vegetation cover (FVC) is calculated using the vegetation index map (NDVI) obtained in step S2. The pixel dichotomy model method is adopted. It is assumed that each pixel consists of two parts: vegetation and soil, and NDVI is the linear weighted average of the NDVI of these two parts. The FVC calculation formula is: FVC = (NDVI - NDVIsoil) / (NDVIveg - NDVIsoil), where NDVIsoil is the NDVI value of a pure soil pixel, and NDVIveg is the NDVI value of a pure vegetation pixel. According to the NDVI histogram, the NDVI value near the minimum NDVI (e.g., 5% cumulative frequency) is taken as NDVIsoil, and the NDVI value near the maximum NDVI (e.g., 95% cumulative frequency) is taken as NDVIveg. For pixels with NDVI less than NDVIsoil, FVC is set to 0; for pixels with NDVI greater than NDVIveg, FVC is set to 1. The entire NDVI map is traversed to obtain the fractional vegetation cover map.

[0126] According to the topographic relief map obtained in step S2, the topographic correction coefficient (Ktopo) is calculated. The topographic relief reflects the roughness of the earth's surface and will affect the estimation of the fractional vegetation cover. The empirical formula is used: Ktopo = 1 / (1 + a * Relief), where Relief is the topographic relief (unit: meter), and a is an empirical coefficient with a value of 0.1. This formula indicates that the greater the topographic relief, the smaller the topographic correction coefficient, and the greater the correction amplitude for the fractional vegetation cover. The entire topographic relief map is traversed to obtain the topographic correction coefficient map.

[0127] The fractional vegetation cover map (FVC) and the topographic correction coefficient map (Ktopo) obtained in the previous two steps are used for fractional vegetation cover correction. The correction formula is: FVC_corrected = FVC * Ktopo. This formula indicates that in areas with greater topographic relief, the fractional vegetation cover will be reduced to reflect the impact of topography on vegetation growth. The entire image is traversed to obtain the corrected fractional vegetation cover map (FVC_corrected).

[0128] Using the corrected fractional vegetation cover map (FVC_corrected) obtained in the previous step and the net radiation map (Rn) calculated in step S3, the ratio of soil heat flux to net radiation (G / Rn) is estimated. The empirical formula proposed by Su (2002) is used: G / Rn = FVC_corrected * Γc + (1 - FVC_corrected) * Γs, where Γc is the G / Rn value under complete vegetation cover, with a value of 0.05; Γs is the G / Rn value under bare soil conditions, with a value of 0.315. This formula indicates that G / Rn decreases with the increase in fractional vegetation cover. The entire image is traversed to obtain the G / Rn ratio map.

[0129] Calculate the soil heat flux (G) based on the G / Rn ratio map obtained in the previous step and the net radiation map (Rn) calculated in step S3. The calculation formula is: G = (G / Rn) * Rn. Traverse the entire image to obtain the soil heat flux map, with the unit of W / m 2 .

[0130] Preferably, the calculation of the apparent thermal inertia in step S3 and the estimation of the sensible heat flux are specifically as follows:

[0131] Extract the temperature data of different time phases from the thermal infrared image to obtain a temperature data set;

[0132] Calculate the apparent thermal inertia based on the meteorological interpolation data and the temperature data set to obtain the apparent thermal inertia map;

[0133] Estimate the sensible heat flux based on the canopy temperature map, the meteorological interpolation data, and the vegetation index map to obtain the sensible heat flux map.

[0134] In the embodiment of the present invention, select at least two image data of different time phases from the thermal infrared image obtained in step S12. To ensure the calculation accuracy, select two time phases with a large temperature difference in a day, such as 10:00 am (strong solar radiation) and 14:00 pm (the highest surface temperature). Extract the thermal infrared images of these two time phases respectively, and ensure that they have been geometrically corrected and radiometrically corrected, and have the same spatial resolution. Store the brightness temperature data (unit: Kelvin) of these two time phases as two independent raster data sets to form the temperature data set.

[0135] Calculate the apparent thermal inertia (ATI) based on the meteorological interpolation data (solar radiation) obtained in step S14 and the temperature data set (brightness temperatures of two time phases) obtained in the previous step. The calculation formula of ATI is: ATI = (1 - α) * S * cos(θ) / ΔT, where α is the surface albedo, S is the solar constant (take 1367 W / m 2 ), θ is the solar altitude angle, and ΔT is the temperature difference between the two time phases. The surface albedo (α) is calculated using the multispectral image obtained in step S12, and the broadband albedo conversion formula proposed by Liang (2001) is adopted: α = 0.356ρblue + 0.130ρred + 0.373ρNIR + 0.085ρSWIR1 + 0.056 * ρSWIR2, where ρblue, ρred, ρNIR, ρSWIR1, and ρSWIR2 are the reflectances of the blue, red, near-infrared, and two shortwave infrared bands respectively (assuming that the UAV multispectral camera has these bands). The solar altitude angle (θ) uses the solar altitude angle map calculated in step S2 and takes the average value of the two time phases. ΔT is the absolute value of the difference in brightness temperatures between the two time phases. Traverse the entire image to obtain the apparent thermal inertia map.

[0136] Based on the canopy temperature map obtained in step S2, the meteorological interpolation data (wind speed, air temperature) obtained in step S14, and the vegetation index map (NDVI) obtained in step S2, the sensible heat flux (H) is estimated. Using the aerodynamic method, the calculation formula for H is: H = ρ * Cp * (Tc - Ta) / ra, where ρ is the air density, Cp is the specific heat capacity at constant pressure of air (taking 1004 J / kg / K), Tc is the canopy temperature, Ta is the air temperature, and ra is the aerodynamic resistance. The air density (ρ) is calculated based on the air pressure and temperature in the meteorological interpolation data. The air temperature (Ta) is directly obtained from the meteorological interpolation data. The aerodynamic resistance (ra) is related to the wind speed and surface roughness, and the surface roughness is estimated using NDVI. The formula proposed by Brutsaert (1982) is used: ra = ln[(z - d) / z0m] * ln[(z - d) / z0h] / (k^2 * u), where z is the reference height (usually taken as 2 meters), d is the zero-plane displacement height, z0m is the momentum roughness length, z0h is the heat roughness length, k is the von Karman constant (taking 0.41), and u is the wind speed at the reference height. d, z0m, and z0h are estimated based on NDVI using empirical formulas. Traverse the entire image to obtain the sensible heat flux map, with the unit of W / m 2 。

[0137] Preferably, the evapotranspiration inversion described in step S3 is specifically as follows:

[0138] Calculate the latent heat flux based on the net radiation map, the apparent thermal inertia map, and the sensible heat flux map to obtain the latent heat flux map;

[0139] Calculate the latent heat of vaporization based on the air temperature in the meteorological interpolation data to obtain the latent heat of vaporization map;

[0140] Conduct initial evapotranspiration calculation based on the latent heat flux map and the latent heat of vaporization map to obtain the initial evapotranspiration map;

[0141] Use the apparent thermal inertia map to assist in selecting cold pixels, and use the initial evapotranspiration map to assist in selecting hot pixels to obtain the cold pixel set and the hot pixel set;

[0142] Calculate the average latent heat flux of cold pixels and the average sensible heat flux of hot pixels, conduct evapotranspiration ratio factor calculation to obtain the evapotranspiration ratio factor map;

[0143] Calculate the average evapotranspiration amount of cold pixels, and conduct final evapotranspiration calculation based on the evapotranspiration ratio factor map to obtain the evapotranspiration distribution map.

[0144] In the embodiments of the present invention, according to the net radiation map (Rn), soil heat flux map (G), and sensible heat flux map (H) calculated in step S3, the latent heat flux (LE) is calculated. According to the energy balance equation: Rn = G + H + LE, the formula for calculating the latent heat flux is: LE = Rn - G - H. Traverse the entire image to obtain the latent heat flux map, with the unit of W / m 2 .

[0145] According to the air temperature (Ta) in the meteorological interpolation data obtained in step S14, the latent heat of vaporization (λ) of water is calculated. The latent heat of vaporization refers to the energy required for a unit mass of water to change from the liquid state to the gaseous state. The empirical formula is used: λ = 2.501 - 0.002361 * Ta, where Ta is in degrees Celsius and the unit of λ is MJ / kg. Traverse the temperature layer in the entire meteorological interpolation data to obtain the latent heat of vaporization map.

[0146] According to the latent heat flux map (LE) and the latent heat of vaporization map (λ) obtained in the above two steps, the initial evapotranspiration (ETi) is calculated. Evapotranspiration refers to the amount of water evaporated and transpired per unit area per unit time, usually in millimeters per hour (mm / h). The calculation formula is: ETi = LE / (λ * ρw) * 3600, where ρw is the density of water (taking 1000 kg / m 3 ), and 3600 is the coefficient for converting seconds to hours. Traverse the entire image to obtain the initial evapotranspiration map, with the unit of mm / h.

[0147] Use the apparent thermal inertia map (ATI) calculated in step S3 to assist in selecting cold pixels, and use the initial evapotranspiration map (ETi) obtained in the previous step to assist in selecting hot pixels. Cold pixels refer to pixels that are fully irrigated and the evapotranspiration reaches the potential evapotranspiration, and hot pixels refer to pixels that are completely dry and the evapotranspiration is zero. First, perform a histogram statistics on the ATI map, and select pixels with higher ATI values (for example, the top 5%) as candidate cold pixels. Then, among the candidate cold pixels, further screen pixels with higher NDVI values (for example, greater than 0.7) as the final set of cold pixels. For hot pixels, perform a histogram statistics on the ETi map, and select pixels with ETi values close to zero (for example, the bottom 5%) as candidate hot pixels. Then, among the candidate hot pixels, further screen pixels with lower NDVI values (for example, less than 0.2) as the final set of hot pixels.

[0148] Based on the cold pixel set and hot pixel set obtained in the previous step, calculate the evapotranspiration fraction (EF). First, calculate the average latent heat flux of cold pixels (LEcold) and the average sensible heat flux of hot pixels (Hhot). Then, calculate the instantaneous evapotranspiration fraction: EF = LE / (Rn - G) = LE / (LE + H). For each pixel, calculate using the average latent heat flux of cold pixels and the average sensible heat flux of hot pixels: EF = (Rn - G - Hhot) / (LEcold - Hhot). Traverse the entire image to obtain the evapotranspiration fraction map.

[0149] Based on the evapotranspiration fraction map (EF) obtained in the previous step and the average evapotranspiration of cold pixels (ETcold), calculate the final evapotranspiration (ET). First, calculate the average initial evapotranspiration of cold pixels (ETcold) as an estimate of potential evapotranspiration. Then, according to the EF map, calculate the final evapotranspiration for each pixel: ET = EF * ETcold. Traverse the entire image to obtain the evapotranspiration distribution map with the unit of mm / h.

[0150] Preferably, step S4 includes the following steps:

[0151] Step S41: Obtain historical irrigation data; perform sample data preprocessing on the historical irrigation data, the initial moisture stress judgment map, and the evapotranspiration distribution map to obtain a sample data set;

[0152] Step S42: Perform multi-dimensional feature extraction on the sample data set to obtain a phenological feature data set;

[0153] Step S43: Select and train a model based on the feature data set to obtain a water shortage diagnosis model;

[0154] Step S44: Input the feature data set into the water shortage diagnosis model to generate a water shortage diagnosis index, and obtain the water shortage diagnosis index.

[0155] In the embodiments of the present invention, historical irrigation data for at least the past year in the research area is collected, including irrigation dates, irrigation amounts (unit: millimeters), and irrigation methods (such as sprinkler irrigation, drip irrigation). These data can be obtained from farm management records or irrigation systems. The historical irrigation data is spatiotemporally matched with the preliminary moisture stress map obtained in step S2 and the evapotranspiration distribution map obtained in step S3. For each irrigation event, its corresponding date and location are found, and the preliminary moisture stress level and evapotranspiration amount at that location are extracted. At the same time, soil moisture data (obtainable through in-situ measurement or remote sensing inversion) and crop yield data (measured when the crops are mature) during the same period are collected. These data are integrated with the preliminary moisture stress map and the evapotranspiration distribution map to construct a sample dataset. Each row of the sample dataset represents a sample and contains the following columns: irrigation date, irrigation amount, preliminary moisture stress level, evapotranspiration amount, soil moisture, and crop yield. The data is cleaned to remove missing values and outliers. The numerical data is standardized so that its mean is 0 and its standard deviation is 1.

[0156] Feature extraction is performed on the sample dataset obtained in the previous step to construct a phenological feature dataset. In addition to directly using the preliminary moisture stress level and evapotranspiration amount as features, the following features are also extracted:

[0157] 1. Cumulative evapotranspiration: Calculate the cumulative evapotranspiration within a period of time (such as one week) before each irrigation.

[0158] 2. Moisture deficit index: According to the preliminary moisture stress level, define a continuous moisture deficit index. For example: severe stress = 1, moderate stress = 0.75, mild stress = 0.5, sufficient moisture = 0.25.

[0159] 3. Temporal feature: Convert the irrigation date to the day of the year as a temporal feature.

[0160] 4. Historical irrigation feature: Extract the irrigation amount of the previous irrigation and the time interval from the previous irrigation to the current irrigation.

[0161] 5. Statistical values of meteorological features: Centered on the irrigation event, select three days before and after, and calculate the average temperature, maximum temperature, minimum temperature, average humidity, average wind speed, and total radiation amount within these seven days.

[0162] 6. Topographic features: Extract the slope, aspect, and topographic undulation obtained in step S2.

[0163] These features are combined with the original features (preliminary moisture stress level, evapotranspiration amount) to form a feature dataset.

[0164] Based on the feature dataset obtained in the previous step, a random forest regression model is selected as the water shortage diagnosis model. The random forest model consists of multiple decision trees, and each decision tree is trained based on a random subset of the feature dataset. The model output is the average of the prediction results of all decision trees. The feature dataset is divided into a training set (70%), a validation set (15%), and a test set (15%). The random forest model is trained using the training set, and the model parameters are adjusted, such as the number of trees (set to 100), the maximum depth (set to 5), the minimum number of samples in the leaf nodes (set to 2), etc. The performance of the model is evaluated using the validation set, and the optimal parameter combination is selected. The evaluation metrics include the root mean square error (RMSE) and the coefficient of determination (R 2 ). The generalization ability of the model is evaluated using the test set. Finally, a trained water shortage diagnosis model is obtained.

[0165] The feature dataset obtained in the previous step is input into the trained random forest model to predict the water shortage diagnosis index (WDI). The model outputs a continuous value between 0 and 1, indicating the severity of crop water shortage. 0 indicates sufficient moisture, and 1 indicates extreme water shortage. Multiply the WDI by 100 to convert it into a percentage form to more intuitively represent the degree of water shortage. Traverse the entire study area to obtain the WDI value of each pixel and generate a water shortage diagnosis index map.

[0166] As an example of the present invention, refer to Figure 2 shown. In this example, step S5 includes:

[0167] Step S51: Set the WDI threshold and classify it according to the water shortage diagnosis index and meteorological interpolation data to obtain a WDI classification map;

[0168] In the embodiment of the present invention, according to the water shortage diagnosis index (WDI) obtained in step S4 and the meteorological interpolation data (especially the rainfall forecast for the next few days) obtained in step S14, the WDI threshold is set and classified. Combining the water demand characteristics of local crops (for example, corn is more sensitive to moisture during the jointing stage) and historical irrigation experience, the following thresholds are determined:

[0169] If there is more than 20 mm of rainfall in the forecast for the next 3 days, then:

[0170] WDI > 85%: Severe water shortage, irrigation is required.

[0171] 65% < WDI ≤ 85%: Moderate water shortage.

[0172] 45% < WDI ≤ 65%: Mild water shortage.

[0173] WDI ≤ 45%: Sufficient moisture.

[0174] If the forecast rainfall for the next 3 days is less than 20 mm, then:

[0175] WDI > 80%: Severe water shortage, irrigation is required.

[0176] 60% < WDI ≤ 80%: Moderate water shortage.

[0177] 40% < WDI ≤ 60%: Mild water shortage.

[0178] WDI ≤ 40%: Sufficient water.

[0179] According to the above thresholds, the WDI map is divided into four levels: severe water shortage, moderate water shortage, mild water shortage, and sufficient water, and the WDI classification map is obtained.

[0180] Step S52: Divide the irrigation management areas according to the WDI classification map to obtain the irrigation management area division map;

[0181] In the embodiment of the present invention, according to the WDI classification map obtained in the previous step, the irrigation management areas are divided. The region growing algorithm is used to merge the pixels with the same WDI level and spatially adjacent into the same irrigation management area. Set the minimum management area area threshold (for example, 100 square meters), and the scattered areas smaller than this threshold will be merged into the adjacent larger areas. Considering the actual layout of the irrigation system (such as the coverage range of the sprinkler unit, the length of the drip irrigation belt), the management area boundary is smoothed to make it more in line with the actual irrigation operation. Output the irrigation management area division map, and each management area has a unique ID.

[0182] Step S53: Calculate the irrigation amount according to the phenological characteristic data and the irrigation management area division map to obtain the irrigation amount distribution table;

[0183] In the embodiment of the present invention, according to the phenological characteristic data set obtained in step S4 and the irrigation management area division map obtained in the previous step, the irrigation amount of each irrigation management area is calculated. The method based on the crop water requirement is adopted. First, according to the meteorological interpolation data (temperature, humidity, wind speed, solar radiation), the reference crop evapotranspiration (ET0) is calculated using the FAO Penman-Monteith formula. Then, according to the crop type and growth stage (which can be obtained from the phenological characteristic data set), the crop coefficient (Kc) is determined. Calculate the actual crop water requirement (ETc) of each management area: ETc = Kc * ET0. Considering the effective rainfall (Pe), Pe is calculated according to the rainfall in the meteorological interpolation data and the empirical formula. Finally, calculate the irrigation amount (IR) of each management area: IR = ETc - Pe. If the calculated IR is less than 0, then IR is set to 0. Store the ID of each management area and the corresponding irrigation amount (unit: millimeter) in the irrigation amount distribution table.

[0184] Step S54: Assign the irrigation amount information in the irrigation amount distribution table to the grid pixels according to the spatial positions in the irrigation management zoning map to obtain an irrigation prescription map;

[0185] In the embodiment of the present invention, an irrigation prescription map is generated based on the irrigation management zoning map and the irrigation amount distribution table obtained in the previous step. First, a blank grid image with the same size as the study area is created, and the pixel values are initialized to 0. Then, according to the irrigation management zoning map, the pixel positions corresponding to each management area are found. According to the irrigation amount distribution table, the irrigation amount value of each management area is assigned to the corresponding pixels. For example, if the irrigation amount of management area 1 is 30 mm, then the values of all pixels with ID 1 in the irrigation management zoning map are set to 30. Finally, an irrigation prescription map is obtained, and the pixel value of this map represents the irrigation amount (unit: mm) at each position. The irrigation prescription map is saved in the GeoTIFF format for importing into the irrigation control system.

[0186] Preferably, step S53 is specifically:

[0187] Calculate the potential evapotranspiration to obtain a potential evapotranspiration map according to the meteorological interpolation data;

[0188] Calculate the crop coefficient to obtain a crop coefficient map according to the phenological characteristic data and the vegetation index map;

[0189] Estimate the actual evapotranspiration to obtain an actual evapotranspiration map according to the evapotranspiration distribution map, the crop coefficient map and the potential evapotranspiration map;

[0190] Calculate the soil water holding capacity correction coefficient to obtain a soil water holding capacity correction coefficient map according to the apparent thermal inertia map;

[0191] Conduct a preliminary calculation of the irrigation amount to obtain a preliminary irrigation amount map according to the actual evapotranspiration map;

[0192] Use the soil water holding capacity correction coefficient map to correct the preliminary irrigation amount map to obtain a corrected irrigation amount map;

[0193] Conduct an irrigation amount distribution according to the corrected irrigation amount map and the irrigation management zoning map to obtain an irrigation amount distribution table.

[0194] In the embodiment of the present invention, the potential evapotranspiration (ET0) is calculated according to the meteorological interpolation data (temperature, humidity, wind speed, solar radiation) obtained in step S14. The FAO Penman-Monteith formula is adopted, which comprehensively considers the energy balance and aerodynamic factors. The specific formula is:

[0195] ET0 = (0.408 * Δ * (Rn - G) + γ * (900 / (T + 273)) * u2 * (es - ea)) / (Δ + γ * (1 + 0.34 * u2));

[0196] Among them, ET0 is the potential evapotranspiration (mm / day), Δ is the slope of the saturation vapor pressure-temperature curve (kPa / °C), Rn is the net radiation (MJ / m 2 / day), G is the soil heat flux (MJ / m 2 / day, which can be ignored for daily-scale calculation), γ is the psychrometric constant (kPa / °C), T is the average air temperature (°C), u2 is the wind speed at 2 m height (m / s), es is the saturation vapor pressure (kPa), and ea is the actual vapor pressure (kPa). All meteorological parameters are obtained from meteorological interpolation data and unit conversions are performed. Traversing the entire study area, a potential evapotranspiration map is obtained, with the unit of mm / day.

[0197] Based on the phenological characteristic dataset obtained in step S4 and the vegetation index map (NDVI) obtained in step S2, the crop coefficient (Kc) is calculated. The crop coefficient reflects the water demand differences of different crops at different growth stages. Using the piecewise linear function method, the crop growth process is divided into four stages: the initial growth stage, the rapid growth stage, the mid-growth stage, and the maturity stage. According to the phenological characteristic data, the growth stage of the current crop is determined. For each stage, the Kc value is determined based on the NDVI value. For example:

[0198] Initial growth stage: Kc = Kc_ini + (NDVI - NDVI_ini) * (Kc_mid - Kc_ini) / (NDVI_mid - NDVI_ini);

[0199] Rapid growth stage: Kc = Kc_mid;

[0200] Mid-growth stage: Kc = Kc_mid;

[0201] Maturity stage: Kc = Kc_mid + (NDVI - NDVI_mid) * (Kc_end - Kc_mid) / (NDVI_end - NDVI_mid);

[0202] Among them, Kc_ini, Kc_mid, and Kc_end are the crop coefficient values for the initial growth stage, mid-growth stage, and maturity stage respectively, and NDVI_ini, NDVI_mid, and NDVI_end are the NDVI thresholds for these three stages. These parameters are determined according to the local crop types and planting experience. Traversing the entire NDVI map and the phenological characteristic dataset, a crop coefficient map is obtained.

[0203] Estimate the actual evapotranspiration (ET) based on the evapotranspiration distribution map obtained in step S3, the crop coefficient map calculated in step S53, and the potential evapotranspiration map. The calculation formula is: ET = Kc * ET0. This formula assumes that the actual evapotranspiration is equal to the product of the crop coefficient and the potential evapotranspiration. The calculation result generates an actual evapotranspiration map with the unit of mm / day.

[0204] Calculate the soil water holding capacity correction coefficient (Ksw) based on the apparent thermal inertia map (ATI) obtained in step S3. ATI reflects the ability of the soil to resist temperature changes and is closely related to soil water content. The linear function method is used: Ksw = a * ATI + b, where a and b are empirical coefficients determined according to the local soil type and measured data. For example, for sandy soil, a = 0.002 and b = 0.5; for clayey soil, a = 0.001 and b = 0.7. The higher the ATI value, the stronger the soil water holding capacity and the larger the Ksw value. Traverse the entire ATI map to obtain the soil water holding capacity correction coefficient map.

[0205] Conduct a preliminary calculation of the irrigation amount based on the actual evapotranspiration map to obtain a preliminary irrigation amount map. The actual evapotranspiration map is the evapotranspiration distribution map obtained in step S3. Preliminary calculation of irrigation amount: Crop water deficit = Actual evapotranspiration amount. Traverse the entire image to obtain the preliminary calculation of the irrigation amount and obtain a preliminary irrigation amount map.

[0206] Use the soil water holding capacity correction coefficient map (Ksw) obtained in the previous step to correct the preliminary irrigation amount map. The correction formula is: IR_corrected = IR_preliminary * Ksw, where IR_preliminary is the preliminary irrigation amount and IR_corrected is the corrected irrigation amount. This formula indicates that in areas with stronger soil water holding capacity, the irrigation amount can be appropriately reduced. Traverse the entire image to obtain the corrected irrigation amount map.

[0207] Based on the corrected irrigation amount map obtained in the previous step and the irrigation management zoning map obtained in step S52, allocate the irrigation amount. For each irrigation management area, calculate the average value of the corrected irrigation amounts of all pixels within its area as the irrigation amount of this management area. Store the ID of each management area and the corresponding irrigation amount (unit: millimeter) in the irrigation amount allocation table.

[0208] Therefore, from any perspective, the embodiments should be regarded as exemplary and non-restrictive. The scope of the present invention is defined by the appended claims rather than the above description. Therefore, all changes falling within the meaning and scope of the equivalent elements of the application documents are intended to be encompassed within the present invention.

[0209] The above are only specific embodiments of the present invention, enabling those skilled in the art to understand or implement the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features invented herein.

Claims

1. A method for inverting irrigation prescription maps based on UAV spectral data, characterized in that: The following steps are involved: Step S1: collect drone data in the farmland area to obtain multispectral images, thermal infrared images, DEM data and meteorological interpolation data; Step S2: Calculate micro-relief parameters based on multispectral images, DEM data and UAV GPS data, and construct a terrain-spectral coupling model to obtain a terrain relief map and a terrain-spectral coupling model; use the terrain-spectral coupling model to analyze the multidimensional state of vegetation and obtain a vegetation index map, a canopy temperature map and a preliminary water stress map; Step S3: Calculate net radiation using meteorological interpolation data to obtain a net radiation map; estimate soil heat flux based on the net radiation map, vegetation index map, and terrain relief map to obtain a soil heat flux map; calculate apparent thermal inertia based on the thermal infrared image and canopy temperature map, and estimate sensible heat flux to obtain an apparent thermal inertia map and a sensible heat flux map; perform evapotranspiration inversion based on the soil heat flux map, apparent thermal inertia map, and sensible heat flux map to obtain an evapotranspiration distribution map; Step S4: extracting phenological characteristic data of the water stress preliminary judgment map and the evapotranspiration distribution map, and determining the water shortage diagnosis index according to the phenological characteristic data; Step S5: Calculate irrigation management and draw up an irrigation prescription map based on the phenological characteristic data and the water shortage diagnosis index.

2. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: Step S1 includes the following steps: Step S11: at least 10 ground control points are evenly distributed in the farmland, and the coordinates of the control points are measured using an RTK GNSS receiver, with horizontal and vertical accuracy better than 2 cm. The UAV route is planned to be an east-west strip flight, with the altitude set to 100 meters, a heading overlap of 80%, and a lateral overlap of 70%. The UAV data is collected from the farmland area by the UAV, and meteorological data is obtained synchronously to obtain original spectral data, original digital elevation model data, UAV GPS data, and original meteorological station data; Step S12: performing spectral data preprocessing on the UAV spectral data to obtain a multispectral image and a thermal infrared image; Step S13: performing DEM data preprocessing on the digital elevation model data to obtain DEM data; Step S14: Perform meteorological data interpolation processing on the original meteorological station data to obtain meteorological interpolation data.

3. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: The micro-relief parameter calculation described in step S2 and the construction of the terrain-spectrum coupling model are specifically as follows: Micro-relief parameters were calculated for the DEM data. According to the elevation difference in the horizontal and vertical directions, the angle between the slope normal vector and the due north direction was calculated in degrees, ranging from 0 to 360 degrees, where 0 degrees represented due north, 90 degrees represented due east, 180 degrees represented due south, and 270 degrees represented due west. A 5x5 pixel window was defined, and the difference between the maximum and minimum elevation values ​​in the window was calculated as the terrain relief value of the central pixel of the window in meters, and the slope map, aspect map, and terrain relief map were obtained. Calculate solar geometry parameters based on UAV GPS data, slope map and aspect map to obtain solar altitude angle map, solar azimuth angle map and incident angle map; Calculate the direct radiation correction coefficient according to the solar altitude angle diagram, solar azimuth angle diagram and incident angle diagram to obtain the total direct radiation correction coefficient diagram; Calculate the scattered radiation correction coefficient according to the slope map and the aspect map to obtain the scattered radiation correction coefficient map; Calculate the proximity effect correction coefficient of the multispectral image to obtain a proximity effect correction coefficient map; According to the total direct radiation correction coefficient map, the scattered radiation correction coefficient map and the proximity effect correction coefficient map, a sub-band radiation transfer model is constructed to obtain a terrain-spectrum coupling model.

4. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: The multi-dimensional state analysis of vegetation described in step S2 is specifically as follows: The spectral index is calculated using multispectral images to obtain the initial spectral index; the spectral index is corrected using the terrain-spectral coupling model to obtain the vegetation index map and the moisture index map; The canopy temperature is calculated using thermal infrared images and vegetation index maps to obtain a canopy temperature map; A preliminary judgment of water stress is made based on the vegetation index map, moisture index map and canopy temperature map to obtain a preliminary judgment map of water stress.

5. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: The soil heat flux estimation described in step S3 is specifically: Calculate vegetation coverage according to the vegetation index map to obtain a vegetation coverage map; Calculate the terrain correction coefficient according to the terrain relief map to obtain a terrain correction coefficient map; Perform vegetation coverage correction calculation on the vegetation coverage map and the terrain correction coefficient map to obtain a corrected vegetation coverage map; The ratio of soil heat flux to net radiation is estimated using the corrected vegetation coverage map and net radiation map, and the G / Rn ratio map is obtained, where G is the soil heat flux and Rn is the net radiation value. The soil heat flux is calculated according to the G / Rn ratio diagram to obtain the soil heat flux diagram.

6. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: The calculation of apparent thermal inertia and estimation of sensible heat flux described in step S3 are specifically as follows: Extract temperature data of different phases of thermal infrared images to obtain temperature data sets; Apparent thermal inertia is calculated based on meteorological interpolation data and temperature data set to obtain the apparent thermal inertia diagram; The sensible heat flux is estimated based on the canopy temperature map, meteorological interpolation data and vegetation index map to obtain the sensible heat flux map.

7. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: The evapotranspiration inversion described in step S3 is specifically as follows: The latent heat flux is calculated based on the net radiation diagram, the apparent thermal inertia diagram and the sensible heat flux diagram to obtain the latent heat flux diagram; The latent heat of vaporization is calculated based on the air temperature in the meteorological interpolation data to obtain a latent heat of vaporization map; The initial evapotranspiration is calculated according to the latent heat flux diagram and the vaporization latent heat diagram to obtain the initial evapotranspiration diagram; The apparent thermal inertia map is used to assist in selecting cold pixels, and the initial evapotranspiration map is used to assist in selecting hot pixels, thereby obtaining a cold pixel set and a hot pixel set; Calculate the average latent heat flux of cold pixels and the average sensible heat flux of hot pixels, calculate the evapotranspiration scale factor, and obtain the evapotranspiration scale factor map; The average evapotranspiration of cold pixels is calculated, and the final evapotranspiration calculation is performed according to the evapotranspiration scale factor map to obtain the evapotranspiration distribution map.

8. The irrigation prescription map inversion method based on UAV spectral data according to claim 1 is characterized in that: Step S4 includes the following steps: Step S41: acquiring historical irrigation data; performing sample data preprocessing on the historical irrigation data, the water stress preliminary judgment map and the evapotranspiration distribution map to obtain a sample data set; Step S42: extracting multidimensional features from the sample data set to obtain a phenological feature data set; Step S43: performing model selection and training according to the feature data set to obtain a water shortage diagnosis model; Step S44: inputting the characteristic data set into the water shortage diagnosis model to generate a water shortage diagnosis index to obtain the water shortage diagnosis index.

9. The irrigation prescription map inversion method based on UAV spectral data according to claim 1, characterized in that: Step S5 includes the following steps: Step S51: performing WDI threshold setting and classification according to the water shortage diagnosis index and meteorological interpolation data to obtain a WDI classification diagram; Step S52: Divide the irrigation management areas according to the WDI classification diagram to obtain an irrigation management zoning map; Step S53: Calculate the irrigation amount according to the phenological characteristic data and the irrigation management zoning map to obtain an irrigation amount allocation table; Step S54: assigning grid pixels to the irrigation quantity information in the irrigation quantity allocation table according to the spatial position of the irrigation management zoning map to obtain an irrigation prescription map.

10. The irrigation prescription map inversion method based on UAV spectral data according to claim 9, characterized in that: Step S53 is specifically as follows: Potential evapotranspiration is calculated based on meteorological interpolation data to obtain a potential evapotranspiration map; The crop coefficient is calculated based on the phenological characteristic data and the vegetation index map to obtain the crop coefficient map; According to the evapotranspiration distribution map, crop coefficient map and potential evapotranspiration map, the actual evapotranspiration is estimated to obtain the actual evapotranspiration map; The soil water holding capacity correction coefficient is calculated according to the apparent thermal inertia diagram to obtain the soil water holding capacity correction coefficient diagram; The irrigation amount is preliminarily calculated based on the actual evapotranspiration map to obtain a preliminary irrigation amount map; The soil water holding capacity correction coefficient map is used to correct the irrigation quantity on the preliminary irrigation quantity map to obtain a corrected irrigation quantity map; the irrigation quantity is allocated according to the corrected irrigation quantity map and the irrigation management zoning map to obtain an irrigation quantity allocation table.

Images