Study on the Variations in Water Storage in Lake Qinghai Based on Multi-Source Satellite Data

Abstract

Performing research on the variation in lake water on the Qinghai–Tibet Plateau (QTP) can give the area’s ecological environmental preservation a scientific foundation. In this paper, we first created a high-precision dataset of lake water level variation every 10 days, from July 2002 to December 2022, using multi-source altimetry satellite SGDR data (Envisat RA-2, SARAL, Jason-1/2, and Sentinel-3A/3B SRAL), which integrated the methods of atmospheric path delay correction, waveform re-tracking, outlier detection, position reduction using a height difference model, and inter-satellite deviation adjustment. Then, using Landsat-5 Thematic Mapper, Landsat-7 Enhanced Thematic Mapper, and Landsat 8 Operational Land Imager data, an averaged area series of Lake Qinghai (LQ) from September to November, each year from 2002 to 2019, was produced. The functional connection between the water level and the area was determined by fitting the water level–area series data, and the lake area time series, of LQ. Using the high-precision lake water level series, the fitted lake surface area time series, and the water storage variation equation, the water storage variation time series of LQ was thus calculated every 10 days, from July 2002 to December 2022. When the hydrological gauge data from the Xiashe station and data from the worldwide inland lake water level database are used as references, the standard deviations of the LQ water level time series are 0.0676 m and 0.1201 m, respectively. The results show that the water storage of LQ increases by 11.022 × 10⁹ m³ from July 2002 to December 2022, with a growth rate of 5.3766 × 10⁸ m³/a. The growth rate from January 2005 to January 2015 is 4.4850 × 10⁸ m³/a, and from January 2015 to December 2022, the growth rate is 8.9206 × 10⁸ m³/a. Therefore, the increased rate of water storage in LQ over the last 8 years has been substantially higher than in the previous 10 years.

1. Introduction

The Qinghai–Tibet Plateau (QTP), as the third pole [1,2,3] of the earth, is crucial for understanding global climate change and geodynamics. The unique and delicate ecology that has developed on the QTP as a result of its specific natural circumstances is sensitive and vulnerable to climate change. Lakes are one of the most significant surface water types occurring on the QTP; they are very vulnerable to climate change and serve as important carriers of information about both worldwide climate change and regional responses to it [4]. Changes in lake water levels are an evident signal of climate dynamics [5,6,7,8,9]. Understanding the variations in lake water levels, and their area and water storage, is of great significance for hydrological and climatic studies. Lake Qinghai (LQ), located northeast of the QTP, is China’s biggest saltwater lake [10]. Its special location makes LQ particularly sensitive to climate change [11], and it is also an excellent location for studying climate response, water balance, and hydrological processes [12]. Monitoring the variations in LQ’s area, water level, and water storage, and understanding its changing patterns, can provide a scientific basis for ecological environmental conservation in the Lake Qinghai basin (LQB).

Remote sensing satellite technology has advanced considerably in recent years. With their advantages of a wide observation range and high efficiency, altimeter satellites, optical remote sensing satellites, and high-resolution synthetic aperture radar (SAR) imaging satellites have gradually become important means of keeping track of the variations in lake water levels and water storage for long-term monitoring, especially in remote areas, and are commonly utilized in the study of lakes on the QTP, e.g., [13,14,15,16,17]. The Ice, Cloud, and Land Elevation Satellite (ICESat) mission is the world’s first spaceborne laser altimetry system, with its main payload being the Geoscience Laser Altimetry System (GLAS), led by Texas State University and developed in collaboration with NASA and other scientific research institutions. The lake level data obtained by the GLAS on the ICESat mission can be precise to the decimeter level [13,14]. The ICESat satellite also has a short footprint radius (70 m [18]), which allows it to gather more important information and is why it is frequently employed in plateau lake water level monitoring [19,20,21,22]. Traditional nadir pulse-limited radar altimeters are not suitable for surveying lakes with small lake surfaces due to their large footprint; however, due to the large size of LQ, using traditional nadir pulse-limited altimeters for observation allows for the retention of more accurate observation results. In 2008, Jiang et al. [23] used Envisat altimeter data from 2002 to 2006 and an EGM96 quasi-geoid model, combined with waveform re-tracking technology to correct ranging and Gaussian filtering and remove altimetry noise, and studied water level changes in LQ. Other new sorts of altimeter satellite, such as SAR or SAR Interferometry (SARIn, e.g., CryoSat-2) altimeters, can also acquire useful data for monitoring lake level variations on the QTP [15,16,24]. Zhao et al. [25] 2007 studied the water level changes in LQ from 2002 to 2015 based on Cryosat-2/SIRAL (2010–2015) waveform re-tracking data and Envisat/RA-2 data. However, the short temporary coverage of single-mission altimeter satellites in lake monitoring impedes the study of various trends of water level over a long timescale. Certain satellites suffer from a long repetition time, such as ICESat/ICESat-2, which has a revisit duration of 91 days, whereas CryoSat-2 has a revisit period of 369 days, resulting in sparse temporary sampling. Because of this, they are ineffective in detecting periodic lake level changes [26]. Time variation series of lake water levels with higher sampling time resolutions can be established using multi-mission altimeter satellites. However, due to the large area of LQ, there are some problems with data fusion for different mission satellites, such as the disunity of the reference ellipsoid data, the disunity of geophysical correction models, and the large deviation in the transit trajectory position. Moreover, the QTP has unique climatic circumstances and complicated terrain, the typical waveform re-tracking algorithm does not function well in the plateau area, and the altimeter data in the lake area are fewer, and their quality unstable.

The water level–area formula can be used to determine changes in lake water storage, e.g., [27]. Moreover, changes in lake water storage can also be estimated using a water storage balance model, e.g., [28]. Both water level and area data must be accurate for using the water level–area relationship to estimate water storage variation, so that the inversion results can be reliable. Weather elements such as cloud cover and rainfall are major influencing factors in lake area extraction. In 2021, Hao and Luo [29] verified the accuracy of the Landsat-8 Operational Land Imager (OLI) using SAR images with higher spatial resolution. Based on the synthetic images with the least cloud cover in the high water period from June to September in 2017 and 2018, the normalized difference water index (NDWI) method was employed in the extraction of water bodies from Landsat-8 OLI and Sentinel-2 images in the same period. When estimating the area of LQ in this study, the two sub-lakes are stated individually, and only the main lake area is utilized to compute the water storage, making the computation more rigorous. In 2020, Du et al. [28] estimated the water storage change datasets for LQ from 1956 to 2017 by establishing a water balance model. In their study, rainfall was estimated using a Tyson polygon and the measured data from four weather stations. The surface runoff into the lake was calculated using the diameter flow of the Buha River. Groundwater recharge was directly regarded as a constant because it is difficult to observe. Because there are few hydrometeorological stations in the LQB, the accuracy of the obtained parameters (precipitation, temperature, etc.) cannot be guaranteed. It is difficult to obtain data on groundwater and lake evaporation involved in the water balance model. Although the assessment of water storage is inaccurate, it can also reflect the general trend of water storage change in LQ to a certain extent, which has a certain degree of reference value.

The lakes’ water levels on the QTP have been altered significantly over the last two decades because of worldwide climate change [30,31]. The water volume in LQ is also rapidly rising. To study the climate reaction in the LQ region, precise variations in water storage in LQ over the last twenty years must be determined. Using multi-source altimetry satellite (Envisat RA-2, SARAL, Jason-1/2, and Sentinel-3A/3B SRAL) data, this paper constructs a high-precision water level variation time series for LQ from July 2002 to December 2022, which integrates atmospheric path delay correction, a waveform re-tracking algorithm, outlier detection, position reduction using a height difference model, and inter-satellite deviation adjustment; additionally, it evaluates accuracy using hydrological gauge data from the Xiashe station and data from the international inland lake water level database. Based on the optical remote sensing images from Landsat-5 Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper (ETM), and Landsat 8 OLI, a high-precision averaged area series for LQ from September to November each year from 2002 to 2019 was constructed. By fitting the water level–area series, the functional relationship between the water level and the area was established. According to the retrieved lake water level change time series, the water storage change time series in LQ was calculated every 10 days, from July 2002 to December 2022, and analyzed accordingly.

2. Materials and Methods

As shown in Figure 1, Lake Qinghai (99°36′–100°46′E, 36°32′–37°15′N) is located on the QTP’s northeastern edge. The lake has a surface size of 4321 km² (2010), an average depth of 21 m (its greatest depth is 27 m), and a drainage area of 29,661 km² [32]. LQ originated from a stratigraphic fault depression, and its strata development is relatively complete.

The LQB is a relatively closed inland basin, with a spindle-shaped boundary outline and a northwest–southeast overall trend. The northwest terrain is high and the southeast terrain is low. The LQB has the temperate continental semi-arid climate of the QTP, with short summers and long winters, a large temperature difference between day and night, abundant sunshine, and strong winds. The lake effect and topography both influence the weather in the basin. The temperature in the lake region is greater than in the nearby mountains, and the average annual temperature in the southeast is higher than that in the northwest. The yearly mean air temperature (1951 to 2007) is 1.2 °C [33], the overall average annual temperature ranges from −1.1 to 4.0 °C, the greatest monthly average temperature is 11 °C, and the lowest is −13.5 °C. The typical annual precipitation varies from 291 mm to 579 mm. The water system in the LQB is unevenly distributed. There are many rivers in the west and north. The rivers with large flows include the Buha River [34], Shaliu River, and Hallgeir River. The runoff of these three rivers exceeds 75% of the runoff flowing into LQ. There are only a few tributaries in the east and south, such as the Daotang River and Heima River, which are mostly seasonal rivers, and their runoff is much smaller than that of rivers in the northwest.

2.1. Materials

2.1.1. Hydrological Gauge Data

The hydrological gauge station of Xiashe is situated to the southeast of LQ, with geographical positions of 100°29′24″E and 36°35′09″N. The hydrological gauge data were provided by Qinghai Provincial Water Conservancy Bureau (website is http://www.qhsl.gov.cn/ (accessed on 19 August 2021)). The gauge data observation time was from 2002 to 2019. The data format included the daily averaged water level, the monthly averaged water level, and the annually averaged water level. The data elevation data included the 1985 national height data of China, and settlement correction and data unification were carried out. The average daily water level data of LQ were utilized in this study to assess the accuracy of the lake water level time series produced using multi-mission altimeter satellite data.

2.1.2. Altimetry Satellite Data

The Level-2 20 Hz sensor geophysical data record (SGDR) from Envisat RA-2 (May 2002–December 2011), Level-2 20 Hz SGDR from off Jason-1 (February 2009–June 2013) and Jason-2 (October 2016–September 2019), Level-2 40 Hz SGDR data from SARAL (March 2013–December 2022), and SR-2 LAN 20 Hz SGDR altimetry satellite data from Sentinel-3A SRAL (April 2018–23 January 2023) and Sentinel-3B SRAL (June 2019–18 January 2023) were used to study the water level changes in LQ.

The Envisat RA-2 (RA2_MWS_2P) level-2 product and the Sentinel-3A/B SRAL altimeter data are distributed by the European Space Agency (https://earth.esa.int/ (accessed on 1 February 2023)). The Envisat satellite, successfully launched on 1 March 2002, had a repetition period of 35 days; it passes 094, 479, and 552 in phase B (14 May 2002–22 October 2010) and 238 and 567 in phase C (26 October 2010–8 April 2012), passing through LQ. A total of 142 trajectory data from Envisat RA-2 were used. The Sentinel-3A/3B track parameters (height, inclination, and coverage) were almost identical to those of Envisat RA-2. The Sentinel-3A/3B satellites had a repetition period of 27 days. In the study, 234 trajectory data from Sentinel-3A (with pass number: 067,332) and Sentinel-3B (with pass number: 136,332) that passed through LQ in this period were used.

The Jason-1/2 altimeter data were published by AVISO (Archiving, Validation, and Interpretation of Satellite Oceanographic data: https://www.aviso.altimetry.fr/ (accessed on 3 February 2021)). For Jason-1, which had a repetition period of 10 days, pass 129 passed through LQ in stage B (10 February 2009–3 March 2012). Those passes passed through LQ in stage C (7 May 2012–21 June 2013) were not fixed. A total of 129 trajectory data from the Jason-1 mission were used. Pass 129 of the Jason-2 mission, which had a repetition period of 10 days, passed through LQ in stage B (13 October 2016–17 May 2017), and the passes in stage C (11 July 2017–18 July 2018) and stage D (25 July 2018–1 October 2019) that passed through LQ were not fixed. A total of 74 trajectory data from the Jason-2 mission were used.

SARAL (Satellite with Argos and AltiKa) is a mission that was launched on 25 February 2013 by the Centre National d’Etudes Spatiales (CNES) and the Indian Space Research Organization. The SARAL satellite uses a Ka-band radar altimeter, which effectively reduces the influence of ionospheric changes on measurement. In the first 35 cycles, passes 94, 479, and 552 passed through LQ with a repetition period of 35 days. SARAL started to change to a drift orbit with an inclination of 98° in July 2016, and then, the passes passing through LQ were not fixed. A total of 323 tracks passing through LQ were obtained. The periodic water level variations in a lake can be obtained using data from the exact repeat mission (ERM) orbit stage, and the time resolution of lake level variations can be enhanced through the addition of data from the geodetic mission (GM) orbit stage, as well as more information about lake water level variations.

2.1.3. Landsat Data

Surface reflectance data from Landsat-5 TM, Landsat-7 ETM, and Landsat 8 OLI (http://www.gscloud.cn/ (accessed on 25 December 2021)) were utilized to investigate LQ’s water body area fluctuations. The Landsat data had a geographic resolution of 30 m and a revisit duration of 16 days. The Landsat-5 TM and Landsat-7 ETM data needed to be spliced by two scenes (strip number 133 and row numbers 34 and 35, respectively) to cover LQ; the Landsat 8 OLI data can be a scene image, with strip number 133 and line number 34. The area of LQ varies with variation in the water level, showing obvious seasonal variation characteristics. Moreover, the area variation is stable from September to November every year, so remote sensing images from September to November were selected to extract the lake boundary. To ensure image quality, the remote sensing images with the least cloud cover over LQ and in the period from September to November were chosen as far as was practicable. The Landsat data usage details are displayed in

Table 1. Landsat data usage details.

Sensor	Time	Cloud Cover (%)
Landsat-5 TM	11 October 2002	12.93, 10.25
	12 September 2003	13.54, 13.39
	20 September 2006	--, --
	28 September 2009	0.41, 1.3
Landsat-7 ETM	8 October 2004	0.06, 18.17
	9 September 2005	0.81, --
	15 September 2007	3.05, 33.19
	4 November 2008	0.59, 49.35
	7 September 2010	0.81, 26.23
	28 October 2011	9.43, 2.45
Landsat 8 OLI	9 October 2013	0.12
	13 November 2014	1.55
	15 October 2015	1.35
	17 October 2016	1.44
	4 October 2017	3.93
	21 September 2018	1.67
	11 November 2019	2.8

.

2.1.4. European Center for Mesoscale Weather Forecasting Data

The data on surface pressure and the vertical integral of water vapor that were released by the European Center for Mesoscale Weather Forecasting (ECMWF: https://www.ecmwf.int/ (accessed on 3 February 2023)) were used to recalculate dry tropospheric delay correction (DTC) and wet tropospheric delay correction (WTC), respectively. The latest version of the ERA-5 global reanalysis product released by ECMWF contains 0.25° × 0.25° grid data with a 1 h temporal resolution.

2.1.5. Lake Level Data from Other Databases

DAHITI (Database for Hydrological Time Series of Inland Waters, https://dahiti.dgfi.tum.de/ (accessed on 14 January 2023)) was established by the Technical University of Munich in 2013 [35]. The water level variation data for LQ were estimated using re-tracked data from Envisat (v3), SARAL (SGDR-F), Sentinel-3A (NTC), and Sentinel-3A (STC); the data cover LQ water level fluctuations from July 2002 to December 2022, and lack information from October 2010 to March 2013.

HTRWLD (high-temporal-resolution water level data: https://data.tpdc.ac.cn/ (accessed on 14 January 2023)) v2.0 is a lake water level, area, and water storage variation dataset (1976–2020) published by the China National QTP Scientific Data Center. The area data for lake variations on the Tibetan Plateau during 1976–2020 (>1 km²) [36,37] were used in the study. The datasets were provided by the National Tibetan Platform Data Center (http://data.tpdc.ac.cn (accessed on 14 January 2023)).

Hydroweb is a data access center for recording the variations in water level, area, and water storage of large rivers and lakes in the world. The data download website is https://hydroweb.theia-land.fr/hydroweb/ (accessed on 14 January 2023), and the format is text format [38]. The website, which was created by LEGOS/GOHS in France, provides information on the water levels, areas, and volumes of nearly 150 inland lakes. Water level and area are time series formed via the fusion of data from multi-mission remote sensing satellites, and water storage is calculated based on water level and area.

2.2. Methods

The data processing of LQ water storage variation is shown in Figure 2. It is mainly divided into four parts: the study of water variation based on multi-mission altimeter data; the study of lake area capture based on Landsat data; water level–area relationship fitting; and water storage estimation based on water level and area series. These are introduced separately in the following:

2.2.1. The Study of Water Variation Based on Multi-Mission Altimeter Data

• Creating a mask file

Landsat images with less cloud cover and that covered LQ were selected in 2005. The water body index method was employed in the extraction of the water body, and a binary map was made. The water body value was 1 and the land value was 0. The binary image was corroded inward by 100 pixels (about 1.5 km) to avoid near-shore contamination of the nadir satellite point. Then, the raster data were vectorized, and the mask file of the water level at the nadir satellite point in LQ was extracted.

2.

Calculation of lake’s water level

According to the satellite altimetry theory [26], the calculation formula for lake water level is

$$\mathrm{L}=\mathrm{A}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}-\left(\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}+\Delta \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\right)-\mathrm{R}\mathrm{C}-\mathrm{g}\mathrm{e}\mathrm{o}\mathrm{i}\mathrm{d}$$

(1)

where $\mathrm{L}$ represents the lake’s water level, and $\mathrm{A}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}$ indicates the orbital height of a satellite. $\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}$ is the altimeter-measured range of the sea surface; $\mathrm{R}\mathrm{C}$ is the distance correction for waveform re-tracking; and $\mathrm{g}\mathrm{e}\mathrm{o}\mathrm{i}\mathrm{d}$ is the geoid’s height with respect to the reference ellipsoid, this paper used the EGM2008 model results of the SGDR products. $\Delta \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}$ represents instrumental errors and various environmental corrections. Environmental corrections include ionospheric delay corrections, DTC, WTC, solid tide corrections, and pole tide corrections. Compared with $\Delta \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}$ errors, ocean tide, inverted barometer correction, thermal expansion, and wind accumulation impact have less influence, so the corrections are disregarded [26].

3.

Atmospheric path delay correction

To obtain a consistent and high-precision geophysical correction model among different satellites, this paper recalculated the DTC of each altimeter satellite (up to September 2022) using the reanalysis numerical model of ECMWF (ERA-5) [26,39,40,41]. As shown in Figure 3, compared with the values of the product, it has better continuity; for example, Figure 3e,f show that reprocessing results in smoother results than those of the product data.

The tropospheric delay for GPS receivers in the zenith direction can reach 2.3 m, and at a low-elevation angle, can reach 20 m [42,43]; it is affected by pressure, temperature, and humidity along the signal transformation route [42]. The correction radiometer was contaminated by land in inland lake locations, the WTC showed discontinuity, and some satellites did not have correction results at certain times. Therefore, the WTC was recalculated using the ERA-5 model (up to September 2022), and the formula was:

$$∆R_{vap}\approx \frac{{\beta }_{vap}^{\prime }}{T_{eff}}{\int }_0^R{\rho }_{vap}\left(z\right)dz$$

(2)

where ${\rho }_{vap}\left(z\right)$ is the density of water vapor at z; R is the maximum height considered for integration [44]; and the integral value is derived from the ERA-5 model. In the mid-latitude region, ${\beta }_{vap}^{\prime }/T_{eff}$ = 6.19 (cm³ g⁻¹) is based on the empirical value. The reprocessed results are shown in Figure 4, which shows better continuity after reprocessing, such as in Figure 4b,e,f, etc.

4.

Waveform re-tracking

Waveform re-tracking is the key to obtaining high-precision water level data for LQ. The difference between the leading-edge gate and the fixed nominal tracking gate is multiplied by a coefficient to achieve waveform re-tracking. At present, there are many waveform re-tracking methods [45,46,47,48], with the empirical method being more appropriate for inland lakes than the functional method [41]. Wang et al. [41] examined the waveform re-tracking methods of ICE, Threshold-20, and Threshold-50, and demonstrated that the ICE re-tracking method is more dependable for data on the plateau lake surface, generates more effective measured values, and has greater accuracy [26]. Therefore, the ICE waveform re-tracking algorithm was used in this study, and the waveform re-tracking gate was interpolated using the position with an amplitude of 30%.

5.

Abnormal value elimination and averaged lake water level calculation

When computing the average lake water level for each pass, outliers were removed. A sliding median filter algorithm and 2σ were employed in this research to find outliers [49]. The steps were as follows:

(1) The sliding median absolute deviation (MAD) method was applied repeatedly to determine outliers for all periods until no more outliers were found; the median filter window was 6 months. That is, observations were discarded as outliers, which were higher than three times the local MAD from the local median during a 6-month-long time window. MAD is defined as [26]:

$$MAD=median\left(\right|{LL}_i-median\left({LL}_i\right)\left|\right)$$

(3)

where $i=1,2,\dots ,N$, and $N$ is the total number of observations in the window.

(2) The 2σ criterion was used repeatedly for the observation of each satellite in each cycle, that is, data from the same satellite with median errors exceeding 2 times in one cycle were determined as outliers and eliminated.

Most of the outliers in the lake level observations were eliminated in steps (1–2). However, there were still a few significant outliers in the Sentinel 3A/B SRAL data, which were eliminated manually. The outlier rejection rates for the Envisat RA-2, Jason-1, Jason-2, SARAL, Sentinel 3A SRAL, and Sentinel 3B SRAL satellites were 43.31%, 27.70%, 26.88%, 44.92%, 40.81%, and 51.36%, respectively, among which Jason-2 was the lowest. Taking SARAL as an example, after these two steps of outlier rejection, SARAL observations tended to be stable, as shown in Figure 5. After the processing of steps (1–2), the averaged standard deviation (STD) was 0.0716 m, and 95.88% of the remaining 97 cycles of SARAL satellite data had a standard deviation of less than 0.15 m. There were 387 mean effective water level observation records in each cycle. After outlier rejections, the mean lake level was estimated per cycle to construct the series, as shown in Figure 6a. Figure 6b–g represent histogram distributions of the standard deviation of each cycle of the lake water level measured by each satellite. Figure 6b–g show that the expectations of μ and the standard deviation σ distributions of all satellites were stable, except for the bipolar distribution of the standard deviation of Sentinel 3A SRAL. Following the standard deviation distribution of Sentinel 3A SRAL data, the gross error in the Sentinel 3A SRAL data was further eliminated by taking other satellites as references. According to the statistics, the STDs of the Envisat RA-2, Jason-1, Jason-2, Sentinel 3A SRAL, and Sentinel 3B SRAL satellites were 0.1129 m, 0.0767 m, 0.1607 m, 0.4136 m, and 0.2391 m, respectively.

6.

Lake orthometric height difference model and position reduction

Because of the wide size of LQ, it is essential to incorporate an orthometric height difference model of the lake water level for position reduction. Therefore, we used the trajectory of each transit in the ERM stage to establish an orthometric height difference model of the lake surface. The steps for establishing the model were as follows:

(1) ERM trajectories with good quality were selected (for which the number of points was more than 30, the STD of the lake level was less than 0.2 m, and the height difference between two ends of the trajectory in the lake was less than 0.5 m), and the relative height difference at each trajectory point was calculated.

(2) Each trajectory of ERMs of the same satellite was averaged to obtain the orthometric height difference model value of the lake where this trajectory was located.

(3) The orthometric height difference model was calculated by decreasing offsets at the intersection sites of distinct trajectories.

Because there are few points in each ERM trajectory of Jason series satellites, this paper used Envisat RA-2 and Sentinel-3A/B SRAL repeated trajectories to establish the orthometric height difference model of lake water level, as shown in Figure 7. We selected the trajectory marked as “Pass No. 479” in Figure 7 to verify the accuracy, and used the repeated trajectory of the SARAL satellite for 1–35 periods to verify the accuracy using the height difference established in the same steps. The result, shown in Figure 8, shows that the STD of the difference between the two satellites passing through the same place was 0.034 m.

(4) The series of lake water level heights at each transit was interpolated to the middle position of the trajectory of the lake surface according to the height difference model, and the average value was taken to determine the lake water level at each pass.

The statistics show the number of different satellites that crossed LQ, and the numbers of Envisat RA-2, Jason-1, Jason-2, Sentinel 3A/B SRAL, and SARAL satellites were 138, 58, 38, 101, 85, and 310, respectively.

(5) We selected 5 points randomly from the center of the lake, the east of the lake, the west of the lake, the south of the lake, and the north of the lake; they were (100.2015, 36.8298), (100.4724, 36.7359), (99.9060, 36.9530), (100.0620, 36.7006), and (100.2754, 37.0880). The distances between these 5 points and the Xiashe station were about 37.38 km, 16.763 km, 66.137 km, 40.268 km, and 59.024 km, respectively. Because a small quantity of data were placed outside the model region, position reduction could be performed, and their own heights were used.

7.

Deviation adjustment

Due to the existence of objective influencing factors, such as instrument and equipment errors, in satellites, the inconsistent reference ellipsoid data, orbit inclination, and revisit period, there are systematic deviations among altimeters in the process of multi-mission altimeter data fusion, and different altimeters data cannot be connected accurately. There are still some systematic deviations in the water level time series after eliminating gross errors and outliers [50]. Multi-mission satellite data generate a continuous time series through deviation correction, which is important for the construction of lake water levels in continuous time series [51]. The steps for deviation adjustment in this paper were as follows:

(1) Calculate the time series of the deviation of the common period observations of two satellites:

$${∆}_j={LL}_j^A-{LL}_j^B$$

(4)

where ${∆}_j$ is the deviation between two satellites at time j; ${LL}_j^A$ indicates the lake water level as observed by satellite A at time j; and ${LL}_j^B$ represents the lake water level as observed by satellite B at the same time j. Due to the existence of system deviation, ${∆}_j$ is not 0. Calculate the deviation value using observations of two satellites during a common period. The common period is the maximum time interval between the observations of the satellite.

(2) Eliminate results greater than 2 times the median error.

(3) Calculate the mean value.

In this paper, the overlapping measurement period between two adjacent altimeters was calculated. The overlapping observation period of the Envisat RA2 and Jason-1 satellites was from February 2009 to December 2011. The overlapping observation period of the Jason-1 and SARAL satellites was from March 2013 to June 2013. The overlapping observation period of the SARAL and Jason-2 satellites was from October 2016 to September 2019. The overlapping observation period of the SARAL and S3A satellites was from April 2018 to December 2022. The overlapping observation period of the SARAL and S3B satellites was from June 2019 to December 2022. The overlapping observation period of the S3A and S3B satellites was from June 2019 to January 2023. Using the overlapping observation period data from two satellites, the deviation between satellites was estimated.

Table 2. Statistics of deviations between different satellites (number of used points/deviations (m)).

	Env/Jas1	Jas1/SARAL	SARAL/Jas2	SARAL/S3A	SARAL/S3B	S3A/S3B
Middle	34/0.4457	3/0.0744	50/−1.1519	79/0.0105	71/0.4266	72/0.5577
East	34/0.4532	3/0.0532	50/−1.1347	79/0.0250	73/0.4552	72/0.5569
West	34/0.4418	3/0.0855	50/−1.1608	79/0.0030	71/0.4191	72/0.5582
North	34/0.4463	3/0.0727	50/−1.1505	79/0.0117	71/0.4278	72/0.5577
South	34/0.4499	3/0.0908	50/−1.1651	79/−6.5 × 10−4	71/0.4156	72/0.5584
Maximum time interval (days)	10	10	10	5	5	5

shows the statistical results. The maximum time interval in Table 2 refers to the observations that can be considered common period observations of different satellites within this time interval, and it can be calculated via substitution into Equation (4). When the number of common period observations was considerable, we set the maximum time intervals to 5 days, and when the number of common period observations was limited, we set the maximum time intervals to 10 days. We adjusted the deviations between Jason-1 and Envisat RA2, SARAL and Jason-1, SARAL and Jason-2, SARAL and Sentinel 3A, and SARAL and Sentinel 3B, using the results in Table 2. The difference between Sentinel 3A SRAL and Sentinel 3B SRAL was used for verification. Under normal circumstances, the sum of deviations between SARAL/Sentinel 3A and SARAL/Sentinel 3B should be the same as the deviation between Sentinel 3A and Sentinel 3B, but there were many abnormal values of water levels for Sentinel 3A/B in Qinghai, resulting in a deviation of 0.138 cm between them.

8.

Gaussian filtering

When the deviations were adjusted, Gaussian filtering was employed to smooth the data and eliminate the noise signal. The Gaussian low-pass filter’s response function was:

$$h\left(\psi \right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\psi }^2}{2{\sigma }^2}}$$

(5)

where $\psi$ is the interval between points, and $\sigma$ is the window of the Gaussian filter. All water levels were convolved with the response function in the filter window $\sigma$ to obtain the filtered water level heights. The larger the size of $\sigma$, the better the Gaussian filtering smoothness. However, if $\sigma$ is too large, it may eliminate the required non-noise signal and reduce the data resolution. In this paper, half a year was used as the filtering window for Gaussian filtering.

2.2.2. The Study of Lake Area Capture Based on Landsat Data

• Landsat image data preprocessing

Because of the effects of spatio-temporal, spectral, and radiative resolution, remote sensing image data contains faults in the data-collecting process, resulting in picture quality degradation. Therefore, before image analysis and processing, it is necessary to preprocess the original image. The preprocessing process mainly includes the following four steps: strip restoration, radiometric calibration, atmospheric correction, and image mosaic

(1) Strip repair

The Landsat 7 satellite scanner failed in 2003, and data were lost in 7 bands of Landsat 7. Landsat 7 ETM images were repaired using the method provided by the website https://www.ixxin.cn/2016/12/03/landsat-7qutiaodaiyuchuli/ (accessed on 25 December 2021).

(2) Radiometric calibration

A satellite transforms electromagnetic wave signals into electrical signals through sensors, stores the information in the original remote sensing image, and then, quantifies it into a dimensionless gray value DN. Such results lack physical significance. To obtain the physical quantities needed for related research (such as reflectivity, temperature, radiance, etc.), radiometric calibration is needed. Using the calibration tool in ENVI 5.3, the image’s gray value is converted into emissivity, and the conversion formula was as follows:

$$\mathrm{L}= \mathrm{Gain}\ast \mathrm{D}\mathrm{N}+\mathrm{O}\mathrm{f}\mathrm{f}\mathrm{s}\mathrm{e}\mathrm{t}$$

(6)

In the formula, $\mathrm{L}$ represents emissivity, $\mathrm{Gain}$ is the increment of emissivity, and $\mathrm{O}\mathrm{f}\mathrm{f}\mathrm{s}\mathrm{e}\mathrm{t}$ is the offset of emissivity. The calibration steps are as follows:

① Set the calibration type: radiance (radiance value);

② Perform storage sequence selection: BLI (stored in band sequence);

③ Perform output type selection: floating-point type;

④ Set the zoom factor: click “Apply FLAASH Setting”;

(3) Apply atmospheric corrections.

When the electromagnetic signals reflected by different ground objects pass through the atmosphere and clouds, they will decay to different degrees, and the gray value in the image will change, so it is impossible to accurately judge ground objects. Therefore, before water extraction, it is necessary to correct the atmosphere of the original image after radiometric calibration. In this process, FLAASH, a tool provided by ENVI 5.3, is used for atmospheric correction. The basic process is as follows:

① Input the radiometric calibration image;

② Set the sensor and image information;

③ Select the parameters of the atmospheric model, the Mid-Latitude Summer model;

④ Select the aerosol model (choose the countryside here);

⑤ Choose 2-Band (K-T) aerosol retrieval;

⑥ In multispectral settings, select Over-Land Retrieval Standard (660:2100) as default.

(4) Create the image mosaic

During image processing, it was found that a single image of Landsat 5 TM or Landsat 7 ETM could not cover LQ, and images with 133/34 and 133/35 serial numbers needed to be spliced to completely cover LQ. In this paper, the image was mosaic, created with the help of Seamless Mosaic (seamless mosaic), which comes with ENVI 5.3. The operation is as follows:

① Overlay the two images in sequence, set “Data Ignore Value” to 0, and set the automatically generated edge line;

② Based on the image with a large area of LQ, select histogram matching to uniformly process the two images;

③ Connect the feathers of the edges to make them naturally blend into the new image.

2.

Extraction of lake surface area

Because only the water surface of LQ needed to be extracted, the ISO clustering unsupervised classification method was adopted, and the number of clusters was set to 6. The unsupervised classified water body was divided into two categories (LQ and other ground objects), and the classification effect is shown in Figure 9. This classification method does not need prior information on the study area, and only needs to set the number of classifications. Because of the relatively large scope of LQ, it was easy to find the target water body. Then, the raster image was vectorized and the main lake area connected with LQ was derived. Finally, the vector boundary was visually interpreted in combination with the original image.

2.2.3. Water Level–Area Relationship Fitting

Considering the cold winter, cloud cover, and other factors, the extraction accuracy for the LQ area was affected. Landsat images from September to November of each year between 2002 and 2019 were utilized in the extraction of the water area of LQ, and the lake’s water level–area relationship was fitted using averaged water level and area data from September to November of each year. To make the computation easier, it was assumed that the water level and area satisfied the following fitting relationship:

$$S=\mathrm{a}{h}^2+\mathrm{b}h+\mathrm{c}$$

(7)

where $S$ is the average lake area; $h$ is the average water level from September to November every year; and $\mathrm{a}$, $\mathrm{b}$, and $\mathrm{c}$ are the parameters. It is worth noting that $h$ is the averaged data of water level from September to November every year in the fitting relationship. Therefore, the average water level in Equation (7) indicates the average area $S$ of the lake from September to November in a year.

2.2.4. Water Storage Estimation Based on Water Level and Area Series

The water level data provided by the hydrological station or multi-source altimeter satellite, and the water area extracted from the Landsat image, were utilized to compute the water volume change in LQ. The area of LQ extends with rising water levels and contracts with falling water levels. Therefore, the lake area $S\left(h\right)$ is an expression that monotonously increases with the change in water level. The expression of lake water quantity change was obtained using the definite integral principle [52]:

$$\Delta V={\int }_{h_0}^{h_1} S(h)dh$$

(8)

where $\Delta V$ indicates the variation in lake water volume from the time of water level $h_0$ to the time of water level $h_1$, and $h_0$ and $h_1$, respectively, indicate the water level of LQ at times $t_0$ and $t_1$. Considering the maintenance and failure of instruments, the gap period between satellites, and other factors, it is impossible to realize continuous observations of water level and area. According to Equation (7), which is used to calculate the water level–area fitting relationship, a lake area with missing time can be fitted using the water level observation value. It can be seen from the remote sensing images of LQ that the lake surface is almost oval. According to previous study findings [53], the area of rising water level in LQ every month (the maximum fluctuation is about 2 cm), and even every year (the maximum fluctuation is about 80 cm), can be ignored relative to the lake area. Combining Equation (7) with Equation (8), the formula of water storage change in LQ can be further expressed as [27,54]:

$$\Delta V={\int }_{h_0}^{h_1} \left(a{h}^2+bh+c\right)dh=\frac{1}{3}\left(h_1-h_0\right)\left[S_0+S_1+a{h}_0{h}_1+\frac{1}{2}b\left(h_0+h_1\right)+c\right]$$

(9)

where $\Delta V$ is the variation in water volume when water levels vary from $h_0$ to $h_1$, assuming that the initial water reserve is 0, and the water reserve at time t is the cumulative variation in water volume from the initial time to that time.

3. Results and Analysis

3.1. Time Series of Lake Water Level in LQ

Continuous water level change findings in LQ were acquired in the east, west, south, north, and center points of the lake using the approach proposed in Section 2.2.1. and multi-mission satellites. After Gaussian filtering processing with the time window set to half a year, the results were resampled at intervals of 10 days. Figure 10a shows the resampled results, and Figure 10b shows an expanded version of the resampled results in Figure 10a’s box.

To validate the inversion of LQ water level data in this research, the findings were compared to the measured data from the Xiashe station and the data on LQ in the DAHITI database. For the convenience of comparison, the data on LQ in the DAHITI database were additionally treated via Gaussian filtering and resampling using the same approach. We used data collected within two days as a common period to compare and generate statistics.

Table 3. Comparison of water level products (m).

Products	Max	Min	Mean	STD	RMS	Distance (km)
East in situ gauge (639 points)	0.2356	−0.1772	−0.0057	0.0696	0.0698	16.763
West in situ gauge	0.3219	−0.0918	0.0778	0.0673	0.1028	66.137
North in situ gauge	0.2849	−0.1256	0.0447	0.0678	0.0811	59.024
South in situ gauge	0.3372	−0.0778	0.0915	0.0673	0.1136	40.268
Middle in situ gauge	0.2898	−0.1211	0.0491	0.0677	0.0835	37.38
Mean value 1 in situ gauge	0.2925	−0.1187	0.0515	0.0676	0.0849
DAHITI data—east (145 points)	0.4154	−0.3347	0.0258	0.1242	0.1264
DAHITI data—west	0.3186	−0.3170	−0.0541	0.1178	0.1292
DAHITI data—south	0.3027	−0.3330	−0.0680	0.1176	0.1356
DAHITI data—north	0.3569	−0.3868	−0.0229	0.1216	0.1233
DAHITI data—north	0.3569	−0.3868	−0.0229	0.1216	0.1233
DAHITI data—middle	0.3518	−0.3918	−0.0273	0.1215	0.1241
DAHITI data—mean value 1	0.3491	−0.3526	−0.0293	0.1201	0.1232

¹ The final results of this study.

shows statistics on the difference between the results of this paper and the measured data from the Xiashe station, the difference between DAHITI data and the results of this paper, and the difference between DAHITI data and the measured data from the Xiashe station. Table 3 demonstrates that compared with the water level of the Xiashe station, the STDs of the calculated results of this paper range from 0.0673 m to 0.0696 m. The size of the root means squares (RMS) varies with the location distribution, and the minimum is 0.0698 m, located at the eastern point. The two largest RMS values are 0.1136 and 0.1028 m at the southern and western points, respectively, and the RMS of the other two points at the northern and middle points are 0.0811 and 0.0835 m, respectively. The number of points compared with the DAHITI data is 145. The STD range of the difference between them is 0.1176~0.1242 m, and the maximum RMS is 0.1242 m, obtained at the southern point in Qinghai. The minimum value of RMS is 0.1233 m, calculated at the northern point. At the same time, this paper also generated STD and RMS statistics on the difference between the measured data from the Xiashe station and those from DAHITI, which are 0.1264 m and 0.1305 m, respectively. The average value of five computed time series of LQ water level variation was used as the final result in this article. The STDs of the discrepancies between the final result and the measured data from the Xiashe station and DAHITI are 0.0676 m and 0.1201 m, respectively, as shown in Table 3.

Figure 11 shows a comparison between the calculation results of each satellite in the middle position and the measured water level data of the hydrological gauge station. Figure 12 shows a comparison between the final results of this paper and the hydrological gauge data from Xiashe and DAHITI. Figure 12 shows that the interannual fluctuation features of LQ water levels are visible regardless of the observed water level data from the hydrological gauge station or satellite altimetry data, and there will be “peak” or “trough” characteristics every other year, with good continuity. The water levels determined by the hydrological gauge station coincide better with the water levels retrieved by altimeters in the rising stage of water levels than in the falling stage of water levels every year, which may be related to the annual freezing of LQ. From the research of Li Xiaofeng et al. [55], it can be seen that LQ starts to freeze in middle and late December of each year and completely freezes in early January of the following year. At the end of March, the temperature rises, and the lake ice gradually melted until it disappears completely in early April. Thus, the difference between the water level observation data from the Xiashe station and the lake water levels retrieved in this paper was compared with the results of the northern hemisphere lake ice thickness dataset [56] (provided by the National Tibetan Platform Data Center (http://data.tpdc.ac.cn (accessed on 14 January 2023)), as shown in Figure 13. Figure 13 shows a relationship between residual error and ice thickness.

3.2. Analysis of Water Level Variation in LQ

It can be seen from Figure 12 that from July 2002 to December 2022, the overall water level of LQ displayed a trend of first decreasing, and then, rising, during which the overall water level of LQ increased by 3.78 m. Before March 2005, the water level of LQ displayed a downward trend, reaching its lowest value in March 2005; additionally, its lowest value decreased by 0.236 m compared with July 2002, with a growth rate of −8.85 cm/a.

There was a “sudden change” between June and October in 2005, and the maximum water level difference was as high as 48 cm. From August 2005 to June, 2009, the lake water level was in a state of equilibrium and in fluctuation, with a slight upward trend, and the change in the lake water level was only about 16 cm over more than three years. From June 2009 to August 2012, the water level increased from 3193.41 m to 3194.39 m, with a total increase of 0.98 m and an annual growth rate of 30.5 cm/a. From August 2012 to August 2016, the water level was in a relatively balanced state, and the water level rose from 3194.39 m to 3194.67 m, with an annual growth rate of 7 cm/a. From August 2016 to December 2022, the water level of LQ was in a state of rapid growth, rising from 3194.67 m to 3196.82 m, with a growth rate of 33.5 cm/a.

3.3. Lake Surface Area Extraction and Fitting

Using the ArcGIS platform, lake boundary extraction was performed on Landsat images with fewer clouds from September to November for each year from 2002 to 2019. Only the area of the main LQ was calculated, excluding the two small lakes (Gahai and Erhai), which have been separated from the main lake area for a long time. To confirm the extraction accuracy of the Landsat images, the results from 2014 and 2017 were compared with the dataset integrated by Hao and Luo [29]. The results show that the coincidence rate between them is over 99%. In addition, the extracted area of the Landsat image was compared to the model data of LQ provided by the Legos official website. Because the model data include sub-lakes separated from LQ, such as Gahai and Erhai, the calculated average area from September to November is larger than that determined in the experimental results of this paper, as seen in Figure 14.

The area series of LQ extracted from the Landsat image is consistent with the changing trend of Hydroweb data, with a correlation of 0.9855 and an RMS of 55.608 km². With time, the absolute value of the difference between them is becoming smaller and smaller, which shows that the increasing speed of the main lake area is gradually speeding up or the percentage of the sub-lake area is gradually slowing down. All of the validation results suggest that Landsat images perform well in extracting water bodies from the Qinghai Lake area; hence, the multi-year area data obtained using Landsat images can be considered reliable.

Since elements such as cold winter and cloud cover have a substantial influence on the accuracy of LQ area extraction, Landsat images from September to November in the years ranging from 2002 to 2019 were utilized in this work. The relationship between the lake water levels of the Xiashe station and the area series was fitted. The fitting relationship is as follows:

$$S=0.43362806h^2-2678.20906561h+4134769.27$$

(10)

It should be noted that the fitting relationship between the water level and the area is based on the average data from September to November, and the average water level of lake $h$ in Formula (10) indicates the average area of lake $S$ in that year. For the determination coefficient, $R^2=0.95 (p<0.01)$. The area fitted by Formula (10) was compared to the water area determined using remote sensing images to validate the fitting effect. The fitted area is similar to the Landsat water area, and the maximum and minimum percentage differences are 0.28% and 0.0044%, respectively, indicating that the fitted result of water level and area is good, as seen in Figure 15.

3.4. Analysis of LQ Area Variation

Based on Landsat images from September to November every year from 2002 to 2019, this paper extracted the water area of LQ using ISO clustering unsupervised classification. The average area of LQ in these 18 years is 4296.36 km². As shown in Figure 14, from 2002 to 2019, the lake’s area grows by 285.18 km², with a growth rate of 15.8 km²/a. From 2002 to 2019, the changing trend of lakes is consistent with the water level, showing an overall trend of decreasing, and then, increasing. From 2002 to 2004, the lake area decreases continuously, with a growth rate of −7.35 km²/a, and the lake area decreases by 14.71 km² in two years. From September to October 2004, the lake area is the smallest at 4190.46 km², which is 105.9 km² lower than the average area over the 18 years, and is consistent with the lowest water level in March 2005. In 2019, the lake area is the largest at 4490.35 km², which is 193.99 km² higher than the average area, for many years. During this period, the water level in LQ was rising at a high speed.

4. Discussion

The final result obtained from satellite altimetry was substituted into Equation (9) to obtain the water storage variation time series, and the water storage time series growth in LQ from May 2002 to December 2022 was obtained via cumulative calculation, as shown in Figure 16. From July 2002 to December 2022, the water volume of LQ increases by 11.022 × 10⁹ m³, with an annual growth rate of 5.3766 × 10⁸ m³. The overall trend increases in fluctuation. From July 2002 to January 2005, the water storage of LQ displays a decreasing trend, with a total decrease of 5.9949 × 10⁸ m³. After this period, the water storage of LQ increases in fluctuation, and by January 2010, the water storage increases by 1.5566 × 10⁹ m³ compared with July 2002. From January 2005 to January 2010, the water storage increases by 2.1561 × 10⁹ m³, so the increased speed from January 2005 to January 2010 is 4.312 × 10⁸ m³/a. The rising rate is 4.6578 × 10⁸ m³/a from January 2010 to January 2015, which is comparable with the conclusion estimated by Zhang et al. (2018): the rising rate of LQ from 2001 to 2016 was 4.5 × 10⁸ m³/a [27]. From January 2005 to January 2015, the rising rate is 4.4850 × 10⁸ m³/a. From August 2016, the water level shows a rapid growth trend. Figure 16 shows that the rising rate from January 2015 to December 2022 is 8.9206 × 10⁸ m³/a. The increased rate of lake water storage in LQ in the last 20 years tends to increase.

The influencing factors of LQ water volume variation include lake runoff, lake surface precipitation, lake surface evaporation, and other factors (such as human activities, groundwater, changes in lake bottom topography, etc.). Due to a lack of data on lake bottom topography, the variation in lake water volume analyzed in this work is the relative variation in lake storage compared to July 2002, and the influence of lake bottom topography fluctuations was not evaluated. This study used the water volume variation model from the literature, which aids in computation but also introduces some errors. This is because the variation in water volume is connected to the topography of the freshly flooded region, as well as the rise and fall of water levels. Yet, because LQ has a tight elliptical shape and we used high-frequency sampling time intervals, the error of water volume variation is within an acceptable range. The following work can be combined with the terrain data of the region to conduct appropriate analyses.

The contribution of glacier to water volume variation in the LQB is extremely tiny due to the limited number of glacial areas, and the influence of human activities on the change in water volume in the LQ is, likewise, quite minimal [57]. Therefore, meteorological and hydrological factors such as precipitation, lake runoff, and evaporation directly govern the variation in LQ water volume [58]. Precipitation and inflow runoff are the LQ’s recharge sources, whereas water surface evaporation is the LQ’s source of water volume expenditure. The monthly cumulative precipitation and annual average precipitation changes in LQ from 2002 to 2019 were estimated using GPM (Global Precision Measurement, https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGDF_06 (accessed on 25 December 2021)) precipitation data, as illustrated in Figure 17. The mean annual precipitation on the lake surface is 342.38 mm, with a growth rate of 3.69 mm/a (the red line in Figure 17b). The average annual precipitation on the lake surface in 2003 is the lowest, and this occurs one year earlier than the lowest point of the LQ water level in 2004, indicating that the response of the water level to precipitation has a certain lag; this is consistent with the research results of Li et al. [59]. There was minimal variation in temperature rise over the research period. Temperature rises typically lead to an increase in evaporation capacity. However, because evaporation is affected by other hydrological and meteorological factors (such as radiation, wind speed, and underlying surface parameters), an increase in the average maximum temperature indirectly inhibits an increase in annual evaporation via the thermal cycle, and the effect of increased evaporation capacity caused by a temperature rise may be masked; this can result in a negative relationship between temperature and evaporation capacity [60]. In the future, the contribution of various influencing variables to variation in lake water volume can be quantitatively analyzed in combination with the lake water balance model.

5. Conclusions

Using multi-mission altimeter data (Envisat RA-2, Jason-1/2, Sentinel 3A/B SRAL, and SARAL SGDR data), this paper presents a method of integrating atmospheric path delay correction, a waveform re-tracking algorithm, outlier detection, position reduction based on a height difference model, and inter-satellite deviation adjustment to construct a high-precision lake water level time variation series over 20.5 years, from July 2002 to December 2022, in LQ. The accuracy of the proposed method was validated using hydrological gauge data from the Xiashe station, as well as data from the International Lake Database. The STDs are 0.0676 m and 0.1201 m, respectively, based on water gauge data and International Lake Database statistics. Therefore, this approach has good monitoring accuracy and may be expanded to monitor water storage variations in lakes that are comparable to LQ but lack tidal data.

Reliable water volume variation data represent an important basic connection in local ecological environment studies. Combined with the optical remote sensing images from Landsat-5 TM, Landsat-7 ETM, and Landsat 8 OLI, and a formula for water volume change, the water storage change time series in LQ was estimated every 10 days, from July 2002 to December 2022. The results show that: (1) The water level of LQ has a dynamic upward trend from July 2002 to December 2022, with a total increase of 3.78 m; from July 2002 to March 2005, the growth rate of the water level is −8.85 cm/a; from August 2005 to June 2009, the lake water level is in a state of equilibrium in fluctuation; from June 2009 to August 2012, the annual growth rate is 30.5 cm/a; from August 2012 to August 2016, the water level is in a relatively balanced state; and from August 2016 to December 2022, the water level of LQ is in a state of rapid growth, with a growth rate of 33.5 cm/a. (2) From July 2002 to December 2022, the water storage of LQ increases by 11.022 × 10⁹ m³, with an annual growth rate of 5.3766 × 10⁸ m³; the water storage of LQ displays a downward trend from July 2002 to January 2005; from January 2005 to January 2015, the rising rate is 4.4850 × 10⁸ m³/a; from August 2016, the water level shown a rapid growth trend; and the rising rate is 8.9206 × 10⁸ m³/a from January 2015 to December 2022. The growth rate of water storage in LQ has been growing during the last 20 years, and it has grown dramatically over the last 8 years as compared to the preceding 10 years.

This paper only studied the dynamic variations in water level, area, and water storage in LQ over the last 20 years from the perspective of acquired data, without further analysis of water storage balance and climate factors. Moreover, it does not reveal the drivers of lake water storage variation or the proportion of various factors’ involvement in the change process. In this paper, the precise variation estimation of water storage in LQ can offer reliable data information for water balance study in this basin, and provide basic data for the sustainable growth of water resources and ecological environmental protection in this area.