Tidal Flat Extraction and Change Analysis Based on the RF-W Model: A Case Study of Jiaozhou Bay, East China

Coastal tidal flats are important ecological resources. As the dividing line between marine and terrestrial ecosystems, tidal flats provide a large number of ecosystem services. However, with the excessive development of coastal areas, tidal flat resources have been drastically reduced, leading to the deterioration of coastal ecosystems. There is an urgent need to acquire accurate information on the changes in tidal flat resources. This research proposes a tidal flat extraction model (RF-W model) that combines the random forest (RF) method and waterline method, which aims to improve the accuracy of tidal flat extraction. This method can effectively eliminate the shortcomings of the RF method in determining the lower boundary of tidal flats and those of the waterline method in distinguishing river channels and tidal flats. The tidal flat extraction of Qingdao Jiaozhou Bay in 2020 is performed as an example of the model. The results show that the user’s and producer’s accuracies of the RF-W model were both the highest, indicating that the improved model can accurately extract tidal flat information. Then, we used the RF-W model to extract tidal flat information for Jiaozhou Bay in seven periods (1990, 1995, 2000, 2005, 2010, 2015, and 2020) and to study the spatiotemporal changes in the tidal flats and influencing factors from 1990 to 2020. The tidal flat area of Jiaozhou Bay showed an overall downward trend before 2015, and the area decreased by 21.9 km2, with a reduction in the rate of approximately 1.1%/year. After 2015, the tidal flat area rebounded slightly. The overall change in Jiaozhou Bay showed reclamation and expansion toward the sea. The reduction in the sand content of the rivers entering the sea, reclamation and cultivation, and land development were the main factors contributing to the reduction in the tidal flat area in Jiaozhou Bay. In addition, sea level rise due to climate warming is a long-term potential factor.


Introduction
Tidal flats refer to the tidal flooding zone between the high tide line and the low tide line in the coastal zone, and they are mostly distributed in estuaries [1]. As one of the most widely distributed ecosystems in coastal areas, tidal flats are the boundary between marine and terrestrial ecosystems and provide many ecosystem service functions [2]. For example, tidal flats support biodiversity and provide a habitat for a variety of animals and plants [3]; as a natural barrier for marine and terrestrial ecosystems, they can absorb multiple pollutants [4] and are an important part of the "blue carbon" ecosystem, playing a role in carbon sequestration [5]. However, with the rapid development of coastal areas and the continuous movement of city centers to coastal areas, tidal flats are sharply shrinking, disrupting the balance of the ecosystem and making tidal flats one of the most vulnerable areas on the planet. From 1989 to 2014, the area of tidal flats in the northern gulf near the Mong Cai watershed, Vietnam, decreased by 208 km 2 due to human activities and natural factors [6]. From 1990 to 2017, the area of tidal flats in the Yangtze River Delta decreased by 324 km 2 at a rate of approximately 21.7 km 2 /year [7]. In the past 50 years, the tidal flats along the Yellow Sea have been reduced by more than 50%, resulting in increased pollution, algae blooms, and loss of many animals and plants [8]. Therefore, there is an urgent need to accurately measure the geographic location and area of undeveloped tidal flats and to understand the extent and causes of changes in developed tidal flats to protect and adjust land use planning [9,10].
Due to the poor accessibility of tidal flat regions, field measurement data are scarce. With the rapid development of remote sensing technology, the richness of remote sensing datasets, and the diversity of data processing methods, it is possible to accurately extract tidal flat information [11,12]. Many studies have used edge detection or threshold segmentation to extract water lines and then obtain high tide lines and low tide lines to determine the tidal flat range and study its temporal and spatial changes [13][14][15]. For example, Qiu et al. (2019) used the normalized difference vegetation index (NDVI) combined with a threshold segmentation method to obtain the water and vegetation lines to extract tidal flats in the Yangtze River Delta and analyze their changes from 1975 to 2017 [7]. Sagar et al. (2017) used a large-scale tide model combined with the normalized difference water index (NDWI) index to survey and map the coastal beaches of Australia [16]. Behling et al. (2018) used modified NDWI (mNDWI) time-series data to monitor changes in the coastline and tidal flats of the Gulf of Namibia from 1984 to 2014 [17]. Based on Landsat imagery, Wang et al. (2018) used the mNDWI and other vegetation indexes (NDVI and the enhanced vegetation index, EVI) to extract water and vegetation lines, to obtain an annual average water frequency map to determine the tidal flat range, and to generate a map of the coastal tidal flat changes in eastern China from 1986 to 2016 [18]. Traditional extraction methods based on waterlines are affected by tidal conditions, residual surface water, and other environmental parameters, resulting in inaccurate waterline extraction [13]. Moreover, because the intertidal zone contains some salt-tolerant vegetation and tidal flat streams, the waterline method cannot distinguish vegetation from tidal flats well, resulting in a decrease in the accuracy of tidal flat extraction.
The intertidal zone is completely exposed at low tide, and it is continuously submerged at high tide. The tidal flat information extracted by remote sensing classification methods, such as visual interpretation [19], supervised classification [20], and object orientation [21], represents only the instantaneous state when an image is acquired. This phenomenon places high requirements on the time and quality of image acquisition, which hinders the rapid extraction and accurate analysis of tidal flat information. With the rapid development of machine learning and the continuous improvement in the spatial and temporal resolution of remote sensing images, an increasing number of studies have appeared in the fields of automatic feature extraction and image classification. As a machine learning algorithm, the random forest (RF) classification method can automate the feature extraction process, reduce work costs, and improve efficiency [10]. For example, Dubeau et al. (2017) used the RF classification method to classify the Dabs' wetland in Ethiopia, with an overall classification accuracy of 94.4% [22]. Gu et al. (2019) used the RF classification method to extract the information for the Wei-Ku oasis wetland on the edge of the Tarim River, with an overall accuracy of 90.1% [23]. Hou et al. (2020) used GF-1 images and combined out-of-bag (OOB) data to perform feature reduction and finally classified alpine wetlands with an RF classifier, with an overall classification accuracy of 90.1% [24]. Although the RF classification method has been widely used, the results of tidal flat information extraction are greatly affected by the tidal conditions at the time of data acquisition, which leads to shortcomings in the application of this method.
This study combines the advantages of the RF method and the waterline method in tidal flat extraction and proposes an improved RF-W tidal flat extraction model. First, multisource remote sensing images are used to extract tidal flat information in the study area through the RF classification method. Second, the enhanced water boundary index (EWBI) of each image is calculated, and the waterline is extracted through threshold segmentation and edge detection of the EWBI information from the image. A digital shoreline analysis system (DSAS) is used to generate a high tide point set and a low tide Remote Sens. 2021, 13, 1436 3 of 18 point set to determine the tidal flat range. Third, the lower boundary of the tidal flats obtained by the waterline method is used to correct the tidal flat information obtained by the RF classification method to obtain more accurate tidal flat information. Finally, the RF-W model proposed in this paper is used to extract the tidal flats in the Qingdao Jiaozhou Bay coastal zone at 5-year intervals (1990-2020) and to study their temporal and spatial changes and influencing factors. This research will provide methods and data support for sustainable management and ecological assessment of coastal zones.

Study Area and Data
Jiaozhou Bay is located on the coast of the Yellow Sea on the south coast of China's Shandong Peninsula (36 • 01 N~36 • 15 N, 120 • 03 E~120 • 25 E). Jiaozhou Bay is the foundation of the formation and development of Qingdao and an important carrier of economic development; it is known as the "Mother Bay of Qingdao" (Figure 1). The terrain along Jiaozhou Bay gradually decreases from southeast to northwest, with mountains in the south and southwest, plains in the north and northwest, and the Laoshan Mountains in the east. Jiaozhou Bay is a typical semienclosed bay with a total area of 446 km 2 , of which the wetland area is approximately 177.6 km 2 [25]. More than 10 rivers flow along the coast of the bay. Among them, the Dagu River, Moshui River, and Baisha River deposit large amounts of sediment in their estuaries, resulting in abundant tidal flat resources in Jiaozhou Bay, and these tidal flats are mainly distributed along the northern and northwestern coasts of Jiaozhou Bay (Figure 1). The tidal flats in this study mainly refer to the intertidal zone, which is the area between the average high tide level and low tide level of the spring tide. mentation and edge detection of the EWBI information from the image. A digital shor analysis system (DSAS) is used to generate a high tide point set and a low tide poin to determine the tidal flat range. Third, the lower boundary of the tidal flats obtaine the waterline method is used to correct the tidal flat information obtained by the RF sification method to obtain more accurate tidal flat information. Finally, the RF-W m proposed in this paper is used to extract the tidal flats in the Qingdao Jiaozhou Bay co zone at 5-year intervals (1990-2020) and to study their temporal and spatial changes influencing factors. This research will provide methods and data support for sustain management and ecological assessment of coastal zones.

Study Area and Data
Jiaozhou Bay is located on the coast of the Yellow Sea on the south coast of Ch Shandong Peninsula (36°01′N∼36°15′N, 120°03′E∼120°25′E). Jiaozhou Bay is the fou tion of the formation and development of Qingdao and an important carrier of econ development; it is known as the "Mother Bay of Qingdao" (Figure 1). The terrain a Jiaozhou Bay gradually decreases from southeast to northwest, with mountains in south and southwest, plains in the north and northwest, and the Laoshan Mountai the east. Jiaozhou Bay is a typical semienclosed bay with a total area of 446 km 2 , of w the wetland area is approximately 177.6 km 2 [25]. More than 10 rivers flow along the of the bay. Among them, the Dagu River, Moshui River, and Baisha River deposit amounts of sediment in their estuaries, resulting in abundant tidal flat resourc Jiaozhou Bay, and these tidal flats are mainly distributed along the northern and n western coasts of Jiaozhou Bay ( Figure 1). The tidal flats in this study mainly refer t intertidal zone, which is the area between the average high tide level and low tide lev the spring tide. The remote sensing data used in the research were mainly Sentinel-2 images Landsat images. From 1990 to 2020, a total of 102 images were obtained to extract tida information ( Figure 2). Sentinel-2 images (L1C) from the official website of the Euro Space Agency (ESA) (https://scihub.copernicus.eu/dhus/#/home; accessed on 12 Sep ber 2020) and the United States Geological Survey (USGS) (https://glovis.usgs.gov/ The remote sensing data used in the research were mainly Sentinel-2 images and Landsat images. From 1990 to 2020, a total of 102 images were obtained to extract tidal flat information ( Figure 2). Sentinel-2 images (L1C) from the official website of the European Space Agency (ESA) (https://scihub.copernicus.eu/dhus/#/home; accessed on 12 September 2020) and the United States Geological Survey (USGS) (https://glovis.usgs.gov/app; accessed on 12 September 2020) cover 13 bands from visible light to shortwave infrared, with a maximum spatial resolution of 10 m and a revisit period of 5 days. Landsat images were obtained from USGS, and the images included those captured by the OLI sensors and TM sensors. The OLI sensor included 9 multispectral bands and 1 panchromatic band, and the spatial resolution after fusion can reach 15 m; the TM sensor included 6 multispectral bands with 30 m spatial resolution and the revisit period was 16 days for both sensors. The selected images had limited cloud coverage and high quality in the bay region. In addition, several Chinese GF-2 remote sensing images in 2020 were also included, with a spatial resolution of up to 1 m, which assisted in the collection of classified samples together with the high-resolution images in Google Earth. accessed on 12 September 2020) cover 13 bands from visible light to shortwav with a maximum spatial resolution of 10 m and a revisit period of 5 days. Land were obtained from USGS, and the images included those captured by the O and TM sensors. The OLI sensor included 9 multispectral bands and 1 panchrom and the spatial resolution after fusion can reach 15 m; the TM sensor include spectral bands with 30 m spatial resolution and the revisit period was 16 day sensors. The selected images had limited cloud coverage and high quality in gion. In addition, several Chinese GF-2 remote sensing images in 2020 were also with a spatial resolution of up to 1 m, which assisted in the collection of classifie together with the high-resolution images in Google Earth. The atmospheric correction needs to be performed on the acquired image true surface reflectance data. Sentinel-2 (L1C) images were processed by the 2.8.0 algorithm [26] published by the ESA and converted into L2A images, an images were resampled to 10 m resolution through SNAP software and conv ENVI standard format. The Sen2Cor algorithm was a plug-in released by the E mospheric correction of Sentinel series images. The core of the algorithm was pheric radiation transmission model libRadtran [27]. This model was suitable atmospheric conditions, including the radiation, irradiance, and photochemical sun and part of Earth's spectrum. The Landsat image was corrected by the FL mospheric correction model in ENVI [28], and finally, the image to be classi Jiaozhou Bay area was obtained by cropping.

Technical Process of Tidal Flats Extraction
When the RF classification method was used to extract tidal flats, the low were not accurately determined, which affects the accuracy of the extracted ti formation. The waterline method had certain advantages in determining the low The atmospheric correction needs to be performed on the acquired images to obtain true surface reflectance data. Sentinel-2 (L1C) images were processed by the Sen2Cor-2.8.0 algorithm [26] published by the ESA and converted into L2A images, and then the images were resampled to 10 m resolution through SNAP software and converted into ENVI standard format. The Sen2Cor algorithm was a plug-in released by the ESA for atmospheric correction of Sentinel series images. The core of the algorithm was the atmospheric radiation transmission model libRadtran [27]. This model was suitable for various atmospheric conditions, including the radiation, irradiance, and photochemical flux of the sun and part of Earth's spectrum. The Landsat image was corrected by the FLAASH atmospheric correction model in ENVI [28], and finally, the image to be classified in the Jiaozhou Bay area was obtained by cropping.

Technical Process of Tidal Flats Extraction
When the RF classification method was used to extract tidal flats, the low tide lines were not accurately determined, which affects the accuracy of the extracted tidal flat information. The waterline method had certain advantages in determining the lower boundary of tidal flats [29]. This research proposed a tidal flat extraction model (RF-W model) that combined the RF method and waterline method to achieve rapid and highprecision extraction of tidal flats. The specific technical process is shown in Figure 3. (1) Determine the optimal segmentation parameters of the images to segment the images and select samples based on field survey data and high-resolution images; (2) build an RF classification feature database and select suitable features; (3) extract tidal flat information through the RF algorithm; (4) based on time series data, use the waterline method to extract the high tide line and low tide line, and finally determine the lower boundary of the tidal flats; (5) use the lower boundary of the tidal flats determined by the waterline method to modify the RF classification algorithm to extract the tidal flats and generate a thematic map of tidal flats. In this technical process, the extraction of tidal flats in Jiaozhou Bay in 2020 was taken as an example to introduce the specific extraction process. the optimal segmentation parameters of the images to segment the images and select samples based on field survey data and high-resolution images; 2) build an RF classification feature database and select suitable features; 3) extract tidal flat information through the RF algorithm; 4) based on time series data, use the waterline method to extract the high tide line and low tide line, and finally determine the lower boundary of the tidal flats; 5) use the lower boundary of the tidal flats determined by the waterline method to modify the RF classification algorithm to extract the tidal flats and generate a thematic map of tidal flats. In this technical process, the extraction of tidal flats in Jiaozhou Bay in 2020 was taken as an example to introduce the specific extraction process.

Image Segmentation
A multiscale algorithm combined with a spectral difference segmentation model was used to segment the image. The multiscale image segmentation algorithm was a region merging segmentation algorithm based on the minimum heterogeneity [30]. When segmenting, we must first analyze the spectral features and shape features of the objects in

Image Segmentation
A multiscale algorithm combined with a spectral difference segmentation model was used to segment the image. The multiscale image segmentation algorithm was a region merging segmentation algorithm based on the minimum heterogeneity [30]. When segmenting, we must first analyze the spectral features and shape features of the objects in the image, then set the segmentation parameters by a certain proportion according to the importance of these features, gradually merge the homogeneous pixels from bottom to top, and reduce the area of heterogeneous regions. Spectral difference segmentation was based on multiscale segmentation, using the spectral reflectance difference between different features to optimize the results of multiscale segmentation to achieve the best segmentation effect [31]. Therefore, when performing image segmentation, the most important aspect is to find the optimal segmentation scale and segmentation parameters.
This study was based on the estimation of the scale parameter (ESP) plug-in in eCognition 9.0 software to calculate the optimal segmentation scale of the image to be classified [32]. This tool can calculate the impact of different segmentation scales on the local variance (LV) of the image object, expressed as the rate of change in LV (ROC-LV). When ROC-LV has an extreme value, the segmentation scale corresponding to this point is the local optimal segmentation scale. Figure 4 shows that there are multiple local optimal segmentation scales. A segmentation scale that is too small will fragment the tidal flats too much, and a segmentation scale that is too large will cause other land types to be included in the tidal flat segmentation results. Therefore, combined with the completeness of the tidal flat information in the segmentation results, the segmentation scale was set to 90; the shape factor and the spectral factor were 0.2 and 0.8, respectively; and the smoothness and compactness were both 0.5 to split the image. After the multiscale segmentation was completed, the result of spectral difference segmentation was optimized to obtain the final segmentation result. top, and reduce the area of heterogeneous regions. Spectral difference segmentation was based on multiscale segmentation, using the spectral reflectance difference between different features to optimize the results of multiscale segmentation to achieve the best segmentation effect [31]. Therefore, when performing image segmentation, the most important aspect is to find the optimal segmentation scale and segmentation parameters. This study was based on the estimation of the scale parameter (ESP) plug-in in eCognition 9.0 software to calculate the optimal segmentation scale of the image to be classified [32]. This tool can calculate the impact of different segmentation scales on the local variance (LV) of the image object, expressed as the rate of change in LV (ROC-LV). When ROC-LV has an extreme value, the segmentation scale corresponding to this point is the local optimal segmentation scale. Figure 4 shows that there are multiple local optimal segmentation scales. A segmentation scale that is too small will fragment the tidal flats too much, and a segmentation scale that is too large will cause other land types to be included in the tidal flat segmentation results. Therefore, combined with the completeness of the tidal flat information in the segmentation results, the segmentation scale was set to 90; the shape factor and the spectral factor were 0.2 and 0.8, respectively; and the smoothness and compactness were both 0.5 to split the image. After the multiscale segmentation was completed, the result of spectral difference segmentation was optimized to obtain the final segmentation result.

Build a Classification Feature Library
Before performing RF classification, a classification feature database must be constructed. To accurately distinguish between water bodies and tidal flats, a series of remote sensing indexes were calculated using ENVI 5.3 software to build a classification feature database, including the NDWI, automatic water body extraction index (AWElsh) [33], high-resolution water index (HRWI) [34], water index based on linear discriminant analysis (WI2015) [35] and an outcropping tidal flat extraction index (TFEI) [36]. Due to its high salt content, salt-tolerant vegetation often grows in tidal flats, and these vegetation types can be distinguished by the NDVI. Additionally, the sliding window size was set to 3*3, the window moving distance was set to 1, and the quantized gray level was set to 8 to calculate the gray level co-occurrence matrix to obtain 8 texture features of the image. The brightness index (BI), greenness index (GI), and wetness index (WI) can be obtained

Build a Classification Feature Library
Before performing RF classification, a classification feature database must be constructed. To accurately distinguish between water bodies and tidal flats, a series of remote sensing indexes were calculated using ENVI 5.3 software to build a classification feature database, including the NDWI, automatic water body extraction index (AWElsh) [33], high-resolution water index (HRWI) [34], water index based on linear discriminant analysis (WI2015) [35] and an outcropping tidal flat extraction index (TFEI) [36]. Due to its high salt content, salt-tolerant vegetation often grows in tidal flats, and these vegetation types can be distinguished by the NDVI. Additionally, the sliding window size was set to 3 × 3, the window moving distance was set to 1, and the quantized gray level was set to 8 to Remote Sens. 2021, 13, 1436 7 of 18 calculate the gray level co-occurrence matrix to obtain 8 texture features of the image. The brightness index (BI), greenness index (GI), and wetness index (WI) can be obtained through tasseled cap transformation [37,38]. The specific classification features and their calculation equation or description are shown in Table 1.
High-Resolution Water Index HRWI [34]  Note: B is the blue band of the images, G is the green band, R is the red band, NIR is the near-infrared band, SWIR 1 is the first shortwave infrared band, and SWIR 2 is the second shortwave infrared band.

RF Classification
The RF classification algorithm is a data-driven nonparametric machine learning algorithm composed of multiple classifications and regression tree (CART) decision trees [39]. By sampling a large number of training samples with replacement N times, N training sample sets were obtained, and m classification features (m ≤ M) were randomly selected from the total features M in each sample set. Then, the nodes of the decision tree were divided by complete splitting, and N decision trees were constructed. In the prediction stage, after each tree independently completes the prediction, a vote was made for the category of a new sample, and the prediction of the new sample was finally completed. Among many machine learning algorithms, the RF algorithm exhibited fast classification speed and high accuracy and was insensitive to noise, and suitable for a variety of classification scenarios.
A major advantage of the RF classification algorithm is that it can perform feature dimensionality reduction, reduce data redundancy, and improve classification efficiency. During each random sampling process, 1/3 of the data were not selected, and these unselected data were called OOB data. OOB data were used to calculate the feature variable importance score (VIS). The higher the VIS value is, the greater the contribution of the feature during classification. For any feature M r , its VIS value is calculated by Equation (1).
In Equation (1), errB 1 represents the OOB error value calculated using the OOB data for each decision tree in the RF classification algorithm. errB 2 means that for any feature, with other features remaining unchanged, noise interference was applied to the feature variable, and the OOB error value was recalculated. In Equation (2), P is the total number of OOB data, and Q is the number of correct classification results obtained by inputting OOB data into the random forest classifier. Figure 5 shows the VIS value calculated from the OOB data for each feature. WI2015 had the largest contribution to the extraction of tidal flat information, with a VIS value exceeding 0.05, and the contribution of texture information was the smallest. Using the forward sequence difference method, each feature was put into the RF classifier for classification according to the VIS value, and the overall classification accuracy was calculated, as shown in Figure 5. With the continuous increase in the number of features, the overall accuracy (OA) continues to increase, and with the addition of NDVI features, the OA tends to stabilize ( Figure 5). Therefore, considering the time cost, 6 feature variables (WI2015, SWIR 2 , SWIR 1 , NIR/B, WI, and GI) with VIS values exceeding 0.04 and NDVI were selected as the input features of the RF classification algorithm.
Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 19 Figure 5 shows the VIS value calculated from the OOB data for each feature. WI2015 had the largest contribution to the extraction of tidal flat information, with a VIS value exceeding 0.05, and the contribution of texture information was the smallest. Using the forward sequence difference method, each feature was put into the RF classifier for classification according to the VIS value, and the overall classification accuracy was calculated, as shown in Figure 5. With the continuous increase in the number of features, the overall accuracy (OA) continues to increase, and with the addition of NDVI features, the OA tends to stabilize ( Figure 5). Therefore, considering the time cost, 6 feature variables (WI2015, SWIR2, SWIR1, NIR/B, WI, and GI) with VIS values exceeding 0.04 and NDVI were selected as the input features of the RF classification algorithm.

Field Investigation and Reference Sample Selection
According to the existing land cover research and tidal flat-related research [40], combined with field investigation and research needs, the land cover types in the study area were divided into tidal flats and nontidal flats, and the nontidal flats were divided into land, water, salt-tolerant vegetation, salt pans, and offshore ponds ( Table 2). The highresolution satellite images from Google Earth and GF-2, as shown in Figure 6, indicated that all land cover types can be accurately identified. Taking land classification in 2020 as an example, referring to the high-resolution remote sensing images from Google Earth and field surveys, 2036 sample points were uniformly selected on the Sentinel-2 images in the study area-the quantity distribution of each type is shown in Table 2. Sample points were evenly distributed in the study area, 70% of which were used for classifier training and 30% for accuracy evaluation.

Field Investigation and Reference Sample Selection
According to the existing land cover research and tidal flat-related research [40], combined with field investigation and research needs, the land cover types in the study area were divided into tidal flats and nontidal flats, and the nontidal flats were divided into land, water, salt-tolerant vegetation, salt pans, and offshore ponds ( Table 2). The highresolution satellite images from Google Earth and GF-2, as shown in Figure 6, indicated that all land cover types can be accurately identified. Taking land classification in 2020 as an example, referring to the high-resolution remote sensing images from Google Earth and field surveys, 2036 sample points were uniformly selected on the Sentinel-2 images in the study area-the quantity distribution of each type is shown in Table 2. Sample points were evenly distributed in the study area, 70% of which were used for classifier training and 30% for accuracy evaluation.

Tidal Flat Extraction by the RF Method
The RF classification code was written in PyCharm integrated development ment based on libraries such as GDAL, Numpy, and Sklearn in Python 3.6. In t of the sample points were randomly selected as the training set, and the remai were used as the verification set. Since the RF classifier was composed of multiple trees, increasing the number of decision trees can improve the classification accu when the number of trees was very large, the computational cost will be high. T to improve accuracy and save costs, the number of decision trees was set to 110 number of classification features was set to 7. RF classification was then perform the classification results were imported into ENVI 5.3 software for post-classifica cessing, merging of adjacent homogeneous areas by cluster analysis, and modi misclassified parts through visual judgment to obtain the final classification inf for tidal flats and various nontidal flats. Tidal flat information was derived for th ing research on tidal flat extraction, and other nontidal flat information was change analysis.

Waterline Extraction Method
Due to the influence of tidal conditions, residual surface water, and reg strictions, some traditional edge detection and threshold segmentation method tracting waterlines were less effective in extracting waterlines in Jiaozhou Ba other areas [41][42][43]. Therefore, this study proposes an EWBI for waterline extr minimize this inaccuracy. As shown in Figure 7, the EWBI index proposed by this was compared with the land surface water index (LSWI) and mNDWI indexes, w better identify the water body boundary and improve the accuracy of waterline e = − + 0.02 = −

Tidal Flat Extraction by the RF Method
The RF classification code was written in PyCharm integrated development environment based on libraries such as GDAL, Numpy, and Sklearn in Python 3.6. In total, 70% of the sample points were randomly selected as the training set, and the remaining 30% were used as the verification set. Since the RF classifier was composed of multiple decision trees, increasing the number of decision trees can improve the classification accuracy, but when the number of trees was very large, the computational cost will be high. Therefore, to improve accuracy and save costs, the number of decision trees was set to 110, and the number of classification features was set to 7. RF classification was then performed. Then, the classification results were imported into ENVI 5.3 software for post-classification processing, merging of adjacent homogeneous areas by cluster analysis, and modifying the misclassified parts through visual judgment to obtain the final classification information for tidal flats and various nontidal flats. Tidal flat information was derived for the following research on tidal flat extraction, and other nontidal flat information was used for change analysis.

Waterline Extraction Method
Due to the influence of tidal conditions, residual surface water, and regional restrictions, some traditional edge detection and threshold segmentation methods for extracting waterlines were less effective in extracting waterlines in Jiaozhou Bay than in other areas [41][42][43]. Therefore, this study proposes an EWBI for waterline extraction to minimize this inaccuracy. As shown in Figure 7, the EWBI index proposed by this research was compared with the land surface water index (LSWI) and mNDWI indexes, which can better identify the water body boundary and improve the accuracy of waterline extraction.
G is the green band of the remote sensing images, NIR is the near-infrared band, and SWIR 1 is the shortwave infrared band 1. First, the EWBI was proposed based on the difference in spectral reflectan tidal flats and nontidal flats along Jiaozhou Bay, as shown in equation (3). Thi highlight the boundary of the waters along Jiaozhou Bay and yield an accura neous waterline. Second, the EWBI image was segmented with 0 as the thresh erate a binary image of land and water. Then, the Canny operator was use whether each pixel was an edge pixel to extract the waterline. Finally, the ex terline was revised through visual interpretation. To extract tidal flat inform outermost low tide waterline was selected as the lower boundary of the tidal fl innermost high tide waterline was selected as the upper boundary of the tidal Due to changes in daily tide levels, weather, and other factors, the tidal m the same area may vary, causing different water lines to cross each other. To the accurate tidal flat range, the DSAS plug-in in ArcGIS software was used to cross-section every 100 m along the coastline of the study area. The intersec cross-section and the waterline was considered the tidal point, which repres stantaneous tide level at the time of satellite overpass. The distance betwee point and the baseline was calculated for each section. The nearest and fart from the baseline were selected to form a high tide point set and low tide po spectively, and all high tide points and low tide points were connected to det First, the EWBI was proposed based on the difference in spectral reflectance between tidal flats and nontidal flats along Jiaozhou Bay, as shown in Equation (3). This index can highlight the boundary of the waters along Jiaozhou Bay and yield an accurate instantaneous waterline. Second, the EWBI image was segmented with 0 as the threshold to generate a binary image of land and water. Then, the Canny operator was used to detect whether each pixel was an edge pixel to extract the waterline. Finally, the extracted waterline was revised through visual interpretation. To extract tidal flat information, the outermost low tide waterline was selected as the lower boundary of the tidal flats, and the innermost high tide waterline was selected as the upper boundary of the tidal flats [14].

Remote Sens. 2021, 13, x FOR PEER REVIEW
Due to changes in daily tide levels, weather, and other factors, the tidal movement in the same area may vary, causing different water lines to cross each other. To determine the accurate tidal flat range, the DSAS plug-in in ArcGIS software was used to generate a cross-section every 100 m along the coastline of the study area. The intersection of the cross-section and the waterline was considered the tidal point, which represents the instantaneous tide level at the time of satellite overpass. The distance between each tide point and the baseline was calculated for each section. The nearest and farthest points from the baseline were selected to form a high tide point set and low tide point set, respectively, and all high tide points and low tide points were connected to determine the annual high tide line and low tide line. Finally, the upper and lower boundaries of the tidal flats were determined in this study (Figure 8). Fifteen Sentinel-2 images in the study area in 2020 were selected for the terline extraction experiments. First, the EWBI values were calculated to e terlines of 15 images through methods such as threshold segmentation, ed and visual interpretation, and the results are shown in Figure 8(a). Second, was used to generate a cross-section every 100 m for a total of 990 cross-s section generated 15 tide points along the waterline, for a total of 14,850 tid distance between each tide point and the baseline was calculated in each se farthest point from the baseline and the nearest point were filtered out, co high tide point set and low tide point set, as shown in Figure 8(b). Finally, al points and low tide points were connected to generate the annual high tide tide line, as shown in Figure 8(c). In this way, the area between the final high lines in Figure 8(c) was the tidal flat area of Jiaozhou Bay in 2020 obtained by method.

RF-W Model Extraction
The results of tidal flat extraction through the RF method and the wat showed that the extracted tidal flat range was roughly the same, but the differences in the details (Figure 9a). Compared with the waterline method, extracted by the RF method were more detailed, allowing the effective dist rivers, tidal flats, and creeks and the accurate extraction of tidal flats (Figu waterline method used multiple images to extract low tide point sets to d Fifteen Sentinel-2 images in the study area in 2020 were selected for the tidal flat waterline extraction experiments. First, the EWBI values were calculated to extract the waterlines of 15 images through methods such as threshold segmentation, edge detection, and visual interpretation, and the results are shown in Figure 8a. Second, the DSAS tool was used to generate a cross-section every 100 m for a total of 990 cross-sections. Each section generated 15 tide points along the waterline, for a total of 14,850 tide points. The distance between each tide point and the baseline was calculated in each section, and the farthest point from the baseline and the nearest point were filtered out, constituting the high tide point set and low tide point set, as shown in Figure 8b. Finally, all the high tide points and low tide points were connected to generate the annual high tide line and low tide line, as shown in Figure 8c. In this way, the area between the final high and low tide lines in Figure 8c was the tidal flat area of Jiaozhou Bay in 2020 obtained by the waterline method.

RF-W Model Extraction
The results of tidal flat extraction through the RF method and the waterline method showed that the extracted tidal flat range was roughly the same, but there were some differences in the details (Figure 9a). Compared with the waterline method, the tidal flats extracted by the RF method were more detailed, allowing the effective distinguishing of rivers, tidal flats, and creeks and the accurate extraction of tidal flats (Figure 9a(i)). The waterline method used multiple images to extract low tide point sets to determine the lower boundary of the tidal flats, which was a more accurate process than the RF method (Figure 9a(ii)). Therefore, we combined the two methods, adopted the shaping element tool in ArcGIS, and revised the extraction results of the RF method based on the lower boundary of the tidal flat extracted by the waterline method. First, select the shaping element tool in the editor toolbar to select a point in the extraction result of the RF method, then select the trace element tool to identify the lower boundary extracted by the waterline method, and finally, return to the starting point to form a closed curve to obtain the final tidal flats extraction result (Figure 9b).

Accuracy Evaluation
To verify the accuracy of the tidal flat information extracted in this resear sion matrix was used for accuracy evaluation [44]. The 30% sample dataset use racy verification described in section 3.1.4 was divided into two sets of tidal fla tidal flats, including 120 tidal flat verification points and 480 nontidal flat points. The confusion matrix of the three methods was calculated to obtain t traction accuracy, as shown in Tables 3, 4, and 5. The tables show that both accuracy and producer's accuracy of the RF-W model were the highest, which and 5.3% higher than those of the RF method alone, and the user's accurac ducer's accuracy of the waterline method were the lowest. These results indica RF-W model can effectively improve the accuracy of tidal flat information ext

Accuracy Evaluation
To verify the accuracy of the tidal flat information extracted in this research, a confusion matrix was used for accuracy evaluation [44]. The 30% sample dataset used for accuracy verification described in Section 3.1.4 was divided into two sets of tidal flats and nontidal flats, including 120 tidal flat verification points and 480 nontidal flat verification points. The confusion matrix of the three methods was calculated to obtain the final extraction accuracy, as shown in Tables 3-5. The tables show that both the user's accuracy and producer's accuracy of the RF-W model were the highest, which were 3.3% and 5.3% higher than those of the RF method alone, and the user's accuracy and producer's accuracy of the waterline method were the lowest. These results indicated that the RF-W model can effectively improve the accuracy of tidal flat information extraction.

Overall Changes in Tidal Flats
The RF-W model proposed in this study was used to extract the tidal flats of Jiaozhou Bay in 1990Bay in , 1995Bay in , 2000Bay in , 2005Bay in , 2010

Overall Changes in Tidal Flats
The RF-W model proposed in this study was used to extract the tidal flats of Jiaozhou Bay in 1990Bay in , 1995Bay in , 2000Bay in , 2005Bay in , 2010Bay in , 2015, and 2020 ( Figure 10). The tidal flats of Jiaozhou Bay are mainly distributed along the western, northern, and northeastern coasts of the gulf. Since 1990, the overall tidal flat area has shown a downward trend (Figure 10h

Tidal Flat Changes in Focus Areas
The area near Banqiaofang on the northeastern shore of Jiaozhou Bay, located to the east of the Moshui River, is a focus area for land reclamation in Qingdao (Figure 10a(i)). Beginning in 1995, the coastal area of Banqiaofang began to reclaim the sea area, resulting in a reduction in tidal flat area by approximately 2.9 km 2 (Figure 11a). In 1997, due to the construction of the Jiaozhou Bay Expressway, the shoreline became straighter and moved toward the sea by a large margin. Large areas of tidal flats disappeared and turned into roads and residential areas.

Tidal Flat Changes in Focus Areas
The area near Banqiaofang on the northeastern shore of Jiaozhou Bay, located to the east of the Moshui River, is a focus area for land reclamation in Qingdao (Figure 10a(i)). Beginning in 1995, the coastal area of Banqiaofang began to reclaim the sea area, resulting in a reduction in tidal flat area by approximately 2.9 km 2 (Figure 11a). In 1997, due to the construction of the Jiaozhou Bay Expressway, the shoreline became straighter and moved toward the sea by a large margin. Large areas of tidal flats disappeared and turned into roads and residential areas. The southwestern shore of Jiaozhou Bay, located in the area south of the D is adjacent to the western coast of Huangdao and is another focus area for lan tion in Qingdao (Figure 10a(ii)). In the past 30 years, Qingdao has made vigor to develop Huangdao District, and the original agricultural and fishery econom gradually replaced by sectors such as industry, commerce, and port-relate Large-scale land reclamation and the construction of factories have caused th area to decrease. In 2010, due to large-scale urbanization, the surface water p decreased, and the requirements for flood storage and discharge greatly inc meet the flood storage requirements, tidal flat reclamation was carried out on th bank of the Yanghe River, and the Yuejin River flood control dike was cons shown in Figure 11b, this dike was basically completed by the end of 2013, re tidal flat area by approximately 5.3 km 2 .
The southern bank of Jiaozhou Bay is located in the northeastern sea area munity of Xiaoshitou and is also a focus area for reclaiming and building po 10a(iii)). With the establishment of the Qingdao economic and technological de zone, Huangdao District's industry has developed rapidly, and the demand fo increased, resulting in frequent seaport building activities. As shown in Figu reclamation project started in the northeastern sea area of Xiaoshitou in 1992. B sea reclamation project was basically completed, and the port construction p mostly completed by 2006. The project construction led to a reduction in tidal approximately 1.9 km 2 . The southwestern shore of Jiaozhou Bay, located in the area south of the Dagu River, is adjacent to the western coast of Huangdao and is another focus area for land reclamation in Qingdao (Figure 10a(ii)). In the past 30 years, Qingdao has made vigorous efforts to develop Huangdao District, and the original agricultural and fishery economy has been gradually replaced by sectors such as industry, commerce, and port-related activity. Large-scale land reclamation and the construction of factories have caused the tidal flat area to decrease. In 2010, due to large-scale urbanization, the surface water permeability decreased, and the requirements for flood storage and discharge greatly increased. To meet the flood storage requirements, tidal flat reclamation was carried out on the northern bank of the Yanghe River, and the Yuejin River flood control dike was constructed. As shown in Figure 11b, this dike was basically completed by the end of 2013, reducing the tidal flat area by approximately 5.3 km 2 .

Analysis of the Causes of Tidal Flat Changes
The southern bank of Jiaozhou Bay is located in the northeastern sea area of the community of Xiaoshitou and is also a focus area for reclaiming and building ports (Figure 10a(iii)). With the establishment of the Qingdao economic and technological development zone, Huangdao District's industry has developed rapidly, and the demand for ports has increased, resulting in frequent seaport building activities. As shown in Figure 11c, the reclamation project started in the northeastern sea area of Xiaoshitou in 1992. By 1998, the sea reclamation project was basically completed, and the port construction project was mostly completed by 2006. The project construction led to a reduction in tidal flat area by approximately 1.9 km 2 .

Direct Cause
Coastal tidal flats are formed by rivers carrying sediment to the sea and nearby coast via an estuary. The tidal flat area depends on the relationship between sediment deposition and the hydrodynamic environment and is also affected by human factors. For Jiaozhou Bay, the interference of human activities dominates. Jiaozhou Bay is a semiclosed bay with superior natural conditions and rich fishery resources, providing people with an ideal living environment. Therefore, the tidal flats in this area are extremely susceptible to interference by human activities. The estuary of the Jiaozhou Bay River is mainly concentrated in the Dagu River estuary and Yanghe River Estuary in the northwest, and the Moshui River estuary and Baisha River estuary in the northeast (Figure 1). The tidal flat area in this region accounts for 80% of the total area of tidal flats in Jiaozhou Bay. After 2000, the regulation of the Dagu River, Yanghe River, Moshui River, and Baisha River began to increase. Reservoir dams were built upstream, and the interception of sediment reduced the sediment flux from the river to the sea, which affected the evolution of estuary beaches. The regulation of runoff by dams also affected the evolution of estuary beaches to a certain extent [45,46].
Because of the benefits of fish and salt in the bay and the convenience of boats providing good development conditions, sea reclamation and land reclamation are also important factors affecting the changes in the tidal flats of Jiaozhou Bay [47,48]. In the 1990s, with the significant improvement in China's economic development, the rate of reclamation in Jiaozhou Bay accelerated, and a large area of tidal flats was used for marine aquaculture. The area used for marine aquaculture increased rapidly from 11 km 2 in 1986 to 33 km 2 in 2000, accounting for 40% of the entire reclamation area [49]. Additionally, transportation also began to develop rapidly: The construction of the highway around Jiaozhou Bay made the eastern coastline straighter and moved the coastline toward the sea by a large margin, and a large number of beaches disappeared. With the vigorous development of Huangdao District after 2000, most reclamation activities in Jiaozhou Bay were concentrated in the Huangdao Qianwan Port and Free Trade Zone and Xuejiadao Coastal Shipyard. Ports and industries have become new factors affecting reclamation. With the promulgation and implementation of the Sea Area Management Measures in 2002, sea area use rights were determined by collecting sea area usage fees, reclamation activities were restricted by law, and the rate of sea reclamation in Jiaozhou Bay decreased. Therefore, the human activities along Jiaozhou Bay have had an important impact on the coastal changes.

Underlying Forces
Government control is also an important factor affecting the changes to the tidal flats in Jiaozhou Bay. Since Qingdao was listed as one of 14 open coastal cities in 1984, the government has encouraged economic development in this area, which has led to a rapid increase in population and has resulted in over-cultivation, weakened upstream soil fixation capacities, increased soil erosion, and increased river sediment loads. These factors caused gradual silting at the downstream entrance to the sea and increased the tidal flat area [50,51]. Additionally, large-scale reclamation activities transformed a large area of tidal flats into aquaculture ponds, which led to a decrease in the tidal flat area. In September 2014, the Qingdao Municipal Government formally implemented local regulations such as the "Regulations on the Protection of Jiaozhou Bay of Qingdao City" and implemented five major projects, including river improvement, dye source improvement, return to the sea, ecological restoration, and greenway construction around the bay. The construction of wetland parks and wetland protection areas led to a slight rebound in the tidal flat area of Jiaozhou Bay after 2015.
In addition, with the warming of the climate and the intensification of the greenhouse effect, the sea level has gradually risen, which has also become one of the factors affecting the changes in the tidal flats of Jiaozhou Bay. According to the data from the China Sea Level Bulletin of the State Oceanic Administration, the sea level in Jiaozhou Bay will continue to fluctuate and rise, and it is estimated that by 2050, the sea level of Shandong Province will have risen by 5-19 cm. Due to sea level rise, the beach surface of the silted bank will gradually become steeper, and the upward convex shape will become larger. The erosion of the beach surface will gradually slow, and the upward concave shape will tend to become straight. The change in the beach surface shape is bound to have an impact on the development of the beach, resulting in a decrease in the area of the beach [52]. Furthermore, the weakening of the hydrodynamic force, the reduction in tidal capacity, the reduction in water exchange and the reduction in seawater carrying capacity due to reclamation will cause an increase in the tidal flat area [53,54]. Although the influence of these factors is a long-term and slow process, their potential impacts on the evolution of tidal flats cannot be ignored.

The Damage Caused by Tidal Flat Changes
The ever-increasing reclamation activities have directly led to a reduction in tidal flat resources in Jiaozhou Bay, causing incalculable harm to the coastal environment and ecosystem. The tidal flats along Jiaozhou Bay are key habitats for fish spawning and migration, sources of food for migratory birds, and areas of animal and plant growth. Large-scale reclamation activities have converted a large number of tidal flat wetlands along Jiaozhou Bay, resulting in a large reduction in the biological species in the tidal flat wetlands and the endangerment or local disappearance of ecological service functions [55]. On the coast of Jiaozhou Bay, activities such as aquaculture, port docks, and port industries have increased the discharge of pollutants into the sea, while the reduction in tidal flats has weakened the self-purification capacity of the sea area. The pollution of the shore water environment and sedimentary environment has continued to intensify As Jiaozhou Bay is a semi-closed harbor, this kind of hazard will become more prominent.

Conclusions
With the excessive development of coastal areas, tidal flat resources have been continuously reduced, and there is an urgent need to accurately understand their distribution and the factors driving the change. This study proposes an improved model (RF-W model) for tidal flat extraction, which combines the RF method and the waterline method, and carries out model verification through the extraction of tidal flats in Jiaozhou Bay in 2020. The results show that the accuracy is distinctly better than those of the RF method and waterline method alone and demonstrate that the RF-W model can be applied for the extraction of large-scale tidal flat information and further spatiotemporal change analysis.
Information on tidal flats in the Jiaozhou Bay from 1990 to 2020 was extracted by the RF-W method, and the results showed a downward trend over the past 30 years. The total tidal flat area decreased by approximately 20 km 2 , of which the area decreased by nearly 17 km 2 from 1995 to 2010, accounting for 85% of the total decrease. The areas with large changes were mainly located on the southern bank of the Dagu River estuary and the areas near the two banks of the Moshui River estuary. The overall change mainly reflected tidal flat reclamation and expansion toward the sea and changes in the sediment content of rivers entering the sea. Human activities, such as reclamation, land development, and government policies, are the main factors affecting the tidal flat changes in Jiaozhou Bay. Over the long term, glacier melting and sea level rise due to climate warming and the greenhouse effect, as well as changes in hydrodynamics and tidal absorption due to reclamation, are also potential factors leading to the changes in the tidal flats in Jiaozhou Bay.