Long-Term and Bimonthly Estimation of Lake Water Extent Using Google Earth Engine and Landsat Data

Abstract

Lakes are one of the most important parts of the terrestrial hydrosphere. The long-term series of lake area dynamic data with high spatial-temporal resolution is of great significance to the study of global change of the water environment. Satellite observations (such as Landsat) have provided images since the 1970s, but there were challenges for the construction of long-term sequences of lake area on a monthly temporal scale. We proposed a temporal-spatial interpolation and rule-based (TSIRB) approach on the Google Earth Engine, which aims to achieve automatic water extraction and bimonthly sequence construction of lake area. There are three main steps of this method which include bimonthly image sequence construction, automatic water extraction, and anomaly rectification. We applied the TSIRB method to five typical lakes (covering salt lakes, river lagoons, and plateau alpine lakes), and constructed the bimonthly surface water dataset (BSWD) from 1987 to 2020. The accuracy assessment that was based on a confusion matrix and random sampling showed that the average overall accuracy (OA) of water extraction was 96.6%, and the average Kappa was 0.90. The BSWD sequence was compared with the lake water level observation data, and the results show that the BSWD data is closely correlated with the water level observation sequence, with correlation coefficient greater than 0.87. The BSWD improves the hollows in the global surface water (GSW) monthly data and has advantages in the temporal continuity of surface water data. The BSWD can provide a 30-m-scale and bimonthly series of surface water for more than 30 years, which shows good value for the long-term dynamic monitoring of lakes, especially in areas that are lacking in situ surveying data.

1. Introduction

As an important part of the terrestrial hydrosphere, lakes not only provide the water resources that are needed by human activities and the ecological environment, but also play an indispensable role in the hydrological cycle [1,2]. Lakes have the function of maintaining the stability of the regional ecological environment, such as agricultural production, fishery resources, flood prevention, and disaster reduction [3]. At the same time, lakes are very sensitive to human activities and climate change [4]. In recent decades, with the increase of human activities and the intensified global climate change, the expansion and shrinkage of lakes have occurred from time to time around the world [5,6]. This change has produced a series of ecological and social problems, such as the deterioration of lake water quality, the degradation of ecosystems, the reduction of biodiversity, and the inundation of arable land, which have severely affected the livelihoods of local people [7]. Therefore, long-term lake observation data with high temporal and spatial resolution is of great significance for revealing the influence of natural factors and human activities on the sustainable development, utilization, and protection of lake waters.

Remote sensing technology enables earth observation from space. The development of remote sensing technology has advanced human understanding of the Earth surface to a new stage, and at the same time has brought convenience to the study of large-scale lake water dynamic monitoring and parameter inversion [8]. Periodically-acquired and remotely sensed imagery enable us to monitor the dynamics of the surface environment in a timely and rapid manner. Remote sensing and GIS technology are widely used in water resource monitoring applications, including flood disaster/loss assessment and management [9,10], surface water dynamic monitoring [11,12], water quality assessment and monitoring [13], and water-related infection diseases [14] etc.

Surface water dynamic monitoring can utilize sensors with different temporal, spatial, and spectral resolutions, such as MODIS, AVHRR, Landsat TM\ETM+\OLI, Sentinel-2 MSI, etc. [15,16,17,18]. Among them, Landsat series satellites have become one of the most widely used optical sensors in surface water and other environmental monitoring research due to their long continuous monitoring time (greater than 30 years) and relatively high spatial resolution (30 m).

General methods for extracting water bodies from optical remote sensing images can be divided into four types: (a) thematic classification methods [19,20,21], (b) linear decomposition methods [22], (c) single-band\band index thresholding methods [23,24], and (d) rule-based methods [25]. The band index threshold methods were widely used in automatic water extraction applications due to their ease of use and low computational time and resource consumption. McFeeters [26] was the first to propose the Normalized Difference Water Index (NDWI) for TM images using the green and near-infrared bands based on the construction principle of the Normalized Difference Vegetation Index (NDVI). Xu proposed modified Normalized Difference Water Index (mNDWI) [27], which improved the separation between water body and built-up area; Feyisa et al. proposed a multi-band index AWEI to improve the separation among water body, shadows and other dark targets [24]. Li et al., proposed BDWI to improve mNDWI by bringing local background information [28]. A combination of different band indexes and thresholding rules were also proposed to improve accuracies [29].

Remote sensing information processing and analysis requires a lot of computing and storage resources, but the traditional workstation-based operation mode can only process a small amount of data. The rise of cloud computing platforms such as the Google Earth Engine (GEE) makes it possible to process massive remote sensing images rapidly. Many scholars have carried out automatic extraction of water bodies at global and regional scales based on the GEE platform [30,31,32,33,34]. One of the most influential is the global surface water (GSW) dataset that was constructed by Pekel et al. [33]. It is worth noting that most of existing research has focused on monitoring the multi-year changes of surface water, while the time series mapping of surface water with high temporal and spatial resolution still faces challenges. The GSW dataset provides monthly and 30 m resolution surface water data, but there are a lot of voids in the data, which limit the direct use of this data for time series analysis.

The main objectives of this paper were: (1) proposed a new method of TSIRB (temporal-spatial interpolation and rule-based approach) for surface water mapping based on the GEE platform and Landsat data, which aims to realize long-term bimonthly time series mapping of lake water at the 30-m scale; (2) selected five typical lake cases in different regions of the world, extracted the bimonthly surface water dataset (BSWD) in the past 34 years, and analyzed the annual and inter-annual variation characteristics of the case lakes; (3) evaluated the accuracy of the TSIRB method for water extraction, and analyzed the statistical correlation characteristics of the obtained BSWD data with in situ and remotely sensed water level data; and (4) discussed the advantages and shortcomings of the TSIRB method for producing long-term bimonthly lake extent and area sequences.

2. Materials and Methods

2.1. Case Lakes

A total of five test lakes from around the world, which feature different climates, water types, and landscapes, were chosen to examine the robustness and general applicability of the TSIRB method. The five case lakes (Figure 1 and

Table 1. Summary of descriptions and characteristics of the case lakes.

Case Lakes	Geo-Locations	Area (km2) *	Type	Characteristics	Climates (Köppen Key) [47]
Chiquita Lake	30°~30°55′S, 62°~63°W	4610	Saline lake	Salt lake in South America; vegetation abundant in the north; precipitation in December to next March.	Humid subtropical climate (Cfa)
Urmia Lake	37°46′′N, 45°20′E	3273	Saline lake	High salinity; low water depth; salt crust around lake shore.	Cold semi-arid climate (Bsk)
Dongting Lake	28°30′~30°20′N, 111°40′~113°10′E	1542	Freshwater lake, river lagoon	The second largest freshwater lake in China; a lot of vegetation around; precipitation concentrated in July to September.	Humid subtropical climate (Cfa)
Poyang Lake	28°22′~29°45′N, 115°47′~116°45′E	2879	Freshwater lake, river lagoon	The largest freshwater lake in China; a lot of vegetation around; precipitation concentrated in July to September.	Humid subtropical climate (Cfa)
Qinghai Lake	36°32′~37°15′N, 99°36′~100°16′E	4504	Saline lake, plateau lake	The largest saline lake in China; precipitation concentrated in June to September; freezes from December to March	Monsoon-influenced subarctic climate (Dwc)

* Average area of lakes in 2020.

) include (a) Chiquita Lake, (b) Urmia Lake, (c) Dongting Lake, (d) Poyang Lake, and (e) Qinghai Lake. These lakes cover several typical lake types: freshwater lake, saline lake, river lagoon, and plateau lake (detailed in Table 1).

Chiquita Lake (30°~30°55′S, 62°~63°W) is located about 150 km northeast of Córdoba, Argentina [35], in the Paraná Basin of South America. The lake is affected by a humid subtropical climate, and precipitation is mainly concentrated from December to March (summer in the southern hemisphere) [36]. There is abundant vegetation around the lake [37].

Urmia Lake (37°46′N, 45°20′E) is an endorheic salt lake in Iran. At its greatest extent, it was the largest lake in the Middle East and the sixth-largest saltwater lake on Earth [38]. The lake is a shallow lake (with an average depth of 1.5 m in 2019) [39], and the lake area varies greatly among different seasons. The lake is affected by a cold semi-arid climate, with an average annual precipitation of 341 mm, an average annual temperature of 11.2 °C, and an average annual cumulative evaporation of 1200 mm. The temperature in the lake area varies greatly between seasons, with the lowest temperature reaching −23 °C in winter and the highest being 39 °C in summer [40].

Poyang Lake (28°22′~29°45′N, 115°47′~116°45′E) and Dongting Lake (28°30′~30°20′N, 111°40′~113°10′E) are located in the middle and lower reaches of the Yangtze River Basin in China [41,42], which are the first and second largest freshwater lakes in China. The lake belongs to Humid subtropical climate [43,44], and the rainfall in the lake area varies significantly seasonally. The flood season is from July to September, during which the lake surface area increases rapidly [45,46]. Aquatic vegetation is widely distributed around the lake.

Qinghai Lake (36°32′~37°15′N, 99°36′~100°16′E), located in the northeastern part of the Qinghai-Tibet Plateau, is the largest saline lake in China [48]. The lake belongs to the monsoon-influenced subarctic climate and is an inner-flow lake that mainly relies on precipitation to supplement its water volume. The average annual precipitation is 395 mm, and the precipitation is concentrated from June to September [49]. The yearly average temperature is below 0 °C. The freezing period begins around November every year, a stable ice surface begins to form in December, and begins to thaw in March or April of the following year [50].

2.2. Platform and Data

2.2.1. GEE Platform

The open Google Earth Engine (GEE) platform was used as the computing environment in this study. GEE is a cloud-based planetary-level geospatial analysis platform. It provides massive remote sensing image resources and programming interfaces, which greatly facilitates remote sensing method research and application at global and regional scales. GEE relies on Google’s powerful cloud computing capabilities to accelerate the processing of remote sensing data [51], which is beneficial for large-scale mapping and long-term sequence analysis based on multisource and high dimensional imagery [33,52,53,54].

2.2.2. Landsat Imagery

We used the multispectral imagery of Landsat5 TM, Landsat7 ETM+, and Landsat8 OLI in this study. The Landsat5 satellite was launched in 1984, but there were no Landsat5 images before 1987 in many parts of the world; we chose 1987 as the starting year for lake monitoring around the world. In addition, due to the contamination of clouds and shadows, the proportion of effective pixels is low in images with large cloud cover, then images with cloud cover greater than 50% were filtered out in this study.

The Landsat Surface Reflectance (SR) data [55] was used for water extraction. Landsat Surface Reflectance products include Quality Assessment (QA) bands which identify pixels that exhibit adverse instrument, atmospheric, or surficial conditions. In this study, the QA bands were used to screen out pixels that might be affected by instrument artifacts or subject to cloud and cloud shadow contamination.

2.2.3. The Lake Water Max-Extent Mask

The max-extent refers to the largest extent of surface water bodies in the lake area over the past 30 years. In this study, the max-extent data was used to define the outer boundary of the lake water body, which could reduce the amount of data processing in non-lake areas. In addition, the use of max-extent data can reduce the confusion of water body extraction that is caused by shadows and urban dark targets and improve the accuracy of water body extraction during non-flood seasons.

In this study, we used the maximum water extent product (max_extent) from the Global Surface Water (GSW) dataset that was provided by the European Commission’s Joint Research Centre (JRC).

2.2.4. Lake Water Level Data

Lake water level series data were used to validate the results of lake area sequences because of the strong correlations between the lake water level and the lake area. In the five case lakes, the water level data came from two different sources. The water level data of Poyang Lake, Dongting Lake, and Qinghai Lake were obtained through site measurements. The data were extracted from the “Hydrology Yearbook of People’s Republic of China”, the time series was 2006–2018, and the time resolution was days. Among them, the water level data for Poyang Lake, Dongting Lake, and Qinghai Lake were measured in Duchang Station, Chenglingji Station, and Xiashe Station, respectively [56,57].

Water level data of the other two lakes were obtained from the Global Lake and Reservoir Monitoring website (G-REALM) (https://ipad.fas.usda.gov/cropexplorer/global_reservoir, last accessed on 30 August 2021), ranging from 1992 to 2020 with monthly temporal resolution.

2.3. Methods

This paper proposes a temporal-spatial interpolation and rule-based (TSIRB) method for automatic water extraction, based on the GEE platform and utilizing Landsat image data. The method consists of three main steps (Figure 2): bimonthly image sequence construction (BISC), automatic water extraction (AWE), and time series anomaly detection and reconstruction (TSADR). Using this method, the problem of void holes in Landsat’s monthly surface water data (such as GSW) can be alleviated, and the automatic construction of bimonthly surface water dataset (BSWD) can be realized.

2.3.1. Bimonthly Image Sequence Construction (BISC)

Optical remote sensing images are always contaminated by clouds, which leads to a large number of void holes in water extraction products. The problem of void holes must be solved in order to obtain the time series of surface water which is continuous in space and time. First, this paper sets the time granularity of the time series to be bimonthly, which can reduce the proportion of void pixels compared with the monthly series. Through cloud and shadow masking, we combined all the valid pixels within two months into one image.

However, in some areas, the heavy cloud coverage in the rainy months makes most of the pixels invalid, which leaves a large number of void holes in the bimonthly composite image (Figure 3b). In this study, we assumed that at the same pixel location and the land cover types are similar for the same month in adjacent years. Then, the current void holes can be replaced by valid pixels of the same month in adjacent years, so as to realize the construction of a bimonthly-scale valid image sequence that is continuous in time and space.

Applying this idea, this paper proposes a method to construct a spatiotemporally continuous bimonthly image sequence. The pixels in each bimonthly composite image are cloud- and shadow-free. The method mainly consists of two main steps.

The first step is the construction of preliminary bimonthly synthetic sequence imagery. This includes cloud and shadow masking, and a bi-monthly composite of valid data. Cloud and cloud shadow masks use QA indicators from Landsat’s surface reflectance data, and areas that are masked by clouds and shadows are defined as invalid. All the masked images of the corresponding months were combined into an image set, and then the median synthesis algorithm was used to obtain the initial bimonthly image imP_year_month. The imP_year_month is a representative of all the valid observations within the two-month timeframe. It should be noted that, at a pixel, if all the original images before compositing are invalid values, the composite result will still be invalid values. Due to the limitation of the relatively low observation frequency of Landsat and the influence of clouds and shadows, there are still many void values in the preliminary bimonthly composite images.

The second step is void filling of the bimonthly image sequences. Here, we assumed that if a pixel was water in the May 2010 imagery, then there is a high probability that it was also water in the May 2009 and May 2011 images. Based on this assumption, for invalid values in the preliminary image sequence, we used valid observations from adjacent years to replace the invalid values for the current cell. There are two ways to retrieve bimonthly valid values: forward current year and backward current year. In order to facilitate the processing of the start year and the end year, we used 2010 as the time boundary. Images after 2010 are processed by year-forward retrieval, and images in 2010 and before are processed by year-backward retrieval to obtain valid values and perform composition. After the void filling operation in the second step, we obtained a continuous bi-monthly image sequence that was composed of all valid observation pixels (Figure 3c).

2.3.2. Automatic Water Extraction (AWE)

After preparing the bimonthly composite images, the automatic extraction of water bodies is carried out using a multi-rule-based classification method. The spectral index makes the pixel of interest stand out from the image background, and the binary thresholding method is very suitable for large-scale and time-series image analysis due to its ease of implementation. Different water body indices (NDWI, mNDWI, AWEI) are widely used in surface water mapping. However, due to the influence of complex interference factors under different temporal and spatial conditions, the overall accuracy of the band index thresholding method is not high enough when it is applied to large-area water body extraction or time series mapping. These disturbance factors mainly include shadow, ice, snow, vegetation, and so on. In this study, a multi-rule-based classification method was used to achieve automatic water extraction (

Table 2. Spectral indexes and image characteristics that were used for water delineation from backgrounds using Landsat images.

Features	Equations	Thresholding	Objectives
AWEI_sh	blue + 2.5 × green − 1.5 × (nir + swir1) − 0.25 × swir2	Empirical: >−0.005	Water mapping
Brightness	(nir + red + swir1)/3	Empirical: varied	Ice and snow noise
mNDWI	(green − swir1)/(green + swir1)	mNDWI > NDVI or mNDWI > EVI	Vegetation noise
NDVI	(nir − red)/(nir + red)
EVI	2.5 × (nir − red)/(nir + 6 × red − 7.5 × blue + 1)
Max_extent	GSW Max_extent	Max_extent = 1	Suppress terrestrial noise in complex environments

), and a long-term time series of bimonthly surface water data was obtained.

(1) First, we apply AWEI_sh [24] for water body extraction and obtain a preliminary result of surface water extraction. AWEI_sh is an index that is proposed to improve the accuracy of water extraction in shadowed or other dark surface areas [24]. Studies have shown that AWEI_sh has good regional adaptability and can realize automatic extraction of water bodies in large areas [58]. Referring to the studies of others [24,59], combined with the spectral analysis results of this study, we adopt AWEI_sh > −0.005 as the decision rule for water bodies. In most cases where the images are clear, using the AWEI_sh decision rule can achieve good water extraction results (Figure 4).

(2) Water extraction in winter images is often disturbed by ice and snow. In this paper, a brightness index is used to assist in water extraction. The reflectivity of snow in each waveband is higher than that of water [60], but its spectral curve shape is similar to that of water, and both snow and water have small calculated AWEI_sh values. In the case of winter images with snow and ice, using AWEI_sh alone to extract water bodies will result in a large error in the results. Considering that snow has higher reflectance in red, near-infrared, and short-wave infrared than water, this study proposes a brightness index, which is the average of the reflectance values of red, near-infrared (nir), and the first short-wave infrared band (swir1). Brightness values that are greater than a certain threshold can be considered non-water bodies. This rule is set according to the time of the freezing season in the northern and southern hemispheres.

(3) In order to reduce the impact of surrounding vegetation noise on water extraction, this paper uses water index (mNDWI) and vegetation index (NDVI and EVI) to assist in the classification. The water index highlights water body information, while the vegetation index highlights vegetation information. If an image pixel’s water index value is lower than the vegetation index value, the pixel should be classified as vegetation rather than water. Studies by others [61,62] also support this hypothesis. In this study, only image pixels satisfying the condition ‘(mNDWI > NDVI) or (mNDWI > EVI)’ were classified as water bodies, and others were classified as non-water bodies.

(4) In order to suppress terrestrial noise in complex environments, this paper uses the maximum water body range data for masking. In winter, the surface water of lakes in high latitudes and plateau areas will freeze, and the automatic extraction of water bodies is greatly affected by ice and snow. In the case of snow cover on the lake surface and around the lake, the snow cover on the land part with lower brightness is often misclassified as a water body. To deal with this type of disturbance, in this study, we can reasonably assume that water bodies are less likely to appear in pixels outside the maximum water body during the snow-covered months. We selected the max_extent data of GSW as the mask to refine the results of water extraction in winter (Figure 5). In the flood season, the inundation range of the water body is wide, and we do not select the max_extent for masking, which can avoid the omission of water bodies in the flood season.

2.3.3. Accuracy Evaluation of Water Extraction

We used the confusion matrix method to evaluate the accuracy of automatic water extraction. The data series processed in this study spanned more than 30 years, resulting in a large amount of processed data. We selected the water body extraction data in 2010 and 2010 as the representative of all water extraction results to carry out accuracy evaluation. Verification points were obtained by random sampling. A total of 200 verification points were randomly selected from each bimonthly image for each lake, and a total of 12,000 verification points were sampled. The reference attributes of verification points were obtained by visual interpretation to ensure the accuracy of attribute acquisition.

There are several accuracy evaluation indicators for the confusion matrix method. We selected the most representative Kappa coefficient and overall accuracy (OA) for the accuracy evaluation. The Kappa coefficient represents the general accuracy of the method, and the OA represents the percentage of correctly classified pixels [63]. The OA and Kappa are calculated as the following formulas:

$$\Delta =\left(TP+FP\right)\times \left(TP+FN\right)+\left(TN+FN\right)\times \left(TN+FP\right)$$

(1)

$$kappa=\left(N\times \left(TN+TP\right) \right)/\left(N^2-\Delta \right)$$

(2)

$$OA=\left( TN+TP\right)/N$$

(3)

In the formulas: $TP$ is the number of the true positive points; $TN$ is the number of the true negative points; $FP$ is the number of the false positive points; $FN$ is the number of the false negative points; and $N$ is the total number of verification points.

2.3.4. Time Series Anomaly Detection and Reconstruction (TSADR)

Using the above AWE method, a bimonthly surface water image sequence can be obtained for each lake. By vectorizing the binarized water extraction images, the time series of the lake area data can be obtained. Although the AWE step can handle most of the disturbances, there may be a small amount of anomalous data in the lake area time series due to snow on the lake area in winter and bright reflection in the shallow water area of the salty lakes. This study adopted the time series processing method to detect and repair sequence anomalies, then the reconstructed bimonthly lake area sequence data was obtained.

1.

Anomaly Detection

For anomaly detection in lake area sequences, the moving average method (MA) was used. The MA method calculates the moving average by adding or subtracting old and new data sequentially, so as to eliminate accidental factors and explore the time development trend of things [64]. The MA method is often used for outlier detection of hydrological feature sequences [65,66]. In this method, the moving average ($y_i$) at each time point is calculated as the average of the value at that time point and the k observations on both sides [67]. Hydrological elements such as surface water will be affected by the values of the previous and subsequent periods, so we chose one year before and after as the sliding period (k is set to 6) and calculated $y_i$. The moving average difference ratio ($z_i$) was calculated at each time point and the standard deviation $\sigma$ of the $z_i$ sequence data was calculated; with $3\sigma$ as the threshold, when $z_i$ falls outside the threshold, the sequence value at this time point is judged to be outlier. Repeat the above steps until no outliers can be detected. Figure 6a shows the result of anomaly detection in lake area sequence of Qinghai lake.

$$y_i=\left(x_{i-k}+\dots +x_i+\dots +x_{i+k}\right)/\left(2k+1\right)$$

(4)

$$z_i=\left(y_i-x_i\right)/\left(y_i\right)$$

(5)

$$\overline{z}={\displaystyle \sum }_{i=1}^n{z}_i/n$$

(6)

$$\sigma =\sqrt{{\displaystyle \sum }_{i=1}^n{\left(z_i-\overline{z}\right)}^2/n}$$

(7)

In the equations: $x_i$ is the value of the original sequence at time i; $y_i$ is the moving average at time i; $z_i$ is the moving average difference ratio at time i; $\overline{z}$ is the average value of the moving average difference ratio sequence; and $\sigma$ is the standard deviation of the moving average difference ratio $z_i$ series.

2.

Time Series Reconstruction

For the reconstruction of the lake area sequences, the weighted moving average method (WMA) was used. Compared with other series data interpolation methods, the WMA method is simple and easy to use, and can quickly process non-stationary series data [68]. The lake water surface area series data generally fluctuates in a 12-month cycle. The lake area value of the current month is affected by the values of the previous and next months (intra-periodic influence), and by the value of the same month of the previous and the next year as well (periodic influence). When carrying out anomaly restoration of lake sequences, it is necessary to comprehensively consider intra-period and periodic impacts. For the intra-periodic effect, we choose 3 as the moving average window size and calculated the intra-periodic moving average (IMA). For periodic effects, we choose 2 as the window size and calculated the periodic moving average (PMA).

$$IM{A}_i=avergae\left(n{v}_{i-3}+n{v}_{i-2}+n{v}_{i-1}+n{v}_{i+1}+n{v}_{i+2}+n{v}_{i+3}\right)$$

(8)

$$PM{A}_i=avergae\left(n{v}_{i-12}+n{v}_{i-6}+n{v}_{i+6}+n{v}_{i+12}\right)$$

(9)

In the formulas: $IM{A}_i$ is the intra-period moving average at time i; $PM{A}_i$ is the periodic moving average at time i; and $n{v}_i$ is the sequence value after removing outliers at time i.

The lake surface area sequence generally exhibits obvious periodicity, and there are always minimum and maximum values in a full cycle. For the restoration of the max area month, the IMA value often underestimates the area of the current month; while for the restoration of the minimum value month, the IMA value often overestimates the area of the current month. In order to improve the accuracy of sequence repair, different formulas are used for extreme months and non-extreme months in this research. For each lake, we first calculate the multi-year monthly average of the lake area series to obtain the extreme month information. Then, for the abnormal point in the extreme month (Em), the repair value is calculated according to the Formula (7). For the abnormal points in non-extreme months (NEm), Formula (8) is used to calculate the repair value. Figure 6b shows the results of outlier repair and sequence reconstruction of the Qinghai Lake area sequence.

$$R{v}_{Em}=IM{A}_i/3+2\ast PM{A}_i/3$$

(10)

$$R{v}_{NEm}=IM{A}_i/2+PM{A}_i/2$$

(11)

In the formulas: $R{v}_{Em}$ represents the repaired value of the extreme month (at time i) and $R{v}_{NEm}$ represents the repaired values of the non-extreme month (at time i).

2.3.5. Lake Area Series Trend Analysis

In order to further analyze the variation characteristics of the lake area time series, the commonly used linear regression model and the Mann–Kendall test were employed in this research. The trend of lake areas of five case lakes from 1987 to 2020 were determined by linear regression while the Mann–Kendall test was conducted to detect whether the slope for the regression line is significant.

3. Results

3.1. Accuracy of Water Extraction

The accuracy evaluation results showed that the average value of the Kappa coefficient and OA were 0.90 and 96.6%, respectively, which indicated that the overall reliability of TSIRB method on water extraction was good. The specific results are shown in

Table 3. Results of Kappa coefficients for accuracy evaluation of surface water extractions by TSIRB.

Month	2010					2020					Average
	Chiquita	Urmia	Dongting	Poyang	Qinghai	Chiquita	Urmia	Dongting	Poyang	Qinghai
1	0.89	0.89	0.88	0.78	1.00	0.97	0.80	0.79	0.76	0.75	0.85
3	0.96	0.96	0.76	0.89	1.00	0.97	0.77	0.76	0.81	- *	0.87
5	0.93	0.93	0.81	0.92	0.96	1.00	0.93	0.82	0.92	0.93	0.92
7	0.97	0.97	0.96	0.97	0.97	0.93	0.93	0.88	0.84	1.00	0.94
9	0.93	0.93	0.96	0.86	0.93	0.90	1.00	0.92	0.83	1.00	0.93
11	1.00	1.00	0.69	0.88	0.89	0.97	0.86	0.65	1.00	1.00	0.89
Average	0.95	0.95	0.84	0.88	0.96	0.96	0.88	0.80	0.86	0.94	0.90

* Qinghai Lake was covered with a lot of snow in March 2020 and was detected as an abnormal point. The data of this month were not included in the accuracy evaluation results.

and

Table 4. Results of overall accuracies (OA) of surface water extractions by TSIRB.

Month	2010 (%)					2020 (%)					Average
	Chiquita	Urmia	Dongting	Poyang	Qinghai	Chiquita	Urmia	Dongting	Poyang	Qinghai
1	95.0	95.0	99.0	91.7	100.0	98.3	90.0	98.0	90.9	88.3	94.6
3	98.3	98.3	96.0	95.0	100.0	98.3	88.3	96.0	91.7	- *	95.8
5	96.7	96.7	96.0	96.5	98.3	100.0	96.7	98.0	96.7	96.7	97.2
7	98.3	98.3	99.0	98.3	98.3	96.7	96.7	97.0	93.3	100.0	97.6
9	96.7	96.7	99.0	93.3	96.7	95.0	100.0	98.0	91.7	100.0	96.7
11	100.0	100.0	96.0	96.7	95.0	98.3	93.3	96.0	100.0	100.0	97.5
Average	97.5	97.5	97.5	95.3	98.1	97.8	94.2	97.2	94.0	97.0	96.6

* Qinghai Lake was covered with a lot of snow in March 2020 and was detected as an abnormal point. The data of this month were not included in the accuracy evaluation results.

. The highest value of OA is 100% and the lowest value is 88.3%; the highest value of Kappa coefficient is 1.00 and the lowest value is 0.65.

There are certain differences in the accuracy of water extraction results for different types of lakes. The average Kappa coefficient of river lagoons (Poyang Lake and Dongting Lake) is lower; the Kappa coefficient of the salty lake and plateau alpine lake (Chiquita Lake, Urmia Lake and Qinghai Lake) is higher.

3.2. Variation Characteristics of the Case Lake Area Sequence

3.2.1. Characteristics of Interannual Variation

The Landsat imagery of five case lakes were processed by the TSIRB method, and 34-year bimonthly time series data of lake area were obtained. The sequence visualization of the five case lakes are shown in Figure 7. The interannual variation trends for the five case lakes were different. The area series of Lake Chiquita, Lake Urmia, Dongting Lake, and Poyang Lake all showed significant decreasing trends (p < 0.01), among which Urmia Lake decreased the fastest, with a change of 20.13 km²/bimonthly, followed by Chiquita Lake, with a change of 8.39 km²/bimonthly. The area reduction rates of Dongting Lake and Poyang Lake were relatively close, with changes of 1.55 and 3.91 km²/bimonthly, respectively. Qinghai Lake showed a different variation trend from the other four case lakes. Its area experienced significant expanding (p < 0.01) in the past 34 years with a relatively slow change rate of 0.64 km²/bimonthly.

3.2.2. Characteristics of Intra-Annual Variation

The long-term bimonthly time series of lake water area can also be used to analyze the seasonal variation characteristics of lakes. Figure 8 plots the peak, median, and valley values of the lake area for each month of the five case lakes. For Dongting Lake and Poyang Lake, the water area showed strong seasonal variation. The lake area gradually expanded from January to July and shrunk from July to November. The area peak in July was two to three times the area in January. Dongting Lake and Poyang Lake were the most representative lakes in the Yangtze River Basin in Asia, and their seasonal trends were similar.

Chiquita Lake, Urmia Lake, and Qinghai Lake also showed seasonal characteristics, although not as prominent as Dongting and Poyang Lake. The maximum value of the water area of Chiquita Lake appeared in May, and the annual change curve showed a single peak shape. The bimonthly variation of the area of Urmia Lake showed a bimodal shape and the peaks of the lake area appeared in May and November, respectively, and the lowest value appeared in September. For Qinghai Lake, the intra-annual variation was completely opposite to that of Urmia Lake, showing a shrinking trend from January to May, expanding to the largest from July to September, and then starting to shrink.

4. Discussion

In order to further discuss the effectiveness of the proposed TSIRB method and the obtained BSWD data, we compared the results of this study with those of previous studies and carried out the correlation analysis between the lake area series and the lake water level series of the case lakes. We also discussed the uncertainties of using the TSIRB method to produce lake area sequence in this section.

4.1. Comparison with the Results of Other Studies

4.1.1. Comparing Different Results on Single Case Lake

Through literature analysis, we collected the results of other studies [69,70,71,72,73,74,75] in four case lakes (Urmia, Dongting, Poyang, and Qinghai Lake), including eight study cases (see

Table 5. Summary of studies by other researchers for case lakes.

Lakes	Image Sources	Methods	Time Span	Temporal Resolution	Correlations	References
Urmia	Landsat	Unsupervised classification	1984–2011	year	0.99	Kabiri et al. [69]
	Landsat	NDWI	2000–2013	year	0.95	Rokni et al. [70]
Dongting	Landsat	AWEI	2002–2013	season	0.93	Li et al. [72]
	Sentinel-1A	thresholding segmentation	2016	month	0.97	Huth et al. [73]
Poyang	Landsat	AWEI	2002–2013	season	0.99	Li et al. [72]
	Sentinel-1A	SWI thresholding	2015–2016	month	0.99	Tian et al. [71]
Qinghai	Landsat	mNDWI and Classification	1987–2016	year	0.94	Tang et al. [74]
	Landsat	mNDWI	1999–2009	year	0.93	Zhu et al. [75]

for details). Statistical analysis showed that BSWD had a strong positive correlation with other research results. The lowest correlation coefficient was 0.93 and the highest value was 0.99. This proves the reliability of the results of this paper from one side.

Plotting the BSWD with area data of other studies on a single graph (Figure 9) makes it easy to compare data from different sources. The BSWD sequence showed similar fluctuation characteristics and changing trends with other research results. Over the comparable time period, BSWD lake areas were highly consistent with areas of most of the other studies. The result of Tang et al. on Qinghai Lake was generally higher than our results. The reason is that they included the range of Qinghai Lake and the two surrounding small lakes to calculate the area, but this research and Zhu et al. calculated the area of Qinghai Lake only.

In addition, it should be noted that BSWD provides complete bimonthly time series of lake area for 34 years with 30 m spatial resolution, the spatial and temporal resolution show good potential for long term lake monitoring applications.

4.1.2. Comparison with GSW MONTHLY Data

GSW is currently the only dataset that provides global long-term monthly surface water data at 30 m resolution [33]. For the five case lakes, the GSW data and the BSWD data were jointly plotted in Figure 10. As the GSW dataset did not handle invalid pixels in the water extraction step, and there were a lot of void values (lake area is 0) in GWS lake area sequences. Among them, Urmia Lake had the most void values (165), and the void value rate reached as high as 40.4%. Except for Chiquita Lake, the void rate of GSW monthly data for the other four case lakes was above 30%. However, the BSWD data was continuous in time and there was no void value month for each case lake.

In comparison with other research results [69,70,71,72,73,74,75], it was found that the long-term change trend of BSWD time series was relatively accurate. While the time series of the GSW monthly data was difficult to be used for trend analysis directly, for example for Qinghai Lake, the seasonal fluctuation of GSW monthly series is very large and does not match the observation data [74,75].

4.2. Correlation Analysis between BSWD Lake Area and Water Level Sequence

The lake surface area generally changes with the lake water level, and there is a certain correlation between them. We collected water level data of some time periods for case lakes and plotted the lake water level and BSWD area data together in one graph (Figure 11). The variation trend of the BSWD area of the case lakes are similar to the variation trend of the water level in the same period. For the five case lakes, the correlation coefficients between the BSWD area and the water level in the same period are all higher than 0.87, and the highest value is 0.95. They all passed the significance test at p < 0.01. These also confirm the reliability of the TSIRB method and the obtained BSWD data from one aspect.

4.3. Uncertainties and Shortcomings

Using the TSIRB method to extract water bodies, lush aquatic vegetation, bright shallow water, and surface snow may introduce uncertainties.

4.3.1. Influence of Aquatic Vegetation

Aquatic vegetation may affect the TSIRB method on water extraction. The TRISB method utilizes the rules of water index (mNDWI) and vegetation index (NDVI and EVI) to suppress the interference of vegetation on water extraction and can deal with the case of general sparse vegetation. However, during the peak vegetation season (summer and autumn), dense aquatic vegetation may appear on the water surface and boundary areas of some lakes. Analyzing images with strong vegetation interference, the extracted water body area would be smaller than the actual area.

4.3.2. Influence of Shallow Water with High Brightness

Shallow water with high brightness may affect the automatic extraction of water bodies. The TSIRB method extracts water bodies based on the rules of various band indices which are affected by the reflectivity of the ground objects themselves. In the boundary area of saline lakes, the water is shallow and clean which makes the spectral characteristics of the shallow water similar to that of the surrounding salt tidal flats. The brightness threshold that is constructed in TSIRB can suppress the influence of high-brightness shallow water to a certain extent, but it is still difficult to accurately distinguish between high-brightness shallow water and surrounding tidal flats. When applying the TSIRB method, the brightness threshold we set was relatively small, which generally results in the omission of high-brightness shallow water areas.

4.3.3. Influence of Lake Icing and Snow Cover

Lake icing and surface snow may introduce uncertainties to the TSIRB method. Although we use a brightness value metric to reduce the effect of ice and snow on water extraction, the generalized threshold setting is difficult to achieve high accuracy. For the case where the lake surface is covered with snow, it is difficult to determine the boundary of the lake water even by visual interpretation [76]. In addition, for time series processing, it is difficult to set a common threshold to apply to all the images. When applying the TSIRB method, the brightness threshold we set is relatively small, which may lead to the possibility of missing the ice surface under the snow cover, making the extraction result smaller than the actual area.

4.3.4. Uncertainty in Bimonthly Image Sequence Construction

For imagery that is influenced by heavy clouds, we replaced invalid holes with valid pixels in the same month from adjacent years in the TSIRB method. This is generally applicable in the central area of the lake. It should be noted that if invalid cells occur in the boundary area of the lake, it is possible that the land cover status of the pixel location varies between adjacent years due to climatic (precipitation/evaporation) differences. The classification of boundary pixels might be mistaken especially under extreme drought or extreme flood conditions.

4.3.5. Uncertainty in Time Series Anomaly Detection and Reconstruction

The uncertainty of the time series anomaly detection and reconstruction steps of TSIRB mainly comes from the extreme conditions of hydrological elements, such as extreme floods or extreme droughts. The bimonthly area data in extreme cases may appear outside the range of three standard deviations and thus be defined as outliers and removed. In this case, in the interpolation step based on the periodic rules, the interpolation value of TSIRB would be larger in the extreme drought cases and smaller in the extreme flood cases than the original actual value. In practice, however, such extreme situations are rare.

5. Conclusions

Long-term monitoring of lake water is an important part of global change research. However, obtaining long-term lake observation data with high temporal and spatial resolution still face challenges. In this study, we developed a novel surface water mapping method—TSIRB, to construct long-term and bimonthly surface water dataset (BSWD) at a 30-m scale.

Application tests were carried out on five typical lakes in different parts of the world, and 34-year BSWD products for the case lakes were obtained. The evaluation results of water extraction accuracy show that for different types of lakes in different seasons, the average values of Kappa coefficient and OA were 0.90 and 96.6%, respectively, which shows a relative high level of accuracy. The BSWD products show strong positive correlations with the lake water level sequences (correlation coefficient 0.87~0.95) and show good agreement with the results of other independent studies on single lakes (correlation coefficients 0.93–0.99). Compared with the GSW monthly water body products, the BSWD products improved the void values in the monthly sequence data, show stronger temporal continuity, and can obtain complete and valid lake water sequences.

The TSIRB method is implemented on the GEE platform and can be fully automated for water detection. This method has good potential to be extended to long-term lake monitoring applications in other regions of the world. BSWD is a long-term continuous lake water product, which is of great significance to the study of global lake water environments, especially in areas where ground observation data is lacking.