A water pollution supervision method and system based on remote sensing images

Abstract

The application relates to the technical field of remote sensing image processing and environment monitoring, and discloses a water pollution supervision method and system based on remote sensing images, which comprises the following steps: acquiring multispectral image data, and performing image preprocessing, water body boundary extraction and spectral feature extraction to obtain a water body spectral feature set; performing pollution analysis and processing according to the water body spectral feature set to obtain a pollution classification result; if the pollution classification result shows an abnormal spectral mode, acquiring historical time series data, and combining the historical time series data with a water body exclusive sub-image to perform dynamic trend analysis and determine pollution diffusion dynamic characteristics; performing pollution hotspot area identification according to the pollution diffusion dynamic characteristics to obtain pollution hotspot distribution characteristic data; and performing supervision decision analysis on the pollution hotspot distribution characteristic data to obtain a water pollution supervision result. The method can realize accurate identification, dynamic monitoring and intelligent early warning of lake water pollution, and provides a scientific basis for water resource protection.

Description

A method and system for water pollution monitoring based on remote sensing images Technical Field

[0001] This invention relates to the fields of remote sensing image processing and environmental monitoring technology, and in particular to a method and system for water pollution monitoring based on remote sensing images. Background Technology

[0002] Currently, with the acceleration of industrialization, water pollution is becoming increasingly serious, making real-time, large-scale monitoring of lakes, rivers, and other water bodies a top priority in environmental protection. Satellite remote sensing technology, due to its wide coverage, is gradually becoming an important tool for water environment monitoring. By analyzing the spectral characteristics of remote sensing images, key water quality parameters such as chlorophyll, suspended solids, and dissolved organic matter can be retrieved, thereby enabling an assessment of water pollution levels.

[0003] In existing technologies, remote sensing images are typically processed manually or semi-automatically, using traditional edge detection operators or simple thresholding methods to extract water areas. However, due to the complexity of the aquatic environment, often affected by cloud cover, uneven lighting, and interference from shoreline vegetation, the quality of the original images often falls short of the requirements for high-precision analysis. Although existing image enhancement techniques can improve image contrast and clarity to some extent, traditional methods often struggle to account for the correlation between different bands when dealing with multispectral data, easily leading to distortion of spectral information.

[0004] Therefore, existing technologies suffer from low automation and insufficient recognition accuracy. Summary of the Invention

[0005] This invention provides a water pollution monitoring method and system based on remote sensing images to solve the problems of low automation and insufficient recognition accuracy in existing technologies.

[0006] In a first aspect, to address the aforementioned technical problems, the present invention provides a water pollution monitoring method based on remote sensing images, comprising:

[0007] Acquire multispectral image data and perform image preprocessing on the multispectral image data to obtain a first corrected image;

[0008] Based on the first corrected image, water body boundary extraction processing is performed to determine the distribution boundary of lake water bodies;

[0009] Based on the distribution boundary of the lake water body, a water-specific sub-image is cropped from the first corrected image, and spectral feature extraction processing is performed on the water-specific sub-image to obtain a water body spectral feature set;

[0010] Based on the spectral feature set of the water body, pollution analysis and processing are performed to obtain pollution classification results;

[0011] If the pollution classification result shows an abnormal spectral pattern, then historical time-series data is acquired, and dynamic trend analysis is performed by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion.

[0012] Based on the pollution diffusion dynamic characteristics, the pollution hotspot area identification processing is performed on the water body-specific sub-image to obtain pollution hotspot distribution characteristic data;

[0013] Based on the data on the distribution characteristics of the pollution hotspots, regional statistics and matching analysis were performed, and the water pollution supervision results were obtained through data integration and processing.

[0014] Secondly, the present invention provides a water pollution monitoring system based on remote sensing images, comprising:

[0015] An image preprocessing module is used to acquire multispectral image data and perform image preprocessing on the multispectral image data to obtain a first corrected image;

[0016] The boundary extraction module is used to perform water body boundary extraction processing based on the first corrected image to determine the distribution boundary of the lake water body;

[0017] The feature extraction module is used to crop a water-specific sub-image from the first corrected image according to the distribution boundary of the lake water body, and to perform spectral feature extraction processing on the water-specific sub-image to obtain a water body spectral feature set;

[0018] The pollution analysis module is used to perform pollution analysis processing based on the spectral feature set of the water body to obtain pollution classification results;

[0019] The dynamic trend module is used to acquire historical time-series data if the pollution classification result shows an abnormal spectral pattern, and to perform dynamic trend analysis by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion.

[0020] The hotspot identification module is used to identify pollution hotspot areas in the water body-specific sub-image based on the pollution diffusion dynamic characteristics, and obtain pollution hotspot distribution feature data.

[0021] The decision analysis module is used to perform regional statistics and matching analysis on the pollution hotspot distribution characteristic data, and obtain water pollution supervision results through data integration and processing.

[0022] Compared with the prior art, the present invention has the following beneficial effects:

[0023] (1) This invention effectively eliminates the interference caused by atmospheric scattering, thin clouds and sensor response differences by performing radiometric calibration and atmospheric correction on multispectral image data and extracting water body boundaries by combining morphological filtering, thus ensuring the authenticity of the physical relationship between spectral response feature values ​​and water body substance concentration. At the same time, by detecting invalid pixels and interpolating to repair them, the problem of data loss caused by noise or occlusion is solved, which significantly improves the accuracy and robustness of water pollution inversion in complex environments.

[0024] (2) This invention breaks through the limitations of traditional single-point static monitoring by introducing historical time series data after identifying abnormal spectral patterns and performing spatiotemporal data alignment and trend fitting analysis with the current water body image. This method can accurately calculate the dynamic characteristics of pollution diffusion (including diffusion direction and speed), thereby realizing dynamic prediction of pollution spread trend, providing environmental management departments with time-dimensional early warning information, and avoiding the lag in supervision.

[0025] (3) This invention identifies hotspot areas by performing density clustering on pixels with high diffusion risk and generates structured water pollution supervision results based on priority ranking, thus realizing the automated transformation from "image interpretation" to "data decision-making". This method can automatically filter isolated noise interference, accurately lock the key areas with the most serious pollution and the highest diffusion risk, and directly output executable decision data, which greatly improves the response speed and resource allocation efficiency of water environment supervision. Attached Figure Description

[0026] Figure 1 is a schematic flowchart of a water pollution monitoring method based on remote sensing images provided in the first embodiment of the present invention;

[0027] Figure 2 is a schematic diagram of a water pollution monitoring system based on remote sensing images provided in the second embodiment of the present invention. Detailed Implementation

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

[0029] Referring to Figure 1, the first embodiment of the present invention provides a water pollution monitoring method based on remote sensing images, comprising the following steps:

[0030] S11, acquire multispectral image data, and perform image preprocessing on the multispectral image data to obtain a first corrected image;

[0031] S12, Based on the first corrected image, perform water body boundary extraction processing to determine the distribution boundary of the lake water body;

[0032] S13, Based on the distribution boundary of the lake water body, a water body-specific sub-image is cropped from the first corrected image, and the water body-specific sub-image is subjected to spectral feature extraction processing to obtain a water body spectral feature set;

[0033] S14. Based on the water body spectral feature set, perform pollution analysis and processing to obtain pollution classification results;

[0034] S15, if the pollution classification result shows an abnormal spectral pattern, then acquire historical time series data, and combine the historical time series data with the water body-specific sub-image to perform dynamic trend analysis to determine the dynamic characteristics of pollution diffusion.

[0035] S16, Based on the pollution diffusion dynamic characteristics, the pollution hotspot area identification processing is performed on the water body-specific sub-image to obtain pollution hotspot distribution characteristic data;

[0036] S17. Based on the distribution characteristics data of the pollution hotspots, regional statistics and matching analysis are performed, and the water pollution supervision results are obtained through data integration and processing.

[0037] In step S11, the multispectral image data is preprocessed to obtain a first corrected image, including:

[0038] Acquire multispectral image data, and perform median filtering and noise reduction on the multispectral image data to obtain a first processed image;

[0039] The first processed image is subjected to radiometric calibration and atmospheric correction to obtain an initial corrected image;

[0040] Based on the initial corrected image, a geometric fine correction operation is performed to obtain the first corrected image.

[0041] It should be noted that the acquisition of multispectral image data is achieved through a multispectral sensor mounted on a satellite platform. This sensor can simultaneously acquire radiation signals in the visible light, near-infrared, and short-wave infrared bands. Median filtering is applied to the multispectral image data to eliminate salt-and-pepper noise caused by sensor electronic noise. In this embodiment, a two-dimensional median filtering method is used. This method traverses every pixel in the image through a sliding window, sorts the grayscale values ​​of all pixels within the window, and selects the median as the output value of the center pixel. This operation effectively removes isolated noise points while preserving the edge features of the image.

[0042] It is worth noting that the size of the sliding window (e.g., 3x3) is determined based on statistical analysis of the typical noise speckle sizes in historical image data. Specifically, a large number of original image samples containing noise are collected, and the diameter distribution of noise specks is statistically analyzed using connected component analysis. The specific steps of the connected component analysis method are as follows: first, the noisy image is binarized using a fixed threshold; then, the 8-connectedness algorithm is used to traverse the image to mark all independent pixel connected regions; finally, the equivalent circle diameter of each connected region is calculated. The 95th percentile of the diameter distribution is selected as the reference size. Set the side length of the sliding window to be no less than... The smallest odd integer, this setting ensures that the vast majority of noise points are effectively covered and smoothed.

[0043] It should be noted that performing radiometric calibration and atmospheric correction on the first processed image is a crucial step in restoring the physical authenticity of the image. Firstly, radiometric calibration is performed using the linear transformation formula... The digital observations (DN) of the image are converted into radiance L at the sensor entrance pupil, where and The absolute calibration coefficients are those used when the sensor leaves the factory. The digital observation (DN) of the image is a standard technical term in the field of remote sensing, referring to the raw quantized value obtained by the remote sensing sensor when acquiring data, which is acquired by the multispectral sensor mounted on the satellite platform as described in step S11. Subsequently, atmospheric correction is performed using Dark Object Subtraction (DOS). This method assumes the presence of water bodies or shadowed areas with extremely low reflectivity (i.e., dark objects) in the image, whose received radiance is mainly composed of atmospheric radiation. By subtracting the radiance of the dark object from the radiance of each pixel in the entire image, the influence of atmospheric scattering is eliminated, thereby calculating the true reflectance of the Earth's surface and obtaining the initial corrected image.

[0044] It is worth noting that the determination of the radiance value of the dark target is based on the statistical characteristics of the image histogram. Specifically, the radiance histogram of the image in a specific band (such as the blue light band) is calculated, and the gray value at which the cumulative probability reaches 0.1% is selected as the benchmark value for the dark target. This statistical method avoids the arbitrariness of manually selecting dark targets.

[0045] It should be noted that the geometric fine-tuning operation based on the initial calibrated image aims to eliminate geometric distortions caused by changes in sensor attitude. This embodiment employs a polynomial calibration method based on ground control points (GCPs). The system utilizes a pre-stored high-precision georeferenced map, extracts pairs of corresponding control points using the SIFT algorithm, and constructs a bivariate quadratic polynomial transformation model.

[0046] It is worth noting that the construction process of the pre-stored high-precision geographic reference base map includes two stages: the generation of the digital elevation model (DEM) and orthorectification. The construction method and training process of the digital elevation model are as follows: First, a depth estimation neural network based on a ResNet-50 encoder and a feature pyramid decoder is constructed. The feature pyramid decoder consists of multiple cascaded upsampling modules. Each module uses bilinear interpolation to upsample the feature map by a factor of 2, and then fuses the upsampled feature map with the corresponding spatial resolution feature map in the encoder through channel concatenation. Subsequently, a 3x3 convolutional layer is used for feature smoothing and channel adjustment. Historical stereo satellite image pairs of the target area are collected as network input, and airborne LiDAR point cloud data of the same area are collected and rasterized into high-precision elevation maps as ground truth labels to construct a training dataset (the airborne LiDAR data and stereo satellite imagery can be obtained from public platforms such as USGS or ESA). A weighted combination of pixel-level L1 norm and image gradient loss is used as the loss function, and the Adam optimizer is used to backpropagate and train the network until the loss converges, resulting in a trained depth estimation model. The currently acquired high-resolution optical satellite imagery is input into the trained depth estimation model, and the digital elevation model is generated as the output. Subsequently, the digital elevation model is used to perform pixel-by-pixel differential correction on the optical satellite imagery, and the corrected images are then seamlessly mosaicked and color balanced, finally cropped into a standard-format geographic base map.

[0047] It should be noted that the specific steps for extracting corresponding control point pairs using the SIFT algorithm are as follows: First, a Gaussian difference pyramid is constructed to detect extreme points in scale space, and the key point positions are accurately located by fitting a three-dimensional quadratic function; then, a principal direction is assigned to the key point based on the gradient direction distribution of the neighboring pixels, generating a 128-dimensional feature descriptor; finally, Euclidean distance is used as a similarity metric to match the feature descriptors of the initial corrected image and the high-precision geographic reference base map to obtain corresponding control point pairs.

[0048] In one implementation, the specific steps for constructing a bivariate quadratic polynomial transformation model are as follows: establishing the transformation equation, and ,in For geographic reference base map coordinates, For the initial calibration image coordinates, and These are the coefficients to be determined. Using the extracted pairs of control points with the same name, the least squares method is employed to solve the system of equations, yielding the optimal coefficients. and This allows us to determine the transformation model.

[0049] It should be noted that this model performs resampling and spatial transformation of the image's pixel coordinates, specifically through bilinear interpolation. First, the transformation model is used to calculate the corresponding floating-point coordinates of each pixel in the output image within the input image. Then find the four integer coordinate pixels closest to the floating-point coordinate; finally, according to the distance weight between the floating-point coordinate and these four integer coordinates, perform a weighted average of the gray values ​​of the four pixels to obtain the gray value of the output pixel, so that it is accurately aligned to the standard geographic coordinate system (such as WGS84), and obtain the first corrected image.

[0050] For example, a multispectral image shows significant salt-and-pepper noise in a lake area. The system uses a 3x3 window for median filtering, effectively suppressing the noise. Subsequently, the DN value is converted to radiance based on the gain coefficient in the metadata, and atmospheric correction is performed by subtracting the dark pixel values ​​at 0.1% of the histogram. At this point, the spectral curve of the water body recovers its typical low reflectance characteristics. Finally, geometric correction is performed using 20 ground control points to control the geographic location error of the image to within 0.5 pixels, generating the first corrected image. The geographic location error control to within 0.5 pixels is an absolute error calculated based on the spatial resolution of the image (e.g., 10 meters / pixel), meaning the maximum positioning error does not exceed 5 meters.

[0051] In step S12, based on the first corrected image, water body boundary extraction processing is performed to determine the lake water body distribution boundary, including:

[0052] The first corrected image is subjected to multi-band weighted fusion processing to obtain the second processed image;

[0053] The second processed image is subjected to pixel-level feature classification processing to obtain a third feature image;

[0054] Based on the third feature image, morphological filtering is performed to obtain the fourth distribution image;

[0055] Edge detection and smoothing are performed on the fourth distribution image to determine the distribution boundary of the lake water body.

[0056] It should be noted that multi-band weighted fusion processing is performed on the first corrected image to enhance the contrast between water and non-water areas. In this embodiment, a weighted variant formula of the Normalized Difference Water Index (NDWI) is used for calculation.

[0057]

[0058] in, and These represent the reflectance in the green light band and the near-infrared band, respectively. and , where is the weighting coefficient. Since water has high reflectivity in the green band and strong absorption in the near-infrared band, this fusion can maximize the water features to obtain the second processed image; the above formula focuses on enhancing water contrast and avoids the denominator sensitivity problem of standard NDWI at the mixed pixels.

[0059] It is worth noting that the weighting coefficients and The determination of the target water body is based on the discriminative analysis of historical spectral data. Specifically, historical spectral samples of the water body and typical background objects in the area are collected, and a Fisher discrimination criterion function is constructed. The construction process involves first calculating the mean vectors of the water body sample set and the background object sample set, respectively. and intraclass scatter matrix Then calculate the overall within-class scatter matrix. and inter-class scatter matrix Finally, the objective function is established.

[0060]

[0061] The optimal projection vector is obtained by solving for the maximum value of the function using the Lagrange multiplier method. The normalized components of the vector are used as the preset values ​​for the weight coefficients.

[0062] It should be noted that the pixel-level feature classification processing of the second processed image is implemented using a deep fully convolutional neural network based on the U-Net architecture. The construction process of the deep fully convolutional neural network is as follows: a network structure including a contraction path and an expansion path is built. The contraction path consists of four downsampling modules, each containing two 3x3 convolutional layers (with a stride of 1 and using the ReLU activation function) and a 2x2 max pooling layer (with a stride of 2); the expansion path consists of four upsampling modules, each containing a 2x2 deconvolutional layer, a feature concatenation layer (used to fuse features from the same layer of the contraction path), and two 3x3 convolutional layers; the output layer uses a 1x1 convolutional layer with a Sigmoid activation function to output a pixel-level probability map. The training process involves collecting historical water body remote sensing images and labeling them to generate a binary mask as the training dataset; using a weighted combination of binary cross entropy and Dice coefficients as the loss function; and employing the Adam optimizer with an initial learning rate of 0.001 to perform backpropagation iterative training on the network until the loss function converges on the validation set, resulting in a trained deep fully convolutional neural network. The second processed image is then input into this network, and a binary threshold is applied to the output probability map. Pixels with probabilities higher than the threshold are labeled as foreground, and the rest as background, resulting in the binarized third feature image.

[0063] It is worth noting that the binarization threshold is determined adaptively using the maximum inter-class variance method. This method iterates through all possible probability thresholds, calculates the inter-class variance between the foreground and background, and selects the threshold with the largest variance as the segmentation point for this classification, thereby adapting to image characteristics under different lighting conditions.

[0064] It should be noted that morphological filtering based on the third feature image aims to remove isolated noise points in the classification results and fill voids within the water body. This embodiment employs a combination of morphological closing and opening operations. First, a closing operation is performed, specifically by using a preset structuring element to dilate the image, i.e., taking the maximum pixel value within the structuring element's coverage area to connect adjacent broken areas, followed by erosion, i.e., taking the minimum pixel value within the structuring element's coverage area to restore boundary positions and fill small voids. Subsequently, an opening operation is performed, specifically by using the same structuring element to erode the image to eliminate small isolated noise spots, followed by dilation to restore the shape of the main regions, resulting in the fourth distribution image with better spatial connectivity.

[0065] It is worth noting that the determination of the structuring element size in the morphological operation is based on statistical analysis of the sizes of typical noise artifacts in the classification images. Specifically, the average equivalent diameter of connected components marked as "noise" in historical classification results is statistically analyzed using a connected component labeling method. And set the size of the structuring element to be no less than The smallest odd integer is used to ensure that noise can be effectively filtered out without disrupting the actual water body boundary.

[0066] It should be noted that the edge detection and smoothing process for the fourth distribution image involves first extracting edge pixels from the binary image using the Canny operator, and then applying a B-spline curve interpolation algorithm to fit the extracted discrete edge points. This operation eliminates the jagged effect caused by pixel meshing, generating a continuous and smooth closed curve, which represents the boundary of the lake's water distribution.

[0067] For example, in the weighted fused image, the water area exhibits a bright feature. The trained deep fully convolutional neural network classifies it as a binary mask, but the mask contains several noise points with a diameter of approximately 3 pixels. After performing morphological opening operations using 5x5 structuring elements, the noise points are completely removed. Finally, the Canny operator extracts the edges, and a smooth, jagged lake boundary line is obtained through B-spline fitting.

[0068] In step S13, the water body-specific sub-image is subjected to spectral feature extraction processing to obtain a water body spectral feature set, including:

[0069] Multi-band reflectance data were extracted from the water body-specific sub-image to obtain a band reflectance dataset;

[0070] For the aforementioned band reflection dataset, invalid pixel detection and interpolation repair are performed to obtain a reflection feature set;

[0071] The reflection feature set is normalized to obtain the water body spectral feature set.

[0072] It should be noted that the multi-band reflectance data extraction of the water-specific sub-image refers to using the boundary determined in the previous step as a mask to crop out a pixel matrix containing only the water area from the first corrected image, and reading the reflectance values ​​of each pixel in the matrix in all bands to form the band reflectance dataset.

[0073] It should be noted that the invalid pixel detection and interpolation repair processing on the aforementioned band reflectance dataset aims to remove and repair abnormal data caused by cloud cover, sensor defects, and water glare. The detection process is implemented using a range threshold discrimination method. If the reflectance value of a pixel in any band exceeds a preset effective physical range or deviates from the local mean by more than a preset statistical threshold, it is marked as an "invalid pixel." Subsequently, an inverse distance weighted interpolation algorithm is used to weight and replace the invalid pixel using the effective pixel values ​​within a preset neighborhood. The calculation formula for the inverse distance weighted interpolation algorithm is as follows:

[0074]

[0075] Where Z is the interpolation result. For the neighboring region Reflectance of each effective pixel For weight (when When the interpolation point coincides with a valid pixel, the value of that valid pixel is directly taken as the interpolation result to avoid division by zero errors. ), For invalid pixels and the first Euclidean distance of effective pixels, It is the power exponent (usually taken as 2). This represents the number of valid pixels in the neighborhood. This process yields the complete set of reflection features.

[0076] It is worth noting that the determination of the preset effective physical range is based on a joint analysis of the theoretical boundary of the radiative transfer physical model and the dynamic range of the sensor. Specifically, the theoretical physical domain of surface reflectance (i.e., 0 to 1) is used as the basic interval, and corrected by combining the statistical values ​​of the sensor's saturation response under historical extreme lighting conditions to cover possible sensor overexposure. The determination of the statistical threshold is based on the 3-Sigma criterion. Specifically, the mean reflectance of all pixels within the current water body area is calculated. and standard deviation , will the interval Values ​​outside of these ranges are defined as outliers. The determination of the neighborhood radius in the inverse distance weighted interpolation algorithm is based on spatial autocorrelation analysis of the water spectrum. The range at which the semivariogram reaches the sill value is selected as the reference value for the neighborhood radius.

[0077] It should be noted that data normalization of the reflection feature set is performed to eliminate differences in dimensions and numerical ranges across different bands, facilitating the convergence of subsequent algorithms. This embodiment employs a minimum-maximum normalization method, linearly mapping the reflectance data of each band to the [0,1] interval. The calculation formula is as follows:

[0078]

[0079] The processed dataset is the water body spectral feature set.

[0080] For example, in the extracted band reflectance dataset, a pixel was detected with a reflectance of 0.25 in the blue light band, a value significantly higher than that of normal water. The system determined this pixel to be invalid based on the 3-Sigma criterion. Subsequently, using the reflectance of valid pixels within a 5x5 radius around the pixel, an estimated value of 0.08 was calculated using an inverse distance-weighted interpolation algorithm and replaced with this value. Finally, all data was normalized to generate a standardized feature set.

[0081] In step S14, pollution analysis is performed based on the water body spectral feature set to obtain pollution classification results, including:

[0082] The spectral feature set of the water body is subjected to band feature decomposition processing to obtain a band feature dataset;

[0083] Spectral response feature values ​​are extracted from the band feature dataset. If the spectral response feature values ​​meet the preset anomaly discrimination conditions, anomaly marking processing is performed to obtain an anomaly marking dataset.

[0084] Based on the anomaly-marked dataset, pollution concentration levels are classified to obtain a preliminary pollution classification list;

[0085] The preliminary pollution classification list is subjected to spatial continuity verification processing to obtain the pollution classification results.

[0086] It should be noted that the band feature decomposition processing of the water body spectral feature set is achieved using principal component analysis (PCA). This process first involves removing the mean center from the original multi-band spectral data matrix X, and then calculating its covariance matrix C using the following formula:

[0087]

[0088] in Let C be the sample size. Next, eigenvalue decomposition is performed on the covariance matrix C, i.e., the characteristic equation is solved. Obtain eigenvalues The eigenvalues ​​and their corresponding eigenvectors are arranged in descending order, and the eigenvectors corresponding to the first two eigenvalues ​​are selected to form the projection matrix V. Finally, the linear transformation formula is used... The original high-dimensional spectral data is projected onto an orthogonal low-dimensional feature space, and the first principal component PC1 and the second principal component PC2 are extracted. These principal components constitute the band feature dataset. This operation can effectively remove redundant information between bands and concentrate contamination features on a few components.

[0089] It should be noted that the spectral response feature values ​​extracted from the band feature dataset are obtained by selecting the value of the first principal component (PC1), which has the strongest correlation with water turbidity, as the metric. If the spectral response feature value meets the preset anomaly discrimination condition, i.e., the value is greater than the preset pollution threshold, it indicates that the spectral features of the pixel deviate from the baseline pattern of clean water, and the system marks the corresponding pixel as "potentially polluted," thus obtaining the anomaly-marked dataset.

[0090] It is worth noting that the predetermined pollution threshold is determined based on regression analysis of historical water quality monitoring data and geostationary satellite data. The historical water quality monitoring data covers at least three years of data from both dry and wet seasons, with a sample size of no less than 100 sets to ensure the representativeness of the regression model. Specifically, historical measured suspended solids concentration data for the water area are collected as the dependent variable, and the corresponding remote sensing image PC1 value is used as the independent variable to construct a least-squares linear regression model. Based on the pollutant concentration limits specified in the national water environment quality standards, the corresponding PC1 characteristic threshold is derived through this regression model and used as a preset discrimination criterion.

[0091] It should be noted that the pollution concentration level classification based on the aforementioned anomaly-marked dataset employs a multi-level threshold segmentation method. The system maps the spectral response feature value of each anomaly pixel to a "mild," "moderate," or "severe" pollution level according to a preset classification standard. The specific logic of this mapping is as follows: if the feature value... satisfy If so, it is mapped to "light pollution"; if the following conditions are met... If so, it is mapped to "moderate pollution"; if If so, it is mapped to "severe pollution". This mapping process generates the preliminary pollution classification list.

[0092] It is worth noting that the preset grading standard (i.e., threshold) The determination of the threshold is based on cluster analysis of historical abnormal spectral data. Specifically, the PC1 values ​​of pixels that have been historically identified as contaminated are collected, and the K-means clustering algorithm is used to divide them into three clusters. The midpoint of the cluster center value of these three clusters is selected as the segmentation threshold for each level.

[0093] It should be noted that the spatial continuity verification process performed on the preliminary pollution classification list aims to eliminate random salt-and-pepper noise in the classification results and ensure the spatial coherence of the polluted areas. This embodiment employs a sliding window majority voting method. The system defines a window of a preset size and traverses each pixel in the classification list; it counts the number of pixels for each pollution level within the window. If the category of the central pixel does not match the category with the highest frequency (mode) within the window, the category of the central pixel is corrected to the mode category. This process effectively removes isolated misclassification points, resulting in the final pollution classification result.

[0094] It is worth noting that the size of the sliding window (e.g., 5x5) is determined based on spatial autocorrelation analysis of water pollution distribution. Specifically, the semi-variogram of historical pollution images is calculated, the range at which the variogram reaches the sill value is selected as the characteristic scale of spatial correlation, and the window side length is set to the largest odd integer not exceeding this range to ensure that the pixels within the window have statistical spatial dependence.

[0095] For example, after performing PCA on the feature set, the system extracts the first principal component (PC1) value of a certain pixel, which is 0.45. Based on the threshold (0.3) determined by historical regression analysis, this value meets the anomaly detection criteria. The system initially labels it as "moderately contaminated". In the spatial continuity verification, it is found that 20 pixels in the 5x5 area around this pixel are also labeled as "moderately contaminated", which satisfies the majority principle, thus confirming the validity of the classification result.

[0096] In step S15, if the pollution classification result shows an abnormal spectral pattern, historical time-series data is acquired, and dynamic trend analysis is performed by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion, including:

[0097] If the pollution classification result shows an abnormal spectral pattern, then historical time-series data is acquired, and the water body-specific sub-image and the historical time-series data are spatiotemporally aligned and integrated to obtain a time-series regional monitoring dataset.

[0098] Trend fitting analysis was performed on the time-series regional monitoring dataset to obtain the dynamic change pattern of pollution.

[0099] Based on the pollution dynamic change pattern, boundary distance changes are calculated to determine the pollution diffusion dynamic characteristics.

[0100] It should be noted that acquiring historical time-series data refers to retrieving archived spectral feature images of the water body from a pre-built historical monitoring database at several past time points (e.g., four times within the past 24 hours). The historical monitoring database is constructed by periodically acquiring satellite images of the target water body, performing image preprocessing in step S11 and spectral feature extraction in step S13, and storing the generated feature images along with their acquisition timestamps and geographic coordinate information in the time-series database. The spatiotemporal data alignment and integration of the water body-specific sub-images and the historical time-series data is achieved using a bilinear interpolation resampling method based on geographic coordinates. The specific mapping process of this method is as follows: first, based on the geographic coordinate projection relationship, the integer coordinate pixels in the current water body-specific sub-image are calculated. The corresponding floating-point coordinates in the historical image coordinate system Then, identify the four nearest integer neighboring pixels to the floating-point coordinate; finally, based on the distance weights between the floating-point coordinate and these four neighboring pixels, perform a weighted average of the spectral values ​​of the four pixels, and assign the result to the current coordinate. This allows for the construction of a spatially aligned three-dimensional array (X,Y,Time), which is the time-seriesd regional monitoring dataset.

[0101] It should be noted that the trend fitting analysis of the time-series regional monitoring dataset involves extracting the spectral response feature value sequence of each spatial pixel location over time and performing linear regression analysis using the least squares method. Specifically, a linear equation is established. , where y is the spectral eigenvalue, t is time, k is the slope, and b is the intercept. This is achieved by minimizing the sum of squared errors. Solve for the optimal slope k. In the formula, Represents the total number of historical points in time. Representing the Spectral characteristic observations at each time point Representing the The timestamps at each time point are given, and k and b are the slope and intercept parameters to be solved, respectively. The slope k represents the rate of change of the pollution level at that location over time. If k is positive and its value exceeds a preset trend threshold, it is determined that the pollution at that location is intensifying. The slope changes of all pixels constitute the dynamic change pattern of the pollution.

[0102] It is worth noting that the preset trend threshold is determined based on statistical analysis of historical spectral fluctuations in normal water bodies. Specifically, historical spectral data of the water body during periods without pollution events are collected, the slope distribution of its change over time is calculated, and the 99th percentile of this distribution is selected as the critical value for determining pollution aggravation.

[0103] It should be noted that the boundary distance change calculation based on the pollution dynamic change pattern aims to quantify the diffusion speed of the polluted area. The system first binarizes the polluted area at the current and previous times based on the pollution dynamic change pattern and historical classification results, and uses the Canny operator to extract the edge pixel sets of both, defining them as the current boundary line and the previous boundary line, respectively. Then, a chamfered distance transformation algorithm is used to calculate the shortest Euclidean distance from each point on the current boundary to the boundary contour at the previous time. The specific calculation process of the chamfered distance transformation is as follows: A distance matrix of the same size as the image is initialized, setting the pixel value corresponding to the boundary line at the previous time to 0, and setting the rest to infinity; a preset distance mask (such as a 3-4 mask) is used to scan the matrix twice. The first scan is from the upper left corner to the lower right corner, updating the distance value of each pixel to the minimum distance between it and the pixels in the mask-covered neighborhood; the second scan is from the lower right corner to the upper left corner and updated again, finally obtaining a distance map containing the shortest distance information for the entire image; the value at the coordinates of the current boundary line is read from this distance map, which is the shortest Euclidean distance. Finally, by combining the time difference between the two moments, the magnitude of the diffusion velocity is calculated, and the normal direction from the boundary of the previous moment to the current boundary is taken as the diffusion direction. This velocity magnitude and diffusion direction together constitute the dynamic characteristics of the pollution diffusion; wherein, the time difference is in hours, and the diffusion velocity is in kilometers per hour.

[0104] It is worth noting that the determination of the preset distance mask (e.g., 3-4 mask) is based on the approximation error analysis of the integer distance transformation method. Specifically, the maximum relative error between the distance field generated by different mask templates (e.g., 3-4 template, 5-7-11 template) and the standard Euclidean distance field is calculated, and the template that can control the maximum relative error within the tolerance range (e.g., 6%) and has the lowest computational cost is selected as the preset mask.

[0105] For example, the system detects contamination in the current image. After retrieving and aligning historical data from the past three times, it finds that the spectral index of the center point of a certain area shows a linear upward trend over time, with a calculated regression slope k of +0.05 / hour, exceeding the preset trend threshold of 0.01 / hour. Comparing the contamination boundary between the current and previous moments (1 hour apart), the shortest distance eastward from the boundary is calculated to be 0.5 kilometers. Based on this, it is determined that the contamination in this area is spreading eastward at a speed of 0.5 kilometers per hour.

[0106] In step S16, based on the pollution diffusion dynamic characteristics, the water body-specific sub-image is processed to identify pollution hotspot areas, resulting in pollution hotspot distribution characteristic data, including:

[0107] Based on the pollution diffusion dynamic characteristics, high diffusion risk pixels are extracted to obtain the target pixel set;

[0108] Density clustering is performed on the target pixel set to determine the center position and coverage area;

[0109] Based on the central location and coverage area, vectorized data encapsulation is performed to obtain the pollution hotspot distribution characteristic data.

[0110] It should be noted that the extraction of high-risk pixels is achieved through multi-condition threshold filtering. The system filters out pixels that simultaneously meet the conditions of "high current pollution concentration" and "fast diffusion speed". Specifically, it iterates through all pixels in the water-specific sub-image and retains only those pixels that were classified as "severely polluted" in step S14 and whose diffusion speed calculated in step S15 exceeds a preset speed threshold, forming the target pixel set.

[0111] It is worth noting that the velocity threshold was determined based on statistical analysis of historical hydrodynamic data for the water area. Specifically, flow velocity data from historically significant pollution diffusion events in the water area were collected, the 80th percentile of the flow velocity distribution was calculated, and this quantile value was used as the critical velocity threshold for determining high diffusion risk.

[0112] It should be noted that the density clustering processing of the target pixel set is implemented using the DBSCAN algorithm. This algorithm does not require pre-specifying the number of clusters and can automatically discover clusters of arbitrary shapes and remove noise points. The algorithm mainly relies on two parameters: neighborhood radius (Eps) and minimum number of points (MinPts). The system takes the spatial coordinates of each pixel as input and clusters density-connected pixels into an independent hotspot region. For each identified cluster, the arithmetic mean of the coordinates of all pixels within the cluster is calculated as the center position; and the Graham Scan algorithm is used to calculate the minimum convex polygon surrounding all pixels within the cluster, which is used as the coverage area. The specific calculation process of the Graham Scan algorithm is as follows: first, the point with the smallest ordinate (or the smallest abscissa if the ordinates are the same) within the cluster is selected as the base point. Then, the remaining points are arranged according to... Sort the polar angles from smallest to largest; then initialize a stack and... The first two sorted points are pushed onto the stack. Then, the remaining points are traversed, and the turning relationship between the current point and the top two points of the stack is determined (by calculating the cross product). If it is not a left turn (the cross product is less than or equal to 0), the top element of the stack is popped. This process continues until a left turn relationship is satisfied or only two points remain in the stack. The current point is then pushed onto the stack. Finally, the sequence of points retained in the stack is the vertex sequence of the convex polygon.

[0113] It is worth noting that the parameters Eps and MinPts of the DBSCAN algorithm are determined based on K-distance graph analysis. Specifically, the distances from all data points to their Kth nearest neighbors are calculated, and the sorted distance curves are plotted. The distance corresponding to the "elbow" point of the curve is selected as the preset value of Eps, and MinPts is set to 4 to ensure the stability of the cluster. The elbow point is determined by calculating the maximum value of the second derivative of the distance between adjacent points, corresponding to the position with the greatest curvature of the curve.

[0114] It should be noted that vectorized data encapsulation involves constructing the calculated center coordinates, coverage area, and average pollution index of the region into a data object conforming to the GeoJSON standard format. Specifically, the construction process involves creating a root object of type "Feature"; within the root object, creating a "geometry" sub-object, setting its "type" attribute to "Polygon" and its "coordinates" attribute to a closed array composed of vertex latitude and longitude sequences output by the Graham scan algorithm; and within the root object, creating a "properties" sub-object, assigning the center latitude and longitude to the "center_lon" and "center_lat" fields respectively, and assigning the average pollution index of the region to the "pollution_index" field. This encapsulation process yields the pollution hotspot distribution characteristic data, facilitating subsequent reading and display by the system.

[0115] For example, the system filters out high-risk pixels with a diffusion speed greater than 0.3 km / h and belonging to the category of heavily polluted areas. The DBSCAN algorithm clusters these discrete points into two independent hotspot regions. The centroid of region A is calculated to be (Lon: 113.5, Lat: 22.8), and the coverage area obtained by the Graham scan method is a convex polygon with an area of ​​1.2 square kilometers. The system encapsulates this information into a GeoJSON object, identifying the region as a "high-risk hotspot".

[0116] In step S17, regional statistics and matching analysis are performed on the pollution hotspot distribution characteristic data, and water pollution supervision results are obtained after data integration and processing, including:

[0117] Regional statistical analysis was performed on the pollution hotspot distribution characteristic data to obtain pollution concentration distribution statistical information;

[0118] The pollution concentration distribution statistics are matched and analyzed with preset alarm triggering conditions to obtain a list of key areas;

[0119] The list of key areas is sorted by priority to obtain a sorted sequence of areas;

[0120] Based on the sorted regional sequence, data integration and format conversion are performed to obtain the water pollution monitoring results.

[0121] It should be noted that the regional statistical analysis of the pollution hotspot distribution data involves traversing all pixels within the coverage area of ​​each independent hotspot region (polygon). The arithmetic mean of the pollution indices of all pixels within the region is calculated as the average pollution concentration; the maximum value is selected as the maximum concentration; the total number of pixels contained in the region is counted and multiplied by the geographical area of ​​a single pixel (determined by spatial resolution, e.g., 100 square meters) to obtain the total area of ​​the region. These indicators are then summarized to generate the pollution concentration distribution statistics.

[0122] It should be noted that the matching analysis between the pollution concentration distribution statistics and the preset alarm triggering conditions is implemented using conditional logic judgment. The preset alarm triggering conditions define a combination of thresholds at different levels. Specifically, if the average concentration is greater than the first-level concentration threshold and the total area is greater than the first-level area threshold, it is marked as "Level 1 Alarm (Red)"; if the average concentration is greater than the second-level concentration threshold but less than the first-level threshold, it is marked as "Level 2 Alarm (Orange)". The system traverses all hotspot areas, adds areas that meet any alarm condition to the key area list, and records their alarm level.

[0123] It is worth noting that the threshold settings for the alarm triggering conditions are based on a digital mapping of the limits for various pollutants in the National Surface Water Environmental Quality Standard. Specifically, the Class V water quality limit specified in the standard is set as the first-level concentration threshold, and the Class IV water quality limit is set as the second-level concentration threshold. The area threshold is determined based on the statistical analysis of emergency response records for historical sudden water pollution events, selecting the smallest polluted area at the time of activation of the first-level response in history as the first-level area threshold.

[0124] It should be noted that the priority ranking of the key areas is intended to determine the order of remediation. This embodiment uses a multi-attribute comprehensive scoring method to calculate the risk score for each area. The scoring formula is as follows: ,in These represent the average concentration, total area, and diffusion rate after max-min normalization, respectively. These are weighting coefficients. The system calculates them... The regions are sorted in descending order to obtain the sorted region sequence.

[0125] It is worth noting that the weighting coefficients The determination is based on the objective weighting of historical pollution data using the entropy weight method. Specifically, this involves collecting... Concentration, area, and velocity data from historical pollution events were used to construct... The original evaluation matrix; the matrix is ​​normalized to eliminate the influence of dimensions; the first evaluation matrix is ​​calculated. The first item under the indicator The proportion of characteristics of each event Using the formula,

[0126]

[0127] Calculate the first Information entropy of the indicator Finally, according to the formula,

[0128]

[0129] Calculate the weight of each indicator. The greater the degree of variation of an indicator, the more information it provides, and the greater the weight it is assigned, thus avoiding the subjectivity of manual weighting.

[0130] It should be noted that the historical pollution events include typical events such as industrial discharge and algal blooms, with a sample size of no less than 50, covering different seasons and hydrological conditions. The data integration and format conversion process involves constructing the sorted regional sequence into a structured data object conforming to the JSON standard. Specifically, the construction process is as follows: create a root object containing an "AlertList" array; traverse the sorted regional sequence, generating a sub-object for each region, which contains "Rank" (sorting position), "Location" (center latitude and longitude), "Level" (alarm level), "Metrics" (containing concentration, area, and velocity values), and "Action" (suggested measures matched according to the level, such as 'interception' or 'continuous monitoring'); fill the sub-objects into the array in order to obtain the water pollution monitoring results.

[0131] For example, the system identifies two key regions. Region A (normalized concentration 0.8, area 0.9, velocity 0.5) is calculated using the entropy weight method (weights set to...). The risk score for region A was 0.78; the risk score for region B (normalized concentration 0.6, area 0.3, velocity 0.2) was 0.40. The system ranked region A first and generated a JSON data packet. , It is sent to the regulatory platform.

[0132] In summary, this invention constructs a closed-loop monitoring system that integrates high-precision preprocessing of multispectral images and water body boundary extraction, quantitative pollution analysis based on spectral response feature values, dynamic trend prediction and intelligent identification of hotspot areas by fusing historical time-series data. This system deeply integrates the physical mechanism of remote sensing inversion with spatiotemporal data mining technology, solving the technical problems of static monitoring lag, unclear pollution boundary identification, and lack of dynamic early warning capabilities in existing technologies. It achieves all-weather, full-process precise monitoring and scientific decision support for lake water pollution.

[0133] Referring to Figure 2, a second embodiment of the present invention provides a water pollution monitoring system based on remote sensing images, comprising:

[0134] An image preprocessing module is used to acquire multispectral image data and perform image preprocessing on the multispectral image data to obtain a first corrected image;

[0135] The boundary extraction module is used to perform water body boundary extraction processing based on the first corrected image to determine the distribution boundary of the lake water body;

[0136] The feature extraction module is used to crop a water-specific sub-image from the first corrected image according to the distribution boundary of the lake water body, and to perform spectral feature extraction processing on the water-specific sub-image to obtain a water body spectral feature set;

[0137] The pollution analysis module is used to perform pollution analysis processing based on the spectral feature set of the water body to obtain pollution classification results;

[0138] The dynamic trend module is used to acquire historical time-series data if the pollution classification result shows an abnormal spectral pattern, and to perform dynamic trend analysis by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion.

[0139] The hotspot identification module is used to identify pollution hotspot areas in the water body-specific sub-image based on the pollution diffusion dynamic characteristics, and obtain pollution hotspot distribution feature data.

[0140] The decision analysis module is used to perform regional statistics and matching analysis on the pollution hotspot distribution characteristic data, and obtain water pollution supervision results through data integration and processing.

[0141] It should be noted that the water pollution monitoring system based on remote sensing images provided in this embodiment of the invention is used to execute all the process steps of the water pollution monitoring method based on remote sensing images in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0142] This invention also provides an electronic device. The electronic device includes a processor, a memory, and a computer program stored in the memory and executable on the processor, such as a water pollution monitoring program based on remote sensing images. When the processor executes the computer program, it implements the steps in the various embodiments of the water pollution monitoring method based on remote sensing images described above, such as step S11 shown in FIG1. ​​Alternatively, when the processor executes the computer program, it implements the functions of each module / unit in the various system embodiments described above, such as an image preprocessing module.

[0143] For example, the computer program may be divided into one or more modules / units, which are stored in the memory and executed by the processor to complete the present invention. The one or more modules / units may be a series of computer program instruction segments capable of performing a specific function, which describe the execution process of the computer program in the electronic device.

[0144] The electronic device may be a desktop computer, laptop, handheld computer, or smart tablet, etc. The electronic device may include, but is not limited to, a processor and memory. Those skilled in the art will understand that the above components are merely examples of electronic devices and do not constitute a limitation on the electronic device. It may include more or fewer components than described above, or combine certain components, or different components. For example, the electronic device may also include input / output devices, network access devices, buses, etc.

[0145] The processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor can be a microprocessor or any conventional processor. The processor is the control center of the electronic device, connecting all parts of the electronic device via various interfaces and lines.

[0146] The memory can be used to store the computer programs and / or modules. The processor implements various functions of the electronic device by running or executing the computer programs and / or modules stored in the memory and by calling data stored in the memory. The memory may mainly include a program storage area and a data storage area. The program storage area may store the operating system, at least one application program required for a function (such as sound playback function, image playback function, etc.), etc.; the data storage area may store data created according to the use of the mobile phone (such as audio data, phonebook, etc.). In addition, the memory may include high-speed random access memory, and may also include non-volatile memory, such as hard disk, memory, plug-in hard disk, smart media card (SMC), secure digital (SD) card, flash card, at least one disk storage device, flash memory device, or other volatile solid-state storage device.

[0147] If the modules / units integrated into the electronic device are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the methods of the above embodiments of the present invention can also be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program includes computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. The computer-readable medium can include: any entity or system capable of carrying the computer program code, recording media, USB flash drives, portable hard drives, magnetic disks, optical disks, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the content included in the computer-readable medium can be appropriately added or removed according to the requirements of legislation and patent practice in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practice, computer-readable media do not include electrical carrier signals and telecommunication signals.

[0148] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the system embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0149] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for monitoring water pollution based on remote sensing images, characterized in that, include: Multispectral image data is acquired and preprocessed to obtain a first corrected image. Based on the first corrected image, water body boundary extraction is performed to determine the lake's water body distribution boundary. Based on the lake's water body distribution boundary, a water-specific sub-image is cropped from the first corrected image, and spectral feature extraction is performed on the water-specific sub-image to obtain a water body spectral feature set. Based on the water body spectral feature set, pollution analysis is performed to obtain pollution classification results. If the pollution classification result shows an abnormal spectral pattern, then historical time-series data is acquired, and dynamic trend analysis is performed by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion. Based on the pollution diffusion dynamic characteristics, pollution hotspot area identification processing is performed on the water body-specific sub-image to obtain pollution hotspot distribution feature data; regional statistics and matching analysis are performed on the pollution hotspot distribution feature data, and water pollution monitoring results are obtained after data integration processing; wherein, the step of extracting water body boundaries based on the first corrected image to determine the lake water body distribution boundary includes: performing multi-band weighted fusion processing on the first corrected image to obtain a second processed image; performing pixel-level feature classification processing on the second processed image to obtain a third feature image; performing morphological filtering processing on the third feature image to obtain a fourth distribution image; performing edge detection and smoothing processing on the fourth distribution image to determine the lake water body distribution boundary; wherein, performing multi-band weighted fusion processing on the first corrected image to obtain a second processed image includes: using a weighted variant formula of the normalized differential water body index. Perform calculations; where, and These represent the reflectance in the green light band and the near-infrared band, respectively, with weighting coefficients. and The determination of the optimal projection vector is based on the discrimination analysis of historical spectral data of the target water area. By collecting historical spectral samples of water bodies and typical background objects in the area, the Fisher discrimination criterion function is constructed to solve for the optimal projection vector, and its normalized component is used as the weight coefficient.

2. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, The step of preprocessing the multispectral image data to obtain a first corrected image includes: acquiring multispectral image data and performing median filtering denoising on the multispectral image data to obtain a first processed image; performing radiometric calibration and atmospheric correction on the first processed image to obtain an initial corrected image; and performing geometric fine correction based on the initial corrected image to obtain the first corrected image.

3. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, The step of extracting spectral features from the water body-specific sub-image to obtain a water body spectral feature set includes: extracting multi-band reflectance data from the water body-specific sub-image to obtain a band reflectance dataset; performing invalid pixel detection and interpolation repair processing on the band reflectance dataset to obtain a reflectance feature set; and performing data normalization processing on the reflectance feature set to obtain the water body spectral feature set.

4. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, The step of performing pollution analysis processing based on the water body spectral feature set to obtain pollution classification results includes: performing band feature decomposition processing on the water body spectral feature set to obtain a band feature dataset; extracting spectral response feature values ​​from the band feature dataset; if the spectral response feature values ​​meet preset anomaly discrimination conditions, performing anomaly labeling processing to obtain anomaly label dataset; performing pollution concentration level classification processing based on the anomaly label dataset to obtain a preliminary pollution classification list; and performing spatial continuity verification processing on the preliminary pollution classification list to obtain the pollution classification results.

5. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, If the pollution classification result shows an abnormal spectral pattern, then historical time-series data is acquired, and dynamic trend analysis is performed on the historical time-series data and the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion. This includes: if the pollution classification result shows an abnormal spectral pattern, then historical time-series data is acquired, and the water body-specific sub-image and the historical time-series data are spatiotemporally aligned and integrated to obtain a time-series regional monitoring dataset; trend fitting analysis is performed on the time-series regional monitoring dataset to obtain the pollution dynamic change pattern; and boundary distance changes are calculated based on the pollution dynamic change pattern to determine the dynamic characteristics of pollution diffusion.

6. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, The step of identifying pollution hotspot areas in the water body-specific sub-image based on the pollution diffusion dynamic characteristics to obtain pollution hotspot distribution feature data includes: extracting high-diffusion-risk pixels based on the pollution diffusion dynamic characteristics to obtain a target pixel set; performing density clustering on the target pixel set to determine the center location and coverage area; and performing vectorized data encapsulation based on the center location and coverage area to obtain the pollution hotspot distribution feature data.

7. The water pollution monitoring method based on remote sensing images according to claim 1, characterized in that, The process of performing regional statistical analysis and matching on the pollution hotspot distribution characteristic data, and then integrating the data to obtain water pollution supervision results, includes: performing regional statistical analysis on the pollution hotspot distribution characteristic data to obtain pollution concentration distribution statistics; matching the pollution concentration distribution statistics with preset alarm triggering conditions to obtain a list of key areas; prioritizing the list of key areas to obtain a sorted area sequence; and integrating and formatting the data according to the sorted area sequence to obtain the water pollution supervision results.

8. A water pollution monitoring system based on remote sensing images, used to implement the water pollution monitoring method based on remote sensing images as described in any one of claims 1-7, characterized in that, include: An image preprocessing module is used to acquire multispectral image data and perform image preprocessing on the multispectral image data to obtain a first corrected image; a boundary extraction module is used to perform water body boundary extraction processing based on the first corrected image to determine the distribution boundary of the lake water body; a feature extraction module is used to crop water body-specific sub-images from the first corrected image based on the lake water body distribution boundary, and perform spectral feature extraction processing on the water body-specific sub-images to obtain a water body spectral feature set; a pollution analysis module is used to perform pollution analysis processing based on the water body spectral feature set to obtain pollution classification results. The dynamic trend module is used to acquire historical time-series data if the pollution classification result shows an abnormal spectral pattern, and to perform dynamic trend analysis by combining the historical time-series data with the water body-specific sub-image to determine the dynamic characteristics of pollution diffusion. The hotspot identification module is used to identify pollution hotspot areas in the water body-specific sub-image based on the pollution diffusion dynamic characteristics, and obtain pollution hotspot distribution characteristic data; the decision analysis module is used to perform regional statistics and matching analysis on the pollution hotspot distribution characteristic data, and obtain water pollution supervision results through data integration processing.

Images