- Society of Exploration Geophysicists
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.
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 (Figure 1). Subtracting the International Geomagnetic Reference Field (IGRF) , we obtain when , applying the first-order Taylor's formula of the square root (Blakely, 1995). The relation between the projection of the anomalous magnetic field , with , on the earth field (supposed constant within the survey area) and the magnetization of the sources, located in the domain with , is described by (Li and Oldenburg, 2000): (1) where the ⋅ indicates the dot (row-column) product. The integrand of equation 1 can be written in the form , where (2) We assume the magnetization is parallel to the inducing field, i.e., , with and m the intensity of magnetization. By projecting on , from equation 1, we have a 3D convolution (Bhattacharyya and Navolio, 1975): (3) where is the projection of the anomalous magnetic field onto the geomagnetic one, , with T denoting transposed.
The single-layer model (Tsokas and Papazachos, 1992) assumes that the sources belong to only one layer, with the top at and bottom at . We refer to as the burial depth and h as the depth extent of the layer. The intensity of magnetization does not depend on inside the layer and is zero outside, so the intensity function is . Finally, all of the measurements are taken on the plane . Under these assumptions, equation 3 becomes (4) where (5) is the so-called impulse response function (IRF).
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 , and taking the measurements on a regular grid on a plane parallel to this layer, the inverse problem is expressed by the discrete convolution product: (6) where the unknown of the problem is the intensity of magnetization of the prism at and where is the discrete IRF, i.e., the total magnetic field produced at point by a prism of unitary magnetization at .
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 to ), considerably increasing processing efficiency.
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.
The maximum likelihood approach, in the case of additive Gaussian noise, is equivalent to the minimization of the least-squares functional: (7) where ‖⋅‖ is the -norm and is the convolution defined by equation 4. In a least-squares sense, the functional is a measure of the goodness of fit of the model predictions to the actual data.
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 : (8) where indicates the reflection of k with respect to the origin, i.e., . The classical algorithms are gradient based; they start from an image that is modified iteratively in the direction of .
Projected Landweber method
Generally, the iterative projected Landweber method (PLM) is used to minimize the functional on the nonnegative closed cone (denoted by C). Being a gradient method, it has the form (9) where is the projection on 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): (10) where the multiplication of arrays is in the Hadamard sense, i.e., pixel by pixel. According to SGM, we look for a decomposition of the gradient in the form (11) where and ; and are not uniquely defined. Then the first Karush-Kuhn-Tucker condition can be written as a fixed-point equation: (12) If we introduce the definitions and (13) then, by inserting and in the expression for (equation 8), the gradient becomes (14) Therefore, the gradient decomposition of equation 11 can be obtained by the following functions: (15)
Next, we apply the successive approximation method to the fixed-point equation 12 to obtain the modified ISRA algorithm: (16) This iterative algorithm can be implemented using the convolution theorem, exploiting the speed of FFT.
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 (17) where comes from the negative log of the prior and α is a positive multiplier called regularization parameter (Bertero and Boccacci, 1998).
This procedure is called regularization of the inverse problem, and is called the regularizing functional, prior, or penalty. More generally, a priori information is conveyed into the inversion algorithm through the choice of and associated parameters. The most common form of , based on the Tikhonov regularization method, is given by (18) which ensures the solution has the minimum energy. This kind of regularization is used by Tsokas and Papazachos (1992), even if it is presented as a Wiener filter for inverting magnetic data for archaeological prospecting.
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) (19) as in the case of the total-variation (TV) functional (Vogel, 1997), for which (20) The parameter β 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), (21) applied to the square norm of the gradient leads to another suitable functional, already developed in a geophysical framework by Portniaguine and Zhdanov (1999) and called the minimum gradient support (MGS). Generally, for any function ψ, the gradient of takes the form (22) where is the first derivative 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 (23) To avoid numerical instability in the case of TV and MGS regularization, the iterative implementation of this algorithm requires computing the step length τ. We apply a dichotomic search of the minimum of (24) for (Vogel, 2002) and find the optimal τ.
Regularized modified ISRA
Regularization of the modified ISRA takes the form (25) in which and are the result of applying SGM to the gradient of . For example, splitting the gradient for the TV functional is (26)because we are considering nonnegative f and the function is increasing.
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 for each given image to deblur. The problem arises because we must know to compute the optimal regularization parameter. The value α is a weight that balances between the minimization of the discrepancy and that of the functionals. In the simulation case, we choose α 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 (27) we stop the iterative process when σ becomes smaller than the estimated standard deviation of the noise in the acquired data.
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 , as shown by equation 5. By fixing h, we have a one-parameter family of IRFs, and the deconvolution problem becomes semiblind. Hence, we consider k depending on the parameter , and we write in place of . We propose a semiblind iterative deconvolution that is a modification of the iterative blind deconvolution introduced in Bertero et al. (1998) and is based on regularized versions of PLM.
Initially, we compute a discrete set of admissible IRFs, corresponding to a given set of burial depth values: (28) where and are the minimum and maximum values of the burial depth . Each external iteration of the semiblind approach consists of two steps. For the first step, represented by the upper rectangle in Figure 2, we use , to denote the output of iteration j. This pair enters in the object box. The deconvolution method implemented in this box is the iterative PLM with TV regularization. By keeping as the IRF and as the initial model, we perform one inner iteration. In this way, we have an update of the object.
Then the pair , enters the IRF box, where the iterative method corresponding to PLM with Tikhonov regularization is implemented. Here, the roles of object and IRF are exchanged: is keep fixed and the deconvolution is applied to the IRF, using as the initial guess. In this box, five iterations are performed, and the output of the first step is the pair , .
In the second step of the iteration (lower rectangle, Figure 2), if , then no operation is performed and we set . If , then we take as the IRF in the set A such that (29) i.e., the one with a minimum distance from . In other words, the projection becomes effective every 250 iterations. However, the IRF changes before the next projection as the iterations proceed (first step box of Figure 2). This number of iterations before projection (250) is chosen for this particular data set after trial and error. However, the number could differ for data sets, so we must select it.
To initialize the procedure, we take as the IRF corresponding to an estimate of , and . The number of internal iterations of the two boxes is optimized by means of numerical simulations, and the number of global iterations is determined by the minimal error (in simulations) or by the discrepancy principle.
The algorithm guarantees that the IRF is physically meaningful; moreover, it automatically estimates the burial depth of the layer.
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 (30) where j is the iteration number. It is a measure of the distance between the solution at the iteration and the true object (Figure 3b).
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 , and its burial depth is (top of the object). The object extends on the horizontal plane about in width and in length (Figure 3b) and is zero padded for matching the number of samples of the survey area.
We divide the subsoil into a set of homogeneous vertical prisms, considering only a single layer of cells. The survey area covers about , and the sampling spacing in both x- and y-direction is , resulting in a grid. The sampling points are located over the buried prism centers.
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 -grid of the object. Note that 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 corresponding to about error in the measurement, a reasonable value for modeling the read-out noise in the acquisition process. In Figure 3c we show the blurred and noisy image, which takes values from −10.4 to 23.2 nT and has negative values in 57% of pixels.
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 , starting from 0.5 to 4.0 m. This guarantees that the IRF used in the deconvolution for the unknown object has a form very close to that of an actual IRF.
Figure 8 shows the result using for the IRF and , for the object. The relative reconstruction error is 40% after 2327 iterations, where it reaches the minimum. We obtain the correct value of the burial depth, .
Finally, we applied the inversion methodologies to a real data set (Figure 9a). The area investigated is located in Burnum, Croatia, and covers . The sample spacing is . We adopt the semiblind deconvolution approach because we do not know the IRF exactly, and we have only an approximate estimate of the burial depth. If the solution belongs to a -thick single layer, we consider the IRF given by a layer at burial depth. The class of admissible IRFs is given by a set computed for burial depth ranging from every . We choose for the IRF and , for the object. Obviously, in this case we cannot compute the relative error, so the iterations are stopped with the discrepancy principle (equation 27), considering . The algorithm chooses the IRF corresponding to a layer buried at (top of the layer) after 398 iterations.
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.
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.
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 , a plausible value for Burnum. However, the estimated depth extent is inexact, so the burial depth of the layer can be biased.
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.
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.