Abstract
A new technique is proposed to process 2D apparent resistivity datasets, in order to obtain a fast and contrasted resistivity image, useful for a rapid data check in field or as a starting model to constrain the inversion procedure. In the past some modifications to the back-projection algorithm, as well as the use of filtering techniques for the sensitivity matrix were proposed. An implementation of this technique is proposed here, considering a two-step approach. Initially a damped least squares solution is obtained after a full matrix inversion of the linearized geoelectrical problem. Furthermore, on the basis of the results, a subsequent filtering algorithm is applied to the Jacobian matrix, aiming at reducing smoothness, and the linearized damped least square inversion is repeated to get the final result. This fast imaging technique aims at increasing the resistivity contrasts and practically, since it does not require a parameter set optimization, it can be used to easily obtain fast and preliminary results. The proposed technique is tested on synthetic data, the objective of which is to find the optimal parameter set. Finally, two field cases are discussed and the comparison between back-projection and inversion is shown.
1. Introduction
Besides the well-known inversion techniques there are other techniques that aim to obtain approximate representations of the subsoil in shorter times compared to those required by the “classical” iterative methods.
These methodologies give approximated images of the resistivity patterns in the subsoil, in which zones of high or low resistivity are approximately represented. However the resulting resistivity values are generally very different from the real ones and the resistivity gradient is generally much lower than the corresponding gradient obtained by the inverted models.
A probabilistic methodology, called probability tomography, was developed by Patella [1] for the SP method and then extended to resistivity measures [2], in which the interpretative model obtained is an image reconstruction of the most probable location of electrical charges induced by the primary source over buried resistivity discontinuities.
In previous papers [3–5] a filtered back-projection resistivity technique was described that was based on a convolution of the experimental data with a filter function calculated on the basis of the sensitivity factors of each voxel on the various resistivity data. The calculation of the influence factors in a discrete set of representative points (centers of the voxels, to be utilized for the back-projection) is based on the integral of the influence density function on the volume of each voxel [6].
2. The Back-Projection Resistivity Technique (BPRT)
The back-projection resistivity technique (BPRT) can be applied to a set of apparent resistivity measures to quickly obtain an approximate image of the resistivity distribution of the investigated volume. This technique is based on the consideration that a resistivity perturbation in a point element (voxel) of a bounded region produces a change in voltage—thus an apparent resistivity anomaly—at the surface of the region, according to a sensitivity coefficient. The value of the coefficient is dependent on the position of the voxel considered in respect of both the current and the voltage dipoles, in agreement with the sensitivity theorem of Geselowitz [7]. This consideration suggests that it is possible to correlate all the measured resistivity values, weighted by the appropriate sensitivity coefficients, to each voxel of the investigated volume [8] and to estimate the resistivity value of each cell of the model using a weighted summation of the apparent resistivity measurements. This method allows us to identify zones of “high” and “low” resistivity in the subsoil [9]. The main disadvantage of these techniques is that the resistivity values of the resulting image can greatly differ from the actual subsurface resistivities.
Barber et al. [10] developed a method of back-projection of measures along equipotential lines. This approximation simplifies the reconstruction process by assuming a linear relationship between electrical measurements and the resistivity distribution. This results in a fast single step image reconstruction algorithm. The image is formed by the back-projection of the boundary measurements normalized to the case of uniform resistivity. Barber and Brown [11] modified the reconstruction by normalizing both sides of the back-projection equation using a reference data set. The consequence of this strategy, however, is that the reconstruction produces differential images rather than static images and thus only the normalized change in conductivity can be determined.
Barber [12] derived a new formulation for the back-projection reconstruction algorithm, in which the measurement data are first transformed and then back-projected. The filtration stage of the “filtered back-projection” method is achieved by premultiplication of the data measured and is followed by the back-projection operation.
Kotre [13, 14] used a linearized sensitivity matrix to carry out a similar data normalization strategy. However, rather than normalizing this matrix with the uniform reference data set, Kotre [13] normalized each column of the matrix (corresponding to each individual cell) with the sum of the sensitivity coefficients for all apparent resistivity measures.
Let us consider a set of measures carried out on the surface of a homogeneous medium characterized by a background resistivity . All data are referenced on the surface by means of the positions of the four electrodes used. The investigated volume can be discretized in small voxels. If the th voxel has an anomalous resistivity , all the apparent resistivity measures are influenced. In particular, the value of the th resistivity measure will be . As long as is small, it can be assumed that the relation is linear, so that .
The proportionality can then be defined by a sensitivity coefficient that correlates the th voxel to the th measure: By considering the resistivity anomaly in the single element in terms of logarithmic difference between the perturbed value of the element and the background resistivity, we establish that perturbation can be expressed as where .
By taking the superposition principle into due consideration, the relative change in the apparent resistivity is given by the sum of the contributions of all the voxels [13, 15]: The perturbation of the apparent resistivity in respect of a homogeneous body can then be expressed as Using the matrix notation we can write where is the vector () of the relative differences of apparent resistivity, is the vector () of the relative resistivity changes in the model elements, and is the matrix of the following terms: The BPRT is based on the estimation of the resistivity of each cell of the model by means of the approximate inversion of (5).
3. Procedure Proposed
The procedure proposed in this work consists of four steps:(1)evaluation of sensitivity matrix ,(2)inversion of matrix using a damped LSQR solution,(3)recalculation of a filtered Jacobian matrix obtained by means of a correlation filter,(4)inversion of the filtered sensitivity matrix.
3.1. Evaluation of the Sensitivity Matrix
In the case of a homogeneous half-space, can easily be calculated by estimating terms by integrating the sensitivity function, derived from the Frechet derivate [16], over each voxel.
In fact, if the pole-pole array is considered with a current electrode and a potential electrode , a small change in resistivity within an elemental volume having coordinates causes a potential variation in given by [17] in which the term within the integral is the Frechet 3D derivative. The Frechet 3D derivative for a four-electrode array can be obtained by adding together the contribution of each potential electrode-current electrode pair.
Equation (7) has the same form as the equation found by Roy and Apparao [6] to calculate the potential difference created by a small volume element centered at with a pole-pole array on the surface of a homogenous earth.
In this way only depends on the positions of the electrodes and on the discretization model.
Between the electrodes in near surface zones, the sensitivity function can be negative [18]. So when the size of the cell is smaller than one-half of the unit electrode spacing, the terms of can oscillate between positive and negative factors with very large absolute values near the electrode. As a consequence, the BPRT becomes instable near the electrodes.
To partially overcome this problem [5] the normalization of the terms of matrix is made by dividing each term by the sum of the absolute values of all the coefficients related to the th measure.
Consider In this way the effects of the oscillation of the sensitivity function near the electrodes are strongly attenuated without there being a significant loss of resolution.
3.2. Inversion of the Sensitivity Matrix
Inversion of the sensitivity matrix in (4) constitutes the basis of the BPRT when estimating . However, the analytical solution of (5) is not trivial and the process is often instable. Approximate solutions [5, 13, 18] have been proposed that substitute the inverted matrix with the transposed matrix . Figure 1 shows an example of back-projection obtained using . Apparent resistivities have been calculated by simulating a resistive square body (1000 m) buried in a homogeneous medium (100 m). It stands to reason that back-projected images are also smoothed. Furthermore, the estimated tomographic resistivity contrasts are generally very low in respect of the real ones. Since matrix is generally ill-conditioned and is not square, a more accurate solution for the inversion adopted here is a damped least-squares solution [19–21]: where is the initial homogeneous background resistivity vector, is a damping factor, is the biggest value of the main diagonal values of the matrix [], and is a unity diagonal matrix of the same size as []. In this way the square matrix [] is regularized to obtain an approximate solution.
For a nonzero value of an inverse can be calculated, although its condition number and hence the stability of the inversion depend on the value of damping factor : for large values the image obtained can result as being too smooth and blurred; however, if is too small, inversion can produce very noisy images.
Figure 2 shows the results of the one-step damped least-square inversion of synthetic data calculated for a model having two square bodies with resistivity equal to 1000 m buried at different depths in a homogeneous medium of 100 m resistivity. Synthetic tests showed that for , the back-projected images were generally noised and furthermore, artifacts are present due to the instability of the process. As increases noise attenuates, while the resistivity gradient decreases. Therefore, if damping is too high () the image becomes too smooth and blurred. All the synthetic tests carried out with 2D data suggested that a good compromise for 2D inversions is to choose .
(a)
(b)
(c)
3.3. Recalculation of the Sensitivity Matrix Using a Correlation Filter (Gradient Enhancing)
The subsequent step of the imaging technique proposed is a recalculation of the sensitivity matrix . This is substituted by a matrix in which each term is calculated by applying a correlation filter: The aim of filtering is to attenuate term linking the th voxel and the th resistivity datum when the first-step resistivity value of the th voxel and the apparent resistivity value of the th datum show a low correlation.
In order to correctly estimate a correlation value between data sets of different kinds, both sets of experimental data (apparent resistivities and first-step resistivities) are previously normalized using the following normalizing function: where is the datum (it can be or ) and and are, respectively, the minimum and the maximum values.
Figure 3 shows the pattern of the normalizing function applied to a set of data measures and to the first-step resistivity values. In this way ranges from −1 (for a minimum resistivity) to 1 (for a maximum).
The filtering function is then defined as follows: where is a correction parameter that controls the grade of attenuation of terms .
In this way varies from (if the correlation between and is the lowest) to 1 (if the correlation is the highest).
Figure 4 shows the pattern of the filtering function for different values of parameter . The lower the correlation, the lower is the value of (the attenuation being higher). Furthermore the higher the value of , the higher is the gradient of attenuation.
3.4. Inversion of the Filtered Matrix B
Finally each sensitivity coefficient is recalculated by means of (10) and using the damped least-squares inversion, now considering the filtered matrix instead of : In Figure 5 the influence of varying correction parameter on the back-projection is shown. The sharpness of the new tomographic image depends on the value of correction parameter : generally, if ≤ 1, the image obtained is generally too smooth. As increases the anomalous bodies become more contrasted. However, when is greater than 10, images become noised and the anomalies break down in several parts.
(a)
(b)
(c)
All the synthetic tests carried out with 2D data suggested that a good compromise for 2D inversions is to choose .
In spite of the long presentation, the four steps are in fact very quick to carry out when using an automatic procedure with a simple notebook.
4. Results of Some Tests on Synthetic Models
The procedure proposed was tested on synthetic data obtained by simulating pseudosections of different arrays carried out on several simple 2D discrete models. All the synthetic tests stressed that the estimated shapes, depths, and resistivity contrasts of the anomalous bodies depend on the kind of array used. Results seem to suggest that the best estimation is obtained by using the Wenner-Schlumberger array. Hereafter, only the results obtained using this latter array will be discussed.
For each model the minimum electrode distance was fixed at 1 m and apparent resistivity measures were calculated considering that order varied between 1 and 10 and that the potential electrode distance MN was equal to 1, 2, and 3 m.
Several tests were carried out by changing the dimensions, the position, and the depth of the targets and by varying the values of the above defined parameters.
The results obtained by varying the depth (from 1 to 3 m) of a square target (3 m × 3 m) are shown in Figure 6. In a homogeneous medium ( = 100 m), a resistive (1000 m, Figure 6(a)) or conductive (10 m, Figure 6(b)) squared prism is placed under the profile at varying depths. The resolution and the contrast of the obtained image obviously decrease with depth. As a consequence the anomalous body is still recognizable until a depth of about 2.5 m.
(a)
(b)
Figure 7 shows the tomographic images obtained starting from square-sized targets in which the resistivity of the targets varies from 2 m to 5000 m. In cases of very low resistivity the tomographic shape does not fit the real shape. The correlation between the tomographic images and the models increases with the resistivity of the target. Specifically, the difference between the real top of the target and the estimated one decreases.
Finally, the results of a more complex model are shown. Apparent resistivity measures were simulated starting from a model consisting of four prisms of different shapes and resistivities, buried in a homogeneous medium. Figure 8 shows a comparison between back-projection using transposed matrix (Figure 8(a)), the image obtained by the technique here proposed using a damped least-squares solution of the filtered sensitivity matrix (Figure 8(b)), and two inversions carried out using the RES2DINV software, respectively, with a standard least-squares constrain (Figure 8(c)) and with a robust data constrain algorithm (Figure 8(d)). Results of back-projections are strongly smoothed and dependent on the lateral position of the anomalous body in respect of the section. Furthermore, the resistivity contrast is very low if compared to the real one. Instead, the damped least-squares solution of the filtered sensitivity matrix allows a less smoothed and more accurate representation of the targets, which show higher resistivity contrasts and more correct shapes, especially if referred to the resistive targets. The resistivity range showed is considerably wider than that of the back-projection, and it is near the range of values obtained in the inversions. Nonuniqueness involved in standard least-squares constrain procedure implies that inverted models can present artifacts linked to the choice of inversion parameters. The image obtained by the damped least-squares solution of the filtered sensitivity matrix is instead less influenced by the choice of parameters and does not produce significant artifacts.
(a)
(b)
(c)
(d)
5. Test on Experimental Measures
The reliability of the proposed method was tested in several archaeological and environmental cases. The results of two of these are discussed below.
Figure 9 shows a test on a set of geoelectrical measures carried out on the plane of Himera [22] to detect buried archeological structures. A 2D dipole-dipole resistivity survey was carried out using 64 electrodes, spaced 0.5 m apart, and with order from 1 to 9.
(a)
(b)
(c)
ERT results showed clear resistive localized anomalies, which are connected with wall-like relics. This was confirmed by an archeological excavation near the profile, in which the structure of a building in line with geoelectrical anomalies was located. In fact, the survey direction was perpendicular to that of the buried archaeological structures. Furthermore, the low background resistivity values shown by the inversion were interpreted as representing river deposits having different moisture characteristics.
The comparison between the fast imaging and the robust inversion shows that the structures detected by the fast imaging technique have very similar shapes to those of the inversion.
The second case regards an area on the outskirts of Marsala (Western Sicily) where many sandstone quarries are present [23], some of which are abandoned. In this area the main outcropping formation is the Calcarenite di Marsala, which has been exploited for a long time as building stone.
A geophysical model of the subsoil was designed using detailed geological information and the geometrical shapes of some galleries that were physically investigated. This allowed us to compare synthetic and experimental data. In the survey here presented, the filtered Jacobian damped LSQR model shows three resistive anomalies corresponding to three buried tunnels; one of which was actually inspected. In Figure 10 the grey lines represent the positions and the shapes of the galleries. The fast imaging model (Figure 10(a)) shows three resistive anomalies in correspondence with the galleries, even though it shows a lower resistivity contrast and smoother boundaries compared with the inversion model (Figure 10(b)).
(a)
(b)
6. Conclusions
An in-depth discussion on this topic involves taking into consideration not only the type of synthetic and/or real (sharp or smooth boundary) model to be used, but also the equivalence problem, which is always relevant. However, the simulations and field data here presented stress that the filtered Jacobian damped LSQR imaging can be a valid tool when looking for a fast approach to the tomographic representation of geoelectrical data.
The modification of the back-projection equation proposed above indeed leads to tomographic images that have higher resistivity contrasts than those obtained using standard back-projection methods. Furthermore, filtering of the sensitivity matrix with the correlation filter proposed and tested above reduces smoothness in the final tomographic representation.
The values chosen for the parameters of the filtered Jacobian damped LSQR imaging (a damping parameter value 0.03 in the first matrix inversion and a correction parameter in the filtering procedure) worked sufficiently well with both synthetic and field data. In all those cases where sharp-boundary models were sufficiently representative (i.e., tunnels, voids, archaeological remains, etc.) the procedure here presented seemed to be efficient enough to give reliable images of the subsoil.
Further simulations on noised data (not presented in this paper) showed that back-projected images are not particularly influenced by noise. In fact, they tend to highlight only the main structures, while masking small anomalies that are probably caused by noise.
In conclusion, the quickness of the technique, together with its almost total independence from choice of parameters and from data noise, suggests that proposed procedure can be a useful as well as a fast technique of approximate image reconstruction. It can therefore be carried out directly in the field for a first and fast preliminary result. Furthermore, these could be used as a starting model, instead of the homogenous model or the pseudosection, to constrain the inversion procedure.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.