Gravity-field recovery of the Earth using reconstruction of spherical harmonic coefficients up to specified degree and order requires proper data sampling based on Shannon-Nyquist rate. Since, these coefficients are globally significant, the sampling must be done uniformly on the Earth, which it takes much time and expense to collect and process data. Many studies have been done in the field of sampling analysis of spherical harmonics [1,2]. Sneeuw [2] showed a lack of Nyquist sampling rate can cause aliasing of second type in gravity-field modeling. Recently based on Compressive Sensing (CS) theorem the sampling rate can be substantially reduced and a signal can be approximated in sparse sense with fewer sampled data that has main role in reconstruction. In this case, the desired signal can be reconstructed, using only some base functions, which are most strongly correlated with the problem. Therefore, based on this strategy, the base functions posing the best solution to the problem will be selected and the sampling rate for regional gravity field modeling will be decreased significantly. When we say a signal is m-sparse, it means that there are at most  nonzero components in the signal. In this case, only m coefficients of the signal have large magnitude, and others are zero, or have very small values. Here, the desired signal can be reconstructed with its large components without loss of more information. The zero-norm of a vector which is defined as , specifies the sparsity-level of a signal. Sparse approximation has been discussed in many studies [4,5,6,7,8,9]. The basic idea proposed by Mallat and Zhang [4] is called matching pursuit (MP), which is an iterative sparse approximation method to reconstruct a signal under specified conditions by replacing a complex sparse problem with a simple optimized solution. Pati et al. [6] modified this algorithm into orthogonal matching pursuit (OMP), which is used for non-orthogonal dictionaries and converges faster than MP. The regularized orthogonal matching pursuit (ROMP) algorithm popularized by Needell and Vershynin [8] is an iterative sparse approximation method where at each iteration m nonzero components of unknown parameters that most closely resemble the properties of the desired signal are selected. Needell and Tropp [9] refined the ROMP algorithm with compressive sampling matching pursuit (CoSaMP), which identifies locations of the large energy of a signal at each iteration. All these algorithms try to find column vectors in the design matrix that most strongly correlate with the desired signal. It is also assumed that the design matrix is well-posed and prior knowledge of the sparsity-level of the signal is clear. Usually, in practical application an ill-posed problem may be encountered, also the sparsity-level of the signal is not exactly clear, which make it difficult to use conventional iterative methods of CS. In this paper we present a new dynamic algorithm called Stabilized Orthogonal Matching Pursuit (SOMP) for gravity-field recovery of the earth using sparse approximation of geopotential spherical harmonic coefficients, which is compatible with the ill-posed problem and can determine the sparsity-level of the signal, properly. Numerical result of the calculated spherical harmonics coefficients up to degree and order 36 shows that the algorithm is able to reconstruct the Earth's gravity-field with precision in mind the number of samples is 50% lower than the Nyquist rate.