Pollution source positioning method and system combining unmanned aerial vehicle remote sensing and AI image recognition

Abstract

The application provides a pollution source positioning method and system combining unmanned aerial vehicle remote sensing and AI image recognition, relates to the technical field of pollution source positioning, and first generates a multidimensional cross image stream containing a ground surface reflection image, an atmospheric transmission image and a vegetation response image and being associated and bound through multidimensional collaborative imaging of an unmanned aerial vehicle remote sensing device; then uses an AI image recognition model to cross verify and extract pollution characteristics from the multidimensional cross image stream to obtain a cross verification characteristic set; then constructs a pollution characteristic reverse tracing path model based on time-space evolution associated data of the cross verification characteristic set; then associates and maps a starting associated node of the pollution characteristic reverse tracing path model with an imaging space coordinate to lock an initial action area; and finally adjusts an area boundary through dynamic evolution trajectory associated data to output a pollution source positioning result, thereby significantly improving the accuracy and efficiency of pollution source positioning.

Description

Technical Field

[0001] This invention relates to the field of pollution source location technology, and more specifically, to a pollution source location method and system that combines UAV remote sensing and AI image recognition. Background Technology

[0002] In the field of environmental protection, accurate and timely location of pollution sources is crucial for effective pollution control and ecological protection. Traditional methods for locating pollution sources have many limitations.

[0003] On the one hand, ground-based monitoring methods have limited coverage, making it difficult to cover large areas, especially in remote or topographically complex regions where monitoring data has many blind spots and cannot comprehensively obtain pollution information. Moreover, ground-based monitoring often only obtains localized, single-dimensional data, making it difficult to grasp the spread of pollution from a macroscopic perspective.

[0004] On the other hand, while a single remote sensing technology can acquire information over a large area, the data obtained by different remote sensing methods can be one-sided. For example, ordinary optical remote sensing mainly reflects surface reflection information, and it is difficult to accurately capture the diffusion patterns of pollutants in the atmosphere and the response of vegetation after pollution. On the other hand, relying solely on atmospheric or vegetation remote sensing lacks the acquisition of information directly related to surface pollution sources, making it difficult to accurately locate pollution sources.

[0005] Furthermore, existing methods lack deep integration and intelligent processing in data processing and analysis, making it difficult to extract effective pollution characteristics from complex environmental data and accurately trace their sources. Summary of the Invention

[0006] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for locating pollution sources by combining unmanned aerial vehicle remote sensing and AI image recognition, the method comprising: Multi-dimensional collaborative imaging is performed using UAV remote sensing equipment to generate a multi-dimensional cross-image stream of the target monitoring area. The multi-dimensional cross-image stream includes surface reflection images, atmospheric transmission images, and vegetation response images. The surface reflection images, atmospheric transmission images, and vegetation response images are associated and bound according to the same spatial coordinates and time nodes. An AI image recognition model is used to extract pollution features from the multi-dimensional cross-image stream, generating a cross-verification feature set that includes optical features of pollutants, diffusion morphology correlation features, and environmental response correlation features. Based on the spatiotemporal evolution correlation data of the cross-verification feature set, a pollution feature reverse source tracing path model is constructed. The pollution feature reverse source tracing path model derives the starting correlation node of pollution diffusion in reverse based on the feature evolution trajectory. The initial associated node of the pollution feature reverse source tracing path model is associated and mapped with the imaging spatial coordinates of the UAV remote sensing equipment to lock the initial area of ​​pollution spread. The spatial boundary of the initial action area is adjusted by the dynamic evolution trajectory association data of the cross-verification feature set, and the location result of the pollution source in the target monitoring area is output.

[0007] Furthermore, embodiments of the present invention also provide a pollution source localization system combining UAV remote sensing and AI image recognition, characterized in that it includes: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to perform the above-described pollution source localization method combining UAV remote sensing and AI image recognition by executing the machine-executable instructions.

[0008] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions, the machine-executable instructions being stored in a computer-readable storage medium, a processor of a pollution source localization system combining UAV remote sensing and AI image recognition reading the machine-executable instructions from the computer-readable storage medium, the processor executing the machine-executable instructions, causing the pollution source localization system combining UAV remote sensing and AI image recognition to perform the above-described pollution source localization method combining UAV remote sensing and AI image recognition.

[0009] Based on the above, multi-dimensional collaborative imaging using UAV remote sensing equipment generates multi-dimensional cross-image streams, encompassing information from various aspects such as surface reflection, atmospheric transmission, and vegetation response. AI image recognition models are then used to cross-verify and extract pollution features from these multi-dimensional cross-image streams. This allows for the extraction of optical characteristics, diffusion morphology correlation features, and environmental response correlation features of pollutants from different dimensions of data, forming a cross-verification feature set, which significantly improves the accuracy and reliability of pollution feature extraction. Based on the spatiotemporal evolution correlation data of the cross-verification feature set, a pollution feature reverse source tracing path model is constructed. Using the feature evolution trajectory as a basis, the starting correlation node of pollution diffusion is deduced in reverse. The initial correlation node is mapped and associated with the imaging spatial coordinates to lock the initial area of ​​action. The spatial boundary is adjusted through dynamic evolution trajectory correlation data, ultimately outputting accurate pollution source location results, significantly improving the accuracy and efficiency of pollution source location. Attached Figure Description

[0010] Figure 1 This is a schematic diagram of the execution flow of the pollution source localization method combining UAV remote sensing and AI image recognition provided in an embodiment of the present invention.

[0011] Figure 2This is a schematic diagram of exemplary hardware and software components of a pollution source localization system combining UAV remote sensing and AI image recognition, provided in an embodiment of the present invention. Detailed Implementation

[0012] Figure 1 This is a flowchart illustrating a pollution source localization method combining UAV remote sensing and AI image recognition, provided in one embodiment of the present invention. A detailed description follows.

[0013] Step S110: Perform multi-dimensional collaborative imaging using UAV remote sensing equipment to generate a multi-dimensional cross-image stream of the target monitoring area. The multi-dimensional cross-image stream includes surface reflection images, atmospheric transmission images, and vegetation response images. The surface reflection images, atmospheric transmission images, and vegetation response images are associated and bound according to the same spatial coordinates and time nodes.

[0014] In this embodiment, the target monitoring area is set as an ecologically sensitive area surrounding an industrial park, where pollutant emissions from industrial production processes may occur. The UAV remote sensing equipment is equipped with a multispectral imaging sensor, a thermal infrared imaging sensor, and a hyperspectral imaging sensor, which are used to collect image data from different dimensions.

[0015] Before executing the imaging mission, the UAV plans its flight path via a ground control station, employing a regional coverage flight path design to ensure a comprehensive scan of the target monitoring area. During flight, the UAV acquires its own spatial coordinates in real time using a Global Positioning System (GPS) and binds this coordinate information to the timestamp of the imaging moment. The timestamp accuracy is down to the millisecond level to ensure the time synchronization of the image data. Surface reflection images are acquired by a multispectral imaging sensor in the visible and near-infrared bands. The number of bands on the sensor is set according to the monitoring requirements, with each band corresponding to a specific wavelength range to capture the surface's reflection characteristics of different wavelengths of light. Atmospheric transmission images are acquired by a thermal infrared imaging sensor, operating in a specific infrared band range, reflecting the temperature distribution in the atmosphere and the transmission of pollutants to thermal radiation. Vegetation response images are acquired by a hyperspectral imaging sensor, which contains multiple continuous spectral bands to precisely capture the spectral response characteristics of vegetation.

[0016] When generating a multi-dimensional cross-image stream, three types of images acquired at the same time and spatial location are associated and bound together. Specifically, a spatial coordinate matching algorithm is used to map the pixels of the three images to a unified geographic coordinate system, and image registration technology is used to eliminate spatial offsets between different sensors. Simultaneously, based on timestamp information, image frames acquired at similar times are associated to form a spatiotemporally synchronized multi-dimensional cross-image stream. For example, if the coordinates of a pixel in the surface reflection image are (Xa, Ya) and the acquisition time is Ta, pixels with the same coordinates (Xa, Ya) and acquisition times within a preset time range near Ta are searched in the atmospheric transmission image and vegetation response image and associated to form an image stream unit containing multi-dimensional information.

[0017] Step S120: Use an AI image recognition model to extract pollution features from the multi-dimensional cross-image stream, and generate a cross-verification feature set that includes optical features of pollutants, diffusion morphology correlation features, and environmental response correlation features.

[0018] In this embodiment, the AI ​​image recognition model adopts a multimodal fusion deep learning architecture, specifically designed to process multi-dimensional cross-image stream data. The input of the model is the multi-dimensional cross-image stream generated in step S110, and the output is a cross-verification feature set containing multiple pollution features. In the model processing flow, the input multi-dimensional cross-image stream is first preprocessed. Preprocessing includes image denoising, size normalization, and pixel value standardization. Image denoising uses a wavelet transform-based denoising algorithm, which removes the high-frequency components corresponding to noise by decomposing the different frequency components of the image; size normalization adjusts the images of the three different modalities to the same spatial resolution, and uses bilinear interpolation to scale the images; pixel value standardization converts the pixel values ​​of each image to a specific numerical range, for example, by subtracting the mean and dividing by the standard deviation, to eliminate the dimensional differences between images of different modalities. After preprocessing, the model enters the feature extraction stage, and independent feature extraction subnetworks are designed for surface reflection images, atmospheric transmission images, and vegetation response images. The feature extraction subnetwork for surface reflectance images employs a convolutional neural network structure, extracting deep spectral and spatial features through multiple convolutional and pooling layers. The feature extraction subnetwork for atmospheric transmission images uses an encoder-decoder structure, focusing on capturing morphological and structural information within the image. The feature extraction subnetwork for vegetation response images combines spectral angle analysis to extract the spectral response features of vegetation. Subsequently, a cross-validation module fuses the extracted multimodal features. This module uses an attention mechanism to assign weights to features of different modalities based on their importance, and then concatenates the weighted features to form preliminary fused features. Finally, a feature selection algorithm removes redundant features, retaining key pollution-related features to generate a cross-validation feature set.

[0019] Step S121: Input the multi-dimensional cross-image stream into the multi-dimensional feature separation module of the AI ​​image recognition model, and perform feature channel separation according to the modal attributes of the surface reflection image, atmospheric transmission image and vegetation response image to obtain the modal feature channel corresponding to the surface reflection image, the modal feature channel corresponding to the atmospheric transmission image and the modal feature channel corresponding to the vegetation response image.

[0020] In this embodiment, the multi-dimensional feature separation module consists of three parallel feature extraction branches, each corresponding to a specific modality of image data. When the surface reflection image is input into this module, it enters the first feature extraction branch. This branch contains multiple convolutional layers and pooling layers. The convolutional layers use convolutional kernels of a specific size and extract features by sliding them across the image. Each convolutional layer is followed by a non-linear activation function to enhance the network's expressive power. The pooling layers use max pooling to downsample the feature maps output by the convolutional layers, reducing feature dimensionality while retaining important feature information. After multiple convolutional and pooling operations, the modal feature channel corresponding to the surface reflection image is obtained. This feature channel is a three-dimensional feature tensor containing deep feature information of the surface reflection image. Similarly, the atmospheric transmission image is input into the second feature extraction branch. The network structure of this branch is similar to that of the surface reflection image feature extraction branch, but parameters such as the number of convolutional kernels and the depth of convolutional layers are adjusted according to the characteristics of the atmospheric transmission image to better extract the modal features of the atmospheric transmission image, ultimately obtaining the modal feature channel corresponding to the atmospheric transmission image. The vegetation response image is input into the third feature extraction branch. This branch, targeting the hyperspectral characteristics of the vegetation response image, employs more convolutional layers and smaller convolutional kernels to extract refined spectral features. After processing, the modal feature channels corresponding to the vegetation response image are obtained. These three modal feature channels correspond to three different modalities of image data.

[0021] Step S122: Perform optical attribute analysis on the modal feature channels corresponding to the surface reflection image, and extract the optical parameters in the modal feature channels corresponding to the surface reflection image that differ from the conventional surface reflection characteristics. The optical parameters in the modal feature channels corresponding to the surface reflection image that differ from the conventional surface reflection characteristics constitute the basic elements of the optical characteristics of pollutants.

[0022] Step S1221: Obtain the surface reflectance reference data of the target monitoring area under uncontaminated conditions. The surface reflectance reference data includes the conventional reflectance spectral parameters and reflectance intensity parameters of different surface types in the target monitoring area.

[0023] In this embodiment, the surface reflectance reference data of the target monitoring area under uncontaminated conditions is obtained through historical data collection and field surveys. First, remote sensing image data of the target monitoring area over the past few years are collected, and images from periods without pollution events are selected as reference images. These reference images are processed to extract reflectance spectral information and reflectance intensity information for different surface types (such as soil, vegetation, and water bodies). Simultaneously, multiple representative sampling points are selected within the target monitoring area for field measurements. Portable spectrometers are used to collect reflectance spectral data for different surface types in multiple bands. Multiple measurements are performed at each sampling point, and the average value is taken as the reflectance spectral parameter for that sampling point. For the reflectance intensity parameter, statistical analysis is performed on the reflectance intensity values ​​of different surface type areas in the reference images to calculate their average value and standard deviation, thus determining the range of conventional reflectance intensity. The collected reflectance spectral parameters and reflectance intensity parameters are categorized and organized according to surface type to establish a surface reflectance reference database. This database contains the conventional optical characteristic parameters of various surface types within the target monitoring area under uncontaminated conditions.

[0024] Step S1222: The modal feature channels corresponding to the surface reflection image are split into pixels to obtain the reflection spectrum data and reflection intensity data of each pixel.

[0025] In this embodiment, the modal feature channel corresponding to the surface reflection image is a three-dimensional feature tensor with dimensions (H, W, C), where H and W represent the height and width of the feature map, corresponding to the number of pixel rows and columns of the image, respectively, and C represents the number of feature channels, corresponding to different spectral bands. When splitting by pixel, each spatial position (i, j) of the feature tensor is traversed, where the value of i ranges from 0 to H-1, and the value of j ranges from 0 to W-1. For each pixel (i, j), its values ​​on all C feature channels are extracted. These values ​​form a vector of length C, which is the reflection spectrum data of that pixel. This reflection spectrum data reflects the reflection characteristics of the pixel in different spectral bands. Simultaneously, according to application requirements, a specific feature channel is selected as the source of reflection intensity data. For example, a certain band in the visible light band is selected, and the value of the pixel on that channel is extracted as the reflection intensity data. Through the above method, the entire modal feature channel is split into reflection spectrum data and reflection intensity data of H×W pixels, with each pixel corresponding to a set of independent optical parameters.

[0026] Step S1223: Compare the reflectance spectral data of each pixel with the corresponding conventional reflectance spectral parameters of the land surface type in the land surface reflectance reference data to identify pixels with differences in spectral curve shape, peak position, and valley position.

[0027] In this embodiment, for the reflectance spectral data of each pixel, the land surface type is first determined based on its spatial coordinates. By matching the pixel's coordinates with the land surface type distribution map of the target monitoring area, the corresponding land surface type, such as farmland, forest, or water body, is determined. Then, the conventional reflectance spectral parameters for this land surface type are obtained from the land surface reflectance reference data. A spectral angle matching algorithm is used to compare the pixel's reflectance spectral data with the conventional reflectance spectral parameters, measuring their similarity by calculating the angle between the two spectral vectors. Specifically, the two spectral vectors are normalized to a magnitude of 1; their dot product is calculated, and the dot product is equal to the cosine of the angle between the two vectors; the angle is obtained using the inverse cosine function. When the calculated angle is greater than a preset spectral angle threshold, it indicates that the spectral curve shape, peak position, and valley position of the pixel are significantly different from the conventional reflectance spectrum, and the pixel is marked as a spectral difference pixel. Through this method, all pixels are compared one by one to filter out pixels with spectral differences.

[0028] Step S1224: Compare the reflectance intensity data of pixels with differences in spectral curve shape, peak position, and valley position with the conventional reflectance intensity parameters of the corresponding land surface type, extract pixels whose reflectance intensity deviates from the conventional range, and form a set of suspected polluted pixels.

[0029] In this embodiment, for the spectral difference pixels identified in step S1223, their reflectance intensity data is extracted. The range of normal reflectance intensity parameters for the corresponding surface type is obtained from the surface reflectance reference data. This range is typically expressed as the mean of the reflectance intensity for that surface type plus or minus a certain standard deviation. The reflectance intensity data of the pixel is compared with this range. If the reflectance intensity data is less than the lower limit of the range or greater than the upper limit, the reflectance intensity of the pixel is considered to deviate from the normal range. All pixels whose reflectance intensity deviates from the normal range are collected to form a suspected pollution pixel set. The pixels in this suspected pollution pixel set differ from normal surface reflectance characteristics in both spectral morphology and reflectance intensity, and may be affected by pollutants.

[0030] Step S1225: Perform spatial clustering on the pixels in the suspected contaminated pixel set, and group pixels that are spatially adjacent and have consistent spectral and intensity differences into the same pixel cluster.

[0031] In this embodiment, a density clustering algorithm is used to spatially cluster the set of suspected contaminated pixels. First, two parameters are set: neighborhood radius and minimum sample size. The neighborhood radius defines the spatial range surrounding a pixel, with other pixels within this range considered neighbors. The minimum sample size specifies the minimum number of pixels required to form a cluster. For each pixel in the suspected contaminated pixel set, the number of neighbors within its neighborhood radius is checked. If the number of neighbors is greater than or equal to the minimum sample size, the pixel is marked as a core point, and a cluster is grown from this point. The core point and all pixels within its neighborhood are added to the cluster. This process is repeated for newly added pixels, checking if any new pixels can be added to their neighborhoods, until no new pixels can be added to the cluster. Pixels with fewer neighbors than the minimum sample size but located within the neighborhood of the core point are also added to the cluster. In this way, pixels that are spatially adjacent and have consistent spectral and intensity differences are grouped into the same pixel cluster, with each cluster representing a potential contaminated area.

[0032] Step S1226: Extract the average reflectance spectrum curve of each pixel cluster. The average reflectance spectrum curve is the result of averaging the reflectance spectrum data of all pixels in the cluster, representing the overall spectral data characteristics of the pixel cluster.

[0033] In this embodiment, for each pixel cluster, the reflectance spectral data of all pixels within it are collected. The reflectance spectral data of each pixel is a vector of length C, corresponding to C spectral bands. The arithmetic mean of the reflectance values ​​of all pixels in each spectral band is calculated to obtain the average reflectance value for that band. The average reflectance values ​​of all bands are arranged in band order to form a new vector of length C, which is the average reflectance spectral curve of the pixel cluster. The average reflectance spectral curve can represent the overall spectral characteristics of the pixel cluster, eliminating noise interference from individual pixels and better reflecting the spectral characteristics of the region represented by the cluster. For example, a pixel cluster contains multiple pixels, each with reflectance values ​​in various bands. By averaging the reflectance values ​​in each band, the resulting average reflectance spectral curve can show the reflectance characteristics of the cluster in different bands.

[0034] Step S1227: Calculate the average reflection intensity value of each pixel cluster. This average reflection intensity value is the statistical average of the reflection intensity data of all pixels in the cluster, reflecting the overall reflection intensity level of the pixel cluster.

[0035] In this embodiment, for each pixel cluster, the reflectance intensity data of all pixels within it is collected. Reflectance intensity data represents the reflectance value of each pixel in a specific wavelength band. The reflectance intensity data is then statistically averaged, i.e., the sum is divided by the number of pixels to obtain the average reflectance intensity value of the pixel cluster. The average reflectance intensity value reflects the overall reflectance intensity level of the pixel cluster. Comparing it with conventional reflectance intensity parameters allows for a direct visual indication of any anomalies in the cluster's reflectance intensity. For example, if a pixel cluster contains multiple pixels, the reflectance intensity values ​​of these pixels are summed and then divided by the number of pixels to obtain the average reflectance intensity value of the cluster. If this average reflectance intensity value deviates significantly from the average reflectance intensity of the corresponding surface type, it may indicate the presence of pollution.

[0036] Step S1228: Analyze the degree of difference between the average reflectance spectrum curve of each pixel cluster and the conventional reflectance spectrum parameters, and determine the distribution range and magnitude of the difference bands through spectral similarity analysis.

[0037] In this embodiment, for the average reflectance spectral curve of each pixel cluster, a spectral similarity analysis is performed with the conventional reflectance spectral parameters of the corresponding surface type. The similarity between the two is calculated using the spectral correlation coefficient method. The calculation process of the spectral correlation coefficient is as follows: first, the covariance of the average reflectance spectral curve and the conventional reflectance spectral parameters is calculated, and then divided by the product of their standard deviations to obtain the correlation coefficient. The correlation coefficient ranges from -1 to 1. The closer it is to 1, the higher the similarity between the two spectral curves; the closer it is to -1, the lower the similarity. When the correlation coefficient is less than a preset similarity threshold, it is considered that there is a significant difference between the average reflectance spectral curve and the conventional reflectance spectral parameters. For spectral curves with differences, the distribution range of the difference bands is determined by comparing the reflectance values ​​of the average reflectance spectral curve and the conventional reflectance spectral parameters band by band, i.e., the band intervals where the reflectance value difference exceeds the preset difference threshold. For each difference band, the average value of the absolute values ​​of the reflectance difference between the average reflectance spectral curve and the conventional reflectance spectral parameters in that band is calculated as the difference amplitude, which reflects the degree of difference between the two spectral curves in the difference band.

[0038] Step S1229: Extract the change in reflection intensity data corresponding to the difference band. This change in reflection intensity data is the difference between the average reflection intensity value and the conventional reflection intensity parameter within the difference band.

[0039] In this embodiment, for the difference band determined in step S1228, the average reflectance intensity value of the pixel cluster within that band and the conventional reflectance intensity parameter of the corresponding surface type are extracted. The difference between the two is calculated, i.e., the average reflectance intensity value minus the conventional reflectance intensity parameter, to obtain the change in reflectance intensity data corresponding to the difference band. This change in reflectance intensity data reflects the degree of deviation of the reflectance intensity of the pixel cluster from the conventional situation within the difference band. If the change is positive, it indicates that the reflectance intensity of that band is higher than the conventional level; if it is negative, it indicates that it is lower than the conventional level. The magnitude of the change in reflectance intensity data can be used as an indicator to judge the degree of pollution; the larger the change, the more severe the pollution may be.

[0040] Step S12210: Integrate the average reflectance spectrum curve, average reflectance intensity value, differential band distribution range and differential amplitude, and differential band reflectance intensity data change of each pixel cluster to form the basic elements constituting the optical characteristics of pollutants.

[0041] In this embodiment, for each pixel cluster, the average reflectance spectrum curve obtained in step S1226, the average reflectance intensity value obtained in step S1227, the difference band distribution range and difference amplitude obtained in step S1228, and the change in reflection intensity data of the difference band obtained in step S1229 are integrated. The above data is organized according to a certain structure, for example, forming a feature vector, which includes the values ​​of each band of the average reflectance spectrum curve, the average reflectance intensity value, the start and end wavelengths of the difference bands, the difference amplitude, and the change in reflection intensity data. The integrated data corresponding to each pixel cluster constitutes the basic elements of the optical characteristics of pollutants, and these elements can comprehensively reflect the optical properties of pollutants in the region represented by the cluster.

[0042] Step S123: Perform morphological structure analysis on the modal feature channels corresponding to the atmospheric transmission image, identify image regions in the image that are continuously distributed and have a diffusion trend, extract the contour morphology, distribution density and extension related data of the image region, and the contour morphology, distribution density and extension related data of the image region form the constituent content of diffusion morphology related features.

[0043] Step S1231: Perform image enhancement processing on the modal feature channels corresponding to the atmospheric transmission image, and use a region growing algorithm to segment the enhanced image into regions. Use pixels with gray value differences within the preset gray value tolerance range and spatially adjacent as seed points to gradually grow into continuous image regions.

[0044] In this embodiment, image enhancement processing is first performed on the modal feature channels corresponding to the atmospheric transmission image. An adaptive histogram equalization algorithm is used to divide the image into multiple sub-regions, and histogram equalization is performed on the histogram of each sub-region to improve the local contrast of the image and make the details in the image clearer, especially the features of low-contrast areas. The enhanced image is then segmented into regions using a region growing algorithm. First, seed points are selected based on the gray value distribution of the image. Pixels with gray values ​​within the gray value range of the suspected pollution area are found in the image, and these pixels are spatially independent. For each seed point, the pixels in its 8-neighborhood are checked, and the gray value difference between the neighboring pixels and the seed point is calculated. If the gray value difference is within a preset gray value tolerance range, the neighboring pixels are merged into the region where the seed point is located, and the neighboring pixels are used as new seed points for further growth. This process is repeated until no new pixels can be merged into the region, thus forming a continuous image region.

[0045] Step S1232: Perform morphological analysis on each segmented continuous image region, calculate the morphological parameters of area, perimeter, and circularity of each image region, and preliminarily screen out image regions with irregular shapes.

[0046] In this embodiment, morphological analysis is performed on each segmented continuous image region. The area of ​​each region is calculated by counting the number of pixels within the region and multiplying this number by the actual area represented by a single pixel. The perimeter is calculated by extracting the boundary pixels of the region, connecting these boundary pixels sequentially to form a closed curve, and then calculating the length of this curve. Circularity is calculated using the formula 4π × area / perimeter². The closer the circularity is to 1, the closer the region is to a circle; conversely, the lower the circularity, the more irregular the shape. Based on these morphological parameters, a circularity threshold is set. When the circularity of an image region is less than this threshold, it is considered irregular in shape and is initially filtered out, as irregularly shaped regions are more likely to be associated with pollution diffusion.

[0047] Step S1233: Extract the coordinates of the boundary contour points of each filtered image region, and construct the contour morphology model of the image region through a curve fitting algorithm. The contour morphology model describes the external shape and edge undulation features of the image region.

[0048] In this embodiment, for the selected irregularly shaped image regions, an edge detection algorithm is used to extract the coordinates of their boundary contour points. The edge detection algorithm first performs Gaussian filtering on the image to remove noise, then calculates the gradient magnitude and direction, next performs non-maximum suppression on the gradient magnitude to retain local maxima, and finally obtains the boundary contour points of the region through double threshold detection and edge connection. The coordinates of these boundary contour points are recorded sequentially to form a coordinate sequence. A B-spline curve fitting algorithm is used to fit this coordinate sequence. B-spline curves have good smoothness and local controllability, and can accurately describe the shape of the boundary contour. By adjusting the control vertices and order of the B-spline curve, the fitted curve is made as close as possible to the boundary contour points, constructing a contour morphology model of the image region. This contour morphology model can describe the external shape and edge undulation features of the image region.

[0049] Step S1234: Calculate the distribution density of pixels in each image region. The distribution density is obtained by the ratio of the number of pixels in the image region to the area of ​​the image region, reflecting the density of material distribution in the image region.

[0050] In this embodiment, for each image region, the number of pixels N it contains is counted. The area S of the image region is calculated in steps S1232, which is the number of pixels in the region multiplied by the actual area represented by a single pixel. The distribution density D is calculated using the formula D=N / S, where N is the number of pixels in the region and S is the area of ​​the region. The distribution density reflects the density of pixels in the image region, that is, the density of material distribution; the higher the distribution density value, the denser the material distribution in the region.

[0051] Step S1235: Analyze the distribution density change gradient of each image region, extract multiple density sampling lines from the center of the image region to the boundary, calculate the density value change at each point on the sampling line, and determine the direction and magnitude of the density gradient. The density gradient direction points in the direction of decreasing density.

[0052] In this embodiment, the center coordinates of each image region are first determined by calculating the average of the coordinates of all pixels within the region. Starting from the center coordinates, multiple sampling lines (e.g., eight sampling lines) are drawn evenly towards the region boundary, each pointing in a different direction. On each sampling line, multiple sampling points are selected at equal intervals from the center to the boundary. For each sampling point, the distribution density value at its location is calculated. This distribution density value is obtained by dividing the number of pixels within a local window centered on the sampling point by the window area. The window size is set according to the image resolution and region size. The density difference between adjacent sampling points is calculated. The direction of the density gradient is determined by the sign of the density difference; the density gradient points in the direction of decreasing density, i.e., from high-density regions to low-density regions. The magnitude of the density gradient is the ratio of the absolute value of the density difference to the distance between the sampling points.

[0053] Step S1236: Determine the extension trend of the image region based on the density gradient direction. If the density gradients in multiple directions show an outward diffusion characteristic, then the image region has a diffusion trend.

[0054] In this embodiment, for each image region, the density gradient directions of all sampling lines are collected. If the density gradient direction of most sampling lines points from the center of the region to the boundary, i.e., exhibits an outward diffusion characteristic, then the image region is determined to have a diffusion trend. For example, among the eight sampling lines, if the density gradient direction of six or more sampling lines points from the center to the boundary, it indicates that the density of the region gradually decreases from the center to the boundary, showing an outward diffusion trend.

[0055] Step S1237: Extract the extension direction parameters of the image region with diffusion trend. The extension direction parameters are obtained by the statistical average of each density gradient direction, which represents the overall diffusion direction of the image region.

[0056] In this embodiment, for image regions exhibiting a diffusion trend, the density gradient directions of all sampling lines are converted into directional angles, ranging from 0 to 360 degrees. The arithmetic mean of these directional angles is calculated to obtain the average directional angle, which is the extension direction parameter. The extension direction parameter represents the overall diffusion direction of the image region. For example, the average value of the density gradient directional angles of multiple sampling lines is a certain angle, which represents the main diffusion direction of the image region.

[0057] Step S1238: Calculate the relevant parameters of the extension amplitude of the image region. By comparing the area change of the same image region at different time points, obtain the area change rate of the image region per unit time. By comparing the boundary movement distance of the same image region at different time points, obtain the boundary movement rate of the image region per unit time. Combine the area change rate and the boundary movement rate to describe the extension amplitude of the image region.

[0058] In this embodiment, area data of the same image region at different time points are acquired. The area change is calculated, which is the area at the later time point minus the area at the earlier time point, and then divided by the time interval to obtain the area change rate per unit time. Simultaneously, the coordinates of the boundary contour points of the image regions at the two time points are extracted. For each feature point on the boundary, its movement distance between the two time points is calculated. The average of the movement distances of all feature points is taken as the boundary movement distance, and then divided by the time interval to obtain the boundary movement rate per unit time. The area change rate and the boundary movement rate are combined as a parameter describing the extent of image region extension.

[0059] Step S1239: Integrate the contour morphology model, distribution density data, density gradient information, extension direction parameters, and extension amplitude data of each image region with a diffusion trend to form diffusion morphology association features.

[0060] In this embodiment, for each image region with a diffusion trend, the contour morphology model obtained in step S1233, the distribution density data obtained in step S1234, the density gradient information (including gradient direction and magnitude) obtained in step S1235, the extension direction parameters obtained in step S1237, and the extension amplitude data (area change rate and boundary movement rate) obtained in step S1238 are integrated. The above data is organized into a structured dataset, such as a data object containing multiple fields, each field corresponding to a parameter, forming a diffusion morphology association feature.

[0061] Step S124: Detect response characteristics of the modal feature channels corresponding to the vegetation response image, capture the features of color change and texture variation of vegetation pixels in the image that differ from the normal vegetation response, and construct the constituent elements of environmental response association features based on the features of color change and texture variation of vegetation pixels in the image that differ from the normal vegetation response.

[0062] In this embodiment, response characteristic detection is performed on the modal feature channels corresponding to the vegetation response image. First, vegetation pixels are extracted from the image using a vegetation index, such as the Normalized Difference Vegetation Index (NDVI), which is calculated based on the reflectance values ​​of the near-infrared and red light bands. Pixels with an NDVI value greater than a preset threshold are identified as vegetation pixels. For the extracted vegetation pixels, color change detection is performed. The vegetation response image is converted to the HSV color space, and the hue, saturation, and brightness components of each vegetation pixel are extracted and compared with the vegetation color reference data under unpolluted conditions. The difference value of each component is calculated. When the difference value exceeds a preset threshold, the vegetation pixel is considered to have a color change. Simultaneously, texture variation detection is performed. The gray-level co-occurrence matrix method is used to calculate the texture features of the vegetation pixels, such as contrast, energy, entropy, and correlation. The calculated texture features are compared with the texture reference features of normal vegetation. When the difference exceeds a preset threshold, texture variation is determined to exist. Vegetation pixels with color variations or texture variations are extracted, and these features together constitute the components of environmental response-related features.

[0063] Step S125: Activate the cross-verification module of the AI ​​image recognition model to identify the spatial location range corresponding to the basic elements of the optical characteristics of pollutants and the spatial location range corresponding to the constituent content of the diffusion morphology-related characteristics, and confirm the overlapping area of ​​the two within the same spatial coordinate range; within the overlapping area, extract the basic element parameters of the optical characteristics of pollutants and the constituent content parameters of the diffusion morphology-related characteristics respectively, analyze the synergistic change relationship between the two in terms of intensity, distribution or morphology, and obtain the correlation intensity data.

[0064] In this embodiment, the cross-verification module of the AI ​​image recognition model first obtains the basic elements of the optical characteristics of pollutants and the constituent content of the diffusion morphology-related features. For the basic elements of the optical characteristics of pollutants, their spatial location range is determined based on the spatial coordinates of their corresponding pixel clusters. This spatial location range is typically a polygonal region composed of the boundary pixel coordinates of the pixel clusters. Similarly, for the constituent content of the diffusion morphology-related features, their spatial location range is determined based on their contour morphology model. The two spatial location ranges are spatially superimposed and their intersection area is calculated. This intersection area is the overlapping area of ​​the two within the same spatial coordinate range. If the area of ​​the overlapping area exceeds a preset threshold in proportion to the total area of ​​the two spatial location ranges, then a spatial correlation is considered to exist between the two. Within the overlapping area, the basic element parameters of the optical characteristics of pollutants are extracted, such as the average reflection intensity value and the change in reflection intensity of the differential band; simultaneously, the constituent content parameters of the diffusion morphology-related features are extracted, such as the distribution density and the extension amplitude. The synergistic relationship between these parameters is analyzed. For example, when the average reflection intensity value of the optical characteristics of pollutants increases, does the distribution density of the diffusion morphology-related features also increase accordingly, or does the extension amplitude show the same trend? The correlation coefficient between the two parameters is calculated using correlation analysis, and the magnitude of the correlation coefficient is the correlation strength data.

[0065] Step S126: Perform a time-series correlation analysis on the basic elements of the optical characteristics of pollutants and the constituent elements of the environmental response correlation characteristics corresponding to the above-mentioned overlapping areas to verify the consistency of the evolution of the basic elements of the optical characteristics of pollutants, the constituent elements of the diffusion morphology correlation characteristics, and the constituent elements of the environmental response correlation characteristics at the same time point.

[0066] In this embodiment, data on the basic elements of the optical characteristics of pollutants, the components of diffusion morphology-related characteristics, and the components of environmental response-related characteristics corresponding to the overlapping areas at multiple time points are collected. For each time point, key parameters of these three characteristics are extracted, such as the average reflectance intensity value of the optical characteristics of pollutants, the area change rate of diffusion morphology-related characteristics, and the color difference value of environmental response-related characteristics. These parameters are arranged in chronological order to form three parameter sequences. A time-series correlation analysis is performed on these three parameter sequences to calculate their synchronicity and correlation. Synchronicity analysis is achieved by comparing whether the changing trends of the three parameter sequences are consistent. For example, when the average reflectance intensity value of the optical characteristics of pollutants begins to increase at a certain time point, do the area change rate of diffusion morphology-related characteristics and the color difference value of environmental response-related characteristics also begin to increase at the same time point? Correlation analysis calculates the pairwise correlation coefficients between the three parameter sequences to assess their linear correlation. Evolutionary consistency data includes synchronicity indices and correlation coefficients. When both the synchronicity indices and correlation coefficients exceed preset thresholds, it indicates that the three characteristics have good evolutionary consistency at the same time point.

[0067] Step S127: Based on the spatial location comparison results and temporal correlation analysis results, assign weights to the basic elements of the optical characteristics of pollutants, the components of the diffusion morphology correlation characteristics, and the components of the environmental response correlation characteristics. The basic elements of the optical characteristics of pollutants are assigned the highest weight, and the components of the diffusion morphology correlation characteristics and the components of the environmental response correlation characteristics are assigned corresponding weights.

[0068] In this embodiment, the spatial location comparison result is represented by the area ratio of overlapping regions, and the temporal correlation analysis result is measured by the correlation coefficient in the evolutionary consistency data. The weight allocation rules are set as follows: the basic elements of the optical characteristics of pollutants are considered core features and assigned the highest weight. The weights of the components of the diffusion morphology correlation features are determined based on the area ratio of overlapping regions in the spatial location comparison results and the correlation coefficient in the temporal correlation analysis. The calculation formula is: Weight_D = a × (area ratio of overlapping regions) + b × (correlation coefficient), where a and b are pre-set coefficients based on feature importance. The weight calculation method for the components of the environmental response correlation features is similar: Weight_E = c × (area ratio of overlapping regions) + d × (correlation coefficient), where c and d are pre-set coefficients. Then, the weights are normalized so that the sum of the three weights is 1, ensuring the rationality of the weight allocation.

[0069] Step S128: The basic elements of the optical characteristics of pollutants, the constituent elements of the diffusion morphology-related characteristics, and the constituent elements of the environmental response-related characteristics are fused to form a preliminary fused feature set. Redundant features are removed from the preliminary fused feature set, and the feature elements that are directly related to the diffusion of pollution and have unique identification characteristics are retained.

[0070] In this embodiment, the fusion process employs a weighted concatenation method. The basic elements of the pollutant's optical characteristics, the components of the diffusion morphology-related characteristics, and the components of the environmental response-related characteristics are each converted into feature vectors. Then, each feature vector is weighted according to the weights assigned in step S127. Finally, the weighted feature vectors are concatenated along the feature dimensions to form a preliminary fused feature set. Redundant features are removed from the preliminary fused feature set using the variance inflation factor method. The variance inflation factor value for each feature is calculated. The variance inflation factor value measures the degree of multicollinearity between the feature and other features; a larger variance inflation factor value indicates more severe multicollinearity. A variance inflation factor threshold is set. When the variance inflation factor value of a feature exceeds this threshold, it is considered a redundant feature and removed, retaining only the feature elements directly related to pollution diffusion and possessing unique identifying characteristics.

[0071] Step S129: Organize the feature elements of optical features, morphological correlation features, and environmental response correlation features after removing redundancy according to the categories of optical features, morphological correlation features, and environmental response correlation features to form a cross-verification feature set.

[0072] In this embodiment, after redundant feature removal, the remaining feature elements are optical features, morphological features, and environmental response features. These feature elements are then categorized according to their respective feature types: optical features are grouped into one category, morphological features into another, and environmental response features into yet another. Each category contains multiple specific feature parameters; for example, optical features include average reflectance spectrum curves and average reflectance intensity values; morphological features include contour morphology models and distribution density data; and environmental response features include color change features and texture variation features. These categorized feature elements are then combined to form a cross-validation feature set.

[0073] Step S130: Based on the spatiotemporal evolution correlation data of the cross-verification feature set, construct a pollution feature reverse source tracing path model. The pollution feature reverse source tracing path model derives the starting correlation node of pollution diffusion based on the feature evolution trajectory.

[0074] In this embodiment, a reverse source tracing path model for pollution features is constructed based on spatiotemporal evolution correlation data of cross-verification feature sets. First, state parameters of each feature element at different time points are extracted from the cross-verification feature set. These state parameters include quantifiable descriptive attributes such as feature intensity, distribution range, and morphological details. The state parameters of each feature element are sorted according to the chronological order of time points, forming an independent evolution sequence for each feature element. All independent evolution sequences together constitute a feature evolution trajectory set. The evolution trend of each feature element in the feature evolution trajectory set is analyzed to identify evolution directions characterized by increasing feature intensity, expanding distribution range, and rich morphological details, thus determining the forward evolution path of pollution diffusion. Using the endpoint of the forward evolution path as the starting point for reverse source tracing, the state parameters of each feature element at historical time points are traced in reverse chronological order to construct a reverse evolution trajectory. The correlation of the state parameters of each feature element in the reverse evolution trajectory is verified to confirm whether there is spatial overlap or adjacency between the state parameters of different feature elements, and to verify whether the trends of the state parameters of each feature element over time are synchronized at the time points, eliminating abnormal trajectory points. Based on the validated reverse evolution trajectory, key time nodes where the state parameters of each feature element tend to stabilize and are interconnected are extracted. The magnitude of feature changes between key time nodes is analyzed, and key time nodes that meet preset change criteria are identified as transition nodes. Using the transition nodes as a benchmark, the state parameters of each feature element are traced in reverse order until the state parameters of all feature elements converge within the same spatial coordinate range and the magnitude of change tends to level off. The combination of feature state parameters corresponding to this spatial coordinate range is taken as the starting correlation node for pollution diffusion. Based on the starting correlation node, transition nodes, and feature state parameters of each key time node, a pollution feature reverse source tracing path model with the reverse evolution trajectory as its core is constructed. This pollution feature reverse source tracing path model includes feature transmission relationships and evolutionary correlation data between each node.

[0075] Step S131: Extract the state parameters of each feature element in the cross-verification feature set at different time nodes. The quantifiable descriptive attributes of the feature intensity, distribution range, and morphological details are the state parameters of each feature element.

[0076] In this embodiment, the state parameters of each feature element at different time nodes are extracted from the cross-verification feature set. The cross-verification feature set includes multiple feature elements such as the optical characteristics of pollutants, the correlation features of diffusion patterns, and the correlation features of environmental responses. For each feature element, there are corresponding state descriptions at different time nodes, and these state descriptions are reflected by quantifiable parameters. For example, the state parameters of the optical characteristics of pollutants include the average reflection intensity value, the change amount of the reflection intensity in different wavelength bands, etc.; the state parameters of the correlation features of diffusion patterns include the distribution density, the extension amplitude, etc.; the state parameters of the correlation features of environmental responses include the color difference value, the texture variation degree, etc. These parameters can quantitatively describe the properties such as the intensity, distribution range, and morphological details of the feature elements at different time nodes.

[0077] Step S132: Sort the state parameters of each feature element in the order of the time nodes to form an independent evolution sequence for each feature element, and the independent evolution sequences of all feature elements together constitute the feature evolution trajectory set.

[0078] In this embodiment, for each feature element, the state parameters at different time nodes are arranged in the order of the time nodes to form an ordered sequence, that is, an independent evolution sequence. For example, if the state parameters of a certain feature element at time nodes T1, T2, T3 (T1 < T2 < T3) are P1, P2, P3 respectively, then its independent evolution sequence is [P1, P2, P3]. Combine the independent evolution sequences of all feature elements together to constitute the feature evolution trajectory set. This feature evolution trajectory set contains the trajectory information of all feature elements changing with time.

[0079] Step S133: Analyze the evolution trends of each feature element in the feature evolution trajectory set, identify the evolution directions in which the feature intensity increases, the distribution range expands, and the morphological details become rich, and determine the positive evolution path of pollution diffusion.

[0080] Step S1331: Arrange each feature element in the feature evolution trajectory set in the order of time nodes to form a feature element time series, and each feature element time series contains the state parameters of this feature element at all captured time points.

[0081] In this embodiment, each feature element in the feature evolution trajectory set has state parameters at multiple time nodes. Arrange the above state parameters in the order of the time nodes to form a feature element time series. Each time series is an ordered set of state parameters, containing the state information of this feature element at all captured time points. For example, if the state parameters of a certain feature element at time points t1, t2, t3... tn are s1, s2, s3... sn respectively, then its time series is [s1, s2, s3... sn].

[0082] Step S1332: Extract the intensity parameter sequence from the time series of each feature element. The intensity parameter sequence is the feature intensity data arranged in chronological order. By comparing the intensity parameter values ​​at adjacent time points, determine whether the intensity change trend is increasing, decreasing or remaining stable.

[0083] In this embodiment, intensity parameters are extracted from the time series of each feature element. These intensity parameters are quantitative indicators describing the strength of the feature element, such as the average reflectance intensity value of the optical characteristics of pollutants and the color difference value of environmental response-related features. These intensity parameters are arranged in chronological order to form an intensity parameter sequence. By comparing the intensity parameter values ​​at adjacent time points, the difference between the intensity parameter value at a later time point and the intensity parameter value at a previous time point is calculated. If the difference is positive and remains positive, the intensity trend is increasing; if the difference is negative and remains negative, the intensity trend is decreasing; if the difference fluctuates around zero, the intensity trend is stable.

[0084] Step S1333: Extract the distribution range parameter sequence from the time series of each feature element. The distribution range parameter sequence contains a set of feature coverage spatial coordinates arranged in chronological order. By comparing the spatial coordinate sets of adjacent time points, calculate the expansion area or volume change of the distribution range and determine whether the trend of the distribution range change is expansion, contraction or stabilization.

[0085] In this embodiment, distribution range parameters are extracted from the time series of each feature element. These parameters describe the spatial coverage of the feature element and are typically represented as a set of spatial coordinates. These distribution range parameters are arranged chronologically to form a distribution range parameter sequence. By comparing the spatial coordinate sets of adjacent time points, the area or volume difference between the distribution range at the later time point and the distribution range at the previous time point is calculated. If the difference is positive, it indicates that the distribution range is expanding; if the difference is negative, it indicates that the distribution range is shrinking; if the difference is close to zero, it indicates that the distribution range remains stable.

[0086] Step S1334: Extract the morphological detail parameter sequence from the time series of each feature element. The morphological detail parameter sequence contains data on feature contour details and internal structure details arranged in chronological order. By comparing the morphological detail parameters at adjacent time points, determine whether the trend of morphological detail changes is enrichment, simplification, or stability.

[0087] In this embodiment, morphological detail parameters are extracted from the time series of each feature element. These parameters describe detailed information such as the contour shape and internal structure of the feature element, such as the contour morphology model parameters and texture feature parameters of diffusion morphology associated features. These morphological detail parameters are arranged in chronological order to form a morphological detail parameter sequence. By comparing the morphological detail parameters at adjacent time points, the changes in the complexity of the feature contour and the refinement of the internal structure are analyzed. If the contour complexity increases and the internal structure becomes more refined, the trend of morphological detail change is enrichment; if the contour complexity decreases and the internal structure becomes simpler, the trend of morphological detail change is simplification; if the morphological detail parameters do not change significantly, they remain stable.

[0088] Step S1335: Select characteristic elements whose intensity parameter sequence shows an increasing trend, whose distribution range parameter sequence shows an expanding trend, and whose morphological detail parameter sequence shows an enriching trend. Characteristic elements whose intensity parameter sequence shows an increasing trend, whose distribution range parameter sequence shows an expanding trend, and whose morphological detail parameter sequence shows an enriching trend are characteristic elements directly related to pollution diffusion.

[0089] In this embodiment, based on the analysis results of steps S1332, S1333 and S1334, feature elements that simultaneously satisfy the following conditions are selected: the intensity parameter sequence shows an increasing trend, the distribution range parameter sequence shows an expanding trend, and the morphological detail parameter sequence shows an enriching trend. The changing trends of the above feature elements are consistent with the process of pollutant diffusion in the environment. That is, as time goes by, the concentration (intensity) of pollutants increases, the diffusion range expands, and the morphological structure becomes more complex. Therefore, they are considered to be feature elements directly related to pollution diffusion.

[0090] Step S1336: Perform a synchronization analysis on the time series of characteristic elements whose intensity parameter sequence shows an increasing trend, whose distribution range parameter sequence shows an expanding trend, and whose morphological detail parameter sequence shows an enriching trend, to determine whether the starting time points of the increasing intensity, expanding distribution range, and enriching morphological details of each characteristic element are consistent, and whether the rate of change has a synergistic relationship.

[0091] In this embodiment, for the selected characteristic elements directly related to pollution diffusion, their time series are analyzed for synchronicity. First, the starting time points for each characteristic element's increasing intensity, expanding distribution range, and richer morphological details are determined—that is, the points at which the trend begins to change. These starting time points are compared for consistency; if the starting time points of each characteristic element are within a preset time error range, they are considered consistent. Then, the rate of change for each characteristic element is calculated, such as the rate of increase in intensity, the rate of expansion of distribution range, and the rate of enrichment of morphological details. The synergistic relationship between these rates of change is analyzed, for example, whether the magnitude and trend of the rate of change are consistent, or whether they increase or decrease simultaneously.

[0092] Step S1337: Based on the synchronicity analysis results, feature elements with consistent starting time points and synergistic change rates are grouped into the same evolutionary group. Each evolutionary group represents a set of interrelated feature evolution trends during the pollution diffusion process.

[0093] In this embodiment, based on the synchronicity analysis results, feature elements with consistent starting times and synergistic change rates are grouped into the same evolutionary group. Feature elements within the same evolutionary group exhibit similar evolutionary behaviors during pollution diffusion; they are interconnected and mutually influential, collectively reflecting a certain aspect of pollution diffusion characteristics. For example, an evolutionary group may include both optical characteristics of pollutants and diffusion morphology-related characteristics. Their starting times are consistent, and their intensity increase rate and distribution range expansion rate show synchronous changes, indicating that these two feature elements are closely related during pollution diffusion.

[0094] Step S1338: Extract the evolution rate data of each characteristic element in each evolution group. The evolution rate data includes the intensity increase rate, the distribution range expansion rate, and the morphological detail enrichment rate. Analyze the magnitude and trend of the intensity increase rate, distribution range expansion rate, and morphological detail enrichment rate in each evolution group to assess the importance of the evolution group in the pollution diffusion process.

[0095] In this embodiment, for each evolutionary group, the evolution rate data of each characteristic element is extracted, including the intensity increase rate, the distribution range expansion rate, and the morphological detail enrichment rate. The intensity increase rate is calculated by the change in intensity parameters per unit time; the distribution range expansion rate is calculated by the change in the area or volume of the distribution range per unit time; and the morphological detail enrichment rate is calculated by the change in morphological detail parameters per unit time. The magnitude of these evolution rates is analyzed; a larger value indicates a more significant change in the corresponding aspect of the evolutionary group. Simultaneously, the trend of the evolution rate is analyzed to determine whether it is accelerating, decelerating, or remaining stable. Based on the magnitude and trend of the evolution rate, the importance of the evolutionary group in the pollution diffusion process is assessed; evolutionary groups with large and continuously increasing evolution rates are generally more important.

[0096] Step S1339: Sort the evolution groups according to the magnitude of the comprehensive evolution rate, and select the evolution groups whose comprehensive evolution rate meets the preset rate standard and whose coverage meets the preset range standard as the main evolution group. The main evolution group dominates the overall evolution process of pollution diffusion.

[0097] In this embodiment, the comprehensive evolution rate of each evolutionary group is calculated. The comprehensive evolution rate is a weighted sum of the intensity increase rate, the distribution range expansion rate, and the morphological detail enrichment rate, with the weights set according to the importance of each rate. All evolutionary groups are sorted according to the magnitude of the comprehensive evolution rate, and evolutionary groups with a comprehensive evolution rate greater than a preset rate standard are selected. Simultaneously, the coverage of these evolutionary groups is checked to see if it meets a preset range standard, which is evaluated using the distribution range parameter of the characteristic elements within the evolutionary group. Evolutionary groups that simultaneously meet both the comprehensive evolution rate standard and the coverage standard are selected as the primary evolutionary group. This primary evolutionary group plays a dominant role in the pollution diffusion process, and its evolution trend represents the overall evolution process of pollution diffusion.

[0098] Step S13310: Based on the time series, change trend and comprehensive evolution rate of each characteristic element in the main evolution group, construct a continuous evolution path from the initial state to the current state. This continuous evolution path serves as the positive evolution path for pollution diffusion.

[0099] In this embodiment, a forward evolution path for pollution diffusion is constructed based on the time series, changing trends, and overall evolution rate of each characteristic element in the main evolutionary group. First, the initial and current states of each characteristic element in the main evolutionary group are determined. The initial state is the state parameter corresponding to the first time node in the time series, and the current state is the state parameter corresponding to the last time node. Then, based on the changing trends and overall evolution rate of each characteristic element, continuous change curves from the initial state to the current state are fitted. These curves together constitute the forward evolution path for pollution diffusion. This forward evolution path describes the process by which pollutants, starting from the initial state, continuously increase in intensity, expand in distribution range, and enrich in morphological details over time.

[0100] Step S134: Using the end point of the forward evolution path as the starting point of the reverse tracing, trace the state parameters of each feature element in reverse order of historical time nodes to construct the reverse evolution trajectory.

[0101] In this embodiment, the endpoint of the forward evolution path, i.e., the feature state parameters at the current time node, is used as the starting point for reverse tracing. The state parameters of each feature element are traced sequentially from back to front, across historical time nodes. For example, if the time node order of the forward evolution path is T1, T2, T3, T4 (T1 being the initial time, T4 being the current time), then the time node order of the reverse tracing is T4, T3, T2, T1. For each historical time node, the state parameters of each feature element are extracted, and these state parameters are arranged in reverse order to construct a reverse evolution trajectory. The reverse evolution trajectory reflects the change process of pollution features from the current state to the historical state.

[0102] Step S135: Perform correlation verification on the state parameters of each feature element in the reverse evolution trajectory. Based on the spatial distribution range of each feature element at the same historical time node, confirm whether there is spatial overlap or adjacency relationship between the state parameters of different feature elements, and verify whether the trend of the state parameters of each feature element changes with time is synchronized at the time node. Eliminate abnormal trajectory points that are not spatially overlapping and not synchronized in time.

[0103] In this embodiment, the correlation of the state parameters of each feature element in the reverse evolution trajectory is verified. First, for the same historical time point, the spatial distribution range parameters of each feature element are extracted, and the overlapping area or adjacent distance of the spatial distribution ranges of different feature elements is calculated. If the overlapping area is greater than a preset area threshold or the adjacent distance is less than a preset distance threshold, it is considered that different feature elements have an overlapping or adjacent relationship in spatial location. Second, it is verified whether the trend of the change of the state parameters of each feature element over time is synchronized at the time point, that is, whether the change direction and magnitude of the state parameters of each feature element are consistent at adjacent time points. For trajectory points that do not overlap spatially and are not synchronized temporally, they are regarded as abnormal trajectory points and removed to ensure the accuracy and reliability of the reverse evolution trajectory.

[0104] Step S136: Based on the verified reverse evolution trajectory, extract the key time nodes where the state parameters of each feature element tend to be stable and interrelated, and the important transition stages in the pollution diffusion process correspond to each key time node.

[0105] In this embodiment, the changes in the state parameters of each feature element are analyzed based on the verified reverse evolution trajectory. When the change in the state parameter of a feature element is less than a preset stability threshold and remains stable at multiple consecutive time points, the state parameter of that feature element is considered to be stabilizing. Simultaneously, the correlation between different feature elements is examined. When multiple feature elements all tend to stabilize at the same time point and are interconnected, that time point is determined as a critical time point. Critical time points correspond to important transitional stages in the pollution diffusion process, such as the time point when pollution diffusion transitions from a rapid diffusion stage to a stable stage.

[0106] Step S137: Analyze the characteristic change range between key time nodes, and lock the key time nodes whose characteristic change range meets the preset change standard as transition nodes. The dividing point where pollution diffusion changes from the initial state to the rapid diffusion state is the time node.

[0107] In this embodiment, the change range of each characteristic element's state parameter between adjacent key time nodes is calculated. The change range is measured by subtracting the absolute value of the state parameter value of the previous time node from the value of the state parameter at the later time node. The change range is compared with a preset change standard. When the change range is greater than the preset change standard, it indicates that the characteristic element has changed significantly within that time period. Key time nodes whose characteristic change range meets the preset change standard are identified as transition nodes. These transition nodes are the dividing points where pollution diffusion transitions from the initial state to a rapid diffusion state. Before this node, pollution diffusion is relatively slow; after this node, pollution begins to diffuse rapidly.

[0108] Step S138: Using the transition node as a reference, continue to trace the state parameters of each feature element in reverse order until the state parameters of all feature elements converge in the same spatial coordinate range and the change amplitude tends to be gradual.

[0109] Step S1381: Take the time point corresponding to the transition node as the reference time point for reverse tracing, and extract the state parameters of all feature elements at the reference time point. The intensity parameter, distribution range parameter, and morphological parameter are the state parameters of all feature elements at the reference time point.

[0110] In this embodiment, the time point corresponding to the transition node is determined as the reference time point for reverse tracing. At this reference time point, the state parameters of all feature elements are extracted, including intensity parameters (such as average reflection intensity value), distribution range parameters (such as spatial coordinate set), and morphological parameters (such as contour morphology model parameters), etc. The above state parameters are used as the starting data for reverse tracing.

[0111] Step S1382: Starting from the reference time point, extract the state parameters of all feature elements of the previous time point in reverse chronological order to form a parameter sequence for reverse tracing.

[0112] In this embodiment, starting from the reference time point, the state parameters of all feature elements are extracted sequentially in reverse chronological order from the previous time point, the two time points before that, and so on. The above state parameters are arranged in chronological order to form a reverse tracing parameter sequence, which reflects the state changes of feature elements from the reference time point to earlier time points.

[0113] Step S1383: Compare the changes in the state parameters of each feature element between the reference time point and the previous time point, and calculate the intensity change difference, distribution range offset, and morphological difference.

[0114] In this embodiment, for each time node in the reverse-tracing parameter sequence, the change in the state parameters of each feature element is calculated by comparing it with that of the reference time point. The intensity change difference is the intensity parameter value of the previous time node minus the intensity parameter value of the reference time point; the distribution range offset is obtained by calculating the offset distance between the center coordinates of the distribution range of the previous time node and the reference time point; the morphological difference is obtained by comparing the degree of difference between the morphological parameters (such as the control vertex coordinates of the contour morphological model) of the previous time node and the reference time point, for example, by calculating the sum of the distances of the corresponding control vertices.

[0115] Step S1384: Analyze the changing trend of each characteristic element, and determine whether the change is gradually increasing, gradually decreasing, or remaining stable. If the change is gradually decreasing, it indicates that the characteristic element is converging towards the initial state.

[0116] In this embodiment, the changing trends of each feature element during the reverse tracing process are analyzed. As time progresses in reverse chronological order, the numerical changes in intensity variation, distribution range offset, and morphological difference are observed. If the values ​​of these changes gradually decrease, it indicates that the state parameters of the feature elements are getting closer and closer to the state parameters at the reference time point, meaning that the feature elements are converging towards their initial state.

[0117] Step S1385: Extract the distribution range parameters of each feature element in the reverse traced parameter sequence, and calculate the intersection area of ​​the distribution range of all feature elements at each time node. The intersection area is the spatial coordinate range covered by all feature elements.

[0118] In this embodiment, for each time point in the reverse-tracing parameter sequence, the distribution range parameters of each feature element are extracted, i.e., the set of spatial coordinates. The intersection region of the distribution ranges of all feature elements is calculated; the intersection region is the spatial coordinate range simultaneously covered by all feature elements. This is obtained by calculating the intersection of multiple spatial coordinate sets. For example, for the distribution ranges A and B of two feature elements, their intersection region is A∩B.

[0119] Step S1386: Analyze the changes of the intersection region over time in reverse order, calculate the area and boundary coordinates of the intersection region at each time node, and determine whether the intersection region gradually expands, gradually shrinks, or remains stable.

[0120] In this embodiment, the area and boundary coordinates of the intersection region at each time point are calculated. The area is obtained by multiplying the number of pixels in the intersection region by the actual area of ​​a single pixel; the boundary coordinates are obtained by extracting the coordinates of the boundary pixels of the intersection region. The changes of these areas and boundary coordinates over time are analyzed. If the area of ​​the intersection region gradually decreases and the boundary coordinates gradually concentrate towards a certain central position, it indicates that the intersection region is gradually shrinking; if the area gradually increases, it indicates that the intersection region is gradually expanding; if the area and boundary coordinates do not change significantly, it indicates that the intersection region remains stable.

[0121] Step S1387: If the intersection area gradually shrinks and the boundary tends to be concentrated, then continue to trace the characteristic element state parameters of the previous time node in reverse order, and repeatedly calculate the intersection area and change.

[0122] In this embodiment, if the intersection region shows a trend of gradually shrinking and the boundary tending to concentrate, it indicates that the distribution range of each feature element is converging towards a common central region, and it is necessary to continue to trace the feature element state parameters of the previous time node in reverse order. Repeat the operations of steps S1385 and S1386 to calculate the new intersection region and change amount until the trend of change of the intersection region changes.

[0123] Step S1388: When tracing back to a certain time node, if the changes in the state parameters of all feature elements meet the preset stability criteria, and the area of ​​the intersection region no longer shrinks significantly and the boundary remains stable, then stop tracing and extract the state parameters of all feature elements at the corresponding stop time node. Confirm that the intensity parameters, distribution range parameters, and morphological parameters of each feature element are in a stable state, and the intersection region of all distribution range parameters is clear.

[0124] In this embodiment, when tracing back to a certain time point, it is checked whether the changes in the state parameters of all feature elements are less than a preset stability threshold, and whether the area of ​​the intersection region fluctuates within a preset area stability range, and whether the boundary coordinates remain stable. If these conditions are met, tracing stops, and this time point is the stop time point. The state parameters of all feature elements at the stop time point are extracted, confirming that the intensity parameters, distribution range parameters, and morphological parameters of each feature element are all in a stable state, and that the intersection region of all distribution range parameters is clear, i.e., the intersection region has clear boundaries and a definite area.

[0125] Step S1389: Determine the intersection area corresponding to the stop time node as the spatial coordinate range where all characteristic element state parameters converge. The spatial coordinate range where all characteristic element state parameters converge is the central spatial area of ​​the initial stage of pollution diffusion.

[0126] In this embodiment, the intersection region corresponding to the stopping time node is defined as the spatial coordinate range where the state parameters of all feature elements converge. This spatial coordinate range is the area jointly covered by all feature elements at the stopping time node. Since the state parameters of each feature element have tended to stabilize at this time, and the intersection region no longer shrinks significantly, this area is considered to be the central spatial region in the initial stage of pollution diffusion, that is, the core area where pollutants are initially released and diffused.

[0127] Step S139: Use the combination of characteristic state parameters corresponding to the spatial coordinate range as the starting correlation node for pollution diffusion. The starting correlation node includes the optical characteristic parameters of pollutants in the initial state, the diffusion morphology correlation characteristic parameters, and the environmental response correlation characteristic parameters.

[0128] In this embodiment, the combination of characteristic state parameters corresponding to the spatial coordinate range determined in step S1389 is used as the starting correlation node for pollution diffusion. The starting correlation node includes parameters of the optical characteristics of pollutants within the spatial coordinate range (such as average reflectance spectrum curve, average reflectance intensity value, etc.), parameters of diffusion morphology correlation characteristics (such as distribution density, contour morphology model, etc.), and parameters of environmental response correlation characteristics (such as color change characteristics, texture variation characteristics, etc.). The above parameters reflect the characteristic information of the initial state of pollution diffusion and are key nodes for pollution source tracing.

[0129] Step S1310: Based on the characteristic state parameters of the initial associated node, transition node and each key time node, construct a pollution feature reverse source tracing path model with the reverse evolution trajectory as the core. The pollution feature reverse source tracing path model includes the feature transmission relationship and evolution association data between the initial associated node, transition node and each key time node in the pollution feature reverse source tracing path model.

[0130] In this embodiment, a reverse source tracing path model for pollution characteristics is constructed based on the characteristic state parameters of the initial associated node, transition nodes, and each key time node. This model uses a reverse evolution trajectory as its core, connecting each node in chronological order to form a path tracing back from the current state to the initial associated node. The model includes the feature transmission relationships between nodes, i.e., how the characteristic state parameters of the previous node affect the characteristic state parameters of the next node; and evolutionary correlation data, such as the magnitude and rate of feature changes between nodes. Through this model, the process of pollution characteristics evolving from the initial associated node through each key time node and transition node to the current state can be demonstrated.

[0131] Step S140: Associate the starting node of the pollution feature reverse source tracing path model with the imaging spatial coordinates of the UAV remote sensing device to lock the initial area of ​​pollution spread.

[0132] In this embodiment, the spatial coordinate range corresponding to the initial associated node of the pollution feature reverse source tracing path model is mapped to the imaging spatial coordinates of the UAV remote sensing device. The spatial coordinate range of the initial associated node is defined in the feature coordinate system and needs to be transformed to the imaging spatial coordinate system of the UAV remote sensing device, which is usually based on the geographic coordinates (longitude and latitude) of the Global Positioning System. Through a coordinate transformation algorithm, the spatial coordinate range in the feature coordinate system is converted into a geographic coordinate range, obtaining the geographic coordinates of the initial area of ​​pollution diffusion. This geographic coordinate range is then identified as the area where pollution diffusion initially begins.

[0133] Step S150: Adjust the spatial boundary of the initial action area by using the dynamic evolution trajectory association data of the cross-verification feature set, and output the location result of the pollution source in the target monitoring area.

[0134] In this embodiment, the spatial boundary of the initial area of ​​action is adjusted by linking the dynamic evolution trajectory data of the cross-verification feature set. First, dynamic evolution trajectory data of the pollutant's optical characteristics, diffusion morphology correlation characteristics, and environmental response correlation characteristics at different time points are extracted from the cross-verification feature set, including changes in the spatial distribution range, intensity, and morphology of the features. These dynamic evolution trajectory data are analyzed to determine the expansion direction and rate of the features. Based on this information, the spatial boundary of the initial area of ​​action is adjusted to better reflect the actual situation of pollution diffusion. Finally, the adjusted spatial boundary parameters are integrated to output the location result of the pollution source in the target monitoring area, which includes information such as the center coordinates, coverage area, and boundary contour of the pollution source.

[0135] Step S151: Extract the spatial distribution range data of the optical features of pollutants in the cross-verification feature set at different time points. The spatial coordinate set and boundary contour parameters of the feature coverage are the spatial distribution range data of the optical features of pollutants at different time points.

[0136] In this embodiment, spatial distribution data of pollutant optical features at different time points are extracted from the cross-verification feature set. This data includes the set of spatial coordinates of the feature coverage, i.e., the spatial coordinates of all pixels covered by the feature at that time point; and boundary contour parameters, i.e., parameters describing the shape of the feature boundary, such as the coordinate sequence of boundary pixels and fitting parameters of the contour curve. These data allow us to understand the spatial distribution of pollutant optical features at different time points.

[0137] Step S152: Analyze the spatial distribution range data of the optical characteristics of pollutants over time to determine the expansion direction and expansion rate of the spatial distribution of the characteristics. The expansion direction is determined by the boundary contour offset trend of adjacent time nodes, and the expansion rate is determined by the amount of change in the spatial distribution range per unit time.

[0138] In this embodiment, the spatial distribution data of the optical characteristics of pollutants at different time points are analyzed, and the changes in the spatial distribution range between adjacent time points are calculated. The expansion direction is determined by comparing the boundary contours of adjacent time points to determine the offset trend of the boundary contours. For example, the expansion direction is determined by calculating the centroid offset direction of the boundary contours. The expansion rate is determined by calculating the area change or boundary movement distance of the spatial distribution range per unit time. The area change is the area of ​​the later time point minus the area of ​​the previous time point, divided by the time interval. The boundary movement distance is the average movement distance of corresponding points on the boundary of adjacent time points, divided by the time interval.

[0139] Step S153: Based on the expansion direction and expansion rate, reversely deduce the original distribution range of the optical features of the pollutants within the initial action area. The original distribution range is the spatial coverage area when the features have not expanded significantly.

[0140] Step S1531: Convert the extension direction data of the optical features of pollutants into spatial direction indicators. Each extension direction corresponds to a spatial direction, representing the spatial extension trend of the optical features of pollutants in that spatial direction.

[0141] In this embodiment, the expansion direction data of the optical features of pollutants are converted into spatial direction indicators. For example, the expansion direction is represented by angles, with 0 degrees representing true north, 90 degrees representing true east, and so on. Each expansion direction corresponds to a specific spatial direction, which represents the extension trend of the optical features of pollutants in that spatial direction, that is, the features mainly expand in that direction.

[0142] Step S1532: Starting from the current spatial distribution range boundary of the optical characteristics of pollutants, extend in the opposite direction of the expansion direction. The range of the reverse extension is estimated based on the historical distribution range change trend in this direction. The starting points of all the reverse extension lines converge to form the initial feature coverage boundary estimated based on the expansion direction and expansion rate. The initial feature coverage boundary estimated based on the expansion direction and expansion rate is the preliminary original distribution range boundary.

[0143] In this embodiment, starting from the current spatial distribution range boundary of the pollutant's optical features, a reverse extension is performed in the opposite direction for each expansion direction. The distance of the reverse extension is estimated based on the historical distribution range change trend in that direction. For example, the total expansion distance in that direction is calculated based on the expansion rate and expansion time, and then this total expansion distance is extended from the current boundary in the opposite direction. Connecting the starting points of the reverse extension lines of all expansion directions forms a preliminary initial feature coverage boundary, which is the preliminary boundary of the original distribution range.

[0144] Step S1533: Extract the spatial distribution range data of the optical characteristics of pollutants at different time points, arrange them in reverse chronological order, and obtain the distribution range sequence from the current state to the earlier state.

[0145] In this embodiment, spatial distribution data of the optical characteristics of pollutants at multiple time points are extracted. This data is then arranged in reverse chronological order, starting from the current time point and sequentially arranging the distribution data of the previous time point, the two time points before that, and so on, forming a distribution range sequence. This sequence reflects the spatial distribution changes of the optical characteristics of pollutants from the current state to earlier states.

[0146] Step S1534: Analyze the variation range of distribution range between adjacent time nodes in the distribution range sequence, and calculate the difference in distribution range between each time node and the previous time node.

[0147] In this embodiment, for each time node in the distribution range sequence, the difference in distribution range between it and the previous time node (i.e., an earlier time node) is calculated. The difference can be measured by calculating the area difference or boundary offset distance between the two distribution ranges. The area difference is the area of ​​the previous time node minus the area of ​​the current time node, and the boundary offset distance is the average distance between corresponding points on the boundary of the two time nodes.

[0148] Step S1535: Select continuous time node groups whose distribution range difference is less than the preset difference threshold as continuous time node groups with gradual changes in feature distribution range.

[0149] In this embodiment, a preset difference threshold is set. When the difference in distribution range between adjacent time nodes in the distribution range sequence is less than this threshold, the change in distribution range is considered to be gradual. Time node groups in which the difference in distribution range between multiple consecutive time nodes is less than the preset difference threshold are selected. These time node groups are considered to be continuous time node groups with gradual changes in the characteristic distribution range, indicating that the spatial distribution range of the optical characteristics of pollutants changes little and is in a relatively stable state during this time period.

[0150] Step S1536: Extract the spatial distribution range data corresponding to the earliest time node in the continuous time node group where the feature distribution range changes smoothly, and use it as the spatial range data of the early stable distribution of the optical features of the pollutant.

[0151] In this embodiment, among a group of consecutive time nodes where the feature distribution range changes gradually, the earliest time node is selected, and the spatial distribution range data corresponding to that time node is extracted. This data is used as the spatial range data of the early stable distribution of the optical features of pollutants. This spatial range data represents the stable distribution range of the optical features of pollutants before they begin to expand significantly.

[0152] Step S1537: Spatially compare the initial feature coverage boundary estimated based on the expansion direction and expansion rate with the spatial distribution range data corresponding to the earliest time node in the stable phase, and calculate the spatial overlap and boundary deviation between the two.

[0153] In this embodiment, the initially estimated feature coverage boundary obtained in step S1532 is spatially compared with the spatial range data of the early stable distribution obtained in step S1536. The spatial overlap is calculated as the ratio of the area of ​​the overlapping region of the two ranges to the total area of ​​the two ranges; the boundary deviation is calculated as the average distance between corresponding points on the boundaries of the two ranges. The accuracy of the initially estimated initial coverage boundary is evaluated by the spatial overlap and boundary deviation.

[0154] Step S1538: Based on the spatial overlap and boundary deviation, adjust the initial feature coverage boundary estimated based on the expansion direction and expansion rate so that the adjusted initial feature coverage boundary estimated based on the expansion direction and expansion rate matches the spatial distribution range data corresponding to the earliest time node in the stable phase to the greatest extent.

[0155] In this embodiment, the initially estimated feature coverage boundary is adjusted based on spatial overlap and boundary deviation. If the spatial overlap is low or the boundary deviation is large, the initially estimated coverage boundary is moved or its shape is adjusted to bring it closer to the spatially stable distribution data from earlier periods, until the spatial overlap is maximized and the boundary deviation is minimized, achieving maximum fit.

[0156] Step S1539: The spatial region enclosed by the initial coverage boundary of the features, which is preliminarily estimated based on the expansion direction and expansion rate, is determined as the original distribution range of the optical features of pollutants within the initial action area. The original distribution range of the optical features of pollutants within the initial action area retains the spatial coverage state when the features have not expanded significantly.

[0157] In this embodiment, the spatial region enclosed by the adjusted preliminary estimated feature initial coverage boundary is determined as the original distribution range of the pollutant's optical features within the initial action area. This original distribution range preserves the spatial coverage state of the pollutant's optical features before significant expansion, i.e., the initial spatial range of the pollutant's distribution in the environment.

[0158] Step S154: Extract the initial morphological parameters of the diffusion morphology associated features in the cross-verification feature set. The initial morphological parameters include the initial contour shape of the feature, the location of the central region, and the distribution density gradient.

[0159] In this embodiment, initial morphological parameters of diffusion morphology-related features are extracted from the cross-verification feature set. Initial morphological parameters refer to the morphological description parameters of diffusion morphology-related features in the initial stage, including the initial contour shape, described by parameters of the contour morphology model (such as control vertex coordinates, curve order, etc.); the position of the central region, i.e., the center coordinates of the feature distribution; and the distribution density gradient, i.e., the rate of change of feature density from the central region to the boundary, reflecting the uniformity of feature distribution and diffusion trend.

[0160] Step S155: Spatial match the initial morphological parameters with the original distribution range derived in reverse, and adjust the boundary of the original distribution range so that the boundary precisely matches the position of the central region and the initial contour shape in the initial morphological parameters.

[0161] In this embodiment, the initial morphological parameters of the diffusion morphology associated features are spatially matched with the original distribution range derived through reverse derivation. First, the position of the central region in the initial morphological parameters is aligned with the center coordinates of the original distribution range to ensure that their center positions are consistent. Then, the boundary of the original distribution range is adjusted according to the initial contour shape to match the boundary shape of the original distribution range with the initial contour shape. By fine-tuning the pixel coordinates on the boundary of the original distribution range, the boundary is made to precisely match the initial contour shape, improving the accuracy of the original distribution range.

[0162] Step S156: Extract the initial response region data of the environmental response correlation features in the cross-verification feature set. The initial response region data is the spatial coordinate range when the environment first shows a pollution-related response.

[0163] In this embodiment, initial response region data of environmental response correlation features are extracted from the cross-verification feature set. Environmental response correlation features reflect the response of the environment (such as vegetation) to pollution. Initial response region data refers to the spatial coordinate range when the environment first shows a pollution-related response, that is, the area where response features such as vegetation color change and texture variation first appear. This area is usually located near the pollution source and is the initial area of ​​pollution impact.

[0164] Step S157: Overlay the initial response area data with the adjusted original distribution range to determine the overlapping area and the difference area between the initial response area data and the adjusted original distribution range. The overlapping area between the initial response area data and the adjusted original distribution range is the basic part of the initial action area, and the difference area between the initial response area data and the adjusted original distribution range is the edge part affected by environmental factors.

[0165] In this embodiment, the initial response region data and the adjusted original distribution range are spatially overlaid to calculate their intersection and difference. The intersection region, also known as the overlapping region, represents both the original distribution range of the pollutant's optical characteristics and the area where the environment first exhibits a pollution response; it forms the foundation of the initial action region. The difference region includes portions of the initial response region data that are not within the original distribution range and portions of the original distribution range that are not within the initial response region data. These differing regions may be marginal areas influenced by environmental factors (such as wind direction, topography, etc.).

[0166] Step S158: Based on the spatial coordinates and boundary parameters of the overlapping area between the initial response area data and the adjusted original distribution range, and combined with the distribution density gradient of the diffusion morphology correlation characteristics, further refine the boundary outline of the initial action area so that the boundary delineates the pollution-affected area and the edge-related area.

[0167] In this embodiment, the boundary contour of the initial affected area is further refined based on the spatial coordinates and boundary parameters of the overlapping region, combined with the distribution density gradient of the diffusion morphology correlation characteristics. The distribution density gradient reflects the change in the concentration of pollutants, and areas with large density gradients are usually the boundary areas of pollution impact. According to the changes in the distribution density gradient, the boundary contour of the overlapping region is adjusted, and the area with more significant pollution impact (the inner part of the area with large density gradient) is divided into the pollution impact area, while the outer part of the area with small density gradient is divided into the edge correlation area, so that the boundary can accurately distinguish these two areas.

[0168] Step S159: Integrate all adjusted spatial boundary parameters to form a set of positioning parameters including the coordinates of the central region, the range of the edge region, and the boundary contour equation.

[0169] For example, step S1591: Collect the spatial coordinate data of the center region of the adjusted initial action area. The center coordinates, vertex coordinates and axis coordinates of the center region are the spatial coordinate data of the center region of the adjusted initial action area. The center coordinates, vertex coordinates and axis coordinates of the center region are measured by a spatial positioning algorithm.

[0170] In this embodiment, spatial coordinate data of the central region of the adjusted initial action area are collected, including center coordinates, i.e., the latitude and longitude coordinates of the geometric center of the central region; vertex coordinates, i.e., the latitude and longitude coordinates of the corner points of the boundary of the central region; and axis coordinates, i.e., the endpoint coordinates of the major and minor axes of the central region. The above coordinate data are measured by a spatial positioning algorithm to ensure its accuracy.

[0171] Step S1592: Extract the boundary contour parameters of the central region. The boundary contour parameters include the coordinates of each feature point on the contour, the curvature change of the contour curve, and the tangent direction. Construct the boundary contour equation of the central region based on the coordinates of each feature point on the contour, the curvature change of the contour curve, and the tangent direction.

[0172] In this embodiment, boundary contour parameters of the central region are extracted, including the coordinates of each feature point on the contour. These feature points are representative points on the boundary, such as inflection points and extreme points. The curvature change of the contour curve, i.e., the curvature value of the boundary curve at each feature point, is also extracted. The tangent direction, i.e., the tangent angle of the boundary curve at each feature point, is also extracted. Based on these parameters, a curve fitting algorithm is used to construct the boundary contour equation of the central region. This boundary contour equation can accurately describe the boundary shape of the central region.

[0173] Step S1593: Collect spatial description data of the edge region range. The boundary distance between the edge region and the center region, the outer boundary coordinates of the edge region, and the area parameter of the edge region are the spatial description data of the edge region range.

[0174] In this embodiment, spatial description data of the edge region range is collected, including the boundary distance between the edge region and the central region, i.e., the average distance from the inner boundary of the edge region to the boundary of the central region; the outer boundary coordinates of the edge region, i.e., the coordinates of the feature points on the outer boundary of the edge region; and the area parameter of the edge region, i.e., the total area of ​​the edge region. The above data are used to describe the spatial range of the edge region and its positional relationship with the central region.

[0175] Step S1594: Perform curve fitting on the outer boundary coordinates of the edge region to generate the outer boundary contour equation of the edge region. The outer boundary contour equation and the boundary contour equation of the central region together constitute complete boundary description data.

[0176] In this embodiment, curve fitting is performed on the outer boundary coordinates of the edge region using the same curve fitting algorithm as the boundary contour equation of the central region to generate the outer boundary contour equation of the edge region. The outer boundary contour equation and the boundary contour equation of the central region together constitute the complete boundary description data of the initial action area, which can comprehensively describe the spatial shape of the initial action area.

[0177] Step S1595: Extract the feature parameters related to the spatial boundary from the cross-verification feature set. The boundary intensity threshold of the optical characteristics of pollutants, the boundary density threshold of the diffusion morphology associated characteristics, and the boundary response threshold of the environmental response associated characteristics are the feature parameters related to the spatial boundary from the cross-verification feature set.

[0178] In this embodiment, feature parameters related to spatial boundaries are extracted from the cross-verification feature set, including the boundary intensity threshold of the optical characteristics of pollutants, i.e., the critical value of optical intensity that distinguishes polluted areas from non-polluted areas; the boundary density threshold of the diffusion morphology associated features, i.e., the critical value of density that distinguishes polluted diffusion areas from non-diffusion areas; and the boundary response threshold of the environmental response associated features, i.e., the critical value of response degree that distinguishes environmental response areas from non-response areas. The above thresholds are used to help determine the boundary of the initial action area.

[0179] Step S1596: Compare the boundary intensity threshold of the optical characteristics of pollutants, the boundary density threshold of the diffusion morphology associated characteristics, and the boundary response threshold of the environmental response associated characteristics with the actual measured values ​​of the corresponding characteristics in the central region and the edge region, respectively. Based on the threshold comparison results, determine the boundary division between the central region and the edge region, so that the boundary divides different regions whose characteristic intensity meets the threshold standard.

[0180] In this embodiment, the extracted boundary strength threshold, boundary density threshold, and boundary response threshold are compared with the actual measured values ​​of the corresponding features in the central and edge regions, respectively. For the central region, the actual measured values ​​of its internal features are ensured to be higher than the corresponding thresholds; for the edge region, the actual measured values ​​of its internal features are near or slightly lower than the thresholds. Based on the comparison results, the boundary division between the central and edge regions is adjusted so that the boundary can accurately divide different regions whose feature strengths meet the threshold criteria.

[0181] Step S1597: Calculate the spatial relationship parameters between the central region and the edge region. The overlap ratio, relative orientation and distance distribution of the central region and the edge region are the spatial relationship parameters between the central region and the edge region. The overlap ratio, relative orientation and distance distribution of the central region and the edge region are used to describe the spatial structure of the pollution source location.

[0182] In this embodiment, the spatial positional relationship parameters of the central region and the edge region are calculated, including the overlap ratio, which is the ratio of the area of ​​the overlapping part of the central region and the edge region to the total area of ​​the two; the relative orientation, which is the direction of the edge region relative to the central region, such as east, south, west, north, etc.; and the distance distribution, which is the statistical distribution of the distances from each point in the edge region to the center of the central region. The above parameters are used to describe the spatial structure of the pollution source location and reflect the relative position and spatial distribution characteristics of the central region and the edge region.

[0183] Step S1598: Integrate the coordinate data of the central region, the boundary contour equation of the central region, the range data of the edge region, the outer boundary contour equation of the edge region, the feature threshold correlation data, and the spatial positional relationship parameters to form a preliminary parameter set.

[0184] In this embodiment, the coordinate data of the central region, the boundary contour equation of the central region, the range data of the edge region, the outer boundary contour equation of the edge region, the feature threshold correlation data, and the spatial position relationship parameters collected and calculated in steps S1591 to S1597 are integrated to form a preliminary parameter set, which includes all the key parameters describing the spatial boundary of the initial action area.

[0185] Step S1599: Perform logical consistency verification on the preliminary parameter set, and sort all the verified parameters according to the categories of central region parameters, edge region parameters, boundary contour equations, feature association parameters and spatial relationship parameters to form a positioning parameter set.

[0186] In this embodiment, a logical consistency check is performed on the initial parameter set to examine whether there are any contradictions or unreasonable aspects among the parameters, such as whether the boundary contour equation is closed or whether the feature threshold matches the actual feature measurement value. Problems discovered during the check are corrected to ensure the logical consistency of the parameter set. Then, all the checked parameters are sorted and organized according to the categories of central region parameters, edge region parameters, boundary contour equation, feature association parameters, and spatial relationship parameters to form the final positioning parameter set.

[0187] Step S1510: Based on the set of positioning parameters, generate a spatial location description of the pollution source in the target monitoring area. This spatial location description includes the center coordinates, coverage area, and boundary contour of the pollution source, and finally outputs a complete positioning result.

[0188] In this embodiment, the center coordinates of the pollution source are determined based on the center region coordinate data in the positioning parameter set; the coverage area of ​​the pollution source is determined based on the range data of the center region and the edge region; and the boundary contour of the pollution source is drawn based on the boundary contour equation. Integrating the above information generates a spatial location description of the pollution source in the target monitoring area. This description indicates the specific location, size, and shape of the pollution source in the target monitoring area, ultimately outputting a complete positioning result.

[0189] In one exemplary embodiment, a pollution source localization system combining UAV remote sensing and AI image recognition is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2As shown, this pollution source localization system combining UAV remote sensing and AI image recognition includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a pollution source localization method combining UAV remote sensing and AI image recognition. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of a pollution source positioning system that combines drone remote sensing and AI image recognition, or an external keyboard, touchpad, or mouse, etc.

[0190] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for locating pollution sources combining UAV remote sensing and AI image recognition, characterized in that, The method includes: Multi-dimensional collaborative imaging is performed using UAV remote sensing equipment to generate a multi-dimensional cross-image stream of the target monitoring area. The multi-dimensional cross-image stream includes surface reflection images, atmospheric transmission images, and vegetation response images. The surface reflection images, atmospheric transmission images, and vegetation response images are associated and bound according to the same spatial coordinates and time nodes. An AI image recognition model is used to extract pollution features from the multi-dimensional cross-image stream, generating a cross-verification feature set that includes optical features of pollutants, diffusion morphology correlation features, and environmental response correlation features. Based on the spatiotemporal evolution correlation data of the cross-verification feature set, a pollution feature reverse source tracing path model is constructed. The pollution feature reverse source tracing path model derives the starting correlation node of pollution diffusion in reverse based on the feature evolution trajectory. The initial associated node of the pollution feature reverse source tracing path model is associated and mapped with the imaging spatial coordinates of the UAV remote sensing equipment to lock the initial area of ​​pollution spread. The spatial boundary of the initial action area is adjusted by the dynamic evolution trajectory association data of the cross-verification feature set, and the location result of the pollution source in the target monitoring area is output.

2. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 1, characterized in that, The AI ​​image recognition model is used to extract pollution features from the multi-dimensional cross-image stream, generating a cross-verification feature set that includes optical features of pollutants, diffusion morphology correlation features, and environmental response correlation features, including: The multidimensional cross-image stream is input into the multidimensional feature separation module of the AI ​​image recognition model. Feature channels are separated according to the modal attributes of the surface reflection image, atmospheric transmission image and vegetation response image to obtain the modal feature channels corresponding to the surface reflection image, the modal feature channels corresponding to the atmospheric transmission image and the modal feature channels corresponding to the vegetation response image. Optical property analysis is performed on the modal feature channels corresponding to the surface reflection image to extract optical parameters that differ from the conventional surface reflection characteristics. These optical parameters that differ from the conventional surface reflection characteristics constitute the basic elements of the optical characteristics of pollutants. Morphological structure analysis is performed on the modal feature channels corresponding to atmospheric transmission images to identify image regions that are continuously distributed and have a diffusion trend. The contour morphology, distribution density, and extension correlation data of the image region are extracted. The contour morphology, distribution density, and extension correlation data of the image region form the composition of diffusion morphological correlation features. The response characteristics of the modal feature channels corresponding to the vegetation response image are detected, and the features that differ from the normal vegetation response in terms of color change and texture variation of vegetation pixels in the image are captured. The features that differ from the normal vegetation response in terms of color change and texture variation of vegetation pixels in the image are used to construct the constituent elements of the environmental response association features. The cross-verification module of the AI ​​image recognition model is activated to identify the spatial location range corresponding to the basic elements of the optical characteristics of pollutants and the spatial location range corresponding to the constituent content of the diffusion morphology-related characteristics, and to confirm the overlapping area of ​​the two within the same spatial coordinate range; within the overlapping area, the basic element parameters of the optical characteristics of pollutants and the constituent content parameters of the diffusion morphology-related characteristics are extracted respectively, and the synergistic change relationship between the two in terms of intensity, distribution or morphology is analyzed to obtain the correlation intensity data. A temporal correlation analysis was conducted on the basic elements of the optical characteristics of pollutants corresponding to the overlapping regions and the constituent elements of the environmental response correlation characteristics to verify the consistency of the evolution of the basic elements of the optical characteristics of pollutants, the constituent elements of the diffusion morphology correlation characteristics, and the constituent elements of the environmental response correlation characteristics at the same time point. Based on the spatial location comparison results and temporal correlation analysis results, the basic elements of the optical characteristics of pollutants, the constituent contents of the diffusion morphology-related characteristics, and the constituent elements of the environmental response-related characteristics are weighted. The basic elements of the optical characteristics of pollutants are assigned the highest weight, and the constituent contents of the diffusion morphology-related characteristics and the constituent elements of the environmental response-related characteristics are assigned corresponding weights. The basic elements of the optical characteristics of pollutants, the components of the diffusion morphology-related characteristics, and the components of the environmental response-related characteristics are fused together to form a preliminary fused feature set. Redundant features are removed from the preliminary fused feature set, and the feature elements that are directly related to pollution diffusion and have unique identification characteristics are retained. The feature elements of optical features, morphological correlation features, and environmental response correlation features after removing redundancy are organized according to the categories of optical features, morphological correlation features, and environmental response correlation features to form a cross-verification feature set.

3. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 1, characterized in that, The pollution feature reverse source tracing path model is constructed based on the spatiotemporal evolution correlation data of the cross-verification feature set, including: Extract the state parameters of each feature element in the cross-verification feature set at different time nodes. The quantifiable descriptive attributes of the feature intensity, distribution range, and morphological details are the state parameters of each feature element. The state parameters of each feature element are sorted according to the order of time nodes to form an independent evolution sequence of each feature element. The independent evolution sequences of all feature elements together constitute a set of feature evolution trajectories. Analyze the evolution trend of each feature element in the feature evolution trajectory set, identify the evolution direction of increasing feature intensity, expanding distribution range and rich morphological details, and determine the positive evolution path of pollution diffusion; Using the end point of the forward evolution path as the starting point of the reverse tracing, the state parameters of each feature element at historical time nodes are traced in reverse order according to the time nodes to construct the reverse evolution trajectory. The correlation of state parameters of each feature element in the reverse evolution trajectory is verified. Based on the spatial distribution range of each feature element at the same historical time node, it is confirmed whether there is spatial overlap or adjacency relationship of the state parameters of different feature elements. It is also verified whether the trend of the state parameters of each feature element changes with time is synchronized at the time node, and abnormal trajectory points that are not spatially overlapping and not synchronized in time are eliminated. Based on the verified reverse evolution trajectory, key time nodes where the state parameters of each feature element tend to be stable and interrelated are extracted, and important transition stages in the pollution diffusion process correspond to each key time node. Analyze the characteristic changes between key time points, identify key time points whose characteristic changes meet the preset change criteria as transition points, and define the time point as the dividing point where pollution diffusion changes from the initial state to the rapid diffusion state. Using the transition node as a reference, continue to trace the state parameters of each feature element in reverse order until the state parameters of all feature elements converge in the same spatial coordinate range and the change rate tends to be flat. The combination of characteristic state parameters corresponding to this spatial coordinate range is used as the starting correlation node for pollution diffusion. The starting correlation node includes the optical characteristic parameters of pollutants in the initial state, the diffusion morphology correlation characteristic parameters, and the environmental response correlation characteristic parameters. Based on the characteristic state parameters of the initial associated node, transition node and each key time node, a pollution feature reverse source tracing path model with reverse evolution trajectory as the core is constructed. The pollution feature reverse source tracing path model includes the feature transmission relationship and evolution association data between the initial associated node, transition node and each key time node in the pollution feature reverse source tracing path model.

4. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 1, characterized in that, The process of adjusting the spatial boundary of the initial action area using the dynamic evolution trajectory association data of the cross-verification feature set, and outputting the location result of the pollution source in the target monitoring area, includes: Extract the spatial distribution range data of the optical features of pollutants in the cross-verification feature set at different time points. The set of spatial coordinates and boundary contour parameters of the feature coverage are the spatial distribution range data of the optical features of pollutants at different time points. By analyzing the spatial distribution range data of the optical characteristics of pollutants over time, the expansion direction and expansion rate of the spatial distribution of the characteristics are determined. The expansion direction is determined by the boundary profile offset trend of adjacent time nodes, and the expansion rate is determined by the amount of change in the spatial distribution range per unit time. Based on the direction and rate of expansion, the original distribution range of the optical features of pollutants in the initial area of ​​action is deduced in reverse. The original distribution range is the spatial coverage area when the features have not expanded significantly. Extract the initial morphological parameters of the diffusion morphology associated features in the cross-verification feature set. The initial morphological parameters include the initial contour shape of the feature, the location of the central region, and the distribution density gradient. Spatial matching is performed between the initial morphological parameters and the original distribution range derived in reverse, and the boundary of the original distribution range is adjusted so that the boundary precisely matches the position of the central region and the initial contour shape in the initial morphological parameters. Extract the initial response region data of the environmental response correlation features in the cross-verification feature set. The initial response region data is the spatial coordinate range when the environment first shows a pollution-related response. By overlaying the initial response area data with the adjusted original distribution range, the overlapping area and the difference area between the initial response area data and the adjusted original distribution range are determined. The overlapping area between the initial response area data and the adjusted original distribution range is the basic part of the initial action area, and the difference area between the initial response area data and the adjusted original distribution range is the edge part affected by environmental factors. Based on the spatial coordinates and boundary parameters of the overlapping area between the initial response area data and the adjusted original distribution range, combined with the distribution density gradient of diffusion morphology correlation characteristics, the boundary contour of the initial action area is further refined, so that the boundary delineates the pollution-affected area and the edge-related area. Integrate all adjusted spatial boundary parameters to form a set of positioning parameters that includes the coordinates of the central region, the range of the edge region, and the boundary contour equation; Based on the set of positioning parameters, a spatial location description of the pollution source in the target monitoring area is generated. This spatial location description includes the center coordinates, coverage area, and boundary outline of the pollution source, and finally outputs a complete positioning result.

5. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 2, characterized in that, The process of performing optical property analysis on the modal feature channels corresponding to the surface reflection image, and extracting optical parameters in the modal feature channels that differ from the conventional surface reflection characteristics, includes: Acquire reference data of surface reflectance in the uncontaminated state of the target monitoring area. This reference data includes conventional reflectance spectral parameters and reflectance intensity parameters of different surface types within the target monitoring area. The modal feature channels corresponding to the surface reflection image are split into pixels to obtain the reflection spectrum data and reflection intensity data of each pixel; Each pixel's reflectance spectral data is compared with the corresponding surface type's conventional reflectance spectral parameters in the surface reflectance reference data to identify pixels with differences in spectral curve shape, peak position, and valley position. The reflectance intensity data of pixels with differences in spectral curve shape, peak position, and valley position are compared with the conventional reflectance intensity parameters of the corresponding land surface type. Pixels with reflectance intensity deviating from the conventional range are extracted, and the pixels with reflectance intensity deviating from the conventional range constitute a set of suspected polluted pixels. Spatial clustering is performed on the pixels in the suspected contaminated pixel set, and pixels that are spatially adjacent and have consistent spectral and intensity differences are grouped into the same pixel cluster; Extract the average reflectance spectrum curve of each pixel cluster. This average reflectance spectrum curve is the result of averaging the reflectance spectrum data of all pixels in the cluster, and represents the overall spectral data characteristics of the pixel cluster. Calculate the average reflection intensity value for each pixel cluster. This average reflection intensity value is the statistical average of the reflection intensity data of all pixels in the cluster, reflecting the overall reflection intensity level of the pixel cluster. The difference between the average reflectance spectrum curve of each pixel cluster and the conventional reflectance spectrum parameters is analyzed, and the distribution range and magnitude of the difference bands are determined by spectral similarity analysis. Extract the change in reflection intensity data corresponding to the difference band. This change in reflection intensity data is the difference between the average reflection intensity value and the conventional reflection intensity parameter within the difference band. The average reflectance spectrum curve, average reflectance intensity value, differential band distribution range and differential amplitude, and differential band reflectance intensity data of each pixel cluster are integrated to form the basic elements constituting the optical characteristics of pollutants.

6. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 3, characterized in that, The analysis focuses on the evolutionary trends of each feature element in the set of analytical feature evolution trajectories, identifying evolutionary directions characterized by increasing feature intensity, expanding distribution range, and rich morphological details, and determining the positive evolutionary path of pollution diffusion, including: Each feature element in the feature evolution trajectory set is arranged in order of time nodes to form a feature element time series. Each feature element time series contains the state parameters of the feature element at all capture time points. Extract the intensity parameter sequence from the time series of each feature element. The intensity parameter sequence is the feature intensity data arranged in chronological order. By comparing the intensity parameter values ​​at adjacent time points, it can be determined whether the intensity change trend is increasing, decreasing or remaining stable. Extract the distribution range parameter sequence from the time series of each feature element. The distribution range parameter sequence contains a set of feature coverage spatial coordinates arranged in chronological order. By comparing the spatial coordinate sets of adjacent time points, calculate the expansion area or volume change of the distribution range and determine whether the trend of the distribution range change is expansion, contraction or stabilization. Extract the morphological detail parameter sequence from the time series of each feature element. The morphological detail parameter sequence contains data on feature contour details and internal structure details arranged in chronological order. By comparing the morphological detail parameters at adjacent time points, determine whether the trend of morphological detail changes is enrichment, simplification, or stability. Feature elements showing an increasing trend in intensity parameter sequences, an expanding trend in distribution range parameter sequences, and an enrichment trend in morphological detail parameter sequences were selected. These feature elements are directly related to pollution diffusion. Synchronicity analysis was performed on the time series of characteristic elements showing an increasing trend in intensity parameters, an expanding trend in distribution range parameters, and an enrichment trend in morphological detail parameters to determine whether the starting time points of the increasing intensity, expanding distribution range, and enrichment of morphological details of each characteristic element were consistent, and whether the rate of change had a synergistic relationship. Based on the results of the synchronization analysis, feature elements with the same starting time and a synergistic relationship in their rate of change are grouped into the same evolutionary group. Each evolutionary group represents a set of interrelated feature evolution trends during the pollution diffusion process. Evolution rate data of each characteristic element in each evolutionary group were extracted. The evolution rate data included the intensity increase rate, the distribution range expansion rate, and the morphological detail enrichment rate. The magnitude and trend of the intensity increase rate, distribution range expansion rate, and morphological detail enrichment rate in each evolutionary group were analyzed to assess the importance of the evolutionary group in the pollution diffusion process. The evolution groups are sorted according to their overall evolution rate. The evolution groups whose overall evolution rate meets the preset rate standard and whose coverage meets the preset range standard are selected as the main evolution groups. The main evolution groups dominate the overall evolution process of pollution diffusion. Based on the time series, changing trends and comprehensive evolution rate of each characteristic element in the main evolutionary group, a continuous evolution path from the initial state to the current state is constructed, which serves as the positive evolution path for pollution diffusion.

7. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 4, characterized in that, The method of inversely deducing the original distribution range of the optical characteristics of pollutants within the initial area of ​​action based on the direction and rate of expansion includes: The data on the expansion direction of the optical features of pollutants are converted into spatial direction indicators. Each expansion direction corresponds to a spatial direction, representing the spatial extension trend of the optical features of pollutants in that spatial direction. Starting from the current spatial distribution range boundary of the optical characteristics of pollutants, the reverse extension is carried out in the opposite direction of the expansion direction. The range of the reverse extension is estimated based on the historical distribution range change trend in this direction. The starting points of all the reverse extension lines converge to form the initial feature coverage boundary estimated based on the expansion direction and expansion rate. The initial feature coverage boundary estimated based on the expansion direction and expansion rate is the preliminary original distribution range boundary. The spatial distribution range data of the optical characteristics of pollutants at different time points are extracted and arranged in reverse chronological order to obtain the distribution range sequence from the current state to the earlier state; Analyze the magnitude of the change in distribution range between adjacent time nodes in the distribution range sequence, and calculate the difference in distribution range between each time node and the previous time node; Select continuous time node groups whose distribution range difference is less than a preset difference threshold as continuous time node groups with gradual changes in feature distribution range. Extract the spatial distribution range data corresponding to the earliest time node in the continuous time node group where the feature distribution range changes smoothly, and use it as the spatial range data of the early stable distribution of the optical features of the pollutant. The initial feature coverage boundary, which is initially estimated based on the expansion direction and expansion rate, is spatially compared with the spatial distribution range data corresponding to the earliest time node in the stable phase, and the spatial overlap and boundary deviation between the two are calculated. Based on spatial overlap and boundary deviation, the initial feature coverage boundary estimated based on the expansion direction and expansion rate is adjusted so that the adjusted initial feature coverage boundary estimated based on the expansion direction and expansion rate matches the spatial distribution range data corresponding to the earliest time node in the stable phase to the greatest extent. The spatial region enclosed by the initial coverage boundary of the features, which is preliminarily estimated based on the expansion direction and expansion rate, is determined as the original distribution range of the optical features of pollutants within the initial action area. The original distribution range of the optical features of pollutants within the initial action area retains the spatial coverage state when the features have not expanded significantly.

8. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 2, characterized in that, The modal feature channels corresponding to the atmospheric transmission image are subjected to morphological structure analysis to identify image regions in the atmospheric transmission image that are continuously distributed and have a diffusion trend. The contour morphology, distribution density, and extension-related data of these image regions are extracted, including: Image enhancement processing is performed on the modal feature channels corresponding to the atmospheric transmission image. The enhanced image is segmented into regions using a region growing algorithm. Pixels with gray value differences within a preset gray value tolerance range and spatially adjacent are used as seed points to gradually grow into continuous image regions. Morphological analysis was performed on each segmented continuous image region to calculate the morphological parameters of area, perimeter, and circularity of each image region, and image regions with irregular shapes were preliminarily screened out. Extract the boundary contour point coordinates of each filtered image region, and construct the contour morphology model of the image region through a curve fitting algorithm. The contour morphology model describes the external shape and edge undulation features of the image region. The distribution density of pixels in each image region is calculated. The distribution density is obtained by the ratio of the number of pixels in the image region to the area of ​​the image region, which reflects the density of the material distribution in the image region. Analyze the density change gradient of each image region, extract multiple density sampling lines from the center of the image region to the boundary, calculate the density value change at each point on the sampling line, and determine the direction and magnitude of the density gradient. The density gradient direction points in the direction of decreasing density. The expansion trend of an image region is determined based on the direction of the density gradient. If the density gradients in multiple directions show an outward diffusion characteristic, then the image region has a diffusion trend. The extension direction parameters of image regions with diffusion trends are extracted. The extension direction parameters are obtained by the statistical average of each density gradient direction, which represent the overall diffusion direction of the image region. The parameters related to the extension of the image region are calculated. By comparing the area change of the same image region at different time points, the area change rate of the image region per unit time is obtained. By comparing the boundary movement distance of the same image region at different time points, the boundary movement rate of the image region per unit time is obtained. The extension of the image region is described by combining the area change rate and the boundary movement rate. The contour morphology model, distribution density data, density gradient information, extension direction parameters, and extension amplitude data of each image region with a diffusion trend are integrated to form diffusion morphology association features.

9. The pollution source localization method combining UAV remote sensing and AI image recognition according to claim 3, characterized in that, The process of using transition nodes as a reference to continue tracing the state parameters of each feature element in reverse order until the state parameters of all feature elements converge within the same spatial coordinate range and the magnitude of change tends to level off includes: The time point corresponding to the transition node is used as the reference time point for reverse tracing. The state parameters of all feature elements at the reference time point are extracted. The intensity parameter, distribution range parameter, and morphological parameter are the state parameters of all feature elements at the reference time point. Starting from the reference time point, the state parameters of all feature elements of the previous time point are extracted in reverse chronological order to form a parameter sequence for reverse tracing. Compare the changes in the state parameters of each characteristic element between the baseline time point and the previous time point, and calculate the intensity change difference, distribution range offset, and morphological difference. Analyze the changing trends of each characteristic element to determine whether the change is gradually increasing, gradually decreasing, or remaining stable. If the change is gradually decreasing, it indicates that the characteristic element is converging towards its initial state. Extract the distribution range parameters of each feature element in the parameter sequence of reverse tracing, and calculate the intersection area of ​​the distribution range of all feature elements at each time node. The intersection area is the spatial coordinate range covered by all feature elements. Analyze the changes of the intersection region over time in reverse order, calculate the area and boundary coordinates of the intersection region at each time point, and determine whether the intersection region gradually expands, shrinks, or remains stable. If the intersection area gradually shrinks and the boundary tends to be concentrated, then continue to trace the characteristic element state parameters of the previous time node in reverse order, and repeatedly calculate the intersection area and change. When tracing back to a certain time point, if the changes in the state parameters of all feature elements meet the preset stability criteria, and the area of ​​the intersection region no longer shrinks significantly and the boundary remains stable, then the tracing stops, and the state parameters of all feature elements at the corresponding stop time point are extracted to confirm that the intensity parameters, distribution range parameters, and morphological parameters of each feature element are in a stable state, and the intersection region of all distribution range parameters is clear. The intersection area corresponding to the stopping time node is determined as the spatial coordinate range where all characteristic element state parameters converge. The spatial coordinate range where all characteristic element state parameters converge is the central spatial region of the initial stage of pollution diffusion.

10. A pollution source localization system combining UAV remote sensing and AI image recognition, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the pollution source localization method combining UAV remote sensing and AI image recognition as described in any one of claims 1 to 9 by executing the machine-executable instructions.

Images