An alternative 3 D inversion method for magnetic anomalies with depth resolution

This paper presents a new method to invert magnetic anomaly data in a variety of non-complex contexts when a priori information about the sources is not available. The region containing magnetic sources is discretized into a set of homogeneously magnetized rectangular prisms, polarized along a common direction. The magnetization distribution is calculated by solving an underdetermined linear system, and is accomplished through the simultaneous minimization of the norm of the solution and the misfit between the observed and the calculated field. Our algorithm makes use of a dipolar approximation to compute the magnetic field of the rectangular blocks. We show how this approximation, in conjunction with other correction factors, presents numerous advantages in terms of computing speed and depth resolution, and does not affect significantly the success of the inversion. The algorithm is tested on both synthetic and real magnetic datasets. Mailing address: Dr. Alessandro Pignatelli, Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: pignatelli@ingv.it


Introduction
It is well known how inversion applied to surface potential field data is a mathematical ill-posed problem which suffers from ambiguity of solutions (Blakely, 1995).This is due to the inherent nonuniqueness of the problem and the inaccuracy and finite number of the potential field observations.
A common strategy to approach the inverse problem is to divide the subsoil into homogeneous magnetized rectangular blocks.The direction of magnetization is supposed to be known in advance and equal for each block.The unknown parameters to be solved are the module of magnetizations.These parameters can be solved by means of a linear system with a number of equations equal to the data available.Unfortunately the magnetic stations are generally lower than the blocks, so the linear system is undetermined and the solution can be obtained only by simultaneously minimizing an objective function and the misfit between the observed magnetic field and the field produced by the magnetized blocks (Hoerl and Kennard, 1970;Last and Kubik, 1983;Silva and Hohmann, 1983;Guillen and Menichetti, 1984;Li and Oldenburg, 1996).The main disadvantage of this strategy is that the definition of the objective function depends on a priori knowledge of the source and this makes algorithms non general, so that a unique approach for inverting magnetic data cannot be used in every situation.However the use of the objective function is necessary in order to reduce the instability of the solutions.
The algorithm proposed is based on the hypothesis that no information about the source is available.Thus a general and appropriate objec-tive function has to be chosen.In this work we show that the use of the norm minimization of solution as objective function is a good choice.
However the inversion procedure using the minimization of norm is influenced by the rapid decrease of the signal amplitude as a function of depth, so the solutions tend to be located only in the upper layers of the earth space model (Fedi and Rapolla, 1999).So an inversion based on a simultaneous minimization of the misfit of the data and the Euclidean norm of the solution causes the lack of depth resolution.To overcome this problem it is possible to use a weighting function (Li and Oldenburg, 1996;Portnianguine and Zhdanov, 2002).Alternatively, Fedi and Rapolla (1999) have proposed to use a multi-observation level dataset, arguing that this kind of data provide depth resolution.However, this statement is still a subject of some debate (Oldenburg and Li, 2003).An additional limitation of inversion methods is related to the long computing time involved in the calculation of the magnetic field of a set of rectangular prisms, which can be implemented using, for instance, the formulation deduced by Sharma (1986).
This paper describes a new magnetic inversion method in which the magnetic field of a rectangular prisms is approximated by the field of a magnetic dipole.This approximation does not affect significantly the success of the inversion method.Furthermore, this approach makes it possible, on the one hand, the choice of a weighting function valid for every magnetic source geometry and, on the other, to improve the calculations in terms of computing time.After a description of our inversion technique, we show its efficiency on both synthetic and real magnetic data sets.

Description of the magnetic inversion method
Let us assume we have an anomaly data set y of N observations y=(y 1 , y 2 , ..., y N ).
(2.1) From the electromagnetic theory, we can obtain the value of the magnetic field along a generic direction f produced by a magnetization distri-bution J(r) as follows: (2.2) where rl is the position vector of the observation point, t is a unit vector in the direction of magnetization and V is the volume containing all magnetic sources.Equation (2.2) can be simplified if a uniform magnetization direction, t, for all the distribution is considered.In fact, in most real cases, a constant direction parallel to the magnetic ambient field is chosen.This approximation is reasonable when the most important component of the total magnetization vector is the induced one or when the remanent magnetization is parallel to the earth's magnetic field.Despite this simplification, calculation of this integral can only be achieved analytically if the geometry of volume V, and the distribution of magnetization J(r) are simple enough.
Theoretically, the purpose of magnetic inversion is to calculate a magnetization distribution from the observed magnetic anomaly data.However, this is not possible because an analytic solution of the integral eq.(2.2) in the parameter J(r) does not exist.For that reason, the common approach consists of dividing the region containing magnetic sources into elementary rectangular prisms cells in which magnetization is assumed to be uniform.Thus, eq. ( 2.2) can be rewritten in the following way: (2.3) where J i is the module of the magnetization of the ith prism, M is the number of blocks, rj is the position vector of the jth observation point and xi,j is the discrete version of the Green's function associated with eq.(2.2), which depends not only on the block geometry, but also on the relative position between the ith cell and the jth observation point.
The discretization process produces a linear system in which the unknown parameters are the magnetizations Ji of each cell.Generally, systems of this kind are underdetermined and overconstrained (Jackson, 1972), so that opti- # mization techniques become necessary in order to evaluate the blocks magnetization.We use the Levemberg and Marquardt technique (Press et al., 1992) minimizing simultaneously the norm of the solution and the misfit between the observed and the computed field.By means of this iterative technique it is possible to introduce the positivity constraint by setting to zero the negative magnetization values.The iterative process is stopped when there are no significant improvements in data convergence.The effect of the minimization of the norm is the concentration of the solution in the upper layers of the cells instead of around the real depth due to the decay of the kernel amplitude with depth.Our method is designed to minimize this problem.
Let us multiply and divide each term of the sum by a weighting function where z ij is the distance between the ith dipole and the altitude of jth station.
Replacing Ji/ w(zij) by J r i and w(zij)xij by Xij, the linear system to be solved becomes (2.5)By solving this new linear system with the iterative method we obtain the set J r i.The initial guess model is the one with all the magnetizations null.
Through this procedure, a geometrical distribution of weighted magnetizations for every block is achieved.However, it is not convenient to recover Ji by simply multiplying J r i by w(z).In fact, due to simultaneous minimization of the misfit between the observed data and the computed field together with the euclidean norm of solution assuming only positive values for magnetization modules, we obtain a partial fit of the data during the inversion procedure.
The main task of this work has been to find an appropriate weighting function and a method for correcting magnetization values, as explained in detail in the following paragraphs.

The weighting function
Equation (2.3) can be used to calculate the magnetic field produced by a given set of magnetized blocks.Usually, the geometrical factors xi,j are computed using the analytical formula outlined by Sharma (1986), which calculates the magnetic anomaly produced by a uniformly magnetized prismatic block.Unfortunately, this approach has the great disadvantage of being CPU-time consuming.
To reduce computational time, we have replaced Sharma's analytical expression with the field produced by magnetic dipole located at the geometrical center of each elementary cell.The dipolar approximation is valid if the distance between the observation point and the source is much higher than the source size, a requirement that is ensured if a relatively small cell size is chosen.
The real advantage of the dipole approximation, however, is that we can define the weighting function as the dipolar field decay (2.6) where z ij is the vertical distance between the ith dipole and the altitude of jth station.

Correction of the magnetization values
The values of J r i obtained by inverting eq.(2.5) by means of Levemberg and Marquardt technique (Press et al., 1992) are the obtained final parameters which best fit simultaneously the anomaly field and the euclidean norm of the solution.
In order to minimize only the misfit of the data, let's write the final magnetizations J i f in the following way: (2.7) where Ai are corrective factors to be determined.We assume that these factors depend only on the depth to the dipole and not on its horizontal position, that is , , ..., , , ..., where zi is the depth of the jth dipole measured from the horizontal plane z=0.
In order to compute the magnetization of each block, we just need the set of parameters: A 1(z1), A2(z2), ..., Ah(zh), where h is the number of layers along the z-axis chosen to perform the inversion.To compute the right factors, it is necessary a new data misfit minimization in which the following equation must be satisfied: (2.9) Now, the number of unknown parameters, Ai(zi), is lower than the number of data and the problem of undertermination has been overcome.So we do not need the minimization of the norm together with the minimization of the misfit of the data.After this second inversion, we can obtain the magnetization value of the ith block by applying . (2.10)

A synthetic case
The proposed inversion method has been tested using a synthetic and a real data sets.For the synthetic test we have constructed two magnetic sources with a complex geometry located at different depths, with the aim of reproducing actual geological structures (fig.1a).The magnetization intensity is of 2.5 A/m for the shallow source and 5.0 A/m for the deep one.
The magnetic anomaly produced by these structures has been calculated on a 21×21 km 2 regular grid with a cell size of 1 km, at an altitude of 100 m (fig. 1b).In addition, to simulate a measured anomaly field, a Gaussian noise with 0 mean and standard deviation of 1 nT was added to the calculated anomaly.The inversion was then performed on a 12×12×5 km lattice , , ..., , , ..., with a 300×300×250 m elementary block cell, which has been selected in order to not match exactly the synthetic sources.On the left hand side of fig.2, the result of the inversion previous to the correction of magnetization values is shown.The various slices represent horizontal cross sections of the lattice at different depths.
The position of the synthetic sources is outlined using thick black lines.The magnetization value correction was then applied according to eq.
(2.10), and the new calculated source distribution is shown on the right hand side of fig. 2. This graphical representation is particularly suitable to easily compare both results at the same depths.Without the correction factors, the solution appears more elongated on the vertical axis, whereas with our method a more reliable location of the sources, along the vertical axis, and a better characterization of the magnetic properties is achieved.

A real case
To test our method on a real case, we have chosen a data set acquired during a ground-based magnetic survey conducted for environmental investigations (Marchetti et al., 1998).In this case study, the parameters of the sources and their exact location are well known in advance, therefore it is an excellent test for our purposes.A set of 12 empty 55-gallons steel drums were buried at a depth of 4.5 m (fig.3a).The magnetic anomaly field produced by the drums was measured at 0.75 m above the ground and the interpolated grid with a cell size of 1 m is shown in fig.3b.
The inversion was then applied using a 22× ×25×10 m lattice with a 1×1×0.5 m elementary block cell.Without any correction, the solution of the inversion is shown in fig. 4 (left hand side).The obtained magnetic sources show an elongated pattern along the vertical axis which is similar to the one observed in the synthetic case.The application of our algorithm improves the source characterization, as shown on right hand side of fig. 4, where the solution is more compact.The obtained source magnetization is about 3 A/m, which is several orders of magnitude lower than the magnetization of steel (Ravat, 1996;Marchetti et al., 1998).This apparent mismatch is due to the peculiar shape of the source, with non-magnetic or empty space inside the drums.Therefore, the estimated magnetization value M Z would represent a bulk magnetization for all the source volume, as given by (3.1) where M(x, y, z) represent the distribution of magnetization in the volume V occupied by the steel drums.

Conclusions
We have proposed a new method for the 3D inversion of magnetic anomalies when a priori knowledge on the sources and a multi-observation level dataset are not available.We have shown that it is possible to estimate the location and magnetization intensity of magnetic sources by means of a discretization of the subsoil into elementary parallelepipeds whose magnetic effect is approximated by a dipolar field.The main advantage of this approach, besides a great improvement in computing speed, is that it allows the choice of an analytic weighting function through which the source is obtained with good depth resolution.Compared with other methods, this weighting function does not depend on the geometry of the elementary cells,   and therefore is valid for every discretized region.Finally, an adjustment on the estimation of the unknown parameters is proposed in order to further improve our knowledge about the location and magnetization intensity of the sources.We have showed the results on a synthetic magnetic field dataset.

Fig. 1a ,Fig. 2 .
Fig. 1a,b.a) Synthetic magnetic source; b) total-field anomaly produced at an altitude of 100 m above the surface.Grid spacing: 1 km.Fig.2. Results of the inversion using synthetic data.The column on the left-hand side shows the solution without the magnetization values correction, whereas the right-hand side reports the results calculated after the magnetization correction.The slices represent horizontal cross sections of the lattice at various depths.The inversion was performed on a 12×12×5 km lattice with 300×300×250 m elementary cells.The thick black lines outline the cross sections of the true source.

Fig
Fig. 3a,b.a) Steel drums buried at a depth of about 4.5 m; b) magnetic anomaly field produced at an altitude of 0.75 m above the surface.Grid spacing: 1 m.

Fig. 4 .
Fig. 4. Results of the inversion using real data.The column on the left-hand side shows the solution without the magnetization values correction, whereas the right-hand side reports the results calculated after the magnetization correction.The slices represent horizontal cross sections of the lattice at various depths.The inversion was performed on a 22×25×10 m lattice with 1×1×0.5 m elementary cells.