## 3 Simulation Experiments
### Determining the Earth's External Gravitational Field  
The GP distribution \(\hat{\phi}_{T}(x,y,z)\) over the TS sphere enables us to estimate the gravitational field outside solid Earth using the least-squares approach. The Earth's GP, \(V\), at an exterior point \((r,\theta,\lambda)\) can be expanded into spherical harmonics [2]:  
\[V(r,\theta,\lambda)=\frac{GM}{a}\sum_{n=0}^{N_{max}}\sum_{m=0}^{n}\left(\frac{ a}{r}\right)^{n+1}\left(\overline{C}_{nm}\cos m\lambda+\overline{S}_{nm}\sin m \lambda\right)\overline{P}_{nm}(\cos\theta). \tag{12}\]  
In Equation (12), \(r\) denotes the geocentric radius, while \((r,\theta,\lambda)\) symbolizes a three-dimensional position in the Earth-centered, Earth-fixed (ECEF) reference frame. The variables \(\theta\) and \(\lambda\) represent spherical co-latitude and longitude, respectively. \(GM\) stands for the geocentric gravitational constant, and \(a\) is the semi-major axis of the reference ellipsoid. The coefficients \(\overline{C}_{nm}\) and \(\overline{S}_{nm}\) are the fully normalized GP coefficients describing the Earth's external gravitational field. \(\overline{P}_{nm}\) refers to the fully normalized associated  
\begin{table}
\begin{tabular}{l l} \hline
**Influence Factor** & **(Residual) Error Magnitude in \(\Delta f\overline{f}_{G}\)** \\ \hline Ionospheric correction residual & \(e_{ion}\sim 10^{-19}\) \\ Tidal correction residual & \(e_{tide}\sim 10^{-18}\) \\ Position error & \(e_{pos}<10^{-18}\) \\ Velocity error & \(e_{vel}\sim 10^{-17}\) \\ Gravitational potential error & \(e_{pot}<10^{-18}\) \\ Clock error & \(e_{clk}\sim 10^{-13},10^{-15},10^{-17}\) \\ \hline \end{tabular}
\end{table}
Table 2: Magnitudes of different error sources in determining the gravitational potential difference between the target satellite (TS) and a geosynchronous equatorial orbit satellite (GS). They are transformed to relative frequency.  
Legendre functions of degree \(n\) and order \(m\). Lastly, \(N_{max}\) indicates the maximum degree of the harmonic expansion.  
In the context of the linear observation as per Equation (12), the statistical and functional models that enable gravitational field recovery from GP distribution observations are established using a standard Gauss-Markov model [56]:  
\[\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{\epsilon},\quad E\{\mathbf{y}\}=\mathbf{A}\mathbf{x},\quad D\{\mathbf{y }\}=\sigma_{0}^{2}\mathbf{Q}=\sigma_{0}^{2}\mathbf{P}^{-1}, \tag{13}\]  
In this formulation, \(\mathbf{y}\) denotes the vector of GP observations, while \(\mathbf{A}\) is the design matrix. The vector \(\mathbf{x}\) symbolizes the GP coefficients \(\overline{C}_{nm}\), \(\overline{S}_{nm}\), which are yet to be estimated. The term \(\mathbf{\epsilon}\) corresponds to the vector of observational errors. \(D\{\mathbf{y}\}\) represents the error variance-covariance matrix, \(\mathbf{P}\) is the weight matrix, \(\mathbf{Q}\) signifies the inverse of the weight matrix, and \(\sigma_{0}^{2}\) is the variance component.  
Based on the above-described data processing method, we estimated three recovered Earth Gravitational Models (REGMs) with the degree and order of 200 from GP values distributed over the TS sphere in the three experiment cases. Since the error distribution is unknown in GP observations, we set the weight matrix \(\mathbf{P}\) as a unit matrix. We consider that the noises in GP observations are white noises. The results are illustrated and discussed in Section 4.