Improving the Inversion Accuracy of Terrestrial Water Storage Anomaly by Combining GNSS and LSTM Algorithm and Its Application in Mainland China

: Densely distributed Global Navigation Satellite System (GNSS) stations can invert the terrestrial water storage anomaly (TWSA) with high precision. However, the uneven distribution of GNSS stations greatly limits the application of TWSA inversion. The purpose of this study was to compensate for the spatial coverage of GNSS stations by simulating the vertical deformation in unobserved grids. First, a new deep learning weight loading inversion model (DWLIM) was constructed by combining the long short-term memory (LSTM) algorithm, inverse distance weight, and the crustal load model. DWLIM is beneﬁcial for improving the inversion accuracy of TWSA based on the GNSS vertical displacement. Second, the DWLIM-based and traditional GNSS-derived TWSA methods were utilized to derive TWSA over mainland China. Furthermore, the TWSA results were compared with the TWSA solutions of the Gravity Recovery and Climate Experiment (GRACE) and Global Land Data Assimilation System (GLDAS) model. The results indicate that the maximum Pearson’s correlation coefﬁcient ( PCC ), Nash–Sutcliffe efﬁciency ( NSE ) coefﬁcient, and root mean square error ( RMSE ) equal 0.81, 0.61, and 2.18 cm, respectively. The accuracy of DWLIM was higher than that of the traditional GNSS inversion method according to PCC , NSE , and RMSE , which were increased by 67.11, 128.15, and 22.75%. The inversion strategy of DWLIM can effectively improve the accuracy of TWSA inversion in regions with unevenly distributed GNSS stations. Third, this study investigated the variation characteristics of TWSA based on DWLIM in 10 river basins over mainland China. The analysis shows that the TWSA amplitudes of Songhua and Liaohe River basins are signiﬁcantly higher than those of the other basins. Moreover, TWSA sequences in each river basin contain annual seasonal signals, and the wave peaks of TWSA estimates emerge between June and July. Overall, DWLIM provides a useful measure to derive TWSA in regions where GNSS stations are uneven or sparse.


Introduction
Terrestrial water storage (TWS) comprises all of the water stored on the crustal surface and underground, including snow, glaciers, soil water, groundwater, runoff, and biological water components, which is an essential part of the water cycle system [1,2]. However, the TWS is extraordinarily limited, only accounting for 3.47% of the total global water resources [3]. The TWS provides an essential function for industry, agriculture, and human Remote Sens. 2022, 14, 535 3 of 18 to harsh geo-climatic conditions, which dramatically limits the application of GNSS for TWSA inversion [3]. Developing methods to accurately derive TWSA based on sparse GNSS arrays has become a research hotspot.
Unlike previous studies, this study proposes a new deep learning weight loading inversion model (DWLIM) by combining the long short-term memory (LSTM) algorithm, inverse distance weight method, and crustal loading model. Moreover, TWSA was derived for mainland China from 2011 to 2020 using DWLIM, GRACE, and GLDAS. The TWSA results were calculated based on DWLIM, and its variation characteristics were investigated in 10 river basins within China. The organization of this study is as follows: Section 2 describes the materials and methods in this study, and Section 3 presents the TWSA results based on DWLIM, including the inversion of TWSA and validation of DWLIM. Section 4 discusses the variation characteristics of TWSA in the river basins, and this section also analyzes the difference among the TWSA results. Finally, the primary findings of this study are summarized in Section 5.

GNSS Datasets
This study utilized GNSS vertical deformation sequences provided by CMONOC, and the distribution of the GNSS stations is shown in Figure 1b. The period of each station is not consistent due to the difference in the station establishment time, and the periods of GNSS arrays are shown in Figure 1a. The study period was chosen as 2011-2020 to ensure the completeness of vertical deformation sequences. There were 263 original GNSS stations after removing 6 stations with large period differences, which are shown by the red shadow in Figure 1a. The GNSS observation sequences were calculated using observation, navigation, precision ephemeris, and table files. Furthermore, the daily coordinate solution file was obtained based on GAMIT/GLOBK 10.4, and its specific solution strategy is shown in Table 1 [42]. The GNSS vertical sequences were preprocessed by removing observed outliers that were three times larger than the standard error and system sequence errors caused by earthquakes or antenna replacement.
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 20 ago, which makes it possible to obtain the crustal deformation over mainland China [38,39]. The GNSS datasets provided by CMONOC have been widely utilized to analyze crustal deformation and surface loading [22,40,41]. However, the distribution of GNSS stations is uneven due to harsh geo-climatic conditions, which dramatically limits the application of GNSS for TWSA inversion [3]. Developing methods to accurately derive TWSA based on sparse GNSS arrays has become a research hotspot.
Unlike previous studies, this study proposes a new deep learning weight loading inversion model (DWLIM) by combining the long short-term memory (LSTM) algorithm, inverse distance weight method, and crustal loading model. Moreover, TWSA was derived for mainland China from 2011 to 2020 using DWLIM, GRACE, and GLDAS. The TWSA results were calculated based on DWLIM, and its variation characteristics were investigated in 10 river basins within China. The organization of this study is as follows: Section 2 describes the materials and methods in this study, and Section 3 presents the TWSA results based on DWLIM, including the inversion of TWSA and validation of DWLIM. Section 4 discusses the variation characteristics of TWSA in the river basins, and this section also analyzes the difference among the TWSA results. Finally, the primary findings of this study are summarized in Section 5.

GNSS Datasets
This study utilized GNSS vertical deformation sequences provided by CMONOC, and the distribution of the GNSS stations is shown in Figure 1b. The period of each station is not consistent due to the difference in the station establishment time, and the periods of GNSS arrays are shown in Figure 1a. The study period was chosen as 2011-2020 to ensure the completeness of vertical deformation sequences. There were 263 original GNSS stations after removing 6 stations with large period differences, which are shown by the red shadow in Figure 1a. The GNSS observation sequences were calculated using observation, navigation, precision ephemeris, and table files. Furthermore, the daily coordinate solution file was obtained based on GAMIT/GLOBK 10.4, and its specific solution strategy is shown in Table 1 [42]. The GNSS vertical sequences were preprocessed by removing observed outliers that were three times larger than the standard error and system sequence errors caused by earthquakes or antenna replacement.  The primary mission of GRACE satellites is to monitor spatiotemporal variations in the Earth's gravity field on a global scale. Specifically, the gravity field anomaly is not only related to the Earth's rotation but also affected by geophysical phenomena, such as earthquakes, glacial equilibrium adjustments, and oceanic and hydrological changes [29,43]. The gravity variations in GRACE inversions are generally attributed to the large-scale hydrological migration in mainland China. To verify the reliability of DWLIM, this study employed GRACE Mascon (GRACE-M) to compare its results with the DWLIM outcomes. However, the difference in solution strategies causes considerable uncertainty in the single GRACE-M solution. This study utilized the GRACE-M products obtained from 2011 to 2020 provided by the Center for Space Research (CSR) and the Jet Propulsion Laboratory (JPL) of NASA. The TWSA datasets in mainland China were extracted from the boundary files, and the mean value of the two products was considered the final GRACE-M result. Moreover, we did not add additional smoothing, empirical destriping, filtering, or a scaling factor. To compare DWLIM and GLDAS, the mean datasets of GRACE-M were corrected by first-order terms.  [46]. The LSTM model is trained by constructing memory storage units and using a temporal backpropagation algorithm. This algorithm can solve the problem of gradient disappearance in RNN, and it has no long-term dependence. The standard LSTM model mainly consists of the following: each step t with its corresponding input sequence X: x 1 , x 2 , . . . , x t , the input gate i t , the forget gate f t , and the output gate o t . The memory unit c t Remote Sens. 2022, 14, 535 5 of 18 can control the memory and forget the data through different gates, and it is calculated as follows [46].
The memory unit c j t of the j LSTM with unit time t can be expressed as follows [46].
When the memory cell is updated, the current hidden layer h j t can be calculated [46].
where W denotes the weight matrix of the input process; U denotes the state transfer weight matrix, which is an S-shaped function; tanh denotes the hyperbolic tangent function; σ denotes the sigmoid function; h t denotes the hidden state vector of the output; and c t denotes the new matrix after updating. The three types of gates jointly control the information entering and leaving the memory cell. The input gate regulates the new information entering the memory cell. The forget gate controls how much information is kept in the memory cell, and the output gate defines how much information can be output. The gate structure of LSTM causes the information in the time series to form a balanced long-and short-term dependence for multiple regression purposes. The original sequences were decomposed into n feature signals based on MEEMD due to the few input sequences in this study. The decomposition process is described as follows [47].
where F denotes the original feature sequence, IMF 1 -IMF n denote the n modal components obtained by decomposing the original sequence, and noi w denotes the Gaussian white noise added by MEEMD to be decomposed in the decomposition process. The geophysical parameters show similar characteristics over a small-scale region. There is a homologous amplitude of crustal deformation where the grid is adjacent to the GNSS station. Therefore, the distance between the grid and GNSS station is considered by using the algorithm of inverse distance weight. The application of inverse distance weight contains three steps. Firstly, the figure center of the grid is regarded as the location coordinates for calculating the distance. Secondly, the distance between the simulated grid and the control GNSS stations is calculated. Finally, we assign the weight to each simulated sequence based on the algorithm of inverse distance weight. The simulated formula is as follows.
where D g denotes the simulated crustal deformation in the unobserved grid by DWLIM; d j represents the distance between the center of the grid and the control station;

The Crustal Load-Deformation Model
The upper part of the continental crust can be considered an elastic layer; it will cause the elastic response of the surface to settle or rebound when the mass of the Earth's surface changes. This deformation is also called crustal load-deformation. Crustal loaddeformation occurs not only in the vertical direction but also in the horizontal direction. Crustal load-deformation is more sensitive in the vertical direction than that in the horizontal direction. The relationship between crustal loading and crustal load-deformation can be established by the Green function [48], which is calculated as follows.
where θ denotes the angular radius from the center of the disk; P n denotes the Legendre polynomials; G denotes Newton universal gravitational constant, which is equal to 6.67 × 10 −11 N × m 2 /kg 2 ; R denotes the radius of the Earth; h n denotes the loading Love number; and g denotes the acceleration of gravity. DWLIM utilizes hydrological deformation as the input data, and it combines the crustal loading inversion model to obtain the TWSA in the study region. In the crustal loading model, the obtained solutions are regularized using a curvature smoothing algorithm, and the solutions are added as constraints in the solution matrix. In other words, the least-squares problem is minimized to estimate the daily terrestrial water storage variability for each segment of time studied [49].
where U green denotes the coefficient matrix of the Green function obtained by Equation (9); σ denotes the standard deviation of the hydrological load-deformation sequence; d denotes the hydrological load-deformation time series, including the simulated U grid and U GNSS ; L denotes the Laplace operator; and β denotes the smoothing factor.

Construction of DWLIM
Broadly speaking, the hydrology and atmosphere on the surface exert stress on the continental crust. At the same time, the crust will produce corresponding elastic deformation when the stress is less than the elasticity of the crustal rocks [50,51]. Fortunately, GNSS can accurately observe crustal deformation with submillimeter accuracy [52]. In recent years, the crustal load-deformation model has been employed to invert the local TWSA in regions where GNSS stations are densely distributed [53][54][55]. However, the distribution of GNSS stations is uneven worldwide due to the influence of geographic conditions [12]. Sparsely distributed GNSS arrays cannot be used to accurately invert TWSA because of the limitation of the disk expansion radius. Therefore, it is one of the keys for accurately deriving TWSA to accurately simulate surface load-deformation in the unobserved regions. In this study, DWLIM was constructed by combining LSTM, the inverse distance weighting method, and the crustal load model. The specific process of DWLIM can be divided into the following five steps. (1) Step I: The study region is divided into 1 • × 1 • grids, and the grids are divided into two situations; specifically, the grids contain or do not contain GNSS stations. This algorithm will proceed to step II if the grid has GNSS stations. Moreover, the grid will be defined as an unobserved grid if it does not contain GNSS stations, and the vertical deformation will be simulated in step III. (2) Step II: The GNSS coordinate solution will be calculated by using observation, precision ephemeris, navigation, and table files based on GAMIT software [42]. The daily coordinates are calculated by the GLOBK software based on baseline data files (h-files), and series outliers and step terms that are three times larger than the standard deviation are removed. (3) Step III: The surface temperature sequence (S T ) and atmosphere pressure sequence (S AP ) are normalized on the grid scale. Furthermore, the normalized results are decomposed using the modified ensemble empirical mode decomposition (MEEMD) method to obtain 2 n feature sequences, including n S T and n S AP feature sequences.
In the unobserved grid, the GNSS vertical deformation sequences are employed as the target sequences, and the 2 n feature sequences are utilized as the input sequences. Then, the LSTM regression method and the inverse distance weight method are employed to simulate the vertical displacement. (4) Step IV: The corrected sequences of atmospheric (NTAL) and non-tidal ocean loading (NTOL) are employed to obtain the hydrologic deformation in all the grids, including the GNSS grids and unobserved grids [56]. (5) Step V: The TWSA results are obtained by combining the Green function and the inversion of the crustal load model with all hydrologic deformation. The flow chart of this study is shown in Figure 2.
two situations; specifically, the grids contain or do not contain GNSS stations. This algorithm will proceed to step II if the grid has GNSS stations. Moreover, the grid will be defined as an unobserved grid if it does not contain GNSS stations, and the vertical deformation will be simulated in step III. (2) Step II: The GNSS coordinate solution will be calculated by using observation, precision ephemeris, navigation, and table files based on GAMIT software [42]. The daily coordinates are calculated by the GLOBK software based on baseline data files (hfiles), and series outliers and step terms that are three times larger than the standard deviation are removed. (3) Step III: The surface temperature sequence (ST) and atmosphere pressure sequence (SAP) are normalized on the grid scale. Furthermore, the normalized results are decomposed using the modified ensemble empirical mode decomposition (MEEMD) method to obtain 2 n feature sequences, including n ST and n SAP feature sequences.
In the unobserved grid, the GNSS vertical deformation sequences are employed as the target sequences, and the 2 n feature sequences are utilized as the input sequences. Then, the LSTM regression method and the inverse distance weight method are employed to simulate the vertical displacement. (4) Step IV: The corrected sequences of atmospheric (NTAL) and non-tidal ocean loading (NTOL) are employed to obtain the hydrologic deformation in all the grids, including the GNSS grids and unobserved grids [56]. (5) Step V: The TWSA results are obtained by combining the Green function and the inversion of the crustal load model with all hydrologic deformation. The flow chart of this study is shown in Figure 2.
where Y and X denote accurate and simulated data, respectively, and Y and X represent the mean value of data. The RMSE can be employed to evaluate the deviation of the inversion results from the actual values. The smaller the value of RMSE, the better the simulation accuracy. The NSE is mainly used to evaluate the performance of the hydrological model, and its value is not larger than 1. The larger the value, the better the hydrological model. When NSE is close to 0, it indicates that the effect of the hydrological model agrees with the Seventy-five GNSS sites were selected in grids where the PCC values between the atmospheric pressure or temperature sequence and the GNSS sequence were greater than 0.5. The 75 GNSS sites were used as control sequences for the regression of LSTM, and 263 GNSS vertical sequences were employed for the validation of regression. It was regressed 74 times when the grid contained control GNSS stations, and it was regressed 75 times when the grid did not contain control GNSS stations. The inverse distance weight was employed to assign weights for 74 or 75 simulations. Furthermore, the GNSS vertical sequences were utilized as the true data to verify the accuracy of regression. The simulated results were contrasted with the GNSS vertical sequence according to the RMSE and PCC. The evaluation results are shown in Figure 3.  It can be seen from Figure 3 that most of the station sequences have RMSE values within 5 mm. The variability of the observation quality among GNSS sequences may lead to large RMSE values for some stations. The statistics of the evaluation index indicate that 68.63% of the RMSE values are within 6 mm. The PCC index was used to evaluate the consistency between the simulated sequences and the in situ measurements; the largest PCC value reaches 0.87, and its mean value is 0.53. Figure 3d shows the mean sequences of simulated results and the true data in China. The features of the annual amplitude are included in the simulated outcomes, and the mean simulated sequence is smoother than the GNSS vertical sequence.

Simulation of Hydrological Load-Deformation
(1) Simulation of crustal deformation It can be seen from Figure 3 that most of the station sequences have RMSE values within 5 mm. The variability of the observation quality among GNSS sequences may lead to large RMSE values for some stations. The statistics of the evaluation index indicate that 68.63% of the RMSE values are within 6 mm. The PCC index was used to evaluate the consistency between the simulated sequences and the in situ measurements; the largest PCC value reaches 0.87, and its mean value is 0.53. Figure 3d shows the mean sequences of simulated results and the true data in China. The features of the annual amplitude are included in the simulated outcomes, and the mean simulated sequence is smoother than the GNSS vertical sequence.

Simulation of Hydrological Load-Deformation
(1) Simulation of crustal deformation In this study, the surface temperature and atmospheric pressure were utilized as the input data for the LSTM algorithm, and the models were established by using the 75 control GNSS vertical sequences. Furthermore, the MEEMD method was employed to decompose the surface temperature and atmospheric pressure sequences into 10 model components, IMF1-IMF10, respectively. The decomposition of the input sequences is shown in Figure 4, and the G456 grid is shown as an example.  The IMF1 components in Figure 4a,b are the normalized original surface temperature and atmospheric pressure sequences, respectively. IMF2-IMF10 are the decomposed feature sequences from high to low frequencies. Specifically, the decomposed results reflect the trend and seasonal and residual terms of the series. In the LSTM regression method, the 10 IMF components were used as the input sequences, and the GNSS vertical deformation sequence was used as output data. Furthermore, the inverse distance weight was used to assign weights to the 75 simulated vertical displacements. The vertical simulated deformation was obtained in the unobserved grids. The distribution between the unknown grids and GNSS stations is shown in Figure 5a, and the simulated results of the unknown grid are shown in Figure 5b  It can be seen from Figure 5a that the GNSS stations (blue points) are unevenly distributed, which cannot achieve the overall coverage of the crust over mainland China. The IMF1 components in Figure 4a,b are the normalized original surface temperature and atmospheric pressure sequences, respectively. IMF2-IMF10 are the decomposed feature sequences from high to low frequencies. Specifically, the decomposed results reflect the trend and seasonal and residual terms of the series. In the LSTM regression method, the 10 IMF components were used as the input sequences, and the GNSS vertical deformation sequence was used as output data. Furthermore, the inverse distance weight was used to assign weights to the 75 simulated vertical displacements. The vertical simulated deformation was obtained in the unobserved grids. The distribution between the unknown grids and GNSS stations is shown in Figure 5a, and the simulated results of the unknown grid are shown in Figure 5b  The IMF1 components in Figure 4a,b are the normalized original surface temperature and atmospheric pressure sequences, respectively. IMF2-IMF10 are the decomposed feature sequences from high to low frequencies. Specifically, the decomposed results reflect the trend and seasonal and residual terms of the series. In the LSTM regression method, the 10 IMF components were used as the input sequences, and the GNSS vertical deformation sequence was used as output data. Furthermore, the inverse distance weight was used to assign weights to the 75 simulated vertical displacements. The vertical simulated deformation was obtained in the unobserved grids. The distribution between the unknown grids and GNSS stations is shown in Figure 5a, and the simulated results of the unknown grid are shown in Figure 5b  It can be seen from Figure 5a that the GNSS stations (blue points) are unevenly distributed, which cannot achieve the overall coverage of the crust over mainland China. Hence, the simulation of vertical deformation in unknown grids (yellow points) is essential. Figure 5b-d presents the simulated outcomes of vertical crustal deformation in unob- It can be seen from Figure 5a that the GNSS stations (blue points) are unevenly distributed, which cannot achieve the overall coverage of the crust over mainland China. Hence, the simulation of vertical deformation in unknown grids (yellow points) is essential. Figure 5b-d presents the simulated outcomes of vertical crustal deformation in unobserved grids using 20 IMF feature components for LSTM regression [60]. The results show that the period term and annual amplitude of the vertical crustal deformation can be well simulated according to this strategy, which provides a reasonable data basis for the inversion of TWSA.
(2) Correction of all deformation sequences This study used the NTAL and NTOL models as correction data to extract the crustal deformation caused by hydrological loading. The two corrected sequences were added to the crustal load-deformation time series in mainland China, including the GNSS vertical deformation and simulated results in the unobserved grids. Furthermore, the annual amplitudes of NATL and NOTL were calculated from 2011 to 2020, as shown in Figure 6a,b, respectively. To evaluate the performance of the correction as a whole, the mean sequence of the vertical deformation and hydrologic displacement was obtained, as shown in Figure 6c.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 20 simulated according to this strategy, which provides a reasonable data basis for the inversion of TWSA. (

2) Correction of all deformation sequences
This study used the NTAL and NTOL models as correction data to extract the crustal deformation caused by hydrological loading. The two corrected sequences were added to the crustal load-deformation time series in mainland China, including the GNSS vertical deformation and simulated results in the unobserved grids. Furthermore, the annual amplitudes of NATL and NOTL were calculated from 2011 to 2020, as shown in Figure 6a,b, respectively. To evaluate the performance of the correction as a whole, the mean sequence of the vertical deformation and hydrologic displacement was obtained, as shown in Figure  6c.  Figure 6 indicates that the mean values of the amplitude of the load-deformation of NATL and NOTL are equal to 3.62 and 0.22 mm, respectively. The raised region of the NTAL annual amplitude is mainly located in northern and eastern China, with a maximum of 5.5 mm. However, the maximum annual amplitude of NTOL is only 1.5 mm, and it is mainly distributed in the eastern coastal regions of China. It can be seen from Figure  6c that there are smaller variations in the amplitude and phase of the corrected sequence. The corrections of NTAL and NTOL provide accurate hydrological load-deformation sequences for DWLIM inversion of TWSA.

Inversion of TWSA Based on DWLIM
The Green function matrix of the point loadings was calculated, and the spatial resolution of the inversion outcomes is 0.25° × 0.25°. The expansion boundary range and β equal 2° and 0.01, respectively. Hydrologic displacement sequences were used as the input data for Equations (9) and (10) to calculate the daily TWSA in mainland China. To verify the accuracy of the DWLIM results, first-order term correction was applied to the inversion results in this study, including TWSA results of DWLIM, GRACE, and GLDAS.  Figure 6 indicates that the mean values of the amplitude of the load-deformation of NATL and NOTL are equal to 3.62 and 0.22 mm, respectively. The raised region of the NTAL annual amplitude is mainly located in northern and eastern China, with a maximum of 5.5 mm. However, the maximum annual amplitude of NTOL is only 1.5 mm, and it is mainly distributed in the eastern coastal regions of China. It can be seen from Figure 6c that there are smaller variations in the amplitude and phase of the corrected sequence. The corrections of NTAL and NTOL provide accurate hydrological load-deformation sequences for DWLIM inversion of TWSA.

Inversion of TWSA Based on DWLIM
The Green function matrix of the point loadings was calculated, and the spatial resolution of the inversion outcomes is 0.25 • × 0.25 • . The expansion boundary range and β equal 2 • and 0.01, respectively. Hydrologic displacement sequences were used as the input data for Equations (9) and (10) to calculate the daily TWSA in mainland China. To verify the accuracy of the DWLIM results, first-order term correction was applied to the inversion results in this study, including TWSA results of DWLIM, GRACE, and GLDAS. The calculated annual amplitudes and mean sequence of the TWSA results are shown in Figure 7a, The calculated annual amplitudes and mean sequence of the TWSA results are shown in Figure 7a,b. It can be seen from Figure 7 that the annual characteristics and amplitudes of TWSA sequences based on DWLIM can be effectively inverted. The raised regions of the annual amplitudes can be calculated using DWLIM, which are located in Yunnan Province, southern Tibet region, and southern North China Plain. Furthermore, this result is consistent with the inversion conclusions of previous studies [8]. To verify the accuracy of DWLIM, the outcomes of this study were compared with the inverted TWSA from the GRACE, GLDAS, and the traditional GNSS-derived method.

Spatial Verification of TWSA Results
In order to compare the accuracy of TWSA based on DWLIM, this study obtained TWSA using the traditional GNSS TWSA inversion method (TRAGNSS), GRACE-M datasets, and the GLDAS hydrological model. The detailed information of these outcomes is summarized in Table 2. Additionally, the annual amplitudes of these TWSA results were calculated, and the results are shown in Figure 8.  It can be seen from Figure 7 that the annual characteristics and amplitudes of TWSA sequences based on DWLIM can be effectively inverted. The raised regions of the annual amplitudes can be calculated using DWLIM, which are located in Yunnan Province, southern Tibet region, and southern North China Plain. Furthermore, this result is consistent with the inversion conclusions of previous studies [8]. To verify the accuracy of DWLIM, the outcomes of this study were compared with the inverted TWSA from the GRACE, GLDAS, and the traditional GNSS-derived method.

Spatial Verification of TWSA Results
In order to compare the accuracy of TWSA based on DWLIM, this study obtained TWSA using the traditional GNSS TWSA inversion method (TRA GNSS ), GRACE-M datasets, and the GLDAS hydrological model. The detailed information of these outcomes is summarized in Table 2. Additionally, the annual amplitudes of these TWSA results were calculated, and the results are shown in Figure 8.  Figure 8 indicates that the DWLIM strategy can effectively invert the raised regions of annual amplitude in mainland China, such as southwestern Yunnan Province, southeast China, and the Qinghai-Tibet region. Overall, the spatial amplitude results of DWLIM are consistent with the outcomes of GRACE and GLDAS. However, the annual amplitude of DWLIM is slightly larger than that of GRACE and GLDAS. The reason is that the influence of crustal deformation is complex, and hydrological displacement cannot be completely extracted using NTAL and NTOL. Specifically, the raised regions of annual amplitude also contain northern Xinjiang and northern Heilongjiang. The spatial distribution of the annual amplitude based on the traditional GNSS-derived TWSA method contains speckle characteristics because of the distance limitation of the radius. Hence, the TWSA results based on the traditional GNSS inversion method can only infer the range around the GNSS stations. This will lead to missed signals in regions with sparse GNSS stations when smoothing, and it greatly limits the application of GNSS for TWSA inversion. Overall, the limitation of the disk radius on the GNSS TWSA inversion can be mitigated by simulating crustal deformation in the unknown grids.  Figure 8 indicates that the DWLIM strategy can effectively invert the raised regions of annual amplitude in mainland China, such as southwestern Yunnan Province, southeast China, and the Qinghai-Tibet region. Overall, the spatial amplitude results of DWLIM are consistent with the outcomes of GRACE and GLDAS. However, the annual amplitude of DWLIM is slightly larger than that of GRACE and GLDAS. The reason is that the influence of crustal deformation is complex, and hydrological displacement cannot be completely extracted using NTAL and NTOL. Specifically, the raised regions of annual amplitude also contain northern Xinjiang and northern Heilongjiang. The spatial distribution of the annual amplitude based on the traditional GNSS-derived TWSA method contains speckle characteristics because of the distance limitation of the radius. Hence, the TWSA results based on the traditional GNSS inversion method can only infer the range around the GNSS stations. This will lead to missed signals in regions with sparse GNSS stations when smoothing, and it greatly limits the application of GNSS for TWSA inversion. Overall, the limitation of the disk radius on the GNSS TWSA inversion can be mitigated by simulating crustal deformation in the unknown grids.

Temporal Verification of the TWSA Results
In order to verify the time series reliability of DWLIM, the DWLIM results were compared with the results of the traditional GNSS TWSA inversion method, GRACE, and GLDAS. To further analyze the relationship between DWLIM inversion results and the results of the other data, cross-wavelet analysis was performed, as shown in Figure 9a-c,

Temporal Verification of the TWSA Results
In order to verify the time series reliability of DWLIM, the DWLIM results were compared with the results of the traditional GNSS TWSA inversion method, GRACE, and GLDAS. To further analyze the relationship between DWLIM inversion results and the results of the other data, cross-wavelet analysis was performed, as shown in Figure 9a-c, respectively. In addition, the mean sequences of DWLIM, traditional GNSS, GRACE, and GLDAS over mainland China are shown in Figure 9d.
It can be seen from Figure 9a-c that the TWSA results of DWLIM are consistent with the TWSA of the traditional GNSS inversion method, GLDAS, and GRACE. In addition, the resonance periods between the DWLIM and the other data are about one year, which is shown by the red strip. DWLIM can effectively derive the annual and semiannual amplitudes of the TWSA sequences, which is consistent with the GRACE and GLDAS results (Figure 9d). However, the annual amplitude of DWLIM is slightly larger than the other TWSA results due to the difference in the observation strategy. Moreover, the corrected crustal deformation sequences also contain other deformation signals, resulting in the inability to separate single hydrological load-deformation sequences. The seasonal feature of the DWLIM results is more pronounced that of the traditional GNSS-derived results. To quantify the advantages of DWLIM over the traditional GNSS TWSA inversion method, this study evaluated the inversion results using PCC, NSE, and RMSE. The results are shown in Figure 10. respectively. In addition, the mean sequences of DWLIM, traditional GNSS, GRACE, and GLDAS over mainland China are shown in Figure 9d. It can be seen from Figure 9a-c that the TWSA results of DWLIM are consistent with the TWSA of the traditional GNSS inversion method, GLDAS, and GRACE. In addition, the resonance periods between the DWLIM and the other data are about one year, which is shown by the red strip. DWLIM can effectively derive the annual and semiannual amplitudes of the TWSA sequences, which is consistent with the GRACE and GLDAS results ( Figure 9d). However, the annual amplitude of DWLIM is slightly larger than the other TWSA results due to the difference in the observation strategy. Moreover, the corrected crustal deformation sequences also contain other deformation signals, resulting in the inability to separate single hydrological load-deformation sequences. The seasonal feature of the DWLIM results is more pronounced that of the traditional GNSS-derived results. To quantify the advantages of DWLIM over the traditional GNSS TWSA inversion method, this study evaluated the inversion results using PCC, NSE, and RMSE. The results are shown in Figure 10.  It can be seen from Figure 9a-c that the TWSA results of DWLIM are consistent with the TWSA of the traditional GNSS inversion method, GLDAS, and GRACE. In addition, the resonance periods between the DWLIM and the other data are about one year, which is shown by the red strip. DWLIM can effectively derive the annual and semiannual amplitudes of the TWSA sequences, which is consistent with the GRACE and GLDAS results (Figure 9d). However, the annual amplitude of DWLIM is slightly larger than the other TWSA results due to the difference in the observation strategy. Moreover, the corrected crustal deformation sequences also contain other deformation signals, resulting in the inability to separate single hydrological load-deformation sequences. The seasonal feature of the DWLIM results is more pronounced that of the traditional GNSS-derived results.
To quantify the advantages of DWLIM over the traditional GNSS TWSA inversion method, this study evaluated the inversion results using PCC, NSE, and RMSE. The results are shown in Figure 10. The results further demonstrate that DWLIM can effectively derive TWSA in regions with sparse GNSS stations. Furthermore, the TWSA of DWLIM is better than the traditional GNSS-derived method in terms of spatial and temporal characteristics.

Comparison with Precipitation over 10 River Basins
It is verified that DWLIM can effectively derive the TWSA in mainland China, and it can detect the raised regions of the TWSA annual amplitude. The crust shows a decreasing trend when the terrestrial water storage load increases. On the contrary, the crust shows an upward rebound trend when the terrestrial water storage load decreases. This study combined monthly precipitation products provided by the China Meteorological Administration (CMA) to analyze the variation characteristics of regional TWSA in mainland China. Furthermore, this study extracted the precipitation and TWSA of 10 river basins in China based on boundary files. TWSA was calculated by DWLIM, and the comparison is shown in Figure 11.  Figure 11 indicates that the annual amplitude of TWSA is generally positively correlated with the annual amplitude of precipitation. The mean precipitation sequences in the Songhua and Liaohe River basins are significantly higher than the others. Correspondingly, the amplitudes of TWSA results are also significantly higher than those in the other basins. The phase relationship between TWSA and precipitation in mainland China shows good consistency. This further indicates the reliability of TWSA in phase for DWLIM inversion in mainland China. However, the sequences of TWSA based on DWLIM and precipitation contain delays on the scale of months due to the time needed for the elastic deformation of TWSA. The results of TWSA and precipitation are consistent with previous studies [8]. The seasonal items of TWSA outcomes are more regular than previous TWSA results. Furthermore, the amplitude performance of TWSA and precipitation can also be used to evaluate the arid situation over the river basins. At the same time, it can also be seen that there is high-frequency noise in the time series of TWSA sequences, which also affects the inversion or prediction of TWSA. It is mainly caused by systematic noise from ionospheric, tropospheric, clock error, and multipath effects during GNSS observations [52,61]. Therefore, we will also focus on the noise classification and removal of GNSS vertical sequences to provide cleaner sequences for TWSA inversion in future research.

Discussion of the Difference between Products
In this study, we utilized DWLIM, GRACE, GLDAS, and the traditional GNSS method to calculate TWSA over mainland China. We compared these TWSA outcomes from the perspectives of spatial amplitude ( Figure 8) and time series (Figure 9). It can be seen from Figure 8 that DWLIM is consistent with GRACE and GLDAS over most regions. However, there are also some differences between DWLIM and other products over certain regions, such as Beijing. The reasons for this can be summarized as follows. First, there are only two available GNSS stations (BJFS and BJSH) over Beijing. Second, vertical crustal deformation in the entirety of the North China Plain is complex and has been greatly influenced by human activity, which can cause inaccuracies in the simulated deformation. Third, the GNSS inversion result is also a little higher than the other products because the loaddeformation contains other components. Therefore, there may be some differences between the results of DWLIM and GRACE and GLDAS in some regions. Furthermore, it can be seen from Figure 8a,b that DWLIM can effectively suppress the speckle effect caused by uneven distribution of GNSS stations. In future research, we will focus on extracting cleaner crustal hydrological load-deformation to increase the accuracy of the inversion results.

Conclusions
The main research results can be summarized in the following three points.
(1) To increase the derived accuracy for TWSA, DWLIM was constructed by combining LSTM, inverse distance weight, and the crustal load-deformation model. First, the study region was divided into 1 • × 1 • grids, and then we determined whether the grid contained GNSS stations. Second, this study selected the surface temperature and atmospheric pressure as input data, and the GNSS vertical sequences were utilized as the output data. Each unobserved grid was simulated 263 times, and the inverse distance weight was used to calculate the weighted sequence. Third, the NTAL and NTOL models were employed to correct vertical deformation over all of the grids to obtain the hydrologic distribution. Finally, all of the corrected sequences were used as the input data for the crustal load model to derive TWSA in mainland China. (2) To verify the accuracy of DWLIM, the TWSA results of DWLIM were compared with the traditional GNSS TWSA inversion, GRACE, and GLDAS results. The results indicate that the annual amplitude distribution of DWLIM is smoother than the traditional GNSS inversion results. The strategy of DWLIM greatly suppresses the effect of a small disk expansion radius. The maximum PCC, NSE, and RMSE of DWLIM results compared with GRACE and GLDAS are equal to 0.81, 0.62, and 2.18 cm, respectively, which are improved by 67.11, 128.15, and 22.75% compared with the traditional GNSS-derived TWSA method, respectively. Overall, the DWLIM can effectively invert the TWSA in regions with an uneven distribution of GNSS stations. (3) This study employed precipitation data to analyze the relationship between TWSA and rainfall. We inverted TWSA based on DWLIM in 10 river basins of mainland China. The results indicate that TWSA is positively correlated with precipitation. The annual amplitudes of precipitation and TWSA in the Songhua River basin and the Liaohe River basin are significantly higher than those in other basins. Furthermore, the wave peaks of precipitation are in good agreement with the peaks of TWSA, which are located in June or July. This result further verifies the reliability of the DWLIM inversion results in terms of phase.
Author Contributions: All authors collaborated to conduct this study; Y.S.: scientific analysis, manuscript writing, and editing. W.Z. and W.Y.: experimental design, project management, and review and editing. A.X. and H.Z.: review and editing. Q.W. and Z.C.: editing. All authors contributed to the article and approved the submitted version. Y.S., W.Z. and W.Y. contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.