Water and fertilizer management control method based on image data processing

Abstract

The application discloses a water and fertilizer management control method based on image data processing, samples each sampling unit of a to-be-controlled land plot, obtains a synchronous nutrient change curve and an extreme weather influence factor, and constructs a soil-crop-climate three-dimensional reference library; performs radiation correction, spectral calibration and real-time correction of spectral drift on a hyperspectral sensor; generates a dynamic collection time table, collects crop spectral data; automatically identifies an optimal sensitive wave band based on a random forest model, extracts spectral waveform features of the optimal sensitive wave band; compares the spectral waveform features with data of a same land plot in a designated period in the reference library in time sequence, calculates a nutrient change rate, constructs a two-dimensional diagnostic matrix of static content+dynamic trend, and obtains a nutrient deficiency type and a nutrient gap amount; constructs a multi-objective optimization function based on the nutrient deficiency type and the nutrient gap amount, calculates target fertilizer application amounts of each sampling unit and a fertilizer formula to obtain a fertilization scheme and executes the fertilization scheme, and realizes transformation from experience fertilization to data-driven precision fertilization.

Description

Technical Field

[0001] This invention belongs to the field of image data processing technology, and in particular relates to a water and fertilizer management and control method based on image data processing. Background Technology

[0002] In agricultural production, fertilization is a crucial step in ensuring crop yield and quality. However, traditional fertilization methods have long suffered from three major drawbacks: low precision, poor efficiency, and heavy pollution, making them unable to meet the needs of sustainable development in modern agriculture.

[0003] Nutrient diagnosis is lagging and crude: Traditional soil and crop nutrient testing relies on manual sampling and laboratory analysis. Not only is the sampling cycle long, but it also often adopts a plot-level average sampling mode. For example, only 1-2 mixed samples are collected from a 10-acre plot, which leads to the problem of local nutrient deficiency not being replenished and local nutrient eutrophic waste. According to statistics, the spatial adaptation error of nutrients under traditional fertilization methods can reach 30%-50%.

[0004] Insufficient reliability of spectral data: In existing UAV spectral detection technology, hyperspectral sensors are easily affected by environmental factors such as temperature, humidity, and light, resulting in spectral drift. For example, for every 5°C change in temperature, the spectral wavelength deviation can reach 0.5-1nm. Furthermore, the lack of standardized radiometric correction and spectral calibration procedures leads to unstable correlation between spectral data and nutrient content, making it difficult to support precise fertilization decisions.

[0005] Fertilization programs lack dynamic adaptability: Traditional fertilization programs are mostly based on empirical formulas for fixed growth periods, such as uniformly applying 20 kg / mu of nitrogen fertilizer during the jointing stage of wheat. They do not take into account the physiological rhythms of crops (such as differences in nutrient requirements during peak photosynthesis periods), historical nutrient change trends (such as the weekly decline in nitrogen in the same period over the past 3 years), and the impact of extreme weather (such as nitrogen leaching loss after heavy rain). As a result, the fertilizer utilization rate is only 30%-40%, while also causing environmental problems such as soil acidification and eutrophication of water bodies.

[0006] To address the aforementioned issues, this invention proposes to provide a precision fertilization control method that integrates multi-source data collaboration, precise spectral detection, and dynamic scheme optimization, thereby enabling a shift from experience-based fertilization to data-driven precision fertilization. Summary of the Invention

[0007] The purpose of this invention is to provide a water and fertilizer management control method based on image data processing, so as to solve the technical problems existing in the background art.

[0008] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0009] A water and fertilizer management control method based on image data processing includes the following steps:

[0010] S1: Divide the plot to be controlled into sampling units, obtain soil data and crop leaf nutrient data based on the sampling units, obtain the nutrient change curves of the same period in the past 3 years and the extreme weather impact factors of each sampling unit, and construct a three-dimensional benchmark library of soil-crop-climate.

[0011] S2: Perform radiometric correction and spectral calibration on the hyperspectral sensor carried by the UAV, and perform real-time correction of spectral drift based on the real-time spectral drift compensation mechanism;

[0012] S3: Based on the optimal acquisition window of crop physiological rhythms and historical data in the three-dimensional benchmark library, automatically generate a dynamic acquisition schedule and acquire crop spectral data;

[0013] S4: Based on the random forest model, the optimal sensitive bands under different crop-soil combinations are automatically identified, and the spectral waveform features of the optimal sensitive bands are extracted by combining continuous wavelet transform.

[0014] S5: Compare the spectral waveform characteristics with the previous period data of the same plot in the benchmark library, calculate the nutrient change rate, construct a two-dimensional diagnostic matrix of static content + dynamic trend, and obtain the nutrient deficiency type and nutrient deficit amount.

[0015] S6: Construct a multi-objective optimization function based on the type of nutrient deficiency and the amount of nutrient deficit. Set three optimization objectives: nutrient deficit filling, fertilizer utilization rate, and environmental risk. Calculate the target fertilizer application amount and fertilizer formula for each sampling unit to obtain the fertilization plan and execute it.

[0016] Preferably, the specific process of acquiring soil data and crop leaf nutrient data based on the sampling unit in step S1 is as follows:

[0017] S11: Each sampling unit adopts a five-point sampling method: at the center point of the unit and four equidistant points from the center point, soil from the 0-20cm topsoil layer is collected using a soil auger. Approximately 200g of soil is collected at each point to obtain the nitrogen, phosphorus, potassium content and pH value of the soil at each point.

[0018] S12: Within each sampling unit, randomly select 10-15 representative plants that are healthy and free from pests and diseases; collect functional leaves from each plant, one leaf from each plant, for a total of 10-15 leaves to form a mixed sample, and obtain the nutrient data of the mixed sample:

[0019] S13: Input the soil data and crop leaf nutrient data of each sampling unit into the database, and associate them with the GPS coordinates and sampling date of the sampling unit.

[0020] Preferably, the specific process for constructing the three-dimensional soil-crop-climate reference library in step S1 is as follows:

[0021] S14: Obtain multi-source basic datasets, including nutrient data from sampling units, extreme weather factor data, and historical management data;

[0022] S15: Perform data cleaning and spatiotemporal alignment on the multi-source basic dataset, and generate nutrient change curves;

[0023] S16: Using the Pearson correlation coefficient method, we analyzed the correlation between extreme weather factors and nutrient changes, the correlation between fertilizer application and crop leaf nutrient content, and the correlation between soil type and nutrient retention capacity, and screened core factors that significantly affect nutrient changes; we added feature tags to each sampling unit, including soil type tags, nutrient sensitivity tags, weather impact tags, and fertilizer suitability tags.

[0024] S17: Construction of a three-dimensional benchmark database of soil-crop-climate, including soil layer: storing soil texture, pH value, organic matter content, and nutrient change curves of the past 3 years for sampling units; crop layer: storing crop varieties, growth period, leaf nutrient change curves, and nutrient-yield response parameters; climate layer: storing extreme weather factors, weather-nutrient influence coefficients, and optimal sampling windows for the same period in the past 3 years.

[0025] Preferably, the specific process of performing radiometric correction and spectral calibration on the hyperspectral sensor carried by the UAV in step S2, and performing real-time spectral drift correction based on the real-time spectral drift compensation mechanism, is as follows:

[0026] S21: Perform radiometric correction: Using the integrating sphere light source as the standard radiation source, turn on the integrating sphere and output five different gradients of standard radiance values ​​in sequence; under each gradient, the sensor collects 10 sets of spectral data and calculates the average radiance of each set of data.

[0027] A radiation response curve was constructed by plotting the standard radiance of the integrating sphere on the x-axis and the measured average radiance of the sensor on the y-axis. A correction equation was obtained by fitting the curve using linear regression.

[0028] L true =a １ *L measured +b ;

[0029] in L true This represents the true radiance. L measured To measure the actual radiance, a １ This is the gain coefficient. b This is the offset coefficient;

[0030] S22: Using a mercury argon lamp as the standard spectral source, the sensor is aligned with the standard light source, the standard light source is turned on, the sensor collects complete emission spectrum data, identifies the position of characteristic peaks in the spectral curve, and records the measured wavelength value of each characteristic peak; the deviation between the measured wavelength and the standard wavelength is compared to perform band-by-band calibration; after calibration, the calibration parameters are saved to the sensor control unit, and data is collected according to the calibrated wavelength coordinates during flight;

[0031] S23: A standard gray diffuse reflectance plate is fixed to the top of the UAV fuselage. Before flight, the reflectance data of the gray plate under the current environment is measured on the ground using a portable spectrometer. This data serves as the reference value for real-time compensation and is denoted as S23. R standard ;

[0032] While the drone flies along a preset route, the hyperspectral sensor simultaneously collects spectral data of the crops in the plot and, at 0.1 seconds, synchronously collects reflectance data of a standard gray board, which is recorded as follows: R measured Real-time calculation of the deviation rate between the measured reflectance of the gray board and the reference value:

[0033] δ =| R measured - R standard | / R standard *100%;

[0034] Drift correction and data correction:

[0035] If the deviation rate δ If the value is ≤2%, it is determined that there is no significant spectral drift, and the current correction parameters remain unchanged;

[0036] If the deviation rate δ If the value is greater than 2%, a spectral shift is determined, and a correction coefficient k is calculated using a compensation algorithm. R standard / R measured This coefficient is then applied to crop spectral data collected concurrently to achieve drift correction.

[0037] Corrected spectral reflectance R corrected = R crop *k ,in R crop This represents the measured reflectance of the crop.

[0038] Preferably, the specific process of automatically generating a dynamic data acquisition schedule in step S3 based on the optimal acquisition window of crop physiological rhythms and historical data in the three-dimensional benchmark library is as follows:

[0039] S31: Extraction of Crop Physiological Rhythm Features and Identification of Key Time Periods:

[0040] Physiological rhythm parameter acquisition: retrieve the basic physiological parameters of the crop variety of the current plot from the crop layer of the three-dimensional benchmark database, including: peak photosynthetic period, key nodes of the growth period, and avoidance period;

[0041] Prioritization of critical time periods: Based on crop physiological characteristics, the collection time periods are divided into optimal, collectable, and avoidable levels, and different weights are assigned to them, including an optimal weight of 1.0, a collectable weight of 0.6, and an avoidable weight of 0.

[0042] S32: Historical Data Mining and Optimal Window Filtering for 3D Benchmark Library

[0043] Historical data retrieval: Extract historical spectral data and meteorological data for the same period of the past three years from the climate layer and crop layer of the three-dimensional benchmark database;

[0044] Historically optimal window: Data collection periods that meet the following conditions in the same period of the past 3 years are selected as historically optimal windows through data analysis: stable meteorological conditions and ≥80% overlap with the peak photosynthetic period of crops;

[0045] S33: Spatiotemporal matching of circadian rhythm periods with historically optimal windows:

[0046] Time dimension matching: The optimal and harvestable time periods identified by crop physiological rhythms are aligned with the historical optimal windows mined from the three-dimensional benchmark database to filter out candidate time periods;

[0047] Spatial dimension adaptation: Based on the spatial characteristics of sampling units in the soil layer of the three-dimensional benchmark database, the sampling time period is adjusted separately for areas sensitive to nutrient variation;

[0048] S34: Generate a basic timetable with the number of days in the reproductive period as the horizontal axis and the daily candidate time periods as the vertical axis;

[0049] UAV flight parameter linkage configuration: Automatically match UAV flight parameters according to the data collection period in the schedule.

[0050] Preferably, the specific process of automatically identifying the optimal sensitive bands for different crop-soil combinations based on the random forest model in step S4 is as follows:

[0051] S41: Retrieve basic data of the target plot from the three-dimensional benchmark library, including sampling unit information of different crop-soil combination types, and corresponding measured values ​​of soil and crop leaf nutrients; extract full-band spectral data collected by UAV, divide the bands at 1nm intervals, and form a raw spectral dataset with full-band features.

[0052] S42: Random Forest Model Construction and Training Model Parameter Initialization Setting Random Forest Core Parameters: Set the number of decision trees to 500, and the maximum depth of decision trees to None; use the normalized reflectance of 601 bands as input features, and the content of nitrogen, phosphorus, and potassium as output labels;

[0053] S43: Calculate the band feature importance scores for the three nutrients nitrogen, phosphorus, and potassium respectively, and obtain three band-importance scores; Threshold screening: Set a threshold for feature importance scores, and screen out bands with scores higher than the threshold as candidate sensitive bands; Calculate the Pearson correlation coefficient between candidate bands, remove highly redundant bands, and retain bands with strong information complementarity; Output the optimal sensitive band set for different crop-soil combinations.

[0054] Preferably, the specific process of extracting the spectral waveform features of the optimal sensitive band by combining continuous wavelet transform in step S4 is as follows:

[0055] S44: Extract the spectral reflectance curve data of the optimal sensitive band interval selected in step S4 from the full-band spectral data collected by the UAV; and use the wavelet threshold noise reduction method to eliminate random noise in the spectral reflectance curve.

[0056] S45: Parameter configuration and time-frequency domain conversion for continuous wavelet transform;

[0057] S45: Spectral waveform feature quantization based on wavelet coefficients: Extract three types of core waveform features that are strongly correlated with crop nutrients from the time-frequency domain wavelet coefficient map to form a spectral waveform feature set of sensitive bands.

[0058] Preferably, the specific process of step S5 is as follows:

[0059] S51: Retrieve the previous three historical data from the same plot and crop growth period from the three-dimensional benchmark library: including the spectral waveform feature set of the previous three monitoring periods, the corresponding measured values ​​of soil and crop leaf nutrients, fertilization records and meteorological data.

[0060] Extract the spectral waveform feature set of the current period, i.e. the spectral waveform feature set in step S45. Using the GPS coordinates of the sampling unit as the unique identifier, match the current feature set, the feature sets of the previous 3 periods, and the measured nutrient values ​​one by one to form the time-series data chain of each sampling unit.

[0061] S52: Convert the spectral waveform characteristic values ​​of the current cycle and the previous 3 cycles into corresponding nutrient estimation values;

[0062] S53: For each sampling unit, calculate the rate of change of nitrogen, phosphorus, and potassium, mark the abnormal rate intervals, and construct a two-dimensional diagnostic matrix of static content + dynamic trend;

[0063] S54: Construct a 3×3 two-dimensional diagnostic matrix to map the nutrient status of each sampling unit to the corresponding quadrant of the matrix, and accurately distinguish the type of nutrient deficiency;

[0064] S55: Precise Classification of Nutrient Deficiency Types Based on the quadrant results of the two-dimensional diagnostic matrix and combined with the soil type labels of the three-dimensional benchmark library, the nutrient deficiency type and nutrient deficit amount of each sampling unit are finally confirmed.

[0065] Preferably, the specific process of step S6 is as follows:

[0066] S61: Construct a multi-objective optimization function, including nutrient gap filling rate, fertilizer utilization rate and environmental risk index, and set constraints, including fertilizer application amount constraints and fertilizer formula constraints;

[0067] S62: Solve the multi-objective optimization function: Generate an initial population based on the nutrient gap and constraints, calculate fitness, generate the next generation population through selection, crossover, and mutation operations, and converge to obtain the optimal solution, i.e., the optimal solution, after 50 iterations.

[0068] S63: Calculate the target fertilization amount and fertilizer formula based on the optimal solution.

[0069] The beneficial effects of this invention include:

[0070] 1. Improve the accuracy of nutrient diagnosis and achieve precise small-scale adaptation: By constructing a three-dimensional benchmark database of soil-crop-climate through a five-point sampling method, the spatiotemporal correlation of nutrient data of sampling units is realized, accurately identifying nutrient hotspots (such as local high phosphorus areas) and cold spots (such as micro-regional potassium deficiency), reducing the error of nutrient spatial adaptation; the two-dimensional diagnostic matrix of static content + dynamic trend, combined with the nutrient change rate of the same period in the past 3 years, can distinguish between short-term nutrient deficiency (such as rapid nitrogen deficiency after heavy rain) and long-term nutrient imbalance (such as slow consumption of phosphorus for 3 consecutive periods), improving the accuracy of nutrient deficiency type determination and providing a precise basis for differentiated fertilization.

[0071] 2. Ensuring the reliability of spectral data: A three-level correction system—radiative correction, spectral calibration, and real-time spectral drift compensation—effectively eliminates environmental interference and improves the signal-to-noise ratio of spectral data through integrating sphere linear regression correction, mercury-argon lamp band-by-band calibration, and real-time compensation with a standard gray card. A random forest + continuous wavelet transform feature extraction strategy can automatically identify the optimal sensitive bands for different crop-soil combinations.

[0072] 3. Optimize the dynamic adaptability of fertilization programs to improve fertilizer utilization and environmental safety: Based on crop physiological rhythms and the historical optimal window of the three-dimensional benchmark library, a dynamic acquisition schedule is automatically generated to ensure that spectral acquisition is carried out during the most sensitive period of nutrient response. At the same time, the multi-objective optimization function can design formulas differently according to the type of nutrient deficiency, thereby improving fertilizer utilization. The environmental risk index is quantitatively calculated through factors such as nitrogen leaching and phosphorus loss. Combined with weather forecasts, the formula is dynamically adjusted to reduce nitrogen and phosphorus loss and effectively reduce the risk of soil acidification and water pollution. Attached Figure Description

[0073] Figure 1 This is a flowchart illustrating the water and fertilizer management control method based on image data processing according to the present invention.

[0074] Figure 2 This is a schematic diagram of the process for constructing a three-dimensional reference library of soil, crop, and climate according to the present invention.

[0075] Figure 3 This is a schematic diagram of the architecture of the random forest model of the present invention. Detailed Implementation

[0076] The following is in conjunction with the appendix Figures 1-3 The present invention will be further described in detail below:

[0077] Example 1

[0078] See appendix Figure 1 As shown, a water and fertilizer management control method based on image data processing includes the following steps:

[0079] S1: Divide the control plots into sampling units, collect soil samples and crop leaf samples based on the sampling units, obtain soil data including nitrogen, phosphorus, potassium content and pH value, and crop leaf nutrient data including nitrogen, phosphorus, and potassium content, respectively, obtain the nutrient change curves of the same period in the past 3 years and the impact factors of extreme weather for each sampling unit, and construct a three-dimensional benchmark library of soil-crop-climate by combining the plot's historical fertilization records, crop yield data over the years and soil type database.

[0080] S2: Perform radiometric correction and spectral calibration on the hyperspectral sensor carried by the UAV to eliminate instrument errors, and perform real-time correction of spectral drift based on the real-time spectral drift compensation mechanism. Simultaneously collect the reflectance of the standard gray board during the flight of the UAV to correct the spectral drift caused by the sensor due to changes in temperature and humidity in real time.

[0081] S3: Based on the optimal acquisition window, including crop physiological rhythms during photosynthetic peak periods and historical data from a three-dimensional benchmark database, a dynamic acquisition schedule is automatically generated. For example, during the wheat jointing stage, acquisitions are conducted twice daily, from 9:00-10:00 and 15:00-16:00. Multispectral data is acquired based on this dynamic acquisition schedule. A drone equipped with a hyperspectral camera (400-1000nm) collects crop spectral data (spectral reflectance), which reflects nutrient content.

[0082] S4: Based on the random forest model, the optimal sensitive bands under different crop-soil combinations are automatically identified. For example, 680nm and 760nm correspond to nitrogen deficiency in rice, and 550nm and 820nm correspond to phosphorus deficiency in maize. The spectral waveform features of the optimal sensitive bands are extracted by continuous wavelet transform (CWT) to enhance the ability to capture weak nutrient signals.

[0083] S5: Compare the spectral waveform characteristics with the data from the previous three periods of the same plot in the benchmark library to calculate the nutrient change rate, such as the weekly decrease in nitrogen, and construct a two-dimensional diagnostic matrix of static content + dynamic trend to accurately distinguish between short-term nutrient deficiency and long-term nutrient imbalance, and determine the nutrient deficiency type and corresponding nutrient deficit of each sampling unit.

[0084] S6: Construct a multi-objective optimization function based on the type of nutrient deficiency and the amount of nutrient deficit, set three optimization objectives: nutrient deficit filling, fertilizer utilization rate, and environmental risk, and calculate the target fertilizer application amount and fertilizer formula for each sampling unit, such as the ratio of available nitrogen to slow-release nitrogen.

[0085] The specific process of acquiring soil data and crop leaf nutrient data based on the sampling unit in step S1 is as follows:

[0086] S11: Each sampling unit adopts a five-point sampling method: at the center point of the unit and four equidistant points from the center point, soil from the 0-20cm topsoil layer is collected using a soil auger. This depth is the main distribution layer of crop roots. Approximately 200g of soil is collected at each point to obtain the nitrogen, phosphorus, potassium content and pH value of the soil at each point.

[0087] Nitrogen content: The total nitrogen content of the soil was determined by the Kjeldahl method, and the alkaline nitrogen (available nitrogen) content was determined by the alkaline hydrolysis diffusion method. The results are expressed in mg / kg.

[0088] Phosphorus content: The available phosphorus content in the soil was determined by the molybdenum-antimony colorimetric method, and the results are expressed in mg / kg.

[0089] Potassium content: The available potassium content in the soil was determined by flame photometry, and the results are expressed in mg / kg.

[0090] pH value: Soil acidity and alkalinity were determined by potentiometric method (water-to-soil ratio 2.5:1), and the results were retained to one decimal place;

[0091] S12: Within each sampling unit, randomly select 10-15 representative plants that are healthy and free from pests and diseases; collect functional leaves from the plants, such as the flag leaf of wheat, the third leaf from the top of rice, and the middle leaf of the unfolded leaf of corn, collecting one leaf from each plant, for a total of 10-15 leaves to form a mixed sample, and obtain the nutrient data of the mixed sample:

[0092] Nitrogen content: determined by the Kjeldahl method, and the results are expressed as a percentage of dry weight (%).

[0093] Phosphorus content: determined by the vanadium molybdenum yellow colorimetric method, and the results are expressed as percentages (%) of dry weight.

[0094] Potassium content: determined by flame photometry, and the results are expressed as a percentage of dry weight (%).

[0095] S13: Input the soil data (nitrogen, phosphorus, potassium content and pH value) and crop leaf nutrient data (nitrogen, phosphorus, potassium content) of each sampling unit into the database, and associate them with the GPS coordinates, sampling date and other information of the sampling unit to provide basic data support for the subsequent construction of a three-dimensional benchmark database of soil-crop-climate.

[0096] See Figure 2 As shown, the specific process of constructing the three-dimensional soil-crop-climate reference library in step S1 is as follows:

[0097] S14: Obtain the multi-source base dataset:

[0098] Nutrient data collection for sampling units: Soil sample test data (nitrogen, phosphorus, potassium content and pH value) and crop leaf sample test data (nitrogen, phosphorus, potassium content) from the same period of the past 3 years (the growth period consistent with the current sampling time) of each sampling unit are retrieved and classified and labeled according to sampling unit number-year-growth period.

[0099] Extreme weather factor collection: By connecting with publicly available meteorological data or monitoring data from small field weather stations, key extreme weather indicators affecting crop nutrient absorption over the same period in the past three years are extracted, including:

[0100] Extreme temperatures, such as high temperatures >35℃ or low temperatures <5℃ for more than 3 consecutive days;

[0101] Abnormal precipitation, such as rainstorms with a single rainfall of more than 50 mm, or drought with no effective precipitation for more than 7 consecutive days;

[0102] Severe convective weather, such as strong winds and hail; match the corresponding weather events according to the region-year-fertility period of the sampling unit, and mark the time of occurrence, duration and intensity level of the weather.

[0103] Historical management data collection: Collect the historical fertilization records of the plots for the past 3 years, including fertilization time, fertilizer type (fast-acting or slow-release), nutrient ratio, and fertilization amount; retrieve crop yield data (yield per unit area) and quality data (such as grain protein content) for the same period; integrate soil type database information to clarify the basic attributes of each sampling unit, such as soil texture (clay, loam, or sandy soil), soil organic matter content, and soil fertilizer and water retention capacity;

[0104] S15: Preprocessing and spatiotemporal matching of multi-source base datasets:

[0105] Data cleaning: Remove abnormal data, including samples with detection values ​​outside the reasonable range or missing records; normalize and transform nutrient data obtained from different detection methods; quantify and encode extreme weather factors, such as assigning drought levels to 1-5.

[0106] Spatiotemporal alignment: Using the GPS coordinates of the sampling unit as the spatial anchor point and the crop growth period as the temporal anchor point, the nutrient data, extreme weather factors, fertilization records, and yield data of the same sampling unit in the same period of the past 3 years are matched one by one to form a three-dimensional data element of space-time-indicator.

[0107] Nutrient variation curve generation: Using the year as the horizontal axis and nutrient content as the vertical axis, plot the soil nitrogen, phosphorus and potassium content variation curve and the crop leaf nitrogen, phosphorus and potassium content variation curve for each sampling unit. Calculate the interannual nutrient variation rate by the curve slope and mark the extreme weather events or fertilization adjustment behaviors corresponding to the curve fluctuation nodes.

[0108] S16: Perform data association and feature extraction:

[0109] Factor correlation analysis: The Pearson correlation coefficient method was used to analyze the correlation between extreme weather factors and nutrient changes (such as the positive correlation between rainstorms and soil nitrogen leaching), the correlation between fertilizer application and crop leaf nutrient content, and the correlation between soil type and nutrient retention capacity. Core factors that have a significant impact on nutrient changes, i.e. factors with a confidence level ≥ 95%, were screened.

[0110] Feature tag generation: Add feature tags to each sampling unit, including soil type tags, nutrient sensitivity tags, weather impact tags, and fertilizer adaptation tags. Nutrient sensitivity tags include nitrogen-sensitive and phosphorus-inefficient types, weather impact tags include those that are easily washed away by heavy rain, and fertilizer adaptation tags include those that are adapted to slow-release nitrogen.

[0111] S17: Construction of a Three-Dimensional Reference Database for Soil-Crop-Climate

[0112] The database employs a hierarchical architecture, divided into three core layers: Soil layer: stores soil texture, pH value, organic matter content, and nutrient variation curves over the past three years for each sampling unit. Crop layer: stores crop variety, growth stage, leaf nutrient variation curves, and nutrient-yield response parameters. Climate layer: stores extreme weather factors for the same period over the past three years, weather-nutrient influence coefficients, and optimal sampling windows (e.g., 9:00-11:00 AM without extreme weather).

[0113] Data entry and indexing: Standardized multi-source data are entered into the database in a hierarchical structure. A relationship is established with the GPS coordinates of the sampling unit as the unique index, enabling the query function to retrieve the full-dimensional data of soil, crop and climate of the corresponding unit by inputting the coordinates.

[0114] Benchmark threshold setting: Based on the statistical analysis results of 3 years of data, benchmark values ​​and upper and lower limit thresholds of soil and crop leaves at each growth stage are set as reference standards for subsequent fertility diagnosis.

[0115] After completing a planting cycle, replenish the seasonal nutrient data, weather data, fertilization data, and yield data, update the nutrient change curve and response parameters, and remove outdated data (such as historical data older than 5 years) to ensure the timeliness of the benchmark database.

[0116] Example 2

[0117] Based on Example 1, the specific process of performing radiometric correction and spectral calibration on the hyperspectral sensor carried by the UAV in step S2, and performing real-time spectral drift correction based on the real-time spectral drift compensation mechanism, is as follows:

[0118] S21: Perform radiation correction.

[0119] Using an integrating sphere light source as the standard radiation source, with known radiance and accuracy conforming to ISO9001 standards, the hyperspectral sensor was fixed on the optical test platform, ensuring that the sensor lens was directly facing the light outlet of the integrating sphere, with the distance controlled at 50cm. Environmental control: the test was conducted in a laboratory environment with constant temperature (25℃±2℃) and constant humidity (50%±5%RH) to avoid radiation response deviations caused by temperature fluctuations.

[0120] The integrating sphere is activated, and five standard radiance values ​​with different gradients are output sequentially, covering the response range of the sensor's working wavelength band of 400-1000nm. The sensor collects 10 sets of spectral data under each gradient, and the average radiance of each set of data is calculated.

[0121] A radiation response curve was constructed by plotting the standard radiance of the integrating sphere on the x-axis and the measured average radiance of the sensor on the y-axis. A correction equation was obtained by fitting the curve using linear regression.

[0122] L true =a １ *L measured +b ;

[0123] in, L true This represents the true radiance. L measured To measure the actual radiance, a １ This is the gain coefficient. b This is the offset coefficient;

[0124] S22: Perform spectral calibration:

[0125] A mercury-argon lamp was used as the standard spectral source. The characteristic peak wavelengths of its emission spectrum are known, such as 546.07nm, 576.96nm, and 696.54nm. The sensor was aligned with the standard light source to ensure that the characteristic peak light rays were perpendicularly incident on the sensor lens. The sensor sampling frequency was adjusted to match the emission frequency of the light source to avoid signal distortion. The spectral resolution was set to 1nm.

[0126] Wavelength calibration and correction:

[0127] The standard light source is activated, the sensor collects complete emission spectrum data, identifies the positions of characteristic peaks in the spectral curve, and records the measured wavelength value of each characteristic peak.

[0128] By comparing the deviation between the measured wavelength and the standard wavelength, such as the measured value of 545.87nm when the standard peak is 546.07nm (deviation of -0.2nm), a wavelength deviation correction table is constructed to correct the wavelength across the entire sensor band (400-1000nm) band by band, ensuring that the wavelength error of each band is ≤±0.1nm.

[0129] After calibration, save the calibration parameters to the sensor control unit, and collect data according to the calibrated wavelength coordinates during flight;

[0130] S23: Perform dynamic compensation during UAV flight, i.e., real-time correction of spectral drift:

[0131] A diffuse reflectance standard gray plate is fixed to the top of the drone fuselage. The reflectance of this standard gray plate is known, for example, an average reflectance of 50% in the 400-1000nm wavelength band with a deviation ≤±2%. Before flight, the reflectance data of the gray plate under the current environment is measured on the ground using a portable spectrometer. This data serves as the reference value for real-time compensation and is denoted as [reference value]. R standard ;

[0132] The spectral drift compensation algorithm is preloaded, and the compensation trigger frequency is set to be consistent with the sensor sampling frequency, such as 10Hz, which means 10 compensations per second to ensure real-time performance.

[0133] Synchronous data acquisition and comparison: While the drone flies along the preset route, the hyperspectral sensor simultaneously acquires the reflectance data of a standard gray board every 0.1 seconds, while also collecting spectral data of the crops in the plot. This data is recorded as follows: R measured Real-time calculation of the deviation rate between the measured reflectance of the gray board and the reference value:

[0134] δ =| R measured - R standard | / R standard *100%;

[0135] Drift correction and data correction:

[0136] If the deviation rate δ If the value is ≤2%, it is determined that there is no significant spectral drift, and the current correction parameters remain unchanged;

[0137] If the deviation rate δ If the value is greater than 2%, a spectral shift is determined to exist (caused by changes in temperature and humidity). A correction coefficient k is calculated using a compensation algorithm. R standard / R measured This coefficient is then applied to crop spectral data collected concurrently to achieve drift correction.

[0138] Corrected spectral reflectance R corrected = R crop *k ,in R crop This represents the measured reflectance of the crop.

[0139] Step S3, which automatically generates a dynamic data acquisition schedule based on the optimal acquisition window of crop physiological rhythms and historical data in the three-dimensional benchmark library, is as follows:

[0140] S31: Extraction of Crop Physiological Rhythm Features and Identification of Key Time Periods:

[0141] Physiological rhythm parameter acquisition: retrieve the basic physiological parameters of the crop varieties in the current plot from the crop layer of the three-dimensional reference library, including: peak photosynthetic period, which for most crops is 9:00-11:00 am and 15:00-16:00 pm on sunny days. At this time, the crop photosynthetic efficiency is high and the canopy spectrum is most sensitive to changes in nutrients.

[0142] At key growth stages, such as the jointing stage of wheat and the tillering stage of rice, the nutrient requirements of crops at different growth stages vary greatly, and their spectral characteristics differ significantly.

[0143] Avoid peak hours, such as the period of direct sunlight at noon (12:00-14:00), when the spectrum is easily affected by high light intensity; and in the early morning when the dew is still wet, the moisture in the leaves will interfere with the spectral reflectance.

[0144] Prioritization of key time periods: Based on crop physiological characteristics, the collection time periods are divided into the optimal level (peak photosynthesis period, with the highest signal-to-noise ratio of spectral data), the collectable level (non-peak but interference-free periods, such as 8:00-9:00 am), and the avoidable level (strong light, dew, high temperature and high humidity periods, during which collection is prohibited), and different weights are assigned to them, with the optimal level having a weight of 1.0, the collectable level having a weight of 0.6, and the avoidable level having a weight of 0.

[0145] S32: Historical Data Mining and Optimal Window Filtering for 3D Benchmark Library

[0146] Historical data retrieval: Extract historical spectral data and meteorological data (temperature, humidity, light intensity) for the same period in the past 3 years (same growth stage) from the climate layer and crop layer of the 3D baseline database.

[0147] Historical Optimal Window Discovery: Through data analysis, data collection periods that meet the following conditions in the same period of the past 3 years are selected as historical optimal windows: stable meteorological conditions (temperature 15℃-28℃, humidity 40%-60%, no strong winds or rain); and overlap with the peak photosynthetic period of crops ≥80%.

[0148] Based on historical data, a correlation between meteorological factors and spectral accuracy was established to clarify the threshold of the impact of different meteorological conditions on spectral data. For example, when the light intensity is >120000lx, the spectral deviation rate will exceed 8%, providing a quantitative basis for subsequent window selection.

[0149] S33: Spatiotemporal matching of circadian rhythm periods with historically optimal windows:

[0150] Time dimension matching: The optimal and collectable time periods identified by crop physiological rhythms are aligned with the historical optimal windows mined from the three-dimensional benchmark database to screen out candidate time periods that both conform to physiological rhythms and have been verified for collection accuracy in historical data.

[0151] Example: During the jointing stage of wheat, the physiologically optimal time is 9:00-11:00 and 15:00-16:00; historical data shows that there has been no extreme weather during this period in the past 3 years, so these two periods are directly listed as candidate windows.

[0152] Spatial dimension adaptation: Based on the spatial characteristics of sampling units in the soil layer of the three-dimensional benchmark database, for nutrient variation sensitive areas, such as plot edges and water-fertilizer transition zones, the sampling units have been densified to 5m×5m, and the collection time period has been adjusted separately: priority is given to collection during the optimal time period to ensure data accuracy in the densified areas; for ordinary sampling units, the data can be flexibly allocated between the optimal and collectable time periods to balance accuracy and collection efficiency.

[0153] S34: Dynamic Acquisition Schedule Generation and Parameter Configuration:

[0154] Basic timetable framework construction: Using the number of days in the reproductive period as the horizontal axis and the daily candidate time period as the vertical axis, a basic timetable is generated, including the total number of days of data collection, covering the key monitoring cycles of the current reproductive period, such as continuous data collection for 5 days during the wheat jointing stage;

[0155] Daily data collection time slots, such as two optimal time slots per day, with the duration of each collection session calculated based on the plot area;

[0156] Sampling units are allocated with priority given to sensitive encrypted areas during the optimal time period, while ordinary areas are allocated in batches.

[0157] UAV flight parameter linkage configuration: Automatically match UAV flight parameters according to the collection period in the schedule. For example, during the optimal period when the light is sufficient, set the flight altitude to 80m and the speed to 4m / s; during the period when the light is slightly weaker, appropriately reduce the flight speed to 3m / s and extend the sensor exposure time to ensure data quality.

[0158] Example 3

[0159] Based on Example 1 or Example 2, see Figure 3 As shown, the specific process of automatically identifying the optimal sensitive bands under different crop-soil combinations based on the random forest model in step S4 is as follows:

[0160] S41: Sample Data Extraction

[0161] Retrieve basic data of the target plot from the three-dimensional benchmark library: including sampling unit information of different crop-soil combination types, corresponding measured values ​​of soil and crop leaf nutrients (nitrogen, phosphorus, and potassium content), where different crop-soil combination types can be rice-loam or corn-clay.

[0162] Extract full-band spectral data collected by the UAV: ​​cover the complete working band of the sensor (400-1000nm), divide the bands at 1nm intervals, and form an original spectral dataset containing 601 band features.

[0163] Based on the correspondence between spectral band characteristics and nutrient content, a labeled sample set was constructed. Each sample contains 601 band reflectance feature values ​​and 3 nutrient content label values ​​(nitrogen, phosphorus, and potassium content). The number of single crop-soil combination samples is ≥300.

[0164] S42: Random Forest Model Parameter Initialization Settings. Core parameters for random forest: number of decision trees set to 500, maximum decision tree depth set to None (depth is automatically determined by the data), feature sampling method is random sampling. n 1 / 2 One characteristic, n The total number of bands is given, and the sample is collected using sampling with replacement. A multi-output random forest model is constructed: the normalized reflectance of 601 bands is used as the input feature, and the content of three nutrients (nitrogen, phosphorus, and potassium) is used as the output label to achieve simultaneous identification of sensitive bands for the three nutrients in a single modeling process.

[0165] S43: Optimal sensitive band selection based on feature importance.

[0166] Band feature importance calculation: Extract the feature importance score of each band in the random forest model: This score reflects the contribution of a single band to the nutrient content prediction result. The higher the score, the stronger the correlation between the band and the nutrient content.

[0167] The importance scores of the band features for the three nutrients nitrogen, phosphorus and potassium were calculated separately, and three band-importance score lists were generated.

[0168] The calculation process for the importance scores of the band features corresponding to the three nutrients, nitrogen, phosphorus, and potassium, is as follows:

[0169] For a given decision tree, a certain band feature f At the node t When used for splitting, the reduction in impurity (using the Gini coefficient as an example) is calculated using the following formula:

[0170] ;

[0171] in, Pre-split node t Gini impurity, c This refers to the number of nutrient content categories, such as classifying nitrogen content into three categories: deficient, critical, and sufficient. p ti For nodes t The middle belongs to the first i The proportion of samples containing different nutrients; k The number of child nodes generated after a split is the number of nodes; decision trees are typically binary splits. k =2, n v child node tv The number of samples in n parent node t The total number of samples in the sample; Gini ( t v ) is a child node t v Gini impurity.

[0172] The greater the reduction in impurity, the greater the contribution of that band to reducing the predicted impurity of nutrients when it splits at that node.

[0173] Calculation of the importance score of a single band in a single decision tree: The importance score of a band f in a single decision tree is the sum of the reductions in impurity of that band at all split nodes.

[0174] ;

[0175] in T f Use bands in this decision tree f The set of all nodes to be split.

[0176] Calculation of the final importance score for a single nutrient corresponding to a single band in a random forest:

[0177] A random forest contains M decision trees, for a given band. f The final Gini importance score for the corresponding nitrogen (or phosphorus, potassium) nutrient is the average importance score of that band across all decision trees:

[0178] .

[0179] Among them, nutrient These refer to nitrogen, phosphorus, and potassium, respectively. Each of the three nutrients needs to be calculated separately, which means constructing three independent random forest models or evaluating the characteristic contribution of each nutrient separately in a multi-output model.

[0180] Optimal sensitive band selection rules:

[0181] Threshold screening: Set a threshold for feature importance scores, which can be the top 20 percentile, and select bands with scores higher than the threshold as candidate sensitive bands.

[0182] Correlation deduplication: Calculate the Pearson correlation coefficient between candidate bands, remove highly redundant bands (correlation coefficient ≥ 0.9), and retain bands with strong information complementarity.

[0183] Optimal sensitive band set output:

[0184] For different crop-soil combinations, a specific set of optimal sensitive bands is output. See the table below for examples:

[0185]

[0186] The specific process of extracting the spectral waveform features of the optimal sensitive band by combining continuous wavelet transform in step S4 is as follows:

[0187] S44: Extract the spectral reflectance data of the optimal sensitive band range selected in step S4 from the full-band spectral data collected by the UAV:

[0188] If the optimal sensitive wavelengths for nitrogen in paddy soil are 685nm and 755nm, then a complete spectral curve containing a continuous range (e.g., 680-760nm) of these two wavelengths needs to be extracted, rather than two isolated wavelength points, to ensure the continuity of the waveform. The curves are categorized by sampling unit, with each sampling unit corresponding to one spectral reflectance curve for a sensitive wavelength range, denoted as R(λ), where λ is the wavelength.

[0189] Spectral data noise reduction preprocessing:

[0190] A wavelet thresholding denoising method is used to eliminate random noise in the spectral curve, such as sensor electronic noise and ambient stray light interference. Specifically, the following steps are performed: A db4 wavelet basis is selected, and a single-level wavelet decomposition is applied to the spectral curve of the sensitive band to obtain approximation coefficients (low frequencies, corresponding to spectral trends) and detail coefficients (high frequencies, corresponding to noise). A soft threshold is set, and the portion of the detail coefficients below the threshold is set to zero, preserving the high-frequency information reflecting the true spectral characteristics. Wavelet reconstruction is then performed on the processed approximation coefficients and detail coefficients to obtain the denoised smooth spectral curve R'(λ), laying the foundation for subsequent waveform feature extraction.

[0191] S45: Perform parameter configuration and time-frequency domain conversion for continuous wavelet transform (CWT);

[0192] The key to continuous wavelet transform is selecting wavelet basis and scaling parameters that fit the spectral waveform characteristics. Specific configurations are as follows:

[0193] Wavelet basis function: Morlet wavelet is selected, a complex-valued wavelet that combines time-domain and frequency-domain localization properties, and can accurately capture the peak and valley positions of the spectral waveform. Its expression is:

[0194] ;

[0195] in, a For scale parameters, For translation parameters, is the conjugate function of the Morlet wavelet basis. The result of the transformation is a complex value. Its magnitude reflects the energy intensity of the spectral signal at that scale and translation position, and its phase reflects the changing trend of the spectral waveform.

[0196] With scale parameters a Using the vertical axis as the ordinate, wavelength λ as the horizontal axis, and wavelet coefficient modulus as the color depth, a time-frequency domain spectrum of the sensitive wavelength band is generated. The darker the color in the spectrum, the more significant the spectral waveform characteristics at that wavelength position, allowing for intuitive identification of nutrient-related feature locations.

[0197] S45: Spectral waveform feature quantification based on wavelet coefficients: From the time-frequency domain wavelet coefficient spectrum, three types of core waveform features strongly correlated with crop nutrients are extracted to form a "sensitive band waveform feature set", as shown in the table below:

[0198]

[0199] Feature optimization and output:

[0200] Constructing a standardized feature set: Feature normalization involves Z-score normalization of the extracted red-edge features, peak-valley features, and trend features to eliminate dimensional differences between different features. The formula is as follows:

[0201] F norm = ( F−μ ) / σ ;

[0202] in F These are the original eigenvalues. μ This is the mean of the feature across all sampling units. σ The standard deviation is denoted as .

[0203] Principal component analysis (PCA) was used to reduce the dimensionality of the normalized feature set, retaining principal components with a cumulative contribution rate of ≥90%, thus forming the final sensitive band spectral waveform feature set.

[0204] The optimal sensitive band and spectral waveform feature set of each sampling unit are associated with GPS coordinates and stored in a three-dimensional benchmark library. This data is then used directly as the core input feature for constructing the two-dimensional diagnostic matrix in step S5, achieving a seamless connection from band selection to feature extraction.

[0205] The specific process of step S5 is as follows:

[0206] S51: Retrieve the historical data of the previous three periods for the same plot and crop growth stage from the three-dimensional benchmark library: including the spectral waveform feature set (red edge position, red edge slope, peak-to-valley difference, etc.) of the previous three monitoring periods, the corresponding measured values ​​of soil and crop leaf nutrients, fertilization records and meteorological data.

[0207] Extract the spectral waveform feature set of the current period (output of step S4), and match the current feature set, the feature sets of the previous 3 periods, and the measured nutrient values ​​one by one according to the GPS coordinates of the sampling unit to form a time-series data chain for each sampling unit.

[0208] S52: Feature-Nutrient Correlation Mapping is based on the response relationship between spectral waveform features and nutrient content in a three-dimensional benchmark library. The spectral waveform feature values ​​of the current period and the previous three periods are converted into corresponding nutrient estimates (nitrogen, phosphorus, and potassium content). For example, for every 1 nm shift of the red edge position towards longer wavelengths, the corresponding leaf nitrogen content increases by 0.2%; for every 0.05 increase in peak-to-valley difference, the corresponding soil available phosphorus content increases by 5 mg / kg. Combined with the measured nutrient values ​​from the previous three periods, the estimation bias is corrected to obtain the precise time-series sequence of "period-nutrient content" for each sampling unit, denoted as... C t , t =0 indicates the current period. t =1,2,3 represent the first three periods.

[0209] S53: Quantitative Calculation of Nutrient Change Rate:

[0210] Calculate the nutrient change between adjacent periods: Δ C t-1，t = C t - C t-1 ;

[0211] Calculate the nutrient change rate based on the monitoring cycle interval, such as 7 days / period:

[0212] v =Δ C t-1，t / Δ T ;

[0213] Where, Δ T The period interval is in days; a negative rate indicates a decrease in nutrients, while a positive rate indicates an accumulation of nutrients.

[0214] For each sampling unit, the rate of change of nitrogen, phosphorus, and potassium is calculated separately, and abnormal rate intervals are marked. For example, if the weekly decrease of nitrogen is >0.3%, it is determined that the rate of decrease is too fast.

[0215] The following is a two-dimensional diagnostic matrix constructed using static content and dynamic trend analysis:

[0216]

[0217] S54: Construct a 3×3 two-dimensional diagnostic matrix to map the nutrient status of each sampling unit to the corresponding quadrant of the matrix, accurately distinguishing the type of nutrient deficiency:

[0218] Static content / Dynamic trend Stable trend slow downward trend rapid downward trend adequate Nutrients are sufficient and stable Potential long-term imbalance risk Short-term nutrient deficiency precursors critical Critical steady state Long-term nutrient imbalance (slow depletion) Short-term nutrient deficiency (rapid depletion) lack Long-term nutrient imbalance (insufficient historical accumulation) Long-term nutrient imbalance (continuous consumption) Acute short-term nutrient deficiency

[0219] Further refine the judgment based on crop physiological characteristics: for example, if the nitrogen level in wheat drops rapidly during the jointing stage, it is judged as a short-term nutrient deficiency and topdressing is required immediately; if the nitrogen level is at a critical value and drops slowly for more than 3 periods, it is judged as a long-term nutrient imbalance and the base fertilizer ratio needs to be adjusted.

[0220] S55: Identification of Nutrient Deficiency Type and Calculation of Nutrient Shortage Amount:

[0221] Precise nutrient deficiency type classification: Based on the quadrant results of the two-dimensional diagnostic matrix and combined with soil type labels from the three-dimensional benchmark library, such as clay soil having strong fertilizer retention capacity and being less prone to short-term nutrient deficiency, the nutrient deficiency type of each sampling unit was finally confirmed:

[0222] Short-term nutrient deficiency: caused by recent factors, such as leaching from heavy rain or rapid absorption by crops, manifested as a rapid decline in dynamic trend, and critical or lack of static content.

[0223] Long-term nutrient imbalance: caused by long-term management or soil characteristics, such as insufficient base fertilizer and high soil phosphorus fixation rate, manifested as a slow dynamic decline or a stable low level, with static content remaining in the critical or deficient range for a long time.

[0224] Nutrient deficit quantitative calculation:

[0225] Determine the target nutrient content: based on the upper limit of suitable nutrient levels for the current growth stage in the three-dimensional baseline database. C target ;

[0226] Calculate the static nutrient deficit: Δ C static = C target - C 0;

[0227] Overlaying dynamic trend correction: For short-term nutrient deficiencies with a rapid downward trend, additional trend compensation is required. For example, if the weekly decrease is 0.3%, and fertilization is planned 10 days later, the compensation amount is 0.3% × 10 / 7 ≈ 0.43%. For long-term imbalances, the soil utilization rate coefficient needs to be considered. For example, if the phosphorus utilization rate of clay is 60%, the deficit amount needs to be divided by the utilization rate.

[0228] Final nutrient deficit: Δ C final =Δ C static +Δ C trend (Short-term nutrient deficiency)

[0229] or Δ C final =Δ C static / η(Long-term imbalance) η (For soil nutrient utilization rate).

[0230] Step S6 constructs a multi-objective optimization function for different nutrient deficiency types (short-term deficiency or long-term nutrient imbalance). This function aims to achieve precise nutrient replenishment while maximizing fertilizer utilization and minimizing environmental risks. The final output is the target fertilizer application rate and fertilizer formula for each sampling unit. The specific process is as follows:

[0231] S61: Set the nutrient gap fill rate f The formula for calculating 1 is:

[0232] f 1 = (Actual total amount of nutrients supplemented / Total nutrient deficit) × 100%;

[0233] The actual total amount of nutrients supplemented = fertilizer application amount × fertilizer nutrient content × soil utilization rate;

[0234] Total nutrient deficit = Δ C final ;

[0235] The primary optimization objective is: nutrient gap filling rate. f 1. Maximize;

[0236] Fertilizer utilization rate f The formula for calculating 2 is:

[0237] f 2 = (Nutrients absorbed by the crop / Nutrients input) × 100%;

[0238] Wherein, the amount of nutrients absorbed by the crop = the total nutrient deficit × the crop absorption coefficient

[0239] Nutrient input = Fertilizer application rate × Fertilizer nutrient content;

[0240] The second optimization objective is: fertilizer utilization rate. f 2. Maximize;

[0241] Environmental Risk Index f The formula for calculating 3 is:

[0242] f 3= α ×Nitrogen leaching risk+ β × Risk of phosphorus loss+ γ × Risk of soil acidification;

[0243] Among them, nitrogen leaching risk: when the probability of rainfall is greater than 60%, the risk value increases by 0.22 for every 10% increase in the proportion of available nitrogen;

[0244] Phosphorus loss risk: When the slope is greater than 15°, the risk value increases by 0.33 for every 5 kg / mu increase in phosphate fertilizer application.

[0245] α , β, γ These are weighting coefficients, assigned values ​​based on soil type, such as α=0.6 for sandy soil and β=0.5 for clay soil;

[0246] Third optimization objective: Environmental risk index f 3 Minimize.

[0247] Set constraints:

[0248] Fertilizer application rate constraint: 0≤ F i ≤ F max ;

[0249] F i For the first i Nutrient application rate F max This refers to the maximum safe amount of fertilizer to be applied during the crop's growth period. For example, the maximum amount of nitrogen fertilizer to be applied during the jointing stage of wheat is ≤20 kg / mu.

[0250] Fertilizer formulation constraints:

[0251] Short-term nutrient deficiency: fast-acting nutrients account for ≥70% (rapid fertilization), and slow-release nutrients account for ≤30%;

[0252] Long-term nutrient imbalance: slow-release nutrients account for ≥50% (long-term improvement), and fast-acting nutrients account for ≤50%.

[0253] The multi-objective optimization function is set as follows:

[0254] max [ f 1( x ), f 2( x )];

[0255] minf 3( x );

[0256] st0≤ F i ≤ F max And fertilizer formulation constraints;

[0257] x =[ F N , F p , F k ]

[0258] in, xAs decision variables, including nitrogen ( F N ),phosphorus( F p ), potassium ( F k ) Amount of fertilizer applied.

[0259] S62: Solve the multi-objective optimization function:

[0260] Initial population generation: 100 initial fertilization schemes are randomly generated based on the nutrient deficit and constraints;

[0261] Fitness calculation: Substitute into the optimization function to calculate the fitness of each solution. f 1. f 2. f Three values ​​are used to evaluate the merits of the proposed solution.

[0262] Genetic iterative optimization: The next generation of population is generated through selection, crossover, and mutation operations. After 50 iterations, the optimal solution is obtained by convergence.

[0263] S63: Based on the optimal solution, calculate the target fertilizer application rate and fertilizer formula:

[0264] F 总 =Δ C final × S / ( ω × η );

[0265] in, F 总 Δ is the target total fertilizer application amount. C final This represents the final nutrient deficit, where S is the weight of the topsoil, typically taken as 200,000 kg / mu (0-20 cm topsoil). ω This refers to the content of the target nutrient in the fertilizer. For example, urea contains 46% nitrogen. ω =0.46, η This refers to soil nutrient utilization rate.

[0266] Determine the type and ratio of fertilizer based on the type of nutrient deficiency and soil properties:

[0267] For short-term nutrient deficiencies, focus on fast-acting nutrients to quickly replenish crop needs; supplement with a small amount of slow-release nutrients to prevent nutrient deficiency later, such as: nitrogen fertilizer: 70% urea (fast-acting) + 30% slow-release urea; phosphorus fertilizer: 100% superphosphate (fast-acting); potassium fertilizer: 100% potassium chloride (fast-acting).

[0268] For long-term nutrient imbalance: focus on slow-release nutrients to improve soil nutrient reserves; add activators or fertilizer retainers to improve nutrient utilization, such as: nitrogen fertilizer: 60% slow-release urea + 40% urea + nitrification inhibitor; phosphate fertilizer: 50% calcium magnesium phosphate (slow-release) + 50% superphosphate + humic acid activator; potassium fertilizer: 100% slow-release potassium sulfate.

Claims

1. An image data processing-based water and fertilizer management control method, characterized by, The method comprises the following steps: S1: sample unit division is performed on the to-be-controlled plot, soil data and crop leaf nutrient data are obtained based on the sample units, a specified-year same-period nutrient change curve and an extreme weather influence factor of each sample unit are obtained, and a soil-crop-climate three-dimensional reference library is constructed; S2: the hyperspectral sensor carried by the unmanned aerial vehicle is subjected to radiation correction and spectral calibration, and spectral drift real-time correction is performed based on a spectral drift real-time compensation mechanism; S3: based on crop physiological rhythm and an optimal collection window of historical data in the three-dimensional reference library, a dynamic collection time table is automatically generated, and crop spectral data are collected; S4: based on a random forest model, optimal sensitive bands under different crop-soil combinations are automatically identified, and spectral waveform features of the optimal sensitive bands are extracted by using continuous wavelet transform; S5: the spectral waveform features are compared with the same-plot previous specified-period data in the reference library in time sequence, a nutrient change rate is calculated, a two-dimensional diagnosis matrix of static content + dynamic trend is constructed, a nutrient deficiency type and a nutrient gap amount are obtained; S6: a multi-objective optimization function is constructed based on the nutrient deficiency type and the nutrient gap amount, three optimization objectives of nutrient gap filling, fertilizer utilization rate and environmental risk are set, target fertilizer application amounts and fertilizer formulations of each sample unit are calculated, a fertilization scheme is obtained, and the fertilization scheme is executed; The specific process of step S5 is as follows: S51: the previous three periods of historical data of the same plot and the same crop growth period are called from the three-dimensional reference library, including the spectral waveform feature set of the previous three monitoring periods, the corresponding soil and crop leaf nutrient measured values, the fertilizer record and the meteorological data; The spectral waveform feature set of the current period is extracted, the current feature set, the previous three period feature sets and the nutrient measured values are matched one by one according to the sample unit GPS coordinates as the unique identifier, and the time sequence data chain of each sample unit is formed; S52: the spectral waveform feature values of the current period and the previous three periods are respectively converted into corresponding nutrient estimated values; S53: for each sample unit, the change rates of nitrogen, phosphorus and potassium are calculated respectively, the rate abnormal interval is marked, and a two-dimensional diagnosis matrix of static content + dynamic trend is constructed; S54: a 3*3 two-dimensional diagnosis matrix is constructed, the nutrient state of each sample unit is mapped to the corresponding quadrant of the matrix, and the nutrient deficiency type is accurately distinguished; S55: the nutrient deficiency type is accurately classified according to the quadrant result of the two-dimensional diagnosis matrix, and the soil type label of the three-dimensional reference library is combined to finally confirm the nutrient deficiency type and the nutrient gap amount of each sample unit; The two-dimensional diagnosis matrix comprises a static nutrient content dimension and a dynamic change trend dimension, the static nutrient content dimension is divided into three levels, including sufficient, critical and deficiency; and the dynamic change trend dimension is divided into three types based on the nutrient change rate, including a stable trend, a slow downward trend and a rapid downward trend.

2. The water and fertilizer management control method based on image data processing according to claim 1, characterized in that, The specific process of obtaining soil data and crop leaf nutrient data based on the sample units in step S1 is as follows: S11: five-point sampling is performed on each sample unit: 0-20cm plough layer soil is collected by using a soil drill at the center point and four equidistant points of the center point, 200g of soil is collected at each point, and the nitrogen, phosphorus and potassium contents and the pH value of the soil at each point are obtained; S12: In each sampling unit, 10-15 representative plants that grow healthily and have no pests and diseases are randomly selected; the functional leaves of the plants are collected, 1 leaf per plant, a total of 10-15 leaves to form a mixed sample, and the nutrient data of the mixed sample are obtained: S13: The soil data and crop leaf nutrient data of each sampling unit are recorded into the database, and the GPS coordinates of the sampling unit, the information of the sampling date are associated.

3. The water and fertilizer management control method based on image data processing according to claim 2, characterized in that, The specific process of constructing the soil-crop-climate three-dimensional reference library in step S1 is as follows: S14: Obtain a multi-source basic data set, including sampling unit nutrient data, extreme weather factor data and historical management data; S15: Data cleaning, time and space dimension alignment are performed on the multi-source basic data set, and a nutrient change curve is generated; S16: Using Pearson correlation coefficient method, the correlation between extreme weather factor and nutrient change, the correlation between fertilizer amount and crop leaf nutrient content, the correlation between soil type and nutrient retention capacity are analyzed, and the core factors that have significant influence on nutrient change are screened; a feature label is added to each sampling unit, including soil type label, nutrient sensitive type label, weather influence label, and fertilizer adaptation label; S17: Soil-crop-climate three-dimensional reference library construction, including soil layer: storing soil texture, pH value, organic matter content, and nutrient change curve in the past 3 years of sampling unit; crop layer: storing crop variety, growth period, leaf nutrient change curve, and nutrient-yield response parameter; climate layer: storing extreme weather factor in the past 3 years, weather-nutrient influence coefficient, and optimal collection window.

4. The water and fertilizer management control method based on image data processing according to claim 1, characterized in that, The specific process of radiation correction and spectral calibration of the hyperspectral sensor carried by the unmanned aerial vehicle in step S2, and real-time compensation of spectral drift based on spectral drift real-time compensation mechanism is as follows: S21: Radiation correction: Take the integrating sphere light source as the standard radiation source, turn on the integrating sphere, and output 5 different gradient standard radiation brightness values in turn; collect 10 groups of spectral data under each gradient, and calculate the average radiation brightness of each group of data; Take the integrating sphere standard radiation brightness as the abscissa and the measured average radiation brightness of the sensor as the ordinate to construct a radiation response curve, and use linear regression fitting to obtain the correction equation: L true =a 1 *L measured +b ； wherein L true is the true radiance, L measured is the measured radiance, a 1 is a gain coefficient, b is an offset coefficient; S22: Use mercury argon lamp as standard spectrum source, align the sensor with the standard light source, start the standard light source, the sensor collects the complete emission spectrum data, identifies the characteristic peak position in the spectrum curve, records the measured wavelength value of each characteristic peak; compare the deviation between the measured wavelength and the standard wavelength for wave band by wave band calibration; save the correction parameters to the sensor control unit, and collect data according to the corrected wavelength coordinates during flight; S23: Fix a diffuse reflection standard gray board on the top of the UAV body. Before flight, measure the reflectance data of the gray board under the current environment on the ground by the portable spectrometer, which is used as the reference value of real-time compensation, denoted as R standard ; When the UAV flies along the preset route, the hyperspectral sensor synchronously collects the reflectivity data of the standard gray board every 0.1 seconds while collecting the spectral data of the crops in the plot, and the data is recorded as R measured The deviation rate of the measured reflectivity of the gray board from the reference value is calculated in real time. δ =| R measured - R standard | / R standard *100%； Drift correction and data correction: If the deviation rate δ ≤ 2%, it is determined that there is no significant spectral drift, and the current correction parameters are kept unchanged; If the deviation rate δ > 2%, it is determined that there is spectral drift, and a correction coefficient is calculated by a compensation algorithm k = R standard / R measured , and the coefficient is applied to the crop spectrum data collected at the same time to correct the drift. Corrected spectral reflectance R corrected = R crop *k wherein R crop is the measured reflectance of the crop.

5. The water and fertilizer management control method based on image data processing according to claim 1, characterized in that, The specific process of automatically generating a dynamic collection schedule based on crop physiological rhythm and the optimal collection window of historical data in the three-dimensional reference library in step S3 is as follows: S31: Extraction of crop physiological rhythm characteristics and key period calibration: Physiological rhythm parameter collection: retrieve the basic physiological parameters of the current crop variety in the crop layer of the three-dimensional reference library, including: photosynthetic peak period, growth period key node and avoidable period; Key period priority division: According to the physiological characteristics of crops, the collection period is divided into optimal level, available level and avoid level, and different weights are given, including optimal level weight 1.0, available level weight 0.6, and avoid level weight 0; S32: Three-dimensional reference library historical data mining and optimal window screening: Historical data retrieval: Extract the historical spectral collection data of the same period in the past three years and the corresponding meteorological data from the three-dimensional reference library climate layer + crop layer; Historical optimal window: Through data analysis, the collection period that meets the following conditions in the same period in the past three years is selected as the historical optimal window: stable weather conditions, and the coincidence degree with the crop photosynthetic peak period is greater than or equal to 80%; S33: Time and space matching of physiological rhythm period and historical optimal window: Time dimension matching: Align the optimal level and available level periods of crop physiological rhythm calibration with the historical optimal window mined from the three-dimensional reference library, and select the candidate period; Space dimension adaptation: Combined with the spatial characteristics of the sampling unit in the soil layer of the three-dimensional reference library, the collection period of the nutrient variation sensitive area is adjusted separately; S34: Generate a basic time table with growth period days as the horizontal coordinate and daily candidate periods as the vertical coordinate; Unmanned aerial vehicle flight parameter linkage configuration: According to the collection period in the time table, automatically match the unmanned aerial vehicle flight parameters.

6. The water and fertilizer management control method based on image data processing according to claim 1, characterized in that, The specific process of automatically identifying the optimal sensitive waveband under different crop-soil combinations based on the random forest model in step S4 is as follows: S41: Retrieve the basic data of the target plot from the three-dimensional reference library, including the sampling unit information of different crop-soil combination types, the corresponding soil and crop leaf nutrient measured values; extract the full-band spectral data collected by the unmanned aerial vehicle, divide the bands by 1 nm interval, and form the original spectral data set with full-band features; S42: Random forest model construction and training model parameter initialization Set the core parameters of the random forest: the number of decision trees is 500, and the maximum depth of the decision tree is None; use the normalized reflectivity of 601 bands as input features, and use nitrogen, phosphorus and potassium as output labels; S43: Calculate the waveband feature importance score of nitrogen, phosphorus and potassium respectively, and obtain three waveband-importance scores; threshold screening: set the feature importance score threshold, select the wavebands with scores higher than the threshold as candidate sensitive wavebands; calculate the Pearson correlation coefficient between the candidate wavebands, remove highly redundant wavebands, and retain wavebands with strong information complementarity; output the optimal sensitive waveband set for different crop-soil combinations.

7. The water and fertilizer management control method based on image data processing according to claim 6, characterized in that, The specific process of extracting the spectral waveform features of the optimal sensitive waveband in step S4 combined with continuous wavelet transform is as follows: S44: From the full-band spectral data collected by the unmanned aerial vehicle, intercept the spectral reflectance curve data of the optimal sensitive waveband interval selected in step S4; and use the wavelet threshold denoising method to eliminate random noise in the spectral reflectance curve; S45: Parameter configuration and time-frequency domain conversion of continuous wavelet transform; S45: Quantification of spectral waveform features based on wavelet coefficients: From the time-frequency domain wavelet coefficient map, extract three types of core waveform features strongly related to crop nutrients to form a spectral waveform feature set of sensitive wavebands. 8.The water and fertilizer management control method based on image data processing of claim 1, wherein, The specific process of step S6 is as follows: S61: Construct a multi-objective optimization function including nutrient gap filling rate, fertilizer utilization rate and environmental risk index, and set constraint conditions including fertilizer application amount constraint and fertilizer formula constraint; S62: Solve the multi-objective optimization function: generate an initial population according to the nutrient gap amount and the constraint conditions, calculate the fitness, generate the next generation population through selection, crossover and mutation operations, and after 50 iterations, the optimal solution, i.e. the optimal scheme, is obtained by convergence; S63: Calculate the target fertilizer application amount and the fertilizer formula based on the optimal scheme.

Images