The gravity field modeling can be performed in global or local scales utilizing satellite, airborne, terrestrial gravity data or a combination of these observations. One of the common methods in gravity field approximation is to use Spherical Harmonic expansion. Due to the global characteristics of the Spherical Harmonic base functions, a small local signal variation can change all the coefficients in the expansion, and therefore, they are not a suitable choice for regional applications. In order to overcome this problem, there are several groups of base functions, including Radial Basis Functions (RBFs). In gravity field modeling using RBFs, the disturbing potential is represented by a linear combination of an infinite set of RBFs. Hence, any linear functional of the disturbing potential, such as gravity anomaly or gravity disturbance, can be also expressed based on the RBF expansion. Thus, measurable quantities of the Earth's gravity field can be utilized in order to determine the RBF parameters in gravimetric modeling. In this study, system of observation equations is set based on the Radial Multi-Poles of order 2 and free-air gravity anomalies and unknown parameters, including RBF centers, RBF bandwidths (or depths) and scaling coefficients, are determined using a least-squares method. In fact, the Levenberg-Marquardt algorithm is applied as a non-linear regularization method to simultaneously optimize all the RBF parameters. In order to enhance the numerical efficiency of this algorithm, a novel scheme is proposed to initialize and update the regularization parameter. Finally, numerical results obtained from the modified Levenberg-Marquardt algorithm are compared with the ones obtained from the simple form of this algorithm. Applying the proposed modifications to this algorithm, the unknown parameters are determined in a fast procedure and with a significant reduction in the number of iterations. Moreover, these modifications can increase the probability of convergence of the solution to the global minimum.