An ecological environment detection method and system based on multispectral remote sensing fusion

Abstract

The application discloses an ecological environment detection method and system based on multispectral remote sensing fusion, and relates to remote sensing image processing. Multispectral remote sensing image data of a target region is collected. Multiband joint atmospheric correction processing is performed on the multispectral remote sensing image according to the scattering coefficients, atmospheric light values and transmittances of various bands. Foreground and background analysis is performed on the corrected multispectral remote sensing image. The vegetation regions and non-vegetation regions of various band images are fused by using different weight strategies respectively, and multispectral fusion images are generated. The spectral features of the multispectral fusion images are extracted, a standard vegetation spectral feature library is constructed, and the regions deviating from the standard vegetation spectrum are identified by using an anomaly detection algorithm according to the extracted spectral features, so that various vegetation cover rates are obtained. In view of the fact that direct foreground and background division of the multispectral remote sensing image under atmospheric interference conditions leads to low vegetation recognition accuracy, the application performs vegetation division and the like after obtaining a clear image, so that the detection accuracy is improved.

Description

Technical Field

[0001] This application relates to the field of remote sensing image processing, and in particular to an ecological environment monitoring method and system based on multispectral remote sensing fusion. Background Technology

[0002] Multispectral remote sensing technology acquires the reflectance information of ground objects in different electromagnetic bands, enabling effective identification of vegetation types and assessment of vegetation health, thereby achieving quantitative monitoring of ecological environment quality. It demonstrates unique advantages, particularly in large-scale areas such as vegetation cover statistics, vegetation growth assessment, and ecosystem health diagnosis.

[0003] However, in practical applications, the quality of multispectral remote sensing images is often severely affected by atmospheric conditions. Atmospheric scattering and absorption effects such as haze and aerosols can lead to image blurring and reduced contrast, seriously affecting subsequent image analysis and information extraction. Existing multispectral remote sensing image processing methods mainly suffer from the following technical problems:

[0004] First, traditional methods typically employ a separation-then-dehazing process, directly analyzing the foreground and background on blurred images affected by atmospheric interference to attempt to identify vegetated and non-vegetated areas. Due to image degradation caused by atmospheric scattering, the spectral differences between vegetated and non-vegetated areas are weakened, resulting in blurred boundaries. This significantly reduces the accuracy of separation methods based on traditional indicators such as the Normalized Difference Vegetation Index (NDVI). Particularly under hazy weather conditions, vegetation edge information is severely lost, easily leading to misidentification and omission of vegetated areas.

[0005] Secondly, in terms of atmospheric correction, existing methods mostly employ a strategy of processing each band independently, i.e., performing dehazing on each band image separately. This method ignores an important characteristic of multispectral data—the inherent correlation between different bands based on physical laws. According to atmospheric scattering theory, the scattering coefficients of electromagnetic waves of different wavelengths in the atmosphere follow specific physical relationships (such as Angstrom's law), and there is also a definite functional relationship between the atmospheric light value and transmittance of each band. Processing each band independently not only fails to utilize these physical constraints to improve correction accuracy but may also lead to spectral inconsistencies in the corrected multi-band images, affecting subsequent spectral feature analysis.

[0006] For example, relevant technical document CN119418226A discloses an ecological environment detection method and system based on multispectral remote sensing fusion technology. This method includes: acquiring multispectral remote sensing images of the target area; determining the foreground information rate at each pixel level in the remote sensing image of each band; calculating the distortion coefficient of the remote sensing image of each band; obtaining the relative loss coefficient at each pixel level in the remote sensing image of each band; dividing the remote sensing image of each band into superpixel blocks and determining the band independence of each superpixel block in the remote sensing image of each band; screening candidate pixels within each superpixel block, determining atmospheric light values, and performing dehazing on the remote sensing image of each band using the dark channel prior method; and performing ecological environment detection based on the dehazed remote sensing images of all bands. However, this application directly calculates the foreground information rate on blurred images affected by atmospheric interference, leading to inaccurate foreground-background separation; simultaneously, it uses the dark channel prior method to independently dehaze each band, ignoring the inherent connection between different bands based on the physical laws of atmospheric scattering, and cannot guarantee the spectral consistency of the multi-band images after dehazing. Furthermore, the band independence processing approach of this application contradicts the band correlation characteristics of multispectral data, fails to fully utilize the advantages of multispectral remote sensing, and does not consider the spectral response differences between vegetated and non-vegetated areas in the subsequent fusion process, resulting in low vegetation identification accuracy and ecological environment detection reliability under complex atmospheric conditions. Summary of the Invention

[0007] To address the issue of low vegetation identification accuracy caused by directly segmenting foreground and background in multispectral remote sensing images under atmospheric interference conditions, this application provides an ecological environment detection method and system based on multispectral remote sensing fusion. Atmospheric correction is performed through multi-band joint optimization, and vegetation segmentation is performed after obtaining a clear image, thereby improving detection accuracy.

[0008] One aspect of this application provides an ecological environment detection method based on multispectral remote sensing fusion, comprising: S1, acquiring multispectral remote sensing image data of a target area to obtain an original image set containing n bands and their corresponding center wavelengths; S2, calculating the scattering coefficient, atmospheric light value, and transmittance between the bands respectively; performing multi-band joint atmospheric correction processing on the multispectral remote sensing images according to the scattering coefficient, atmospheric light value, and transmittance of each band; S3, performing foreground and background analysis on the corrected multispectral remote sensing images to identify vegetated and non-vegetated areas; S4, fusing the vegetated and non-vegetated areas of each band image using different weighting strategies to generate a multispectral fused image; S5, extracting the spectral features of the multispectral fused image, constructing a standard vegetation spectral feature library, and identifying areas deviating from the standard vegetation spectrum using an anomaly detection algorithm based on the extracted spectral features to obtain the coverage rates of various vegetation types.

[0009] Multispectral remote sensing image data refers to surface reflection or radiation information acquired simultaneously by remote sensing sensors across multiple different electromagnetic wavelength ranges. In this application, the multispectral remote sensing image data includes a set of original images in n bands and their corresponding center wavelengths. Each band corresponds to a specific electromagnetic spectrum range, such as visible light (red, green, blue), near-infrared, and short-wave infrared bands.

[0010] The scattering coefficient is a physical parameter describing the intensity of light scattering by particles in the atmosphere. In this application, the scattering coefficient β(λ) is related to the wavelength λ and follows the relationship... Where α is the Angstrom index, determined according to weather conditions (1.3 for sunny days and 0.8 for hazy days). The scattering coefficient is used to quantify the degree of atmospheric scattering influence on different wavebands.

[0011] Atmospheric light value (A): refers to the contribution of atmospheric scattering to light at infinity, representing the brightness value of a pixel when it is completely blocked by the atmosphere. In this scheme, it is estimated by averaging the top 0.1% of the largest pixel values ​​in a reference band image.

[0012] Transmittance (t): Represents the proportion of light that is not scattered or absorbed after passing through the atmosphere, with a value ranging from 0 to 1. Higher transmittance indicates less atmospheric influence and a clearer image. The calculation formula is as follows: .

[0013] Atmospheric correction is a preprocessing step to eliminate the influence of the atmosphere on remote sensing images. This application employs multi-band joint atmospheric correction. By constructing an optimized objective function E, and considering data fidelity terms, band consistency constraints, and spatial smoothness constraints, the optimal transmittance and atmospheric light values ​​are solved using the alternating direction multiplier method, ultimately restoring clear images for each band. .

[0014] Foreground / background analysis is an image analysis technique that segments remotely sensed images into regions of interest (foreground) and other regions (background). In this application, foreground / background analysis specifically refers to the process of identifying vegetated and non-vegetated regions. This is achieved by calculating the Normalized Difference Vegetation Index (NDVI), determining the optimal segmentation threshold using the Otsu's method, and generating a vegetation mask by combining morphological processing.

[0015] Standard vegetation spectra, typical spectral features extracted from healthy vegetation samples, serve as a reference standard for judging vegetation health status. This application selects the top P% pixels with the highest normalized vegetation index as healthy vegetation samples, and uses K-means clustering to obtain k cluster centers for healthy vegetation types. The mean vector of each cluster... Covariance Matrix A standard vegetation spectral feature library is constructed for the detection of abnormal vegetation and the statistics of vegetation coverage.

[0016] Furthermore, S2, calculate the scattering coefficient, atmospheric light value, and transmittance between wavebands, including: establishing the relationship between the scattering coefficient and wavelength: ,in, Let λ be the scattering coefficient at wavelength λ. Reference wavelength The scattering coefficient at the location, where α is the Angstrom exponent; α is taken as a value according to the weather conditions of the image, 1.3 for sunny days and 0.8 for hazy days;

[0017] Based on the relationship, calculate the scattering coefficient ratio between each band, and generate the scattering coefficient ratio matrix R, including: based on the center wavelengths of the n bands. Initialize an n×n scattering coefficient ratio matrix R; iterate through each position (i, j) of matrix R, where i is the row index and j is the column index, i, j ∈ [1, n]; i and j are positive integers; for each position (i, j). , to perform calculations, where, The center wavelength of the i-th band is Let λ be the center wavelength of the j-th band, and α be the Angstrom exponent.

[0018] Select the band with the longest wavelength as the reference band, and then use the reference band image. The image is divided into m x m blocks. The maximum pixel value of each block is calculated, and the average of the top 0.1% of the maximum pixel values ​​of all blocks is selected as the atmospheric light value. According to atmospheric light values Calculate the transmittance of the reference band. , where ω is an adjustment factor with a value of 0.95; This represents the pixel value of the reference band at coordinates (x, y); based on the transmittance of the reference band... Calculate the atmospheric light values ​​of other bands i using the scattering coefficient ratio matrix R. , transmittance Where ref represents the reference band, Indicates the atmospheric light value of the reference band;

[0019] The Angstrom exponent (α) is an important parameter describing the size distribution characteristics of atmospheric aerosol particles, characterizing the sensitivity of the scattering coefficient to wavelength variations. In this application, the Angstrom exponent is used in the relationship between the scattering coefficient and wavelength. The Angstrom index is an exponential term. This parameter is determined based on the weather conditions at the time the image was captured: under clear weather conditions, α is 1.3, indicating that small particles dominate the atmosphere; under hazy weather conditions, α is 0.8, indicating that the proportion of large particles in the atmosphere increases. The larger the Angstrom index, the more drastic the change in the scattering coefficient with wavelength.

[0020] The scattering coefficient ratio matrix R is an n×n matrix used to describe the relative relationships of scattering coefficients between different wavebands. Each element R(i,j) in the matrix represents the ratio of the scattering coefficient of the i-th waveband to that of the j-th waveband. When i=j, R(i,j)=1 (the ratio for the same waveband is 1); when i≠j, according to the formula... Calculate, where, and These represent the center wavelengths of the i-th and j-th bands, respectively. This matrix establishes a quantitative relationship between the scattering characteristics of different bands, providing a basis for subsequent calculations of atmospheric light values ​​and transmittance in other bands.

[0021] Reference band image The remote sensing image selected is the one with the longest wavelength as the reference. The longest wavelength band is chosen because longer wavelengths are less affected by atmospheric scattering, providing a more stable basis for atmospheric parameter estimation. In this application, the reference band image is divided into m×m image blocks. The atmospheric light value is estimated by calculating the maximum pixel value of each image block and selecting the average of the top 0.1% of the maximum pixel values ​​from all blocks. Transmittance of the reference band After the calculation is completed, the atmospheric light value and transmittance of other bands can be derived using the scattering coefficient ratio matrix R, realizing the transfer of atmospheric parameters between multiple bands.

[0022] In particular, traditional methods estimate atmospheric light values ​​and transmittance independently for each band, which easily leads to a lack of physical correlation between parameters of different bands, and even contradictions. This application uses a formula... and This ensures that atmospheric parameters across all bands follow the same physical model, maintaining inherent consistency among the parameters. This consistency is crucial for subsequent multispectral fusion and spectral feature analysis.

[0023] Furthermore, choosing the longest wavelength band as the reference band is physically reasonable because long-wavelength bands are less affected by atmospheric scattering, resulting in more stable and reliable atmospheric light value estimates. By transferring parameters from the stable reference band to other bands through the scattering coefficient ratio matrix, the potential error amplification problem caused by directly estimating parameters in the heavily scattered short-wavelength band is avoided. Simultaneously, the Angstrom exponent α is adaptively adjusted according to weather conditions (1.3 for sunny days, 0.8 for hazy days), further improving the model's adaptability to different atmospheric conditions.

[0024] Finally, compared to performing complex atmospheric parameter estimation independently for each band, this method only requires a complete parameter estimation (including image segmentation, statistical analysis, etc.) on the reference band, and the parameters of other bands can be obtained through simple matrix operations. This design greatly reduces computational complexity, especially for hyperspectral data containing multiple bands, where the improvement in computational efficiency is even more significant.

[0025] Furthermore, multi-band joint atmospheric correction processing is performed on the multispectral remote sensing images based on the scattering coefficient, atmospheric light value, and transmittance of each band. This includes: constructing a multi-band joint optimization objective function based on the scattering coefficient, atmospheric light value, and transmittance of each band.

[0026] ;in, Let be the observed pixel value of band i at coordinates (x, y). Let be the sharp pixel value to be recovered for band i at coordinates (x, y). Let be the transmittance of band i at coordinates (x, y). The atmospheric light value for band i. This represents the band consistency constraint weight, with a value range of [0.1, 1.0]. The weights for spatial smoothing constraints range from [0.01, 0.1]. and Transmittance Gradient operators in the x and y directions, , n represents the total number of bands;

[0027] The objective function E is solved using the alternating direction multiplier method to obtain the optimized transmittance for each band. and atmospheric light value ;

[0028] Based on the optimized transmittance and atmospheric light value Through formula Clear images were obtained after correction for each band. .

[0029] Specifically, in the objective function of this application, the first data fidelity term ensures that the corrected image conforms to the atmospheric scattering physical model; the second band consistency constraint introduces Angstrom's law into the optimization process, forcing the transmittance of different bands to meet physical laws; and the third spatial smoothness constraint ensures the spatial continuity of the transmittance map.

[0030] The second term of the objective function Angstrom's law of atmospheric scattering Transforming the constraints into linear constraints in the logarithmic domain ensures that the transmittance of each band always satisfies the physical correlation during the optimization process, avoiding inconsistencies between bands that may arise from independent processing.

[0031] Furthermore, traditional methods solve for each band independently, which easily leads to local optima, and the local optima for each band may contradict each other. This invention, through joint optimization, solves for the parameters of all bands simultaneously within the same objective function, ensuring the global optimality of the solution. The use of the Alternating Direction Multiplier Method (ADMM) further guarantees convergence and computational efficiency when dealing with this large-scale optimization problem.

[0032] Furthermore, S3, performs foreground and background analysis on the corrected multispectral remote sensing image to identify vegetated and non-vegetated areas, including: based on the clear images corrected for each band. Select red light band image and near-infrared band images Calculate the normalized vegetation index : ;right Histogram analysis was performed, and the vegetation segmentation threshold was automatically determined using the Otsu's method. Based on vegetation segmentation threshold Generate initial vegetation cover :when hour, , indicating a vegetated area; when hour, , indicating a non-vegetated area;

[0033] Initial vegetation cover Morphological processing is performed, including first performing an opening operation with a radius of 3 pixels to remove noise points, and then performing a closing operation with a radius of 5 pixels to fill the voids inside the vegetation, resulting in an optimized vegetation mask. ;

[0034] Based on the optimized vegetation cover Corrected images for each band Image segmented into vegetation regions Images of non-vegetated areas : ; .

[0035] In particular, atmospheric scattering severely weakens the spectral differences between vegetated and non-vegetated areas in the original multispectral image. The scattering effect caused by haze homogenizes the strong reflectivity of vegetation in the near-infrared band and its strong absorption in the red band, resulting in a compressed dynamic range of NDVI values ​​and a significant reduction in the distinguishability between vegetated and non-vegetated areas. This invention provides a clearer image after correction. The NDVI calculation restores the true differences in vegetation spectral characteristics, enabling the NDVI value to accurately reflect the actual condition of the vegetation and significantly improving the accuracy of vegetation identification.

[0036] Furthermore, atmospheric scattering has a spatial low-pass filtering effect, blurring edge information in images and causing vegetation boundaries to become indistinct. Traditional methods for segmenting blurred images are prone to boundary localization errors, resulting in over-segmentation or under-segmentation of vegetation areas. This invention analyzes clear images after dehazing, restoring edge information and enabling the Otsu's method to accurately find the optimal segmentation threshold between vegetation and non-vegetation areas. This enables precise boundary positioning.

[0037] Finally, the Otsu's method relies on the bimodal characteristics of the histogram to determine the optimal threshold. In hazy images, due to reduced contrast, Histograms often exhibit a single peak or indistinct peaks, leading to the failure of automatic threshold selection or significant bias. The corrected, clear image restores vegetation and non-vegetation areas. The normal distribution of values ​​results in a clear bimodal structure in the histogram, ensuring the stability and reliability of automatic threshold selection.

[0038] Furthermore, in S4, different weighting strategies are applied to fused the vegetated and non-vegetated areas of each band image to generate a multispectral fused image. This includes: for the vegetated area image, calculating the fusion weight based on the sensitivity of each band to the spectral response of vegetation, where the near-infrared band has a weight of 0.4, the infrared band has a weight of 0.2, the green band has a weight of 0.2, and for the remaining (n-3) bands, each band is assigned a weight of 0.2 / (n-3); for the non-vegetated area image, based on the optimized transmittance... Calculate the fusion weights for each band: Based on the fusion weights, the images of vegetated areas and non-vegetated areas are weighted and fused separately to obtain a multispectral fused image.

[0039] In particular, vegetation possesses unique and stable spectral response characteristics: it exhibits strong reflectivity in the near-infrared band (700-1100nm) and strong absorption in the red band (600-700nm). This red-edge feature is a key marker distinguishing vegetation from other ground features. This invention employs a fixed-weight strategy for vegetation areas, assigning the highest weight of 0.4 to the near-infrared band, thus fully highlighting the spectral characteristics of vegetation. This design ensures maximum preservation of vegetation information in the fused image, which is beneficial for subsequent vegetation health status analysis and anomaly detection.

[0040] Non-vegetated areas contain various land cover types (such as soil, water bodies, and buildings), and their spectral characteristics are complex and variable, with no single optimal band. This invention innovatively utilizes transmittance obtained during atmospheric correction. As an indicator of image quality, higher transmittance indicates less atmospheric influence on that wavelength band, resulting in better image quality. This is achieved through the formula... By calculating the fusion weights, adaptive fusion based on image quality is achieved, enabling the selection of the clearest band information for fusion at each pixel location.

[0041] Furthermore, traditional fusion methods require additional calculation of image quality evaluation metrics, increasing computational complexity. This invention cleverly reuses the transmittance map already optimized during atmospheric correction. This not only avoids redundant calculations, but more importantly, transmittance itself is a direct physical representation of image sharpness.

[0042] Finally, for vegetated areas, the fixed weight strategy ensures that all vegetation pixels fully reflect their spectral characteristics, avoiding the weakening of vegetation features due to differences in local atmospheric conditions. For non-vegetated areas, the adaptive weight strategy ensures that the most reliable band information is selected under different atmospheric conditions. This differentiated processing preserves the essential characteristics of ground features while adapting to spatial variations in atmospheric conditions.

[0043] Furthermore, in step S5, spectral features are extracted from the multispectral fusion image to construct a standard vegetation spectral feature library. Based on the extracted spectral features, an anomaly detection algorithm is used to identify regions deviating from the standard vegetation spectrum, resulting in various vegetation coverage rates. This includes: extracting vegetation region pixels from the multispectral fusion image based on the morphologically processed vegetation mask, and combining this with the clear images corrected for each band. Construct a multidimensional spectral feature vector for each vegetation pixel. The multidimensional spectral feature vector includes the pixel values ​​of each band. and the ratio of pixel values ​​in adjacent bands ;

[0044] From the pixels of the vegetation region, the top P% pixels with the highest normalized vegetation index are selected as healthy vegetation samples. K-means clustering is performed on the multidimensional spectral feature vectors of the healthy vegetation samples to obtain k cluster centers of healthy vegetation types. The mean vector of each cluster is calculated. Covariance Matrix , and store them in the feature library as the spectral features of the i-th type of standard vegetation;

[0045] For all vegetation region pixels in the multispectral fusion image, calculate the distance between the corresponding multidimensional spectral feature vector and the spectral features of various standard vegetation types in the feature library. ; Where v(x, y) is the multidimensional spectral feature vector of the vegetation pixel at coordinates (x, y). Let be the mean vector of the i-th type of standard vegetation in the feature library. Let be the covariance matrix of the i-th type of standard vegetation. is the inverse of the covariance matrix, and the superscript T denotes the vector transpose.

[0046] For each vegetation pixel, select the minimum distance. and the corresponding categories, Less than Pixels with good vegetation type are classified as corresponding healthy vegetation types, while those with poor vegetation type are marked as abnormal vegetation areas. Let be the mean standard deviation of each dimension of the feature vector of the i-th type of healthy vegetation;

[0047] The ratio of the number of pixels of each type of healthy vegetation to the number of pixels of each type of abnormal vegetation is calculated to obtain the coverage rate of each type of healthy vegetation and the coverage rate of each type of abnormal vegetation.

[0048] In particular, the multidimensional spectral feature vector contains not only the absolute pixel values ​​of each band. It also includes the ratio of pixel values ​​in adjacent bands. The absolute value reflects the reflectance intensity of vegetation in each band, while the band ratio characterizes the slope change of the spectral curve, especially capturing the red-edge effect unique to vegetation. Extracting these features from clear images avoids the distortion of spectral features by atmospheric scattering, ensuring the authenticity and reliability of the features.

[0049] Furthermore, methods typically use predefined vegetation index thresholds for health assessment, lacking adaptability to specific scenarios. This invention selects the top P% pixels with the highest NDVI as healthy vegetation samples, automatically learning the spectral characteristics of healthy vegetation from the data. K-means clustering further identifies k different healthy vegetation types, each with its specific spectral feature distribution (mean). Covariance This data-driven approach can adapt to spectral differences in different regions, seasons, and vegetation types, making it more universally applicable.

[0050] Finally, this application uses Mahalanobis distance. Unlike Euclidean distance, Mahalanobis distance considers the correlation between different dimensions of the feature (via the covariance matrix Σᵢ) and normalizes the scale differences between different dimensions. This is particularly important for multispectral data, as the reflectance ranges of different bands can vary greatly, and there are correlations between bands. Mahalanobis distance can more accurately measure the degree of anomalies in spectral features.

[0051] Another aspect of this application provides an ecological environment monitoring system based on multispectral remote sensing fusion, comprising: a data acquisition module for acquiring multispectral remote sensing image data of a target area, obtaining an original image set containing n bands and the corresponding center wavelengths; and an atmospheric correction module, comprising: a scattering coefficient calculation unit for establishing the relationship between the scattering coefficient and the wavelength. The Angstrom index α is determined based on the weather conditions in the image, and the scattering coefficient ratios between the n bands are calculated based on the center wavelengths of each band, generating an n×n scattering coefficient ratio matrix R, where the matrix elements... The atmospheric parameter estimation unit is used to select the longest wavelength band as the reference band, divide the reference band image into m×m image blocks, and extract the atmospheric light value Aᵣ. e f, calculate the transmittance of the reference band. The atmospheric light values ​​for other bands were calculated based on the scattering coefficient ratio matrix R. and transmittance The joint optimization unit is used to construct a multi-band joint optimization objective function E that includes data fidelity terms, band consistency constraints, and spatial smoothness constraints. The optimized transmittance of each band is obtained by solving the alternating direction multiplier method. and atmospheric light value And based on the optimized parameters, through the formula Restore clear images in each band;

[0052] The vegetation recognition module is used to calculate the normalized vegetation index based on the corrected red light band and near-infrared band images, determine the vegetation segmentation threshold by the maximum inter-class variance method, generate an initial vegetation mask and perform morphological processing to obtain an optimized vegetation mask, and segment each band image into vegetation region and non-vegetation region.

[0053] The image fusion module is used to fuse vegetated areas using fixed weights of 0.4 for the near-infrared band, 0.2 for the infrared band, 0.2 for the green band, and 0.2 / (n-3) for the remaining bands. For non-vegetated areas, the module uses the optimized transmittance... Calculate adaptive weights The images are fused to generate a multispectral fused image;

[0054] The anomaly detection and statistics module includes: a feature extraction unit, used to extract vegetation pixels from clear images of each band based on vegetation masks, and construct a feature extraction module containing pixel values ​​of each band. Ratio with adjacent bands The multidimensional spectral feature vector; the feature library construction unit, used to select the top P% pixels with the highest normalized vegetation index as healthy vegetation samples, obtain k healthy vegetation types through K-means clustering, and calculate the mean vector of each type. Covariance Matrix As a standard vegetation spectral feature; the anomaly detection unit is used to calculate the Mahalanobis distance between the feature vector v(x, y) of each vegetation pixel and various standard features in the feature library. ,according to The threshold determines whether a pixel belongs to healthy or abnormal vegetation; the coverage statistics unit is used to count the number of pixels of various types of healthy and abnormal vegetation and calculate the corresponding vegetation coverage.

[0055] Compared to existing technologies, the advantages of this application are:

[0056] In existing technologies, direct foreground and background analysis of multispectral remote sensing images under atmospheric interference conditions and independent atmospheric correction of each band fail to fully utilize the physical correlation between bands, resulting in low accuracy in vegetation identification.

[0057] This application establishes a physical relationship model of inter-band scattering coefficients, employs multi-band joint optimization for atmospheric correction, and then performs precise separation of vegetated and non-vegetated areas after obtaining clear images. Differentiated fusion strategies are adopted for different areas, and finally, anomaly detection algorithms are used to achieve accurate monitoring of the ecological environment. This improves the identification accuracy of vegetated areas under atmospheric interference conditions and makes full use of the band correlation of multispectral data to detect the status of healthy vegetation and accurately count the coverage of various types of vegetation. Attached Figure Description

[0058] This application will be further described by way of exemplary embodiments, which will be described in detail with reference to the accompanying drawings. These embodiments are not limiting; in these embodiments, the same reference numerals denote the same structures, wherein:

[0059] Figure 1 This is an exemplary flowchart of an ecological environment monitoring method based on multispectral remote sensing fusion, according to some embodiments of this application;

[0060] Figure 2 This is an exemplary flowchart illustrating the segmentation of remote sensing images according to some embodiments of this application;

[0061] Figure 3 This is an exemplary flowchart illustrating the acquisition of vegetation coverage according to some embodiments of this application. Detailed Implementation

[0062] The methods and systems provided in the embodiments of this application will now be described in detail with reference to the accompanying drawings.

[0063] like Figure 1 As shown, S1 involves acquiring multispectral remote sensing image data of the target area to obtain a raw image set containing n bands. and the corresponding center wavelength ;

[0064] S2, calculate the scattering coefficient, atmospheric light value, and transmittance between wavebands, including:

[0065] Based on atmospheric scattering theory, a power function model of the scattering coefficient and wavelength is established. ,in This represents the scattering coefficient at wavelength λ. Reference wavelength The scattering coefficient at the location is given by α, which is the Angstrom exponent. The value of α is determined based on the weather conditions at the time of image acquisition: 1.3 for sunny conditions and 0.8 for hazy conditions.

[0066] Based on the above relationship, an n×n dimensional scattering coefficient ratio matrix R is constructed to characterize the proportional relationship of scattering coefficients between different wavebands. For each element in the matrix... When i equals j, it is set to 1, representing the ratio within the same band; when i is not equal to j, it is determined by the formula... Calculate, where, and These are the center wavelengths of the i-th and j-th bands, respectively. This matrix establishes the physical basis for the transfer of atmospheric parameters between bands.

[0067] The longest wavelength band is chosen as the reference band because it is least affected by atmospheric scattering, resulting in more stable parameter estimation. The reference band image is then used. The image is divided into m×m blocks. For each block, the maximum pixel value is extracted; these maximum values ​​typically correspond to the region with the strongest atmospheric illumination. All maximum values ​​from all blocks are sorted in descending order, and the top 0.1% of pixel values ​​are selected and their average value is calculated as the atmospheric illumination value for the reference band. .

[0068] Based on the estimated atmospheric light values, the transmittance distribution of the reference band is calculated using the dark channel prior principle. ω=0.95 is an adjustment factor used to retain a small amount of haze to enhance the realism of the image.

[0069] Atmospheric parameters are transferred from the reference band to other bands using the scattering coefficient ratio matrix R. For the i-th band, its atmospheric light value is obtained through... The transmittance was calculated to be obtained through... Calculations show that This represents the ratio of the i-th band to the reference band in the scattering coefficient ratio matrix. This parameter transfer based on a physical model ensures the consistency of parameters across all bands.

[0070] Construct a multi-band joint optimization objective function:

[0071] ;

[0072] in, Let be the observed pixel value of band i at coordinates (x, y). Let be the sharp pixel value to be recovered for band i at coordinates (x, y). Let be the transmittance of band i at coordinates (x, y). The atmospheric light value for band i. This represents the band consistency constraint weight, with a value range of [0.1, 1.0]. The weights for spatial smoothing constraints range from [0.01, 0.1]. and Transmittance Gradient operators in the x and y directions, , n represents the total number of bands;

[0073] The objective function E is solved using the alternating direction multiplier method to obtain the optimized transmittance for each band. and atmospheric light value :

[0074] Transmittance obtained from physical model and atmospheric light value As the initial value for optimization, the iteration termination condition is set to a relative error of less than 10⁻. 4 Or reach the maximum number of iterations of 100. Initialize a clear image. For observation images .

[0075] The first step is to fix the transmittance and atmospheric light values ​​and update the image to a clearer state. For each pixel in each band, the least squares method is used to solve for the data fidelity term of the atmospheric scattering model. Specifically, this refers to the objective function with respect to... Taking the partial derivative and setting it to zero, we get .

[0076] The second step involves fixing the clear image and atmospheric light values, and then updating the transmittance. This step requires considering the influence of three constraints simultaneously. An optimization subproblem regarding transmittance is constructed and solved using gradient descent. After each gradient update, the transmittance value is constrained within the range of [0.1, 1.0] to ensure physical plausibility. The update process requires calculating: the gradient of the data fidelity term with respect to transmittance; the gradient of the band consistency constraint term in the logarithmic domain; and the Laplace operator for the spatial smoothing term.

[0077] The third step involves fixing the clear image and transmittance, and then updating the atmospheric light value. For each band, the optimal atmospheric light value is solved by minimizing the data fidelity term and band consistency constraints. Considering that the atmospheric light value is a global parameter, a weighted average is applied to the contributions of all pixels, with the weights determined by... Decision. Simultaneously, through the proportional relationships between bands. Constraints are applied to adjust atmospheric light values.

[0078] The fourth step is to check the convergence condition. Calculate the relative rate of change of transmittance and atmospheric light value between the current iteration and the previous iteration. If the rate of change of parameters in all bands is less than the set threshold, the algorithm is considered to have converged.

[0079] After each iteration, a consistency correction is performed on the transmittance of all bands. The deviation between the transmittance ratio of each pair of bands (i, j) and the theoretical value is calculated, and a weighted adjustment is made to bring it closer to the prediction of the physical model. The weights are dynamically determined based on the image quality of each band.

[0080] Edge-preserving filtering is applied to the converged transmittance map using a guided filter that uses the original image as a guide to enhance spatial continuity while preserving edge details. The filter radius is adaptively set according to the image resolution.

[0081] After iterative optimization and post-processing, the optimized transmittance for each band was obtained. and atmospheric light value The entire solution process involves alternately optimizing different variables to gradually reduce the objective function value, ultimately obtaining the optimal solution that satisfies the atmospheric scattering physical model while maintaining band consistency and spatial smoothness.

[0082] Based on the optimized atmospheric light value and transmittance Clear images of each band are recovered through inverse operations of the atmospheric scattering model. The formula is applied. Image restoration is performed, with 0.1 in the denominator serving as a protection threshold to prevent division by zero. This process removes the effects of atmospheric scattering and absorption, resulting in clear images corrected for each band. .

[0083] like Figure 2 As shown, S3 performs foreground and background analysis on the corrected multispectral remote sensing image to identify vegetated and non-vegetated areas, including:

[0084] Based on the clear images after correction for each band Select red light band image and near-infrared band images Calculate the normalized vegetation index : ];

[0085] right Histogram analysis was performed, and the vegetation segmentation threshold was automatically determined using the Otsu's method. ;

[0086] Based on vegetation segmentation threshold Generate initial vegetation cover :when hour, , indicating a vegetated area; when hour, , indicating a non-vegetated area;

[0087] Initial vegetation cover Morphological processing is performed, including first performing an opening operation with a radius of 3 pixels to remove noise points, and then performing a closing operation with a radius of 5 pixels to fill the voids inside the vegetation, resulting in an optimized vegetation mask. ;

[0088] Based on the optimized vegetation cover Corrected images for each band Image segmented into vegetation regions Images of non-vegetated areas : ; .

[0089] S4, different weighting strategies are applied to the vegetated and non-vegetated areas of each band image to fuse them, generating a multispectral fused image, including:

[0090] Based on vegetation cover From the clear images after correction in each band Pixels from vegetation regions are extracted. Fixed fusion weights are assigned based on the spectral response characteristics of vegetation: the near-infrared band is assigned the highest weight of 0.4 due to its high sensitivity to vegetation chlorophyll and cell structure; the infrared and green bands are assigned weights of 0.2 each, as they are important indicators of vegetation growth status; the remaining (n-3) bands are evenly distributed with the remaining weight of 0.2.

[0091] For each vegetation pixel location (x, y), the fusion value is calculated by weighted summation: , where i is the index of other bands. This fixed weight strategy ensures that all vegetation pixels can fully reflect their characteristic spectral information.

[0092] For non-vegetated areas, an adaptive fusion strategy based on image quality is employed. The transmittance optimized during atmospheric correction is utilized. As a metric for image quality across all bands, for each non-vegetated pixel location (x, y), the sum of transmittance across all bands is first calculated. Then calculate the normalized weights for each band. .

[0093] Bands with higher transmittance are less affected by atmospheric conditions and have better image quality, thus receiving a greater fusion weight. The fusion value for non-vegetated areas is obtained through weighted summation: This adaptive strategy allows each pixel to select the clearest band information for fusion.

[0094] Based on vegetation cover Merge the results of vegetation area fusion and non-vegetation area fusion: The final result is a multispectral fusion image F that retains both the spectral characteristics of vegetation and has high image quality.

[0095] like Figure 3 As shown in step S5, spectral features of the multispectral fusion image are extracted to construct a standard vegetation spectral feature library. Based on the extracted spectral features, an anomaly detection algorithm is used to identify areas deviating from the standard vegetation spectrum, resulting in various vegetation coverage rates, including:

[0096] Based on morphologically optimized vegetation mask The system iterates through the entire image and extracts the spatial coordinates of all vegetation pixels. For each pixel location (x, y) marked as vegetation, the system extracts the spatial coordinates from the atmospherically corrected, sharp image across n bands. to The pixel value at the corresponding position is read to form the basic spectral vector at that position.

[0097] To enhance the discriminative ability of features, the system further calculates the spectral ratio between adjacent bands. For the i-th band and the (i+1)-th band, the ratio is calculated. A total of (n-1) ratio features are generated. These ratio features reflect the rate of change of the spectral curve between different bands, especially the sharp change between the red and near-infrared bands (i.e., the red edge effect), and are important indicators of vegetation health.

[0098] The n original spectral values ​​are concatenated with (n-1) spectral ratios to construct a comprehensive feature vector v(x, y) of dimension (2n-1). This feature representation method preserves both the absolute intensity information of the spectrum and the shape features of the spectral curve, providing a rich information foundation for subsequent vegetation classification and anomaly detection.

[0099] The system calculates the normalized vegetation index for all vegetation pixels. ,in These are the corrected pixel values ​​for the near-infrared and red light bands, respectively. All vegetation pixels are sorted in descending order of NDVI value to form an ordered list. The pixels at the top P% (typically 20%) of the list are selected as representative samples of healthy vegetation. These high NDVI value pixels usually correspond to vegetation areas with high chlorophyll content and vigorous growth.

[0100] Extract the corresponding (2n-1)-dimensional feature vectors from the selected healthy vegetation pixels to form the training dataset. The dataset size is M×(2n-1), where M is the number of healthy vegetation samples. Apply the K-means clustering algorithm to this dataset to divide the healthy vegetation into k subcategories. The clustering process iteratively optimizes to maximize the feature similarity of samples within the same category and the difference between samples from different categories. For each clustering result, calculate the mean of the feature vectors of all samples in that category. The central features of this class are obtained. Simultaneously, the covariance matrix is ​​calculated. This describes the distribution range of this type of feature and the correlation between its dimensions. The statistical parameters of the k categories are then compared... This data is stored as standard vegetation spectral characteristics to form a feature library. Each category may represent different vegetation types (such as broadleaf forests, coniferous forests, and grasslands) or different growth stages of the same vegetation.

[0101] For each pixel in the vegetation area, calculate the Mahalanobis distance between its feature vector v(x, y) and various standard features in the feature library. Mahalanobis distance By taking into account the scale differences and correlations of each dimension of the features, the degree of anomalousness of spectral features can be measured more accurately.

[0102] For each vegetation pixel, find the minimum Mahalanobis distance. And the corresponding category index. Calculate the average of the standard deviations of each dimension of the category's feature vector. This serves as the benchmark for determining the threshold. When If the pixel is healthy, it is considered to belong to the corresponding healthy vegetation type; otherwise, it is marked as abnormal vegetation. The 3σ criterion is based on the normal distribution assumption and can theoretically cover 99.7% of normal samples. Therefore, pixels outside this range are likely to have abnormal growth or pests and diseases.

[0103] Traverse the entire image and count the number of pixels for each type of healthy vegetation. and the number of pixels of abnormal vegetation Calculate the coverage rate of various types of healthy vegetation: Abnormal vegetation coverage: ,in, The image contains a total number of pixels. A vegetation distribution map is generated, using different colors to mark various types of healthy and abnormal vegetation areas, providing intuitive visualization results for ecological environment assessment. In the early stages of pests and diseases, vegetation may appear normal in the visible light spectrum, but its spectral characteristics have already undergone subtle changes. These changes can be detected early through anomaly detection. This application does not require retraining the model for each newly emerging pest or disease; as long as the spectral characteristics it causes deviate from the normal range, it can be detected. It not only calculates vegetation coverage but also assesses vegetation health, providing a more comprehensive evaluation index for ecological environment quality.

[0104] The foregoing illustrative description of the present application and its embodiments is not restrictive and can be implemented in other specific forms without departing from the spirit or essential characteristics of the present application. The accompanying drawings are only one embodiment of the present application, and the actual structure is not limited thereto. Therefore, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the present application, such designs should fall within the scope of protection of this application. Furthermore, the word "comprising" does not exclude other elements or steps, and the word "a" preceding an element does not exclude the inclusion of "a plurality" of that element. Terms such as "first," "second," etc., are used to indicate names and do not indicate any specific order.

Claims

1. An ecological environment monitoring method based on multispectral remote sensing fusion, characterized in that, include: S1, acquire multispectral remote sensing image data of the target area, and obtain a set of original images containing n bands and their corresponding center wavelengths; S2, calculate the scattering coefficient, atmospheric light value and transmittance between the bands respectively; Multi-band joint atmospheric correction processing is performed on multispectral remote sensing images based on the scattering coefficient, atmospheric light value, and transmittance of each band. S3, perform foreground and background analysis on the corrected multispectral remote sensing image to identify vegetated and non-vegetated areas; S4. Different weighting strategies are applied to the vegetation areas and non-vegetation areas of each band image to fuse them and generate a multispectral fused image. S5. Extract spectral features from the multispectral fusion image, construct a standard vegetation spectral feature library, and identify areas that deviate from the standard vegetation spectrum based on the extracted spectral features using an anomaly detection algorithm to obtain the coverage of various types of vegetation. S2 is used to calculate the scattering coefficient, atmospheric light value, and transmittance between wavebands, including: Establish the relationship between scattering coefficient and wavelength: ,in, Let λ be the scattering coefficient at wavelength λ. Reference wavelength The scattering coefficient at that point, where α is the Angstrom exponent; Based on the relationship, calculate the scattering coefficient ratio between each band and generate the scattering coefficient ratio matrix R. Select the band with the longest wavelength as the reference band, and then use the reference band image. The image is divided into m x m blocks. The maximum pixel value of each block is calculated, and the average of the top 0.1% of the maximum pixel values ​​of all blocks is selected as the atmospheric light value. ; According to atmospheric light value Calculate the transmittance of the reference band. ; Based on the transmittance of the reference band Calculate the atmospheric light values ​​of other bands i using the scattering coefficient ratio matrix R. , transmittance Where ref represents the reference band, Indicates the atmospheric light value of the reference band; This includes multi-band joint atmospheric correction processing of multispectral remote sensing images, including: Based on the scattering coefficient, atmospheric light value, and transmittance of each band, a multi-band joint optimization objective function is constructed: ； in, Let be the observed pixel value of band i at coordinates (x, y). Let be the sharp pixel value to be recovered for band i at coordinates (x, y). Let be the transmittance of band i at coordinates (x, y). The atmospheric light value for band i. As the weight for band consistency constraints, Weights for spatial smoothing constraints; and Transmittance Gradient operators in the x and y directions; n is the total number of bands; The objective function E is solved using the alternating direction multiplier method to obtain the optimized transmittance for each band. and atmospheric light value ; Based on the optimized transmittance and atmospheric light value Through formula Clear images were obtained after correction for each band. .

2. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 1, characterized in that: Generate the scattering coefficient ratio matrix R, including: Based on the center wavelength of n bands Initialize the n×n scattering coefficient ratio matrix R; Iterate through each position (i, j) of matrix R, where i is the row index, j is the column index, i, j ∈ [1, n], and i and j are positive integers; For each position (i, j), when i = j, R(i, j) = 1; when i ≠ j, according to the formula... Calculations are performed, in which, The center wavelength of the i-th band is Let α be the center wavelength of the j-th band, and α be the Angstrom exponent.

3. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 2, characterized in that: Calculate the transmittance of the reference band , Where ω is the adjustment factor; This represents the pixel value of the reference band at coordinates (x, y).

4. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 2, characterized in that: S3 identifies vegetated and non-vegetated areas, including: Based on the clear images after correction for each band Red light band images and near-infrared band images were selected respectively; The normalized vegetation index is calculated for each pixel based on the spectral differences between the red band image and the near-infrared band image. Histogram statistical analysis was performed on the normalized vegetation index, and the optimal segmentation threshold between vegetation and non-vegetation was determined by the maximum inter-class variance method. Each pixel is binarized according to the optimal segmentation threshold to generate an initial vegetation mask. Pixels with a value greater than the optimal segmentation threshold are marked as vegetation areas, and pixels with a value less than or equal to the optimal segmentation threshold are marked as non-vegetation areas. Morphological processing of the initial vegetation mask is performed sequentially using opening and closing operations. The opening operation is used to remove isolated noise points, and the closing operation is used to fill the voids and breaks inside the vegetation, resulting in a vegetation mask with optimized spatial continuity. Based on the morphologically processed vegetation mask, the corrected clear images of each band are divided into two subsets: images of vegetated areas and images of non-vegetated areas.

5. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 3, characterized in that: S4 generates a multispectral fused image, including: For images of vegetated areas, fusion weights are calculated based on the sensitivity of each band to the spectral response of vegetation. The weight of the near-infrared band is 0.4, the weight of the infrared band is 0.2, the weight of the green band is 0.2, and the weight of the remaining (n-3) bands is 0.2 / (n-3). For images of non-vegetated areas, based on the optimized transmittance Calculate the fusion weights for each band: ; Based on the fusion weights, the images of vegetated areas and non-vegetated areas are weighted and fused separately to obtain a multispectral fused image.

6. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 3, characterized in that: S5 yields various vegetation coverage rates, including: Based on the morphologically processed vegetation mask, pixels of vegetation regions are extracted from the multispectral fused image and combined with the clear images after correction for each band. Construct a multidimensional spectral feature vector for each vegetation pixel. The multidimensional spectral feature vector includes the pixel values ​​of each band. and the ratio of pixel values ​​in adjacent bands ; From the pixels of the vegetation region, the top P% pixels with the highest normalized vegetation index are selected as healthy vegetation samples. K-means clustering is performed on the multidimensional spectral feature vectors of the healthy vegetation samples to obtain k cluster centers of healthy vegetation types. The mean vector of each cluster is calculated. Covariance Matrix , and store them in the feature library as the spectral features of the i-th type of standard vegetation; For all vegetation region pixels in the multispectral fusion image, calculate the distance between the corresponding multidimensional spectral feature vector and the spectral features of various standard vegetation types in the feature library. ; For each vegetation pixel, select the minimum distance. and the corresponding categories, Pixels with good vegetation type are classified as corresponding healthy vegetation types, while those with poor vegetation type are marked as abnormal vegetation areas. Let be the mean standard deviation of each dimension of the feature vector of the i-th type of healthy vegetation; The ratio of the number of pixels of each type of healthy vegetation to the number of pixels of each type of abnormal vegetation is calculated to obtain the coverage rate of each type of healthy vegetation and the coverage rate of each type of abnormal vegetation.

7. The ecological environment monitoring method based on multispectral remote sensing fusion according to claim 6, characterized in that: Calculate the distance between the corresponding multidimensional spectral feature vector and the spectral features of various standard vegetation types in the feature library. The following formula is used: Where v(x, y) is the multidimensional spectral feature vector of the vegetation pixel at coordinates (x, y). Let be the mean vector of the i-th type of standard vegetation in the feature library. Let be the covariance matrix of the i-th type of standard vegetation.

8. An ecological environment monitoring system based on multispectral remote sensing fusion, used to implement the method described in any one of claims 1 to 7, characterized in that, include: The data acquisition module collects multispectral remote sensing image data of the target area, and obtains a set of original images containing n bands and their corresponding center wavelengths; The atmospheric correction module establishes a physical relationship model of scattering coefficients between bands, calculates the ratio of scattering coefficients, atmospheric light value and transmittance of each band, constructs a multi-band joint optimization objective function to perform multi-band atmospheric correction processing, and obtains clear images after correction for each band. The vegetation recognition module calculates the vegetation index based on the corrected clear image, identifies vegetated and non-vegetated areas through threshold segmentation and morphological processing, and generates an optimized vegetation mask. The image fusion module partitions the image according to the optimized vegetation mask. It adopts a fixed weight fusion strategy based on spectral response characteristics for vegetation areas and an adaptive weight fusion strategy based on transmittance for non-vegetation areas to generate a multispectral fusion image. The detection module extracts multidimensional spectral features of vegetation areas, constructs a standard spectral feature library for healthy vegetation, performs anomaly detection by calculating Mahalanobis distance, identifies abnormal vegetation that deviates from the standard spectrum, and calculates the coverage rate of various types of healthy and abnormal vegetation.

Images