## 3 Results  
### Equivalent Source Inversion  
#### 3.1.1 Basic Principles of the Inversion  
Changes in the gravity field are not only the results of crustal movements, but they are also affected by density variations in the materials inside the earth [43]. The time-varying gravity field obtained by repeated surface gravity observations has the advantages of a high spatial resolution, observation accuracy, and data reliability, which provides important support for extracting anomalous signals during the earthquake preparation process. Therefore, this field can be inverted using the \"equivalent source\" model to obtain the characteristics of apparent density changes and infer the dynamic density variations in the subsurface media [44].  
The experimental method can be summarized as the tuning and inversion of time-varying gravity field signals. First, the resolution of the source model is assessed and adjusted through a checkerboard test. The source term is discretized into hexahedral grids in a spherical domain [45], and spatiotemporal constraints are introduced. Then, the \"equivalent source\" inversion model is optimized and solved following Akaike Information Criterion and Bayesian Information Criterion (_ABIC_) [46].  
The equation of equivalent source inversion [47] can be expressed as  
\[\Phi=W_{0}|Gm-f(x,y,t)|^{2}+\Phi_{s}(m)+\Phi_{T}(m) \tag{1}\]  
The three terms on the right-hand side of the equation represent the observations, spatial constraints, and temporal constraints, respectively, the details of which depend on the a priori assumptions of the model; \(W_{0}\) is a weight parameter; \(G\) is the kernel function matrix; \(m\) is the source model; and \(f(x,y,t)\) denotes the three-dimensional space-time model, where \(x\), \(y\), and \(t\) represent the latitude, longitude, and time, respectively.  
The spatial constraint term \(\Phi_{s}(m)\) can be expressed as  
\[\Phi_{s}(m)=\iint\left[W_{1}\left(\frac{\partial m}{\partial x}\right)^{2}+W_ {2}\left(\frac{\partial m}{\partial x\partial y}\right)^{2}+W_{3}\left(\frac{ \partial m}{\partial y}\right)^{2}\right]dxdy \tag{2}\]  
\(W_{1}\), \(W_{2}\), and \(W_{3}\) are the weight coefficients (hyperparameters) related to the spatial smoothness of the model in each direction, whose physical implication is that the source model to be solved is expected to have second-order smoothness in horizontal directions of \(X\), \(XY\), and \(Y\).  
The temporal constraint term \(\varPhi_{T}(m)\) can be expressed as  
\[\varPhi_{T}(m)=\iint W_{4}\bigg{(}\frac{\partial m}{\partial t}\bigg{)}^{2}dt \tag{3}\]  
\(W_{4}\) is the weight coefficient (hyperparameter) related to the temporal smoothness of the model, whose physical implication is that the source model to be solved is expected to have second-order smoothness on the temporal scale.  
During the inversion, it is necessary to determine all the hyperparameters in the equation. Therefore, the \(ABIC\) criterion needs to be introduced to optimize the hyperparameters [48]. It is defined as follows:  
\[ABIC=-2\text{log}\,L+2P \tag{4}\]  
where \(L\) is the likelihood function of the model, and \(P\) is the number of hyperparameters. The criterion can be minimized using multi-parameter and nonlinear optimization methods, such as [PERSON]'s method or the simplex method.