Vegetable pest and disease warning method combining low-altitude remote sensing with environmental parameter analysis

Abstract

The application discloses a vegetable disease and pest early warning method fusing low-altitude remote sensing and environmental parameter analysis, and particularly relates to the technical field of agricultural information monitoring; by acquiring a first mapping relationship set F and a second mapping relationship set G, a stress initial value time sequence matrix P and an environmental parameter list set E are constructed, a lag correlation coefficient p between a stress initial value variation sequence and an environmental parameter sequence of each spatial position is calculated, and an environmental lag dominant area is identified; then, a correction coefficient a is calculated based on an environmental parameter ET of a Tth time sequence point, and the stress initial value PT of the Tth time sequence point is fused with the correction coefficient a to obtain a comprehensive early warning index Q; finally, the comprehensive early warning index Q is compared with an early warning threshold Qth, and an early warning signal is output; the method can reduce false positives caused by environmental time sequence lag interference, and improve the early warning accuracy of vegetable diseases and pests.

Description

Technical Field

[0001] This invention relates to the field of agricultural information monitoring technology, specifically to a method for early warning of vegetable diseases and pests that integrates low-altitude remote sensing and environmental parameter analysis. Background Technology

[0002] In vegetable cultivation in river valleys and hilly areas, early warning of diseases such as downy mildew and powdery mildew is easily affected by morning fog, diurnal temperature variations, and local microclimates. Existing technologies fall into two categories: one based on low-altitude remote sensing multispectral images, using vegetation indices to identify stress; however, changes in vegetation indices are highly coupled with environmental fluctuations such as temperature, humidity, and leaf wetness duration, making it difficult to distinguish between actual disease stress and environmental disturbances. The other category relies solely on environmental parameter modeling, which fails to reflect spatial heterogeneity. In particular, existing methods do not consider the temporal lag effect of environmental parameters on vegetation indices, easily leading to false alarms, and lack a technical solution for dynamically decoupling temporal environmental disturbances from actual physiological stress. Summary of the Invention

[0003] The purpose of this invention is to provide a method for early warning of vegetable diseases and pests that integrates low-altitude remote sensing and environmental parameter analysis, so as to overcome the shortcomings of the prior art.

[0004] To achieve the above objectives, the present invention provides the following technical solution: a vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis, comprising: S1, obtain the first mapping relationship set F between vegetation index and the degree of pest and disease stress, and the second mapping relationship set G between environmental parameters and the probability of pest and disease occurrence; S2, acquire multispectral images of the target area at T consecutive time points, extract the spectral reflectance of each sampling point, and calculate the initial stress value of each time point based on the first mapping relationship set F to obtain the initial stress value time series matrix P; S3, synchronously collect environmental parameters for the T consecutive time points to obtain an environmental parameter list set E; S4. Based on the stress initial value time series matrix P and the environmental parameter list set E, calculate the lag correlation coefficient ρ between the stress initial value variation sequence and the environmental parameter sequence at each spatial location, and mark the corresponding spatial location when ρ is greater than the preset threshold ρth. S5. Based on the second mapping relationship set G and the environmental parameter ET of the Tth time series point, calculate the correction coefficient α. For unmarked spatial locations, fuse the initial stress value PT of the Tth time series point with the correction coefficient α to obtain the comprehensive early warning index Q. For marked spatial locations, set Q=0. S6 compares the comprehensive early warning index Q with the early warning threshold Qth. When Q is greater than or equal to Qth, it is determined that pests and diseases have occurred at the corresponding sampling point and an early warning signal is output.

[0005] Preferably, step S1 includes: selecting vegetable sample areas in healthy, mildly stressed, moderately stressed, and severely stressed states, and spatially locating each sample area; based on the spatial location, simultaneously collecting multispectral images, temperature, air humidity, leaf surface wetness duration, and soil moisture content data for each sample area, and obtaining corresponding disease level labeling results; extracting vegetation index features based on the multispectral images, and establishing a correspondence between vegetation indices and the degree of disease and pest stress in combination with the disease level labeling results; and establishing a correspondence between environmental parameters and the probability of disease and pest occurrence using a time-segmented fitting method based on the environmental parameter data, historical disease occurrence records, and corresponding disease level labeling results, thereby generating the first mapping relationship set F and the second mapping relationship set G respectively.

[0006] Preferably, the step of calculating the initial stress values ​​for each time series point based on the first mapping relationship set F to obtain the initial stress value time series matrix P includes: repeatedly conducting aerial surveys of the target area at preset time intervals to acquire multispectral images of each time series point; sequentially performing radiometric calibration, geometric correction, and spatial registration on the multispectral images of each time series point to obtain standardized images under a unified coordinate system; extracting the spectral reflectance of each sampling point in each band based on the standardized images and calculating the corresponding vegetation index; inputting the vegetation index of each sampling point into the first mapping relationship set, calculating the initial stress values ​​corresponding to each time series point, and arranging them according to time order and spatial location to generate the initial stress value time series matrix P.

[0007] Preferably, the step of inputting the vegetation index of each sampling point into the first mapping relationship set and calculating the initial stress value corresponding to each time series point includes: combining multiple vegetation indices corresponding to each sampling point according to a preset weight to construct a vegetation index feature vector of the sampling point; performing interval normalization processing on the vegetation index feature vector and performing feature matching according to the preset stress classification interval in the first mapping relationship set to determine the candidate stress level corresponding to each sampling point; determining the upper bound reference interval and lower bound reference interval adjacent to the candidate stress level in the first mapping relationship set according to the candidate stress level; extracting the interval position parameters of the vegetation index feature vector of the sampling point in the upper bound reference interval and lower bound reference interval, and constructing a piecewise interpolation function accordingly; and using the piecewise interpolation function to perform continuous numerical mapping on the vegetation index feature vector of the sampling point to obtain the initial stress value of each time series point.

[0008] Preferably, the calculation of the lag correlation coefficient ρ between the initial stress value variation sequence and the environmental parameter sequence at each spatial location includes: According to the spatial location, the initial stress values ​​of each sampling point at multiple consecutive time points are extracted from the initial stress value time series matrix, and the corresponding initial stress value variation sequence is constructed; based on the spatial location of the sampling points, the environmental parameter data of the corresponding time points in the environmental parameter list are matched to construct the environmental parameter sequence corresponding to the initial stress value variation sequence. The environmental parameter sequence is shifted hourly according to a preset lag time window, and the correlation value between the shifted environmental parameter sequence and the stress initial value variation sequence is calculated respectively; the maximum value among the correlation values ​​is selected as the lag correlation coefficient ρ of the corresponding spatial location.

[0009] Preferably, the step of calculating the correlation between the shifted environmental parameter sequence and the initial stress value variation sequence includes: performing detrending and amplitude normalization processing on each shifted environmental parameter sequence to obtain a standard environmental parameter subsequence; extracting the effective stress subsequence for the corresponding time period based on the initial stress value variation sequence aligned with the standard environmental parameter subsequence; using the standard environmental parameter subsequence and the effective stress subsequence as input, obtaining the local correlation value at each window position using a sliding window correlation calculation method, and weighting and accumulating each local correlation value to obtain the correlation value corresponding to the current shift amount.

[0010] Preferably, the step of calculating the correction coefficient α based on the second mapping relationship set G and the environmental parameter ET at the Tth time point includes: extracting the environmental parameter ET corresponding to each unmarked spatial location at the Tth time point, and inputting it into the second mapping relationship set for probability mapping to obtain the environmental risk value corresponding to each unmarked spatial location; and calculating the correction coefficient α corresponding to each unmarked spatial location based on the environmental risk value according to a preset nonlinear transformation rule.

[0011] Preferably, the comprehensive early warning index Q is obtained by fusing the initial stress value PT at the Tth time point with the correction coefficient α, including: constraining the correction coefficient within a preset numerical range; fusing the initial stress value corresponding to each unmarked spatial location at the Tth time point with the corresponding correction coefficient to obtain the initial early warning value corresponding to each unmarked spatial location; and then performing neighborhood consistency correction based on the distribution continuity of the initial early warning value in adjacent spatial locations to output the comprehensive early warning index Q corresponding to each unmarked spatial location.

[0012] Preferably, the calculation of the correction coefficient α corresponding to each unmarked spatial location based on the environmental risk value according to a preset nonlinear transformation rule includes: inputting the environmental risk value corresponding to each unmarked spatial location into a preset piecewise nonlinear response function to obtain the corresponding initial correction value; performing environmental sensitivity modulation on the initial correction value according to the temperature, air humidity, and leaf wetness duration corresponding to each unmarked spatial location at the current time point to obtain an intermediate correction value; then performing time-series smoothing constraint on the intermediate correction value according to the environmental risk change amplitude of each unmarked spatial location at multiple consecutive time points to obtain a stable correction value; finally, performing upper and lower limit truncation and monotonic consistency correction on the stable correction value to output the correction coefficient α corresponding to each unmarked spatial location.

[0013] Preferably, the step of fusing the initial stress value and the corresponding correction coefficient corresponding to each unmarked spatial location at the Tth time point includes: extracting the initial stress value and correction coefficient corresponding to each unmarked spatial location at the Tth time point, and normalizing the initial stress value and the correction coefficient to obtain a fusionable input quantity; constructing a nonlinear coupled fusion function based on the normalized initial stress value and correction coefficient, and calculating the initial warning value corresponding to each unmarked spatial location.

[0014] The technical effects and advantages provided by the present invention in the above technical solution are as follows: 1. This invention constructs a first mapping relationship set F and a second mapping relationship set G, and conducts lag correlation analysis using the initial stress value time series matrix P and the environmental parameter list set E at T consecutive time points. This effectively decouples environmental time series lag interference from actual pest and disease physiological stress. Compared to existing methods that rely solely on multispectral imagery or single environmental parameter modeling, this invention introduces a lag correlation coefficient ρ between the initial stress value variation sequence and the environmental parameter sequence, and uses a preset threshold ρth to mark the dominant environmental stress areas. This eliminates spatial locations significantly affected by factors such as morning fog, short-term high humidity, or temperature fluctuations at the source, significantly reducing the probability of false alarms. Especially in complex environments with large diurnal temperature differences and frequent morning fog in river valleys and hilly areas, this invention avoids misjudgments caused by the lag recovery of canopy reflectivity, making the early warning results closer to the actual occurrence of pests and diseases.

[0015] 2. This invention introduces an environmental risk value based on a second mapping relationship set G, and utilizes a piecewise nonlinear response function, environmental sensitivity modulation, temporal smoothing constraints, and a nonlinear coupling fusion function to fuse the initial stress value PT at the T-th time-series point with the correction coefficient α, obtaining a comprehensive early warning index Q. This achieves a synergistic quantitative expression of environmental conditions and vegetation physiological states. Through neighborhood consistency correction and early warning threshold Qth determination rules, the spatial continuity and stability of early warning results are further enhanced, effectively improving the accuracy and reliability of early pest and disease identification. This invention not only reduces unnecessary pesticide application and lowers production costs, but also reduces pesticide residue risks and soil ecological damage, demonstrating significant agricultural production application value and ecological environmental benefits. Attached Figure Description

[0016] 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.

[0017] Figure 1 This is a flowchart of the vegetable pest and disease early warning method that integrates low-altitude remote sensing and environmental parameter analysis according to the present invention.

[0018] Figure 2 This is a flowchart of the method for calculating the lag correlation coefficient of the present invention. Detailed Implementation

[0019] 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.

[0020] For examples, please refer to Figure 1 and Figure 2 As shown in this embodiment, the vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis includes: In this embodiment, step S1 is used to generate a first mapping relationship set F between vegetation index and the degree of pest and disease stress, and a second mapping relationship set G between environmental parameters and the probability of pest and disease occurrence, so as to provide a unified data foundation for subsequent initial stress value calculation and environmental correction.

[0021] Specifically, firstly, vegetable sample areas under healthy, mildly stressed, moderately stressed, and severely stressed conditions were selected within the target vegetable planting area. The sample areas were selected by dividing the field into zones based on topography, planting density, and variety distribution. Within each zone, a representative plant group was selected as a vegetable sample area, with each sample area ranging from 4 to 25 square meters. Spatial positioning was achieved using a combination of GPS and ground markers. First, the latitude and longitude coordinates of the center point of each vegetable sample area were recorded using GPS, and then the boundaries were corrected using ground markers to ensure a one-to-one correspondence between the subsequent multispectral image coverage and each vegetable sample area. Disease severity labeling results were obtained through a combination of manual investigation and laboratory testing. A healthy state corresponds to no disease spots and normal leaf color; a mildly stressed state corresponds to disease spots covering less than 10% of the leaf area; a moderately stressed state corresponds to disease spots covering 10% to 30% of the leaf area; and a severely stressed state corresponds to disease spots covering more than 30% of the leaf area. The above disease level labeling results serve as the benchmark labels for the subsequent construction of the first mapping relationship set F and the second mapping relationship set G.

[0022] After spatial positioning is completed, multispectral images, temperature, air humidity, leaf wetness duration, and soil moisture content data for each vegetable sample area are simultaneously collected based on the spatial positioning, and corresponding disease level labeling results are obtained. The multispectral images are acquired by a low-altitude remote sensing platform equipped with multi-band sensors, including at least the red light band, near-infrared band, and red-edge band; temperature and air humidity are acquired using field micro-meteorological sensors; leaf wetness duration is accumulated using leaf wetness sensors; and soil moisture content is obtained using buried moisture probes. For each vegetable sample area, continuous collection is conducted for no less than 7 days, with no less than 4 collections per day, thus forming an original training sample set containing multispectral images, environmental parameter data, and disease level labeling results. For example, in river valleys and hilly areas where morning fog is frequent, data can be collected at 6:00, 10:00, 14:00, and 18:00 to cover typical environmental change stages such as the impact of morning fog, post-fog recovery, and afternoon high temperatures.

[0023] After obtaining the original training sample set, vegetation index features are extracted based on the multispectral images, and a correspondence between vegetation index and the degree of pest and disease stress is established in combination with the disease level labeling results, thereby generating the first mapping relationship set F.

[0024] In practice, the multispectral images are first subjected to radiometric consistency correction and spatial cropping to extract the image regions corresponding to each vegetable sample area. Then, vegetation index features are calculated based on the reflectance information of each band. The vegetation index features may include at least two of the following: normalized difference vegetation index, red-edge difference vegetation index, and canopy color difference index. To ensure comparability between different vegetation index features, each vegetation index feature is normalized according to the minimum and maximum values ​​of the overall sample. The normalization method is as follows: subtract the minimum value of the current vegetation index feature from the minimum value of the vegetation index feature in all samples, and then divide by the difference between the maximum and minimum values ​​of the vegetation index feature in all samples, so that each vegetation index feature falls within the range of 0 to 1.

[0025] Based on the disease severity labeling results, all samples were divided into four categories: healthy, mildly stressed, moderately stressed, and severely stressed. The center value, fluctuation range, and boundary intervals of each vegetation index feature were statistically analyzed for each category. The first mapping relationship set F was constructed as follows: using the disease severity labeling results as the classification axis, the normalized intervals of the vegetation index features as the input range, and configuring a continuous stress value output interval for each disease severity labeling result. The values ​​for healthy status corresponded to 0 to 0.25, mildly stressed status to 0.25 to 0.5, moderately stressed status to 0.5 to 0.75, and severely stressed status to 0.75 to 1. For samples located in the boundary intervals of two adjacent categories, a linear interpolation method was used to calculate the continuous stress value. That is, based on the position of the current vegetation index feature value between the lower and upper boundary values, its specific position in the corresponding continuous stress value output interval was determined proportionally.

[0026] For example, if the vegetation index characteristics of a vegetable sample area fall within the boundary region between mild and moderate stress, and its location is closer to the moderate stress range, then the corresponding continuous stress value for this sample can be assigned a value of 0.58 or 0.62 to reflect that it is close to moderate stress but has not yet fully entered the moderate stress state. Through the above method, a first mapping relationship set F can be formed, consisting of "vegetation index characteristic range, pest and disease stress level, and continuous stress value".

[0027] After establishing the first mapping relationship set F, based on the environmental parameter data, disease history records, and corresponding disease level labeling results, a time-segmented fitting method is used to establish the correspondence between environmental parameters and the probability of disease occurrence, thereby generating the second mapping relationship set G.

[0028] The specific implementation of the time-segmented fitting method is as follows: First, the sampling time is divided into multiple fixed time periods according to the daily environmental change pattern, preferably into early morning, morning, afternoon, and nighttime periods, for example, 0:00 to 6:00, 6:00 to 12:00, 12:00 to 18:00, and 18:00 to 24:00. Then, within each time period, data on temperature, air humidity, leaf surface wetness duration, and soil moisture content are extracted, and combined with the corresponding disease level labeling results and disease history occurrence records to form a time period sample subset. The disease history occurrence records can consist of field inspection records from the past three planting cycles, disease confirmation dates, and occurrence locations. For each time period sample subset, temperature, air humidity, leaf surface wet duration, and soil moisture content are first divided into multiple parameter intervals. For example, temperature can be divided into intervals of 2 degrees Celsius, air humidity into intervals of 5%, leaf surface wet duration into intervals of 1 hour, and soil moisture content into intervals of 5% volumetric moisture content. Then, the number of samples with disease occurrence under each parameter interval combination is counted against the total number of samples. The probability of disease and pest occurrence corresponding to that parameter interval combination is calculated by dividing the number of disease samples by the total number of samples.

[0029] To improve the continuity and applicability of the second mapping set G, the probability of pest and disease occurrence corresponding to each parameter interval combination is further subjected to curve fitting. The curve fitting adopts a piecewise logical growth fitting method, that is, within each time period, the combination of environmental parameters is used as input, and the probability of pest and disease occurrence is used as output. The fitting curve parameters are determined by minimizing the deviation between the fitted output and the statistical probability. When the air humidity is higher than 90%, the leaf wetness duration is greater than 4 hours, and the temperature is between 18°C ​​and 24°C within a certain time period, the fitted probability of pest and disease occurrence is usually higher than 0.7; when the air humidity is lower than 60%, the leaf wetness duration is less than 1 hour, and the temperature is higher than 30°C, the fitted probability of pest and disease occurrence is usually lower than 0.2. In this way, the second mapping set G can be specifically represented as a set of correspondences between "time period category, combination of environmental parameter intervals, and probability of pest and disease occurrence".

[0030] In this embodiment, both the first mapping set F and the second mapping set G can be stored in the same data table or in different data tables. The first mapping set F includes at least four types of fields: vegetation index feature type, normalized interval, pest and disease stress level category, and continuous stress value interval. The second mapping set G includes at least six types of fields: time period category, temperature interval, air humidity interval, leaf surface wet duration interval, soil moisture content interval, and pest and disease occurrence probability. When the first mapping set F is called in subsequent steps, the corresponding stress value can be directly queried or interpolated based on the input vegetation index features; when the second mapping set G is called, the corresponding pest and disease occurrence probability can be output based on the interval and time period category of the current environmental parameters.

[0031] For example, in the early morning hours in a valley or hilly area, if the air humidity of a vegetable sample area is 93%, the leaf surface is moist for 5 hours, the temperature is 20 degrees Celsius, and the soil moisture content is 28%, then the corresponding combination of environmental parameters can be matched in the second mapping relationship set G, and a higher probability of pest and disease occurrence can be obtained for subsequent environmental correction calculations.

[0032] In this embodiment, step S2 is used to acquire multispectral images of the target area at T consecutive time points, and convert the vegetation index of each sampling point into the initial stress value of each time point based on the first mapping relationship set F. Furthermore, a stress initial value time series matrix P is constructed according to the time sequence and spatial location. The T consecutive time points refer to T sampling times set sequentially at fixed time intervals within the same monitoring period. The preset time interval is determined based on the duration of morning fog, leaf moisture decay time, and the operational capability of the low-altitude remote sensing platform in the target area, and is preferably 30 to 120 minutes. In application scenarios where morning fog is frequent in river valleys and hilly areas, the preset time interval can be set to 60 minutes, and eight time points can be continuously set between 6:00 and 13:00 to enhance the ability to capture the lag changes in canopy reflectivity after the dissipation of morning fog.

[0033] When conducting repeated aerial surveys of the target area at preset time intervals, multiple acquisitions are performed using the same low-altitude remote sensing platform, the same flight path, the same flight altitude, and the same flight speed to ensure comparability between multispectral images at different time points. Specifically, the flight altitude can be set to 40 to 80 meters, and the flight path overlap rate can be set to a forward overlap rate of over 70% and a lateral overlap rate of over 60%. Before each aerial survey, the start time, meteorological conditions, and sensor status information are recorded. To reduce the impact of changes in solar altitude angle on the consistency of multispectral image brightness, automatic exposure locking and fixed gain acquisition are preferably used for repeated aerial surveys within the same monitoring period. After acquisition at each time point, corresponding multispectral images are obtained, each containing information from at least two of the following bands: red, near-infrared, and red-edge.

[0034] Radiometric calibration, geometric correction, and spatial registration were performed sequentially on the multispectral images at each time point to obtain standardized images under a unified coordinate system. The radiometric calibration was implemented as follows: before each repeated aerial survey, a standard reflector with known reflectivity was deployed at the edge of the target area. The measured gray values ​​of the standard reflector in each band were recorded. Based on the correspondence between the known reflectivity of the standard reflector and the measured gray values, the image pixel gray values ​​were converted into ground surface reflectivity. Radiometric calibration was then performed band by band on the multispectral images obtained from the same aerial survey.

[0035] For example, if the known reflectivity of a standard reflector in the near-infrared band is 0.50, and the average pixel value corresponding to the measured grayscale value in this band is 5000, a conversion relationship between near-infrared grayscale value and reflectivity can be established, and the pixel grayscale values ​​within the target area can be uniformly converted to near-infrared reflectivity values. Geometric correction is implemented by calling the position and attitude data recorded by the low-altitude remote sensing platform and combining it with the coordinates of ground control points deployed at the four corners and center of the target area to perform distortion correction and geolocation correction on the multispectral image, ensuring that the location of ground features in the image corresponds to the actual ground location. Spatial registration is implemented by selecting the geometrically corrected multispectral image of the first time-series point as the reference image, establishing a set of corresponding points using stable feature points such as road edges, field boundaries, and ground landmarks between the images of each time-series point, and then using nearest neighbor resampling or bilinear resampling to register the remaining time-series point images to the unified coordinate system of the reference image, ultimately obtaining a standardized image with consistent pixel size, consistent spatial range, and one-to-one correspondence between rows and columns.

[0036] When extracting the spectral reflectance of each sampling point in each band based on the standardized image and calculating the corresponding vegetation index, the pixel position corresponding to each sampling point is first located in each standardized image according to a pre-set sampling grid or based on the field-calibrated sampling point coordinates. To reduce the influence of individual pixel noise on the calculation results, the average value of 3×3 neighboring pixels centered on the sampling point's center pixel can be used as the spectral reflectance of that sampling point in the corresponding band. Subsequently, vegetation indices are calculated based on the vegetation index type used when constructing the first mapping relationship set F, using the spectral reflectance of each sampling point. The vegetation indices preferably include the Normalized Difference Vegetation Index, the Red Edge Difference Vegetation Index, and the Canopy Color Difference Index. The Normalized Difference Vegetation Index (NDVI) can be calculated as follows: the difference between near-infrared reflectance and red reflectance is divided by the sum of near-infrared and red reflectance. The Red Edge Difference Vegetation Index (REDVI) can be calculated as follows: the degree of difference between near-infrared and red edge reflectance characterizes leaf activity changes. The Canopy Color Difference Index (CDI) can be calculated as follows: the combination of red, green, and blue light bands reflects the degree of canopy color shift. For example, when the near-infrared reflectance of a sampling point at a certain time point is 0.62 and the red reflectance is 0.24, the NDVI can be calculated by subtracting 0.24 from 0.62 and then dividing by the sum of 0.62 and 0.24. The resulting value is used to characterize the canopy vitality level at that sampling point.

[0037] When inputting the vegetation indices of each sampling point into the first mapping relationship set F and calculating the initial stress value corresponding to each time series point, the vegetation indices corresponding to each sampling point are first combined according to preset weights to construct the vegetation index feature vector of the sampling point. The construction method of the preset weights is consistent with the training samples in step S1, specifically: based on the four types of vegetable sample areas of healthy, mildly stressed, moderately stressed, and severely stressed, the distinguishing contribution of each vegetation index between different levels of pest and disease stress is calculated; then, it is normalized according to the method of "dividing the distinguishing contribution of each vegetation index by the sum of the distinguishing contributions of all vegetation indices" to obtain the preset weight corresponding to each vegetation index. For example, when the distinguishing contribution of the normalized differential vegetation index, the red edge differential vegetation index, and the canopy color differential index accounts for 0.40, 0.35, and 0.25 of the total contribution, respectively, the preset weights of the three are set to 0.40, 0.35, and 0.25, respectively. When constructing the vegetation index feature vector of the sampling points, the current value of each vegetation index is matched with its corresponding preset weight one by one, and arranged in a fixed order to form a vector, so as to ensure that the input structure of different sampling points and different times is completely consistent.

[0038] After constructing the vegetation index feature vector, the vegetation index feature vector is subjected to interval normalization, and feature matching is performed according to the preset stress grading intervals in the first mapping relationship set to determine the candidate stress level corresponding to each sampling point. The interval normalization process adopts the same normalization benchmark as in step S1, specifically: for each dimension component in the vegetation index feature vector, the current vegetation index value is subtracted from the minimum value of the vegetation index in the training sample, and then divided by the difference between the maximum and minimum values ​​of the vegetation index in the training sample, to transform each dimension component to the interval of 0 to 1. After completing the interval normalization, the stress grading intervals stored in the first mapping relationship set F for healthy, mild stress, moderate stress, and severe stress are read respectively. The stress grading interval includes the lower boundary value, upper boundary value, and center value of each vegetation index under the corresponding pest and disease stress level. During feature matching, it is determined whether the normalized components of each dimension of the current sampling point fall within the stress grading interval corresponding to the severity of each pest or disease stress. If the hit dimension corresponding to a certain pest or disease stress level reaches more than 70% of the total dimensions of the vegetation index, then the pest or disease stress level is determined as a candidate stress level. If two or more pest or disease stress levels simultaneously meet the hit condition, the weighted deviation between the current vegetation index feature vector and the central value vector of each pest or disease stress level is further compared, and the pest or disease stress level with the smallest weighted deviation is determined as a candidate stress level. The calculation logic of the weighted deviation is as follows: the absolute difference between each normalized component and its corresponding central value is multiplied by the preset weight of that dimension, and then summed over all dimensions. This method can maintain the stability of the candidate stress level determination even when there are local fluctuations in each vegetation index.

[0039] Based on the candidate stress levels, upper and lower bound reference intervals adjacent to the candidate stress levels are determined in the first mapping relationship set F. Specifically, when the candidate stress level is mild stress, the lower bound reference interval is the stress level interval corresponding to health, and the upper bound reference interval is the stress level interval corresponding to moderate stress; when the candidate stress level is moderate stress, the lower bound reference interval is the stress level interval corresponding to mild stress, and the upper bound reference interval is the stress level interval corresponding to severe stress. For health and severe stress at the boundary positions, the first mapping relationship set F pre-stores endpoint extension reference intervals, where the lower bound reference interval corresponding to health is the health endpoint extension reference interval, and the upper bound reference interval corresponding to severe stress is the severe stress endpoint extension reference interval, thereby ensuring that any candidate stress level can correspond to a set of upper and lower bound reference intervals. The above-mentioned endpoint extension reference intervals are constructed by using the interval width of adjacent stress level intervals as the extension width, and extending the endpoints of health or severe stress outward by the same width to form a reference interval.

[0040] After determining the upper and lower reference intervals, the interval position parameters of the vegetation index feature vector of the sampling points within these intervals are extracted, and a piecewise interpolation function is constructed accordingly. These interval position parameters characterize the relative position of the current sampling point between adjacent stress grading intervals. The specific calculation logic is as follows: for each dimension of the vegetation index feature vector, the relative distance between the upper boundary of the lower reference interval and the lower boundary of the upper reference interval is first determined. Then, this relative distance is divided by the total distance between the two boundaries to obtain the interval position parameter between 0 and 1. When a dimension is less than the upper boundary of the lower reference interval, the interval position parameter is 0; when a dimension is greater than the lower boundary of the upper reference interval, the interval position parameter is 1. For example, if the normalized vegetation index of a certain dimension at a certain sampling point is 0.46, and the upper boundary of the corresponding lower reference interval is 0.40 and the lower boundary of the upper reference interval is 0.60, then the position parameter of this dimension can be calculated by subtracting 0.40 from 0.46 and then dividing by 0.60 minus 0.40. The result is 0.30, which means that this dimension component is located on the side closer to the lower reference interval between adjacent hierarchical boundaries.

[0041] When constructing the piecewise interpolation function, the continuous stress value output interval corresponding to the candidate stress level is used as the main interval, and the interval position parameters of each dimension and the preset weights of each dimension are used as function inputs. The specific implementation of the piecewise interpolation function is as follows: first, the interval position parameters of each dimension are multiplied by their respective preset weights, and then the weighted interval position parameters are summed to obtain the overall position coefficient; then, based on the starting value and ending value of the continuous stress value output interval corresponding to the candidate stress level, the continuous numerical mapping result of the sampling point is generated by adding the overall position coefficient to the starting value and multiplying it by the interval width. When the candidate stress level is mild stress, its continuous stress value output interval can be 0.25 to 0.50; when the candidate stress level is moderate stress, its continuous stress value output interval can be 0.50 to 0.75. If a sampling point has a candidate stress level of mild stress at a certain time series point and an overall location coefficient of 0.60, its continuous numerical mapping result can be calculated by adding 0.25 to 0.60 and multiplying by 0.25, resulting in 0.40. This value is the initial stress value for that sampling point at that time series point. Through the above piecewise interpolation function, the originally discrete pest and disease stress levels can be converted into continuous initial stress values, thereby improving the sensitivity of subsequent time series analysis to subtle changes.

[0042] After continuously mapping the vegetation index feature vectors of the sampling points using the piecewise interpolation function, the initial stress values ​​for each time series point are obtained. For each time series point, the initial stress values ​​of all sampling points within the target area are backfilled according to their spatial location to form the spatial distribution of the initial stress values ​​for that time series point. Then, the spatial distributions of the initial stress values ​​from the 1st time series point to the Tth time series point are arranged sequentially according to time order to generate the initial stress value time series matrix P. The initial stress value time series matrix P can be represented as a time series set consisting of P1, P2, up to PT, where P1 represents the spatial distribution of the initial stress values ​​at the 1st time series point, and PT represents the spatial distribution of the initial stress values ​​at the Tth time series point. If the target area is divided into a sampling grid of 50 rows and 60 columns, and 8 time series points are continuously monitored, then each time series point corresponds to a 50-row by 60-column spatial distribution of initial stress values, and the initial stress value time series matrix P is composed of 8 spatial distributions of 50 rows by 60 columns. Subsequent steps can directly call the numerical sequence of any spatial location in the initial stress time series matrix P at T consecutive time points to analyze the dynamic response process of that spatial location under the combined effects of environmental changes and pest and disease stress.

[0043] In this embodiment, step S3 is used to synchronously collect environmental parameters at the T consecutive time points to obtain an environmental parameter list set E, thereby providing a time-consistent data foundation for calculating the lag correlation coefficient based on the stress initial value time series matrix P and the environmental parameter list set E.

[0044] Specifically, environmental parameter acquisition devices are set up within the target area according to the spatial distribution of sampling points. These environmental parameters include at least temperature, air humidity, leaf surface wetness duration, and soil moisture content. Temperature and air humidity are acquired in real-time by microclimate sensors, leaf surface wetness duration is recorded by leaf surface wetness sensors using a cumulative time-based method, and soil moisture content is continuously measured by a soil moisture probe buried in the root activity layer. Synchronous acquisition means that the recording time of each time-series environmental parameter is consistent with the acquisition time of the corresponding time-series multispectral image, with the time deviation between the two controlled within 5 minutes. When the environmental parameter acquisition frequency is higher than the multispectral image acquisition frequency, the environmental parameter value closest to the current time-series point is selected, or a time-weighted average of the environmental parameter values ​​at adjacent times before and after the current time-series point is used to determine the environmental parameter at that time-series point, ensuring the correspondence between the environmental parameters and the initial stress value in the time dimension.

[0045] After collecting environmental parameters at each time point, the environmental parameters corresponding to each time point are organized in chronological order to obtain an environmental parameter list set E = (E1, E2, ..., ET), where E1 represents the environmental parameter list corresponding to the first time point, and ET represents the environmental parameter list corresponding to the Tth time point. Each environmental parameter list contains the values ​​of temperature, air humidity, leaf wetness duration, and soil moisture content corresponding to each spatial location within the target area, and maintains a one-to-one correspondence with the spatial location in the initial stress time series matrix P. To improve the accuracy of subsequent calculations, before generating the environmental parameter list set E, outlier removal and missing value completion can be performed on the environmental parameters collected at each time point. Outlier removal can be achieved by judging whether the value exceeds a preset fluctuation limit compared to two adjacent time points, and missing value completion can be achieved by interpolation between adjacent time points at the same spatial location. For example, if the air humidity data for a certain spatial location is missing at the 5th time point, the average of the air humidity values ​​from the 4th and 6th time points can be used as the air humidity value for the 5th time point and written into E5. The environmental parameter list set E formed in this way can accurately reflect the environmental change process of the target area over T consecutive time points and provide a unified input for subsequent lag correlation analysis and comprehensive early warning index calculation.

[0046] In this embodiment, step S4 is used to calculate the lag correlation coefficient ρ between the stress initial value variation sequence and the environmental parameter sequence at each spatial location based on the stress initial value time series matrix P and the environmental parameter list set E, and when ρ is greater than a preset threshold ρth, the corresponding spatial location is marked to identify the spatial location that is significantly affected by the environmental time series lag.

[0047] Specifically, the initial stress values ​​of each sampling point at multiple consecutive time points are first extracted from the initial stress value time series matrix P according to their spatial location. If the initial stress values ​​of a sampling point at the 1st to the Tth time series points are 0.32, 0.35, 0.41, 0.39, and 0.46 respectively, then the initial stress value variation sequence of that sampling point is obtained by constructing the difference between adjacent time series points, that is, subtracting the initial stress value of the previous time series point from the initial stress value of the next time series point, resulting in 0.03, 0.06, -0.02, and 0.07 respectively. For each spatial location in the initial stress value time series matrix P, the corresponding initial stress value variation sequence is constructed in the same way, thereby forming a time series input that can be aligned point by point with the environmental parameter change process.

[0048] After obtaining the initial stress value variation sequence, environmental parameter sequences corresponding to the initial stress value variation sequence are constructed by matching the spatial location of the sampling points with the environmental parameter data of the corresponding time points in the environmental parameter list set E. The environmental parameter sequence is formed by arranging the values ​​of temperature, air humidity, leaf wetness duration, and soil moisture content at multiple consecutive time points corresponding to the sampling point in chronological order, maintaining consistency with the record order in the environmental parameter list set E. After spatial matching, an environmental parameter sequence with a time span consistent with the initial stress value variation sequence is formed for each sampling point.

[0049] In this embodiment, the environmental parameter sequence is shifted hourly according to a preset lag time window to characterize the leading influence of environmental parameter changes on the initial stress value changes. The preset lag time window is determined based on the actual duration of delayed recovery of canopy reflectivity after the dissipation of morning fog in river valleys and hilly areas, preferably covering 1 to 6 time points. When the preset time interval is 60 minutes, the preset lag time window corresponds to a lag range of 1 to 6 hours. The specific implementation of the hourly shift is as follows: keeping the time position of the initial stress value variation sequence unchanged, the environmental parameter sequence is shifted forward by 1 time point, 2 time points, up to 6 time points; after each shift, the overlapping segment aligned with the current shift result is extracted as the data interval for participating in the correlation value calculation under the current shift amount.

[0050] For each set of shifted environmental parameter sequences, detrending and amplitude normalization are performed to obtain standard environmental parameter subsequences. The detrending process is performed separately for temperature, air humidity, leaf wetness duration, and soil moisture content. Specifically, using data points of environmental parameters changing over time within the overlapping segment corresponding to the current shift amount as input, a trend line for each environmental parameter within that segment is established using least-squares linear fitting. Then, the fitted value on the trend line at each time point is subtracted from the original parameter value, resulting in a residual sequence after removing long-term monotonic changes. The amplitude normalization process uses zero-mean unit variance normalization. First, the mean of each environmental parameter within the current overlapping segment is calculated. Then, the mean is subtracted from the parameter value at each time point and divided by the standard deviation of the environmental parameter within the current overlapping segment, thus allowing environmental parameters with different dimensions to participate in subsequent calculations on a unified scale. If the standard deviation of an environmental parameter within the current overlapping segment is less than 0.01, it is considered that its variation is too small to form effective comparative information, and the normalized value of that environmental parameter under the current shift is set to 0. After completing the above processing, the normalized results of each environmental parameter within the current overlapping segment are reorganized into a standard environmental parameter subsequence according to a fixed order of temperature, air humidity, leaf wetness duration, and soil moisture content.

[0051] Based on the stress initial value variation sequence aligned with the standard environmental parameter subsequence in time, the effective stress subsequence for the corresponding time period is extracted. The length of the effective stress subsequence is consistent with the time length of the standard environmental parameter subsequence under the current shift, and corresponds one-to-one with the standard environmental parameter subsequence at each time point. To avoid excessive amplification of local correlation calculations due to occasional spikes in the stress initial value variation sequence, after extracting the effective stress subsequence, outliers exceeding three times the absolute deviation of the overall median of the stress initial value variation sequence at this sampling point are replaced with the mean of the two nearest time points.

[0052] For example, if the effective stress subsequence of a sampling point shows a sudden increase of 0.35 at a certain time point, while the initial stress value variation of most time points at that sampling point is between -0.05 and 0.08, then the 0.35 is identified as an outlier and replaced with the average of the corresponding values ​​at the two adjacent time points. The effective stress subsequence obtained in this way can more stably reflect the stress change process of each sampling point within the current comparison segment.

[0053] Using the standard environmental parameter subsequence and the effective stress subsequence as input, a sliding window correlation calculation method is employed to obtain the local correlation values ​​at each window position. The length of the sliding window is determined based on the length of the current overlapping segment, preferably between 1 / 2 and 2 / 3 of the current overlapping segment length, and not less than 3 time points; when the current overlapping segment length is 6 time points, the sliding window length can be 3 time points. The sliding window moves from front to back in steps of 1 time point, and at each window position, the local correlation between temperature, air humidity, leaf surface wet duration, and soil moisture content and the effective stress subsequence is calculated. The local correlation is calculated as follows: First, the normalized value of a certain environmental parameter within the window and the corresponding effective stress subsequence value are averaged respectively. Then, the deviation of the environmental parameter from its average value at each time point is multiplied by the deviation of the effective stress subsequence from its average value and summed. Finally, the sum is divided by the square root of the product of the squares of the deviations of the two values ​​from their respective average values ​​to obtain the correlation coefficient of the environmental parameter at the current window position.

[0054] Considering the varying degrees of influence of different environmental parameters on the characterization of pest and disease stress, when forming local correlation values, environmental parameters such as temperature, air humidity, leaf wetness duration, and soil moisture content are assigned weights. The construction method for these environmental parameter weights is as follows: based on the original training samples used to establish the second mapping relationship set G in step S1, the distinguishing contribution of each environmental parameter among the four categories of samples (healthy, mildly stressed, moderately stressed, and severely stressed) is calculated, and then normalized by dividing the distinguishing contribution of each environmental parameter by the sum of the distinguishing contributions of all environmental parameters. For example, when the distinguishing contributions of air humidity, leaf wetness duration, temperature, and soil moisture content are 0.35, 0.30, 0.20, and 0.15, respectively, the environmental parameter weights for these four environmental parameters are 0.35, 0.30, 0.20, and 0.15, respectively. Subsequently, at each window location, the absolute value of the correlation coefficients of the four environmental parameters is taken, multiplied by the corresponding environmental parameter weight, and the products are summed to obtain the local correlation value at that window location. The reason for using absolute values ​​is that different environmental parameters may have the same or opposite effects on the variation sequence of the initial stress value, but both reflect the existence of environmental lag effects.

[0055] After obtaining the local correlation values ​​at each window position, the local correlation values ​​are weighted and accumulated to obtain the correlation value corresponding to the current shift. The window weights used for the weighted accumulation are determined based on the fluctuation intensity of environmental parameters and the variation amplitude of the initial stress value within each window. Specifically, firstly, the fluctuation intensity of the standard environmental parameter subsequence within each window position is calculated. The fluctuation intensity is the sum of the average absolute values ​​of the differences between adjacent time series points of the four environmental parameters within the window, multiplied by the corresponding environmental parameter weights. Secondly, the variation amplitude of the initial stress value of the effective stress subsequence within the window position is calculated. The variation amplitude of the initial stress value is the difference between the maximum and minimum values ​​within the window. Then, the fluctuation intensity of the environmental parameters is multiplied by the variation amplitude of the initial stress value to obtain the original window weight for that window position. This weight is then normalized by dividing the original window weight of each window position by the sum of the original window weights of all window positions to obtain the final window weight. Finally, the local correlation values ​​at each window position are multiplied by the corresponding final window weights and summed to obtain the correlation value corresponding to the current shift.

[0056] For example, if there are four window positions for a given shift, with local correlation values ​​of 0.42, 0.61, 0.58, and 0.33 respectively, and corresponding final window weights of 0.20, 0.35, 0.30, and 0.15, then the correlation value corresponding to the current shift can be obtained by multiplying 0.42 by 0.20, 0.61 by 0.35, 0.58 by 0.30, and 0.33 by 0.15, and then summing them. This method can increase the contribution of window positions with significant fluctuations in environmental parameters and significant changes in initial stress values ​​to the overall judgment, while reducing the interference of stationary and random noise sections on lag correlation analysis.

[0057] After performing the above calculations on each shift within the preset lag time window, correlation values ​​corresponding to multiple shifts are obtained, and the maximum value among these correlation values ​​is selected as the lag correlation coefficient ρ for the corresponding spatial location. The lag correlation coefficient ρ characterizes the strength of the most significant lag association between the stress initial value variation sequence and the environmental parameter sequence at that spatial location. The larger ρ is, the more likely the change in the stress initial value at that spatial location is mainly driven by the time-series lag effect of environmental parameter changes preceding stress characterization changes. A preset threshold ρth is used to determine whether a spatial location is significantly affected by environmental lag. The preset threshold ρth is constructed as follows: From the labeled samples obtained in step S1, sample areas confirmed to be free of pests and diseases but experiencing morning fog, dew, or short-term high humidity disturbances are selected as environmental disturbance samples. Sample areas confirmed to have actual pests and diseases are selected as disease samples. The lag correlation coefficient ρ distribution for each spatial location in both types of samples is calculated. The median value between the 75th percentile of the lag correlation coefficient ρ in the environmental disturbance samples and the 25th percentile of the lag correlation coefficient ρ in the disease samples is then taken as the preset threshold ρth. This construction method allows the preset threshold ρth to balance the ability to suppress false environmental alarms and the ability to maintain the presence of actual diseases.

[0058] In a set of monitoring samples of leafy vegetable crops in a valley and hilly area, the preset threshold ρth obtained using the above method can be 0.65. When the lag correlation coefficient ρ of a certain spatial location is greater than 0.65, it indicates that the stress initial value variation sequence of that spatial location has a strong lag correlation with the environmental parameter sequence, and the spatial location should be marked; when the lag correlation coefficient ρ of a certain spatial location is not greater than 0.65, it is considered that the spatial location does not belong to the situation dominated by significant environmental lag, and can be retained for subsequent comprehensive early warning index calculation.

[0059] In this embodiment, step S5 is used to calculate the correction coefficient α based on the second mapping relationship set G and the environmental parameter ET of the Tth time series point, and to fuse the initial stress value PT of the Tth time series point with the correction coefficient α for the unmarked spatial location to obtain the comprehensive early warning index Q. At the same time, Q=0 is set for the marked spatial location, thereby realizing the joint judgment of environmental conditions and physiological stress signals.

[0060] Specifically, firstly, the T-th time-series environmental parameter ET is extracted from the unlabeled spatial locations in step S4. The T-th time-series environmental parameter ET includes at least the temperature, air humidity, leaf wetness duration, and soil moisture content of the corresponding spatial location. Then, the environmental parameters ET corresponding to each unlabeled spatial location are input into the second mapping relationship set G for probability mapping to obtain the environmental risk value corresponding to each unlabeled spatial location. The probability mapping is implemented as follows: first, the sampling period category to which the current spatial location belongs is determined; then, it is determined which set of environmental parameter intervals in the second mapping relationship set G the temperature, air humidity, leaf wetness duration, and soil moisture content fall into; if the input environmental parameter happens to fall into a stored environmental parameter interval combination in the second mapping relationship set G, the probability of pest and disease occurrence corresponding to that environmental parameter interval combination is directly read as the environmental risk value of that spatial location; if the input environmental parameter is between two adjacent environmental parameter interval combinations, the probability of pest and disease occurrence corresponding to the two or four adjacent environmental parameter interval combinations is interpolated and mapped according to the distance ratio of each environmental parameter to the boundary of the adjacent interval to obtain continuous environmental risk values. For example, during the morning period, if the temperature at an unmarked spatial location is 21°C, the air humidity is 92%, the leaf surface is wet for 4 hours, and the soil moisture content is 27%, and the probability of pest and disease occurrence for the corresponding adjacent environmental parameter interval combinations in the second mapping relationship set G is 0.68 and 0.76 respectively, then an environmental risk value of 0.72 can be obtained by proportional interpolation according to the relative position of the spatial location's environmental parameters falling between the two sets of intervals.

[0061] After obtaining the environmental risk value corresponding to each unmarked spatial location, the correction coefficient α corresponding to each unmarked spatial location is calculated based on the environmental risk value according to a preset nonlinear transformation rule.

[0062] The preset nonlinear transformation rule first maps the environmental risk value to an initial correction value through a piecewise nonlinear response function. The piecewise nonlinear response function is constructed by dividing the environmental risk value range into a low-risk range, a medium-risk range, and a high-risk range, where the low-risk range is defined as 0 to 0.30, the medium-risk range is defined as greater than 0.30 and not greater than 0.70, and the high-risk range is defined as greater than 0.70 and not greater than 1. Within the low-risk range, a slow-rising mapping method is used, causing the initial correction value to increase slightly with the environmental risk value. Specifically, it can be calculated by "starting with 0.60 and adding the environmental risk value multiplied by 0.50". Within the medium-risk range, an enhanced mapping method with a slope higher than that of the low-risk range is used, making the initial correction value more sensitive to changes in the environmental risk value. Specifically, it can be calculated by "based on the initial correction value corresponding to the end of the low-risk range, plus the portion of the environmental risk value exceeding 0.30 multiplied by 0.90". Within the high-risk range, a saturated stabilizing mapping method is used to avoid the correction coefficient α from being infinitely amplified due to extreme high humidity environments. Specifically, it can be calculated by "based on the initial correction value corresponding to the end of the medium-risk range, plus the portion of the environmental risk value exceeding 0.70 multiplied by 0.40". For example, when the environmental risk value of an unmarked spatial location is 0.20, its initial correction value can be calculated as "0.60 plus 0.20 multiplied by 0.50", resulting in 0.70; when the environmental risk value is 0.60, its initial correction value can be calculated as "the initial correction value at the end of the low-risk interval plus 0.30 multiplied by 0.90"; when the environmental risk value is 0.85, an increment of 0.15 multiplied by 0.40 is added to the initial correction value at the end of the medium-risk interval. The purpose of using the above piecewise nonlinear response function is to make the adjustment effect of moderate-intensity environmental risk changes on the correction coefficient α more obvious, while keeping the impact of extremely low or extremely high environmental risk values ​​on the correction coefficient α relatively mild.

[0063] After obtaining the initial correction value, the initial correction value is modulated with environmental sensitivity based on the temperature, air humidity, and leaf wetness duration corresponding to each unmarked spatial location at the current time point to obtain an intermediate correction value. The environmental sensitivity modulation is implemented as follows: sensitivity factors are first constructed for temperature, air humidity, and leaf wetness duration. The temperature sensitivity factor is determined according to the suitable temperature range for disease growth: when the temperature is between 18 and 24 degrees Celsius, the temperature sensitivity factor is 1.15; when the temperature is between 12 and 18 degrees Celsius or between 24 and 28 degrees Celsius, the temperature sensitivity factor is 1.00; and when the temperature is below 12 degrees Celsius or above 28 degrees Celsius, the temperature sensitivity factor is 0.85. The sensitivity factor for air humidity was determined based on the degree to which a moist foliage environment promotes pathogen transmission. When the air humidity was greater than or equal to 90%, the sensitivity factor was 1.20; when the air humidity was between 75% and 90%, the sensitivity factor was 1.00; and when the air humidity was less than 75%, the sensitivity factor was 0.80. The sensitivity factor for the duration of foliage moisture was determined based on the supporting effect of continuous moisture on spore germination. When the duration of foliage moisture was greater than or equal to 4 hours, the sensitivity factor was 1.20; when the duration of foliage moisture was between 2 and 4 hours, the sensitivity factor was 1.00; and when the duration of foliage moisture was less than 2 hours, the sensitivity factor was 0.85. Next, the sensitivity factors for temperature, air humidity, and leaf wetness duration are multiplied by their corresponding weights and then summed to form a comprehensive modulation factor. The weights for temperature, air humidity, and leaf wetness duration can be set to 0.25, 0.35, and 0.40, respectively. The initial correction value is then multiplied by this comprehensive modulation factor to obtain an intermediate correction value. For example, if the temperature at an unmarked spatial location at time point T is 20°C, the air humidity is 93%, and the leaf wetness duration is 5 hours, then the sensitivity factors for temperature, air humidity, and leaf wetness duration are 1.15, 1.20, and 1.20, respectively. The comprehensive modulation factor can be obtained by multiplying 1.15 by 0.25, 1.20 by 0.35, and 1.20 by 0.40, and then summing the results. This summation is then used to modulate the initial correction value.

[0064] After obtaining the intermediate correction value, a time-series smoothing constraint is applied to the intermediate correction value based on the environmental risk variation amplitude of each unmarked spatial location at multiple consecutive time series points to obtain a stable correction value. The environmental risk variation amplitude is used to characterize the degree of fluctuation of the environmental risk value of the unmarked spatial location at multiple consecutive time series points. It is calculated as follows: extract the environmental risk value of the unmarked spatial location at the T-th time series point and the two time series points before it, calculate the absolute value of the difference between the environmental risk values ​​between two adjacent time series points, and then average the absolute values. If the change in environmental risk is no greater than 0.10, the current environmental change is considered relatively stable, and the stability correction value is directly taken as the intermediate correction value. If the change in environmental risk is greater than 0.10 but no greater than 0.25, the environmental change is considered to have some fluctuations, and the stability correction value is obtained by multiplying the intermediate correction value by 0.70 and adding the correction coefficient α of the previous time point by 0.30. If the change in environmental risk is greater than 0.25, the environmental change is considered relatively drastic. To suppress excessive correction caused by sudden fluctuations at a single moment, the stability correction value is obtained by multiplying the intermediate correction value by 0.50, adding the correction coefficient α of the previous time point by 0.30, and adding the correction coefficient α of the previous two time points by 0.20. For example, if the environmental risk values ​​of an unmarked spatial location at time point T, the previous time point, and the previous two time points are 0.80, 0.58, and 0.55, respectively, then the magnitude of the environmental risk change is the average of the absolute values ​​of the differences between 0.80 and 0.58 and between 0.58 and 0.55, which is 0.135. In this case, a time-series smoothing constraint method corresponding to moderate fluctuation intensity is used to calculate the stability correction value. This step can effectively avoid abnormal amplification of the environmental risk value at a single time point due to sudden enhancement of morning fog, short-term sprinkler irrigation, or local heat flow.

[0065] After obtaining a stable correction value, the stable correction value is truncated at upper and lower limits and subjected to monotonic consistency correction. The correction coefficient α corresponding to each unmarked spatial location is output, and the correction coefficient is constrained within a preset numerical range. The preset numerical range is preferably defined as 0.50 to 1.50. When the stable correction value is less than 0.50, the correction coefficient α is set to 0.50; when the stable correction value is greater than 1.50, the correction coefficient α is set to 1.50; when the stable correction value is between 0.50 and 1.50, this stable correction value is temporarily taken as the correction coefficient α. The monotonic consistency correction is implemented as follows: The unmarked spatial locations under the same time period category are checked after sorting the environmental risk values ​​from smallest to largest. If a situation arises where the environmental risk value is larger but the correction coefficient α is smaller, the latter's correction coefficient α is adjusted to be no less than the former's correction coefficient α. The adjustment range is 50% to 100% of the difference between the two, preferably 100%, to ensure that the correction coefficient α does not decrease as the environmental risk value increases. For example, if the environmental risk values ​​for two unmarked spatial locations are 0.62 and 0.70, respectively, and the corresponding stable correction values ​​are 1.05 and 1.02, then a monotonic consistency correction is applied to the latter, adjusting its correction coefficient α to 1.05. Through upper and lower limit truncation and monotonic consistency correction, situations where the correction coefficient α is excessively amplified, excessively compressed, or inconsistent with the direction of change in the environmental risk value can be avoided.

[0066] After obtaining the correction coefficient α corresponding to each unlabeled spatial location, the initial stress value PT corresponding to each unlabeled spatial location at the T-th time point is fused with the corresponding correction coefficient α to obtain the initial warning value corresponding to each unlabeled spatial location. Specifically, firstly, the initial stress value PT and the correction coefficient α corresponding to each unlabeled spatial location at the T-th time point are extracted, and the initial stress value PT and the correction coefficient α are normalized to obtain the fused input quantity. Since the initial stress value PT has been mapped to the 0 to 1 interval in step S2, the normalization of the initial stress value PT adopts the identity preservation method, that is, the current initial stress value PT is directly used as the normalized initial stress value. The normalization of the correction coefficient α is performed by linear conversion according to the preset numerical interval, specifically: the current correction coefficient α is subtracted from the lower limit of the preset numerical interval 0.50, and then divided by the difference between the upper limit of the preset numerical interval 1.50 and the lower limit 0.50, thereby converting the correction coefficient α to the 0 to 1 interval. For example, when the correction coefficient α for an unlabeled spatial location is 1.10, its normalized correction coefficient can be calculated by subtracting 0.50 from 1.10 and then dividing by 1.00, resulting in 0.60. This normalization process ensures that the initial stress value PT and the correction coefficient α participate in subsequent fusion calculations on the same numerical scale.

[0067] After obtaining the fusionable input, a nonlinear coupled fusion function is constructed based on the normalized initial stress value PT and the correction coefficient α to calculate the initial warning value corresponding to each unlabeled spatial location. The purpose of constructing the nonlinear coupled fusion function is to ensure that the initial stress value PT and the correction coefficient α are not only linearly added together, but also amplified through coupling terms the spatial locations of "high stress and high environmental risk". In specific implementation, the nonlinear coupled fusion function includes three parts: the first part is the linear term corresponding to the normalized initial stress value PT, the second part is the linear term corresponding to the normalized correction coefficient α, and the third part is the coupling term formed by multiplying the normalized initial stress value PT and the normalized correction coefficient α. The fusion weights corresponding to the three parts can be set to 0.55, 0.20, and 0.25, respectively, and the sum of the three fusion weights is 1. Therefore, the initial warning value can be calculated using the following textual logic: multiply the normalized initial stress value PT by 0.55, add the normalized correction coefficient α by 0.20, and then add the product of the normalized initial stress value PT and the normalized correction coefficient α by 0.25. For example, when the initial stress value PT corresponding to an unmarked spatial location at the T-th time point is 0.68 and the normalized correction coefficient α is 0.60, its initial warning value can be obtained by multiplying 0.68 by 0.55, adding 0.60 by 0.20, adding 0.68 by 0.60, and then multiplying by 0.25. This design ensures that the initial stress value PT dominates the overall judgment, while using the correction coefficient α to enhance or suppress environmental suitability conditions, and significantly improves the initial warning value through coupling terms when both are high.

[0068] After obtaining the initial warning values ​​corresponding to each unmarked spatial location, neighborhood consistency correction is performed based on the continuity of the distribution of the initial warning values ​​in adjacent spatial locations, and the comprehensive warning index Q corresponding to each unmarked spatial location is output. The neighborhood consistency correction uses a 3×3 spatial neighborhood. For each unmarked spatial location, the initial warning values ​​of its center location and its eight surrounding adjacent spatial locations are first extracted; then, adjacent spatial locations with a difference of no more than 0.20 from the initial warning value of the center location are selected as effective neighbors; if the number of effective neighbors is no less than four, the initial warning value of the center location and the initial warning values ​​of each effective neighbor are weighted and averaged according to the method of "the center location has a weight of 0.50, and the remaining effective neighbors are allocated the remaining 0.50 weights according to the inverse distance normalization" to obtain the comprehensive warning index Q after neighborhood consistency correction; if the number of effective neighbors is less than four, the initial warning value of the center location is kept unchanged and is used as the comprehensive warning index Q. The purpose of adopting the above judgment rule is to preserve the continuous spatial distribution characteristics of disease patches and suppress false high values ​​caused by isolated noise points. For example, if the initial warning value of an unlabeled spatial location is 0.72, and the initial warning values ​​of 5 of its 8 neighboring spatial locations are between 0.60 and 0.80, then a neighborhood consistency correction is performed on that spatial location. If only one of its neighboring spatial locations has an initial warning value close to 0.72, while the rest are below 0.30, then that location is considered an isolated anomaly, and no neighborhood diffusion correction is performed. Meanwhile, for spatial locations already labeled in step S4, they are no longer included in the above fusion calculation and neighborhood consistency correction; instead, Q=0 is directly set to explicitly exclude the interference of the environmental lag-dominant region on the final warning output.

[0069] In this embodiment, step S6 is used to compare the comprehensive early warning index Q with the early warning threshold Qth. When Q is greater than or equal to Qth, it is determined that pests and diseases have occurred at the corresponding sampling point and an early warning signal is output.

[0070] Specifically, the comprehensive early warning index Q corresponding to each spatial location within the target area is read point by point, and a unified judgment rule is used for threshold comparison. The early warning threshold Qth is pre-determined based on the sample data of the disease level labeling results completed in step S1. The construction method is as follows: select sample areas where diseases and pests have been confirmed to occur and sample areas where diseases and pests have not occurred, and calculate the corresponding comprehensive early warning index Q distribution for each; then use the boundary value that takes into account both the missed report rate and the false alarm rate as the early warning threshold Qth, preferably taking the comprehensive early warning index Q value corresponding to the maximum sum of the disease sample identification rate and the non-disease sample exclusion rate as the early warning threshold Qth. In a set of leafy vegetable monitoring samples in a river valley and hilly area, the early warning threshold Qth can be set to 0.62, that is, when the comprehensive early warning index Q corresponding to a certain sampling point is 0.62 or above, it is determined that the sampling point has diseases and pests; when the comprehensive early warning index Q corresponding to a certain sampling point is less than 0.62, it is determined that the sampling point has not yet met the disease and pest early warning conditions. To avoid misjudgment caused by a single isolated high value, it can be further stipulated that: only when the comprehensive early warning index Q of the target sampling point reaches the early warning threshold Qth, and the comprehensive early warning index Q of at least two sampling points in its adjacent spatial location simultaneously reaches the early warning threshold Qth, will the target sampling point be finally determined to have pests and diseases, thereby improving the stability of spatial determination.

[0071] After completing the threshold comparison and determining the sampling points where pests and diseases occur, an early warning signal corresponding to the sampling points is output. The early warning signal includes at least the location of the pest or disease occurrence, the early warning level, and the early warning time. The location of the pest or disease occurrence is determined by the spatial coordinates of the corresponding sampling point. The early warning time uses the collection time of the T-th time series point. The early warning level is determined based on the magnitude by which the comprehensive early warning index Q exceeds the early warning threshold Qth. Specifically, when the comprehensive early warning index Q is greater than or equal to the early warning threshold Qth and less than "early warning threshold Qth plus 0.10", a Level 1 early warning signal is output; when the comprehensive early warning index Q is greater than or equal to "early warning threshold Qth plus 0.10" and less than "early warning threshold Qth plus 0.20", a Level 2 early warning signal is output; and when the comprehensive early warning index Q is greater than or equal to "early warning threshold Qth plus 0.20", a Level 3 early warning signal is output. For example, when the comprehensive early warning index Q of a certain sampling point is 0.74 and the early warning threshold Qth is 0.62, since 0.74 is higher than 0.72 but lower than 0.82, the sampling point outputs a secondary early warning signal, and its spatial location and corresponding time are recorded in the early warning results. Finally, all sampling points where pests and diseases have occurred can be marked on the target area base map to form a pest and disease early warning distribution result, which can be used to guide growers to conduct targeted inspections and precise control in high-risk areas of pests and diseases.

[0072] 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 early warning of vegetable diseases and pests integrating low-altitude remote sensing and environmental parameter analysis, characterized by: include: S1, obtain the first mapping relationship set F between vegetation index and the degree of pest and disease stress, and the second mapping relationship set G between environmental parameters and the probability of pest and disease occurrence; S2, acquire multispectral images of the target area at T consecutive time points, extract the spectral reflectance of each sampling point, and calculate the initial stress value of each time point based on the first mapping relationship set F to obtain the initial stress value time series matrix P; S3, synchronously collect environmental parameters for the T consecutive time points to obtain an environmental parameter list set E; S4. Based on the stress initial value time series matrix P and the environmental parameter list set E, calculate the lag correlation coefficient ρ between the stress initial value variation sequence and the environmental parameter sequence at each spatial location, and mark the corresponding spatial location when ρ is greater than the preset threshold ρth. S5. Based on the second mapping relationship set G and the environmental parameter ET of the Tth time series point, calculate the correction coefficient α. For unmarked spatial locations, fuse the initial stress value PT of the Tth time series point with the correction coefficient α to obtain the comprehensive early warning index Q. For marked spatial locations, set Q=0. S6 compares the comprehensive early warning index Q with the early warning threshold Qth. When Q is greater than or equal to Qth, it is determined that pests and diseases have occurred at the corresponding sampling point and an early warning signal is output.

2. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 1, characterized in that: Step S1 includes: selecting vegetable sample areas in healthy, mildly stressed, moderately stressed, and severely stressed states, and spatially locating each sample area; based on the spatial location, simultaneously collecting multispectral images, temperature, air humidity, leaf surface wetness duration, and soil moisture content data for each sample area, and obtaining corresponding disease level labeling results; extracting vegetation index features based on the multispectral images, and establishing a correspondence between vegetation indices and the degree of disease and pest stress based on the disease level labeling results; and establishing a correspondence between environmental parameters and the probability of disease and pest occurrence using a time-segmented fitting method based on the environmental parameter data, historical disease occurrence records, and corresponding disease level labeling results, thereby generating the first mapping relationship set F and the second mapping relationship set G respectively.

3. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 1, characterized in that: The step of calculating the initial stress values ​​for each time series point based on the first mapping relationship set F to obtain the initial stress value time series matrix P includes: repeatedly conducting aerial surveys of the target area at preset time intervals to acquire multispectral images of each time series point; sequentially performing radiometric calibration, geometric correction, and spatial registration on the multispectral images of each time series point to obtain standardized images under a unified coordinate system; extracting the spectral reflectance of each sampling point in each band based on the standardized images and calculating the corresponding vegetation index; inputting the vegetation index of each sampling point into the first mapping relationship set, calculating the initial stress values ​​corresponding to each time series point, and arranging them according to time order and spatial location to generate the initial stress value time series matrix P.

4. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 3, characterized in that: The step of inputting the vegetation indices of each sampling point into the first mapping relationship set and calculating the initial stress value corresponding to each time series point includes: combining multiple vegetation indices corresponding to each sampling point according to preset weights to construct a vegetation index feature vector for the sampling point; performing interval normalization on the vegetation index feature vector and performing feature matching according to the preset stress classification interval in the first mapping relationship set to determine the candidate stress level corresponding to each sampling point; determining the upper and lower bound reference intervals adjacent to the candidate stress levels in the first mapping relationship set based on the candidate stress levels; extracting the interval position parameters of the vegetation index feature vector of the sampling point in the upper and lower bound reference intervals, and constructing a piecewise interpolation function accordingly; and using the piecewise interpolation function to perform continuous numerical mapping on the vegetation index feature vector of the sampling point to obtain the initial stress value for each time series point.

5. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 1, characterized in that: The calculation of the lag correlation coefficient ρ between the initial stress value variation sequence and the environmental parameter sequence at each spatial location includes: According to the spatial location, the initial stress values ​​of each sampling point at multiple consecutive time points are extracted from the initial stress value time series matrix, and the corresponding initial stress value variation sequence is constructed; based on the spatial location of the sampling points, the environmental parameter data of the corresponding time points in the environmental parameter list are matched to construct the environmental parameter sequence corresponding to the initial stress value variation sequence. The environmental parameter sequence is shifted hourly according to a preset lag time window, and the correlation value between the shifted environmental parameter sequence and the stress initial value variation sequence is calculated respectively; the maximum value among the correlation values ​​is selected as the lag correlation coefficient ρ of the corresponding spatial location.

6. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 5, characterized in that: The step of calculating the correlation between the shifted environmental parameter sequence and the initial stress value variation sequence includes: performing detrending and amplitude normalization processing on each shifted environmental parameter sequence to obtain a standard environmental parameter subsequence; extracting the effective stress subsequence for the corresponding time period based on the initial stress value variation sequence aligned with the standard environmental parameter subsequence; using the standard environmental parameter subsequence and the effective stress subsequence as input, obtaining the local correlation value at each window position using a sliding window correlation calculation method, and weighting and accumulating each local correlation value to obtain the correlation value corresponding to the current shift amount.

7. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 1, characterized in that: The step of calculating the correction coefficient α based on the second mapping relationship set G and the environmental parameter ET at the Tth time point includes: extracting the environmental parameter ET corresponding to each unmarked spatial location at the Tth time point, and inputting it into the second mapping relationship set for probability mapping to obtain the environmental risk value corresponding to each unmarked spatial location; and calculating the correction coefficient α corresponding to each unmarked spatial location based on the environmental risk value according to a preset nonlinear transformation rule.

8. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 7, characterized in that: The comprehensive early warning index Q is obtained by fusing the initial stress value PT at the T-th time point with the correction coefficient α. This includes: constraining the correction coefficient within a preset numerical range; fusing the initial stress value corresponding to each unmarked spatial location at the T-th time point with the corresponding correction coefficient to obtain the initial early warning value corresponding to each unmarked spatial location; and then performing neighborhood consistency correction based on the distribution continuity of the initial early warning value in adjacent spatial locations to output the comprehensive early warning index Q corresponding to each unmarked spatial location.

9. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 7, characterized in that: The correction coefficient α corresponding to each unmarked spatial location is calculated based on the environmental risk value according to a preset nonlinear transformation rule. This includes: inputting the environmental risk value corresponding to each unmarked spatial location into a preset piecewise nonlinear response function to obtain the corresponding initial correction value; modulating the initial correction value based on the temperature, air humidity, and leaf wetness duration corresponding to each unmarked spatial location at the current time point to obtain an intermediate correction value; then, according to the environmental risk change amplitude of each unmarked spatial location at multiple consecutive time points, applying a time-series smoothing constraint to the intermediate correction value to obtain a stable correction value; finally, performing upper and lower limit truncation and monotonic consistency correction on the stable correction value to output the correction coefficient α corresponding to each unmarked spatial location.

10. The vegetable pest and disease early warning method integrating low-altitude remote sensing and environmental parameter analysis according to claim 8, characterized in that: The step of fusing the initial stress value and the corresponding correction coefficient corresponding to each unmarked spatial location at the Tth time point includes: extracting the initial stress value and correction coefficient corresponding to each unmarked spatial location at the Tth time point, and normalizing the initial stress value and the correction coefficient to obtain a fusionable input; constructing a nonlinear coupled fusion function based on the normalized initial stress value and correction coefficient, and calculating the initial warning value corresponding to each unmarked spatial location.

Images