A Novel Method for River Bank Detection from Landsat Satellite Data: A Case Study in the Vietnamese Mekong Delta

: River bank (RB) erosion is a global issue a ﬀ ecting livelihoods and properties of millions of people. However, it has not received enough attention in the Vietnamese Mekong Delta (VMD), i.e., the world’s third largest delta, compared to salinity intrusion and ﬂooding. There have been several studies examining RB and coastal erosion in the VMD using remotely sensed satellite data, but the applied methodology was not adequately validated. Therefore, we developed a novel SRBED (Spectral RB Erosion Detection) method, in which the M-AMERL (Modiﬁed Automated Method for Extracting Rivers and Lakes) is proposed, and a new RB change detection algorithm using Landsat data. The results show that NDWI (Normalized Di ﬀ erence Water Index) and MNDWI (Modiﬁed Normalized Di ﬀ erence Water Index) using the M-AMERL algorithm (i.e., NDWI M-AMERL , MNDWI M-AMERL ) perform better than other indices. Furthermore, the NDWI M-AMERL; SMA (i.e., NDWI M-AMERL using the SMA (Spectral Mixture Analysis) algorithm) is the best RB extraction method in the VMD. The NDWI M-AMERL; SMA performs better than the MNDWI, NDVI (Normalized Di ﬀ erence Vegetation Index), and WNDWI (Weighted Normalized Di ﬀ erence Water Index) indices by 35–41%, 70% and 30%, respectively. Moreover, the NDVI index is not recommended for assessing RB changes in the delta. Applying the developed SRBED method and RB change detection algorithm, we estimated a net erosion area of the RB of –1.5 km 2 from 2008 to 2014 in the Tien River from Tan Chau to My Thuan, with a mean erosion width of –2.64 m and maximum erosion widths exceeding 60 m in places. Our advanced method can be applied in other river deltas having similar characteristics, and the results from our study are helpful in future studies in the VMD. than the MNDWI M-AMERL , revealing that this method is the best for RB extraction in the VMD.


Introduction
As the third-largest delta in the world [1], the Vietnamese Mekong Delta (VMD) is home to 17 million people and is the most important food basket of Vietnam and Southeast Asia. It contributes about 52% and more than 85% to the national annual food production and rice export, respectively [2]. In terms of biodiversity, the Mekong is only second to the Amazon [1]. The Mekong River is the world's 8 th largest river with critical natural resources. Specifically, the river supports the basis of the livelihood of local people, i.e., facilitating in-land water navigation while providing fresh-water fisheries and nutrient-rich sediment for agriculture [3]. Like many of the large river deltas worldwide, the VMD has been recently threatened by various environmental pressures, such as sea-level rise, land subsidence, floods, droughts, salinity intrusion, river bed incision, coastal erosion, and river bank (RB) erosion [1,2,[4][5][6][7][8]. The ever-increasing environment poses multiple critical impacts on the sustainability of the livelihoods of millions of VMD inhabitants [9][10][11][12]. Among these, RB erosion has not been given enough attention due to its complex phenomenon involving many physical-chemical processes and data shortage. However, RB erosion incidents have exacerbated continuously in the last decade [1,5,6,13] that critically threaten the residents' livelihoods. Therefore, the assessment of RB erosion with a reliable methodology is of utmost importance for the survivability and development of the VMD.
There are several techniques in studying the development of RB, including extraction and erosion/deposition, among which, bathymetric and topographic field measurements using total station surveying and terrestrial photogrammetry are the most widely used traditional methods [14,15]. Erosion pins and radio frequency identification equipped pebbles have also been widely used in the past decades [16,17]. Nowadays, UAV (Unmanned Aerial Vehicle) and airborne based techniques such as drone and LiDAR (Light Detection and Ranging) employing digital cameras and terrestrial laser scanning are increasingly used in river and coastal morphology [18][19][20]. These new methods are particularly useful in providing highly accurate topography and bathymetry for shallow water bodies. However, they are expensive, time-consuming, and labour-intensive in both pre-and post-processing. Other notable approaches in RB erosion studies are physical experiment and numerical modelling. While physical experiments have several advantages in studying RB erosion/deposition, especially for long-term assessments of dynamically stable channels [21], they are limited by the extensive labour involved. Scaling physical experiments appropriately is also challenging. Being comparatively cheaper, numerical models are increasingly popular [22,23]; nevertheless, they are not always applicable for long-term assessments and large-scale complex studies due to extensive data requirements (e.g., topography, bathymetry, geotechnical characteristics of river banks) and computational demand.
In recent years, many studies [1,5,6,[24][25][26][27][28][29][30][31] have employed remote sensing data for studying RB erosion. Satellite data with relatively reasonable temporal and spatial resolutions are especially suitable for long-term assessment of morphological changes in large-scale studies. Medium spatial resolution satellites are preferable because they can provide reasonable accuracy with high temporal resolution. Among them, Landsat is one of the most successful satellites in history because of its long-term catalogue of observational data of the Earth's surface [32].
The RB can be identified by classifying the satellite pixels into either water or land (i.e., hard classification). A widely used approach of this concept is based on the combination of two or more spectral bands to enhance the discrepancy between water bodies and land (i.e., spectral indices). Some widely used spectral indices based on this characteristic are NDVI (Normalized Difference Vegetation Index) [33], NDWI (Normalized Difference Water Index) [34], MNDWI (Modified Normalized Difference Water Index) [35], WNDWI (Weighted Normalized Difference Water Index) [36], and AWEI (Automated Water Extraction Index) [37]. Another approach to extract water from satellite data are mathematical morphological (MM) techniques based on topological and geometrical concepts. However, the hard classification concept may have considerable errors at mixed water-land pixels due to misclassification. To overcome this limitation, the soft classification concept has been developed, which estimates the percentage of water at mixed water-land pixels to achieve sub-pixel accuracy.
One popular approach in the soft classification concept is pixel unmixing and reconstruction which aims to exploit as much information as possible from a coarse pixel via a spectral mixture analysis (SMA), to reconstruct the fraction of each class at subpixel scale [32].
Although there have been some studies analyzing RB and coastal erosion processes in the VMD using satellite data [1,5,6], the methodology applied in previous studies was not adequately validated and the outcomes were not quantitatively assessed. Specifically, [6] tracked shoreline changes in the VMD using the NDVI which does not directly extract water from land but rather through above-ground vegetation detection [32]. Based on the literature in other river systems, [5] used the MNDWI to investigate the RB erosion in the VMD. However, we speculate that a single-index approach to identify RB is insufficient and unreliable given the complexities of flow dynamics, geometry, and land covers (especially in the case of VMD where aquaculture is situated very close to the river). Therefore, we develop a methodology for RB extraction and RB change detection in the VMD using Landsat data. Our methodology is a promising reference for RB erosion assessment in other large deltas in the globe. First, we assess the performance of five existing spectral indices (i.e., the NDWI, MNDWI, WNDWI, AWEI nsh (AWEI with shadow) and NDVI) in combination with three different segmentation methods (i.e., default thresholding, Otsu method, and MM techniques) for RB extraction. Second, the SMA technique is further applied to the best RB extraction method to enhance the accuracy of RB identification to sub-pixel which was not considered by any previous study in the VMD. We also develop a new algorithm for RB change detection. Our important contribution is an improved algorithm of MM and SMA techniques for RB segmentation methods in the VMD.

Study Area
The VMD, stretching from the Vietnam-Cambodia border to the East and West seas of Vietnam, is located in the lowermost reach of the Mekong River which runs through six countries from China to Vietnam ( Figure 1). The VMD is a flat delta with a mean ground elevation below 1.2 m (relative to the mean sea level at the Ha Tien datum), covering an area of 39,000 km 2 . The citizens living within the delta rely mainly on agro-aquaculture for their livelihoods. Tien and Hau (Vietnamese names of the Mekong and Bassac Rivers, respectively) are two main rivers discharging approximately 300-550 km 3 of water to the seas annually [38] via eight distributaries through a complex and spacious channel and floodplain system. The discharge of the Tien River is approximately fourfold that of the Hau River. The hydrograph in the VMD is distinguished by the flood (from June/July to November/December) and dry (from December/January to May/June of the following year) seasons. The dry season usually causes freshwater shortage and salinity intrusion [10][11][12]39], especially in severe drought years such as the one in 2015-2016 which was the most disastrous drought over the past 90 years [40]. The annual flooding normally causes large-scale inundation, leading to damage to agricultural land, infrastructures, and property. However, annual flood pulses are important for the development of the delta because it provides valuable sediment and fisheries, estimated as US$8-10 billion which is many times greater than the annual damages caused by floods, estimated as US$60-70 million [41]. The annually supplied sediment from the Mekong River is crucial for the sustainability of the VMD to prevent land subsidence due to groundwater extraction and sea-level rise [42,43]. zone in the Tien river is well above 5 m/year with the most severe cases reaching 30-50 m/years in certain places [5]. This study focuses on a 20 km long RB stretch around Cao Lanh in the Tien River ( Figure 1) to verify our developed RB detection model. This area is chosen because the RBs are highly active with the appearance of two dynamic islands and clear Google Earth (GE) images are available to assess the performance of different RB detection methods. We then apply the developed model to examine long-term RB erosion/deposition between 2008 and 2014 along 100 km in the Tien River from Tan Chau to My Thuan.

Landsat
The Landsat mission is a joint program between the USGS (U.S. Geological Survey) and NASA (National Aeronautics and Space Administration) which started in the 1970s. Since then, Landsat data have become increasingly applicable in monitoring applications in the Nile, Yellow, and Yangtze River Deltas e.g., [50][51][52], including ecological monitoring, change detection, biodiversity conservation, and flood and drought monitoring [24,[26][27][28]30], especially when these data became freely accessible in 2009. Raw Landsat images are prone to biases from the effects of the sensor, solar, atmosphere, and topographic features [31]. Therefore, preprocessing is important to minimize undesirable errors from these effects. To reduce our efforts, Landsat Level-2 products which are already processed to surface reflectance are used in this study. Pixels affected by instrument artefacts or subjected to cloud contamination are masked prior to the analysis using the provided quality assessment bands.
In developing a model for RB detection along 20 km around Cao Lanh area (Figure 1 Notably, in flat areas like the VMD having wide floodplains, water level fluctuations may influence RB detection significantly. Therefore, we intentionally selected the Landsat images for which the water level differences between the two selected datasets (i.e., for model development and In recent years, the VMD is increasingly prone to RB and coastal erosion [1,5,6,13], resulting from substantial reductions of the sediment supply from the Mekong River caused by upstream dams [4]. Specifically, the sediment flux has reduced from approximately 160 Mt/year in the normal conditions [44] to 40 Mt/year in 2012-2015 [45]. The situation is even made worse by the accelerated sand mining activities [46][47][48][49]. Sixty-four large hydropower dams were completed by 2017 [41], and more than 100 dams are additionally planned and under construction in the Mekong River basin. The total storage capacity of the existing 64 dams (over 80 km 3 ) accounts for >96% of the annual discharge at Chiang Saen, of which six mainstream dams in the upper Mekong River basin account for >50% [4].
In the VMD, it is speculated that the RB erosion is more fluvial than tidal, and the Tien River has its banks eroded substantially more than the Hau River [5]. The erosion rate in the river-dominated zone in the Tien river is well above 5 m/year with the most severe cases reaching 30-50 m/years in certain places [5]. This study focuses on a 20 km long RB stretch around Cao Lanh in the Tien River ( Figure 1) to verify our developed RB detection model. This area is chosen because the RBs are highly active with the appearance of two dynamic islands and clear Google Earth (GE) images are available to assess the performance of different RB detection methods. We then apply the developed model to examine long-term RB erosion/deposition between 2008 and 2014 along 100 km in the Tien River from Tan Chau to My Thuan.

Landsat
The Landsat mission is a joint program between the USGS (U.S. Geological Survey) and NASA (National Aeronautics and Space Administration) which started in the 1970s. Since then, Landsat data have become increasingly applicable in monitoring applications in the Nile, Yellow, and Yangtze River Deltas e.g., [50][51][52], including ecological monitoring, change detection, biodiversity conservation, and flood and drought monitoring [24,[26][27][28]30], especially when these data became freely accessible Remote Sens. 2020, 12, 3298 5 of 20 in 2009. Raw Landsat images are prone to biases from the effects of the sensor, solar, atmosphere, and topographic features [31]. Therefore, preprocessing is important to minimize undesirable errors from these effects. To reduce our efforts, Landsat Level-2 products which are already processed to surface reflectance are used in this study. Pixels affected by instrument artefacts or subjected to cloud contamination are masked prior to the analysis using the provided quality assessment bands.
In developing a model for RB detection along 20 km around Cao Lanh area (Figure 1), Landsat 8 OLI datasets on 22.02.2014 and 08.03.2019 are used. The developed model is then applied to examine RB erosion/deposition along 100 km from Tan Chau to My Thuan ( Figure 1) between 25.03.2008 and 21.01.2014. Notably, in flat areas like the VMD having wide floodplains, water level fluctuations may influence RB detection significantly. Therefore, we intentionally selected the Landsat images for which the water level differences between the two selected datasets (i.e., for model development and application) are neglected which guarantees reliable assessment of RB changes ( Table 1). The dry season is considered rather than the flood season because we attempt to avoid overestimation of the river-body due to flood water overtopping the RB. Moreover, the water levels at the time of retrieving the images are higher than the in-river floodplain levels to capture the real RB ( Figure 2). Images substantially covered by clouds are combined with other images in the same season to conserve the continuity of the RB for the analysis without notable error. application) are neglected which guarantees reliable assessment of RB changes ( Table 1). The dry season is considered rather than the flood season because we attempt to avoid overestimation of the river-body due to flood water overtopping the RB. Moreover, the water levels at the time of retrieving the images are higher than the in-river floodplain levels to capture the real RB ( Figure 2). Images substantially covered by clouds are combined with other images in the same season to conserve the continuity of the RB for the analysis without notable error.

Reference Dataset
Two GE datasets on 15.02.2014 and 27.01.2019 are used as reference datasets to verify the performance of the different RB detection methods in the present study. The RBs are digitized manually for both years (about 75 km RB per year). This is done with an "eye altitude" (zoom indication in GE) of ≤ 80 m which results in a map scale of roughly 1:400. Owing to the digitization precision of the RB is about 1 mm, a digitization error of ≤ 0.8 m can therefore be expected when comparing the two-years GE images. Where the RB is bordered by tall vegetation, it is furthermore possible that an inclination error is added to the reference dataset, which is related to the off-nadir of GE imagery. In the 2019 field survey, we observed that most of the trees along the main rivers in the for both years (about 75 km RB per year). This is done with an "eye altitude" (zoom indication in GE) of ≤ 80 m which results in a map scale of roughly 1:400. Owing to the digitization precision of the RB is about 1 mm, a digitization error of ≤ 0.8 m can therefore be expected when comparing the two-years GE images. Where the RB is bordered by tall vegetation, it is furthermore possible that an inclination error is added to the reference dataset, which is related to the off-nadir of GE imagery. In the 2019 field survey, we observed that most of the trees along the main rivers in the VMD were approximately 5 m, comparable to an average house. By measuring the shadow of the houses in the GE images in 2014 and 2019, the inclination angle is estimated at approximately 70 • relative to the RB. Therefore, the inclination error is estimated to be ≤ 2 m.
To increase the reliability of comparison between the two GE images, the two datasets are geometrically corrected to one another. More specifically, we used the Horn's quaternion-based ABSOR (Absolute Orientation) algorithm [53] to identify the rotation and translation that best maps one collection of point coordinates to the other. For its application, the RBs are divided into multiple 200m-segments with the marked ground control points (GCPs) ( Figure 3A). As the GE dataset is a mosaic of different aerial photographs, it is important to choose for each segment the nearest GCPs for its geometrical correction to minimize errors from mosaic processes. Well, identifiable fixed points in 2014 and 2019 are chosen as GCPs (e.g., road cross-sections, house corners) which are at a distance ≤550 m (D max ) to a segment, resulting in two or more GCPs per segment ( Figure 3B). We estimated that the maximum residual error (E max ) of GCPs between 2014 and 2019 after geometrical correction is 4.1 m, which is related to the mosaicing process while merging different aerial images into the final GE dataset. In essence, this error emerges from the GCPs lying on different aerial images, hence it cannot be corrected by the ABSOR. On the other hand, the error may come from the identification processes of the GCPs. Near the Mekong River, vegetation, fast-changing land occupation and the lack of a consistent infrastructure network with clear marks make the identification of suitable GCPs challenging. Finally, we estimated the total error (i.e., ≤0.8 m from digitization, ≤2 m from inclination, and ≤4.1 m from geometrical correction) from using GE imagery in RB change detection is <7 m.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 20 houses in the GE images in 2014 and 2019, the inclination angle is estimated at approximately 70° relative to the RB. Therefore, the inclination error is estimated to be ≤ 2 m.
To increase the reliability of comparison between the two GE images, the two datasets are geometrically corrected to one another. More specifically, we used the Horn's quaternion-based ABSOR (Absolute Orientation) algorithm [53] to identify the rotation and translation that best maps one collection of point coordinates to the other. For its application, the RBs are divided into multiple 200m-segments with the marked ground control points (GCPs) ( Figure 3A). As the GE dataset is a mosaic of different aerial photographs, it is important to choose for each segment the nearest GCPs for its geometrical correction to minimize errors from mosaic processes. Well, identifiable fixed points in 2014 and 2019 are chosen as GCPs (e.g., road cross-sections, house corners) which are at a distance ≤550 m (Dmax) to a segment, resulting in two or more GCPs per segment ( Figure 3B). We estimated that the maximum residual error (Emax) of GCPs between 2014 and 2019 after geometrical correction is 4.1 m, which is related to the mosaicing process while merging different aerial images into the final GE dataset. In essence, this error emerges from the GCPs lying on different aerial images, hence it cannot be corrected by the ABSOR. On the other hand, the error may come from the identification processes of the GCPs. Near the Mekong River, vegetation, fast-changing land occupation and the lack of a consistent infrastructure network with clear marks make the identification of suitable GCPs challenging. Finally, we estimated the total error (i.e., ≤0.8 m from digitization, ≤2 m from inclination, and ≤4.1 m from geometrical correction) from using GE imagery in RB change detection is <7 m.

RB Extraction
To extract the RB, we develop a new spectral RB erosion detection (SRBED) method, which can be divided into three main steps ( Figure 4): 1) calculation of the spectral indices using five different band combinations; 2) image segmentation into land and water using three different methodologies; 3) increase the accuracy of the best performing method by applying the SMA algorithm. Not all segmentation methods are compatible with all spectral indices, resulting in 11 different RB detection methods in total ( Table 2).

RB Extraction
To extract the RB, we develop a new spectral RB erosion detection (SRBED) method, which can be divided into three main steps ( Figure 4): 1) calculation of the spectral indices using five different band combinations; 2) image segmentation into land and water using three different methodologies; 3) increase the accuracy of the best performing method by applying the SMA algorithm. Not all segmentation methods are compatible with all spectral indices, resulting in 11 different RB detection methods in total (Table 2).   In this study, we employ five spectral indices: NDVI, NDWI, MNDWI, WNDWI, and AWEInsh, which are the most widely used for RB extraction. The algorithms and thresholds dividing into water and land of each spectral index are shown in Table 3. Table 3. Algorithm and thresholds dividing into water and land of the analyzed spectral indices. The equations display the band combinations (G = green; R = red; NIR = near-infrared; SWIR1,2 = short wave infrared).  In this study, we employ five spectral indices: NDVI, NDWI, MNDWI, WNDWI, and AWEI nsh, which are the most widely used for RB extraction. The algorithms and thresholds dividing into water and land of each spectral index are shown in Table 3. Table 3. Algorithm and thresholds dividing into water and land of the analyzed spectral indices. The equations display the band combinations (G = green; R = red; NIR = near-infrared; SWIR1,2 = short wave infrared).

Segmentation Methods
We apply three segmentation methods to distinguish water from land: default thresholding, Otsu method, and MM technique. Default thresholding refers to a fixed value per spectral index which groups the pixels into land and water. The thresholds for each index are indicated in Table 3. The Otsu method is one of the most successful methods for image thresholding [54]. It minimizes the intra-class variance, or equivalently, maximizes the inter-class variance of the input data values to select the optimal global threshold. The Otsu method uses either bi-level or multilevel thresholding to segment the images [55][56][57], of which the latter performs better than the former when the pictures have complicated features [54].
MM technique is only applied for NDWI and MNDWI as we found it unsuitable for the other spectral indices. It is inspired by the automated method for extracting rivers and lakes (AMERL) algorithm [25] which makes use of the image gradient by the magnitude (Gm) and the direction of the maximum increase of the intensity. Spectral index pixels with low Gm represent areas of relatively homogeneous land cover types, whereas pixels with high Gm may indicate a change between water and land. AMERL detects the RB using image segmentation by the watershed algorithm, which treats the image as a topographic map with the brightness of each pixel representing its height and finds the lines along the ridges. The RB is therefore extracted by considering not only the spectral characteristics of the pixels but also their topological connections.
Direct application of the watershed transform to a gradient image can result in serious over-segmentation due to local irregularities of the gradient [25]. To avoid this over-segmentation, AMERL contains an intermediate step where markers of identical low gradient values are identified in the image. The watershed algorithm identifies no ridge over these makers. [25] proposed thresholds for the indices NDWI, MNDWI and AWEI nsh to identify such markers (e.g., for NDWI: the threshold of land = -0.2; the threshold of the pure water = 0). Pixels with a spectral index value between these two thresholds are considered as mixed water-land pixels on the border of water surfaces. By marking all other pixels, the watershed algorithm finally identifies the RB over only the mixed water pixels.
However, we observed that these thresholds are not applicable in the VMD because, after applying the markers, the remaining unmarked pixels were too sparse to be connected to become a continuous line using segmentation method by the watershed algorithm. For this reason, we proposed a modified AMERL (M-AMERL), which contains the following processes. First, we calculate the Gm of the spectral indices using the Sobel operator, which is a discrete differentiation operator that is widely used to estimate a digital image gradient. The Sobel operator performs better than other gradient operators, such as the Prewitt operator [25]. Second, we identify markers to mask pixels which are not near the RB. For this, Gm is cleaned by assigning outliers (values that are more than three scales of the median deviation) and inland pixels to NaN (undefined values). Inland pixels are defined as pixels which do not have at least one pure water pixel (after the thresholds given by [25] for each index) in their close surroundings (2 pixels around analyzed pixel). Once this step is completed, the natural break over the cleaned Gm is calculated using the Otsu method. By doing this, it is assured that the obtained threshold is lower than Gm of RB pixels which have the highest Gm, and that the threshold is of water-land border pixels rather than due to some random inland pixels. The obtained threshold value is applied over the original Gm to segment the image into water-land border pixels and no-border pixels. No-border and inland pixels were then assigned to the value zero. Third, we segment the pixels into water and land by using the watershed algorithm, resulting in the identification of the RB. Finally, the identified RB pixels were sorted by solving the nearest neighbour algorithm developed by [58] for the MATLAB environment.

SMA Algorithm
The identified RB resulting from the M-AMERL segmentation algorithm is subject to a certain estimation error. Figure 5 shows that the RB is identified as running through the centre of the pixels where, for example, in pixel 1, land occupies most of the pixel. Therefore, the accuracy of the defined RB is low. To further improve the accuracy of the extracted RB position, the SMA is applied on the best RB extraction method identified in previous steps, based on a linear spectral mixture model (LSMM) [59] which is a pixel unmixing and reconstruction approach. It is based on the assumption that the spectral signature of a mixed pixel can be represented as the linear sum of the most representative spectral value of endmembers, weighted by their corresponding fractions, and has been applied to derive fraction images along the shoreline of large rivers [60]. This allows estimating the water fraction α contained in RB pixels and shifting the RB position accordingly.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 20 The obtained threshold value is applied over the original Gm to segment the image into water-land border pixels and no-border pixels. No-border and inland pixels were then assigned to the value zero. Third, we segment the pixels into water and land by using the watershed algorithm, resulting in the identification of the RB. Finally, the identified RB pixels were sorted by solving the nearest neighbour algorithm developed by [58] for the MATLAB environment.

SMA Algorithm
The identified RB resulting from the M-AMERL segmentation algorithm is subject to a certain estimation error. Figure 5 shows that the RB is identified as running through the centre of the pixels where, for example, in pixel 1, land occupies most of the pixel. Therefore, the accuracy of the defined RB is low. To further improve the accuracy of the extracted RB position, the SMA is applied on the best RB extraction method identified in previous steps, based on a linear spectral mixture model (LSMM) [59] which is a pixel unmixing and reconstruction approach. It is based on the assumption that the spectral signature of a mixed pixel can be represented as the linear sum of the most representative spectral value of endmembers, weighted by their corresponding fractions, and has been applied to derive fraction images along the shoreline of large rivers [60]. This allows estimating the water fraction α contained in RB pixels and shifting the RB position accordingly.

Figure 5. Schematic of RB identification resulted from the Modified Automated Method for Extracting
Rivers and Lakes (M-AMERL) segmentation method. Blue represents water; green represents land; the squares represent the satellite image pixels. The red crosses are the centres of the satellite pixels, and the dashed line is the identified RB using the M-AMERL segmentation method. We apply a dynamic nearest neighbour searching algorithm on a square of 5 × 5 pixels to estimate the water fraction in pixel 1.
In this study, we focus on land-water mapping; only water and land are therefore chosen as endmembers for the LSMM. The reflectance value of a mixed pixel (Rmix) is expressed as where Rw and RL are the representative reflectance values of water and land, respectively. The choice of these representative reflectance values per endmembers is crucial in this analysis. We apply a dynamic nearest neighbour searching algorithm on a square of 5 × 5 pixels ( Figure 5) to estimate the water fraction of each pixel. Within this square, pure land and pure water pixels are identified by applying the thresholds proposed by [25] per spectral index. Pixels obstructed by clouds or Figure 5. Schematic of RB identification resulted from the Modified Automated Method for Extracting Rivers and Lakes (M-AMERL) segmentation method. Blue represents water; green represents land; the squares represent the satellite image pixels. The red crosses are the centres of the satellite pixels, and the dashed line is the identified RB using the M-AMERL segmentation method. We apply a dynamic nearest neighbour searching algorithm on a square of 5 × 5 pixels to estimate the water fraction in pixel 1.
In this study, we focus on land-water mapping; only water and land are therefore chosen as endmembers for the LSMM. The reflectance value of a mixed pixel (R mix ) is expressed as where R w and R L are the representative reflectance values of water and land, respectively. The choice of these representative reflectance values per endmembers is crucial in this analysis. We apply a dynamic nearest neighbour searching algorithm on a square of 5 × 5 pixels ( Figure 5) to estimate the water fraction of each pixel. Within this square, pure land and pure water pixels are identified by applying the thresholds proposed by [25] per spectral index. Pixels obstructed by clouds or containing erroneous pixel values are excluded from the analysis. Once the pure water and pure land pixels within the 5 × 5 pixel square are identified, the mean of their spectral index values (i.e., representative reflectance values) is calculated. Finally, the water fraction α (in the range from 0 to 1) is estimated and the ultimate position of the RB in each pixel is shifted accordingly. The principle of this approach is the assumption that the subpixel water/land portion in the RB pixel is composed of a similar water/land cover type in its surrounding pixels. This approach is similar to the one described by [61] which has shown promising results. If α < 0 or α > 1, R mix is smaller or bigger than R L or R w respectively. This means that the analyzed pixel contains only water or only land. In this case, the RB pixel is moved one pixel in the corresponding direction (e.g., α > 1; R mix > R w ; RB pixel is then moved one-pixel inland direction). The SMA is then applied over this new pixel to define the exact RB position, following the aforementioned procedure.

RB Change Detection Algorithm
Once the RBs are reliably extracted using the developed method, temporal RB changes (i.e., erosion or deposition) are estimated using the newly developed algorithm as illustrated in Figure 6 where we describe a step-by-step procedure of change detection between 2014 and 2019. Starting from the fixed main control points (MCPs) such as the river mouth of side channels, we divided the 2014 RB into multiple segments with a length of L seg = 200m ( Figure 6A). In every segment, two land points (LPs) are identified, which serves as the baseline for the estimation of the surface change between the two years ( Figure 6B). Then we estimate the area from the two LPs to the 2014 RB, i.e., the orange shape ( Figure 6C). Similarly, we estimate the area from the two LPs to the 2019 RB, i.e., the blue shape ( Figure 6D-E). Finally, we estimate the surface change between the 2014 and 2019 RBs as the difference between the two estimated areas, i.e., the red shape ( Figure 6F). If the area of the orange shape is greater than that of the blue shape, the RB is eroded; otherwise, the RB is deposited. This newly developed RB change detection algorithm is applied for the GE data (reference data set) and all 11 methods using Landsat data in Table 2. Finally, the best RB extraction method is achieved by comparing the 2014-2019 RB changes among 11 methods using the Landsat and the GE data. Notably, segments obstructed by clouds, cloud shadows or containing erroneous pixel values are excluded from the analysis. If the area of the orange shape is greater than that of the blue shape, the RB is eroded; otherwise, the RB is deposited. Green represents land; light blue represents water. The background image is from 2014.

Performance of the Developed M-AMERL Method in RB Extraction
In selecting the best RB extraction method out of 11 possibilities including M-AMERL (Table 2), the surface changes between the 2014 and 2019 RBs using the Landsat data are compared with those using the GE data by considering the errors of the variance and the mean. Figure 7 shows that the vegetation index NDVI has a much lower accuracy than the water indices, suggesting that water indices should be used in quantifying RB erosion/deposition in the VMD. Using the default thresholding and Otsu segmentation methods, the WNDWI and AWEInsh perform better than the MNDWI and NDWI. It indicates that the WNDWI and AWEInsh are relatively applicable in the VMD. In fact, the WNDWI has been specifically developed for turbid water environments [36] and the AWEInsh was proposed to increase water extraction accuracy in the presence of various sources of If the area of the orange shape is greater than that of the blue shape, the RB is eroded; otherwise, the RB is deposited. Green represents land; light blue represents water. The background image is from 2014.

Performance of the Developed M-AMERL Method in RB Extraction
In selecting the best RB extraction method out of 11 possibilities including M-AMERL (Table 2), the surface changes between the 2014 and 2019 RBs using the Landsat data are compared with those using the GE data by considering the errors of the variance and the mean. Figure 7 shows that the vegetation index NDVI has a much lower accuracy than the water indices, suggesting that water indices should be used in quantifying RB erosion/deposition in the VMD. Using the default thresholding and Otsu segmentation methods, the WNDWI and AWEI nsh perform better than the MNDWI and NDWI. It indicates that the WNDWI and AWEI nsh are relatively applicable in the VMD. In fact, the WNDWI has been specifically developed for turbid water environments [36] and the AWEI nsh was proposed to increase water extraction accuracy in the presence of various sources of environmental noise [37], which are the case in the VMD. For instance, [4] found from field surveys that the suspended sediment concentration in the VMD may reach up to 1132 g/m 3 . Normally, the Otsu segmentation method performs better than the default thresholding in terms of mean estimation error.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 20 environmental noise [37], which are the case in the VMD. For instance, [4] found from field surveys that the suspended sediment concentration in the VMD may reach up to 1132 g/m 3 . Normally, the Otsu segmentation method performs better than the default thresholding in terms of mean estimation error. Figure 7. Performance of 11 RB extraction method compared to the reference Google Earth (GE) data sets in terms of the variance of and mean errors. The NDWI using our developed M-AMERL segmentation algorithm results in the best method of RB extraction in the VMD which is recommended to be used in quantifying RB erosion/deposition.
Notably, the MNDWI and NDWI applying the default and Otsu segmentation methods are not as good as the WNDWI and AWEInsh (Figure 7). In other words, these methods should not be used in RB extraction in the VMD. However, the MNDWI and NDWI indices using the M-AMERL segmentation method, i.e., MNDWIM-AMERL and NDWIM-AMERL respectively, are the best among 11 methods, verifying the potential applicability of M-AMERL. Figure 8A,B shows that the M-AMERL outperforms the Otsu segmentation method, especially in areas having aquaculture ponds close to the RBs (Figure 8 C,D) which are common in the VMD. The Otsu method could not distinguish the RBs from the banks of these ponds properly. Finally, we found that the NDWIM-AMERL is slightly better than the MNDWIM-AMERL, revealing that this method is the best for RB extraction in the VMD. Notably, the MNDWI and NDWI applying the default and Otsu segmentation methods are not as good as the WNDWI and AWEI nsh (Figure 7). In other words, these methods should not be used in RB extraction in the VMD. However, the MNDWI and NDWI indices using the M-AMERL segmentation method, i.e., MNDWI M-AMERL and NDWI M-AMERL respectively, are the best among 11 methods, verifying the potential applicability of M-AMERL. Figure 8A,B shows that the M-AMERL outperforms the Otsu segmentation method, especially in areas having aquaculture ponds close to the RBs ( Figure 8C,D) which are common in the VMD. The Otsu method could not distinguish the RBs from the banks of these ponds properly. Finally, we found that the NDWI M-AMERL is slightly better than the MNDWI M-AMERL , revealing that this method is the best for RB extraction in the VMD.

Advantage of the SMA Algorithm
Studies on RB erosion/deposition in the VMD using satellite data that does not refine the images into sub-pixels may be subject to mix-identify the exact position of the RB in each pixel ( Figure 5). As such, the RB becomes a zigzag line which is not natural and realistic. This is because the algorithm defines the RB at the centre of each pixel ( Figure 5). To overcome this limitation, we apply the SMA on the best RB extraction method, i.e., NDWIM-AMERL, to enhance the accuracy of RB extraction. Figure  9 shows that the RB becomes smoother when applying the SMA algorithm compared to the pure NDWIM-AMERL because the SMA algorithm shifts the RB in each pixel close to the real position. Compared to the pure NDWIM-AMERL, the SMA algorithm (i.e., NDWIM-AMERL; SMA) reduces the mean estimation error and the variance by approximately 20% (Figure 10). Notably, the relative improvement in RB extraction between the NDWIM-AMERL; SMA and the NDWI using the default threshold and Otsu methods considering the errors of the mean and the variance is about 40-50% and 80-90%, respectively.
Moreover, the NDWIM-AMERL; SMA performs better than the MNDWI using the default threshold and Otsu methods by 41% and 35%, respectively, in terms of the mean estimation error (Figure 10). The numbers of the error of the variance are approximately 87% and 78%, respectively. The accuracy of the MNDWIM-AMERL and NDWIM-AMERL is relatively similar, which is lower than the NDWIM-AMERL; SMA by nearly 20%. Similarly, the NDWIM-AMERL; SMA performs better than the WNDWI and AWEInsh indices by 30% and 60% considering the estimation error of the mean and the variance, respectively. Importantly, the NDVI is the poorest index in RB detection in the VMD: the NDWIM-AMERL; SMA is approximately 70% and >90% better than the NDVI index considering the estimation error of the mean and the variance, respectively.

Advantage of the SMA Algorithm
Studies on RB erosion/deposition in the VMD using satellite data that does not refine the images into sub-pixels may be subject to mix-identify the exact position of the RB in each pixel ( Figure 5). As such, the RB becomes a zigzag line which is not natural and realistic. This is because the algorithm defines the RB at the centre of each pixel ( Figure 5). To overcome this limitation, we apply the SMA on the best RB extraction method, i.e., NDWI M-AMERL , to enhance the accuracy of RB extraction. Figure 9 shows that the RB becomes smoother when applying the SMA algorithm compared to the pure NDWI M-AMERL because the SMA algorithm shifts the RB in each pixel close to the real position. Compared to the pure NDWI M-AMERL , the SMA algorithm (i.e., NDWI M-AMERL; SMA ) reduces the mean estimation error and the variance by approximately 20% (Figure 10). Notably, the relative improvement in RB extraction between the NDWI M-AMERL; SMA and the NDWI using the default threshold and Otsu methods considering the errors of the mean and the variance is about 40-50% and 80-90%, respectively.
Moreover, the NDWI M-AMERL; SMA performs better than the MNDWI using the default threshold and Otsu methods by 41% and 35%, respectively, in terms of the mean estimation error ( Figure 10). The numbers of the error of the variance are approximately 87% and 78%, respectively. The accuracy of the MNDWI M-AMERL and NDWI M-AMERL is relatively similar, which is lower than the NDWI M-AMERL; SMA by nearly 20%. Similarly, the NDWI M-AMERL; SMA performs better than the WNDWI and AWE Insh indices by 30% and 60% considering the estimation error of the mean and the variance, respectively. Importantly, the NDVI is the poorest index in RB detection in the VMD: the NDWI M-AMERL; SMA is approximately 70% and >90% better than the NDVI index considering the estimation error of the mean and the variance, respectively.

Applying the Developed SRBED Methodology on Quantifying RB Changes in the River-Dominated Zone in the VMD
To demonstrate the applicability of the developed SRBED methodology (e.g., the newly developed NDWIM-AMERL; SMA and change detection algorithm), we perform an analysis for quantifying the RB changes over a 100 km river segment in the Tien River, spanning from Tan Chau to My Thuan ( Figure 1B). This river segment is river-dominated influenced by tidal water level variations in the dry season. Over this area, the RBs are extracted and compared between 25.03.2008 and 21.01.2014. These dates are used because they have comparable water levels ( Table 1) to eliminate the influence

Applying the Developed SRBED Methodology on Quantifying RB Changes in the River-Dominated Zone in the VMD
To demonstrate the applicability of the developed SRBED methodology (e.g., the newly developed NDWIM-AMERL; SMA and change detection algorithm), we perform an analysis for quantifying the RB changes over a 100 km river segment in the Tien River, spanning from Tan Chau to My Thuan ( Figure 1B). This river segment is river-dominated influenced by tidal water level variations in the dry season. Over this area, the RBs are extracted and compared between 25.03.2008 and 21.01.2014. These dates are used because they have comparable water levels ( Table 1) to eliminate the influence

Applying the Developed SRBED Methodology on Quantifying RB Changes in the River-Dominated Zone in the VMD
To demonstrate the applicability of the developed SRBED methodology (e.g., the newly developed NDWI M-AMERL; SMA and change detection algorithm), we perform an analysis for quantifying the RB changes over a 100 km river segment in the Tien River, spanning from Tan Chau to My Thuan ( Figure 1B). This river segment is river-dominated influenced by tidal water level variations in the dry season. Over this area, the RBs are extracted and compared between 25.03.2008 and 21.01.2014. These dates are used because they have comparable water levels (Table 1) to eliminate the influence of the water level on the position of the real banks. Figure 11 shows that the RBs are on average eroded with the erosion area and erosion width having the same pattern. Moreover, we also checked the quantification manually and found no error. These indicate that our developed algorithms work properly. eroded with the erosion area and erosion width having the same pattern. Moreover, we also checked the quantification manually and found no error. These indicate that our developed algorithms work properly.
Over a short period from 2008 to 2014, the RBs from Tan Chau and My Thuan are on average eroded at considerable rates ( Figure 11). The total net RB erosion area (i.e., erosion subtracting deposition) and mean erosion width are -1.5 km 2 and -2.64 m for this period, which are equivalent to mean annual rates of -0.25 km 2 /year and -0.44 m/year, respectively. Notably, we identify some segments having high RB erosion such as the right and left banks around Tan Chau (from 6 km upstream to 4 km downstream), the head of the islands, around Cao Lanh, and from My Thuan to Sa Dec (about 15 km upstream of My Thuan). The most severe segment is around Sa Dec on the right bank of the Tien River where the river is meandering.

Discussion
Several studies have been conducted to investigate long-term morphological changes in the VMD; however, the RB and shoreline changes were quantified using spectral indices based on either the authors' experience or the literature in other river systems without adequate justification in the context of the VMD. [6] applied the NDVI index to investigate the shoreline evolution of the VMD from 1973 to 2015 based on the Landsat data. Recently, [5] quantified long-term RB and coastal erosion in the entire VMD using Landsat data from 1989 to 2014. They used the MNDWI index because it has been widely used in similar studies in other river systems. We found that SRBED combining NDWI and MNDWI with the M-AMERL segmentation algorithm performs better than other methods among 11 combinations (Table 2; Figures 7 and 8). This is because the M-AMERL algorithm depends not only on the spectral characteristics of the pixels but also on their topological connections. This characteristic allows detecting the RB at relatively correct locations even in noisy environments (e.g., near aquaculture areas which are common in the VMD) where other approaches are not suitable. Incorporating the SMA further improves the performance of RB extraction ( Figure  9). We conclude that the NDWIM-AMERL; SMA is the best RB extraction method in the VMD using the Landsat data. Considering the mean estimation error, the NDWIM-AMERL; SMA is 35-41%, 70% and 30% better than the MNDWI, NDVI and WNDWI indices, respectively ( Figure 10). Notably, the NDVI index has the poorest performance in RB extraction in the VMD. Moreover, the MNDWI and NDWI using the default thresholding and Otsu methods do not yield highly accurate RB positions. Over a short period from 2008 to 2014, the RBs from Tan Chau and My Thuan are on average eroded at considerable rates ( Figure 11). The total net RB erosion area (i.e., erosion subtracting deposition) and mean erosion width are -1.5 km 2 and -2.64 m for this period, which are equivalent to mean annual rates of -0.25 km 2 /year and -0.44 m/year, respectively. Notably, we identify some segments having high RB erosion such as the right and left banks around Tan Chau (from 6 km upstream to 4 km downstream), the head of the islands, around Cao Lanh, and from My Thuan to Sa Dec (about 15 km upstream of My Thuan). The most severe segment is around Sa Dec on the right bank of the Tien River where the river is meandering.

Discussion
Several studies have been conducted to investigate long-term morphological changes in the VMD; however, the RB and shoreline changes were quantified using spectral indices based on either the authors' experience or the literature in other river systems without adequate justification in the context of the VMD. [6] applied the NDVI index to investigate the shoreline evolution of the VMD from 1973 to 2015 based on the Landsat data. Recently, [5] quantified long-term RB and coastal erosion in the entire VMD using Landsat data from 1989 to 2014. They used the MNDWI index because it has been widely used in similar studies in other river systems. We found that SRBED combining NDWI and MNDWI with the M-AMERL segmentation algorithm performs better than other methods among 11 combinations (Table 2; Figures 7 and 8). This is because the M-AMERL algorithm depends not only on the spectral characteristics of the pixels but also on their topological connections. This characteristic allows detecting the RB at relatively correct locations even in noisy environments (e.g., near aquaculture areas which are common in the VMD) where other approaches are not suitable. Incorporating the SMA further improves the performance of RB extraction (Figure 9). We conclude that the NDWI M-AMERL; SMA is the best RB extraction method in the VMD using the Landsat data. Considering the mean estimation error, the NDWI M-AMERL; SMA is 35-41%, 70% and 30% better than the MNDWI, NDVI and WNDWI indices, respectively ( Figure 10). Notably, the NDVI index has the poorest performance in RB extraction in the VMD. Moreover, the MNDWI and NDWI using the default thresholding and Otsu methods do not yield highly accurate RB positions.
Applying our developed SRBED methodology, we quantify RB changes along 100 km in the Tien River from Tan Chau to My Thuan from 2008 to 2014. It is estimated that the RB in this segment is dominated by erosion. The total net RB erosion area is approximately -1.5 km 2 with a mean erosion width of −2.64 m. These values are equivalent to annual rates of −0.25 km 2 /year and −0.44 m/year, respectively. It is revealed that our estimated values in 2008-2014 are higher than the RB erosion estimated by [5] in 1989-2014. It is likely due to the accelerated RB erosion in recent years compared to the past because of increased sand mining [38,[46][47][48][49] and decreased sediment load caused by upstream dams [44,[62][63][64] and a shift in tropical cyclones [65]. Similarly, [4] estimated that the riverbed incision volume in the upper Tien River from Tan Chau to My Thuan in 2014-2017 is approximately three times greater than that in 1998-2008, of which 85% is attributed to sediment reduction (from 166.7 Mt/year in the pre-1992 period to 43.1 Mt/year in 2012-2015-reduced by 74%). They claimed such a significant sediment reduction is caused dominantly by upstream dams, of which more than 50% is associated with six mega-dams in the upper Mekong River. Other drivers, including flow regime alterations due to dams and dyke systems and ship waves from navigation may also contribute to RB erosion in the VMD [5]. Therefore, the sustainability of the VMD requires an integrated river basin management, involving political negotiation among all riparian countries along the Mekong River and holistic and active water resources planning within the VMD.

Conclusions
Similar to the world's large mega deltas, the VMD faces various environmental pressures, among which RB erosion is usually overlooked. As such, RB erosion incidents have constantly exacerbated in the last decade that critically threaten the livelihoods of thousands of local people. Among various methods in RB erosion estimation, remote sensing of satellite data has shown advantages in long-term assessment of morphological changes in large-scale studies, of which Landsat has been widely used. Although there have been some studies analyzing RB and coastal erosion processes in the delta using satellite data, the applied methodology in previous studies was not adequately validated, and the outcomes were not quantitatively assessed. Therefore, we develop the SRBED method and RB change detection algorithm using Landsat data for RB assessment in the VMD, which can be applied in other river deltas in the world having similar characteristics.
The results show an outperformance of water indices (e.g., NDWI) compared to the vegetation index (i.e., NDVI). Using the M-AMERL algorithm, NDWI and MNDWI are better than other indices. We found that the NDWI M-AMERL; SMA (i.e., NDWI index combined with our developed M-AMERL segmentation method accompanied by the SMA algorithm) is the best RB extraction method in the VMD using the Landsat data. The NDWI M-AMERL; SMA is 35-41%, 70% and 30% better than the MNDWI, NDVI, and WNDWI indices, respectively. Moreover, we found an inappropriate application of the NDVI index in RB erosion assessment in the delta. The results of this study provide a useful reference for ongoing studies in the VMD.
Applying the developed SRBED and RB change detection algorithm, we estimated a net erosion area of −1.5 km 2 from 2008 to 2014 in the Tien River from Tan Chau to My Thuan, with a mean erosion width of −2.64 m and maximum erosion widths exceeding 60 m in places. We reveal an accelerated RB erosion in 2008-2014 compared to that in 1989-2014 by [5] which is likely due to increased sand mining activities, upstream dams, and a shift in tropical cyclones.
Influenced by a meso-tidal regime, tidal water level fluctuations during the dry season may strongly affect the accuracy of RB detection in the VMD where wide tidal flats are common. Therefore, future studies should take into consideration the effects of water levels on RB changes before drawing any conclusion. Additionally, a comprehensive investigation of the effects of water levels on RB detection accuracy is recommended.