We have developed a multiscale approach for solving 2D and 3D nonlinear inverse problems of gravity data in estimating the basement topography. The inversion is carried out in two stages in which the long-wavelength features of the basement are first estimated from smoothed gravity data via a stochastic optimization algorithm. The solution of this stage is used as the starting model for a deterministic optimization algorithm to reconstruct the short-wavelength features from the full-spectrum gravity data. The forward problem is capable of handling lateral and vertical variations in the density of sediments. Two cases are considered regarding prior knowledge about the density: (1) The density contrast between sediments at the surface and the underlying basement and its vertical variations are a priori known, and (2) only the density contrast at the surface is known with its vertical gradient to be recovered in the inversion. In the former case, the unknowns of the problem are the depths, whereas in the latter case, they are the depths and density gradients defined individually for each prism. Therefore, the inverse problem is ill-posed and has many local minima. The stochastic optimization algorithm uses a random initial model and estimates a coarse model of the basement topography. By repeating the stochastic inversion, an ensemble of solutions is formed defining an equivalent domain in the model space supposed to be within the neighborhood of the global minimum of which several starting solutions are extracted for the secondary deterministic inversion. The presented methodology has been tested successfully in converging to the global minima in 2D and 3D cases with 50 and 2352 total number of prisms, respectively. Finally, the inversion algorithm is used to calculate the thickness of the sediments in the South Caspian Basin using the EIGEN-6c4 global gravity model.