Spatiotemporal Changes of Coastline over the Yellow River Delta in the Previous 40 Years with Optical and SAR Remote Sensing

: The integration of multi-source, multi-temporal, multi-band optical, and radar remote sensing images to accurately detect, extract, and monitor the long-term dynamic change of coastline is critical for a better understanding of how the coastal environment responds to climate change and human activities. In this study, we present a combination method to produce the spatiotemporal changes of the coastline in the Yellow River Delta (YRD) in 1980–2020 with both optical and Synthetic Aperture Radar (SAR) satellite remote sensing images. According to the measurement results of GPS RTK, this method can obtain a high accuracy of shoreline extraction, with an observation error of 71.4% within one pixel of the image. Then, the inﬂuence of annual water discharge and sediment load on the changes of the coastline is investigated. The results show that there are two signiﬁcant accretion areas in the Qing 8 and Qingshuigou course. The relative high correlation illustrates that the sediment discharge has a great contribution to the change of estuary area. Human activities, climate change, and sea level rise that affect waves and storm surges are also important drivers of coastal morphology to be investigated in the future, in addition to the sediment transport.


Introduction
As one of the most active, complex, and fragile regions in the world, coastal, estuary, and river delta areas that are densely populated and economically developed contribute little to global landmass, where hydrosphere, lithosphere, atmosphere, biosphere, and human society interact very frequently [1]. Therefore, accurate and dynamic monitoring of the changing coastal and delta zone is critical for a better understanding of how the Earth's surface system responds to human activities.
As one of the most important topographic elements on land and sea maps, coastline is named as one of the 27 most important surface features by the International Geographic Data Committee (IGDC) [2]. Many large river deltas in the world are facing the impacts of shoreline change, and the study of delta erosion and expansion has attracted the attention of the scientific community and policymakers [3][4][5][6][7][8][9][10]. The change of coastline has an extremely important impact on the security of ports, the change of coastal ecological environment,

Study Area
Originating from the eastern Qinghai-Tibet Plateau, the Yellow River drains a wide basin of >75,000 km 2 with a total length of 5464 km [39]. The place where the Yellow River meets the Bohai Sea is the YRD located in Dongying City, Shandong Province, projecting in a fan shape, with a total area of 5943 km 2 ( Figure 1). There are ecological tourism areas, harbor industrial areas, high-end industrial areas, efficient ecological agriculture areas, and so on in this area.
YRD in the past 40 years and reveal its dynamic change drivers. Therefore, in this study, we firstly introduce a combination method for extracting the coastline of the YRD with both optical and SAR satellite remote sensing images and then evaluate the accuracy of results. Furthermore, we investigate the influence of annual water discharge and sediment load as well as critical sediment on the changes of the coastline.

Study Area
Originating from the eastern Qinghai-Tibet Plateau, the Yellow River drains a wide basin of >75,000 km 2 with a total length of 5464 km [39]. The place where the Yellow River meets the Bohai Sea is the YRD located in Dongying City, Shandong Province, projecting in a fan shape, with a total area of 5943 km 2 ( Figure 1). There are ecological tourism areas, harbor industrial areas, high-end industrial areas, efficient ecological agriculture areas, and so on in this area. Meanwhile, on both sides of the estuary of the new Yellow River and the old Yellow River, the YRD National Reserve of Shandong Province dominated by wetland types was established in 1992 and declared Ramsar wetland sites in 2013 [27]. The wetland is composed of two units, with the southern part located along the course of the Yellow River and extending out to the Bohai Sea, while the northern part is located at Diaokou River, referred to as the 'Ancient Yellow River'. This site is an almost naturally intact estuary wetland composed of marshes, reed swamps, tidal flats, canals and drainage channels, shallow estuarine waters, and aquaculture ponds at the mouth of the Yellow River Estuary (hereinafter referred to as YRE).
In modern times, there were nine major diversions of the river estuary, three of which occurred after the 1950s. Before 1953, the water flowed from the Tianshuigou Course (60%) and Shenxian'gou Course (40%) to the Bohai Sea. In July 1953, Tianshuigou river groove into the Shenxian'gou into the sea alone. On New Year's Day of 1964, the dike was destroyed and the Yellow River was diverted to the sea by the Diaokou River. After 1972, as silt accumulated, the flow began to branch again. In May 1976, the Yellow River shifted to Qingshuigou Course into the sea. In 1996, artificial diversion made the water flow from Qing 8 Course into the sea. This flow path continues to this day. In this study, the east side Meanwhile, on both sides of the estuary of the new Yellow River and the old Yellow River, the YRD National Reserve of Shandong Province dominated by wetland types was established in 1992 and declared Ramsar wetland sites in 2013 [27]. The wetland is composed of two units, with the southern part located along the course of the Yellow River and extending out to the Bohai Sea, while the northern part is located at Diaokou River, referred to as the 'Ancient Yellow River'. This site is an almost naturally intact estuary wetland composed of marshes, reed swamps, tidal flats, canals and drainage channels, shallow estuarine waters, and aquaculture ponds at the mouth of the Yellow River Estuary (hereinafter referred to as YRE).
In modern times, there were nine major diversions of the river estuary, three of which occurred after the 1950s. Before 1953, the water flowed from the Tianshuigou Course (60%) and Shenxian'gou Course (40%) to the Bohai Sea. In July 1953, Tianshuigou river groove into the Shenxian'gou into the sea alone. On New Year's Day of 1964, the dike was destroyed and the Yellow River was diverted to the sea by the Diaokou River. After 1972, as silt accumulated, the flow began to branch again. In May 1976, the Yellow River shifted to Qingshuigou Course into the sea. In 1996, artificial diversion made the water flow from Qing 8 Course into the sea. This flow path continues to this day. In this study, the east side of the Tiaohe course to the Qingshuigou course is selected as the study area.

Datasets
In this study, 10 key years are selected for coastline analysis on the premise that optical and SAR data sources are available, by taking into account the position change of the Yellow River estuary and the principle of uniform time sampling interval. Besides, the characteristics of water flow and sediment transport in the Yellow River basin vary significantly with the seasons [21]. The satellite remote sensing images as shown in Table 1, including Landsat-3 Multispectral Scanner (MSS), Landsat-5 Thematic Mapper (TM), Landsat-8 Operational Land Imager (OLI), Envisat Advanced Synthetic Aperture Radar (ASAR) Image Mode Precision (IMP) L1 Image, Sentinel-1 Interferometric Wide Swath (IWS) mode, and GaoFen-3 (GF-3) Fine Strip (FSII), are mainly focused on June, July, and August of each year. All images were selected with consideration of both spatial coverage and pixel resolution. The optical and SAR images are freely and readily available from United States Geological Survey (USGS), European Space Agency (ESA), and National Satellite Ocean Application Service (NSOAS). The optical false color composite images and SAR intensity images are shown in Figure 2. In this paper, SAR intensity information and band operation of optical images are used to carry out land and sea segmentation. To compare the differences in the extraction details between the two kinds of images, we investigated the shoreline extraction results of the YRD derived from the Landsat-8 OLI and Sentinel-1 IWS images by taking the example of the year 2017. The results are superposed on GaoFen-2 (GF-2) images with nadir pixel resolution of 3.2 m (see Figure 3). Due to the limitation of the GF-2 data source, the acquisition time was distributed between 2017 and 2020, leading to the actual distribution of coastlines in 2017 slightly different from those of the Landsat-8 and Sentinel-1.
In Figure 3, we can find that there is a significant difference between the shorelines in those nearshore areas where the slope is very flat and sediment is enriched (a, e, f, g) by visual interpretation. The result extracted from Sentinel-1 is closer inland than the Landsat-8. The biggest difference exists in the southern sediment-rich area of Qingshuigou Course in Figure 3f. In SAR images, the backscattering intensity is used to extract the coastline, which is the boundary between land and water. Therefore, we used existing Sentinel-1, GF-3, and Envisat ASAR images with a higher spatial resolution to replace Landsat images from 2004 onwards, except for that of 2013 (see Table 1).

Coastline Detection
As shown in Figure 4, the preprocessing of Landsat images was firstly performed with ENVI ® 5.5, including radiometric calibration and atmospheric correction. As a satellite-derived index derived from the near-infrared (NIR) and Green bands [40] (Equations (1)-(3)), normalized difference water index (NDWI) is usually used to make a preliminary distinction between land and seawater.

Coastline Detection
As shown in Figure 4, the preprocessing of Landsat images was firstly performed with ENVI ® 5.5, including radiometric calibration and atmospheric correction. As a satellitederived index derived from the near-infrared (NIR) and Green bands [40] (Equations (1)-(3)), normalized difference water index (NDWI) is usually used to make a preliminary distinction between land and seawater.  The NIR reflectance is affected by leaf internal structure and leaf dry matter content but not by water content [41]. Although the NDWI index was created for Landsat MSS images [30], it has been successfully used with other sensor systems in applications where the measurement of the extent of open water is needed [42]. In this study, we used Landsat-3 band 4/7, Landsat-5 band 2/5, and Landsat-8 band 3/5 to compute the NDWI values.
As for SAR images, GF-3 was preprocessed with the Pixel Information Expert (PIE ® ) SAR 6.0, while Envisat ASAR and Sentinel-1 were preprocessed with the Sentinel Application Platform (SNAP ® ) 7.0. After preprocessing, binarization processing was then employed. According to a series of morphological operations, the results showed a clear boundary of water and land. Enhanced Lee filter is used primarily to suppress speckle in radar imagery while simultaneously preserving texture information in heterogeneous areas. This filter that uses local statistics (coefficient of variation) within individual filter windows smooths images without removing edges or sharp features while minimizing the loss of radiometric and textural information. In areas containing isolated point targets, this filter preserves the observed value.
The most widely used threshold segmentation method is the fusion of NDWI and image binarization. Thresholding is a conceptually simple segmentation technique by using a histogram to segment images. At present, the optimal threshold is designed to determine and improve the threshold segmentation method. In this study, we used the threshold segmentation method proposed by Otsu [43], which is considered a gold standard in the field of thresholding. Figure 5 shows the flowchart of detailed water edge extraction. In the binary images obtained above, the gray level of pixels in the image was first converted to 0 or 1. The pixels with gray level 0 are represented by black, while the pixels with gray level 1 are represented by white. Morphology defines the white area with a gray level of 1 as the 'highlighted part' and the black part as the 'background'. Then, small water patches were The NIR reflectance is affected by leaf internal structure and leaf dry matter content but not by water content [41]. Although the NDWI index was created for Landsat MSS images [30], it has been successfully used with other sensor systems in applications where the measurement of the extent of open water is needed [42]. In this study, we used Landsat-3 band 4/7, Landsat-5 band 2/5, and Landsat-8 band 3/5 to compute the NDWI values.
As for SAR images, GF-3 was preprocessed with the Pixel Information Expert (PIE ® ) SAR 6.0, while Envisat ASAR and Sentinel-1 were preprocessed with the Sentinel Application Platform (SNAP ® ) 7.0. After preprocessing, binarization processing was then employed. According to a series of morphological operations, the results showed a clear boundary of water and land. Enhanced Lee filter is used primarily to suppress speckle in radar imagery while simultaneously preserving texture information in heterogeneous areas. This filter that uses local statistics (coefficient of variation) within individual filter windows smooths images without removing edges or sharp features while minimizing the loss of radiometric and textural information. In areas containing isolated point targets, this filter preserves the observed value.
The most widely used threshold segmentation method is the fusion of NDWI and image binarization. Thresholding is a conceptually simple segmentation technique by using a histogram to segment images. At present, the optimal threshold is designed to determine and improve the threshold segmentation method. In this study, we used the threshold segmentation method proposed by Otsu [43], which is considered a gold standard in the field of thresholding. Figure 5 shows the flowchart of detailed water edge extraction. In the binary images obtained above, the gray level of pixels in the image was first converted to 0 or 1. The pixels with gray level 0 are represented by black, while the pixels with gray level 1 are represented by white. Morphology defines the white area with a gray level of 1 as the 'highlighted part' and the black part as the 'background'. Then, small water patches were Remote Sens. 2021, 13,1940 8 of 27 removed after binarization. Furthermore, opening and closing operations were used to smooth the coastlines. pansion and corrosion are usually used in pairs to eliminate small particle noise and make the water edge smoother. Dilation operation process is to highlight the effect of target pixels extend outward into the background, and can splice and combine the fractured object for extracting the main body of the object. It also can be used to fill the patches in the region of interest, eliminate speckle noise in the region of interest. The function of erosion operation is to eliminate the highlighted boundary points. It is the process of boundary shrinking inward, aiming to eliminate the small and meaningless target objects at the boundary. Sobel operator method is simple, fast, and the extracted boundary is smooth and continuous. Therefore, Sobel edge detection operator was used to extract water edges in this paper. A good overview of the state of the art for coastline and shoreline derivation from remote sensing images is given by Gens [44], who notes that the term 'coastline' is customarily used in remote sensing-based studies, while the coastal community commonly refers to it as 'shoreline'. The definition of the shoreline theoretically is supposed to represent the linear boundary between the maritime. Though its apparent simplicity, for practical purposes, the dynamic nature of this boundary and its dependence on the temporal and spatial scale at which it is being considered results in the use of a range of shoreline indicators as following [30]. Therefore, the chosen definition should consider the shoreline in both a temporal and spatial sense and take account of the dependence of this variability on the time scale by which it is being investigated.
Accurate shoreline extraction is very important to find out the relationship between shoreline change and sediment load. Tidal variation plays a leading role during high resolution aerial survey, whilst it can be safely neglected for satellite data acquisition due to the coarser resolution of the satellite data. The aspect of time averaging can be critical for some remote sensing techniques [44]. Boak and Turner [14] proposed three groups of shoreline indicators as the 'true' shoreline positions as shown below.  A visually discernible indicator is a feature that can be physically seen, e.g., a previous high-tide line or the wet and dry boundary, usually preferred for photo interpretation, aerial photography, and earth observation data.  A specific tidal datum based shoreline indicator is determined by the intersection of the coastal profile with a specific vertical elevation, defined by the tidal constituents of a particular area, for example, mean high water (MHW) or mean sea level.  An indicator for the application of image processing techniques to extract proxy shoreline features from digital coastal images that are not necessarily visible to the Opening operation is a dilation after a erosion, using the same structuring element for both dilation and erosion. Closing operation is an erosion followed by a dilation. Expansion and corrosion are usually used in pairs to eliminate small particle noise and make the water edge smoother. Dilation operation process is to highlight the effect of target pixels extend outward into the background, and can splice and combine the fractured object for extracting the main body of the object. It also can be used to fill the patches in the region of interest, eliminate speckle noise in the region of interest. The function of erosion operation is to eliminate the highlighted boundary points. It is the process of boundary shrinking inward, aiming to eliminate the small and meaningless target objects at the boundary. Sobel operator method is simple, fast, and the extracted boundary is smooth and continuous. Therefore, Sobel edge detection operator was used to extract water edges in this paper.
A good overview of the state of the art for coastline and shoreline derivation from remote sensing images is given by Gens [44], who notes that the term 'coastline' is customarily used in remote sensing-based studies, while the coastal community commonly refers to it as 'shoreline'. The definition of the shoreline theoretically is supposed to represent the linear boundary between the maritime. Though its apparent simplicity, for practical purposes, the dynamic nature of this boundary and its dependence on the temporal and spatial scale at which it is being considered results in the use of a range of shoreline indicators as following [30]. Therefore, the chosen definition should consider the shoreline in both a temporal and spatial sense and take account of the dependence of this variability on the time scale by which it is being investigated.
Accurate shoreline extraction is very important to find out the relationship between shoreline change and sediment load. Tidal variation plays a leading role during high resolution aerial survey, whilst it can be safely neglected for satellite data acquisition due to the coarser resolution of the satellite data. The aspect of time averaging can be critical for some remote sensing techniques [44]. Boak and Turner [14] proposed three groups of shoreline indicators as the 'true' shoreline positions as shown below.

•
A visually discernible indicator is a feature that can be physically seen, e.g., a previous high-tide line or the wet and dry boundary, usually preferred for photo interpretation, aerial photography, and earth observation data.

•
A specific tidal datum based shoreline indicator is determined by the intersection of the coastal profile with a specific vertical elevation, defined by the tidal constituents of a particular area, for example, mean high water (MHW) or mean sea level.

•
An indicator for the application of image processing techniques to extract proxy shoreline features from digital coastal images that are not necessarily visible to the human eye and the accuracy of all techniques usually rely on the spatial resolution of the images used.
Due to the lack of tidal data for long-term span that could be used for tidal correction to get the mean high water, the second method would not be realizable. Therefore, Kuenzer et al. [27] used the first indicator to derive shorelines of YRD. However, the range of error caused by resolution that lies in the range of half a pixel is much less than the range of shoreline change.
In this study ( Figure 6), we mainly used the third indicator group to derive the waterland boundary as the shoreline. For the early low-resolution images in 1980 and 1985, the shoreline extracted were also inspected by combining the visual interpretation.
human eye and the accuracy of all techniques usually rely on the spatial resolution of the images used.
Due to the lack of tidal data for long-term span that could be used for tidal correction to get the mean high water, the second method would not be realizable. Therefore, Kuenzer et al. [27] used the first indicator to derive shorelines of YRD. However, the range of error caused by resolution that lies in the range of half a pixel is much less than the range of shoreline change.
In this study ( Figure 6), we mainly used the third indicator group to derive the waterland boundary as the shoreline. For the early low-resolution images in 1980 and 1985, the shoreline extracted were also inspected by combining the visual interpretation.
To evaluate the coastline extraction accuracy, we conducted Global Position System real time kinematic (GPS RTK) measurements of water boundary along the Gudong oilfield artificial coastal dam on 11 January 2020. The Gudong coastal dam, with a length of 12 km, was built in the 1980s to prevent the intrusion of seawater and provide the most important security barrier between the Gudong oilfield and Bohai Sea [15]. The highest point of the dam below the seawall is 4.3 m and the lowest point filled with dolosse to sea water is 3.4 m, both of which are higher than the high sea-level of 2.5 m.

Change Rate of Coastline
The rate of change in the past four decades is derived from a time series of shoreline vector data of the YRD with the Digital Shoreline Analysis System (DSAS) version 5.0. As is shown in Figures 7 and 8, end point rate (EPR) and linear regression rate (LRR) are used to describe the change rate of the coastline (Equations (4)-(7)). To evaluate the coastline extraction accuracy, we conducted Global Position System real time kinematic (GPS RTK) measurements of water boundary along the Gudong oilfield artificial coastal dam on 11 January 2020. The Gudong coastal dam, with a length of 12 km, was built in the 1980s to prevent the intrusion of seawater and provide the most important security barrier between the Gudong oilfield and Bohai Sea [15]. The highest point of the dam below the seawall is 4.3 m and the lowest point filled with dolosse to sea water is 3.4 m, both of which are higher than the high sea-level of 2.5 m.

Change Rate of Coastline
The rate of change in the past four decades is derived from a time series of shoreline vector data of the YRD with the Digital Shoreline Analysis System (DSAS) version 5.0. As is shown in Figures 7 and 8, end point rate (EPR) and linear regression rate (LRR) are used to describe the change rate of the coastline (Equations (4)- (7)).
where net shoreline movement (NSM) is equal to distance oldest and youngest shoreline; time span (SP) equals to time between oldest and most recent shoreline. where net shoreline movement (NSM) is equal to distance oldest and youngest shoreline; time span (SP) equals to time between oldest and most recent shoreline.  where net shoreline movement (NSM) is equal to distance oldest and youngest shoreline; time span (SP) equals to time between oldest and most recent shoreline.  A linear regression rate-of-change statistic can be determined by fitting a least-squares regression line to all shoreline points for a transect. As the slope of the line (Equation (6)), the calculation of the LRR is based on accepted statistical concepts [45,46]. The method of linear regression uses all the data, regardless of changes in trend or accuracy.
where x i is the year; y i is the distance from the intersection to the baseline along one transect; x and y are the averages of x i and y i , respectively. In the selection of the baseline position, the trend of all shorelines should be taken into account, especially in areas where the shoreline curvature changes greatly. In addition, the true change of shoreline of each section still needs to be manually identified. When selecting the location of the section, the selection of the tangent line depends on the setting of the baseline, especially the places with very drastic changes for some jumping points of the curve, which may also be difficult to reflect the real changes of the shoreline. In the section of uncertainty analysis of coastline extraction below, we set up a tangent with 50 m as the sampling interval to calculate the distance between the extracted coastline and the reference coastline along the tangent.
EPR and LRR can only roughly reflect the changes of shoreline. In cases where more data are available, EPR only uses the oldest and most recent shoreline, which would ignore additional information, e.g., the trends of accretion and erosion. The method of linear regression is susceptible to the effects of outliers and trends to underestimate the rate of change relative such as EPR [45][46][47]. In order to describe how shorelines change over time in detail, Figure 9 shows the staged 'end point rate' based on the distance of the intersection point along a transect between the adjacent shorelines, which presents the detailed information of shoreline accretion and erosion year by year in the YRD.
In shoreline analysis section, we use Pearson's correlation coefficient and coefficient of determination to analysis that shoreline change. As shown in Formula (8) where x and y are the average of two variables x, y, respectively. As shown in Equation (9), coefficient of determination R 2 is used to measure the percentage of variance in the data that is explained by a regression. As a dimensionless index, the value ranges from 1 to 0, where 1 is perfect fit. R 2 values close to 1.0 imply that the best-fit line explains most of the variation in the dependent variable.
where y, y , and y is measured data, predicted value, and mean of the measured data, respectively.
The observations of sediment load and water discharge are collected at the Yellow River Lijin Hydrometric Station (http://www.mwr.gov.cn accessed on 1 April 2021). The Lijin Station is the last hydrological station before the Yellow River enters the sea. The measurements of water discharge and sediment load from this station are used to investigate the correlation between the amount of water and sediment and the change of shoreline length and area of YRD.
Remote Sens. 2021, 13,1940 12 of 29 In shoreline analysis section, we use Pearson's correlation coefficient and coefficient of determination to analysis that shoreline change. As shown in Formula (8) where x and y are the average of two variables x, y, respectively.
As shown in Equation (9), coefficient of determination R 2 is used to measure the percentage of variance in the data that is explained by a regression. As a dimensionless index, the value ranges from 1 to 0, where 1 is perfect fit. R 2 values close to 1.0 imply that the best-fit line explains most of the variation in the dependent variable.
where y , y′, and y is measured data, predicted value, and mean of the measured data, respectively. The observations of sediment load and water discharge are collected at the Yellow River Lijin Hydrometric Station (http://www.mwr.gov.cn). The Lijin Station is the last hydrological station before the Yellow River enters the sea. The measurements of water discharge and sediment load from this station are used to investigate the correlation between the amount of water and sediment and the change of shoreline length and area of YRD.  Figure 7 for the analysis of shoreline change, the extracted shoreline ranges from the eastern part of Tiaohe Course to the southern part of Qingshuigou Course, where the distance between adjacent transects is set as 1000 m.

As shown in
The rate of change from 1980 to 2020 in Figure 8 shows spatial heterogeneous variability as a whole with average EPR of 87.87 m/a and LRR of 47.82 m/a, respectively. There are significant positive bursts in four areas including Tiaohe, Shenxian'gou, Qing 8, and Qingshuigou course, which are mainly due to artificial land reclamation and deposition of sediment. However, a great negative change with the erosion rate more than −150 m/a in EPR and LRR exists in the northernmost area of YRD near Diaokou river course that probably attributes to winter storm surge erosion.
In Shenxian'gou Course, the maximum change rates of the shoreline are 93.86 m/a in EPR and 52.18 m/a in LRR, which are 2.9 times the average EPR and 3.6 times of the average LRR within transects from no. 55 to no. 66 The staged average annual rate of change of adjacent shorelines in three key areas of the YRD is illustrated in Figure 10. From 1980 to 1985, all the three sections of shorelines were in a state of accretion. From 1985 to 1996, the annual change rate of shoreline in the Qingshuigou Course was as high as 400 m/a due to the high sediment transport, while there was a sharp decline in the period of 1996-1999 as a result of the artificial diversion of the Yellow River. From 1999 to the present, the overall annual change rate in the three areas is less than 100 m/a, which indicates that erosion and accretion alternate and gradually reach a balance, except that of Qing 8 in 2009Qing 8 in -2013 EPR and 52.18 m/a in LRR, which are 2.9 times the average EPR and 3.6 times of the average LRR within transects from no. 55 to no. 66. In Qing 8 Course, the maximum change rate reach 523.21 m/a in EPR and 385.4 m/a in LRR, which are 4.9 times the average EPR and 6.2 times the average LRR within transects from no. 82 to no. 101, respectively. In Qinghsuigou Course, the maximum change rates of the coastline reach 462.54 m/a in EPR and 311.48 m/a in LRR, which are 4.5 times the average EPR and 5.8 times of the average LRR within transects from no. 107 to no. 129. In summary, Qing 8 and Qingshuigou Course are the areas with the most dramatic changes of shoreline in the past 40 years.
The staged average annual rate of change of adjacent shorelines in three key areas of the YRD is illustrated in Figure 10. From 1980 to 1985, all the three sections of shorelines were in a state of accretion. From 1985 to 1996, the annual change rate of shoreline in the Qingshuigou Course was as high as 400 m/a due to the high sediment transport, while there was a sharp decline in the period of 1996-1999 as a result of the artificial diversion of the Yellow River. From 1999 to the present, the overall annual change rate in the three areas is less than 100 m/a, which indicates that erosion and accretion alternate and gradually reach a balance, except that of Qing 8 in 2009-2013.

Area Change, Water Discharge, and Sediment Load
As shown in Figure 11, the staged area changes over the YRD for the last four decades are estimated from the extracted coastline vector data. From 1980 to 1985, the entire YRD was basically in the state of accretion. In the years 1985-1996, the accretion area was mainly located in the Qingshuigou Course due to sufficient sediment accumulation, while the erosion area was located in the north part of the YRD. In 1996, the coastal boundary outline of Gudong Oilfield was initially formed, and the coastline of this area has been stable since 1996. From 1996 to 2003, the change of the mouth of the Yellow River led to the reduction of the area of Qingshuigou channel, and Qing 8 Course gradually expanded to the sea. From 2003 to 2013, the Qing 8 Course further expanded to the sea, mainly in the mouth of the Yellow River to the direction of the lobes on both sides of the river, while the Qingshuigou Course was eroded to a lesser extent. Since 2013, the YRD has been in a state of erosion on the whole, with only a slight expansion at the mouth of the river.
In Figure 12a, we can find that a similar variation trend between water discharge and sediment load, especially before the implementation of water and sediment regulation projects during the period of 1980-2002. According to the official Yellow River Sediment Bulletin (http://www.mwr.gov.cn/sj/, accessed on 13 May 2021), six major water and sediment regulation projects were carried out between 2002 and 2010. However, both water discharge and sediment transport began to decrease after 2002 [48], in which sediment transport was small in the years (2002,2003,2004,2005,2007, and 2010) with large water volume. Artificial hydraulic engineering increases the water discharge, but the total amount of sediment load does not increase significantly, which mainly resulted in the corresponding decrease of the shoreline length in the estuary area.
was basically in the state of accretion. In the years 1985-1996, the accretion area was mainly located in the Qingshuigou Course due to sufficient sediment accumulation, while the erosion area was located in the north part of the YRD. In 1996, the coastal boundary outline of Gudong Oilfield was initially formed, and the coastline of this area has been stable since 1996. From 1996 to 2003, the change of the mouth of the Yellow River led to the reduction of the area of Qingshuigou channel, and Qing 8 Course gradually expanded to the sea. From 2003 to 2013, the Qing 8 Course further expanded to the sea, mainly in the mouth of the Yellow River to the direction of the lobes on both sides of the river, while the Qingshuigou Course was eroded to a lesser extent. Since 2013, the YRD has been in a state of erosion on the whole, with only a slight expansion at the mouth of the river. In Figure 12a, we can find that a similar variation trend between water discharge and sediment load, especially before the implementation of water and sediment regulation projects during the period of 1980-2002. According to the official Yellow River Sediment Bulletin (http://www.mwr.gov.cn/sj/, accessed on 13 May 2021), six major water and sediment regulation projects were carried out between 2002 and 2010. However, both water discharge and sediment transport began to decrease after 2002 [48], in which sediment transport was small in the years (2002,2003,2004,2005,2007, and 2010) with large water volume. Artificial hydraulic engineering increases the water discharge, but the total amount of sediment load does not increase significantly, which mainly resulted in the corresponding decrease of the shoreline length in the estuary area.
In Figure 12b, the total length of coastline continued to increase (from 150 km to 220 km) and reached a peak during the period 1980-2005. From 2005 to 2013, the growth rate of shoreline slowed down and the length of the coastline decreased to 180 km. Since 2013, the coastline length has started to grow rapidly again, reaching 220 km. However, the total In Figure 12b, the total length of coastline continued to increase (from 150 km to 220 km) and reached a peak during the period 1980-2005. From 2005 to 2013, the growth rate of shoreline slowed down and the length of the coastline decreased to 180 km. Since 2013, the coastline length has started to grow rapidly again, reaching 220 km. However, the total area is accompanied by an increase and then a decrease with a period of about 5-10 years, and the highest peak of the total area is about 210-220 square kilometers.
Due to the twists and turns of the shoreline in the YRE, the change of shoreline length and area is not always positively correlated. As shown in Figure 12c, the coastline change is usually accompanied by a positive change in area, as in 1985, 1991, 1996, 1999, and 2009. However, the increase of length does not necessarily correspond to the increase of area, and vice versa, mainly due to the irregularity of the estuarine shoreline shape. Therefore, there are also shoreline changes accompanied by inverse changes in area, such as in 2003 and 2013.
In Figure 13, we can find that there is a strong positive correlation (up to 0.865) between water discharge and sediment load, which means sediment load is correspondingly high in the years with high water discharge. Furthermore, there is also a good positive correlation between the area change of the YRE and the entire area change of the YRD. Besides, both area changes are highly correlated with the water discharge and sediment load. and the highest peak of the total area is about 210-220 square kilometers.
Due to the twists and turns of the shoreline in the YRE, the change of shoreline length and area is not always positively correlated. As shown in Figure 12c, the coastline change is usually accompanied by a positive change in area, as in 1985, 1991, 1996, 1999, and 2009. However, the increase of length does not necessarily correspond to the increase of area, and vice versa, mainly due to the irregularity of the estuarine shoreline shape. Therefore, there are also shoreline changes accompanied by inverse changes in area, such as in 2003 and 2013. In Figure 13, we can find that there is a strong positive correlation (up to 0.865) between water discharge and sediment load, which means sediment load is correspondingly high in the years with high water discharge. Furthermore, there is also a good positive correlation between the area change of the YRE and the entire area change of the YRD. Besides, both area changes are highly correlated with the water discharge and sediment load.

Considerations of Sensors, Image Resolution, and Extration Methods
The sensors that collected images can be divided into two types: multi-spectral and SAR. The difference is mainly in the data preprocessing stage before the binarization step. NDWI that is calculated from two bands of the multispectral image can be considered as a nonlinear function related to reflectance, while filtered and radiometric calibrated amplitude that is generated from the SAR SLC image represents the backscattering intensity corresponding to the ground target closely related to the target medium, water content, and roughness. Then, both NDWI and amplitude images are used for binarization, morphologic operation, and edge detection, in which basically there is no loss of precision involved ( Figure 14).

Considerations of Sensors, Image Resolution, and Extration Methods
The sensors that collected images can be divided into two types: multi-spectral and SAR. The difference is mainly in the data preprocessing stage before the binarization step. NDWI that is calculated from two bands of the multispectral image can be considered as a nonlinear function related to reflectance, while filtered and radiometric calibrated amplitude that is generated from the SAR SLC image represents the backscattering intensity corresponding to the ground target closely related to the target medium, water content, and roughness. Then, both NDWI and amplitude images are used for binarization, morphologic operation, and edge detection, in which basically there is no loss of precision involved ( Figure 14). Image resolution can be considered the biggest source of uncertainty caused by the process of vector to raster conversion. The process of image acquisition by sensor is rasterization. Each pixel in raster data contains only one attribute data value, which is a generalization of various attributes within the pixel. The type of error that occurs when the feature is not captured properly is called a positional error, as opposed to attribute errors where information about the feature capture is inaccurate or false.
Geometric error refers to the positional error caused by the transformation of vector data into raster data, as well as the error of length, area, and topology matching caused by the positional error. The size of geometric error is proportional to the size of the pixel. When the polygon network represented by vector data is approximated by pixels, the problem of topology matching is serious.
Attribute error is related to the sensors that collect images and the image resolution. For example, the ground area corresponding to each pixel in the Landsat-3 MSS image is approximately 80 × 80 m. The attribute value of the pixel is the average reflected reflectivity of the objects. Some objects with high reflectivity will have a great influence on the pixel attribute value, even if their area is small. This will lead to misclassification and loss of some other useful information. Attribute error occurs when images are acquired. Due to the absence of measured data, especially as early as 1980, we cannot evaluate the attribute errors. However, when discussing uncertainty due to these various images used, attribute errors due to sensors and resolutions must be taken into account.
Positional error can be easily calculated with the measured data. According to overlying real vector and raster maps, the amount of correctly classified grids and misclassified grids can be counted. However, in many cases, the measured data cannot be provided. The rasterization error must then be estimated from the grid itself. The error is related to two factors. One is the distance between a random point in the vector map and pixel center in the raster map, and the other is the geometric characteristic of sampling grids. That means the higher the resolution, the smaller the pixel error.
It is not possible to use satellite images with the same resolution over a period of 40 years, especially from different sensors. However, the understanding of the relationship Image resolution can be considered the biggest source of uncertainty caused by the process of vector to raster conversion. The process of image acquisition by sensor is rasterization. Each pixel in raster data contains only one attribute data value, which is a generalization of various attributes within the pixel. The type of error that occurs when the feature is not captured properly is called a positional error, as opposed to attribute errors where information about the feature capture is inaccurate or false.
Geometric error refers to the positional error caused by the transformation of vector data into raster data, as well as the error of length, area, and topology matching caused by the positional error. The size of geometric error is proportional to the size of the pixel. When the polygon network represented by vector data is approximated by pixels, the problem of topology matching is serious.
Attribute error is related to the sensors that collect images and the image resolution. For example, the ground area corresponding to each pixel in the Landsat-3 MSS image is approximately 80 × 80 m. The attribute value of the pixel is the average reflected reflectivity of the objects. Some objects with high reflectivity will have a great influence on the pixel attribute value, even if their area is small. This will lead to misclassification and loss of some other useful information. Attribute error occurs when images are acquired. Due to the absence of measured data, especially as early as 1980, we cannot evaluate the attribute errors. However, when discussing uncertainty due to these various images used, attribute errors due to sensors and resolutions must be taken into account.
Positional error can be easily calculated with the measured data. According to overlying real vector and raster maps, the amount of correctly classified grids and misclassified grids can be counted. However, in many cases, the measured data cannot be provided. The rasterization error must then be estimated from the grid itself. The error is related to two factors. One is the distance between a random point in the vector map and pixel center in the raster map, and the other is the geometric characteristic of sampling grids. That means the higher the resolution, the smaller the pixel error.
It is not possible to use satellite images with the same resolution over a period of 40 years, especially from different sensors. However, the understanding of the relationship between image resolution and uncertainty of shoreline extraction is crucial for weighing the reliability of the results of shoreline changes in different periods.
The method we used to extract the shorelines is the most widely used threshold segmentation method from the fusion of NDWI and image binarization. Thresholding is a conceptually simple segmentation technique by using a histogram to segment images. At present, the optimal threshold is designed to determine and improve the threshold segmentation method. In this study, we used the threshold segmentation method proposed by Otsu [43], which is considered a gold standard in the field of thresholding.

Uncertainty Analysis of Coastline Extraction
With respect to river course or coastline extraction, several methods exist, such as image differencing, unsupervised land-water classification, supervised land-water classification, image segmentation, and manual coastline digitization. Most of them are based on multispectral bands or derived indices, such as the NDWI. A detailed review of the latest developments in coastline/shoreline derivation from remote sensing data has been presented by Gens [44]. However, it should be mentioned that while automatic digital image processing is fast and might seem more elegant, for coastline extraction, so far nothing can yield better results than a thorough human interpretation and manual digitization.
Since there are many types of coastline in the YRD region, it is difficult to use a single tidal station to correct the tide of the whole coastline. Furthermore, there is a lack of long time series of historical tidal observation data, especially in the YRE. Therefore, a dry-wet boundary is a more suitable choice for the coastline [27]. As far as the results of shoreline change are concerned, the annual average rate of change reached 87.87 m/a, and the range of change at the mouth of the Yellow River reached~500 m/a. The remote sensing images used in this paper were acquired from June to August and the tidal coastline migration made little contribution to the coastline change. Considering that the spatial resolution of the optical and SAR images is about 10-30 m except for earlier 80 m Landsat-3 MSS imagery, the accuracy of image extraction will have an impact on the result of coastline change.
By comparing the performance of SAR and optical images in shoreline extraction, we find that SAR images are more suitable for intertidal shoreline detection. Therefore, we used them instead of multispectral images at the corresponding time, when SAR images are available. Visual interpretation and manual identification are necessary in the procedure of automatic shoreline extraction, especially for the low-resolution multispectral images from 1980 and 1985. Due to the extremely frequent and drastic changes in the Yellow River Delta, there is still a lack of long-term continuous GPS time-series measurements. We only used GPS observation to evaluate the accuracy of the results in 2020.
As shown in Figure 15, the coastline extracted from GF-3 SAR on 11 July 2020, is consistent with observations from the GPS RTK to a large extent. Morphological operation was used to make the artificial coastal edge smooth, especially for those long, man-made piers that go deep into the sea, which would probably produce larger errors in the corner area. Figure 16 illustrates that the percentage of errors within one pixel of GF-3 SAR imagery (10 m) is 71.4%, while the percentage within two pixels (20 m) is 91.8%. The outliers more than 20 m come from the corner area as shown in Figure 15d. Therefore, through the evaluation from the GPS RTK measurements and manual visual interpretation, the authors think that the shoreline extraction accuracy is generally controlled at the level of 1-2 pixels.
For Qingshuigou and Qing 8 shorelines, we analyzed the confidence interval of the calculated change rate to investigate the uncertainty of the change rate. As illustrated in Figure 17, there are two transects selected, no. 97 for Qingshuigou course and no. 112 for Qingshuigou course. In Figure 18, 95% confidence belts were plotted using the distance from the intersections of the shorelines and the transect to the baseline. The standard error of the slope with confidence interval (LCI for ordinary linear regression) describes the uncertainty of the LRR.   Figure 15d. Therefore, through the evaluation from the GPS RTK measurements and manual visual interpretation, the authors think that the shoreline extraction accuracy is generally controlled at the level of 1-2 pixels. For Qingshuigou and Qing 8 shorelines, we analyzed the confidence interval of the calculated change rate to investigate the uncertainty of the change rate. As illustrated in Figure 17, there are two transects selected, no. 97 for Qingshuigou course and no. 112 for Qingshuigou course. In Figure 18, 95% confidence belts were plotted using the distance from the intersections of the shorelines and the transect to the baseline. The standard error of the slope with confidence interval (LCI for ordinary linear regression) describes the   In Qingshuigou course (Figure 18a), the band of confidence around the reported rate of change is 320 m/a ± 82 m/a. In other words, we can be 95% confident that the true rate of change is between 238 m/a and 402 m/a, leaving a 5% chance that the true line is outside those boundaries. Due to the artificial diversion of the Yellow River mouth from Qingshuigou to Qing 8 in 1996, two sections of 95% confidence band can be obtained by  In Qingshuigou course (Figure 18a), the band of confidence around the reported rate of change is 320 m/a ± 82 m/a. In other words, we can be 95% confident that the true rate of change is between 238 m/a and 402 m/a, leaving a 5% chance that the true line is outside those boundaries. Due to the artificial diversion of the Yellow River mouth from Qingshuigou to Qing 8 in 1996, two sections of 95% confidence band can be obtained by In Qingshuigou course (Figure 18a), the band of confidence around the reported rate of change is 320 m/a ± 82 m/a. In other words, we can be 95% confident that the true rate of change is between 238 m/a and 402 m/a, leaving a 5% chance that the true line is outside those boundaries. Due to the artificial diversion of the Yellow River mouth from Qingshuigou to Qing 8 in 1996, two sections of 95% confidence band can be obtained by linear fitting in Qing 8 (Figure 18b), that is, 1312 m/a ± 290 m/a (1022 m/a~1602 m/a) and −361 m/a ± 45 m/a (−406 m/a~−316 m/a).

Influence of Water Discharge and Sediment Load
Since the 1950s, human activities have had a more serious impact on the Yellow River basin due to the great economic and social development. Sediment losses due to urbanization, upstream irrigation, diversion, and hydraulic dams have altered the delta's erosion and accretion balance [49]. Delta ecosystems are often developed well beyond sustainable limits, especially in developing and emerging countries [22]. Human activities that impact the river system include dams and reservoirs, flow diversions, and water extractions. By far >45,000 dams over 15 m in height have been registered in the world with the purpose of flood control and water extraction [38].
Sediment transport plays an important role in the shoreline change of the YRD. The water-sediment regulation scheme, beginning in 2002, presented unexpected disturbances on the delta and coastal system. Increasing grain size of suspended sediment and decreasing suspended sediment concentration at the river mouth resulted in a regime shift of sediment transport patterns that enhanced the disequilibrium of the delta [21]. Over the past four decades, sediment load decreased stepwise revealed by the datasets from key gauging stations in the main stream, on account of both natural and anthropogenic impacts [39]. The sediment load of the Yellow River has an effect on many aspects of the YRD, including delta morphology and land subsidence [50,51].
From the overall trend of change, the sediment load of the Yellow River decreases gradually with the anthropological activities. The correlations between the change of shoreline length and water discharge, sediment load are very low, indicating that water discharge or sediment transport has little influence on the change of shoreline. The correlations between the area of the YRD and water discharge, sediment load reached 0.562 and 0.572, respectively, indicating that water discharge and sediment load are vital to area change. The impact of water transport and sediment transport on the YRD is more reflected in the change of area than the change of shoreline length.
Since the 1990s, the Yellow River flow cutoff has become increasingly serious, which brings many problems to the evolution of river course and estuaries [52]. The objective situation of the Estuary of the Yellow River also requires a new understanding of the estuary and the adoption of new countermeasures. At that time, it was one of the strategies of oil and gas development in the coastal area of Shengli Oil Field to exploit oil and gas turning from submerged coastal area to land, by using the silt of the Yellow River to make land. In 1996, owing to the Yellow River shifted to Qing 8 Course from Qingshuigou Couse, and the status of the YRE began to change obviously, with Qing 8 Lobe forming. From 1996 to 1999, there was obvious accretion in Qing 8 Course, with the change rate reaching 116.5 m/a, and obvious erosion occurred in Qingshuigou Course, with the erosion rate reaching 253.2 m/a. Along with the reduction of sediment discharge, there is a downward trend of the change range of the YRE coastline. It is a special time period from 2003 to 2009, Qing 8 experienced a state of obvious accretion, on account of water-sediment modulation. Water-sediment Modulation has played a vital role in regulating the delivery of material from the Yellow River to the sea [53,54]. From 2013 to 2020, the range of coastline change in the YRD further decreased, and there hardly occurs a distinct change in the three regions (Shenxian'gou Course, Qingshuigou, Qing 8), with the range of change within 50 m/a basically.
To further understand the relationship between sediment load and the area of the YRE (Figure 19), we analyzed the relationship between the area change of the YRE and water discharge, sediment load. In Figure 20, we can roughly see that the area at the YRE increases sharply in years with high sediment load (1980)(1981)(1982)(1983)(1984)(1985), and decreases slightly in years with low sediment load (1985)(1986)(1987)(1988)(1989)(1990)(1991). After the Yellow River shifted to Qing 8 (1996Qing 8 ( -2009, the area of erosion (mainly located in Qingshuigou Course) was larger than the area of accretion (mainly located in Qing 8). The correlations between estuary area change and water discharge, sediment load are shown in Figure 13, reaching 0.539 and 0.615, respectively, higher than the previous results of correlation analysis for the whole study area (0.562 and 0.572). The outcome proves that the influence of water discharge and sediment load on the area change of the YRD is mainly concentrated in the YRE.
Remote Sens. 2021, 13,1940 23 of 29 water discharge, sediment load are shown in Figure 13, reaching 0.539 and 0.615, respectively, higher than the previous results of correlation analysis for the whole study area (0.562 and 0.572). The outcome proves that the influence of water discharge and sediment load on the area change of the YRD is mainly concentrated in the YRE.   water discharge, sediment load are shown in Figure 13, reaching 0.539 and 0.615, respectively, higher than the previous results of correlation analysis for the whole study area (0.562 and 0.572). The outcome proves that the influence of water discharge and sediment load on the area change of the YRD is mainly concentrated in the YRE.   We can see that there is a strong correlation between sediment load and water discharge before 2003. From 2003 to 2009, water discharge increases slightly and remains in a relatively stable state, while sediment load decreases obviously due to the water and sediment regulation projects. However, there is still a weak correlation between them from the similar trends. The sediment load from upstream in the Yellow River mouth that does not keep a static state is an important factor affecting delta morphology. It can be found that the increase or decrease of sediment transport is accompanied by the increase or decrease of shoreline length, but there is a certain lag in the occurrence time of the latter. Due to the complicated ocean dynamics from waves and tides, accretion and erosion will gradually reach a balance.
Theoretically, the greater the water discharge, the greater the sediment transport. However, this is not always the case. Human activities have a great impact on the change of sediment in the upper reaches of the Yellow River, such as afforestation activities on the banks of the Yellow River, flood storage, and sediment control activities of hydropower stations, and gradual solidification and deposition of coarse sediment on the riverbed. Therefore, sediment transport does not necessarily increase at the same time as water flow increases.

Critical Sediment Deriven by Human Activities and Climate Change
The optical and SAR remote sensing images selected in this study could not fully reflect the detailed relationship between water, sediment transport, and estuary area with a time interval of 3-6 years. There are other underlying factors that are affecting the YRD. Previous studies usually assume that there is a linear correlation between sediment discharge and coastline change [26,51,52,55,56]. The analysis in Figure 13 shows that the sediment discharge has a great contribution to the change of estuary area with a correlation of 0.615. Therefore, the impulsive delivery of muds and sands indeed transforms the present YRD from a destructive phase to an accretion phase [57]. However, this study shows that the contribution of sediment transport to the change of the estuary is not the only dominant factor.
Due to the coarsening of the surface bed material sediments, the erosion efficiency of the downstream riverbed has reduced. Since 2013-2017 in Figures 9 and 10, continuous scouring and reduced sediment transport have caused the current delta to re-enter a phase of erosion and destruction [57]. The previous study has shown that erosion of the abandoned delta lobes can sustain the supply of sediment to remote offshore depocenters, even after decades of declining river sediment discharge [58]. The artificial water and sediment regulation scheme increases the ratio of water flow and sediment transport, resulting in more sediment being carried away by the current entering the Bohai Sea, and further reducing the impact of sediment transport on the estuary area [59].
Furthermore, human activities intensified gradually from a negligible factor from 7000 to 3000 cal yr BP to a considerable force compared with natural forces over the last 3000 years, and finally dominated the Yellow River's sediment discharge within the last few decades [60,61]. From the 1980s to the present, the anthropogenic influence on the YRE has been increasing, which further reduces the proportion of the impact of sediment transport on the estuary area. The transport of sediment into the sea significantly changes in the geomorphology of the estuary delta, and then the coupling of sediment transport and ocean dynamics has a superposition effect on the topography of the continental shelf [62,63]. As shown in Figure 21, the sediments from the Yellow River are mainly transported outward along two pathways, including the west side of the central Bohai Shoal to the central part of the Bohai Sea and the Bohai Strait to the Yellow Sea, which is closely related to sea surface temperature, flow field, wind speed, and wind direction [64]. Last but not the least, climate change and sea level rise affect waves and storm surges, which are also important drivers of coastal morphology [65,66]. Much of the world's sandy coastline is already eroding, and climate change is likely to exacerbate this as sea level rise is accelerating, in addition to the combined result of a wide range of potentially erosive or accretive factors [9]. In particular, the southern section of the Qingshuigou course has a lower slope and the effect of sea level rise is more serious [37]. In recent years, the impact of sediment discharge on the area change of the YRD is decreasing, so understanding the redistribution of sediments by waves and tides as well as dynamics of extreme weather patterns is critical for human-driven and climate change to deltas.

Conclusions
In this study, we present a combination method to extract the coastline of the Yellow River Delta (YRD) in 1980-2020 with both optical and SAR satellite remote sensing images. Then, we investigate the influence of annual water discharge and sediment load as well as critical sediment on the changes of the coastline. The long-term time series fusion of multi-source optical and SAR remote sensing data provides an accurate quantitative estimate of the coastline geolocation of the YRD in the last 40 years and reveals its dynamic change drivers.
The results indicate that the rate of change shows spatial heterogeneous variability during 40 years. There are four significant accretion areas (Tiaohe, Shenxian'gou, Qing 8, and Qingshuigou course) due to artificial land reclamation and deposition of sediment. In the Qing 8 Course, the maximum change rate reaches 523.21 m/a in end point rate (ERR) and 385.4 m/a in linear regression rate (LRR), while 462.54 m/a in ERR and 311.48 m/a in LRR in the Qingshuigou Course. However, there was a sharp decline in the period of 1996-1999 as a result of the artificial water and sediment regulation scheme. Erosion and accretion gradually reach a balance from 1999 to the present, except for Qing 8 in 2009-2013.
The relative high correlation supports that the sediment discharge has a great contribution to the change of estuary area. However, the influence of water discharge and sediment load on the area change of the YRD is mainly concentrated in the Yellow River Estuary (YRE). Although the artificial hydraulic engineering increases the water discharge, the total amount of sediment load does not increase significantly, resulting in the corresponding decrease of the shoreline length in the estuary area. Therefore, apart from Last but not the least, climate change and sea level rise affect waves and storm surges, which are also important drivers of coastal morphology [65,66]. Much of the world's sandy coastline is already eroding, and climate change is likely to exacerbate this as sea level rise is accelerating, in addition to the combined result of a wide range of potentially erosive or accretive factors [9]. In particular, the southern section of the Qingshuigou course has a lower slope and the effect of sea level rise is more serious [37]. In recent years, the impact of sediment discharge on the area change of the YRD is decreasing, so understanding the redistribution of sediments by waves and tides as well as dynamics of extreme weather patterns is critical for human-driven and climate change to deltas.

Conclusions
In this study, we present a combination method to extract the coastline of the Yellow River Delta (YRD) in 1980-2020 with both optical and SAR satellite remote sensing images. Then, we investigate the influence of annual water discharge and sediment load as well as critical sediment on the changes of the coastline. The long-term time series fusion of multi-source optical and SAR remote sensing data provides an accurate quantitative estimate of the coastline geolocation of the YRD in the last 40 years and reveals its dynamic change drivers.
The results indicate that the rate of change shows spatial heterogeneous variability during 40 years. There are four significant accretion areas (Tiaohe, Shenxian'gou, Qing 8, and Qingshuigou course) due to artificial land reclamation and deposition of sediment. In the Qing 8 Course, the maximum change rate reaches 523.21 m/a in end point rate (ERR) and 385.4 m/a in linear regression rate (LRR), while 462.54 m/a in ERR and 311.48 m/a in LRR in the Qingshuigou Course. However, there was a sharp decline in the period of 1996-1999 as a result of the artificial water and sediment regulation scheme. Erosion and accretion gradually reach a balance from 1999 to the present, except for Qing 8 in 2009-2013.
The relative high correlation supports that the sediment discharge has a great contribution to the change of estuary area. However, the influence of water discharge and sediment load on the area change of the YRD is mainly concentrated in the Yellow River Estuary (YRE). Although the artificial hydraulic engineering increases the water discharge, the total amount of sediment load does not increase significantly, resulting in the corresponding decrease of the shoreline length in the estuary area. Therefore, apart from the contribution of sediment transport to the change of the estuary, human activities, climate change, and