This article describes a robust simultaneous joint inversion scheme for the interpretation of gravity and self-potential data measured along the profile. The developed scheme jointly inverts the two data sets of the causative targets by some geometrically simple idealized bodies (the so-called approximative/interpretive models) in the restricted class of spheres and cylinders. It simultaneously recovers the characteristic inverse parameters of all interpretive bodies (that is, the depth <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$z$ </tex-math></inline-formula> , amplitude coefficient <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A$ </tex-math></inline-formula> , electric dipole moment <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$K$ </tex-math></inline-formula> , and polarization angle <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\theta $ </tex-math></inline-formula> of each body). It employs the Gauss–Newton (GN) method in the space of the logarithmed and nonlogarithmed model parameters ( <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\log (|A|)$ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\log (K)$ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\theta $ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\log (z)$ </tex-math></inline-formula> of each body) (rather than in the space of the model parameters themselves ( <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$A$ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$K$ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$\theta $ </tex-math></inline-formula> , <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$z$ </tex-math></inline-formula> )) to maintain the convergence. The scheme has been successfully validated on a number of noise-free numerical examples, after which the accuracy and stability of the scheme have been carefully assessed on various noisy data. The influence of the scaling factors of the objective functional subjected to minimization on the convergence of the GN method and on the nonuniqueness of the approximative solution that describes and resembles the underlying buried targets has been investigated. It has been found that the developed scheme is very robust and capable of extracting accurate and useful information that is of some significance in mineral exploration. Finally, the scheme has been applied to a real data example from Germany. Careful analysis of this case study suggests new results that are of some value in mining geophysics.