Application of Image Segmentation in Surface Water Extraction of Freshwater Lakes using Radar Data

: Freshwater lakes supply a large amount of inland water resources to sustain local and regional developments. However, some lake systems depend upon great ﬂuctuation in water surface area. Poyang lake, the largest freshwater lake in China, undergoes dramatic seasonal and interannual variations. Timely monitoring of Poyang lake surface provides essential information on variation of water occurrence for its ecosystem conservation. Application of histogram-based image segmentation in radar imagery has been widely used to detect water surface of lakes. Still, it is challenging to select the optimal threshold. Here, we analyze the advantages and disadvantages of a segmentation algorithm, the Otsu Method, from both mathematical and application perspectives. We implement the Otsu Method and provide reusable scripts to automatically select a threshold for surface water extraction using Sentinel-1 synthetic aperture radar (SAR) imagery on Google Earth Engine, a cloud-based platform that accelerates processing of Sentinel-1 data and auto-threshold computation. The optimal thresholds for each January from 2017 to 2020 are − 14.88, − 16.93, − 16.96 and − 16.87 respectively, and the overall accuracy achieves 92% after rectiﬁcation. Furthermore, our study contributes to the update of temporal and spatial variation of Poyang lake, conﬁrming that its surface water area ﬂuctuated annually and tended to shrink both in the center and boundary of the lake on each January from 2017 to 2020.


Introduction
Water is significant for all ecosystems on Earth. The presence of surface water on Earth mainly consists of oceans, lakes and rivers [1]. The extent of lakes accounts for nearly 3% of the surface [2] and is endowed with irreplaceable functions to supply water [3], control flooding [4], sustain species [5] and provide ecosystem services to nations and regions [6] due to the unique role of water in climate [7], biological diversity [8] and human wellbeing [9]. Meanwhile, natural phenomena and human activities affect the variation of water occurrence in response, especially the water dynamics of inland freshwater lakes [10]. Timely monitoring of freshwater lake surface is indispensable for sustainable development [11] and regional and global ecosystem dynamics [12].
Remote sensing, the science and art of detecting objects from a distance, has been the most common approach to monitor and analyze land features for several decades [13]. In imagery, land features are is to maximize conditional probability of the predictionŷ given any dB values (x). According to the Bayes Theorem, this equates to maximizing the multiplication of the conditional probability of x overŷ and the marginal probability ofŷ.
However, the issue with such formulation of the problem is based on the assumption of the prior distribution of water and land as a Gaussian Distribution. However, such an assumption cannot be directly assumed to be correct for universal cases. Moreover, the distribution parameters µ 1 , σ 1 , µ 2 , σ 2 are unknowns. The researcher also needs to identify estimators for these four parameters through the density diagram. Possible solutions for estimation of these unknown parameters can be iterative methods such as Expectation Maximization Method [24], however, it is unstable for two reasons. First, the iteration process is time consuming to reach a satisfied accuracy. Second, it is also likely to be constrained in some local optimum points and thus never reaches global optimal solution [25].
Instead, we propose to use the Otsu Method to solve this thresholding problem. The Otsu Method is an unsupervised method and it was initially designed to select a threshold to separate an object out of its background, through the gray-level histogram of the image [25]. In application, the Otsu Method can be widely extended to work on other density histograms or distributions other than gray-level histogram from images and can also be applied for multi-thresholding problems. The Otsu Method is a better approach for this problem as compared to some conventional solutions because it automatically selects a threshold from two mixed distributions through the density histogram [25]. In addition, the Otsu Method does not require prior knowledge nor assumptions of the distribution of objects [25]. Furthermore, the Otsu Method is equivalent to the K-Means Method but the Otsu Method can provide the global optimal solution, while K -Means Method may be limited to the local optimum point [26]. Although it is computationally complex and heavy because of iterative searching [26], GEE can speed up the Otsu Method with its cloud computing platform. For instance, the Otsu Method has been applied on the cloud-free Landsat TM images for urban land cover detection, which focused on differentiating the urban land and nonurban land region in Haidian District of Beijing, China [27]. This research resulted in an accuracy of 84.83% for the Otsu Method, which was larger than the accuracy of 74% for the conventional postclassification change detection method [27]. Another study used the Otsu Method on the SAR data for the detection of oil spills over sea surfaces, which tried to find a threshold on the radar data to draw the edge of spilled oil film floating over the sea [28]. It examined the Penglai oil field and the Gulf of Dalian, resulting in an error rate of 3.0% on the Penglai oil field and an error rate of 13.0% on the Gulf of Dalian for the Otsu Method [28]. Even though the Otsu Method has already been widely applied in thresholding problems, it has been seldom used for surface water extraction. Furthermore, most previous studies do not provide algorithms and detailed scripts for implementation of the Otsu Method. Thus, we were interested in the application of the Otsu Method for surface water detection and providing reusable code for future implementation.
Therefore, the objectives of the present work are to: 1. Implement the Otsu Method and write reusable scripts to automatically select thresholds for surface water extraction using Sentinel-1 data on Google Earth Engine 2. Analyze the advantages and disadvantages of an unsupervised classifier from both mathematical and application perspectives 3. Contribute to the knowledge base of hydrological variation at Poyang lake by mapping surface water extent of the lake in January 2017, 2018, 2019 and 2020

Study Area
Poyang lake, the largest freshwater lake in China, is located between Nanchang City and Jiujiang City, to the north of Jiangxi Province. The basin crosses from 28 • 22 to 29 • 45 N and 115 • 47 to 116 • 45 E, which belongs to a humid, subtropical monsoon climate zone, with an average annual temperature of 17.5 • C and average annual precipitation of 1665 mm [29]. Poyang lake basin is fed by the Xiu, Gan, Fu, Xin and Rao rivers, while the basin connects to the Yangtze river through an outflow channel at the north end of the lake ( Figure 1). The lake has a surface area of approximately 4000 square kilometers at its summer high-water level [30,31]. Beyond its size, Poyang lake is also significant for several economic and ecological reasons. For instance, Poyang lake's aquatic ecosystems are wintering home to thousands of migratory waterbirds, including the Siberian crane-a critically endangered species whose 4000 surviving individuals spend their winters almost solely in the wetlands around Poyang lake [32]. However, Poyang lake has undergone a series of significant transformations that threaten the variability and critical habitats in the region. While surface water areas have traditionally fluctuated on a seasonal scale-peaking in the summer and receding in the winter-large interannual declines in mean water level have been observed in recent years [33]. The substantial variations of the surface water area and dramatic seasonal water level fluctuations of 8 to 22 m each year are caused by the regional hydrological regime, which is controlled both by the five catchment rivers and the Yangtze River [34]. Additionally, groundwater dynamics are highly affected by the variations in the lake water level, rather than local precipitation, indicating a close hydraulic relationship between groundwater and the lake [35].

Platform and Data
Google Earth Engine (GEE, https://earthengine.google.com) consists of a multipetabyte satellite imagery data catalog colocated with a high-performance, intrinsically parallel cloud computation service. Users can access GEE through an Internet-accessible application programming interface (API) and an associated web-based interactive development environment (IDE) that enables rapid prototyping and visualization of results. This cloud computing platform not only makes it easy to access most of the geospatial datasets but also enables high throughput analysis. There are many examples where environmental scientists empowered their research with help of GEE, such as population mapping [36], cropland mapping [37], extraction of glacial lakes [38] and probabilistic wetland mapping [39].
Sentinel-1 is the first Copernicus Program satellite constellation deployed by the European Space Agency. This space mission is composed of two satellites, Sentinel-1A and Sentinel-1B, carrying a C-band synthetic-aperture radar instrument which collects data in all weather, day or night [40]. Since radar sensors have the advantage in detecting moisture and water because of their ability to penetrate clouds, Sentinel-1 is one of the most common datasets for surface water detection [41] and flood mapping [42].
The winter low-water season of Poyang lake provides important foraging habitat and wintering area for many waterbirds of special concern, including the critically endangered Siberian crane. Because of the importance of water level during this time, we looked at images taken in January over subsequent years. We loaded Sentinel-1 Level-1 IW GRD images from the data catalog of GEE from January 2017-2020 ( Table 1). The imagery acquired on January of 2020 was used to evaluate our Otsu Method implementation on GEE. The others were used to analyze the water area change in January across 4 years from 2017 to 2020.

Otsu Method
In this section, we firstly introduce the main idea of the Otsu Method [25,26] in a general framework and then we discuss how the Otsu Method is applied on this thresholding problem with a radar value density histogram.
Here we use the following notations: where i ∈ C refers to i-th element belonging to the whole set C we are considering. x i is the value for this i-th element and without loss of generality, we can assume that x i are sorted. Explicitly, p i is the probability or density of the element i. It is clear that ∑ i∈C p i = 1.
• we try to split up the set C into two disjoint subclusters of index C 0 , µ j is the center or the mean value of cluster C j , for j = 0, 1: • µ is the center or the mean value of the whole set C: • V j is denoted as the inner-variance of the cluster C j , which is defined as the weighted summation of the squared distance of cluster C j 's each data point from its center µ j , for j = 0, 1: • V 0,1 is denoted as the interclass variance between the cluster C 0 and cluster C 1 [25], which is defined as the weighted summation of the squared distance of each cluster's center µ j from the center of the whole set µ: • V is denoted as the total-variance, which is defined as the weighted summation of the squared distance of all data points from the center of the whole set µ. Furthermore, we can see that V is actually exactly the variance σ 2 C of the set C: The main idea of the Otsu Method is to minimize the summation of the inner-variance V j of all clusters C j , which is called intraclass variance [25]. The inner-variance of a cluster shows the summation of squared distance of each element to the center of the cluster as we defined, and the smaller value of the inner-variance presents the closer distance of each point toward the center of the cluster, which shows a closer relationship or higher similarity that the elements in this cluster share. Therefore, the best separation of the whole set of elements should group the similar elements in the same cluster as optimally as possible. In mathematics, this is equivalent to minimizing the summation of inner-variance inside each cluster. The objective function is formulated as follows: Furthermore, the summation of each cluster's inner-variance and the interclass variance should be equal to the total-variance of the whole set [25], which is a constant for a fixed data set.
Therefore, the previous objective function Equation (1) is equivalent to maximizing the interclass variance V 0,1 : max Now, in applying the Otsu Method on the density histogram, we can have:

•
The set of all possible bin's values on the density histogram as Θ, which is also the hypothesis space for the estimation of the threshold.

•
The density corresponding to the bin with value θ is denoted as p θ and we should have 1 = ∑ θ∈Θ p θ .
• For each bin value θ ∈ Θ, we can put a corresponding index i θ into the indexed set C, where i θ -th bin on the density histogram has a bin value of θ. Therefore, x i θ = θ and p i θ is equal to the density p θ of the bin with value θ from the density histogram.

•
The final prediction of the threshold is denoted asθ ∈ Θ, which corresponds to the index iθ ∈ C. • the two separate clusters formed by a threshold θ is The final estimation of the thresholdθ should be the one based on which the subcluster C * 0 and C * 1 can reach the optimal value of the objective function in Equation (3). Then, we can have: This optimization problem can be solved in at most quadratic polynomial time of the size of the finite set Θ, i.e., the time complexity will be O(|Θ| 2 ). One possible implementation as shown in Algorithm 1 is to iterate through the finite set and record the element in the set that provides the highest value for the objective function. Each inner iteration takes linear time to calculate the objective function.
is the set of bins for the density histogram. Output:θ is the final prediction of the optimal threshold between the two classes. Procedure: θ ← ∞ objective_value ← 0 for each bin i in the input density histogram do C 0 ← ∅ C 1 ← ∅ for each bin j in the input density histogram do if θ j < θ i then C 0 adds j else C 1 adds j end end curr ← Compute objective function value based on Equation (3)

Water Detection
Once the preprocessing procedure was completed by GEE, the histogram of VV band was generated, and the Otsu method was used to search over the thresholds that are represented by the bins in the histogram. The optimal threshold was computed to classify the data, where the partition whose values are smaller than the threshold are labeled as water while the partition whose values are larger than the threshold are labeled as nonwater. In order to reduce the effect of double bounce scattering issues, we defined the label of water as purely open water area, while the label of nonwater included the submerged and emergent aquatic vegetation and land features. The specific implementation of the Otsu Method can be found through the link: source code for the Otsu Method (by Sulong Zhou). The postprocessing procedure that removes noise and improves the quality of the classified output involving mask extraction, majority filtering and boundary clean was carried out in ArcGIS to remove the water bodies not geographically related to Poyang lake, small islands of pixels and odd edge of clusters.

Accuracy Assessment
We denote a point x with its true label y drawn from the true distribution D as (x, y) ∼ D, where the true distribution D is actually unknown. Specifically, x is the radar dB value for a pixel and defined as: Accuracy measures the agreement between a standard assumed to be correct and a classified image of unknown quality [43]. Classification errors occur when a pixel (or feature) belonging to one category is assigned to another category. Errors of omission occur when a feature is left out of the category being evaluated; errors of commission occur when a feature is incorrectly included in the category being evaluated [44]. An error of omission in one category will be counted as an error in commission in another category. Explicitly, for a pixel's dB value x and its true label y from the unknown true distribution D, the error happens whenŷ = y. Therefore, accuracy can be mathematically defined as the follow: Since the true distribution D is unknown, it is not possible to calculate the accuracy through Equation (7). Therefore, we need an estimator to estimate such accuracy. One possible way to estimate is based on the Empirical Distribution [45]. A test set T = {(x i , y i )} n i=1 , which forms an empirical distributionD, is used to approximate the true distribution D, where each element (x i , y i ) is independently and identically (i.i.d) drawn from the true distribution D.

∼ D
The empirical estimator of accuracy Accuracy for classification function fθ can be expressed as: Since (x i , y i ) ∈ T, for ∀i = 1, . . . , n, is i.i.d drawn from the true distribution D, this estimator Accuracy in Equation (8) is an unbiased estimator of the true Accuracy as defined in Equation (7) [46], proved as follows: Proof of Unbias Estimator.
A modified double-blind visual assessment of a random sample of test sites was used to assess classification accuracy. Firstly, a random set of 304 test sites was generated across the region, and the algorithm can be found through the link: source code for random points (by Sulong Zhou). This random set corresponds to the test set T with n = 304 and each element (x i , y i ) is i.i.d drawn from the true distribution D, which is the distribution of dB value and label of the locations in the study area as shown in Figure 1. Next, it was assigned to a team in Nanchang who visited all accessible points from the set of 304. Then based on their experience, knowledge, and observation in both real field settings and Google Earth, they distinguished the visited sites and labeled them as water and nonwater areas. We finally verified the ground truth data by comparison with false color composite Landsat 8 imagery, and rectified 19 labels. Explicitly, this step is to assess the true label y i as water or nonwater for each x i in the test set T. Finally, these labeled test sites {(x i , y i ,ŷ i )} n=304 i=1 , whereŷ i = fθ(x i ), were input as the ground truth information to generate a confusion matrix.

Confusion Matrix
The overall accuracy before test set rectification was 83.88% (Table 2) while the overall accuracy after test set rectification increased to 92.11% (Table 3). The diagonal elements (left to right, top to bottom) in the matrix represent the number of correctly classified pixels of each class, for example, the number of ground truth pixels with a label of water that was actually predicted as water during classification.  In contrast, the cross-diagonal elements represent misclassified pixels. A large loss of accuracy (40 out of 93) occurs at the pixels that are water in ground truth data but are classified into nonwater (Figure 2a). This happens for two reasons. First, many of these points are located at boundary pixels between two classes. Second, most points are isolated from their neighbor clusters. The boundary area has the mixed pixels problem which means both water and nonwater contribute to the observed spectral response of the pixel. In addition, the penetrating ability of C-band is unable to detect water hidden below rocks or vegetation cover where the normalized difference vegetation index (NDVI) is greater than 0.7 [47]. By contrast, there are only nine pixels that are misclassified into water while they are labeled as nonwater in ground truth data (Figure 2c). This is likely because the VV band can be affected by wind so that the wavy water surface will be classified into nonwater because of the diffuse refection. On the other hand, we examined each of those misclassified points and discovered that human errors in test data also compromise the overall accuracy. Twenty-two out of 40 nonwater test points and three out of nine water test points were rectified by using Landsat 8 imagery as a reference. As a result, the overall accuracy increased nearly 9%.
Note that the rectification of the test set T does not influence the training process nor the estimation of the optimal thresholdθ provided by the Otsu Method. Since the Otsu Method is an unsupervised learning algorithm, it does not depend on the label y of data for its training process. This property presents the feature of data corruption tolerance of the Otsu Method. In other words, corruption in the input data's label does not influence the actual training or performance of the Otsu Method. In addition, the test set T is only used for the statistical estimation or evaluation for the performance of the classification based on the estimation of optimal thresholdθ from the Otsu Method.
Furthermore, for factors influencing the accuracy, it is worth noting that the radar dB value with its corresponding label is linear nonseparable data [48]. In other words, there does not exist a hyperplane to clearly separate the dB value corresponding to the label of Water and the dB value corresponding to Nonwater. Because there exists two different regions or pixels i and j, where i = j, such that they have the same dB value x i = x j , but they are actually having different label y i = y j , one region corresponds to water and another region corresponds to nonwater. Such linear nonseparability may decrease the accuracy of this learning algorithm, which is unavoidable because the Otsu Method is trying to use a linear threshold to separate the data. One possible example of such region i and j is shown in the Figure 2b.

Water Area
The optimal threshold for 2020 was selected as −16.87 through the histogram, where the low peak corresponded to water pixels while the high peak corresponded to nonwater pixels (see Figure 3). Similarly, the optimal thresholds for each year from 2017 to 2019 were −14.88, −16.93 and −16.96 respectively. Based on the auto-selected thresholds derived from the Otsu Method, the imagery was classified into water and nonwater regions. As a result, the surface water acreages of Poyang lake from 2017 to 2020 were obtained and are presented in Table 4 and visualized from left to right in the Figure 4. The surface water area decreased by nearly 650 km 2 between 2017 and 2018, then increased by nearly 640 km 2 between 2018 and 2019 and finally decreased by nearly 856 km 2 between 2019 and 2020. This shows that surface water area of Poyang lake decreased with fluctuation, which is consistent with other research on variation of the surface water of Poyang lake during the time period of 1988-2016 [49].
In addition to the significant interannual variation, our results also showed the spatial variation of surface water area. The dry or draw down areas mainly occurred in the center and the boundary of the lake at the same time. The water areas located to the north (connected to Yangtze River) and west (connected to Gan River) accounted for most of the variation, while the water areas located at the east and south maintained much less variation.
The water area variations typically are closely associated to water level variations in Poyang lake basin [50]. Both water area and water level are dominant factors for wetlands in Poyang lake, and thus affect habitat distribution and accessibility. In this case, our classification results that show the spatiotemporal water area variations can provide robust linkages to habitat availability and suggest future research to further quantify this relationship.  , the x-axis represents that backscatter coefficient is calculated in dB scale, the y-axis represents that how many pixels have the same dB value in a bin and the interval of bins is 0.5.

Conclusions
Through this research we mapped the spatio-temporal variation of Poyang lake in January from 2017 to 2020 and showed that the surface water area fluctuated annually and tended to shrink both in the center and boundary of the lake over the past four years. The variation was consistent with related Poyang lake research for earlier decades. Our mapping approach involved a novel implementation of the Otsu Method and processing of Sentinel-1 data in Google Earth Engine. GEE performed well as a powerful cloud computing platform to implement an exhaustive searching algorithm. We provided detailed mathematical explanation to enumerate the advantages and limitations of the Otsu Method that were not clearly indicated in previous remote sensing research. We also demonstrated that the Otsu Method can be an effective classifier for threshold auto-selected algorithms to extract water surface with use of Sentinel-1 data. As a result, the Otsu Method has potential to be applied to other water related studies, such as water extraction applications for other lake regions, water pollutant detection for environmental assessment and aquatic habitat mapping for ecological conservation, using the open access scripts of the threshold algorithm contributed here.
In the future, to reduce the influence of linear nonseparability nature of the data, the 2D Otsu Method [51] can be applied. However, since the Otsu Method is an unsupervised method, we have not compared its performance with supervised learning algorithms. The supervised learning algorithm requires training data that is unavailable for the January of 2020 at this time. In addition, the Otsu Method is affected by the penetration ability of single C-band radar signal so that it is difficult to capture water beneath the vegetation. To advance our research toward mapping aquatic habitat availability, we recommend the comparison between supervised and unsupervised methods by using different series of imagery to discriminate vegetation zones. This next step will allow us to identify and project the spatial distribution of available foraging habitat under varying hydrological conditions for species of concern like Siberian crane and other aquatic organisms.
Funding: This research was funded by the Incubator program of Institute for Regional and International Studies at University of Wisconsin-Madison, grant number AAC3592.