*3.1. Inversion of TWSA Using DWLIM*

3.1.1. Validation of Simulated Crustal Deformation

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.

**Figure 3.** The verification outcomes of simulated crustal deformation. (**a**) The *RMSE* between the simulation and in situ measurement; (**b**) the *PCC* between the derived sequence and true data; (**c**) the Taylor figure of the simulated results; (**d**) the mean simulated sequences in mainland China.

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.


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.

**Figure 4.** The result of normalization and decomposition in the unobserved grids, showing G456 as an example. (**a**) The result of the temperature sequence; (**b**) the result of the atmospheric sequence.

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–d. The G464, G740, and G456 grids are shown as examples.

**Figure 5.** The simulated results in the unobserved grids. (**a**) The distribution of the unobserved grids and GNSS sites; (**b**) the result of G464; (**c**) the result of G740; (**d**) the result of G456.

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.

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.

#### 3.1.3. 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,b.

**Figure 7.** The derived TWSA results based on DWLIM in China. (**a**) The distribution of TWSA annual amplitude in China; (**b**) the mean time series of TWSA in mainland China.

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.
