Data analysis method for monitoring soil pollution spatial distribution based on image fusion

Abstract

The invention discloses a data analysis method for monitoring soil pollution spatial distribution based on image fusion. The method comprises the following steps: step 1, collecting soil pollution monitoring data; 2, preprocessing the soil pollution monitoring data; step 3, performing image fusion calibration; step 4, extracting soil pollutant characteristics; 5, generating a soil health index; and step 6, monitoring the soil environment. Compared with the prior art, the soil pollution monitoring method has the advantages that all-directional monitoring of soil pollution is achieved through the image fusion technology, the limitation of a traditional monitoring method is broken through, the spatial distribution characteristics of pollutants can be comprehensively reflected, and the soil pollution monitoring range and accuracy are effectively improved.

Description

Technical Field

[0001] This invention is a data analysis method for monitoring the spatial distribution of soil pollution based on image fusion, belonging to the field of electronic data processing technology. Background Technology

[0002] Soil pollution is a significant environmental issue that has serious impacts on ecosystems and human health. To effectively monitor and assess soil pollution, scientists have conducted extensive research and developed various methods for analyzing sampling and monitoring data.

[0003] The following technical methods are mainly used for the detection of soil pollution in existing technologies: I. Laboratory Chemical Analysis Techniques: After collecting soil samples, pollutant components and contents are detected through chemical reagent reactions and instrumental analysis (such as chromatography and mass spectrometry). Soil pollution monitoring relies on manual sampling and laboratory analysis to determine the concentration of heavy metals or organic pollutants at specific sampling points. This method has high accuracy in the quantitative detection of pollutants, but its monitoring range is limited and cannot fully reflect the spatial distribution characteristics of pollutants.

[0004] II. In-situ direct detection technology: The sensor directly contacts or detects the soil at close range to obtain pollutant-related signals in real time without laboratory analysis. However, this technology has a limited spatial coverage. The sensor needs to be in direct contact with the soil or detect at close range (the detection radius is usually in the centimeter to meter range). It can only obtain pollution data at points or small areas. If you want to monitor an area of ​​square kilometers, you need to move the sensor point by point, which is extremely inefficient. At the same time, the results are greatly affected by the physical and chemical properties of the soil (such as moisture, organic matter content, and particle size).

[0005] III. Single Remote Sensing or Spectroscopic Detection Technologies: These technologies rely on a single type of sensor (such as optical remote sensing or hyperspectral imaging) to acquire data. While they can achieve macroscopic coverage, they lack information dimensions. Single data sources cannot simultaneously address both spatial and spectral resolution. For example, Landsat™ remote sensing imagery has a spatial resolution of 30 meters (macroscopic but lacking detail, unable to identify small pollution points), while ground-based hyperspectral instruments have high spectral resolution (able to identify pollutants) but can only acquire point data (unable to capture macroscopic distribution). Furthermore, they have poor anti-interference capabilities, and single data sources are easily affected by the external environment. For instance, optical remote sensing imagery (such as visible-near-infrared) is easily obstructed by clouds and vegetation (if the soil is covered by vegetation, soil pollution cannot be directly detected), resulting in low quantitative accuracy. Single data sources are insufficient to establish accurate spectral-pollutant content inversion models. For example, when using only hyperspectral data to invert heavy metal content, the inversion error is significant due to spectral interference from soil background (such as iron oxides).

[0006] In recent years, with the advancement of data processing technology, multidimensional pollutant data integration analysis has gradually become a research hotspot. Researchers have attempted to comprehensively analyze the concentrations of various pollutants in soil, soil physicochemical parameters, and topographic features to reveal the correlations between pollutants and their diffusion patterns. However, due to the lack of a dynamic update mechanism, these methods cannot reflect the diffusion process and changing trends of pollutants in real time.

[0007] Publication No. CN118067960A discloses a soil environmental pollution monitoring system and method based on multi-source data. The system includes a multi-source data monitoring module, a database, a data acquisition device, an internet platform, and a monitoring terminal. The internet platform includes a first partitioning unit, a second partitioning unit, a first pollution prediction module, a second pollution prediction module, and a monitoring unit anomaly judgment module. This system can improve the accuracy of soil monitoring and reduce manpower and material costs. However, the system still suffers from insufficient consideration of environmental factors and soil characteristics during grid partitioning and prediction unit establishment, affecting the accuracy and reliability of predictions. Summary of the Invention

[0008] The technical problem to be solved by this invention is to provide a data analysis method for monitoring the spatial distribution of soil pollution based on image fusion, which addresses the above-mentioned shortcomings. By using image fusion technology, it achieves comprehensive monitoring of soil pollution, breaks through the limitations of traditional monitoring methods, and can fully reflect the spatial distribution characteristics of pollutants, effectively improving the scope and accuracy of soil pollution monitoring.

[0009] To solve the above technical problems, the present invention adopts the following technical solution: A data analysis method for monitoring the spatial distribution of soil pollution based on image fusion includes the following steps: Step 1: Collect soil pollution monitoring data; Step 2: Preprocess soil pollution monitoring data; Step 3: Perform image fusion calibration, integrating preprocessed multi-source data through multi-scale, multi-modal image fusion technology, while simultaneously performing spatial and spectral calibration to generate a fused image with both high spatial and spectral resolution; Step 4: Extract soil pollutant characteristics, extracting spectral characteristics, spatial distribution characteristics, and pollution intensity characteristics of soil pollutants based on the fused and calibrated image data; Step 5: Generate a soil health index, constructing a multi-dimensional soil health evaluation system by comprehensively considering soil pollutant characteristics, soil physicochemical parameters, and ecological function requirements, calculating the soil health index, and achieving a quantitative assessment of soil health status; Step 6: Conduct soil environmental monitoring, establishing a real-time dynamic monitoring system based on the soil health index and pollutant characteristic extraction results to achieve dynamic tracking, risk warning, and monitoring report output for soil pollution.

[0010] Furthermore, step 1 includes the following steps: Step 1.1 Laboratory chemical analysis data: After collecting soil samples, the concentrations of heavy metals and organic pollutants are detected by instruments, and the latitude, longitude, and sampling depth of each sampling point are recorded. Step 1.2 In-situ sensor detection data: Multi-parameter in-situ sensors are evenly deployed in the monitoring area. The sensors collect pollutant concentration signals and soil physicochemical parameters in real time. The collection frequency is set to once per hour. The data is uploaded to the data platform in real time through the wireless transmission module. At the same time, the sensor location coordinates and collection timestamp are recorded. Step 1.3 Remote sensing and spectral data: Simultaneously acquire multi-source remote sensing and spectral data, including 30m spatial resolution optical remote sensing images from Landsat-8 satellite, 10m spatial resolution multispectral images from Sentinel-2 satellite, hyperspectral data with spatial resolution below 1m acquired by a hyperspectral instrument carried by an UAV, and hyperspectral curves of soil samples collected simultaneously at sampling points by a ground-based portable hyperspectral instrument.

[0011] Furthermore, step 2 includes the following steps: Step 2.1, Preprocessing of laboratory chemical analysis data: Step 2.1.1, outlier removal: For pollutant concentration data from multiple tests of the same sample, outlier judgment is performed, and outlier data caused by experimental operation errors or instrument malfunctions are removed. Step 2.1.2, Data Standardization: Convert the concentration data of different pollutants into a standardized form, using the following formula: ,in Let be the concentration of the j-th pollutant at the i-th sampling point. Let be the average concentration of the j-th pollutant. Let be the standard deviation of the concentration of the j-th pollutant, to eliminate the influence of differences in the order of magnitude of different pollutant concentrations on subsequent analyses; Step 2.2, preprocessing of in-situ sensor detection data; Step 2.2.1, Noise Filtering: For the real-time data collected by the sensor, the moving average filtering method is used to filter high-frequency noise caused by sensor fluctuations and external electromagnetic interference. Step 2.2.2, Data Completion: If data is missing due to transmission interruption or temporary sensor failure, and the missing time is ≤4 hours, linear interpolation is used to complete the short-term missing data; if the missing time exceeds 4 hours, Kriging interpolation is used to complete the data by combining the historical data of the sensor and the data of neighboring sensors to ensure data continuity. Step 2.3, remote sensing and spectral data preprocessing; Step 2.3.1, Remote sensing image preprocessing: For Landsat and Sentinel-2 satellite images, radiometric calibration is performed sequentially to convert the DN values ​​received by the sensor into surface reflectance; atmospheric correction is performed using the FLAASH atmospheric correction model to eliminate the influence of atmospheric scattering and absorption on the image; geometric correction is performed using a quadratic polynomial correction model based on the high-precision topographic map of the monitoring area to convert the image coordinate system to the WGS84 coordinate system, with the correction error controlled within 1 pixel. Step 2.3.2, Hyperspectral data preprocessing: For the hyperspectral images acquired by the UAV, radiometric calibration, atmospheric correction, and geometric correction are performed in sequence. Then, the SIFT feature matching algorithm is used to stitch together multiple images and perform radiometric normalization. Spectral data acquired by ground-based portable hyperspectral analyzers and UAV hyperspectral analyzers were smoothed using Savitzky-Golay filtering; baseline correction was performed using an adaptive iterative reweighted penalized least squares method to eliminate spectral baseline drift; and noise reduction was performed to remove abnormal spectral bands caused by instrument noise and soil particle scattering, retaining the effective spectral range of 400-2500 nm. Step 2.4, Data Format and Coordinate System 1: Convert all preprocessed data sources to GeoTIFF and CSV formats, and set the coordinate system of all spatial data to WGS84 coordinate system. Associate attribute data with corresponding spatial coordinates to ensure the consistency of multi-source data in spatial location.

[0012] Furthermore, step 3 includes the following steps: Step 3.1, the specific process of data layer fusion is as follows: Step 3.1.1: Perform multi-scale wavelet decomposition on the images to be fused, divided into first-level decomposition and second-level decomposition. The preprocessed Landsat-8 image is denoted as image A, and the Sentinel-2 image is denoted as image B. Perform wavelet decomposition with the same parameters on both images, using row-column convolution and downsampling to decompose the images into low-frequency approximation components and high-frequency detail components. Step 3.1.2: Reconstruct the high-frequency detail components of the high spatial resolution image and the low-frequency approximation components of the high spectral resolution image to generate the fused image. Step 3.2, Feature Layer Fusion: The remote sensing data is fused with the in-situ sensor data. The specific process is as follows: Step 3.2.1, Spatial matching: Using the fused remote sensing image as a spatial framework, based on the latitude and longitude coordinates of the in-situ sensor, the pollutant concentration and soil physicochemical parameter data collected by the sensor are embedded as feature points into the remote sensing image to establish the spatial relationship between image pixels and sensor data. Step 3.2.2, Interpolation Supplement: Ordinary Kriging interpolation is used. The spatial correlation of pollutant concentration is analyzed by semi-variogram analysis. In-situ sensor data is used as sample points. Combined with the influence of soil physicochemical parameters on the interpolation results, the pollutant concentration and physicochemical parameters in the remote sensing image without sensor deployment are interpolated and estimated to generate continuous spatial distribution maps of pollutant concentration and soil physicochemical parameters. 3.3, Fusion Calibration: Optimization of Spatial and Spectral Accuracy; Spatial calibration: Select 10-20 ground control points within the monitoring area, compare the coordinate deviations of the control points between the fused image and the high-precision topographic map, and use affine transformation to perform spatial calibration on the fused image to ensure that the spatial position error of the fused image is ≤0.5 pixels; Spectral calibration: Using laboratory chemical analysis data as the true value, select 50-100 sampling points, correlate the spectral values ​​of corresponding pixels in the fused image with the pollutant concentration detected in the laboratory, establish a linear regression model of spectral value-pollutant concentration, and adjust the spectral response coefficient of the fused image to make the determination coefficient of the model ≥0.85.

[0013] Furthermore, step 3.1.1 includes the following steps: Step 3.1.1.1, first-level decomposition, extracting large-scale information; Row direction decomposition: For each row of pixels in image A, perform convolution operations using a low-pass filter and a high-pass filter of the db4 wavelet basis function: Input pixel sequence: Suppose the pixel value sequence of a certain row of the image. Where n is the length of the pixel sequence, the pixel value sequence X is padded with zeros to extend it, and the extended sequence is: Ensure that the window covers all pixels; Calculation process of low-frequency approximate components: Low-pass filter length coefficient For the expanded sequence Using the length of the low-pass filter as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding low-pass filter coefficients is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: The window moves one pixel at a time until it covers the entire extended sequence. The final convolution result has a length of n+7. Perform interval sampling to obtain the low-frequency approximate component Lrow in the row direction; Calculation process of high-frequency detail components: High-pass filter length coefficient Where k is the coefficient index, for the expanded sequence Using the high-pass filter length as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding high-pass filter coefficient is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: For the convolution result Perform interval sampling to obtain the high-frequency detail component Hrow in the row direction; Column-direction decomposition: Repeat convolution and downsampling operations on Lrow and Hrow obtained from row decomposition by column: Convolve and downsample Lrow by column: Obtain a low-frequency approximate component L1; Column-wise convolution and downsampling of Hrow yields one horizontal high-frequency detail component Hh1, one vertical high-frequency detail component Hv1, and one diagonal high-frequency detail component Hd1, respectively. Output results: After the first layer decomposition of image A, {L1_A,Hh1_A,Hv1_A,Hd1_A} are obtained; after the first layer decomposition of image B, {L1_B,Hh1_B,Hv1_B,Hd1_B} are obtained. Step 3.1.1.2, Second-level decomposition: Extracting mesoscale details; Using the low-frequency approximation component L1 obtained from the first-level decomposition as input, the row-column decomposition process is repeated to extract finer details: Low-pass / high-pass filtering convolution and downsampling are performed on each row of L1_A to obtain the low-frequency sequence L1_row and the high-frequency sequence H1_row in the row direction; L1_row and H1_row are filtered and downsampled column by column to obtain two layers of low-frequency approximate components L2 and two layers of high-frequency detail components Hh2, Hv2, and Hd2; the final result of the two-layer decomposition of image A is: {L2_A,Hh2_A,Hv2_A,Hd2_A,Hh1_A,Hv1_A,Hd1_A}; the result of the two-layer decomposition of image B is: {L2_B,Hh2_B,Hv2_B,Hd2_B,Hh1_B,Hv1_B,Hd1_B}.

[0014] Furthermore, step 3.1.2 includes the following steps: Step 3.1.2.1: Complete the replacement of high-frequency detail components according to the decomposition level. The wavelet decomposition level corresponds to the scale of spatial detail: The high-frequency detail components (Hh2, Hv2, Hd2) corresponding to the second layer L2 reflect mesoscale details, while the high-frequency detail components (Hh1, Hv1, Hd1) corresponding to the first layer L1 reflect small-scale details. They need to be replaced in order from the upper layer to the lower layer to ensure that spatial details of different scales are accurately supplemented. The high-frequency detail components of the second layer are replaced to supplement the mesoscale details: the three high-frequency detail components (Hh2_A, Hv2_A, Hd2_A) of the second layer of image A are directly replaced with the corresponding high-frequency detail components (Hh2_B, Hv2_B, Hd2_B) of the second layer of image B. Replacement of high-frequency detail components in the first layer for small-scale detail supplementation: Replace the three high-frequency detail components (Hh1_A, Hv1_A, Hd1_A) in the first layer of image A with the corresponding high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) in the first layer of image B. Low-frequency approximation component preservation: The two low-frequency approximation components (L2_A, L1_A) of image A are completely preserved without replacement, ensuring that the fused image can identify the core of the pollutant's spectral characteristics; After replacement, image A forms a fusion decomposition component set: {L2_A,Hh2_B,Hv2_B,Hd2_B,Hh1_B,Hv1_B,Hd1_B}; Step 3.1.2.2, wavelet inverse transform to generate fused image, is the reverse restoration of the wavelet decomposition sequence, first layer then second layer. It employs row-column convolution and upsampling, specifically divided into two layers of inverse transform to progressively reconstruct the image: Second-level inverse transform: generates mesoscale fused low-frequency approximation components; Using the second-level components (L2_A, Hh2_B, Hv2_B, Hd2_B) in the fusion decomposition component set as input, the operation is performed in the order of column inverse decomposition followed by row inverse decomposition: Inverse column decomposition: Column convolution and upsampling are performed on L2_A, Hh2_B, Hv2_B, and Hd2_B respectively. Using a db4 wavelet basis low-pass or high-pass filter, after convolution operation on each column of each component, the number of columns is doubled by upsampling with zero values ​​inserted at intervals. The specific process is as follows: Step L1: Extract the pixel value sequence of each column from column 1 to column 125 according to the column index order. If the component is a low-frequency approximate component, use the low-pass filter coefficient h for convolution. If the component is a high-frequency detail component, use the high-pass filter coefficient g for convolution. Calculate the column sequence X_col and the filter coefficients in a sliding window manner. The window size is equal to the filter length. The convolution result of each window is the sum of the X_col pixel values ​​in the window × the corresponding filter coefficients. After convolution of each column, a convolution sequence with the same length as the original column is generated. Step L2: After completing the single-column convolution of all input components, the convolution result of the low-frequency approximate component of the same column index is summed with the convolution result of the three high-frequency detail components to obtain the intermediate merged sequence of column inverse decomposition, which integrates the low-frequency contour information and the high-frequency detail information to restore the complete signal of the column in the spatial domain. Step L3: Upsampling is performed on the intermediate merged sequence by inserting zero values ​​at intervals. One zero is inserted between every two adjacent values, doubling the length of the intermediate merged sequence, that is, doubling the column dimension size, matching the spatial resolution of the original image. The upsampled sequence is the final result of the inverse decomposition of the column in the column direction, corresponding to the pixel value sequence of the column in the spatial domain. Step L4: Repeat steps L1-L3 for all columns of the input component. After completing the inverse decomposition calculation of all columns, output the intermediate component with doubled column dimension size. The input component of 125×125 pixels is output as an intermediate component of 250×125 pixels after inverse decomposition in the column direction. Row-direction inverse decomposition: For the result after column inverse decomposition, row-direction repeated convolution and upsampling are performed according to the column-direction inverse decomposition process, doubling the number of rows from 125 to 250, and finally generating the first layer of fused low-frequency approximation component L1_fuse with a size of 250×250 pixels. At this time, L1_fuse has integrated the mesoscale spectral information of image A with the mesoscale spatial details of image B. First inverse transform: Generates the final fused image. Using L1_fuse and the first layer high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) as input, repeat the inverse transform process of the second layer inverse transform: Inverse column decomposition: Perform column convolution and upsampling on L1_fuse, Hh1_B, Hv1_B, and Hd1_B to double the number of columns from 250 to 500; Row-direction inverse decomposition: The column inverse decomposition results are subjected to row convolution and upsampling to double the number of rows from 250 to 500, finally generating a fused image with a size of 500×500 pixels. The spatial resolution is consistent with that of image B, and the spectral resolution is consistent with that of image A.

[0015] Furthermore, step 3.2.1 includes the following steps: Step 3.2.1.1, Coordinate System Verification: First, confirm the consistency of the coordinate systems of the two types of input data. The fused image has been unified into the WGS84 coordinate system. The in-situ sensor data has been recorded with WGS84 latitude and longitude coordinates by GPS positioning during deployment. It is necessary to use spatial analysis software to import the sensor coordinates into the spatial framework of the fused image and verify whether the coordinates of each sensor point fall within the corresponding pixel range. Step 3.2.1.2, Pixel-Sensor Data Association: For each in-situ sensor point, extract the spectral feature information of its corresponding pixel in the fused image, and bind the pollutant concentration data and soil physicochemical parameters of the sensor with the extracted spectral information to form a spatial location-spectral feature-quantitative parameter association dataset. Step 3.2.1.3, outlier screening: Spatial outlier detection is performed on the associated dataset to determine whether it is an outlier. Outliers need to be removed or data needs to be re-collected to avoid affecting the accuracy of subsequent interpolation.

[0016] Furthermore, step 3.2.2 includes the following steps: Step 3.2.2.1, Sample point data screening; Basic data extraction: Extract core data from all valid sensors in the monitoring area from the in-situ sensor database, including: spatial coordinates: sensor latitude and longitude, unified in the WGS84 coordinate system, consistent with the fused remote sensing image; target variable values: pollutant concentration and soil physicochemical parameters to be interpolated; auxiliary variable values: original values ​​of soil physicochemical parameters collected by the corresponding sensors, used as covariates for interpolation to correct spatial correlation. Step 3.2.2.2, modeling the spatial variogram; Step 3.2.2.3: Correct the variogram by combining soil physicochemical parameters; Step 3.2.2.4: Interpolation grid generation. Using the fused remote sensing image as the spatial frame, an interpolation grid is generated that corresponds one-to-one with the image pixels, ensuring that the interpolation results can be directly superimposed on the remote sensing image. Grid parameter settings: The grid resolution should be consistent with the remote sensing image. Grid coordinate generation: Based on the latitude and longitude of the upper left and lower right corners of the remote sensing image, calculate the center point coordinates (x_f, y_f) of each grid cell as the point to be predicted, where f is the grid cell number, to ensure that all grid cells cover the entire monitoring area without omissions or overlaps; Step 3.2.2.5, Ordinary Kriging interpolation calculation: For the center point of each grid cell, calculate the estimate of its target variable based on the modified variogram. Step 3.2.2.6: Spatial distribution map generation. The interpolated grid data is converted into a visual spatial distribution map and overlaid with the fused remote sensing image.

[0017] Furthermore, the modeling process of the spatial variogram function is as follows: The formula for the variogram function is: for the target variable In spatial location and The variogram function is defined as the value at a given position: ; The variogram value with a lag distance of g is used to quantify the spatial variability of the target variable; g is the spatial lag distance, which represents the spatial distance between two sample points. When the lag distance is g, the total number of sample point pairs within the monitoring area that meet this distance condition; the index number of sample point pair i, with a value ranging from 1 to... This is used to iterate through all sample point pairs that satisfy the lag distance g. The spatial location of the i-th sample point; The measured value of the target variable at the i-th sample point; The spatial location of the paired sample point that is at a distance g from the i-th sample point; The measured value of the target variable at the paired sample point with a distance of g from the i-th sample point; Based on the pollutant diffusion characteristics of the monitoring area, a suitable theoretical model is selected to fit the experimental variogram: If the pollutants diffuse uniformly, choose the spherical model: When the hysteresis distance g is less than or equal to the range a: When the lag distance g is greater than the range a: If pollutants are blocked by terrain, choose the index model: ;in, denoted as nugget value, reflecting random error; C is the sill value, reflecting spatial structural variation; a is the range, the maximum distance of spatial correlation, beyond which there is no spatial correlation between sample points.

[0018] Furthermore, the specific process of step 3.2.2.5 is as follows: To reduce computational cost and ensure interpolation accuracy, a variable radius neighborhood selection method is used for selecting neighborhood sample points. With the desired center point (x_f, y_f) as the center, the initial search radius is set to the range a of the mutation function; Within the search radius, sample points are selected. If the number of sample points is less than 10, the radius is gradually expanded by 10% of a each time until the number of sample points is greater than or equal to 10. If the number of sample points is greater than 30, the 30 sample points closest to the point to be predicted are retained to avoid interference from edge sample points. Based on the variogram, a weight matrix for sample points to prediction points is constructed. This integrates the spatial location information of sample points with auxiliary information on soil physicochemical parameters. The aim is to adapt to the complex heterogeneity of the soil environment within the monitoring area and satisfy two constraints: Unbiasedness constraint: ; Minimum variance constraint: The weights are solved using the Kriging equations, which are in the following form: ; in, Must meet , The variogram value between the i-th and j-th neighboring sample points is calculated using a theoretical model modified with soil physicochemical parameters, reflecting the spatial correlation of the target variable between the two points. The Lagrange multiplier is an auxiliary parameter set to satisfy the unbiasedness constraint of interpolation. The spatial location of the point to be predicted The i-th neighboring sample point and the point to be predicted The variogram values ​​between the sample points and the points to be predicted quantify the spatial correlation between them. The number of neighborhood sample points participating in the interpolation calculation is determined by the variable radius neighborhood selection method: with the point to be predicted as the center, the initial search radius is set to the range 'a' of the variogram, and the neighborhood sample points are selected... One effective neighborhood sample point, The value range is 10-30 to ensure a balance between interpolation accuracy and computational efficiency. The weight of each sample point is obtained by solving the system of equations through matrix inversion. ; Calculate the estimated value of the point to be predicted based on the measured values ​​and weights of the sample points: , Points to be predicted The interpolation results of the target variable, i.e. the estimated values ​​of the corresponding points in the grid cells, are used to generate a continuous spatial distribution map; The measured value of the target variable at the i-th sample point comes from the effective monitoring data of the in-situ sensor.

[0019] The present invention adopts the above technical solution and has the following technical effects compared with the prior art: 1. This invention achieves comprehensive monitoring of soil pollution through image fusion technology, breaking through the limitations of traditional monitoring methods. It can comprehensively reflect the spatial distribution characteristics of pollutants and effectively improve the scope and accuracy of soil pollution monitoring. 2. This invention employs advanced image processing algorithms to conduct in-depth analysis of soil environmental factors and characteristics, optimizes grid division and prediction unit establishment, improves prediction accuracy and reliability, and overcomes the problem of insufficient consideration of environmental factors in existing technologies; 3. This invention utilizes image fusion technology to achieve real-time monitoring of soil pollution status and establishes a dynamic update mechanism, which can promptly reflect the diffusion process and changing trend of pollutants, meeting the needs of modern environmental protection supervision for real-time monitoring of soil pollution status. 4. This invention optimizes the drawing and acquisition of pollution contour lines through an innovative image fusion method, improves the accuracy of pollution area and volume calculation, and overcomes the problem of unreasonable pollution contour line processing in the prior art. 5. This invention effectively integrates multi-source data and utilizes advanced image processing technology to achieve effective utilization and integration of data from different sensors and monitoring systems, significantly improving the comprehensiveness and accuracy of soil pollution monitoring. Attached Figure Description

[0020] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the accompanying drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. In all the drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn to scale.

[0021] Figure 1 This is a flowchart of the data analysis method for monitoring the spatial distribution of soil pollution based on image fusion according to the present invention; Figure 2 This is a flowchart of step 2 of the present invention, which involves preprocessing soil pollution monitoring data. Figure 3 This is a flowchart of step 3, image fusion calibration, of the present invention; Figure 4 This is a flowchart of step 4 of the present invention, which involves extracting soil pollutant characteristics. Figure 5 This is a flowchart of step 3.1, data layer fusion, of the present invention; Figure 6 This is a flowchart of step 3.1.1 of the present invention, which involves multi-scale wavelet decomposition of the images to be fused. Figure 7 This is a flowchart of step 3.1.2 of the present invention, which involves reconstructing the image to generate a fused image; Figure 8 This is a flowchart of step 3.2, feature layer fusion, of the present invention; Figure 9 This is a flowchart illustrating step 3.2.1 of the present invention, which involves establishing spatial association between image pixel and sensor data through spatial matching. Figure 10 This is a flowchart of step 3.2.2 of the present invention, which uses ordinary Kriging interpolation. Detailed Implementation

[0022] Examples, such as Figure 1 As shown, a data analysis method for monitoring the spatial distribution of soil pollution based on image fusion includes: Step 1: Collect soil pollution monitoring data. To ensure the comprehensiveness and accuracy of the data for subsequent analysis, this step requires systematically collecting multi-dimensional and multi-source soil pollution-related data: Step 1.1 Laboratory chemical analysis data; Sampling principles: A combination of grid-based sampling and key-area densification is adopted. Within the monitoring area, sampling units are divided according to a basic grid of 1km×1km, and basic sampling is carried out at the grid intersections. For potentially high-pollution areas such as industrial sites, the vicinity of chemical industrial parks, and sewage irrigation areas, the sampling grid is densified to 200m×200m to ensure accurate capture of pollution hotspots.

[0023] Detection indicators: After collecting soil samples, the concentrations of heavy metals (such as cadmium, mercury, arsenic, lead, and chromium) and organic pollutants (such as polycyclic aromatic hydrocarbons, pesticide residues, and volatile organic compounds) were detected by instruments such as high performance liquid chromatography-mass spectrometry (HPLC-MS) and inductively coupled plasma mass spectrometry (ICP-MS). At the same time, the spatial attribute information of each sampling point, such as latitude and longitude and sampling depth (0-20cm topsoil, 20-60cm middle soil, and below 60cm deep soil), was recorded.

[0024] Step 1.2 In-situ sensor detection data; Sensor Deployment: Multi-parameter in-situ sensors are evenly deployed within the monitoring area. The sensor spacing is set to 500m-1km according to the monitoring accuracy requirements, and can be shortened to 100-200m in heavily polluted areas. Sensor types include soil heavy metal ion selective electrode sensors, organic pollutant optical sensors, soil moisture sensors, organic matter content sensors, and soil particle size sensors.

[0025] Data Acquisition: The sensor collects pollutant concentration signals and soil physicochemical parameters in real time. The pollutant concentration signals include heavy metal ion activity and characteristic spectral intensity of organic pollutants. The soil physicochemical parameters include soil moisture, organic matter content, and particle size. The acquisition frequency is set to once per hour. The data is uploaded to the data platform in real time through a wireless transmission module, and the sensor location coordinates and acquisition timestamp are recorded at the same time.

[0026] Step 1.3 Remote sensing and spectral data; Data types: Simultaneous acquisition of multi-source remote sensing and spectral data, including 30m spatial resolution optical remote sensing images from the Landsat-8 satellite for macro-area coverage; 10m spatial resolution multispectral images from the Sentinel-2 satellite for supplementing mesoscale details, such as the edges of small pollution spots at the 10m level, micro-topography around sampling points, and the boundary between industrial sites and farmland; hyperspectral data with spatial resolution below 1m acquired by a UAV-mounted hyperspectral instrument for refined detection of key areas; and hyperspectral curves of soil samples simultaneously collected at sampling points by a ground-based portable hyperspectral instrument for establishing a spectral-pollutant content correlation model.

[0027] Data acquisition time: Select cloudless or lightly cloudy weather conditions to acquire remote sensing data, ensuring that the images are not obscured by clouds; keep the acquisition time of UAV hyperspectral data consistent with that of ground sampling and in-situ sensor data, with a maximum time difference of no more than 24 hours, to avoid reduced data correlation due to time differences.

[0028] While ensuring comprehensive monitoring, the sampling density in highly polluted areas is increased by more than 5 times, which can accurately capture small pollution spots at the 10m level and avoid missing pollution hotspots due to sparse sampling.

[0029] Multi-source data spatiotemporal synchronization eliminates correlation deviations caused by time misalignment between different data sources, providing a high-quality data foundation for subsequent fusion analysis, reducing the proportion of invalid data, and improving overall analysis efficiency.

[0030] like Figure 2 As shown, step 2 involves preprocessing soil pollution monitoring data; Step 2.1, Preprocessing of laboratory chemical analysis data: Step 2.1.1, Outlier Removal: Outlier analysis is performed on the pollutant concentration data from multiple tests of the same sample. Outliers caused by experimental errors or instrument malfunctions are removed, as detailed below: Arrange the datasets in ascending order: y1 is the minimum value and ym is the maximum value, both of which are potential outliers; Calculate basic statistics: Calculate the sample mean ; Calculate the sample standard deviation ; Calculate the Grubbs statistic G: If the minimum value y1 is suspected to be an outlier: If the maximum value ym is suspected to be an outlier: Based on the significance level α and sample size m, with α = 0.05, the critical value G is obtained by consulting the Grubbs critical value table. α,m .

[0031] If the calculated G > G α,m If G ≤ G, then at a significance level of α, this value is an outlier and can be removed; α,m If the value is not an outlier, it cannot be considered an outlier and must be retained. After removing the confirmed outliers, repeat the above steps for the remaining data until there are no outliers.

[0032] If the pollutant concentration data at a certain sampling point exceeds the background value of the soil in that area by more than 10 times, it is necessary to combine the on-site investigation (such as whether there is a pollution source) to determine whether it is truly high pollution data, so as to avoid erroneous rejection.

[0033] Step 2.1.2, Data Standardization: Convert the concentration data of different pollutants into a standardized form, using the following formula: ,in Let be the concentration of the j-th pollutant at the i-th sampling point. Let be the average concentration of the j-th pollutant. Let be the standard deviation of the concentration of the j-th pollutant, to eliminate the influence of differences in the concentration magnitudes of different pollutants on subsequent analyses.

[0034] Outliers are removed by using the Grubbs method, and the authenticity of data with values ​​more than 10 times the background value is determined by on-site investigation, thus avoiding the accidental deletion of real high-pollution data. Standardized formulas are used to eliminate differences in the concentration levels of different pollutants, so that the concentration data of heavy metals such as cadmium and mercury and organic pollutants such as polycyclic aromatic hydrocarbons can be directly used in comprehensive analysis, solving the analytical bias problem caused by different units and numerical ranges of pollutant concentrations in traditional methods.

[0035] Step 2.2, preprocessing of in-situ sensor detection data; Step 2.2.1, Noise Filtering: For the real-time data collected by the sensor, the moving average filtering method is adopted, and the window size is set to 5, that is, the average value of the data at 5 consecutive time points is taken to filter out high-frequency noise caused by sensor fluctuations and external electromagnetic interference. For parameters that change slowly, such as soil moisture and organic matter content, the window size can be appropriately increased to 10 to further smooth the data.

[0036] Step 2.2.2, Data Completion: If data is missing due to transmission interruption or temporary sensor failure, and the missing time is ≤4 hours, linear interpolation is used to complete the short-term missing data; if the missing time exceeds 4 hours, Kriging interpolation is used to complete the data by combining the historical data of the sensor and the data of neighboring sensors to ensure data continuity.

[0037] Moving average filtering effectively filters high-frequency noise. Linear interpolation (≤4 hours missing) and Kriging interpolation (>4 hours missing) are combined with historical data from the same period and surrounding areas to complete the data, ensuring that the continuity of sensor data reaches ≥95%, avoiding monitoring gaps caused by data interruptions, and providing continuous and stable point data support for feature layer fusion.

[0038] Step 2.3, remote sensing and spectral data preprocessing; Step 2.3.1, Remote sensing image preprocessing: For Landsat and Sentinel-2 satellite images, radiometric calibration is performed sequentially to convert the DN values ​​received by the sensor into surface reflectance; atmospheric correction is performed using the FLAASH atmospheric correction model to eliminate the influence of atmospheric scattering and absorption on the image; geometric correction is performed using a quadratic polynomial correction model based on the high-precision topographic map of the monitoring area to convert the image coordinate system to the WGS84 coordinate system, and the correction error is controlled within 1 pixel.

[0039] Step 2.3.2, Hyperspectral data preprocessing: For the hyperspectral images acquired by the UAV, radiometric calibration, atmospheric correction, and geometric correction are performed in sequence. In addition, multiple images need to be stitched together using the SIFT feature matching algorithm, and radiometric normalization is performed to eliminate radiometric differences caused by different flight sorties and lighting conditions.

[0040] Spectral data acquired by ground-based portable hyperspectral analyzers and UAV hyperspectral analyzers were smoothed using Savitzky-Golay filtering with a window size of 7 points and a polynomial order of 2. Baseline correction was performed using an adaptive iterative reweighted penalized least squares method to eliminate spectral baseline drift. Noise reduction was performed to remove abnormal spectral bands caused by instrument noise and soil particle scattering, retaining the effective spectral range of 400-2500 nm.

[0041] Operations such as FLAASH atmospheric correction, quadratic polynomial geometric correction (error ≤ 1 pixel), and Savitzky-Golay filtering (7-point window, second-order polynomial) eliminate interference from atmospheric scattering, coordinate bias, and spectral drift, improving the signal-to-noise ratio of the effective band (400-2500nm) of hyperspectral data by more than 30%, ensuring that the spectral and spatial accuracy of remote sensing data meets the fusion requirements.

[0042] Step 2.4, Data Format and Coordinate System 1: Convert all preprocessed data sources (laboratory data, in-situ sensor data, remote sensing and spectral data) into GeoTIFF format (spatial data) and CSV format (attribute data), and set the coordinate system of all spatial data to WGS84 coordinate system. Associate the attribute data with the corresponding spatial coordinates (latitude and longitude) to ensure the consistency of multi-source data in spatial location.

[0043] like Figure 3 As shown, step 3 involves image fusion calibration. This is achieved by integrating preprocessed multi-source data using multi-scale, multi-modal image fusion technology, while simultaneously performing spatial and spectral calibration to generate a fused image with both high spatial and spectral resolution. The specific process is as follows: Step 3.1, Data Layer Fusion: Multi-source remote sensing data fusion; A wavelet transform fusion algorithm is used to fuse Landsat images, with high spatial resolution Sentinel-2 imagery or UAV hyperspectral imagery as the spatial reference and high spectral resolution ground hyperspectral data as the spectral reference. The specific process is as follows: like Figure 5As shown in step 3.1.1, multi-scale wavelet decomposition is performed on the images involved in the fusion, which is divided into a first-level decomposition and a second-level decomposition. The preprocessed Landsat-8 image is denoted as image A and the Sentinel-2 image is denoted as image B. Wavelet decomposition with the same parameters is performed on each image. Row-column convolution and downsampling are used to decompose the image into low-frequency approximation components and high-frequency detail components. The low-frequency approximation components reflect the overall contour of the image, and the high-frequency detail components reflect the spatial details of the image.

[0044] like Figure 6 As shown, step 3.1.1.1, the first layer decomposition, extracts large-scale information; Line direction decomposition: For each row of pixels in image A, convolution is performed using a low-pass filter and a high-pass filter based on the db4 wavelet basis function. By low-pass filtering and downsampling the pixel value sequence of each row, the low-frequency approximate component Lrow in the row direction is obtained, representing the overall trend or smoothing information of the pixels in that row. Similarly, by high-pass filtering and downsampling the pixel value sequence of each row, the high-frequency detail component Hrow in the row direction is obtained, reflecting detail information, edge features, and rapidly changing signals in the row direction. Input pixel sequence: Suppose the pixel value sequence of a certain row of the image. Where n is the length of the pixel sequence, the pixel value sequence X is padded with zeros to avoid calculation errors at boundary pixels. The expanded sequence is: Ensure that the window covers all pixels; Calculation process of low-frequency approximate components: Low-pass filter length coefficient For the expanded sequence Using the length of the low-pass filter as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding low-pass filter coefficients is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: The window moves one pixel at a time until it covers the entire extended sequence. The final convolution result has a length of n+7. Perform interval sampling and take the value at the even index position to obtain the row direction low frequency approximate component Lrow, with a length of 1 / 2 of the original sequence (rounded up).

[0045] Calculation process of high-frequency detail components: High-pass filter length coefficient Where k is the coefficient index (0-7), for the expanded sequence Using the high-pass filter length as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding high-pass filter coefficient is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: For the convolution result Perform interval sampling, take the value at even index positions to obtain the high-frequency detail component Hrow in the row direction, with a length of 1 / 2 (rounded up) of the original sequence.

[0046] Column-direction decomposition: Repeat convolution and downsampling operations on Lrow and Hrow obtained from row decomposition by column: Convolve and downsample Lrow by column: Obtain a low-frequency approximate component L1; Convolution and downsampling of Hrow by column yields one layer of horizontal high-frequency detail components Hh1, one layer of vertical high-frequency detail components Hv1, and one layer of diagonal high-frequency detail components Hd1. Hh1 reflects horizontal details, such as the edge of a horizontal road; Hv1 reflects vertical details, such as the boundary of a vertical river; and Hd1 reflects diagonal details, such as diagonally distributed pollution points. Output results: After the first layer decomposition of image A, we get {L1_A,Hh1_A,Hv1_A,Hd1_A}; after the first layer decomposition of image B, we get {L1_B,Hh1_B,Hv1_B,Hd1_B}.

[0047] Step 3.1.1.2, Second-level decomposition: Extracting mesoscale details; Using the low-frequency approximation component L1 obtained from the first-level decomposition as input, the row-column decomposition process is repeated to extract finer details: Perform low-pass / high-pass filtering convolution and downsampling on each row of L1_A (250×250 pixels) to obtain the low-frequency sequence L1_row and the high-frequency sequence H1_row in the row direction; Filter and downsample L1_row and H1_row column by column: Two low-frequency approximation components L2 (size 125×125 pixels, reflecting a more macroscopic outline of the image, such as the land use type distribution of the entire monitoring area) were obtained. Two high-frequency detail components, Hh2, Hv2, and Hd2 (125×125 pixels in size, reflecting mesoscale details, such as the edges of soil pollution patches in the range of 30-50m), were obtained.

[0048] Output results: The final result of the 2-layer decomposition of image A is: {L2_A,Hh2_A,Hv2_A,Hd2_A,Hh1_A,Hv1_A,Hd1_A}; The result of the 2-layer decomposition of image B is: {L2_B,Hh2_B,Hv2_B,Hd2_B,Hh1_B,Hv1_B,Hd1_B}.

[0049] Through multi-scale wavelet decomposition, the contours and details of remote sensing images at different resolutions can be effectively separated, providing a high-quality decomposition basis for subsequent fusion reconstruction (combining the low-frequency approximation components of high spectral resolution with the high-frequency detail components of high spatial resolution), and finally generating a fused image of soil pollution monitoring that has both macroscopic coverage and precise details.

[0050] like Figure 5 As shown in step 3.1.2, the high-frequency detail components of the high spatial resolution image and the low-frequency approximation components of the high spectral resolution image are reconstructed to generate a fused image. This fused image retains both high spectral resolution, enabling the identification of pollutant characteristic spectra, and high spatial resolution, allowing the presentation of details of small-area pollution points. like Figure 7 As shown in step 3.1.2.1, the high-frequency detail components are replaced according to the decomposition level. The wavelet decomposition level corresponds to the scale of spatial detail: The high-frequency detail components (Hh2, Hv2, Hd2) corresponding to the second layer L2 reflect mesoscale details (such as the outline of pollution patches in the range of 30-50m), while the high-frequency detail components (Hh1, Hv1, Hd1) corresponding to the first layer L1 reflect small-scale details (such as the edge of pollution points and the location of sampling points in the range of 10-30m). The replacement should be completed in order from high to low layers to ensure that spatial details of different scales are accurately supplemented. Image A is a high-spectral-resolution image, with the advantage of spectral information represented by low-frequency approximation components. Its high-frequency detail components include instrument noise and soil particle scattering interference, resulting in low spatial resolution and inability to reflect small-scale spatial details. Image B is a high spatial resolution image, with the advantage of refined spatial information in high-frequency detail components, which can accurately locate small-scale pollution areas. Its low-frequency approximation components are invalid outlines. Replacing the inferior high-frequency detail components of image A with the superior high-frequency detail components of image B can achieve the integration effect of removing the inferior and retaining the superior. The high-frequency detail components of the second layer are replaced to supplement the mesoscale details: the three high-frequency detail components (Hh2_A, Hv2_A, Hd2_A) of the second layer of image A are directly replaced with the corresponding high-frequency detail components (Hh2_B, Hv2_B, Hd2_B) of the second layer of image B. Replacement of high-frequency detail components in the first layer for small-scale detail supplementation: Replace the three high-frequency detail components (Hh1_A, Hv1_A, Hd1_A) in the first layer of image A with the corresponding high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) in the first layer of image B. Low-frequency approximation component preservation: The two low-frequency approximation components (L2_A, L1_A) of image A are completely preserved without replacement, ensuring that the fused image can identify the core of the pollutant's spectral characteristics; After replacement, image A forms a fusion decomposition component set: {L2_A, Hh2_B, Hv2_B, Hd2_B, Hh1_B, Hv1_B, Hd1_B}.

[0051] Step 3.1.2.2, wavelet inverse transform to generate fused image, is the reverse restoration of the wavelet decomposition sequence, first layer, then second layer. Wavelet inverse transform is the process of restoring the fused decomposition component set from the frequency domain to the spatial domain image, ensuring that spatial details are gradually supplemented from macroscopic to microscopic, avoiding image reconstruction distortion due to scale misalignment, and ultimately generating a high-fidelity fused image that meets the needs of soil pollution monitoring. Row-column convolution and upsampling are used, specifically divided into two layers of inverse transform to progressively reconstruct the image. Second-level inverse transform: generates mesoscale fused low-frequency approximation components; Using the second-level components (L2_A, Hh2_B, Hv2_B, Hd2_B) in the fusion decomposition component set as input, the operation is performed in the order of column inverse decomposition followed by row inverse decomposition: Inverse column decomposition: Column convolution and upsampling are performed on L2_A (125×125 pixels), Hh2_B, Hv2_B, and Hd2_B respectively. Using a db4 wavelet basis low-pass or high-pass filter, after convolution operation on each column of each component, the number of columns is doubled (from 125 to 250) by inserting zero values ​​at intervals during upsampling. The specific process is as follows: Step L1: Extract the pixel value sequence of each column from column 1 to column 125 according to the column index order. If the component is a low-frequency approximate component, use the low-pass filter coefficient h for convolution. If the component is a high-frequency detail component, use the high-pass filter coefficient g for convolution. Calculate the column sequence X_col and the filter coefficients in a sliding window manner. The window size is equal to the filter length. The convolution result of each window is the sum of the X_col pixel values ​​in the window × the corresponding filter coefficients. After convolution of each column, a convolution sequence with the same length as the original column is generated. Step L2: After completing the single-column convolution of all input components, the convolution result of the low-frequency approximate component of the same column index is summed with the convolution result of the three high-frequency detail components to obtain the intermediate merged sequence of column inverse decomposition, which integrates the low-frequency contour information and the high-frequency detail information to restore the complete signal of the column in the spatial domain. Step L3: Upsampling is performed on the intermediate merged sequence by inserting zero values ​​at intervals. One zero is inserted between every two adjacent values, doubling the length of the intermediate merged sequence, that is, doubling the column dimension size, matching the spatial resolution of the original image. The upsampled sequence is the final result of the inverse decomposition of the column in the column direction, corresponding to the pixel value sequence of the column in the spatial domain. Step L4: Repeat steps L1-L3 for all columns of the input component. After completing the inverse decomposition calculation of all columns, output the intermediate component with doubled column dimension size. The input component of 125×125 pixels is output as an intermediate component of 250×125 pixels after inverse decomposition in the column direction.

[0052] Row-direction inverse decomposition: For the results after column inverse decomposition, row-direction repeated convolution and upsampling are performed according to the column-direction inverse decomposition process, doubling the number of rows from 125 to 250, and finally generating the first layer of fused low-frequency approximation component L1_fuse with a size of 250×250 pixels. At this time, L1_fuse has integrated the mesoscale spectral information of image A with the mesoscale spatial details of image B.

[0053] First inverse transform: Generates the final fused image. Using L1_fuse and the first layer high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) as input, repeat the inverse transform process of the second layer inverse transform: Inverse column decomposition: Perform column convolution and upsampling on L1_fuse (250×250 pixels), Hh1_B, Hv1_B, and Hd1_B to double the number of columns from 250 to 500.

[0054] Row-direction inverse decomposition: The column inverse decomposition results are subjected to row convolution and upsampling to double the number of rows from 250 to 500, finally generating a fused image with a size of 500×500 pixels. The spatial resolution is consistent with image B, such as 10m; the spectral resolution is consistent with image A, such as nm level.

[0055] The fused image generated by the reconstruction process perfectly solves the shortcomings of traditional single images. High spectral resolution is preserved: the fused image inherits the spectral core of image A, which can accurately identify pollutant characteristics. For example, it can clearly capture the "double absorption peaks of polycyclic aromatic hydrocarbons at 1780nm and 2200nm", distinguish the spectral differences between ordinary organic matter and organic pollutants, and the error of pollutant concentration inversion is ≤10%, far lower than the 30% error of single remote sensing images. High spatial resolution is improved: the fused image supplements the high-frequency details of image B, which can present small-scale pollution points. For example, it can identify small-area pollution spots at the scale of 10m (such as temporary chemical waste dumping sites), while traditional hyperspectral images (such as ground hyperspectral instruments) can only acquire point data and cannot present their spatial distribution. Traditional remote sensing images (such as Landsat-8) with a resolution of 30m will mix the pollution spot with the surrounding soil and cannot identify it.

[0056] The fused image generated by data layer fusion not only solves the problem that traditional Landsat-8 imagery (30m resolution) cannot identify small-area pollution points, but also makes up for the limitation that ground hyperspectral instruments can only acquire point data. It can clearly capture the characteristic spectra of pollutants, and the error of pollutant concentration inversion is ≤10%, which is far lower than the 30% error of single remote sensing imagery. It provides high-precision image data for subsequent pollutant type identification and concentration inversion, enabling pollution monitoring to be upgraded from qualitative to quantitative.

[0057] Step 3.2, Feature layer fusion: Fusion of remote sensing data and in-situ sensor data; like Figure 8 As shown, step 3.2.1, spatial matching: using the fused remote sensing image as a spatial framework, and based on the latitude and longitude coordinates of the in-situ sensor, the pollutant concentration and soil physicochemical parameter data collected by the sensor are embedded as feature points into the remote sensing image to establish the spatial association between image pixels and sensor data. The specific process is as follows: like Figure 9 As shown, step 3.2.1.1, coordinate system verification: First, confirm the consistency of the coordinate systems of the two types of input data. The fused image has been unified into the WGS84 coordinate system. The in-situ sensor data has been recorded with WGS84 latitude and longitude coordinates via GPS positioning during deployment (positioning accuracy ≤1m). It is necessary to use spatial analysis software such as ArcGIS or ENVI to import the sensor coordinates into the spatial framework of the fused image and verify whether the coordinates of each sensor point fall within the corresponding pixel range. A deviation of ±0.5 pixels is allowed. If it exceeds the range, the sensor coordinates need to be recalibrated. Step 3.2.1.2, Pixel-Sensor Data Association: For each in-situ sensor point, extract the spectral feature information of its corresponding pixel from the fused image. Simultaneously, bind the sensor's pollutant concentration data and soil physicochemical parameters with the extracted spectral information to form a spatial location-spectral feature-quantitative parameter associated dataset. For example: Sensor point A (latitude and longitude: 116.32°E, 39.95°N) corresponds to the 200th row and 300th column pixel of the fused image. The cadmium characteristic band (2345nm) of this pixel is extracted, and the reflectance is 0.32. At the same time, the actual measured cadmium concentration of the sensor is 0.9mg / kg, soil moisture is 16%, and organic matter content is 22g / kg.

[0058] Step 3.2.1.3, Outlier Detection: Spatial outlier detection is performed on the associated dataset. If the spectral value of the image pixel corresponding to a certain sensor point differs too much from the spectral value of the surrounding pixels, such as exceeding 3 times the standard deviation; or if the measured concentration of the sensor differs significantly from the concentration of the surrounding sensors, such as exceeding 5 times the regional background value; it is necessary to combine on-site investigation, such as whether the sensor is faulty or whether there is a local pollution source; to determine whether it is an outlier. Outliers need to be removed or data needs to be re-collected to avoid affecting the accuracy of subsequent interpolation.

[0059] like Figure 8 As shown in step 3.2.2, interpolation supplementation: Ordinary Kriging interpolation is adopted, which differs from traditional linear interpolation and inverse distance weighted interpolation. Its advantage lies in its ability to combine the influence of soil physicochemical parameters on pollutant distribution, analyze the spatial correlation of pollutant concentration through semi-variogram analysis, and improve interpolation accuracy. It is especially suitable for spatially structured data such as soil pollution. This method uses in-situ sensor data as sample points, combines the influence of soil physicochemical parameters on the interpolation results, and performs interpolation estimation of pollutant concentration and physicochemical parameters in areas of remote sensing image where no sensors are deployed. It generates continuous spatial distribution maps of pollutant concentration and soil physicochemical parameters, realizing feature layer fusion of remote sensing image and in-situ data. The specific process is as follows: like Figure 10 As shown, step 3.2.2.1 involves filtering sample point data; Basic data extraction: Extract core data from all valid sensors within the monitoring area from the in-situ sensor database, including: Spatial coordinates: sensor latitude and longitude, unified in the WGS84 coordinate system, consistent with the fused remote sensing imagery; Target variable values: pollutant concentration and soil physicochemical parameters to be interpolated; Auxiliary variable values: These correspond to the original values ​​of soil physicochemical parameters collected by the sensors and are used as covariates for interpolation to correct spatial correlation.

[0060] Step 3.2.2.2, modeling the spatial variogram; The formula for the variogram function is: for the target variable In spatial location and The variogram function is defined as the value at a given position: ; The variogram value with a lag distance of g is used to quantify the spatial variability of the target variable; g is the spatial lag distance, which represents the spatial distance between two sample points. When the lag distance is g, the total number of sample point pairs within the monitoring area that meet this distance condition; the index number of sample point pair i, with a value ranging from 1 to... This is used to iterate through all sample point pairs that satisfy the lag distance g. The spatial location of the i-th sample point; The measured value of the target variable at the i-th sample point; The spatial location of the paired sample point that is at a distance g from the i-th sample point; The measured value of the target variable at the paired sample point at a distance of g from the i-th sample point.

[0061] Based on the pollutant diffusion characteristics of the monitoring area, such as industrial site pollution spreading in concentric circles and agricultural pollution spreading in strips, a suitable theoretical model is selected to fit the experimental variogram function: If pollutants are dispersed evenly, such as in plains or farmland, a spherical model is preferred. When the hysteresis distance g is less than or equal to the range a: ; When the hysteresis distance g is greater than the range a: ; If pollutants are blocked by terrain, such as mountains, choose the index model: ; in, C is the nugget value, reflecting random errors, such as differences in sensor accuracy; C is the sill value, reflecting spatial structural variations, such as the spatial patterns of pollutant diffusion; a is the range, the maximum distance of spatial correlation, beyond which there is no spatial correlation between sample points.

[0062] Step 3.2.2.3: Correct the variogram by combining soil physicochemical parameters; Zonal variation analysis: The monitoring area is divided into different sub-regions according to the key thresholds of soil physicochemical parameters, for example: Soil moisture is categorized into high humidity zone (>25%), medium humidity zone (15%-25%), and low humidity zone (<15%). Based on organic matter content, the zones are divided into high organic matter zone (>30g / kg), medium organic matter zone (15-30g / kg), and low organic matter zone (<15g / kg). Calculate the experimental variogram for each subregion separately. If the nugget effect ratio of the variogram in a certain subregion is... A value below 0.25 indicates strong spatial correlation of the target variables within the region, requiring separate modeling to improve interpolation accuracy.

[0063] Using standardized soil physicochemical parameters as covariates, the sill value C of the variogram is corrected using the following formula: ,in By combining the sill values ​​of the variogram after covariate correction, the spatial structural variation of the target variable can be quantified more accurately. The initial sill value obtained by fitting the experimental variogram without combining covariates; Φ is the number of soil physicochemical parameters involved in the correction. For example, when soil moisture and organic matter content are corrected simultaneously, Φ=2. If soil particle size is added, then Φ=3. Let be the weight of the j-th parameter.

[0064] Step 3.2.2.4: Interpolation mesh generation; Using the fused remote sensing image as a spatial framework, an interpolation grid is generated that corresponds one-to-one with the image pixels, ensuring that the interpolation results can be directly superimposed on the remote sensing image. Grid parameter settings: The grid resolution should be consistent with the remote sensing image. Grid coordinate generation: Based on the latitude and longitude of the upper left and lower right corners of the remote sensing image, calculate the center point coordinates (x_f, y_f) of each grid cell as the point to be predicted, where f is the grid cell number, to ensure that all grid cells cover the entire monitoring area without omissions or overlaps.

[0065] Step 3.2.2.5, Ordinary Kriging interpolation calculation: For the center point of each grid cell, the estimated target variable is calculated based on the corrected variogram. The specific steps are as follows: To reduce computational cost and ensure interpolation accuracy, a variable radius neighborhood selection method is used for selecting neighborhood sample points. With the desired center point (x_f, y_f) as the center, the initial search radius is set to the range a of the mutation function; Within the search radius, sample points are selected. If the number of sample points is less than 10, the radius is gradually expanded by 10% of a each time until the number of sample points is greater than or equal to 10. If the number of sample points is greater than 30, the 30 sample points closest to the point to be predicted are retained to avoid interference from edge sample points.

[0066] Based on the variogram, a weight matrix for sample points to prediction points is constructed. This integrates the spatial location information of sample points with auxiliary information on soil physicochemical parameters. The aim is to adapt to the complex heterogeneity of the soil environment within the monitoring area and satisfy two constraints: Unbiasedness constraint: ; Minimum variance constraint: The weights are solved using the Kriging equations, which are in the following form: ; in, Must meet , The variogram value between the i-th and j-th neighboring sample points is calculated using a theoretical model modified with soil physicochemical parameters, reflecting the spatial correlation of the target variable between the two points. The Lagrange multiplier is an auxiliary parameter set to satisfy the unbiasedness constraint of interpolation. The spatial location of the point to be predicted The i-th neighboring sample point and the point to be predicted The variogram values ​​between the sample points and the points to be predicted quantify the spatial correlation between them. The number of neighborhood sample points participating in the interpolation calculation is determined by the variable radius neighborhood selection method: with the point to be predicted as the center, the initial search radius is set to the range 'a' of the variogram, and the neighborhood sample points are selected... One effective neighborhood sample point, The value range is 10-30 to ensure a balance between interpolation accuracy and computational efficiency.

[0067] The weight of each sample point is obtained by solving the system of equations through matrix inversion. .

[0068] Calculate the estimated value of the point to be predicted based on the measured values ​​and weights of the sample points: , Points to be predicted The interpolation results of the target variable, i.e. the estimated values ​​of the corresponding points in the grid cells, are used to generate a continuous spatial distribution map; The measured value of the target variable at the i-th sample point comes from the effective monitoring data of the in-situ sensor.

[0069] Step 3.2.2.6: Spatial distribution map generation. The interpolated grid data is converted into a visual spatial distribution map, which is then overlaid with the fused remote sensing image, as detailed below: Data format conversion: Convert the coordinates of the mesh cells generated in step 3.2.2.4 to the estimated values ​​calculated in step 3.2.2.5. Convert the data to GeoTIFF format to ensure a perfect match of spatial coordinates; Color mapping: Setting the color gradient based on the numerical range of the target variable. Pollutant concentration: A blue-green-yellow-red gradient is used, with blue representing low concentration and red representing high concentration, and color nodes corresponding to risk screening values ​​and control values ​​are marked; Soil physicochemical parameters: A light blue to dark blue gradient was used, with light blue representing low humidity and dark blue representing high humidity; Map enhancement: Add elements such as coordinate system, scale, legend, north arrow, and monitoring area name to generate the final spatial distribution map of pollutant concentration and spatial distribution map of soil physicochemical parameters.

[0070] Feature layer fusion breaks through the point-area data barrier: it resolves the contradiction between the limited measurement range of in-situ sensors and the insufficient quantitative accuracy of remote sensing images. By interpolating, point data is expanded into area data, achieving full-area coverage and pixel-level precise quantification. It considers the influence of soil physicochemical parameters: physicochemical parameters such as soil moisture and organic matter content are included as auxiliary variables in the interpolation, making the results more consistent with the actual distribution patterns of soil pollution. For example, the concentration of heavy metals is higher in soils with high organic matter around industrial areas. It is deeply integrated with data layer fusion: the spatial distribution map generated by interpolation has the same spatial resolution and coordinate system as the fused image, providing a unified data foundation for subsequent fusion calibration (spectral-concentration correlation) and pollutant feature extraction (such as pollution area segmentation), avoiding process breaks caused by differences in data format or accuracy.

[0071] 3.3, Fusion Calibration: Optimization of Spatial and Spectral Accuracy; Spatial calibration: Select 10-20 ground control points (such as road intersections, landmark buildings, and sampling points with known coordinates) in the monitoring area, compare the coordinate deviation of the control points in the fused image with those in the high-precision topographic map (accuracy ≤ 1:10,000), and use Affine transformation to perform spatial calibration on the fused image to ensure that the spatial position error of the fused image is ≤ 0.5 pixels.

[0072] Spectral calibration: Using laboratory chemical analysis data as the true value, select 50-100 sampling points, correlate the spectral values ​​of corresponding pixels in the fused image with the pollutant concentration detected in the laboratory, establish a linear regression model of spectral value-pollutant concentration, and adjust the spectral response coefficient of the fused image to make the determination coefficient of the model ≥0.85, so as to ensure that the spectral information of the fused image can accurately reflect the pollutant concentration, and complete the spectral calibration.

[0073] Step 4: Extract soil pollutant characteristics. Based on the fused and calibrated image data, the spectral characteristics, spatial distribution characteristics, and pollution intensity characteristics of soil pollutants are extracted through spectral analysis and spatial pattern recognition. The specific process is as follows: Step 4.1, Extraction of spectral features of pollutants: Feature band screening: For different types of pollutants (such as heavy metals and organic pollutants), analyze their characteristic absorption peaks in the hyperspectral range. For example, the heavy metal cadmium has a characteristic absorption peak near 2345 nm, and polycyclic aromatic hydrocarbons have characteristic absorption peaks near 1780 nm and 2200 nm. Use the continuous projection algorithm (SPA) or genetic algorithm (GA) to screen out 5-10 characteristic bands that are most sensitive to pollutant concentration from the spectral band of 400-2500 nm, and eliminate interference from irrelevant bands.

[0074] Spectral characteristic parameter calculation: Based on the screened characteristic bands, spectral characteristic parameters are calculated, including the depth, width, and area of ​​characteristic absorption peaks, as well as the first and second derivatives of spectral reflectance. These parameters are used as spectral characteristic indicators of pollutants for subsequent pollutant type identification and concentration inversion.

[0075] 4.2 Extraction of Spatial Distribution Characteristics of Pollutants: Contaminated area segmentation: An object-oriented image segmentation algorithm is adopted, which uses the spectral information (reflectance of characteristic bands), texture information (contrast and correlation extracted from gray-level co-occurrence matrix), and spatial information (pixel adjacency relationship) of the fused image as the segmentation basis. The segmentation scale is set (5-10 for key areas and 20-30 for macro areas), shape factor (0.3-0.5), and compactness factor (0.5-0.7) to segment the fused image into different contaminated units (such as high-contaminated areas, medium-contaminated areas, low-contaminated areas, and uncontaminated areas).

[0076] Spatial pattern recognition: Spatial autocorrelation analysis (Moran's I index) is used to determine the spatial clustering of pollution units: if the Moran's I index is positive and significant (P<0.05), it indicates that the pollutants are clustered, such as pollution clumps around industrial sites; if it is negative and significant, it indicates that the pollutants are discretely distributed, such as non-point source pollution in agricultural areas; at the same time, hotspot analysis is used to identify pollution hotspots (high concentration clusters) and coldspots (low concentration clusters) to clarify the spatial distribution pattern of pollutants.

[0077] 4.3, Pollution Intensity Feature Extraction: Pollution Level Classification: Based on the "Soil Environmental Quality Standard for Agricultural Land Soil Pollution Risk Control (GB15618-2018)" and the "Soil Environmental Quality Standard for Construction Land Soil Pollution Risk Control (GB36600-2018)," and combined with the land use type (agricultural land, construction land) of the monitoring area, the pollutant concentration is divided into four levels: No pollution: concentration ≤ risk screening value; Low pollution: risk screening value < concentration ≤ mild risk control value; Medium pollution: mild risk control value < concentration ≤ severe risk control value; High pollution: concentration > severe risk control value.

[0078] Pollution intensity quantification: Calculate the pollution intensity index for each pollution unit using the following formula: ,in Let j be the concentration of the i-th pollutant in the i-th pollution unit. Let j be the risk screening value for the j-th pollutant. denoted by , represents the weight of the j-th pollutant (determined based on the pollutant's toxicity coefficient; for example, mercury has a higher weight than cadmium, and polychlorinated biphenyls (PCBs) have a higher weight than common pesticides among organic pollutants), and n represents the number of pollutant types. This index quantifies the severity of pollution in different pollution units.

[0079] Step 5: Generate the soil health index. By integrating soil pollutant characteristics, soil physicochemical parameters, and ecological function requirements, a multi-dimensional soil health evaluation system is constructed, and the soil health index is calculated to achieve a quantitative assessment of soil health status.

[0080] Step 6: Conduct soil environmental monitoring. Based on the soil health index and pollutant characteristic extraction results, establish a real-time dynamic monitoring system to achieve dynamic tracking, risk warning, and monitoring report output of soil pollution.

[0081] The description of this invention is given for illustrative and descriptive purposes only and is not intended to be exhaustive or to limit the invention to the forms disclosed. Many modifications and variations will be apparent to those skilled in the art. The embodiments were chosen and described in order to better illustrate the principles and practical application of the invention and to enable those skilled in the art to understand the invention and design various embodiments with various modifications suitable for a particular purpose.

Claims

1. A data analysis method for monitoring the spatial distribution of soil pollution based on image fusion, characterized by: Includes the following steps: Step 1: Collect soil pollution monitoring data; Step 2: Preprocess soil pollution monitoring data; Step 3: Perform image fusion calibration. Through multi-scale and multi-modal image fusion technology, integrate the preprocessed multi-source data, and simultaneously perform spatial and spectral calibration to generate a fused image with both high spatial and high spectral resolution. Step 4: Extract soil pollutant features. Based on the fused and calibrated image data, extract the spectral features, spatial distribution features, and pollution intensity features of soil pollutants. Step 5: Generate the soil health index. By integrating soil pollutant characteristics, soil physicochemical parameters and ecological function requirements, construct a multi-dimensional soil health evaluation system, calculate the soil health index, and achieve a quantitative assessment of soil health status. Step 6: Conduct soil environmental monitoring. Based on the soil health index and pollutant characteristic extraction results, establish a real-time dynamic monitoring system to achieve dynamic tracking, risk warning, and monitoring report output of soil pollution.

2. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 1, characterized in that: Step 1 includes the following steps: Step 1.1 Laboratory chemical analysis data: After collecting soil samples, the concentrations of heavy metals and organic pollutants are detected by instruments, and the latitude, longitude, and sampling depth of each sampling point are recorded. Step 1.2 In-situ sensor detection data: Multi-parameter in-situ sensors are evenly deployed in the monitoring area. The sensors collect pollutant concentration signals and soil physicochemical parameters in real time. The collection frequency is set to once per hour. The data is uploaded to the data platform in real time through the wireless transmission module. At the same time, the sensor location coordinates and collection timestamp are recorded. Step 1.3 Remote sensing and spectral data: Simultaneously acquire multi-source remote sensing and spectral data, including 30m spatial resolution optical remote sensing images from Landsat-8 satellite, 10m spatial resolution multispectral images from Sentinel-2 satellite, hyperspectral data with spatial resolution below 1m acquired by a hyperspectral instrument carried by an UAV, and hyperspectral curves of soil samples collected simultaneously at sampling points by a ground-based portable hyperspectral instrument.

3. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 2, characterized in that: Step 2 includes the following steps: Step 2.1, Preprocessing of laboratory chemical analysis data: Step 2.1.1, outlier removal: For pollutant concentration data from multiple tests of the same sample, outlier judgment is performed, and outlier data caused by experimental operation errors or instrument malfunctions are removed. Step 2.1.2, Data Standardization: Convert the concentration data of different pollutants into a standardized form, using the following formula: ,in Let be the concentration of the j-th pollutant at the i-th sampling point. Let be the average concentration of the j-th pollutant. Let be the standard deviation of the concentration of the j-th pollutant, to eliminate the influence of differences in the order of magnitude of different pollutant concentrations on subsequent analyses; Step 2.2, preprocessing of in-situ sensor detection data; Step 2.2.1, Noise Filtering: For the real-time data collected by the sensor, the moving average filtering method is used to filter high-frequency noise caused by sensor fluctuations and external electromagnetic interference. Step 2.2.2, Data Completion: If data is missing due to transmission interruption or temporary sensor failure, and the missing time is ≤4 hours, linear interpolation is used to complete the short-term missing data; if the missing time exceeds 4 hours, Kriging interpolation is used to complete the data by combining the historical data of the sensor and the data of neighboring sensors to ensure data continuity. Step 2.3, remote sensing and spectral data preprocessing; Step 2.3.1, Remote sensing image preprocessing: For Landsat and Sentinel-2 satellite images, radiometric calibration is performed sequentially to convert the DN values ​​received by the sensor into surface reflectance; atmospheric correction is performed using the FLAASH atmospheric correction model to eliminate the influence of atmospheric scattering and absorption on the image; geometric correction is performed using a quadratic polynomial correction model based on the high-precision topographic map of the monitoring area to convert the image coordinate system to the WGS84 coordinate system, with the correction error controlled within 1 pixel. Step 2.3.2, Hyperspectral data preprocessing: For the hyperspectral images acquired by the UAV, radiometric calibration, atmospheric correction, and geometric correction are performed in sequence. Then, the SIFT feature matching algorithm is used to stitch together multiple images and perform radiometric normalization. Spectral data acquired by ground-based portable hyperspectral analyzers and UAV hyperspectral analyzers were smoothed using Savitzky-Golay filtering; baseline correction was performed using an adaptive iterative reweighted penalized least squares method to eliminate spectral baseline drift; and noise reduction was performed to remove abnormal spectral bands caused by instrument noise and soil particle scattering, retaining the effective spectral range of 400-2500 nm. Step 2.4, Data Format and Coordinate System 1: Convert all preprocessed data sources to GeoTIFF and CSV formats, and set the coordinate system of all spatial data to WGS84 coordinate system. Associate attribute data with corresponding spatial coordinates to ensure the consistency of multi-source data in spatial location.

4. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 2, characterized in that: Step 3 includes the following steps: Step 3.1, the specific process of data layer fusion is as follows: Step 3.1.1: Perform multi-scale wavelet decomposition on the images involved in the fusion, which is divided into first-level decomposition and second-level decomposition. The preprocessed Landsat-8 image is denoted as image A and the Sentinel-2 image is denoted as image B. Perform wavelet decomposition with the same parameters on each image, and use row-column convolution and downsampling to decompose the images into low-frequency approximate components and high-frequency detail components. Step 3.1.2: Reconstruct the high-frequency detail components of the high spatial resolution image and the low-frequency approximation components of the high spectral resolution image to generate a fused image; Step 3.2, Feature Layer Fusion: The remote sensing data is fused with the in-situ sensor data. The specific process is as follows: Step 3.2.1, Spatial matching: Using the fused remote sensing image as a spatial framework, based on the latitude and longitude coordinates of the in-situ sensor, the pollutant concentration and soil physicochemical parameter data collected by the sensor are embedded as feature points into the remote sensing image to establish the spatial relationship between image pixels and sensor data. Step 3.2.2, Interpolation Supplement: Ordinary Kriging interpolation is used. The spatial correlation of pollutant concentration is analyzed by semi-variogram analysis. In-situ sensor data is used as sample points. Combined with the influence of soil physicochemical parameters on the interpolation results, the pollutant concentration and physicochemical parameters in the remote sensing image without sensor deployment are interpolated and estimated to generate continuous spatial distribution maps of pollutant concentration and soil physicochemical parameters. 3.3, Fusion Calibration: Optimization of Spatial and Spectral Accuracy; Spatial calibration: Select 10-20 ground control points within the monitoring area, compare the coordinate deviations of the control points between the fused image and the high-precision topographic map, and use affine transformation to perform spatial calibration on the fused image to ensure that the spatial position error of the fused image is ≤0.5 pixels; Spectral calibration: Using laboratory chemical analysis data as the true value, select 50-100 sampling points, correlate the spectral values ​​of corresponding pixels in the fused image with the pollutant concentration detected in the laboratory, establish a linear regression model of spectral value-pollutant concentration, and adjust the spectral response coefficient of the fused image to make the determination coefficient of the model ≥0.

85.

5. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 4, characterized in that: Step 3.1.1 includes the following steps: Step 3.1.1.1, first-level decomposition, extracting large-scale information; Line direction decomposition: For each row of pixels in image A, perform convolution operations using a low-pass filter and a high-pass filter based on the db4 wavelet basis function: Input pixel sequence: Suppose the pixel value sequence of a certain row of the image. Where n is the length of the pixel sequence, the pixel value sequence X is padded with zeros to extend it, and the extended sequence is: Ensure that the window covers all pixels; Calculation process of low-frequency approximate components: Low-pass filter length coefficient For the expanded sequence Using the length of the low-pass filter as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding low-pass filter coefficients is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: The window moves one pixel at a time until it covers the entire extended sequence. The final convolution result has a length of n+7. Perform interval sampling to obtain the low-frequency approximate component Lrow in the row direction; Calculation process of high-frequency detail components: High-pass filter length coefficient Where k is the coefficient index, for the expanded sequence Using the high-pass filter length as the window size, the convolution is performed by sliding the window from left to right. Within each window, the sum of the pixel value multiplied by the corresponding high-pass filter coefficient is calculated, thus obtaining the convolution result. For the pixel at the starting position t of the window, the convolution result is: For the convolution result Perform interval sampling to obtain the high-frequency detail component Hrow in the row direction; Column-direction decomposition: Repeat convolution and downsampling operations on Lrow and Hrow obtained from row decomposition by column: Convolve and downsample Lrow by column: Obtain a low-frequency approximate component L1; Column-wise convolution and downsampling of Hrow: respectively, yields 1 layer of horizontal high-frequency detail component Hh1, 1 layer of vertical high-frequency detail component Hv1, and 1 layer of diagonal high-frequency detail component Hd1; Results output: After the first layer decomposition of image A, we get {L1_A,Hh1_A,Hv1_A,Hd1_A}; After the first layer decomposition of image B, we get {L1_B,Hh1_B,Hv1_B,Hd1_B}. Step 3.1.1.2, Second-level decomposition: Extracting mesoscale details; Using the low-frequency approximation component L1 obtained from the first-level decomposition as input, the row-column decomposition process is repeated to extract finer details: Perform low-pass / high-pass filtering convolution and downsampling on each row of L1_A to obtain the low-frequency sequence L1_row and the high-frequency sequence H1_row in the row direction; Filter and downsample L1_row and H1_row column by column: Two layers of low-frequency approximation components L2 and two layers of high-frequency detail components Hh2, Hv2, and Hd2 are obtained; Output results: The final result of the 2-layer decomposition of image A is: {L2_A,Hh2_A,Hv2_A,Hd2_A,Hh1_A,Hv1_A,Hd1_A}; The result of the 2-layer decomposition of image B is: {L2_B,Hh2_B,Hv2_B,Hd2_B,Hh1_B,Hv1_B,Hd1_B}.

6. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 5, characterized in that: Step 3.1.2 includes the following steps: Step 3.1.2.1: Complete the replacement of high-frequency detail components according to the decomposition level. The wavelet decomposition level corresponds to the scale of spatial detail: The high-frequency detail components (Hh2, Hv2, Hd2) corresponding to the second layer L2 reflect mesoscale details, while the high-frequency detail components (Hh1, Hv1, Hd1) corresponding to the first layer L1 reflect small-scale details. They need to be replaced in order from the upper layer to the lower layer to ensure that spatial details of different scales are accurately supplemented. The high-frequency detail components of the second layer are replaced to supplement the mesoscale details: the three high-frequency detail components (Hh2_A, Hv2_A, Hd2_A) of the second layer of image A are directly replaced with the corresponding high-frequency detail components (Hh2_B, Hv2_B, Hd2_B) of the second layer of image B. Replacement of high-frequency detail components in the first layer for small-scale detail supplementation: Replace the three high-frequency detail components (Hh1_A, Hv1_A, Hd1_A) in the first layer of image A with the corresponding high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) in the first layer of image B. Low-frequency approximation component preservation: The two low-frequency approximation components (L2_A, L1_A) of image A are completely preserved without replacement, ensuring that the fused image can identify the core of the pollutant's spectral characteristics; After replacement, image A forms a fusion decomposition component set: {L2_A,Hh2_B,Hv2_B,Hd2_B,Hh1_B,Hv1_B,Hd1_B}; Step 3.1.2.2, wavelet inverse transform to generate fused image, is the reverse restoration of the wavelet decomposition sequence, first layer then second layer. It employs row-column convolution and upsampling, specifically divided into two layers of inverse transform to progressively reconstruct the image: Second-level inverse transform: generates mesoscale fused low-frequency approximation components; Using the second-level components (L2_A, Hh2_B, Hv2_B, Hd2_B) in the fusion decomposition component set as input, the operation is performed in the order of column inverse decomposition followed by row inverse decomposition: Inverse column decomposition: Column convolution and upsampling are performed on L2_A, Hh2_B, Hv2_B, and Hd2_B respectively. Using a db4 wavelet basis low-pass or high-pass filter, after convolution operation on each column of each component, the number of columns is doubled by upsampling with zero values ​​inserted at intervals. The specific process is as follows: Step L1: Extract the pixel value sequence of each column from column 1 to column 125 according to the column index order. If the component is a low-frequency approximate component, use the low-pass filter coefficient h for convolution. If the component is a high-frequency detail component, use the high-pass filter coefficient g for convolution. Calculate the column sequence X_col and the filter coefficients in a sliding window manner. The window size is equal to the filter length. The convolution result of each window is the sum of the X_col pixel values ​​in the window × the corresponding filter coefficients. After convolution of each column, a convolution sequence with the same length as the original column is generated. Step L2: After completing the single-column convolution of all input components, the convolution result of the low-frequency approximate component of the same column index is summed with the convolution result of the three high-frequency detail components to obtain the intermediate merged sequence of column inverse decomposition, which integrates the low-frequency contour information and the high-frequency detail information to restore the complete signal of the column in the spatial domain. Step L3: Upsampling is performed on the intermediate merged sequence by inserting zero values ​​at intervals. One zero is inserted between every two adjacent values, doubling the length of the intermediate merged sequence, that is, doubling the column dimension size, matching the spatial resolution of the original image. The upsampled sequence is the final result of the inverse decomposition of the column in the column direction, corresponding to the pixel value sequence of the column in the spatial domain. Step L4: Repeat steps L1-L3 for all columns of the input component. After completing the inverse decomposition calculation of all columns, output the intermediate component with doubled column dimension size. The input component of 125×125 pixels is output as an intermediate component of 250×125 pixels after inverse decomposition in the column direction. Row-direction inverse decomposition: For the result after column inverse decomposition, row-direction repeated convolution and upsampling are performed according to the column-direction inverse decomposition process, doubling the number of rows from 125 to 250, and finally generating the first layer of fused low-frequency approximation component L1_fuse with a size of 250×250 pixels. At this time, L1_fuse has integrated the mesoscale spectral information of image A with the mesoscale spatial details of image B. First inverse transform: Generates the final fused image. Using L1_fuse and the first layer high-frequency detail components (Hh1_B, Hv1_B, Hd1_B) as input, repeat the inverse transform process of the second layer inverse transform: Inverse column decomposition: Perform column convolution and upsampling on L1_fuse, Hh1_B, Hv1_B, and Hd1_B to double the number of columns from 250 to 500; Row-direction inverse decomposition: The column inverse decomposition results are subjected to row convolution and upsampling to double the number of rows from 250 to 500, finally generating a fused image with a size of 500×500 pixels. The spatial resolution is consistent with that of image B, and the spectral resolution is consistent with that of image A.

7. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 4, characterized in that: Step 3.2.1 includes the following steps: Step 3.2.1.1, Coordinate System Verification: First, confirm the consistency of the coordinate systems of the two types of input data. The fused image has been unified into the WGS84 coordinate system. The in-situ sensor data has been recorded with WGS84 latitude and longitude coordinates by GPS positioning during deployment. It is necessary to use spatial analysis software to import the sensor coordinates into the spatial framework of the fused image and verify whether the coordinates of each sensor point fall within the corresponding pixel range. Step 3.2.1.2, Pixel-Sensor Data Association: For each in-situ sensor point, extract the spectral feature information of its corresponding pixel in the fused image, and bind the pollutant concentration data and soil physicochemical parameters of the sensor with the extracted spectral information to form a spatial location-spectral feature-quantitative parameter association dataset. Step 3.2.1.3, outlier screening: Spatial outlier detection is performed on the associated dataset to determine whether it is an outlier. Outliers need to be removed or data needs to be re-collected to avoid affecting the accuracy of subsequent interpolation.

8. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 4, characterized in that: Step 3.2.2 includes the following steps: Step 3.2.2.1, Sample point data screening; Basic data extraction: Extract core data from all valid sensors within the monitoring area from the in-situ sensor database, including: Spatial coordinates: sensor latitude and longitude, unified in the WGS84 coordinate system, consistent with the fused remote sensing imagery; Target variable values: pollutant concentration and soil physicochemical parameters to be interpolated; Auxiliary variable values: corresponding to the original values ​​of soil physicochemical parameters collected by the sensors, used as covariates for interpolation to correct spatial correlation; Step 3.2.2.2, modeling the spatial variogram; Step 3.2.2.3: Correct the variogram by combining soil physicochemical parameters; Step 3.2.2.4: Interpolation grid generation. Using the fused remote sensing image as the spatial frame, an interpolation grid is generated that corresponds one-to-one with the image pixels, ensuring that the interpolation results can be directly superimposed on the remote sensing image. Grid parameter settings: The grid resolution should be consistent with the remote sensing image. Grid coordinate generation: Based on the latitude and longitude of the upper left and lower right corners of the remote sensing image, calculate the center point coordinates (x_f, y_f) of each grid cell as the point to be predicted, where f is the grid cell number, to ensure that all grid cells cover the entire monitoring area without omissions or overlaps; Step 3.2.2.5, Ordinary Kriging interpolation calculation: For the center point of each grid cell, calculate the estimate of its target variable based on the modified variogram. Step 3.2.2.6: Spatial distribution map generation. The interpolated grid data is converted into a visual spatial distribution map and overlaid with the fused remote sensing image.

9. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 8, characterized in that: The spatial variogram modeling process is as follows: The formula for the variogram function is: for the target variable In spatial location and The variogram function is defined as the value at a given position: ; The variogram value with a lag distance of g is used to quantify the spatial variability of the target variable; g is the spatial lag distance, which represents the spatial distance between two sample points. When the lag distance is g, the total number of sample point pairs within the monitoring area that meet this distance condition; the index number of sample point pair i, with a value ranging from 1 to... This is used to iterate through all sample point pairs that satisfy the lag distance g. The spatial location of the i-th sample point; The measured value of the target variable at the i-th sample point; The spatial location of the paired sample point that is at a distance g from the i-th sample point; The measured value of the target variable at the paired sample point with a distance of g from the i-th sample point; Based on the pollutant diffusion characteristics of the monitoring area, a suitable theoretical model is selected to fit the experimental variogram: If the pollutants diffuse uniformly, choose the spherical model: When the hysteresis distance g is less than or equal to the range a: ; When the hysteresis distance g is greater than the range a: ; If pollutants are blocked by terrain, choose the index model: ; in, denoted as nugget value, reflecting random error; C is the sill value, reflecting spatial structural variation; a is the range, the maximum distance of spatial correlation, beyond which there is no spatial correlation between sample points.

10. The data analysis method for monitoring the spatial distribution of soil pollution based on image fusion as described in claim 8, characterized in that: The specific process of step 3.2.2.5 is as follows: To reduce computational cost and ensure interpolation accuracy, a variable radius neighborhood selection method is used for selecting neighborhood sample points. With the desired center point (x_f, y_f) as the center, the initial search radius is set to the range a of the mutation function; Within the search radius, sample points are selected. If the number of sample points is less than 10, the radius is gradually expanded by 10% of a each time until the number of sample points is greater than or equal to 10. If the number of sample points is greater than 30, the 30 sample points closest to the point to be predicted are retained to avoid interference from edge sample points. Based on the variogram, a weight matrix for sample points to prediction points is constructed. This integrates the spatial location information of sample points with auxiliary information on soil physicochemical parameters. The aim is to adapt to the complex heterogeneity of the soil environment within the monitoring area and satisfy two constraints: Unbiasedness constraint: ; Minimum variance constraint: The weights are solved using the Kriging equations, which are in the following form: ； in, Must meet , The variogram value between the i-th and j-th neighboring sample points is calculated using a theoretical model modified with soil physicochemical parameters, reflecting the spatial correlation of the target variable between the two points. The Lagrange multiplier is an auxiliary parameter set to satisfy the unbiasedness constraint of interpolation. The spatial location of the point to be predicted The i-th neighboring sample point and the point to be predicted The variogram values ​​between the sample points and the points to be predicted quantify the spatial correlation between them. The number of neighborhood sample points participating in the interpolation calculation is determined by the variable radius neighborhood selection method: with the point to be predicted as the center, the initial search radius is set to the range 'a' of the variogram, and the neighborhood sample points are selected... One effective neighborhood sample point, The value range is 10-30 to ensure a balance between interpolation accuracy and computational efficiency. The weight of each sample point is obtained by solving the system of equations through matrix inversion. ; Calculate the estimated value of the point to be predicted based on the measured values ​​and weights of the sample points: , Points to be predicted The interpolation results of the target variable, i.e. the estimated values ​​of the corresponding points in the grid cells, are used to generate a continuous spatial distribution map; The measured value of the target variable at the i-th sample point comes from the effective monitoring data of the in-situ sensor.

Images