Sunflower field weeding agent variable spraying optimization method based on multispectral recognition

Abstract

This invention discloses a method for optimizing variable herbicide spraying in sunflower fields based on multispectral recognition, belonging to the field of intelligent agriculture technology. The method includes: acquiring raw multispectral image data of the sunflower canopy at multiple consecutive time points, generating a time-spectral joint feature trajectory for each spatial pixel; performing dynamic time-spectral feature decoupling on the time-spectral joint feature trajectory of each spatial pixel to identify key spectral inflection points characterizing weed stress response; calculating the cumulative spectral deviation of each spatial pixel based on the spectral response differences within the time interval before and after the key spectral inflection point; within each heterogeneous spraying operation unit, inverting and deducing the herbicide demand potential corresponding to that unit based on the cumulative spectral deviation of all spatial pixels it contains; and generating a spraying amount control command for each heterogeneous spraying operation unit based on the herbicide demand potential, which drives the variable spraying execution mechanism.

Description

Technical Field

[0001] This invention relates to the field of intelligent agriculture technology, specifically to a method for optimizing the variable spraying of herbicides in sunflower fields based on multispectral recognition. Background Technology

[0002] Weed management is a particularly prominent issue in large-scale sunflower cultivation. Traditional extensive, uniform herbicide application not only wastes pesticides and pollutes the environment but also easily leads to insufficient or excessive application in certain areas, affecting control efficacy and crop safety. To achieve precise variable-rate spraying, the industry has attempted to use multispectral remote sensing technology to identify weed distribution in the field and generate spray prescription maps accordingly. Existing solutions typically acquire multispectral images at a single time point, distinguishing between crops and weeds by calculating vegetation indices or performing supervised classification, and then dividing the operational area and allocating spray amounts based on weed density maps. This analysis model based on static spectral snapshots has two main drawbacks. First, the spectral characteristics of a single time phase only reflect the canopy state at the moment of acquisition and cannot capture the dynamic evolution of weed stress effects. The spectral response patterns of sunflowers and weeds exhibit temporal differences at different stress stages. Existing technologies struggle to extract the dynamic characteristics of stress responses from temporal spectral changes, resulting in insufficient sensitivity to early or intermittent stress and a tendency to misinterpret spectral fluctuations caused by non-stress factors as weed stress signals. Secondly, the division of operational units typically relies on spatial clustering based on spectral features or fixed grid partitioning. This clustering process fails to incorporate temporal stress evolution information, making it difficult for the divided regions to accurately reflect the spatial variability and patch heterogeneity of weed stress intensity in the field. This results in a spatial mismatch between spraying dosage allocation and actual needs. To address these issues, effectively decoupling the temporal characteristics representing the dynamic response to weed stress from multispectral data and incorporating the cumulative effects of stress evolution into the spatial division of operational units and pesticide dosage decision-making processes has become a key technical challenge for improving the accuracy of variable-rate spraying in sunflower fields. Summary of the Invention

[0003] The purpose of this invention is to provide a method that can dynamically analyze the temporal evolution characteristics of weed stress and optimize variable spraying operations accordingly, so as to solve the problems that existing static spectral analysis is difficult to identify the dynamic evolution of stress and the mismatch between the division of operation units and the spatial distribution of stress.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a method for optimizing the variable spraying of herbicides in sunflower fields based on multispectral recognition. This method acquires raw multispectral image data of the sunflower canopy at multiple consecutive time points and generates a time-spectral joint feature trajectory for each spatial pixel. As one technical solution of the present invention, the reflectance values ​​of each spatial pixel under multiple preset bands are extracted from the raw multispectral image data to form the initial spectral feature vector of that pixel; the initial spectral feature vectors of the same spatial pixel at multiple consecutive time points are stacked in chronological order to construct the initial time-spectral feature matrix of that pixel; principal component dimensionality reduction is performed on the initial time-spectral feature matrix to extract the top few principal component components whose cumulative contribution rate meets a preset standard; the score sequence of the top few principal component components is used as the final time-spectral joint feature trajectory of that spatial pixel. Preferably, the multiple preset bands include at least the green light band, the red light band, and the near-infrared band.

[0005] This invention performs dynamic temporal-spectral feature decoupling on the temporal-spectral joint feature trajectory of each spatial pixel to identify key spectral inflection points characterizing the weed stress response. As a preferred embodiment of this invention, the identification operation specifically involves: performing first-order difference calculation on the temporal-spectral joint feature trajectory to obtain a sequence of feature value changes over time; locating local extrema from the rate of change sequence, where each extrema is a neighborhood where the rate of change changes from positive to negative or vice versa; using the time point corresponding to the local extrema as a candidate inflection point and extracting the change in feature value amplitude before and after the candidate inflection point; and determining the candidate inflection point as a key spectral inflection point when the change in feature value amplitude exceeds a preset weed stress response threshold. Preferably, this weed stress response threshold is dynamically set based on the statistical distribution differences in the change in feature value amplitude between healthy and weed-stressed sunflower canopies, thereby improving the adaptability and accuracy of stress identification.

[0006] Based on the identification of key spectral inflection points, the cumulative spectral deviation of each spatial pixel is calculated according to the difference in spectral response within the time interval before and after the key spectral inflection point. In a specific embodiment of the invention, the time-spectral joint feature trajectory is divided into a feature subsequence before the inflection point and a feature subsequence after the inflection point, using the key spectral inflection point as the boundary; the mean center of the feature subsequence before the inflection point is calculated as the reference spectral reference vector; the Euclidean distance between the feature vector at each time point in the feature subsequence after the inflection point and the reference spectral reference vector is calculated to generate a distance change sequence; the distance change sequence is integrated over time and accumulated, and the accumulated result is used as the cumulative spectral deviation. Preferably, when integrating the distance change sequence over time, the trapezoidal integral method is used to improve the calculation accuracy of the cumulative amount, making the quantification of the deviation more accurate.

[0007] Furthermore, spatial heterogeneity clustering of all spatial pixels in the field is performed based on the cumulative spectral deviation, dynamically dividing the area into multiple heterogeneous spraying operation units. The clustering process is as follows: the cumulative spectral deviation of all spatial pixels is used as a clustering feature and input into an adaptive spatial constraint clustering model; in this model, a spatial adjacency weight is set between adjacent spatial pixels, which is inversely proportional to the Euclidean distance between pixels; iterative clustering operations are performed with the optimization objective of minimizing the variance of the cumulative spectral deviation within the same cluster and maximizing the variance of the cumulative spectral deviation between different clusters; when the iteration converges or reaches the preset number of iterations, the clustering result is output, and the spatial region corresponding to each cluster is a heterogeneous spraying operation unit. This approach takes into account both spectral response differences and spatial continuity, making the divided operation units more consistent with the actual distribution patterns of weeds in the field.

[0008] Within each heterogeneous spraying operation unit, the herbicide demand potential corresponding to that unit is inverted and deduced based on the cumulative spectral deviation of all spatial pixels contained therein. Specifically, the statistical distribution parameters of the cumulative spectral deviation of all spatial pixels within the heterogeneous spraying operation unit are obtained, including the median and interquartile range. The median is used as the central intensity of weed stress in that unit, and the interquartile range is used as the spatial variability of weed stress in that unit. A preset stress-response inversion mapping relationship is invoked to map both the central intensity of weed stress and the spatial variability of weed stress to an initial demand intensity value. The initial demand intensity value is compared with a preset minimum herbicide activation threshold, and the larger of the two values ​​is taken as the final herbicide demand potential. This method not only considers the average level of weed stress but also integrates the spatial heterogeneity information of stress, making the assessment of demand potential more in line with actual control needs.

[0009] Based on the obtained herbicide demand potential, a spraying rate control command is generated for each heterogeneous spraying operation unit. This command drives the variable spraying actuator. The command generation process includes: inputting the herbicide demand potential of each heterogeneous spraying operation unit into the spraying rate decision model, which outputs a theoretical spraying flow rate value linearly related to the potential value; reading the maximum rated flow rate and minimum stable flow rate of a single nozzle of the variable spraying actuator as flow rate constraint boundaries; limiting the theoretical spraying flow rate value within these flow rate constraint boundaries; outputting the maximum rated flow rate when the theoretical spraying flow rate value exceeds the maximum rated flow rate, and outputting a zero flow rate command when it is lower than the minimum stable flow rate; and generating the duty cycle parameter of the corresponding pulse width modulation signal based on the limited flow rate value. This duty cycle parameter constitutes the spraying rate control command. This control mechanism achieves a precise response of the spraying rate to the degree of weed stress while ensuring the safe and stable operation of the actuator.

[0010] As a further optimization of the present invention, after generating the spraying dosage control command, a feedback calibration operation is also performed to form a closed-loop control. This operation includes: acquiring multispectral image data of the sunflower canopy after the spraying operation; generating the cumulative spectral deviation after spraying for each heterogeneous spraying unit based on the feedback multispectral image data; calculating the attenuation ratio between the cumulative spectral deviation after spraying and the cumulative spectral deviation before spraying for each heterogeneous spraying unit; and when the attenuation ratio is lower than a preset response qualification threshold, incremental compensation adjustment is performed on the herbicide demand potential of the heterogeneous spraying unit. Specifically, the incremental compensation adjustment is performed as follows: calculating the difference between the response qualification threshold and the attenuation ratio, using this difference as a compensation intensity coefficient; multiplying the compensation intensity coefficient by a preset compensation step size factor to generate a compensation increment value; adding the compensation increment value to the current herbicide demand potential to generate a corrected herbicide demand potential; and re-inputting the corrected herbicide demand potential into the spraying dosage decision model to generate a compensated spraying dosage control command for performing a secondary spraying operation. This feedback calibration mechanism can effectively correct insufficient spraying caused by environmental factors or model deviations, ensuring weed control effectiveness.

[0011] This invention also possesses cyclic monitoring and adaptive adjustment capabilities. After generating the spraying dosage control command, the variable spraying operation of the entire field is cyclically monitored. Specifically, at preset time intervals, all steps from generating the time-spectral joint feature trajectory to generating the spraying dosage control command are repeatedly executed. The difference between the spraying dosage control command generated in the current cycle and that generated in the previous cycle is compared. When the difference is lower than a preset convergence threshold, or when the number of cycles reaches a preset upper limit, subsequent spraying control operations are stopped. Otherwise, the spraying dosage control command generated in the current cycle is sent to the variable spraying execution mechanism for a new round of spraying operations. This method achieves continuous tracking of weed occurrence dynamics in sunflower fields and dynamic adjustment of spraying strategies, avoiding over-application while ensuring control efficacy.

[0012] The technical effects and advantages provided by the present invention in the above technical solution are as follows: By performing dynamic time-spectral feature decoupling on the time-spectral joint feature trajectory formed by each spatial pixel at multiple consecutive time points, key spectral inflection points characterizing weed stress response are identified, effectively separating weed stress signals from spectral fluctuations caused by crop phenological changes and environmental noise. This method does not rely on the absolute spectral value or static vegetation index threshold at a single time point, but rather captures the moment when the rate of change in spectral features undergoes a directional inflection, achieving precise location of stress trigger points. Therefore, the spectral information upon which spraying decisions are based is no longer an isolated instantaneous state, but a temporal event with clear physiological stress implications. This allows subsequent calculations of deviation and demand intensity to truly reflect the cumulative effect of weed stress, overcoming the problem of overestimation or underestimation of stress caused by the inability to distinguish spectral change inducements in traditional methods. Based on the difference in spectral response within the time interval before and after the key spectral inflection point, the cumulative spectral deviation is calculated and used as the unique feature for clustering spatial pixels across the entire field. By setting spatial adjacency weights inversely proportional to the actual spatial distance between pixels for spatial heterogeneous clustering, heterogeneous spraying operation units can be dynamically divided. This process directly maps the temporal evolution of weed stress intensity to the spatial division of the work area, ensuring that each work unit not only has an inherent consistency in the degree of spectral deviation but also maintains spatial continuity. Within each unit, the cumulative stress effect is similar, resulting in homogeneous application requirements; however, differences in stress levels between units naturally create heterogeneous boundaries. Based on this, the herbicide demand potential is inverted according to the statistical distribution of the cumulative spectral deviation within each unit, generating spraying dosage control instructions adapted to the actual stress intensity of each unit. This ensures that the pesticide distribution accurately matches the true spatial differentiation pattern of weed stress in the field, fundamentally solving the problems of spatial mismatch and ineffective spraying caused by static zoning and fixed spraying dosage settings. Attached Figure Description

[0013] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0014] Figure 1 This is a flowchart of a method for optimizing the variable spraying of herbicides in sunflower fields based on multispectral recognition. Figure 2 This is a flowchart of the time-spectral joint feature trajectory generation process; Figure 3 This is a flowchart for identifying key spectral inflection points in weed stress response; Figure 4 This is a flowchart for calculating the cumulative spectral deviation; Figure 5This is a flowchart of the heterogeneous spraying operation unit division and herbicide demand potential inversion; Figure 6 This is a flowchart of the process for generating spraying dosage control instructions, feedback calibration, and cyclic monitoring. Figure 7 It is the time series variation curve of the first principal component score of spatial pixels; Figure 8 This is a schematic diagram of the sequence of spatial pixel spectral deviation distance changes after the key spectral inflection point; Figure 9 This is a spatial distribution map of heterogeneous spraying operation units in a sunflower field; Figure 10 This is a schematic diagram showing the mapping relationship between the potential demand for herbicides and the theoretical application rate. Detailed Implementation

[0015] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0016] See Figure 1 This invention provides a method for optimizing variable herbicide spraying in sunflower fields based on multispectral recognition. The method includes acquiring raw multispectral image data of the sunflower canopy at multiple consecutive time points, generating a time-spectral joint feature trajectory for each spatial pixel; performing dynamic time-spectral feature decoupling on the time-spectral joint feature trajectory of each spatial pixel to identify key spectral inflection points characterizing weed stress response; calculating the cumulative spectral deviation of each spatial pixel based on the spectral response differences within the time interval before and after the key spectral inflection point; performing spatial heterogeneity clustering of all spatial pixels in the field based on the cumulative spectral deviation, dynamically dividing the field into multiple heterogeneous spraying operation units; within each heterogeneous spraying operation unit, inverting and deducing the herbicide demand potential corresponding to that unit based on the cumulative spectral deviation of all spatial pixels it contains; and generating a spraying amount control command for each heterogeneous spraying operation unit based on the herbicide demand potential, which drives the variable spraying execution mechanism.

[0017] Example 1:

[0018] In specific implementation, please refer to Figure 2The process of generating the temporal-spectral joint feature trajectory for each spatial pixel further includes the following implementation: From the acquired raw multispectral image data, for each spatial pixel, reflectance values ​​under multiple preset bands are extracted. The multiple preset bands include at least the green band, red band, and near-infrared band. In some embodiments, in addition to the green band, red band, and near-infrared band, at least one of the red-edge band and short-wave infrared band may be added, but the green band, red band, and near-infrared band must be included. The extracted green band reflectance value, red band reflectance value, and near-infrared band reflectance value of a single spatial pixel at a single time point together constitute the initial spectral feature vector of that spatial pixel. For the same spatial pixel, its initial spectral feature vectors at multiple consecutive time points are stacked in chronological order to construct the initial temporal-spectral feature matrix of that spatial pixel. In the specific implementation, the initial spectral feature vector collected at the first time point is used as the first row vector of the initial time-spectral feature matrix, and the initial spectral feature vector collected at the second time point is used as the second row vector of the initial time-spectral feature matrix. The same arrangement is followed until the initial spectral feature vector of the last time point is used as the last row vector of the initial time-spectral feature matrix, thus forming a matrix with the number of rows equal to the total number of consecutive time points and the number of columns equal to the preset number of bands.

[0019] Principal component analysis (PCA) is performed on the constructed initial time-spectral feature matrix to extract the top principal components whose cumulative contribution rate meets a preset standard. In this implementation, the preset standard is set to a cumulative contribution rate greater than or equal to 0.90. Based on the spectral reflectance variation characteristics of sunflower canopy in multiple bands of multispectral data, the top principal components can explain most of the original spectral variance information. Setting the cumulative contribution rate threshold to 0.90 can effectively reduce the feature dimension while retaining sufficient details of spectral dynamic changes. It is understood that the preset standard can be adjusted under other spatial resolutions and temporal sampling intervals, but 0.90 is used as the selected threshold in this embodiment. The specific process of PCA dimensionality reduction includes: first, centering each column of the initial time-spectral feature matrix to obtain a centered initial time-spectral feature matrix; then calculating the covariance matrix of the centered initial time-spectral feature matrix; solving for all eigenvalues ​​of the covariance matrix and the unit eigenvector corresponding to each eigenvalue; arranging all eigenvalues ​​in descending order and calculating the first few eigenvalues... The ratio of the sum of individual eigenvalues ​​to the total sum of all eigenvalues, when this ratio first becomes greater than or equal to 0.90, The value of is the number of principal component components to be retained; before extraction The unit eigenvectors corresponding to each eigenvalue are arranged in columns to form a transformation matrix.

[0020] After principal component dimensionality reduction, the initial temporal-spectral feature matrix of the spatial pixel is mapped to a low-dimensional representation, yielding a score sequence of the first few principal component components. This score sequence serves as the final temporal-spectral joint feature trajectory of the spatial pixel. The formula for calculating the score sequence is as follows:

[0021] in, Indicates the first The initial time-spectral feature matrix of each spatial pixel, the number of rows of which is equal to the total number of consecutive time points. The number of columns is the preset number of bands. ; Indicates that by the above The transformation matrix, composed of the unit eigenvectors corresponding to the eigenvalues, arranged column-wise, has the following number of rows: The number of columns is ; This represents the score matrix obtained after mapping, with rows. The number of columns is , Each line corresponds to a point in time. Principal component score vectors, arranged in chronological order. The sequence of principal component score vectors is the temporal-spectral joint feature trajectory of the spatial pixels.

[0022] In practical implementation, the total number of consecutive time points The number of imaging operations is determined by the time series length of the multispectral images acquired, for example, continuous multispectral imaging of the sunflower canopy at fixed time intervals. The value, A positive integer. Preset number of bands. Not less than 3. Transformation matrix It is calculated independently based on the initial time-spectral feature matrix of each spatial pixel. The transformation matrix corresponding to different spatial pixels can be different, thereby realizing adaptive dimensionality reduction for each spatial pixel and completely preserving the morphological features of its respective spectral trajectory.

[0023] See Figure 7 In the figure, the horizontal axis represents consecutive time points, and the vertical axis represents the first principal component score value of the spatial pixel. The legend identifies three curves for spatial pixel 1, spatial pixel 2, and spatial pixel 3. According to the content of Example 1, the curves in the figure correspond to the first principal component score sequence of three different spatial pixels at multiple consecutive time points, reflecting the dynamic change characteristics of the principal components in the time-spectral joint feature trajectory of each spatial pixel.

[0024] Specifically, the first principal component score of spatial pixel 1 showed a gradual upward trend overall, and remained at a relatively stable high value after time point 150, indicating that the spectral characteristics of this spatial pixel underwent significant and continuous changes during the monitoring period, possibly corresponding to an area where the sunflower canopy was under strong weed stress. The score sequence of spatial pixel 3 showed obvious fluctuations, with a larger increase than that of spatial pixel 1. The score rapidly climbed to its highest point in the time point range of 170 to 240, showing significant nonlinear dynamic changes in the spectral response of this spatial pixel, suggesting that the weed stress response in this area was relatively severe. The score of spatial pixel 2 fluctuated around zero overall, without showing a significant trend change, indicating that the spectral characteristics of the location corresponding to this spatial pixel were relatively stable during the monitoring period, and did not show obvious weed stress characteristics.

[0025] The trend and fluctuation characteristics of the above curves are consistent with the description of the time-spectral joint feature trajectory obtained by dimensionality reduction through principal component analysis in Example 1. The score sequence of multiple time points shown in the figure is the basic data for the identification of key spectral inflection points and the calculation of spectral deviation accumulation in subsequent examples. It reflects the dynamic process of the spectral characteristics of different spatial pixels changing over time, and provides key time-series spectral information for realizing spatial heterogeneous clustering and variable herbicide spraying.

[0026] Example 2:

[0027] In specific implementation, please refer to Figure 3 The process of identifying key spectral inflection points characterizing weed stress responses further includes the following implementation: First-order difference calculation is performed on the time-spectral joint feature trajectory to obtain a sequence of eigenvalue changes over time. In this embodiment, the time-spectral joint feature trajectory is represented as a score sequence of the first few principal component components. For a single spatial pixel, the time-spectral joint feature trajectory contains k-dimensional principal component score vectors at T time points. Among these, the first principal component score sequence retains the most original spectral variance information, and is used as the object for first-order difference calculation. Let the first principal component score sequence be represented as... ,in Indicates a time index. , Indicates the first The first principal component score at each time point. The process of performing first-order difference calculations is defined by the following expression:

[0028] in, Indicates the first The rate of change of the characteristic value at each time point The range of values ​​is ; Indicates the first The first principal component score at each time point; Indicates the first The first principal component score at each time point. Through the above point-by-point differencing operation, a length of [length missing] is obtained. The sequence of eigenvalue change rates, where each element corresponds to the change in eigenvalue between two adjacent time points.

[0029] In practical implementation, local extrema are located from the rate of change sequence. A local extrema is defined as the location where the rate of change changes from positive to negative or vice versa. This is done by iterating through each element in the eigenvalue rate of change sequence, for the first element... rate of change ,when and At that time, the first Each time point is determined as a local extremum where the rate of change changes from positive to negative; when and At that time, the first Each time point is identified as a local extremum where the rate of change turns from negative to positive. All time points that satisfy the above conditions are marked to form a candidate inflection point set.

[0030] In practice, the time point corresponding to each local extreme point is taken as a candidate inflection point, and the change in the amplitude of the eigenvalues ​​before and after the candidate inflection point is extracted. The change in the amplitude of the eigenvalues ​​before and after the candidate inflection point is extracted in the following way: taking the candidate inflection point as the benchmark, backtracking to before the candidate inflection point. Before calculating candidate inflection points at each time point. The mean of the first principal component scores at each time point is denoted as . Extending beyond the candidate inflection point After calculating candidate inflection points at each time point. The mean of the first principal component scores at each time point is denoted as . The change in the amplitude of the eigenvalues ​​before and after the candidate inflection point is calculated and expressed as:

[0031] in, This represents the magnitude change of the eigenvalue before and after the candidate inflection point; Indicates after the candidate inflection point The mean of the first principal component scores at each time point; Indicates before the candidate inflection point The mean of the first principal component scores at each time point; This indicates the half-width of the average time window. The value is set to 3. The reason for setting the value to 3 is that the change process of the spectral response of sunflower canopy under weed stress usually appears on a medium time scale. The window width of 3 time points can smooth out short-term random noise interference, while retaining the trend change characteristics of stress response, and will not drown out effective inflection point information due to excessively wide windows.

[0032] In practical implementation, the change in the amplitude of the eigenvalues ​​before and after the candidate inflection point will be... Compare with the preset weed stress response threshold. When When the preset weed stress response threshold is exceeded, the candidate inflection point is identified as a key spectral inflection point. The preset weed stress response threshold is dynamically set based on the statistical distribution differences in the amplitude changes of eigenvalues ​​between healthy and weed-stressed sunflower canopies. The specific process of dynamic setting includes: pre-collecting and labeling multispectral image data of healthy and weed-stressed sunflower canopy samples at different growth stages; calculating the amplitude changes of eigenvalues ​​for all spatial pixels in the healthy canopy sample set to form a healthy amplitude change set; calculating the amplitude changes of eigenvalues ​​for all spatial pixels in the weed-stressed canopy sample set to form a stress amplitude change set; and calculating the mean of the healthy amplitude change set. and standard deviation ; Calculate the mean of the set of stress amplitude changes and standard deviation Preset weed stress response threshold The calculation formula is dynamically set according to Fisher's linear discriminant criterion:

[0033] in, This represents the preset threshold for weed stress response; This represents the mean of the set of changes in health status. The standard deviation of the set of health range variations; This represents the mean of the set of changes in the magnitude of stress. This represents the standard deviation of the set of stress amplitude changes. Using the above method, the preset weed stress response threshold is adaptively adjusted according to the sunflower growth stage and the actual spectral response distribution characteristics in the field, setting different discrimination thresholds at different growth stages.

[0034] In practical implementation, for each candidate inflection point identified in the temporal-spectral joint feature trajectory of a spatial pixel, the aforementioned feature value amplitude change extraction and threshold comparison operations are performed to ensure that the above-mentioned feature value amplitude change is satisfied. Candidate inflection points are identified as key spectral inflection points characterizing the weed stress response. A single spatial pixel may contain zero, one, or multiple key spectral inflection points in its temporal-spectral joint feature trajectory.

[0035] Example 3:

[0036] In specific implementation, please refer to Figure 4 The process of calculating the cumulative spectral deviation of each spatial pixel further includes the following implementation: Using key spectral inflection points as boundaries, the time-spectral joint feature trajectory is divided into a feature subsequence before the inflection point and a feature subsequence after the inflection point. When multiple key spectral inflection points exist in the time-spectral joint feature trajectory of a spatial pixel, the key spectral inflection point with the earliest time index is selected as the dividing criterion. Let the score sequence of the first principal component in the time-spectral joint feature trajectory be... Several time points were identified as key spectral inflection points, and the time index was adjusted from... arrive The subsequence formed by the first principal component score is used as the feature subsequence before the inflection point, and the time index is from arrive The subsequence formed by the scores of the first principal components is used as the feature subsequence after the inflection point. If there are zero key spectral inflection points in the time-spectral joint feature trajectory, the cumulative spectral deviation of the spatial pixel is directly assigned to zero, and no subsequent calculation steps are performed.

[0037] In practice, the mean center of the feature subsequence before the inflection point is calculated, and the calculated mean center is used as the reference spectral baseline vector. The feature subsequence before the inflection point contains... The scores of the first principal components, represented by the reference spectral baseline vector, are as follows:

[0038] in, This represents the reference spectral baseline vector, whose value is equal to the arithmetic mean of the scores of all first principal components in the feature subsequence before the inflection point. This represents the time index value corresponding to the key spectral inflection point. greater than 1 and less than Integers; Indicates time index The corresponding first principal component score; For time index variables, The value is from arrive Integers.

[0039] It is understandable that the feature subsequence before the inflection point corresponds to the relatively stable state of the sunflower canopy before it is affected by weed stress. Using the mean center of this subsequence to represent the baseline spectral characteristics of the canopy when it is not under significant stress can reflect the inherent spectral attributes of the spatial pixel itself.

[0040] In practice, the Euclidean distance between the feature vector at each time point in the feature subsequence after the inflection point and the reference spectral baseline vector is calculated to generate a distance change sequence. The feature subsequence after the inflection point contains... For each time point in the feature subsequence after the inflection point, calculate the first principal component score and the reference spectral baseline vector at that time point. The Euclidean distance between them. Since the score of the first principal component is a one-dimensional scalar, the Euclidean distance simplifies to the absolute value of the difference between the two. The distance change sequence of the [missing information]... The calculation method for each element is expressed as follows:

[0041] in, Indicates the distance change sequence of the th The Euclidean distance values ​​corresponding to each element; For time indexing after the inflection point, The value is from arrive Integers; Indicates the time index is The first principal component score; This represents the reference spectral reference vector. Through the above point-by-point calculations, a length of [length missing] is obtained. The distance variation sequence, where each element reflects the degree of deviation of the canopy spectral characteristics from the baseline state at the corresponding time point.

[0042] In practice, the distance variation sequence is integrated over time, and the accumulated result is used as the cumulative spectral deviation. The trapezoidal integration method is used to improve the accuracy of the cumulative calculation. The trapezoidal integration method approximates the interval formed by each pair of adjacent elements in the distance variation sequence as a trapezoid, calculates the area of ​​the trapezoid, and sums the areas of all trapezoids. The formula for calculating the cumulative spectral deviation is expressed as:

[0043] in, Indicates the first Cumulative spectral deviation of each spatial pixel; Indicates the distance change sequence of the th The Euclidean distance values ​​corresponding to each element; Indicates the distance change sequence of the th The Euclidean distance values ​​corresponding to each element; This indicates the time sampling interval between adjacent time points. The value is determined by the time sampling frequency of the multispectral image acquisition device. When a fixed time interval is used for acquisition, It is a constant positive real number; The range of summation is from arrive The trapezoidal integral method constructs a trapezoidal region using two adjacent distance values. Compared to the rectangular integral method, it can more accurately approximate the true area under the distance curve, improving the numerical accuracy of the cumulative spectral deviation when the length of the distance variation sequence is limited.

[0044] See Figure 8 The horizontal axis in the figure represents the time point number after the key spectral inflection point. The value ranges from 0 to 190, and the vertical axis represents the Euclidean distance at the corresponding time point. The values ​​range from approximately 0 to 3.5. The legend identifies the distance variation sequence curves of the three spatial pixels, with "pixel 1" represented by a solid line, "pixel 2" by a dashed line, and "pixel 3" by a dotted line.

[0045] According to the description of Example 3, the distance change sequence This reflects the degree of deviation of the first principal component score at each time point after the inflection point from the mean center before the inflection point. The three curves in the figure as a whole show the changes over time. The gradually increasing trend indicates that the degree to which the canopy spectral characteristics of the three spatial pixels deviate from the baseline state gradually intensifies over time, reflecting the cumulative development of weed stress effects over time.

[0046] Specifically, the distance change curve of pixel 2 is the highest overall, with the maximum Euclidean distance exceeding 3.0, indicating that the spectral characteristics of this pixel deviate the most after the critical spectral inflection point, and the stress level is relatively severe; the distance change curve of pixel 1 is in the middle, with the maximum distance value being about 2.3; the curve of pixel 3 is the lowest, with the maximum distance being about 1.5, indicating that its spectral deviation is the smallest and the stress response is relatively weak.

[0047] In addition, the three curves at an early time point (approximately to All showed a fluctuating upward trend, but the increase was limited, and then from to Within the interval, the slope of the curve increases significantly, indicating an enhanced rate of spectral deviation accumulation. Combined with the trapezoidal integral method used in Example 3 to integrate the distance change sequence over time, the growth trend of the curve shown in this figure directly affects the magnitude of the spectral deviation accumulation, thus having a decisive impact on subsequent spatial heterogeneity clustering and herbicide demand potential projection.

[0048] Furthermore, the area between the distance curves of pixel 1 and pixel 2 is shaded in the figure to highlight the difference in their spectral deviations, reflecting the heterogeneity between different spatial pixels. This spatial heterogeneity provides a data foundation for the adaptive spatial constraint clustering based on the cumulative amount of spectral deviation in the subsequent embodiment 4.

[0049] In summary, this figure fully demonstrates that after calculating the key spectral inflection point in Example 3, the spectral deviation characteristics of different spatial pixel time series are revealed by the accumulation of time integrals of distance change sequences, which reveals the heterogeneity of canopy spectral response over time. This provides key input data for dynamically dividing heterogeneous spraying operation units and accurately inverting the potential demand for herbicides.

[0050] Example 4:

[0051] In practical implementation, the process of dynamically dividing multiple heterogeneous spraying operation units further includes the following implementation methods: The cumulative spectral deviation of all spatial pixels is used as a clustering feature and input into an adaptive spatial constraint clustering model. The core framework of the adaptive spatial constraint clustering model is an iterative optimization framework that combines spatial constraints and discriminative objectives; see [reference needed]. Figure 5 The input consists of the cumulative spectral deviation of each spatial pixel within the entire field and the two-dimensional spatial coordinates of each pixel. The output is the label of the heterogeneous spraying operation unit to which each spatial pixel belongs. The optimization process of the adaptive spatially constrained clustering model is driven by the following objective function:

[0052] in, This represents the total energy function value for cluster optimization. The adaptive spatially constrained clustering model minimizes... To determine the clustering results; This indicates the preset number of heterogeneous spraying operation units. The value is preset based on the field area and the minimum operating width of the variable spraying actuator. It is an integer greater than or equal to 2; For clustering index variables, The value ranges from 1 to Integers; Indicates being divided into the first A set of spatial pixels in a cluster; For spatial cell index variables, Representing spatial pixels Belongs to clustering ; Indicates the first Cumulative spectral deviation of each spatial pixel; Indicates the first The mean of the spectral deviations of each cluster center. ; Indicates the first The number of spatial pixels contained in each cluster; This represents the global mean of the cumulative spectral deviation of all spatial pixels in the entire field. This represents the inter-class separation adjustment coefficient. The range of values ​​is an interval. The specific values ​​were selected by using five-fold cross-validation on historical multispectral data to maximize the silhouette coefficient. This represents the spatial smoothing penalty coefficient. The range of values ​​is an interval. The specific values ​​were selected from historical multispectral data through five-fold cross-validation to achieve the highest consistency between the clustering results and the actual weed distribution boundaries in the field. This represents a set of spatially adjacent pairs of pixels, with the adjacency relationship defined using four-neighbor or eight-neighbor domains. Represents a set A pair of adjacent spatial pixels in the image; Representing spatial pixels With spatial pixels Spatial adjacency weights between them; As an indicator function, when spatial pixels Clustering tags With spatial pixels Clustering tags When they are not the same, When the cluster labels are the same, .

[0053] In practical implementation, spatial adjacency weight Calculated based on the Euclidean distance between adjacent spatial pixels. It is inversely proportional to the Euclidean distance between pixels, specifically set as follows: ,in, Representing spatial pixels With spatial pixels The Euclidean distance between two-dimensional spatial coordinates is used. Adding 1 to the denominator avoids the denominator being meaningless when the distance is zero, and also ensures that the weight is close to the maximum value of 1 when adjacent spatial pixels are closely adjacent. Through this weight setting, the closer the spatial pixel pair is, the greater the penalty it receives when the cluster labels are inconsistent, thus prompting spatially continuous regions with similar spectral deviation characteristics to be classified into the same heterogeneous spraying operation unit.

[0054] In practical implementation, the specific iterative computation steps of the adaptive spatial constraint clustering model are as follows. The first step is to calculate the cumulative spectral deviation of all spatial pixels within the entire field. Perform K-means clustering initialization to obtain initial cluster labels and cluster centers. The second step is to fix the cluster labels and cluster centers, and then calculate the total energy function under the current clustering configuration. The third step is to iterate through each spatial cell, attempting to change the spatial cell from its current cluster label to another cluster label, and then calculate the change in the energy function after the change. ;like If so, then accept the label change and update the corresponding cluster centers; if If the label remains unchanged, the algorithm continues with the third step. The fourth step iterates through all spatial pixels until no label changes occur, or the preset number of iterations is reached. The preset number of iterations is set to 100. This value is based on the fact that the energy function value typically stabilizes after dozens of iterations, and 100 iterations are sufficient to ensure that the algorithm converges under most field conditions.

[0055] In practice, when the iteration converges or reaches the preset number of iterations, the clustering results are output. Each spatial cell is assigned a unique clustering label. Spatial cells with the same clustering label form a connected operational region in space, which is a heterogeneous spraying operational unit. Each heterogeneous spraying operational unit is separable in terms of both cumulative spectral deviation and spatial location, thus reflecting the spatial heterogeneity distribution formed by different weed stress levels within the sunflower field.

[0056] In practical implementation, the process of inverting and deducing the herbicide demand potential corresponding to the unit further includes the following implementation methods. For each heterogeneous spraying operation unit, the statistical distribution parameters of the cumulative spectral deviation of all spatial pixels within the heterogeneous spraying operation unit are obtained. The statistical distribution parameters include the median and interquartile range. The median is used as the weed stress center intensity of the heterogeneous spraying operation unit, which reflects the typicality of weed stress within the unit; the interquartile range is used as the spatial variability of weed stress within the heterogeneous spraying operation unit, which reflects the unevenness of stress levels within the unit.

[0057] In practical implementation, a pre-defined stress-response inversion mapping relationship is invoked to map the central intensity of weed stress and the spatial variability of weed stress together to the initial demand intensity value. The pre-defined stress-response inversion mapping relationship is constructed as follows: During historical sunflower planting seasons, multispectral image data of sunflower canopies under different weed coverage and distribution uniformity are simultaneously collected through field trials. The median and interquartile range of the cumulative spectral deviation for each experimental plot are calculated. At the same time, the recommended herbicide dosage required for optimal weed control in each experimental plot is determined through manual calibration. Using the median and interquartile range of the cumulative spectral deviation as two-dimensional input variables and the recommended herbicide dosage as the output target, a two-input single-output support vector regression model is trained. The kernel function of the support vector regression model adopts a radial basis function, and the scaling parameter and regularization constant of the kernel function are determined through grid search and ten-fold cross-validation. The trained support vector regression model constitutes a pre-defined stress-response inversion mapping relationship. In practical applications, this mapping relationship directly receives the central intensity of weed stress and the degree of spatial variation of weed stress, and outputs the initial demand intensity value.

[0058] It can be understood that the initial demand intensity value represents the theoretically appropriate relative dosage of herbicide to be applied under given weed stress center intensity and spatial variability. After obtaining the initial demand intensity value, it is compared with the preset minimum herbicide activation threshold. The preset minimum herbicide activation threshold is pre-calibrated based on the multispectral response characteristics corresponding to the economic damage threshold of weeds in sunflower fields. Specifically, the calibration method is as follows: in a control area with extremely low weed density and no need for chemical control, the spectral deviation of spatial pixels from the cumulative distribution range is measured, and the 90th percentile value of this distribution range is taken as the minimum herbicide activation threshold. When the initial demand intensity value is lower than the minimum activation threshold for herbicides, it indicates that the weed stress level of the heterogeneous spraying operation unit has not reached the level that requires spraying. In this case, the larger value between the preset minimum activation threshold for herbicides and the initial demand intensity value is taken as the final herbicide demand potential, thereby avoiding ineffective spraying in areas with excessively low weed pressure. When the initial demand intensity value is greater than or equal to the minimum activation threshold for herbicides, the initial demand intensity value is the final herbicide demand potential.

[0059] See Figure 9 The figure shows the spatial distribution of heterogeneous spraying operation units based on spatial coordinates, with the horizontal axis representing the X coordinate (in meters) and the vertical axis representing the Y coordinate (in meters). Three heterogeneous spraying operation units are identified by different symbols: heterogeneous spraying operation unit 1 (circled), heterogeneous spraying operation unit 2 (squared), and heterogeneous spraying operation unit 3 (triangled). Each symbol represents the geographical location of a spatial pixel, reflecting the corresponding heterogeneous spraying operation unit to which that pixel belongs.

[0060] The spatial boundary areas of the three heterogeneous spraying operation units are marked by dashed rectangles in the figure. The size and position of the dashed rectangles correspond to the spatial range covered by each operation unit. Heterogeneous spraying operation unit 1 is mainly distributed in the spatial coordinate range of approximately 10 to 40 meters X and approximately 20 to 50 meters Y, showing a relatively concentrated and spatially continuous distribution. Heterogeneous spraying operation unit 3 is mainly distributed in the region of approximately 30 to 70 meters X and approximately 30 to 70 meters Y, exhibiting characteristics of an intermediate transition zone. Heterogeneous spraying operation unit 2 covers a spatial range of approximately 50 to 90 meters X and approximately 55 to 85 meters Y, with its location slightly to the upper right.

[0061] The three heterogeneous spraying operation units in the figure are spatially adjacent but have clear boundaries, and the spatial pixels within each unit are relatively clustered. This demonstrates that the adaptive spatially constrained clustering model based on cumulative spectral deviation and spatial adjacency weights has effectively divided the field into heterogeneous segments. This segmentation result reflects the spatial heterogeneity of weed stress within the sunflower field, which is beneficial for implementing differentiated herbicide variable spraying strategies for each heterogeneous spraying operation unit.

[0062] Example 5:

[0063] In specific implementation, please refer to Figure 6 The process of generating spraying rate control instructions for each heterogeneous spraying operation unit further includes the following implementation: The herbicide demand potential of each heterogeneous spraying operation unit is input into the spraying rate decision model, and the spraying rate decision model outputs a theoretical spraying flow rate value linearly correlated with the herbicide demand potential. The core framework of the spraying rate decision model is a linear mapping function, which is predetermined through field calibration experiments. Specifically, multiple sets of known herbicide demand potentials are set for different heterogeneous spraying operation units. For each set of herbicide demand potentials, the optimal spraying flow rate value is found by gradually adjusting the spraying flow rate and evaluating the weed control efficacy and sunflower phytotoxicity. A univariate linear regression is performed on the herbicide demand potentials of all calibration samples and the optimal spraying flow rate value to obtain the linear correlation coefficient. The coefficient of determination of the regression equation must meet the requirement of being greater than or equal to 0.90 to ensure the reliability of the model. The input of the spraying rate decision model is the herbicide demand potential, and the output is the theoretical spraying flow rate value.

[0064] In practice, the maximum rated flow rate and minimum stable flow rate of a single nozzle of the variable spraying actuator are read and used as flow constraint boundaries. The maximum rated flow rate and minimum stable flow rate of a single nozzle are obtained from the product technical manual of the variable spraying actuator, reflecting the maximum flow rate that a single nozzle can output under rated operating conditions and the minimum flow rate that can maintain stable atomization performance, respectively. The theoretical spraying flow rate value is limited within the flow constraint boundaries. The specific rules for limiting the flow rate are as follows: when the theoretical spraying flow rate value exceeds the maximum rated flow rate, the maximum rated flow rate is output as the limited flow rate value; when the theoretical spraying flow rate value is lower than the minimum stable flow rate, a zero flow command is output; when the theoretical spraying flow rate value is between the minimum stable flow rate and the maximum rated flow rate, the theoretical spraying flow rate value is output as the limited flow rate value. Outputting a zero flow command indicates that the weed stress level of this heterogeneous spraying operation unit is too low, and herbicide spraying is not required.

[0065] In practical implementation, based on the flow rate value after limiting, a duty cycle parameter for the corresponding pulse width modulation signal is generated. This duty cycle parameter constitutes the spray volume control command. The duty cycle parameter is calculated as follows: assuming the pulse duty cycle corresponding to the continuously fully open state of the variable spray actuator solenoid valve is 100%, and the pulse duty cycle corresponding to the continuously fully closed state is zero; when the flow rate value after limiting is zero, the duty cycle parameter is directly set to zero; when the flow rate value after limiting is not zero, the duty cycle parameter is expressed as the ratio of the limited flow rate value to the maximum rated flow rate of a single nozzle multiplied by 100%. This duty cycle parameter is encoded into a digital control signal recognizable by the variable spray actuator and sent to the controller of the variable spray actuator.

[0066] In practice, after generating the spraying dosage control command, a feedback calibration operation is performed. Multispectral image data of the sunflower canopy is collected after the spraying operation. The data collection time is set to a preset number of days after the spraying operation is completed. This preset number of days is determined based on the herbicide's onset rate in the sunflower field, generally 3 to 7 days after spraying. Based on the feedback multispectral image data, the cumulative spectral deviation after spraying for each heterogeneous spraying unit is generated using the same method as generating the cumulative spectral deviation before spraying. The attenuation ratio of the cumulative spectral deviation after spraying for each heterogeneous spraying unit to the cumulative spectral deviation before spraying is calculated. The formula for the attenuation ratio is expressed as:

[0067] in, Indicates the attenuation ratio. The range of is the interval of real numbers; This indicates the cumulative spectral deviation before the spraying operation was performed. Determined by the median of the cumulative spectral deviation of all spatial pixels within the heterogeneous spraying operation unit; This indicates the cumulative deviation of the spectrum after spraying. Similarly, the spectral deviation was determined by the median of the cumulative spectral deviation of all spatial pixels within the heterogeneous spraying unit after spraying. When the herbicide effectively suppresses weed stress, the spectral deviation of the sunflower canopy tends to be alleviated, and the cumulative spectral deviation after spraying decreases compared to the cumulative spectral deviation before spraying, with a reduction rate. It is a positive number; when the attenuation ratio The higher the value, the more significant the relief from weed stress.

[0068] In practice, the calculated attenuation ratio is compared with a preset acceptable response threshold. The preset acceptable response threshold is determined based on the weed control efficacy acceptance standard for sunflower fields, and is set at 0.65. This threshold is based on the principle that in field weed control efficacy trials, when the attenuation of the sunflower canopy spectrum deviates from the cumulative amount by more than 65%, the stress of weeds on the sunflower has decreased to an economically acceptable level where further re-spraying is unnecessary. When the attenuation ratio is lower than the preset acceptable response threshold, incremental compensation adjustments are made to the herbicide demand potential of this heterogeneous spraying operation unit. The incremental compensation adjustment process includes: calculating the difference between the response qualification threshold and the attenuation ratio, and using this difference as the compensation intensity coefficient; multiplying the compensation intensity coefficient by a preset compensation step size factor to generate the compensation increment value; setting the compensation step size factor to 0.25, which is based on historical experience in supplementary spraying control, adopting a strategy that the single compensation amount does not exceed one-quarter of the theoretical maximum adjustment range of the demand potential, to avoid the risk of herbicide damage due to excessive one-time compensation; adding the compensation increment value to the current herbicide demand potential to generate the corrected herbicide demand potential; and re-inputting the corrected herbicide demand potential into the spraying amount decision model to generate the compensated spraying amount control instruction for executing the secondary supplementary spraying operation.

[0069] In practice, after generating the spraying rate control command, the variable spraying operation of the entire field is monitored cyclically. The cyclic monitoring is implemented by repeatedly executing all steps from generating the time-spectral joint feature trajectory of each spatial pixel to generating the spraying rate control command for each heterogeneous spraying operation unit at preset time intervals. The preset time interval is set as a growth cycle based on the sunflower growth rate and weed recovery rate, with a value of 7 days. The difference between the spraying rate control command generated in this cycle and the spraying rate control command generated in the previous cycle is compared. The difference is calculated by, for each heterogeneous spraying operation unit, calculating the absolute difference between the limited flow rate value generated in this cycle and the limited flow rate value generated in the previous cycle, and taking the maximum absolute difference among all heterogeneous spraying operation units as the difference. When the difference is lower than the preset convergence threshold, or when the number of cycles reaches the preset upper limit, the subsequent spraying control operation is stopped; the preset convergence threshold is set to 5% of the minimum stable flow of the variable spraying actuator, and the preset upper limit is set to 5 times; otherwise, the spraying volume control command generated in this cycle is sent to the variable spraying actuator to carry out a new round of spraying operation.

[0070] See Figure 10 In the figure, the horizontal axis represents the herbicide demand potential, ranging from 0 to 10, and the vertical axis represents the theoretical spraying flow rate (in L / min), reflecting the spraying flow rate calculated for different herbicide demand potentials. The figure uses scatter dots “○” to indicate the theoretical spraying flow rate values ​​of multiple sample points, reflecting the mapping result of the spraying rate decision model to the input demand potential; scatter dots “×” indicate spraying flow rate values ​​exceeding the flow rate constraint boundaries, showing cases where some theoretical flow rates exceed the maximum rated flow rate. The dashed lines in the figure represent the maximum rated flow rate (approximately 7.3 L / min) and minimum stable flow rate (approximately 0.35 L / min) of a single nozzle, serving as the upper and lower limits of the spraying flow rate, respectively. The solid lines are fitted straight lines, representing the linear mapping relationship of the spraying rate decision model, meeting the requirements of univariate linear regression, with a determination coefficient greater than or equal to 0.90, ensuring the reliability of the model.

[0071] As shown in the graph, the theoretical spraying flow rate increases linearly as the herbicide demand potential increases from 0 to approximately 10, indicating that the spraying flow rate decision model can effectively adjust the spraying flow rate according to the demand potential. When the spraying flow rate approaches the maximum rated flow rate, some theoretical spraying flow rates exceed the maximum rated flow rate boundary, indicated by an "×" mark. This suggests that in actual implementation, flow rate limiting is required, limiting the spraying flow rate to the maximum rated flow rate to avoid overloading the equipment. When the theoretical spraying flow rate is lower than the minimum stable flow rate, a zero flow rate command is output, indicating that the weed stress level in the corresponding area is too low, and herbicide spraying is unnecessary.

[0072] The spraying rate decision model shown in this figure is the linear mapping function described in Example 5. The linear regression coefficients obtained through field calibration experiments, combined with the single-nozzle flow rate constraint boundary, enable the generation of spraying flow rate control commands for each heterogeneous spraying operation unit. The key numerical ranges and curve trends in the figure clearly verify the applicability of the spraying rate decision model within different demand potential ranges and the rationality of the amplitude limiting treatment.

[0073] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition, characterized in that, include: Obtain raw multispectral image data of sunflower canopy at multiple consecutive time points, and generate temporal-spectral joint feature trajectory for each spatial pixel; A dynamic temporal-spectral feature decoupling operation is performed on the temporal-spectral joint feature trajectory of each spatial pixel to identify key spectral inflection points characterizing weed stress response; Based on the difference in spectral response within the time interval before and after the key spectral inflection point, calculate the cumulative spectral deviation of each spatial pixel; Based on the cumulative spectral deviation, spatial heterogeneity clustering of spatial pixels in the entire field is performed to dynamically divide multiple heterogeneous spraying operation units. Within each of the heterogeneous spraying operation units, the herbicide demand potential corresponding to that unit is inverted and deduced based on the cumulative spectral deviation of all spatial pixels contained therein. Based on the herbicide demand potential, a spraying volume control instruction is generated for each of the heterogeneous spraying operation units, and the spraying volume control instruction is used to drive the variable spraying actuator.

2. The method for optimizing herbicide application in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The generation of the temporal-spectral joint feature trajectory for each spatial pixel includes: The reflectance values ​​of each spatial pixel in multiple preset bands are extracted from the original multispectral image data to form the initial spectral feature vector of the pixel. The initial spectral feature vectors of the same spatial pixel at multiple consecutive time points are stacked in chronological order to construct the initial time-spectral feature matrix of the pixel. Perform principal component dimensionality reduction on the initial time-spectral feature matrix to extract the first few principal component components whose cumulative contribution rate meets the preset standard; The score sequence of the first few principal component components is used as the final time-spectral joint feature trajectory of the spatial pixel.

3. The method for optimizing herbicide application in sunflower fields based on multispectral recognition according to claim 2, characterized in that, The multiple preset bands include at least the green light band, the red light band, and the near-infrared band.

4. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The identification of key spectral inflection points characterizing weed stress response includes: The first-order difference calculation is performed on the time-spectral joint feature trajectory to obtain the sequence of the rate of change of feature values ​​with time; Locate local extreme points from the rate of change sequence, where the local extreme point is a location in the vicinity where the rate of change changes from positive to negative or from negative to positive. The time point corresponding to the local extreme point is taken as the candidate inflection point, and the change in the amplitude of the feature value before and after the candidate inflection point is extracted. When the change in the amplitude of the characteristic value exceeds the preset weed stress response threshold, the candidate inflection point is determined as the key spectral inflection point.

5. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 4, characterized in that, The preset weed stress response threshold is dynamically set based on the statistical distribution differences in the magnitude of characteristic value changes between healthy sunflower canopies and weed-stressed canopies.

6. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The calculation of the cumulative spectral deviation for each spatial pixel includes: Using the key spectral inflection point as the boundary, the time-spectral joint feature trajectory is divided into a feature subsequence before the inflection point and a feature subsequence after the inflection point; Calculate the mean center of the feature subsequence before the inflection point, and use it as the reference spectral baseline vector; Calculate the Euclidean distance between the feature vector at each time point in the feature subsequence after the inflection point and the reference spectral reference vector to generate a distance change sequence; The distance change sequence is integrated over time, and the accumulated result is used as the cumulative spectral deviation.

7. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 6, characterized in that, When performing time integration and accumulation on the distance change sequence, the trapezoidal integration method is used to improve the calculation accuracy of the accumulated amount.

8. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The dynamic division into multiple heterogeneous spraying operation units includes: using the cumulative spectral deviation of all spatial pixels as clustering features and inputting them into an adaptive spatial constraint clustering model; In the adaptive spatial constraint clustering model, a spatial adjacency weight is set between adjacent spatial cells, and the spatial adjacency weight is inversely proportional to the Euclidean distance between cells; The optimization objective is to minimize the variance of the cumulative spectral deviation within the same cluster and maximize the variance of the cumulative spectral deviation between different clusters, and to perform iterative clustering operations. When the iteration converges or reaches the preset number of iterations, the clustering results are output, and the spatial region corresponding to each cluster is a heterogeneous spraying operation unit.

9. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The inversion deduction of the herbicide demand potential corresponding to this unit includes: Obtain the statistical distribution parameters of the cumulative spectral deviation of all spatial pixels within the heterogeneous spraying operation unit, wherein the statistical distribution parameters include the median and interquartile range; The median is used as the central intensity of weed stress in the unit, and the interquartile range is used as the degree of spatial variability of weed stress in the unit. The preset stress-response inversion mapping relationship is invoked to map the weed stress center intensity and the weed stress spatial variability into an initial demand intensity value. The initial demand intensity value is compared with the preset minimum activation threshold for herbicides, and the larger of the two values ​​is taken as the final demand potential for herbicides.

10. The method for optimizing variable spraying of herbicides in sunflower fields based on multispectral recognition according to claim 1, characterized in that, The generation of spray volume control instructions for each of the heterogeneous spraying operation units includes: The herbicide demand potential of each heterogeneous spraying operation unit is input into the spraying rate decision model, which outputs a theoretical spraying flow rate that is linearly related to the potential value. The maximum rated flow and minimum stable flow of a single nozzle of the variable spraying actuator are read as the flow constraint boundary; The theoretical spraying flow rate is limited within the flow constraint boundary. When the theoretical spraying flow rate exceeds the maximum rated flow rate, the maximum rated flow rate is output. When it is lower than the minimum stable flow rate, a zero flow command is output. Based on the flow rate value after amplitude limiting, a duty cycle parameter corresponding to the pulse width modulation signal is generated, and this duty cycle parameter constitutes the spraying volume control command.

Images