A wetland information extraction method based on multi-source remote sensing data

Abstract

The present application relates to a kind of wetland information extraction method based on multi-source remote sensing data, first to image pre-processing, using optical image extraction spectral feature, index feature, kernel index feature, principal component segmentation feature, texture feature, using SAR image extraction SAR feature, texture feature, using DEM data extraction elevation, slope feature, and all features are overlaid with multi-temporal phase;Then based on Mahalanobis distance measure interclass separability, to carry out feature sensitivity experiment and feature optimization experiment, obtain optimal feature set;Finally using the obtained feature set, using integrated learning classification method, reach optimal classification effect.The present application designs method to integrate multiple types of data, truly and comprehensively express ground information;Multiple time phase data are overlaid, retain time sequence characteristics;Introduce principal component segmentation feature, fully obtain the geometric shape information of ground;Using the feature optimization method based on Mahalanobis distance, effectively retain the feature.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image classification technology, specifically to a method for extracting wetland information based on multi-source remote sensing data. Background Technology

[0002] Wetlands, known as the "kidneys of the Earth," provide habitat for one-third of the world's endangered species. They not only reduce flood risk, protect coastlines, conserve soil and water, filter sediments, remove pollution, and purify water sources, but also possess aesthetic and recreational value and serve as an important indicator of environmental health. Furthermore, wetlands foster fisheries and animal husbandry, bringing countless natural benefits to the local economy.

[0003] Although wetlands share many similar characteristics, they exhibit considerable variability in size, location, and hydrology. Furthermore, they often constitute transitional zones between clearly defined terrestrial and aquatic areas. This fact makes defining these natural resources difficult and explains the numerous different definitions of "wetland."

[0004] my country boasts 180 million meters of mainland coastline and 140 million meters of island coastline, spanning temperate and subtropical zones, and possesses abundant wetland resources. These resources are distributed across various regions including: Northeast wetlands, the middle and lower reaches of the Yellow River, the middle and lower reaches of the Yangtze River, the northern coast of Hangzhou Bay, areas south of Hangzhou Bay, the Yunnan-Guizhou Plateau, the semi-arid wetlands of Inner Mongolia and Xinjiang, and the alpine wetlands of the Qinghai-Tibet Plateau. Sixty-four wetlands within my country are listed as Wetlands of International Importance, such as the Chinese Sturgeon Nature Reserve in the Yangtze River Estuary of Shanghai and the Zhalong National Nature Reserve in Heilongjiang Province. However, relative to sea-level rise, rising temperatures, changes in precipitation, storms, and other climate changes, as well as human activities such as land irrigation, groundwater extraction, and drainage leading to drought, salinization, eutrophication, pollution, and overuse, wetland protection faces severe challenges. Therefore, employing new technologies for rapid wetland monitoring is of great significance for wetland conservation.

[0005] As early as the 1950s, Chinese scholars began conducting field surveys and resource inventories of wetlands. While traditional field surveys may be the most accurate and reliable method, many wetlands are located in remote areas with vegetation, topography, and hydrological conditions unfavorable for on-site investigation. Therefore, this method is time-consuming and costly. Furthermore, on-site investigations often cannot achieve large-scale, real-time statistics on wetland conditions, making it easier to miss some areas. The emergence of remote sensing technology has effectively solved the time and cost problems associated with traditional wetland inventory methods.

[0006] In the late 1970s, the creation of wetland land cover and land use change maps based on remote sensing technology emerged in my country. In the mid-to-late 1990s, research on remote sensing wetlands gradually expanded. By the 21st century, the technology for remote sensing wetland detection had matured. The achievements of remote sensing in large-scale wetland monitoring, wetland biodiversity monitoring, and biomass estimation have provided a foundation and ideas for future research.

[0007] Compared to traditional field surveys, remote sensing technology can solve the problems of complex environments and high costs, enabling large-scale, short-term repeated measurements, providing information on changes at different scales around wetlands, and can be easily integrated into geographic information systems for multi-faceted analysis. In addition, the combination of many remote sensing data types also brings benefits to high-precision wetland mapping.

[0008] However, some characteristics of wetlands and existing classification systems pose challenges to remote sensing in wetland information extraction. (1) Some wetlands exhibit high intraclass variability and low interclass variability, resulting in spectral information and backscattering information similar to those of many wetland categories. Furthermore, the same wetland type has different names in different classification systems, making classification difficult. (2) Water levels around wetland plants change seasonally and are also affected by snowmelt, precipitation, and human activities. (3) Wetland boundaries are indistinct, transition zones are narrow, and overlap with other wetlands often leads to insufficient training samples. (4) Detection of small wetlands requires high-resolution remote sensing imagery. At lower resolutions, small wetlands are often mistaken for other land features and difficult to detect.

[0009] Based on the characteristics of rapid wetland evolution and previous research findings, it is clear that timely and high-precision monitoring of wetland changes is necessary. However, most current research focuses on exploring high-precision algorithms, with little attention paid to long-term spatiotemporal dynamic monitoring, which is detrimental to the study of the driving mechanisms of wetland change. Furthermore, in the exploration of wetland evolution, most studies lack comprehensive analysis that is closely linked to real-world applications. Summary of the Invention

[0010] The purpose of this invention is to provide a wetland information extraction method based on multi-source remote sensing data. It employs a novel design strategy and exhibits superior performance in terms of computational efficiency, information extraction accuracy, and robustness to noise.

[0011] To achieve the above functions, this invention designs a wetland information extraction method based on multi-source remote sensing data. For the target area, multi-temporal remote sensing images of the target area containing various preset categories of land features are collected. For the multi-temporal remote sensing images, steps A-F are executed to complete the classification of each pixel in the multi-temporal remote sensing images, thereby achieving land feature classification of the target area.

[0012] Step A: Based on the phenological characteristics and distribution features of each land use type in the wetland classification system, classify the land features in the target area, collect multi-temporal remote sensing images of the target area containing each preset category of land features, including optical multispectral images, SAR remote sensing images, and DEM data, and then proceed to Step B.

[0013] Step B: Preprocess the multispectral image and SAR remote sensing image respectively, and fuse the preprocessed multispectral image, SAR remote sensing image and DEM data to obtain a three-temporal median composite image, i.e., overlaid image X, and then proceed to step C.

[0014] Step C: For the optical band portion of the superimposed image X, extract multi-temporal index features, which include sub-features NDVI, MNDWI, and NDBBI.

[0015] Extract the texture features of the superimposed image X, including the texture features of the optical band portion of the superimposed image X, which includes sub-features ASM, contrast, entropy, and correlation; and the texture features of the microwave band portion of the superimposed image X, which includes sub-features entropy, intertia, gearys, gradient, direction, and edgy.

[0016] For the optical band portion of the overlay image X, kernel index features are calculated, including sub-features kNDVI, kNDBBI, and kMNDWI. Principal component analysis is used to obtain the principal components of the optical band portion of the overlay image X. Object-oriented segmentation is performed on the first principal component to obtain principal component segmentation features. At this point, the extraction of all types of features is completed, resulting in a feature set W1. The overlay image X is superimposed with the feature set W1 to obtain the feature fusion image Y, and then the process proceeds to step D.

[0017] Step D: For the feature fusion image Y, calculate the Mahalanobis distance value of the sample points of each preset category of land cover in the feature fusion image Y to represent the separability between different preset category combinations of land cover; then proceed to step E; in the following steps, Y is replaced with other images, and the Mahalanobis distance value in other cases is calculated using the same principle;

[0018] Step E: Combine the optical band portion of the superimposed image X with multi-temporal index features, kernel index features, texture features of the optical band portion, and texture features of the microwave band portion to obtain feature maps after each combination. Use a gradient boosting tree classifier to classify each pixel in each feature map after each combination, calculate the overall classification accuracy and Mahalanobis distance value, and define the combined features whose classification accuracy and Mahalanobis distance value are both higher than the preset threshold as sensitive features.

[0019] The feature fusion image Y obtained in step C is fed into a gradient boosting tree classifier for training. Simultaneously, the explain algorithm in the GEE platform is used to calculate the feature importance index of all sub-features in the feature fusion image Y. Then, the optical band portion of the overlay image X is filtered out, and all remaining sub-features are sorted in descending order of feature importance index. Based on the optical portion of the overlay image X, new features are added one by one in descending order. After each feature addition, the overall classification accuracy and Mahalanobis distance are calculated, resulting in curves showing the changes in overall classification accuracy and Mahalanobis distance with the number of features. By analyzing the trends of the curves, the first N features are determined as necessary features. Finally, it is checked whether any sensitive features have not been selected among the remaining features; if so, they are added to obtain the final feature set W2, and then the process proceeds to step F.

[0020] Step F: Input the final feature set W2 into the gradient boosting tree classifier to classify each pixel in the final feature set W2, obtain the preset category to which each pixel belongs, and thus realize the land cover classification of the target area.

[0021] As a preferred technical solution of the present invention: the preprocessing of multispectral images in step B includes cloud removal processing, which removes areas with cloud cover greater than 5%, and performs atmospheric correction, date filtering, and cropping according to the vector boundary of the study area; the preprocessing of SAR remote sensing images includes atmospheric correction, cropping according to the vector boundary of the study area, and fusing upper and lower orbit data, and performing date filtering and pattern filtering on the fused SAR remote sensing images.

[0022] As a preferred embodiment of the present invention, the specific steps of step C are as follows:

[0023] Step C1: For the optical band portion of the superimposed image X, calculate the sub-features NDVI, MNDWI, and NDBBI of the multi-temporal exponential features as follows:

[0024]

[0025] In the formula, NIR is the near-infrared band in the optical band portion of the superimposed image X, Red is the red light band, Green is the green light band, MIR is the mid-infrared band, and SWIR2 is the short-wave infrared band.

[0026] Step C2: Extract the sub-features ASM, contrast, entropy, and correlation of the texture features in the optical band portion of the overlay image X, and the sub-features entropy, intertia, gearys, gradient, direction, and edgy of the texture features in the microwave band portion of the overlay image X, as follows:

[0027] Let the pixel value of the pixel at positions (x,y) and (x+d,y+l) in the overlay image X be (i,j), the number of pixels with pixel value (i,j) in the overlay image X be N(i,j), and the total number of pixels in the overlay image X be N. glcm Then the probability P(i,j) of pixel value (i,j) appearing in the superimposed image X is:

[0028]

[0029] The ASM expression for the energy feature of pixel value (i,j) is:

[0030]

[0031] The contrast feature expression for pixel value (i,j) is:

[0032]

[0033] The entropy expression for pixel value (i,j) is:

[0034]

[0035] The autocorrelation expression for pixel value (i,j) is:

[0036]

[0037] In the formula,

[0038] Set the matrix size to 6, and use the GEE built-in algorithm glcmTexture to calculate the number of times a pixel with pixel value X is adjacent to a pixel with pixel value Y in a preset direction and distance. Count the number of adjacent pixels to obtain the feature intertia.

[0039] Define a 9×9 kernel with a focal position of (-4, -4). Its horizontal and vertical coordinates represent the offsets from the left and top of the kernel, respectively. The weight of the kernel is 0 at the center pixel and 1 at the other pixels. Using this kernel, for each pixel, take 16 pixels in the neighborhood of each pixel and convert them into a set of bands to obtain the filtered image. Divide the original image and the filtered image, raise the power, and sum them to obtain the feature gemary.

[0040] The gradients in the x and y directions of the SAR remote sensing image are calculated using the built-in algorithm of GEE. The feature gradient is calculated as follows:

[0041]

[0042] The characteristic direction is calculated according to the following formula:

[0043] direction = ytan -1 x

[0044] The image is filtered using a 3×3 Laplacian 8 edge detection normalization kernel to obtain the sharpening feature edgy.

[0045] Step C3: For the optical band portion of the superimposed image X, calculate the kernel exponent features. The kernel exponent features include sub-features kNDVI, kNDBBI, and kMNDWI. The calculation of kNDVI is as follows:

[0046] Based on NDVI, using the radial basis function (RBF) regeneration kernel, we have the following equation:

[0047]

[0048] In the formula, nir and r are the reflectivities of the near-infrared channel and the red channel, respectively, and k is the radial basis function (RBF) kernel function:

[0049] k(a,b)=exp(-(ab) 2 / (2σ 2 ))

[0050] Substituting into the RBF kernel formula, the exponent can be simplified to:

[0051]

[0052] In the formula, σ is the distance between the near-infrared spectrum and the red spectrum;

[0053] Similarly, the kNDBBI formula can be obtained as follows:

[0054]

[0055] In the formula, SWIR2 is the shortwave infrared band; g is the green light band.

[0056] The formula for kMNDWI is as follows:

[0057]

[0058] Where mir represents the mid-infrared band;

[0059] Step C4: Using principal component analysis, obtain the principal components of the optical band portion of the overlay image X. Perform object-oriented segmentation on the first principal component to obtain segmented images in units of pixel groups. Take the average pixel value of each segmented image pixel group to obtain the principal component segmentation features. At this point, the extraction of all types of features is completed, and the feature set W1 is obtained. Overlay the overlay image X with the feature set W1 to obtain the feature fusion image Y.

[0060] As a preferred embodiment of the present invention: the Mahalanobis distance value in step D, i.e., the feature separability of two different category feature sets, is calculated as follows:

[0061]

[0062] In the formula, μ is the mean distance of a certain category of sample set, and μ = (μ1, μ2, μ3, ..., μ n ) T x is a multivariate, that is, the mean distance of the sample sets of other classes, x = (x 1 ,x 2 ,…,x n ) T S is the covariance matrix of x.

[0063] As a preferred embodiment of the present invention, the calculation steps for the overall classification accuracy in step E are as follows:

[0064] Step E1: Calculate the confusion matrix as follows:

[0065]

[0066]

[0067] Element X in the confusion matrix ij This represents the number of pixels that are actually of category j but are classified as category i. If i = j, then this represents the number of pixels whose classified category matches the actual category. n is the total number of categories.

[0068] Step E2: Calculate the user classification accuracy UA, which is the percentage of pixels whose classified category matches the actual category among all pixels. The expression is as follows:

[0069]

[0070] In the formula, X ni The number of pixels that fall into the nth category of the ground feature detection point are correctly classified, which is represented by the row in the confusion matrix: class n; This represents the number of pixels classified into class n, which is the sum of all elements in the row of class n in the confusion matrix;

[0071] Step E3: Calculate the producer classification accuracy PA, i.e., calculate the omission rate of each category, as shown in the following formula:

[0072]

[0073] In the formula, X in This represents the number of pixels correctly classified into class n, which is the column in the confusion matrix: class n. This represents the total number of pixels that actually belong to class n, which is the sum of all elements in class n of the confusion matrix;

[0074] Step E4: Calculate the overall classification accuracy OA, which is the ratio of the total number of correctly classified pixels to the total number of all pixels. The expression is as follows:

[0075]

[0076] Step E5: Calculate the Kappa coefficient, which represents the rate of reduction in error for random classification. Its expression is as follows:

[0077]

[0078] In the formula, n represents the number of categories, N represents the number of test samples for one of the categories, and X ii X represents the number of pixels that were correctly classified. i+ X represents the number of pixels that actually belong to other classes but are classified as class i. +i This represents the number of pixels that actually belong to class i but are classified as other classes.

[0079] As a preferred technical solution of the present invention: In step F, in the gradient boosting tree classifier, the training dataset D is defined as... The boosting tree model can then be represented as:

[0080]

[0081] Among them, g m (x) represents the m-th decision tree, and M represents the number of decision trees.

[0082] Define the initial lifting tree f0(x) = 0. According to the forward step-by-step method, the model for the m-th step is:

[0083] f m (x)=f m-1 (x)+g m (x)

[0084] In m iterations, a base learner g is found. m (x) minimizes the loss and ultimately yields the regression boosting tree function model.

[0085] The gradient boosting tree algorithm uses the negative gradient of the loss function at the current model level as an approximation of the residual of the boosting tree model, and the gradient value is continuous. The expression for the negative gradient of its loss function is as follows:

[0086]

[0087] Beneficial effects: Compared with the prior art, the advantages of the present invention include:

[0088] This invention designs a wetland information extraction method based on multi-source remote sensing data. This method integrates multiple data types to comprehensively and accurately represent ground information; it overlays multiple temporal data to preserve temporal characteristics; it introduces principal component segmentation features to fully obtain the geometric shape information of the ground; and it employs a feature optimization method based on Mahalanobis distance to reasonably retain effective features. This method can comprehensively acquire ground features while avoiding feature redundancy, balancing computational cost and classification accuracy, thus achieving high-precision wetland classification. Attached Figure Description

[0089] Figure 1 This is a flowchart illustrating the wetland information extraction method based on multi-source remote sensing data provided in an embodiment of the present invention.

[0090] Figure 2 This is the result of drawing samples of the study area provided in the embodiments of the present invention;

[0091] Figure 3 This is the visualization result of the multi-temporal index features extracted in the embodiments of the present invention;

[0092] Figure 4 This is the visualization result of the texture features extracted in the embodiments of the present invention;

[0093] Figure 5 This is the visualization result of the multi-temporal kernel index features extracted in the embodiments of the present invention;

[0094] Figure 6 This is the visualization result of the principal component segmentation features extracted in the embodiments of the present invention;

[0095] Figure 7 This describes the variation of Mahalanobis distance for different feature combinations in the feature sensitivity experiment of this invention embodiment;

[0096] Figure 8 This describes the changes in OA for different feature combinations in the feature sensitivity experiment of this invention embodiment;

[0097] Figures 9(a)-9(d) This describes the variation of Mahalanobis distance values ​​under different numbers of features in this embodiment of the invention.

[0098] Figure 10 This describes the changes in OA under different feature quantities in the embodiments of the present invention;

[0099] Figure 11 These are curves showing the variation of OA with classifier parameters under different numbers of features in embodiments of the present invention;

[0100] Figure 12 This is a classification diagram of single-phase and multi-phase methods in the Yangtze River Estuary study area according to embodiments of the present invention;

[0101] Figure 13 This is a classification diagram of different ensemble learning methods in the Yangtze River Estuary research area according to embodiments of the present invention. Detailed Implementation

[0102] The present invention will be further described below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and should not be used to limit the scope of protection of the present invention.

[0103] Previous wetland information extraction studies have faced two main challenges. Firstly, the traditional method of processing images one by one using specialized offline software is often slow, costly, and unsuitable for large-scale data synthesis and computation. Secondly, given the complex and ever-changing nature of wetlands, data from a single time period is insufficient to represent the complete wetland situation. Cloud processing using the GEE platform can significantly improve computational efficiency. Furthermore, using overlays of three months' worth of optical, microwave, and DEM data for collaborative classification not only complements information from different images, fully expressing ground characteristics from multiple perspectives, but also comprehensively represents temporal information, reducing classification errors. Additionally, traditional feature processing methods often extract a small number of features and directly stack them into the classifier. However, using weighted judgments based on feature importance indices and Mahalanobis distance to optimize feature selection avoids feature redundancy, reduces computational costs, and improves classification efficiency and accuracy.

[0104] In terms of feature extraction, the exponential bands are placed in a Gaussian kernel space to form high-dimensional kernel exponential features, which are more conducive to classification. Utilizing object-oriented principles, the multi-source data fusion image is segmented into pixel-based segments, effectively preserving the geometry of the ground. Based on these ideas, a multi-source, multi-temporal feature extraction and optimization method for wetland information extraction is proposed.

[0105] The present invention provides a wetland information extraction method based on multi-source remote sensing data, referring to... Figure 1 Using the Yangtze River Estuary study area as the target region, multi-temporal remote sensing images of land features including various preset categories were collected. For these multi-temporal remote sensing images, steps A-F were performed to classify each pixel, thereby achieving land feature classification in the target region.

[0106] Step A: Refer to Figure 2 Based on the phenological characteristics and distribution features of each land use type in the wetland classification system, the land features in the target area are classified, and multi-temporal remote sensing images of the target area containing each preset category of land features are collected. The multi-temporal remote sensing images include optical multispectral images, SAR remote sensing images, and DEM data. Then proceed to step B.

[0107] Step B: Preprocess the multispectral image and SAR remote sensing image respectively, and fuse the preprocessed multispectral image, SAR remote sensing image and DEM data to obtain a three-temporal median composite image, i.e., overlaid image X, and then proceed to step C.

[0108] Step C: For the optical band portion of the overlaid image X, extract multi-temporal index features, including sub-features: Normalized Difference Vegetation Index (NDVI), Modified Normalized Difference Water Index (MNDWI), and Normalized Difference Bare Land and Built-up Area Index (NDBBI). The visualization results of each multi-temporal index feature are shown in [reference needed]. Figure 3 ;

[0109] Texture features are extracted from the overlay image X, including texture features of the optical band portion of the overlay image X, which includes sub-features ASM, contrast, entropy, and correlation; and texture features of the microwave band portion of the overlay image X, which includes sub-features entropy, intertia, gearys, gradient, direction, and edgy. The visualization results of each texture feature are shown in [reference needed]. Figure 4 ;

[0110] For the optical band portion of the overlaid image X, kernel index features are calculated. These kernel index features include sub-features kNDVI, kNDBBI, and kMNDWI. The visualization results of each kernel index feature are shown in [reference needed]. Figure 5For the optical band portion of the overlaid image X, principal component analysis is used to obtain the principal components of the optical band portion of the overlaid image X. Object-oriented segmentation is then performed on the first principal component to obtain principal component segmentation features. This completes the extraction of all types of features, resulting in a feature set W1. The overlaid image X is then superimposed with the feature set W1 to obtain the feature fused image Y, and then proceed to step D. The visualization results of the principal component segmentation features are shown in [reference needed]. Figure 6 ;

[0111] Step D: For the feature fusion image Y, calculate the Mahalanobis distance value of the sample points of each preset category of land cover in the feature fusion image Y to represent the separability between different preset category combinations of land cover; then proceed to step E; in the following steps, Y is replaced with other images, and the Mahalanobis distance value in other cases is calculated using the same principle;

[0112] Step E: Combine the optical band portion of the superimposed image X with multi-temporal index features, kernel index features, texture features of the optical band portion, and texture features of the microwave band portion to obtain feature maps after each combination. Use a gradient boosting tree classifier to classify each pixel in each feature map after each combination, calculate the overall classification accuracy and Mahalanobis distance value, and define the combined features whose classification accuracy and Mahalanobis distance value are both higher than the preset threshold as sensitive features.

[0113] The feature fusion image Y obtained in step C is fed into a gradient boosting tree classifier for training. Simultaneously, the explain algorithm in the GEE platform is used to calculate the feature importance index of all sub-features in the feature fusion image Y. Then, the optical band portion of the overlay image X is filtered out, and all remaining sub-features are sorted in descending order of feature importance index. Based on the optical portion of the overlay image X, new features are added one by one in descending order. After each feature addition, the overall classification accuracy and Mahalanobis distance are calculated, resulting in curves showing the changes in overall classification accuracy and Mahalanobis distance with the number of features. By analyzing the trends of the curves, and assuming that most sensitive features are included and the feature dimension is small, the first N features are determined as necessary features. Finally, it is checked whether any sensitive features were not selected from the remaining features; if so, they are added to obtain the final feature set W2, and then the process proceeds to step F.

[0114] Step F: Input the final feature set W2 into the gradient boosting tree classifier to classify each pixel in the final feature set W2, obtain the preset category to which each pixel belongs, and thus realize the land cover classification of the Yangtze River Estuary study area.

[0115] In this embodiment, taking the Yangtze River Estuary as the study area, the effectiveness of the wetland information extraction method based on multi-source remote sensing data is verified as follows:

[0116] 1 Experimental Setup

[0117] (1) Feature sensitivity analysis

[0118] Different features are superimposed and combined with the original optical spectral features, and then fed into the classifier for training. The Mahalanobis distance value and the overall classification accuracy are calculated to evaluate the inter-class separability and classification accuracy. The sensitivity of features is evaluated through the above two indicators.

[0119] Among them, Mahalanobis distance is obtained by substituting the mean values ​​of various samples; the feature combinations are divided into 9 categories, namely: spectral features plus spectral texture features, spectral features plus DEM terrain features, spectral features plus microwave features, spectral features plus exponential features, spectral features plus kernel exponential features, spectral features plus principal component segmentation features, spectral features plus microwave texture features, multi-temporal spectral features, and single-temporal spectral features. Among them, the single-temporal spectral features are features extracted from the synthetic effects in August 2020.

[0120] (2) Feature selection based on Mahalanobis distance

[0121] Mahalanobis distance was used for feature selection experiments. First, all features were ranked by importance, and the initial features were set as the original spectral features from multiple time periods. Then, features were added one by one according to their importance index from high to low for training. The Mahalanobis distance change curve and the overall classification accuracy change curve were calculated, and features with reasonable dimensions were selected based on the results of the feature sensitivity experiment.

[0122] While keeping other conditions unchanged, the accuracy of feature sets of different dimensions is calculated. 50% is randomly selected as training samples and 50% as validation samples, with no overlap between the two. The overall classification accuracy is calculated. The above steps are repeated 30-50 times to calculate the average value of the overall classification accuracy. A line graph of the overall classification accuracy as a function of the number of features is then plotted, and the error is marked.

[0123] (3) Classification parameter optimization

[0124] Parameter optimization is based on the GBDT classifier, with the parameter being the number of trees. The number of trees is set to be between 10 and 120, and a value is taken every 10 units to obtain a graph showing the overall classification accuracy as a function of parameter values ​​over 30 iterations.

[0125] (4) Comparison between single-phase and multi-phase

[0126] The results obtained by combining four different features, one single-phase and one multi-phase, are compared. Single-phase spectrum means that the input feature is only the spectral band of a single phase of the optical image, while multi-phase spectrum means that the spectral bands of three phases are superimposed. Single-phase feature means that the superimposed spectral bands of the existing feature set are removed, while multi-phase feature means that the multi-phase feature set is selected after feature selection.

[0127] (5) Classifier comparison

[0128] The comparison methods used to verify the effectiveness of the model include:

[0129] Two basic learning tools: Decision Tree and Classification and Regression Tree (CART);

[0130] Two ensemble learning methods: Gradient Boosting, where the base learners belong to the Boosting ensemble learning framework of CART regression trees; and Random Forest, where the base learners are multiple decision trees constructed from multiple different training samples and different features.

[0131] (6) Evaluation indicators

[0132] The classification results were quantitatively evaluated by statistically analyzing and comparing overall accuracy (OA), user accuracy (UA), producer accuracy (PA), and the Kappa coefficient (Kappa). For all classification algorithms used, the evaluation metric was the average of the results from 30-50 independent runs with randomly initialized training samples.

[0133] 2 Experimental Results

[0134] (1) Feature sensitivity analysis

[0135] Figure 7 Figure 8 The graphs show the variation curves of Mahalanobis distance and the overall accuracy results. Generally, the larger the Mahalanobis distance value, the higher the overall classification accuracy and the better the classification effect. However, although the Mahalanobis distance variation curves show the obvious superiority of the combination of multi-temporal spectral features and multi-temporal spectral texture features, this feature combination did not yield good results in the overall accuracy calculation results after averaging ten calculations.

[0136] By observing the variation values ​​of Mahalanobis distance, it was found that the combination of multi-temporal spectral features and multi-temporal index features also showed good results, especially in distinguishing between impermeable surfaces and Spartina alterniflora, and between impermeable surfaces and reeds. It was even superior to multi-temporal spectral texture features. The overall accuracy obtained indicates that multi-temporal index features have high sensitivity and are suitable for the classification of wetland ecosystems in this region. Furthermore, the fusion of multi-temporal spectral features and microwave features brought complementary information and also showed good inter-class separability and overall classification accuracy results, especially between several difficult-to-distinguish categories such as impermeable surfaces and bare tidal flats, impermeable surfaces and reeds, cultivated land and bare tidal flats, Spartina alterniflora and bare tidal flats, and bare tidal flats and ponds.

[0137] In summary, multi-phase exponential features and multi-phase microwave features have the highest feature sensitivity and should be given priority in subsequent feature selection. Using multi-phase feature superposition has a significant effect improvement over using single-phase features. In subsequent experiments, a comprehensive comparative verification experiment will be conducted to illustrate the differences between single-phase and multi-phase features.

[0138] (2) Evaluation of high-quality feature set

[0139] As shown in Figure 9, Mahalanobis distance generally exhibits a positive correlation with the number of features. When the number of features is less than 44, the Mahalanobis distance increases rapidly with the number of features, with a large slope. However, when the number of features is greater than 44, the Mahalanobis distance value shows a slow upward trend. Even if there are sudden increases in the separability between individual classes, it is not enough to affect the overall separability. In addition, when the number of features reaches 120, the Mahalanobis distance values ​​between some classes show a discontinuous increase, such as between impermeable surfaces and Spartina alterniflora, impermeable surfaces and ponds, impermeable surfaces and water bodies, cultivated land and water bodies, cultivated land and Spartina alterniflora, ponds and reeds, as well as the distinguishability between Rubus tricuspidata and other classes, all of which show a significant increase.

[0140] Figure 10 The graph shows the change in OA with the number of features. It can be seen that when the feature dimension is less than 58, the number of features has a significant impact on OA, showing a positive and rapid growth trend. However, when the feature dimension is greater than 58, the OA region is stable, fluctuating around 88%. This also verifies that Mahalanobis distance and OA change are roughly positively correlated. When the Mahalanobis distance value is large, the inter-class separability is better, the OA value is relatively large, and the classification effect is better.

[0141] Based on the above analysis, in order to balance computational cost and inter-class separability, weigh the huge load that high-dimensional features bring to the classifier against the relatively small inter-class separability of low-dimensional features, and at the same time ensure that too much redundancy is avoided, this study first selected the top 88 features as necessary bands for final classification according to feature importance, and created a table of Mahalanobis distance variation to verify the rationality of the method.

[0142] In addition, considering the conclusions of the feature sensitivity experiment, it is necessary to ensure that the multi-temporal exponential features and multi-temporal microwave features are added to the final feature set. Among the first 88 features, two multi-temporal exponential features (MNDWI) are missing, so they are included to obtain the final set of 90 feature dimensions, as shown in Table 1.

[0143] Table 1

[0144]

[0145]

[0146] (3) Influence of classification parameters

[0147] According to the curve Figure 11 It can be seen that when the parameter is greater than or equal to 70, the OA value tends to be stable, with an accuracy of around 88.02%. To balance accuracy and computational cost, a parameter of 70 is used for classification.

[0148] (4) The effect of multiple time phases superposition

[0149] The experimental data results comparing single-phase and multi-phase data are shown in Table 2, and the visualization images are as follows: Figure 12 As shown, compared to multi-temporal features, inputting only the original single-temporal spectral features or inputting all other single-temporal features results in more noise, especially in the northeastern marine region where cloud cover severely obscures the image. Furthermore, single-temporal features cannot effectively represent changing land types, leading to lower accuracy. Single-temporal spectral features yielded an OA of 72.73% and a Kappa of 0.67; while multi-temporal spectral features achieved an OA of 83.81% and a Kappa of 0.8, with an OA increase of 11.8% compared to single-temporal features, reflecting the importance of multi-temporal spectral bands. Moreover, single-temporal features, based on existing features, achieved better classification results than using only spectral features, with an overall accuracy of 85.62%, an increase of 12.89% compared to single-temporal spectral features and 1.81% compared to multi-temporal spectral features, demonstrating the effectiveness of the final feature set. Finally, the classification accuracy obtained by using the final multi-temporal feature set reached 88.02%, which is 2.4% higher than the OA and 0.02 higher Kappa than the feature set obtained by removing the other two temporal phases. The experiment shows that the multi-temporal feature set used in this study has richer information than other single-temporal feature sets and single multi-temporal spectral feature sets, and is suitable for the study of wetland information extraction in the Yangtze River Estuary.

[0150] Table 2

[0151]

[0152] (5) Analysis of the performance of different classifiers

[0153] Figure 13Table 3 visualizes the accuracy evaluation results of four different classifiers in the Yangtze River Estuary region, showing the specific classification accuracy of different methods under the current multi-temporal feature set. It was found that the classification accuracy obtained using the Decision Tree (DT) and CART classifiers was relatively low. DT's OA was 77.63% with a Kappa coefficient of 0.72; CART's was 78.34% with a Kappa coefficient of 0.74. Overall, CART's accuracy was 0.71% higher than DT's, with little difference in overall performance. Observing the classification results, it was found that the CART classifier resulted in large areas of misclassified patches, especially in the eastern water body, which was misclassified as bare tidal flats and other wetland vegetation. The DT classifier exhibited significant salt-and-pepper noise, resulting in poor classification quality. Compared to the above two classifiers, Random Forest (RF) and Gradient Boosting Tree (GBDT) achieved better classification results. The RF classifier achieved an OA of 88.09% with a Kappa of 0.85; the GBDT classifier also achieved over 88% accuracy with a Kappa coefficient of 0.85. The classification results also show that both classifiers can effectively preserve the true ground shape features, ensuring pixel purity in the classification results. Although the overall classification accuracy of RF is slightly higher than that of GBDT, and it also shows superior performance in the 2020 classification results, by observing user classification accuracy (UA) and producer classification accuracy (PA), GBDT has a more significant advantage in distinguishing important wetland vegetation (Scirpus triqueter, Spartina alterniflora, and Phragmites australis) in this study area. For example, the PA of Scirpus triqueter under RF is 52.2%, and the UA is 67.18; while GBDT's are 73.09% and 71.27%, respectively, with PA and UA being 20.89% and 4.09% higher than RF, respectively. Similarly, the PA and UA of Spartina alterniflora under GBDT classification are 5.5% and 6.41% higher than RF, respectively; and 5.76% and 1.97% higher for Phragmites australis. Based on the above analysis, the effectiveness of using GBDT as a classifier for wetland information extraction in the Yangtze River Estuary is fully demonstrated.

[0154] Table 3

[0155]

[0156] As can be seen from the above embodiments, this invention is a wetland information extraction method based on multi-source remote sensing data. It couples optical, microwave, and DEM data. First, after preprocessing such as screening, fusion, correction, and cropping, microwave and optical data are superimposed with three temporal data to extract multi-temporal features. Elevation and slope features are extracted using the DEM. Samples are drawn based on the Google Earth Pro platform, with 50% selected as training and validation samples respectively. Then, feature sensitivity analysis, feature importance ranking, and optimization methods are used to determine the feature set. Single-temporal and multi-temporal feature comparison, classifier comparison, and parameter optimization methods are used to determine the classifier. Finally, spatiotemporal variation analysis of the Yangtze River Estuary wetland ecosystem is performed. This invention can comprehensively acquire ground features while avoiding feature redundancy, considering both computational cost and classification accuracy, achieving high-precision wetland classification.

[0157] The embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention.

Claims

1. A method for extracting wetland information based on multi-source remote sensing data, characterized in that, For the target area, acquire multi-temporal remote sensing images of the target area containing various preset categories of land features. For the multi-temporal remote sensing images, perform the following steps A-F to complete the classification of each pixel in the multi-temporal remote sensing images, thereby achieving land feature classification of the target area: Step A: Based on the phenological characteristics and distribution features of each land use type in the wetland classification system, classify the land features in the target area, collect multi-temporal remote sensing images of the target area containing each preset category of land features, including optical multispectral images, SAR remote sensing images, and DEM data, and then proceed to Step B. Step B: Preprocess the multispectral image and SAR remote sensing image respectively, and fuse the preprocessed multispectral image, SAR remote sensing image and DEM data to obtain a three-temporal median composite image, i.e., overlaid image X, and then proceed to step C. Step C: For the optical band portion of the superimposed image X, extract multi-temporal index features, which include sub-features NDVI, MNDWI, and NDBBI. Extract the texture features of the superimposed image X, including the texture features of the optical band portion of the superimposed image X, which includes sub-features ASM, contrast, entropy, and correlation; and the texture features of the microwave band portion of the superimposed image X, which includes sub-features entropy, intertia, gearys, gradient, direction, and edgy. For the optical band portion of the overlay image X, kernel index features are calculated, including sub-features kNDVI, kNDBBI, and kMNDWI. Principal component analysis is used to obtain the principal components of the optical band portion of the overlay image X. Object-oriented segmentation is performed on the first principal component to obtain principal component segmentation features. At this point, the extraction of all types of features is completed, resulting in a feature set W1. The overlay image X is superimposed with the feature set W1 to obtain the feature fusion image Y, and then the process proceeds to step D. Step D: For the feature fusion image Y, calculate the Mahalanobis distance value of the sample points of each preset category of land cover in the feature fusion image Y to represent the separability between different preset category combinations of land cover; then proceed to step E; In the following steps, Y is replaced with other images, and the Mahalanobis distance values ​​for other cases are calculated using the same principle; Step E: Combine the optical band portion of the superimposed image X with multi-temporal index features, kernel index features, texture features of the optical band portion, and texture features of the microwave band portion to obtain feature maps after each combination. Use a gradient boosting tree classifier to classify each pixel in each feature map after each combination, calculate the overall classification accuracy and Mahalanobis distance value, and define the combined features whose classification accuracy and Mahalanobis distance value are both higher than the preset threshold as sensitive features. The feature fusion image Y obtained in step C is fed into a gradient boosting tree classifier for training. Simultaneously, the explain algorithm in the GEE platform is used to calculate the feature importance index of all sub-features in the feature fusion image Y. Then, the optical band portion of the overlay image X is filtered out, and all remaining sub-features are sorted in descending order of feature importance index. Based on the optical portion of the overlay image X, new features are added one by one in descending order. After each feature addition, the overall classification accuracy and Mahalanobis distance are calculated, resulting in curves showing the changes in overall classification accuracy and Mahalanobis distance with the number of features. By analyzing the trends of the curves, the first N features are determined as necessary features. Finally, it is checked whether any sensitive features have not been selected among the remaining features; if so, they are added to obtain the final feature set W2, and then the process proceeds to step F. Step F: Input the final feature set W2 into the gradient boosting tree classifier to classify each pixel in the final feature set W2, obtain the preset category to which each pixel belongs, and thus realize the land cover classification of the target area.

2. The wetland information extraction method based on multi-source remote sensing data according to claim 1, characterized in that, Step B involves preprocessing the multispectral imagery, including cloud removal (removing areas with cloud cover greater than 5%), atmospheric correction, date filtering, and cropping according to the vector boundaries of the study area. Preprocessing the SAR remote sensing imagery includes atmospheric correction, cropping according to the vector boundaries of the study area, and fusing the upper and lower orbital data. The resulting SAR remote sensing imagery is then filtered by date and mode.

3. The wetland information extraction method based on multi-source remote sensing data according to claim 1, characterized in that, The specific steps of step C are as follows: Step C1: For the optical band portion of the superimposed image X, calculate the sub-features NDVI, MNDWI, and NDBBI of the multi-temporal exponential features as follows: In the formula, NIR is the near-infrared band in the optical band portion of the superimposed image X, Red is the red light band, Green is the green light band, MIR is the mid-infrared band, and SWIR2 is the short-wave infrared band. Step C2: Extract the sub-features ASM, contrast, entropy, and correlation of the texture features in the optical band portion of the overlay image X, and the sub-features entropy, intertia, gearys, gradient, direction, and edgy of the texture features in the microwave band portion of the overlay image X, as follows: Let the pixel value of the pixel at positions (x,y) and (x+d,y+l) in the overlay image X be (i,j), the number of pixels with pixel value (i,j) in the overlay image X be N(i,j), and the total number of pixels in the overlay image X be N. glcm Then the probability P(i,j) of pixel value (i,j) appearing in the superimposed image X is: The ASM expression for the energy feature of pixel value (i,j) is: The contrast feature expression for pixel value (i,j) is: The entropy expression for pixel value (i,j) is: The autocorrelation expression for pixel value (i,j) is: In the formula, Set the matrix size to 6, and use the GEE built-in algorithm glcmTexture to calculate the number of times a pixel with pixel value X is adjacent to a pixel with pixel value Y in a preset direction and distance. Count the number of adjacent pixels to obtain the feature intertia. Define a 9×9 kernel with a focal position of (-4, -4). Its horizontal and vertical coordinates represent the offsets from the left and top of the kernel, respectively. The weight of the kernel is 0 at the center pixel and 1 at the other pixels. Using this kernel, for each pixel, take 16 pixels in the neighborhood of each pixel and convert them into a set of bands to obtain the filtered image. Divide the original image and the filtered image, raise the power, and sum them to obtain the feature gemary. The gradients in the x and y directions of the SAR remote sensing image are calculated using the built-in algorithm of GEE. The feature gradient is calculated as follows: The characteristic direction is calculated according to the following formula: direction=surface -1 x The image is filtered using a 3×3 Laplacian 8 edge detection normalization kernel to obtain the sharpening feature edgy. Step C3: For the optical band portion of the superimposed image X, calculate the kernel exponent features. The kernel exponent features include sub-features kNDVI, kNDBBI, and kMNDWI, where kNDVI is calculated as follows: Based on NDVI, using the radial basis function (RBF) regeneration kernel, we have the following equation: In the formula, nir and r are the reflectances of the near-infrared channel and the red channel, respectively, k is the radial basis function (RBF) kernel function, and σ is the distance between the near-infrared spectrum and the red spectrum; Similarly, the kNDBBI formula can be obtained as follows: In the formula, SWIR2 is the shortwave infrared band; g is the green light band; The formula for kMNDWI is as follows: Where mir represents the mid-infrared band; Step C4: Using principal component analysis, obtain the principal components of the optical band portion of the overlay image X. Perform object-oriented segmentation on the first principal component to obtain segmented images in units of pixel groups. Take the average pixel value of each segmented image pixel group to obtain the principal component segmentation features. At this point, the extraction of all types of features is completed, and the feature set W1 is obtained. Overlay the overlay image X with the feature set W1 to obtain the feature fusion image Y.

4. The wetland information extraction method based on multi-source remote sensing data according to claim 3, characterized in that, In step D, the Mahalanobis distance, i.e., the feature separability of two different feature sets, is calculated as follows: In the formula, μ is the mean distance of a certain category of sample set, and μ = (μ1, μ2, μ3, ..., μ n ) T x is a multivariate, that is, the mean distance of the sample sets of other classes, x = (x 1 ,x 2 ,…,x n ) T S is the covariance matrix of x.

5. The wetland information extraction method based on multi-source remote sensing data according to claim 4, characterized in that, The steps for calculating the overall classification accuracy in step E are as follows: Step E1: Calculate the confusion matrix as follows: Element X in the confusion matrix ij This represents the number of pixels that are actually of category j but are classified as category i. If i = j, then this represents the number of pixels whose classified category matches the actual category. n is the total number of categories. Step E2: Calculate the user classification accuracy UA, which is the percentage of pixels whose classified category matches the actual category among all pixels. The expression is as follows: In the formula, X ni The number of pixels that fall into the nth category of the ground feature detection point are correctly classified, which is represented by the row in the confusion matrix: class n; This represents the number of pixels classified into class n, which is the sum of all elements in the row of class n in the confusion matrix; Step E3: Calculate the producer classification accuracy PA, i.e., calculate the omission rate of each category, as shown in the following formula: In the formula, X in This represents the number of pixels correctly classified into class n, i.e., the column in the confusion matrix: class n; This represents the total number of pixels that actually belong to class n, which is the sum of all elements in class n of the confusion matrix; Step E4: Calculate the overall classification accuracy OA, which is the ratio of the total number of correctly classified pixels to the total number of all pixels. The expression is as follows: Step E5: Calculate the Kappa coefficient, which represents the rate of reduction in error for random classification. Its expression is as follows: In the formula, n represents the number of categories, N represents the number of test samples for one of the categories, and X ii X represents the number of pixels that were correctly classified. i+ X represents the number of pixels that actually belong to other classes but are classified as class i. +i This represents the number of pixels that actually belong to class i but are classified as other classes.

6. The wetland information extraction method based on multi-source remote sensing data according to claim 5, characterized in that, In step F, within the gradient boosting tree classifier, the training dataset is defined. The boosting tree model can then be represented as: Among them, g m (x) represents the m-th decision tree, and M represents the number of decision trees; Define the initial lifting tree f0(x) = 0. According to the forward step-by-step method, the model for the m-th step is: f m (x)=f m-1 (x)+g m (x) In m iterations, a base learner g is found. m (x) minimizes the loss, ultimately yielding the regression boosting tree function model; The gradient boosting tree algorithm uses the negative gradient of the loss function at the current model level as an approximation of the residual of the boosting tree model, and the gradient value is continuous. The expression for the negative gradient of its loss function is as follows:

Images