Next Article in Journal
Spatiotemporal Variations and Uncertainty in Crop Residue Burning Emissions over North China Plain: Implication for Atmospheric CO2 Simulation
Next Article in Special Issue
AEM in Norway: A Review of the Coverage, Applications and the State of Technology
Previous Article in Journal
The Identification and Analysis of Gas-Related Volcanic Features within Chang’e-5 Landing Region
Previous Article in Special Issue
3D Wavelet Finite-Element Modeling of Frequency-Domain Airborne EM Data Based on B-Spline Wavelet on the Interval Using Potentials
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

1D Stochastic Inversion of Airborne Time-Domain Electromagnetic Data with Realistic Prior and Accounting for the Forward Modeling Error

1
Department of Civil and Environmental Engineering and Architecture (DICAAR), University of Cagliari, 09123 Cagliari, Italy
2
Geological Survey of Denmark and Greenland (GEUS), 8000 Aarhus, Denmark
3
Department of Geoscience, Aarhus University, 8000 Aarhus, Denmark
*
Author to whom correspondence should be addressed.
Submission received: 2 September 2021 / Accepted: 24 September 2021 / Published: 28 September 2021
(This article belongs to the Special Issue Airborne Electromagnetic Surveys)

Abstract

:
Airborne electromagnetic surveys may consist of hundreds of thousands of soundings. In most cases, this makes 3D inversions unfeasible even when the subsurface is characterized by a high level of heterogeneity. Instead, approaches based on 1D forwards are routinely used because of their computational efficiency. However, it is relatively easy to fit 3D responses with 1D forward modelling and retrieve apparently well-resolved conductivity models. However, those detailed features may simply be caused by fitting the modelling error connected to the approximate forward. In addition, it is, in practice, difficult to identify this kind of artifacts as the modeling error is correlated. The present study demonstrates how to assess the modelling error introduced by the 1D approximation and how to include this additional piece of information into a probabilistic inversion. Not surprisingly, it turns out that this simple modification provides not only much better reconstructions of the targets but, maybe, more importantly, guarantees a correct estimation of the corresponding reliability.

Graphical Abstract

1. Introduction

Airborne time-domain electromagnetics (ATEM) was originally developed for mineral prospection [1,2,3,4] and is now a very consolidated methodology. Its first application probably dates back to 1948 [5], but since then, many different helicopter and fixed-wing acquisition systems have become commercially available. Amongst the most internationally known examples are: EQUATOR [6,7,8], Xcite [9,10], VTEM [11,12], SkyTEM [13], AeroTEM [14], MULTIPULSE [15], MEGATEM/GEOTEM [16,17], TEMPEST [18,19], and Spectrem [20].
In the last decades, technological advancements made it possible to move from the simple mineral target detection to more detailed groundwater mapping [21,22] and geological modeling applications [23,24].
Together with the enhanced capabilities of the instrumentation, the strategies to deal with the collected data have also been largely improved. For example, similarly to what happens in seismics (e.g., [25]), the stacking of the recorded transient curves—by means of moving windows whose widths depend on the time-gate [26,27,28]—can be used to preserve high lateral resolution at shallow, while increasing the reliability of the data at depth; this is, indeed, compatible, at the same time: (i) with the higher signal-to-noise at the early-gates, and also (ii) with the larger spatial footprint at the late-gates.
Moreover, novel strategies to effectively incorporate prior information into the translation of the observed data into conductivity models have led to the development of, for example, spatially constrained inversion schemes [29,30,31,32] (again, similarly to approaches utilized for other, very different, kinds of data [33,34,35,36,37]). This capability of enforcing spatial coherence allowed the reconstruction of (pseudo-)3D conductivity distributions even by means of simple 1D forward modeling [38,39].
Clearly, approaches based on 1D forward modeling are particularly convenient for their computational performances [40]. Recent strategies exploit the neural networks potential, allowing almost real-time inversions with no significant quality reductions [41]. In addition, in the attempt of retrieving, not merely the conductivity distribution, but, rather, immediately useful pieces of information about the targets and, at the same time, in the effort of supplying reliable assessment of the associated uncertainty, probabilistic petrophysical inversions are becoming more and more common [42].
Unfortunately, as mentioned, all these approaches rely on 1D approximation of the physical phenomena. Very often, the error introduced by using this rough description of the reality is neglected, and when, in the most conservative workflows, it is not, the modeling error is accounted for by arbitrarily increasing the data fitting threshold [43,44]. Of course, disregarding this further noise component during 1D inversions could potentially be extremely dangerous as fitting the 3D observations with 1D forward modeling might produce reconstructions appearing very detailed, and, seemingly, characterized by high resolution, but that, in reality, are affected by correlated noise artifacts [45]. The necessity of tackling three-dimensionality issues in effective and pragmatic ways is even more crucial nowadays as large datasets are (and are going to be) collected for (ultra-)shallow mapping [46] where even sedimentary environments can be characterized by extremely high heterogeneity; in such a case, 2D/3D artifacts, due to improper modeling error assessment, would pervade the conductivity (and/or petrophysical) reconstructions without being detrimental in terms of the data fitting levels [47]. Consequently, the results could nicely fit (in a 1D sense) the (3D) observed data and potentially seem very accurate, but actually be permeated by pitfalls.
In principle, the most straightforward way to deal with these dimensionality problems would consist of using fully 3D forward modeling algorithms. Indeed, much research is currently devoted to practical 3D inversion approaches [1,39,48,49,50,51,52,53,54]. However, 3D approaches are still too computationally demanding for being routinely used [55], especially within probabilistic frameworks [56].
In the present research, we show that the 1D modeling error can be significant and lead to erroneous/misleading inversion results. In addition, we discuss how such modeling error can be quantified based on the available prior information about the system under investigation. Once quantified, this resulting noise model could be incorporated, in principle, in both deterministic and probabilistic inversions to correctly account for the 1D approximation. Finally, we demonstrate how using such a correlated noise model positively impacts even the probabilistic inversion of ATEM data.
We test the proposed approach on synthetic datasets characterized by increasing complexity levels. Through them, we demonstrate the effectiveness of the proposed approach: (i) in reconstructing features that are indeed in the 3D data, and (ii) in providing the crucial, associated uncertainties. For each of the considered Tests, the results obtained via: (i) a “standard” deterministic 1D Occam’s inversion, (ii) a stochastic inversion with realistic prior, and (iii) the same stochastic inversion, but now taking into account the modeling error consistent with the used prior are compared.

2. Methodology

The calculation of the response of a physical system is always characterized by some level of modeling error. For example, even in the case of very sophisticated 3D forward modeling tools used to calculate the d B / d t responses of an electrical conductivity distribution caused by the excitation induced by a “perfectly” described ATEM system (is it really possible to characterize a priori an acquisition system? [8,13,44]), the parameterization used and the size of the discretized domain might affect the retrieved response. In this respect, every time we use a 1D forward modeling approximation [43] for the inversion of ATEM data (inherently 3D), we introduce some errors that need to be handled; neglecting this source of additional uncertainty would inevitably lead to artifacts, paving the road to successive geological misinterpretations. In the best scenario, the magnitude of the modeling error will be negligible compared to uncertainty in the data, and can be ignored [57]. More often, modeling errors will be significant, and must be taken into account.
Here, we briefly recall a formal framework in which the modeling error is described through a multivariate normal correlated probability distribution, which can be naturally used when the measurement errors are also Gaussian [57]. In probabilistic formulations of inverse problems, the goal is to retrieve the posterior probability density function p ( m | d ) measuring the probability of having the model m compatible with the measurements d . In accordance with Bayes’ theorem, p ( m | d ) is proportional to the product between the prior probability density function of the model parameters p ( m ) , and the conditional probability density function p ( d | m ) . Hence, p ( m | d ) p ( d | m ) p ( m ) , with p ( d | m ) connecting the measured data and the model parameters, and, in the specific case of a Gaussian noise distribution, can be written as
p ( d | m ) = k d exp ( 1 2 ( d F ( m ) ) T W d T W d ( d F ( m ) ) ) ,
where: (i) k d is just a normalization factor, (ii) F is the forward modeling operator used during the inversion, and (iii) W d is related to the data covariance C , and, often, by assuming mutually independent data, can be considered equal to W d = diag ( σ d ) 1 =   C 1 / 2 where the i -th component of the vector σ d is the standard deviation of the i -th data component (in the specific case of the ATEM data, [ σ d ] i is the standard deviation of the d B / d t value associated with the i -th time-gate).
If the model parameters are also assumed to follow a Gaussian distribution, then the prior information about the solution can be formalized as follows: p ( m ) = k m exp ( 1 2 ( m m 0 ) T W m T W m ( m m 0 ) ) , in which: (i) k m is another normalization factor and (ii) the Gaussian is centered on the reference model m 0 . In this specific case of Gaussian distributions assumed both for the data noise and the model parameters, we obtain that the maximizer of the probability p ( m | d ) is also the minimizer of the regularized inversion objective functional W d ( d F ( m ) ) L 2 2 + W m ( m m 0 ) L 2 2 . In fact, this is a very well-known result (e.g., [36,58,59]), to some extent, reconciling probabilistic and deterministic approaches; in particular, if C m 1 is taken equal to λ 2 L T L with λ being the Tikhonov parameter controlling the relative importance of the regularization term with respect of the data misfit, and L a discrete approximation of the spatial derivative—then, the minimization of the objective functional coincides with the standard Occam’s inversion [60].
However, the approach discussed in the present research loosens several of these ansätze, and, in the following:
  • we do not restrict ourselves to the Gaussian assumption for the model parameters distribution p ( m ) as we are going to consider quite general prior distributions defined through the realizations of those distributions and that will be generated via a geologically informed procedure;
  • the W d =   C 1 / 2 will not consist uniquely of the component attributable to the noise in the observations, but it will also include a term incorporating the modeling error. In particular, the modeling error will be assumed to be consistent with a Gaussian probability density N ( d Δ ,   C Δ ) defined by the mean d Δ and the covariance C Δ . Hence, the p ( d | m ) in Equation (1) will have now the following expression [47]
    p ( d | m ) = k d exp ( 1 2 ( d + d Δ F ( m ) ) T W Δ T W Δ ( d + d Δ F ( m ) ) ) ,
    with W Δ = ( C + C Δ ) 1 / 2 .
By construction, as it will be detailed in what follows, these new terms d Δ and C Δ will also depend on the prior geological knowledge available about the investigated area.
It is important to stress that despite us demonstrating the effects of taking into account the modeling error within a probabilistic framework, the W Δ = ( C + C Δ ) 1 / 2 can actually be, in a very immediate way, incorporated into a deterministic framework as well. However, the evident advantage of using a stochastic approach is that we can naturally incorporate complex prior information (rather than enforcing simple—e.g., smooth or sharp—constraints) and those pieces of complex information need to be available in any case since they are used for the assessment of the modeling error.
Moreover, given the arbitrariness we can benefit from by defining the prior distribution directly via its samples, in general, we could potentially have the maximum flexibility and even use very powerful strategies as, for example, those based on Multiple-point statistics (MPS) approaches [61]; in those cases, the prior geological information can be formalized by means of the so-called Training Image (TI) that is, basically, representing the conceptual geological model of the expected target subsurface; MPS algorithms can generate samples of the prior that are statistically stationary with respect to the original TI and that can be used as detailed in the following Subsection. In the present research, however, we use other geostatistical strategies to populate the prior (and consistently estimate the associated modeling error). They will be described in the following Section 3 “Results”.

2.1. Estimation of Gaussian Correlated Modeling Errors

Here, we follow the strategy detailed in [47] to (i) simulate and (ii) quantify modeling errors caused by using a 1D forward as opposed to a full 3D forward for simulating ATEM data.
Firstly, a sample of the underling probability distribution representing the modeling errors is generated. This is done by generating a relatively large sample of N Δ realizations of the prior distribution p ( m ) as M = [ m 1 , m 2 , ,   m N Δ ] . The forward response is then calculated by using the approximate 1D forward model, F a p p , and the (assumed) exact 3D forward model, F e x . This provides a set of ‘approximate’ and ‘exact’ data in the form of D a p p = [ F a p p ( m 1 ) , F a p p ( m 2 ) , ,   F a p p ( m N Δ ) ] and D e x = [ F e x ( m 1 ) , F e x ( m 2 ) , ,   F e x ( m N Δ ) ] . Hence, the difference between the approximate and the exact forward models represents a realization, [ D d i f f ] i = [ D e x D a p p ] i , of the modeling error associated to the specific i -th sample of the prior. As will also be discussed in the following, to feed the 1D forward F a p p , unidimensional conductivity models have been extracted from the original 3D realizations of the prior in the locations just below the acquisition system positions.
Assuming that the modeling error can be characterized by a multivariate Gaussian distribution N ( d Δ ,   C Δ ) , the mean and covariance can be trivially computed from the samples D d i f f .
Finally, the assessment of the 1D modeling error can be plugged into Equation (2).
The number of realizations needed, N Δ , and the validity of the Gaussian assumption on the modeling errors will be addressed below.
In the present research, as the best approximation F e x , it has been considered an implementation of the forward modeling discussed in detail in [62]. The simulations have been preceded by a validation phase, in which the 3D forward modeling results have been compared against known solutions assumed to be exact. In our specific case, we performed preliminary tests of the 3D forward against semi-analytical solutions for unidimensional conductivity distributions; the chosen 3D simulation settings have led to mismatches of a few percentage points (generally around 5%). Consequently, as it will be shown later in the paper, the modeling error inherited from the 3D modeling has been assumed negligible with respect to the other noise sources and it has not been further considered in our analysis.

2.2. Inversion Strategies

Concerning the deterministic inversion results discussed in the following, as mentioned, we use an Occam’s inversion scheme [63] in which each individual d B / d t sounding is inverted independently from the adjacent ones and, therefore, the roughness operator in the regularization terms acts only vertically. To be fair, it is true that this specific kind of prior information, formalized by the stabilizer, is not in accordance with the investigated models (as it will be clear in the descriptions of the tests, characterized by abrupt conductivity changes). Nevertheless, the deterministic algorithm retrieves (smooth) models whose 1D responses are in almost perfect agreement with the inverted 3D d B / d t   data. The level of data fitting used for the inversion of the noise-free data is, for all tests, 0.01%.
The stochastic inversion consists of an independent extended Metropolis algorithm, which is analogous to the standard extended Metropolis algorithm [64], but with an infinitely long step-length; this choice implies that every new model proposal from the prior is a completely independent realization. Accordingly, its implementation consists of a slight modification of the SIPPI toolbox [65,66], making the proposed samples independently drawn from the prior ensemble. The 105 realizations of the prior are generated in advance since a subset of them needs to be used for the modeling error assessment; the appropriateness of this choice, and the convergency properties will be discussed in the section “Discussion”.

3. Results

In this section, we perform two synthetic tests of increasing complexity to investigate the effects of including in the inversion process: (i) the proper prior information, and (ii) an estimation of the modeling error. In both cases, we compare our results against a solution provided by a more standard 1D deterministic inversion.

3.1. Test 1: 3D Conductivity Distribution with Homogeneous Layers

The first test (Test 1—Figure 1) consists of a conductivity distribution mimicking possible glacial geological settings typical, for example, of Denmark [24,67,68] and characterized by an intricated network of paleovalleys.
Figure 2 compares, for the sounding locations considered in Figure 1, the 3D responses calculated for the entire 3D model against the corresponding 1D measurements that would be obtained by considering exclusively the 1D portion of the original distribution just below the acquisition position. The acquisition specifications are those of a typical VTEM system [11,12] (e.g., each sounding consists of 54 measurements). Not surprisingly, the 1D responses (Figure 2b) are characterized by abrupt lateral changes associated with the lateral variations of the conductivity model (Figure 2a), whereas, in accordance with the physics, the 3D d B / d t   data are much smoother (Figure 2c). Here, we treat the 3D calculated responses as the ‘observed data’ and invert them by means of both deterministic and probabilistic inversion methods.
In all cases we assume negligible measurement errors during the inversion, as the main purpose is to focus on the modeling error effects. In this respect, it is important to highlight the level that the modeling error can reach: the difference in the 1D and 3D responses (blue line in Figure 2b) can be as big as ~20%, and, very seldom, smaller than several percentage points. Therefore, the size of this mismatch should make accounting for the modeling error unavoidable; on the contrary, as mentioned before, the common practice is to tackle it by discretionally increase the measurement noise. By properly including the modeling error in the inversion, the assumed measurement error could be potentially reduced.

3.1.1. Deterministic Occam’s Inversion

Figure 3 shows the result of 1D deterministic inversion—implementing an Occam’s regularization strategy. Clearly, the deterministic result fits the observations extremely well (blue line in Figure 3b) and is capable to retrieve the major features of the true model whereas, in some cases, infers deceptive discontinuous reconstructions of the true interfaces (e.g., ~1600 to ~2000 m). It is worth mentioning that the data sensitivity to the model parameters drops below the first reconstructed conductive interface; hence, eventually, all the conductivity variation below that retrieved deep interface would not be considered reliable by skilled interpreters. In this respect, the deterministic result cannot be considered, at least from a practical point of view, much different from the probabilistic result obtained without accounting for the modeling error (and discussed in depth in the following Subsection). Nevertheless, it is true that the retrieved (laterally discontinuous) features (caused by fitting the coherent modeling error—as it will be clear later) might be challenging to be correctly deciphered and might lead to erroneous conclusions.

3.1.2. Stochastic Inversion without Modeling Error Assessment

The stochastic approach allows incorporation of (in principle) arbitrarily complex prior information, as long as realizations from the prior can be generated. In the present research, we generated independent realizations of p ( m ) representing buried valley structures (similar to those in Figure 1) by means of a Fast Fourier Transform Moving Average (FFT-MA) strategy [69] providing unconditional realizations of a Gaussian random field of the interfaces’ locations; in particular, the corresponding mean values are chosen to be uniformly distributed: (i) between 20 and 30 m for the shallowest interface, and (ii) between 65 and 85 m for the deepest, whereas the associated covariance models are characterized by a standard deviations and ranges, respectively, of: (i) 5 m and 100 m, for the first interface, and (ii) 80 m and 500 m for the second.
More details about possible ways of constructing the samples can be found in [65,66] including, for example, some details about the above-mentioned MPS strategies.
Whilst the geometry is varying, the conductivity values of each layer are kept constant realization-by-realization. For clarity, two examples of prior distribution realizations can be seen in Figure 4. It is evident that the parameters defining the realizations of the prior have been selected to be in agreement with our expectations about the geological structures we are dealing with.
It is worth highlighting that, since the inversion is one-dimensional, the actual (1D) models used to feed the inversion algorithm are the individual columns (and the corresponding 1D response) of each 3D realization (e.g., in Figure 4). Therefore, each single 1D prior model can be described by five parameters (the depths to the interfaces and the three conductivities).
Undoubtedly, the stochastic inversion, with such an informative prior, is quite facilitated and is basically reduced to the inference of the locations of the interfaces. As a matter of fact, as Figure 5 shows, by using the Metropolis algorithm [64,70], we can reconstruct quite satisfactorily the first interface at about 25–30 m depth, whereas we generally underestimate the depth of the second layer (except for the last ~1500 m on the right). However, what is most disturbing is that, even in this simple case, the misleading features reconstructed by the stochastic inversion (here, performed without accounting for the modeling error) appear to be almost certain.

3.1.3. Stochastic Inversion Incorporating the 1D Modeling Error

If instead, we use a subset of the 3D prior samples (Figure 4) to calculate their actual 3D responses and compare them with the data calculated, this time, by means of the 1D forward modeling applied to the 1D conductivity vertical profile in the center of each selected prior realization, we can estimate the appropriate mean d Δ and covariance C Δ , and incorporate them into the 1D stochastic inversion scheme. The results of the application of this new scheme to Test 1 is shown in Figure 6. Now, the result almost perfectly matches the true model. When the mean depth of the interfaces does not fit the true conductivity change—for example, near the steep lateral variations around 1000 m and 1400 m—the correct location of the boundaries still lies within the uncertainty bands.
Since we are dealing with 3D conductivity models, it is probably more appropriate to visualize the consequences of the different inversion schemes over several flight lines (Figure 7) across the Test 1 model. If we examine Figure 7b,c, the same conclusions drawn for an individual vertical section are clearly valid also for the other acquisition lines (and, consistently, also for the lateral intra-line resolution—kindly, compare Figure 7e,g).
A direct assessment of the performances of the different inversion schemes in terms of data fitting can be performed by looking at Figure 8b; in general, the relative mismatch between each of the channels of the observed 3D responses and those calculated (via a 3D forward) from the 3D conductivity distribution obtained via the stochastic inversion (with modeling error appraisal) lays within 10 % . It is worth noticing that the relatively poor data fitting on the right side of Figure 8b is caused by the higher level of heterogeneity of that end of the model, beyond the surveyed area. In fact, Figure 7d shows the rapidly varying morphology for X > 2500 m in correspondence of the section considered in Figure 8 ( Y ~ 500 m); that conductivity variation is not in accordance with the inevitable laterally homogeneous extension of the reconstructed solution (Figure 7g). This mismatch affects, not surprisingly, mainly the late-gate measurements. The same does not happen on the other end of the section, where the data fitting is particularly good (Figure 8b); in this case, the reason is that, differently from before, the true model continues largely unchanged towards low X values, at least for Y ~ 500 m (Figure 7d). On the other hand, the 3D responses generated by the conductivity volume retrieved by the 1D deterministic inversion (Figure 8a) demonstrates, once more, that fitting the data with a 1D forward (as in Figure 3b) does not necessarily guarantee that the corresponding 3D calculated data are in agreement with the (3D) observations; indeed, in Figure 8a, the relative misfit between the 3D calculated responses and the measurements ranges approximately between 30 % .

3.2. Test 2: 3D Conductivity Distribution with Heterogeneous Layers

In Test 2, we apply the same strategy to a more elaborated 3D conductivity model. Test 2’s model consists of three layers with similar geometries as in Test 1, but, now, characterized by heterogeneous conductivity values (Figure 9). For sake of completeness, analogously to what has been done for Test 1, also for Test 2, we show the 3D responses (Figure 10c) along the central profile (Figure 10a) of the 3D conductivity model (Figure 9a). In Figure 10b, we display the 1D response as they would result by considering each column of the true conductivity model as independent (for comparison, kindly, see Figure 2b concerning Test 1). Clearly, for Test 2, the importance of taking into account the modeling error is even more evident: the error we introduce when interpreting the 3D model response in terms of 1D data is never below 4% (blue line in Figure 10b).
Figure 11a presents the result of Occam’s inversion of the 3D data (Figure 10c). Again, from Figure 11b, it is clear that it is not difficult to perfectly fit the 3D data with 1D responses, even with a model barely capable to get the very major features of the conductivity model to be inferred.
Also for this more complex test, a 1D stochastic inversion can be performed by following the previously discussed Metropolis approach, making use of precalculated samples. The 105 realizations of the prior distribution are very similar to those shown in Figure 4 for Test 1; the only difference is that, consistently with the new model to be reconstructed, also for the prior realizations, the conductivity of each layer is allowed to vary. Figure 12 shows the mean map of the posterior distribution obtained without taking into consideration the modeling error; as for Figure 5 (about Test 1), also for Test 2, the retrieved mean model (Figure 12a) does not capture the lateral variation of the top of the deepest layer and smooths the large majority of the incisions out. Perhaps worse than that is the fact that, as shown by the standard deviations of the conductivity of each layer (plotted in Figure 12b), the leveled-out reconstruction is given for (almost) undeniable.
On the contrary, if we take into account the modeling error, the mean map derived from the 1D stochastic inversion provides a quite satisfactory reconstruction of the true model with all its complex morphology (Figure 13a), and when the inferred mean conductivity values do not reflect the true distribution, they are associated with high levels of uncertainty (for example at depth ~100 m and X ~1300 m in Figure 13b).

4. Discussion

Are 105 samples from the prior enough to guarantee the convergence of the stochastic inversions? What about the subset of 500 prior 3D realizations used for the assessment of the modeling error for the geological settings considered? Is the retrieved modeling error Gaussian as it should be to be able to use the proposed inversion scheme? In this section, we try to answer all these legitimate questions.

4.1. About the Numerosity of the Prior Samples for the Convergency of the Stochastic Inversion

Since the propositional scheme of the adopted Metropolis approach is based on a finite number of precalculated samples of the prior, it is important to establish if the abundance of those samples is sufficient to guarantee the convergence. In this respect, we run several inversions characterized by an increasing number of prior’s samples. From Figure 14f–h, it is evident that, at least for the simple problem of Test 1, convergence is reached with a numerosity of the prior samples of a few thousand. Hence, an abundance of two orders of magnitude higher should reasonably be compatible with our goal. A similar conclusion can be deduced by considering the evolution of the correlation coefficient between the depth of the deepest interface inferred and the true one, as shown in Figure 15.

4.2. About the Numerosity of the Prior’s Samples for the Estimation of the Modeling Error

Still considering Test 1’s dataset, we try to set up a strategy for the evaluation of the minimum number of realizations of the prior to be considered for an effective assessment of the modeling error. For all the tests performed in the present research, a maximum number of 500 models (similar to those in Figure 4) have been used to calculate the difference D d i f f and, in turn, the parameters defining the Gaussian probability density N ( d Δ ,   C Δ ) describing the modeling error. If we consider the behavior of the mean vector d Δ calculated for an increasing number of samples number N Δ , we can plot the results as in Figure 16a; from that figure, it appears evident that the assessment of d Δ reaches significant stability after considering a few hundred realizations of the prior.
A similar conclusion can be even more directly deduced by checking Figure 16b, in which the mean of the difference between d Δ and its best estimation d Δ ( 500 )   (based on N Δ = 500 realizations) is plotted against the increasing numerosity N Δ : approximately after considering 400 samples, d Δ does not show significant variations.
Regarding the evolution of the covariance matrix C Δ , we can, again, study how it changes over the number of prior samples. Figure 17 demonstrates, once more, that for the considered case, N Δ = 500 should guarantee a reasonable estimation of the modeling error. In particular, Figure 17b shows that the estimation of C Δ has already reached convergence after considering ~400 samples.

4.3. About the Gaussianity of the Modeling Error

As mentioned before, the proposed approach is based on the assumption that the modeling error is actually Gaussian. In this section, we demonstrate that this working hypothesis is largely met. In this respect, each subpanel of Figure 18 shows one of the 54 histograms (one for each time-gate) of the corresponding component of the vector D d i f f of the difference between the 3D and 1D responses calculated for the 500 realizations used for the modeling error assessment. In particular, the red lines represent the Gaussian that is better fitting the experimental histograms, whereas the black lines show the Gaussian profile as inferred from the distribution N ( d Δ ,   C Δ ) . It is worth noticing the excellent agreement between the black and red curves, but, more than that, the fact that for the large majority of the time-gates the modeling error is indeed compatible with Gaussian distributions. Surprisingly, the histograms for the late channels (in particular those for i going from ~49 to 54—i.e., basically the last row of Figure 18) show some sort of bimodal behavior whose possible justification is not evident. A possible reason might be connected with the finite size of the 3D simulation domain. However, this highly hypothetical guess will need to be investigated and verified.

5. Conclusions

Throughout this paper, two Tests are used to demonstrate the capabilities (and limitations) of stochastic inversion with geologically informed prior distribution to properly reconstruct the complex target conductivity distributions and the associated uncertainties.
In particular, we show how crucial the correct quantification of the modeling error (based on the prior choices) is. In our specific cases, we deal with ATEM data (mimicking the collection of VTEM measurements over geologies recalling glacial sedimentary environments typical, for example, of some regions in Northern Europe) inverted via a 1D stochastic approach with realistic prior and with the assessment of the corresponding modeling error. The conclusion is that, even in the case of stochastic approaches (already improving the results, for example, when compared with 1D deterministic approaches relying merely on the regularization term for including the prior knowledge about the targets), neglecting the fact that we are using a brute 1D approximation, instead of a more sophisticated 3D one, leads to unrealistically low levels of uncertainty in correspondence of those features that might turn out to be just artifacts. On the other hand, the assessment (and utilization) of the modeling error allows a more effective reconstruction of the true models and their associated reliability levels.
This research shows how the initial working hypotheses—concerning the minimum possible number of realizations defining the prior distribution and the numerosity of the smallest subset useful for an effective estimation of the modeling error—can be validated. Moreover, also the ansatz regarding the Gaussianity of the modeling error is largely verified.
Thus, the discussed workflow (presented and tested before on other kinds of data [47]) paves the way to the implementation of 1D probabilistic inversions of ATEM measurements capable to incorporate the complex pieces of geological information available and overcomes many of the difficulties connected with the utilization of efficient 1D approximations.
It is evident that the inclusion of the modeling error does not slow down the already available (and extremely fast) algorithms for 1D inversions (it does not matter if stochastic or deterministic). Hence, this approach will possibly remain useful also when fully 3D stochastic inversions will be practical; in fact, it will not be possible to consider any forward modeling tool perfect and, consequently, accounting for the modeling error will be, most likely, always beneficial.
It is also worth being highlighted that, at least for the investigated Tests, simply a few hundreds of 3D forward simulations are needed to retrieve a robust assessment of d Δ and C Δ , and, actually, the same estimation for d Δ and C Δ can be used, in principle, in any survey characterized by similar conditions (i.e. similar prior). Hence, in those cases, the efforts for the calculation of d Δ and C Δ would impact merely the first survey and clearly would not increase with the size of the survey. On the contrary, a full 3D inversion requires at least a few tens of iterations (i.e. 3D calculations) for each sounding location; this easily results in thousands of expensive 3D forward simulations. From these considerations, it clearly appears how convenient the proposed approach is. On the other hand, it is probably true that severely 3D targets will be, in any case, poorly reconstructed by using 1D approaches; however, the proper inclusion of the modeling error will be always useful to correctly estimate the high uncertainty of the 1D reconstruction of 3D inclusions, whereas not taking into account the modeling error will, most likely, lead to wrong solutions that look (incorrectly) certain.

Author Contributions

Conceptualization, P.B., G.V. and T.M.H.; methodology, P.B., G.V. and T.M.H.; software, P.B., G.V. and T.M.H.; validation, P.B., G.V. and T.M.H.; formal analysis, P.B., G.V. and T.M.H.; investigation, P.B., G.V. and T.M.H.; resources, G.V.; data curation, P.B.; writing—original draft preparation, P.B. and G.V.; writing—review and editing, P.B., G.V. and T.M.H.; visualization, P.B. and G.V.; supervision, G.V. and T.M.H.; project administration, G.V. and T.M.H.; funding acquisition, G.V. and T.M.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially funded: by the Independent Research Fund Denmark, grant number 7017-00160B; by the Italian Ministry of University and Research, through the initiative PON-RI 2014-2020, Asse I "Capitale Umano"–Azione I.1 "Dottorati innovativi con caratterizzazione industriale Ciclo XXXIII"–project: "GEOPROBARE: stochastic inversion of time-domain electromagnetic data" (CUP: F22J17000090007); by the University of Padova, through the research programme STARS@UNIPD–project: CHANGED – CHAracteriziNG pEatlands from Drones; by the Fondazione di Sardegna, through the initiative Progetti Biennali 2019–project: “GEO-CUBE” (CUP: F72F20000250007).

Acknowledgments

Thanks are due to Antonello Sanna, Roberto Ricciu and Giuseppe Desogus of the University of Cagliari (DICAAR) for providing most of the computer resources necessary to perform these tests.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhdanov, M.S.; Alfouzan, F.A.; Cox, L.; Alotaibi, A.; Alyousif, M.; Sunwall, D.; Endo, M. Large-scale 3D modeling and in-version of multiphysics airborne geophysical data: A case study from the Arabian Shield, Saudi Arabia. Minerals 2018, 8, 271. [Google Scholar] [CrossRef] [Green Version]
  2. Witherly, K. The quest for the Holy Grail in mining geophysics: A review of the development and application of airborne EM systems over the last 50 years. Geophysics 2000, 19, 270–274. [Google Scholar] [CrossRef]
  3. Alfouzan, F.A.; Alotaibi, A.M.; Cox, L.H.; Zhdanov, M.S. Spectral induced polarization survey with distributed array system for mineral exploration: Case study in Saudi Arabia. Minerals 2020, 10, 769. [Google Scholar] [CrossRef]
  4. Cudahy, T. Mineral mapping for exploration: An australian journey of evolving spectral sensing technologies and industry collaboration. Geosciences 2016, 6, 52. [Google Scholar] [CrossRef] [Green Version]
  5. Fountain, D. Airborne electromagnetic systems—50 years of development. Explor. Geophys. 1998, 29, 1–11. [Google Scholar] [CrossRef]
  6. Karshakov, E.V.; Podmogov, Y.G.; Kertsman, V.M.; Moilanen, J. Combined frequency domain and time domain airborne data for environmental and engineering challenges. J. Environ. Eng. Geophys. 2017, 22, 1. [Google Scholar] [CrossRef]
  7. Moilanen, E.; Karshakov, E.; Volkovitsky, A. Time domain helicopter EM system Equator: Resolution, sensitivity, universality. In 6th International AEM Conference & Exhibition; European Association of Geoscientists & Engineers: Houten, The Netherlands, 2013. [Google Scholar]
  8. Volkovitsky, A.; Karshakov, E. Airborne EM systems variety: What is the difference? In 6th International AEM Conference & Exhibition; European Association of Geoscientists & Engineers: Houten, The Netherlands, 2013. [Google Scholar]
  9. Combrinck, M.; Wright, R. Achieving accurate interpretation results from full-waveform streamed data AEM surveys. ASEG Ext. Abstr. 2016, 2016, 1–4. [Google Scholar] [CrossRef] [Green Version]
  10. Viezzoli, A.; Dauti, F.; Wijns, C. Robust scanning of AEM data for IP effects. Explor. Geophys. 2021, 52, 563–574. [Google Scholar] [CrossRef]
  11. Legault, J.M.; Izarra, C.; Prikhodko, A.; Zhao, S.; Saadawi, E.M. Helicopter EM (ZTEM–VTEM) survey results over the Nuqrah copper–lead–zinc–gold SEDEX massive sulphide deposit in the Western Arabian Shield, Kingdom of Saudi Arabia. Explor. Geophys. 2015, 46, 36–48. [Google Scholar] [CrossRef]
  12. Kwan, K.; Prikhodko, A.; Legault, J.M.; Plastow, G.C.; Kapetas, J.; Druecker, M. VTEM airborne EM, aeromagnetic and gamma-ray spectrometric data over the Cerro Quema high sulphidation epithermal gold deposits, Panama. Explor. Geophys. 2016, 47, 179–190. [Google Scholar] [CrossRef] [Green Version]
  13. Sørensen, K.I.; Auken, E. SkyTEM-A new high-resolution helicopter transient electromagnetic system. Explor. Geophys. 2004, 35, 191–199. [Google Scholar] [CrossRef]
  14. Boyko, W.; Paterson, N.R.; Kwan, K. AeroTEM characteristics and field results. Geophysics 2001, 20, 1130–1138. [Google Scholar] [CrossRef]
  15. Chen, T.; Hodges, G.; Miles, P. MULTIPULSE–high resolution and high power in one TDEM system. Explor. Geophys. 2015, 46, 49–57. [Google Scholar] [CrossRef]
  16. Smith, R.; Lee, T.J. The moments of the impulse response: A new paradigm for the interpretation of transient electromagnetic data. Geophysics 2002, 67, 1095–1103. [Google Scholar] [CrossRef]
  17. Annan, A.P.; Lookwood, R. An application of airborne geotem in australian conditions. Explor. Geophys. 1991, 22, 5–12. [Google Scholar] [CrossRef]
  18. Macnae, J. Improving the accuracy of shallow depth determinations in AEM sounding. Explor. Geophys. 2004, 35, 203–207. [Google Scholar] [CrossRef]
  19. Peters, G.; Street, G.; Kahimise, I.; Hutchins, D. Regional TEMPEST survey in north-east Namibia. Explor. Geophys. 2015, 46, 27–35. [Google Scholar] [CrossRef]
  20. Leggatt, P.B.; Klinkert, P.S.; Hage, T.B. The Spectrem airborne electromagnetic system—Further developments. Geophysics 2000, 65, 1976–1982. [Google Scholar] [CrossRef]
  21. Sapia, V.; Viezzoli, A.; Jørgensen, F.; Oldenborger, A.G.; Marchetti, M. The impact on geological and hydrogeological mapping results of moving from ground to airborne TEM. J. Environ. Eng. Geophys. 2014, 19, 53–66. [Google Scholar] [CrossRef] [Green Version]
  22. Siemon, B.; Christiansen, A.V.; Auken, E. A review of helicopter-borne electromagnetic methods for groundwater exploration. Near Surf. Geophys. 2009, 7, 629–646. [Google Scholar] [CrossRef] [Green Version]
  23. Siemon, B.; Seht, M.I.-V.; Frank, S. Airborne electromagnetic and radiometric peat thickness mapping of a bog in Northwest Germany (Ahlen-Falkenberger Moor). Remote Sens. 2020, 12, 203. [Google Scholar] [CrossRef] [Green Version]
  24. Høyer, A.-S.; Jørgensen, F.; Sandersen, P.B.; Viezzoli, A.; Møller, I. 3D geological modelling of a complex buried-valley network delineated from borehole and AEM data. J. Appl. Geophys. 2015, 122, 94–102. [Google Scholar] [CrossRef] [Green Version]
  25. Vignoli, G.; Gervasio, I.; Brancatelli, G.; Boaga, J.; Della Vedova, B.; Cassiani, G. Frequency-dependent multi-offset phase analysis of surface waves: An example of high-resolution characterization of a riparian aquifer. Geophys. Prospect. 2015, 64, 102–111. [Google Scholar] [CrossRef]
  26. Liu, R.; Guo, R.; Liu, J.; Liu, Z. An efficient footprint-guided compact finite element algorithm for 3-D airborne electro-magnetic modeling. IEEE Geosci. Remote. Sens. Lett. 2019, 16, 1809–1813. [Google Scholar] [CrossRef]
  27. Reninger, P.-A.; Martelet, G.; Perrin, J.; Dumont, M. Processing methodology for regional AEM surveys and local implications. Explor. Geophys. 2019, 51, 143–154. [Google Scholar] [CrossRef]
  28. Dzikunoo, E.A.; Vignoli, G.; Jørgensen, F.; Yidana, S.M.; Banoeng-Yakubo, B. New regional stratigraphic insights from a 3D geological model of the Nasia sub-basin, Ghana, developed for hydrogeological purposes and based on reprocessed B-field data originally collected for mineral exploration. Solid Earth 2020, 11, 349–361. [Google Scholar] [CrossRef] [Green Version]
  29. Christiansen, V.A.; Auken, E. Optimizing a layered and laterally constrained 2D inversion of resistivity data using Broyden’s update and 1D derivatives. J. Appl. Geophys. 2004, 56, 247–261. [Google Scholar] [CrossRef]
  30. Vignoli, G.; Fiandaca, G.; Christiansen, A.V.; Kirkegaard, C.; Auken, E. Sharp spatially constrained inversion with applications to transient electromagnetic data. Geophys. Prospect. 2014, 63, 243–255. [Google Scholar] [CrossRef]
  31. Brodie, R.; Sambridge, M. A holistic approach to inversion of frequency-domain airborne EM data. Geophysics 2006, 71, G301–G312. [Google Scholar] [CrossRef]
  32. Brodie, R.; Sambridge, M. A holistic approach to inversion of time-domain airborne EM. ASEG Ext. Abstr. 2006, 2006, 1–4. [Google Scholar] [CrossRef]
  33. Socco, L.V.; Boiero, D.; Foti, S.; Wisén, R. Laterally constrained inversion of ground roll from seismic reflection records. Geophysics 2009, 74, G35–G45. [Google Scholar] [CrossRef] [Green Version]
  34. Vignoli, G.; Deiana, R.; Cassiani, G. Focused inversion of vertical radar profile (VRP) traveltime data. Geophysics 2012, 77, H9–H18. [Google Scholar] [CrossRef]
  35. de Groot-Hedlin, C.; Constable, S. Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics 1990, 55, 1613–1624. [Google Scholar] [CrossRef]
  36. Vignoli, G.; Guillemoteau, J.; Barreto, J.; Rossi, M. Reconstruction, with tunable sparsity levels, of shear wave velocity profiles from surface wave data. Geophys. J. Int. 2021, 225, 1935–1951. [Google Scholar] [CrossRef]
  37. Zhdanov, M.; Vignoli, G.; Ueda, T. Sharp boundary inversion in crosswell travel-time tomography. J. Geophys. Eng. 2006, 3, 122–134. [Google Scholar] [CrossRef] [Green Version]
  38. Viezzoli, A.; Auken, E.; Munday, T. Spatially constrained inversion for quasi 3D modelling of airborne electromagnetic data–an application for environmental assessment in the Lower Murray Region of South Australia. Explor. Geophys. 2009, 40, 173–183. [Google Scholar] [CrossRef] [Green Version]
  39. Ley-Cooper, A.Y.; Viezzoli, A.; Guillemoteau, J.; Vignoli, G.; Macnae, J.; Cox, L.; Munday, T. Airborne electromagnetic mod-elling options and their consequences in target definition. Explor. Geophys. 2015, 46, 74–84. [Google Scholar] [CrossRef]
  40. Christiansen, A.V.; Auken, E.; Kirkegaard, C.; Schamper, C.; Vignoli, G. An efficient hybrid scheme for fast and accurate inversion of airborne transient electromagnetic data. Explor. Geophys. 2016, 47, 323–330. [Google Scholar] [CrossRef] [Green Version]
  41. Bai, P.; Vignoli, G.; Viezzoli, A.; Nevalainen, J.; Vacca, G. (Quasi-) real-time inversion of airborne time-domain electro-magnetic data via artificial neural network. Remote Sens. 2020, 12, 3440. [Google Scholar] [CrossRef]
  42. Hansen, T.M. Probabilistic inverse problems using machine learning-applied to inversion of airborne EM data. Geophys. Res. Lett. 2021. submitted. [Google Scholar] [CrossRef]
  43. Auken, E.; Christiansen, A.V.; Kirkegaard, C.; Fiandaca, G.; Schamper, C.; Behroozmand, A.A.; Binley, A.; Nielsen, E.; Effersø, F.; Christensen, N.B.; et al. An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data. Explor. Geophys. 2015, 46, 223–235. [Google Scholar] [CrossRef] [Green Version]
  44. Vest Christiansen, A.V.; Auken, E.; Viezzoli, A. Quantification of modeling errors in airborne TEM caused by inac-curate system description. Geophysics 2011, 76, F43–F52. [Google Scholar] [CrossRef] [Green Version]
  45. Oldenburg, D.W.; Heagy, L.J.; Kang, S.; Cockett, R. 3D electromagnetic modelling and inversion: A case for open source. Explor. Geophys. 2019, 51, 25–37. [Google Scholar] [CrossRef] [Green Version]
  46. Auken, E.; Foged, N.; Larsen, J.J.; Lassen, K.V.T.; Maurya, P.K.; Dath, S.M.; Eiskjær, T.T. tTEM—A towed transient electro-magnetic system for detailed 3D imaging of the top 70 m of, the subsurface. Geophysics 2019, 84, E13–E22. [Google Scholar] [CrossRef]
  47. Hansen, T.M.; Cordua, K.S.; Jacobsen, B.H.; Mosegaard, K. Accounting for imperfect forward modeling in geophysical inverse problems—Exemplified for crosshole tomography. Geophysics 2014, 79, H1–H21. [Google Scholar] [CrossRef] [Green Version]
  48. Cuma, M.; Cox, L.; Zhdanov, M. Paradigm change in interpretation of AEM data by using a large-scale parallel 3D inversion. In Near Surface Geoscience—First Conference on Geophysics for Mineral Exploration and Mining; European Association of Geoscientists & Engineers: Houten, The Netherlands, 2016. [Google Scholar]
  49. Yang, D.; Oldenburg, D.W. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit. Geophysics 2012, 77, B23–B34. [Google Scholar] [CrossRef]
  50. Liu, Y.; Yin, C. 3D inversion for multipulse airborne transient electromagnetic data. Geophysics 2016, 81, E401–E408. [Google Scholar] [CrossRef]
  51. Haber, E.; Ascher, U.M.; Oldenburg, D.W. Inversion of 3D electromagnetic data in frequency and time domain using an inexact all-at-once approach. Geophysics 2004, 69, 1216–1228. [Google Scholar] [CrossRef]
  52. McMillan, M.S.; Schwarzbach, C.; Haber, E.; Oldenburg, D.W. 3D parametric hybrid inversion of time-domain airborne electromagnetic data. Geophysics 2015, 80, K25–K36. [Google Scholar] [CrossRef] [Green Version]
  53. Cox, L.H.; Wilson, G.A.; Zhdanov, M. 3D inversion of airborne electromagnetic data using a moving footprint. Explor. Geophys. 2010, 41, 250–259. [Google Scholar] [CrossRef]
  54. Zhdanov, M. Foundations of Geophysical Electromagnetic Theory and Methods; Elsevier: Amsterdam, The Netherlands, 2018. [Google Scholar] [CrossRef]
  55. Viezzoli, A.; Munday, T.; Auken, E.; Christiansen, A.V. Accurate quasi 3D versus practical full 3D inversion of AEM data–the Bookpurnong case study. Preview 2010, 149, 23–31. [Google Scholar] [CrossRef] [Green Version]
  56. Hauser, J.; Gunning, J.; Annetts, D. Probabilistic inversion of airborne electromagnetic data under spatial constraints. Geophysics 2015, 80, E135–E146. [Google Scholar] [CrossRef]
  57. Tarantola, A. Inverse Problem Theory and Methods for Model Parameter Estimation; SIAM: Philadelphia, PA, USA, 2005. [Google Scholar] [CrossRef] [Green Version]
  58. Zhdanov, M.S. Geophysical Inverse Theory and Regularization Problems; Elsevier: Amsterdam, The Netherlands, 2002. [Google Scholar]
  59. Tarantola, A.; Valette, B. Generalized nonlinear inverse problems solved using the least squares criterion. Rev. Geophys. 1982, 20, 219–232. [Google Scholar] [CrossRef]
  60. Constable, S.; Parker, R.L.; Constable, C. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 1987, 52, 289–300. [Google Scholar] [CrossRef]
  61. Høyer, A.S.; Vignoli, G.; Hansen, T.M.; Vu, L.T.; Keefer, D.A.; Jørgensen, F. Multiple-point statistical simulation for hydrogeological models: 3-D training image development and conditioning strategies. Hydrol. Earth Syst. Sci. 2017, 21, 6069–6089. [Google Scholar] [CrossRef] [Green Version]
  62. Haber, E. Computational Methods in Geophysical Electromagnetics; SIAM: Philadelphia, PA, USA, 2014. [Google Scholar] [CrossRef]
  63. Vallée, M.A.; Smith, R.S. Application of Occam’s inversion to airborne time-domain electromagnetics. Lead. Edge 2009, 28, 284–287. [Google Scholar] [CrossRef]
  64. Mosegaard, K.; Tarantola, A. Monte Carlo sampling of solutions to inverse problems. J. Geophys. Res. Space Phys. 1995, 100, 12431–12447. [Google Scholar] [CrossRef]
  65. Hansen, T.M.; Cordua, K.S.; Looms, M.C.; Mosegaard, K. SIPPI: A Matlab toolbox for sampling the solution to inverse problems with complex prior information: Part 1—Methodology. Comput. Geosci. 2013, 52, 470–480. [Google Scholar] [CrossRef] [Green Version]
  66. Hansen, T.M.; Cordua, K.S.; Looms, M.C.; Mosegaard, K. SIPPI: A Matlab toolbox for sampling the solution to inverse problems with complex prior information: Part 2—Application to crosshole GPR tomography. Comput. Geosci. 2013, 52, 481–492. [Google Scholar] [CrossRef] [Green Version]
  67. Jørgensen, F.; Sandersen, P.B. Buried and open tunnel valleys in Denmark—Erosion beneath multiple ice sheets. Quat. Sci. Rev. 2006, 25, 1339–1363. [Google Scholar] [CrossRef]
  68. Kehew, A.E.; Piotrowski, J.A.; Jørgensen, F. Tunnel valleys: Concepts and controversies—A review. Earth-Sci. Rev. 2012, 113, 33–58. [Google Scholar] [CrossRef]
  69. Le Ravalec, M.; Noetinger, B.; Hu, L.Y. The FFT Moving Average (FFT-MA) Generator: An efficient numerical method for generating and conditioning gaussian simulations. Math. Geol. 2000, 32, 701–723. [Google Scholar] [CrossRef]
  70. Hansen, T.M.; Minsley, B.J. Inversion of airborne EM data with an explicit choice of prior model. Geophys. J. Int. 2019, 218, 1348–1366. [Google Scholar] [CrossRef]
Figure 1. Conductivity distribution for Test 1: (a) 3D view of the model consisting of a sequence of three homogeneous layers with varying thicknesses; (b) Vertical section of the model in panel (a), along the survey line highlighted by the red dots indicating the locations of the ATEM soundings; (c) Plain view at 100 m depth (dash green line in panel (b)).
Figure 1. Conductivity distribution for Test 1: (a) 3D view of the model consisting of a sequence of three homogeneous layers with varying thicknesses; (b) Vertical section of the model in panel (a), along the survey line highlighted by the red dots indicating the locations of the ATEM soundings; (c) Plain view at 100 m depth (dash green line in panel (b)).
Remotesensing 13 03881 g001
Figure 2. Comparison of the 3D and 1D responses for the conductivity model of Test 1: (a) Vertical section of the 3D conductivity model in Figure 1a (it is a portion of the section in Figure 1b); (b) 1D responses calculated for the 1D portions of the original model in Figure 1a taken at the locations of the ATEM soundings (red dots in Figure 1b) – the blue line represents, sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right, in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 1a).
Figure 2. Comparison of the 3D and 1D responses for the conductivity model of Test 1: (a) Vertical section of the 3D conductivity model in Figure 1a (it is a portion of the section in Figure 1b); (b) 1D responses calculated for the 1D portions of the original model in Figure 1a taken at the locations of the ATEM soundings (red dots in Figure 1b) – the blue line represents, sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right, in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 1a).
Remotesensing 13 03881 g002
Figure 3. 1D deterministic inversion of the 3D data (Figure 2c) associated to the conductivity model of Test 1 (Figure 1a): (a) The solution of the 1D deterministic inversion – the black lines show the interfaces of the original conductivity model to be reconstructed; (b) The 1D responses resulting from the conductivity model in panel (a)–the blue line represents sounding-by-sounding, the relative misfit between the 3D and 1D responses (the corresponding axis is on the right in blue).
Figure 3. 1D deterministic inversion of the 3D data (Figure 2c) associated to the conductivity model of Test 1 (Figure 1a): (a) The solution of the 1D deterministic inversion – the black lines show the interfaces of the original conductivity model to be reconstructed; (b) The 1D responses resulting from the conductivity model in panel (a)–the blue line represents sounding-by-sounding, the relative misfit between the 3D and 1D responses (the corresponding axis is on the right in blue).
Remotesensing 13 03881 g003
Figure 4. (a,b) show two examples of possible, distinct, realizations of the prior distribution used for the stochastic inversion of Test 1’s dataset.
Figure 4. (a,b) show two examples of possible, distinct, realizations of the prior distribution used for the stochastic inversion of Test 1’s dataset.
Remotesensing 13 03881 g004
Figure 5. 1D stochastic inversion of the 3D data (Figure 2c) associated to the conductivity model of Test 1 (Figure 2a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used (Figure 4), but no modeling error has been taken into account. The orange bands represent the reconstruction uncertainty defined by the mean and standard deviation values (red points and vertical bars) deduced by the retrieved realizations of the posterior p ( m | d ) . As for Figure 3a, the black lines show the locations of the interfaces to be reconstructed.
Figure 5. 1D stochastic inversion of the 3D data (Figure 2c) associated to the conductivity model of Test 1 (Figure 2a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used (Figure 4), but no modeling error has been taken into account. The orange bands represent the reconstruction uncertainty defined by the mean and standard deviation values (red points and vertical bars) deduced by the retrieved realizations of the posterior p ( m | d ) . As for Figure 3a, the black lines show the locations of the interfaces to be reconstructed.
Remotesensing 13 03881 g005
Figure 6. 1D stochastic inversion of the 3D data (Figure 2c) associated with the conductivity model of Test 1 (Figure 1a). In this case, realistic prior is used (Figure 4) together with the modeling error assessment. As in Figure 5, the orange bands represent the uncertainty (the associated mean and standard deviations are plotted as vertical red bars) calculated by the retrieved realizations of the posterior p ( m | d ) . The black lines show the locations of the interfaces to be reconstructed (Figure 2a).
Figure 6. 1D stochastic inversion of the 3D data (Figure 2c) associated with the conductivity model of Test 1 (Figure 1a). In this case, realistic prior is used (Figure 4) together with the modeling error assessment. As in Figure 5, the orange bands represent the uncertainty (the associated mean and standard deviations are plotted as vertical red bars) calculated by the retrieved realizations of the posterior p ( m | d ) . The black lines show the locations of the interfaces to be reconstructed (Figure 2a).
Remotesensing 13 03881 g006
Figure 7. The inversion results of the 3D data from (a) Test 1’s conductivity model (Figure 1) and obtained with: (b) Occam’s deterministic strategy, and (c) The stochastic strategy incorporating the modeling error. The locations of the inverted soundings are shown as red dots in (d) on top of the plain view of the topography of the second interface of the true model–shown in (e). In (f), the 3D view of the results from (b) is plotted, whereas, in (g), it is possible to see a similar 3D view but now based on the stochastic results in (c). For clarity, in panels (eg), only the sounding locations of the central flight line are shown (as red dots).
Figure 7. The inversion results of the 3D data from (a) Test 1’s conductivity model (Figure 1) and obtained with: (b) Occam’s deterministic strategy, and (c) The stochastic strategy incorporating the modeling error. The locations of the inverted soundings are shown as red dots in (d) on top of the plain view of the topography of the second interface of the true model–shown in (e). In (f), the 3D view of the results from (b) is plotted, whereas, in (g), it is possible to see a similar 3D view but now based on the stochastic results in (c). For clarity, in panels (eg), only the sounding locations of the central flight line are shown (as red dots).
Remotesensing 13 03881 g007
Figure 8. Comparison of the percentage misfit between the observed 3D responses and those calculated from: (a) The 3D conductivity model retrieved by the 1D deterministic approach (Figure 7b,f); (b) The 3D conductivity model obtained via the stochastic inversion accounting for the modeling error (Figure 7c,g). Each of the 54 plotted time-gates (channels) is depicted with a different color.
Figure 8. Comparison of the percentage misfit between the observed 3D responses and those calculated from: (a) The 3D conductivity model retrieved by the 1D deterministic approach (Figure 7b,f); (b) The 3D conductivity model obtained via the stochastic inversion accounting for the modeling error (Figure 7c,g). Each of the 54 plotted time-gates (channels) is depicted with a different color.
Remotesensing 13 03881 g008
Figure 9. Conductivity distribution for Test 2: (a) 3D view of the model consisting of a sequence of three heterogeneous layers with varying thicknesses; (b) Vertical section of the model in panel (a), along the survey line highlighted by the red dots indicating the locations of the ATEM soundings; (c) Plain view at 100 m depth (dash green line in panel (b)).
Figure 9. Conductivity distribution for Test 2: (a) 3D view of the model consisting of a sequence of three heterogeneous layers with varying thicknesses; (b) Vertical section of the model in panel (a), along the survey line highlighted by the red dots indicating the locations of the ATEM soundings; (c) Plain view at 100 m depth (dash green line in panel (b)).
Remotesensing 13 03881 g009
Figure 10. Comparison of the 3D and 1D responses for the conductivity model of Test 2: (a) Vertical section of the 3D conductivity model in Figure 9a (it is a portion of the section in Figure 9b); (b) 1D responses calculated for the 1D portions of the original model in Figure 9a at the locations of the ATEM soundings (red dots in Figure 9b)–the blue line represents sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 9a).
Figure 10. Comparison of the 3D and 1D responses for the conductivity model of Test 2: (a) Vertical section of the 3D conductivity model in Figure 9a (it is a portion of the section in Figure 9b); (b) 1D responses calculated for the 1D portions of the original model in Figure 9a at the locations of the ATEM soundings (red dots in Figure 9b)–the blue line represents sounding-by-sounding, the relative misfit between the 1D and 3D responses (the corresponding axis is on the right in blue); (c) 3D responses measured at the same location in panel (b), but, here, calculated for the entire 3D conductivity model (Figure 9a).
Remotesensing 13 03881 g010
Figure 11. 1D deterministic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 9a): (a) The solution of the 1D deterministic inversion–the black lines show the interfaces of the original conductivity model to be reconstructed; (b) The 1D responses resulting from the conductivity model in panel (a)—the blue line represents, sounding-by-sounding, the relative misfit between the 3D and 1D responses (the corresponding axis is on the right in blue).
Figure 11. 1D deterministic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 9a): (a) The solution of the 1D deterministic inversion–the black lines show the interfaces of the original conductivity model to be reconstructed; (b) The 1D responses resulting from the conductivity model in panel (a)—the blue line represents, sounding-by-sounding, the relative misfit between the 3D and 1D responses (the corresponding axis is on the right in blue).
Remotesensing 13 03881 g011
Figure 12. (a) Mean map of the 1D stochastic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used, but no modeling error has been taken into account. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with the mean and standard deviation vertical profile superimposed for each vertical conductivity profile.
Figure 12. (a) Mean map of the 1D stochastic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used, but no modeling error has been taken into account. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with the mean and standard deviation vertical profile superimposed for each vertical conductivity profile.
Remotesensing 13 03881 g012
Figure 13. (a) Mean map of the 1D stochastic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used together with the modeling error assessment. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with, superimposed, the mean and standard deviation vertical profiles.
Figure 13. (a) Mean map of the 1D stochastic inversion of the 3D data (Figure 10c) associated to the conductivity model of Test 2 (Figure 10a shows one vertical section of that 3D conductivity model). In this case, realistic prior is used together with the modeling error assessment. The black lines show the locations of the interfaces to be reconstructed. (b) The same reconstructed conductivity distribution as in panel (a), but now with, superimposed, the mean and standard deviation vertical profiles.
Remotesensing 13 03881 g013
Figure 14. Result of the extended Metropolis inversion as a function of the numerosity of the considered precalculated prior samples–the title of each panel (ah) reports that numerosity. Clearly, here, we are incorporating the information about the modeling error into the inversion of Test 1’s data.
Figure 14. Result of the extended Metropolis inversion as a function of the numerosity of the considered precalculated prior samples–the title of each panel (ah) reports that numerosity. Clearly, here, we are incorporating the information about the modeling error into the inversion of Test 1’s data.
Remotesensing 13 03881 g014
Figure 15. Correlation coefficient between the retrieved and the true depth of the deepest layer as a function of the number of precalculated samples of the prior. Clearly, as for Figure 14, we are considering the stochastic inversion with modeling error assessment applied to Test 1’s dataset.
Figure 15. Correlation coefficient between the retrieved and the true depth of the deepest layer as a function of the number of precalculated samples of the prior. Clearly, as for Figure 14, we are considering the stochastic inversion with modeling error assessment applied to Test 1’s dataset.
Remotesensing 13 03881 g015
Figure 16. Mean of the modeling error as the number of the prior samples used for its estimation increases: (a) the mean d Δ vectors calculated for a number of prior samples ranging, for example, from 1 to 100 are represented by the solid yellow lines, whereas, considering another example, for a number of prior samples between 401 and 500, the same results for the d Δ vectors are plotted in dark red. (b) shows the behavior of the mean misfit between the estimated d Δ and our best assessment d Δ ( 500 ) (based on 500 realizations).
Figure 16. Mean of the modeling error as the number of the prior samples used for its estimation increases: (a) the mean d Δ vectors calculated for a number of prior samples ranging, for example, from 1 to 100 are represented by the solid yellow lines, whereas, considering another example, for a number of prior samples between 401 and 500, the same results for the d Δ vectors are plotted in dark red. (b) shows the behavior of the mean misfit between the estimated d Δ and our best assessment d Δ ( 500 ) (based on 500 realizations).
Remotesensing 13 03881 g016
Figure 17. Covariance C Δ of the modeling error as the number of prior samples used for its estimation increases: (a) different C Δ ’s calculated for the number of prior samples indicated by the title of each subpanel (from N Δ = 10 for the panel on the top-left to N Δ = 500 for our best estimation on the bottom-right corner of panel (a)). (b) shows the behavior of the maximum misfit between all the components of the matrix C Δ and the corresponding C Δ ( 500 ) based on our maximum number of realizations ( N Δ = 500).
Figure 17. Covariance C Δ of the modeling error as the number of prior samples used for its estimation increases: (a) different C Δ ’s calculated for the number of prior samples indicated by the title of each subpanel (from N Δ = 10 for the panel on the top-left to N Δ = 500 for our best estimation on the bottom-right corner of panel (a)). (b) shows the behavior of the maximum misfit between all the components of the matrix C Δ and the corresponding C Δ ( 500 ) based on our maximum number of realizations ( N Δ = 500).
Remotesensing 13 03881 g017
Figure 18. Each panel j (with j , time-gate’s index, varying from 1 to 54) shows the histogram of the corresponding component   [ D d i f f ] j of the difference between the 3D and 1D responses for the 500 samples of the prior used for the modeling error estimation. The solid black line is the Gaussian curve better fitting the histogram, whereas the red line is the Gaussian derived by N ( d Δ ,   C Δ ) .
Figure 18. Each panel j (with j , time-gate’s index, varying from 1 to 54) shows the histogram of the corresponding component   [ D d i f f ] j of the difference between the 3D and 1D responses for the 500 samples of the prior used for the modeling error estimation. The solid black line is the Gaussian curve better fitting the histogram, whereas the red line is the Gaussian derived by N ( d Δ ,   C Δ ) .
Remotesensing 13 03881 g018
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bai, P.; Vignoli, G.; Hansen, T.M. 1D Stochastic Inversion of Airborne Time-Domain Electromagnetic Data with Realistic Prior and Accounting for the Forward Modeling Error. Remote Sens. 2021, 13, 3881. https://doi.org/10.3390/rs13193881

AMA Style

Bai P, Vignoli G, Hansen TM. 1D Stochastic Inversion of Airborne Time-Domain Electromagnetic Data with Realistic Prior and Accounting for the Forward Modeling Error. Remote Sensing. 2021; 13(19):3881. https://doi.org/10.3390/rs13193881

Chicago/Turabian Style

Bai, Peng, Giulio Vignoli, and Thomas Mejer Hansen. 2021. "1D Stochastic Inversion of Airborne Time-Domain Electromagnetic Data with Realistic Prior and Accounting for the Forward Modeling Error" Remote Sensing 13, no. 19: 3881. https://doi.org/10.3390/rs13193881

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop