Early diagnosis method and system for crop nutrient stress based on time-series remote sensing images

Abstract

The present application provides a crop nutrient stress early diagnosis method and system based on time sequence remote sensing image, belonging to the technical field of agricultural remote sensing and crop growth monitoring, comprising: loading a multi-temporal remote sensing image sequence covering the whole growth period of the target field, calculating the time sequence spectral vegetation index pixel by pixel after preprocessing; loading a pre-constructed phenology knowledge graph, retrieving the corresponding normal growth curve template according to the crop variety and geographical location of the target field; constructing a time sequence differential attention network; inputting the spatio-temporal spectral data cube and the normal growth curve template into the time sequence differential attention network, outputting the nutrient stress probability value; determining the nutrient stress grade of the target field according to the nutrient stress probability value, and generating a stress grade spatial distribution map. The present application solves the problem that the existing method is difficult to separate the crop's own phenological change from the spectral difference caused by nutrient stress in time sequence remote sensing data, resulting in the inability to achieve high-precision early diagnosis before the stress symptoms appear.

Description

Technical Field

[0001] This invention relates to the field of agricultural remote sensing and crop growth monitoring technology, specifically to a method and system for early diagnosis of crop nutrient stress based on time-series remote sensing images. Background Technology

[0002] Currently, precision agriculture places increasingly higher demands on the timeliness and accuracy of crop nutrient management. Early identification of nutrient stress states such as nitrogen, phosphorus, and potassium in crops allows for precise variable-rate fertilization before irreversible yield reductions occur. This is crucial for ensuring grain yields, reducing fertilizer application, and minimizing agricultural non-point source pollution. While traditional field sampling and laboratory chemical analysis methods are highly accurate, they suffer from destructiveness, latency, and sparse sampling points, making them unsuitable for large-scale, high-time-sensitivity monitoring. Therefore, crop growth and stress diagnosis technologies based on non-destructive remote sensing have emerged.

[0003] Existing technologies utilize single-period remote sensing imagery to construct statistical relationships between specific vegetation indices, such as the Normalized Difference Vegetation Index (NDVI) or red-edge parameters, and crop nitrogen content for nutrient diagnosis. However, static spectral information from a single period is highly susceptible to interference from other factors such as crop variety, growth stage, planting density, and soil background, leading to confusion of nutrient stress signals and poor universality of diagnostic models across different fields or growth stages. Furthermore, stress is a dynamic process, with a delay between initial physiological changes and spectral response changes; single-period images often only capture stress after symptoms appear, failing to provide true early warning. Another approach utilizes hyperspectral sensors to acquire continuous narrow-band information, which can capture more subtle spectral changes, but the high cost of data acquisition and the complexity of massive data processing limit its widespread application at the field scale.

[0004] Therefore, how to effectively separate the dynamic spectral differences caused by crop phenological changes and nutrient stress from easily accessible time-series remote sensing data, and construct a method that is applicable to different growth environments and can achieve high-precision early diagnosis before stress symptoms appear, is a key problem that urgently needs to be solved in this field. Summary of the Invention

[0005] This invention provides a method and system for early diagnosis of crop nutrient stress based on time-series remote sensing images, which solves the problem that existing methods are unable to separate the phenological changes of crops themselves from the spectral differences caused by nutrient stress from time-series remote sensing data, thus making it impossible to achieve high-precision early diagnosis before stress symptoms appear.

[0006] In view of the above problems, the present invention provides a method and system for early diagnosis of crop nutrient stress based on time-series remote sensing images.

[0007] In a first aspect, the present invention provides a method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery, comprising:

[0008] Load multi-temporal remote sensing image sequences covering the entire growth period of the target field, and calculate the temporal spectral vegetation index pixel by pixel after preprocessing to obtain a spatiotemporal spectral data cube. Load the pre-built phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field; A temporal differential attention network is constructed, wherein the temporal differential attention network includes a phenological gating unit and a differential comparison unit; The spatiotemporal spectrum data cube and the normal growth curve template are input into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve, and outputs the nutrient stress probability value after self-attention encoding. The nutrient stress probability value is used as the output of the temporal differential attention network. The nutrient stress level of the target field is determined based on the nutrient stress probability value, and a spatial distribution map of the stress level is generated.

[0009] Secondly, the present invention provides an early diagnosis system for crop nutrient stress based on time-series remote sensing imagery, comprising: The spatiotemporal spectral data cube construction module is used to load multi-temporal remote sensing image sequences covering the entire growth period of the target field, and calculate the temporal spectral vegetation index pixel by pixel after preprocessing to obtain the spatiotemporal spectral data cube. The normal growth curve template retrieval module is used to load a pre-built phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field. A temporal differential attention network construction module is used to construct a temporal differential attention network, wherein the temporal differential attention network includes a phenological gating unit and a differential comparison unit; The nutrient stress probability calculation module is used to input the spatiotemporal spectrum data cube and the normal growth curve template into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve, and outputs the nutrient stress probability value after self-attention encoding. The stress level distribution map generation module is used to use the nutrient stress probability value as the output of the temporal differential attention network, determine the nutrient stress level of the target field based on the nutrient stress probability value, and generate a stress level spatial distribution map.

[0010] One or more technical solutions provided in this invention have at least the following technical effects or advantages: The technical solution of this invention first loads a multi-temporal remote sensing image sequence covering the entire growth period of the target field. After preprocessing, four spectral vegetation indices are calculated pixel by pixel, and stacked along the time dimension to form a four-dimensional spatiotemporal spectral data cube containing spatial, temporal, and spectral index dimensions. This integrates the scattered multi-phase remote sensing images into a structured, pixel-by-pixel indexable temporal spectral dataset, ensuring that each spatial pixel obtains a complete multi-dimensional feature curve recording the dynamic spectral changes throughout the entire growth period, providing a unified data foundation for subsequent differential comparative analysis.

[0011] Furthermore, a pre-constructed phenological knowledge graph is loaded, and the corresponding normal growth curve template is retrieved based on the crop variety and geographical location of the target field. A personalized normal growth reference benchmark adapted to its variety characteristics and regional climate conditions is provided for each target field, eliminating the inherent influence of variety differences and geographical environment differences on the spectral curve, so that subsequent differential comparisons can focus on the real abnormal deviations caused by nutrient stress.

[0012] Furthermore, a temporal differential attention network comprising phenological gating units and differential contrast units is constructed. The phenological gating mechanism guides the network to focus on stress-sensitive critical growth periods, multi-level differential features capture the cumulative effects and trend reversals of stress, and physical constraint loss ensures that the prediction results conform to radiative transfer laws, enabling the network to distinguish between normal crop phenological changes and spectral anomalies caused by nutrient stress.

[0013] Furthermore, the spatiotemporal spectral data cube and the normal growth curve template are jointly input into a temporal differential attention network. After self-attention encoding, the fully connected classification head outputs a four-class discrete probability distribution. This achieves an end-to-end mapping from temporal remote sensing data to pixel-level stress probabilities, and the self-attention mechanism effectively captures the transmission and amplification relationship of stress differences in the temporal dimension.

[0014] Finally, the stress level corresponding to the maximum probability value of the discrete probability distribution output by the network is taken as the diagnostic result for each pixel, and rasterized according to the pixel spatial coordinates to generate a spatial distribution map of stress levels. The time series analysis results abstracted by the network are transformed into intuitive field zoning prescription maps, providing operational unit-level decision-making basis for precision variable topdressing.

[0015] In summary, the technical solution of this invention can achieve high-precision early identification of nutrient stress before the symptoms of stress appear, providing a prescription basis at the operational unit level for precise variable topdressing, reducing the amount of chemical fertilizer applied while ensuring crop yield, thus achieving both economic and ecological benefits. Attached Figure Description

[0016] Figure 1This is a flowchart illustrating the method for early diagnosis of crop nutrient stress based on time-series remote sensing images provided by the present invention.

[0017] Figure 2 This is a schematic diagram of the temporal differential attention network in the crop nutrient stress early diagnosis method based on temporal remote sensing images provided by the present invention.

[0018] Figure 3 This is a schematic diagram of the structure of the crop nutrient stress early diagnosis system based on time-series remote sensing images provided by the present invention.

[0019] In the attached diagram, the labels representing each component are as follows: The module includes: 11 Spatiotemporal spectral data cube construction module; 12 Normal growth curve template retrieval module; 13 Temporal differential attention network construction module; 14 Nutrient stress probability calculation module; and 15 Stress level distribution map generation module. Detailed Implementation

[0020] This invention provides a method and system for early diagnosis of crop nutrient stress based on time-series remote sensing images. This solves the problem that existing methods are unable to separate the spectral differences caused by crop phenological changes and nutrient stress from time-series remote sensing data, thus making it impossible to achieve high-precision early diagnosis before stress symptoms appear.

[0021] It should be noted that the terms "comprising" and "having" are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or server that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or modules that are not explicitly listed or that are inherent to these processes, methods, products, or devices.

[0022] Example 1, as Figure 1 As shown, this invention provides a method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery, the method comprising: S100: Load multi-temporal remote sensing image sequences covering the entire growth period of the target field, calculate the temporal spectral vegetation index pixel by pixel after preprocessing, and obtain a spatiotemporal spectral data cube.

[0023] This step integrates scattered, multi-period remote sensing images into a structured time-series dataset, providing a unified and comparable input for subsequent network models. Coverage throughout the entire growth period ensures that the complete dynamic process of crops from normal growth to stress and then to stress recovery is recorded, while the pixel-by-pixel processing method guarantees that each field unit at the spatial scale has independent and continuous temporal spectral characteristics, which is the data foundation for achieving pixel-level early stress diagnosis.

[0024] In this step, the first step is to load a multi-temporal remote sensing image sequence covering the entire growth period of the target field, including: Acquire satellite imagery that includes the red-edge band; Each image undergoes radiometric calibration, atmospheric correction, and cloud mask preprocessing to obtain a preprocessed image. Based on the preprocessed image, a spectral vegetation index is calculated pixel by pixel. The spectral vegetation index includes at least the normalized difference red edge index, the red edge chlorophyll index, the improved chlorophyll uptake ratio index, and the normalized difference vegetation index.

[0025] In this step, the first step is to acquire a satellite image sequence covering the entire growth period of the target field, including the red-edge bands. The satellite image source is, for example, the Sentinel-2 satellite, whose multispectral imager has three red-edge bands: B5 (705nm), B6 ​​(740nm), and B7 (783nm), with a spatial resolution of 20 meters and a revisit period of 5 days, meeting the spatiotemporal resolution requirements for time-series monitoring of crops throughout their entire growth period. The image sequence spans from one week before sowing to one week after the expected harvest, ensuring complete coverage of the spectral changes throughout the entire growth period.

[0026] For example, taking winter wheat field A as the monitoring target, its sowing date was October 15, 2024, and the expected harvest date was June 10, 2025. Sentinel-2 L1C level images covering all fields A with cloud cover below 40% from October 8, 2024 to June 17, 2025 were downloaded from the Copernicus Open Access Center, resulting in a total of 32 valid images.

[0027] Furthermore, each acquired satellite image undergoes radiometric calibration, atmospheric correction, and cloud masking preprocessing sequentially. Radiometric calibration converts the image's digital quantization values ​​into physically meaningful apparent radiance or apparent reflectance; atmospheric correction eliminates interference from atmospheric molecules and aerosols on electromagnetic wave transmission, obtaining the true surface reflectance; cloud masking utilizes the image's built-in cloud classification bands or QA bands to identify and remove pixels affected by clouds, cloud shadows, and snow, preventing outliers from interfering with subsequent analysis. Finally, a reliable surface reflectance image sequence is obtained scene-by-scene and pixel-by-pixel. All of the above preprocessing can be performed using the Sen2Cor processor.

[0028] For example, atmospheric correction is performed on each of the 32 Sentinel-2 L1C images using the Sen2Cor processor to output L2A-level surface reflectance products. At the same time, the SCL scene classification band of each image is read, and pixels with a classification value of 8 (probability in clouds), 9 (probability of cloud height), 10 (probability of thin cirrus clouds), and 11 (probability of snow) are marked as invalid pixels and removed by masking in subsequent calculations.

[0029] Furthermore, multiple spectral vegetation indices are calculated pixel-by-pixel on each preprocessed image to construct a temporal spectral feature vector. The selected vegetation indices include at least: Normalized Difference Red Edge Index (NDRE), Red Edge Chlorophyll Index (CIre), Modified Chlorophyll Uptake Ratio Index (MCARI), and Normalized Difference Vegetation Index (NDVI). Each index characterizes crop chlorophyll content, canopy structure, and photosynthetic activity from different spectral dimensions. The red edge-related indices are more sensitive to chlorophyll changes and early stress, while NDVI is used to characterize overall greenness trends. The combination of multiple indices provides complementary spectral feature information for the network model, thereby improving the accuracy of nutrient stress diagnosis.

[0030] The formulas for calculating NDRE are: NDRE=(RNIR-RRE2) / (RNIR+RRE2); the formula for calculating the red-edged chlorophyll index CIre is: CIre=(RRE3 / RRE1)-1; the formula for calculating the improved chlorophyll absorption ratio index MCARI is: MCARI=[(RRE2-RRE1)-0.2×(RRE2-RRED)]×(RRE2 / RRE1); and the formula for calculating the normalized difference vegetation index NDVI is: NDVI=(RNIR-RRED) / (RNIR+RRED).

[0031] In the above physical quantities, RRED represents the surface reflectance of a pixel in the red band, RRE1 represents the surface reflectance of a pixel at the beginning of the red edge, RRE2 represents the surface reflectance of a pixel in the middle of the red edge, RRE3 represents the surface reflectance of a pixel at the end of the red edge, and RNIR represents the surface reflectance of a pixel in the near-infrared band. These correspond to B4, B5, B6, B7, and B8 in the following examples, respectively.

[0032] For example, taking the target pixel with coordinates (x1, y1) in field A as an example, on the pre-processed image on March 15, 2025, which is the jointing stage, the surface reflectance values ​​of each band are: B4=0.045, B5=0.072, B6=0.135, B7=0.162, B8=0.285. The calculations yield: NDRE = (B8-B6) / (B8+B6) = (0.285-0.135) / (0.285+0.135) = 0.357; CIre = (B7 / B5)-1 = (0.162 / 0.072)-1 = 1.250; MCARI = [(B6-B5)-0.2×(B6-B4)]×(B6 / B5) = [(0.135-0.072)-0.2×(0.135-0.045)]×(0.135 / 0.072) = 0.081; NDVI = (B8-B4) / (B8+B4) = (0.285-0.045) / (0.285+0.045) = 0.727. The four index values ​​of the above time phases are combined into a vector [0.357, 1.250, 0.081, 0.727] and stored in the time series feature sequence.

[0033] Finally, based on the preprocessed multi-temporal image sequence, each spatial location is traversed pixel by pixel, and the four spectral vegetation indices NDVI, NDRE, CIre, and MCARI for each corresponding pixel in each image are calculated. The calculation results of each temporal phase for the same pixel are arranged in chronological order to form a temporal spectral index matrix. After all pixels at all spatial locations have undergone the above processing, the temporal spectral index matrices of all pixels are combined along the spatial dimension to obtain a spatiotemporal spectral data cube covering the entire growth period of the target field.

[0034] The spatiotemporal spectral data cube is a four-dimensional data structure, with each dimension defined as follows: Spatial dimension, i.e., rows × columns: corresponding to the geographic spatial range of the target field, with each grid point representing a pixel, and the resolution being the spatial resolution of the input remote sensing image, such as 20 meters for Sentinel-2; Temporal dimension, containing all effective observation phases throughout the entire growth period, with the number of phases depending on the satellite revisit cycle and the number of effective images behind cloud masking, typically ranging from 20 to 50 phases; Spectral index dimension: containing four spectral vegetation indices calculated per pixel and per phase.

[0035] For example, taking winter wheat field A as an example, its spatial range contains M rows × N columns of effective pixels. For each pixel, the surface reflectance values ​​of five bands (B4, B5, B6, B7, and B8) are extracted one by one from all 32 preprocessed Sentinel-2 images, and four indices are calculated respectively. Taking the target pixel with coordinates (x1, y1) as an example, a set of four-dimensional index vectors is obtained for each of the 32 temporal phases. After being arranged chronologically, a 32-row × 4-column temporal spectral index matrix is ​​formed, that is, 32 temporal phases × 4 index types. The 32×4 matrices of all M×N pixels in field A are superimposed along the spatial dimensions to finally obtain a four-dimensional spatiotemporal spectral data cube with dimensions of M×N×32×4. Among them, the first dimension is the spatial coordinates in the row direction, the second dimension is the spatial coordinates in the column direction, the third dimension is time, and the fourth dimension is the spectral index type.

[0036] If, after cloud masking, the number of valid images of the target field throughout its entire growth period is lower than the preset minimum temporal threshold (e.g., less than 5 images), then the field is deemed not to meet the conditions for temporal analysis. A message indicating insufficient data quality is generated, and it is suggested that synthetic images from nearby temporal phases or ground sampling data from the current year be used for supplementary analysis.

[0037] It should be noted that the above values ​​are for illustrative purposes only and do not constitute a limitation on the present invention.

[0038] In summary, this step integrates scattered multi-period remote sensing images into a structured, indexable time-series spectral dataset, providing a unified data foundation for subsequent network models to conduct pixel-by-pixel and full-growth-cycle nutrient stress difference analysis.

[0039] S200: Load the pre-built phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field.

[0040] This step provides a personalized reference benchmark for normal crop growth for subsequent differential analysis. Because different crop varieties, such as winter wheat and summer maize, have vastly different phenological rhythms and spectral response patterns, even the same variety can exhibit significant differences in growth period length and growth curve morphology across different climatic zones. By accurately retrieving normal growth curve templates that match the target field variety and geographical location using a phenological knowledge graph, the interference of variety and regional differences on stress diagnosis can be eliminated.

[0041] In this step, a phenological knowledge graph is pre-constructed, including: The data collected from field trials over many years includes measured values ​​of chlorophyll content at key growth stages under different crop varieties and nitrogen gradients, as well as synchronously acquired spectral vegetation indices. Based on the measured chlorophyll content, samples with normal nutrient status were selected. Based on samples with normal nutrient status and prior knowledge of crop cultivation, the normal dynamic range and its rate of change confidence interval of the spectral vegetation index under each crop variety-ecological zone-growth period combination are labeled. The normal dynamic range includes an upper limit and a lower limit. The labeled knowledge is stored in a graph data structure, where the node attributes include the mean and standard deviation of the spectral vegetation index for the corresponding growth period, and the edge attributes include the upper and lower limits of the allowable rate of change of the spectral vegetation index between adjacent growth periods.

[0042] Specifically, field test sites were first set up in several typical ecological zones, and different nitrogen gradient treatments were set up for different crop varieties, including a control treatment with sufficient nitrogen supply and stress treatments with nitrogen deficiency at each gradient. Field observation data were collected for key growth periods for many consecutive years.

[0043] The field trial observation data includes: measured values ​​of chlorophyll content at key growth stages, such as the greening stage, jointing stage, booting stage, heading stage, grain-filling stage, and maturity stage of winter wheat, which can be obtained through SPAD chlorophyll meter or leaf chemical analysis method; and canopy spectral vegetation indices obtained simultaneously, which can be collected by ground spectrometer or UAV equipped with multispectral sensor to calculate indices such as NDVI, NDRE, CIre, and MCARI.

[0044] For example, three typical ecological sites were selected in the main winter wheat producing area of ​​the North China Plain, including Shijiazhuang in Hebei Province, Dezhou in Shandong Province, and Xinxiang in Henan Province. Five main winter wheat varieties—wheat A, wheat B, wheat C, wheat D, and wheat E—were used, with five nitrogen gradient treatments for each variety: N0 (no nitrogen application); N1 (75 kg / hm²); N2 (150 kg / hm²); N3 (225 kg / hm²); and N4 (300 kg / hm²). N4 served as a sufficient nitrogen control based on the locally recommended nitrogen application rate. For three consecutive years, leaf SPAD values ​​were measured at the greening, jointing, booting, heading, grain-filling, and maturity stages, and canopy multispectral images were simultaneously collected. A total of 4050 sets of valid sample data were obtained (3 stations × 5 varieties × 5 gradients × 6 growth stages × 3 replicates × 3 years). The above leaf SPAD values ​​are the measured chlorophyll content values ​​obtained using a SPAD chlorophyll meter.

[0045] Secondly, based on the measured chlorophyll content, samples with normal nutrient status were selected from all experimental samples. The specific criteria were as follows: Statistical analysis, such as analysis of variance, was performed on the measured chlorophyll content of each nitrogen gradient treatment under the same variety, ecological zone, and growth stage. Treatment combinations with the highest chlorophyll content and no significant difference from the adequate nitrogen control were selected, and their corresponding samples were marked as samples with normal nutrient status. For samples in the adequate nitrogen control treatment, if their growth was significantly better than other gradients and their chlorophyll content reached the local reference upper limit for that variety, they were directly marked as normal samples.

[0046] For example, using data from the wheat jointing stage (A) at the Shijiazhuang site, an analysis of variance was performed on the SPAD values ​​of the five nitrogen gradient treatments (N0 to N4). The results showed that the mean SPAD value after treatment N3 was 52.1, with no significant difference between this and the mean SPAD value after treatment N4 (53.4), p>0.05 (mean significance level α=0.05). Both were significantly higher than treatments N0, N1, and N2. Therefore, all replicates under treatments N3 and N4 were labeled as samples with normal nutrient status, resulting in a total of 2 gradients × 3 replicates = 6 groups of normal sample data. Treatments N0 to N2 were excluded because their chlorophyll content was significantly lower than normal levels.

[0047] Secondly, based on the selected samples with normal nutrient status, and according to prior knowledge of crop cultivation, the spectral vegetation indices of each normal sample were statistically analyzed using a three-level grouping dimension of crop variety, ecozone, and growth stage. The statistics included the mean and standard deviation of each spectral vegetation index within each growth stage, as well as the mean and standard deviation of the rate of change of each index between adjacent growth stages. Based on this, upper and lower limits of the normal dynamic range, as well as upper and lower limits of the allowable rate of change, were set.

[0048] The rules for setting the upper and lower limits of the normal dynamic range are as follows: Lower limit = Mean - Confidence coefficient × Standard deviation; Upper limit = Mean + Confidence coefficient × Standard deviation. The rules for setting the allowable upper and lower limits of the rate of change are as follows: Allowable upper limit of rate of change = Mean of rate of change + Confidence coefficient × Standard deviation of rate of change; Allowable lower limit of rate of change = Mean of rate of change - Confidence coefficient × Standard deviation of rate of change.

[0049] The preferred range for the confidence coefficient is 1.5 to 2.5, which can be set differently according to the characteristics of the variety: for varieties with long growth period and robust growth, the confidence coefficient value is set to 2.0; for varieties with short growth period and rapid growth, the confidence coefficient value is set to 1.5, so as to appropriately narrow the normal range and improve the sensitivity of stress detection.

[0050] For example, taking the NDRE index of wheat at the jointing stage (A) at the Shijiazhuang site, the mean NDRE of the normal sample is 0.452, and the standard deviation is 0.028. Setting a confidence coefficient of 2.0, the upper limit of the normal dynamic range of NDRE at the jointing stage is 0.452 + 2.0 × 0.028 = 0.508, and the lower limit is 0.452 - 2.0 × 0.028 = 0.396. Similarly, statistically analyzing the rate of change of NDRE from the jointing stage to the booting stage, the mean rate of change of NDRE in the normal sample during this stage is +0.038 / ten-day period, and the standard deviation is 0.012 / ten-day period. Therefore, the upper limit of the allowable rate of change is +0.038 + 2.0 × 0.012 = +0.062 / ten-day period, and the lower limit is +0.038 - 2.0 × 0.012 = +0.014 / ten-day period.

[0051] Finally, the labeled knowledge is stored and indexed using a graph data structure. Nodes in the graph are uniquely identified by variety-ecoregion-growth stage. Node attributes include the mean, standard deviation, upper limit, and lower limit of the normal dynamic range for each spectral vegetation index during that growth stage. Edges connect nodes in adjacent growth stages, and edge attributes include the mean, standard deviation, upper limit, and lower limit of the rate of change of each spectral vegetation index between the two growth stages. The top-level index dimensions of the graph are crop variety and ecoregion. During a query, the corresponding subgraph can be quickly located based on the input variety and geographical location, obtaining the normal growth curve template information carried by each node and edge throughout the entire growth stage.

[0052] For example, a node is stored in the graph, identified as [Variety = Wheat A, Ecological Region = Northern North China Plain, Growth Stage = Jointing Stage]. The NDRE attribute values ​​of the above node are: mean = 0.452, standard deviation = 0.028, lower limit of normal = 0.396, upper limit of normal = 0.508; the NDVI attribute values ​​are: mean = 0.785, standard deviation = 0.035, lower limit of normal = 0.715, upper limit of normal = 0.855; the attribute values ​​of other indices are stored similarly. The edge attribute pointing to [Growth Stage = Heading Stage] is: NDRE rate of change mean = +0.038 / ten-day, standard deviation = 0.012 / ten-day, lower limit of allowable = +0.014 / ten-day, upper limit of allowable = +0.062 / ten-day. When the target field is identified as wheat variety A, located near Shijiazhuang, Hebei, the graph substructure can be automatically matched, and the normal range and rate of change threshold of the index for the entire sequence from jointing stage to maturity can be output to form a complete normal growth curve template.

[0053] In this invention, when a perfectly matching normal growth curve template cannot be retrieved from the phenological knowledge graph for the target field, a transfer learning strategy is employed: Search for alternative templates that belong to the same agricultural climate zone as the target field and have the highest genetic similarity to the crop variety; Using the limited-temporal remote sensing image data of the target field from the seedling stage to the present, the alternative template is subjected to unsupervised domain adaptive calibration. By minimizing the feature distribution difference between the alternative template and the acquired limited-temporal remote sensing image data, the adapted normal growth curve template is obtained.

[0054] Specifically, when the combination of crop variety and geographical location of the target field cannot find a perfectly matching normal growth curve template in the phenological knowledge graph, a strategy of tiered relaxation of matching conditions is adopted to search for alternative templates. First, templates belonging to the same agricultural climate zone and containing the same crop variety as the target field are searched; if no match is found, then templates of varieties with the highest genetic similarity to the target variety within the same agricultural climate zone are searched. The agricultural climate zone refers to a homogeneous region defined based on agricultural climate indicators such as accumulated temperature, precipitation, and light intensity. Genetic similarity is comprehensively evaluated based on indicators such as pedigree, growth period length, and plant type characteristics.

[0055] For example, target field B is located in northern Jiangsu Province and is planted with winter wheat variety F. The phenological knowledge graph does not yet include a normal growth curve template for wheat F in this region. First, a search for templates for wheat F within the same agroclimatic zone was conducted, but no match was found. Then, the variety with the highest genetic similarity to wheat F in the Huang-Huai-Hai region was searched. Wheat A is included in the pedigree of wheat F, and the difference in growth period between the two is only 3 days. Both have a semi-compact plant type and the highest genetic similarity score. Therefore, a normal growth curve template for wheat A in southern Huang-Huai-Hai region was retrieved as a substitute template. The substitute template includes normal expected values ​​and rate of change thresholds for NDVI, NDRE, CIre, and MCARI at 32 time points.

[0056] Next, acquire the limited-temporal remote sensing image data accumulated for the target field from the seedling emergence period to the current diagnosis time. The limited-temporal phase refers to the effective remote sensing image observation phase retained after preprocessing and cloud masking within the time span from the seedling emergence period to the current time.

[0057] For example, the winter wheat emergence date for field B was November 5, 2024, and the current diagnosis time was March 15, 2025. A total of 14 Sentinel-2 images were acquired, and after preprocessing and cloud masking, 10 valid images were retained. The four spectral vegetation indices NDVI, NDRE, CIre, and MCARI for the 10 valid temporal phases were calculated pixel by pixel to obtain a finite-temporal spectral index matrix with a dimension of 10×4 for each pixel in the field.

[0058] Finally, using the limited temporal remote sensing image data already acquired for the target field, unsupervised domain adaptive calibration is performed on the alternative template to eliminate the feature distribution deviation between the alternative template and the target field caused by differences in variety characteristics, local climate, and soil conditions.

[0059] The specific method is as follows: take the first T' finite temporal phases of the replacement template as the source domain features, and take the field mean sequence of the finite temporal image data of the target field as the target domain features. By minimizing the feature distribution difference between the two, the replacement template is transformed linearly or nonlinearly to obtain the adapted normal growth curve template.

[0060] The preferred value of T' is greater than or equal to 5. When T' is less than 5, the finite time-phase data is insufficient to reliably estimate the mean and standard deviation of the characteristic distribution. In this case, adaptive calibration is not performed, and the alternative template is directly used as the template for the normal growth curve after adaptation. It is marked as a low-confidence template adaptation, and the template confidence level is added to the subsequent stress diagnosis results.

[0061] The difference in characteristic distribution can be measured by the maximum mean difference distance or the correlation alignment distance. The optimization objective is to minimize the distribution difference between the transformed replacement template and the target field data in a limited time period, while maintaining the curve shape and physical rationality of the replacement template throughout the entire growth period.

[0062] The preferred transformation method is piecewise linear scaling plus translation. Specifically, a scaling factor *a* and a translation factor *b* are applied to each dimension of the spectral vegetation index of the replacement template. The scaling factor *a* = standard deviation of the target domain spectral vegetation index / standard deviation of the source domain spectral vegetation index. Here, the target domain standard deviation refers to the statistical standard deviation of the spectral vegetation index values ​​of the target field at each effective observation time within a finite time phase, and the source domain standard deviation refers to the statistical standard deviation of the spectral vegetation index values ​​of the replacement template at each corresponding time within the same finite time phase, with a constraint range of [0.8, 1.2]. The translation factor *b* = mean of the target domain spectral vegetation index - scaling factor × mean of the source domain. Here, the target domain mean refers to the arithmetic mean of the spectral vegetation index values ​​of the target field at each effective observation time within a finite time phase, and the source domain mean refers to the arithmetic mean of the spectral vegetation index values ​​of the replacement template at each corresponding time within the same finite time phase, with a constraint range of [-0.1, +0.1], to ensure that the adapted template does not undergo excessive distortion.

[0063] For example, taking the NDRE index of target field B as an example, the mean NDRE of the substitute template in the first 10 time phases is 0.385, while the mean NDRE of the target field B in the first 10 time phases is 0.412, with a systematic deviation of +0.027. Using the correlation alignment method, the scaling factor a = 0.068 / 0.062 ≈ 1.097 is calculated, where the standard deviation of the target domain (0.068) is the statistical standard deviation of the NDRE values ​​of the target field B in the first 10 time phases, and the standard deviation of the source domain (0.062) is the statistical standard deviation of the NDRE values ​​of the substitute template in the first 10 time phases; the translation factor b = 0.412 - 1.097 × 0.385 ≈ -0.010. The transformation parameters are applied to all 32 time phases of the substitute template, i.e., multiplied, to obtain the adapted template. The mean NDRE of the first 10 time phases of the adapted template becomes 0.412, which matches the actual data of the target field. The same adaptive calibration was performed on the three indices NDVI, CIre, and MCARI to obtain a complete normal growth curve template after adaptation, with a dimension of 32×4, for use by the subsequent temporal difference attention network.

[0064] It should be noted that the above values ​​are for illustrative purposes only and do not constitute a limitation on the present invention.

[0065] In summary, this step provides each target field with a personalized normal growth reference benchmark adapted to its varietal characteristics and regional climatic conditions. This benchmark includes not only the normal dynamic range of spectral vegetation indices for each key growth stage but also the allowable thresholds for the rate of index change between adjacent growth stages. This effectively eliminates the inherent influence of varietal and geographical differences on the spectral curves, allowing subsequent differential comparison networks to focus their analysis on abnormal deviations caused by nutrient stress, rather than normal fluctuations caused by varietal or regional phenological differences.

[0066] S300: Construct a temporal differential attention network, which includes a phenological gating unit and a differential comparison unit.

[0067] like Figure 2 As shown, in this step, a temporal differential attention network is first constructed. The input layer of the temporal differential attention network receives the temporal spectral index matrix of a single pixel in the spatiotemporal spectral data cube, with a dimension of T×C; and the normal growth curve template, with the same dimension of T×C. The two input data are linearly projected onto a d-dimensional embedding space and superimposed with learnable temporal position codes to preserve temporal information.

[0068] For example, taking the target pixel with coordinates (x1, y1) in winter wheat field A as an example, the temporal spectral index matrix of this pixel has a dimension of 32×4, and the normal growth curve template is also 32×4. After linear projection, the two data streams are mapped to 32×64 embedding vector matrices respectively. The temporal position encoding adopts a sine-cosine position encoding scheme, generating a 64-dimensional position vector for each of the 32 time phases, and adding them element by element into the embedding vector.

[0069] Furthermore, a phenological gating unit is added. The phenological gating unit is used to receive the actual curve embedding sequence and the template curve embedding sequence after superimposed position encoding, and generate a phase-dimensional gating weight vector based on the key growth period window information corresponding to the variety-ecoregion in the phenological knowledge graph.

[0070] The critical growth window refers to the growth stage where nutrient stress is most likely to manifest and has the greatest impact on yield, such as the jointing to grain-filling stage of winter wheat. Phases within this window are assigned higher gating weights, while phases outside the window, such as the overwintering period and late maturity, are assigned lower weights. The gating weight vector is multiplied with the embedded sequence phase by phase to achieve selective weighting of the input features.

[0071] The specific implementation of the phenological gating unit is as follows: First, the start and end times of each growth period stored in the phenological knowledge graph are mapped to specific temporal index ranges; then, initial weights are generated according to the importance level of the growth period to which the temporal index belongs, and the importance level is determined by prior knowledge of crop cultivation; finally, the weights are fine-tuned by learnable scaling and translation parameters, and the final gating weights are output after activation by Sigmoid.

[0072] The gating weights are set as follows: the initial weights for the four key growth stages—jointing, booting, heading, and grain-filling—are set to 1.0; the initial weights for the greening and maturity stages are set to 0.5; and the initial weights for the overwintering and emergence stages are set to 0.1. The gating weights can be adjusted using scaling and translation parameters. The initial values ​​for scaling and translation are set to 1.0 and 0.0, respectively, and can be automatically optimized during subsequent training.

[0073] For example, among the 32 time phases mentioned above, based on the growth period records of wheat A in the northern North China Plain in the phenological knowledge graph, phases 8 to 12 correspond to the jointing stage, phases 13 to 16 correspond to the booting stage, phases 17 to 20 correspond to the heading stage, and phases 21 to 25 correspond to the grain-filling stage. The initial gating weights for these time phases are set to 1.0; phases 4 to 7 correspond to the greening stage, with a weight of 0.5; phases 1 to 3 correspond to the overwintering stage, with a weight of 0.1; and phases 26 to 32 correspond to the maturity stage, with a weight of 0.5. The gating weights are calculated using the Sigmoid function, i.e., final gating weight = Sigmoid((initial weight + translation parameter) × scaling parameter). Each time phase vector of the actual curve embedding sequence is multiplied by the corresponding gating coefficient to enhance the characteristics of key growth periods and suppress non-key periods.

[0074] Furthermore, a differential comparison unit is added. The differential comparison unit is used to receive the actual curve embedding sequence and the template curve embedding sequence after phenological gating weighting. First, the difference embedding vector between the two in each time phase is calculated. The formula is: difference embedding vector = actual curve embedding vector of time phase t - template curve embedding vector of time phase t. The difference vector sequence constitutes a difference feature matrix with dimension T×d.

[0075] Subsequently, the difference feature matrix is ​​input into the self-attention encoding layer. This layer consists of a multi-head self-attention mechanism and a stacked feedforward network. By calculating the attention weights between the difference vectors of each time phase, it captures the temporal dependencies and cumulative effects of stress differences. Specifically, the self-attention mechanism enables the network to focus on the process of difference propagation and amplification from the initial time phase to subsequent time phases. For example, slight differences appearing at the jointing stage accumulate continuously during the heading and grain-filling stages, thus integrating scattered time-phase differences into a global stress response pattern.

[0076] After the attention encoding layer, a global average pooling layer and a fully connected layer are connected to compress the variable-length temporal difference features into a fixed-length feature vector. Finally, the stress probability values ​​of each category are output through the Softmax layer, or the continuous stress probability values ​​are output through the Sigmoid layer.

[0077] For example, after phenological gating weighting, the actual curve embedding vector and template curve embedding vector from phase 8 to phase 25, from the jointing stage to the grouting stage, are retained with a relatively large amplitude. The difference is calculated phase by phase to obtain the difference embedding vector sequence. Taking the NDRE dimension as an example, the expected template NDRE value for phase 8 is 0.435, the actual value is 0.418, and the difference is -0.017; the expected template NDRE value for phase 16 is 0.482, the actual value is 0.398, and the difference is -0.084. These gradually expanding negative differences are captured by the self-attention layer, and the attention weights are concentrated around phase 16, indicating a significant stress response in this phase. After multi-head self-attention with 8 heads, the differential feature sequence is encoded into a 32×64 context-aware feature matrix, compressed into a 64-dimensional vector by global average pooling, mapped to a 4-dimensional output by a fully connected layer, and then obtained by Softmax to obtain the stress probability distribution [normal: 0.10, mild: 0.65, moderate: 0.20, severe: 0.05].

[0078] Finally, the output of the temporal difference attention network is the nutrient stress probability value of that pixel. The nutrient stress probability value can be output in the form of a discrete probability distribution, that is, the probability of each stress level. The discrete form is convenient for generating a spatial distribution map of stress levels with clear classification.

[0079] For example, after the pixel undergoes forward inference through a temporal differential attention network, the discrete probability distribution output by the Softmax layer is [normal: 0.10, mild nitrogen deficiency: 0.65, moderate nitrogen deficiency: 0.20, severe nitrogen deficiency: 0.05]. The maximum probability value is 0.65, corresponding to the category of mild nitrogen deficiency. Furthermore, 0.65 satisfies the condition that the probability range for this category is 0.50 to 1.00, thus the pixel is classified as having mild nitrogen deficiency, with a diagnostic confidence level of 0.65.

[0080] The method provided by this invention also includes: Load the physical simulation sample set generated in advance using the radiative transfer model and the non-differentiable physical constraint rules therein; The non-differentiable physical constraint rules are encoded as physical constraint loss terms, which, together with the classification loss and regression loss, constitute a joint loss function used to train the temporal difference attention network.

[0081] This includes pre-generating a physical simulation sample set using a radiative transfer model and extracting non-differentiable physical constraint rules, including: The nutrient stress dynamic parameters measured in field trials are used as the a priori range. These nutrient stress dynamic parameters include at least chlorophyll content, leaf area index, leaf tilt angle, and soil background reflectance. Multiple sets of nutrient stress dynamic parameters are generated by sampling methods. Each set of parameters is input into the radiative transfer model to obtain a corresponding set of simulated values ​​of canopy reflectance spectrum and vegetation index. Multiple sets of samples are generated repeatedly. Among them, the extracted non-differentiable physical constraint rules include at least: the monotonic relationship between vegetation index and chlorophyll content, and the upper limit of the spectral feature shift amplitude.

[0082] In this step, the physical simulation sample set refers to the physiological parameter-reflectivity-vegetation index data pairs generated using the radiative transfer model within a known range of physiological parameters, used to expand the training data. Non-differentiable physical constraint rules are intrinsic laws extracted from this physical model that cannot be automatically learned by a neural network, such as the monotonically decreasing relationship between chlorophyll content and a specific vegetation index.

[0083] Specifically, the dynamic parameters of nutrient stress measured in field experiments are used as the prior range for parameter sampling. These dynamic parameters include at least: chlorophyll content (Cab), leaf area index (LAI), leaf tilt distribution (LAD), and soil background reflectance (ρsoil). The measured variation range and distribution characteristics of each parameter under normal and stress conditions are obtained through field experiments, serving as reasonable boundaries for the input parameters of the radiative transfer model.

[0084] The aforementioned field trials included: setting up experimental sites in representative ecological zones of the target crop varieties; setting up nitrogen gradient treatments covering a sufficient nitrogen supply control and multiple nitrogen deficiency stress gradients for different crop varieties; and simultaneously collecting three types of observation data at key growth stages of the crop's entire growth period, such as the greening stage, jointing stage, booting stage, heading stage, grain-filling stage, and maturity stage of winter wheat: leaf biochemical parameters represented by chlorophyll content (Cab) measured by SPAD chlorophyll meter and spectrophotometry; canopy structure parameters represented by leaf area index (LAI) and leaf tilt angle distribution (LAD) measured by canopy analyzer and protractor; and soil background reflectance (ρsoil) measured in bare soil areas by ground object spectrometer. Through continuous observation over many years, at multiple sites, and for multiple varieties, the measured variation range and distribution characteristics of nutrient stress dynamic parameters at each growth stage and at each nitrogen level were obtained.

[0085] For example, based on field trial data of winter wheat sample A from the North China Plain, the measured ranges of various parameters throughout the entire growth period were statistically obtained: chlorophyll content (Cab) was 15-55 μg / cm², leaf area index (LAI) was 0.5-6.0, average leaf inclination angle (ALIA) was 30°-65°, and soil background reflectance (ρsoil) in the red light band was 0.08-0.25, corresponding to a range from moist dark soil to dry light soil. These ranges cover the complete dynamic interval from severe nitrogen deficiency stress to normal growth with sufficient nitrogen supply.

[0086] Secondly, a large number of parameter combinations are generated in the multidimensional parameter space through sampling methods. Each parameter combination is input into the radiative transfer model to simulate the corresponding canopy reflectance spectrum, and the simulated values ​​of various vegetation indices are calculated accordingly. Latin hypercube sampling can be used to ensure uniform coverage in the parameter space. The generation process is repeated to accumulate multiple sets of samples to form a physical simulation sample set.

[0087] Among them, the radiative transfer model is a physical model describing the transmission process of solar radiation between the vegetation canopy and the atmospheric medium. By inputting leaf biochemical parameters, such as chlorophyll content, dry matter content, and water content, as well as canopy structure parameters, such as leaf area index and leaf tilt angle distribution, and soil background reflectance, it simulates and outputs the directional reflectance spectrum of the canopy top. Commonly used radiative transfer models include the PROSAIL model, the SAIL model, and the PROSPECT model. The PROSAIL model is formed by coupling the leaf optical characteristic model PROSPECT with the canopy biaxial reflectance model SAIL, and it is widely used in the field of agricultural remote sensing for vegetation biochemical parameter inversion and spectral characteristic analysis.

[0088] Latin hypercube sampling is a hierarchical random sampling method, commonly used in existing multidimensional parameter space sampling techniques. Its basic principle is to divide the range of each parameter's value into several non-overlapping intervals with equal probability, then randomly select a sample point within each interval, and finally randomly combine the sampling results of different parameters to generate a multidimensional parameter combination set. Compared to simple Monte Carlo random sampling, Latin hypercube sampling can achieve more uniform coverage of the entire parameter space with a smaller number of samples, effectively avoiding the problem of sample points clustering or sparse in local areas, thus improving the representativeness of parameter combinations and simulation efficiency.

[0089] It should be noted that the methods used above are all methods disclosed in the prior art and are preferred methods of this solution. Those skilled in the art may also use other methods according to the actual situation, and this invention does not limit them.

[0090] For example, within a defined parameter range, 10,000 parameter combinations are generated using the Latin hypercube sampling method. Each combination contains random values ​​for four parameters: Cab, LAI, ALIA, and ρsoil, with the values ​​uniformly distributed within the measured range of each parameter. A parameter combination [Cab=32μg / cm², LAI=3.5, ALIA=45°, ρsoil=0.15] is input into the PROSAIL model to simulate the canopy reflectance spectrum curve from 400 to 2500 nm. The reflectance values ​​for the corresponding band of Sentinel-2 are extracted from the simulated spectrum, and NDVI=0.782, NDRE=0.418, CIre=1.310, and MCARI=0.095 are calculated. The above simulation samples are recorded as [input parameter combination, output reflectance spectrum, output index value]. All 10,000 combinations are processed repeatedly to form a physical simulation sample set.

[0091] Next, non-differentiable physical constraint rules are extracted from the physical simulation sample set. Non-differentiable physical constraint rules refer to physical relationship constraints that objectively exist in the radiation transfer process but cannot be directly expressed by differentiable functions.

[0092] The aforementioned non-differentiable physical constraint rules include at least the following two categories: The first category is the monotonic relationship constraint between vegetation index and chlorophyll content. That is, under the condition that other parameters are fixed, the increase of chlorophyll content (Cab) will inevitably lead to a monotonic increase in red-edge vegetation indices such as NDRE and CIre. If the stress probability trend predicted by the network conflicts with this monotonic relationship, it indicates that the prediction result violates physical laws.

[0093] For example, in the physical simulation sample set, a subset with LAI=3.5, ALI=45°, and ρsoil=0.15 was selected. It was observed that as Cab increased from 15 μg / cm² to 55 μg / cm², NDRE monotonically increased from 0.285 to 0.512, and CIre monotonically increased from 0.720 to 1.850. Based on this, a monotonic constraint rule was extracted: under the same canopy structure conditions, both NDRE and CIre should be monotonically increasing functions of chlorophyll content (Cab).

[0094] The second type: Upper limit constraint on the amplitude of spectral feature shifts, meaning that under given canopy structure conditions, there is a physical upper limit to the amplitude of spectral index changes caused by changes in chlorophyll content. For example, when LAI is fixed, the change in NDRE caused by extreme changes in Cab from complete nitrogen deficiency to complete nitrogen abundance will not exceed a certain threshold. If the differential features predicted by the network exceed the threshold, it indicates that the prediction results are outside the physical feasible region.

[0095] The threshold is determined as follows: Different combinations of canopy structure parameters are iterated through in the physical simulation sample set. Within each layer, the maximum change in vegetation index corresponding to the change of Cab from the minimum to the maximum prior value is extracted, and this maximum change is recorded as the upper limit threshold for the current LAI level. This process is repeated for all discrete LAI levels to establish a lookup table relationship between LAI and the upper limit threshold. In the physical constraint loss term, if the differential feature value output by the network exceeds the upper limit threshold listed in the lookup table for the corresponding LAI level, a penalty is applied to the excess portion.

[0096] For example, based on statistical analysis of 10,000 sets of PROSAIL model simulation samples, under the conditions of LAI=3.5, ALI=45°, and ρsoil=0.15, when Cab changes from 15 μg / cm2 to 55 μg / cm2, NDRE monotonically increases from 0.285 to 0.512, with a maximum change of 0.227; CIre monotonically increases from 0.720 to 1.850, with a maximum change of 1.130. Based on this, an upper limit threshold constraint is set: for the canopy with LAI≈3.5, the NDRE stress deviation per unit reproductive period must not exceed 0.227, and the CIre stress deviation must not exceed 1.130. In subsequent training, if the difference feature of a certain pixel output by the temporal difference attention network exceeds this threshold, the physical constraint loss term penalizes the sample.

[0097] Finally, the extracted non-differentiable physical constraint rules are encoded into physical constraint loss terms, which, together with the classification loss and regression loss, constitute a joint loss function for training the temporal difference attention network.

[0098] Preferably, the classification loss term is used to measure the difference between the nutrient stress probability value output by the network and the true stress level label, and is calculated using the cross-entropy loss function. The true stress level label is derived from the stress level labels corresponding to different nitrogen gradient treatments in the field experiment.

[0099] For example, if the true stress level of a training sample is mild nitrogen deficiency, the label is coded as category 1, corresponding to the nitrogen gradient N2, and the predicted probability distribution output by the network through the Softmax layer is [normal: 0.10, mild: 0.65, moderate: 0.20, severe: 0.05], then the probability of this sample being predicted as mild nitrogen deficiency is 0.65. The cross-entropy loss is calculated as: classification loss term = -ln(0.65) = 0.4308.

[0100] Preferably, the regression loss term is used to measure the numerical error between the stress probability value output by the network and the continuous stress level code, and is calculated using the mean squared error loss function. The stress level coding rules are: normal = 0, mild = 0.33, moderate = 0.67, and severe = 1.0. The predicted continuous stress probability value output by the network is obtained by weighting the continuous codes of each category in the discrete probability distribution.

[0101] For example, if the true stress code for the same training sample is 0.33, which is mild nitrogen deficiency, the discrete probability distribution output by the network through the Softmax layer is [normal: 0.10, mild: 0.65, moderate: 0.20, severe: 0.05]. The predicted continuous value is 0.10×0+0.65×0.33+0.20×0.67+0.05×1.0≈0.399, and the regression loss term is (0.399-0.33)2≈0.0047.

[0102] The physical constraint loss term includes at least two parts: a monotonicity constraint term and an amplitude constraint term. The monotonicity constraint term penalizes predictions from temporal difference attention networks that violate the physical law that vegetation indices should monotonically increase with chlorophyll content. Therefore, for the same pixel's predictions at two adjacent time phases t and t+1, if NDRE decreases from time phase t to time phase t+1, it indicates a decrease in chlorophyll content and a corresponding increase in stress, with the stress probability change value being positive. When the direction of the stress probability change value output by the network is inconsistent with the direction expected by the aforementioned physical law, it is judged as a violation of the monotonic relationship, and a penalty is imposed on the prediction result. The formula for calculating the monotonicity constraint term is: Monotonicity constraint term = max(0, stress probability change value - normalized NDRE change value)² In the above formula, the change in coercion probability = Pt+1-Pt, where Pt and Pt+1 are the coercion probability values ​​output by the temporal differential attention network for the same pixel at time phase t and time phase t+1, respectively. The range of the change in coercion probability is [-1,1], with positive values ​​indicating an increase in coercion probability and negative values ​​indicating a decrease in coercion probability.

[0103] The normalized NDRE change value = |NDRE change value| / the maximum possible range of NDRE change value. The maximum possible range of NDRE change value is obtained statistically from the simulation sample set of the radiative transfer model and stored in a lookup table. The normalized NDRE change value ranges from [0,1], and the larger the value, the more drastic the change in NDRE.

[0104] For example, the data of a target pixel in a winter wheat field at two adjacent time phases are as follows: During the jointing stage, the NDRE value at time phase t is 0.452, and the stress probability Pt output by the temporal differential attention network is 0.60. During the late jointing stage, the NDRE value at time phase t+1 is 0.398, and the stress probability Pt+1 output by the network is 0.35. Calculating the normalized NDRE change: With the current LAI of 3.5, the maximum possible change in NDRE value obtained through a lookup table is 0.227. The normalized NDRE change = |0.398 - 0.452| / 0.227 = 0.054 / 0.227 ≈ 0.238. The stress probability change = 0.35 - 0.60 = -0.25. A negative stress probability change indicates that the network predicts a decrease in the stress probability.

[0105] Analysis of the physical laws governing the above process: NDRE decreases from 0.452 to 0.398, indicating a decrease in chlorophyll content and an intensification of stress, with a positive expected change in stress probability. However, the network output a change in stress probability of -0.25, showing a decrease instead of an increase, with the prediction direction completely opposite to the physical laws. Substituting into the monotonicity constraint formula: monotonicity constraint term = max(0, -0.488)² = 0. Since the change in stress probability is negative, the difference between the change in stress probability and the normalized change in NDRE must also be negative. After being zeroed by the max function, it becomes 0, and no penalty is imposed. Therefore, this indicates that the above formula only imposes a penalty when the direction of the change in stress probability is consistent with the physical laws, but its positive change amplitude exceeds the reasonable range indicated by the normalized change in NDRE. For cases where the direction of the change in stress probability is completely wrong, the above formula does not impose a penalty, and the constraint of wrong direction is corrected by the classification loss and regression loss through label supervision.

[0106] The amplitude constraint term is used to penalize situations where the amplitude of the difference feature output by the differential contrast unit in the temporal differential attention network exceeds the physical upper limit under the current leaf area index (LAI) of that pixel. If the amplitude of the difference feature between the actual curve output by the network and the normal growth curve template exceeds this upper limit, it indicates that the prediction result has exceeded the physical feasible region and is an abnormal prediction that violates the laws of radiative transfer. Therefore, a penalty should be imposed on this prediction result.

[0107] The formula for calculating the amplitude constraint term is: Amplitude constraint term = max(0, |Difference Feature Value| - Maximum Allowable Deviation)². In this formula, the difference feature value is the output of the difference contrast unit in the temporal difference attention network, representing the deviation of the actual temporal spectral curve from the normal growth curve template in a certain spectral vegetation index dimension. It is obtained by mapping the feature vector after self-attention encoding through a fully connected layer. The maximum allowable deviation is the maximum allowable deviation of the vegetation index under the current LAI condition for that pixel, obtained by searching a pre-constructed lookup table. The lookup table is statistically generated from the radiative transfer model simulation sample set and records the maximum change in each vegetation index caused by the change from the minimum to the maximum prior value of chlorophyll at different discrete LAI levels.

[0108] For example, taking the prediction of the same pixel in the late stage of growth at time t+1 as an example, the current LAI of this pixel is 3.5. The NDRE difference feature value output by the differential contrast unit in the temporal differential attention network is 0.310. By looking up the table, the maximum allowable deviation of NDRE under the condition of LAI=3.5 is 0.227. Substituting into the amplitude constraint formula, the amplitude constraint term = max(0,|0.310|-0.227)2 = 0.0069, which means that the difference feature output by the network exceeds the physical upper limit of 0.083, and the square of the excess part is 0.0069.

[0109] If the NDRE difference feature value of the differential contrast unit output for the same pixel in another training iteration is 0.150, the maximum allowable deviation is still 0.227 according to the table. At this time, the difference feature amplitude is within the physical upper limit, the amplitude constraint term = 0, no penalty is applied, indicating that the difference feature amplitude output by the network does not exceed the physical feasible region, and the prediction result conforms to the radiative transfer law.

[0110] The combination of the monotonic constraint term and the amplitude constraint term is the physical constraint loss term, i.e., physical constraint loss term = monotonic constraint term + amplitude constraint term.

[0111] The joint loss function is expressed as: Joint Loss Function = Classification Loss Term + α × Regression Loss Term + β × Physical Constraint Loss Term. In the above expression, α and β are weighting coefficients used to balance the contributions of each loss term. The setting rules are as follows: α ranges from 0.1 to 1.0, adjusted according to the granularity of the stress level classification; the finer the classification, such as 5 levels or higher, the larger the value of α, with 0.5 being the preferred value; β ranges from 0.01 to 0.5, adjusted according to the consistency between the physical simulation sample set and the measured data; the higher the consistency, the larger the value of β, with 0.1 being the preferred value. During training, the physical constraint loss term penalizes predictions that violate monotonic relationships and exceed the amplitude upper limit, guiding the network to learn stress diagnosis logic that conforms to physical laws.

[0112] For example, based on the above values, the losses for the current training samples are summarized and the joint loss function is calculated. Classification loss term = 0.4308, regression loss term = 0.0047, monotonicity constraint term = 0, amplitude constraint term = 0.0069, physical constraint loss term = monotonicity constraint term + amplitude constraint term = 0 + 0.0069 = 0.0069, weight coefficients: α = 0.5, β = 0.1, joint loss function = 0.4308 + 0.5 × 0.0047 + 0.1 × 0.0069 = 0.43384. The above joint loss value updates the weight parameters of each layer of the temporal difference attention network through the backpropagation algorithm, completing one training iteration. In batch training, the average of all samples within the batch is used for gradient calculation.

[0113] Preferably, the validation set performance metrics no longer improve as a criterion for model convergence, i.e., an early stopping mechanism. Specifically, during training, after each training epoch, the nutrient stress level classification accuracy is calculated on an independently partitioned validation set. If the validation set accuracy does not reach a new highest value within M consecutive training epochs, early stopping is triggered, the model is deemed to have converged, and the model parameters at the point of highest validation set accuracy are saved as the final model. The preferred value for M is 15 epochs, which strikes a balance between avoiding underfitting due to premature stopping and overfitting due to late stopping.

[0114] For example, taking the training of a nutrient stress diagnostic model for winter wheat in the North China Plain as an example, 4050 samples obtained from field trials were divided into training, validation, and test sets in an 8:1:1 ratio. During the training process, the validation set accuracy reached its highest value of 91.8% in the 65th round. Subsequently, the validation set accuracy fluctuated between 90.1% and 91.5% in the 66th to 80th rounds (a total of 15 rounds), never exceeding 91.8%, triggering the early stopping mechanism. The model was determined to have converged in the 65th round, training was stopped, and the model parameters for that round were saved.

[0115] It should be noted that the above values ​​are for illustrative purposes only and do not constitute a limitation on the present invention.

[0116] In summary, this step fundamentally eliminates the confusion caused by crop variety, geographical environment, and normal phenological changes in the diagnostic results; the physical constraint loss explicitly introduces the non-differentiable physical laws in the radiative transfer model into the training process, guiding the network to learn stress diagnostic logic that conforms to physical laws, thus avoiding the prediction results that may be produced by a purely data-driven model that violate physical common sense.

[0117] S400: The spatiotemporal spectrum data cube and the normal growth curve template are input into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve and outputs the nutrient stress probability value after self-attention encoding.

[0118] This step organically integrates the data foundation, knowledge reference, and network model built in the previous steps, achieving an end-to-end mapping from raw time-series remote sensing data to pixel-level stress probabilities. Through a three-stage processing of gated weighting, differential calculation, and self-attention encoding, the network can focus on abnormal deviations of the actual spectrum from the normal template during critical growth periods, thereby distinguishing between spectral differences caused by normal phenological changes and those caused by nutrient stress, and outputting a reliable stress probability estimate before stress symptoms become apparent.

[0119] Step S400 in the method provided by the present invention includes: Read the time windows of seedling stage, jointing stage, heading stage, flowering stage, grain filling stage and maturity stage from the phenological knowledge graph, and generate an initial gating vector based on the time windows; Among them, the position located within the window from the jointing stage to the heading stage is assigned a first weight value, and the position located within the window from the seedling stage to the maturity stage is assigned a second weight value, and the first weight value is greater than the second weight value. During the training of the temporal difference attention network, the parameters are updated starting from the initial gate vector, and regularization constraints are applied to limit the deviation of the updated gate vector from the initial gate vector.

[0120] In this step, the time windows for each key growth stage corresponding to the crop variety-ecoregion combination are first retrieved from the phenological knowledge graph loaded in the previous step, based on the crop variety and ecoregion of the target field. The time window defines the start and end phase index range of each growth stage, including at least six key growth stages: seedling stage, jointing stage, heading stage, flowering stage, grain-filling stage, and maturity stage. The time windows for each growth stage are determined by statistical analysis of multi-year field trial observation data, based on the average growth period length and accumulated temperature requirement of the variety in the corresponding ecoregion.

[0121] For example, taking wheat A in the northern ecological zone of the North China Plain as an example, the entire growth period covers 32 effective observation phases. The time windows of each growth period stored in the phenological knowledge graph are as follows: seedling stage corresponds to phases 1 to 4, jointing stage corresponds to phases 5 to 9, heading stage corresponds to phases 10 to 13, flowering stage corresponds to phases 14 to 17, grain filling stage corresponds to phases 18 to 25, and maturity stage corresponds to phases 26 to 32.

[0122] Next, based on the time windows of each reproductive period read in step one, the corresponding initial gating weights are assigned to each phase position of the entire reproductive period, generating an initial gating vector of dimension T×1, where T is the number of effective phases in the entire reproductive period.

[0123] The gating weights are assigned as follows: Phases within the jointing to heading stage window are assigned a first weight value; phases within the seedling and maturity stages window are assigned a second weight value, with the first weight value being greater than the second weight value. Phases within the flowering and grain-filling stages window are assigned a third weight value, which falls between the first and second weight values.

[0124] The specific rules for setting the weight values ​​are as follows: The first weight value is preferably 1.0, corresponding to the critical stage of crop transition from vegetative growth to reproductive growth. This is the period when nutrient stress has the most sensitive impact on yield and should be given the highest attention. The second weight value is preferably 0.1, corresponding to the early and late stages of crop growth. The spectral vegetation index is mainly dominated by seedling coverage and normal senescence processes, and stress signals are easily confused, so the attention should be reduced. The third weight value is preferably 0.7, corresponding to the flowering, pollination, and grain filling stages. Stress still has a significant impact on yield, but the spectral response pattern is more complex than that from jointing to heading, so it should be given a moderate to high attention.

[0125] For example, an initial gating vector is generated based on the time window read in step one. Phases 1 to 4 are located in the seedling stage window and are assigned a second weight value of 0.1; phases 5 to 9 are located in the jointing stage window and are assigned a first weight value of 1.0; phases 10 to 13 are located in the heading stage window and are assigned a first weight value of 1.0; phases 14 to 17 are located in the flowering stage window and are assigned a third weight value of 0.7; phases 18 to 25 are located in the grain-filling stage window and are assigned a third weight value of 0.7; and phases 26 to 32 are located in the maturity stage window and are assigned a second weight value of 0.1. The initial gating vector is then [0.1,0.1,0.1,0.1,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.7,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1].

[0126] Finally, during the training of the temporal difference attention network, starting from the initial gating vector, the gating weights are set as learnable parameters and updated along with other network parameters through backpropagation. To prevent the gating parameters from deviating excessively from the initial values ​​set based on phenological prior knowledge during training, a regularization constraint is applied to the updated gating vector to limit the magnitude of deviation between the updated gating vector and the initial gating vector.

[0127] The specific form of the regularization constraint is to add a gating regularization term to the joint loss function. The calculation formula is: Gating regularization term = Regularization strength coefficient × ||Updated gating vector - Initial gating vector||². The rules for setting the regularization strength coefficient are as follows: the value ranges from 0.01 to 0.5. When there is sufficient observation data of this variety-ecoregion in the phenological knowledge graph, such as a sample size of more than 100 groups, the value is smaller, such as 0.05, allowing for a larger adaptive adjustment space for the gating parameter; when the observation data is limited, such as a sample size of less than 30 groups, the value is larger, such as 0.3, making the gating parameter more dependent on prior phenological knowledge and avoiding overfitting.

[0128] For example, in one training iteration, the current gate vector is updated to 0.85 at phase 6 of the jointing stage and to 0.22 at phase 2 of the seedling stage. The gate regularization term = 0.1 × [(0.85 - 1.0)² + (0.22 - 0.1)²] is calculated by multiplying the sum of squared deviations for all 32 phase positions by the regularization strength coefficient = 0.1 to obtain the regularization penalty value.

[0129] Step S400 in the method provided by the present invention further includes: The difference values ​​at each observation time are obtained by subtracting the actual spectral vegetation index at each observation time in the spatiotemporal spectral data cube from the normal spectral vegetation index at the same time in the normal growth curve template, and then arranging them in chronological order to form a difference vector. The difference vector is accumulated over time to obtain the cumulative difference. The first-order difference and the second-order difference between the time series of the actual spectral vegetation index and the time series of the normal spectral vegetation index are calculated respectively, and the difference series of the corresponding orders are subtracted to obtain the first-order difference and the second-order difference. The difference vector, the cumulative difference, the first-order difference, and the second-order difference are concatenated to form the difference feature output.

[0130] In this step, the actual spectral vegetation index value of each observation time for the target pixel in the spatiotemporal spectral data cube is subtracted from the corresponding index value at the same time in the normal growth curve template, item by item, to obtain the difference value at that time. The difference values ​​of all observation times are arranged in chronological order to form the difference vector for that pixel. The dimension of the difference vector is T×C, where T is the number of effective observation phases throughout the entire growth period, and C is the number of spectral vegetation index species. Specifically, in the normal growth curve template, the expected values ​​for NDRE (NDRE) in phase 8 are 0.435, CIre (CIre) are 1.370, MCARI (MCARI) are 0.089, and NDVI (NDVI) are 0.762.

[0131] For example, taking the NDRE index of a pixel at the jointing stage, phase t=8 as an example, the actual spectral vegetation index value is 0.418, and the expected NDRE value at the same time in the normal growth curve template is 0.435. Then the difference value of NDRE at that time is 0.418-0.435=-0.017. Traversing all 32 phases and 4 indices in the same way, the difference vector of phase 8 is [-0.017,-0.12,-0.008,-0.035]. Arranging the difference vectors of the 32 phases in chronological order forms a 32×4 difference matrix, that is, the difference vector sequence.

[0132] Secondly, the difference vectors obtained above are accumulated along the time dimension to calculate the cumulative difference at each moment. The cumulative difference reflects the persistence and amplification trend of the stress deviation over time, and can distinguish between temporary fluctuations and persistent stress.

[0133] For example, the difference vector for phase 8 in the previous example is [-0.017, -0.12, -0.008, -0.035], and the difference vector for phase 9 is [-0.022, -0.15, -0.012, -0.042]. The cumulative difference for phase 8 is [-0.017, -0.12, -0.008, -0.035]; the cumulative difference for phase 9 is the sum of phases 8 and 9, which is [-0.039, -0.27, -0.020, -0.077]. By phase 32, the cumulative difference in NDRE reaches -0.583, indicating that the actual NDRE curve remains below the normal template throughout the entire reproductive period, with a significant cumulative deviation.

[0134] Next, the first-order and second-order differences are calculated for the actual spectral vegetation index time series and the normal spectral vegetation index time series, respectively. The difference series of the corresponding orders are subtracted one time phase at a time to obtain the first-order difference and the second-order difference.

[0135] The first-order difference reflects the temporal rate of change of the spectral vegetation index, and is calculated using the formula: First-order difference = Current temporal spectral vegetation index - Previous temporal spectral vegetation index. The first-order difference is the difference between the first-order differences of the actual sequence and the template sequence. The second-order difference reflects the acceleration of the rate of change of the spectral vegetation index, i.e., the turning point characteristic of the changing trend, and is calculated using the formula: Second-order difference = Current temporal first-order difference - Previous temporal first-order difference. The second-order difference is the difference between the second-order differences of the actual sequence and the template sequence.

[0136] For example, consider the NDRE index of the aforementioned pixels in phases 15 and 16. The actual NDRE in phase 15 is 0.426, and in phase 16 it is 0.398. The actual first-order difference is 0.398 - 0.426 = -0.028. The template NDRE in phase 15 is 0.465, and in phase 16 it is 0.482. The template first-order difference is 0.482 - 0.465 = +0.017. The first-order difference is -0.028 - (+0.017) = -0.045, indicating that the actual NDRE decreases at a rate of 0.028 per ten-day period, while the normal template increases at a rate of 0.017 per ten-day period. The two change in opposite directions, and the first-order difference is significant.

[0137] Furthermore, the second-order difference is calculated. In phase 14, the actual NDRE is 0.441; in phase 15, the actual NDRE is 0.426; the actual first-order difference is 0.426 - 0.441 = -0.015. In phase 16, the actual first-order difference is -0.028. The actual second-order difference is -0.028 - (-0.015) = -0.013, indicating that the rate of NDRE decrease is accelerating. On the template side, in phase 14, the template NDRE is 0.452; in phase 15, the template NDRE is 0.465; the template first-order difference is +0.013; in phase 16, the template first-order difference is +0.017. The template second-order difference is +0.017 - (+0.013) = +0.004, indicating that the rate of increase of the normal template is slightly faster. The second-order difference = -0.013 - (+0.004) = -0.017, reflecting the structural divergence between the actual sequence and the template sequence in terms of their changing trends. The actual curve is decreasing at an accelerated rate while the template curve is rising steadily.

[0138] Finally, the obtained difference vector, cumulative difference, first-order difference, and second-order difference are concatenated along the feature dimension to form the complete difference feature matrix of the pixel, which serves as the output of the difference contrast unit for subsequent self-attention encoding layers to model temporal dependencies. The concatenated difference feature dimension is T×(C×4), meaning that each temporal region contains four sets of features: original difference, cumulative difference, first-order difference, and second-order difference. Each set of features includes C vegetation indices.

[0139] For example, at time phase 16, the concatenated difference feature vector for the aforementioned pixel is as follows: original difference [-0.084, -0.45, -0.031, -0.098], cumulative difference [-0.312, -1.68, -0.112, -0.395], first-order difference [-0.045, -0.22, -0.018, -0.052], and second-order difference [-0.017, -0.09, -0.008, -0.022]. These four sets of vectors are concatenated into a 16-dimensional feature vector. All 32 time phases are concatenated to form a 32×16 difference feature matrix, which is input into a self-attention encoding layer. A multi-head self-attention mechanism captures the correlation patterns of differences at each order in the temporal dimension. Finally, after global average pooling and a fully connected layer, the nutrient stress probability value of the pixel is output.

[0140] Step S400 in the method provided by the present invention further includes: The weighted features output by the phenology gating unit and the difference features output by the differential comparison unit are concatenated along the channel dimension to obtain the fused features; The fused features are input into a multilayer self-attention encoder, and the output of the self-attention encoder is fed into a fully connected classification head to calculate the nutrient stress probability value.

[0141] In this step, the weighted features output by the phenological gating unit and the difference features output by the differential comparison unit are first concatenated along the channel dimension to obtain the fused features. The weighted features output by the phenological gating unit are the actual curve embedding sequence after phase-by-phase weighting by the gating vector, with a dimension of T×d, where T is the number of effective phases throughout the entire growth period and d is the embedding dimension. The difference features output by the differential comparison unit are a multi-order difference feature matrix between the actual curve and the template curve, with a dimension of T×D, where D is the total dimension of all order difference features. After concatenating the two types of features along the channel dimension, a fused feature matrix with a dimension of T×(d+D) is obtained.

[0142] For example, the weighted feature matrix output by the phenological gating unit has a dimension of 32×64, i.e., 32 time phases and 64-dimensional embedding vectors, while the difference feature matrix output by the differential comparison unit has a dimension of 32×16, i.e., 32 time phases and 16-dimensional difference features. Concatenating the two along the channel dimension yields a fusion feature matrix of dimension 32×80. In this matrix, the first 64 dimensions of the 80-dimensional fusion feature vector for phase 8 are the gated weighted embedding representation, and the last 16 dimensions are the concatenated values ​​of the original difference, cumulative difference, first-order difference, and second-order difference for that phase. Gating weighting amplifies the feature response of the jointing phase, with a gating weight of approximately 1.0, while the difference features provide details of the deviation of the actual spectrum from the template for that phase. The concatenation of the two along the channel dimension allows the subsequent self-attention encoder to simultaneously utilize the gated enhanced temporal representation and the finely quantized difference information for stress pattern learning.

[0143] Next, the obtained fusion feature matrix is ​​input into a multilayer self-attention encoder. This multilayer self-attention encoder consists of stacked Transformer encoders, with each layer containing a multi-head self-attention sublayer and a feedforward network sublayer. Residual connections and layer normalization are applied after each sublayer. The self-attention encoder captures the long-range dependencies of the fusion features in the temporal dimension by calculating attention weights between different temporal positions. This enables the network to associate gated weighted features and differential features from different reproductive stages, identifying dynamic patterns of stress transmission, accumulation, and amplification from the initial temporal phase to subsequent temporal phases.

[0144] The preferred number of layers in a multi-layer self-attention encoder is 2 to 6. The layer number setting rule is as follows: when the number of phases in the entire reproductive period is less than 20, L is 2 to 3 layers to avoid overfitting due to excessive model complexity; when the number of phases is between 20 and 50, L is 4 layers; when the number of phases exceeds 50, L is 5 to 6 layers to fully capture complex dependencies in long time series. The preferred number of heads in the multi-head self-attention encoder is 8, with the dimension of each head being the total embedding dimension divided by the number of heads. The hidden layer dimension of the feedforward network is 4 times the total embedding dimension, and the GELU activation function is used.

[0145] For example, the aforementioned pixel has 32 phases throughout its entire reproductive period, and 4 layers are selected according to the setting rule L. A 32×80 fused feature matrix is ​​input into the first layer of the self-attention encoder. An attention weight matrix with a dimension of 32×32 is calculated using 8-head self-attention. In network inference after training convergence, the self-attention weight matrix shows a high attention weight between phase 16 and phases 8 to 15, indicating that the network has learned the dynamic process of stress differences accumulating from the jointing stage and reaching a significant level in the late spikelet stage. After layer-by-layer abstraction by the 4-layer self-attention encoder, a context-aware feature matrix with a dimension of 32×80 is finally output, where the feature vector of each phase has fused the dependency information of the entire reproductive period.

[0146] Finally, the output of the self-attention encoder is fed into the fully connected classification head to calculate the nutrient stress probability value. The fully connected classification head first performs global average pooling on the temporal feature matrix output by the self-attention encoder, compressing it into a fixed-length feature vector with dimension d+D along the time dimension. Then, it passes through one or two fully connected layers to map the feature vector to the stress level category space. Finally, it outputs the probability distribution of each stress level after Softmax normalization, or outputs the continuous stress probability value after the Sigmoid function.

[0147] The structure setting rules for the fully connected classification head are as follows: if the stress level is divided into discrete categories, such as normal, mild, moderate and severe, then the number of neurons in the output layer is the number of categories, and the activation function is Softmax.

[0148] For example, the 32×80 feature matrix output from the 4-layer self-attention encoder is subjected to global average pooling, and the feature vectors of 32 time phases are averaged along the time dimension to obtain an 80-dimensional fixed-length feature vector. This vector is then fed into a fully connected classification head, first passing through a fully connected layer with 32 neurons and the ReLU activation function, then mapped to a 4-dimensional vector by the output layer. After Softmax normalization, the output stress probability distribution is [normal: 0.10, mild: 0.65, moderate: 0.20, severe: 0.05]. The probability of the highest probability category, mild nitrogen deficiency (0.65), is taken as the stress level diagnosis result for that pixel.

[0149] It should be noted that the above values ​​are for illustrative purposes only and do not constitute a limitation on the present invention.

[0150] In summary, this step guides the network to focus its analysis on stress-sensitive stages such as jointing to grain-filling stages using a phenological gating mechanism, while suppressing interference from non-critical stages such as overwintering and maturity. It comprehensively characterizes the anomalous patterns of the actual spectrum relative to the normal template using multi-level difference features, and captures the transmission and amplification relationships of differences over time using self-attention encoding, integrating scattered temporal phase differences into a global stress response pattern. This achieves an end-to-end mapping from time-series remote sensing data to pixel-level stress probabilities, enabling the network to output reliable stress probability estimates before stress symptoms manifest.

[0151] S500: The nutrient stress probability value is used as the output of the temporal differential attention network. The nutrient stress level of the target field is determined based on the nutrient stress probability value, and a spatial distribution map of the stress level is generated.

[0152] In this step, the discrete probability distribution output by the temporal differential attention network is used as the basis for nutrient stress diagnosis of each pixel. The stress level category corresponding to the maximum probability value of each pixel is taken as the nutrient stress level. The pixel spatial coordinates are used for raster rendering to generate a spatial distribution map of stress level of the target field.

[0153] Specifically, the output layer of the temporal difference attention network, after being normalized by the Softmax function, outputs a probability distribution belonging to each nutrient stress level. The stress levels are divided into four categories: normal, mild nitrogen deficiency, moderate nitrogen deficiency, and severe nitrogen deficiency, corresponding to an output dimension of 4. The sum of the probability values ​​for the four categories is 1. The stress level corresponding to the highest probability value is taken as the diagnostic result and used as the diagnostic confidence level.

[0154] Pixels with confidence levels below a preset threshold are marked as low-confidence areas on the distribution map, prompting users to pay attention or perform supplementary sampling for verification. The confidence threshold is set as follows: the preferred value is 0.5. When the maximum probability value of a pixel is below 0.5, it indicates that the network's diagnostic result for that pixel has high uncertainty and may be affected by factors such as residual cloud shadows, mixed pixels, or boundary effects, requiring special identification.

[0155] The spatial distribution map of stress levels has the same resolution as the input remote sensing image, with different color schemes for different stress levels: normal areas are represented in green, mild nitrogen deficiency areas in yellow-green, moderate nitrogen deficiency areas in orange, severe nitrogen deficiency areas in red, and low-confidence areas filled with gray diagonal lines. The spatial distribution map is output in GeoTIFF format and overlaid on the field boundary vector map, providing operational unit-level prescription basis for precise variable topdressing.

[0156] For example, taking winter wheat field A as an example, all valid pixels in the field are input into a temporal difference attention network one by one, and a discrete probability distribution is output for each pixel. The output probability distribution for the pixel with coordinates (x1, y1) is [normal: 0.10, mild nitrogen deficiency: 0.65, moderate nitrogen deficiency: 0.20, severe nitrogen deficiency: 0.05]. The maximum probability of 0.65 corresponds to the mild nitrogen deficiency category, and the confidence level of 0.65 is greater than the threshold of 0.5, so it is judged as mild nitrogen deficiency and rendered as yellow-green. The output probability distribution for the pixel with coordinates (x2, y2) is [normal: 0.08, mild nitrogen deficiency: 0.15, moderate nitrogen deficiency: 0.32, severe nitrogen deficiency: 0.45]. The maximum probability of 0.45 corresponds to the severe nitrogen deficiency category, but the confidence level of 0.45 is lower than the threshold of 0.5. In the distribution map, it is marked with a red diagonal line as low confidence - suspected severe nitrogen deficiency. After processing all pixels in the same way, a spatial distribution map of stress levels with a spatial resolution of 20 meters is generated. The distribution map is then output to the variable fertilizer application system in GeoTIFF format to guide the application of corresponding amounts of nitrogen fertilizer to different areas.

[0157] It should be noted that the above values ​​are for illustrative purposes only and do not constitute a limitation on the present invention.

[0158] In summary, this step uses the four-category probability distribution after Softmax normalization as the diagnostic basis. By taking the stress level corresponding to the maximum probability value and verifying whether it meets a confidence threshold of not less than 0.5, it ensures that the diagnostic result of each pixel has a clear category classification and a quantifiable confidence level. The final output is a spatial distribution map of stress levels with a spatial resolution of 20 meters. This transforms the temporal spectral difference analysis results abstracted by the network into an intuitive field zoning prescription map that can guide variable fertilization operations, realizing an end-to-end closed loop from remote sensing image sequences to field-level precision nitrogen management decisions.

[0159] In summary, the stress level spatial distribution map output by this invention has the same spatial resolution as the input image, enabling high-precision early identification of nutrient stress before stress symptoms appear. This provides operational unit-level prescription basis for precise variable topdressing, reducing fertilizer application while ensuring crop yield, thus achieving both economic and ecological benefits.

[0160] Example 2, as Figure 3 As shown, this invention provides an early diagnosis system for crop nutrient stress based on time-series remote sensing imagery, the system comprising: The spatiotemporal spectral data cube construction module 11 is used to load multi-temporal remote sensing image sequences covering the entire growth period of the target field, and calculate the temporal spectral vegetation index pixel by pixel after preprocessing to obtain the spatiotemporal spectral data cube.

[0161] The loaded multi-temporal remote sensing image sequence covering the entire growth period of the target field includes: Acquire satellite imagery that includes the red-edge band; Each image undergoes radiometric calibration, atmospheric correction, and cloud mask preprocessing to obtain a preprocessed image. Based on the preprocessed image, a spectral vegetation index is calculated pixel by pixel. The spectral vegetation index includes at least the normalized difference red edge index, the red edge chlorophyll index, the improved chlorophyll uptake ratio index, and the normalized difference vegetation index.

[0162] The normal growth curve template retrieval module 12 is used to load a pre-constructed phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field.

[0163] Among them, the pre-construction of a phenological knowledge graph includes: The data collected from field trials over many years includes measured values ​​of chlorophyll content at key growth stages under different crop varieties and nitrogen gradients, as well as synchronously acquired spectral vegetation indices. Based on the measured chlorophyll content, samples with normal nutrient status were selected. Based on samples with normal nutrient status and prior knowledge of crop cultivation, the normal dynamic range and its rate of change confidence interval of the spectral vegetation index under each crop variety-ecological zone-growth period combination are labeled. The normal dynamic range includes an upper limit and a lower limit. The labeled knowledge is stored in a graph data structure, where the node attributes include the mean and standard deviation of the spectral vegetation index for the corresponding growth period, and the edge attributes include the upper and lower limits of the allowable rate of change of the spectral vegetation index between adjacent growth periods.

[0164] The temporal differential attention network construction module 13 is used to construct a temporal differential attention network, which includes a phenological gating unit and a differential comparison unit.

[0165] This system also includes: Load the physical simulation sample set generated in advance using the radiative transfer model and the non-differentiable physical constraint rules therein; The non-differentiable physical constraint rules are encoded as physical constraint loss terms, which, together with the classification loss and regression loss, constitute a joint loss function used to train the temporal difference attention network.

[0166] This includes pre-generating a physical simulation sample set using a radiative transfer model and extracting non-differentiable physical constraint rules, including: The nutrient stress dynamic parameters measured in field trials are used as the a priori range. These nutrient stress dynamic parameters include at least chlorophyll content, leaf area index, leaf tilt angle, and soil background reflectance. Multiple sets of nutrient stress dynamic parameters are generated by sampling methods. Each set of parameters is input into the radiative transfer model to obtain a corresponding set of simulated values ​​of canopy reflectance spectrum and vegetation index. Multiple sets of samples are generated repeatedly. Among them, the extracted non-differentiable physical constraint rules include at least: the monotonic relationship between vegetation index and chlorophyll content, and the upper limit of the spectral feature shift amplitude.

[0167] The nutrient stress probability calculation module 14 is used to input the spatiotemporal spectrum data cube and the normal growth curve template into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve and outputs the nutrient stress probability value after self-attention encoding.

[0168] The phenological gating unit weights the input features based on key reproductive period time windows in the phenological knowledge graph, including: Read the time windows of seedling stage, jointing stage, heading stage, flowering stage, grain filling stage and maturity stage from the phenological knowledge graph, and generate an initial gating vector based on the time windows; Among them, the position located within the window from the jointing stage to the heading stage is assigned a first weight value, and the position located within the window from the seedling stage to the maturity stage is assigned a second weight value, and the first weight value is greater than the second weight value. During the training of the temporal difference attention network, the parameters are updated starting from the initial gate vector, and regularization constraints are applied to limit the deviation of the updated gate vector from the initial gate vector.

[0169] The differential comparison unit calculates the difference characteristics between the actual curve and the template curve, including: The difference values ​​at each observation time are obtained by subtracting the actual spectral vegetation index at each observation time in the spatiotemporal spectral data cube from the normal spectral vegetation index at the same time in the normal growth curve template, and then arranging them in chronological order to form a difference vector. The difference vector is accumulated over time to obtain the cumulative difference. The first-order difference and the second-order difference between the time series of the actual spectral vegetation index and the time series of the normal spectral vegetation index are calculated respectively, and the difference series of the corresponding orders are subtracted to obtain the first-order difference and the second-order difference. The difference vector, the cumulative difference, the first-order difference, and the second-order difference are concatenated to form the difference feature output.

[0170] The self-attention encoding outputs a nutrient stress probability value, including: The weighted features output by the phenology gating unit and the difference features output by the differential comparison unit are concatenated along the channel dimension to obtain the fused features; The fused features are input into a multilayer self-attention encoder, and the output of the self-attention encoder is fed into a fully connected classification head to calculate the nutrient stress probability value.

[0171] The stress level distribution map generation module 15 is used to use the nutrient stress probability value as the output of the temporal differential attention network, determine the nutrient stress level of the target field based on the nutrient stress probability value, and generate a stress level spatial distribution map.

[0172] In this system, when a perfectly matching normal growth curve template cannot be retrieved from the phenological knowledge graph for the target field, a transfer learning strategy is employed: Search for alternative templates that belong to the same agricultural climate zone as the target field and have the highest genetic similarity to the crop variety; Using the limited-temporal remote sensing image data of the target field from the seedling stage to the present, the alternative template is subjected to unsupervised domain adaptive calibration. By minimizing the feature distribution difference between the alternative template and the acquired limited-temporal remote sensing image data, the adapted normal growth curve template is obtained.

[0173] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0174] This specification and accompanying drawings are merely illustrative examples of the invention and are intended to cover any and all modifications, variations, combinations, or equivalents within the scope of the invention. Clearly, those skilled in the art can make various alterations and modifications to the invention without departing from its scope. Therefore, if such modifications and modifications fall within the scope of the invention and its equivalents, the invention is intended to include these modifications and modifications.

Claims

1. A method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery, characterized in that, include: Load multi-temporal remote sensing image sequences covering the entire growth period of the target field, and calculate the temporal spectral vegetation index pixel by pixel after preprocessing to obtain a spatiotemporal spectral data cube. Load the pre-built phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field; A temporal differential attention network is constructed, wherein the temporal differential attention network includes a phenological gating unit and a differential comparison unit; The spatiotemporal spectrum data cube and the normal growth curve template are input into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve, and outputs the nutrient stress probability value after self-attention encoding. The nutrient stress probability value is used as the output of the temporal differential attention network. The nutrient stress level of the target field is determined based on the nutrient stress probability value, and a spatial distribution map of the stress level is generated.

2. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, Load a multi-temporal remote sensing image sequence covering the entire growth period of the target field, including: Acquire satellite imagery that includes the red-edge band; Each image undergoes radiometric calibration, atmospheric correction, and cloud mask preprocessing to obtain a preprocessed image. Based on the preprocessed image, a spectral vegetation index is calculated pixel by pixel. The spectral vegetation index includes at least the normalized difference red edge index, the red edge chlorophyll index, the improved chlorophyll uptake ratio index, and the normalized difference vegetation index.

3. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, Pre-constructing a phenological knowledge graph, including: The data collected from field trials over many years includes measured values ​​of chlorophyll content at key growth stages under different crop varieties and nitrogen gradients, as well as synchronously acquired spectral vegetation indices. Based on the measured chlorophyll content, samples with normal nutrient status were selected. Based on samples with normal nutrient status and prior knowledge of crop cultivation, the normal dynamic range and its rate of change confidence interval of the spectral vegetation index under each crop variety-ecological zone-growth period combination are labeled. The normal dynamic range includes an upper limit and a lower limit. The labeled knowledge is stored in a graph data structure, where the node attributes include the mean and standard deviation of the spectral vegetation index for the corresponding growth period, and the edge attributes include the upper and lower limits of the allowable rate of change of the spectral vegetation index between adjacent growth periods.

4. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, Also includes: Load the physical simulation sample set generated in advance using the radiative transfer model and the non-differentiable physical constraint rules therein; The non-differentiable physical constraint rules are encoded as physical constraint loss terms, which, together with the classification loss and regression loss, constitute a joint loss function used to train the temporal difference attention network.

5. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 4, characterized in that, A physical simulation sample set is generated in advance using the radiative transfer model, and non-differentiable physical constraint rules are extracted, including: The nutrient stress dynamic parameters measured in field trials are used as the a priori range. These nutrient stress dynamic parameters include at least chlorophyll content, leaf area index, leaf tilt angle, and soil background reflectance. Multiple sets of nutrient stress dynamic parameters are generated by sampling methods. Each set of parameters is input into the radiative transfer model to obtain a corresponding set of simulated values ​​of canopy reflectance spectrum and vegetation index. Multiple sets of samples are generated repeatedly. Among them, the extracted non-differentiable physical constraint rules include at least: the monotonic relationship between vegetation index and chlorophyll content, and the upper limit of the spectral feature shift amplitude.

6. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, The phenological gating unit weights the input features based on key reproductive period time windows in the phenological knowledge graph, including: Read the time windows of seedling stage, jointing stage, heading stage, flowering stage, grain filling stage and maturity stage from the phenological knowledge graph, and generate an initial gating vector based on the time windows; Among them, the position located within the window from the jointing stage to the heading stage is assigned a first weight value, and the position located within the window from the seedling stage to the maturity stage is assigned a second weight value, and the first weight value is greater than the second weight value. During the training of the temporal difference attention network, the parameters are updated starting from the initial gate vector, and regularization constraints are applied to limit the deviation of the updated gate vector from the initial gate vector.

7. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, The differential comparison unit calculates the difference characteristics between the actual curve and the template curve, including: The difference values ​​at each observation time are obtained by subtracting the actual spectral vegetation index at each observation time in the spatiotemporal spectral data cube from the normal spectral vegetation index at the same time in the normal growth curve template, and then arranging them in chronological order to form a difference vector. The difference vector is accumulated over time to obtain the cumulative difference. The first-order difference and the second-order difference between the time series of the actual spectral vegetation index and the time series of the normal spectral vegetation index are calculated respectively, and the difference series of the corresponding orders are subtracted to obtain the first-order difference and the second-order difference. The difference vector, the cumulative difference, the first-order difference, and the second-order difference are concatenated to form the difference feature output.

8. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, The self-attention encoding outputs a nutrient stress probability value, including: The weighted features output by the phenology gating unit and the difference features output by the differential comparison unit are concatenated along the channel dimension to obtain the fused features; The fused features are input into a multilayer self-attention encoder, and the output of the self-attention encoder is fed into a fully connected classification head to calculate the nutrient stress probability value.

9. The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery according to claim 1, characterized in that, When the target field cannot find a perfectly matching normal growth curve template from the phenological knowledge graph, a transfer learning strategy is adopted: Search for alternative templates that belong to the same agricultural climate zone as the target field and have the highest genetic similarity to the crop variety; Using the limited-temporal remote sensing image data of the target field from the seedling stage to the present, the alternative template is subjected to unsupervised domain adaptive calibration. By minimizing the feature distribution difference between the alternative template and the acquired limited-temporal remote sensing image data, the adapted normal growth curve template is obtained.

10. A crop nutrient stress early diagnosis system based on time-series remote sensing imagery, characterized in that, The method for early diagnosis of crop nutrient stress based on time-series remote sensing imagery as described in any one of claims 1-9 includes: The spatiotemporal spectral data cube construction module is used to load multi-temporal remote sensing image sequences covering the entire growth period of the target field, and calculate the temporal spectral vegetation index pixel by pixel after preprocessing to obtain the spatiotemporal spectral data cube. The normal growth curve template retrieval module is used to load a pre-built phenological knowledge graph and retrieve the corresponding normal growth curve template based on the crop variety and geographical location of the target field. A temporal differential attention network construction module is used to construct a temporal differential attention network, wherein the temporal differential attention network includes a phenological gating unit and a differential comparison unit; The nutrient stress probability calculation module is used to input the spatiotemporal spectrum data cube and the normal growth curve template into the temporal differential attention network. The phenological gating unit weights the input features according to the key growth period window in the phenological knowledge graph. The differential comparison unit calculates the difference features between the actual curve and the template curve, and outputs the nutrient stress probability value after self-attention encoding. The stress level distribution map generation module is used to use the nutrient stress probability value as the output of the temporal differential attention network, determine the nutrient stress level of the target field based on the nutrient stress probability value, and generate a stress level spatial distribution map.

Images