Inverted Algorithm of Terrestrial Water-Storage Anomalies Based on Machine Learning Combined with Load Model and Its Application in Southwest China

: Dense Global Position System (GPS) arrays can be used to invert the terrestrial water-storage anomaly (TWSA) with higher accuracy. However, the uneven distribution of GPS stations greatly limits the application of GPS to derive the TWSA. Aiming to solve this problem, we grid the GPS array using regression to raise the reliability of TWSA inversion. First, the study uses the random forest (RF) model to simulate crustal deformation in unobserved grids. Meanwhile, the new Machine-Learning Loading-Inverted Method (MLLIM) is constructed based on the traditional GPS derived method to raise the truthfulness of TWSA inversion. Second, this research selects southwest China as the study region, the MLLIM and traditional GPS inversion methods are used to derive the TWSA, and the inverted results are contrasted with datasets of the Gravity Recovery and Climate Experiment (GRACE) Mascon and the Global Land Data Assimilation System (GLDAS) model. The comparison shows that values of Pearson Correlation Coefﬁcient (PCC) between the MLLIM and GRACE and GRACE Follow-On (GRACE-FO) are equal to 0.91 and 0.88, respectively; and the values of R-squared ( R 2 ) are equal to 0.76 and 0.65, respectively; the values of PCC and R 2 between MLLIM and GLDAS solutions are equal to 0.79 and 0.65. Compared with the traditional GPS inversion, the MLLIM improves PCC and R 2 by 8.85% and 7.99% on average, which indicates that the MLLIM can improve the accuracy of TWSA inversion more than the traditional GPS method. Third, this study applies the MLLIM to invert the TWSA in each province of southwest China and combines the precipitation to analyze the change of TWSA in each province. The results are as follows: (1) The spatial distribution of TWSA and precipitation is coincident, which is highlighted in southwest Yunnan and southeast Guangxi; (2) this study compares TWSA of MLLIM with GRACE and GLDAS solutions in each province, which indicates that the maximum value of PCC is as high as 0.86 and 0.94, respectively, which indicates the MLLIM can be used to invert the TWSA in the regions with sparse GPS stations. The TWSA based on the MLLIM can be used to ﬁll the vacancy between GRACE and GRACE-FO.


Introduction
Terrestrial water is crucial for supporting industrial, agricultural, and human life in arid and semiarid areas, while accounting for only 3.5% of global water resources. In China, the uneven spatiotemporal distribution of terrestrial water has always been a serious problem for the survival and development of human society [1,2]. In particular, due to China based on the MLLIM. We also divide the TWSA into the five provinces, namely Yunnan, Sichuan, Guangxi, Guizhou, and Chongqing. The precipitation is also considered to analyze the changes of TWSA in the five provinces. The organization of the study is shown as follows: Section 2 introduces the datasets and the methods; Section 3 verifies the reliability of MLLIM and shows the results; Section 4 applies MLLIM to derive the TWSA in each province; finally, Section 5 summarizes the primary findings of this study.

GPS Solutions
This study uses the sequences of GPS provided by CMONOC and the distribution of the 57 GPS sites shown in Figure 1. These daily solutions of CORS are calculated by the files of observation, navigation, and precise ephemeris. The table files are used to correct the errors of the measurement, which includes the collection of the ionosphere, troposphere, absolute antenna phase center, and ocean tide. The GAMIT is used to calculate the baseline among the GPS stations; meanwhile, the GLOBK is used to calculate the coordinated location based on the baselines [40]. However, the obtained sequences contain lots of leap and are anonymous due to the earthquake and instrument factors, so we correct these to acquire the cleaner crustal deformation. The time period of the GPS station is 2011-2019 to avoid the influence of missing data for the TWSA inversion. From Figure 1, there are substantial GPS stations over the Yunnan province; however, the spatial distribution of GPS stations is sparse in eastern Sichuan, Chongqing, Guizhou, and Guangxi.

GRACE Solutions
The spherical harmonic coefficient of gravity field that observed by the GRACE and the GRACE-FO can invert the TWSA effectively [12]. The formula of the inversion is shown as: where the ∆h denotes the height of equivalent water for the TWSA; A denotes the Earth radius that equals 6371.39 km. ρ e represents the density of the Earth that equals 5.5 × 10 3 kg/m 3 ; ρ w denotes the density of the water that equals 10 3 kg/m 3 ; h l and k l represent the Love numbers; W l denotes the kernel function of Gaussian smoothing; P l,m denotes the Legendre functions; ∆C lm and ∆S lm represent the spherical harmonic coefficient of Earth gravity field. The GRACE Mascon (0.5 • × 0.5 • ) is supplied by the Center for Space Research (CSR) and Jet Propulsion Laboratory (JPL) 2001-2019 [41]. Subsequently, the TWSA over southwest China is extracted and the mean result of the TWSA is used as the final GRACE solution, as follows:

Products of GLDAS and ERA 5
To simulate the crustal deformation (U) over the unobserved grids, the sequences of the temperature and the atmospheric pressure are used as the input data for regression. The sequences of surface temperature and atmospheric pressure in the study region are obtained from the GLDAS V2.2 and ERA 5, respectively. The GLDAS V 2.2 can be searched in the Goddard Earth Sciences Data and Information Services Center (GES DISC), and it contains two production streams, one with GRACE data assimilation outputs, and the early production stream without the GRACE DA (data assimilation). Specifically, the datasets of GLDAS V2. 2 Figure 2 shows the sequences of vertical deformation (GPS), temperature, and atmospheric pressure, taking the CORS stations located in Panzhihua of Sichuan (SCPZ), Xiaojin of Sichuan (SCXJ), Mojiang of Yunnan (YNMJ), and Mile of Yunnan (YNML) as examples. The distribution of the four GPS sites is shown in Figure 1b with the red points.
This study uses the sequences of GPS provided by CMONOC and the distribution of the 57 GPS sites shown in Figure 1. These daily solutions of CORS are calculated by the files of observation, navigation, and precise ephemeris. The table files are used to correct the errors of the measurement, which includes the collection of the ionosphere, troposphere, absolute antenna phase center, and ocean tide. The GAMIT is used to calculate the baseline among the GPS stations; meanwhile, the GLOBK is used to calculate the coordinated location based on the baselines [40]. However, the obtained sequences contain lots of leap and are anonymous due to the earthquake and instrument factors, so we correct these to acquire the cleaner crustal deformation. The time period of the GPS station is 2011-2019 to avoid the influence of missing data for the TWSA inversion. From Figure 1, there are substantial GPS stations over the Yunnan province; however, the spatial distribution of GPS stations is sparse in eastern Sichuan, Chongqing, Guizhou, and Guangxi.   Figure 2 indicates that the period of vertical deformation, temperature, and atmospheric pressure are close to one year. There is a negative correlation between the sequences of vertical deformation and atmospheric pressure. In addition, there is a certain phase difference between the sequence of atmospheric pressure and the sequence of vertical deformation. Overall, the periodicity of temperature and atmospheric pressure sequences Remote Sens. 2021, 13, 3358 5 of 19 are more obvious than the GPS. The reason is that the influential factors of GPS are more complex than the sequences of temperature and atmospheric pressure.  Figure 2 indicates that the period of vertical deformation, temperature, and atmospheric pressure are close to one year. There is a negative correlation between the sequences of vertical deformation and atmospheric pressure. In addition, there is a certain phase difference between the sequence of atmospheric pressure and the sequence of vertical deformation. Overall, the periodicity of temperature and atmospheric pressure sequences are more obvious than the GPS. The reason is that the influential factors of GPS are more complex than the sequences of temperature and atmospheric pressure.

Random forest
The RF is a learning method for regression, classification, and other purposes proposed by Breiman in 2001 [42]. It constructs a large number of random, de-correlated decision trees and averages them, which contains many advantages for multiple linear regression [43]. Specifically, the factors can be included or excluded based on the applicability of the data and the RF can combine the predictors of the continuous classification. In the task of regression, fewer parameters of the model need to be defined and there is a little risk in the overfitting. Furthermore, the variable important score can be automatically calculated to evaluate the contribution of the single predictor in the final model. The sequences of surface temperature and atmospheric pressure are used as the input sequences; meanwhile, the GPS vertical sequences are employed as the output sequences. Thus, each unobserved grid is simulated N times by these GPS stations and the average of the N sequences is regarded as the final result of the grid's deformation [43], shown as follows:

Random Forest
The RF is a learning method for regression, classification, and other purposes proposed by Breiman in 2001 [42]. It constructs a large number of random, de-correlated decision trees and averages them, which contains many advantages for multiple linear regression [43]. Specifically, the factors can be included or excluded based on the applicability of the data and the RF can combine the predictors of the continuous classification. In the task of regression, fewer parameters of the model need to be defined and there is a little risk in the overfitting. Furthermore, the variable important score can be automatically calculated to evaluate the contribution of the single predictor in the final model. The sequences of surface temperature and atmospheric pressure are used as the input sequences; meanwhile, the GPS vertical sequences are employed as the output sequences. Thus, each unobserved grid is simulated N times by these GPS stations and the average of the N sequences is regarded as the final result of the grid's deformation [43], shown as follows: where U grid (x) represents the final result of crustal deformation in the unobserved grid; g t (x) denotes each tree's simulated result of RF; T represents the number of spanning trees for RF; N denotes the number of the GPS stations. The load-deformation will be caused when the load of the Earth changes, due to surface water, snow, and ice. Fortunately, the Green's function can establish the relationship between the mass of loading and the crustal deformation [22]. The structural displacement is caused by the internal tectonic locomotion of the Earth, which is exhausted in the linear crustal movement at the horizontal direction (N, E) [44]. The seasonal non-structural displacement is caused by the seasonal crustal surface loadings [45]. However, the crustal load-deformation is sensitive at vertical direction (U), and its amplitude is about 2-3 times the horizontal direction [20]. The vertical load-deformation can be described by Greens' function, shown as follows [46]: where the θ denotes the angular radius from the center to the disk; P n denotes the Legendre polynomials; G represents the constant of Newton's universal gravitation; R denotes the radius of the Earth; h n represents the Love's number; g denotes the constant of the acceleration for the gravity; The deduction of the Γ n is shown as follows [46]: where the θ denotes the angular radius from the center to the disk; P denotes the Legendre polynomials; the expression of the function is as follows when the n equals 0 [46]: Figure 3 displays the relationship between the load-deformation and distance when the disks load different mass, radius, and thickness. It can be seen that the load response in the near field is very obvious, while it is only half of the response in the center of the disk when the distance from the disk's center equals radio. The load response will disappear when the distance from the center of the disk is about five times, i.e., the GPS station can only reflect the change of the TWSA in a limited range [23]. Hence, it is essential to simulate the vertical deformation sequences in the unobserved grids for the inversion of TWSA. The MLLIM derives the 0.25° × 0.25° TWSA using all the corrected crustal sequences, which contain the 57 vertical GPS sequences and the 87 simulated sequences. The curvature smoothing algorithm is used to regularize the obtained solution, which is attached to The MLLIM derives the 0.25 • × 0.25 • TWSA using all the corrected crustal sequences, which contain the 57 vertical GPS sequences and the 87 simulated sequences. The curvature smoothing algorithm is used to regularize the obtained solution, which is attached to the solution matrix as a set of additional constraints [47]. Specifically, the suppression leastsquare problem is minimized to estimate the daily change of land water storage in each period of inversion [23], as follows: where G represents the Green's function calculated by Formula (4); σ denotes the standard deviation of the vertical sequence; b represents crustal displacement, which contains the vertical GPS sequences and the simulated crustal displacement; L denotes the Laplacian operator; β denotes the smoothing factor.

Construction of the MLLIM
The model of crustal loading-deformation was put forward by Argus in 2014 [23], which can invert the seasonal load by vertical deformation. This method is different from the traditional measurements, including well depth, hydrologic model, satellite altimetry, and gravity [48]. The theories of the crustal load-deformation are used to invert the TWSA by the vertical GPS sequences, which can increase the physical character of the GPS solution [39]. Meanwhile, the inverted result can be combined with the GRACE to fill the missing data and provides a scale factor for GRACE solutions [39,49]. However, the inverted result cannot indicate the spatial features in the regions where the distribution of GPS stations is sparse, which greatly limits the application of the inversion [50,51]. Nevertheless, it is critical for the derivation of TWSA to simulate the crustal deformation in the unobserved regions [52]. The algorithm of MLLIM can be divided into the following steps.
Step I: we search the GPS station in the 1 • × 1 • grid; in the event of there being a GPS station in the grid, we will go to Step II.
Step II: the MLLIM will preprocess the GPS solution, removing the leap and exception values of the GPS vertical sequences. However, if the grid does not contain the GPS station, we will move to Step III.
Step III: the RF is used to simulate the deformation using the sequences of temperature [53] and the atmospheric pressure [54]. Specifically, the known 57 GPS stations are taken as the controlled sequences and every grid is simulated 57 times by the RF, then, the mean value of the 57 simulated sequences is regarded as the final simulated result.
Step IV: the crustal displacement caused by atmosphere and non-tidal ocean are corrected by MERRA and Ocean Model for Circulation and Tides (OMCT), respectively. Furthermore, all the corrected sequences (the GPS and the simulated result) are taken as the input data and use the model of the crustal load-deformation to invert the TWSA in the study region [55].
Step V: this study also derives the TWSA by the traditional method of GPS, GRACE Mascon, and GLDAS solutions to demonstrate the reliability of the MLLIM. The technological process of the MLLIM is shown in Figure 4. thermore, all the corrected sequences (the GPS and the simulated result) are taken as the input data and use the model of the crustal load-deformation to invert the TWSA in the study region [55].
Step V: this study also derives the TWSA by the traditional method of GPS, GRACE Mascon, and GLDAS solutions to demonstrate the reliability of the MLLIM. The technological process of the MLLIM is shown in Figure 4.

Evaluation index
This study estimates the accuracy between the inverted and other products (i.e., GRACE Mascon and GLDAS solutions) based on the following statistics metrics: rootmean-square error (RMSE) [56], Pearson's correlation coefficient (PCC) [57], and Rsquared (R 2 ) [58], respectively, as follows: where the Y and S denote true and simulated data, respectively, and the Y and S represent the mean value of data. The RMSE describes a degree of dispersion between actual and simulated sequences. If the value of RMSE is small, it implies that the simulated result is stable in amplitude. The PCC is a linearly dependent coefficient, and it reflects the degree of linear correlation between two sequences. The value of PCC is between −1 and 1, the closer the absolute value is to 1, the stronger the correlation. It can mainly extract the seasonal relationship between sequences. Meanwhile, R 2 is also called the coefficient of

Evaluation Index
This study estimates the accuracy between the inverted and other products (i.e., GRACE Mascon and GLDAS solutions) based on the following statistics metrics: root-mean-square error (RMSE) [56], Pearson's correlation coefficient (PCC) [57], and R-squared (R 2 ) [58], respectively, as follows: where the Y and S denote true and simulated data, respectively, and the Y and S represent the mean value of data. The RMSE describes a degree of dispersion between actual and simulated sequences. If the value of RMSE is small, it implies that the simulated result is stable in amplitude. The PCC is a linearly dependent coefficient, and it reflects the degree of linear correlation between two sequences. The value of PCC is between −1 and 1, the closer the absolute value is to 1, the stronger the correlation. It can mainly extract the seasonal relationship between sequences. Meanwhile, R 2 is also called the coefficient of determination, which can explain the variation degree according to the independent variables. The value of R 2 is between 0 and 1, and the larger it is, the better the fitting effect.

Inversion of the TWSA Based on the MLLIM
The crustal deformation can be divided into the horizontal (N, E) and vertical (U) directions. There are mainly the characteristics of linearity in horizontal deformation, while the character of the vertical deformation is nonlinear [59]. The bedrock near the surface will vibrate up and down periodically when the surface temperature changes due to the thermal barrier shrinkage phenomenon [60]. Correspondingly, the different altitudes among the GPS stations lead to the atmospheric pressure of the region also being different, which affects the amplitude of the vertical deformation [61,62]. Hence, the temperature and the atmospheric pressure are used as the input data for the RF in this study and the GPS vertical sequences are used as the output data. Aiming the error for single-station regression, the vertical sequences of the unknown grids are predicted N times by the known GPS sequences. Furthermore, the average of the N sequences is regarded as the final sequence in the unobserved grid.
To verify the availability of the regression method, the sequences of the temperature and atmospheric pressure of the grids (57 GPS stations) are used as input data, and the vertical GPS sequence is used as the output data. Each GPS vertical sequence is regressed 56 times, which excludes the station itself. The indexes of PCC and RMSE are used to evaluate this regression, and the result is shown in Figure 5a,b, respectively. To display the effectiveness of the regression method, the mean values of the vertical simulation and the real sequences of 57 stations are obtained, as shown in Figure 5c.
The crustal deformation can be divided into the horizontal (N, E) and vertical (U) directions. There are mainly the characteristics of linearity in horizontal deformation, while the character of the vertical deformation is nonlinear [59]. The bedrock near the surface will vibrate up and down periodically when the surface temperature changes due to the thermal barrier shrinkage phenomenon [60]. Correspondingly, the different altitudes among the GPS stations lead to the atmospheric pressure of the region also being different, which affects the amplitude of the vertical deformation [61,62]. Hence, the temperature and the atmospheric pressure are used as the input data for the RF in this study and the GPS vertical sequences are used as the output data. Aiming the error for singlestation regression, the vertical sequences of the unknown grids are predicted N times by the known GPS sequences. Furthermore, the average of the N sequences is regarded as the final sequence in the unobserved grid.
To verify the availability of the regression method, the sequences of the temperature and atmospheric pressure of the grids (57 GPS stations) are used as input data, and the vertical GPS sequence is used as the output data. Each GPS vertical sequence is regressed 56 times, which excludes the station itself. The indexes of PCC and RMSE are used to evaluate this regression, and the result is shown in Figure 5a,b, respectively. To display the effectiveness of the regression method, the mean values of the vertical simulation and the real sequences of 57 stations are obtained, as shown in Figure 5c.    Figure 5 shows that the crustal deformation can be simulated effectively by this simulated method. According to the statistics in Figure 5a, 84.21% values of PCC are more than 0.5, with the maximum value of 0.79. The dispersion degree between the simulated sequences and the real sequences performs well, which can be determined by RMSE (Figure 5b). The 68.42% values of RMSE are less than 5 mm, and the value of RMSE is related to the quality of the GPS sequence. To verify the reliability of the simulation of crustal deformation, the mean simulated, and true sequences are obtained, respectively. Figure 5c indicates that the MLLIM can effectively simulate the periodic term and amplitude of the sequences with the PCC of 0.75, with RMSE of 3.45 mm, respectively. Thus, the GPS vertical deformation sequence can be simulated by this method in the unobserved grid.
Specifically, each grid is derived 57 times and the mean sequence is regarded as the final simulated result, which is used as the input data for the crustal load model.
To invert the TWSA based on the crustal load-displacement, the 57 vertical sequences observed by GPS and the 87 simulated sequences are employed as the training data for the simulated model. The coordinates of the simulated sequence are defined as the center of the grid (Figure 6a) and the effect of the simulated sequences is shown in Figure 6b,c.
sequences and the real sequences performs well, which can be determined by RMSE (Figure 5b). The 68.42% values of RMSE are less than 5 mm, and the value of RMSE is related to the quality of the GPS sequence. To verify the reliability of the simulation of crustal deformation, the mean simulated, and true sequences are obtained, respectively. Figure  5c indicates that the MLLIM can effectively simulate the periodic term and amplitude of the sequences with the PCC of 0.75, with RMSE of 3.45 mm, respectively. Thus, the GPS vertical deformation sequence can be simulated by this method in the unobserved grid. Specifically, each grid is derived 57 times and the mean sequence is regarded as the final simulated result, which is used as the input data for the crustal load model.
To invert the TWSA based on the crustal load-displacement, the 57 vertical sequences observed by GPS and the 87 simulated sequences are employed as the training data for the simulated model. The coordinates of the simulated sequence are defined as the center of the grid (Figure 6a) and the effect of the simulated sequences is shown in Figure 6b,c. It can be seen from Figure 6 that the annual and the semi-annual items can be effectively simulated in the unobserved grids. Then the models of MERRA and OMCT are used to correct the loading of atmospheric [33] and non-tidal ocean [63], respectively. First, we calculate the Green's function over southwest China with the spatial resolution and radius of 0.25° and 13.90 km [24]. Meanwhile, the values of the expanding border and β equal 2° and 0.01. The TWSA in southwest China from 2011 to 2019 can be derived based on Formula (5) and the result of the TWSA is shown in Figure 7. It can be seen from Figure 6 that the annual and the semi-annual items can be effectively simulated in the unobserved grids. Then the models of MERRA and OMCT are used to correct the loading of atmospheric [33] and non-tidal ocean [63], respectively. First, we calculate the Green's function over southwest China with the spatial resolution and radius of 0.25 • and 13.90 km [24].  Figure 7 indicates that the TWSA in the study region can be inverted by the MLLIM. In addition, the inverted TWSA sequence contains the annual and semi-annual items, which is shown in Figure 7c. Collectively, the total water storage of Yunnan, Sichuan, and Guangxi is higher than other provinces. Meanwhile, the annual amplitude of TWSA in  Figure 7 indicates that the TWSA in the study region can be inverted by the MLLIM. In addition, the inverted TWSA sequence contains the annual and semi-annual items, which is shown in Figure 7c. Collectively, the total water storage of Yunnan, Sichuan, and Guangxi is higher than other provinces. Meanwhile, the annual amplitude of TWSA in Yunnan province is higher than that in the surrounding provinces significantly, with an amplitude of 120 mm (Figure 7b). To prove the reliability of the inverted result, the result based on the MLLIM is compared with the traditional GPS inverted method, GRACE, and hydrologic model.

Comparison with Traditional Method of GPS Inversion
To contrast with the traditional inverted method by GPS, the 57 corrected vertical GPS sequences are used as the input data. Meanwhile, the Green's function and the model of load-deformation are combined to invert the TWSA in northwest China. In this test, the derived parameters are consistent with the first inversion. Then the TWSA signals (0.25 • × 0.25 • ) are calculated by the traditional GPS inverted method between 2011 and 2019. Correspondingly, this study calculates the annual, semi-annual, and periodic items in the grids, respectively, the result is shown in Figure 8. As shown in Figure 8, the inverted results of the MLLIM are obviously different from those of traditional GPS inversion. The subgraphs in Figure 8 can be divided into three contrast groups: (I) is a/d, (II) is b/e, and (III) is c/f. Groups (I)-(III) denote the annual, semiannual, and periodic comparison, respectively. As shown in Group (I), since there are fewer GPS stations in Guangxi province than other provinces, some features are ignored due to smoothing when using traditional methods derive TWSA. However, the characteristics of annual amplitude can be clearly reflected in Figure 8a, and this phenomenon is more obvious in Group (II), which can effectively suppress the loss of signal features caused by interpolation. By analyzing the period time of the sequence (Group (III)), there are obvious differences between Guangxi and Guizhou. The MLLIM can effectively retrieve the spatiotemporal characteristics of the annual term in Guangxi province. To further judge the accuracy of the two inversion methods, TWSA results obtained by the two methods are compared with the GRACE inversion and hydrological model (GLDAS V2.2).

Comparison with GRACE
GRACE and GRACE-FO satellites can detect the time-varying gravity field change As shown in Figure 8, the inverted results of the MLLIM are obviously different from those of traditional GPS inversion. The subgraphs in Figure 8 can be divided into three contrast groups: (I) is a/d, (II) is b/e, and (III) is c/f. Groups (I)-(III) denote the annual, semi-annual, and periodic comparison, respectively. As shown in Group (I), since there are fewer GPS stations in Guangxi province than other provinces, some features are ignored due to smoothing when using traditional methods derive TWSA. However, the characteristics of annual amplitude can be clearly reflected in Figure 8a, and this phenomenon is more obvious in Group (II), which can effectively suppress the loss of signal features caused by interpolation. By analyzing the period time of the sequence (Group (III)), there are obvious differences between Guangxi and Guizhou. The MLLIM can effectively retrieve the spatiotemporal characteristics of the annual term in Guangxi province. To further judge the accuracy of the two inversion methods, TWSA results obtained by the two methods are compared with the GRACE inversion and hydrological model (GLDAS V2.2).

Comparison with GRACE
GRACE and GRACE-FO satellites can detect the time-varying gravity field change and retrieve the TWSA with high accuracy. To verify the reliability of the inversion, this study uses the Mascon data results issued by JPL and Center for Space Research at University of Texas, Austin (CSR) institutions and calculates their mean values as the final results in southwest China. Because of there is a gap between GRACE and GRACE-FO satellite, thus, the GRACE Mascon products are divided into the periods of GRACE and GRACE-FO. Meanwhile, to facilitate the comparison with the inverted results of GRACE, the GPS inversion is processed by monthly average. The results are shown in Figure 9.  Figure 9 shows that the results of GRACE Mascon are consistent with the results of MLLIM in the annual and semi-annual amplitude. From Figure 9a, the largest annual amplitude is 120 mm located in Yunnan province, and the annual amplitude of Guangxi province is also prominent. In marked contrast (Figures 8a and 9a), the annual signals can be detected accurately, while the annual signals are ignored in the traditional method (Figure 8d). Figure 9b indicates that the detection effect is prominent for the semi-annual   Figure 9a, the largest annual amplitude is 120 mm located in Yunnan province, and the annual amplitude of Guangxi province is also prominent. In marked contrast (Figures 8a and 9a), the annual signals can be detected accurately, while the annual signals are ignored in the traditional method (Figure 8d). Figure 9b indicates that the detection effect is prominent for the semi-annual signal of the sequence. The raised dots of the semi-annual signal are located in Yunnan, Guizhou, and central Sichuan. Due to the different measurement methods between the GPS and GRACE, there are some differences in the cyclic item between Figures 8c and 9c. However, the overall order of the magnitude is consistent. Then the TWSA results of GPS and GRACE are fitted. The correlations between the GPS and GRACE are 0.91 and 0.88 for the GRACE and GRACE-FO, respectively. In addition, the values of R 2 equal 0.76 and 0.65. Hence, compared with the traditional GPS inverted method, the MLLIM improves the PCC and R 2 by 7.98% and 9.30% on average.

Comparison with the GLDAS Solutions
To better verify the MLLIM, the TWSA of GLDAS V2.2 is used to compare with the inversion of this study. The spatial resolution of GLDAS V2.2 is 0.5 • × 0.5 • , and its temporal resolution is daily. The results of the comparison are shown in Figure 10. Figure 10 indicates that the GLDAS solution is consistent with the experimental results in annual, semi-annual, and period items. The above conclusions show that the TWSA can be retrieved accurately based on MLLIM in the regions where the distribution of GPS stations is uneven. From Figure 10b, the raised dots of the annual amplitudes are detected based on the MLLIM and they are distributed in Yunnan, western Sichuan, and eastern Guangxi provinces. It can be seen from Figure 10c that it is more consistent with the results of this paper (Figure 8b

Comparison with the GLDAS solutions
To better verify the MLLIM, the TWSA of GLDAS V2.2 is used to compare with the inversion of this study. The spatial resolution of GLDAS V2.2 is 0.5° × 0.5°, and its temporal resolution is daily. The results of the comparison are shown in Figure 10.

Discussion
In Section 3, the values of PCC between the inversion and GRACE and GLDAS are 0.90 and 0.79, respectively. Meanwhile, the annual, semi-annual and periodic signals of TWSA sequences can be detected effectively. The crust will produce downward vertical deformation when the water mass increases; conversely, the crust will produce upward rebound displacement. To analyze the relationship between the TWSA and precipitation in southwest China, this study uses the dataset of precipitation provided by the China Meteorological Administration (CMA). The spatial resolution of precipitation is 0.5 • × 0.5 • , and its temporal resolution is monthly. The compared result is shown in Figure 11. 0.90 and 0.79, respectively. Meanwhile, the annual, semi-annual and periodic signals of TWSA sequences can be detected effectively. The crust will produce downward vertical deformation when the water mass increases; conversely, the crust will produce upward rebound displacement. To analyze the relationship between the TWSA and precipitation in southwest China, this study uses the dataset of precipitation provided by the China Meteorological Administration (CMA). The spatial resolution of precipitation is 0.5° × 0.5°, and its temporal resolution is monthly. The compared result is shown in Figure 11.  It can be seen from Figure 11 that the inverted result of TWSA has a good spatiotemporal consistency with GRACE, GLDAS solutions, and precipitation. However, Figure 10c indicates that the GLDAS solutions in Chongqing and Guizhou are abnormal; the GLDAS solutions of Chongqing and Guizhou are not compared. From Figure 11a-e, the space distribution of precipitation exists in the obvious raised regions in southwest Yunnan, central Sichuan, and southeast Guangxi. It is consistent with the inversion of TWSA in this paper ( Figure 8). It can be seen from Figure 11f-j that: (1) the inverted result based on the MLLIM is almost consistent with GRACE and GLDAS, and the phase of the inverted result is better than that of GLDAS solutions. However, due to the difference from the observations, there is a phase difference of about 2 months with the GRACE; (2) the amplitude of the inverted result is slightly larger than the GRACE and GLDAS, because there will be some residual signals reserved in the GPS sequences; (3) The average PCC of the MLLIM with GRACE and GLDAS equals 0.77 and 0.86, respectively. The results show that the TWSA can be inverted precisely by the MLLIM.
However, due to the characters of GPS, there are various noise components in the crustal deformation sequences observed by GPS [64]. It makes the inverted TWSA se-quences contain a large number of noise components. Therefore, the noise components of GPS sequence will be better removed in the follow-up studies. There is a certain phase difference with GRACE due to the characteristics of GPS inversion TWSA. Therefore, the follow-up studies will process the phase difference to achieve better fitting with the GRACE solution.

Conclusions
It is essential for the inversion of TWSA to obtain the crustal deformation effectively, which determines the reliability of the inversion. To improve the accuracy of TWSA inversion by GPS, the MLLIM is proposed by combining the RF and the model of crustal load-deformation in this study. The following useful conclusions are drawn.
(1) To increase the inverted accuracy for TWSA, the MLLIM is constructed using RF and the crustal load model. To simulate the crustal displacement in unobserved grids, the sequences of temperature and atmospheric pressure are used as the input data, and the GPS vertical sequence is used as the output data. Thus, the unobserved grids' crustal deformation can be simulated based on the MLLIM. Furthermore, all the corrected crustal sequences are employed as the input data for the crustal load model; specifically, the corrected sequences include the GPS vertical sequences and the simulated sequences. (2) To verify the reliability of the MLLIM, the results of traditional GPS inverted method (without adding the simulated crustal deformation), GRACE, and GLDAS are used to compare with the TWSA based on the MLLIM. The compared results show that this inverted method can detect the raised region of annual amplitude in Yunnan, central Sichuan, and Guangxi provinces, and is consistent with the results of GRACE Mascon and GLDAS. However, the traditional GPS inverted method loses this characteristic signals due to the fitting of inversion (Figure 8d-f). From the comparison with GRACE, the values of PCC with GRACE and GRACE-FO equal 0.91 and 0.88, and the values of R 2 equal 0.76 and 0.65, respectively. The MLLIM improves PCC and R 2 by 7.98% and 9.30% than traditional method, respectively. From the comparison with GLDAS, the values of PCC and R 2 equal 0.79 and 0.64, respectively, it improves PCC and R 2 by 9.72% and 6.67% compared with the traditional method, respectively. On the whole, compared with the traditional GPS inverted method, the MLLIM effectively improves the accuracy of TWSA inversion. (3) The data of precipitation is combined to analyze the relationship between the TWSA and the rainfall. We apply the MLLIM to derive the TWSA in the five provinces, and the results indicate that the precipitation in Yunnan and Guangxi is significantly higher than other provinces. Meanwhile, the spatial distribution characteristics of TWSA is consistent with the precipitation in southwest China. From the comparison with GRACE and GLDAS solutions, the maximum PCC values equal 0.86 and 0.94, respectively. The above conclusions show that the TWSA can be inverted accurately based on the MLLIM in the region where the distribution of GPS stations is uneven.
Author Contributions: All authors collaborated to conduct this study; Y.S.: scientific analysis, manuscript writing, and editing. W.Z. and W.Y.: experiment design, project management, review and editing. A.X. and H.Z.: review and editing. S.Y. and K.S.: processing of GRACE. 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.