The Change Pattern and Its Dominant Driving Factors of Wetlands in the Yellow River Delta Based on Sentinel-2 Images

: There were signiﬁcant differences in the dominant driving factors of the change process of different types of wetlands in the Yellow River delta. In addition, to our knowledge, the optimal classiﬁcation feature sets with the Random Forest algorithm for wetlands in the Yellow River delta were least explored. In this paper, the wetland information in the study area was extracted based on a Random Forest algorithm with de-feature variable redundancy, and then the change process of wetland and its dominant factors from 2015 to 2021 was monitored and analyzed using the Geodetector and gravity center model. The results showed that (1) the optimal variable sets composed of red edge indexes based on the Random Forest algorithm had the highest classiﬁcation accuracy, with the overall accuracy and Kappa coefﬁcient of 95.75% and 0.93. (2) During 2015–2021, a large area of natural wetland in the Yellow River delta was transformed into an artiﬁcial wetland. The wetlands showed an overall development direction of “northwest–southeast” along the Yellow River. (3) The interaction between vegetation coverage and accumulated temperature had the largest explanatory power of the change in the natural wetland area. The interaction between solar radiation and DEM had the largest explanatory power for the change in the artiﬁcial wetland area. The research results could better provide decisions for wetland protection and restoration in the Yellow River delta.


Introduction
Wetlands, an important ecosystem, can provide critical functions and values in water supply, biodiversity maintenance, climate regulation, and water purification. Under the double stress of global climate change and intensified human activity, the wetlands showed a degradation trend, with the shrink area and declined function. Effective monitoring of wetland change processes can help reveal and understand the causes of degradation, which is conducive to better restoring ecological functions and promoting wetland ecological health. Therefore, monitoring the change patterns of wetlands and determining their dominant driving factors are essential for regional sustainable development [1]. The Yellow River delta (YRD) refers to the alluvial plain formed by the sediments carried by the Yellow River in the depression of the Bohai Sea, which has the most complete temperate wetland ecosystem in China. It was officially listed as an important international wetland in 2013. At present, the YRD faces three serious ecological challenges: the severe formalization of different periods are unclear. The objectives of this paper are as follows: (1) to determine the optimal feature variable sets for wetlands classification in the YRD; (2) to extract the information of wetlands in the YRD based on the RF algorithm and the optimal feature variable sets; (3) to monitor and reveal the change patterns of the wetland ecosystem and its driving mechanisms in the YRD based on the Geodetector and gravity center model, which provided the basis for the protection and rational utilization of wetland resources in the YRD [20].

Study Area
The Yellow River delta, located in the Yellow River estuary of Dongying City, Shandong Province, China (117°31′-119°18′E, 36°55′-38°16′N), has the youngest and most complete wetland ecosystem, with an area of approximately 7033.405 km 2 ( Figure 1). In this study region, a large area of mud flats and wetlands is formed under river-sea integration actions. The soil types in the study area are mainly composed of coastal saline soil and tidal soil. The terrain is flat, with an average altitude less than 15 m. It belongs to a temperate monsoon climate with distinct seasons and sufficient light. Since the 1980s, a decrease in sediments in the Yellow River Basin has led to serious seawater recharge and coastal erosion, which is conducive to the formation of salinization.

Data Source and Preprocessing
To achieve the accurate extraction of wetland information in the YRD, this study selected a Sentinel-2 remote sensing image with a higher spatial resolution (10 m) and more available bands (13 bands) (United States Geological Survey, https://earthexplorer.usgs.gov/, accessed on 13 September 2015 and 8 September 2021) in August and September as the data source. The hydrological characteristics of the YRD in August and September were relatively stable in the growing seasons, which could better distinguish

Data Source and Preprocessing
To achieve the accurate extraction of wetland information in the YRD, this study selected a Sentinel-2 remote sensing image with a higher spatial resolution (10 m) and more available bands (13 bands) (United States Geological Survey, https://earthexplorer.usgs. gov/, accessed on 13 September 2015 and 8 September 2021) in August and September as the data source. The hydrological characteristics of the YRD in August and September were relatively stable in the growing seasons, which could better distinguish different vegetation types. The red-edge index derived from the red-edge bands of Sentinel-2 images had a large contribution rate in the process of extracting and classifying wetland information in the YRD [21]. The Sen2cor atmospheric correction model of SNAP software was used to elimi-Remote Sens. 2022, 14, 4388 4 of 24 nate the noises of atmospheric and aerosols. In addition, ENVI 5.5 was applied to conduct the process of BRDF (Bidirectional Reflectance Distribution Function) correction [22]. The final datasets after atmospheric correction and BRDF correction were at a spatial resolution of 10 m. The driving factors, including vegetation coverage (FVC), precipitation, soil type, and Gross Domestic Product (GDP), were obtained from the Resource and Environmental Science and Data Centre (https://www.resdc.cn/, accessed on 15 September 2021).

Methods
In this paper, the overall technical flowchart is shown in Figure 2. different vegetation types. The red-edge index derived from the red-edge bands of Sentinel-2 images had a large contribution rate in the process of extracting and classifying wetland information in the YRD [21]. The Sen2cor atmospheric correction model of SNAP software was used to eliminate the noises of atmospheric and aerosols. In addition, ENVI 5.5 was applied to conduct the process of BRDF (Bidirectional Reflectance Distribution Function) correction [22]. The final datasets after atmospheric correction and BRDF correction were at a spatial resolution of 10 m. The driving factors, including vegetation coverage (FVC), precipitation, soil type, and Gross Domestic Product (GDP), were obtained from the Resource and Environmental Science and Data Centre (https://www.resdc.cn/, accessed on 15 September 2021).

Methods
In this paper, the overall technical flowchart is shown in Figure 2.

Classification System of Wetlands
Wetlands were composed of natural or artificial, permanent, or temporary swamps, wet sources, mire lands or water areas, water bodies with static or flowing fresh, brackish, or salt water, and waters with a depth of not more than 6 m at low tide [23]. Combined with the actual situation of the YRD, the wetland classification system (Table 1) in this paper was formulated based on the national wetland classification standard. To accurately extract wetland information in the YRD, this study described the geometric and image features of ground objects based on Sentinel-2 images and high-resolution Google Earth images.

Classification System of Wetlands
Wetlands were composed of natural or artificial, permanent, or temporary swamps, wet sources, mire lands or water areas, water bodies with static or flowing fresh, brackish, or salt water, and waters with a depth of not more than 6 m at low tide [23]. Combined with the actual situation of the YRD, the wetland classification system (Table 1) in this paper was formulated based on the national wetland classification standard. To accurately extract wetland information in the YRD, this study described the geometric and image features of ground objects based on Sentinel-2 images and high-resolution Google Earth images.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.  Shallow sea Green, striped, distributed in coastal areas.
Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.

Remote Sensing
Image Geometric Feature

Natural wetland
River Green or yellow, slender striped.

Swamp
Dark green, irregular shape, largely distributed in nature reserves.
Mud flat Yellow, striped, distributed in coastal areas.
Shallow sea Green, striped, distributed in coastal areas.
Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.
Dark green, distributed along roads and rivers.

Construction land
Remote Sens. 2022, 14, x FOR PEER REVIEW 5 of 25 Table 1. Classification system of wetlands in the YRD.

Remote Sensing
Image Geometric Feature

Natural wetland
River Green or yellow, slender striped.

Swamp
Dark green, irregular shape, largely distributed in nature reserves.
Mud flat Yellow, striped, distributed in coastal areas.
Shallow sea Green, striped, distributed in coastal areas.
Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.
Pink and gray hybrid, regular in shape, distributed inland.

Remote Sensing
Image Geometric Feature

Natural wetland
River Green or yellow, slender striped.

Swamp
Dark green, irregular shape, largely distributed in nature reserves.
Mud flat Yellow, striped, distributed in coastal areas.
Shallow sea Green, striped, distributed in coastal areas.
Artificial wetland

Salt pans
Green or gray-white, grid-like, coastal distribution.

Paddy field
Bright green, striped, distributed along the river.

Reservoirs pond
Reservoir pond: Green, regular in shape, scattered in the study area; Loose pond: Green, striped, distributed in coastal areas.

Non-wetland
Dryland Green, regular shape, with evident texture features.

Woodland
Dark green, distributed along roads and rivers.

Construction land
Pink and gray hybrid, regular in shape, distributed inland.
Unused land Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.
Soil yellow, scattered distribution.

Description of Feature Variables
Based on Sentinel-2 images, this study extracted 35 feature variables, including spectral feature, texture feature [24][25][26], vegetation index, water index, red-edge index, and others [27], which are shown in Table 2.
Modified Simple Ratio Index red-edge MSRre B 2 is the blue band, B 3 is the green band, B 4 is the red band, B 5 is the vegetation red edge1 band, B 6 is the vegetation red edge2 band, B 7 is the vegetation red edge3 band, B 8 is near-infrared band, B 8a is narrow near-infrared band, B 11 is short-wave infrared1 band, B 12 is short-wave infrared2 band.

Texture Information Extraction
(1) Principal component analysis Principal component analysis (PCA) is a dimension reduction mathematical statistical method [30]. With the help of the orthogonal transformation, the component-related original random variables were transformed into new component-independent random variables, named the principal component [31,32]. In this paper, PCA was applied to eliminate the mutual influences among the bands of Sentinel-2 images, reduce the workload of texture information extraction, and improve work efficiency with the tools of spatial analyst of ArcGIS 10.5.
(2) Gray level co-occurrence matrix A gray level co-occurrence matrix (GLCM) can be used to extract the texture information of Sentinel-2 images, which was proposed by Haralick in 1973 [33]. First, the GLCM of the image was calculated, and then second-order statistical characteristics describing the texture were obtained based on this matrix, which has strong adaptability and robustness [34,35]. GLCM has been widely used in the texture extraction of SAR, multispectral, high-resolution images. Texture is a degree-dependent phenomenon caused by the interaction of the hue primitive space of an image [36]. The equation is as follows: where, # is the number of elements in the dataset; f (x 1 , y 1 ) = i represents the pixel gray of image position (x 1 , y 1 ); d represents the distance between two pixels, θ represents the direction angle between two pixels; d ranges from 1 to 4; θ refers to the values of 0 • , 45 • , 90 • , and 135 • .

Random Forest Algorithm
The random forest algorithm plays an important role in feature selection and dimension reduction, and is an algorithm for generating multiple decision trees through ensemble learning based on the CART algorithm invented by Breiman et al. [37,38]. It allows the fusion of high-dimensional data from multiple sources and has a strong tolerance for missing values and outliers. It is suitable for dealing with high-dimensional complex datasets and can automatically determine the importance of feature variables. A classification accuracy model sequence a is obtained by K-round training, which forms a multi-classification model system {h 1 (X), h 2 (X), . . . , h k (X)}. The final classification result can be obtained by a simple majority voting decision. The final classification decision was as follows: where H(x) represents the combination classification model, h i (x) represents a single decision tree classification model, Y represents an output variable,

Confusion Matrix
The confusion matrix is also called the error matrix, which has been widely used to compare the degree of confusion between classification results and actual measured values. This study applied the overall accuracy (OA), Kappa coefficient (Kappa), producer's accuracy (PA), and user's accuracy (UA) to evaluate the accuracy of each classification [40]. where N is the total number of samples, k is the total number of categories, N ii is the number of samples correctly classified, N i+ and N +i are the number of actual samples and predicted samples of the i-th category, respectively.

Wetland Information Extraction and Classification
In the classification process, there is a certain information redundancy among the 35 feature variables, which can increase the calculation complexity and reduce the accuracy in the process of classification. Therefore, it is necessary to optimize the feature variables before classification. The RF algorithm can not only be used to realize the classification of wetlands but also play an important role in feature selection and dimensional reduction. In this paper, the out-of-bag (OOB) data error of RF was used to calculate the importance of the feature variables and to determine the optimal feature variables. The expression of the importance of feature variables V X j was as follows: where, n is the number of decision trees in random forests, e t is the OOB error of each decision tree calculated by OOB, X j is the random assignment of the j-th feature variable of the OOB data, e j t is the OOB error calculated by the new value X j . The change of variable X j leads to an increase in OOB error and a decrease in accuracy.
According to the formula above, the importance score of 35 feature variables in the YRD in 2021 was calculated and is then shown in Figure 3. The GLCM_Variance ranked first, with an importance score of 0.1317, followed by the MNDWI, with an importance score of 0.1084. The importance scores of NDWI and GLCM_Mean were the same at 0.1046. The texture features for Homogeneity and Dissimilarity had smaller importance scores of 0.005 and 0.0005, respectively. According to the comprehensive feature set, the contribution rate of feature variables in the classification process of wetland in the YRD was ranked as red-edge index > vegetation and water index > spectral feature > texture feature. The contribution rate of the red-edge index was about 37.5%, while that of vegetation and water index was 25%. The contribution rates of spectral features and texture features were smaller at 20% and 17.5%, respectively.
The correlation matrix of the 35 feature variables in 2021 is shown in Figure 4. GLCM_variance was extremely strongly correlated with B11, the correlation coefficient was 0.84, and the correlation coefficient of MNDWI with B11 was 0.83; therefore, B11 should be removed [41]. The correlation coefficient thresholds are shown in Table 3. Similar to this process, 15 feature variables were removed. Finally, 20 feature variables were selected to construct the optimal feature variable sets for the classification of wetlands based on Sentinel-2 images. Compared to 2021, there was a certain degree of change in the types and spatial distribution patterns of wetlands in 2015, which led to the fluctuation of the importance scores of feature variables. Depending on the selection principles, the optimal feature variable sets for wetlands classification of the YRD in 2015 and 2021 are obtained in Table 4. There were slight differences in the feature variables between 2015 and 2021. However, the contribution rates of the feature variable set had the same ranking, which showed that the red-edge index had the greatest contribution rate in the process of wetland information extraction and classification.  The correlation matrix of the 35 feature variables in 2021 is shown in Figure 4. GLCM_variance was extremely strongly correlated with B11, the correlation coefficient was 0.84, and the correlation coefficient of MNDWI with B11 was 0.83; therefore, B11 should be removed [41]. The correlation coefficient thresholds are shown in Table 3. Similar to this process, 15 feature variables were removed. Finally, 20 feature variables were selected to construct the optimal feature variable sets for the classification of wetlands based on Sentinel-2 images. Compared to 2021, there was a certain degree of change in the types and spatial distribution patterns of wetlands in 2015, which led to the fluctuation of the importance scores of feature variables. Depending on the selection principles, the optimal feature variable sets for wetlands classification of the YRD in 2015 and 2021 are obtained in Table 4. There were slight differences in the feature variables between 2015 and 2021. However, the contribution rates of the feature variable set had the same ranking, which showed that the red-edge index had the greatest contribution rate in the process of wetland information extraction and classification. Importance score

Feature variation
Important score Ranking of important score   As shown in Table 5, six classification schemes were designed in this paper. Scheme A only used the spectral features of Sentinel-2 images (except B1, B9, B10). Scheme B was   As shown in Table 5, six classification schemes were designed in this paper. Scheme A only used the spectral features of Sentinel-2 images (except B1, B9, B10). Scheme B was composed of spectral features, vegetation index, and water index. Scheme C was used to select the spectral features of the Sentinel-2 images and the red-edge index. Scheme D was composed of spectral and texture features, while scheme E was a combination of all feature variables. Scheme F was composed of the optimal feature variable sets determined by the RF algorithm. The main purposes of setting different classification schemes are: (1) To investigate the contribution rate of different feature variables to wetland information extraction and determine the importance ranking of different feature variables; (2) To explore the optimal method to improve the extraction and classification of wetland information by comparing the accuracy among different classification schemes.

Accuracy Evaluation Analysis
It can be seen from Table 6 that the OA and Kappa of scheme A based on the spectral features were 82.91% and 0.82, respectively. After adding the vegetation index, water index, and red-edge index to scheme A, the classification accuracy of schemes B and C was slightly improved, with OA of 83.87% and 84.21%, and Kappa of 0.83 and 0.82, respectively. However, the classification accuracy of scheme D with texture features and scheme F with all feature variables decreased, with the OA of 81.74% and 81.54% and Kappa values of 0.80 and 0.79, respectively. This indicated that the vegetation index, water index, and red-edge index could improve the accuracy of wetland information extraction and classification, but the texture features derived from Sentinel-2 images with a spatial resolution of 10 m were not conducive to improving the accuracy of wetland information extraction. The OA and Kappa of classification scheme F were enhanced to 95.75% and 0.93, respectively, indicating that classification scheme F had the best applicability for wetland classification in the YRD. Single type of PA and UA could show that vegetation index and the red-edge index were conductive to improve the accuracy of wetland information classification. In schemes B and C, the accuracy of swamp, paddy field, and dryland were enhanced after adding the above two feature sets, and the contribution rate of the red-edge index was higher than that of the vegetation index. Texture features improved the accuracy of mud flat, shallow sea, reservoirs, ponds, and construction land, which indicated that the texture features could contribute to the classification of objects with evident texture information. When all the feature variables were involved in the classification, classification accuracy and Kappa were reduced due to information redundancy. In classification scheme F, the feature variables with high correlation were removed, which avoided the interference of overlapping redundant information on wetland information extraction and classification. The classification results are shown in Figure 5.

Wetland Information Extraction and Classification Results
Utilizing the optimal feature variable sets and RF algorithm, the spatial distribution of wetlands in 2015 and 2021 was obtained. As shown in Figure 6a, in 2015, the area of the natural wetland was 986.99 km 2 , while that of the artificial wetland was 1334.33 km 2 . In the natural wetland, the mud flat had the largest area of 432.90 km 2 , accounting for 6.16 % of the total area, which was mainly distributed in the northern Hekou district, the coastal area of Kenli district, and eastern Dongying district. The shallow sea had the second largest area of 353.47 km 2 , accounting for 5.03% of the total area, which was mostly located in the northern Hekou district, coastal areas in the Kenli district, and eastern Dongying district. The area of swamp was 129.79 km 2 , accounting for 1.85 %, which was mostly distributed in the estuary of the Yellow River in the northeast of Kenli district. The river had the smallest area of 70.82 km 2 , accounting for only 1.01% of the total area. In the artificial wetland, the reservoirs and pond had the largest area of 926.30 km 2 , accounting for 13.19%, followed by salt pans, with an area of 274.66 km 2 , accounting for 3.91%, which was mainly distributed in the northern Hekou district, Kenli district, Dongying district, and eastern Guangrao county. The paddy field area was 133.37 km 2 , accounting for only 1.90% of the total area, which was largely distributed in the central and northeastern Kenli district. accuracy and Kappa were reduced due to information redundancy. In classification scheme F, the feature variables with high correlation were removed, which avoided the interference of overlapping redundant information on wetland information extraction and classification. The classification results are shown in Figure 5.

Wetland Information Extraction and Classification Results
Utilizing the optimal feature variable sets and RF algorithm, the spatial distribution of wetlands in 2015 and 2021 was obtained. As shown in Figure 6a, in 2015, the area of the natural wetland was 986.99 km 2 , while that of the artificial wetland was 1334.33 km 2 . In the natural wetland, the mud flat had the largest area of 432.90 km 2 , accounting for 6.16 % of the total area, which was mainly distributed in the northern Hekou district, the coastal area of Kenli district, and eastern Dongying district. The shallow sea had the second largest area of 353.47 km 2 , accounting for 5.03% of the total area, which was mostly located in the northern Hekou district, coastal areas in the Kenli district, and eastern Dongying district. The area of swamp was 129.79 km 2 , accounting for 1.85 %, which was mostly distributed in the estuary of the Yellow River in the northeast of Kenli district. The river had the smallest area of 70.82 km 2 , accounting for only 1.01% of the total area. In the artificial wetland, the reservoirs and pond had the largest area of 926.30 km 2 , accounting for 13.19%, followed by salt pans, with an area of 274.66 km 2 , accounting for 3.91%, which was mainly distributed in the northern Hekou district, Kenli district, Dongying district, and eastern Guangrao county. The paddy field area was 133.37 km 2 , accounting for only 1.90% of the total area, which was largely distributed in the central and northeastern Kenli district.
As shown in Figure 6b, in 2021, the area of the natural wetland was 1239.93 km 2 , while that of the artificial wetland was 1573.22 km 2 . In the natural wetland, the mud flat had the largest area of 658.23 km 2 , accounting for 9.36% of the total area, which was mostly distributed in the northern Hekou district, the coastal area of Kenli district, and the eastern Dongying district. The river had the second largest area of 290.97 km 2 , accounting for 4.14%, which was composed of the Yellow River and the Guangli River. The area of swamp was 155.18 km 2 , accounting for 2.21%, which was mostly distributed in the northeast of the YRD Nature Reserve in the Kenli district. The shallow sea had the smallest area of 135.55 km 2 , accounting for only 1.93%, which was mostly distributed along the coastline. In the artificial wetland, reservoirs, and ponds had the largest area of 1286.21 km 2 , accounting for 18.29% of the total area. The salt pans had the second largest area of 151.22 km 2 , accounting for 2.15% of the total area, which was mostly distributed in coastal areas. The paddy field had the smallest area of 135.79 km 2 , accounting for only 1.93% of the total area, which was mainly distributed in the rivers and lakes of the Kenli district.

Spatial distribution of wetland changes
The change patterns of the primary classification of wetlands in the YRD during 2015-2021 are shown in Figure 7. The YRD had the largest area of land cover as unchanged (5871.75 km 2 ), followed by that of the new-born area, with an area of 674.82333 km 2 . The internal conversion area of wetland occupied 288.98 km 2 and the extinction group had the smallest area of 239.59 km 2 . The extinction area of wetlands had the smallest area of 239.59 km 2 . Among them, the unchanged area, internal conversion area, and extinction area covered the largest area in the Kenli district with 1826.29 km 2 , 118.33 km 2 and 83.12 km 2 , respectively. In the Hekou district, the unchanged area, internal conversion area, and extinction area covered an area of 1769.66 km 2 , 97.90 km 2 and 65.46 km 2 , respectively. However, the new-born wetland had the largest area of 242.69 km 2 . In Guangrao county, the unchanged area, the internal conversion area, extinction area, and new-born area covered the smallest area of 276.83 km 2 , 11.44 km 2 , 5.09 km 2 and 23.70 km 2 . Table 7 is the comparison table of wetland change types in Figure 7. As shown in Figure 6b, in 2021, the area of the natural wetland was 1239.93 km 2 , while that of the artificial wetland was 1573.22 km 2 . In the natural wetland, the mud flat had the largest area of 658.23 km 2 , accounting for 9.36% of the total area, which was mostly distributed in the northern Hekou district, the coastal area of Kenli district, and the eastern Dongying district. The river had the second largest area of 290.97 km 2 , accounting for 4.14%, which was composed of the Yellow River and the Guangli River. The area of swamp was 155.18 km 2 , accounting for 2.21%, which was mostly distributed in the northeast of the YRD Nature Reserve in the Kenli district. The shallow sea had the smallest area of 135.55 km 2 , accounting for only 1.93%, which was mostly distributed along the coastline. In the artificial wetland, reservoirs, and ponds had the largest area of 1286.21 km 2 , accounting for 18.29% of the total area. The salt pans had the second largest area of 151.22 km 2 , accounting for 2.15% of the total area, which was mostly distributed in coastal areas. The paddy field had the smallest area of 135.79 km 2 , accounting for only 1.93% of the total area, which was mainly distributed in the rivers and lakes of the Kenli district.

Spatial Distribution of Wetland Changes
The change patterns of the primary classification of wetlands in the YRD during 2015-2021 are shown in Figure 7. The YRD had the largest area of land cover as unchanged (5871.75 km 2 ), followed by that of the new-born area, with an area of 674.82333 km 2 . The internal conversion area of wetland occupied 288.98 km 2 and the extinction group had the smallest area of 239.59 km 2 . The extinction area of wetlands had the smallest area of 239.59 km 2 . Among them, the unchanged area, internal conversion area, and extinction area covered the largest area in the Kenli district with 1826.29 km 2 , 118.33 km 2 and 83.12 km 2 , respectively. In the Hekou district, the unchanged area, internal conversion area, and extinction area covered an area of 1769.66 km 2 , 97.90 km 2 and 65.46 km 2 , respectively. However, the new-born wetland had the largest area of 242.69 km 2 . In Guangrao county, the unchanged area, the internal conversion area, extinction area, and new-born area covered the smallest area of 276.83 km 2 , 11.44 km 2 , 5.09 km 2 and 23.70 km 2 . Table 7 is the comparison table of wetland change types in Figure 7. Internal conversion area; Ⅲ Extinction area; Ⅳ New-born area. Table 7. Comparison of primary classification change patterns.

Number
Type Meaning Ⅰ Unchanged area Land cover did not change Ⅱ Internal conversion area Conversion between natural/constructed wetlands Ⅲ Extinction area Conversion of natural/constructed wetlands to non-wetlands Ⅳ New-born area Conversion of non-wetlands to natural / constructed wetlands  The change patterns of the secondary classification of wetlands during 2015-2021 are shown in Figure 8. The zone of unclassified land converted to reservoirs and ponds had the largest area of 142.29 km 2 . The zone of salt pans converted to reservoirs and ponds had the second largest area of 139.79 km 2 . The conversion area from shallow sea to mud flat covered an area of 129.89 km 2 , while that of the conversion area from mud flat to shallow sea was only 1.79 km 2 . For the internal conversion area of wetland, the zone of salt pans converted to reservoirs and ponds in the Hekou district had the largest area of 70.62 km 2 , accounting for 50.52 % of this conversion-type area, followed by that of the Kenli district, with an area of 28.41 km 2 , accounting for 20.32% of this conversion-type area. The zone of shallow sea converted to mud flat in the Kenli district had the largest area of 68.94 km 2 , accounting for 53.08% of this conversion-type area, followed by that of the Hekou district, with an area of 54.48 km 2 , accounting for 41.94%. In the extinction area of wetland, the zone of mud flat converted to construction land in the Dongying district had the largest area of 55.80 km 2 , accounting for 34.91%, followed by that of the Hekou district, accounting for 30.87%. The swamp zone converted to dryland in the Kenli district had the second largest area of 43.32 km 2 , accounting for the largest proportion of 61.56%. In the new-born wetland area, the zone of unused land converted to reservoirs and ponds had the largest area, while that of the Hekou district and Kenli district accounted for 46.13% and 28.85%, respectively. The dryland zone converted to paddy field had the second largest area of 87.00 km 2 , while that of Kenli district and Lijin county accounted for 41.46% and 23.66 %, respectively. Secondary classification coding of land cover and wetland types in the YRD is shown in Table 8.
As shown in Figure 9a,b, the area of the new-born paddy field that derived from dryland was 111.63 km 2 , accounting for 77.93%. The extinction area of the paddy field was 88.31 km 2 , of which 49.06% converted to dryland. The area of the new-born shallow sea was 26.99 km 2 , while that of the extinction area of the shallow sea was 51.87 km 2 . The areas of new-born salt pans and reservoirs and ponds were 22.35 km 2 and 316.27 km 2 , and the conversion area percentages from unused land were 41.13% and 44.99%, respectively. The areas of extinction salt pans and reservoirs and reservoirs and ponds were 18.69 km 2 and 42.63 km 2 , which mostly changed into construction land, with conversion ratios of 51.35% and 59.60%, respectively. As shown in Figure 9c, 50.45% of paddy fields, 77.55% of salt pans, 82.35% of rivers, 47.97% of mud flats, and 46.14 % of swamps were converted to reservoirs and ponds, and 62.97% of shallow seas were converted to mud flats. Notes: The first two digits and the last two digits represent the land cover types shown in Table 8, and the third digit 0 represents the conversion process. For example, 01001 represents the unchanged area of paddy fields, and 01003 represents paddy fields converted to salt pan areas. new-born wetland area, the zone of unused land converted to reservoirs and ponds had the largest area, while that of the Hekou district and Kenli district accounted for 46.13% and 28.85%, respectively. The dryland zone converted to paddy field had the second largest area of 87.00 km 2 , while that of Kenli district and Lijin county accounted for 41.46% and 23.66 %, respectively. Secondary classification coding of land cover and wetland types in the YRD is shown in Table 8.  conversion area percentages from unused land were 41.13% and 44.99%, respectively. The areas of extinction salt pans and reservoirs and reservoirs and ponds were 18.69 km 2 and 42.63 km 2 , which mostly changed into construction land, with conversion ratios of 51.35% and 59.60%, respectively. As shown in Figure 9c, 50.45% of paddy fields, 77.55% of salt pans, 82.35% of rivers, 47.97% of mud flats, and 46.14 % of swamps were converted to reservoirs and ponds, and 62.97% of shallow seas were converted to mud flats.

Gravity Center Migrations for Different Wetlands
The migration trajectory of the gravity center for wetlands could reflect the spatial heterogeneity and bias of the change patterns of different wetland types in varying periods [42][43][44]. As shown in Figure 10, the gravity center of each wetland type in the YRD was distributed along the Yellow River Basin from 2015 to 2021. Among them, the gravity center of the paddy field migrated from northeast to southwest, while that of the shallow sea and river moved from southeast to northwest. The gravity center of reservoirs and pond migrated from northwest to southeast, while that of mud flat moved from south to north. The long half axis of the standard deviation ellipse was significantly larger than the short half axis, indicating that the change direction of the wetlands in the YRD was more remarkable from 2015 to 2021. The increasing trend of wetlands in the southwest and northeast was significantly greater than that in the other directions. The migration directions of gravity centers of different wetland types in the YRD from 2015 to 2021 are shown in Table 9.
sea and river moved from southeast to northwest. The gravity center of reservoirs and pond migrated from northwest to southeast, while that of mud flat moved from south to north. The long half axis of the standard deviation ellipse was significantly larger than the short half axis, indicating that the change direction of the wetlands in the YRD was more remarkable from 2015 to 2021. The increasing trend of wetlands in the southwest and northeast was significantly greater than that in the other directions. The migration directions of gravity centers of different wetland types in the YRD from 2015 to 2021 are shown in Table 9.  The spatial distribution of various wetlands in the YRD was affected by both natural and human factors. Correlation analysis was applied to measure the correlations between  The spatial distribution of various wetlands in the YRD was affected by both natural and human factors. Correlation analysis was applied to measure the correlations between two or more associated variables [45,46]. Based on the method for the simple linear regression analysis, the correlation between the wetland area and the mean value of influencing factors was analyzed under the confidential level of 95%. The correlation coefficient R is shown in Table 10. As shown in Figure 11a, there was high correlation between wetland area and FVC, with r = −0.961. The primary reason was that zones with higher FVC were mostly dominated by dryland and woodland, which were mainly distributed in the northwest of Kenli district on both sides of the Yellow River Basin and inland, the southern part of Hekou district, Lijin county, Dongying district, and the western part of Guangrao county. The wetland was mostly distributed in the northern part of Hekou district, Kenli district, and the eastern part of Dongying district. The wetlands were mainly composed of Suaeda salsa, which is a typical salt-accumulating halophyte with lower vegetation coverage. As shown in Figure 11b, the correlation between wetland area and precipitation was moderate, with r = 0.850. The wetland area increased with the increase in precipitation. The primary reason was that more precipitation would increase soil moisture content and river runoff supply, and promote the expansion of wetland area. As shown in Figure 11c, the correlation between wetland area and coastal saline soil area was high, with r = 0.960. The reason was that aquaculture reservoirs, ponds, and salt pans were mainly concentrated in coastal areas with severe salinization. As shown in Figure 11d, there was a moderate correlation between wetland area and GDP, with r = 0.847. The primary reason was that, since 2009, economic development has been beneficial to the rational development and utilization of wetland resources, which would protect wetland resources.
As shown in Figure 11a, there was high correlation between wetland area and FVC, with r = −0.961. The primary reason was that zones with higher FVC were mostly dominated by dryland and woodland, which were mainly distributed in the northwest of Kenli district on both sides of the Yellow River Basin and inland, the southern part of Hekou district, Lijin county, Dongying district, and the western part of Guangrao county. The wetland was mostly distributed in the northern part of Hekou district, Kenli district, and the eastern part of Dongying district. The wetlands were mainly composed of Suaeda salsa, which is a typical salt-accumulating halophyte with lower vegetation coverage. As shown in Figure 11b, the correlation between wetland area and precipitation was moderate, with r = 0.850. The wetland area increased with the increase in precipitation. The primary reason was that more precipitation would increase soil moisture content and river runoff supply, and promote the expansion of wetland area. As shown in Figure 11c, the correlation between wetland area and coastal saline soil area was high, with r = 0.960. The reason was that aquaculture reservoirs, ponds, and salt pans were mainly concentrated in coastal areas with severe salinization. As shown in Figure 11d, there was a moderate correlation between wetland area and GDP, with r = 0.847. The primary reason was that, since 2009, economic development has been beneficial to the rational development and utilization of wetland resources, which would protect wetland resources.

Interaction Factors
The combined impacts (Figure 12c, d) of the change of wetlands in YRD had a higher q-value than individual impacts (Figure 12a, b) on the change of wetlands. Approximately 36 pairwise combinations with three different types of interactions were formed between the 9 driving factors. Most drivers exhibited non-linear enhancement when interacting with the change of wetlands, followed by bi-factor enhancement unique-variable weakening interactions. As shown in Figure 12c, solar radiation, when interacting with accumulated temperature, showed bi-factor enhancement and explained 86.3% of the change in the natural wetland. Bi-factor enhancement occurred when the combined influence (solar radiation + accumulated temperature = 1.253) of the two variables was greater than individual (solar radiation=0.448 and accumulated temperature=0.805) and interaction influences (solar radiation and accumulated temperature = 0.863) on the response variable. Non-linear enhancement occurred when the interaction influence (FVC and temperature = 0.797) was greater than that of the individual (FVC=0.098 and temperature=0.648) and combined (FVC + temperature = 0.746) variables. Unique variable weakening occurred when individual variable (precipitation=0.531) influence was greater than either combined (precipitation + slope = 0.594) or interaction (precipitation + slope = 0.512) effects on a response variable [47].
As shown in Figure 12d, FVC interaction with DEM and accumulated temperature caused bi-factor enhancement and explained 89.6% and 82.0% of the change in artificial wetlands, respectively. Socioeconomic variable interactions, such as GDP density and

Interaction Factors
The combined impacts (Figure 12c, d) of the change of wetlands in YRD had a higher q-value than individual impacts (Figure 12a, b) on the change of wetlands. Approximately 36 pairwise combinations with three different types of interactions were formed between the 9 driving factors. Most drivers exhibited non-linear enhancement when interacting with the change of wetlands, followed by bi-factor enhancement unique-variable weakening interactions. As shown in Figure 12c, solar radiation, when interacting with accumulated temperature, showed bi-factor enhancement and explained 86.3% of the change in the natural wetland. Bi-factor enhancement occurred when the combined influence (solar radiation + accumulated temperature = 1.253) of the two variables was greater than individual (solar radiation = 0.448 and accumulated temperature = 0.805) and interaction influences (solar radiation and accumulated temperature = 0.863) on the response variable. Nonlinear enhancement occurred when the interaction influence (FVC and temperature = 0.797) was greater than that of the individual (FVC=0.098 and temperature=0.648) and com-bined (FVC + temperature = 0.746) variables. Unique variable weakening occurred when individual variable (precipitation = 0.531) influence was greater than either combined (precipitation + slope = 0.594) or interaction (precipitation + slope = 0.512) effects on a response variable [47].

Advantages of the Classification Method of Wetlands
Most previous studies have realized the classification of wetlands by extracting many feature variables. However, classification accuracy and efficiency were greatly affected by many feature variables. Therefore, it was crucial to optimize the feature variables according to local conditions before classification. The RF algorithm had greater advantages in determining the optimal subset of feature variables because it could automatically judge As shown in Figure 12d, FVC interaction with DEM and accumulated temperature caused bi-factor enhancement and explained 89.6% and 82.0% of the change in artificial wetlands, respectively. Socioeconomic variable interactions, such as GDP density and population density, explained 69.7% of the change in artificial wetlands. Climate and socioeconomic variable interactions, such as GDP density and accumulated temperature, population density, and precipitation, explained 86.2% and 76.8% of the change in artificial wetland, respectively. On the other hand, DEM showed non-linear enhancement when interacting with slope because of their weak individual influence on the change of artificial wetlands.
It can be seen that solar radiation and accumulated temperature had larger explanatory power on the change of natural wetlands, while FVC and DEM contributed greatly to the changes of artificial wetlands [48][49][50].

Advantages of the Classification Method of Wetlands
Most previous studies have realized the classification of wetlands by extracting many feature variables. However, classification accuracy and efficiency were greatly affected by many feature variables. Therefore, it was crucial to optimize the feature variables according to local conditions before classification. The RF algorithm had greater advantages in determining the optimal subset of feature variables because it could automatically judge the importance of variables and give the dependence among variables. In this paper, the optimal feature variable sets in the YRD were determined utilizing the RF method, which could provide important references for relative investigations of wetland in the Yellow River delta. Then, six classification schemes were designed and tested, and it was found that scheme F had the best classification accuracy. This was because scheme F eliminated the feature variables with a low contribution rate to classification through the RF method, so that the remaining feature variables had low information redundancy, which was conductive to enhance the classification accuracy. The classification accuracy was reduced when all feature variables were utilized to extract wetland information based on the RF method. This was consistent with the results of Song et al. [51]. Scheme C had better classification accuracies than schemes E and D because there was low redundancy between spectral features and the red-edge index. In addition, the red-edge index was a better indicator to reflect the grown condition of vegetation, which was conductive to enhancing the classification accuracy of wetlands [52]. The overall classification accuracy based on the RF method and the optimal feature variable sets was higher than other methods (SVM classification [53], supervised classification [54]), which could provide technical references for related research in the YRD and other regions.

Causes of Wetland Change Pattern
The wetlands in the YRD showed coastal and riverside distribution characteristics. The reason was that a large area of land extending to the ocean had formed mud flats and shallow waters. In addition, reservoirs, ponds, and salt pans were mostly developed in coastal or riverside regions, with severe salinization that were not conducive to the growth of ordinary crops. To ensure economic income, marine aquaculture reservoirs, ponds, and salt pans are vigorously developed in the region.
During 2015-2021, the new-born and distinction wetlands mainly occurred in the coastal areas. The reasons were that the northern Hekou district and the eastern coast of Kenli district had large areas of shallow waters, mud flats, reservoirs, ponds, and other wetland types, and the coastal areas of the YRD had dramatic alternation of the land-sea boundary. Due to the continuous water and sediment transport of the Yellow River in the northeast of the Kenli district, the delta showed an expansion trend toward the ocean. In recent years, more unused land has been converted to paddy fields, rivers, reservoirs, ponds, and salt pans to improve the regional economy in coastal areas.

Conclusions
Based on Sentinel-2 images and RF methods, the optimal feature variable sets for wetland classification in the YRD were determined, and then the spatial distribution of wetlands in 2015 and 2021 were obtained. Finally, the change pattern and driving mechanism of wetlands during 2015-2021 were analyzed based on the gravity center model and geographical detector. The major conclusions were as follows: (1) Vegetation index and water index have a positive impact on wetland information extraction and classification accuracy. In addition, the red-edge index made greater contributions to enhance classification accuracy.
(2) The Random Forest algorithm could be used to optimize the feature variables, and remove the redundancy among feature variables. Scheme F, based on optimal feature variable sets, had the best classification accuracy of wetlands, which could provide references for the investigation of wetlands in the YRD. (3) During 2015-2021, a large area of natural wetland in the YRD was transformed into an artificial wetland. The wetlands showed an overall development direction of "northwest-southeast" along the Yellow River. (4) Wetland changes in the YRD were affected by both natural and human factors. The interaction between FVC and accumulated temperature had the largest explanatory power of the change in the natural wetland area. The interaction between solar radiation and DEM had the largest explanatory power for the change in the artificial wetland area.
In this paper, the study period (2015-2021) was not long enough to reveal the evolutionary process of the wetland in the YRD. In further studies, longer time series remote sensing images would be applied to investigate the spatio-temporal change patterns of wetlands. Moreover, the landscape pattern index should be adopted to analyze the spatial structure characteristics of different types of wetlands.