## 3 Methodology
Additionally, [PERSON] considered latitude differences and used a one-dimensional Fast Fourier Transform (1D FFT) for global gravity field calculations, making the theoretical framework more rigorous:  
\[\Delta g_{\phi_{p}}(\lambda_{p}) =\frac{\gamma_{0}\Delta\phi\Delta\lambda}{4\pi}\sum\limits_{\phi_ {q}=\phi_{1}}^{\phi_{p}}\sum\limits_{\lambda_{q}=\lambda_{1}}^{\lambda_{q}}H ^{\prime}\bigl{(}\Delta\lambda_{qp}\bigr{)}\times\bigl{(}\xi_{\cos}\cos\alpha _{qp}+\eta_{\cos}\sin\alpha_{qp}\bigr{)} \tag{8}\] \[=\frac{\gamma_{0}\Delta\phi\Delta\lambda}{4\pi}F_{1}^{-1}\times \Biggl{\{}\sum\limits_{\phi_{q}=\phi_{1}}^{\phi_{p}}\bigl{[}F_{1}\left(H^{ \prime}\bigl{(}\Delta\lambda_{qp}\bigr{)}\cos\alpha_{qp}\right)F_{1}(\xi_{\cos })\] \[+F_{1}\bigl{(}H^{\prime}\bigl{(}\Delta\lambda_{qp}\bigr{)}\sin \alpha_{qp}\bigr{)}F_{1}(\eta_{\cos})\bigr{\}}\]  
where \(\xi_{\cos}=\xi\cos\phi\), \(\eta_{\cos}=\eta\cos\phi\), and \(\Delta\phi\) and \(\Delta\lambda\) represent the grid spacing in the latitude and longitude directions, respectively. However, when the spherical distance \(\psi_{pq}=0\), the kernel function will be singular. Therefore, the innermost zone effect must be considered, which can be calculated using Equation (9):  
\[\eta(i,j)=-\frac{1}{Rcos\varphi_{i}}\frac{\left(N_{j+1}-N_{j}\right)}{\left( \lambda_{j+1}-\lambda_{j}\right)} \tag{3}\]0  
\[\eta(i,j)=-\frac{1}{Rcos\varphi_{i}}\frac{\left(N_{j+1}-N_{j}\right)}{\left( \lambda_{j+1}-\lambda_{j}\right)} \tag{3}\]1  
where \(s_{0}\) is the innermost zone, \(\xi_{y}=\frac{\partial\xi}{\partial y}\), and \(\eta_{x}=\frac{\partial\eta}{\partial x}\). It should be noted that the grid DOV calculated from the simulated data is only available for the marine region. However, since the gravity of the land area also affects the long-wave component of the marine gravity field, the DOV resolved by the EGM2008 model is used to fill the missing values in the land area before calculating the marine gravity anomaly. It can be seen that we did not use the remove-restore method to restore the long-wave component of the marine gravity field for the residual gravity anomalies [19], but instead used the total vertical deflection to invert the gravity anomalies. While both methods aim to recover the long-wave component of the gravity field, the remove-restore method is unable to effectively suppress noise in the simulation scenarios. On the other hand, using total vertical deflection for gravity anomaly inversion allows for a better assessment of the potential of utilizing spaceborne GNSS-R interferometric altimetry in simulation gravity field inversion. Lastly, the gravity anomaly model EGM2008, released by the National Geospatial-Intelligence Agency (NGA), is employed as the true gravity anomaly values. These reference values are then compared with the gravity anomaly results calculated using the simulated data in a month. The EGM2008 model provides up to 2159 degrees and 2190 orders of spherical harmonic coefficients, which can be used to calculate gravity anomaly results with a spatial resolution of up to 5\({}^{\prime}\times 5^{\prime}\) using spherical harmonic functions, as shown in \[\eta(i,j)=-\frac{1}{Rcos\varphi_{i}}\frac{\left(N_{j+1}-N_{j}\right)}{\left( \lambda_{j+1}-\lambda_{j}\right)} \tag{3}\]2 By using different maximum spherical harmonic coefficients for expansion orders,we can obtain gravity anomaly results with different grid resolutions and compare them with the simulated gravity anomaly results at corresponding resolutions. The relationship between the order of maximum spherical harmonic coefficients and spatial resolution can be described by Equation (11):  
\[\eta(i,j)=-\frac{1}{Rcos\varphi_{i}}\frac{\left(N_{j+1}-N_{j}\right)}{\left( \lambda_{j+1}-\lambda_{j}\right)} \tag{3}\]3  
where \(n\) is the maximum degree of the spherical harmonic coefficients used, and \(\Omega\) is the corresponding resolution of the model.