3.2.1. 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.

**Table 2.** Statistical parameters of DWLIM, traditional GNSS TWSA inversion results, GRACE, and GLDAS.


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.** The spatial distribution of TWSA annual amplitude in mainland China. (**a**) The result of DWLIM; (**b**) the result of traditional inversion method based on GNSS; (**c**) the result of GRACE; (**d**) the result of GLDAS.

3.2.2. 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.

**Figure 9.** The analysis and time series of TWSA results. (**a**) Wavelet analysis between DWLIM and traditional GNSS TWSA inversion method; (**b**) wavelet analysis between DWLIM and GLDAS; (**c**) wavelet analysis between DWLIM and GRACE; (**d**) the mean time series of the TWSA results in mainland China.

**Figure 10.** The heat figure of the evaluation index based on DWLIM. (**a**) The value of *PCC* in the GRACE period; (**b**) the value of *NSE* in the GRACE period; (**c**) the value of *RMSE* in the GRACE period; (**d**) the value of *PCC* in the GRACE-FO period; (**e**) the value of *NSE* in the GRACE-FO period; (**f**) the value of *RMSE* in the GRACE-FO period.

It can be seen from Figure 10 that, based on the evaluation indexes *PCC*, *NSE*, and *RMSE*, the TWSA results based on DWLIM are superior to the traditional GNSS-derived results. For the period of the GRACE mission (2011–2017), the maximum *PCC*, *NSE*, and *RMSE* indicators of DWLIM inversion results reach 0.81, 0.62, and 2.18 cm, respectively. For the period of the GRACE-FO mission (2018–2020), the maximum *PCC*, *NSE*, and *RMSE* of DWLIM inversion results reach 0.71, 0.49, and 2.4 cm. The results show that the TWSA results of DWLIM are more consistent with the GLDAS results, which is attributed to the monthly scale resolution of GRACE, leading to signal loss. Further statistics from the data show that the DWLIM results improve the *PCC*, *NSE*, and *RMSE* by 67.11, 128.15, and 22.75% on average compared to the traditional GNSS inversion method, respectively. 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.
