Abstract: The polar gaps in satellite geodesy lead to poor estimability of low order spherical harmonic coefficients of the Earth’s gravity field model. The Slepian localized spectral analysis method is introduced in this paper and its mathematical characteristics are analyzed. Grünbaum operator is introduced into the solution process in order to improve the computational efficiency and stability. Then, the practical formulas of Slepian method for recovering the gravity field with satellite gravity gradiometry data are educed and numerical analysis are done. The results indicate that, Slepian function is orthogonal in both the sphere and the belt, and its spectral energy distribution characteristic is in accord with the orbit characteristic. Also, the accuracy of Slepian method in Slepian domain is better than that of Legendre solution in spherical harmonic domain. However, the contribution of Slepian method to the improvement of geoid is limited.