The inverse problem of recovering the Earth's density distribution from satellite data of the first or second derivative of the gravitational potential at orbit height is discussed. This problem is exponentially ill-posed. In this paper a multiscale regularization technique using scaling functions and wavelets constructed for the corresponding integro-differential equations is introduced and its numerical applications are discussed. In the numerical part the second radial derivative of the gravitational potential at 200 km orbit height is calculated on a point grid out of the NASA/GSFC/NIMA Earth Geopotential Model (EGM96). Those simulated derived data out of SGG satellite measurements are taken for convolutions with the introduced scaling functions yielding a multiresolution analysis of harmonic density variations in the Earth's crust.