Inversion of electromagnetic induction data using a novel wavelet-based and scale-dependent regularization term
Wouter Deleersnyder, Benjamin Maveau, Thomas Hermans, David Dudal
ssubmitted to
Geophys. J. Int.
Inversion of electromagnetic induction data using a novelwavelet-based and scale-dependent regularization term
Wouter Deleersnyder , , Benjamin Maveau , Thomas Hermans , David Dudal , KU Leuven Campus Kortrijk - KULAK, Department of Physics, Etienne Sabbelaan 53, 8500 Kortrijk, Belgium Ghent University, Department of Geology, Krijgslaan 281 - S8, 9000 Gent, Belgium Ghent University, Department of Physics and Astronomy, Ghent, Krijgslaan 281 - S9, 9000 Gent, Belgium
SUMMARY
The inversion of electromagnetic induction data to a conductivity profile is an ill-posedproblem. Regularization improves the stability of the inversion and, based on Occam’srazor principle, a smoothing constraint is typically used. However, the conductivity pro-files are not always expected to be smooth. Here, we develop a new inversion schemein which we transform the model to the wavelet space and impose a sparsity constraint.This sparsity constrained inversion scheme will minimize an objective function with aleast-squares data misfit and a sparsity measure of the model in the wavelet domain.A model in the wavelet domain has both temporal as spatial resolution, and penalizingsmall-scale coefficients effectively reduces the complexity of the model. Depending onthe expected conductivity profile, an optimal wavelet basis function can be chosen. Thescheme is thus adaptive. Finally, we apply this new scheme on a frequency domain elec-tromagnetic sounding (FDEM) dataset, but the scheme could equally apply to any other1D geophysical method.
Key words:
Inverse theory – Wavelet transform – Electromagnetic theory – Controlledsource electromagnetics (CSEM) a r X i v : . [ phy s i c s . g e o - ph ] S e p Wouter Deleersnyder
Electromagnetic induction (EMI) surveys aim to image the electrical properties of the subsurfacenon-invasively. Via Archie’s empirical petrophysical law (Archie 1942), the electrical properties arerelated to the soil characteristics. Archie’s law relates the water conductivity ( σ w ), the porosity ( θ ), thesaturation of the pores ( S ) and experimental coefficients m and n (they are geometric factors) to theconductivity of the soil ( σ ): σ = σ w θ m S n . EMI methods are used in archaeological prospection (Saey et al. 2012), soil contamination (Pettersson& Nobes 2003), detection of electric anomalies (such as unexploded ordnance (Fern´andez et al. 2010)),agriculture (Jadoon et al. 2015) and saltwater intrusion (Paepen et al. 2020).EMI instrumentation can often be used with different intercoil spacings or heights, resulting in datawith a sensitivity varying with depth, allowing to invert the data and obtain a conductivity profile of thesubsurface. Geophysical inversion is an ill-posed problem, because the solution is usually non-unique,which is generally cured with regularization techniques. Additionally, a slightly different dataset willgenerally yield an entirely different result. Deterministic methods impose one or two constraints, suchas the smoothness constraint imposed by the traditional Tikhonov regularization (Tikhonov 1943).The smoothness constraint is not always adapted to the subsurface structure (Linde et al. 2015). Whilemany alternatives exist, it is not always easy to find the best adapted scheme. In earlier work on EMIinversion, Tikhonov regularization, truncated singular value decomposition (Deidda et al. 2014) orsparsity based inversion (Tantum et al. 2012) have beencommonly applied. In our work, a new regu-larization scheme is developed that is more adaptive than many other schemes as it can both representsmooth and sharp models. The regularization term is a sparsity constraint of the model parameters inwavelet domain. The wavelet basis determines the type of constraint (blocky, smooth or in between)that is imposed and hence this is an adaptive inversion procedure.The literature about full-waveform inversion is richer and faces similar problems. Full-waveforminversion is concerned with fitting the model parameters, velocities and impedances, to the recordedseismic or ground penetrating radar data. The ill-posedness of the problem similarly requires a suit-able regularization scheme in order to obtain a realistic model. The non-uniqueness is tackled typicallywith Tikhonov regularization (Virieux & Operto 2009). More recently, total variation regularization isapplied to this field, including blocky regularization schemes with several focusing functions (Guitton2012), such as the (cid:96) -norm.The high non-linearity of full-waveform inversion results typically in a multimodal objective function, nversion with wavelet-based and scale-dependent regularization (cid:96) -norm in the space domain. In theirmethodology, a finite difference method was used to solve the forward problem. Due to the slow for-ward modelling, only about ten iterations could be considered in the inverse problem. In this study,we develop a different inversion scheme, that uses a scale-dependent penalizing term (a comparison isprovided in Section 3.1). Inspired by the technicalities of the wavelet transform (cfr. Theorem 2.1 inSection 2.3) and the pragmatic idea that small scale details should be penalized more heavily than thecoarser main structure of the conductivity profile, each length scale of the conductivity profile is penal-ized with a different weight. Furthermore, our scheme uses a faster (approximate) forward modellingoperator (see Section 2.1). This allows us to compute much more iterations (in fact, computationalpower is not an issue anymore), at the cost of an additional modelling error due to the approximationsin the forward operator. In Section 3.3, we demonstrate that this modelling error poses not much ofan issue. Our results show mainly two advantages. First, the scheme is more adaptive as it can recoverboth blocky and smooth conductivity profiles. Secondly, it can recover high amplitude anomalies in aglobally smooth conductivity profile.In this work, we first present the methods and the technicalities of the inversion scheme, includingforward modelling, the inversion procedure and the basics of wavelet theory. In particular, we reviewthe properties of wavelet basis and discuss why they are good candidates as a family of regularizationfunctionals for geophysical inversion. In Sections 3.1-3.3, we apply our regularization scheme onsynthetic data and compare it with the commonly employed Tikhonov regularization as a benchmark.Finally, we apply our scheme to a real dataset of the Gontrode forest, Belgium, a context at which an1D inversion scheme is realistic (Bobe et al. 2020). Wouter Deleersnyder
The forward model describes the soil’s response d ∈ C n d × to the magnetic dipole placed at height h above the surface and a soil with specific parameter distribution m ∈ R n m × , measured at n d intercoilspacings s . The parameter distribution m is referred to as model or model parameters and containsthe electrical conductivities in Siemens per meter for a specific (usually equidistant) parametrizationof the subsurface. Finding the operator F , such that d = F ( m ) (1)holds, is called the forward problem . Wait (1951) derived such a model, taking into account the loop-loop couplings between eddy currents and the dampening of the electromagnetic fields, for the 1Dapproximation (horizontally stratified earth) in the quasi-stationary Maxwell regime. The model in-volves recursive relations with integration over Bessel functions, which have to be executed numeri-cally. This computationally demanding task slows down an inversion procedure in which the forwardproblem needs to be calculated at every iteration. Wait later derived a linear model (Wait 1962), whichis extensively used in the work of McNeill (McNeill 1980) and many others (e.g. (Simpson et al. 2009)or (De Smedt et al. 2011) ), known as the LIN approximation. The LIN approximation is linear andneglects the couplings and electromagnetic dampening. It is only valid at Low Induction Numbers (cid:63) . Amore recently proposed model by (Maveau et al. 2020), based on the LIN approach, takes into accountthe electromagnetic dampening and provides accurate results under more lenient conditions, especiallyin more conductive (e.g. high salinity) backgrounds. Using an approximate forward operator F approx introduces a biased modelling error κ , we write F ( m ) = F approx ( m ) + κ . (2)Another method to solve the forward problem is to use finite elements methods (FEM) or finitevolume methods (FVM), such as provided for example in SimPEG (Cockett et al. 2015). This type ofmodelling is considered exact if the mesh is sufficiently fine. For cylindrically symmetric set-ups (i.e.horizontally stratified conductivity profile and a vertical dipole), the computation is relatively fast. In (cid:63) The= induction number B is defined as the intercoil distance s divided by the skin depth δ , i.e. B = ωµ σs , where ω is the angularfrequency of the dipole, µ is the magnetic permeability of the vacuum and σ the electric conductivity of the medium. The meaningof sufficiently low induction numbers is that if the skin depth is much larger than the path the electromagnetic field has to traverse, thedampening can be neglected. Indeed, the path length is restricted by the fall-off of the magnetic dipole, which has the same order as theintercoil distance s . nversion with wavelet-based and scale-dependent regularization κ and the forward operator F . Due to the non-linearity of many forward operators, the inverse problem is typically solved as an opti-mization problem. A suitable parameter distribution of the subsurface is determined by a minimum ofan objective function φ . The objective function takes the parameters distribution or model parametersas input and expresses how well the parameters fit the data. In minimum structure inversion proce-dures, an additional measure of model complexity, a regularization or model misfit term, is added tothe objective function. The minimum is usually obtained via gradient descent-like methods.The objective function in minimum structure inversion is a combination of a data misfit functional φ d , a model misfit functional φ m and a regularization parameter λ . The latter balances the relative im-portance of both misfit terms. The data misfit functional φ d measures how well the model parameters m fit the data d and is defined as a least-squares term φ d = 12 || D d ( d − F ( m )) || , (3)where D d is a diagonal matrix whose elements are usually the reciprocals of the estimated noisestandard deviation (Strutz 2010). This weighting matrix penalizes the data misfit for data points withminimal measurement error less than data points with larger measurement error. The model misfitfunctional is a regularization term that is introduced to handle the ill-posedness of the problem. In-deed, geophysical inversion problems are ill-conditioned and the solution is typically not unique. Atraditional method for minimum structure inversion is to apply a smoothing constraint, such as withTikhonov regularization (Tikhonov 1943). Tikhonov regularization improves the stability of the inver-sion, but promotes smooth solutions and hence cannot recover blocky structures. Other regularizationterms exist, such as the (cid:96) -norm inversion in original model space (Farquharson 2007), the minimum Wouter Deleersnyder gradient support functional (Portniaguine & Zhdanov 1999) or a sequential inversion (Guillemoteauet al. 2016).In this paper, we develop a different approach. Suppose that there exists a basis in which the truemodel parameters m , known to possess minimum structure, are represented in a sparse fashion in x ∈ R n x × . Then, x = W m , (4)where W is the basis transformation. Such a basis transformation, in combination with a sparsitypromoting measure, yields an appropriate model misfit, because a complex structure in the modelparameters m will have many non-zero entries in the x -representation and will therefore be heavilypenalized by the sparsity promoting measure. The (cid:96) -norm is a well-known sparsity promoting mea-sure. We will rely on gradient decent-like methods and since the (cid:96) -norm is not differentiable at zero,we will use the perturbed (cid:96) -norm measure of Ekblom (Ekblom 1987) µ Ekblom ( x ) = (cid:112) x + (cid:15), (5)which is very similar to the (cid:96) -norm for small (cid:15) and which is also convex. The latter is a welcomeproperty for the optimization problem. The model misfit in terms of x in our inversion scheme is thus φ m ( x ) = (cid:88) j (cid:113) x j + (cid:15). (6)In summary, we have min x φ = min x ( φ d + λφ m ) (7) = min x (cid:18) || D d ( d − F approx (cid:0) W − x (cid:1) || + λ n x (cid:88) j =1 (cid:113) x j + (cid:15) , (8)where the regularization parameter λ is still undetermined. A regularization parameter can be toosmall, that is the case when the data misfit φ d is smaller than the actual noise on the data, i.e. whenwe have overfitting. On the other hand, with a too large regularization parameter, the parameter dis-tribution will exhibit too little structure. Hansen (2010) discusses the common automatic methods toestimate an optimal regularization parameter. We will rely on the L-curve criterion. This method re-quires us to solve the inverse problem for multiple regularization parameters and plots the model misfitin terms of the data misfit. It is expected that the outcomes of each inversion will appear in an L-shape:there will be branches in which a change in the regularization parameter will yield a significant changein either the data or model misfit. There will be a corner in which the change in the regularization pa- nversion with wavelet-based and scale-dependent regularization ( φ d , φ m ) -plane closest to the origin and thus with both misfits minimized. The corner is usuallymore pronounced and easier to detect in log-log plane. In general, there is no guarantee that the curvewill exhibit an L-shape. Local corners pose difficulties for the automatic corner detection. We rely onthe adaptive pruning algorithm (Hansen et al. 2007), because it is a robust method for the appropriatecorner selection. The basis transformation W transforms a minimum structure model m to a sparse representation x inthe wavelet domain. The wavelet transform has both spatial and temporal resolution and allows to rep-resent the conductivity profile in a sparse fashion. By studying the wavelet transform, we can analyseboth the location of the peaks (spatial resolution) and periodicities at different frequencies (temporalresolution). The sparsity of the transform is guaranteed by a theorem presented later in this text. Mostof the material in this section can be found in (Strang & Nguyen 1996) and (Mallat 1989).While the wavelet transform has a complete and independent theory based on the definition ofmultiresolution analysis (Daubechies 1992), we sketch a more intuitive summary with analogies toFourier series. In Fourier series, a π periodic signal in the continuous time domain is represented interms of basis functions { e ikt } k ∈ Z with Fourier coefficients f k = (cid:104) e ikt , f (cid:105) (Mallat 1999). Indeed, wewrite f ( t ) = (cid:88) k ∈ Z f k e ikt . (9)Wavelet theory constructs a set of basis functions, consisting of scaling functions ϕ and wavelet func-tions ψ . Analogously to Fourier series, the signal is expanded in a set of basis functions. Unlike Fourierseries, these basis functions have compact support (i.e. they only exist on a finite interval) and f doesnot need to be periodic. Finding such bases and determining their properties are the major matter of in-terest of wavelet theory. Instead of considering sines and cosines at different frequencies, the followingtranslations and dilations of the wavelet and scaling functions span the function space: ϕ n,k ( t ) = 2 n/ ϕ (2 n t − k ) or ψ n,k ( t ) = 2 n/ ψ (2 n t − k ) , (10)where n ∈ Z is the dilation parameter that makes the basis function’s compact support larger or smaller(temporal resolution). The k -parameter describes translations along the t -axis (spatial resolution). Werefer to ψ n or ϕ n as a set of basis functions at resolution level n and these sets span the spaces W n and V n respectively. Assume that the bases are orthonormal (this is not required, biorthogonal wavelets Wouter Deleersnyder can also be used). The scaling function is more fundamental than the wavelet function, the waveletfunction at resolution level n is a linear combination of scaling functions at resolution level n + 1 .As with Fourier series, we can decompose an arbitrary function f ( t ) ∈ ( W n ⊕ V n ) as a series inboth scaling and wavelet functions: f ( t ) = (cid:88) k ∈ Z v nk ϕ nk ( t ) , f ( t ) = (cid:88) k ∈ Z w nk ψ nk ( t ) , (11)where v nk = (cid:104) φ nk , f (cid:105) and w nk = (cid:104) ψ nk , f (cid:105) are scaling and wavelet coefficients respectively. In thediscrete wavelet transform, the function f ( t ) and wavelets are discretely sampled. We drop the time-dependence in this notation and we refer to it as a signal instead of a function.A one level (discrete) wavelet transform is the decomposition of the signal f n ∈ V n in its compo-nents f n − ∈ V n − and g n − ∈ W n − . f n = (cid:88) k v nk ϕ nk = (cid:88) k v n − ,k ϕ n − ,k (cid:124) (cid:123)(cid:122) (cid:125) f n − + (cid:88) k w n − ,k ψ n − ,k (cid:124) (cid:123)(cid:122) (cid:125) g n − (12)where v n − ,k are the scaling coefficients and w n − ,k the wavelet coefficients at resolution level n − .The signal g n − describes the details that were present in the signal on scale n , but have disappearedfrom the coarser ( n − -scale f n − . A higher level wavelet transform is obtained by repeating theone level wavelet transform on f n − . This process can be interpreted in the framework of filter banks,where the wavelet coefficients are obtained after passing through a high pass filter (contains detailsof signal) and a subsequent downsampling procedure. A low pass filter and the subsequent downsam-pling generate the scaling coefficients. These scaling functions are then again passed through the highand low pass filter etc. The high and low pass filter coefficients are a result from specific conditions,for example that they allow for perfect reconstruction to the model space, i.e. such that there exists aninverse procedure to go from the coefficients in wavelet space to the original signal. For every waveletbasis function, there exists a unique set of low pass and high pass filters. Their relation is not obvi-ous, for which we refer to Mallat (1999). However, the filter bank interpretation allows for a fast andreliable computation of the wavelet transform with complexity O ( N ) , in which the explicit form ofthe basis function is not required. This computation, known as the Fast Wavelet Transform (FWT), isanalogous to the Fast Fourier Transform.In our inversion approach, we require an explicit matrix W ∈ R n x × n m for the basis transforma- nversion with wavelet-based and scale-dependent regularization n m wavelet transforms of a Dirac train † W = W I = W [ δ δ · · · δ n m ] = [ W δ W δ · · · W δ n m ] , (13)where W δ i is computed with an available implementation of the FWT, such as the PyWavelets pack-age (Lee et al. 2006) in Python. We distinguish the lengths of the signal in the different domains, morespecifically n x ≥ n m . This results from signal extension, which is an analogous situation as withconvolutions ‡ . With finite signals, one must think about how to handle the boundaries. It is impos-sible to correlate a basis function with a part of a signal when that basis function has wider supportthan the length of that part of the signal. The signal is therefore by default symmetrically or period-ically extended (e.g. in Matlab’s implementation (Misiti et al. 2009)). Another strategy to limit theboundary distortions, is to compute a maximum useful level of decomposition. One often considersthe maximum level N max N max = (cid:22) log (cid:18) n m F − (cid:19)(cid:23) , (14)where (cid:98)·(cid:99) is the floor-operator and F is the number of filter coefficients, which scales with the supportof the wavelet.Finally, some properties about wavelets are introduced to better understand the approximatingabilities of a specific wavelet. The number of vanishing moments is the most decisive property of awavelet. A wavelet has p vanishing moments when (cid:90) t k ψdt = 0 for k = 0 , · · · , p − . (15)The number of vanishing moments is related with the compact support of the wavelet: an orthonormalwavelet with p vanishing moments has at least a support of size p − (Daubechies 1988). Waveletswith minimal compact support for a given p are called Daubechies wavelets (db p ). Further, if f ( t ) is p times differentiable, its wavelet coefficients decay (Strang & Nguyen 1996): Theorem 2.1.
Decay of the Wavelet: If f ( t ) is p times differentiable, its wavelet coefficients decaylike − jp : | w j,k | ≤ C − jp || f ( p ) ( t ) || . (16)Hence, wavelet coefficients decay with increasing resolution level. Inversely, discontinuities or sin-gularities in the signal f yield large wavelet coefficients corresponding to the same spatial location † A Dirac train δ i is a vector with zeros, and with a spike (a one) at the i -th index. ‡ In the filter bank interpretation, the relation between v nk and v n +1 ,k is a convolution with a downsampled and time reversed impulseresponse. Wouter Deleersnyder . . . (a) Scaling function ϕ − (b) Wavelet function ψ . . . . . × − (c) Signal f Signal f Sparse f . . . . × − (d) 1 level DWT . . . (e) 2 level DWT . . . . (f) 3 level DWT .
00 0 .
25 0 .
50 0 .
75 1 . (g) Signal g − (h) 4 level DWT of g − − v , − ϕ , − + v , ϕ , + v , ϕ , g (i) Reconstruction of signal g Figure 1.
Example of the wavelet transform with Daubechies wavelets with two vanishing moments (db2) ontwo discrete signals f and g . Vertical lines in figures (d)-(f) show the boundary between scaling coefficients andwavelet coefficients at different resolution levels. on all resolution levels. This theorem, therefore, reveals an important trade-off: higher vanishing mo-ments improve the approximating abilities of a wavelet (which is interesting in applications as it yieldssparser representations) but involves wavelets with larger compact support. The larger the compactsupport, the more difficult to meet the conditions of the theorem: the signal f must then be piecewisesmooth on a larger interval. Daubechies wavelets are by definition the wavelets with best approxi-mating abilities: they are for a given number of vanishing moments, the orthogonal wavelet with thesmallest compact support.The Daubechies wavelet with two vanishing moments (db2) is shown in Fig 1. Their scaling andwavelet functions in 1(a) and 1(b) have a support width of p − or 3. Figs 1(d)-(f) show respectivelythe one, two and three level discrete wavelet transform of signal f in Fig 1(c). Signal f is a discretelysampled conductivity profile obtained from (Hermans & Irving 2017). It exhibits minimum-structurein the sense that neighbouring samples have similar electrical conductivities. The samples can beviewed as scaling coefficients v n at resolution level n . The one level DWT (with by default symmetric nversion with wavelet-based and scale-dependent regularization Space domain:Initial guess m (0) log -space: m (0) , log = log m (0) Wavelet- log domain: x log , (0) = W m (0) , log Inversion: Solve φ ( x ) with L-BFGS-B methodWavelet- log domain: x out,log log -space: m out,log = W − x out,log Space domain: m out = 10 ( m out,log ) Figure 2.
The inversion scheme in log -space domain signal extension) transforms these scaling coefficients to scaling coefficients v n − and wavelet coef-ficients w n − at a coarser resolution level n − . The former describes the signal at a coarser scale,while the latter coefficients contain the details. In Fig 1(d), the scaling coefficients, the 20 coefficientsin the first half of the figure, describe the original signal with less detail (and are also rescaled with afactor √ , see Equation (10)). One can immediately notice that for such minimum-structure signalsthe wavelet coefficients, in the second half of the figure, are small. When the small wavelet coefficientsare ignored (set to zero), one gets a sparser representation of the original signal. The back transforma-tion to model space (sparse signal f or dashed line in Fig 1(c)) is still be very similar to the originalsignal f . This is obtained by setting all wavelet coefficients equal to zero if it is smaller than 5% ofthe largest scaling coefficient. This common practice in wavelet compression yields a representationwith all scaling coefficients and only one wavelet coefficient. To go from the one level to the two levelDWT, the same procedure is applied on the scaling coefficients v n − (the wavelet coefficients w n − remain unaltered). The scaling coefficients v n − again describe the original signal at an even coarserscale and w n − describe the details that are lost since v n − . Fig 1(f) shows the three level DWT. Notethat the number of coefficients in wavelet space is slightly larger than in model space. This is due thesignal extension at the boundaries. In general n m ≤ n x, ≤ n x, ≤ n x, , where n m is the length of the signal f in model space and n x,i the i-th level DWT of signal f . Waveletswith larger support widths yield wavelet representations with more coefficients.A wavelet with p vanishing moments is orthogonal to polynomials of degree p − (see Eq. (15)).Hence, the db2 wavelet is orthogonal to linear functions. This guarantees that wavelet coefficients willbe zero for linear pieces in the signal f . The greater the number of vanishing moments, the more com-plex structures can be represented in a sparse fashion. This clearly shows that wavelets with highervanishing moments have better approximating abilities. In Figs 1(j)-(i) we consider a discretely sam-2 Wouter Deleersnyder pled linear function g and its maximum four level wavelet transform. One immediately notices that allwavelet coefficients (right-side) are zero. In the scaling coefficients, one can still recognize the original(rescaled) structure. The minimum-structure linear function can consequently be represented exactly with only three coefficients in the wavelet representation. To increase intuition, we return to Equation(9) in which we expand a signal into a basis. It is not obvious by eye that a linear combination ofthe irregular shaped db2 scaling function yields the original function. In this example, there are onlythree terms: three scaling functions with support width 3. In Fig 1(i), we show each scaling functionmultiplied by its scaling coefficient and in the correct location (depends on the translation parameter k ). If we add the signals, we get a linear function in the original interval [0 , . Note that we usedsmooth signal extension here, as a result of which we effectively avoided boundary distortions. Finally, we introduce what we call scale-dependent regularization. The problem with the regularizationterm in equation (7) is that it gives each coefficients in x the same weight. However, as we know thatsome components are more likely to be zero than others, we account for this in the objective function.The first coefficient will never be zero, because it would mean that the integrated value over theconductivity profile would be zero. Wavelet coefficients at smaller scales are expected to be zero, sincethese correspond with neighbouring model parameters having equal conductivities. Our wavelet basis,however, exhibits better minimum structure when the high resolution wavelet coefficients x i ∈ W − are sparse. Theorem 2.1 makes this more formal for any type of wavelet with p vanishing moments,recall that if f is locally smooth and p times differentiable, the theorem states that at scale j + 1 the wavelet coefficients, localized where f is smooth and p times differentiable, are approximatelysmaller than those on scale j by a factor of p . This result can be used to define a new regularizationterm (where we assume no signal extension): φ m ( x ) = 1 E (cid:32) µ ( x ) + 2 (cid:88) i =3 µ ( x i )+2 (cid:88) i =5 µ ( x i ) + 2 (cid:88) i =9 µ ( x i ) + · · · (cid:33) , (17)where we have made the sparsity constraint for high resolution coefficients more stringent and µ is theEkblom measure. Note that we have dropped the number of vanishing moments p in the exponential asin Theorem 2.1, because this would be an excessive penalization of the small-scale wavelet coefficientsfor a large number of vanishing moments p (the results for j are in general better than for jp ). nversion with wavelet-based and scale-dependent regularization E ) of coefficients tobe able to better compare between different parametrizations. Solving the inverse problem is basically minimizing an objective function. Constructing a suitablefunction, however, requires some thought and multiple assumptions must be made. Our inversionscheme is illustrated in Fig 2. An iterative optimization method starts with an initial model or guess.In principle, this model can be based on prior knowledge of the geological context (e.g. groundtruthknowledge obtained in the vicinity of the surveying site) or by relying on the apparent conductivity,although the initial model can likewise be generated randomly. The latter option is less suitable formultimodal objective functions.Following the generation of the initial model, it is transformed into the logarithmic domain. Thisis a common action in geophysical inversion schemes to ensure positive parameters, such as for theelectrical conductivity. Of course, we will conduct a back-transform to the original domain later in theinversion scheme, where all values from the logarithmic domain are transformed to the R + domain.The model in the logarithmic domain is subsequently transformed into wavelet domain by thediscrete wavelet transform. This involves a few preliminary choices that must be made before thetransform can be applied. Each decision affects the objective function and hence potentially affectsits minimum. First of all, one must decide which wavelet to use. Since in this paper only one dimen-sion is considered, we choose wavelets from the Daubechies family for their approximating abilities,given the trade-off between compact support width and the number of vanishing moments. A seconddecision is the level of the wavelet transform. In principle, we require a sparse representation of ageologically realistic model, because we work with a sparsity promoting measure. We have alreadyargued in the previous section that for a signal of length N , it is not necessarily better to apply an N level DWT, due to boundary distortions. Nevertheless, we choose N = log ( n m ) which turns out towork equally well. To minimize those boundary distortions, the choice which signal extension modeto use plays a role. In a geological context, we generally do not expect the model parameters to beperiodic, so we opt for symmetrization or a smooth signal extension. Finally note that the shape of thewavelet will affect the result of the inversion, especially if we opt for a strong degree of regularization.Accordingly, if a blocky model is expected, it is better to opt for blocky wavelets which have sharperedges, such as the db1 wavelet. If soft boundaries are expected between the layers, then smootherwavelets are more appropriate. These wavelets usually have a larger number of vanishing moments4 Wouter Deleersnyder . . . . . . . . . d (m) . . . . . σ ( S m − ) Wavelet: db1 modelconst. λ , λ = 1.6e-05 j λ , λ = d (m) . . . . . σ ( S m − ) Wavelet: db6
True modelconst. λ , λ = 1.0e-07 j λ , λ =3.2e-07 Figure 3.
Comparison of traditional vs. scale-dependent wavelet-based regularization schemes with the LINapproximation as forward operator. and have a stabilizing effect due to their larger compact support.There are various optimization methods available, classified into trust-region, line-search and non-gradient methods. In this paper we exclusively use the BroydenFletcherGoldfarbShanno (BFGS) al-gorithm. This algorithm is implemented in Python (SciPy) (Jones et al. 2001) with the line search al-gorithm that meets the Wolfe conditions (Mor´e & Thuente 1994). This method is a Newton’s methodand approximates the Hessian at every iteration. In practice, the choice of tolerances (stop criteria)will affect the actual outcome. In general, the tolerance on the objective function is set to − and − for the gradient is sufficient (for further details, see (Jones et al. 2001)). In the first examples, we demonstrate the advantage of the scale-dependence of our novel regulariza-tion scheme. Whereas other schemes (such as (Liu et al. 2017)) penalize each wavelet coefficient withequal weight, our scheme penalizes each wavelet coefficient proportionally with their resolution level n . To make a valid comparison between the two wavelet-based regularization schemes, we first elim-inate all other factors in the inversion process that affects the outcome. Such factors are, for example,the presence of numerous local minima in the objective function, the presence of a modelling error κ in the forward operator or a limited maximal number of iterations due to limited computationalresources. We eliminate the modelling error by using the same forward operator for the synthetic data nversion with wavelet-based and scale-dependent regularization j λ ) and the traditional sparsity inversion with equalregularization (const. λ ). The L-curve criterion is used for the estimation of the optimal regularizationparameter. While the convexity of the L-curve is not guaranteed, the L-curves displayed the typicalL-shape (not shown).The synthetic data is generated for a vertical dipole in both components (horizontal coplanar HCPand perpendicular PERP) for 20 equidistant intercoil distances from 1 to 20 meters at a height of 0.1meters above the soil and where we have added 5% Gaussian multiplicative noise. We choose to fit32 model parameters. The result of the scale-dependent sparsity inversion outperforms the traditionalwavelet-based inversion. The scale-dependent scheme follows the true model almost exactly, whilethe latter has too much structure at the small scales. This example illustrates the advantage of scale-dependent regularization over traditional wavelet-based inversion with db1 wavelets.The second example shows a conductivity profile from a borehole logging from an alluvial aquifer(Hermans & Irving 2017). It has ten times smaller conductivity values and therefore the variation inthe magnetic field data will be more subtle, especially compared with the noise level of 5%. In thiscase, we opt for a more realistic undetermined problem in which we fit 64 model parameters with thedb6 wavelet to data from a vertical dipole for both components and 5 equidistant intercoil distancesfrom 1 to 10 meters. The results in Fig 3 show for both schemes a wrong location of the main peak at2.5m depth. However, we observe that the scale-dependent regularization has a much narrower peakand whose amplitude is closer to the exact values indicating that the scale-dependent regularizationseems also better suited for other wavelet bases. The misidentification of the peak is related to thecombination of the noise level with the challenging nature of the inverse problem. The inversion ofthe noise-free data set (shown in the Appendix B in Fig A2) yields a correct location of the peak (for6 Wouter Deleersnyder . . . . . . . . . d (m) . . . . . σ ( S m − ) (a) modeldb1, λ = ‘ , λ = . . . . . . . . . d (m) . . . . . . σ ( S m − ) (b) j λ , λ = modeldb1db2db3 . . . . . . . . . d (m) . . . . . . σ ( S m − ) (c) j λ , λ = modeldb4db5db6 . . . . . . . . . d (m) . . . . . . σ ( S m − ) (d) j λ , λ = modeldb7db8db9 Figure 4.
Models after inversion with no modelling error ( κ = 0 ) and the adaptive wavelet-based regularizationwith various Daubechies wavelets vs. Tikhonov regularization. both regularization methods). d (m) . . . . . σ ( S m − )
5% noise modeldb1, λ = λ = ‘ , λ = Figure 5.
Models after inversion with no modelling error ( κ = 0 ) and the adaptive wavelet-based regularizationwith the db1 and db6-wavelet vs. Tikhonov regularization. nversion with wavelet-based and scale-dependent regularization d (m) . . . . . . . σ ( S m − )
32 parameters
True modeldb1, λ = ‘ , λ = d (m) . . . . . . . σ ( S m − )
64 parameters
True modeldb2, λ = ‘ , λ = d (m) . . . . . σ ( S m − )
64 parameters
True modeldb3, λ = ‘ , λ = d (m) . . . . . σ ( S m − )
128 parameters
True modeldb6, λ = ‘ , λ = Figure 6.
Models after inversion with no modelling error ( κ = 0 ) and the adaptive wavelet-based regularizationwith various Daubechies vs. Tikhonov regularization. In this section, we illustrate our proposed regularization scheme and take the traditional scheme withsmoothing Tikhonov regularization as benchmark. Our adaptive scheme can recover blocky structuresand sharp boundaries, which is lacking in minimum-structure inversions with smoothing constraints.To illustrate the advantage of wavelet space over standard model space, we perform a numerical teston the three-layered model described in the previous section, with the only difference that we fit 64model parameters instead of 32. The results after inversion are shown in Fig 4. From the shapes ofthe wavelet, it is expected that the db1-wavelet yields the best (blocky) inversion. We therefore com-pare the results of the db1-wavelet with smoothing inversion in Fig 4(a). The wavelet-based inversionfollows the model closely, while the smoothing inversion locates the peak correctly, but with an under-estimated value for the electrical conductivity as expected. The smooth inversion also creates a highconductivity artefact close to the surface. In addition, the interface between the two different layers isnot identified. Figs 4(b) to 4(d) show inversion results for wavelets with increasing vanishing moment(and increasing regularity). To facilitate result comparison, the regularization parameter was fixed to − for all inversions (The L-curve criterion selected very similar regularization parameters for each8 Wouter Deleersnyder wavelet). As expected from the shape of the wavelet, the db2-wavelet is capable of both showing sharpand relatively smooth transitions. For other wavelets, the inversion results are relatively identical andsmoother, but with a larger peak than commonly found with Tikhonov regularization. So far, we couldnot determine a method to select an optimal wavelet from the data, as this remains strongly dependenton the expected geological variations. On the other hand, one can impose the characteristics of theinversion scheme as desired. For blocky inversion, one chooses a low number of vanishing moments,while for smooth inversion, a larger vanishing moment can be opted for. Note that in general the max-imal electrical conductivity value is recovered more exactly with wavelet inversion than for smoothTikhonov inversion (i.e. large number of vanishing moments) and that the results are more symmet-rical near the actual conductive layer. The only disadvantage is the extra structure at larger depths.Those are oscillations which are a result of the shape of the wavelet (the more vanishing moments, themore oscillations as can be seen from Fig A1 in Appendix A).The scale-dependent wavelet based inversion is advantageous for recovering smoother inversionas well. As in the previous section, we compare the inversion with the conductivity profile from thealluvial aquifer. For each inversion, the L-curve criterion determines the optimal regularization param-eter which is different for each wavelet. The results for Tikhonov, db1 and db6 wavelet inversion areshown in Fig 5. The results for db1 and Tikhonov are similar. Both schemes recover peaks at the samelocations with more or less the same values for electrical conductivity. As with the previous test withthis profile, the wavelets with more vanishing moments can recover the peaks better. The case withthe correct location of the peak in the absence of experimental noise is shown in Fig A3 in Appendix B.The shape of the wavelet (and scaling function) affects the inversion result, and choosing thewavelet should be based on prior knowledge. This can be understood from the theoretical framework,namely by penalizing large coefficients of the wavelet representation of a model, the outcome is alteredwith exactly the shape of the wavelet or scaling function. Therefore, the shape is often recognized inthe inversion result. Consider another example with a softer interface between two different layers. InFig 6, such a conductivity profile from (Hermans et al. 2012) is shown from De Panne, Flanders, wherea saltwater lens is lying on a clay layer. We show the inversion results for 4 different wavelets: db1,db2, db3 and db6. In general, better results are obtained when more model parameters are consideredfor more vanishing moments. For each case, the Tikhonov smoothing inversion is also shown in thecorresponding number of model parameters. For all those cases, the Tikhonov regularization cannotfit the actual peak, while the wavelet inversion can, independently from the chosen wavelet. We againobserve a decreasing blockyness for larger vanishing moments and an artificial second peak for db3 nversion with wavelet-based and scale-dependent regularization
In the previous section in the absence of unmodelled errors ( || κ || = 0) , our proposed inversion schemesuccessfully recovered blocky structures. Next, we test the effectiveness of our scheme with an approx-imated forward operator F approx instead of the exact computation F , for the purpose of computationalefficiency. The synthetic data is generated from the operator F to mimic a true survey (including 5%multiplicative Gaussian noise), hence we examine our scheme in the presence of biased error κ . Usingthe scheme for both the LIN approximation and the damped model allows comparing the effect of thedegree of unmodelled errors on the result of the inversion, especially in a conductive setting where theLIN assumption breaks down.The magnetic field data generated from the exact model F , the LIN approximation F LIN and thedamped model F damped are shown in Fig 7. From this data, it is clear that the noise is biased andlarger for the LIN approximation than for the damped model, i.e. their error is || κ LIN || = 0 . and || κ damped || = 0 . respectively. As expected, we get || κ damped || < || κ LIN || . In terms of data-misfit(3), the unbiased experimental noise equals 0.000149, the biased unmodelled error in damped andexperimental noise equals 0.0002664, while for the LIN approximation this equals 0.00996. Note thatthe error and its norm will differ for different conductivity profiles.The L-curve criterion can easily locate vertical branch due to the experimental random noise, butwe can not use this principle to estimate a regularization parameter such that the data misfit meets thevalue of the sum of the experimental and biased noise. Furthermore, we cannot rely on the discrep-ancy principle, since the biased noise κ is unknown. The optimal estimated regularization parametersfollowing the L-curve criterion yield the estimates in the range of the experimental noise. Their datamisfits φ d for the LIN and damped model are 0.000217 and 0.000143 respectively. The damped modelhas a misfit approximately equal to the experimental noise level, while for LIN this is slightly larger.The results are shown in Fig 8. The conductivity profile obtained via the damped model has only aslightly deviant electrical conductivity, while for the LIN approximation slightly too much structurewas introduced. However, both models recover the profile relatively well. Surprisingly, the smoothinversion with the LIN approximation yields better results. This example illustrates the importance of0 Wouter Deleersnyder
Intercoil distance s (m) . . . . . M a g n e t i c fi e l d r a t i o FF damped F LIN
Figure 7.
Magnetic field ratio generated via the different forward operators from the conductivity profile in Fig8. a good approximation for wavelet-based inversion, yet, our method seems relatively robust.
In this section, we demonstrate the wavelet-based inversion scheme on real FDEM field data, obtainedfrom (Bobe et al. 2020). The measurement location is in Gontrode forest, Belgium, known for its hori-zontal stratigraphy of different soil layers (i.e. no 2D or 3D effects are expected). For the sounding, theDUALEM 421S was used with six receiver coils (three HCP coils at 1,2 and 4 metres and three PRPcoils at 1.1,2.1 and 4.1 metres) with an operating frequency of 9.0 kHz. The sounding was performedat zero, 30 cm and 60 cm height. The magnetic field ratios are shown in Fig 9.Bobe et al. (2020) conducted a visual inspection of the subsurface for the validation of the results,by drilling a hole with an Edelman auger. The top 0.15 meters consists of organic rich forest material . . . . . . . . d (m) . . . . . . σ ( S m − ) Wavelet: db1 modeldb1, F LINdb1, F damped . . . . . . . . d (m) . . . . . σ ( S m − ) Smooth inversion model ‘ , F LIN ‘ , F damped Figure 8.
Effect of forward operator (and modelling error κ ) on the inversion result. nversion with wavelet-based and scale-dependent regularization . . . . . . . Intercoil distance s M a g n e t i c fi e l d r a t i o ( % ) Field data
HCP, 0 mHCP, 0.3 mHCP, 0.6 mPRP, 0 mPRP, 0.3 mPRP, 0.6 m
Figure 9.
FDEM field data obtained with the DUALEM 421S, from Bobe et al. (2020). with many air-filled voids. Between 0.15 and 0.8 metres, a clay-rich sand was observed. The claycontent decreased with further depth. At 1.5 metres, the groundwater table was reached.The result of the inversion procedure with the damped model as forward operator is shown fordb1,db2,db3 and db6 wavelets and for smooth inversion in Fig 10. The results of the inversion withother wavelets exhibited similar characteristics. We have used an equidistant parameterization of 128model parameters up to 10 m depth and only the outcome of the first three meters is shown (becauseof the little sensitivity in the lower depth region). The results can be compared to the outcome in figure7 (red line) in Bobe et al. (2020), where with FDEM data alone, the clay-rich sand was not recovered.It shows a steady increase in electrical conductivity from 10 mS/m at the surface to almost 100 mS/mat two meters depth, followed by a constant plateau.We expect the first 0.15 metres to be highly resistive (air voids are an insulator). Clay is known to . . . . . . . d (m) − − σ ( m S m − ) ‘ db6 . . . . . . . d (m) σ ( m S m − ) db1db2db3 Figure 10.
Inversion results from Tikhonov regularization ( (cid:96) ) and wavelet-based regularization withDaubechies (db) wavelets. Estimation of regularization parameter via L-curve criterion. Wouter Deleersnyder have a relatively high bulk electrical conductivity value, which explains the peak in the conductivityprofiles. Sand, on the other hand, has a low bulk electrical conductivity. Groundwater contains morefree electrical charges and therefore, the electrical conductivity usually increases in the saturated zone.This behaviour is recovered in almost all wavelet-based inversion results, however a peak as observedin db1 and db2 is not necessarily expected. Based on the lithological profile, db3 seems to providethe more realistic results. The main difference with smooth inversion from Bobe et al. (2020) is theincrease in electrical conductivity due to the water table, which is not present in the smooth inversionand inversion with db6 and the unrealistic presence of the low conductivity between 15 and 50 cmdepth. .The wavelet-based regularization scheme yields results with similar characteristics, yet the resultsare different. This is the faith of deterministic inversion, however the wavelet-based inversion schemewould also suit in stochastic inversion methods.
We have introduced an improved inversion scheme for EMI surveys that can be extended to any other1D geophysical method. It involves a new model misfit or regularization term based on the wavelettransform and scale-dependent weighting which can easily be combined with the existing frameworkof deterministic inversion (gradient-based optimization methods, L-curve criterion for optimal regu-larization parameter).Our results show that the scale-dependent wavelet-based inversion scheme with sparsity constraintis more adaptive than Tikhonov regularization since it can recover both blocky and smooth conductiv-ity profiles, by adequately choosing the number of vanishing moments of the wavelet. Furthermore,our scheme can recover high amplitude anomalies in combination with globally smooth profiles whatis generally not possible, as discussed in Section 3.2. The scale-dependency of our scheme allows touse of wavelets with few vanishing moments and is, therefore, an improvement with respect to exist-ing wavelet-based regularization schemes. The adaptive nature of the inversion method allows for highflexibility because the shape of the wavelet can be exploited to generate multiple representations of theinverse model. Depending on available prior knowledge, the final result can be chosen or the resultscan be interpreted together. Although the choice os the wavelet is case-specific, we have observed fewvariations in the recovered model above db6. Together with the computational speed-up, our methodthus offers a viable alternative to the existing methodologies, as discussed in the Introduction. nversion with wavelet-based and scale-dependent regularization
ACKNOWLEDGMENTS
The authors thank Christin Bobe for making available the data from the field work in Gontrode, Bel-gium and Marieke Paepen for the conversations on the inversion of those data. The research leadingto these results have gratefully received funding from FWO (Fund for Scientific Research, Flanders,grant 1113020N).
REFERENCES
Archie, G. E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics,
Transactions of the AIME , (01), 54–62.Bobe, C., Hanssens, D., Hermans, T., & Van De Vijver, E., 2020. Efficient probabilistic joint inversion ofdirect current resistivity and small-loop electromagnetic data, Algorithms , (6), 144.Bunks, C., Saleck, F. M., Zaleski, S., & Chavent, G., 1995. Multiscale seismic waveform inversion, Geo-physics , (5), 1457–1473.Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., & Oldenburg, D. W., 2015. SimPEG: An open sourceframework for simulation and gradient based parameter estimation in geophysical applications, Computers& Geosciences , , 142–154.Daubechies, I., 1988. Orthonormal bases of compactly supported wavelets, Communications on pure andapplied mathematics , (7), 909–996.Daubechies, I., 1992. Ten lectures on wavelets , vol. 61, Siam.De Smedt, P., Van Meirvenne, M., Meerschman, E., Saey, T., Bats, M., Court-Picon, M., De Reu, J., Zwert-vaegher, A., Antrop, M., Bourgeois, J., et al., 2011. Reconstructing palaeochannel morphology with a mobilemulticoil electromagnetic induction sensor,
Geomorphology , (3-4), 136–141.Deidda, G. P., Fenu, C., & Rodriguez, G., 2014. Regularized solution of a nonlinear problem in electromagneticsounding, Inverse Problems , (12), 125014.Ekblom, H., 1987. The l1-estimate as limiting case of an lp-or huber-estimate, in Statistical data analysisbased on the L1-norm and related methods: 31/08/1987-04/09/1987 , pp. 109–116, Elsevier.Farquharson, C. G., 2007. Constructing piecewise-constant models in multidimensional minimum-structureinversions,
Geophysics , (1), K1–K9.Fern´andez, J. P., Shubitidze, F., Shamatava, I., Barrowes, B. E., & O’Neill, K., 2010. Realistic subsurfaceanomaly discrimination using electromagnetic induction and an svm classifier, EURASIP Journal on Ad-vances in Signal Processing , (1), 305890. Wouter Deleersnyder
Guillemoteau, J., Simon, F.-X., L¨uck, E., & Tronicke, J., 2016. 1d sequential inversion of portable multi-configuration electromagnetic induction data,
Near Surface Geophysics , (5), 423–432.Guitton, A., 2012. Blocky regularization schemes for full-waveform inversion, Geophysical Prospecting , (5), 870–884.Hansen, P. C., 2010. Discrete inverse problems: insight and algorithms , vol. 7, Siam.Hansen, P. C., Jensen, T. K., & Rodriguez, G., 2007. An adaptive pruning algorithm for the discrete l-curvecriterion,
Journal of computational and applied mathematics , (2), 483–492.Hermans, T. & Irving, J., 2017. Facies discrimination with electrical resistivity tomography using a proba-bilistic methodology: effect of sensitivity and regularisation, Near Surface Geophysics , (1), 13–25.Hermans, T., Vandenbohede, A., Lebbe, L., Martin, R., Kemna, A., Beaujean, J., & Nguyen, F., 2012. Imag-ing artificial salt water infiltration using electrical resistivity tomography constrained by geostatistical data, Journal of Hydrology , , 168–180.Jadoon, K. Z., McCabe, M. F., & Moghadas, D., 2015. Application of electromagnetic induction to monitorchanges in soil electrical conductivity profiles in arid agriculture, in First Conference on Proximal SensingSupporting Precision Agriculture , vol. 2015, pp. 1–5, European Association of Geoscientists & Engineers.Jones, E., Oliphant, T., Peterson, P., et al., 2001. SciPy: Open source scientific tools for Python, [Online;accessed 12-01-2019].Lee, G., Wasilewski, F., Gommers, R., Wohlfahrt, K., O’Leary, A., & Nahrstaedt, H., 2006. Pywavelets–wavelet transforms in python.Linde, N., Renard, P., Mukerji, T., & Caers, J., 2015. Geological realism in hydrogeological and geophysicalinverse modeling: A review,
Advances in Water Resources , , 86–101.Liu, Y., Farquharson, C. G., Yin, C., & Baranwal, V. C., 2017. Wavelet-based 3-d inversion for frequency-domain airborne em data, Geophysical Journal International , (1), 1–15.Mallat, S., 1999. A wavelet tour of signal processing , Elsevier.Mallat, S. G., 1989. A theory for multiresolution signal decomposition: the wavelet representation,
IEEETransactions on Pattern Analysis & Machine Intelligence , (7), 674–693.Maveau, B., Delrue, S., & Dudal, D., 2020. A damped forward EMI model for a horizontally stratified earth, Exploration Geophysics , pp. 1–12.McNeill, J., 1980. Electromagnetic terrain conductivity measurement at low induction numbers.Misiti, M., Misiti, Y., Oppenheim, G., & Poggi, J., 2009. Matlab wavelet toolbox tm 4 user’s guide,
TheMathWorks, Inc. Natick, Massachusetts. 153p .Mor´e, J. J. & Thuente, D. J., 1994. Line search algorithms with guaranteed sufficient decrease,
ACM Transac-tions on Mathematical Software (TOMS) , (3), 286–307.Paepen, M., Hanssens, D., Smedt, P. D., Walraevens, K., & Hermans, T., 2020. Combining resistivity andfrequency domain electromagnetic methods to investigate submarine groundwater discharge in the littoralzone, Hydrology and Earth System Sciences , (7), 3539–3555.Pettersson, J. K. & Nobes, D. C., 2003. Environmental geophysics at scott base: ground penetrating radar nversion with wavelet-based and scale-dependent regularization and electromagnetic induction as tools for mapping contaminated ground at antarctic research bases, Coldregions science and technology , (2), 187–195.Portniaguine, O. & Zhdanov, M. S., 1999. Focusing geophysical inversion images, Geophysics , (3), 874–887.Saey, T., De Smedt, P., Meerschman, E., Islam, M. M., Meeuws, F., Van De Vijver, E., Lehouck, A., &Van Meirvenne, M., 2012. Electrical conductivity depth modelling with a multireceiver emi sensor forprospecting archaeological features, Archaeological Prospection , (1), 21–30.Simpson, D., Van Meirvenne, M., Saey, T., Vermeersch, H., Bourgeois, J., Lehouck, A., Cockx, L., & Vitha-rana, U. W., 2009. Evaluating the multiple coil configurations of the em38dd and dualem-21s sensors todetect archaeological anomalies, Archaeological Prospection , (2), 91–102.Strang, G. & Nguyen, T., 1996. Wavelets and filter banks , SIAM.Strutz, T., 2010.
Data fitting and uncertainty: A practical introduction to weighted least squares and beyond ,Vieweg and Teubner.Tantum, S. L., Scott, W. R., Morton, K. D., Collins, L. M., & Torrione, P. A., 2012. Target classificationand identification using sparse model representations of frequency-domain electromagnetic induction sensordata,
IEEE transactions on geoscience and remote sensing , (5), 2689–2706.Tikhonov, A. N., 1943. On the stability of inverse problems, in Dokl. Akad. Nauk SSSR , vol. 39, pp. 195–198.Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics,
Geophysics , (6), WCC1–WCC26.Wait, J. R., 1951. The magnetic dipole over the horizontally stratified earth, Canadian Journal of Physics , (6), 577–592.Wait, J. R., 1962. A note on the electromagnetic response of a stratified earth, Geophysics , (3), 382–385. APPENDIX A: DAUBECHIES WAVELETS
Scaling functions of the wavelets used in this paper in Fig A1
APPENDIX B: RESULTS WITHOUT EXPERIMENTAL NOISE
Identical inverse problem as in Section 3.1 and 3.2 but without any multiplicative noise in Figs A2 andA3 respectively.
This paper has been produced using the Blackwell Scientific Publications GJI L A TEX2e class file. Wouter Deleersnyder . . . . . . . . . . . . db1 . . . . . . . . . . db2 . . . db3 . . . db4 . . . db5 . . . db6 . . . . . . . . . db7 . . . . . . . − . . . . db8 − . . . . db9 − . . . . db10 − . . . . db11 − . . . . db12 Figure A1.
Scaling functions of the wavelets used in this paper. d (m) . . . . . . . σ ( S m − ) Wavelet: db6
True modelconst. λ , λ = 1.0e-07 j λ , λ =1.0e-07 Figure A2.
Identical inverse problem as in Section 3.1 but without any multiplicative noise. nversion with wavelet-based and scale-dependent regularization d (m) . . . . σ ( S m − )