Subpixel Surface Water Extraction (SSWE) Using Landsat 8 OLI Data

Surface water extraction from remote sensing imagery has been a very active research topic in recent years, as this problem is essential for monitoring the environment, ecosystems, climate, and so on. In order to extract surface water accurately, we developed a new subpixel surface water extraction (SSWE) method, which includes three steps. Firstly, a new all bands water index (ABWI) was developed for pure water pixel extraction. Secondly, the mixed water–land pixels were extracted by a morphological dilation operation. Thirdly, the water fractions within the mixed water–land pixels were estimated by local multiple endmember spectral mixture analysis (MESMA). The proposed ABWI and SSWE have been evaluated by using three data sets collected by the Landsat 8 Operational Land Imager (OLI). Results show that the accuracy of ABWI is higher than that of the normalized difference water index (NDWI). According to the obtained surface water maps, the proposed SSWE shows better performance than the automated subpixel water mapping method (ASWM). Specifically, the root-mean-square error (RMSE) obtained by our SSWE for the data sets considered in experiments is 0.117, which is better than that obtained by ASWM (0.143). In conclusion, the SSWE can be used to extract surface water with high accuracy, especially in areas with optically complex aquatic environments.


Introduction
Surface water provides water resources for various human uses [1] and is critical to floral and faunal habitat [2]. With the ever-increasing human population and accelerated social development, persistent stress is exerted on surface water resources [3], as surface water is essential for many applications, such as sustained water resource management [4][5][6][7]. However, conventional in situ measurements of the extent of surface water are limited in terms of spatial coverage and frequency [8].
In the past years, remote sensing has been an efficient technique for the extraction of surface water. It can provide data at various temporal scales, with large spatial coverage, low cost, among other benefits. Specifically, many satellite sensors with different spatial, temporal, and spectral resolution have been applied for the extraction of surface water, such as the Moderate-Resolution Imaging Spectroradiometer (MODIS) [9,10], Worldview [11], Satellite Pour l'Observation de la Terre (SPOT) [12], mixed water-land pixels by a morphological dilation operation; and (3) developing a local MESMA approach to dynamically extract multiple water endmembers and identify the water fractions in mixed water-land pixels.
The remainder of this work is organized as follows. Section 2 introduces three test sites (and their reference maps) that have been considered for test and validation purposes. Section 3 introduces the proposed SSWE method, as well as the adopted accuracy metrics. Section 4 presents the results obtained by the proposed ABWI and SSWE in the task of surface water extraction. Comparisons with other methods are also included in this section. Section 5 discusses the results and hints at plausible future research lines. Finally, Section 6 concludes the paper, with some remarks.

Test Sites
The three test sites considered in experiments are located in Guangdong Province, China. They comprise different surface water types, landscapes, and mixed water pixels, including rivers, reservoirs, and estuaries that differ in turbidity, depth, and the level of water pollution. These sites were deliberately selected so that the subscenes contain water bodies with diverse spectra. Figure 1 illustrates the location of the three test sites and their false color images. The first data set ( Figure 1b) is the Dongjiang River, which is a complex aquatic environment under an urban background, polluted by domestic sewage and waste water from manufactories. The second data set ( Figure 1c) is the Dongbao Estuary, which contains diverse water bodies, including turbid rivers with urban surfaces, shallow mulberry fish ponds, tidal creeks, and the ocean. The last test site (Figure 1d) is the Tiegang Reservoir, characterized by the presence of clear water with many branches. Water-land mixed pixels exist in all these test sites, among which the number of water-land mixed pixels in Dongbao Estuary is the most significant, owing to the presence of many mulberry fish ponds.
Water 2018, 10, x FOR PEER REVIEW 3 of 14 MESMA approach to dynamically extract multiple water endmembers and identify the water fractions in mixed water-land pixels. The remainder of this work is organized as follows. Section 2 introduces three test sites (and their reference maps) that have been considered for test and validation purposes. Section 3 introduces the proposed SSWE method, as well as the adopted accuracy metrics. Section 4 presents the results obtained by the proposed ABWI and SSWE in the task of surface water extraction. Comparisons with other methods are also included in this section. Section 5 discusses the results and hints at plausible future research lines. Finally, Section 6 concludes the paper, with some remarks.

Test Sites
The three test sites considered in experiments are located in Guangdong Province, China. They comprise different surface water types, landscapes, and mixed water pixels, including rivers, reservoirs, and estuaries that differ in turbidity, depth, and the level of water pollution. These sites were deliberately selected so that the subscenes contain water bodies with diverse spectra. Figure 1 illustrates the location of the three test sites and their false color images. The first data set ( Figure 1b) is the Dongjiang River, which is a complex aquatic environment under an urban background, polluted by domestic sewage and waste water from manufactories. The second data set ( Figure 1c) is the Dongbao Estuary, which contains diverse water bodies, including turbid rivers with urban surfaces, shallow mulberry fish ponds, tidal creeks, and the ocean. The last test site (Figure 1d) is the Tiegang Reservoir, characterized by the presence of clear water with many branches. Water-land mixed pixels exist in all these test sites, among which the number of water-land mixed pixels in Dongbao Estuary is the most significant, owing to the presence of many mulberry fish ponds.

Landsat 8 OLI Data
For the three considered test sites, orthorectified and terrain-corrected Landsat 8 OLI images in GeoTIFF format were obtained from the United States Geological Survey (USGS) EarthExplorer [48], all of them of product type Level 1 T. All Landsat 8 OLI images were defined in Universal Transverse Mercator (UTM) map projection and georeferenced with a precision better than 0.4 pixels [49]. Additional details, including the acquisition date, cloud cover, and scene ID of the OLI images are given in Table 1. All Landsat images were clipped to 230 × 230 pixels because the sizes of Landsat images were much larger than that of the reference data. All subsets were cloud-free to decrease the influence of clouds on the data. These images, acquired on specific dates, were selected because high-resolution images acquired on the same dates for reference can be found in Google Earth™. The acquisition dates of all the Landsat 8 OLI images were in the dry season. On the acquisition dates of the images used for the three test sites, days from the last rainfall event were 6, 4, and 6, respectively. The Landsat 8 OLI image at Dongbao Estuary was acquired during low tide.

Reference Data
For the three test sites, the fine spatial resolution Google Earth™ images (<1 m) were used for reference. The acquisition dates of the reference data and the Landsat 8 OLI images were exactly the same. From Google Earth, we can obtain the acquisition dates of the Google Earth™ images, however, the acquisition time of these images were unavailable. The short time span between the Google Earth™ images and the Landsat 8 OLI images was inevitable. Generally, the surface water boundaries change little during a short time period, though the tide may influence surface water extent at Dongbao Estuary, where the changes in water level were less than 0.4 m from 9 a.m. to 12 p.m. on the image acquisition date. The Landsat 8 OLI images were registered to the Google Earth™ images using a first-order polynomial with 10 ground control points (GCPs) manually selected from each image in ENVI v.5.3 [50]. The GCPs were located mainly at road junctions, and the corresponding root-mean-square error (RMSE) values for the three test sites were 0.21, 0.26, and 0.28 pixels, respectively. The "true" boundaries of surface water bodies were obtained via manual digitization on-screen from the reference data. Then, the "true" surface waterbody maps were overlaid on the 30 m grids of the Landsat data to generate the "true" water fraction images at the 30 m resolution, which were used to assess the accuracy of the different water extraction methods considered in this work.

Image Preprocessing
The raw Landsat 8 OLI images are in the form of digital numbers (DNs), which often need to be converted into apparent surface reflectance before processing. In this work, radiometric calibration and atmospheric correction were applied to all used images. Radiometric calibration converted the original DNs into top of atmosphere (TOA) radiances R λ [W/(m 2 ·sr·µm)]. For a given wavelength λ, the TOA radiance R λ was derived from the following equation: where G λ is the gain factor and O λ is the offset, which were provided in the Level 1 Product Generation System (LPGS) metadata file. Atmospheric correction was applied to the radiometrically calibrated images using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module in ENVI v.5.3 [50]. The atmospheric parameters were obtained from a look-up table based on a seasonal-latitude surface temperature module [51]. Aerosol optical depth (AOD) data was obtained from the Level 2.0 AERONET data [52]. Though many approaches existed to correct for atmospheric effects with high accuracy, FLAASH was selected because it was commonly used in many recent surface water extraction researches based on Landsat 8 OLI data [16,42].

Subpixel Surface Water Extraction
The subpixel surface water extraction (SSWE) method includes three major steps. Firstly, an all bands water index (ABWI) was developed to classify the pure water and land classes. In this step, the water-land mixed pixels were classified as land. Secondly, we used a morphological dilation operation to identify the mixed water-land pixels. Thirdly, a local MESMA procedure was developed for unmixing the mixed pixels. The details of the SSWE are given below.

Extraction of Pure Water Pixels
Previous studies have revealed the spectral signatures of major land cover types [22,53]. It was found that the water shows higher reflectance of visible bands than infrared bands in most cases. Few non-water surface possesses this similar spectral characteristic. Therefore, in order to maximize separability of water and non-water pixels, we designed a new water index taking full advantage of the multispectral bands of Landsat 8 OLI sensor. Then, the new method was called all bands water index (ABWI), which can be obtained from the Landsat 8 OLI images as follows: where ρ 1 represents the reflectance of the coastal band, ρ 2 represents the reflectance of the blue band, ρ 3 represents the reflectance of the green band, ρ 4 represents the reflectance of the red band, ρ 5 represents the reflectance of the NIR band, ρ 6 represents the reflectance of the SWIR1 band, and ρ 7 represents the reflectance of the SWIR2 band. The pure water pixels can be extracted automatically based on the histogram of ABWI using the method of Xie et al. [42]. However, this method may not always achieve the highest accuracy of pure water pixel extraction. In order to achieve the highest possible water extraction accuracy, the performance of the developed ABWI was tested in the optimal threshold. Similar to Feyisa et al. [22], multiple thresholds were considered to determine the optimal threshold for ABWI. The corresponding commission errors and omission errors of each threshold were plotted on a graph, on which the intersection point was regarded as the optimal, in that it approximates the minimum number of total errors. Specifically, 1500 pure water pixels and 1500 non-water pixels were randomly selected from the reference data to determine the optimal threshold, and the rest of the samples were used for accuracy assessment for each test site. The optimal threshold of ABWI for the three test sites were 0.66, 0.59, and 0.47, respectively.

Mixed Water-Land Pixels Extraction
Given the fact that the mixed waterland pixels are generally located in the transitions from water to land, these pixels are expected to be found on the boundaries of water and land classes. Therefore, these pixels can be determined based on their spatial relationship to the pure water pixels. A morphological dilation operation was used to exploit this spatial characteristic to extract the mixed pixels. Conceptual graph of extraction of mixed water-land pixels and selection of water endmembers is shown in Figure 2. This is a conceptual graph, and Google Earth™ data were used only as reference data, with the OLI data being the only reflectance input data for the SSWE method. As shown in Figure 2a, a grid is superimposed on a high spatial resolution image obtained through Google Earth™. The mixed water-land pixels were determined based on the water and land classification result shown in Figure 2b. For example, as shown in Figure 2c, among the 8-pixel neighborhood connected to a pure water pixel located at (i, j), all the land pixels were considered to be mixed pixels. In this work, a 3 × 3 moving window was used to search all the mixed water-land pixels.

Local Multiple Endmember Spectral Mixture Analysis
After the determination of the mixed water-land pixels, the image can be classified into three classes: water, land, and water-land mixed pixels. As shown in Figure 2d, the main confusion is precisely in the water-land mixed pixels. SMA can then be used to estimate water fractions in these pixels. The commonly used linear constrained spectral mixture model assumes that the spectral signature of a mixed pixel can be represented as the linear sum of a set of spectrally pure endmembers, weighted by their corresponding fractions [54]. That is, for a given wavelength, λ: where is the observed reflectance at the mixed pixel i, N is the number of endmembers, is the fraction of the jth endmember at pixel i, is the reflectance of the jth endmember, and is a residual term indicating the disagreement between the observed and modeled reflectance. The fractions are typically nonnegative and satisfy the following constraint: Subsequently, for the pixel i, the fractions can be estimated by minimizing the root-mean-square error (RMSE) of the residuals across all bands. The RMSE can be obtained as follows: where T is the number of bands, is the residual at pixel i with band k. For the surface water, there are more severe spectral variabilities than for other land cover types. This is because the spectral features of surface water are influenced by depth and the aquatic

Local Multiple Endmember Spectral Mixture Analysis
After the determination of the mixed water-land pixels, the image can be classified into three classes: water, land, and water-land mixed pixels. As shown in Figure 2d, the main confusion is precisely in the water-land mixed pixels. SMA can then be used to estimate water fractions in these pixels. The commonly used linear constrained spectral mixture model assumes that the spectral signature of a mixed pixel can be represented as the linear sum of a set of spectrally pure endmembers, weighted by their corresponding fractions [54]. That is, for a given wavelength, λ: where R iλ is the observed reflectance at the mixed pixel i, N is the number of endmembers, f ji is the fraction of the jth endmember at pixel i, R jλ is the reflectance of the jth endmember, and e iλ is a residual term indicating the disagreement between the observed and modeled reflectance. The fractions are typically nonnegative and satisfy the following constraint: Subsequently, for the pixel i, the fractions can be estimated by minimizing the root-mean-square error (RMSE) of the residuals across all bands. The RMSE can be obtained as follows: where T is the number of bands, e ik is the residual at pixel i with band k.
For the surface water, there are more severe spectral variabilities than for other land cover types. This is because the spectral features of surface water are influenced by depth and the aquatic environment. As shown in Figure 3, the spectral features of 13 water pixels varied greatly from pixel to pixel. In order to tackle the high variability of water endmembers, we used a local spatial constraint to dynamically extract the water endmembers. Specifically, the pure water pixels in the nearest spatial neighborhood were considered as the water endmembers for a mixed water-land pixel. Through experimental comparison, the optimal window size was set to 3 × 3. In this work, we used an 8-pixel connected neighborhood system to identify the water endmembers. As shown in Figure 2d, in the 8-pixel neighborhood connected to a water-land mixed pixel located at (k, l), all the pure water pixels were considered to be the water endmembers. neighborhood were considered as the water endmembers for a mixed water-land pixel. Through experimental comparison, the optimal window size was set to 3 × 3. In this work, we used an 8-pixel connected neighborhood system to identify the water endmembers. As shown in Figure 2d, in the 8pixel neighborhood connected to a water-land mixed pixel located at (k, l), all the pure water pixels were considered to be the water endmembers. The local MESMA was implemented as follows. Firstly, a spectral endmember library was created by extracting three categories of land endmembers (i.e., vegetation, soil, impervious surfaces), and was derived from Landsat 8 OLI images by using the pixel purity index (PPI) in ENVI v.5.3 [50]. Pixels identified as extreme by the PPI were visually inspected, and the spectra of those which could be classified with confidence as one of the above three land cover components were collected. Countbased endmember selection (CoB), endmember average RMSE (EAR), and minimum average spectral angle (MASA) [45] were used to choose the most suitable endmembers. Secondly, for each mixed water-land pixel, we dynamically selected the water endmembers from Landsat 8 OLI images using a 3 × 3 moving window, where the pure water pixels in the nearest spatial neighborhood were selected as water endmembers for the given mixed pixel. The numbers of land cover endmembers for the three test sites are presented in Table 2. Lastly, we adopted multiple endmember spectral unmixing (based on the linear spectral mixture model) to unmix the mixed water-land pixels. Here, combinations of land and water endmembers were considered, including vegetation-water-shade (V-W-shade), soil-water-shade (S-W-shade), impervious-surfaces-water-shade (I-W-shade), V-S-W-shade, V-I-W-shade, S-I-W-shade, and V-S-I-W-shade. The shade endmember was included in each endmember combination to account for the variation in illumination. Similar to previous studies [47], the following selection criteria were adopted: the fraction values should be between −0.05 and 1.05; the shade fraction value should be below 0.8; and the RMSE should be below 0.025. Among the eligible combinations, the one with the lowest RMSE was considered to be the best endmember combination. If no endmember combination met the RMSE constraint, the pixel was neglected, and the water fraction of the pixel was set to 0. The local MESMA was implemented as follows. Firstly, a spectral endmember library was created by extracting three categories of land endmembers (i.e., vegetation, soil, impervious surfaces), and was derived from Landsat 8 OLI images by using the pixel purity index (PPI) in ENVI v.5.3 [50]. Pixels identified as extreme by the PPI were visually inspected, and the spectra of those which could be classified with confidence as one of the above three land cover components were collected. Count-based endmember selection (CoB), endmember average RMSE (EAR), and minimum average spectral angle (MASA) [45] were used to choose the most suitable endmembers. Secondly, for each mixed water-land pixel, we dynamically selected the water endmembers from Landsat 8 OLI images using a 3 × 3 moving window, where the pure water pixels in the nearest spatial neighborhood were selected as water endmembers for the given mixed pixel. The numbers of land cover endmembers for the three test sites are presented in Table 2. Lastly, we adopted multiple endmember spectral unmixing (based on the linear spectral mixture model) to unmix the mixed water-land pixels. Here, combinations of land and water endmembers were considered, including vegetation-water-shade (V-W-shade), soil-water-shade (S-W-shade), impervious-surfaces-water-shade (I-W-shade), V-S-W-shade, V-I-W-shade, S-I-W-shade, and V-S-I-W-shade. The shade endmember was included in each endmember combination to account for the variation in illumination. Similar to previous studies [47], the following selection criteria were adopted: the fraction values should be between −0.05 and 1.05; the shade fraction value should be below 0.8; and the RMSE should be below 0.025. Among the eligible combinations, the one with the lowest RMSE was considered to be the best endmember combination. If no endmember combination met the RMSE constraint, the pixel was neglected, and the water fraction of the pixel was set to 0.

Accuracy Assessment
The kappa coefficient (KC) and total error (TE) were calculated to assess the classification accuracy of the developed ABWI. For comparison, the accuracies of NDWI were also obtained. A preliminary test of different water indices, which included the normalized difference water index (NDWI) [23], the modified NDWI (MNDWI) [24], and the automated water extraction index (AWEI) [22], showed that all indices, except the NDWI, had poor performance at our test sites. Therefore, only NDWI was considered for comparison with the ABWI. The ABWI and NDWI were all tested in the optimal threshold.
Root-mean-square error (RMSE) and systematic error (SE) were used to evaluate the accuracy of the water fraction estimated by the proposed SSWE method. The RMSE quantifies the relative error of the estimated land cover fraction, and the SE is used to evaluate the overall trend of the over-or underestimation. For comparison, the accuracy of the automated subpixel water mapping method (ASWM) was also obtained. In order to conduct a comprehensive assessment of the SSWE method, all the pixels, except the 3000 training samples, were used for validation. As the whole sample size was 62,900 in each test site, the validation sample size was 59,900. Here, the water fractions in the pure water pixels and pure land pixels were set to 1 and 0, respectively. The RMSE and SE can be obtained as follows: where f i (re f ) and f i (esti) are the water fraction of the ith pixel of the reference image and estimated image, respectively. N is the number of pixels in the whole image. Figure 4 shows the pure water pixel extraction results obtained by the proposed ABWI for the three test sites. For comparison, the results obtained by the NDWI are also reported. The visual interpretation indicates that the ABWI shows good performance in extracting pure water pixels. The ABWI is more effective for distinguishing water and mixed water-land pixels than the NDWI (Region 1 and Region 3). Particularly, in the second test set with the Dongbao Estuary, the new index (ABWI) is better at suppressing mud flat and other non-water surface (Region 2).

Extraction of Pure Water Pixels
The accuracies for ABWI and NDWI over the three test sites are shown in Table 3, from which it can be clearly observed that the ABWI achieves better accuracy than the NDWI at all test sites. Overall, the means of the KC and TE achieved by ABWI over the three test sites are 0.957 and 7.56%, respectively, which improves the results provided by NDWI (KC = 0.934, TE = 10.26%). Especially at Dongbao Estuary (with significant presence of mixed water-land pixels and mud flats), the ABWI achieves the largest increments of KC (from 0.912 to 0.945) and the greatest decrements of TE (from 11.67% to 8.16%). Since the NDWI considers only a few bands, the developed all bands water index (ABWI) performs better than the NDWI in the extraction of pure water pixels.   Figure 5 shows the subpixel surface water extraction maps obtained by the proposed SSWE for the three test sites. For comparison, the results obtained by the ASWM method are also reported. Several observations can be obtained according to visual inspection of the obtained results. First of all, it is remarkable that, for the mixed water-land pixels, the SSWE properly extracts the fractions of water in these pixels, which shows clear transitions from water to non-water (Region 1 and Region 3). Furthermore, the SSWE shows very good performance in highly mixed areas, such as the estuary (Region 2). These areas are generally difficult to characterize with the ASWM, mainly due to the complex spectra present in those areas. Finally, it can be observed that both methods fail to identify too narrow rivers and too small waterbodies (width less than 2 pixels). In order to further evaluate the performance of the proposed SSWE, the root-mean-square error (RMSE) and systematic error (SE) were obtained for the three test sites. For comparative purposes, the results obtained from ASWM are also reported. Table 4 indicates that the SSWE reduces RMSE, as compared to the ASWM, for all the three test sites. Especially at Dongbao Estuary (with significant presence of mixed water-land pixels and mud flats), the SSWE achieves the greatest decrements of RMSE (from 0.214 to 0.163).  Figure 5 shows the subpixel surface water extraction maps obtained by the proposed SSWE for the three test sites. For comparison, the results obtained by the ASWM method are also reported. Several observations can be obtained according to visual inspection of the obtained results. First of all, it is remarkable that, for the mixed water-land pixels, the SSWE properly extracts the fractions of water in these pixels, which shows clear transitions from water to non-water (Region 1 and Region 3). Furthermore, the SSWE shows very good performance in highly mixed areas, such as the estuary (Region 2). These areas are generally difficult to characterize with the ASWM, mainly due to the complex spectra present in those areas. Finally, it can be observed that both methods fail to identify too narrow rivers and too small waterbodies (width less than 2 pixels). In order to further evaluate the performance of the proposed SSWE, the root-mean-square error (RMSE) and systematic error (SE) were obtained for the three test sites. For comparative purposes, the results obtained from ASWM are also reported. Table 4 indicates that the SSWE reduces RMSE, as compared to the ASWM, for all the three test sites. Especially at Dongbao Estuary (with significant presence of mixed water-land pixels and mud flats), the SSWE achieves the greatest decrements of RMSE (from 0.214 to 0.163). Overall, the mean of the RMSE and SE achieved by SSWE over the three test sites is 0.117 and −0.005, respectively, which improves the results provided by ASWM (RMSE = 0.143, SE = 0.008). Negative SE values are observed for SSWE over all the test sites, indicating that the overall water fraction is underestimated.

Discussion
Surface water extraction is a complex problem affected by many factors, especially in estuaries. Since estuaries form a transition zone between river environments and maritime environments, the aquatic environments are complex and highly variable. Mud flats and suspended particulate sediments in estuaries show similar reflectance spectra with their adjacent water pixels, making the water boundaries optical-confused. Meanwhile, there is no distinct rule in the distribution of Figure 5. Comparison of surface water extraction results of the subpixel surface water extraction (SSWE) and the automated subpixel water mapping method (ASWM) for the three test sites. The color range indicates the range of water fractions (i.e., blue indicates pure water pixels and red indicates pixels with low water fractions). Dark indicates non-water pixels.

Discussion
Surface water extraction is a complex problem affected by many factors, especially in estuaries. Since estuaries form a transition zone between river environments and maritime environments, the aquatic environments are complex and highly variable. Mud flats and suspended particulate sediments in estuaries show similar reflectance spectra with their adjacent water pixels, making the water boundaries optical-confused. Meanwhile, there is no distinct rule in the distribution of suspended particulate sediments because of the complicated hydrodynamic conditions in estuaries, which result in high spectral variability of surface water. The results indicate that the SSWE shows better performance than the ASWM in Dongbao Estuary. It can be explained as follows. Firstly, the ABWI using information from all the bands is more effective to distinguish water and non-water pixels than the NDWI, which is used in the ASWM. Secondly, it is difficult for ASWM to extract water-land mixed pixels based on only the spectral information. It seems that SSWE, based on both spectral and spatial information, is a better choice to extract water-land interface mixed pixels. Lastly, the ASWM uses a fixed combination of endmembers and it is limited when mapping land cover with high spectral variability. Thus, the applications of ASWM are limited by high spectral variability of surface water in Dongbao Estuary. As a result, the SSWE integrating the advantages of ABWI and local MESMA appears to be more appropriate in areas that include diverse water types and where the reflectance spectra of surface water are highly variable.
Despite a number of surface water extraction methods reported in the literature, limited research has been undertaken on accuracy assessment using high spatial resolution reference images acquired on the same day with Landsat 8 OLI data. This is particularly important when a subpixel water extraction method is developed. Given that water conditions are sometimes changing rapidly, short or no time span between reference data and Landsat 8 OLI data is essential for producing the "true" water map, especially when the bank is not artificial. For instance, in coastal areas, the surface water changes with the tide. If the reference image is not acquired at the same time as Landsat 8 OLI data, a bias in the surface water boundaries could arise, since tides influence the water level. However, the short time span between reference data and Landsat 8 OLI data is inevitable. In order to minimize the bias in the surface water boundaries and guarantee a convincing validation, the SSWE is evaluated using reference data acquired on the same day with Landsat 8 OLI data.
In the task of surface water extraction from a large data set of images, an automated process is preferred. Though only the second step of SSWE is automatically implemented in this work, the three-step strategy in SSWE has the potential to perform automated water extraction. First of all, the pure water pixels can be extracted automatically based on the histogram of ABWI using the method of Xie et al. [42]. This method was tested in the ABWI image histogram of our three test sites to extract pure water pixels, and the corresponding kappa coefficients were a little less than that of the optimal threshold. Furthermore, in the third step of SSWE, the land endmembers are able to be extracted automatically from the image with the land mask produced from the first two steps by using the methods proposed by Somers et al. [55]. Besides, in situ data is not required in the SSWE method. Thus, the whole procedure of SSWE can be automated. Nevertheless, the performance of the automated process needs to be evaluated.
Although our new water extraction method has been tested in three test sites with different surface water types and achieves promising performance, further studies should be conducted in order to improve the robustness of the proposed method. At this stage, the problem of too narrow rivers and too small waterbodies (width less than 2 pixels) has not yet been taken into account. Despite that methods reported in the literature [28,36] can partially detect these pixels, extraction of these pixels is still a challenge. To accurately identify water fraction in these pixels, it is crucial to collect the most representative water endmember. However, because of little or no pure water pixels adjacent to these pixels, the selection of the most representative water endmember remains a difficulty. In the future, we will extend the capacity of the SSWE to extract these waterbodies. The ABWI will be also tested in other water types for water extraction purpose. Furthermore, since the SSWE was designed on the basis of the land cover reflectance and local spatial constraints, its application can be extended to satellite image types with similar spectral bands. The SSWE will be tested in additional sites and other satellite image types with similar spectral bands, such as Sentinel-2 MSI images, to validate its performance and robustness.

Conclusions
In this study, we developed a new subpixel surface water extraction (SSWE) method for Landsat 8 OLI data that is able to improve the extraction of surface water in mixed water-land pixels by means of subpixel analysis. Specifically, we first developed an all bands water index (ABWI) to extract the pure water pixels. The ABWI takes full advantage of the multispectral bands and mines more information than the previous water indices. Results indicated that the ABWI performed reasonably well in pure water pixel extraction. Then, we developed a local multiple endmember spectral mixture analysis (MESMA) to identify water fractions at mixed water-land pixels. By incorporating a local MESMA, the proposed SSWE dynamically characterizes the high spectral variability of the water endmembers and complex land classes. The proposed method considered the nearest pure water pixels as water endmembers, and further allowed the number and type of endmembers to vary on a per-pixel basis. Our experimental results with three different test sites in China indicated that the newly developed SSWE exhibited very good performance, with an average root-mean-square error of 0.117. A thorough comparison with the automated subpixel water mapping method (ASWM) indicated that our SSWE outperformed the method in most cases. Especially at Dongbao Estuary (with significant presence of mixed water-land pixels and mud flats), the SSWE achieved the greatest decrements of RMSE. Therefore, the SSWE can be used to extract surface water with high accuracy, especially in areas with optically complex aquatic environments.
Author Contributions: L.X. designed and conducted the overall research, experiment and analysis, and prepared the manuscript. R.D. and J.L. led this research and reviewed the manuscript. X.L. supported the reference data collection. Y.Q. helped perform the analysis with constructive discussions. Y.L. supported the Landsat data collection. Y.L. joined the experiment. All authors wrote the paper.