Intelligent recognition and early warning system for forest diseases and insect pests based on deep learning

Abstract

This invention discloses a deep learning-based intelligent identification and early warning system for forest pests and diseases, relating to the field of forest pest and disease monitoring technology. It constructs a forest observation matrix by acquiring multimodal data from forest areas, extracts suspected pest and disease feature maps, and constructs a forest risk association map by combining spatial clustering relationships. Furthermore, it utilizes a graph-constrained time-series prediction model to generate a lesion evolution probability sequence and an expansion direction vector. Based on this, it calculates a graded early warning index and outputs the early warning results by comprehensively considering the risk growth rate, cumulative risk value, and directional consistency index. This invention enables early identification and dynamic prediction of forest pests and diseases, improves the accuracy and timeliness of early warnings, and provides a basis for prevention and control decisions.

Description

Technical Field

[0001] This invention relates to the field of forest pest and disease monitoring technology, specifically to a deep learning-based intelligent identification and early warning system for forest pests and diseases. Background Technology

[0002] Current methods for identifying forest pests and diseases largely rely on visible light image classification or manual inspections. These methods typically focus on obvious symptoms such as large leaf spots and overall tree discoloration, making it difficult to identify fine-grained early anomalies formed in shady areas at forest edges, in foggy and humid regions, and during the incubation period of insect eggs. Especially during the early stages of pine wilt nematode transmission, the early stages of longhorn beetle borer infestation, and the stages of localized fungal infection, the canopy texture shows less variation, temperature differences are small, and spatial distribution is scattered. Trees are easily affected by wind, branch shadows, and topographic microclimate disturbances, leading to misjudgments of small lesions and difficulty in depicting cross-temporal expansion trends. Existing methods are insufficient for early identification and accurate warning of forest pests and diseases. Summary of the Invention

[0003] The purpose of this invention is to provide a deep learning-based intelligent identification and early warning system for forest pests and diseases, in order to address the shortcomings of the prior art.

[0004] To achieve the above objectives, the present invention provides the following technical solution: a deep learning-based intelligent identification and early warning system for forest pests and diseases, comprising: The acquisition and registration module acquires visible light canopy images, thermal infrared temperature maps, and corresponding environmental parameters of the target forest area, performs spatiotemporal registration, and generates a forest observation matrix. Feature extraction and fusion module: Input the forest observation matrix into a dual-branch deep network to extract canopy texture features, temperature difference anomalous features and humidity-heat coupling features, and fuse them to generate a suspected feature map of pests and diseases; Risk association modeling module: Based on the spatial clustering relationship of the suspected pest and disease feature map, construct a forest risk association map and extract the diffusion edge weights between adjacent sample trees; Temporal evolution prediction module: Input the diffusion edge weights and the suspected feature map of pests and diseases into the temporal prediction model to generate a lesion evolution probability sequence and an expansion direction vector; Early warning output module: Based on the lesion evolution probability sequence and expansion direction vector, calculate the graded early warning index of the target forest area, and output early warning results including pest and disease type, risk level, spread path and treatment area.

[0005] Preferably, the process of generating the forest observation matrix includes: Radiometric correction and distortion correction were performed on visible light canopy images to extract standardized canopy texture images; The thermal infrared temperature map is subjected to noise suppression and temperature field reconstruction to generate a spatially continuous temperature distribution map; A time synchronization sequence is constructed based on the corresponding location environmental parameters, and the standardized canopy texture image and the temperature distribution map are timestamped. The geographic coordinate mapping model was used to spatially register the time-aligned multi-source data, and the data were fused at a unified grid scale to generate a forest observation matrix that includes canopy texture features, temperature field features and environmental parameter features.

[0006] Preferably, the process of generating suspected pest and disease feature maps includes: The forest observation matrix is ​​divided into a canopy texture sub-matrix and a temperature environment sub-matrix according to the feature channels, and then input into a structurally symmetric dual-branch convolutional neural network for feature encoding to obtain the canopy texture feature tensor and the temperature difference response feature tensor. Based on the temperature difference response feature tensor and environmental parameter features, a humid-thermal coupling mapping function is constructed, and a humid-thermal coupling feature tensor is generated by channel-by-channel weighted calculation. The canopy texture feature tensor and the humid-thermal coupling feature tensor are spatially aligned and channel-separated, and the responses of each channel are adaptively enhanced through an attention weight allocation strategy to obtain a fused feature tensor. Threshold mapping and spatial smoothing are performed on the fused feature tensor to output a suspected feature map of pests and diseases.

[0007] Preferably, the process of constructing the forest risk association map includes: Based on the suspected feature map of pests and diseases, suspected feature values ​​of each tree node are extracted, and an initial adjacency matrix is ​​constructed according to spatial proximity. The suspected feature values ​​are then assigned to each node as graph signals. A graph Fourier transform is performed on the graph signals to construct a graph Laplace operator and optimize the adjacency weights by minimizing the graph signal smoothness index, resulting in a first association weight matrix. A transmission cost matrix is ​​constructed based on the first association weight matrix and the node spatial coordinates. Using the normalized result of the suspected feature map of pests and diseases as the probability distribution, the optimal transmission model with entropy regularization is solved to obtain the probability transmission matrix between nodes. The probability transmission matrix is ​​then weighted and fused with the first association weight matrix to generate a forest risk association map.

[0008] Preferably, the weights between adjacent sample nodes in the fusion result are used as diffusion edge weights.

[0009] Preferably, the process of inputting the diffusion edge weights and the suspected pest and disease feature map into the time series prediction model to generate the lesion evolution probability sequence and the expansion direction vector includes: The suspected feature maps of pests and diseases are stacked sequentially according to a preset time window, and the risk evolution tensor of the sample tree node is constructed by combining the diffusion edge weights. The risk evolution tensor of the sample tree node is input into the graph-constrained time series prediction model, and the historical risk propagation intensity between adjacent sample tree nodes is recursively encoded using diffusion edge weights to obtain the node time series evolution characteristics. Based on the temporal evolution characteristics of the nodes and the diffusion edge weights, a directional response decoding function is constructed to calculate the lesion diffusion probability value and corresponding directional component of each wood node in future consecutive moments; The lesion diffusion probability value and the directional component are subjected to time consistency correction and spatial vector synthesis to generate a lesion evolution probability sequence and an expansion direction vector.

[0010] Preferably, the process of calculating the graded early warning index of the target forest area includes: Based on the lesion evolution probability sequence, the risk growth rate and cumulative risk value of each sample tree node in future continuous time are extracted, and the consistency index of risk propagation direction is calculated in combination with the expansion direction vector. A multidimensional risk assessment function is constructed based on the risk growth rate, cumulative risk value and risk propagation direction consistency index. Different risk factors are weighted and fused by setting weight coefficients to obtain the comprehensive risk score of the sample tree node. Based on the comprehensive risk score of the sample tree node and its spatial distribution characteristics, the quantile division method is used to determine the multi-level early warning threshold interval, and the risk level of the target forest area is divided. According to the spatial continuous area corresponding to each risk level and the dominant direction of the expansion direction vector, the graded early warning index of the target forest area is calculated.

[0011] Preferably, the graded early warning index is denoted as S, and the value of S ranges from 0 to 1. The final early warning level is divided according to the graded early warning index: when 0≤S<0.30, it is defined as a low early warning level; when 0.30≤S<0.55, it is defined as a medium early warning level; when 0.55≤S<0.75, it is defined as a high early warning level; and when 0.75≤S≤1, it is defined as an emergency early warning level.

[0012] The technical effects and advantages provided by the present invention in the above technical solution are as follows: 1. This invention deeply integrates suspected pest and disease feature maps, diffusion edge weights, and temporal evolution features to construct a multi-dimensional risk assessment and propagation modeling method for complex forest environments. Compared to traditional methods that rely solely on single-moment images or simple statistical indicators, this invention introduces graph-constrained temporal prediction and optimal transmission and diffusion modeling, which can simultaneously characterize the spatial aggregation characteristics and temporal evolutionary trends of pests and diseases. This effectively solves the problems of strong concealment of early lesions and unclear diffusion paths, thereby significantly improving the foresight and accuracy of pest and disease identification in forest areas, especially achieving stable identification even in the stage of weak anomalies.

[0013] 2. This invention constructs a multi-dimensional fusion assessment system integrating risk growth rate, cumulative risk value, and directional consistency indicators. It also calculates a graded early warning index by combining spatially continuous regions and expansion direction vectors. This enables the early warning results to not only quantitatively express risk but also reflect the direction of diffusion and the potential scope of impact. Compared to existing technical solutions that only output risk levels, this invention can simultaneously provide diffusion paths and treatment areas, achieving a shift from "identification" to "decision support." This improves the targeting and timeliness of forest pest and disease control, reduces the cost of manual inspections, and enhances adaptability to complex terrain and diverse climatic conditions. Attached Figure Description

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

[0015] Figure 1 This is a flowchart of the module of the intelligent identification and early warning system for forest pests and diseases based on deep learning, according to the present invention.

[0016] Figure 2 This is a flowchart illustrating the method for generating the lesion evolution probability sequence and expansion direction vector of the present invention. Detailed Implementation

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

[0018] For examples, please refer to Figure 1 and Figure 2 As shown in this embodiment, the intelligent identification and early warning system for forest pests and diseases based on deep learning includes: Acquisition and Registration Module: Acquires visible light canopy images, thermal infrared temperature maps, and corresponding location environmental parameters of the target forest area, performs spatiotemporal registration, and generates a forest observation matrix.

[0019] In one embodiment, after acquiring the visible canopy image, radiometric correction is performed first. Specifically, three calibration plates with known reflectivities (0.18, 0.50, and 0.90) are pre-positioned within the imaging area. During each acquisition of a visible canopy image, the three calibration plates are simultaneously photographed at the edge of the image, and the original grayscale values ​​of the red, green, and blue channels are read. The gain and bias coefficients of each color channel are then calculated using the least squares method, converting the original grayscale values ​​into apparent reflectivity. The conversion relationship is: Apparent reflectivity = Gain coefficient × Original grayscale value + Bias coefficient. After radiometric correction, visible canopy images acquired at different times have a unified brightness reference, eliminating brightness drift caused by changes in solar altitude angle and thin cloud cover.

[0020] After radiometric correction, distortion correction is performed on the visible light canopy image. A joint model of radial and tangential distortion is used for distortion correction. Let the coordinates of the distorted pixels be... The corrected pixel coordinates are The principal point coordinates of the image are The radial distance is The radial correction term is determined by the first, second, and third order radial distortion coefficients, while the tangential correction term is determined by two tangential distortion coefficients. These coefficients are obtained using multiple sets of checkerboard calibration images, and then all pixels are mapped to the correction coordinates. Fifteen sets of calibration images are used, with shooting angles covering frontal, left-tilt, right-tilt, top-down, and bottom-up views to ensure the stability of the obtained distortion parameters.

[0021] After distortion correction, the forest canopy region is separated from the visible light canopy image. Specifically, the vegetation response value is first calculated. Here, R, G, and B represent the apparent reflectance of the red, green, and blue channels, respectively. Then, the segmentation threshold is obtained using the maximum inter-class variance thresholding method. All candidate thresholds are traversed, and the inter-class variance between the target class and the background class is calculated. The threshold with the largest inter-class variance is taken as the canopy segmentation threshold. Then, an opening operation and a closing operation are performed using a 5×5 pixel structuring element to remove isolated noise and fill small holes. Finally, connected regions with an area of ​​less than 300 pixels are deleted, and only the main canopy region is retained.

[0022] After canopy separation, the canopy region is standardized. Specifically, the image is cropped using the canopy's bounding rectangle, resampled to 512×512 pixels, and the pixel values ​​are normalized to the 0-1 range to obtain a standardized canopy texture image. This standardized canopy texture image serves as the temporal alignment object and spatial registration reference in subsequent steps.

[0023] In one embodiment, the thermal infrared temperature map is first calibrated. Specifically, before acquisition, calibration is performed using two blackbody temperature points, set to 20 degrees Celsius and 45 degrees Celsius respectively. The raw response values ​​of the thermal infrared camera at these two temperature points are recorded. A temperature conversion coefficient is calculated using a linear relationship, converting the raw thermal infrared response values ​​into absolute temperature values. To reduce temperature errors caused by differences in tree species, the emissivity is set according to the dominant tree species in the target forest area; 0.96 is used for coniferous forests, and 0.98 for broadleaf forests; when the proportion of dominant tree species is less than 0.70%, a uniform value of 0.97 is used.

[0024] After temperature calibration, the thermal infrared temperature map is cropped to the same range as the canopy area corresponding to the standardized canopy texture image, ensuring that the thermal infrared temperature map and the standardized canopy texture image describe the same batch of tree canopies. Noise suppression is then performed. Noise suppression consists of two parts: temporal domain suppression and spatial domain suppression. Temporal domain suppression uses the median temperature method across three consecutive frames, taking the median temperature value for the same pixel location across three consecutive frames to eliminate instantaneous thermal noise. Spatial domain suppression uses bilateral filtering, with a spatial weight standard deviation of 3 pixels and a temperature weight standard deviation of 0.8 degrees Celsius, thereby removing local shot noise while preserving temperature boundaries.

[0025] After noise suppression, anomalous pixels are identified and temperature field reconstruction is performed. An anomalous pixel is defined as a pixel whose average temperature difference with its eight-yard neighborhood is greater than 3.0 degrees Celsius, and this difference appears in three consecutive frames. For anomalous and missing pixels, inverse distance-weighted interpolation is used for compensation. The interpolation window is 7×7 pixels, and the weights are defined as follows: ,in Let be the Euclidean distance from the pixel to be compensated to the i-th valid pixel. The temperatures of all valid pixels within the window are summed according to the aforementioned weights to obtain the compensated temperature value. After compensation, a spatially continuous temperature distribution map is formed. This temperature distribution map, along with the standardized canopy texture image, undergoes time alignment processing in subsequent steps.

[0026] In one embodiment, the corresponding location environmental parameters include ambient temperature, relative humidity, wind speed, and rainfall. Ambient temperature and relative humidity are obtained from a forest meteorological data acquisition device, wind speed is the average wind speed at a height of 2 meters above the ground, and rainfall is the cumulative rainfall over the past 10 minutes. First, the four types of environmental parameters are uniformly resampled to a 60-second time interval to form the original sequence of environmental parameters. If a certain environmental parameter is missing for no more than two sampling intervals, linear interpolation is used to complete it; if the missing duration exceeds two sampling intervals, the standardized canopy texture image and temperature distribution map corresponding to that time period are discarded.

[0027] After generating the complete original sequence of environmental parameters, these parameters were normalized. Normalization employed a sliding time window standardization method with a window length of 24 hours. Using the mean and standard deviation within the time window as benchmarks, environmental temperature, relative humidity, wind speed, and rainfall were converted into dimensionless values. This process reduces the impact of the diurnal cycle on different batches of image data.

[0028] When aligning timestamps, the acquisition time of the standardized canopy texture image is used as the reference time. Search for the closest temperature acquisition time in the temperature distribution map time series, and search for the surrounding area in the original environmental parameter series. The two sampling times before and after the standardization canopy texture image are used. If the time difference between the standardized canopy texture image and the temperature distribution map is no greater than 1 second, they are directly paired; if the time difference is greater than 1 second but no greater than 30 seconds, linear interpolation is used to generate the reference time. The corresponding temperature distribution map; if the time difference is greater than 30 seconds, the data set is discarded. Environmental parameters at the reference time... The values ​​at each location are also obtained using linear interpolation. After the above processing, a time-synchronized sequence is obtained, consisting of a standardized canopy texture image, a temperature distribution map, and environmental parameters at the same time. This time-synchronized sequence serves as the input for spatial registration and feature fusion in the next step.

[0029] In one implementation, the geographic coordinate mapping model is constructed using a quadratic polynomial coordinate mapping method. Specifically, at least six ground control points are selected within the target forest area. The planar coordinates of each ground control point are obtained using a differential positioning method, and the corresponding pixel positions are manually marked on a standardized canopy texture image and temperature distribution map. Let the image pixel coordinates be (x, y) and the geographic planar coordinates be (X, Y), then the mapping relationship is written as: ; .in, For the quadratic polynomial coordinate mapping coefficients in the X direction, Let i = 1, 2, ..., 6 be the quadratic polynomial coordinate mapping coefficients in the Y direction. A system of equations is established using all ground control points, and the coefficients are solved using the least squares method. After obtaining the coefficients, the temperature distribution map is mapped to the coordinate reference of the standardized canopy texture image. Registration accuracy is evaluated using the root mean square error (RMSE), defined as the square root of the average of the squares of the mapping errors of all control points. When the RMS error exceeds 0.5 meters, control points are reselected and the mapping coefficients are resolved.

[0030] After spatial registration, the registered standardized canopy texture image and temperature distribution map are divided into a unified grid. The unified grid size is 0.25 m × 0.25 m. For each grid cell, canopy texture features, temperature field features, and environmental parameter features are extracted. Canopy texture features include gray-level mean, local contrast, and texture entropy. The gray-level mean is the average gray-level value of all pixels within the grid cell; the local contrast is the square root of the average of the squares of the differences between the gray-level values ​​and the gray-level mean of the grid cell; the texture entropy is calculated based on the gray-level probability distribution. Temperature field features include the cell average temperature, temperature gradient amplitude, and relative temperature anomalies. The temperature gradient amplitude is obtained by taking the square root of the sum of the squares of the horizontal and vertical temperature differences; the relative temperature anomaly is defined as the average temperature of the current grid cell minus the average temperature within a 2-meter radius centered on the center of that grid cell. Environmental parameter features include the normalized values ​​of ambient temperature, relative humidity, wind speed, and rainfall at the reference time.

[0031] Subsequently, the three canopy texture features, three temperature field features, and four environmental parameter features of each grid cell are concatenated in a fixed order to form a one-dimensional feature sequence of length 10. Then, according to the spatial arrangement order of the grid cells, all one-dimensional feature sequences are stacked row by row to obtain the forest observation matrix. Each row in the forest observation matrix corresponds to a grid cell, and each column corresponds to a predefined feature, which can be directly used as input data for the subsequent two-branch deep network.

[0032] Feature extraction and fusion module: The forest observation matrix is ​​input into a two-branch deep network to extract canopy texture features, temperature difference abnormal features and humidity coupling features, and fuse them to generate a suspected feature map of pests and diseases.

[0033] In one embodiment, the forest observation matrix is ​​composed of a feature sequence after uniform grid division, and its size is defined as a three-dimensional matrix with height H, width W, and 10 channels. The 10 channels are arranged in a fixed order. The first 3 channels are, in order, the grayscale mean, local contrast, and texture entropy, forming the canopy texture sub-matrix; the last 7 channels are, in order, the unit average temperature, temperature gradient amplitude, relative temperature anomaly, normalized ambient temperature, normalized ambient relative humidity, normalized wind speed, and normalized rainfall, forming the temperature environment sub-matrix. Since the normalized ambient temperature, normalized ambient relative humidity, normalized wind speed, and normalized rainfall correspond to the entire observation area at the same time, these four environmental parameters can be copied to all grid positions in the spatial dimension, thereby ensuring that the temperature environment sub-matrix and the canopy texture sub-matrix have the same spatial size.

[0034] The symmetrical bi-branch convolutional neural network consists of a canopy texture coding branch and a temperature environment coding branch. The hierarchical structure, kernel size, stride, and downsampling method of the two branches are consistent, and only the network parameters are trained independently. Each branch includes an input mapping layer, a first coding layer, a second coding layer, and a third coding layer. The input mapping layer employs a single 1×1 convolution operation, projecting the input channels uniformly into 16 feature channels. The first encoding layer consists of two convolution operations with a 3×3 kernel, a stride of 1, and padding of 1. After each convolution, batch normalization and linear rectified function activation are performed sequentially, resulting in 32 feature channels. The second encoding layer, based on the output of the first encoding layer, first performs a maximum value downsampling with a 2×2 window and a stride of 2, then performs two convolution operations with a 3×3 kernel, outputting 64 feature channels. The third encoding layer, based on the output of the second encoding layer, again performs a maximum value downsampling with a 2×2 window and a stride of 2, then performs two convolution operations with a 3×3 kernel, outputting 128 feature channels. The linear rectified function is defined as: outputting the original value when the input value is greater than 0, and outputting 0 when the input value is less than or equal to 0.

[0035] The input to the canopy texture encoding branch is the canopy texture submatrix. After the above three layers of encoding, the first-scale canopy texture feature tensor, the second-scale canopy texture feature tensor, and the third-scale canopy texture feature tensor are obtained respectively. These three tensors together constitute the multi-scale canopy texture feature tensor. The input to the temperature environment encoding branch is the temperature environment submatrix. Under the same hierarchical structure, the first-scale temperature environment feature tensor, the second-scale temperature environment feature tensor, and the third-scale temperature environment feature tensor are obtained. In order to further transform the temperature environment features into a more sensitive temperature difference abnormality expression for pests and diseases, local difference response calculation is performed on the temperature environment feature tensor at each scale. Specifically, for the temperature-related channels in each scale, the feature value of the current grid position is subtracted from the mean feature value of its surrounding 5×5 neighborhood to obtain the local temperature difference response value. Then, the local temperature difference response value is superimposed with the original temperature environment feature by channel and then compressed by a convolution with a kernel size of 1×1 to form the temperature difference response result at that scale. The temperature difference response results at the three scales together constitute the multi-scale temperature difference response feature tensor. Therefore, the canopy texture submatrix and temperature environment submatrix obtained in the previous step, after being encoded by a structurally symmetric bi-branch convolutional neural network, respectively form the multi-scale canopy texture feature tensor and the multi-scale temperature difference response feature tensor for the next step of humidity-heat coupling calculation.

[0036] In one embodiment, the heat-humidity coupling mapping function is used to characterize the amplifying effect of high temperature, high humidity, rainfall accumulation, and low wind speed conditions on the probability of pest and disease occurrence. Since different tree species, slope locations, and pest and disease types have varying sensitivities to environmental changes, the heat-humidity coupling mapping function is constructed using an "environmental parameter projection - weight generation - channel-by-channel weighting" approach.

[0037] First, scale matching is performed on the environmental parameter features. Specifically, normalized values ​​of ambient temperature, relative humidity, wind speed, and rainfall are spatially replicated to create an environmental parameter map. This map is then averaged and pooled according to the spatial dimensions of the first, second, and third scales, yielding environmental parameter feature maps at each scale. Next, a 1×1 convolution operation is performed on each scale's feature map, mapping the four environmental parameter channels to a number of projected channels matching the number of temperature difference response feature channels at that scale, thus obtaining the environmental parameter projection tensor.

[0038] Secondly, the thermo-humidity coupling weights are calculated based on the multi-scale temperature difference response feature tensor and the environmental parameter projection tensor. For any channel c and any grid position (i,j) at any scale, let the temperature difference response value at that position be A(c,i,j) and the environmental parameter projection value be B(c,i,j). Then, the thermo-humidity coupling mapping function is defined as follows: first, add A(c,i,j) and B(c,i,j), then add the trainable bias value of that channel, and then input the result into a compression function, outputting a weight value between 0 and 1. The compression function is a logarithmic probability compression function, calculated as "1 divided by 1 plus the negative exponent of the natural constant". To amplify the impact of humidity and rainfall on the early stages of lesion expansion, the normalized values ​​of relative humidity (0.35), rainfall (0.30), temperature (0.25), and wind speed (-0.10) were assigned initial weights during the environmental parameter projection process. The negative sign indicates that increased wind speed inhibits the probability of pest and disease expansion. These initial weights can be further adjusted during the training phase.

[0039] After obtaining the thermo-hygroscopic coupling weights, a weighted calculation is performed on each channel of the multi-scale temperature difference response feature tensor. For any scale, any channel, and any grid location, the temperature difference response value is multiplied by the corresponding thermo-hygroscopic coupling weight to obtain the thermo-hygroscopic coupling response value. Then, the thermo-hygroscopic coupling response values ​​of all channels are arranged in their original spatial order to form the thermo-hygroscopic coupling feature tensor at that scale. The thermo-hygroscopic coupling feature tensors at the three scales together constitute the complete thermo-hygroscopic coupling feature tensor. Thus, the multi-scale temperature difference response feature tensor obtained in the previous step is further transformed into a thermo-hygroscopic coupling feature tensor that includes environmental constraints, providing input for spatial alignment and fusion in the next step.

[0040] In one embodiment, since the multi-scale canopy texture feature tensor and the humid-thermal coupling feature tensor contain three resolutions—original scale, half-scale, and quarter-scale—spatial alignment is performed before fusion. Specifically, using the half-scale as the unified fusion scale, the original scale feature tensor is converted to the half-scale using average pooling with a window size of 2×2 and a stride of 2. The quarter-scale feature tensor is then enlarged to the half-scale using bilinear interpolation. Subsequently, a 1×1 convolution operation is performed on each scale-converted canopy texture feature tensor, uniformly compressing the number of channels to 64. The same convolution compression process is performed on each scale-converted humid-thermal coupling feature tensor, also resulting in 64 channels. The canopy texture features from the three sources are concatenated along the channel direction to obtain a 192-channel canopy texture aggregate tensor; the humid-thermal coupling features from the three sources are also concatenated along the channel direction to obtain a 192-channel humid-thermal coupling aggregate tensor; finally, both are concatenated along the channel direction to form an initial fusion tensor with 384 channels.

[0041] To ensure that channels in the initial fusion tensor that truly reflect signs of pests and diseases receive higher weights, an attention weight allocation strategy is applied to the initial fusion tensor. This strategy comprises two consecutive steps: channel weight calculation and spatial weight calculation. Specifically, the channel weight calculation is performed as follows: global average pooling and global max pooling are applied to each channel of the initial fusion tensor, resulting in two one-dimensional vectors of length 384. These two vectors are then added and input into two fully connected layers. The first layer compresses the dimension from 384 to 96, and the second layer restores it to 384. A linear rectified function is applied after the first layer, and a log-odds compression function is applied after the second layer, resulting in 384 channel weight values. Each channel weight value is multiplied by all spatial location feature values ​​of the corresponding channel to obtain a channel-weighted tensor. The specific steps for calculating spatial weights are as follows: calculate the average and maximum values ​​of the channel-weighted tensor along the channel dimension, stitch the two images together along the channel direction to form a 2-channel feature map, and then input it into a convolution operation with a kernel size of 7×7 and a stride of 1. After the convolution result is processed by the log-odds compression function, the spatial weight map is obtained. Finally, multiply the spatial weight map with the channel-weighted tensor position by position to form the enhanced fused feature tensor.

[0042] To ensure that the fused feature tensor retains both canopy texture boundaries and temperature anomaly regions, it undergoes a final convolution with a 3×3 kernel before output, along with batch normalization and linear rectified function activation, to ensure consistent feature distributions from different sources. The resulting fused feature tensor inherits both texture variation information from the multi-scale canopy texture feature tensor and environment-driven anomaly information from the thermo-hygroscopic coupling feature tensor, and is directly used as input for the next step of threshold mapping and spatial smoothing.

[0043] In one implementation, the fused feature tensor is first mapped to a single-channel pest and disease probability map. Specifically, a 1×1 convolution operation is performed on the fused feature tensor to compress the number of channels from the current value to 1, resulting in a single-channel response map. Then, a log-probability compression function is applied to the single-channel response map, causing each grid position to output a probability value between 0 and 1, forming an initial pest and disease probability map. The larger the value at each grid position in the initial pest and disease probability map, the higher the probability that the position belongs to a suspected pest and disease area.

[0044] After obtaining the initial pest and disease probability map, threshold mapping is performed. To balance the recall rate of weak lesions and the ability to suppress false positives, a dual-threshold hysteresis discrimination method is adopted for threshold mapping. The high threshold is obtained by traversing the validation samples. Specifically, candidate high thresholds are selected sequentially between 0.50 and 0.90 with a step size of 0.02. For each candidate high threshold, the precision and recall are calculated, and then the harmonic mean of the two is calculated. The candidate high threshold with the largest harmonic mean is taken as the final high threshold. The low threshold is defined as 0.75 times the final high threshold. During threshold mapping, grid positions in the initial pest and disease probability map that are greater than or equal to the high threshold are marked as strongly suspected positions, grid positions that are less than the low threshold are marked as non-suspected positions, and grid positions between the low and high thresholds are marked as pending positions. If a pending position is connected to any strongly suspected position within an eight-yard neighborhood, the pending position is retained as a suspected position; otherwise, it is classified as a non-suspected position. Through the above processing, extended regions with weak edges but adjacent to core lesions can be preserved.

[0045] After thresholding, spatial smoothing is performed on the suspected location distribution results to eliminate isolated noise and enhance spatial continuity. Spatial smoothing is achieved using anisotropic diffusion guided by canopy texture. Specifically: First, a canopy texture guidance map is constructed based on the gradient magnitudes of the multi-scale canopy texture feature tensors at a unified fusion scale. The gradient magnitudes are obtained by taking the square root of the sum of the squares of the horizontal and vertical differences. Then, the suspected location distribution results after thresholding are iteratively updated five times. In each iteration, the new value of the current grid position is determined by the difference between its four adjacent positions (upper, lower, left, and right), multiplied by the conduction coefficient corresponding to the canopy texture guidance map. The conduction coefficient is defined as an exponential function with the natural constant as the base and the square of the negative gradient ratio as the exponent, where the denominator constant in the gradient ratio is 0.15, and the step size for each iteration is 0.20. This setting reduces the diffusion rate at locations with obvious canopy texture boundaries and increases the diffusion rate at locations with gentle canopy texture changes, thus filling local fracture areas without damaging lesion boundaries. After the diffusion process is completed, isolated connected regions with an area of ​​less than 4 grid cells are deleted, and the remaining connected regions are retained as the final output to obtain a suspected feature map of pests and diseases with spatial continuity.

[0046] In one embodiment, the symmetrical bi-branch convolutional neural network completes parameter learning during the offline training phase. Training samples are obtained by labeling on-site pest and disease survey results with corresponding historical monitoring images, marking pest and disease locations as 1 and non-pest and disease locations as 0. During training, the forest observation matrix is ​​input into the network to obtain an initial pest and disease probability map, which is then compared with the manually labeled map for loss calculation. The loss function consists of pixel-level binary cross-entropy loss and region overlap loss. The pixel-level binary cross-entropy loss constrains the classification accuracy of individual grid locations, while the region overlap loss constrains the shape consistency of the overall lesion region. The weight ratio of the two is set to 0.6 to 0.4. Parameter updates employ an adaptive moment estimation optimization method, with an initial learning rate of 0.001, a batch size of 8, and 120 iterations. When the validation loss decreases by less than 0.001 in 10 consecutive training iterations, the learning rate is reduced to 0.5 of its original value. The above training method can simultaneously learn the bi-branch convolution parameters, the channel weight parameters in the wet-heat coupling mapping function, and the fully connected mapping parameters in the attention weight allocation strategy.

[0047] The high threshold in the dual-threshold hysteresis discrimination is determined independently based on validation samples after network training, while the low threshold is calculated based on the high threshold at a fixed ratio. The anisotropic diffusion processing guided by canopy texture is only performed after obtaining the initial pest and disease probability map and does not participate in the reverse update of network parameters, thereby ensuring clear boundaries between the training and post-processing processes and facilitating reproduction.

[0048] Risk association modeling module: Based on the spatial clustering relationship of the suspected pest and disease feature map, construct a forest risk association map and extract the diffusion edge weights between adjacent sample trees.

[0049] In one implementation, sample tree nodes are first extracted from the suspected pest and disease feature map. Specifically, the response value of each grid position in the suspected pest and disease feature map is recorded as a suspected pest and disease probability between 0 and 1. An adaptive threshold method is used to extract high-risk candidate regions. The adaptive threshold is defined as the mean of all grid response values ​​in the current suspected pest and disease feature map plus 0.8 times the standard deviation. Grid positions greater than or equal to the adaptive threshold are marked as candidate positions, and adjacent candidate positions are merged into candidate patches using an eight-code neighborhood connectivity rule. If the area of ​​a candidate patch is less than 3 grid cells, it is deleted to eliminate scattered noise. For each retained candidate patch, its geometric center coordinates are calculated, and the geometric center is matched with the nearest neighbor coordinates of the canopy center obtained in the previous stage. When the Euclidean distance between the geometric center and a certain canopy center is not greater than 1.2 times the average radius of the canopy, the candidate patch is determined to belong to that canopy, and a sample tree node is generated accordingly. If multiple candidate patches exist within the canopy of the same tree, the candidate patch with the largest average response value shall be taken as the disease and pest characterization area of ​​that sample tree node.

[0050] After the sample tree nodes are determined, the suspected feature values ​​of each sample tree node are extracted. Specifically, based on the pest and disease characterization area to which the sample tree node belongs, the average and maximum values ​​of the suspected pest and disease probability within the area, as well as the average gradient amplitude, are calculated. These are then weighted and summed with weights of 0.5, 0.3, and 0.2 to obtain the suspected feature value of that sample tree node. This weighting ensures that the regional average risk accounts for the majority, while local peak risk and boundary variation are used as supplementary parameters. Thus, each sample tree node obtains a corresponding suspected feature value, which is recorded as the node signal amplitude.

[0051] After obtaining the spatial coordinates and potential eigenvalues ​​of all sample tree nodes, an initial adjacency matrix is ​​constructed. Specifically, each sample tree node is treated as a node in the graph, and the planar Euclidean distance between any two sample tree nodes is calculated. For each sample tree node, only the four closest sample tree nodes are retained as candidate connection objects, thus forming local spatial proximity relationships. The reason for selecting four candidate connection objects is that the direct influence of a single tree at the canopy scale in a forest area is usually concentrated in its nearest neighbors; too many connections will introduce long-distance spurious correlations. If two sample tree nodes are candidate connection objects for each other, an undirected connection is established between them. For node pairs with established connections, initial weights are assigned according to the Gaussian distance decay function. The specific calculation formula is: the initial weight between nodes is equal to the function value with the natural constant as the base and the square of the distance between nodes divided by the square of the distance scale parameter, taking the negative value as the exponent. The distance scale parameter is the average distance of all connected nodes. For node pairs without established connections, the initial weight is recorded as 0. The initial adjacency matrix is ​​formed by the initial weights of all node pairs.

[0052] After constructing the initial adjacency matrix, the suspected feature values ​​of each sample tree node are arranged in node order to form a graph signal column vector. Each element in the graph signal column vector corresponds to the pest and disease risk amplitude of a sample tree node. This graph signal column vector serves as the input for the next step of graph Fourier transform and graph signal smoothness index calculation.

[0053] In one implementation, a degree matrix is ​​first constructed based on the initial adjacency matrix. This degree matrix is ​​a diagonal matrix, where the i-th diagonal element is equal to the sum of all weights in the i-th row of the initial adjacency matrix. Then, a graph Laplacian operator is constructed based on the degree matrix and the initial adjacency matrix, specifically defined as the degree matrix minus the initial adjacency matrix. To eliminate the dimensionality effect caused by inconsistent node connection numbers, the graph Laplacian operator is further normalized to obtain a normalized graph Laplacian operator. The normalized graph Laplacian operator is constructed by first raising each diagonal element of the degree matrix to the power of -1 / 2, and then multiplying it by the graph Laplacian operator on both sides. After this processing, sample nodes with higher and lower degrees are comparable in frequency domain analysis.

[0054] Subsequently, eigenvalue decomposition is performed on the normalized graph Laplace operator, yielding a set of eigenvalues ​​arranged in ascending order and a corresponding set of orthogonal eigenvectors. These orthogonal eigenvectors constitute the graph Fourier basis. Projecting the graph signal column vectors onto the graph Fourier basis yields the graph frequency domain coefficients; this process is the graph Fourier transform. If the graph frequency domain coefficients are mainly concentrated in the low-frequency part, it indicates that the changes in the suspected pest and disease feature values ​​of spatially adjacent sample tree nodes are gradual, suggesting that the current adjacency relationship is relatively consistent with the actual diffusion relationship. If the graph frequency domain coefficients have a large proportion in the high-frequency part, it indicates that the initial adjacency matrix cannot fully explain the spatial clustering changes in the suspected pest and disease feature map, and further adjustment of the adjacency weights is needed.

[0055] To achieve adjacency weight inversion, a graph signal smoothness index is defined. This index is a scalar value obtained by left-multiplying the graph signal column vector by the graph Laplacian operator and then right-multiplying it by the graph signal column vector. Its physical meaning is the weighted sum of the differences in graph signals between all connected node pairs. A smaller graph signal smoothness index indicates smoother node signals corresponding to connected edges, and the constructed graph structure better conforms to the assumption of continuity in the spatial diffusion of pests and diseases. Based on this, an adjacency weight optimization objective function is established, consisting of three parts: the first part is the graph signal smoothness index, used to ensure that interconnected sample tree nodes have similar pest and disease risk amplitudes; the second part is a difference constraint term between the adjacency weights and the initial weights, used to limit the optimization result from deviating from the fact of spatial proximity; and the third part is a total adjacency weight constraint term, used to avoid all edge weights increasing simultaneously. The weight coefficients of the first, second, and third parts are 1.0, 0.4, and 0.1, respectively.

[0056] In the solution process, the initial adjacency matrix is ​​used as the initial value, and the weights of each connecting edge are iteratively updated using the projective gradient descent method. In each iteration, the gradient of each connecting edge weight is first calculated according to the objective function, and then the edge weight is updated along the negative gradient direction. The updated edge weights are then constrained to be non-negative, while maintaining the symmetry of the adjacency matrix. To ensure convergence stability, the learning step size is set to 0.01, and the maximum number of iterations is set to 200. The calculation terminates when the change in the objective function between two consecutive iterations is less than 0.0001. After the above process, the inverted optimized adjacency matrix is ​​obtained, which is denoted as the first association weight matrix. The first association weight matrix retains the spatial proximity constraint and incorporates the smoothing pattern of the suspected pest and disease feature map in the graph frequency domain, which can more accurately represent the static association strength between adjacent sample tree nodes. The first association weight matrix serves as an important input for constructing the transmission cost matrix in the next step.

[0057] In one implementation, an optimal transmission model with entropy regularization is used to characterize the directionality and distribution ratio of pest and disease risk spreading from highly suspected sample tree nodes to adjacent sample tree nodes. First, a transmission cost matrix is ​​constructed. For any two sample tree nodes i and j, the transmission cost is obtained by a linear combination of three parts: the first part is the normalized result of the Euclidean distance between the nodes as a percentage of the maximum distance of all connected nodes; the second part is the reciprocal normalized result of the edge weights corresponding to the first association weight matrix, where a larger edge weight indicates a stronger risk association between nodes, thus resulting in a lower transmission cost; the third part is a slope consistency correction term, where if node j is located within a 60-degree angle between node j and node i in the downstream direction of the prevailing wind, the correction term is 0; otherwise, it is 0.15. The combined weights of the three parts are 0.5, 0.4, and 0.1, respectively. The transmission cost matrix is ​​constructed from the transmission costs of all node pairs. For node pairs not directly connected in the initial adjacency matrix, a large constant of 1.5 is assigned to the transmission cost to suppress direct transmission from non-nearest neighbor nodes.

[0058] Subsequently, a probability distribution is constructed. Specifically, the total risk value is obtained by summing the suspected feature values ​​of all sample nodes. Then, the risk quality ratio of each sample node is obtained by dividing its suspected feature value by the total risk value. The source probability distribution is constructed from the risk quality ratios of all nodes. The target probability distribution is constructed in the same way, but to reflect the potential for infection, a very small smoothing constant of 0.001 is added to the suspected feature value of each sample node before normalization to avoid some low-response sample nodes having a target probability of 0. Therefore, the optimal transmission model with entropy regularization is transformed into: finding a non-negative transmission matrix that minimizes the weighted sum of the total transmission cost and the entropy regularization term, while satisfying the edge constraints of the source and target probability distributions. The entropy regularization term is defined as the negative value of the sum of the product of each element in the transmission matrix and its logarithmic value, used to encourage a smooth distribution of transmission results and avoid sharp concentrations only on a few edges.

[0059] To solve for the optimal transmission model with entropy regularization, an alternating scaling iterative method is adopted. First, a kernel matrix is ​​constructed based on the transmission cost matrix and the regularization coefficient. Each element in the kernel matrix is ​​equal to the exponent of a function with the natural constant as the base and the corresponding transmission cost value divided by the regularization coefficient, taking the negative value as the exponent. The regularization coefficient is set to 0.05. Then, two scaling vectors are initialized, both initially as vectors of all 1s. In each iteration, the first scaling vector is updated using the source probability distribution and the kernel matrix, and the second scaling vector is updated using the target probability distribution and the kernel matrix. After the two sets of scaling vectors converge, the first scaling vector, the kernel matrix, and the second scaling vector are multiplied sequentially to obtain the final transmission matrix, i.e., the probability transmission matrix. To ensure the stability of the results, the maximum number of iterations is set to 300; iteration stops when the average difference between the elements of the transmission matrix obtained from two consecutive iterations is less than 0.00001. In the obtained probability transmission matrix, the element (i,j) represents the proportion of pest and disease risk transmission allocated from sample tree node i to sample tree node j. Because this matrix is ​​simultaneously constrained by spatial distance, the first association weight, and the direction correction term, it not only represents the magnitude of the risk but also the direction of risk diffusion. In the next step, the probability transfer matrix is ​​fused together with the first association weight matrix to generate the final forest risk association map.

[0060] In one implementation, to balance the static spatial correlation strength obtained from graph signal processing and the dynamic directional propagation strength obtained from optimal transmission, the first correlation weight matrix and the probability transfer matrix are normalized before weighted fusion. Specifically, all non-zero elements in the first correlation weight matrix are linearly normalized from minimum to maximum value, falling within the 0-1 range; similarly, all elements in the probability transfer matrix are linearly normalized from minimum to maximum value. After normalization, the fusion weight is defined as: the fusion weight equals 0.6 times the normalized first correlation weight plus 0.4 times the normalized probability transfer value. The weight of the first correlation weight matrix is ​​set to 0.6 to ensure that the actual spatial continuity between adjacent sample tree nodes dominates the final graph structure; the weight of the probability transfer matrix is ​​set to 0.4 to supplement the expression of the flow trend of pest and disease risk from high-risk nodes to low-risk nodes.

[0061] Considering the directionality of the probability transfer matrix and the symmetric nature of the first association weight matrix, directed edges are retained when constructing the forest risk association graph. For any two sample tree nodes i and j, if they are connected in the initial adjacency matrix, the fusion weight from i to j and the fusion weight from j to i are calculated, and these two weights can be different. If the fusion weight in a certain direction is less than 0.12, the diffusion trend in that direction is considered insufficient, and the directed edge is not retained in the forest risk association graph. The 0.12 threshold is determined through historical sample traversal. Specifically, candidate thresholds are selected sequentially with a step size of 0.01 between 0.05 and 0.20, and the consistency rate between the risk diffusion path and the measured lesion expansion direction under different thresholds is calculated. The candidate threshold with the highest consistency rate is taken as the final threshold. After threshold screening, the retained directed connections and their corresponding fusion weights together constitute the forest risk association graph.

[0062] Diffusion edge weights are extracted from the forest risk association graph. For any retained directed edge pointing from sample tree node i to sample tree node j, its fusion weight is directly defined as the diffusion edge weight between adjacent sample trees. To further enhance the physical meaning of the edge weight, it can be multiplied by the product of the suspected feature value of node i and the infection sensitivity factor of node j, where the infection sensitivity factor is obtained by weighting the normalized value of the relative humidity and the normalized value of the temperature gradient amplitude at the location of node j by 0.5 and 0.5 respectively. After this correction, the diffusion edge weight not only reflects the connection strength between nodes, but also the intensity of the upstream risk source and the downstream infection conditions. All diffusion edge weights are stored in the order of the start and end points of the edge, serving as the basic input for characterizing the lesion evolution probability sequence and the expansion direction vector in subsequent time series prediction.

[0063] In one implementation, the core function of the graph signal smoothness index is to map the spatial clustering patterns in the suspected pest and disease feature map onto the graph structure. If adjacent sample tree nodes show similar responses in the suspected pest and disease feature map, the graph signal smoothness index decreases, indicating that the connection edge more closely matches the actual propagation path. Therefore, the first association weight matrix obtained by minimizing the graph signal smoothness index is essentially a structured expression of the spatial continuity of the suspected pest and disease feature map. The optimal transmission model with entropy regularization further introduces a quantitative allocation relationship of "where the risk quality flows from and to," thus enabling the final forest risk association map to possess both correlation and directionality. A threshold of 0.12 is used to suppress low-confidence edges, a candidate region area threshold of 3 grid cells is used to remove small noise points, and the nearest neighbor number of 4 is used to limit the local connection radius. These parameters can all be repeatedly obtained from sample data, demonstrating a clear construction process and feasibility.

[0064] Temporal evolution prediction module: Input the diffusion edge weights and the suspected feature map of pests and diseases into the temporal prediction model to generate a lesion evolution probability sequence and an expansion direction vector.

[0065] In one implementation, suspected pest and disease feature images are first acquired at a fixed sampling period for multiple consecutive time points. The sampling period is 1 day, and the preset time window length is 7 consecutive time points, i.e., the previous 6 days and the current day constitute a complete time window. If more than one suspected pest and disease feature image is missing in any time window, the time window is discarded; if only one image is missing, it is filled in using suspected pest and disease feature images from adjacent time points using a linear interpolation method.

[0066] In the suspected pest and disease feature map at each time point, the risk response value of the canopy area where each sample tree node is located is read according to the aforementioned sample tree node extraction results. Specifically, a circular area corresponding to the average radius of the sample tree canopy is taken, centered on the canopy center corresponding to the sample tree node. The average probability of suspected pest and disease at all grid positions within this circular area is calculated, and this average is defined as the node risk value of the sample tree node at that time. The node risk values ​​of all sample tree nodes within the time window are arranged in chronological order to form a node risk time series. Subsequently, to ensure the comparability of batches of data from different forest areas, time-window normalization is performed on the node risk time series of each sample tree node. The normalization method is: the current node risk value minus the minimum value of the node within the current time window, then divided by the difference between the maximum and minimum values ​​plus 0.0001, so that the normalized node risk value falls within the range of 0 to 1.

[0067] After obtaining the normalized node risk time series, diffusion edge weights are introduced into the time window structure. Specifically, for any given time, the diffusion edge weights of all directed edges in the forest risk association graph are read, keeping the starting and ending order of the edges unchanged. Then, for each sample tree node, its weighted input value from adjacent upstream sample tree nodes is calculated. This weighted input value is defined as the sum of all diffusion edge weights pointing to that sample tree node multiplied by the node risk value of the upstream sample tree node at the same time. This weighted input value reflects the external propagation effect exerted by adjacent sample tree nodes on the target sample tree node at the current time. Then, the normalized node risk value of the current sample tree node, the weighted input value at the corresponding time, and the sum of the incoming and outgoing edge weights of that sample tree node are concatenated in a fixed order to form a node state vector for a single time moment. The node state vectors of all sample tree nodes at all times are stacked in the order of "time-node-feature" to form the sample tree node risk evolution tensor.

[0068] In one specific implementation, assuming a time window length of 7, a sample tree node count of N, and a state vector length of 4 for each node, the size of the sample tree node risk evolution tensor is 7×N×4. In this sample tree node risk evolution tensor, the first dimension represents the temporal order, the second dimension represents the sample tree node order, and the third dimension represents the risk state feature order. Thus, the spatial risk changes in the suspected pest and disease feature map are converted into a temporal structure input that can be used for subsequent recursive prediction, and the diffusion edge weights are explicitly written into the tensor as features constraining the propagation relationship.

[0069] In one embodiment, the graph-constrained temporal prediction model is constructed using a two-level cascaded structure of "graph constraint propagation aggregation layer - gated recursive coding layer". Graph constraint propagation aggregation is performed first, followed by temporal recursive coding, to ensure that the propagation relationship between adjacent sample nodes takes precedence in the temporal state update.

[0070] The input to the graph-constrained propagation aggregation layer is the total node state vector corresponding to the risk evolution tensor of the sample tree node at a certain moment. For any sample tree node, first read the propagation weights of all directed edges pointing to that node, then sum the weighted node state vectors of each upstream sample tree node to obtain the neighborhood propagation vector. To avoid excessively large edge weights leading to abnormally concentrated influence of individual upstream nodes, the propagation edge weights are first normalized. The normalization method is as follows: add all propagation edge weights pointing to the same target sample tree node to obtain the denominator, then divide each propagation edge weight by the denominator and add 0.0001. Then, use the normalized propagation edge weights as weighting coefficients to sum the node state vectors of the upstream sample tree nodes element by element. If a sample tree node has no upstream connecting edges, its neighborhood propagation vector is a zero vector. After completing the above processing, the node state vector of the sample tree node itself and the neighborhood propagation vector are concatenated along the feature dimension, and then subjected to a linear mapping and hyperbolic tangent function activation to obtain the propagation encoding vector of the sample tree node at the current moment. The output dimension of the linear mapping is set to 16 to enhance the coupling expression between features from different sources.

[0071] The gated recursive coding layer is used to extract the cumulative impact of historical risk propagation intensity along the time direction. For any tree node, its propagation coding vectors at each time step are sequentially input into the gated recursive unit according to the time window order. The gated recursive unit includes three computational stages: update gate, retention gate, and candidate state. The update gate is used to determine what proportion of the propagation coding vector at the current time step should be written into the new hidden state; the retention gate is used to determine how much information should be retained in the hidden state at the previous time step; and the candidate state is used to generate new candidate expressions based on the current input and historical states. The specific calculation process is as follows: first, the propagation coding vector at the current time step and the hidden state at the previous time step are multiplied by the trainable weight matrix and added together, then a bias is added and input into the log-probability compression function to obtain the update gate and the retention gate; then, the propagation coding vector at the current time step and the hidden state at the previous time step after being filtered by the retention gate are input together into the hyperbolic tangent function to obtain the candidate state; finally, the hidden state at the current time step is obtained by multiplying the update gate by the hidden state at the previous time step, adding 1 and subtracting the update gate, and then multiplying by the candidate state. The initial hidden state is a vector of all zeros, and the dimension of the hidden state is 32.

[0072] To enhance sensitivity to sudden lesion growth, historical risk propagation intensity is introduced into the gating recursive coding process. The historical risk propagation intensity is defined as follows: Let the neighborhood propagation vector of sample node i at time t be... The intensity of historical risk transmission is defined as ;in, This represents the historical risk propagation intensity of sample node i at time t, used to characterize the magnitude of risk propagation changes. Let i represent the neighborhood propagation vector of sample node i at time t-1. This represents the diffusion edge weight from sample node j to sample node i, reflecting the risk propagation strength from node j to node i; Let represent the set of all adjacent nodes pointing to sample node i, i.e., the in-neighbor set. If the historical risk propagation intensity is greater than 0.18, it indicates that the external propagation input received by the sample node has significantly increased at the current moment, and it is necessary to increase the impact of the current input on the hidden state update; in this case, the update gate output value is increased by 0.10 and then restricted to the interval between 0 and 1. The 0.18 threshold is obtained through historical sample statistics. Specifically, the historical risk propagation intensity distribution corresponding to the day before the rapid expansion of the actual lesion is statistically analyzed, and its 60th percentile is used as the threshold. Through the above processing, the graph-constrained time series prediction model can not only extract the general risk recursive law, but also highlight the abnormal propagation enhancement stage.

[0073] After recursively encoding each sample tree node for seven consecutive time steps, the hidden state sequence at each time step is retained and defined as the node temporal evolution feature. The node temporal evolution feature includes both the historical risk change information of the sample tree node itself and the propagation influence from neighboring sample tree nodes under the diffusion edge weight constraint, providing a basis for the next step to calculate the future lesion diffusion probability value and directional component.

[0074] In one embodiment, the directional response decoding function is used to convert the temporal evolution features of nodes into lesion spread probability values ​​and directional components for consecutive future time periods. The future prediction duration is taken as three consecutive time periods, i.e., predicting the evolution results for the first, second, and third future days. The directional response decoding function consists of a probability decoding part and a directional decoding part. The two parts share the same set of node temporal evolution feature inputs, but their output targets are different.

[0075] The specific implementation of the probabilistic decoding part is as follows: The hidden state of each sample node at the last moment of the current time window is taken and concatenated with the weighted average of the hidden states of the same node at the last three moments to form the probabilistic prediction input vector. The weights of the last three moments are set to 0.2, 0.3, and 0.5 according to the principle of time proximity. Subsequently, a two-layer linear mapping is performed on the probabilistic prediction input vector. The first layer has an output dimension of 24, and the activation function is a linear rectified function; the second layer has an output dimension of 3, corresponding to the next three consecutive moments. The output of the second layer is then input into a log-odds compression function to obtain the lesion spread probability values ​​for the next three consecutive moments. All lesion spread probability values ​​fall within the range of 0 to 1; the larger the value, the higher the probability that the corresponding sample node will experience lesion expansion or increased risk at that future moment.

[0076] The specific method for directional decoding is as follows: First, taking each sample tree node as the center, read the diffusion weights and corresponding spatial azimuths of all directed edges pointing from that node to adjacent sample tree nodes. The spatial azimuths are defined using a Cartesian coordinate system, with due east as 0 degrees and increasing counterclockwise. Then, for each directed edge emanating from that sample tree node, calculate the directional response value. The directional response value is defined as the diffusion weight of the directed edge multiplied by the probability of lesion spread at a future time for the target sample tree node, and then multiplied by the directional consistency coefficient. The directional consistency coefficient describes the degree of consistency between the current edge direction and the historical expansion direction. The historical expansion direction is calculated using the risk center displacement direction of the sample tree node in the previous two time moments. If the risk center cannot be calculated at a certain time, the most recent valid time moment is used instead. The directional consistency coefficient is defined as a value between 0 and 1 obtained by translating and scaling the cosine of the angle between the current edge azimuth and the historical expansion direction; the smaller the angle, the larger the directional consistency coefficient.

[0077] After obtaining the directional response values ​​of all adjacent directions, the lateral and longitudinal directional components are calculated separately. The lateral directional component is equal to the sum of all directional response values ​​multiplied by their respective azimuth cosine values; the longitudinal directional component is equal to the sum of all directional response values ​​multiplied by their respective azimuth sine values. This calculation is performed for the next three future days to obtain the directional component sequence for each sample node at consecutive future times. To suppress directional oscillations caused by minimal directional responses, when the square root of the sum of the squares of the lateral and longitudinal directional components is less than 0.05, the directional component at that moment is uniformly set to a zero vector. The 0.05 threshold is determined by the average amplitude of directional components without significant expansion phases in the training samples. This processing avoids outputting false expansion directions at times when risk growth is not significant.

[0078] In one implementation, the probability values ​​of lesion spread at consecutive future times are first corrected for temporal consistency. Since the spread of pests and diseases typically involves a slow accumulation followed by acceleration, if the probability value of lesion spread on a particular tree node on the second or third future day is significantly higher than the previous day, while the probability on the previous day was too low, a sudden jump that does not conform to the biological transmission process is likely to occur. Therefore, the temporal consistency correction rule is defined as follows: if the difference between the probability value of lesion spread on the second future day and the probability value of lesion spread on the first future day is greater than 0.35, then the probability value of lesion spread on the second future day is adjusted to the probability value of lesion spread on the first future day plus 0.35; if the difference between the probability value of lesion spread on the third future day and the probability value of lesion spread on the second future day is greater than 0.35, then the same adjustment is made. The 0.35 threshold is derived from the 90th percentile of the daily increase in the probability of lesion spread in historical samples and is used to limit abnormal sudden increases. After the burst restriction is completed, a one-dimensional time smoothing filter with weights of 0.2, 0.3 and 0.5 is used to smooth the lesion spread probability values ​​on the next 1 day, the next 2 days and the next 3 days, to obtain the final lesion evolution probability sequence.

[0079] Subsequently, temporal consistency correction is performed on the directional components. Specifically, the directional angles between future days 1 and 2, and between future days 2 and 3, are calculated. If the directional angle between any two adjacent future times is greater than 75 degrees, and the directional amplitudes at both times are greater than 0.08, the directional change is considered too drastic and requires correction. The correction method is as follows: the direction at the time with the larger directional amplitude remains unchanged, and the direction at the other time is adjusted to the weighted average of the two, with the weights determined by the proportion of their respective amplitudes. If the directional amplitude at a certain time is less than 0.08, its direction is directly adjusted to the direction of the previous valid time. The 75-degree threshold is used to limit large-angle reversals in biological diffusion paths within a short period, and the 0.08 threshold is used to distinguish between valid directions and weak noise directions.

[0080] After time consistency correction, spatial vector synthesis is performed on the directional components. Specifically, at each future time step, the lateral and longitudinal directional components of all sample tree nodes are multiplied by their corresponding lesion spread probability values ​​and then weighted to highlight the directional contribution of high-risk sample tree nodes. Then, the directional components within a local neighborhood are spatially smoothed. The neighborhood is a set of sample trees with a radius of 5 meters centered on the sample tree node. The smoothing weight is the product of the Gaussian decay weight of the distance between nodes and the normalized value of the spread edge weight. After smoothing, the lateral and longitudinal directional components of each sample tree node at each future time step are recalculated and sequentially combined into a two-dimensional vector, thus obtaining the expansion directional vector. The direction of the expansion directional vector is determined by both the lateral and longitudinal components, and the length of the expansion directional vector represents the strength of the spread trend.

[0081] Finally, the lesion spread probability values ​​for each sample tree node on the next three days (day 1, day 2, and day 3) are arranged in chronological order to form a lesion evolution probability sequence; the two-dimensional direction vectors corresponding to the same time are arranged in chronological order to form an expansion direction vector sequence. If an overall trend result for the entire target forest area is required, the mean of the lesion evolution probability sequences of all sample tree nodes is taken, and all expansion direction vectors are weighted and summed according to the lesion spread probability values ​​to obtain a comprehensive lesion evolution probability sequence and a comprehensive expansion direction vector at the forest area scale.

[0082] In one embodiment, the training samples of the graph-constrained time-series prediction model consist of suspected pest and disease feature maps at consecutive time points and subsequent measured changes in lesions. During training, the risk evolution tensor of the sample tree nodes within the current time window is used as input, and the node risk growth labels and directional displacement labels actually observed at the next three consecutive time points are used as supervision targets. The lesion spread probability value is constrained by a combination of binary cross-entropy loss and mean squared error loss, where the binary cross-entropy loss is used to distinguish whether expansion has occurred, and the mean squared error loss is used to constrain the expansion probability magnitude; the directional component is constrained by a combination of cosine similarity loss and mean squared error loss, where the cosine similarity loss is used to constrain the consistency of the directional angle, and the mean squared error loss is used to constrain the directional magnitude. The weights of the four losses are 0.35, 0.25, 0.25, and 0.15, respectively. Parameter updates use an adaptive moment estimation optimization method, with an initial learning rate of 0.001, a batch size of 4 time windows, and 150 training rounds; when the verification loss does not decrease for 8 consecutive rounds, the learning rate is reduced to 0.5 of the original value. The linear mapping weights in the graph constraint propagation aggregation layer, the weight matrix in the gated recursive unit, and the linear mapping weights in the directional response decoding function can be determined using the above method.

[0083] The historical risk transmission intensity threshold of 0.18, the lesion spread probability surge threshold of 0.35, the directional amplitude thresholds of 0.05 and 0.08, and the directional angle threshold of 75 degrees can all be obtained through repeated statistical analysis of historical monitoring samples, with clear data sources and determination processes. Therefore, the lesion evolution probability sequence generated in this step can quantify the risk change trend in future continuous timeframes, and the generated expansion direction vector can characterize the dominant direction and strength of risk expansion. Both can be directly used as inputs for subsequent calculation of the graded early warning index and delineation of treatment areas.

[0084] Early warning output module: Based on the lesion evolution probability sequence and expansion direction vector, calculate the graded early warning index of the target forest area, and output early warning results including pest and disease type, risk level, spread path and treatment area.

[0085] In one implementation, for each sample node, the probability sequence of lesion evolution at consecutive future time points is read, with the probability values ​​for day 1, day 2, and day 3 being P1, P2, and P3, respectively. The risk growth rate is defined as the weighted sum of the probability differences between adjacent time points, specifically calculated as follows: The aforementioned weighting settings give greater weight to the later growth trend, highlighting the stage of accelerated lesion expansion.

[0086] The cumulative risk value is defined as the weighted average of the probability values ​​at consecutive future time points, and is calculated as follows: This definition emphasizes the cumulative impact of risk in subsequent time periods.

[0087] After obtaining the risk growth rate and cumulative risk value, the consistency index of risk propagation direction is calculated. Specifically, the expansion direction vector of the sample node at each future time point is read, and the angle between adjacent time points is calculated. This angle is then converted into a consistency value, calculated by linearly mapping the cosine of the angle to the interval between 0 and 1. Finally, the consistency values ​​at each time point are averaged to obtain the risk propagation direction consistency index. This index reflects whether the expansion path is stable and continuous.

[0088] In one implementation, a multidimensional risk assessment function is constructed using risk growth rate, cumulative risk value, and risk propagation direction consistency index as input variables. To ensure the comparability of variables with different dimensions, the three indicators are first normalized by linear normalization from minimum to maximum value. The multidimensional risk assessment function is defined as: Comprehensive Risk Score = 0.35 × Normalized Risk Growth Rate + 0.45 × Normalized Cumulative Risk Value + 0.20 × (1 - Normalized Direction Consistency Index). This setting prioritizes the cumulative risk value while considering the diffusion uncertainty caused by the instability of growth trends and directions.

[0089] To enhance the ability to identify extreme risk nodes, when the normalized cumulative risk value is greater than 0.85 and the normalized risk growth rate is greater than 0.60, an additional 0.10 is added to the comprehensive risk score, but the final result is limited to the range of 0 to 1. The thresholds of 0.85 and 0.60 are determined by the statistical distribution of historical lesion rapid expansion samples.

[0090] In one implementation, the comprehensive risk scores of all sample tree nodes are sorted, and a quantile-based method is used to determine multi-level early warning threshold intervals. Specifically, the comprehensive risk scores are sorted from smallest to largest, and the 50th, 75th, and 90th percentiles are used as the first, second, and third thresholds, respectively. Based on these thresholds, the sample tree nodes are divided into four risk levels: a comprehensive risk score less than the first threshold is defined as low risk; a score between the first and second thresholds is defined as medium risk; a score between the second and third thresholds is defined as high risk; and a score greater than the third threshold is defined as extremely high risk. After completing the node-level division, spatial distribution characteristics are further considered. Specifically, the eight-code neighborhood connectivity rule is used to cluster sample tree nodes of the same risk level, forming multiple spatially continuous regions. If the number of nodes in a spatially continuous region is less than three, its risk level is lowered by one level to eliminate isolated outliers.

[0091] In one implementation, a regional early warning index is calculated for each spatially continuous region. Specifically, the average of the comprehensive risk scores of all sample tree nodes within the region is first calculated as the basic risk value for the region; then, the weighted composite result of the expansion direction vectors within the region is calculated, with the weights being the comprehensive risk scores of each node; the amplitude of the composite direction vector is used as the diffusion intensity index and normalized.

[0092] The graded early warning index is defined as: Graded Early Warning Index = 0.6 × Regional Basic Risk Value + 0.4 × Diffusion Intensity Index. The graded early warning index is denoted as S, and the value of S ranges from 0 to 1. Subsequently, the final early warning level is determined based on the graded early warning index: when 0 ≤ S < 0.30, it is defined as a low early warning level; when 0.30 ≤ S < 0.55, it is defined as a medium early warning level; when 0.55 ≤ S < 0.75, it is defined as a high early warning level; and when 0.75 ≤ S ≤ 1, it is defined as an emergency early warning level. The above thresholds are determined through statistical results of different levels of pest and disease occurrence areas in historical monitoring data. When generating early warning results, the following content is output for each spatially continuous area: the pest and disease type is determined through the classification results of the aforementioned identification model; the risk level is determined according to the corresponding level of the graded early warning index; the diffusion path is obtained by extending the regional expansion direction vector along the spatial coordinates; the treatment area is obtained by expanding the spatially continuous area outwards by a certain distance along the expansion direction, the outward expansion distance being defined as the diffusion intensity index multiplied by the maximum influence radius, where the maximum influence radius is 10 meters.

[0093] Finally, the results from all spatially contiguous areas are summarized to form a complete early warning result for the target forest area. This result includes information on risk intensity, expansion trend, and spatial extent, and can be directly used for decision-making regarding pest and disease control in forest areas.

[0094] 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 deep learning-based forest disease and pest intelligent identification and early warning system, characterized in that: include: The acquisition and registration module acquires visible light canopy images, thermal infrared temperature maps, and corresponding environmental parameters of the target forest area, performs spatiotemporal registration, and generates a forest observation matrix. Feature extraction and fusion module: Input the forest observation matrix into a dual-branch deep network to extract canopy texture features, temperature difference anomalous features and humidity-heat coupling features, and fuse them to generate a suspected feature map of pests and diseases; Risk association modeling module: Based on the spatial clustering relationship of the suspected pest and disease feature map, construct a forest risk association map and extract the diffusion edge weights between adjacent sample trees; Temporal evolution prediction module: Input the diffusion edge weights and the suspected feature map of pests and diseases into the temporal prediction model to generate a lesion evolution probability sequence and an expansion direction vector; Early warning output module: Based on the lesion evolution probability sequence and expansion direction vector, calculate the graded early warning index of the target forest area, and output early warning results including pest and disease type, risk level, spread path and treatment area.

2. The deep learning-based forest disease and pest intelligent identification and early warning system according to claim 1, characterized in that: The process of generating the forest observation matrix includes: Radiometric correction and distortion correction were performed on visible light canopy images to extract standardized canopy texture images; The thermal infrared temperature map is subjected to noise suppression and temperature field reconstruction to generate a spatially continuous temperature distribution map; A time synchronization sequence is constructed based on the corresponding location environmental parameters, and the standardized canopy texture image and the temperature distribution map are timestamped. The geographic coordinate mapping model was used to spatially register the time-aligned multi-source data, and the data were fused at a unified grid scale to generate a forest observation matrix that includes canopy texture features, temperature field features and environmental parameter features.

3. The deep learning-based forest disease and pest intelligent identification and early warning system according to claim 1, characterized in that: The process of generating suspected feature maps of pests and diseases includes: The forest observation matrix is ​​divided into a canopy texture sub-matrix and a temperature environment sub-matrix according to the feature channels, and then input into a structurally symmetric dual-branch convolutional neural network for feature encoding to obtain the canopy texture feature tensor and the temperature difference response feature tensor. Based on the temperature difference response feature tensor and environmental parameter features, a humid-thermal coupling mapping function is constructed, and a humid-thermal coupling feature tensor is generated by channel-by-channel weighted calculation. The canopy texture feature tensor and the humid-thermal coupling feature tensor are spatially aligned and channel-separated, and the responses of each channel are adaptively enhanced through an attention weight allocation strategy to obtain a fused feature tensor. Threshold mapping and spatial smoothing are performed on the fused feature tensor to output a suspected feature map of pests and diseases.

4. The intelligent identification and early warning system for forest pests and diseases based on deep learning according to claim 1, characterized in that: The process of constructing the forest risk association map includes: Based on the suspected feature map of pests and diseases, suspected feature values ​​of each tree node are extracted, and an initial adjacency matrix is ​​constructed according to spatial proximity. The suspected feature values ​​are then assigned to each node as graph signals. A graph Fourier transform is performed on the graph signals to construct a graph Laplace operator and optimize the adjacency weights by minimizing the graph signal smoothness index, resulting in a first association weight matrix. A transmission cost matrix is ​​constructed based on the first association weight matrix and the node spatial coordinates. Using the normalized result of the suspected feature map of pests and diseases as the probability distribution, the optimal transmission model with entropy regularization is solved to obtain the probability transmission matrix between nodes. The probability transmission matrix is ​​then weighted and fused with the first association weight matrix to generate a forest risk association map.

5. The intelligent identification and early warning system for forest pests and diseases based on deep learning according to claim 4, characterized in that: The weights between adjacent sample nodes in the fusion result are used as the diffusion edge weights.

6. The intelligent identification and early warning system for forest pests and diseases based on deep learning according to claim 1, characterized in that: The process of inputting the diffusion edge weights and the suspected feature map of pests and diseases into the time series prediction model to generate the lesion evolution probability sequence and the expansion direction vector includes: The suspected feature maps of pests and diseases are stacked sequentially according to a preset time window, and the risk evolution tensor of the sample tree node is constructed by combining the diffusion edge weights. The risk evolution tensor of the sample tree node is input into the graph-constrained time series prediction model, and the historical risk propagation intensity between adjacent sample tree nodes is recursively encoded using diffusion edge weights to obtain the node time series evolution characteristics. Based on the temporal evolution characteristics of the nodes and the diffusion edge weights, a directional response decoding function is constructed to calculate the lesion diffusion probability value and corresponding directional component of each wood node in future consecutive moments; The lesion diffusion probability value and the directional component are subjected to time consistency correction and spatial vector synthesis to generate a lesion evolution probability sequence and an expansion direction vector.

7. The intelligent identification and early warning system for forest pests and diseases based on deep learning according to claim 1, characterized in that: The process of calculating the graded early warning index for the target forest area includes: Based on the lesion evolution probability sequence, the risk growth rate and cumulative risk value of each sample tree node in future continuous time are extracted, and the consistency index of risk propagation direction is calculated in combination with the expansion direction vector. A multidimensional risk assessment function is constructed based on the risk growth rate, cumulative risk value and risk propagation direction consistency index. Different risk factors are weighted and fused by setting weight coefficients to obtain the comprehensive risk score of the sample tree node. Based on the comprehensive risk score of the sample tree node and its spatial distribution characteristics, the quantile division method is used to determine the multi-level early warning threshold interval, and the risk level of the target forest area is divided. According to the spatial continuous area corresponding to each risk level and the dominant direction of the expansion direction vector, the graded early warning index of the target forest area is calculated.

8. The intelligent identification and early warning system for forest pests and diseases based on deep learning according to claim 7, characterized in that: The graded early warning index is denoted as S, and the value of S ranges from 0 to 1. The final early warning level is divided according to the graded early warning index: when 0≤S<0.30, it is defined as a low early warning level; when 0.30≤S<0.55, it is defined as a medium early warning level; when 0.55≤S<0.75, it is defined as a high early warning level; and when 0.75≤S≤1, it is defined as an emergency early warning level.

Images