# Geophysics

- Society of Exploration Geophysicists

## Abstract

In archaeological magnetic prospecting, most targets can be modeled by a single layer of constant burial depth and thickness. With this assumption, recovery of the magnetization distribution of the buried layer from magnetic surface measurements is a 2D deconvolution problem. Because this problem is ill posed, it requires regularization techniques to be solved. In analogy with image reconstruction, the solution showing the resolved subsoil features can be considered a focused version of the blurred and noisy magnetic image. Exploiting image deconvolution tools, two iterative reconstruction methods are applied to minimize the least-squares functional: the standard projected Landweber method and a proposed modification of the iterative space reconstruction algorithm. Different regularization functionals inject a priori information in the optimization problem, and the split-gradient method modifies the algorithms. Numerical simulations in the case of perfect knowledge of the impulse response functions demonstrate that the edge-preserving, total-variation functionals give the best results. An iterative semiblind deconvolution method to estimate the burial depth of the source layer was used with a real data set to test the effectiveness of the method.

## INTRODUCTION

One of the most difficult problems in potential-field prospecting is to extract information about the burial depth of the sources causing the measured anomalies. However, in some cases, the survey focuses more on delineating the horizontal position, lateral extent, and edges of the sources rather than their burial depths. For example, lineaments extrapolated from magnetic maps are used for small-scale regional geologic interpretation. Another interesting case is the use of magnetics in archaeological exploration (Bozzo et al., 1994; Herwanger et al., 2000). In fact, magnetic prospecting is used routinely to map buried stone/brick foundation structures, roads, or graves and to outline the locations of kilns, hearths, and ferrous objects at archaeological sites (Wynn, 1986 and references therein).

The main problem in archaeological geophysics is to convert the geophysical results into an image that resembles a plan view of the buried antiquities, i.e., to have a final image that nonexperts can read easily and that looks like an aerial photograph or a drawing of the ruins as if an excavation had taken place. Archaeological interpretation of magnetic data requires detecting the signal from the background noise and reconstructing images of the subsoil features, showing location, shape, and dimensions of the archaeological targets in a form easily read by nongeophysicists.

To achieve these goals, additional data processing generally is needed, with emphasis on image enhancement and data filtering or inversion (Jeng et al., 2003; Tsivouraki and Tsokas, 2007; Tassis et al., 2008). Treating the buried bodies as a set of vertical-sided rectangular prisms is a technique frequently used in 2D and 3D modeling of buried archaeological structures. In a simplified approach (the single layer), the prisms are assumed buried at equal depths (Tsokas and Papazachos, 1992) and characterized by homogeneous induced magnetization or, in the case of remanent magnetization, along a known direction. This single-layer model is effective in situations such as foundation structures, defense walls, roads, and graves. In cases when remanent magnetization is relevant but unknown (as for buried kilns or ferrous bodies), usually target localization can be successfully reached without further image enhancement (Bozzo et al., 1994). If the effect of remanent magnetization is dominant, the model is no more valid. In fact, the linearization discussed below is performed by supposing that the anomaly vector is parallel to the inducing field.

Traditionally, the inverse problem has been solved in the frequency domain by the apparent-susceptibility mapping method or in the space domain by 2D deconvolution. Apparent-susceptibiity mapping (e.g., Urquhart and Strangway, 1985) consists of reduction to the pole of the total field, downward continuation, and recovery of the distribution of the apparent magnetization by dividing by the shape factor of the prism. The 2D deconvolution approach consists of the convolution of the total field with an optimum filter operator (Tsokas and Papazachos, 1992) or a stochastic approach (Aykroyd et al., 2001; Glazunov and Efimova, 2003). In both cases, some low-pass filtering must be applied to suppress high-frequency noise.

Major drawbacks of traditional methods arise because, in practice, the single-layer model can only roughly approximate the actual distribution of magnetization in the subsoil. The effects from a possible variable burial depth and magnetization of the archaeological targets, their geometric complexity, and the presence of other magnetization contrasts (e.g., geologic noise) also must be taken into account. Moreover, required a priori estimates of burial depth and depth extent of sources may be biased. All of these effects may produce anomalous and complex magnetization distributions in solutions based on the simplified assumptions. Nevertheless, the single-layer model is effective, simple, and robust because in archaeological prospecting we often look for ruins at very shallow depth and in rather small areas, so the simplified assumptions are reasonable.

These traditional methods may be improved if we insert more suitable a priori information about the target and the possibility of an automated estimate of burial depth and depth extent of the layer when they are unknown. Therefore, we develop an approach that, taking advantage of iterative image-restoration tools, includes such a priori information in the algorithm in the form of constraints. Adopting the semiblind deconvolution method, we can obtain an estimate of burial depth and depth extent of the layer.

## THE SINGLE-LAYER MODEL

The quantity usually measured in magnetic surveys is the modulus of total geomagnetic field *m* the intensity of magnetization. By projecting

The single-layer model (Tsokas and Papazachos, 1992) assumes that the sources belong to only one layer, with the top at *h* as the depth extent of the layer. The intensity of magnetization does not depend on

Hence, if the object to be retrieved belongs to only one layer constituted of homogeneous vertical-sided rectangular prisms buried at depth *h* and centered at positions

The term *ε* represents additive noise, defined as instrumental and cultural noise plus the effect of deviations from the simple model assumptions. In magnetic prospecting, the measured magnetic field is acquired by a magnetometer (typically proton precession, fluxgate, or Cesium vapor) that produces read-out noise, which is described by an additive Gaussian noise. “Geologic” noise cannot be modeled because it depends on the object to be reconstructed and is an interpretive task. The most pronounced component comes from microundulation of the relief of the ground surface where the readings are taken. Another source of geologic noise is the plowing track, which create severe coherent noise. This kind of noise is difficult to discriminate from the effect of linear subsurface structures; however, one can keep notes during the survey and then omit the features by directional filtering.

Assuming a single layer, the IRF is space invariant on every horizontal plane, changing only with the variation of burial depth and depth extent. In our case, knowledge of the IRF is equivalent to knowledge of the burial depth of the layer. The number of magnetic data points is often smaller than the desired number of prisms that compose the single layer of unknown susceptibility sources; but when using the fast Fourier transform (FFT), we must match the number of data points and unknowns. So we divide the single layer into a number of prisms equal to the amount of magnetic data and sample the IRF accordingly.

Major computational advantages of the convolution approach are smaller array size and ability to use the FFT algorithm (reducing the computational complexity from

## THE OPTIMIZATION PROBLEM

The inverse problem is ill conditioned; hence, the presence of noise completely corrupts the result of the inversion of detected data, usually called the generalized solution, in such a way that it has no physical meaning. By exploiting image deconvolution tools, we solve this problem using certain iterative algorithms derived from the maximum likelihood approach (Vogel, 2002). These algorithms allow us to reconstruct the solution iteratively, starting from the lower frequencies and adding the higher ones step by step. Knowing the variance of the noise, we can stop the iterations before adding to the solution the high frequencies corresponding to noise amplification. Moreover, these algorithms are very flexible tools, giving us the opportunity to customize them to add a priori constraints on the solutions. This can be achieved by suitable regularization functionals.

To perform the inversion, we focus on two different iterative methods adapted for IRF that, instead of classical point-spread functions, can have negative values.

### Least-squares methods

The maximum likelihood approach, in the case of additive Gaussian noise, is equivalent to the minimization of the least-squares functional:

Our deconvolution problem is a constrained minimization issue because we search for nonnegative solutions. The susceptibility of the buried object is always considered positive (Li and Oldenburg, 1996). To construct the iterative algorithms, one needs to compute the gradient of *k* with respect to the origin, i.e.,

#### Projected Landweber method

Generally, the iterative projected Landweber method (PLM) is used to minimize the functional *C*). Being a gradient method, it has the form
*C* and *τ* is the step length. This method alternates a step in the descent direction and a projection of *f* on *C* (Bertero and Boccacci, 1998). The projection on the positive cone is obtained at each iteration by setting all negative pixels to zero. The convergence of the algorithm is proved in Eicke (1992). However, because the least-squares solution is corrupted completely by the amplification of noise, the important feature of the algorithm is the so-called *semiconvergence* (Engl et al., 1996): the iterates first approach the correct solution, then go away. Therefore, regularization can be imposed by stopping the iterations early. The drawback is that, in real cases, it is not easy to estimate the number of iterations providing the optimal solution.

#### Modified iterative space reconstruction algorithm (ISRA)

In our problem, the solution is nonnegative, but the IRF and the data can take negative values. In such a situation, we cannot use the iterative space-reconstruction algorithm of Daube-Witherspoon and Muehllehner (1986), whose convergence is proved by De Pierro (1987). However, by using the split-gradient method (SGM) (Lanteri et al., 2001), one can introduce a modified ISRA that can be used.

For this purpose, we write the Karush-Kuhn-Tucker conditions for the constrained minima of the convex functional (equation 7):

Next, we apply the successive approximation method to the fixed-point equation 12 to obtain the modified ISRA algorithm:

### Regularization methods

A comparison with the image deconvolution problem offers several techniques to retrieve the optimal solution. Thanks to regularization theory, we can take into account some a priori information about the solution. Adopting a Bayesian approach, we consider the constrained minimization of a functional of the form
*α* is a positive multiplier called regularization parameter (Bertero and Boccacci, 1998).

This procedure is called regularization of the inverse problem, and

In many archaeological surveys, we can assume that we are seeking a plan view, so we can use an edge-preserving regularizing functional. Generally it has the form (Vogel, 1997)
*β* is introduced to avoid the singularity of the gradient at the origin. Alternative choices of *ψ* lead to other edge-preserving functionals (Vogel, 1997). For example, the Geman-McClure prior (Geman and McClure, 1987),
*ψ*, the gradient of *ψ*. These two edge-preserving methods allow the inversion procedure to preserve the blocky structures present in the images.

#### Regularized projected Landweber algorithm

The regularized PLM takes the form
*τ*. We apply a dichotomic search of the minimum of
*τ*.

#### Regularized modified ISRA

Regularization of the modified ISRA takes the form
*f* and the function

#### Estimate of *α* and stopping criterion

The choice of the regularization parameter *α* is crucial to obtain a solution that is physically admissible. In the Tikhonov regularization, there is an optimum *α* for which the solution has the minimum distance from the true image *α* is a weight that balances between the minimization of the discrepancy and that of the *α* so the relative reconstruction error, defined by equation 30, reaches the minimum possible value during the iterative process.

For the real data set, we cannot use the relative reconstruction error, so we adopt the discrepancy principle for stopping the iterations and an interactive process to choose *α*. Defining
*σ* becomes smaller than the estimated standard deviation of the noise in the acquired data.

### Semiblind deconvolution

In deconvolution methods, one assumes that the IRF is known exactly. However, in practical applications, the burial depth and the depth extent of the layer are unknown. An estimate of the depth of archaeological targets can be reached using prior knowledge of the site or other geophysical techniques. For example, ground-penetrating radar or geoelectrical prospecting can provide insight about the depth of the buried bodies. However, it is often hard to compute this parameter.

If the magnetization is approximatively constant within the layer, then the IRF depends mainly on *h*, we have a one-parameter family of IRFs, and the deconvolution problem becomes semiblind. Hence, we consider *k* depending on the parameter

Initially, we compute a discrete set of admissible IRFs, corresponding to a given set of burial depth values:
*j*. This pair enters in the object box. The deconvolution method implemented in this box is the iterative PLM with TV regularization. By keeping

Then the pair

In the second step of the iteration (lower rectangle, Figure 2), if

To initialize the procedure, we take as

The algorithm guarantees that the IRF is physically meaningful; moreover, it automatically estimates the burial depth of the layer.

## NUMERICAL EXPERIMENTS

We have tested the inversion scheme using synthetic models. We performed our numerical experiments using both Landweber and modified ISRA, focusing on three objectives:

1) To test the inversion scheme, exploiting the performances of the two iterative methods

2) To compare the performances of different kinds of regularization

3) To test the robustness of the method (using a semiblind deconvolution approach) when the estimate of the IRF (and hence of the burial depth and depth extent of the targets) is biased

The relative reconstruction error, computable only in the case of simulations, is defined by
*j* is the iteration number. It is a measure of the distance between the solution at the

### Inversion tests with exact IRF

We simulated a buried archaeological target that may represent a set of partially collapsed walls. The object depth extension is

We divide the subsoil into a set of homogeneous vertical prisms, considering only a single layer of cells. The survey area covers about *x*- and *y*-direction is

The total magnetic field anomaly and the IRF are computed following Bhattacharyya (1964). All simulations are performed with in-house code implemented in Python using the Scipy software package (www.scipy.org).

We first describe the IRF used in our simulations (Figure 3a). Assuming the north direction parallel to the *x*-axis, we compute a discrete IRF corresponding to a geomagnetic field with inclination and declination of 60° and 0° north, respectively, and sampling points on the same *f* is given in amperes/meter as intensity of magnetization, and the resulting IRF is expressed in nanoteslas, the unit of magnetic anomaly intensity. The magnetic field is obtained by convolving the simulated object with the IRF.

Next, the magnetic field is perturbed by additive white Gaussian noise, with

We perform the inversion by means of the nonregularized PLM and the nonregularized ISRA algorithms (Figures 4a and 5a). For the regularized algorithms (Figures 4b–4d, 5b–d, 6, and 7), we first choose the best value of *α* and *β*, minimizing the relative reconstruction error. Table 1 gives the reconstruction error defined in equation 30 for different values of *α*; the best result is about 32%. Figure 4b–d shows the reconstructions provided by regularized PLM, and Figure 5b–d shows the reconstructions provided by regularized ISRA. Relative reconstruction errors as a function of the number of iterations, computed with optimal parameters, are shown in Figure 6a and b. The modified ISRA algorithm for Tikhonov and TV regularization has a very flat minimum. Thanks to the flatness of the minimum, generally we can stop the iterations before the minimum, obtaining practically the same reconstructions.

### Burial depth and depth extent errors in IRF estimation

In the single-layer approximation of the magnetic problem, the IRF depends on the depth extent and burial depth of the object. The semiblind deconvolution method minimizes reconstruction errors by using a rough estimate of the IRF. The true object is the same used in the exact IRF simulation (Figure 3b). We start from an estimate of the IRF in the context of an archaeological survey by assuming a 1-m-depth-extent and 1.5-m-burial-depth single layer. Class A consists of IRFs corresponding to different burial depths of the top of the layer, one every

Figure 8 shows the result using

### Field example

Finally, we applied the inversion methodologies to a real data set (Figure 9a). The area investigated is located in Burnum, Croatia, and covers

## DISCUSSION

### Exact IRF

The best reconstructions provided by nonregularized PLM and modified ISRA are similar, but they do not fit the true solution well. In fact, the structure is too sharpened and the values of magnetization are not very close to those of the true object.

Better reconstructions can be obtained using the regularized algorithms. The regularized inversions produce reconstruction errors slightly smaller than those of the nonregularized ones. However, their reconstructions are much better from the visual point of view. Both methods are quite slow with respect to other algorithms (PLM and ISRA provide nonegative solutions and for this reason can be superior, for instance, to conjugate gradients, even if the conjugate gradients algorithm is much faster). However, slowness is irrelevant because in real applications we must stop the iterations at a certain point; the fact that the solution varies little, iteration after iteration, simplifies the stopping process. The solution provided by Tikhonov regularization is close to the nonregularized one in the PLM case, but in the ISRA case the reconstruction has less energy and hence is more similar to the true solution.

The TV functional clearly demonstrates the best performance in shape similarity and magnetization values. In particular, TV enhances the blocky structures that characterize the buried object. By inspecting the recovered objects, it is clear that the result obtained by TV is closest to the shape of the true object. The MGS functional in this case behaves similarly to the TV functional, providing a more focused image with respect to the Tikhonov and nonregularized results.

A further comparison of the different regularization methods is provided in Figure 7a and b. Profiles A and B (Figure 7c) reveal the behavior of different a priori information on the solution. The edge-preserving TV and MGS functionals provide a better approximation of the true shape of the object than Tikhonov regularization. However, in this case, TV achieves the best result, showing a flat structure in the two profiles.

### Semiblind deconvolution

The IRF selected by the algorithm is given by a 1.5-m-thick, 1.0-m-deep layer, corresponding to the IRF used to generate the synthetic data. In this case, the method can retrieve the true IRF, enhancing the inversion results and reducing the relative reconstruction error. The resulting deconvolved image (Figure 8) is quite similar in shape and in terms of boundary detection to the results obtained with exact IRF. The magnetizations are also quite close to those retrieved with exact IRF. However, this result requires estimating a number of parameters (*β*, *α*, and the number of iterations between projections onto the admissible set of the IRFs), which may be difficult in a real case where the values of the magnetization could be biased. This suggests that if the aim of inversion is oriented mostly to trace the position and edges of the buried objects, then the single-layer model also works well in the case of incorrectly estimating burial depth or depth extent.

### Field example

The resulting deconvolved image shows a set of linear structures that fit well in the Burnum case, where buried walls have been found (Figure 9b). This image enhances and delineates the structures only slightly visible in the acquired data set, confirming the simulation experiments. In particular, some orthogonal lineaments can be associated with the buried walls. However, the acquired data are quite noisy, so some noise is present in the inversion. The recovered IRF gives a top-of-layer burial depth of

## CONCLUSIONS

Assuming that buried structures may be modeled by a single layer of constant burial depth and depth extent, interpretation of magnetic data over 2D arrays is equivalent to the problem of image reconstruction. Our methodology differs from standard filtering/inversion techniques for two reasons. First, we use a constrained minimization that accounts for appropriate a priori information. Second, with the semiblind deconvolution method, we can treat the case of unknown impulse response function. The two iterative reconstruction methods, PLM and modified ISRA, work well for a single-layer approximation in inverse magnetic prospecting. However, PLM and modified ISRA do not take into account remanent magnetization, which could affect the buried material and possibly buried structures of different depths and thicknesses.

Tests comparing standard Tikhonov and other edge-preserving TV regularizing functionals reveal that the other functionals are the best choice of a priori information if the target is represented by piecewise constant (blocky) magnetization distribution (e.g., foundation structures, roads, graves). In the case of perfect knowledge of the IRF, the algorithms behave as in standard image restoration problem, even if, in the magnetic problem, the IRF takes negative values. The TV functional clearly outperforms the other algorithms. The FFT improves speed and reduces the size of the matrices; however, the rate of convergence can be ameliorated by applying acceleration techniques.

By applying our semiblind deconvolution method to numerical simulations, we verify the robustness of these algorithms with respect to errors in the magnetic IRF as normally happens in real cases. The semiblind approach damps the errors caused by the initial estimated IRF and allows a semiautomated estimate of the unknown burial depth and depth extent of the layer. Hence, the semiblind method is an important tool for archaeogeophysical prospecting. We performed a quantitative analysis by computing the reconstruction error in the cases we investigate. Next, we applied the semiblind deconvolution method to a real data set. The result provides an enhanced image of the structures buried in the subsoil, illustrating the effectiveness of our approach.

## ACKNOWLEDGMENTS

We thank Barbara Frezza (Università degli Studi di Siena-Area di Archeologia Medievale) for providing the data set acquired in Burnum. The authors are grateful to Roman Pasteka, René E. Chavez, an anonymous reviewer, and the associate editor for their helpful suggestions that improved the manuscript.

- Received July 31, 2008.
- Revision received January 8, 2009.