Wetland ecosystem health evaluation method based on multi-source remote sensing data

Abstract

The invention provides a wetland ecosystem health evaluation method based on multi-source remote sensing data, and belongs to the technical field of wetland ecosystems, and the method comprises the steps: carrying out the dense dark pixel and Kalman filtering cooperative atmospheric correction of a multi-spectral remote sensing image to obtain a surface reflectance image, and achieving the precise segmentation of a wetland landscape through a graph cut theory, a spectral unmixing model based on an improved Gaussian kernel is utilized to embed physical constraints to invert water quality parameters, a laser radar canopy height priori constraint three-dimensional radiation transmission model is combined to invert vegetation parameters, and a mixed pixel decomposition result is optimized through a Markov random field. A water quality, vegetation and landscape three-dimensional comprehensive health evaluation system is constructed, degradation dominant factors are identified, and the technical problem that it is difficult for multi-source remote sensing data to cooperatively invert wetland ecosystem multi-dimensional health state parameters is solved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of wetland ecological system, and in particular, relates to a wetland ecological system health evaluation method based on multi-source remote sensing data. BACKGROUND

[0002] Wetland ecological system health evaluation needs to comprehensively consider parameters in multiple dimensions of water quality, vegetation and landscape. Traditional methods mainly rely on multispectral remote sensing images to extract single type parameters, and use empirical band ratio models to retrieve chlorophyll concentration or vegetation index, which lacks the ability to accurately analyze complex water body environment and mixed pixel scenes. In current wetland monitoring practices, due to the lack of vertical structure information in multispectral data and the lack of spectral characteristics in laser radar data, it is difficult to simultaneously obtain multi-dimensional information such as water quality concentration distribution, vegetation physiological parameters and landscape spatial pattern by using a single data source, resulting in one-sidedness of the health evaluation results. The existing methods usually use simple linear unmixing or hard classification strategy when dealing with mixed pixels, which fails to fully utilize the complementary advantages of multi-source data for collaborative inversion, and lacks physical constraint mechanisms specific to the ecological environment of wetlands, resulting in significant decrease in inversion accuracy under conditions of high concentration of suspended solids and complex canopy structure. That is, there is a technical problem in the prior art that multi-source remote sensing data is difficult to collaboratively retrieve multi-dimensional health state parameters of wetland ecological system. SUMMARY

[0003] Therefore, the present application provides a wetland ecological system health evaluation method based on multi-source remote sensing data, which can solve the technical problem in the prior art that multi-source remote sensing data is difficult to collaboratively retrieve multi-dimensional health state parameters of wetland ecological system.

[0004] The present application is implemented as follows: The present application provides a wetland ecological system health evaluation method based on multi-source remote sensing data, which acquires multispectral remote sensing images and laser radar point cloud data of a target wetland area, performs radiation calibration and atmospheric correction on the multispectral remote sensing images to obtain a surface reflectance image, segments the surface reflectance image into four landscape types of open water, emergent vegetation, submerged plants and bare mudflat based on the graph cut theory to obtain a wetland landscape segmentation result image, extracts open water regions and vegetation regions from the wetland landscape segmentation result image, inputs a water quality parameter unmixing model to obtain chlorophyll a concentration distribution, suspended solids concentration distribution and dissolved organic matter concentration distribution, combines the laser radar point cloud data to retrieve leaf area index distribution and chlorophyll content distribution through a three-dimensional radiation transfer model, judges the proportion of mixed pixels in the wetland landscape segmentation result image, uses a total least squares method and a Markov random field model to obtain an optimized abundance coefficient distribution, inputs a wetland health evaluation model to calculate comprehensive health index distribution, water quality health sub-index distribution, vegetation health sub-index distribution and landscape health sub-index distribution, divides wetland health grade distribution and identifies dominant factor types.

[0005] The atmospheric correction adopts a dense dark pixel method to extract an initial value of an atmospheric parameter, combines ground observation data of a sun photometer to obtain an optimal atmospheric correction parameter through Kalman filtering dynamic updating, and introduces a water vapor continuous absorption correction term to calculate a surface reflectance image.

[0006] The dense dark pixel method selects a pixel with a vegetation coverage greater than 80% and a normalized vegetation index greater than 0.7 as a dark pixel candidate set, and calculates the average value of the 5% pixels with the lowest apparent reflectance in the blue light band and the red light band as an estimated value of the atmospheric path radiation to inversely calculate an initial value of the aerosol optical thickness.

[0007] The Kalman filtering takes the initial value of the aerosol optical thickness as a system state variable, establishes a state transition equation, takes the initial value of the aerosol optical thickness obtained by the dense dark pixel method and the observation value of the ground observation data of the sun photometer as measurement variables, and iteratively obtains a time sequence of the regional aerosol optical thickness through a prediction step and an update step.

[0008] The water vapor continuous absorption correction term calculates the water vapor optical thickness of each band according to the atmospheric water vapor column concentration and the water vapor absorption cross section, adopts an exponential decay function to describe the influence of water vapor continuous absorption on radiation transmission to obtain a water vapor absorption transmittance, and multiplies the water vapor absorption transmittance with Rayleigh scattering transmittance and aerosol scattering transmittance to obtain the atmospheric correction equation.

[0009] The graph cut theory takes the pixels of the surface reflectance image as undirected graph nodes, calculates the normalized spectral distance between adjacent pixels as the edge weight, sets four types of terminal nodes, solves the minimum value of the energy function through a minimum cut maximum flow algorithm, and decomposes the multi-label problem into a series of binary label expansion sub-problems through an alpha expansion algorithm and iterates to convergence.

[0010] The normalized spectral distance extracts the reflectance vectors of adjacent pixels in all spectral bands, calculates the Euclidean distance of the two reflectance vectors, and obtains the normalized spectral distance as the edge weight by dividing the Euclidean distance by the maximum modulus length of the reflectance vector.

[0011] The water quality parameter unmixing model is a spectral unmixing model based on an improved Gaussian kernel, the improved Gaussian kernel introduces four types of physical constraint conditions including a water quality concentration range constraint, a spectral physical constraint, an absorption peak position constraint and a scattering angle distribution constraint, increases a constraint term on the basis of the traditional Gaussian kernel, and integrates the constraint term into the kernel function optimization target through a Lagrange multiplier method.

[0012] The water quality parameter unmixing model is a five-layer feedforward architecture, the input layer receives principal component feature vectors, ultraviolet band reflectance and linear fluorescence height index, the first hidden layer contains 128 neurons and uses an improved Gaussian kernel as an activation function, the second hidden layer contains 64 neurons, the third hidden layer contains 32 neurons, and the output layer contains 3 neurons.

[0013] Before inputting the water quality parameter unmixing model, the multispectral band reflectance data of the open water area is subjected to principal component analysis dimension reduction, the eigenvalues and eigenvectors of the covariance matrix are calculated, the first three eigenvectors with a cumulative variance contribution rate of 95% are selected according to the descending order of the eigenvalues, and the principal component feature vectors are obtained.

[0014] The three-dimensional radiative transfer model discretizes the wetland vegetation canopy into a three-dimensional voxel grid, each three-dimensional voxel grid is assigned with leaf area density, leaf inclination distribution and leaf optical property parameters, the Monte Carlo ray tracing method is used to simulate the propagation path of direct sunlight and sky scattered light in the vegetation canopy, and the energy distribution of a large number of light rays is counted to obtain the simulated canopy reflectance.

[0015] The three-dimensional radiative transfer model establishes a lookup table, Latin hypercube sampling is used to generate parameter combinations in the parameter space, the three-dimensional radiative transfer model is run to calculate the corresponding multispectral reflectance, and a mapping lookup table of parameter combinations and multispectral reflectance is established.

[0016] The three-dimensional radiative transfer model establishes a lookup table, Latin hypercube sampling is used to generate parameter combinations in the parameter space, the three-dimensional radiative transfer model is run to calculate the corresponding multispectral reflectance, and a mapping lookup table of parameter combinations and multispectral reflectance is established.

[0017] The full-constraint least squares method performs linear spectral mixture decomposition on the mixed pixel, solves the constraint condition that the abundance coefficients are non-negative and the sum is 1, uses the interior point method to solve the constraint optimization problem, iteratively calculates the optimal abundance coefficient to minimize the mean square error between the simulated reflectance and the measured reflectance, and obtains the abundance coefficient distribution of each land cover type.

[0018] The Markov random field model defines the neighborhood of a pixel as its eight-connected adjacent pixels, establishes an energy function containing a data fidelity term and a smoothing regularization term, the data fidelity term measures the deviation of the optimized abundance coefficient from the abundance coefficient distribution of each land cover type, and the smoothing regularization term imposes a penalty on the difference between the abundance coefficients of adjacent pixels, and the graph cut algorithm or the iterative conditional mode algorithm is used to minimize the energy function.

[0019] The wetland health evaluation model is used to divide the wetland health grade distribution, and the contribution rate distribution of the water quality health index distribution, the vegetation health index distribution and the landscape health index distribution to the reduction of the comprehensive health index is calculated, and the type of the index with the maximum contribution rate is identified as the dominant factor type.

[0020] The present application realizes the synchronous and accurate extraction of water quality parameters and vegetation parameters by constructing a collaborative inversion framework of multispectral remote sensing and laser radar data, deeply fusing spectral information and vertical structure information. In view of the defects of the traditional method, the present application introduces a spectral unmixing model based on improved Gaussian kernel to embed water quality physical constraint conditions, which significantly improves the accuracy of multi-component concentration inversion in complex water environment. At the same time, the present application uses the crown height prior information extracted by laser radar to constrain the parameter search space of the three-dimensional radiation transfer model, effectively solving the ill-posed problem in vegetation parameter inversion, and through the Markov random field model to optimize the mixed pixel decomposition result, enhancing the spatial continuity of the landscape pattern. The present application establishes a comprehensive evaluation system of water quality health, vegetation health and landscape health, which can accurately identify the dominant driving factor of the degradation area, and solves the technical problem that the existing technology is difficult to collaboratively invert the multi-dimensional health state parameters of the wetland ecosystem by multi-source remote sensing data. BRIEF DESCRIPTION OF DRAWINGS

[0021] Figure 1 The figure is a spatial distribution map of the surface reflectance image after atmospheric correction in the embodiment.

[0022] Figure 2 The figure is a spatial distribution map of the wetland landscape segmentation result image in the embodiment.

[0023] Figure 3 The figure is a spatial distribution map of the water quality parameter in the embodiment.

[0024] Figure 4 The figure is a spatial distribution map of the vegetation parameter in the embodiment.

[0025] Figure 5 The figure is a spatial distribution map of the wetland health grade distribution and the dominant factor type in the embodiment. DETAILED DESCRIPTION

[0026] In order to make the purpose, technical scheme and advantages of the embodiments of the present application more clear, the technical scheme in the embodiments of the present application is described clearly and completely below.

[0027] The present application provides a wetland ecosystem health evaluation method of multi-source remote sensing data, comprising:

[0028] S01, acquire the multispectral remote sensing image and the laser radar point cloud data of the target wetland area, perform radiation calibration on the multispectral remote sensing image to convert it into an apparent reflectance image, extract the initial value of an atmospheric parameter by using a dense dark object method, combine the ground observation data of a sun photometer to obtain optimal atmospheric correction parameters by using Kalman filtering for dynamic updating, introduce a water vapor continuous absorption correction term to calculate a ground reflectance image;

[0029] S02, based on the graph cut theory, take the pixels of the ground reflectance image as nodes of an undirected graph, calculate the normalized spectral distance between adjacent pixels as an edge weight, set four types of terminal nodes, i.e., open water, emergent vegetation, submerged plants and bare mudflat, solve the minimum value of an energy function by using a minimum cut maximum flow algorithm to obtain a wetland landscape segmentation result image;

[0030] S03, extract an open water area from the wetland landscape segmentation result image, perform principal component analysis on the multispectral band reflectance data of the open water area to obtain a principal component feature vector, extract the ultraviolet band reflectance and the fluorescence linear height index, and input them into a water quality parameter unmixing model to obtain the concentration distribution of chlorophyll a, the concentration distribution of suspended matter and the concentration distribution of dissolved organic matter; extract a vegetation area from the wetland landscape segmentation result image, use a three-dimensional radiation transfer model to establish a lookup table of ground reflectance and leaf area index and chlorophyll content, combine the prior information of canopy height extracted from the laser radar point cloud data to constrain the parameter search space and obtain the distribution of leaf area index and chlorophyll content by inversion;

[0031] S04, judge whether the proportion of mixed pixels in the wetland landscape segmentation result image is greater than 30%, when the proportion of mixed pixels is greater than 30%, extract pure endmember spectra from a high-resolution remote sensing image to establish an endmember spectral library, use a total constrained least squares method to perform linear spectral mixture decomposition on the mixed pixels to obtain the abundance coefficient distribution of each land cover type, introduce a Markov random field model to perform spatial neighborhood smoothing optimization to obtain the optimized abundance coefficient distribution, and when the proportion of mixed pixels is less than or equal to 30%, set the abundance coefficient distribution of each land cover type as the class purity distribution of the wetland landscape segmentation result image;

[0032] S05, inputting the chlorophyll-a concentration distribution, the suspended matter concentration distribution, the dissolved organic matter concentration distribution, the leaf area index distribution, the chlorophyll content distribution and the optimized abundance coefficient distribution into a wetland health evaluation model, calculating and outputting a comprehensive health index distribution, a water quality health sub-index distribution, a vegetation health sub-index distribution and a landscape health sub-index distribution, dividing a wetland health grade distribution according to the comprehensive health index distribution, judging whether there is a degraded state or a severely degraded state of a pixel region, when there is, calculating a contribution rate distribution of the water quality health sub-index distribution, the vegetation health sub-index distribution and the landscape health sub-index distribution to the decrease of the comprehensive health index, identifying the sub-index type with the largest contribution rate as a dominant factor type, and outputting the wetland health grade distribution and the dominant factor type.

[0033] The processing steps of the dense dark pixel method specifically include: selecting pixels with vegetation coverage greater than 80% and normalized difference vegetation index greater than 0.7 in the apparent reflectance image as a dark pixel candidate set, counting the apparent reflectance distribution of the dark pixel candidate set in the blue light band and the red light band, taking the average value of the 5% pixels with the lowest apparent reflectance in the apparent reflectance distribution as an atmospheric path radiation estimation value, and inversely calculating an aerosol optical depth initial value according to the linear relationship between the atmospheric path radiation estimation value and the aerosol optical depth, the aerosol optical depth initial value being the atmospheric parameter initial value.

[0034] The step of dynamically updating the atmospheric parameter initial value by the Kalman filter specifically includes: taking the aerosol optical depth initial value as a system state variable, establishing a state transition equation to describe the time evolution trend of the aerosol optical depth initial value, taking the aerosol optical depth initial value obtained by the dense dark pixel method and the observation value of the sun photometer ground observation data as measurement variables, calculating a prior state estimation and a prior error covariance through a prediction step, calculating a Kalman gain and correcting the prior state estimation through an update step, iteratively executing a prediction update cycle to obtain a regional aerosol optical depth time series, and extracting an optimal value of the aerosol optical depth at the current time from the regional aerosol optical depth time series as the optimal atmospheric correction parameter.

[0035] The calculation method of the water vapor continuous absorption correction term is: calculating the water vapor optical thickness of each band according to the atmospheric water vapor column concentration and the water vapor absorption cross section, using an exponential decay function to describe the influence of water vapor continuous absorption on radiation transmission to obtain a water vapor absorption transmittance, multiplying the water vapor absorption transmittance with Rayleigh scattering transmittance and aerosol scattering transmittance, and substituting the product into an atmospheric correction equation to convert the apparent reflectance image into the surface reflectance image through the atmospheric correction equation.

[0036] The energy function of the wetland landscape segmentation algorithm constructed by the graph cut theory is composed of a data item and a smoothing item. The data item quantifies the Mahalanobis distance between the spectral vector of a pixel and the center spectral vector of each class. The smoothing item imposes a penalty on the assignment of different labels to adjacent pixels, and the penalty weight of the smoothing item is inversely proportional to the edge weight. The Boykov-Kolmogorov maximum flow algorithm is used to solve the graph cut problem, and the minimum cut set is iteratively calculated through the augmented path search and residual graph update. The pixel label configuration corresponding to the minimum cut set is the global minimum solution of the energy function. For multi-label segmentation of the four types of terminal nodes, the alpha expansion algorithm is used to decompose the multi-label problem into a series of binary label expansion sub-problems. In each iteration, a label is selected as an alpha label to attempt to expand its area. The graph cut solution accepts or rejects the expansion, and the wetland landscape segmentation result image is obtained by iteration to convergence.

[0037] The normalized spectral distance is calculated by extracting the reflectance vectors of adjacent pixels in all spectral bands, calculating the Euclidean distance between the two reflectance vectors, and dividing the Euclidean distance by the maximum modulus of the reflectance vector to obtain the normalized spectral distance. The normalized spectral distance ranges from 0 to 2, and the smaller the value, the higher the spectral similarity. The normalized spectral distance is the edge weight.

[0038] The method of principal component analysis dimension reduction is as follows: covariance matrix calculation is performed on the multi-spectral band reflectance data of the open water area, the eigenvalues and eigenvectors of the covariance matrix are solved, the first three eigenvectors with a cumulative variance contribution rate of 95% are selected according to the descending order of eigenvalues, and the multi-spectral band reflectance data is projected onto the first three eigenvectors to obtain the principal component eigenvectors.

[0039] The ultraviolet band reflectance refers to the surface reflectance of ultraviolet light, and the center wavelength of the ultraviolet light band is between 350 nm and 400 nm. The fluorescence linear height index is defined as the difference between the reflectance of the fluorescence peak band and the baseline of the adjacent band. The fluorescence peak band is located near 685 nm, and the adjacent bands are located at 665 nm and 705 nm, respectively.

[0040] The water quality parameter unmixing model is a spectral unmixing model based on an improved Gaussian kernel. The improved Gaussian kernel introduces four types of physical constraints, including water quality concentration range constraints, spectral physical constraints, absorption peak position constraints, and scattering angle distribution constraints. The improved Gaussian kernel adds constraint terms to the traditional Gaussian kernel, including reasonable range constraints of chlorophyll-a concentration, reasonable range constraints of suspended matter concentration, non-negativity constraints of absorption coefficients of each component, and physical range constraints of backscattering ratio. The constraint terms are integrated into the kernel function optimization objective through the Lagrange multiplier method.

[0041] The structure of the water quality parameter unmixing model is a five-layer feedforward architecture. The input layer receives five variables, including the principal component feature vector, the ultraviolet band reflectivity, and the fluorescence linear height index. The first hidden layer contains 128 neurons and uses the improved Gaussian kernel as the activation function. The second hidden layer contains 64 neurons and uses the rectified linear unit activation function. The third hidden layer contains 32 neurons and uses the rectified linear unit activation function. The output layer contains three neurons corresponding to the chlorophyll-a concentration distribution, the suspended matter concentration distribution, and the dissolved organic matter concentration distribution, respectively. The water quality parameter unmixing model adopts a residual connection structure to directly connect the input layer to the second hidden layer and adopts a batch normalization layer to accelerate training convergence.

[0042] The steps of establishing the training data set of the water quality parameter unmixing model specifically include: collecting wetland water remote sensing image samples under different seasons and different weather conditions, synchronously collecting water surface samples for laboratory analysis to obtain measured values of chlorophyll-a concentration, suspended matter concentration, and dissolved organic matter concentration, performing atmospheric correction processing on the wetland water remote sensing image samples according to step S01 to obtain ground reflectance samples, extracting the multispectral reflectivity, ultraviolet band reflectivity, and fluorescence linear height index of water pixels in the ground reflectance samples, grouping the multispectral reflectivity, ultraviolet band reflectivity, and fluorescence linear height index into a feature vector, and pairing the feature vector with the measured values of chlorophyll-a concentration, suspended matter concentration, and dissolved organic matter concentration to form training samples, dividing the training set and the validation set in a ratio of 8 to 2, and the training set contains no less than 5000 groups of training samples.

[0043] The steps of training the water quality parameter unmixing model specifically include: updating the weight parameters of the water quality parameter unmixing model using the Adam optimization algorithm, setting the initial learning rate to 0.001, decaying the initial learning rate to 0.8 times the original value every 20 training periods, setting the batch size to 64, setting the total number of training periods to 200, using the weighted sum of the mean square error and the physical constraint violation penalty term as the loss function, applying a quadratic penalty to the prediction of concentration values that exceed the reasonable range, and setting the penalty weight coefficient to 10. In the training process, the performance of the water quality parameter unmixing model is evaluated on the validation set every 5 periods, the training is terminated in advance when the validation set loss does not decrease for 10 consecutive periods, and the model weight with the minimum validation set loss is selected as the final model parameter.

[0044] The improved Gaussian kernel-based spectral unmixing mechanism significantly improves the generalization ability and inversion accuracy of the model by introducing water quality physical constraints. The traditional Gaussian kernel only considers the statistical distribution of input features, which can easily produce prediction results that do not conform to physical laws in sparse areas of training data. The improved Gaussian kernel embeds prior knowledge such as reasonable range constraints of chlorophyll-a concentration, reasonable range constraints of suspended sediment concentration, non-negative constraints of component absorption coefficients, and physical range constraints of backscattering ratio in the kernel function, so that the water quality parameter unmixing model automatically avoids unreasonable parameters during parameter space search, reducing the uncertainty of the solution space of the ill-posed inversion problem. During the activation process of the first hidden layer, the improved Gaussian kernel-based spectral unmixing mechanism performs constraint checking on the output value of each neuron. When the output value corresponds to a water quality parameter that exceeds the physically reasonable range, a soft constraint penalty function is used to adjust the output value of the neuron, so that the output value returns to the feasible region. This physically embedded design enables the water quality parameter unmixing model to effectively distinguish the spectral contributions of different components in complex water environments, especially in the presence of high concentrations of suspended sediment and dissolved organic matter, avoiding the confusion caused by spectral overlap in traditional band ratio methods. Experimental verification shows that after introducing the improved Gaussian kernel-based spectral unmixing mechanism, the root mean square error of the chlorophyll-a concentration distribution inversion is reduced from 50% to 18%, and the root mean square error of the suspended sediment concentration distribution inversion is reduced from 45% to 22%. The water quality parameter unmixing model still maintains good prediction performance for lake and wetland samples outside the training data, proving that the improved Gaussian kernel-based spectral unmixing mechanism significantly enhances the physical consistency and transfer generalization ability of the water quality parameter unmixing model.

[0045] The method for simulating the multiple scattering process of the vegetation canopy in the three-dimensional radiation transfer model is as follows: the wetland vegetation canopy is discretized into a three-dimensional voxel grid, each voxel grid is assigned with leaf area density, leaf inclination distribution, and leaf optical property parameters, and the Monte Carlo ray tracing method is used to simulate the propagation path of direct sunlight and sky scattered light in the vegetation canopy. When each light ray propagates in the three-dimensional voxel grid, the scattering probability is calculated according to the extinction coefficient, and the new propagation direction is determined according to the phase function when scattering occurs. The light ray undergoes multiple reflections and transmissions between leaves, stems, and water surfaces, records the light energy that finally escapes to the sensor direction, and obtains the simulated value of canopy albedo by statistical energy distribution of a large number of light rays.

[0046] The lookup table establishment method is as follows: the value range of the leaf area index is set to 0.5-8.0, the value range of the chlorophyll content is set to 20-80 micrograms per square centimeter, and the value range of the canopy height is set to 0.5-5.0 meters; 1000 sets of parameter combinations are generated in a parameter space by using Latin hypercube sampling; the three-dimensional radiation transfer model is run for each set of the parameter combinations to calculate a corresponding multi-spectral reflectance; a mapping lookup table of the parameter combinations and the multi-spectral reflectance is established; and in the inversion, the measured reflectance of the surface reflectance image is compared with the simulated canopy reflectance values in the lookup table one by one, and the parameter combination with the minimum error is selected as the inversion result.

[0047] The method for constraining the parameter search space by the canopy height prior information is as follows: a canopy height model is extracted from the laser radar point cloud data; for each vegetation pixel, a measured value of the canopy height thereof is obtained; the parameter combinations in the lookup table, for which the deviation between the canopy height and the measured value of the canopy height is greater than 0.5 meters, are removed; and the measured reflectance matching is only performed in the remaining parameter combinations. The method can reduce the parameter search space to 30% of the original size, and significantly improve the inversion efficiency and accuracy.

[0048] The calculation method of the mixed pixel proportion is as follows: the class label uncertainty index of each pixel in the wetland landscape segmentation result image is counted; when the class label uncertainty index is greater than a threshold value 0.3, the pixel is determined as a mixed pixel; and the mixed pixel proportion is obtained by dividing the number of mixed pixels by the total number of pixels.

[0049] The mathematical expression of the full-constrained least square method for linear spectral mixture decomposition of the mixed pixels of the surface reflectance image is as follows: the reflectance of a mixed pixel is equal to the linear combination of the reflectance of each endmember and the abundance coefficient thereof; the solution of the abundance coefficient needs to satisfy the constraint conditions that all the abundance coefficients are non-negative and the sum is 1; the interior point method is used to solve the constraint optimization problem; the optimal abundance coefficient is iteratively calculated to minimize the mean square error between the simulated reflectance and the measured reflectance; and the optimal abundance coefficient is the abundance coefficient distribution of each land cover type.

[0050] The method for spatial neighborhood smoothing optimization of the abundance coefficient distribution of each land cover type by the Markov random field model is as follows: the neighborhood of a pixel is defined as its eight-connected neighboring pixels; an energy function is established to include a data fidelity term and a smoothing regularization term; the data fidelity term measures the deviation between the optimized abundance coefficient and the abundance coefficient distribution of each land cover type; the smoothing regularization term imposes a penalty on the difference between the abundance coefficients of adjacent pixels; the graph cut algorithm or the iterative conditional mode algorithm is used to minimize the energy function, so that the optimized abundance coefficient distribution with enhanced spatial continuity is obtained.

[0051] The definition of the category purity distribution is that when the mixed pixel ratio is less than or equal to 30%, the pixel category distribution result of the wetland landscape segmentation result image is considered to be basically reliable, the abundance corresponding to the dominant category label of each pixel is directly set to 1, and the abundance of the remaining categories is set to 0 to form the category purity distribution as the replacement input of the abundance coefficient distribution of each feature type.

[0052] The calculation method of the wetland health evaluation model is that the comprehensive health index distribution is equal to the weighted sum of the water quality health sub-index distribution, the vegetation health sub-index distribution and the landscape health sub-index distribution, and the weight coefficients are determined according to the wetland type and ecological function. For the wetland with water quality purification as the main function, the weight of the water quality health sub-index distribution is set to 0.5, the weight of the vegetation health sub-index distribution is set to 0.3, and the weight of the landscape health sub-index distribution is set to 0.2. The water quality health sub-index distribution is calculated by comparing the chlorophyll a concentration distribution, the suspended matter concentration distribution and the dissolved organic matter concentration distribution with the water quality standard threshold value. When the concentration is lower than the water quality standard threshold value, the sub-index is close to 1, and when the concentration is much higher than the water quality standard threshold value, the sub-index is close to 0. The vegetation health sub-index distribution is calculated by comparing the leaf area index distribution and the chlorophyll content distribution with the healthy vegetation benchmark value. When the parameters are close to the healthy vegetation benchmark value, the sub-index is close to 1, and when the parameters deviate significantly from the healthy vegetation benchmark value, the sub-index decreases. The landscape health sub-index distribution is calculated by the diversity index and the connectivity index of the optimized abundance coefficient distribution. When the abundance distribution is uniform and the spatial connectivity is high, the sub-index is higher, and when a single type dominates or the fragmentation is serious, the sub-index decreases.

[0053] The division rule of the wetland health grade distribution is that when the comprehensive health index of the pixel in the comprehensive health index distribution is greater than 0.8, the pixel is marked as healthy state, when the comprehensive health index is between 0.6 and 0.8, the pixel is marked as sub-healthy state, when the comprehensive health index is between 0.4 and 0.6, the pixel is marked as degradation state, and when the comprehensive health index is less than 0.4, the pixel is marked as serious degradation state.

[0054] The calculation method of the contribution rate distribution is to calculate the contribution rate of the water quality health sub-index distribution, the vegetation health sub-index distribution and the landscape health sub-index distribution to the decrease of the comprehensive health index distribution respectively. The contribution rate is equal to the product of each sub-index and its weight coefficient divided by the theoretical maximum value of the comprehensive health index distribution. The ecological element corresponding to the sub-index type with the maximum contribution rate is the dominant factor type.

[0055] Optionally, the application also provides a computer-implemented way to form a wetland ecosystem health evaluation system of multi-source remote sensing data, wherein a readable storage medium is arranged in the computer, and program instructions are stored in the readable storage medium, and the program instructions execute the wetland ecosystem health evaluation method of multi-source remote sensing data when running in the computer.

[0056] The specific implementation of the above steps is described in detail below.

[0057] The specific implementation of step S01 is to first obtain multispectral remote sensing images and lidar point cloud data of the target wetland area, wherein the multispectral remote sensing images should include visible light bands, near-infrared bands and ultraviolet bands, and the lidar point cloud data is used to extract canopy height information subsequently. The purpose of this step is to provide a raw data basis for subsequent processing. Radiometric calibration is performed on the obtained multispectral remote sensing images to convert the digital quantization values of the images into apparent reflectance images. This process calculates the top-of-atmosphere reflectance using sensor calibration coefficients and solar irradiance parameters. Radiometric calibration eliminates the influence of sensor response differences on reflectance calculation. The initial value of the atmospheric parameter is extracted using the dense dark object method. This method is based on the assumption that the apparent reflectance of dark objects in the blue and red bands should be close to zero. In the apparent reflectance image, the pixels with a vegetation coverage greater than 80% and a normalized vegetation index greater than 0.7 are selected as the dark object candidate set. These pixels are usually areas covered by dense vegetation, and the surface reflectance is extremely low. The apparent reflectance distribution of the dark object candidate set in the blue and red bands is statistically analyzed. The average value of the apparent reflectance of the 5% pixels with the lowest apparent reflectance in the apparent reflectance distribution is taken as the estimated value of the atmospheric path radiation. According to the linear relationship between the atmospheric path radiation and the aerosol optical thickness, the initial value of the aerosol optical thickness is calculated. The initial value of the aerosol optical thickness is dynamically updated using Kalman filtering combined with solar photometer ground observation data. The initial value of the aerosol optical thickness is taken as the system state variable, and a state transition equation is established to describe the time evolution trend of the aerosol optical thickness. This equation assumes that the aerosol optical thickness changes slowly in a short time. The initial value of the aerosol optical thickness obtained by the dense dark object method and the observation value of the solar photometer ground observation data are taken as the measurement variables. The prior state estimation and prior error covariance are calculated through the prediction step. The Kalman gain is calculated through the update step and the prior state estimation is corrected. The Kalman gain is dynamically adjusted according to the relative size of the prediction error and the measurement error. The prediction update cycle is iteratively executed to obtain the time series of the regional aerosol optical thickness. The optimal value of the aerosol optical thickness at the current time is extracted from the time series as the optimal atmospheric correction parameter. This method combines the advantages of remote sensing inversion and ground observation, and improves the estimation accuracy of the atmospheric parameter. A water vapor continuous absorption correction term is introduced to calculate the surface reflectance image. According to the atmospheric water vapor column concentration and the water vapor absorption cross section, the water vapor optical thickness of each band is calculated. The influence of water vapor continuous absorption on radiation transmission is described by an exponential decay function to obtain the water vapor absorption transmittance. The water vapor absorption transmittance is multiplied by the Rayleigh scattering transmittance and the aerosol scattering transmittance, and then substituted into the atmospheric correction equation. The apparent reflectance image is converted into the surface reflectance image through the atmospheric correction equation. The introduction of the water vapor continuous absorption correction term can correct the continuous absorption effect of water vapor in the near-infrared band, and improve the surface reflectance inversion accuracy. The output of this step is the surface reflectance image after accurate atmospheric correction, which provides high-quality basic data for subsequent wetland landscape segmentation and parameter inversion.

[0058] The specific implementation of step S02 is that the pixels of the surface reflectance image are taken as nodes of an undirected graph based on the graph cut theory, and the undirected graph is a data structure in which the nodes represent the pixels and the edges represent the relationship between the pixels. The purpose of this step is to convert the wetland landscape segmentation problem into a graph optimization problem. The normalized spectral distance between adjacent pixels is calculated as the edge weight. The reflectance vectors of adjacent pixels in all spectral bands are extracted, the Euclidean distance between the two reflectance vectors is calculated, the Euclidean distance measures the degree of similarity of the spectral characteristics of the two pixels, the normalized spectral distance is obtained by dividing the Euclidean distance by the maximum modulus of the reflectance vector, and the normalized processing eliminates the influence of the range difference of the reflectance value. The normalized spectral distance has a value range of 0 to 2, and the smaller the value, the higher the spectral similarity. This distance is the edge weight. Four types of terminal nodes, open water, emergent vegetation, submerged plants and bare mudflat, are set, and the terminal nodes represent seed pixels of known classes, which are determined by manual annotation or prior knowledge. These nodes have fixed class labels during the optimization process. An energy function of the graph cut theory is constructed, which consists of a data term and a smoothing term. The data term quantifies the Mahalanobis distance between the spectral vector of a pixel and the center spectral vector of each class. The Mahalanobis distance takes into account the covariance relationship between the bands and can more accurately measure the spectral similarity. The smoothing term imposes a penalty on adjacent pixels assigned different labels. The penalty weight of the smoothing term is inversely proportional to the edge weight. Adjacent pixels with high spectral similarity will be penalized more when they are assigned different classes, which encourages spatial continuity. The Boykov-Kolmogorov maximum flow algorithm is used to solve the graph cut problem. This algorithm calculates the minimum cut set through augmented path search and residual graph update iteration. The pixel label configuration corresponding to the minimum cut set is the global minimum solution of the energy function. This algorithm can obtain the global optimal solution in polynomial time, avoiding the local optimal trap. For multi-label segmentation of the four types of terminal nodes, the alpha expansion algorithm is used to decompose the multi-label problem into a series of binary label expansion sub-problems. Each iteration selects a label as the alpha label to try to expand its area. The graph cut is used to accept or reject the expansion, and the iteration is converged to obtain the wetland landscape segmentation result image. This algorithm ensures efficient solution of the multi-label segmentation problem, and the output segmentation result image not only guarantees spectral consistency but also has good spatial continuity.

[0059] The specific implementation of step S03 is to extract the open water area from the wetland landscape segmentation result image, the pixel class label of the area is open water, and the multispectral band reflectivity data of these pixels is extracted. The purpose of this step is to perform parameter inversion for water body and vegetation area respectively. Principal component analysis dimension reduction is performed on the multispectral band reflectivity data of the open water area, the covariance matrix of the multispectral band reflectivity data is calculated, the eigenvalues and eigenvectors of the covariance matrix are solved, the first three eigenvectors with a cumulative variance contribution rate of 95% are selected according to the descending order of eigenvalues. Principal component analysis can project high-dimensional spectral data into low-dimensional space while retaining the main information, reducing data redundancy and computational complexity. The multispectral band reflectivity data is projected into the first three eigenvectors to obtain the principal component eigenvectors. The ultraviolet band reflectivity and the fluorescence linear height index are extracted. The ultraviolet band reflectivity refers to the surface reflectivity of the ultraviolet band, and the center wavelength of the ultraviolet band is between 350 nm and 400 nm. The fluorescence linear height index is defined as the difference between the fluorescence peak band reflectivity and the baseline of the adjacent band reflectivity. The fluorescence peak band is located near 685 nm, and the adjacent bands are located at 665 nm and 705 nm, respectively. The fluorescence linear height index can effectively represent the chlorophyll fluorescence characteristics, which is sensitive to water body chlorophyll concentration. The principal component eigenvectors, ultraviolet band reflectivity and fluorescence linear height index are input into the water quality parameter unmixing model, which is a spectral unmixing model based on improved Gaussian kernel. The model adopts a five-layer feedforward architecture. The input layer receives five input variables. The first hidden layer contains 128 neurons and uses an improved Gaussian kernel as the activation function. The improved Gaussian kernel introduces four types of physical constraints: water concentration range constraint, spectral physical constraint, absorption peak position constraint and scattering angle distribution constraint. The second hidden layer contains 64 neurons using a rectified linear unit activation function. The third hidden layer contains 32 neurons using a rectified linear unit activation function. The output layer contains three neurons corresponding to chlorophyll concentration distribution, suspended matter concentration distribution and dissolved organic matter concentration distribution, respectively. The model uses a residual connection structure to connect the input layer directly to the second hidden layer, and uses a batch normalization layer to accelerate training convergence. During the activation process in the first hidden layer, the output value of each neuron is checked for constraints. When the output value of the water quality parameter exceeds the physically reasonable range, the output value is adjusted to return to the feasible region by using a soft constraint penalty function. The chlorophyll Concentration distribution, suspended matter concentration distribution and dissolved organic matter concentration distribution. The vegetation area including the pixels of emergent vegetation and submerged plant categories is extracted from the wetland landscape segmentation result image, a lookup table of surface reflectance and leaf area index, chlorophyll content is established by using a three-dimensional radiation transfer model, the three-dimensional radiation transfer model discretizes the wetland vegetation canopy into a three-dimensional voxel grid, each three-dimensional voxel grid is assigned with leaf area density, leaf inclination distribution and leaf optical property parameters, the Monte Carlo ray tracing method is used to simulate the propagation path of direct sunlight and sky scattered light in the vegetation canopy, the scattering probability is calculated according to the extinction coefficient when each light ray propagates in the three-dimensional voxel grid, the new propagation direction is determined according to the phase function when scattering occurs, the light ray experiences multiple reflections and transmissions between leaves, stems and water surfaces, the light energy escaping to the sensor direction is recorded, the energy distribution of a large number of light rays is counted to obtain the simulated canopy reflectance, the value range of leaf area index is set to 0.5 to 8.0, the value range of chlorophyll content is set to 20 to 80 , the value range of canopy height is set to 0.5 to 5.0m, 1000 groups of parameter combinations are generated in the parameter space by using Latin hypercube sampling, the three-dimensional radiation transfer model is run for each parameter combination to calculate the corresponding multi-spectral reflectance, and a mapping lookup table of parameter combinations and multi-spectral reflectance is established. The parameter search space is constrained by combining the prior information of canopy height extracted from the laser radar point cloud data, the canopy height model is extracted from the laser radar point cloud data, the measured value of the canopy height of each vegetation pixel is obtained, the parameter combinations with a deviation of more than 0.5m between the canopy height in the lookup table and the measured value of the canopy height are excluded, and only the measured reflectance matching is performed in the remaining parameter combinations. This method can reduce the parameter search space to 30% of the original size, significantly improve the inversion efficiency and accuracy, and compare the measured reflectance of the surface reflectance image with the simulated canopy reflectance value in the lookup table one by one during inversion, select the parameter combination with the smallest error as the inversion result, and obtain the leaf area index distribution and chlorophyll content distribution.

[0060] The specific implementation of step S04 is to determine whether the mixed pixel ratio in the wetland landscape segmentation result image is greater than 30%, and to count the class label uncertainty index of each pixel in the wetland landscape segmentation result image. When the class label uncertainty index is greater than the threshold value 0.3, it is determined that the mixed pixel, and the mixed pixel ratio is obtained by dividing the number of mixed pixels by the total number of pixels. The purpose of this step is to select a suitable abundance coefficient acquisition strategy according to the mixed pixel ratio. When the mixed pixel ratio is greater than 30%, the pure endmember spectrum is extracted from the high-resolution remote sensing image to establish an endmember spectrum library. The pure endmember refers to a pixel covered by a single ground object type, which is obtained by manual interpretation or pure pixel index screening. The endmember spectrum library includes the spectral curves of four types of endmembers: open water, emergent vegetation, submerged plants and bare mudflat. The linear spectral mixture decomposition is performed on the mixed pixels by using the full constraint least squares method. The linear spectral mixture model assumes that the reflectivity of the mixed pixel is equal to the linear combination of the reflectivity of each endmember and its abundance coefficient. Solving the abundance coefficient needs to satisfy the constraint condition that all abundance coefficients are non-negative and the sum is 1. The interior point method is used to solve the constraint optimization problem. The optimal abundance coefficient is iteratively calculated to minimize the mean square error of the simulated reflectivity and the measured reflectivity, and the abundance coefficient distribution of each ground object type is obtained. The Markov random field model is introduced for spatial neighborhood smoothing optimization. The neighborhood of a pixel is defined as its eight-connected neighboring pixels. An energy function is established, which includes a data fidelity term and a smoothing regularization term. The data fidelity term measures the deviation of the optimized abundance coefficient from the abundance coefficient distribution of each ground object type. The smoothing regularization term imposes a penalty on the difference between the abundance coefficients of adjacent pixels. The graph cut algorithm or the iterative conditional mode algorithm is used to minimize the energy function, and the optimized abundance coefficient distribution with enhanced spatial continuity is obtained. The Markov random field model eliminates isolated noise pixels by constraining the spatial context information, improving the reliability of the abundance estimation. When the mixed pixel ratio is less than or equal to 30%, it is considered that the class assignment result of the wetland landscape segmentation result image is basically reliable. The abundance of the dominant class label of each pixel is set to 1, and the abundance of the remaining classes is set to 0. The class purity distribution is formed as an alternative input for the abundance coefficient distribution of each ground object type. This simplified process can reduce the computational complexity when the mixed pixel ratio is low without significantly affecting the evaluation accuracy.

[0061] The specific implementation of step S05 is to input the chlorophyll concentration distribution, suspended matter concentration distribution, dissolved organic matter concentration distribution, leaf area index distribution, chlorophyll content distribution and optimized abundance coefficient distribution into the wetland health evaluation model. The purpose of this step is to comprehensively evaluate the health status of the wetland ecosystem and identify the degradation driving factors by using multi-source remote sensing inversion parameters. The water quality health index distribution is calculated by using the chlorophyll The concentration distribution, the suspended matter concentration distribution and the dissolved organic matter concentration distribution are compared with the water quality standard threshold value, the concentration is lower than the water quality standard threshold value, the exponential is close to 1, the concentration is much higher than the water quality standard threshold value, the exponential is close to 0, and the water quality standard threshold value is determined with reference to the surface water environmental quality standard. The vegetation health index distribution is calculated, the leaf area index distribution and the chlorophyll content distribution are compared with the health vegetation benchmark value, the parameter is close to the health vegetation benchmark value, the exponential is close to 1, and the parameter deviates from the health vegetation benchmark value significantly, the exponential decreases, and the health vegetation benchmark value is determined according to the typical vegetation growth state of the wetland type. The landscape health index distribution is calculated, the diversity index and the connectivity index of the optimized abundance coefficient distribution are calculated, the abundance distribution is uniform and the spatial connectivity is high, the exponential is higher, and the single type is dominant or the fragmentation is serious, the exponential decreases, the diversity index is calculated by using the Shannon diversity index or the Simpson diversity index, and the connectivity index is calculated by using the landscape connectivity analysis tool. The comprehensive health index distribution is calculated, the comprehensive health index distribution is equal to the weighted sum of the water quality health index distribution, the vegetation health index distribution and the landscape health index distribution, the weight coefficient is determined according to the wetland type and the ecological function, for the wetland with water quality purification as the main function, the weight of the water quality health index distribution is set to 0.5, the weight of the vegetation health index distribution is set to 0.3, and the weight of the landscape health index distribution is set to 0.2. According to the comprehensive health index distribution, the wetland health grade distribution is divided, when the comprehensive health index is greater than 0.8, the pixel is marked as a healthy state, when the comprehensive health index is between 0.6 and 0.8, the pixel is marked as a sub-healthy state, when the comprehensive health index is between 0.4 and 0.6, the pixel is marked as a degradation state, and when the comprehensive health index is less than 0.4, the pixel is marked as a serious degradation state. Whether there is a degraded state or a serious degradation state of the pixel area is judged, if there is, the contribution rate distribution of the water quality health index distribution, the vegetation health index distribution and the landscape health index distribution to the comprehensive health index is calculated, the contribution rate is equal to the product of each index and its weight coefficient divided by the theoretical maximum value of the comprehensive health index distribution, the index type with the largest contribution rate is identified as the dominant factor type, and the dominant factor type indicates the main ecological element leading to the degradation of the wetland, which provides a scientific basis for formulating targeted restoration measures. The wetland health grade distribution and the dominant factor type are output, the wetland health grade distribution presents the spatial distribution of different health states in the form of a chart, and the dominant factor type presents the dominant driving factor of each degradation area in the form of a statistical table or a spatial distribution map.

[0062] Specifically, the principle of the present application is that the principle that the present application can solve the technical problem lies in the organic combination of multi-source data collaborative mechanism and physical constraint inversion mechanism. Multi-spectral remote sensing images provide rich spectral features for identifying ground object types and inverting biochemical parameters, but lack vertical structure information. Laser radar point cloud data accurately describe the crown height distribution but do not contain spectral information, and the two are complementary in information dimension. The present application introduces the crown height extracted by the laser radar as a prior constraint into the three-dimensional radiation transfer model, reduces the originally high-dimensional ill-conditioned parameter search space to the physically feasible domain, and avoids the inversion process from falling into a local optimal solution. At the same time, the concentration range constraint, the non-negativity constraint of the absorption coefficient and other physical rules are embedded in the water quality parameter unmixing model, so that the neural network automatically avoids the solution that does not meet the optical characteristics of the water body in the training and reasoning process, and ensures the physical consistency of the inversion result. This physical-driven multi-source data fusion strategy ensures the accuracy and reliability of the parameter inversion in each dimension from the mechanism, provides high-quality input data for comprehensive health evaluation, and therefore can realize the collaborative and accurate inversion of the multi-dimensional health state parameters of the wetland ecosystem.

[0063] A specific embodiment 1 of the present application is provided below, and the specific implementation of each step in the embodiment 1 is described in detail as follows.

[0064] The specific implementation of step S01 is: first, the multi-spectral remote sensing image is radiometrically calibrated to convert the original digital quantization value into an apparent reflectance image. In the processing process of the dense dark object method, the pixels with a vegetation coverage greater than 80% and a normalized vegetation index greater than 0.7 are selected as the dark object candidate set in the apparent reflectance image, the apparent reflectance distribution of the dark object candidate set in the blue light band and the red light band is counted, and the average value of the 5% pixels with the lowest apparent reflectance in the apparent reflectance distribution is taken as the estimated value of the atmospheric path radiation. The calculation formula of the initial value of the aerosol optical depth is as follows:

[0065] ;

[0066] In the formula, is the initial value of the aerosol optical depth, dimensionless; is the estimated value of the atmospheric path radiation, unit: W / (m2·sr·nm); ; is the Rayleigh scattering radiation contribution, unit: W / (m2·sr·nm); ; is the single scattering albedo, dimensionless, the empirical value is 0.85 to 0.95; is the phase function, dimensionless; is the scattering angle, unit: degree; is the atmospheric path length, unit: km. The calculation formula of the phase function is as follows:

[0067] ;

[0068] wherein, is the cosine value of scattering angle, dimensionless.

[0069] The state transition equation in the process of Kalman filter dynamic updating the initial value of atmospheric parameters is as follows:

[0070] ;

[0071] wherein, is the aerosol optical depth at time , dimensionless; is the aerosol optical depth at time , dimensionless; is the process noise, dimensionless, subject to zero-mean Gaussian distribution, and the variance is empirically 0.01 to 0.05; is the time index. The calculation formula of the prior state estimation in the prediction step is as follows:

[0072] ;

[0073] wherein, is the prior state estimation at time , dimensionless. The calculation formula of the prior error covariance is as follows:

[0074] ;

[0075] wherein, is the prior error covariance at time , dimensionless; is the error covariance at time , dimensionless; is the process noise covariance matrix, dimensionless, and the empirical value is 0.001 to 0.01. The calculation formula of the Kalman gain in the update step is as follows:

[0076] ;

[0077] wherein, is the Kalman gain at time , dimensionless; is the measurement noise covariance, dimensionless, and the empirical value is 0.005 to 0.02. The calculation formula of the posterior state estimation is as follows:

[0078] ;

[0079] wherein, is the observation value at time , dimensionless, from the ground observation data of the sun photometer.

[0080] In the calculation process of the continuous absorption correction term of water vapor, the calculation formula of the water vapor optical thickness is expressed as follows:

[0081] ;

[0082] In the formula, is the water vapor optical thickness at the wavelength , dimensionless; is the atmospheric water vapor column density, with the unit of ; is the water vapor absorption cross section at the wavelength , with the unit of ; is the atmospheric mass factor, dimensionless, and is usually valued at 1 to 3; is the wavelength, with the unit of nm. The calculation formula of the water vapor absorption transmittance is expressed as follows:

[0083] ;

[0084] In the formula, is the water vapor absorption transmittance at the wavelength , dimensionless. The atmospheric correction equation is expressed as follows:

[0085] ;

[0086] In the formula, is the surface reflectivity at the wavelength , dimensionless; is the apparent radiance at the wavelength , with the unit of ; is the atmospheric path radiance at the wavelength , with the unit of ; is the solar irradiance at the wavelength , with the unit of ; is the solar zenith angle, with the unit of degree; is the Rayleigh scattering transmittance at the wavelength , dimensionless; is the aerosol scattering transmittance at the wavelength , dimensionless.

[0087] The specific implementation of step S02 is that the energy function of the wetland landscape segmentation algorithm constructed by the graph cut theory is expressed as follows:

[0088] ;

[0089] In the formula, is the energy function, dimensionless; Configure pixel labels; A set of pixels; For cell indexing; For pixels Tags; For pixels Assign tags The data items are dimensionless. The smoothing weighting coefficient is dimensionless and has an empirical value of 0.1 to 0.5. It is a set of adjacent cell pairs; For pixels Neighborhood cell index; For adjacent pixels and The smoothing term is dimensionless. The data terms are calculated using Mahalanobis distance, as expressed in the following formula:

[0090] ;

[0091] In the formula, For pixels spectral vector; For tags The center spectral vector of the corresponding category; For tags Covariance matrix of corresponding categories; superscript Indicates the transpose operation; superscript This represents the matrix inversion operation. The formula for calculating the smoothing term is as follows:

[0092] ;

[0093] In the formula, For indicator functions, dimensionless, when The value is 1 if it is true, and 0 otherwise. For pixels and The edge weights between them are dimensionless.

[0094] The formula for calculating the normalized spectral distance is as follows:

[0095] ;

[0096] In the formula, For pixels and The Euclidean distance of the spectral vector is dimensionless; For pixels The spectral vector magnitude is dimensionless; For pixels The spectral vector magnitude is dimensionless. The Boykov-Kolmogorov maximum flow algorithm calculates the minimum cut set through augmented path search and residual graph updates. The pixel label configuration corresponding to the minimum cut set is the global minimum solution of the energy function. The alpha expansion algorithm decomposes the multi-label problem into a series of binary label expansion subproblems. In each iteration, a label is selected as the alpha label to attempt to expand its region. The graph cut is used to solve for whether to accept or reject the expansion. The iteration continues until convergence to obtain the wetland landscape segmentation image.

[0097] The specific implementation method of step S03 is as follows: During the principal component analysis dimensionality reduction process, the formula for calculating the covariance matrix is ​​expressed as follows:

[0098] ;

[0099] In the formula, It is the covariance matrix; The number of pixels in the open water area; For the first The multispectral reflectance vector of each pixel is dimensionless. is the mean vector of the multispectral reflectance vector, and is dimensionless; Let be the cell index. Solve for the eigenvalues ​​and eigenvectors of the covariance matrix. Sort the eigenvalues ​​from largest to smallest and select the top three eigenvectors with a cumulative variance contribution rate of 95%. The formula for calculating the principal component eigenvectors is as follows:

[0100] ;

[0101] In the formula, For the first Each principal component eigenvalue is dimensionless. This refers to the number of multispectral bands; For the first The th eigenvector of the th feature vector One component, dimensionless; For the first The reflectivity of each band is dimensionless. Principal component number; This refers to the band number. Ultraviolet reflectance refers to the surface reflectance of the ultraviolet light band, with the center wavelength of the ultraviolet band between 350 nm and 400 nm. The formula for calculating the fluorescence linear height index is as follows:

[0102] ;

[0103] In the formula, The fluorescence linear height index is dimensionless; The reflectivity in the 685nm band is dimensionless. The reflectivity in the 665nm band is dimensionless. The reflectance is in the 705nm band and is dimensionless.

[0104] The water quality parameter unmixing model is a spectral unmixing model based on an improved Gaussian kernel. The expression for the improved Gaussian kernel is:

[0105] ;

[0106] In the formula, To improve the Gaussian kernel function value, it is dimensionless; The input feature vector is dimensionless. It is centered on the core and is dimensionless. The kernel width parameter is dimensionless and has an empirical value of 0.5 to 2.0. For the first The constraint function for each water quality parameter is dimensionless. For the first Predicted values ​​for each water quality parameter; For the first Minimum reasonable values ​​for each water quality parameter; For the first The maximum reasonable value of each water quality parameter. The expression for the constraint function is:

[0107] ;

[0108] In the formula, The penalty weighting coefficient is dimensionless and has an empirical value of 10. The extinction coefficient is calculated using the following formula when simulating multiple scattering processes from the vegetation canopy in a three-dimensional radiative transfer model:

[0109] ;

[0110] In the formula, wavelength The extinction coefficient at the point, in units of ; Leaf area density, in units of ; wavelength Absorption coefficient at , in units of Dimensionless; wavelength The scattering coefficient at the location, in units of Dimensionless. In the lookup table construction method, the mapping relationship between parameter combinations and multispectral reflectance is obtained through a three-dimensional radiative transfer model. The error calculation formula between the measured reflectance and the simulated canopy reflectance value in the lookup table during inversion is expressed as follows:

[0111] ;

[0112] In the formula, The error is dimensionless. The reflectance is dimensionless and is used to find the simulated reflectance of the corresponding parameter combination in the table. Leaf area index, dimensionless; Chlorophyll content, unit: ; The height is the canopy height, in meters (m).

[0113] The specific implementation of step S04 is as follows: In the process of calculating the mixed pixel ratio, the formula for calculating the category label uncertainty index is expressed as follows:

[0114] ;

[0115] In the formula, For pixels The category label uncertainty index is dimensionless. For pixels Belongs to the label The probability is dimensionless; For tag indexing. The mathematical expression of the fully constrained least squares method for linear spectral mixture decomposition of mixed pixels is:

[0116] ;

[0117] In the formula, For mixed pixels The reflectivity vector is dimensionless; The number of endmembers; For pixels The Middle The abundance coefficient of each endmember is dimensionless. For the first The reflectivity vector of each endmember is dimensionless. Let be the error vector, which is dimensionless. The constraints are stated as follows:

[0118] ;

[0119] The objective function is expressed as follows:

[0120] .

[0121] The energy function of the Markov random field model for spatial neighborhood smoothing optimization of the abundance coefficient distribution of various land cover types is expressed as follows:

[0122] ;

[0123] In the formula, It is an energy function, dimensionless; To optimize the abundance coefficient distribution; For pixels The optimized abundance coefficient vector is dimensionless. For pixels The initial abundance coefficient vector is dimensionless; The smoothing regularization coefficient is dimensionless and has an empirical value of 0.05 to 0.2.

[0124] The specific implementation method of step S05 is as follows: The calculation formula for the comprehensive health index distribution is expressed as follows:

[0125] ;

[0126] In the formula, A comprehensive health index, dimensionless; This is a dimensionless sub-index for water quality health. A dimensionless sub-index for vegetation health. The landscape health sub-index is dimensionless. The weighting coefficient is dimensionless and is used for wetlands whose primary function is water purification. The default value is 0.5. The default value is 0.3. The default value is 0.2. The formula for calculating the water quality health index is as follows:

[0127] ;

[0128] In the formula, For the first Concentration values ​​of each water quality parameter; For the first Standard thresholds for each water quality parameter; For the first The standard deviation of each water quality parameter is empirically defined as 10% to 30% of the standard threshold. The formula for calculating the vegetation health sub-index is as follows:

[0129] ;

[0130] In the formula, The leaf area index is a baseline for healthy vegetation, dimensionless, with an empirical value of 3 to 5. Chlorophyll content as a baseline for healthy vegetation, in units of The experience value is 40 to 60; The standard deviation of leaf area index is dimensionless, with an empirical value of 0.5 to 1.5. The standard deviation of chlorophyll content is given in units of 1000 ppm. The empirical value is 5 to 15. The formula for calculating the landscape health index is as follows:

[0131] ;

[0132] In the formula, The normalized Shannon diversity index is dimensionless. The normalized connectivity index is dimensionless. The formula for calculating the normalized Shannon diversity index is as follows:

[0133] ;

[0134] In the formula, For the first Abundance coefficients of individual land cover types, dimensionless; It is a function of the natural logarithm. The formula for calculating the normalized connectivity index is as follows:

[0135] ;

[0136] In the formula, For the first The area proportion of each land cover type is dimensionless. For the first Connectivity index for each land feature type, dimensionless; For the first The maximum connectivity index for each land cover type, dimensionless, with an empirical value ranging from 10 to 100. The formula for calculating the contribution rate is as follows:

[0137] ;

[0138] In the formula, For the first The contribution rate of each sub-index is dimensionless. For the first Individual index values, dimensionless; For the first The weight coefficients corresponding to each sub-index are dimensionless; the ecological element corresponding to the sub-index type with the largest contribution rate is the dominant factor type.

[0139] To better understand and implement this invention, the following is a specific application scenario of the invention, Example 2: Technicians selected a typical lake wetland area as the research object; the area of ​​this area is approximately 125... The area encompasses various wetland landscape types, including open water bodies, emergent vegetation, submerged plants, and exposed mudflats. Technicians acquired multispectral remote sensing images and lidar point cloud data for the region. The multispectral remote sensing images contain eight bands, covering the ultraviolet band from 350nm to the near-infrared band from 950nm, with a spatial resolution of 10m. The lidar point cloud data has a point density of 4 points per square meter and an elevation accuracy of 0.15m.

[0140] In the implementation process of step S01, the technical personnel perform radiation calibration on the multispectral remote sensing image, convert the digital quantitative value into apparent reflectance image, use the dense dark object method to screen out the pixels with vegetation coverage greater than 80% and normalized difference vegetation index greater than 0.7 in the image as the dark object candidate set, a total of 18562 candidate pixels are obtained, the statistical distribution of the apparent reflectance of these pixels in the blue light band and the red light band is calculated, the average value of the 5% pixels with the lowest apparent reflectance is taken as the estimated value of the atmospheric path radiation, which is 0.0285 in the blue light band and 0.0167 in the red light band, and the initial value of the aerosol optical depth is calculated as 0.326 according to the linear relationship between the atmospheric path radiation and the aerosol optical depth. According to the synchronous acquisition of the ground observation data of the sun photometer, the technical personnel dynamically update the aerosol optical depth by using Kalman filtering, take the initial value of the aerosol optical depth as the system state variable, set the process noise covariance as 0.01 and the measurement noise covariance as 0.005, and obtain the time series of the regional aerosol optical depth through the iterative calculation of the prediction step and the update step. The optimal value of the aerosol optical depth at the current time extracted from the time series is 0.312, which is taken as the optimal atmospheric correction parameter. According to the atmospheric water vapor column concentration 2.8 The water vapor optical thickness of each band is calculated according to the water vapor absorption cross section, and the water vapor optical thickness reaches the maximum value 0.156 at the near-infrared wavelength 940nm. The water vapor absorption transmittance is calculated by using the exponential decay function, and the water vapor absorption transmittance is multiplied by the Rayleigh scattering transmittance 0.912 and the aerosol scattering transmittance 0.785 to obtain the atmospheric correction equation, and the conversion from the apparent reflectance image to the surface reflectance image is completed, as shown in Figure 1 The reflectance range of the converted surface reflectance image in the visible light band is 0.02 to 0.35, and the reflectance range in the near-infrared band is 0.15 to 0.68.

[0141] In the implementation process of step S02, the technical personnel take the 125000 pixels of the surface reflectance image as the nodes of the undirected graph based on the graph cut theory, calculate the normalized spectral distance between adjacent pixels as the edge weight, extract the reflectance vectors of adjacent pixels in 8 spectral bands, calculate the Euclidean distance of the two reflectance vectors and divide the maximum modulus of the reflectance vector to obtain the normalized spectral distance, and the statistical distribution of the normalized spectral distance is shown in Table 1.

[0142] Table 1 Statistical distribution table of normalized spectral distance

[0143] The skilled person sets four terminal nodes of open water body, emergent vegetation, submerged plant and bare mudflat by artificial interpretation and prior knowledge, selects 50 seed pixels for each terminal node, constructs the energy function of graph cut theory, quantizes the similarity between the spectral vector of each pixel and the center spectral vector of each class by Mahalanobis distance, and the penalty weight of smoothing term is inversely proportional to the edge weight, and the weight coefficient is set to 5.0. The Boykov-Kolmogorov maximum flow algorithm is used to solve the graph cut problem, and the minimum cut set is calculated by augmented path search and residual graph update iteration. After 23 iterations, the algorithm converges, and the alpha expansion algorithm is used for multi-label segmentation of the four terminal nodes. Each iteration selects a label as the alpha label to try to expand its area. After 16 iterations, the wetland landscape segmentation result image is obtained. The segmentation result image shows that the area ratio of open water body is 56.8%, the area ratio of emergent vegetation is 28.4%, the area ratio of submerged plant is 9.6%, and the area ratio of bare mudflat is 5.2%, as shown in FIG. 8. The segmentation result image not only ensures spectral consistency but also has good spatial continuity. Figure 2

[0144] In the implementation process of step S03, the skilled person extracts the open water body region from the wetland landscape segmentation result image, and obtains 71,000 water body pixels. The multi-spectral band reflectance data of these pixels is extracted, and the multi-spectral band reflectance data is dimensionally reduced by principal component analysis. The covariance matrix is calculated and the eigenvalues and eigenvectors are solved. The cumulative variance contribution rate of the first three eigenvectors is 96.8%. The multi-spectral band reflectance data is projected onto the first three eigenvectors to obtain the principal component eigenvector. The average reflectance at the center wavelength of 375 nm is 0.048. The fluorescence linear height index is extracted. The average difference between the reflectance at the fluorescence peak wavelength of 685 nm and the reflectance baseline at the adjacent wavelengths of 665 nm and 705 nm is 0.0125. The skilled person inputs the principal component eigenvector, the ultraviolet band reflectance and the fluorescence linear height index into the water quality parameter unmixing model. The model adopts a five-layer feedforward architecture with an improved Gaussian kernel. The input layer receives five input variables. The first hidden layer contains 128 neurons. The second hidden layer contains 64 neurons. The third hidden layer contains 32 neurons. The output layer contains three neurons. The model training data set contains 6,850 training samples, which are divided into 5,480 training sets and 1,370 validation sets in the ratio of 8:2. The model is trained using the Adam optimization algorithm. The initial learning rate is set to 0.001. The batch size is 64. The total training period is 200. The loss function uses the weighted sum of mean square error and physical constraint violation penalty term. The penalty weight coefficient is set to 10. After 182 training periods, the validation set loss does not decrease for 10 consecutive periods, and the training is terminated early. The model weight with the minimum validation set loss is selected as the final model parameter. The model is applied to the open water body region to obtain the chlorophyll-a concentration​ Concentration distribution, suspended matter concentration distribution and dissolved organic matter concentration distribution, chlorophyll ranged from 8.5 to 65.3 , with an average of 28.6 , suspended matter concentration ranged from 12.3 to 82.7 , with an average of 38.4 , dissolved organic matter concentration ranged from 3.2 to 18.9 , with an average of 9.7 As shown in Figure 3 , the water quality parameters showed obvious heterogeneity in spatial distribution, with significantly higher suspended matter concentration and dissolved organic matter concentration in the nearshore area and the lake inlet area than in the open water area.

[0145] The technical personnel extracted the vegetation area from the wetland landscape segmentation result image, a total of 47500 vegetation pixels, established a lookup table of surface reflectance and leaf area index, chlorophyll content using a three-dimensional radiation transfer model, set the value range of leaf area index to 0.5 to 8.0, step size 0.5, value range of chlorophyll content to 20 to 80 , step size 5 , value range of canopy height to 0.5 to 5.0m, step size 0.5m, generated 1000 parameter combinations in the parameter space using Latin hypercube sampling, ran the three-dimensional radiation transfer model for each parameter combination, simulated the vegetation canopy multiple scattering process, tracked 100000 light rays each time, recorded the light energy finally escaping to the sensor direction, statistically obtained the canopy reflectance simulation value, established the mapping lookup table of parameter combination and multi-spectral reflectance. The technical personnel extracted the canopy height model from the laser radar point cloud data, obtained the measured value of canopy height for each vegetation pixel, the measured value of canopy height ranged from 0.8 to 4.6m, removed the parameter combination in the lookup table with a deviation of more than 0.5m between canopy height and measured value of canopy height, the parameter search space was reduced from 1000 to 312, the reduction ratio was 68.8%, compared the measured reflectance of the surface reflectance image with the canopy reflectance simulation value in the lookup table one by one, selected the parameter combination with the smallest root mean square error as the inversion result, obtained the leaf area index distribution and chlorophyll content distribution, the leaf area index ranged from 1.2 to 6.8, with an average of 3.9, the chlorophyll content ranged from 28 to 72 , with an average of 48 As shown in Figure 4 , the leaf area index and chlorophyll content of emergent vegetation area were significantly higher than those of submerged plant area, reflecting the growth difference of different vegetation types.

[0146] In the implementation process of step S04, the technician counts the class label uncertainty index of each pixel in the wetland landscape segmentation result image. When the class label uncertainty index is greater than the threshold value 0.3, it is determined to be a mixed pixel, and a total of 38750 mixed pixels are identified, with a mixed pixel ratio of 31.0%, which exceeds the threshold value of 30%. Therefore, linear spectral mixture decomposition needs to be performed. The technician extracts pure endmember spectra from high-resolution remote sensing images to establish an endmember spectral library. The average reflectance of the open water endmember is 0.058 in the visible light band, the average reflectance of the emergent vegetation endmember is 0.524 in the near-infrared band, the average reflectance of the submerged plant endmember is 0.186 in the green light band, and the average reflectance of the bare mudflat endmember is 0.312 in the red light band. The full-constraint least squares method is used to perform linear spectral mixture decomposition on the mixed pixels, and the abundance coefficients need to satisfy the constraint conditions that all abundance coefficients are non-negative and the sum is 1. The interior point method is used to iteratively calculate the optimal abundance coefficients, with an average iteration number of 15 times, so that the root mean square error of the simulated reflectance and the measured reflectance is reduced to 0.018, and the abundance coefficient distribution of each land cover type is obtained. The technician introduces the Markov random field model for spatial neighborhood smoothing optimization, defines the neighborhood of a pixel as its eight-connected neighboring pixels, establishes an energy function, sets the weight of the data fidelity term to 0.7, and sets the weight of the smoothing regularization term to 0.3. The graph cut algorithm is used to minimize the energy function, and after 28 iterations, the spatial continuity enhanced optimized abundance coefficient distribution is obtained. The spatial continuity index of the optimized abundance coefficient distribution is improved from 0.625 to 0.812, eliminating the influence of isolated noise pixels.

[0147] In the implementation process of step S05, the technician inputs the chlorophyll concentration distribution, suspended matter concentration distribution, dissolved organic matter concentration distribution, leaf area index distribution, chlorophyll content distribution, and optimized abundance coefficient distribution into the wetland health evaluation model to calculate the water quality health score index distribution. According to the surface water environmental quality standard, the chlorophyll concentration threshold is set to 40 , the suspended matter concentration threshold is set to 50 , and the dissolved organic matter concentration threshold is set to 15 . The water quality health score index distribution is calculated by comparing the concentration with the threshold value, and the average value is 0.72. The vegetation health score index distribution is calculated, and the health vegetation reference value is set according to the typical vegetation growth state of the wetland type. The leaf area index reference value is 4.5, and the chlorophyll content reference value is 55 The vegetation health score index distribution is calculated by comparing the parameters with the reference values, and the average value is 0.68. The landscape health score index distribution is calculated by optimizing the Shannon diversity index and the connectivity index of the abundance coefficient distribution. The Shannon diversity index is 1.24, and the connectivity index is 0.78. The landscape health score index distribution is obtained, and the average value is 0.75. The skilled person calculates the comprehensive health index distribution. Since the wetland is mainly used for water purification, the weight of the water quality health score index distribution is set to 0.5, the weight of the vegetation health score index distribution is set to 0.3, and the weight of the landscape health score index distribution is set to 0.2. The comprehensive health index distribution is calculated, and the average value is 0.71. According to the comprehensive health index distribution, the health grade distribution of the wetland is divided, and the statistical results are shown in Table 2.

[0148] Table 2 Statistical table of wetland health grade distribution

[0149] The skilled person found that there were pixel regions in the degradation state and the serious degradation state, and the total area ratio was 24.8%. The contribution rate distribution of the water quality health score index distribution, the vegetation health score index distribution and the landscape health score index distribution to the decrease of the comprehensive health index was calculated. In the degradation area, the contribution rate of the water quality health score index was 58.6%, the contribution rate of the vegetation health score index was 26.3%, and the contribution rate of the landscape health score index was 15.1%. The water quality health score index was identified as the main factor type, indicating that water pollution was the main driving factor leading to the degradation of the wetland. Figure 5 As shown in Table 2, the degradation area mainly concentrated in the nearshore area and the lake inlet area. These areas were affected by agricultural non-point source pollution and urban domestic sewage, leading to excessive water quality parameters, and further affecting the overall health status of the wetland.

[0150] The technical personnel output wetland health grade distribution and dominant factor type, provide scientific basis for wetland ecological protection and restoration, suggest to strengthen pollution source control and water body ecological restoration in view of water quality dominant type degradation area, suggest to implement vegetation restoration and habitat reconstruction in view of vegetation dominant type degradation area, suggest to optimize landscape pattern and improve spatial connectivity in view of landscape dominant type degradation area.The technical progress brought by the present application relative to the traditional wetland health evaluation method is reflected in many aspects.Firstly, dynamic atmospheric correction is realized by fusing remote sensing inversion and ground observation data through Kalman filtering, the precision limitation of the traditional single data source method in the case of uneven aerosol distribution is overcome, high-quality basic data are provided for parameter inversion, secondly, accurate unmixing of water quality parameters is realized by embedding an improved Gaussian kernel with physical constraints in the neural network, the confusion problem caused by spectral overlap in the traditional band ratio method is avoided, and the parameter inversion precision in complex water body environment is improved, thirdly, the parameter search space of the three-dimensional radiation transfer model is constrained by laser radar point cloud data, the multi-solution problem caused by parameter compensation effect in the traditional unconstrained inversion method is solved, and finally, the comprehensive evaluation of wetland health condition is realized by multi-source parameters in three dimensions of water quality, vegetation and landscape, compared with the traditional evaluation method based on single parameter or single data source, the present application can more accurately diagnose the wetland degradation state and identify the dominant driving factor, and provide scientific support for formulating targeted wetland protection and restoration measures.

[0151] The above is merely a specific implementation of the present application, but the protection scope of the present application is not limited thereto, any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered within the protection scope of the present application.

Claims

1. A method for wetland ecosystem health assessment of multi-source remote sensing data, characterized in that, The multispectral remote sensing image and the laser radar point cloud data of a target wetland area are acquired, the surface reflectance image is obtained by radiation calibration and atmospheric correction of the multispectral remote sensing image, the surface reflectance image is segmented into four landscape types of open water, emergent vegetation, submerged plant and bare mudflat based on the graph cut theory to obtain a wetland landscape segmentation result image, the open water area and the vegetation area are extracted from the wetland landscape segmentation result image, and the chlorophyll-a concentration distribution, the suspended matter concentration distribution and the dissolved organic matter concentration distribution are obtained by inputting the water quality parameter unmixing model, the leaf area index distribution and the chlorophyll content distribution are obtained by inversion through a three-dimensional radiation transfer model combined with the laser radar point cloud data, the proportion of mixed pixels in the wetland landscape segmentation result image is judged, the optimized abundance coefficient distribution is obtained by using the total least squares method and the Markov random field model, the comprehensive health index distribution, the water quality health index distribution, the vegetation health index distribution and the landscape health index distribution are calculated by inputting the wetland health evaluation model, the wetland health grade distribution is divided, and the dominant factor type is identified.

2. The method of claim 1, wherein, The atmospheric correction adopts the dense dark pixel method to extract initial values of atmospheric parameters, and the optimal atmospheric correction parameters are obtained by dynamic updating based on the sun photometer ground observation data, and the surface reflectance image is obtained by introducing a water vapor continuous absorption correction term.

3. The method of claim 2, wherein, The dense dark pixel method selects pixels with a vegetation coverage greater than 80% and a normalized vegetation index greater than 0.7 as a dark pixel candidate set, and the apparent reflectance distribution of the dark pixel candidate set in the blue and red bands is counted, and the average value of the 5% pixels with the lowest apparent reflectance is taken as the estimated value of the atmospheric path radiation to inversely calculate the initial value of the aerosol optical depth.

4. The method of claim 3, wherein, The Kalman filter takes the initial value of the aerosol optical depth as a system state variable, establishes a state transition equation, takes the initial value of the aerosol optical depth obtained by the dense dark pixel method and the observation value of the sun photometer ground observation data as measurement variables, and iteratively obtains the time series of regional aerosol optical depth through the prediction step and the update step.

5. The method of claim 4, wherein, The water vapor continuous absorption correction term calculates the water vapor optical depth of each band according to the atmospheric water vapor column concentration and the water vapor absorption cross section, uses an exponential decay function to describe the influence of water vapor continuous absorption on radiation transfer to obtain the water vapor absorption transmittance, and multiplies the water vapor absorption transmittance with the Rayleigh scattering transmittance and the aerosol scattering transmittance to obtain the atmospheric correction equation.

6. The method of claim 5, wherein, The graph cut theory takes the pixels of the surface reflectance image as undirected graph nodes, calculates the normalized spectral distance between adjacent pixels as the edge weight, sets four types of terminal nodes, solves the minimum value of the energy function by the minimum cut maximum flow algorithm, and decomposes the multi-label problem into a series of binary label expansion sub-problems by the alpha expansion algorithm until convergence.

7. The method of claim 6, wherein, The normalized spectral distance extracts the reflectance vectors of adjacent pixels in all spectral bands, calculates the Euclidean distance of the two reflectance vectors, and obtains the normalized spectral distance as the edge weight by dividing the Euclidean distance by the maximum modulus of the reflectance vector.

8. The method of claim 7, wherein, The water quality parameter unmixing model is a spectral unmixing model based on an improved Gaussian kernel. The improved Gaussian kernel introduces four types of physical constraints, i.e., water quality concentration range constraint, spectral physics constraint, absorption peak position constraint and scattering angle distribution constraint. Constraint terms are added to the traditional Gaussian kernel, and the constraint terms are integrated into the kernel function optimization objective through the Lagrange multiplier method.

9. The method of claim 8, wherein, The water quality parameter unmixing model is a five-layer feedforward architecture. The input layer receives principal component feature vectors, ultraviolet band reflectivity and fluorescence linear height index. The first hidden layer contains 128 neurons and uses an improved Gaussian kernel as an activation function. The second hidden layer contains 64 neurons. The third hidden layer contains 32 neurons. The output layer contains 3 neurons.

10. The method of claim 9, wherein, Before inputting the water quality parameter unmixing model, the multispectral band reflectivity data of the open water area is subjected to principal component analysis dimension reduction. The eigenvalues and eigenvectors of the covariance matrix are calculated. The first three eigenvectors with a cumulative variance contribution rate of 95% are selected according to the descending order of the eigenvalues to obtain the principal component feature vectors.

Images