A variable fertilization prescription image generation method combined with periodicity and phenology

Abstract

The application relates to a variable fertilization prescription image generation method combining cycles and phenology, and relates to the technical field of fertilization control. The application is aimed at solving the problems of poor adaptability of an existing fertilization prescription image and low accuracy of a given fertilization decision. The application comprises the following steps: calling a phenology historical accumulation vegetation state vectorization module to generate a fertilization prescription image, counting the results of each partition of the fertilization prescription image and outputting a final fertilization prescription image; the input of the phenology historical accumulation vegetation state vectorization module comprises the following: a time sequence list of unmanned aerial vehicle multi-spectrum data, a ground wireless sensor data table, a crop phenology history table, a fertilization history table, a basic fertilization amount, a cluster number, a minimum fertilization amount and a maximum fertilization amount; the output is a prescription image generation output; the accurate fertilization amount of each pixel is obtained based on the fertilization prescription image generation output, and a final fertilization prescription image is obtained. The application is used for obtaining a fertilization prescription image.

Description

Technical Field

[0001] This invention relates to the field of fertilization control technology, and in particular to a method for generating variable fertilization prescription images that combines periodicity and phenology. Background Technology

[0002] Food security is a core strategic issue in my country's agricultural development. The application efficiency of chemical fertilizers such as nitrogen, phosphorus, and potassium directly affects crop yield and the degree of agricultural non-point source pollution. In traditional field production, the amount of chemical fertilizers applied far exceeds the actual needs of crops. Excessive fertilization not only wastes resources but also leads to a chain of environmental problems such as soil acidification and nitrate pollution in groundwater. Variable-rate fertilization technology is an effective way to solve these problems. Based on the differences in crop growth and soil nutrient distribution within a field, it generates differentiated fertilization requirements for different spatial locations, forming a spatially differentiated, pixel-by-pixel fertilization prescription image. This image is then executed by drones or precision fertilization equipment according to the generated prescription image. The fertilization prescription image is a single-band raster image aligned with the field, and the value of each pixel is a floating-point fertilization amount to be applied at the corresponding spatial location.

[0003] Current methods for generating variable fertilization prescription images mainly rely on the following technical means: 1. Fixed threshold zoning generation method: Based on experience, the field is divided into several fixed management zones, and a fertilization prescription image with uniform pixel values ​​is generated in each zone. However, this method relies on manual subjective division of management zones, which makes it difficult to reflect the heterogeneity of crop growth that changes dynamically with the seasons within the field. Furthermore, it cannot utilize multi-temporal data to perceive the changing trends of the same plot at different stages of the fertilization cycle. Once the zone boundaries are determined, they are not updated for the current season, resulting in the generated fertilization prescription images showing large areas of uniform value plateaus within the field, with poor pixel-by-pixel adaptability, thus failing to achieve precise fertilization. 2. Single-period UAV data index generation method: This method collects one period of UAV multispectral data to calculate vegetation indices such as NDVI, and generates fertilization level images pixel by pixel based on index levels. However, single-phase data reflects the instantaneous canopy state, which is easily affected by light and meteorological conditions at the time of collection, and cannot integrate historical information on crop phenology. Moreover, this method does not consider the key influences of soil moisture, electrical conductivity, and nitrogen metabolism indicators on fertilization effects, resulting in low accuracy of fertilization decisions based on fertilization prescription maps. Third, the static soil sampling and testing method obtains soil nutrient distribution data through field sampling and laboratory testing, and then uses spatial interpolation to generate prescription images as fertilization base maps. However, the number of sampling points is limited, and the spatial interpolation error is large, which leads to artifacts at the pixel level in the generated prescription images, making it difficult for fertilization prescription maps to accurately guide fertilization decisions. Meanwhile, the long nutrient testing cycle cannot reflect the dynamic changes in nutrients during the crop's growth period. Therefore, it lacks dynamic response capability to the differentiated needs of different fertilization windows (such as topdressing at the jointing stage and supplemental fertilization at the grain-filling stage), resulting in the generated fertilization prescription map being unable to adapt to the needs of different crop growth stages, leading to low accuracy in fertilization decisions based on the fertilization prescription map. Furthermore, soil test data does not include information on crop canopy growth, which leads to insufficient pixel perception of areas such as lodging, lesions, and growth stagnation, making it difficult to identify abnormal crop growth areas in the field, further resulting in low accuracy in fertilization decisions based on the fertilization prescription map. IV. A simple multi-source image overlay generation method directly combines or weights-averages UAV multispectral images and ground sensor outputs to create a fertilization prescription image. However, there are significant misalignments between UAV images and ground sensor images in terms of acquisition time and spatial resolution. Simple superposition will directly introduce the misalignment error into the prescription image generation process, resulting in misalignment deviations in the fertilization prescription image at the pixel level. At the same time, this method does not use the crop fertilization cycle to anchor the time-series data, treating the contribution weight of data from different phenological stages equally and ignoring the special value of data near the critical fertilization window. The cumulative impact of fertilization history on current growth is ignored, causing the generated prescription image to fail to reflect the compensation needs of areas with historical fertilization deficiencies at the pixel level. This can easily lead to local over-fertilization or under-fertilization, ultimately resulting in low accuracy of the fertilization decisions given by the fertilization prescription map. Summary of the Invention

[0004] This invention addresses the problems of poor adaptability of existing fertilizer prescription images and low accuracy of the resulting fertilization decisions by proposing a variable fertilizer prescription image generation method that combines periodicity and phenology.

[0005] A method for generating variable fertilization prescription images that combines periodicity and phenology, specifically:

[0006] The phenological history accumulation and vegetation state vectorization module JLBQModel is called to generate fertilization prescription images, the results of each partition of the fertilization prescription images are statistically analyzed, and the final fertilization prescription image is output.

[0007] The input JLBQInput of the phenological history accumulation vegetation state vectorization module includes: UAV multispectral data time series table UAVList, ground wireless sensor data table SensorList, crop phenological history table PhenocalList, fertilization history table FertHistList, base fertilizer amount BaseAmount, cluster number JLBQClusterK, minimum fertilizer amount PrescMin, and maximum fertilizer amount PrescMax; the output of JLBQModel is the prescription image generation output JLBQOutput;

[0008] The fertilizer prescription image PrescImg is the pixel fertilizer amount image JLBQPrescImg field of the JLBQOutput generated from the prescription image;

[0009] The process of statistically analyzing the results of each partition of the fertilizer prescription image and outputting the final fertilizer prescription image involves: obtaining the precise amount of fertilizer for each pixel based on the total amount of fertilizer applied to each region and the total number of pixels in each region of the fertilizer prescription image PrescImg, thus obtaining the final fertilizer prescription image.

[0010] Furthermore, the phenological history-accumulated vegetation state vectorization module JLBQModel is obtained in the following way:

[0011] S1. Construct UAVList, SensorList, PhenocalList, and FertHistList and their global control parameters;

[0012] S2. Use UAVList, PhenocalList, SensorList and global control parameters to obtain the air-to-ground joint generation input ZQMZInput, and establish a periodically anchored air-to-ground timing joint generation module ZQMZModel. The input of ZQMZModel is ZQMZInput, and the output is the air-to-ground joint generation output ZQMZOutput.

[0013] S3. Use UAVList, SensorList, PhenocalList, FertHistList and global control parameters to obtain the phenological accumulation input WHSJInput, construct the phenological history accumulation vegetation state vectorization module WHSJModel, the input of WHSJModel is WHSJInput, and WHSJModel calls ZQMZModel to obtain the phenological accumulation output WHSJOutput.

[0014] S4. Use WHSJOutput and global control parameters to obtain the prescription image generation input JLBQInput, construct the cluster label weighted prescription image generation module JLBQModel, the input of JLBQModel is JLBQInput, and JLBQModel calls WHSJModel to obtain the prescription image generation output JLBQOutput.

[0015] Furthermore, the UAVList, SensorList, PhenocalList, and FertHistList, along with the global control parameters, are specifically as follows:

[0016] Each record in UAVList contains the following fields: UAVDate: the data collection date number calculated from the sowing date of the current season; UAVImg: the multispectral data collected by the UAV; UAVBandList: a list of band names;

[0017] Each record in SensorList contains the following fields: SensorDate: the observation date number calculated from the sowing date of the current season; SensorID: the sensor node number; SensorLon: the node longitude; SensorLat: the node latitude; SensorSM: soil moisture; SensorEC: soil electrical conductivity; SensorN: nitrogen metabolism index;

[0018] Each record in PhenocalList contains the following fields: PhenStage: Phenological stage number; PhenDateStart: Phenological stage start date; PhenDateEnd: Phenological stage end date; PhenAnchorDate: Center date of the fertilization window for the current phenological stage; PhenFertWindow: Whether the current phenological stage contains a fertilization window (PhenFertWindow=1 indicates the current phenological stage contains a fertilization window, otherwise PhenFertWindow=0 indicates the current phenological stage does not contain a fertilization window); PhenWeight: Phenological weight for the current phenological stage.

[0019] Each record in FertHistList contains the following fields: FertDate: Fertilization date number; FertAmount: Fertilization amount; FertType: Fertilization type number;

[0020] The global control parameters include: FieldHeight (number of cell rows in the field height direction), FieldWidth (number of cell columns in the field width direction), PrescReso (prescription cell resolution), BaseAmount (basic fertilizer application amount), Sensor spatial search radius ZQMZRadius (sensor spatial search radius), Periodic anchoring time sequence Gaussian width ZQMZPeriodSigma (period time sequence width), Number of clusters JLBQClusterK (number of clusters), Minimum fertilizer application PrescMin (minimum fertilizer application amount), and Maximum fertilizer application PrescMax (maximum fertilizer application amount).

[0021] Furthermore, in S2, the air-to-ground joint generation input ZQMZInput is obtained using UAVList, PhenocalList, SensorList, and global control parameters. A periodically anchored air-to-ground timing joint generation module ZQMZModel is then established. The input of ZQMZModel is ZQMZInput, and the output is the air-to-ground joint generation output ZQMZOutput. Specifically:

[0022] S201. Establish a periodic anchored air-to-ground time series joint generation module ZQMZModel. The input of ZQMZModel is the air-to-ground joint generation input ZQMZInput.

[0023] The ZQMZInput includes: a UAVList record, the PhenAnchorDate of the phenological stage corresponding to the current UAVList record, a subset of SensorList within the time window of the current phenological stage, ZQMZRadius, and ZQMZPeriodSigma;

[0024] S202. Add the ZQMZFeatMap field to the air-ground joint generation output ZQMZOutput, and initialize the ZQMZFeatMap field of the air-ground joint generation output ZQMZOutput to a 6-dimensional array of all zeros with FieldHeight row and FieldWidth column.

[0025] S203. Initialize row traversal counter ZQMZCounter1=1;

[0026] S204. Initialize column traversal counter ZQMZCounter2=1;

[0027] S205. Extract the 5-dimensional pixel value of UAVImg from the ZQMZCounter1 row and ZQMZCounter2 column. Use the extracted 5-dimensional pixel value to obtain the normalized vegetation index temporary value ZQMZNdvi and the red-edge vegetation index temporary value ZQMZNdre. Specifically:

[0028] ZQMZNdvi=(ZQMZTemp2-ZQMZTemp1) / (ZQMZTemp2+ZQMZTemp1+0.001);

[0029] ZQMZNdre=(ZQMZTemp3-ZQMZTemp2) / (ZQMZTemp3+ZQMZTemp2+0.001);

[0030] Among them, ZQMZTemp1, ZQMZTemp2 and ZQMZTemp3 are all temporary values, which are the R-band, NIR-band and RE-band values ​​in the 5-dimensional pixel values, respectively. ZQMZNdvi is the temporary value of the normalized vegetation index, and ZQMZNdre is the temporary value of the red edge vegetation index.

[0031] S206. Based on the row and column numbers of the current pixel (ZQMZCounter1, ZQMZCounter2) and the geographical range of the field, calculate the true geographical latitude and longitude corresponding to the current pixel through linear mapping to obtain the temporary longitude ZQMZCurLon and the temporary latitude ZQMZCurLat.

[0032] S207. Initialize the weighted and cumulative soil moisture ZQMZWSumSM=0, the weighted and cumulative soil electrical conductivity ZQMZWSumEC=0, the weighted and cumulative nitrogen metabolism index ZQMZWSumN=0, and the total weighted and cumulative ZQMZWSum=0; initialize the sensor weight counter ZQMZCounter3=1.

[0033] S208. Use the ZQMZCounter3 record of SensorList, ZQMZCurLon, and ZQMZCurLat to obtain the distance ZQMZDist between the sensor and the current pixel. If ZQMZDist is greater than ZQMZRadius, go to S210; otherwise, obtain the temporal Gaussian weight ZQMZTempW, update ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZTempW and ZQMZDist, and then execute S209.

[0034] S209. Set ZQMZCounter3 = ZQMZCounter3 + 1; if ZQMZCounter3 is less than or equal to SensorTotal, go to S208; otherwise, go to S210.

[0035] Where SensorTotal is the total number of observation records in SensorList;

[0036] S210. If ZQMZWSum is less than 0.001, then set the current cell soil moisture value after spatiotemporal weighted interpolation to ZQMZEstSM=0, the current cell soil electrical conductivity value after spatiotemporal weighted interpolation to ZQMZEstEC=0, and the current cell nitrogen metabolism index value after spatiotemporal weighted interpolation to ZQMZEstN=0; otherwise, let ZQMZEstSM=ZQMZWSumSM / ZQMZWSum; ZQMZEstEC=ZQMZWSumEC / ZQMZWSum; let ZQMZEstN=ZQMZWSumN / ZQMZWSum.

[0037] S211, Calculate the sensor-enhanced nitrogen metabolism production: ZQMZEstNAdj=ZQMZEstN×(0.7+0.3×Min(ZQMZEstEC / 2.0,1.0));

[0038] Wherein, Min represents the smaller of the two values;

[0039] S212. Assign values ​​to each dimension of the ZQMZFeatMap in the ZQMZCounter1 row and ZQMZCounter2 column: 1st dimension = ZQMZNdvi; 2nd dimension = ZQMZNdre; 3rd dimension = ZQMZTempW; 4th dimension = ZQMZEstSM; 5th dimension = ZQMZEstEC; 6th dimension = ZQMZEstNAdj;

[0040] Where ZQMZTempW is the temporal Gaussian weight;

[0041] S213. Set ZQMZCounter2 = ZQMZCounter2 + 1; if ZQMZCounter2 is less than or equal to FieldWidth, go to S205; otherwise, go to S214.

[0042] S214. Set ZQMZCounter1 = ZQMZCounter1 + 1; if ZQMZCounter1 is less than or equal to FieldHeight, go to S204; otherwise, go to S215.

[0043] S215. Output ZQMZOutput as the result of ZQMZModel.

[0044] Furthermore, the distance ZQMZDist between the sensor and the current pixel in S208 is specifically as follows:

[0045] ZQMZDist = ABS(ZQMZCurLon - SensorLon of the 3rd record in SensorList) × FieldWidth + ABS(ZQMZCurLat - SensorLat of the 3rd record in SensorList) × FieldHeight;

[0046] Where ABS is the absolute value;

[0047] The specific steps for updating ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZDist are as follows:

[0048] ZQMZWSumSM = ZQMZWSumSM + ZQMZW × SensorSM of the 3rd ZQMZCounter item in SensorList;

[0049] ZQMZWSumEC = ZQMZWSumEC + ZQMZW × SensorEC of the 3rd ZQMZCounter item in SensorList;

[0050] ZQMZWSumN = ZQMZWSumN + ZQMZW × SensorN of the 3rd ZQMZCounter item in SensorList;

[0051] ZQMZWSum=ZQMZWSum+ZQMZW;

[0052] ZQMZW=ZQMZTempW×ZQMZDistW;

[0053] ZQMZDistW=1.0 / (ZQMZDist+1.5);

[0054] ZQMZTempW = exp(-((SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput) × (SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput)) / (2.0 × ZQMZPeriodSigma × ZQMZPeriodSigma));

[0055] Where exp is an exponential function with base e, and ZQMZW is the combined weight.

[0056] Furthermore, in step S3, the phenological accumulation input WHSJInput is obtained using UAVList, SensorList, PhenocalList, FertHistList, and global control parameters. A phenological history accumulation vegetation state vectorization module WHSJModel is then constructed. The input to WHSJModel is WHSJInput, and WHSJModel calls ZQMZModel to obtain the phenological accumulation output WHSJOutput. Specifically:

[0057] S301. Obtain the phenological accumulation input WHSJInput, and construct the phenological history accumulation vegetation state vectorization module WHSJModel. The input of WHSJModel is the phenological accumulation input WHSJInput.

[0058] The phenological accumulation input WHSJInput includes: all records in UAVList, all records in SensorList, all records in PhenocalList, all records in FertHistList, field height (number of rows), field width (number of columns), sensor spatial search radius (ZQMZRadius), and periodic anchoring time sequence Gaussian width (ZQMZPeriodSigma).

[0059] S302. Use the pixel-by-pixel state vector image WHSJStateMap and the last phase joint fusion feature image WHSJLastFeat to form the phenological accumulation output WHSJOutput, and initialize WHSJOutput; then establish the phenological weighted cumulative total image WHSJWtotalMap=FieldHeight row FieldWidth column all zero two-dimensional array;

[0060] Initialize the phenological accumulation output WHSJOutput as a WHSJStateMap = an 8-dimensional array of all zeros with FieldHeight rows and FieldWidth columns;

[0061] Initialize the phenological accumulation output WHSJOutput into a 6-dimensional array of all zeros, consisting of the last phase joint fusion feature image WHSJLastFeat = FieldHeight rows and FieldWidth columns;

[0062] S303. Establish a phenological stage traversal counter WHSJCounter1=1;

[0063] S304. Retrieve the WHSJCounter1 record from PhenocalList and store it in temporary storage WHSJTemp1; retrieve the record from UAVList whose UAVDate is closest to the PhenAnchorDate of WHSJTemp1 and store it in WHSJTemp2; retrieve all records from SensorList whose SensorDate is within the range of PhenDateStart to PhenDateEnd of WHSJTemp1 and store them in the sensor subset WHSJSensorSub.

[0064] S305. Fill the parameters of UAVImg of WHSJTemp2, PhenAnchorDate of WHSJTemp1, WHSJSensorSub, ZQMZRadius and ZQMZPeriodSigma into ZQMZInput, call ZQMZModel, obtain the fusion feature image ZQMZOutput of the current stage and store it in the temporary storage WHSJTemp3.

[0065] S306. Initialize the cell row counter WHSJCounter2=1; Initialize the cell column counter WHSJCounter3=1;

[0066] S307. Extract the 6-dimensional values ​​of ZQMZFeatMap in the WHSJCounter2 row and WHSJCounter3 column of WHSJTemp3, and denote the 1st dimension data, 2nd dimension data, and 4th to 6th dimension data as WHSJCurNdvi, WHSJCurNdre, WHSJCurSM, WHSJCurEC and WHSJCurN respectively.

[0067] S308. Use PhenWeight of WHSJTemp1 as the phenological weight WHSJPhenW for the current stage, update each dimension of WHSJStateMap with WHSJCounter2 and WHSJCounter3, and simultaneously assign WHSJPhenW to the corresponding row and column position of WHSJWtotalMap.

[0068] S309. If WHSJCounter1 equals PhenTotal, then set the value of WHSJLastFeat in the WHSJCounter2 row and WHSJCounter3 column to the value of ZQMZFeatMap of WHSJTemp3 in the WHSJCounter2 row and WHSJCounter3 column, and then execute S310; otherwise, execute S310 directly.

[0069] Wherein, PhenTotal is the number of phenological stages in PhenocalList;

[0070] S310. Set WHSJCounter3 = WHSJCounter3 + 1. If WHSJCounter3 is less than or equal to FieldWidth, go to S307. Otherwise, set WHSJCounter2 = WHSJCounter2 + 1. If WHSJCounter2 is less than or equal to FieldHeight, go to S307. Otherwise, go to S311.

[0071] S311. Set WHSJCounter1 = WHSJCounter1 + 1; if WHSJCounter1 is less than or equal to PhenTotal, go to S304; otherwise, go to S312.

[0072] S312. Normalize the first to fifth dimensions of WHSJStateMap according to the corresponding position values ​​of WHSJWtotalMap:

[0073] Traverse each pixel location and read the cumulative total of phenological weights WHSJWtotalMap corresponding to each location; if the WHSJWtotalMap of the current location is less than 0.001, then set the values ​​of the first to fifth dimensions of the current pixel's WHSJStateMap to 0; otherwise, divide the values ​​of the first to fifth dimensions of the current pixel's WHSJStateMap by the total weight value of the WHSJWtotalMap at the same location to obtain the normalized feature value after weighted average.

[0074] S313. Establish a fertilizer application history response counter WHSJCounter4=1; establish a fertilizer application attenuation cumulative amount WHSJFertAccum=0;

[0075] S314. Retrieve the 4th record of FertHistList (WHSJCounter4), take the last phenological stage PhenDateEnd as the reference date WHSJRefDate, and calculate the number of days since today: WHSJDeltaDays = WHSJRefDate - FertDate of the 4th record of FertHistList (WHSJCounter4). If WHSJDeltaDays is less than 0, set WHSJDeltaDays = 0; otherwise, keep the original value unchanged. Let WHSJFertAccum = WHSJFertAccum + FertAmount × exp(-WHSJDeltaDays / 30.0) of the 4th record of FertHistList (WHSJCounter4). Let WHSJCounter4 = WHSJCounter4 + 1. If WHSJCounter4 is less than or equal to FertTotal, go to S314; otherwise, go to S315.

[0076] Where FertTotal is the total number of fertilization events in FertHistList;

[0077] S315. Assign the 6th dimension of all pixels in WHSJStateMap to WHSJFertAccum / 200.0; if the 6th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 7th dimension of all pixels in WHSJStateMap to UAVTotal / Max(PhenTotal,1); if the 7th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 8th dimension of all pixels in WHSJStateMap to tanh(WHSJFertAccum / 100.0), and then update the 1st to 5th dimensions of WHSJStateMap in WHSJOutput: WHSJStateMap1st to 5th dimensions = arctan(WHSJStateMap1st to 5th dimensions × 2.0) / (π / 2); output WHSJOutput as the result of WHSJModel;

[0078] Where Max is the larger of the two values, tanh is the hyperbolic tangent function, arctan is the arctangent function, π is pi, and UAVTotal is the number of data periods in UAVList.

[0079] Furthermore, in S308, WHSJPhenW is used to update the data of each dimension of WHSJStateMap in the WHSJCounter2 row and WHSJCounter3 column, specifically as follows:

[0080] Updated first-dimensional data = original first-dimensional data + WHSJPhenW × WHSJCurNdvi;

[0081] Updated second-dimensional data = Unupdated second-dimensional data + WHSJPhenW × WHSJCurNdre;

[0082] Updated 3rd dimension data = Unupdated 3rd dimension data + WHSJPhenW × WHSJCurSM;

[0083] Updated fourth dimension data = original fourth dimension data + WHSJPhenW × WHSJCurEC;

[0084] The updated 5th dimension data = the original 5th dimension data + WHSJPhenW × WHSJCurN.

[0085] Further, in step S4, the prescription image generation input JLBQInput is obtained using WHSJOutput and global control parameters, and a clustered label weighted prescription image generation module JLBQModel is constructed. The input of JLBQModel is JLBQInput, and JLBQModel calls WHSJModel to obtain the prescription image generation output JLBQOutput, specifically as follows:

[0086] S401. Construct a clustered label weighted prescription image generation module JLBQModel, where the input of JLBQModel is JLBQInput;

[0087] The prescription image generation input JLBQInput includes: the entire contents of WHSJOutput, BaseAmount, JLBQClusterK, PrescMin, and PrescMax;

[0088] S402. Using UAVList, SensorList, PhenocalList, FertHistList and the global constructor WHSJInput, call WHSJModel to obtain the temporary storage JLBQTemp1 of WHSJOutput;

[0089] S403. Expand the WHSJStateMap of JLBQTemp1 into a two-dimensional feature matrix JLBQFeatMat by rows and columns;

[0090] S404. Perform K-means clustering on JLBQFeatMat based on JLBQClusterK to obtain the cluster label image JLBQLabelMap and the cluster center matrix JLBQCenters for each pixel.

[0091] S405. Calculate the mean of the first dimension of data for each cluster center, sort the cluster numbers according to the mean of the first dimension of data from smallest to largest, and obtain the list of cluster numbers in ascending order JLBQSortedClust.

[0092] S406. Create a fertilizer weight list JLBQWeightList; create a weight assignment counter JLBQCounter1=1; set the minimum fertilizer weight JLBQWMin=0.60 and the total fertilizer weight span JLBQWRange=0.70;

[0093] S407, the cluster weight of the JLBQSortedClust[JLBQCounter1] in JLBQWeightList = JLBQWMin + JLBQWRange × (JLBQCounter1-1) / (JLBQClusterK-1);

[0094] Set JLBQCounter1 = JLBQCounter1 + 1; if JLBQCounter1 is less than or equal to JLBQClusterK, go to S407, otherwise go to S408;

[0095] Where JLBQSortedClust[JLBQCounter1] is the value of the JLBQCounter1 element in JLBQSortedClust;

[0096] S408. Calculate the first-dimensional mean value of the entire field JLBQTemp1, JLBQNDVIMean = Sum of the first-dimensional values ​​of all cells / (FieldHeight × FieldWidth);

[0097] S409. Initialize the prescription image and generate the pixel fertilizer amount image JLBQOutput JLBQPrescImg field = FieldHeight row FieldWidth column single-band all-zero floating-point raster image; establish the prescription row counter JLBQCounter2=1;

[0098] S410. Set the prescription column counter JLBQCounter3=1;

[0099] S411. Using the WHSJLastFeat data of JLBQTemp1 in the first dimension of the JLBQCounter2 row and JLBQCounter3 column and the value of JLBQLabelMap in the JLBQCounter2 row and JLBQCounter3 column, obtain the prescription image value JLBQPrescrip of the current pixel, and assign the pixel value of JLBQPrescImg in the JLBQCounter2 row and JLBQCounter3 column to JLBQPrescrip;

[0100] S412, Set JLBQCounter3 = JLBQCounter3 + 1; If JLBQCounter3 is less than or equal to FieldWidth, go to S411; Otherwise, Set JLBQCounter2 = JLBQCounter2 + 1; If JLBQCounter2 is less than or equal to FieldHeight, go to S410; Otherwise, go to S413.

[0101] S413. Output JLBQOutput as the result of JLBQModel.

[0102] Further, in S411, the prescription image value JLBQPrescrip of the current pixel is obtained by using the WHSJLastFeat data of JLBQTemp1 in the first dimension of the JLBQCounter2 row and the JLBQCounter3 column, and the value of JLBQLabelMap in the JLBQCounter2 row and the JLBQCounter3 column.

[0103] Obtain the WHSJLastFeat of the latest NDVI value JLBQCurNdvi=JLBQTemp1 of the current pixel in the first dimension of the JLBQCounter2 row and JLBQCounter3 column; calculate the base amount spatial modulation JLBQBase=BaseAmount×(1.0-0.25×JLBQCurNdvi / JLBQNDVIMean); if JLBQBase is less than PrescMin×0.5, then set JLBQBase=PrescMin×0.5, otherwise keep JLBQBase unchanged;

[0104] Get the cluster label JLBQCurLabel=JLBQLabelMap of the current cell in the JLBQCounter2 row and JLBQCounter3 column; calculate the prescription image value of the current cell JLBQPrescrip=JLBQBase×JLBQWeightList[JLBQCurLabel];

[0105] Where JLBQWeightList[JLBQCurLabel] is the element value of JLBQCurLabel in JLBQWeightList.

[0106] If JLBQPrescrip is less than PrescMin, then set JLBQPrescrip=PrescMin; if JLBQPrescrip is greater than PrescMax, then set JLBQPrescrip=PrescMax; otherwise, keep JLBQPrescrip unchanged.

[0107] Furthermore, the step of calling the phenological history accumulation vegetation state vectorization module JLBQModel to generate a fertilization prescription image, statistically analyzing the results of each partition of the fertilization prescription image, and outputting the final fertilization prescription image specifically involves:

[0108] A1. Call JLBQModel, input JLBQInput into JLBQModel, and obtain JLBQOutput; assign the pixel fertilizer application image JLBQPrescImg of JLBQOutput to the final fertilizer prescription image PrescImg;

[0109] A2. Create a cumulative fertilizer application list for each cluster region, PrescStatList = JLBQCluster, containing K zero-element lists; create a pixel count list for each cluster, PrescPixList = JLBQCluster, containing K zero-element lists.

[0110] A3. Initialize the output row traversal counter OutCounter1=1; Initialize the output column traversal counter OutCounter2=1;

[0111] A4. Take the cell value of PrescImg in the OutCounter1 row and OutCounter2 column and store it in the temporary storage OutTemp1; take the cluster label OutTemp2 at the corresponding position of JLBQLabelMap;

[0112] Let PrescStatList[OutTemp2] = PrescStatList[OutTemp2] + OutTemp1;

[0113] Set PrescPixList[OutTemp2] = PrescPixList[OutTemp2] + 1;

[0114] Where PrescStatList[OutTemp2] is the element value corresponding to OutTemp2 in PrescStatList;

[0115] A5. Set OutCounter2 = OutCounter2 + 1; if OutCounter2 is less than or equal to FieldWidth, go to A4; otherwise, set OutCounter1 = OutCounter1 + 1; if OutCounter1 is less than or equal to FieldHeight, go to A3; otherwise, go to A6.

[0116] A6. Set the partitioned statistics output counter StatCounter=0;

[0117] A7. Calculate the average fertilizer application for the StatCounter-th cluster: PrescStatAvg = PrescStatList[StatCounter] / Max(PrescPixList[StatCounter], 1); Output the number of pixels in the StatCounter-th cluster: PrescPixList[StatCounter] and the average fertilizer application: PrescStatAvg; Set StatCounter = StatCounter + 1; If StatCounter is less than or equal to JLBQClusterK-1, go to A7; otherwise, go to A8.

[0118] Where PrescStatList[StatCounter] is the value of the StatCounter-th element in PrescStatList, and Max(PrescPixList[StatCounter],1) is to get the larger value between PrescPixList[StatCounter] and 1;

[0119] A8. Output PrescImg as the final fertilizer prescription image.

[0120] The beneficial effects of this invention are as follows:

[0121] This invention provides a method for generating variable fertilization prescription images that combines periodicity and phenology. It constructs a periodic anchored air-ground time-series joint generation module and a phenological history-accumulated vegetation state vectorization module. Based on these, a clustering label-weighted prescription image generation module is established. The periodic anchored air-ground time-series joint generation module uses the fertilization phenological window as the time anchor point to perform Gaussian time-series weighted spatial joint processing on UAV multispectral data and ground wireless sensor data, eliminating spatiotemporal alignment bias and forming a pixel-by-pixel fused feature image centered on the key fertilization period. The phenological history-accumulated vegetation state vectorization module... The module iteratively calls the phenological stage weighted vegetation index and fertilization history decay response to accumulate each stage, forming an 8-dimensional state vector for each pixel location. Then, the cluster label weighted prescription image generation module performs a fixed number of unsupervised clusters on the state vectors of all field pixels and labels them in ascending order of cluster mean. Combined with the spatial basic fertilizer amount, the module calculates and generates a floating-point fertilizer prescription image pixel by pixel. The prescription image is a single-band raster image, and the value of each pixel is a continuous floating-point fertilizer amount at the corresponding spatial location. This realizes the pixel-by-pixel programmatic generation of multimodal variable fertilizer prescription images for open fields during the crop growth period. The periodic anchored air-ground time-series joint generation module and the phenological history accumulation vegetation state vectorization module of this invention do not rely on single-period static data and manual zoning rules. Instead, they use the fertilization phenological window period as the time anchor point, utilize Gaussian decay weights to prioritize the adoption of UAV and sensor data near the critical period, and accumulate vegetation state vectors with weighted differences in phenological stages across periods. This is unaffected by light disturbances at the moment of a single collection and incorporates the 30-day decay response of fertilization history into the state vector. It integrates operational phenological history accumulation information and fertilization history information, while also considering the key influences of soil moisture, electrical conductivity, and nitrogen metabolism indicators on fertilization effects, as well as the needs of different operations and different growth stages. This improves the ability to perceive abnormal areas of crop production in the field, avoids the problem of local over-fertilization or under-fertilization, and thus improves the accuracy of fertilization prescription map decision-making. On the other hand, the clustering label weighted prescription image generation module does not use manual thresholding to delineate management areas. Instead, it performs a fixed number of K-means clustering on the 8-dimensional state vectors of all field pixels, assigns fertilization weight coefficients in ascending order of the mean of each cluster's state vector, and calculates and generates a floating-point prescription image that takes into account spatial heterogeneity and growth differences on a pixel-by-pixel basis using the real-time NDVI modulation of the phenological fusion image. Each pixel in the generated prescription image has a continuous floating-point fertilization value, which can be directly used as a pixel-by-pixel variable fertilization instruction image for drones or precision fertilization equipment, improving the ability to adapt to differences and the accuracy of the fertilization prescription map, thereby achieving precision fertilization. Attached Figure Description

[0122] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0123] Specific implementation method one: as follows Figure 1 As shown, the specific process of the variable fertilization prescription image generation method combining periodicity and phenology in this embodiment is as follows:

[0124] S1. Input the UAV multispectral data time series table UAVList, the ground wireless sensor data table SensorList, the crop phenological history table PhenocalList, and the fertilization history table FertHistList, as well as the global control parameters. Obtain the data period number UAVTotal from UAVList, the total number of observation records SensorTotal from SensorList, the number of phenological stages PhenTotal from PhenocalList, and the total number of fertilization events FertTotal from FertHistList. Then, create the fertilization prescription image PrescImg, specifically:

[0125] S101. Input the UAV multispectral data time series table UAVList. Each record in UAVList corresponds to one period of UAV flight data acquisition results and includes the following fields:

[0126] UAVDate: Data collection date number (days, calculated from the sowing date of the current season);

[0127] UAVImg: Multispectral data, with dimensions of FieldHeight rows × FieldWidth columns × UAVBandNum bands;

[0128] Where FieldHeight is the number of pixel rows in the field height direction, FieldWidth is the number of pixel columns in the field width direction, and UAVBandNum is the total number of multispectral bands of the UAV.

[0129] UAVBandList: A list of band names, containing five bands in sequence: R (red), G (green), B (blue), RE (red edge), and NIR (near-infrared).

[0130] In this invention, the date numbers refer to the day following the sowing date of the current season;

[0131] S102. Input the ground wireless sensor data table SensorList. Each record in SensorList corresponds to one observation of a sensor node at a certain time and contains the following fields:

[0132] SensorDate: Observation date number (days);

[0133] SensorID: Sensor node number;

[0134] SensorLon: Node longitude;

[0135] SensorLat: Node latitude;

[0136] SensorSM: Soil Moisture (percentage);

[0137] SensorEC: Soil electrical conductivity (milliseconds per centimeter);

[0138] SensorN: Nitrogen metabolism index (relative units, 0 to 1);

[0139] In this step, nitrogen metabolism indicators are obtained by normalizing the data acquired by the soil ammonium nitrogen sensor; node longitude and node latitude are in decimal form;

[0140] S103. Input the crop phenological history table PhenocalList. Each record in PhenocalList corresponds to one phenological stage within the current growing season and contains the following fields:

[0141] PhenStage: Phenological stage number (1 to 5, in order: seedling emergence-seedling establishment stage, vegetative growth stage, vigorous growth-construction stage, reproductive initiation stage, and grain filling-prematurity stage);

[0142] PhenDateStart: Date number (in days) at the start of the phenological phase;

[0143] PhenDateEnd: Date number (in days) at the end of the phenological phase;

[0144] PhenAnchorDate: The center date (in days) of the fertilization window for the current phenological stage. If there is no fertilization window for this stage, it is the same as PhenDateStart.

[0145] PhenFertWindow: Whether the current phenological stage contains a fertilization window (1 for yes, 0 for no);

[0146] PhenWeight: The phenological weight of the current phenological stage (0 to 1), with a higher weight for important fertilization stages;

[0147] In this step, when the phenological stage is the seedling emergence-seedling establishment period, the phenological weight of the current phenological stage is 0.5; when the phenological stage is the vegetative growth period, the phenological weight of the current phenological stage is 0.88.

[0148] When the phenological stage is the vigorous growth-construction period, the phenological weight of the current phenological stage is 1.00;

[0149] When the phenological stage is the reproductive initiation period, the phenological weight of the current phenological stage is 0.82;

[0150] When the phenological stage is from grain filling to early maturity, the phenological weight of the current phenological stage is 0.48.

[0151] S104. Input the fertilization history table FertHistList. Each record in FertHistList corresponds to one historical fertilization event and includes the following fields:

[0152] FertDate: Fertilization date (days);

[0153] FertAmount: Fertilizer application rate (kg per hectare);

[0154] FertType: Fertilizer type number;

[0155] S105. Input global control parameters, specifically:

[0156] Input the number of cell rows in the field height direction, FieldHeight. The default value of FieldHeight is 200.

[0157] Enter the number of cell columns in the field width direction, FieldWidth. The default value for FieldWidth is 200.

[0158] Enter the prescription pixel resolution PrescReso (meters). The default value for PrescReso is 10.

[0159] Enter the base fertilizer application amount BaseAmount (kilograms per hectare). The default value for BaseAmount is 120.

[0160] Input the sensor spatial search radius ZQMZRadius (number of pixels), the default value of ZQMZRadius is 5;

[0161] Input the Gaussian width of the period anchoring sequence ZQMZPeriodSigma (days), the default value of ZQMZPeriodSigma is 7;

[0162] Input the number of clusters, JLBQClusterK. The default value for JLBQClusterK is 5.

[0163] Enter the minimum fertilizer application amount PrescMin (kg per hectare). The default value for PrescMin is 60.

[0164] Enter the maximum amount of fertilizer to apply (PrescMax, kilograms per hectare). The default value for PrescMax is 200.

[0165] S106. Get the number of data periods in UAVList (UAVTotal); get the total number of observation records in SensorList (SensorTotal); get the number of phenological stages in PhenocalList (PhenTotal); get the total number of fertilization events in FertHistList (FertTotal);

[0166] S107. Create a single-band floating-point raster image with PrescImg = FieldHeight row and FieldWidth column for fertilizer prescription image;

[0167] In the fertilizer prescription image PrescImg, the initial value of each pixel is 0, and each pixel represents the amount of fertilizer applied at the corresponding spatial location.

[0168] S2. Using UAVList, PhenocalList, SensorList, and global control parameters, obtain the air-to-ground joint generation input ZQMZInput, and establish a periodically anchored air-to-ground timing joint generation module ZQMZModel. The input of ZQMZModel is the air-to-ground joint generation input ZQMZInput, and the output of ZQMZModel is the air-to-ground joint generation output ZQMZOutput. Specifically:

[0169] S201. Use UAVList, PhenocalList, SensorList and global control parameters to obtain the air-to-ground joint generation input ZQMZInput, and establish a periodically anchored air-to-ground timing joint generation module ZQMZModel. The input of ZQMZModel is the air-to-ground joint generation input ZQMZInput.

[0170] The air-ground joint generation input ZQMZInput includes: a UAVList record, the PhenAnchorDate of the phenological stage corresponding to the current UAVList record, a subset of SensorList within the time window of the current phenological stage, and ZQMZRadius and ZQMZPeriodSigma parameters;

[0171] S202. Add a ZQMZFeatMap field to the air-ground joint generation output ZQMZOutput. Initialize the ZQMZFeatMap field of the air-ground joint generation output ZQMZOutput to a 6-dimensional array of all zeros with FieldHeight rows and FieldWidth columns (i.e., an image that stores 6-dimensional feature vectors pixel by pixel).

[0172] Among them, ZQMZFeatMap is a fused feature image of FieldHeight×FieldWidth×6 dimensions, where each pixel position in the image stores a 6-dimensional floating-point vector;

[0173] S203. Initialize row traversal counter ZQMZCounter1=1;

[0174] S204. Initialize column traversal counter ZQMZCounter2=1;

[0175] S205. Extract the 5-dimensional pixel value of UAVImg from the ZQMZCounter1 row and ZQMZCounter2 column. Use the extracted 5-dimensional pixel value to obtain the normalized vegetation index temporary value ZQMZNdvi and the red-edge vegetation index temporary value ZQMZNdre. Specifically:

[0176] First, the R-band, NIR-band, and RE-band values ​​extracted from the 5-dimensional pixel values ​​are stored in temporary storage values ​​ZQMZTemp1, ZQMZTemp2, and ZQMZTemp3, respectively.

[0177] Then, the normalized vegetation index transient ZQMZNdvi is calculated using ZQMZTemp1 and ZQMZTemp2, specifically as follows:

[0178] ZQMZNdvi=(ZQMZTemp2-ZQMZTemp1) / (ZQMZTemp2+ZQMZTemp1+0.001);

[0179] Then, the red-edge vegetation index transient ZQMZNdre is calculated using ZQMZTemp2 and ZQMZTemp3, specifically as follows:

[0180] ZQMZNdre=(ZQMZTemp3-ZQMZTemp2) / (ZQMZTemp3+ZQMZTemp2+0.001);

[0181] If ZQMZNdvi is greater than 1.0, then ZQMZNdvi = 1.0; if ZQMZNdvi is less than -1.0, then ZQMZNdvi = -1.0; if ZQMZNdvi is less than or equal to 1.0 and greater than -1.0, then the original value remains unchanged.

[0182] S206. Based on the row and column numbers of the current pixel (ZQMZCounter1, ZQMZCounter2) and the geographical range of the field, calculate the true geographical latitude and longitude corresponding to the current pixel through linear mapping to obtain the temporary longitude ZQMZCurLon and the temporary latitude ZQMZCurLat.

[0183] S207. Initialize the weighted and cumulative soil moisture ZQMZWSumSM=0, the weighted and cumulative soil electrical conductivity ZQMZWSumEC=0, the weighted and cumulative nitrogen metabolism index ZQMZWSumN=0, and the total weighted and cumulative ZQMZWSum=0; initialize the sensor weight counter ZQMZCounter3=1.

[0184] S208. Retrieve the ZQMZCounter3 record from SensorList. Use the ZQMZCounter3 record from SensorList, the longitude temporary storage ZQMZCurLon, and the latitude temporary storage to obtain the distance ZQMZDist between the sensor and the current pixel. If ZQMZDist is greater than ZQMZRadius, proceed to S210; otherwise, update ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZDist, and then execute S209.

[0185] The method of obtaining the distance ZQMZDist between the sensor and the current pixel using the ZQMZCounter3 record of SensorList, the longitude temporary storage ZQMZCurLon, and the latitude temporary storage is as follows:

[0186] ZQMZDist = ABS(ZQMZCurLon - SensorLon of the 3rd record in SensorList) × FieldWidth + ABS(ZQMZCurLat - SensorLat of the 3rd record in SensorList) × FieldHeight;

[0187] Where ABS represents the absolute value; the specific steps of updating ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZDist are as follows:

[0188] First, calculate the time-series Gaussian weights ZQMZTempW:

[0189] ZQMZTempW = exp(-((SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput) × (SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput)) / (2.0 × ZQMZPeriodSigma × ZQMZPeriodSigma));

[0190] Where exp is an exponential function with the natural constant e as its base;

[0191] Then, the distance decay weight ZQMZDistW is calculated using ZQMZDist, specifically as follows:

[0192] ZQMZDistW = 1.0 / (ZQMZDist + 1.5)

[0193] Finally, the comprehensive weight ZQMZW is obtained using ZQMZDistW and ZQMZTempW, and ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum are updated using ZQMZW, specifically as follows:

[0194] ZQMZW=ZQMZTempW×ZQMZDistW;

[0195] ZQMZWSumSM = ZQMZWSumSM + ZQMZW × SensorSM of the 3rd ZQMZCounter item in SensorList;

[0196] ZQMZWSumEC = ZQMZWSumEC + ZQMZW × SensorEC of the 3rd ZQMZCounter item in SensorList;

[0197] ZQMZWSumN = ZQMZWSumN + ZQMZW × SensorN of the 3rd ZQMZCounter item in SensorList;

[0198] ZQMZWSum=ZQMZWSum+ZQMZW;

[0199] S209. Set ZQMZCounter3 = ZQMZCounter3 + 1; if ZQMZCounter3 is less than or equal to SensorTotal, go to S208; otherwise, go to S210.

[0200] S210. If ZQMZWSum is less than 0.001, then set the current cell soil moisture value after spatiotemporal weighted interpolation to ZQMZEstSM=0, the current cell soil electrical conductivity value after spatiotemporal weighted interpolation to ZQMZEstEC=0, and the current cell nitrogen metabolism index value after spatiotemporal weighted interpolation to ZQMZEstN=0; otherwise, let ZQMZEstSM=ZQMZWSumSM / ZQMZWSum; let ZQMZEstEC=ZQMZWSumEC / ZQMZWSum; let ZQMZEstN=ZQMZWSumN / ZQMZWSum.

[0201] S211, Calculate the sensor-enhanced nitrogen metabolism production: ZQMZEstNAdj=ZQMZEstN×(0.7+0.3×Min(ZQMZEstEC / 2.0,1.0));

[0202] Wherein, Min represents the smaller of the two values;

[0203] S212. Assign values ​​to each dimension of the 6-dimensional vector in the ZQMZFeatMap at the ZQMZCounter1 row and ZQMZCounter2 column: 1st dimension = ZQMZNdvi; 2nd dimension = ZQMZNdre; 3rd dimension = ZQMZTempW (take the temporal weight of the current sensor, set to 0 if ZQMZWSum is less than 0.001); 4th dimension = ZQMZEstSM; 5th dimension = ZQMZEstEC; 6th dimension = ZQMZEstNAdj;

[0204] S213. Set ZQMZCounter2 = ZQMZCounter2 + 1; if ZQMZCounter2 is less than or equal to FieldWidth, go to S205; otherwise, go to S214.

[0205] S214. Set ZQMZCounter1 = ZQMZCounter1 + 1; if ZQMZCounter1 is less than or equal to FieldHeight, go to S204; otherwise, go to S215.

[0206] S215. Output ZQMZOutput as the result of ZQMZModel;

[0207] In this step, the periodic anchored air-to-ground time-series joint generation module ZQMZModel uses the fertilization phenology anchor date as a reference to perform Gaussian time-series weighted spatial joint processing on UAV data and ground sensor data, generating a fused feature image pixel by pixel to eliminate the spatiotemporal alignment deviation of air-to-ground data.

[0208] S3. Using UAVList, SensorList, PhenocalList, FertHistList, and global control parameters, obtain the phenological accumulation input WHSJInput, and construct the phenological history accumulation vegetation state vectorization module WHSJModel. The input of WHSJModel is the phenological accumulation input WHSJInput. WHSJModel calls ZQMZModel to obtain the phenological accumulation output WHSJOutput, specifically:

[0209] S301. Use UAVList, SensorList, PhenocalList, FertHistList and global control parameters to obtain the phenological accumulation input WHSJInput, and construct the phenological history accumulation vegetation state vectorization module WHSJModel. The input of WHSJModel is the phenological accumulation input WHSJInput.

[0210] The phenological accumulation input WHSJInput includes: all records in UAVList, all records in SensorList, all records in PhenocalList, all records in FertHistList, field height (number of rows), field width (number of columns), sensor spatial search radius (ZQMZRadius), and periodic anchoring time sequence Gaussian width (ZQMZPeriodSigma).

[0211] S302. Initialize the phenological accumulation output WHSJOutput, which includes: a pixel-by-pixel state vector image WHSJStateMap and a last-period joint fusion feature image WHSJLastFeat; initialize the phenological accumulation output WHSJOutput as an 8-dimensional array of all zeros (i.e., an image storing 8-dimensional state vectors pixel by pixel); initialize the last-period joint fusion feature image WHSJOutput as a 6-dimensional array of all zeros (WHSJLastFeat = FieldHeight = FieldWidth = FieldWidth); then establish the phenological weight cumulative total image WHSJWtotalMap = a 2-dimensional array of all zeros (WHSJWtotalMap = FieldHeight = FieldWidth = FieldWidth).

[0212] S303. Establish a phenological stage traversal counter WHSJCounter1=1;

[0213] S304. Retrieve the WHSJCounter1 record from PhenocalList and store it in temporary storage WHSJTemp1; retrieve the record from UAVList whose UAVDate is closest to the PhenAnchorDate of WHSJTemp1 and store it in WHSJTemp2; retrieve all records from SensorList whose SensorDate is within the range of PhenDateStart to PhenDateEnd of WHSJTemp1 and store them in the sensor subset WHSJSensorSub.

[0214] S305. Fill the parameters of UAVImg of WHSJTemp2, PhenAnchorDate of WHSJTemp1, WHSJSensorSub, ZQMZRadius and ZQMZPeriodSigma into ZQMZInput, call ZQMZModel, obtain the fusion feature image ZQMZOutput of the current stage and store it in the temporary storage WHSJTemp3.

[0215] S306. Initialize the cell row counter WHSJCounter2=1; Initialize the cell column counter WHSJCounter3=1;

[0216] S307. Extract the 6-dimensional values ​​of ZQMZFeatMap in the WHSJCounter2 row and WHSJCounter3 column of WHSJTemp3. Record the data of the first dimension (NDVI) as WHSJCurNdvi, the data of the second dimension (NDRE) as WHSJCurNdre, and the data of the fourth to sixth dimensions (SM, EC, N) as WHSJCurSM, WHSJCurEC, and WHSJCurN respectively.

[0217] S308. Use PhenWeight of WHSJTemp1 as the phenological weight WHSJPhenW for the current stage. Update each dimension of WHSJStateMap in the WHSJCounter2 row and WHSJCounter3 column using WHSJPhenW. Specifically: Updated 1st dimension = Updated 1st dimension + WHSJPhenW × WHSJCurNdvi; Updated 2nd dimension = Updated 2nd dimension + WHSJPhenW × WHSJCurNdre; Updated 3rd dimension = Updated 3rd dimension + WHSJPhenW × WHSJCurSM; Updated 4th dimension = Updated 4th dimension + WHSJPhenW × WHSJCurEC; Updated 5th dimension = Updated 5th dimension + WHSJPhenW × WHSJCurN. Simultaneously, assign WHSJPhenW to the corresponding row and column positions of WHSJWtotalMap.

[0218] S309. If WHSJCounter1 equals PhenTotal, then set the value of WHSJLastFeat in the WHSJCounter2 row and WHSJCounter3 column to the value of ZQMZFeatMap of WHSJTemp3 in the WHSJCounter2 row and WHSJCounter3 column, and then execute S310; otherwise, execute S310 directly.

[0219] S310. Set WHSJCounter3 = WHSJCounter3 + 1. If WHSJCounter3 is less than or equal to FieldWidth, go to S307. Otherwise, set WHSJCounter2 = WHSJCounter2 + 1. If WHSJCounter2 is less than or equal to FieldHeight, go to S307. Otherwise, go to S311.

[0220] S311. Set WHSJCounter1 = WHSJCounter1 + 1; if WHSJCounter1 is less than or equal to PhenTotal, go to S304; otherwise, go to S312.

[0221] S312. Normalize the first to fifth dimensions of WHSJStateMap according to the corresponding position values ​​of WHSJWtotalMap:

[0222] Iterate through each pixel location and read the cumulative total phenological weight WHSJWtotalMap value corresponding to each location;

[0223] If the current position WHSJWtotalMap is less than 0.001, to avoid division by zero, the values ​​of the first to fifth dimensions of the current cell WHSJStateMap are set to 0; if WHSJWtotalMap is greater than or equal to 0.001, the values ​​of the first to fifth dimensions of the current cell WHSJStateMap are divided by the total weight value of WHSJWtotalMap at the same position to obtain the normalized feature value after weighted average.

[0224] S313. Establish a fertilizer application history response counter WHSJCounter4=1; establish a fertilizer application attenuation cumulative amount WHSJFertAccum=0;

[0225] S314. Retrieve the 4th record of FertHistList (WHSJCounter4), take the last phenological stage PhenDateEnd as the reference date WHSJRefDate, and calculate the number of days since today: WHSJDeltaDays = WHSJRefDate - FertDate of the 4th record of FertHistList (WHSJCounter4). If WHSJDeltaDays is less than 0, set WHSJDeltaDays = 0; otherwise, keep the original value unchanged. Let WHSJFertAccum = WHSJFertAccum + FertAmount × exp(-WHSJDeltaDays / 30.0) of the 4th record of FertHistList (WHSJCounter4). Let WHSJCounter4 = WHSJCounter4 + 1. If WHSJCounter4 is less than or equal to FertTotal, go to S314; otherwise, go to S315.

[0226] S315. Assign the 6th dimension of all pixels in WHSJStateMap to WHSJFertAccum / 200.0; if the 6th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 7th dimension of all pixels in WHSJStateMap to UAVTotal / Max(PhenTotal,1); if the 7th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 8th dimension of all pixels in WHSJStateMap to tanh(WHSJFertAccum / 100.0), and then update the 1st to 5th dimensions of WHSJStateMap in WHSJOutput: WHSJStateMap1st to 5th dimensions = arctan(WHSJStateMap1st to 5th dimensions × 2.0) / (π / 2); output WHSJOutput as the result of WHSJModel;

[0227] Where Max is the larger of the two values, tanh is the hyperbolic tangent function, arctan is the arctangent function, and π is pi.

[0228] S4. Using WHSJOutput and global control parameters, obtain the prescription image generation input JLBQInput, and construct the cluster label weighted prescription image generation module JLBQModel. The input of JLBQModel is JLBQInput. JLBQModel calls the phenological history accumulation vegetation state vectorization module WHSJModel to obtain the prescription image generation output JLBQOutput. Specifically:

[0229] S401. Obtain the prescription image generation input JLBQInput using WHSJOutput and global control parameters, and construct the cluster label weighted prescription image generation module JLBQModel. The input of JLBQModel is JLBQInput.

[0230] The prescription image generation input JLBQInput includes: the entire contents of WHSJOutput, BaseAmount, JLBQClusterK, PrescMin, and PrescMax;

[0231] S402. Using UAVList, SensorList, PhenocalList, FertHistList and the global constructor WHSJInput, call WHSJModel to obtain the temporary storage JLBQTemp1 of WHSJOutput;

[0232] S403. Expand the WHSJStateMap of JLBQTemp1 into a two-dimensional feature matrix JLBQFeatMat with dimensions of (FieldHeight×FieldWidth) rows × 8 columns.

[0233] S404. Perform K-means clustering on JLBQFeatMat (K-Means algorithm, number of clusters JLBQClusterK, maximum number of iterations is 100) to obtain the cluster label image JLBQLabelMap for each pixel (dimension is FieldHeight rows × FieldWidth columns, and the value of each pixel is an integer from 0 to JLBQClusterK-1) and the cluster center matrix JLBQCenters (JLBQClusterK rows × 8 columns).

[0234] S405. Calculate the mean of the first dimension data (normalized phenological accumulation NDVI) of each cluster center, sort the cluster numbers according to the mean of the first dimension data from smallest to largest, and obtain the ascending list of cluster numbers JLBQSortedClust (the cluster with the smallest mean NDVI is ranked first, corresponding to the area with the worst crop growth, and labeled as the partition that needs more fertilizer).

[0235] S406. Create a fertilizer weight list JLBQWeightList; create a weight assignment counter JLBQCounter1=1; set the minimum fertilizer weight JLBQWMin=0.60 and the total fertilizer weight span JLBQWRange=0.70;

[0236] S407, the cluster weight of the JLBQSortedClust[JLBQCounter1] in JLBQWeightList = JLBQWMin + JLBQWRange × (JLBQCounter1-1) / (JLBQClusterK-1);

[0237] Set JLBQCounter1 = JLBQCounter1 + 1; if JLBQCounter1 is less than or equal to JLBQClusterK, go to S407, otherwise go to S408;

[0238] Where JLBQSortedClust[JLBQCounter1] is the value of the JLBQCounter1 element in JLBQSortedClust;

[0239] S408. Calculate the mean value of the entire field in the first dimension (last NDVI image) of WHSJLastFeat JLBQTemp1: JLBQNDVIMean = Sum of the first dimension values ​​of all pixels / (FieldHeight × FieldWidth); If JLBQNDVIMean is less than 0.01, set JLBQNDVIMean = 0.01; otherwise, keep the JLBQNDVIMean value unchanged.

[0240] S409. Initialize the prescription image and generate the pixel fertilizer amount image JLBQOutput JLBQPrescImg field = FieldHeight row FieldWidth column single-band all-zero floating-point raster image (each pixel represents the fertilizer amount at the corresponding spatial location); establish the prescription row counter JLBQCounter2=1;

[0241] S410. Set the prescription column counter JLBQCounter3=1;

[0242] S411. Take the latest NDVI value of the current pixel, JLBQCurNdvi=JLBQTemp1, and its WHSJLastFeat data in the first dimension of the JLBQCounter2 row and JLBQCounter3 column; calculate the base spatial modulation amount JLBQBase=BaseAmount×(1.0-0.25×JLBQCurNdvi / JLBQNDVIMean); if JLBQBase is less than PrescMin×0.5, set JLBQBase=PrescMin×0.5, otherwise keep JLBQBase unchanged; take the cluster label of the current pixel, JLBQCurLabel=JLBQLabelMap, in the first dimension of the JLBQCounter2 row and JLBQCounter3 column; The value in the 3rd column of JLBQCounter in row ter2; calculate the prescription image value of the current cell JLBQPrescrip = JLBQBase × JLBQWeightList[JLBQCurLabel]; if JLBQPrescrip is less than PrescMin, set JLBQPrescrip = PrescMin; if JLBQPrescrip is greater than PrescMax, set JLBQPrescrip = PrescMax; otherwise, keep JLBQPrescrip unchanged; assign the cell value of JLBQPrescImg in the 3rd column of JLBQCounter in row 2 to JLBQPrescrip;

[0243] Where JLBQWeightList[JLBQCurLabel] is the element value of JLBQCurLabel in JLBQWeightList.

[0244] S412, Set JLBQCounter3 = JLBQCounter3 + 1; If JLBQCounter3 is less than or equal to FieldWidth, go to S411; Otherwise, Set JLBQCounter2 = JLBQCounter2 + 1; If JLBQCounter2 is less than or equal to FieldHeight, go to S410; Otherwise, go to S413.

[0245] S413. Output JLBQOutput as the result of JLBQModel.

[0246] In this step, the phenological history accumulation vegetation state vectorization module WHSJModel iteratively calls ZQMZModel to generate air-ground joint data for each phenological stage, and accumulates vegetation states with phenological weight differences stage by stage. At the same time, the fertilization history decay response is introduced to form a multi-dimensional comprehensive state vector for each spatial location on a pixel-by-pixel basis.

[0247] S5. Call JLBQModel to generate a fertilizer prescription image, analyze the results of each partition of the fertilizer prescription image, and output the final fertilizer prescription image. Specifically:

[0248] S501. Construct JLBQInput, call JLBQModel, input JLBQInput into JLBQModel, and obtain JLBQOutput; assign JLBQPrescImg of JLBQOutput to the final fertilizer prescription image PrescImg;

[0249] JLBQInput includes: UAVList, SensorList, PhenocalList, FertHistList, BaseAmount, JLBQClusterK, PrescMin, and PrescMax;

[0250] The value of each pixel in PrescImg represents the amount of fertilizer applied at the corresponding spatial location;

[0251] S502. Create a cumulative fertilizer application list for each cluster region, PrescStatList = JLBQCluster, containing K all-zero elements; create a pixel count list for each cluster, PrescPixList = JLBQCluster, containing K all-zero elements.

[0252] S503. Initialize the output row traversal counter OutCounter1=1; Initialize the output column traversal counter OutCounter2=1;

[0253] S504. Take the cell value of PrescImg in the OutCounter1 row and OutCounter2 column and store it in the temporary storage OutTemp1; take the cluster label OutTemp2 at the corresponding position of JLBQLabelMap;

[0254] Let PrescStatList[OutTemp2] = PrescStatList[OutTemp2] + OutTemp1;

[0255] Set PrescPixList[OutTemp2] = PrescPixList[OutTemp2] + 1;

[0256] Where PrescStatList[OutTemp2] is the element value corresponding to OutTemp2 in PrescStatList;

[0257] S505, Set OutCounter2 = OutCounter2 + 1; If OutCounter2 is less than or equal to FieldWidth, go to S504; otherwise, set OutCounter1 = OutCounter1 + 1; If OutCounter1 is less than or equal to FieldHeight, go to S503; otherwise, go to S506.

[0258] S506. Establish a partitioned statistical output counter StatCounter=0;

[0259] S507. Calculate the average fertilizer application for the StatCounter-th cluster: PrescStatAvg = PrescStatList[StatCounter] / Max(PrescPixList[StatCounter], 1); Output the number of pixels in the StatCounter-th cluster: PrescPixList[StatCounter] and the average fertilizer application: PrescStatAvg; Set StatCounter = StatCounter + 1; If StatCounter is less than or equal to JLBQClusterK-1, go to S507; otherwise, go to S508.

[0260] Where PrescStatList[StatCounter] is the value of the StatCounter-th element in PrescStatList, and Max(PrescPixList[StatCounter],1) is to get the larger value between PrescPixList[StatCounter] and 1;

[0261] S508. Output PrescImg as the final generated fertilizer prescription image, stored in a floating-point single-band raster image format. The value of each pixel is the fertilizer amount (kg per hectare) at the corresponding spatial location. The pixel resolution is PrescReso meters, and the coordinate system is consistent with the first phase of UAVList data. Each pixel in the prescription image has a continuous floating-point fertilizer amount value, which can be directly used as a pixel-by-pixel variable fertilizer instruction image for UAVs or precision fertilization equipment.

Claims

1. A method for generating variable fertilization prescription images combining periodicity and phenology, characterized in that: The specific process of the method is as follows: The phenological history accumulation and vegetation state vectorization module JLBQModel is called to generate fertilization prescription images, the results of each partition of the fertilization prescription images are statistically analyzed, and the final fertilization prescription image is output. The input JLBQInput of the phenological history accumulation vegetation state vectorization module includes: UAV multispectral data time series table UAVList, ground wireless sensor data table SensorList, crop phenological history table PhenocalList, fertilization history table FertHistList, base fertilizer amount BaseAmount, cluster number JLBQClusterK, minimum fertilizer amount PrescMin, and maximum fertilizer amount PrescMax; the output of JLBQModel is the prescription image generation output JLBQOutput; The fertilizer prescription image PrescImg is the pixel fertilizer amount image JLBQPrescImg field of the JLBQOutput generated from the prescription image; The process of statistically analyzing the results of each partition of the fertilizer prescription image and outputting the final fertilizer prescription image specifically involves: obtaining the precise amount of fertilizer for each pixel based on the total amount of fertilizer applied to each region and the total number of pixels in each region of the fertilizer prescription image PrescImg, thereby obtaining the final fertilizer prescription image. The phenological history accumulation vegetation state vectorization module JLBQModel is obtained through the following method: S1. Construct UAVList, SensorList, PhenocalList, and FertHistList along with global control parameters, specifically: Each record in UAVList contains the following fields: UAVDate: the data collection date number calculated from the sowing date of the current season; UAVImg: the multispectral data collected by the UAV; UAVBandList: a list of band names; Each record in SensorList contains the following fields: SensorDate: the observation date number calculated from the sowing date of the current season; SensorID: the sensor node number; SensorLon: the node longitude; SensorLat: the node latitude; SensorSM: soil moisture; SensorEC: soil electrical conductivity; SensorN: nitrogen metabolism index; Each record in PhenocalList contains the following fields: PhenStage: Phenological stage number; PhenDateStart: Phenological stage start date; PhenDateEnd: Phenological stage end date; PhenAnchorDate: Center date of the fertilization window for the current phenological stage; PhenFertWindow: Whether the current phenological stage contains a fertilization window (PhenFertWindow=1 indicates the current phenological stage contains a fertilization window, otherwise PhenFertWindow=0 indicates the current phenological stage does not contain a fertilization window); PhenWeight: Phenological weight for the current phenological stage. Each record in FertHistList contains the following fields: FertDate: Fertilization date number; FertAmount: Fertilization amount; FertType: Fertilization type number; The global control parameters include: FieldHeight (number of cell rows in the field height direction), FieldWidth (number of cell columns in the field width direction), PrescReso (prescription cell resolution), BaseAmount (basic fertilizer application amount), Sensor spatial search radius ZQMZRadius (period anchoring time sequence Gaussian width ZQMZPeriodSigma), Number of clusters JLBQClusterK, Minimum fertilizer application PrescMin, and Maximum fertilizer application PrescMax. S2. Use UAVList, PhenocalList, SensorList and global control parameters to obtain the air-to-ground joint generation input ZQMZInput, and establish a periodically anchored air-to-ground timing joint generation module ZQMZModel. The input of ZQMZModel is ZQMZInput, and the output is the air-to-ground joint generation output ZQMZOutput. S3. Use UAVList, SensorList, PhenocalList, FertHistList and global control parameters to obtain the phenological accumulation input WHSJInput, construct the phenological history accumulation vegetation state vectorization module WHSJModel, the input of WHSJModel is WHSJInput, and WHSJModel calls ZQMZModel to obtain the phenological accumulation output WHSJOutput. S4. Use WHSJOutput and global control parameters to obtain the prescription image generation input JLBQInput, construct the cluster label weighted prescription image generation module JLBQModel, the input of JLBQModel is JLBQInput, and JLBQModel calls WHSJModel to obtain the prescription image generation output JLBQOutput.

2. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 1, characterized in that: In step S2, the air-to-ground joint generation input ZQMZInput is obtained using UAVList, PhenocalList, SensorList, and global control parameters. A periodically anchored air-to-ground timing joint generation module ZQMZModel is then established. The input to ZQMZModel is ZQMZInput, and the output is the air-to-ground joint generation output ZQMZOutput. Specifically: S201. Establish a periodic anchored air-to-ground time series joint generation module ZQMZModel. The input of ZQMZModel is the air-to-ground joint generation input ZQMZInput. The ZQMZInput includes: a UAVList record, the PhenAnchorDate of the phenological stage corresponding to the current UAVList record, a subset of SensorList within the time window of the current phenological stage, ZQMZRadius, and ZQMZPeriodSigma; S202. Add the ZQMZFeatMap field to the air-ground joint generation output ZQMZOutput, and initialize the ZQMZFeatMap field of the air-ground joint generation output ZQMZOutput to a 6-dimensional array of all zeros with FieldHeight row and FieldWidth column. S203. Initialize row traversal counter ZQMZCounter1=1; S204. Initialize column traversal counter ZQMZCounter2=1; S205. Extract the 5-dimensional pixel value of UAVImg from the ZQMZCounter1 row and ZQMZCounter2 column. Use the extracted 5-dimensional pixel value to obtain the normalized vegetation index temporary value ZQMZNdvi and the red-edge vegetation index temporary value ZQMZNdre. Specifically: ZQMZNdvi=(ZQMZTemp2-ZQMZTemp1) / (ZQMZTemp2+ZQMZTemp1+0.001); ZQMZNdre=(ZQMZTemp3-ZQMZTemp2) / (ZQMZTemp3+ZQMZTemp2+0.001); Among them, ZQMZTemp1, ZQMZTemp2 and ZQMZTemp3 are all temporary values, which are the R-band, NIR-band and RE-band values ​​in the 5-dimensional pixel values, respectively. ZQMZNdvi is the temporary value of the normalized vegetation index, and ZQMZNdre is the temporary value of the red edge vegetation index. S206. Based on the row and column number of the current pixel and the geographical range of the field, calculate the real geographical latitude and longitude corresponding to the current pixel through linear mapping to obtain the temporary longitude ZQMZCurLon and the temporary latitude ZQMZCurLat. S207. Initialize the weighted and cumulative soil moisture ZQMZWSumSM=0, the weighted and cumulative soil electrical conductivity ZQMZWSumEC=0, the weighted and cumulative nitrogen metabolism index ZQMZWSumN=0, and the total weighted and cumulative ZQMZWSum=0; initialize the sensor weight counter ZQMZCounter3=1. S208. Use the ZQMZCounter3 record of SensorList, ZQMZCurLon, and ZQMZCurLat to obtain the distance ZQMZDist between the sensor and the current pixel. If ZQMZDist is greater than ZQMZRadius, go to S210; otherwise, obtain the temporal Gaussian weight ZQMZTempW, update ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZTempW and ZQMZDist, and then execute S209. S209. Set ZQMZCounter3 = ZQMZCounter3 + 1; if ZQMZCounter3 is less than or equal to SensorTotal, go to S208; otherwise, go to S210. Where SensorTotal is the total number of observation records in SensorList; S210. If ZQMZWSum is less than 0.001, then set the current cell soil moisture value after spatiotemporal weighted interpolation to ZQMZEstSM=0, the current cell soil electrical conductivity value after spatiotemporal weighted interpolation to ZQMZEstEC=0, and the current cell nitrogen metabolism index value after spatiotemporal weighted interpolation to ZQMZEstN=0; otherwise, let ZQMZEstSM=ZQMZWSumSM / ZQMZWSum; ZQMZEstEC=ZQMZWSumEC / ZQMZWSum; let ZQMZEstN=ZQMZWSumN / ZQMZWSum. S211, Calculate the sensor-enhanced nitrogen metabolism production: ZQMZEstNAdj=ZQMZEstN×(0.7+0.3×Min(ZQMZEstEC / 2.0,1.0)); Wherein, Min represents the smaller of the two values; S212. Assign values ​​to each dimension of the ZQMZFeatMap in the ZQMZCounter1 row and ZQMZCounter2 column: 1st dimension = ZQMZNdvi; 2nd dimension = ZQMZNdre; 3rd dimension = ZQMZTempW; 4th dimension = ZQMZEstSM; 5th dimension = ZQMZEstEC; 6th dimension = ZQMZEstNAdj; Where ZQMZTempW is the temporal Gaussian weight; S213. Set ZQMZCounter2 = ZQMZCounter2 + 1; if ZQMZCounter2 is less than or equal to FieldWidth, go to S205; otherwise, go to S214. S214. Set ZQMZCounter1 = ZQMZCounter1 + 1; if ZQMZCounter1 is less than or equal to FieldHeight, go to S204; otherwise, go to S215. S215. Output ZQMZOutput as the result of ZQMZModel.

3. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 2, characterized in that: The distance ZQMZDist between the sensor and the current pixel in S208 is specifically as follows: ZQMZDist = ABS(ZQMZCurLon - SensorLon of the 3rd record in SensorList) × FieldWidth + ABS(ZQMZCurLat - SensorLat of the 3rd record in SensorList) × FieldHeight; Where ABS is the absolute value; Update ZQMZWSumSM, ZQMZWSumEC, ZQMZWSumN, and ZQMZWSum using ZQMZDist, specifically as follows: ZQMZWSumSM = ZQMZWSumSM + ZQMZW × SensorSM of the 3rd ZQMZCounter item in SensorList; ZQMZWSumEC = ZQMZWSumEC + ZQMZW × SensorEC of the 3rd ZQMZCounter item in SensorList; ZQMZWSumN = ZQMZWSumN + ZQMZW × SensorN of the 3rd ZQMZCounter item in SensorList; ZQMZWSum=ZQMZWSum+ZQMZW; ZQMZW=ZQMZTempW×ZQMZDistW; ZQMZDistW=1.0 / (ZQMZDist+1.5); ZQMZTempW = exp(-((SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput) × (SensorDate of the 3rd record in SensorList - PhenAnchorDate of ZQMZInput)) / (2.0 × ZQMZPeriodSigma × ZQMZPeriodSigma)); Where exp is an exponential function with base e, and ZQMZW is the combined weight.

4. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 3, characterized in that: In step S3, the phenological accumulation input WHSJInput is obtained using UAVList, SensorList, PhenocalList, FertHistList, and global control parameters. A phenological history accumulation vegetation state vectorization module WHSJModel is then constructed. The input to WHSJModel is WHSJInput, and WHSJModel calls ZQMZModel to obtain the phenological accumulation output WHSJOutput. Specifically: S301. Obtain the phenological accumulation input WHSJInput, and construct the phenological history accumulation vegetation state vectorization module WHSJModel. The input of WHSJModel is the phenological accumulation input WHSJInput. The phenological accumulation input WHSJInput includes: all records in UAVList, all records in SensorList, all records in PhenocalList, all records in FertHistList, field height (number of rows), field width (number of columns), sensor spatial search radius (ZQMZRadius), and periodic anchoring time sequence Gaussian width (ZQMZPeriodSigma). S302. Use the pixel-by-pixel state vector image WHSJStateMap and the last phase joint fusion feature image WHSJLastFeat to form the phenological accumulation output WHSJOutput, and initialize WHSJOutput; then establish the phenological weighted cumulative total image WHSJWtotalMap=FieldHeight row FieldWidth column all zero two-dimensional array; Initialize the phenological accumulation output WHSJOutput as a WHSJStateMap = an 8-dimensional array of all zeros with FieldHeight rows and FieldWidth columns; Initialize the phenological accumulation output WHSJOutput into a 6-dimensional array of all zeros, consisting of the last phase joint fusion feature image WHSJLastFeat = FieldHeight rows and FieldWidth columns; S303. Establish a phenological stage traversal counter WHSJCounter1=1; S304. Retrieve the WHSJCounter1 record from PhenocalList and store it in temporary storage WHSJTemp1; retrieve the record from UAVList whose UAVDate is closest to the PhenAnchorDate of WHSJTemp1 and store it in WHSJTemp2; retrieve all records from SensorList whose SensorDate is within the range of PhenDateStart to PhenDateEnd of WHSJTemp1 and store them in the sensor subset WHSJSensorSub. S305. Fill the parameters of UAVImg of WHSJTemp2, PhenAnchorDate of WHSJTemp1, WHSJSensorSub, ZQMZRadius and ZQMZPeriodSigma into ZQMZInput, call ZQMZModel, obtain the fusion feature image ZQMZOutput of the current stage and store it in the temporary storage WHSJTemp3. S306. Initialize the cell row counter WHSJCounter2=1; Initialize the cell column counter WHSJCounter3=1; S307. Extract the 6-dimensional values ​​of ZQMZFeatMap in the WHSJCounter2 row and WHSJCounter3 column of WHSJTemp3, and denote the 1st dimension data, 2nd dimension data, and 4th to 6th dimension data as WHSJCurNdvi, WHSJCurNdre, WHSJCurSM, WHSJCurEC and WHSJCurN respectively. S308. Use PhenWeight of WHSJTemp1 as the phenological weight WHSJPhenW for the current stage, update each dimension of WHSJStateMap with WHSJCounter2 and WHSJCounter3, and simultaneously assign WHSJPhenW to the corresponding row and column position of WHSJWtotalMap. S309. If WHSJCounter1 equals PhenTotal, then set the value of WHSJLastFeat in the WHSJCounter2 row and WHSJCounter3 column to the value of ZQMZFeatMap of WHSJTemp3 in the WHSJCounter2 row and WHSJCounter3 column, and then execute S310; otherwise, execute S310 directly. Wherein, PhenTotal is the number of phenological stages in PhenocalList; S310. Set WHSJCounter3 = WHSJCounter3 + 1. If WHSJCounter3 is less than or equal to FieldWidth, go to S307. Otherwise, set WHSJCounter2 = WHSJCounter2 + 1. If WHSJCounter2 is less than or equal to FieldHeight, go to S307. Otherwise, go to S311. S311. Set WHSJCounter1 = WHSJCounter1 + 1; if WHSJCounter1 is less than or equal to PhenTotal, go to S304; otherwise, go to S312. S312. Normalize the first to fifth dimensions of WHSJStateMap according to the corresponding position values ​​of WHSJWtotalMap: Traverse each pixel location and read the cumulative total of phenological weights WHSJWtotalMap corresponding to each location; if the WHSJWtotalMap of the current location is less than 0.001, then set the values ​​of the first to fifth dimensions of the current pixel's WHSJStateMap to 0; otherwise, divide the values ​​of the first to fifth dimensions of the current pixel's WHSJStateMap by the total weight value of the WHSJWtotalMap at the same location to obtain the normalized feature value after weighted average. S313. Establish a fertilizer application history response counter WHSJCounter4=1; establish a fertilizer application attenuation cumulative amount WHSJFertAccum=0; S314. Retrieve the 4th record of FertHistList (WHSJCounter4), take the last phenological stage PhenDateEnd as the reference date WHSJRefDate, and calculate the number of days since today: WHSJDeltaDays = WHSJRefDate - FertDate of the 4th record of FertHistList (WHSJCounter4). If WHSJDeltaDays is less than 0, set WHSJDeltaDays = 0; otherwise, keep the original value unchanged. Let WHSJFertAccum = WHSJFertAccum + FertAmount × exp(-WHSJDeltaDays / 30.0) of the 4th record of FertHistList (WHSJCounter4). Let WHSJCounter4 = WHSJCounter4 + 1. If WHSJCounter4 is less than or equal to FertTotal, go to S314; otherwise, go to S315. Where FertTotal is the total number of fertilization events in FertHistList; S315. Assign the 6th dimension of all pixels in WHSJStateMap to WHSJFertAccum / 200.0; if the 6th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 7th dimension of all pixels in WHSJStateMap to UAVTotal / Max(PhenTotal,1); if the 7th dimension is greater than 1.0, set it to 1.0, otherwise keep the original value unchanged; assign the 8th dimension of all pixels in WHSJStateMap to tanh(WHSJFertAccum / 100.0), and then update the 1st to 5th dimensions of WHSJStateMap in WHSJOutput: WHSJStateMap1st to 5th dimensions = arctan(WHSJStateMap1st to 5th dimensions × 2.0) / (π / 2); output WHSJOutput as the result of WHSJModel; Where Max is the larger of the two values, tanh is the hyperbolic tangent function, arctan is the arctangent function, π is pi, and UAVTotal is the number of data periods in UAVList.

5. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 4, characterized in that: In step S308, WHSJPhenW is used to update the data of each dimension of WHSJStateMap in the WHSJCounter2 row and WHSJCounter3 column, specifically as follows: Updated first-dimensional data = original first-dimensional data + WHSJPhenW × WHSJCurNdvi; Updated second-dimensional data = Unupdated second-dimensional data + WHSJPhenW × WHSJCurNdre; Updated 3rd dimension data = Unupdated 3rd dimension data + WHSJPhenW × WHSJCurSM; Updated fourth dimension data = original fourth dimension data + WHSJPhenW × WHSJCurEC; The updated 5th dimension data = the original 5th dimension data + WHSJPhenW × WHSJCurN.

6. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 5, characterized in that: In step S4, the prescription image generation input JLBQInput is obtained using WHSJOutput and global control parameters. A clustered label-weighted prescription image generation module JLBQModel is then constructed. The input to JLBQModel is JLBQInput. JLBQModel calls WHSJModel to obtain the prescription image generation output JLBQOutput. Specifically: S401. Construct a clustered label weighted prescription image generation module JLBQModel, where the input of JLBQModel is JLBQInput; The prescription image generation input JLBQInput includes: the entire contents of WHSJOutput, BaseAmount, JLBQClusterK, PrescMin, and PrescMax; S402. Using UAVList, SensorList, PhenocalList, FertHistList and the global constructor WHSJInput, call WHSJModel to obtain the temporary storage JLBQTemp1 of WHSJOutput; S403. Expand the WHSJStateMap of JLBQTemp1 into a two-dimensional feature matrix JLBQFeatMat by rows and columns; S404. Perform K-means clustering on JLBQFeatMat based on JLBQClusterK to obtain the cluster label image JLBQLabelMap and the cluster center matrix JLBQCenters for each pixel. S405. Calculate the mean of the first dimension of data for each cluster center, sort the cluster numbers according to the mean of the first dimension of data from smallest to largest, and obtain the list of cluster numbers in ascending order JLBQSortedClust. S406. Create a fertilizer weight list JLBQWeightList; create a weight assignment counter JLBQCounter1=1; set the minimum fertilizer weight JLBQWMin=0.60 and the total fertilizer weight span JLBQWRange=0.70; S407, the cluster weight of the JLBQSortedClust[JLBQCounter1] in JLBQWeightList = JLBQWMin + JLBQWRange × (JLBQCounter1-1) / (JLBQClusterK-1); Set JLBQCounter1 = JLBQCounter1 + 1; if JLBQCounter1 is less than or equal to JLBQClusterK, go to S407, otherwise go to S408; Where JLBQSortedClust[JLBQCounter1] is the value of the JLBQCounter1 element in JLBQSortedClust; S408. Calculate the first-dimensional mean value of the entire field JLBQTemp1, JLBQNDVIMean = Sum of the first-dimensional values ​​of all cells / (FieldHeight × FieldWidth); S409. Initialize the prescription image and generate the pixel fertilizer amount image JLBQOutput JLBQPrescImg field = FieldHeight row FieldWidth column single-band all-zero floating-point raster image; establish the prescription row counter JLBQCounter2=1; S410. Set the prescription column counter JLBQCounter3=1; S411. Using the WHSJLastFeat data of JLBQTemp1 in the first dimension of the JLBQCounter2 row and JLBQCounter3 column and the value of JLBQLabelMap in the JLBQCounter2 row and JLBQCounter3 column, obtain the prescription image value JLBQPrescrip of the current pixel, and assign the pixel value of JLBQPrescImg in the JLBQCounter2 row and JLBQCounter3 column to JLBQPrescrip; S412, Set JLBQCounter3 = JLBQCounter3 + 1; If JLBQCounter3 is less than or equal to FieldWidth, go to S411; Otherwise, Set JLBQCounter2 = JLBQCounter2 + 1; If JLBQCounter2 is less than or equal to FieldHeight, go to S410; Otherwise, go to S413. S413. Output JLBQOutput as the result of JLBQModel.

7. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 6, characterized in that: In S411, the prescription image value JLBQPrescrip of the current pixel is obtained by using the first dimension data of WHSJLastFeat in the JLBQCounter2 row and JLBQCounter3 column of JLBQTemp1 and the value of JLBQLabelMap in the JLBQCounter2 row and JLBQCounter3 column. Obtain the WHSJLastFeat of the latest NDVI value JLBQCurNdvi=JLBQTemp1 of the current pixel in the first dimension of the JLBQCounter2 row and JLBQCounter3 column; calculate the base amount spatial modulation JLBQBase=BaseAmount×(1.0-0.25×JLBQCurNdvi / JLBQNDVIMean); if JLBQBase is less than PrescMin×0.5, then set JLBQBase=PrescMin×0.5, otherwise keep JLBQBase unchanged; Get the cluster label JLBQCurLabel=JLBQLabelMap of the current cell in the JLBQCounter2 row and JLBQCounter3 column; calculate the prescription image value of the current cell JLBQPrescrip=JLBQBase×JLBQWeightList[JLBQCurLabel]; Where JLBQWeightList[JLBQCurLabel] is the element value of JLBQCurLabel in JLBQWeightList. If JLBQPrescrip is less than PrescMin, then set JLBQPrescrip=PrescMin; if JLBQPrescrip is greater than PrescMax, then set JLBQPrescrip=PrescMax; otherwise, keep JLBQPrescrip unchanged.

8. The method for generating variable fertilization prescription images combining periodicity and phenology according to claim 7, characterized in that: The process involves calling the phenological history accumulation vegetation state vectorization module JLBQModel to generate a fertilizer prescription image, statistically analyzing the results of each partition of the fertilizer prescription image, and outputting the final fertilizer prescription image. Specifically: A1. Call JLBQModel, input JLBQInput into JLBQModel, and obtain JLBQOutput; assign the pixel fertilizer application image JLBQPrescImg of JLBQOutput to the final fertilizer prescription image PrescImg; A2. Create a cumulative fertilizer application list for each cluster region, PrescStatList = JLBQCluster, containing K zero-element lists; create a pixel count list for each cluster, PrescPixList = JLBQCluster, containing K zero-element lists. A3. Initialize the output row traversal counter OutCounter1=1; Initialize the output column traversal counter OutCounter2=1; A4. Take the cell value of PrescImg in the OutCounter1 row and OutCounter2 column and store it in the temporary storage OutTemp1; take the cluster label OutTemp2 at the corresponding position of JLBQLabelMap; Let PrescStatList[OutTemp2] = PrescStatList[OutTemp2] + OutTemp1; Set PrescPixList[OutTemp2] = PrescPixList[OutTemp2] + 1; Where PrescStatList[OutTemp2] is the element value corresponding to OutTemp2 in PrescStatList; A5. Set OutCounter2 = OutCounter2 + 1; if OutCounter2 is less than or equal to FieldWidth, go to A4; otherwise, set OutCounter1 = OutCounter1 + 1; if OutCounter1 is less than or equal to FieldHeight, go to A3; otherwise, go to A6. A6. Set the partitioned statistics output counter StatCounter=0; A7. Calculate the average fertilizer application for the StatCounter-th cluster: PrescStatAvg = PrescStatList[StatCounter] / Max(PrescPixList[StatCounter], 1); Output the number of pixels in the StatCounter-th cluster: PrescPixList[StatCounter] and the average fertilizer application: PrescStatAvg; Set StatCounter = StatCounter + 1; If StatCounter is less than or equal to JLBQClusterK-1, go to A7; otherwise, go to A8. Where PrescStatList[StatCounter] is the value of the StatCounter-th element in PrescStatList, and Max(PrescPixList[StatCounter],1) is to get the larger value between PrescPixList[StatCounter] and 1; A8. Output PrescImg as the final fertilizer prescription image.

Images