Search

We have developed an iterative, efficient, and stable 3D magnetic inversion method to estimate the depth to the basement of a sedimentary basin, the volume of which is approximated by a grid of rectangular prisms with tops at an arbitrary surface and thicknesses estimated iteratively. To compute the thicknesses corrections at the kth iteration, we first compute an approximate sparse sensitivity matrix in which the nonnull elements of a given line are obtained by differentiating the fitting function at the observation position (which occupies the center of a small moving data window) relative to the thicknesses of the prisms situated inside the window. By moving the data window around, all nonnull elements of all lines are computed. This sparse matrix is used to calculate an approximate Gauss-Newton gradient and an approximate diagonal Gauss-Newton Hessian matrix. Finally, by dividing each element of the gradient vector by the corresponding diagonal element of the approximate Hessian matrix, we obtain the corrections of all parameters at a given iteration. The solution is stable because of (1) the small parameter corrections inherent to the method, combined with the initial estimate at the surface and (2) the application of a moving average operator to the solution. Any magnetization orientation, except horizontal (and parallel to any vertical side of the prisms), yields good results. In the case of magnetization inclinations of approximately 45°, the rectangular prisms may induce spurious ripples on the surface of the estimated relief, which can be eliminated by simultaneously reducing the horizontal dimensions of the prisms and applying an appropriate moving average to the solution. Application of our method to real data in a rift area discloses a confined basin with steep and asymmetric borders, confirming the practical utility of the method.

We have developed a new and robust method (in the sense of it being applicable to a wide range of situations) to estimate the regularization parameter [Formula: see text] in a regularized inverse problem. For each tentative value of [Formula: see text], we perturb the observations with [Formula: see text] sequences of pseudorandom noise and we track down the instability effect on the solutions. Then, we define a quantitative measure [Formula: see text] of the solution instability consisting of the largest value among the Chebyshev norms of the vectors obtained by the differences between all pairs of the perturbed solutions. Despite being quantitative, [Formula: see text] cannot be used directly to estimate the best value of [Formula: see text] (the smallest value that stabilizes the solution) because, in practice, instability may depend on the particular and specific interests of the interpreter. Then, we determine that the interpreter, at each iteration of a bisection method, visually compares, in the ([Formula: see text], [Formula: see text], [Formula: see text]) space, the pair [Formula: see text] and [Formula: see text] of the solutions most distant from each other and associated with the current [Formula: see text]. From this comparison, the interpreter decides if the current [Formula: see text] produces stable solutions. Because the bisection method can be applied only to monotonic functions (or segments of monotonic functions) and because [Formula: see text] has a theoretical monotonic behavior that can be corrupted, in practice by a poor experiment design, the set of values of [Formula: see text] can be used as a quality control of the experiments in the proposed bisection method to estimate the best value of [Formula: see text]. Because the premises necessary to apply the proposed method are very weak, the method is robust in the sense of having broad applicability. We have determined part of this potential by applying the proposed method to gravity, seismic, and magnetotelluric synthetic data, using two different interpretation models and different types of pseudorandom noise.

We have developed a novel approach to compute, in an efficient and versatile way, the gravity anomaly produced by an arbitrary, discrete 3D distribution of density contrast. The method allows adjustable precision and is particularly suited for the interpretation of sedimentary basins. Because the gravity field decays with the square of the distance, we use a discrete Green’s operator that may be made much smaller than the whole study area. For irregularly positioned observations, this discrete Green’s operator must be computed just at the first iteration, and because each of its horizontal layers presents a center of symmetry, only one-eighth of its total elements need to be calculated. Furthermore, for gridded data on a plane, this operator develops translation symmetry for the whole region of interest and has to be computed just once for a single arbitrary observation position. The gravity anomaly is obtained as the product of this small operator by any arbitrary density contrast distribution in a convolution-like operation. We use the proposed modeling to estimate the basement relief of a [Formula: see text] basin with density contrast varying along [Formula: see text] only using a smaller Green’s operator at all iterations. The median of the absolute differences between relief estimates, using the small and a large operator (the latter covering the whole basin) has been approximately 9 m for a 3.6 km deep basin. We also successfully inverted the anomaly of a simulated basin with density contrast varying laterally and vertically, and a real anomaly produced by a steeply dipping basement. The proposed modeling is very fast. For 10,000 observations gridded on a plane, the inversion using the proposed approach for irregularly spaced data is two orders of magnitude faster than using an analytically derived fitting, and this ratio increases enormously with the number of observations.

We have developed a gravity inversion method to estimate a discontinuous basement relief based on an extended version of Bott’s method that allows variable density contrasts between sediments and basement, optimizes the modulus of the solution correction at each iteration, and provides for solution stabilization. Initially, we obtain a linear approximation stabilized by the total variation functional that correctly maps the horizontal positions of the existing high-angle faults but produced poor estimates of the basin depths. Subsequent iterations update the depth estimates toward the correct values, at the same time preserving the correct fault horizontal positions. Additionally, we stabilize each solution correction by the smoothness constraint without inverting any matrix. The method was substantially more efficient than the nonlinear method, which solves a system of linear equations by the conjugate gradient method at each iteration. For 3000 parameters, it is almost four times faster than the conjugate gradient, and this ratio increased with the number of parameters.

We have combined modeling and gravity inversion to estimate an overall smooth but locally discontinuous basement relief of a sedimentary basin. Through an initial global smoothness solution, the method indicates fault positions that may be accepted or modified. Then, the interpreter specifies the fault dip and a tentative fault slip, so the fault plane is fixed on the estimated relief. The segments of the current interpreted relief that are not accepted faults are estimated through smoothness inversion, and the response of the current estimated relief is computed. The slip of the fault being incorporated in the estimated relief is then modified interactively until the solution presents no side lobes about the fault extremes. The process is repeated for all faults. The iteration stops when the tentative solution becomes compatible with the geologic knowledge of the area and produces an acceptable data fit. The method was applied to data produced by a simulated graben defined by step faults dipping 60º and by inclined and arcuate terraces. The solution virtually coincided with the true source. The method’s ability in testing geologic hypotheses was demonstrated in gravity profiles across the Büyük Menderes Valley in west Turkey and the San Jacinto graben in southern California.

Gravity interpretation of large sedimentary basins requires a large number of observations and of parameters estimations, so the development of very efficient gravity inversion methods applicable to such environments is of utmost importance. We found that the highly efficient Bott’s method can be cast in the framework of Gauss-Newton’s method for minimizing nonlinear functions. Also, we extended Bott’s method to (1) be applicable to sedimentary basins with density contrast between the sediments and the basement decreasing with depth according to different laws, (2) optimize the modulus of the parameter correction vector, (3) stabilize the solution by imposing that it be smooth, and (4) quit the iteration according to an objective and effective criterion. To keep up the efficiency of Bott’s method, stable solutions were obtained by applying a moving average to the usually unstable solution obtained at the last iteration. Solutions as efficacious as those obtained by the usual nonlinear method were produced by the present extension at a much smaller computation time. Nonlinear methods, with different implementations for solving a linear system of equations, required between one and two orders of magnitude more time to produce similar solutions using between 2000 and 3000 parameters and observations, for example.

We found that minimum [Formula: see text]-norm and smoothness-constrained continuous solutions of the linear inverse problem of potential field data are harmonic and biharmonic, respectively. In the case of a discrete distribution, the minimum [Formula: see text]-norm and smoothness-constrained solutions become biased toward being harmonic or biharmonic, respectively. As a result, the estimated discrete distribution of density or magnetization contrast tends to be smooth and to satisfy the maximum principle, which forces the solution maxima and minima to lie on any boundary of the discretized region. The above findings were illustrated with 2D numerical examples. The harmonic or biharmonic bias is brought forth when the strengths of the minimum [Formula: see text]-norm or the smoothness constraint are enhanced (relative to all other constraints) by approximating the continuous case (a large number of discretizing cells relative to the number of independent observations) and/or by using a regularizing parameter. We discovered that, by inspecting the rearranged normal equations, it is possible to qualify three different possibilities of designing nonharmonic estimators. Then we found that all three possibilities have in fact already been implemented in the literature, reinforcing, in this way, the validity of the theoretical demonstration.

We have developed a new method that drastically reduces the number of the source location estimates in Euler deconvolution to only one per anomaly. Our method employs the analytical estimators of the base level and of the horizontal and vertical source positions in Euler deconvolution as a function of the [Formula: see text]- and [Formula: see text]-coordinates of the observations. By assuming any tentative structural index (defining the geometry of the sources), our method automatically locates plateaus, on the maps of the horizontal coordinate estimates, indicating consistent estimates that are very close to the true corresponding coordinates. These plateaus are located in the neighborhood of the highest values of the anomaly and show a contrasting behavior with those estimates that form inclined planes at the anomaly borders. The plateaus are automatically located on the maps of the horizontal coordinate estimates by fitting a first-degree polynomial to these estimates in a moving-window scheme spanning all estimates. The positions where the angular coefficient estimates are closest to zero identify the plateaus of the horizontal coordinate estimates. The sample means of these horizontal coordinate estimates are the best horizontal location estimates. After mapping each plateau, our method takes as the best structural index the one that yields the minimum correlation between the total-field anomaly and the estimated base level over each plateau. By using the estimated structural index for each plateau, our approach extracts the vertical coordinate estimates over the corresponding plateau. The sample means of these estimates are the best depth location estimates in our method. When applied to synthetic data, our method yielded good results if the bodies produce weak- and mid-interfering anomalies. A test on real data over intrusions in the Goiás Alkaline Province, Brazil, retrieved sphere-like sources suggesting 3D bodies.

An aeromagnetic data filter was tested as a tool in enhancing the magnetic responses of deep sources associated with marginal basins. Our approach consists in amplifying the responses of magnetic sources which extend through the continental slope and continental rise. This amplification compensates for the magnetic response attenuation over the continental slope as a result of (1) the seaward increase of the distance between the source and the magnetometer and (2) the seaward decrease in the volume of sources within the continental crust. To amplify the source response, we multiply the total-field anomaly at a given location by the exponential of the bathymetric depth at the same point. Tests on synthetic and real data showed that the filter was able to recover the locations of linear magnetic responses not clearly seen in the total-field-anomaly map. When compared with the analytic signal amplitude and the total horizontal derivative filters, the proposed filter proved to be complementary in enhancing lineaments with different orientations. When applied in conjunction with these filters, the proposed filter produced a synergistic effect, revealing lineaments that were not present in individual applications of the filters. The application of the filter to an offshore Brazilian margin comprising the Camamu-Almada and Jequitinhonha marginal basins highlighted northwest–southeast and northeast–southwest lineaments not clearly seen in the total-field data. The northeast–southwest lineaments were interpreted as related to fracture zones that may cut across the continental margin, controlling basement structures and fluid migration paths.

We have developed an iterative scheme for inverting gravity data produced by salt bodies with density contrasts relative to the sediments varying from positive to negative, crossing, in this way, the nil zone. Our inversion method estimates a 3D density-contrast distribution, through a piecewise constant function defined on a user-specified grid of cells. It consists of two nested iterative loops. The outer loop uses an adaptive learning strategy that starts with a coarse grid of cells, a set of first-guess geometric elements (axes and points) and the corresponding assigned density contrasts. From the second iteration on, this strategy refines the grid and automatically creates a new set of geometric elements (points only) and associated density contrasts. Each geometric element operates as the first-guess skeletal outline of a section of the salt body to be imaged. The inner loop estimates the 3D density-contrast distribution for the grid of cells and for the set of geometric elements defined in the outer loop. The outer loop allows for easy incorporation of prior geologic information about the lithologic units and automatic evolution of the prior information. The inner loop forces the estimated density contrast of each cell to be close either to a null or to a non-null prespecified value. The iteration stops when the geometries of the estimated salt bodies are invariant along successive iterations. We apply our method to synthetic gravity data produced by a homogeneous salt body embedded in heterogeneous sediments. We tested two geologic hypotheses about the real gravity data from Galveston Island salt dome, USA. In the first, the estimated salt body attains a maximum bottom depth of 5 km, whereas in the second hypothesis, it is shallower and discloses an overhang. Both solutions fit the data and are feasible geologically, so both hypotheses are acceptable.

We applied the mathematical basis of the total variation (TV) regularization to analyze the physicogeologic meaning of the TV method and compared it with previous gravity inversion methods (weighted smoothness and entropic Regularization) to estimate discontinuous basements. In the second part, we analyze the physicogeologic meaning of the TV method and compare it with previous gravity inversion methods (weighted smoothness and entropic regularization) to estimate discontinuous basements. Presenting a mathematical review of these methods, we show that minimizing the TV stabilizing function favors discontinuous solutions because a smooth solution, to honor the data, must oscillate, and the presence of these oscillations increases the value of the TV stabilizing function. These three methods are applied to synthetic data produced by a simulated 2D graben bordered by step faults. TV regularization and weighted smoothness are also applied to the real anomaly of Steptoe Valley, Nevada, U.S.A. In all applications, the three methods perform similarly. TV regularization, however, has the advantage, compared with weighted smoothness, of requiring no a priori information about the maximum depth of the basin. As compared with entropic regularization, TV regularization is much simpler to use because it requires, in general, the tuning of just one regularization parameter.

We have developed an inversion approach that estimates the basement relief of a fault-bounded sedimentary basin. The sedimentary pack is approximated by a grid of 3D or 2D vertical prisms juxtaposed in the horizontal directions of a right-handed coordinate system. The prisms’ thicknesses represent the depths to the basement and are the parameters to be estimated from the gravity data. To obtain depth-to-basement estimates, we introduce the total variation (TV) regularization as a stabilizing function. This approach lets us estimate a nonsmooth basement relief because it does not penalize sharp features of the solution. We have deduced a compact matrix form of the gradient vector and the Hessian matrix of the approximation to the TV function that allows a regularized Gauss-Newton minimization approach. Because the Hessian matrix of the approximation to the TV function is ill conditioned, we have modified this Hessian matrix to improve its condition and to accelerate the convergence of the Gauss-Newton algorithm. Tests conducted with synthetic data show that the inversion method can delineate discontinuous basements presenting large slips or sequences of small-slip step faults. Tests on field data from the Almada Basin, Brazil, and from the San Jacinto Graben, California, U.S.A., confirm the potential of the method in detecting and locating in-depth normal faults in the basement relief of a sedimentary basin.

We have developed a gravity interpretation method for estimating the discontinuous basement relief of a sedimentary basin. The density contrast between the basement and the sediments is assumed to be known, and it could be either constant or vary monotonically with depth. The interpretation model consists of a set of vertical, juxtaposed prisms, whose thicknesses are the parameters to be estimated. We used the entropic regularization that combines the minimization of the first-order entropy measure with the maximization of the zeroth-order entropy measure of the solution vector. We validated the method by applying it to synthetic data produced by a simulated basin bordered by high-angle step faults; we obtained a good definition of the relief, particularly of the discontinuities. We also applied the method to a profile across the Büyük Menderes Valley in West Turkey and obtained a solution exhibiting a gravity fault with large slip on the northern border of the valley. When applied to the interpretation of a discontinuous basement relief, the method has a better performance than the global smoothness method. It is comparable to the weighted smoothness method, but it does not require the a priori knowledge about the maximum basin depth.

We have developed a gravity-inversion method for simultaneously estimating the 3D basement relief of a sedimentary basin and the parameters defining a presumed parabolic decay of the density contrast with depth in a sedimentary pack, assuming prior knowledge about the basement depth at a few points. The sedimentary pack is approximated by a grid of 3D vertical prisms juxtaposed in both horizontal directions of a right-handed coordinate system. The prisms’ thicknesses represent the depths to the basement and are the parameters to be estimated from the gravity data. To estimate the parameters defining the parabolic decay of the density contrast with depth and to produce stable depth-to-basement estimates, we imposed smoothness on the basement depths and proximity between estimated and known depths at boreholes. We applied our method to synthetic data from a simulated complex 3D basement relief with two sedimentary sections having distinct parabolic laws describing the density-contrast variation with depth. The results provide good estimates of the true parameters of the parabolic law of density-contrast decay with depth and of the basement relief. Inverting the gravity data from the onshore and part of the shallow offshore Almada Basin on Brazil’s northeastern coast shows good correlation with known structural features.

A new apparent-magnetization mapping method on the horizontal plane combines minimization of first-order entropy with maximization of zeroth-order entropy of the estimated magnetization. The interpretation model is a grid of vertical, juxtaposed prisms in both horizontal directions. To estimate the magnetization of the prisms, assume that the top and bottom of the magnetic sources are horizontal. Minimization of the first-order entropy favors solutions with sharp borders, and the maximization of zeroth-order entropy prevents the tendency of the estimated source to become a single prism with large magnetization. Thus, a judicious combination of both constraints can lead to solutions characterized by regions with virtually constant magnetizations separated by sharp discontinuities. This is applied to synthetic data from simulated intrusive bodies in sediments that have horizontal tops. By comparing the results with those obtained with the common Tikhonov regularization (smoothness constraint) method, it is shown that both methods produce good and equivalent locations of the central positions of the sources. However, entropic regularization delineates the boundaries of the bodies with greater detail. Both the proposed and the smoothness constraints are applied to real anomaly data over a magnetic skarn in Butte Valley, Nevada, U.S.A. Entropic regularization produced an estimated magnetization distribution with sharper boundaries, smaller volume, and higher apparent magnetization as compared with results produced by incorporating the smoothness constraint.

We have developed a gravity inversion method to estimate a 3D density-contrast distribution producing strongly interfering gravity anomalies. The interpretation model consists of a grid of 3D vertical, juxtaposed prisms in the horizontal and vertical directions. Iteratively, our approach estimates the 3D density-contrast distribution that fits the observed anomaly within the measurement errors and favors compact gravity sources closest to prespecified geometric elements such as axes and points. This method retrieves the geometry of multiple gravity sources whose density contrasts (positive and negative values) are prescribed by the interpreter through the geometric element. At the first iteration, we set an initial interpretation model and specify the first-guess geometric elements and their target density contrasts. Each geometric element operates as the first-guess skeletal outline of the entire homogeneous gravity source or any of its homogeneous parts to be reconstructed. From the second iteration on, our method automatically redefines a new set of geometric elements, the associated target density contrasts, and a new interpretation model whose number of prisms increases with the iteration. The iteration stops when the geometries of the estimated sources are invariant along successive iterations. Tests on synthetic data from geometrically complex bodies and on field data collected over a mafic-ultramafic body and a volcanogenic sedimentary sequence located in the Tocantins Province, Brazil, confirmed the potential of our method in producing a sharp image of multiple and adjacent bodies.

Wave propagation in a layered medium when the wavelength is much greater than each layer thickness (low frequency) produces a response equivalent to that of wave propagation in an equivalent single-layer medium. This equivalent medium is transversely isotropic with symmetry about a vertical axis (VTI), and the elastic parameters are computed with the Backus averaging technique. Conversely, when the wavelength is comparable to each layer thickness (high frequency), the directional dependence of the phase velocity in the transmission response also can be simulated by replacing the layered medium with a single homogeneous medium with properties derived from a time average. It then can be treated approximately as a VTI medium. To compute the medium parameters, a method based on fitting the traveltime parameters is used. We investigated the relationship between Thomsen’s anisotropic parameters [Formula: see text] and [Formula: see text] computed for the equivalent medium in the low-frequency limit and for the homogenized medium in the high-frequency limit. In our experiments, we used a medium in which layers of only two isotropic materials alternate repeatedly. For the high-frequency limit, we obtained solutions for PP- and SS-wave propagation.

We present a new apparent-density mapping method on the horizontal plane that combines the minimization of the first-order entropy with the maximization of the zeroth-order entropy of the estimated density contrasts. The interpretation model consists of a grid of vertical, juxtaposed prisms in both horizontal directions. We assume that the top and the bottom of the gravity sources are flat and horizontal and estimate the prisms’ density contrasts. The minimization of the first-order entropy favors solutions presenting sharp borders, and the maximization of the zeroth-order entropy prevents the tendency of the source estimate to become a single prism. Thus, a judicious combination of both constraints may lead to solutions characterized by regions with virtually constant estimated density contrasts separated by sharp discontinuities. We apply our method to synthetic data from simulated intrusive bodies in sediments that present flat and horizontal tops. By comparing our results with those obtained with the smoothness constraint, we show that both methods produce good and equivalent locations of the sources’ central positions. However, the entropic regularization delineates the boundaries of the bodies with greater resolution, even in the case of 100-m-wide bodies separated by a distance as small as [Formula: see text]. Both the proposed and the global smoothness constraints are applied to real anomalies from the eastern Alps and from the Matsitama intrusive complex, northeastern Botswana. In the first case, the entropic regularization delineates two sources, with a horizontal and nearly flat top being consistent with the known geologic information. In the second case, both constraints produce virtually the same estimate, indicating, in agreement with results of synthetic tests, that the tops of the sources are neither flat nor horizontal.

We demonstrate the potential of gravity data to detect and to locate in-depth subtle normal faults in the basement relief of a sedimentary basin. This demonstration is accomplished by inverting the gravity data with the constraint that the estimated basement relief presents local abrupt faults and is smooth elsewhere. We inverted the gravity data from the onshore Almada Basin in northeastern Brazil, and we mapped several normal faults whose locations and plane geometries were already known from seismic imaging. The inversion method delineated well both the discontinuities with small or large slips and a sequence of step faults. Using synthetic data, we performed a systematic search of normal fault slips versus fault displacement depths to map the fault-detectable region in this space. This mapping helps to assess the ability of gravity inversion to detect normal faults. Mapping shows that normal faults with small [Formula: see text], medium (about [Formula: see text]), and large (about [Formula: see text]) vertical slips can be detected if the maximum midpoint depths of the fault planes are smaller than 1.8, 3.8, and [Formula: see text], respectively.

We present a new semiautomatic gravity interpretation method for estimating a complex interface between two media containing density heterogeneities (referred to as interfering sources) that give rise to a complex and interfering gravity field. The method combines a robust fitting procedure and the constraint that the interface is very smooth near the interfering sources, whose approximate horizontal coordinates are defined by the user. The proposed method differs from the regional-residual separation techniques by using no spectral content assumption about the anomaly produced by the interface to be estimated, i.e., the interface can produce a gravity response containing both low- and high-wavenumber features. As a result, it may be applied to map the relief of a complex interface in a geologic setting containing either shallow or deep-seated interfering sources. Tests conducted with synthetic data show that the method can be of utility in estimating the basement relief of a sedimentary basin in the presence of salt layers and domes or in the presence of mafic intrusions in the basement or in both basement and the sedimentary section. The method was applied to real gravity data from two geologic settings having different kinds of interfering sources and interfaces to be interpreted: (1) the interface between the upper and lower crusts over the Bavali shear zone of southern India and (2) the anorthosite-tonalite interface over the East Bull Lake gabbro-anorthosite complex outcrop in Ontario, Canada.

We present a method to estimate the basement relief as well as the density contrast at the surface and the hyperbolic decaying factor of the density contrast with depth, assuming that the gravity anomaly and the depth to the basement at a few points are known. In both cases, the interpretation model is a set of vertical rectangular 2D prisms whose thicknesses are parameters to be estimated and that represent the depth to the interface separating sediments and basement. The solutions to both problems are stable because of the incorporation of additional prior information about the smoothness of the estimated relief and the depth to the basement at a few locations, presumably provided by boreholes. The method was tested with synthetic gravity anomalies produced by simulated sedimentary basins with smooth relief, providing not only well-resolved estimated relief, but also good estimates for the density contrasts at the surface and for the decaying factors of the density contrast with depth. The method was applied to the Bouguer anomaly from Recôncavo Basin, estimating the surface density contrast and the decaying factor of the density contrast with depth as [Formula: see text] and [Formula: see text], respectively.

We present a method for inverting magnetic data with interfering anomalies produced by multiple complex 2D magnetic sources having arbitrary shapes and known magnetization vectors. Our method is stable and can recover a complex 2D magnetization distribution, leading to a reliable delineation of sectionally homogeneous sources with complex shapes. Our method, although similar to interactive forward modeling, is unique in that it automatically fits the observations and only requires that the interpreter know the outlines of the sources expressed by simple geometric elements such as points and line segments. Each geometric element operates as a skeletal outline of a particular homogeneous section of the magnetic source to be reconstructed. Also, the interpreter can define the geometric elements interactively without worrying about data fitting because data are fit automatically. The examples with synthetic data illustrate the good performance of the method in mapping the complex geometry of magnetic sources. The solution sensitivity to uncertainties in the a priori information shows that to produce good results, the uncertainty on the magnetization intensity of each homogeneous extent of the source should be smaller than 40%. A wrong magnetization vector direction can be detected easily because it often leads to poor data fitting and to estimated sources with abrupt borders. The method is also applied to two sets of real data from the Northwest Ore Body at Iron Mountain Mine, Missouri, and the Hatton-Rockall Basin in the northeast Atlantic Ocean. The estimated magnetization distribution in all tests demonstrates a good correlation of estimated magnetic sources with corresponding known geologic features.

We have developed a new approach for estimating the location and geometry of several density anomalies that give rise to a complex, interfering gravity field. The user interactively defines the assumed outline of the true gravity sources in terms of points and line segments, and the method estimates sources closest to the specified outline to achieve a match between the predicted and observed gravity fields. Each gravity source is assumed to be a homogeneous body with a known density contrast; different density contrasts may be assigned to each source. Tests with synthetic data show that the method can be of use in estimating (1) multiple laterally adjacent and closely situated gravity sources, (2) single gravity sources consisting of several homogeneous compartments with different density contrasts, and (3) two gravity sources with different density contrasts of the same sign, one totally enclosed by the other. The method is also applied to three different sets of field data where the gravity sources belong to the same categories established in the tests with synthetic data. The method produces solutions consistent with the known geologic attributes of the gravity sources, illustrating its potential practicality.

We introduce a new 2D method for inverting potential‐field data with model constraints designed by the interpreter. Our method uses an interpretation model consisting of a source with polygonal cross‐section whose vertices are described by polar coordinates with an origininside the source. With this coordinate system, constraints in an inversion are easier to develop and apply. Our inversion method assumes a known physical property contrast for the source and estimates the radii associated with the polygon vertices for a fixed number of equally spaced angles from 0° to 360°. A wide variety of constraints may be used to stabilize the solutions by introducing information about the source shape. The method recovers stable solutions whose shapes range from almost circular or pear‐shaped to elongated in one or more directions. The convexity constraint applied to the source shape, despite requiring no quantitative information, is more versatile than the other constraints. The convexity constraint efficiently recovers source geometries that are either isometric or elongated in one direction.

We derive the analytical estimators for the horizontal and vertical source positions in 3D Euler deconvolution as a function of the x‐, y‐, and z‐derivatives of the magnetic anomaly within a data window. From these expressions we show that, in the case of noise‐corrupted data, the x‐, y‐, and z‐coordinate estimates computed at the anomaly borders are biased toward the respective horizontal coordinate of the data window center regardless of the true or presumed structural indices and regardless of the magnetization inclination and declination. On the other hand, in the central part of the anomaly, the x‐ and y‐coordinate estimates are very close to the respective source horizontal coordinates regardless of the true or presumed structural indices and regardless of the magnetization inclination and declination. This contrasting behavior of the horizontal coordinate estimates may be used to automatically delineate the region associated with the best solutions. Applying the Euler deconvolution operator inside this region would decrease the dispersion of all position estimates, improving source location precision.

We illustrate the importance of establishing solution uniqueness through mathematical restrictions reflecting a source attribute. We also illustrate the validity and utility of a guideline derived in an accompanying paper for constructing sound gravity inversion methods for the class of sources presenting either homogeneous or depth‐independent density distributions. The two‐part guideline is (1) to introduce a priori information favoring uniqueness, either by assuming that a nonnull density distribution depending only on x and y is confined to the interior of a horizontal slab with known position or by limiting the class of possible solutions to homogeneous, simply connected polygons (or polyhedra) with known density, displaying no fancy shapes and no curling apophyses at their borders, and (2) to introduce information favoring solution stability by estimating only the features of the source which may be resolved by the data. Following the guideline, we apply different methods to gravity data using interpretation models consisting of a grid of cells on the x‐y and x‐z planes. In both cases the estimates are very close to the true synthetic source. The data produced by the distribution varying with x and z are also inverted using the method, which minimizes the norm of the first‐order derivative of the density. This constraint does not reflect a true source attribute but is strong enough to stabilize the solution and to guarantee its uniqueness. Because of the strong bias imposed to the solution, the estimated distribution, although unique and stable, is far from the true source, concentrating most of the anomalous mass at the surface. Finally, we present an alternative method which redistributes the estimated anomalous mass downward. To be effective, this technique requires prior knowledge about the source depth to the top. In addition, the source should not be too small and deep. Although being able to produce good results, this alternative method requires a great dose of the interpreter's art.

To obtain a unique and stable solution to the gravity inverse problem, a priori information reflecting geological attributes of the gravity source must be used. Mathematical conditions to obtain stable solutions are established in Tikhonov's regularization method, where the a priori information is introduced via a stabilizing functional, which may be suitably designed to incorporate some relevant geological information. However, there is no unifying approach establishing general uniqueness conditions for a gravity inverse problem. Rather, there are many theorems, usually establishing just abstract mathematical conditions and making it difficult to devise the type of geological information needed to guarantee a unique solution. In Part I of these companion papers, we show that translating the mathematical uniqueness conditions into geological constraints is an important step not only in establishing the type of geological setting where a particular method may be applied but also in designing new gravity inversion methods. As an example, we analyze three uniqueness theorems in gravimetry restricted to the class of homogeneous bodies with known density and show that the uniqueness conditions established by them are more probably met if the solution is constrained to be a compact body without curled protrusions at their borders. These conditions, together with stabilizing conditions (assuming a simple shape for the source), form a guideline to construct sound gravity inversion methods. A historical review of the gravity interpretation methods shows that several methods implicitly follow this guideline. In Part II we use synthetic examples to illustrate the theoretical results derived in Part I. We also illustrate that the presented guideline is not the only way to design sound inversion methods for the class of homogeneous bodies. We present an alternative approach which produces good results but whose design requires a good dose of the interpreter's art.

Despite being widely used, Euler deconvolution has received little attention in its theoretical aspects. We obtain analytical estimators [Formula: see text] and [Formula: see text] for the horizontal [Formula: see text] and the vertical [Formula: see text] source positions, which are the parameters to be estimated when the structural index is assigned a tentative value. By analyzing the estimators’ properties, we show two things. (1) Scattering of the alternative solutions associated with different data window positions has two components: data noise and the use of a wrong tentative value for the structural index, (2) The differences [Formula: see text] and [Formula: see text], as a function of the data window position, exhibit antisymmetric and symmetric behaviors, respectively, about [Formula: see text] when the magnetization and the geomagnetic field inclinations are 90°, 0°, or 45°. For intermediate inclinations there is a moderate departure from the symmetry properties. This analysis shows that the criterion for determining the structural index as the tentative value producing the smallest solution scattering is theoretically sound but occasionally fails in practice because data noise also contributes to solution scattering. In addition, the lower sensitivity of the averages of estimates of the horizontal position (as compared with the vertical position) occurs because [Formula: see text] displays an antisymmetric behavior with respect to the true horizontal position while the estimates of the vertical position display a symmetric behavior with respect to the same point.

To produce a unique and stable solution in potential‐field interpretation, an inversion method must introduce particular constraints. These constraints will inevitably restrict the type of geological setting where the method may be applied. We present a nonmathematical overview of most stabilizing constraints used in inversion methods. Our purpose is to demonstrate that the inversion results are valuable only if the mathematical stabilizing constraints are translated from the geological setting. We identify five basic types of constraints: (1) lower and upper bounds of parameter estimates; (2) proximity of a parameter estimate to a specified value; (3) proximity between pairs of parameter estimates; (4) concentration of the anomalous source about a geometrical element such as an axis; and (5) source compactness. In practice, if used in isolation, constraints (1), (2), (4), and (5) will not produce geologically meaningful results, regardless of the geological setting of the interpretation area. Constraint (3) may produce geologically meaningful results if the anomalous source has a spatially smooth attribute such as the physical property. We illustrate that constraints 1–4, if used in isolation, cannot delineate the geometry of a simulated sill intruded into a sedimentary basin. The basic constraints may (and should) be combined in inversion to produce geologically meaningful results. We present two examples of effective constraint combination: (1) proximity to a specific value and mass concentration about an axis (used to delineate the thickness variation of a sill intruded in a sedimentary basin) and (2) inequality, proximity of a parameter estimate to a specified value, and proximity between pairs of parameter estimates (used to map a discontinuous basement relief). Usually, the stabilizing constraints are too restrictive to hold at all points of a given geological environment. In this case, we use different constraints in different sub‐areas. Each constraint is based on its compatibility with the actual geology of the subarea.

We present a constraint which incorporates in the geophysical inverse problem a priori information about the source convexity. Two kinds of convexity are considered: directional convexity, defined as the attribute of a body being intersected at most at two points by a straight line with a fixed spatial orientation, and global convexity, defined as the attribute of a body being intersected at most at two points by any arbitrarily oriented straight line. The interpretation model consists of several juxtaposed 2-D prisms whose positions, widths, and physical property are established by the interpreter. The depth to the top and the thickness of each prism are the parameters to be determined. The directional convexity is incorporated by minimizing a functional which expresses the number of times that several lines, parallel to a given direction and spaced at regular intervals, intersect the interpretation model. On the other hand, the global convexity is introduced in an algorithmic way by constraining the depths to the top and to the bottom of each prism to be, respectively, smaller or greater than the average of the depths to the top or to the bottom of the adjacent prisms. The use of the presented convexity contraint is demonstrated using synthetic gravity data produced by an isolated source. We found that convexity constraints (either directional along the horizontal axis or global) alone are insufficient to produce strongly stable solutions, although they reduce substantially the solution instability. Despite the “slight” instability, the directional convexity constraint leads to a reasonable delineation of the source shape and produces a better resolution of its base as compared with the inverse method imposing smoothness constraints on the top and bottom surfaces of the anomalous source. The potential use of the convexity constraint is the possibility of combining it with other constraints to produce geologically meaningful solutions. As an example, the global convexity was combined with the minimum moment of inertia of the anomalous masses about a vertical axis passing through its center of mass. This combination allowed introducing a priori information about a diapir shape. This is important because a maximum spread at the diapir middle portion indicates that it ascended like a viscous bubble in an isotropic medium. On the other hand, a spread located at its top indicates that the host medium is mechanically anisotropic, exhibiting horizontal planes of weakness where the mobile material was forced into. So, this combination of constraints alows introducing a priori information about the diapir emplacement. The combination of global convexity and minimum moment of inertia may also be applied to isolated intrusions, as illustrated with the real gravity anomaly produced by the granitic body of Castelsarrasin, in the Aquitaine Basin, France. The estimated source has its top located between 0.45 and 0.57 km from the surface, a base 6.4 km deep, and center of mass at a depth of 3.5 km, which agree reasonably with a priori geological and geophysical information about the body.

We present a new stable gravity inversion method applied to the mapping of an interface separating two homogeneous media. In contrast with previous similar methods, it does not impose an overall smoothness on the estimated interface to stabilize the solution. The density contrast between the media is assumed to be known. The interpretation model for the upper medium consists of rectangular juxtaposed prisms whose thicknesses represent the depths to the interface and are the parameters to be estimated. The true interface is assumed to be flat everywhere except at faults. To incorporate this attribute into the estimated relief, we developed an iterative process in which three kinds of constraints are imposed on parameters: (1) proximity between values of adjacent parameters, (2) lower and upper bounds to parameters, and (3) proximity between the values of parameters and fixed numerical values. Starting with an initial solution which presents an overall smooth relief, the method enhances initially estimated geometric features of the interface; that is, flat areas will tend to become flatter and steep areas will tend to become steeper. This is accomplished by weighting the constraints, which requires proximity between adjacent parameters. The weights are initialized with values equal to unity and are updated automatically to enhance any discrepancy between adjacent depths that have been detected at the initial solution. Constraints 2 and 3 are used both to compensate for the decrease in solution stability caused by the introduction of small weights and to reinforce flatness at the basin bottom. Constraint 2 imposes that any depth be nonnegative and smaller than an a priori known maximum depth value, whereas constraint 3 imposes that all depths be closest to a value greater than the maximum depth. The trade‐off between these conflicting constraints is attained with a final relief presenting flat bottom and steep borders. The method was tested with a synthetic gravity anomaly produced by a simulated sedimentary cratonic extensional basin whose basement consists of steep edges and a flat bottom. The results showed an improvement in the resolution of the relief, leading to a reliable mapping both of the sharp discontinuities at the borders and of the lateral extent of the base of the basin. Additionally, the method produced excellent estimates for the average dip angles of the basin edges (presumably controlled by normal faults), indicating, in this way, its potential in interpreting data produced by this kind of basin. The method was applied to the Bouguer anomaly from the northern portion of Steptoe Valley, Nevada, delineating an isolated basin with a wider, flat base and relatively straight borders as compared with the estimate imposing overall smoothness on the relief.

We present a new, stable method for interpreting the basement relief of a sedimentary basin which delineates sharp discontinuities in the basement relief and incorporates any law known a priori for the spatial variation of the density contrast. The subsurface region containing the basin is discretized into a grid of juxtaposed elementary prisms whose density contrasts are the parameters to be estimated. Any vertical line must intersect the basement relief only once, and the mass deficiency must be concentrated near the earth’s surface, subject to the observed gravity anomaly being fitted within the experimental errors. In addition, upper and lower bounds on the density contrast of each prism are introduced a priori (one of the bounds being zero), and the method assigns to each elementary prism a density contrast which is close to either bound. The basement relief is therefore delineated by the contact between the prisms with null and nonnull estimated density contrasts, the latter occupying the upper part of the discretized region. The method is stabilized by introducing constraints favoring solutions having the attributes (shared by most sedimentary basins) of being an isolated compact source with lateral borders dipping either vertically or toward the basin center and having horizontal dimensions much greater than its largest vertical dimension. Arbitrary laws of spatial variations of the density contrast, if known a priori, may be incorporated into the problem by assigning suitable values to the nonnull bound of each prism. The proposed method differs from previous stable methods by using no smoothness constraint on the interface to be estimated. As a result, it may be applied not only to intracratonic sag basins where the basement relief is essentially smooth but also to rift basins whose basements present discontinuities caused by faults. The method’s utility in mapping such basements was demonstrated in tests using synthetic data produced by simulated rift basins. The method mapped with good precision a sequence of step faults which are close to each other and present small vertical slips, a feature particularly difficult to detect from gravity data only. The method was also able to map isolated discontinuities with large vertical throw. The method was applied to the gravity data from Reco⁁ncavo basin, Brazil. The results showed close agreement with known geological structures of the basin. It also demonstrated the method’s ability to map a sequence of alternating terraces and structural lows that could not be detected just by inspecting the gravity anomaly. To demostrate the method’s flexibility in incorporating any a priori knowledge about the density contrast variation, it was applied to the Bouguer anomaly over the San Jacinto Graben, California. Two different exponential laws for the decrease of density contrast with depth were used, leading to estimated maximum depths between 2.2 and 2.4 km.

Euler deconvolution has been widely used in automatic aeromagnetic interpretations because it requires no prior knowledge of the source magnetization direction and assumes no particular interpretation model, provided the structural index defining the anomaly falloff rate related to the nature of the magnetic source, is determined in advance. Estimating the correct structural index and electing optimum criteria for selecting candidate solutions are two fundamental requirements for a successful application of this method. We present a new criterion for determining the structural index. This criterion is based on the correlation between the total‐field anomaly and the estimates of an unknown base level. These estimates are obtained for each position of a moving data window along the observed profile and for several tentative values for the structural index. The tentative value for the structural index producing the smallest correlation is the best estimate of the correct structural index. We also propose a new criterion to select the best solutions from a set of previously computed candidate solutions, each one associated with a particular position of the moving data window. A current criterion is to select only those candidates producing a standard deviation for the vertical position of the source smaller than a threshold value. We propose that in addition to this criterion, only those candidates producing the best fit to the known quantities (combinations of anomaly and its gradients) be selected. The proposed modifications to Euler deconvolution can be implemented easily in an automated algorithm for locating the source position. The above results are grounded on a theoretical uniqueness and stability analysis, also presented in this paper, for the joint estimation of the source position, the base level, and the structural index in Euler deconvolution. This analysis also reveals that the vertical position and the structural index of the source cannot be estimated simultaneously because they are linearly dependent; the horizontal position and the structural index, on the other hand, are linearly independent. For a known structural index, estimates of both horizontal and vertical positions are unique and stable regardless of the value of the structural index. If this value is not too small, estimates of the base level for the total field are stable as well. The proposed modifications to Euler deconvolution were tested both on synthetic and real magnetic data. In the case of synthetic data, the proposed criterion always detected the correct structural index and good estimates of the source position were obtained, suggesting the present theoretical analysis may lead to a substantial enhancement in practical applications of Euler deconvolution. In the case of practical data (vertical component anomaly over an iron deposit in the Kursk district, Russia), the estimated structural index (corresponding to a vertical prism) was in accordance with the known geology of the deposit, and the estimates of the depth and horizontal position of the source compared favorably with results reported in the literature.

We present a gravity interpretation method for estimating the relief of an arbitrary interface separating two homogeneous media. The upper medium is discretized into rectangular, juxtaposed prisms whose thicknesses represent the depths to the interface and are the parameters to be estimated from the gravity anomaly. The density contrast of each prism is assumed to be constant and known. To stabilize the inversion, we introduce two kinds of constraints on the depths. The first one requires proximity between the observed and computed depths at isolated points such as those obtained from boreholes (absolute equality constraint). The second one requires that groups of depths approximately follow an established linear relationship among the depths (relative equality constraint). Both kinds of constraints are imposed in the least‐squares sense. We illustrate the method performance by applying it to a synthetic anomaly produced by a simulated basement relief consisting of four narrow and adjacent structural lows. Only two structural lows produced isolated gravity lows. Nonetheless, the whole basement topography was successfully reconstructed with an average error of 4% of the maximum relief amplitude. In this example, the relative constraints established that the thicknesses of adjacent prisms should be as close as possible to each other (overall relief smoothness). As absolute constraints we used point information about the basement depth at five points. In addition to the overall relief smoothness, other relevant geologic information, such as localized relief smoothness, occurring at structural terraces, for example, can be incorporated by assigning different weights to the relative equality constraints. The method was applied to the gravity anomaly of Recôncavo Basin, Brazil, leading to a sharper definition of the basement features relative to previous gravity interpretations of this area.

We present an interpretation method for the gravity anomaly of an arbitrary interface separating two homogeneous media. It consists essentially of a downward continuation of the observed anomaly and the division of the continued anomaly by a scale factor involving the density contrast between the media. The knowledge of the interface depth at isolated points is used to estimate the depth [Formula: see text] of the shallowest point of the interface, the density contrast Δρ between the two media, and the coefficients [Formula: see text] and [Formula: see text] of a first‐order polynomial representing a linear trend to be removed from data. The solutions are stabilized by introducing a damping parameter in the computation of the downward‐continued anomaly by the equivalent layer method. Different from other interface mapping methods using gravity data, the proposed method: (1) takes into account the presence of an undesirable linear trend in data; (2) requires just intervals for both Δρ (rather than the knowledge of its true value) and coefficients [Formula: see text] and [Formula: see text]; and (3) does not require the knowledge of the average interface depth [Formula: see text]. As a result of (3), the proposed method does not call for extensive knowledge of the interface depth to obtain a statistically significant estimate of [Formula: see text]; rather, it is able to use the knowledge of the interface depth at just a few isolated points to estimate [Formula: see text], Δρ, [Formula: see text], and [Formula: see text]. Tests using synthetic data confirm that the method produces good and stable estimates as far as the established premises (smooth interface separating two homogeneous media and, at most, the presence of an unremoved linear trend in data) are not violated. If the density contrast is not uniform, the method may still be applied using Litinsky’s concept of effective density. The method was applied to gravity data from Recôncavo Basin, Brazil, producing good correlations of estimated lows and terraces in the basement with corresponding known geological features.

Current inversion methods estimating a physical property on a grid use two kinds of constraints: (1) absolute constraints requiring that the solution vector be close to a prespecified vector as, for example, proximity to the null vector in the ridge regression, and (2) relative constraints requiring that the elements of the solution vector be close to each other as, for example, smoothness imposed by methods requiring that spatial derivatives of the physical property be continuous. Using just absolute constraints may produce a biased solution conflicting with some true geological attributes of the source. On the other hand, using just relative constraints may be insufficient to stabilize the inversion. We present a stable inversion method incorporating absolute constraints only at points where the physical property is known as, for example, outcroppings and boreholes. Elsewhere, relative constraints are introduced according to a linear relationship for the spatial distribution of the physical property. The proposed method is analogous to an interpolation method applied to the physical property distribution: the interpolating function must satisfy some property, like continuity, and pass through the data. The proposed method is versatile in incorporating not only factual, but also indirect information such as continuity, symmetry, and trends. In addition, it may be applied both to linear and nonlinear problems. The relevant results obtained from applying the method to synthetic potential‐field data were: (1) ability to produce reliable mappings of the basement relief of a sedimentary basin (from gravity data) and of an areal distribution of magnetization, using only smoothness constraints and restricted surface and borehole information; (2) operational simplicity because of the solution insensitivity to the damping parameter. The method is also applied to the Bouguer anomaly over San Jacinto Graben, California, producing results consistent with previous interpretations.

Magnetic interpretations are usually carried out either by assuming induced magnetization and estimating the model geometry, or by presuming a known source spatial orientation to estimate the total magnetization. We present a 3-D magnetic interpretation method that estimates simultaneously the total magnetization direction and the spatial orientation of the source. It is based on the approximation of the anomaly by the series derived from expanding the magnetic potential into multipoles and retaining source moments up to second order. The moments and linear combinations of moments appearing in the series are then inverted from the magnetic anomaly. The total magnetization is assumed constant in direction but not in magnitude. It is also presumed implicitly that the anomalous distribution of magnetization intensity has nonzero values in a finite‐volume region, is far from the observation points, and presents three othogonal planes of symmetry intersecting at the center of the dipole moment. The method is essentially linear and requires no a priori explicit assumption of a fixed geometry for the sources. The method is particularly suited to interpret compact, isolated or disjoint, but spatially correlated sources. If the source satisfies all assumptions presumed by the method, it is possible to obtain accurate, stable estimates of the total dipole moment vector, the position of the center of dipole moment, and the directions of all three principal axes of symmetry. If the source is not far from the observation plane and/or if the total magnetization direction is not constant, it is still possible to obtain accurate and stable estimates of the direction of the mean total magnetization and the projection, on the observation plane, of the center of dipole moment. The method is applied to magnetic data from the Gulf of Guinea Seamount. The estimated magnetic palaeopole is at 50°48′S and 74°54′E which is in good agreement with estimates published by other authors.

We present a 3-D gravity interpretation method based on the inversion of source moments employing the series derived from the gravity anomaly expansion in multipoles and retaining moments up to second order only. It presents the advantages of being essentially linear and allowing a straightforward computer implementation. The method requires neither an explicit assumption about the source geometry nor a priori knowledge about the density contrast distribution, which may even be nonuniform. The method assumes implicitly: (1) that the source fits in the interior of an imaginary sphere whose center coincides with the source center of mass and whose radius is smaller than the depth to the source center of mass relative to the measuring plane, and (2) that the spatial density distribution presents three ortogonal planes of symmetry intersecting at the source center of mass. The first assumption can be met by upward continuing the observed anomaly. When both assumptions are met, the method produces reliable and stable estimates of the total anomalous mass, the coordinates of the center of mass, the three principal axes directions of the anomalous body, and the relative importance among the axes. The method is particularly suited for interpreting compact, isolated or disjoint, but spatially correlated sources. The method is sensitive to an incomplete removal of the regional field (presumably superimposed on the residual anomaly of interest) because the unremoved part of the regional field may strongly affect the estimates of the residual source moments.

Interpolation using only the observations at discrete points is an ill‐posed problem because it admits infinite solutions. Usually, to reduce ambiguity, a priori information about the sample function is introduced. Current interpolation methods in mineral exploration introduce only the constraints of continuity and smoothness of the interpolating function. In interpolating potential‐field anomalies, the constraint that the sampled function is harmonic may be introduced by the equivalent‐layer method (ELM). We compare the performance of the ELM and the minimum curvature method (MCM) in interpolating potential‐field anomalies by applying these methods to synthetic magnetic data simulating an aeromagnetic survey. In the case the anomaly flanks and peak are undersampled, the ELM performs better than the MCM in recovering the anomaly gradients and peak. In the case of elongated linear anomalies, the ELM recovers the exact linear pattern, but the MCM introduces spurious oscillations in the linear pattern. Also, the ELM is able to reduce the data from a survey flown at different heights to a common level. In contrast, the MCM, being a 2-D interpolation method, cannot account for variations in the vertical coordinates of the observation points.

We present an empirical ambiguity analysis method based on a finite number of acceptable solutions that are representative of the ambiguity region. These solutions are submitted to a Q‐mode factor analysis that indicates which parameters are ambiguous and their ambiguity range. We illustrate, with a synthetic nonlinear example, that our method is more effective than singular value decomposition analysis in producing an average trend of the ambiguity region. It requires less restrictive hypotheses and is more robust than analytical methods of ambiguity analysis, in the sense of being applicable to a broader class of problems.

The equivalent layer calculation becomes more efficient by first converting the observed potential data set to a much smaller equivalent data set, thus saving considerable CPU time. This makes the equivalent‐source method of data interpolation very competitive with other traditional gridding techniques that ignore the fact that potential anomalies are harmonic functions. The equivalent data set is obtained by using a least‐squares iterative algorithm at each iteration that solves an underdetermined system fitting all observations selected from previous iterations and the observation with the greatest residual in the preceding iteration. The residuals are obtained by computing a set of “predicted observations” using the estimated parameters at the current iteration and subtracting them from the observations. The use of Cholesky’s decomposition to implement the algorithm leads to an efficient solution update everytime a new datum is processed. In addition, when applied to interpolation problems using equivalent layers, the method is optimized by approximating dot products by the discrete form of an analytic integration that can be evaluated with much less computational effort. Finally, the technique is applied to gravity data in a 2 × 2 degrees area containing 3137 observations, from Equant‐2 marine gravity survey offshore northern Brazil. Only 294 equivalent data are selected and used to interpolate the anomalies, creating a regular grid by using the equivalent‐layer technique. For comparison, the interpolation using the minimum‐curvature method was also obtained, producing equivalent results. The number of equivalent observations is usually one order of magnitude smaller than the total number of observations. As a result, the saving in computer time and memory is at least two orders of magnitude as compared to interpolation by equivalent layer using all observations.

Ambiguity in a geophysical interpretation problem is the possibility of accepting more than one solution caused either by solution nonuniqueness or instability. Nonuniqueness is related to the existence of more than one solution regardless of the precision of observations. On the other hand, instability is related to the acceptance of different solutions producing data fittings within the expected observational errors. We studied the ambiguity in the inversion of well‐logging data using a method based on the analysis of a finite number of acceptable solutions, which are ordered, in the solution space, according to their contributions to the overall ambiguity. The analysis of the parameter variations along these ordered solutions provides an objective way to characterize the most ambiguous parameters. Because this analysis is based on the geometry of an ambiguity region, empirically estimated by a finite number of alternative solutions, it is possible to analyze the ambiguity due not only to errors in the observations, but also to discrepancies between the interpretation model and the true geology. Moreover, the analysis can be applied even in the case of a nonlinear interpretation model. The analysis was performed with recorded data, and compared with the analysis using singular value decomposition, leading to comparable results. Following the determination of the most ambiguous parameters, a reparameterization is possible by grouping these parameters into a single parameter leading to a simpler interpretation model and, therefore, to a drastic reduction in the ambiguity.

Extending the compact gravity inversion technique by incorporating a priori information about the maximum compactness of the anomalous sources along several axes provides versatility. Thus, the method may also incorporate information about limits in the axes lengths or greater concentration of mass along one or more directions. The judicious combination of different constraints on the anomalous mass distribution allows the introduction of several kinds of a priori information about the (arbitrary) shape of the sources. This method is particularly applicable to constant, linear density sources such as mineralizations along faults and intruded sills, dikes, and laccoliths in a sedimentary basin. The correct source density must be known with a maximum uncertainty of 40 percent; otherwise, the inversion produces thicker bodies for densities smaller than the true value and vice‐versa. Because of the limitations of the inverse gravity problem, the proposed technique requires an empirical technique to analyze the sensitivity of solutions to uncertainties in the a priori information. The proposed technique is based on a finite number of acceptable solutions, presumably representative of the ambiguity region. By using standard statistical techniques, each parameter is assigned a coefficient measuring its uncertainty. The known hematite and magnetite ore body shape, in the vicinity of Iron Mountain, MO, was reproduced quite well using this inversion technique.

Apparent density maps are derived from observed residual Bouguer maps under the assumptions that the sources are restricted to a depth interval and that the density distribution is not a function of the vertical coordinate. If the depth interval containing the sources is known, the computed apparent density maps are very close to the true density distributions. Even when the depth interval of the sources is not known, but the sources have constant density, their horizontal extent may be delineated by the half‐maximum contours of the apparent densities. We developed an interpretation method using a family of curves based on source density, depth to top of sources, and source thickness. The main advantage in the use of these families of curves is the broad range of possible interpretations involving the above‐mentioned parameters. Rough semiquantitative interpretations are possible in the absence of any further a priori information; more refined semiquantitative interpretations require additional quantitative knowledge of one parameter and just a lower or upper bound for another parameter. Finally, quantitative interpretations are possible only with additional quantitative knowledge of two of the three parameters involved.

We combine a stabilized reduction‐to‐the‐pole and an upward continuation filter to produce meaningful reduced‐to‐the‐pole fields at low magnetic latitudes. The stabilizing procedure is based on the development, in Taylor’s series, of the theoretical expression for the reduction‐to‐the‐pole filter in the wavenumber domain. The filter instability is caused by the huge filter amplitudes along the magnetization azimuth, which are expressed by the infinite sum of terms close to unity. The stabilizing procedure reduces to simply truncating the infinite series. The upward continuation filter attenuates the high wavenumber component of the noise and allows us to design a stabilized filter closer to the theoretical one. Besides, quantitative interpretations of source depth based on the filtered field are more reliable when using upward continuation as compared with arbitrary low‐pass filters. The proposed filter was applied to synthetic data of a single prism uniformly magnetized along a supposedly known direction, and it produced a reduced‐to‐the‐pole field very close to the theoretical field at pole. We also applied the filter to magnetic data from Dixon Seamount assuming induced magnetization only. We obtained, within the central part of the anomaly, roughly circular contours of the reduced‐to‐the‐pole anomaly due to the nearly circular shape of the Seamount (evidenced by topographic data).

We present an inversion approach that estimates the basement relief of a fault-bounded sedimentary basin. We parameterize the Earth’s subsurface, containing the sedimentary pack, into prismatic cells, with known horizontal dimensions and known density contrast, and estimate their thicknesses, which represent the depths to the basement. To obtain depth-to-basement estimates we use the total-variation regularization as a stabilizing function that allows the estimation of a faulting basement relief. Tests on synthetic and on field data from Almada Basin, Brazil, and from Steptoe Valley, Nevada, confirmed the potential of our method in detecting and locating normal faults in the basement relief of a sedimentary basin.

Standard polynomial fitting methods are inconsistent in their formulation. The regional field is approximated by a polynomial fitted to the observed field. As a result, in addition to the nonuniqueness in the definition of the regional field, the fitted polynomial is strongly influenced by the residual field (observed field minus regional field). We present a regional‐residual separation method for gravity data which uses a robust procedure to determine the coefficients of a polynomial fitted to the observations. Under the hypothesis that the regional can be modeled correctly by the polynomial surface, the proposed method minimizes the influence of the residual field in the fitted surface. The proposed method was applied to real gravity data from Ceará state, Brazil, and produced information on zones of possible crustal thickening and the occurrence of lower‐crustal granulitic rocks thrust into the shallow subsurface.

We present a 3D gravity inversion for estimating the location and geometry of 3D salt bodies. Iteratively, our approach estimates a 3D density-contrast distribution by assuming a piecewise constant function defined on a userspecified grid of cells. To estimate a unique and stable density-contrast distribution, we look for the solution that fits the observed anomaly within the measurement errors and favors compact salt bodies closest to pre-specified geometric elements, such as axes and points. Our inversion method consists of two nested iterative loops. The outer loop uses an adaptive learning strategy that starts with a coarse grid of cells, a set of first-guess geometric elements (axes and points) and the corresponding assigned density contrast. From the second iteration on, this strategy refines the grid and creates a new set of geometric elements (points only) and associated density contrasts. The inner loop estimates the density-contrast distribution for the grid of cells and for the set of geometric elements defined in the outer loop. Each geometric element operates as the firstguess skeletal outline of a particular, homogeneous section of the salt body to be imaged. The algorithm easily incorporates known density contrast information by assigning these values to each geometric element. The iteration stops when the geometries of the estimated salt bodies are invariant along successive iterations. We apply our method to synthetic gravity data obtained from a salt body, whose density contrast with the host rocks varies from positive to negative. We tested two different geologic hypotheses about the geometric elements and the assigned density contrast on the real gravity anomaly from Galveston Island salt dome, offshore Texas, USA. In the first one, we estimated a vertically elongated salt dome with depth to the bottom at 5 km. In the second hypothesis, we estimated a cylindrical salt dome with an associated overhang and depth to the bottom at 4 km.