# Wavelet-based 3-D inversion for frequency-domain airborne EM data

Wavelet-based 3-D inversion for frequency-domain airborne EM data Summary In this paper, we propose a new wavelet-based 3-D inversion method for frequency-domain airborne electromagnetic (FDAEM) data. Instead of inverting the model in the space domain using a smoothing constraint, this new method recovers the model in the wavelet domain based on a sparsity constraint. In the wavelet domain, the model is represented by two types of coefficients, which contain both large- and fine-scale informations of the model, meaning the wavelet-domain inversion has inherent multiresolution. In order to accomplish a sparsity constraint, we minimize an L1-norm measure in the wavelet domain that mostly gives a sparse solution. The final inversion system is solved by an iteratively reweighted least-squares method. We investigate different orders of Daubechies wavelets to accomplish our inversion algorithm, and test them on synthetic frequency-domain AEM data set. The results show that higher order wavelets having larger vanishing moments and regularity can deliver a more stable inversion process and give better local resolution, while the lower order wavelets are simpler and less smooth, and thus capable of recovering sharp discontinuities if the model is simple. At last, we test this new inversion algorithm on a frequency-domain helicopter EM (HEM) field data set acquired in Byneset, Norway. Wavelet-based 3-D inversion of HEM data is compared to L2-norm-based 3-D inversion's result to further investigate the features of the new method. Controlled source electromagnetics (CSEM), Non-linear electromagnetic, Inverse theory, Wavelet transform 1 INTRODUCTION For 3-D electromagnetic (EM) inversion, the geo-electrical model is generally parametrized into regular or unstructured cells. Size of the model are usually not small, so the number of unknowns is much larger than the number of observed data, which always results in an ill-conditioned problem to be solved. In order to make the 3-D inversion stable, some kind of regularization such as Tikhonov regularization needs to be applied (Aster et al.2013). Such regularization methods usually use a smoothing operator to build a relationship between different cells. This technique indeed makes the inversion more stable, but if the regularization is strong, it often smears out small-scale features of the model even where there is sufficient data density available to provide more detailed information. If we use a weak regularization, the resolution may be improved, but the inversion may not be stable. Therefore, it is not easy to get a stable and high-resolution inversion at the same time. For some inversion problems with different site spacing/line spacing, the smoothing regularization may not give a reasonable geo-electrical structure because of the uneven distribution of data. Several techniques such as focusing inversion and L1-norm-based inversion were proposed to improve the resolution of stable gravity and 2-D magnetotelluric (MT) inversion (Portniaguine & Zhdanov 1999; Farquharson 2008). These techniques give the same constraint to the whole inverted model without considering the uneven data distribution. Actually, for the uneven data distribution, the model should be recovered in multiresolution. However, for the model parametrized by regular cell, it is hard to make the model to be adaptive with the uneven data distribution. In seismic, to get different multiscale features of the model according to the ray distribution, several adaptive seismic tomography methods have been proposed in which the cells are irregularly distributed to match with the uneven ray distribution (Abers & Roecker 1991; Spakman & Bijwaard 2001; Zhang & Thurber 2005). These methods work well for recovering multiscale features of the model, but they are based purely on data and do not consider whether the model is appropriately parametrized. The wavelet representation of a signal has an inherent multiscale nature (Daubechies 1992). The wavelet-based multiscale tomography methods are successfully used to recover the multiscale features of the model parametrized by regular cells in seismic applications (Chiao & Kuo 2001; Chiao & Liang 2003; Loris et al.2007; Delost et al.2008; Hung et al.2011; Simons et al.2011; Fang & Zhang 2014). Nittinger & Becken (2016) have also used the sparsity constraint in wavelet domain to accomplish 2-D MT inversion. In fact, before the wavelet is used in seismic tomography for multiscale inversion, it was used in gravity data inversion by Li & Oldenburg (1999) to compress the sensitivity matrix. Later, they made some improvements in wavelet compression-based gravity and magnetic inversions as well (Li & Oldenburg 2003, 2010; Davis & Li 2011). Actually in their research, they already recovered the model in the wavelet domain, but they focused only on the matrix compress and acceleration of the inversion. Besides these applications in linear inversion, Chiao & Kuo (2001) and Chiao & Liang (2003) were the first to propose wavelet-based multiscale inversion for the non-linear tomography inversion problem. In their research, they used the wavelet representation as the model constraints instead of using the common smoothness regularization. The wavelet representation of the model can directly give multiresolution inversion results. For the smoother part of the model, few wavelet coefficients at large scales are significant, while for the part of model with small-scale features and with good data coverage, more wavelet coefficients will be used to sufficiently represent the model. Daubechies et al. (2004) proposed an iterative thresholding method for linear inverse problems using a sparsity constraint. Hung et al. (2010, 2011) proposed a wavelet-based finite-frequency teleseismic tomography method to recover the velocity structure of the central Tibet. Delost et al. (2008) used the wavelet transform for the first-arrival traveltime tomography by designing a bit mask operator to only invert for the wavelet coefficients in specific areas where the resolution is high. Based on work of Daubechies et al. (2004), Loris et al. (2007) and Simons et al. (2011) used the L1-norm of model wavelet coefficients in their inversion scheme to make the solution sparse, which is different from the L2-norm of the model wavelet coefficients used by Chiao & Kuo (2001), Chiao & Liang (2003) and Hung et al. (2010, 2011). Following strategy of Loris et al. (2007) and Nittinger & Becken (2016) accomplished 2-D MT inversion in a sparse model. For the L1-norm inversion of wavelet coefficients, Donoho (2006) proved that it mostly gave a sparse solution and the inversion is a sparsity constraint. Aster et al. (2013) stated that if the model was sparse in nature, the sparsity-constraint inversion could recover it better without a priori smoothing regularization term. Besides the wavelet-based inversion, there are also some other subspace inversion methods. Oldenburg et al. (1993) proposed generalized subspace methods for large-scale inverse problems using basis vectors representation to accelerate the convergence rate. Siripunvaraporn et al. (2005) inverted the MT data using a data-space method which could reduce the inversion space efficiently in comparison to a traditional model-space method. Though these two methods do not have any concept of the multiresolution inversion. In this research, we want to test the wavelet-based multiscale inversion on the EM data that have lower resolution than the seismic data and its sparsity-constraint inversion is less studied. Different from the previous work in seismic, we propose a new wavelet-based inversion method for EM data that uses the model constraints on three directions, which is quite similar to the common smoothing constraint in space domain. We want to further analyse the difference between wavelet-based inversion and traditional Lp-norm inversion method from aspects of both resolution and stability. This new method is derived from a general measure inversion algorithm, and the final equation for this iteratively reweighted least-squares (IRLS) problem is a Gauss–Newton equation solved by conjugate gradient (CG) method. We solve the final equation directly using an IRLS without using the soft threshold method as proposed by Daubechies et al. (2004). We use helicopter EM (HEM) data as an example to test the character of this new inversion method which will be valid in inversion of any frequency-domain airborne EM (FDAEM) data. The AEM survey most of the time has an uneven data distribution. Since the AEM system is compact and has a limited ‘footprint’ (Beamish 2003) for rather larger survey-line spacing, it is hard for conventional 3-D inversion to connect the conductivity structures from two consecutive flight lines. The wavelet-based inversion treats the whole model as a 3-D volume, which main features are recovered. This property can bridge the conductivity structures from two survey lines naturally. For small line spacing, the AEM survey can provide denser data distribution, which can recover fine local structures as well. Based on Liu & Yin's (2013) work, we use the finite-difference method to solve the forward modeling problem and test different kinds of the wavelet basis functions to compare it with other traditional methods on both synthetic and field HEM data set. We further demonstrate the impact of different order of the wavelet basis functions on the inversion result. 2 METHODOLOGY 2.1 General inversion approach Generally, the objective function for an inversion problem is expressed as (Farquharson 2008)   $${\bf \Phi } = {\phi _d} + \lambda {\phi _m}.$$ (1)Here, ϕd and ϕm are the measures of data misfit and model complexity, respectively. The factor λ is the trade-off parameter which balances the contributions of the data misfit and the model complexity term. ϕd can be expressed in a general form as   $${\phi _d} = {\phi _d}({\bf u}),$$ (2)and   $${\bf u} = {{\bf W}_d}({{\bf d}^{\rm obs}} - {{\bf d}^{\rm prd}}),$$ (3)where dobs is the vector of observed data, dprd is the predicted data computed using the model m, and Wd is a diagonal matrix whose elements are the reciprocals of the estimates of the standard deviations of the noise in the observations. ϕm can have the form as   $${\phi _m} = \sum\limits_k {{\alpha _k}} {\phi _k}({{\bf v}_k}),$$ (4)and   $${{\bf v}_k} = {{\bf W}_k}({{\bf m}} - {\bf m}_k^{{\rm{ref}}}).$$ (5)Here, m is the vector of the recovered model parameters, mkref is the reference model, Wk is the smoothing operator matrix which is used to constrain the complexity of the recovered model and αk are the coefficients used to balance the contributions of different kinds of measure of the recovered model. A general form of ϕd and ϕk is   $$\phi ({{\bf x}}) = \sum\limits_j {\rho ({x_j})} ,$$ (6)where xj are the elements of the vector x, which can be either u or vk. There are numerous choices for the specific form of ρ(xj), for example, L2-norm, L1-norm and so on. In the traditional inversion approach, differentiating both sides of eq. (1) with respect to the space-domain model parameters and neglecting the higher order terms of the Taylor-series expansion gives the linear system of equations to be solved at each iteration as   \begin{eqnarray} &&\left[{{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf J}} + {\lambda ^n}\sum\limits_k {{\alpha _k}{{\bf W}}_k^T{{{\bf R}}_k}{{{\bf W}}_k}} } \right]\delta {{\bf m}}\nonumber\\ &&= {{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad +\,\, {\lambda ^n}\sum\limits_k {{\alpha _k}{{\bf W}}_k^T{{{\bf R}}_k}{{{\bf W}}_k}} ({{\bf m}}_k^{{\rm{ref}}} - {{{\bf m}}^{n - 1}}), \end{eqnarray} (7)where the matrix R is a diagonal matrix which has the form as   $${{\bf R}} = {\rm diag}\left\{ {\begin{array}{*{20}{c}} {{{\rho '({x_j})} / {{x_j},}}}&{ \cdots \cdots ,}&{{{\rho '({x_N})} / {{x_N},}}} \end{array}} \right\}.$$ (8) J is the sensitivity matrix in the space-domain model parameters and it can be calculated by an adjoint modeling technique (Liu & Yin 2016), dn−1 is the vector of the data calculated from the model mn−1 which is obtained from the previous iteration and δm is model update. We apply the ‘moving footprint’ method of Cox et al. (2010) to accelerate the forward modeling. 2.2 Wavelet-domain inversion with sparsity constraint In the smoothness constrained inversion, we usually use an L2-norm measure for the model's complexity to get a smooth-edged model or an L1-norm to get a sharper edged model or some other focusing measure to get a somewhat ‘focused’ model with greater boundary resolution. In numerical tests, an L2-norm inversion is typically more stable than an L1-norm inversion, particularly if data density is low. It is hard to combine the advantages of L2- and L1-norms in one single-inversion scheme. The wavelet transform can divide one signal into two groups of coefficients: one is the approximate coefficients and the other is the detail coefficients. This means that this transform has the potential to resolve the model at different levels at the same time according to the data set distribution. This property of wavelet transform motivates us to perform the inversion in the wavelet domain. A space-domain model m can be transformed into wavelet-domain model coefficients $${{\bf \tilde{m}}}$$ using a wavelet transform matrix Ww as   $${{\bf \tilde{m}}} = {{{\bf W}}_w}{{\bf m}}.$$ (9) Since the wavelet transform usually demands that the dimension of the vector should be a power of two, we first extend the three dimensions of m to be power of two, then we apply the 1-D partial wavelet transform scheme of Press et al. (1996) to each dimension of the 3-D model sequentially. This wavelet transform algorithm is a filtering process and the filtering coefficients can be found in Daubechies (1992). According to the chain rule, the sensitivity matrix in the wavelet domain is computed as   $${{\bf \tilde{J}}} = {{\bf JW}}_w^{{\rm{ - 1}}},$$ (10)where $${{\bf W}}_w^{ - 1}$$ is the inverse wavelet transform operator. Since the wavelet coefficients are usually sparse, we use a sparsity constraint in the wavelet domain instead of the smoothing constraint in the space domain (Daubechies et al.2004). According to Donoho (2006), the L1-norm measure can lead to a sparse solution, so we simply use the L1-norm measure on the wavelet coefficients and set k = 1 and α1 = 1 in eq. (4). For the data misfit term, we still use the L2-norm. Then, the linear system (eq. 7) to be solved at each iteration of the inversion becomes   \begin{eqnarray} &&\left[ {{{{{\bf \tilde{J}}}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf \tilde{J}}} + {\lambda ^n}{{\bf R}}} \right]\delta {{\bf \tilde{m}}}\nonumber\\ &&= {{{{\bf \tilde{J}}}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad+\,\, {\lambda ^n}{{\bf R}}({{\bf \tilde{m}}}_{}^{{\rm{ref}}} - {{{{\bf \tilde{m}}}}^{n - 1}}), \end{eqnarray} (11)where the smoothing operator Wk is no longer needed because the wavelet-based representation of the model can give the constraint inherently. Since we use an L2-norm measure for the data misfit term and an L1-norm measure for the wavelet coefficient, Rd is a diagonal matrix with all elements of 2, while for R, we use the specific form ρ(x) proposed by Ekblom (1987) as   $$\rho (x) = {(x_{}^2 + {\varepsilon ^2})^{\frac{1}{2}}},$$ (12)and we get   $${R_{jj}} = {(x_j^2 + {\varepsilon ^2})^{ - \frac{1}{2}}}, {x_j} = \tilde{m}_j^{{\rm{ref}}} - \tilde{m}_j^{n - 1},$$ (13)where j is the index of the element in the vector x, and ε is a very small real number which is used to guarantee the derivative exists at x = 0. For the large-scale 3-D inversion problem of EM data, the sensitivity matrix is usually calculated by an adjoint forward modeling technique, in which the product of sensitivity matrix and a vector or the product of transpose of sensitivity matrix and a vector is calculated instead of calculating sensitivity matrix itself explicitly. In this situation, we can rewrite formula (11) using eq. (10) and assuming that wavelet basis function is orthogonal (i.e. $${{\bf W}}_w^{ - 1}  =$${{\bf W}}_w^T) as   \begin{eqnarray} &&\left[ {{{\bf W}}_w^{}{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf JW}}_w^{ - 1} + {\lambda ^n}{{\bf R}}} \right]\delta {{\bf \tilde{m}}}\nonumber\\ && = {{\bf W}}_w^{}{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad+\,\, {\lambda ^n}{{\bf R}}({{\bf \tilde{m}}}_{}^{{\rm{ref}}} - {{{{\bf \tilde{m}}}}^{n - 1}}), \end{eqnarray} (14)and therefore the adjoint modeling method can be used to calculate the sensitivity matrix implicitly. The matrices R and the sensitivity matrix J, all depend on the model and they will be updated at each iteration resulting it be solved by IRLS technique. At each iteration, the CG method is used to solve eq. (14). After we obtain the update in the wavelet domain, the model in the space domain is updated as   $${{{\bf m}}^n} = {{{\bf m}}^{n - 1}} + s{{\bf W}}_w^{ - 1}\delta {{\bf \tilde{m}}}.$$ (15) Here, s is a scaling factor which is determined by a linear search proposed by Haber (2014). In practice, a hard threshold to zero out small or nearly zero wavelet coefficients are set in each iteration of the inversion. In the beginning of the inversion, a relatively large threshold value is set to recover the main structures of the model, then it is gradually reduced to obtain more details of the model. 2.3 Selection of the wavelet There are many approaches to do the forward and inverse wavelet transform. In this work, we select the Daubechies (db) wavelet to test the resolution and stability of the wavelet-based inversion because it has a good regularity to guarantee a relatively smooth recovery of the signal. The db wavelet basis functions are orthogonal and biorthogonal, so they can be used directly in eq. (14). For a dbN wavelet, the support width is 2N − 1 with vanishing moment of N and regularity of about 0.2 N for large N. Too large support width will cause problems in wavelet transformation at the ends of the grid sequence, while if the support width is too small, then vanishing moment will be too low to focus the energy. Regularity and vanishing moment are important factors for selecting a wavelet. Regularity defines the smoothness of wavelet function and the vanishing moment defines the flatness and the oscillation of the wavelet function. For dbN wavelet, if we choose a smaller N, then wavelet can represent discontinuous signal better. In this case, we can recover the discontinuous boundary better in the model. While if we choose a larger N, we can get a smoother and more focusing representation of the model, so the inversion will be more stable and has higher resolution. However, there is another problem with larger vanishing moment wavelets. The strong oscillation of the wavelet basis function may produce small fake anomalous regions therefore in this research we mainly test the db2, db6 and db10 wavelets, which scaling and wavelet functions are illustrated in Fig. 1. For the transformation of 2n grids in one direction, we do n − 1 level of transformations to better use the multiscale properties of db wavelets. Figure 1. View largeDownload slide Scaling functions and wavelet functions of db2, db6 and db10 wavelets. Figure 1. View largeDownload slide Scaling functions and wavelet functions of db2, db6 and db10 wavelets. 3 SYNTHETIC EXAMPLES 3.1 Synthetic example 1 for resolution analysis We consider the frequency-domain HEM system of the Geology Survey of Norway (NGU) as an example. This system has total set of five of transmitter–receiver coils working at five different frequencies. Three of them are horizontal coplanar transmitter–receiver coils set at frequencies namely 880, 6606 and 34 133 Hz with corresponding coil separations of 6, 6.3 and 4.9 m, respectively. Two of them are vertical coaxial transmitter–receiver coil set at frequencies 980 and 7001 Hz for coil separations of 6 and 6.3 m, respectively. Terrain clearance is generally maintained around 60 m. The synthetic model 1 is divided in 48 × 48 × 21 cells. A central slice of the model in xz-plane is shown in Fig. 2(a). The cell sizes in the x- and y-directions are 20 m, while in the z-direction the thickness of the top layer is 5 m and the thicknesses of the underlying layers increases with a ratio of 1.2. The background half-space resistivity is 100 Ω⋅m. A flat conductive plate of dimensions 360 m × 360 m × 6 m having a resistivity of 20 Ω⋅m is placed at the centre of the model and at 5 m depth. Another 60 m × 360 m × 80 m inclining plate having a resistivity of 10 Ω⋅m is placed at 18 m depth (see Fig. 2a). Two 10 Ω⋅m anomalies of dimensions 20 m × 360 m × 8 m and 40 m × 360 m × 12 m are located at the right- and left-hand sides of the inclining plate, respectively. The survey stations are assumed to be evenly distributed in the central area of the model with a spacing of 40 m along eight flight lines with line spacing of 100 m to give 19 sites on each line and a total of 152 survey locations. One per cent Gaussian noise is added to the synthetic data to add some noise in it. In the first several iterations of our 3-D inversion, a relatively larger threshold of 90 per cent is set to zero out most of the coefficients to invert the main structures of the model and in later iterations, the threshold is gradually reduced to 0 per cent to recover more detailed features. Figure 2. View largeDownload slide Central slice in xz-plane of (a) synthetic model 1, models constructed by different inversion approaches as (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The circles on the top mark the site locations. Figure 2. View largeDownload slide Central slice in xz-plane of (a) synthetic model 1, models constructed by different inversion approaches as (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The circles on the top mark the site locations. We investigate three different db wavelets namely db2, db6 and db10. Inversion results based on the wavelet-domain method is compared with result based on traditional space-domain method. For the space-domain inversion, we set the initial value of the trade-off parameter λ to be 10 000, while for the wavelet-domain inversion we set the initial value of λ to be 1000. If the data misfit reduces very slowly, we set λ to be 10 per cent of its previous value. After it reduces to 0.1, it will be kept constant in the subsequent iterations. For the L1-norm-type inversion, we set the ε in eq. (13) to be 10−4 and keep it constant throughout the inversion. Inversion results for L1-norm, L2-norm, db2, db6 and db10 are shown in Figs 2(b)–(f) for a slice in the xz-plane passing through centre of the model. From the figures, we can see that the wavelet-based inversions ( Figs 2e and f ) have better resolution for the overlying conductive plate than the Lp-norm method in the space domain. The L1-norm inversion recovers quite focusing result for the two small anomalies at the top, especially for the right one (Fig. 2c). Similar to L1-norm inversion, the db2 wavelet can also generate sharper boundaries because its basis function is not smooth. As can be seen from Fig. 2(d), the two small anomalies recovered by db2 have better resolution than those recovered by other wavelet-based inversions. As shown in Fig. 2(d), db2 is not suitable for complex structures (such as the inclining plate in Fig. 2a) because its wavelet functions are relatively simple. Since db10 basis function has a larger vanishing moment than db2 and db6, so it has more focused energy distribution that lets it generate relatively focusing result for complex structures. The recovered inclining plate in Fig. 2( f ) has a similar focusing as in Fig. 2(c). However, in comparison to L1-norm inversion, the db10 basis function inversion resolves both the small anomalies and large complex bodies. This means that using wavelet-based inversion, we can resolve underground targets of different scales. Also, the regularity of db10 is greater than db2 and db6, which defines the smoothness of the basis function, so the inversion based on db10 is more stable and suitable for complex models. Fig. 3 shows that all the inversion methods can converge to same level of the data misfit (i.e. an rms 1), which means the inversions have worked well. However, for the L1-norm-based inversions, the misfit reduces very slowly at the first 2–3 iterations. This is because in the first several iterations, the diagonal elements Rjj of R calculated by eq. (13) is quite large, so we cannot obtain a large update for the model by solving eq. (11). One possibility would be to increase ε to a large number in the first several iterations to get a large model update, and then change ε to a smaller value in the subsequent iterations. Through numerical testing, we found that we cannot significantly accelerate the convergence of the inversion by using this technique, so we just keep ε small throughout the inversion for simplicity. Figure 3. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 1 as shown in Fig. 1. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 3. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 1 as shown in Fig. 1. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. 3.2 Synthetic example 2 for stability analysis To test the stability of wavelet-based inversion, we use the same overall dimension of the model and HEM system parameters as defined for synthetic example 1, but replace the earlier anomalous bodies with newer anomalous bodies. The background half-space resistivity is still 100 Ω⋅m. Fig. 4(a) shows the synthetic model 2 as layers 1–5 at 0, 5, 11, 18 and 27 m depths, respectively. Upper and lower bodies are 21 and 31 m thick with resistivity of 10 and 20 Ω⋅m, respectively. Figure 4. View largeDownload slide 3-D inversion result with good data coverage for synthetic model 2. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 4. View largeDownload slide 3-D inversion result with good data coverage for synthetic model 2. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. In order to test the stability of different inversion approaches, we calculate two synthetic data set with different line spacing. In the first case, survey stations are evenly distributed in the central area of the model with a spacing of 40 m along eight flight lines having line spacing of 100 m to give 19 sites on each line and a total of 152 survey locations. In the second case, we change the line spacing to 180 m, so there are total five lines and 95 survey locations only. One per cent Gaussian noise is added to both the synthetic data. We compare inversion results of db2, db6 and db10 Daubechies wavelets with inversion results from traditional space-domain methods for both the data set. Same selection strategies for λ and ε is used as done for the synthetic model 1. The inversion results of two different synthetic data sets are shown in Figs 4 and 6, respectively, for layers 1−5. Fig. 4 reveals that when there is good data coverage, then all the five inversion methods recover the true model well and all of them converge to the target data misfit of rms 1 (Fig. 5a). L1-norm inversion could not recover the second body well in layer 5 and recovers lower resistivity of the bodies in all the layers. L2-norm and wavelet-based inversions recovers both bodies well in xy-plane for all the layers except db6 for layer 5. Wavelet-based inversions recovers a closer resistivity to the true resistivity for layers 2–4, while the L2- and L1-norm inversions recover much lower resistivity for the layer 3 than those for layers 2 and 4. The smoothness constraint in space domain always forces the model to change smoothly. But wavelet-domain inversions have constraints in multiscale on the model complexity to let it adapt better. If we have a good data density, then we can recover much sharper boundary in a smoother background through wavelet-domain inversions than the inversion methods based on smoothness constraint in space domain. Figure 5. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with good data coverage as shown in Fig. 3. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 5. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with good data coverage as shown in Fig. 3. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Fig. 6 shows the 3-D inversion result for larger flight line spacing data for same synthetic model 2. Since we use much less data than that in Fig. 4, the inversion is much more ill-posed. Figs 6(b) and (f) show although the resolution is not high in this situation, still these two inversion methods (L2-norm and db10) can recover the model better, which means the db10 wavelet-based inversion is as stable as the L2-norm inversion in the space domain. The other three inversion methods can also get a good data misfit (as shown in Fig. 7a), but the recovered models are different from the true model. From definition of norm, L1-norm prefers to generate a sparser model variation than the L2-norm, so L1-norm is not as stable as the L2-norm for seriously ill-posed inversion problems. When there is enough data as in inversion of synthetic model 1 (Fig. 2), the L1-norm-based inversion can give much sharper boundary than L2-norm-based inversion. Comparing Figs 6(d)–(f) reveals when the number N of dbN increases then the inversion becomes more stable and recovers more features. This is because when N is larger, the regularity of the wavelet is also larger and the basis function is smoother, which makes the recovered model smoother and the inversion more stable. We also observe from Fig. 6(c) that the inversion result of db2 wavelet is simpler with clear discontinuity. Theoretically, if the model is simple and has quite clear discontinuities, the db2 wavelet can give sharper boundaries because its basis function is not smooth. Figure 6. View largeDownload slide 3-D inversion result for synthetic model 2 with poor data coverage. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 6. View largeDownload slide 3-D inversion result for synthetic model 2 with poor data coverage. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 7. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with poor data coverage as shown in Fig. 4. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 7. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with poor data coverage as shown in Fig. 4. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. For better understanding of recovery of the model and sensitivity information in the wavelet-domain inversion, we show wavelet coefficients of the right-hand side of eq. (14) (the minus the gradient of the objective function) in Figs 8(b)–(d) for the first layer (at the surface) of the good data coverage case. Fig. 8(a) shows the minus gradient of objective function in space domain. The distribution in Fig. 8(a) resembles the shape of two anomalous bodies that means if we select a suitable step length then we can use this minus gradient to update the model directly. The wavelet representation of the minus gradient of the objective function transformed by db2, db6 and db10 wavelets are shown in Figs 8(b), (c) and (d), respectively. Since we use three times 1-D wavelet to transform it subsequently in three directions (x, y and z), the wavelet coefficients are mainly focused in one corner of the model. They mainly locate at the bottom left corner. Most of the coefficients are close to zero in Figs 8(b) and (c) for db2 and db6, respectively, except in Fig. 8(d) for db10 coefficients which has a more flat distribution than db2 and db6 because it has a larger vanishing moment. Figure 8. View largeDownload slide Minus gradient of objective function from the first layer of synthetic model 2 for good data coverage case. (a) In space domain, (b) in wavelet domain transformed by db2 wavelet, (c) in wavelet domain transformed by db6 wavelet and (d) in wavelet domain transformed by db10 wavelet. Figure 8. View largeDownload slide Minus gradient of objective function from the first layer of synthetic model 2 for good data coverage case. (a) In space domain, (b) in wavelet domain transformed by db2 wavelet, (c) in wavelet domain transformed by db6 wavelet and (d) in wavelet domain transformed by db10 wavelet. 4 FIELD DATA INVERSION In order to further test the resolution and stability of wavelet-based inversion technique in real world, we inverted a field FDAEM data set. This HEM data set was acquired by NGU in Byneset, Norway. Byneset is located in Trondheim municipality in Mid Norway. Fig. 9 shows the Quaternary geology map of Byneset. Most of the area is covered by marine sediments deposited during the deglaciation around 10 000 yr ago. Due to glacio-isostatic rebound, the sediments are now exposed on surface with a ground surface elevation around 160 m a.s.l. (Reite et al.1999). Central parts of Byneset are covered by ravines and about 100 landslide scars (Solberg et al.2016). A recent landslide in 2012 happened due to presence of the quick clay in the area and it is marked by a green curve in the middle of Fig. 9. The bedrock in the area is dominated by chlorite slate in west and phyllites in the east. The phyllite may contain graphite which is an electronically conductive mineral (Solberg et al.2016). Figure 9. View largeDownload slide Quaternary geological map of Byneset with location of HEM flight lines. Figure 9. View largeDownload slide Quaternary geological map of Byneset with location of HEM flight lines. This data set was collected mostly at 100 m line interval along 60 flight lines shown by black lines in Fig. 9. In-phase and quadrature components of induced secondary magnetic field in part per million (ppm) of transmitted primary field were recorded at 0.1 s interval resulting in 3 m data interval. The data were collected at five different frequencies, three of them in coplanar (34 kHz, 6.6 kHz and 880 Hz) and two of them in co-axial (7 kHz and 980 Hz) setting of transmitter and receiver coils. Details of the processing and spatially constrained inversion of the data are discussed by Baranwal et al. (2015, 2017). For inverting this field data set, we build a model with 231 × 261 × 21 cells. The cell sizes in the x- and y-directions are 30 m, while in the z-direction the thickness of the top layer is 5 m and the thicknesses of the underlying layers increases with a ratio of 1.2. The initial resistivity of background half-space is set to be 100 Ω⋅m. Since the 3-D inversion forward modeling speed is not very fast, we extract total 4148 data from the original data set with a surveying point space of about 50 m on the flight lines. The noise floor is set to be 5 ppm plus 1 per cent of the amplitude. Flight line spacing ranges from about 60 to 300 m (Fig. 9). In order to get a stable inversion, we only test db10 wavelet-based inversion and compare its result with the result of L2-norm-based inversion in the space domain. For both of these two inversion methods, we set the initial value of the trade-off parameter λ to be 1000. Data misfit reduces very slowly so we set λ to be 10 per cent of its previous value. After it reduces to 1, it is kept constant in the subsequent iterations. The inversion result of L2-norm and db10 wavelet methods are illustrated in Fig. 10 for top eight layers (at 0, 5, 11, 18, 27, 37, 50 and 65 m depths). We can see that the two inversion methods give overall similar results which has good agreement with the geology as shown in Fig. 9. The whole surveying area is surrounded by high resistivity outcropping bedrock (observed in the field). The central region of the area is conductive due to thick layers of marine deposits. It comprises a thin and moderately resistive (10–100 Ω⋅m) layer at the top (in layers 1 and 2 of Fig. 10) interpreted as leached clay or silty sediments with a thick unleached marine clay (< 10 Ω⋅m) underneath. Unleached marine sediment is interpreted thickest in the central region of Byneset from 3-D inversion of HEM data. Exposed to deep-seated high resistive bedrock is interpreted at different locations in the area. Upper part of the horizontal depth slices in Fig. 10 ( y > 7 027 500 m) recovered by L2-norm and db10 wavelet look very alike, while lower part ( y < 7 027 500 m) recovered by these two inversions are quite different. There is a good data density to ensure a stable inversion in the upper part where the flight line spacing is smaller. While in the lower part, the line spacing is quite large, therefore it is difficult for smoothing constraint in space domain to bridge the structures from two flight lines because of smaller footprint (about 70 m) of the AEM system in comparison to the line spacing (150–300 m). Figs 10(a)–(h) reveal that the recovered model from L2-norm in the lower part is not continuous. The db10 wavelet-based inversion (Figs 10i–p) gives a smoother result in the lower part because inherent constraint of wavelet inversion can recover the main structure correctly. Inversion result around y = 7 028 000 m is not reliable from both the inversions because the data have large noise. Figure 10. View largeDownload slide 3-D inversion result of field data acquired in Byneset, Norway. (a)–(h) are the horizontal depth slices of top eight layers of L2-norm-based 3-D inversion. (i) – (p) are the horizontal depth slices of top eight layers of db10 wavelet-based 3-D inversion. Layers–8 correspond to the depth at 0, 5, 11, 18, 27, 37, 50 and 65 m, respectively. Figure 10. View largeDownload slide 3-D inversion result of field data acquired in Byneset, Norway. (a)–(h) are the horizontal depth slices of top eight layers of L2-norm-based 3-D inversion. (i) – (p) are the horizontal depth slices of top eight layers of db10 wavelet-based 3-D inversion. Layers–8 correspond to the depth at 0, 5, 11, 18, 27, 37, 50 and 65 m, respectively. Plots of the inversion parameters during iterations are shown in Fig. 11. L2-norm-based inversion converges a little faster than the db10 wavelet-based inversion. Because we do not provide an accurate noise floor to the data, therefore the data misfit (rms) of these two methods reduces from initial 9.52 to about 2 only. However, we observe a good fit between field data and predicted data visually from both the inversions. As illustrated in Fig. 11, we have reset the trade-off parameter λ to be 100 after four iterations for db10 inversion. This is done because trade-off parameter λ is allowed to be smaller to get a suitable model update in the beginning of the inversion and later it is reset to a large λ again to ensure a stable inversion. Figure 11. View largeDownload slide Plots of inversion parameters obtained from the inversion for the field data from Byneset as shown in Fig. 9. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 11. View largeDownload slide Plots of inversion parameters obtained from the inversion for the field data from Byneset as shown in Fig. 9. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. To show the vertical resolution of the inverted models obtained from these two inversions of the field data, we illustrate vertical sections at y = 7 030 191 m from both the inversions in Fig. 12. This vertical section passes through 2012 landslide location. Two vertical slices obtained from L2-norm and db10 inversions have similar layered structure. The 2012 landslide took place at around x = 557 000 m, where we see a clear discontinuity of the conductive layers (unleached clay) in both the images. Presence of the quick clay at this location is confirmed by soil samples analysis (Solberg et al.2016). The result shown in Fig. 12 has a good agreement with HEM and ERT inversion results shown in Solberg et al. (2016). The inversion result of db10 wavelet (Fig. 12b) provides a clearer image of unleached clay layer than the L2-norm inversion (Fig. 12a). Due to the smoothing constraint in the space-domain method, the result of L2-norm inversion connects the surface and middle thick conductive layer of unleached clay at several locations. The field data and predicted data illustrated in Figs 12(c) and (d) show that the two inversion methods fit the data to a similar level. Figure 12. View largeDownload slide Vertical sections of inversion result at y = 7 030 191 m from (a) L2-norm-based method and (b) db10 wavelet-based method; (c) and (d) are data fitting for in-phase and quadrature components. Figure 12. View largeDownload slide Vertical sections of inversion result at y = 7 030 191 m from (a) L2-norm-based method and (b) db10 wavelet-based method; (c) and (d) are data fitting for in-phase and quadrature components. In Fig. 13, we illustrate the first (at 0 m depth) and fourth layer (at 18 m depth) inversion results at different iterations to show how the two inversion methods work. The two methods have different convergence process, so it is not suitable to compare the result from the same iteration, therefore we compare the results from the iterations with closer rms values. Since the AEM data are dense enough and one site relates to several nearby sites. The smoothing constraint in L2-norm-based inversion recovers the local structures quickly through the inversion (as shown in Figs 13a–d for the first layer and in Figs 13e–h for the fourth layer). In the wavelet-based inversion, we use a large lambda in the beginning, which enforces a very strong sparsity constraint and only allows few large-scale coefficients, so the inversion recovers the main structure of this model first. Figs 13(i)–(l) (for first layer) and Figs 13(m)–(p) (for fourth layer) show this property of the wavelet inversion by recovering the main structure of the entire model without any local detail in the first several iterations (Figs 13j and k and n and o). More local details appear in later iterations of the inversion (Figs 13l and p). Therefore, the wavelet-based inversion can recover structures of the model at different scales. Figure 13. View largeDownload slide Inversion results of first and fourth layers from two inversion approaches in different iterations. (a)–(d) depict first layer at the surface obtained from L2-norm inversion at rms values 8.51, 5.38, 4.06 and 1.90. (e)–(h) depict fourth layer at 18 m depth obtained from L2-norm-based inversion at rms values 8.51, 5.38, 4.06 and 1.90. (i)–(l) depict first layer at the surface obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. (m)–(p) depict fourth layer at 18 m depth obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. Figure 13. View largeDownload slide Inversion results of first and fourth layers from two inversion approaches in different iterations. (a)–(d) depict first layer at the surface obtained from L2-norm inversion at rms values 8.51, 5.38, 4.06 and 1.90. (e)–(h) depict fourth layer at 18 m depth obtained from L2-norm-based inversion at rms values 8.51, 5.38, 4.06 and 1.90. (i)–(l) depict first layer at the surface obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. (m)–(p) depict fourth layer at 18 m depth obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. 5 CONCLUSIONS We propose a new stable and well-resolved 3-D inversion method based on the wavelet transform and sparsity constraint to invert FDAEM data. This new method takes advantages of multiscale property of the wavelet representation which can be used to recover different scale structure of the model in different iterations. The wavelet-based inversion gives constraints in multiscale on the model complexity, while the smooth constraint in space domain is affected more by nearby cells. Synthetic and field data inversion tests show that convergence is stable and the inversion can produce as good resolution as L1-norm approaches in the space domain for sufficiently dense data. Both the wavelet-based inversion and L1-norm inversion give more focusing result than L2-norm inversion. The sparsity constraint on model complexity in wavelet domain can overcome the problem of sparse data distribution which always results in discontinuous structure, and generates a data-based resolution that adaptive to the data distribution rather than the model-based resolution that depends on the model discretization. The method has the potential to improve inversion resolution for many other EM methods without compromising the stability of the inversion. Acknowledgements We want to thank the reviewers and editors for their suggestions that helped to clarify this paper. This research is funded by the Key National Research Project of China (2016YFC0303100, 2017YFC0601900). This research is also funded by Key Program of National Natural Science Foundation of China (41530320) and Natural Science Foundation for Young Scientist (41404093). NGU is thanked for making the Byneset HEM data available. REFERENCES Abers G.A., Roecker S.W., 1991. Deep structure of an arc-continent collision: earthquake relocation and inversion for upper mantle P and S wave velocities beneath Papua New Guinea. J. Geophys. Res. , 96( B4), 6379– 6401. https://doi.org/10.1029/91JB00145 Aster R.C., Borchers B., Thurber C.H., 2013. Parameter Estimation and Inverse Problems , Academic Press. Baranwal V.C., Dalsegg E., Tønnesen J.F., Rønning J.S., Solberg I.L., Rodionov A., Dretvik H., 2015. Mapping of marine clay layers using airborne EM and ground geophysical surveys at Byneset, Trondheim municipality. Geological Survey of Norway (NGU) Report 2015.006 . Available at: https://www.ngu.no/en/publikasjon/mapping-marine-clay-layers-using-airborne-em-and-ground-geophysical-methods-byneset. Baranwal V.C., Rønning J.S., Solberg I.-L., Dalsegg E., Tønnesen J.F., Long M., 2017, Investigation of a sesnsitive clay landslide area using frequency-domain helicopter-borne EM and ground geophysical methods, in Landslides in Sensitive Clays, Series in Advances in Natural and Technological Hazards Research , eds Thakur V., L’Heureux J.-S., Locat A., p. 46. Springer International Publishing. Beamish D., 2003. Airborne EM footprints, Geophys. Prospect. , 51( 1), 49– 60. https://doi.org/10.1046/j.1365-2478.2003.00353.x Chiao L.Y., Kuo B.Y., 2001. Multiscale seismic tomography, Geophys. J. Int. , 145( 2), 517– 527. https://doi.org/10.1046/j.0956-540x.2001.01403.x Chiao L.Y., Liang W.T., 2003. Multiresolution parameterization for geophysical inverse problems, Geophysics , 68( 1), 199– 209. https://doi.org/10.1190/1.1543207 Cox L.H., Wilson G.A., Zhdanov M.S., 2010. 3D inversion of airborne electromagnetic data using a moving footprint, Explor. Geophys. , 41( 4), 250– 259. https://doi.org/10.1071/EG10003 Daubechies I., 1992. Ten Lectures on Wavelets , Society for Industrial and Applied Mathematics. Daubechies I., Defrise M., De Mol C., 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Comm. Pure appl. Math. , 57( 11), 1413– 1457. https://doi.org/10.1002/cpa.20042 Davis K., Li Y., 2011. Fast solution of geophysical inversion using adaptive mesh, space-filling curves and wavelet compression, Geophys. J. Int. , 185( 1), 157– 166. https://doi.org/10.1111/j.1365-246X.2011.04929.x Delost M., Virieux J., Operto S., 2008. First-arrival traveltime tomography using second generation wavelets, Geophys Prospect , 56( 4), 505– 526. https://doi.org/10.1111/j.1365-2478.2008.00710.x Donoho D.L., 2006. For most large underdetermined systems of linear equations the minimal?1-norm solution is also the sparsest solution, Comm. Pure appl. Math. , 59( 6), 797– 829. https://doi.org/10.1002/cpa.20132 Ekblom H., 1987, The L1-estimate as limiting case of an Lp-or Huber-estimate, in Statistical Data Analysis based on the L1-Norm and Related Methods , ed. Doge Y., pp. 109– 116. Elsevier Science Publ. Co., Inc. Fang H., Zhang H., 2014. Wavelet-based double-difference seismic tomography with sparsity regularization, Geophys. J. Int. , 199( 2), 944– 955. https://doi.org/10.1093/gji/ggu305 Farquharson C.G., 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions, Geophysics , 73( 1), K1– K9. https://doi.org/10.1190/1.2816650 Haber E., 2014, Computational Methods in Geophysical Electromagnetics , SIAM. Hung S.H., Chen W.P., Chiao L.Y., 2011. A data‐adaptive, multiscale approach of finite‐frequency, traveltime tomography with special reference to P and S wave data from central Tibet, J. geophys. Res. , 116, B06307, doi:10.1029/2010JB008190. Hung S.H., Chen W.P., Chiao L.Y., Tseng T.L., 2010. First multi‐scale, finite‐frequency tomography illuminates 3‐D anatomy of the Tibetan plateau, Geophys. Res. Lett. , 41, 7025– 7034, doi:10.1002/2014GL061644. Li Y., Oldenburg D.W., 1999. Rapid construction of equivalent sources using wavelets, in 69th Annual International Meeting, Expanded Abstracts , Vol. 18, pp. 374– 378. SEG, Houston TX. Li Y., Oldenburg D.W., 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method, Geophys. J. Int. , 152( 2), 251– 265. https://doi.org/10.1046/j.1365-246X.2003.01766.x Li Y., Oldenburg D.W., 2010. Rapid construction of equivalent sources using wavelets, Geophysics , 75( 3), L51– L59. https://doi.org/10.1190/1.3378764 Liu Y.H., Yin C.C., 2013. 3D inversion for frequency-domain HEM data, Chin. J. Geophys. , 56, 4278– 4287. Liu Y.H., Yin C.C., 2016. 3D inversion for multipulse airborne transient electromagnetic data, Geophysics , 81( 6), E401– E408. https://doi.org/10.1190/geo2015-0481.1 Loris I., Nolet G., Daubechies I., Dahlen F., 2007. Tomographic inversion using? 1 -norm regularization of wavelet coefficients, Geophys. J. Int. , 170( 1), 359– 370. https://doi.org/10.1111/j.1365-246X.2007.03409.x Nittinger C.G., Becken M., 2016. Inversion of magnetotelluric data in a sparse model domain, Geophys. J. Int. , 206( 2), 1398– 1409. https://doi.org/10.1093/gji/ggw222 Oldenburg D., McGillivray P., Ellis R., 1993. Generalized subspace methods for large-scale inverse problems, Geophys. J. Int. , 114( 1), 12– 20. https://doi.org/10.1111/j.1365-246X.1993.tb01462.x Portniaguine O., Zhdanov M.S., 1999. Focusing geophysical inversion images, Geophysics , 64( 3), 874– 887. https://doi.org/10.1190/1.1444596 Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1996. Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing , Cambridge University Press. Reite A.J., Sveian H., Erichsen E., 1999. Trondheim fra istid til nåtid – landskapshistorie og løsmasser , Geological Survey of Norway (NGU), Gråsteinen 5 (in Norwegian). Available at: https://www.ngu.no/publikasjon/gr%C3%A5steinen-5-trondheim-fra-istid-til-n%C3%A5tid. Simons F.J.et al.  , 2011. Solving or resolving global tomographic models with spherical wavelets, and the scale and sparsity of seismic heterogeneity, Geophys. J. Int. , 187( 2), 969– 988. https://doi.org/10.1111/j.1365-246X.2011.05190.x Siripunvaraporn W., Egbert G., Lenbury Y., Uyeshima M., 2005. Three-dimensional magnetotelluric inversion: data-space method, Phys. Earth planet. Inter. , 150( 1–3), 3– 14. https://doi.org/10.1016/j.pepi.2004.08.023 Solberg I.L., Long M., Baranwal V.C., Gylland A.S., Rønning J.S., 2016. Geophysical and geotechnical studies of geology and sediment properties at a quick-clay landslide site at Esp, Trondheim, Norway, Eng. Geol. , 208, 214– 230. https://doi.org/10.1016/j.enggeo.2016.04.031 Spakman W., Bijwaard H., 2001. Optimization of cell parameterizations for tomographic inverse problems, Pure appl. geophys. , 158( 8), 1401– 1423. https://doi.org/10.1007/PL00001227 Zhang H., Thurber C., 2005. Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: application to Parkfield, California, J. geophys. Res. , 110( B4), B04303, doi:10.1029/2004JB003186. https://doi.org/10.1029/2004JB003186 © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Wavelet-based 3-D inversion for frequency-domain airborne EM data

, Volume 213 (1) – Apr 1, 2018
15 pages

/lp/ou_press/wavelet-based-3-d-inversion-for-frequency-domain-airborne-em-data-jFS4Ae22dD
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx545
Publisher site
See Article on Publisher Site

### Abstract

Summary In this paper, we propose a new wavelet-based 3-D inversion method for frequency-domain airborne electromagnetic (FDAEM) data. Instead of inverting the model in the space domain using a smoothing constraint, this new method recovers the model in the wavelet domain based on a sparsity constraint. In the wavelet domain, the model is represented by two types of coefficients, which contain both large- and fine-scale informations of the model, meaning the wavelet-domain inversion has inherent multiresolution. In order to accomplish a sparsity constraint, we minimize an L1-norm measure in the wavelet domain that mostly gives a sparse solution. The final inversion system is solved by an iteratively reweighted least-squares method. We investigate different orders of Daubechies wavelets to accomplish our inversion algorithm, and test them on synthetic frequency-domain AEM data set. The results show that higher order wavelets having larger vanishing moments and regularity can deliver a more stable inversion process and give better local resolution, while the lower order wavelets are simpler and less smooth, and thus capable of recovering sharp discontinuities if the model is simple. At last, we test this new inversion algorithm on a frequency-domain helicopter EM (HEM) field data set acquired in Byneset, Norway. Wavelet-based 3-D inversion of HEM data is compared to L2-norm-based 3-D inversion's result to further investigate the features of the new method. Controlled source electromagnetics (CSEM), Non-linear electromagnetic, Inverse theory, Wavelet transform 1 INTRODUCTION For 3-D electromagnetic (EM) inversion, the geo-electrical model is generally parametrized into regular or unstructured cells. Size of the model are usually not small, so the number of unknowns is much larger than the number of observed data, which always results in an ill-conditioned problem to be solved. In order to make the 3-D inversion stable, some kind of regularization such as Tikhonov regularization needs to be applied (Aster et al.2013). Such regularization methods usually use a smoothing operator to build a relationship between different cells. This technique indeed makes the inversion more stable, but if the regularization is strong, it often smears out small-scale features of the model even where there is sufficient data density available to provide more detailed information. If we use a weak regularization, the resolution may be improved, but the inversion may not be stable. Therefore, it is not easy to get a stable and high-resolution inversion at the same time. For some inversion problems with different site spacing/line spacing, the smoothing regularization may not give a reasonable geo-electrical structure because of the uneven distribution of data. Several techniques such as focusing inversion and L1-norm-based inversion were proposed to improve the resolution of stable gravity and 2-D magnetotelluric (MT) inversion (Portniaguine & Zhdanov 1999; Farquharson 2008). These techniques give the same constraint to the whole inverted model without considering the uneven data distribution. Actually, for the uneven data distribution, the model should be recovered in multiresolution. However, for the model parametrized by regular cell, it is hard to make the model to be adaptive with the uneven data distribution. In seismic, to get different multiscale features of the model according to the ray distribution, several adaptive seismic tomography methods have been proposed in which the cells are irregularly distributed to match with the uneven ray distribution (Abers & Roecker 1991; Spakman & Bijwaard 2001; Zhang & Thurber 2005). These methods work well for recovering multiscale features of the model, but they are based purely on data and do not consider whether the model is appropriately parametrized. The wavelet representation of a signal has an inherent multiscale nature (Daubechies 1992). The wavelet-based multiscale tomography methods are successfully used to recover the multiscale features of the model parametrized by regular cells in seismic applications (Chiao & Kuo 2001; Chiao & Liang 2003; Loris et al.2007; Delost et al.2008; Hung et al.2011; Simons et al.2011; Fang & Zhang 2014). Nittinger & Becken (2016) have also used the sparsity constraint in wavelet domain to accomplish 2-D MT inversion. In fact, before the wavelet is used in seismic tomography for multiscale inversion, it was used in gravity data inversion by Li & Oldenburg (1999) to compress the sensitivity matrix. Later, they made some improvements in wavelet compression-based gravity and magnetic inversions as well (Li & Oldenburg 2003, 2010; Davis & Li 2011). Actually in their research, they already recovered the model in the wavelet domain, but they focused only on the matrix compress and acceleration of the inversion. Besides these applications in linear inversion, Chiao & Kuo (2001) and Chiao & Liang (2003) were the first to propose wavelet-based multiscale inversion for the non-linear tomography inversion problem. In their research, they used the wavelet representation as the model constraints instead of using the common smoothness regularization. The wavelet representation of the model can directly give multiresolution inversion results. For the smoother part of the model, few wavelet coefficients at large scales are significant, while for the part of model with small-scale features and with good data coverage, more wavelet coefficients will be used to sufficiently represent the model. Daubechies et al. (2004) proposed an iterative thresholding method for linear inverse problems using a sparsity constraint. Hung et al. (2010, 2011) proposed a wavelet-based finite-frequency teleseismic tomography method to recover the velocity structure of the central Tibet. Delost et al. (2008) used the wavelet transform for the first-arrival traveltime tomography by designing a bit mask operator to only invert for the wavelet coefficients in specific areas where the resolution is high. Based on work of Daubechies et al. (2004), Loris et al. (2007) and Simons et al. (2011) used the L1-norm of model wavelet coefficients in their inversion scheme to make the solution sparse, which is different from the L2-norm of the model wavelet coefficients used by Chiao & Kuo (2001), Chiao & Liang (2003) and Hung et al. (2010, 2011). Following strategy of Loris et al. (2007) and Nittinger & Becken (2016) accomplished 2-D MT inversion in a sparse model. For the L1-norm inversion of wavelet coefficients, Donoho (2006) proved that it mostly gave a sparse solution and the inversion is a sparsity constraint. Aster et al. (2013) stated that if the model was sparse in nature, the sparsity-constraint inversion could recover it better without a priori smoothing regularization term. Besides the wavelet-based inversion, there are also some other subspace inversion methods. Oldenburg et al. (1993) proposed generalized subspace methods for large-scale inverse problems using basis vectors representation to accelerate the convergence rate. Siripunvaraporn et al. (2005) inverted the MT data using a data-space method which could reduce the inversion space efficiently in comparison to a traditional model-space method. Though these two methods do not have any concept of the multiresolution inversion. In this research, we want to test the wavelet-based multiscale inversion on the EM data that have lower resolution than the seismic data and its sparsity-constraint inversion is less studied. Different from the previous work in seismic, we propose a new wavelet-based inversion method for EM data that uses the model constraints on three directions, which is quite similar to the common smoothing constraint in space domain. We want to further analyse the difference between wavelet-based inversion and traditional Lp-norm inversion method from aspects of both resolution and stability. This new method is derived from a general measure inversion algorithm, and the final equation for this iteratively reweighted least-squares (IRLS) problem is a Gauss–Newton equation solved by conjugate gradient (CG) method. We solve the final equation directly using an IRLS without using the soft threshold method as proposed by Daubechies et al. (2004). We use helicopter EM (HEM) data as an example to test the character of this new inversion method which will be valid in inversion of any frequency-domain airborne EM (FDAEM) data. The AEM survey most of the time has an uneven data distribution. Since the AEM system is compact and has a limited ‘footprint’ (Beamish 2003) for rather larger survey-line spacing, it is hard for conventional 3-D inversion to connect the conductivity structures from two consecutive flight lines. The wavelet-based inversion treats the whole model as a 3-D volume, which main features are recovered. This property can bridge the conductivity structures from two survey lines naturally. For small line spacing, the AEM survey can provide denser data distribution, which can recover fine local structures as well. Based on Liu & Yin's (2013) work, we use the finite-difference method to solve the forward modeling problem and test different kinds of the wavelet basis functions to compare it with other traditional methods on both synthetic and field HEM data set. We further demonstrate the impact of different order of the wavelet basis functions on the inversion result. 2 METHODOLOGY 2.1 General inversion approach Generally, the objective function for an inversion problem is expressed as (Farquharson 2008)   $${\bf \Phi } = {\phi _d} + \lambda {\phi _m}.$$ (1)Here, ϕd and ϕm are the measures of data misfit and model complexity, respectively. The factor λ is the trade-off parameter which balances the contributions of the data misfit and the model complexity term. ϕd can be expressed in a general form as   $${\phi _d} = {\phi _d}({\bf u}),$$ (2)and   $${\bf u} = {{\bf W}_d}({{\bf d}^{\rm obs}} - {{\bf d}^{\rm prd}}),$$ (3)where dobs is the vector of observed data, dprd is the predicted data computed using the model m, and Wd is a diagonal matrix whose elements are the reciprocals of the estimates of the standard deviations of the noise in the observations. ϕm can have the form as   $${\phi _m} = \sum\limits_k {{\alpha _k}} {\phi _k}({{\bf v}_k}),$$ (4)and   $${{\bf v}_k} = {{\bf W}_k}({{\bf m}} - {\bf m}_k^{{\rm{ref}}}).$$ (5)Here, m is the vector of the recovered model parameters, mkref is the reference model, Wk is the smoothing operator matrix which is used to constrain the complexity of the recovered model and αk are the coefficients used to balance the contributions of different kinds of measure of the recovered model. A general form of ϕd and ϕk is   $$\phi ({{\bf x}}) = \sum\limits_j {\rho ({x_j})} ,$$ (6)where xj are the elements of the vector x, which can be either u or vk. There are numerous choices for the specific form of ρ(xj), for example, L2-norm, L1-norm and so on. In the traditional inversion approach, differentiating both sides of eq. (1) with respect to the space-domain model parameters and neglecting the higher order terms of the Taylor-series expansion gives the linear system of equations to be solved at each iteration as   \begin{eqnarray} &&\left[{{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf J}} + {\lambda ^n}\sum\limits_k {{\alpha _k}{{\bf W}}_k^T{{{\bf R}}_k}{{{\bf W}}_k}} } \right]\delta {{\bf m}}\nonumber\\ &&= {{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad +\,\, {\lambda ^n}\sum\limits_k {{\alpha _k}{{\bf W}}_k^T{{{\bf R}}_k}{{{\bf W}}_k}} ({{\bf m}}_k^{{\rm{ref}}} - {{{\bf m}}^{n - 1}}), \end{eqnarray} (7)where the matrix R is a diagonal matrix which has the form as   $${{\bf R}} = {\rm diag}\left\{ {\begin{array}{*{20}{c}} {{{\rho '({x_j})} / {{x_j},}}}&{ \cdots \cdots ,}&{{{\rho '({x_N})} / {{x_N},}}} \end{array}} \right\}.$$ (8) J is the sensitivity matrix in the space-domain model parameters and it can be calculated by an adjoint modeling technique (Liu & Yin 2016), dn−1 is the vector of the data calculated from the model mn−1 which is obtained from the previous iteration and δm is model update. We apply the ‘moving footprint’ method of Cox et al. (2010) to accelerate the forward modeling. 2.2 Wavelet-domain inversion with sparsity constraint In the smoothness constrained inversion, we usually use an L2-norm measure for the model's complexity to get a smooth-edged model or an L1-norm to get a sharper edged model or some other focusing measure to get a somewhat ‘focused’ model with greater boundary resolution. In numerical tests, an L2-norm inversion is typically more stable than an L1-norm inversion, particularly if data density is low. It is hard to combine the advantages of L2- and L1-norms in one single-inversion scheme. The wavelet transform can divide one signal into two groups of coefficients: one is the approximate coefficients and the other is the detail coefficients. This means that this transform has the potential to resolve the model at different levels at the same time according to the data set distribution. This property of wavelet transform motivates us to perform the inversion in the wavelet domain. A space-domain model m can be transformed into wavelet-domain model coefficients $${{\bf \tilde{m}}}$$ using a wavelet transform matrix Ww as   $${{\bf \tilde{m}}} = {{{\bf W}}_w}{{\bf m}}.$$ (9) Since the wavelet transform usually demands that the dimension of the vector should be a power of two, we first extend the three dimensions of m to be power of two, then we apply the 1-D partial wavelet transform scheme of Press et al. (1996) to each dimension of the 3-D model sequentially. This wavelet transform algorithm is a filtering process and the filtering coefficients can be found in Daubechies (1992). According to the chain rule, the sensitivity matrix in the wavelet domain is computed as   $${{\bf \tilde{J}}} = {{\bf JW}}_w^{{\rm{ - 1}}},$$ (10)where $${{\bf W}}_w^{ - 1}$$ is the inverse wavelet transform operator. Since the wavelet coefficients are usually sparse, we use a sparsity constraint in the wavelet domain instead of the smoothing constraint in the space domain (Daubechies et al.2004). According to Donoho (2006), the L1-norm measure can lead to a sparse solution, so we simply use the L1-norm measure on the wavelet coefficients and set k = 1 and α1 = 1 in eq. (4). For the data misfit term, we still use the L2-norm. Then, the linear system (eq. 7) to be solved at each iteration of the inversion becomes   \begin{eqnarray} &&\left[ {{{{{\bf \tilde{J}}}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf \tilde{J}}} + {\lambda ^n}{{\bf R}}} \right]\delta {{\bf \tilde{m}}}\nonumber\\ &&= {{{{\bf \tilde{J}}}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad+\,\, {\lambda ^n}{{\bf R}}({{\bf \tilde{m}}}_{}^{{\rm{ref}}} - {{{{\bf \tilde{m}}}}^{n - 1}}), \end{eqnarray} (11)where the smoothing operator Wk is no longer needed because the wavelet-based representation of the model can give the constraint inherently. Since we use an L2-norm measure for the data misfit term and an L1-norm measure for the wavelet coefficient, Rd is a diagonal matrix with all elements of 2, while for R, we use the specific form ρ(x) proposed by Ekblom (1987) as   $$\rho (x) = {(x_{}^2 + {\varepsilon ^2})^{\frac{1}{2}}},$$ (12)and we get   $${R_{jj}} = {(x_j^2 + {\varepsilon ^2})^{ - \frac{1}{2}}}, {x_j} = \tilde{m}_j^{{\rm{ref}}} - \tilde{m}_j^{n - 1},$$ (13)where j is the index of the element in the vector x, and ε is a very small real number which is used to guarantee the derivative exists at x = 0. For the large-scale 3-D inversion problem of EM data, the sensitivity matrix is usually calculated by an adjoint forward modeling technique, in which the product of sensitivity matrix and a vector or the product of transpose of sensitivity matrix and a vector is calculated instead of calculating sensitivity matrix itself explicitly. In this situation, we can rewrite formula (11) using eq. (10) and assuming that wavelet basis function is orthogonal (i.e. $${{\bf W}}_w^{ - 1}  =$${{\bf W}}_w^T) as   \begin{eqnarray} &&\left[ {{{\bf W}}_w^{}{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}{{\bf JW}}_w^{ - 1} + {\lambda ^n}{{\bf R}}} \right]\delta {{\bf \tilde{m}}}\nonumber\\ && = {{\bf W}}_w^{}{{{\bf J}}^T}{{\bf W}}_d^T{{{\bf R}}_d}{{{\bf W}}_d}({{{\bf d}}^{{\rm{obs}}}} - {{{\bf d}}^{n - 1}})\nonumber\\ &&\quad+\,\, {\lambda ^n}{{\bf R}}({{\bf \tilde{m}}}_{}^{{\rm{ref}}} - {{{{\bf \tilde{m}}}}^{n - 1}}), \end{eqnarray} (14)and therefore the adjoint modeling method can be used to calculate the sensitivity matrix implicitly. The matrices R and the sensitivity matrix J, all depend on the model and they will be updated at each iteration resulting it be solved by IRLS technique. At each iteration, the CG method is used to solve eq. (14). After we obtain the update in the wavelet domain, the model in the space domain is updated as   $${{{\bf m}}^n} = {{{\bf m}}^{n - 1}} + s{{\bf W}}_w^{ - 1}\delta {{\bf \tilde{m}}}.$$ (15) Here, s is a scaling factor which is determined by a linear search proposed by Haber (2014). In practice, a hard threshold to zero out small or nearly zero wavelet coefficients are set in each iteration of the inversion. In the beginning of the inversion, a relatively large threshold value is set to recover the main structures of the model, then it is gradually reduced to obtain more details of the model. 2.3 Selection of the wavelet There are many approaches to do the forward and inverse wavelet transform. In this work, we select the Daubechies (db) wavelet to test the resolution and stability of the wavelet-based inversion because it has a good regularity to guarantee a relatively smooth recovery of the signal. The db wavelet basis functions are orthogonal and biorthogonal, so they can be used directly in eq. (14). For a dbN wavelet, the support width is 2N − 1 with vanishing moment of N and regularity of about 0.2 N for large N. Too large support width will cause problems in wavelet transformation at the ends of the grid sequence, while if the support width is too small, then vanishing moment will be too low to focus the energy. Regularity and vanishing moment are important factors for selecting a wavelet. Regularity defines the smoothness of wavelet function and the vanishing moment defines the flatness and the oscillation of the wavelet function. For dbN wavelet, if we choose a smaller N, then wavelet can represent discontinuous signal better. In this case, we can recover the discontinuous boundary better in the model. While if we choose a larger N, we can get a smoother and more focusing representation of the model, so the inversion will be more stable and has higher resolution. However, there is another problem with larger vanishing moment wavelets. The strong oscillation of the wavelet basis function may produce small fake anomalous regions therefore in this research we mainly test the db2, db6 and db10 wavelets, which scaling and wavelet functions are illustrated in Fig. 1. For the transformation of 2n grids in one direction, we do n − 1 level of transformations to better use the multiscale properties of db wavelets. Figure 1. View largeDownload slide Scaling functions and wavelet functions of db2, db6 and db10 wavelets. Figure 1. View largeDownload slide Scaling functions and wavelet functions of db2, db6 and db10 wavelets. 3 SYNTHETIC EXAMPLES 3.1 Synthetic example 1 for resolution analysis We consider the frequency-domain HEM system of the Geology Survey of Norway (NGU) as an example. This system has total set of five of transmitter–receiver coils working at five different frequencies. Three of them are horizontal coplanar transmitter–receiver coils set at frequencies namely 880, 6606 and 34 133 Hz with corresponding coil separations of 6, 6.3 and 4.9 m, respectively. Two of them are vertical coaxial transmitter–receiver coil set at frequencies 980 and 7001 Hz for coil separations of 6 and 6.3 m, respectively. Terrain clearance is generally maintained around 60 m. The synthetic model 1 is divided in 48 × 48 × 21 cells. A central slice of the model in xz-plane is shown in Fig. 2(a). The cell sizes in the x- and y-directions are 20 m, while in the z-direction the thickness of the top layer is 5 m and the thicknesses of the underlying layers increases with a ratio of 1.2. The background half-space resistivity is 100 Ω⋅m. A flat conductive plate of dimensions 360 m × 360 m × 6 m having a resistivity of 20 Ω⋅m is placed at the centre of the model and at 5 m depth. Another 60 m × 360 m × 80 m inclining plate having a resistivity of 10 Ω⋅m is placed at 18 m depth (see Fig. 2a). Two 10 Ω⋅m anomalies of dimensions 20 m × 360 m × 8 m and 40 m × 360 m × 12 m are located at the right- and left-hand sides of the inclining plate, respectively. The survey stations are assumed to be evenly distributed in the central area of the model with a spacing of 40 m along eight flight lines with line spacing of 100 m to give 19 sites on each line and a total of 152 survey locations. One per cent Gaussian noise is added to the synthetic data to add some noise in it. In the first several iterations of our 3-D inversion, a relatively larger threshold of 90 per cent is set to zero out most of the coefficients to invert the main structures of the model and in later iterations, the threshold is gradually reduced to 0 per cent to recover more detailed features. Figure 2. View largeDownload slide Central slice in xz-plane of (a) synthetic model 1, models constructed by different inversion approaches as (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The circles on the top mark the site locations. Figure 2. View largeDownload slide Central slice in xz-plane of (a) synthetic model 1, models constructed by different inversion approaches as (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The circles on the top mark the site locations. We investigate three different db wavelets namely db2, db6 and db10. Inversion results based on the wavelet-domain method is compared with result based on traditional space-domain method. For the space-domain inversion, we set the initial value of the trade-off parameter λ to be 10 000, while for the wavelet-domain inversion we set the initial value of λ to be 1000. If the data misfit reduces very slowly, we set λ to be 10 per cent of its previous value. After it reduces to 0.1, it will be kept constant in the subsequent iterations. For the L1-norm-type inversion, we set the ε in eq. (13) to be 10−4 and keep it constant throughout the inversion. Inversion results for L1-norm, L2-norm, db2, db6 and db10 are shown in Figs 2(b)–(f) for a slice in the xz-plane passing through centre of the model. From the figures, we can see that the wavelet-based inversions ( Figs 2e and f ) have better resolution for the overlying conductive plate than the Lp-norm method in the space domain. The L1-norm inversion recovers quite focusing result for the two small anomalies at the top, especially for the right one (Fig. 2c). Similar to L1-norm inversion, the db2 wavelet can also generate sharper boundaries because its basis function is not smooth. As can be seen from Fig. 2(d), the two small anomalies recovered by db2 have better resolution than those recovered by other wavelet-based inversions. As shown in Fig. 2(d), db2 is not suitable for complex structures (such as the inclining plate in Fig. 2a) because its wavelet functions are relatively simple. Since db10 basis function has a larger vanishing moment than db2 and db6, so it has more focused energy distribution that lets it generate relatively focusing result for complex structures. The recovered inclining plate in Fig. 2( f ) has a similar focusing as in Fig. 2(c). However, in comparison to L1-norm inversion, the db10 basis function inversion resolves both the small anomalies and large complex bodies. This means that using wavelet-based inversion, we can resolve underground targets of different scales. Also, the regularity of db10 is greater than db2 and db6, which defines the smoothness of the basis function, so the inversion based on db10 is more stable and suitable for complex models. Fig. 3 shows that all the inversion methods can converge to same level of the data misfit (i.e. an rms 1), which means the inversions have worked well. However, for the L1-norm-based inversions, the misfit reduces very slowly at the first 2–3 iterations. This is because in the first several iterations, the diagonal elements Rjj of R calculated by eq. (13) is quite large, so we cannot obtain a large update for the model by solving eq. (11). One possibility would be to increase ε to a large number in the first several iterations to get a large model update, and then change ε to a smaller value in the subsequent iterations. Through numerical testing, we found that we cannot significantly accelerate the convergence of the inversion by using this technique, so we just keep ε small throughout the inversion for simplicity. Figure 3. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 1 as shown in Fig. 1. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 3. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 1 as shown in Fig. 1. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. 3.2 Synthetic example 2 for stability analysis To test the stability of wavelet-based inversion, we use the same overall dimension of the model and HEM system parameters as defined for synthetic example 1, but replace the earlier anomalous bodies with newer anomalous bodies. The background half-space resistivity is still 100 Ω⋅m. Fig. 4(a) shows the synthetic model 2 as layers 1–5 at 0, 5, 11, 18 and 27 m depths, respectively. Upper and lower bodies are 21 and 31 m thick with resistivity of 10 and 20 Ω⋅m, respectively. Figure 4. View largeDownload slide 3-D inversion result with good data coverage for synthetic model 2. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 4. View largeDownload slide 3-D inversion result with good data coverage for synthetic model 2. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. In order to test the stability of different inversion approaches, we calculate two synthetic data set with different line spacing. In the first case, survey stations are evenly distributed in the central area of the model with a spacing of 40 m along eight flight lines having line spacing of 100 m to give 19 sites on each line and a total of 152 survey locations. In the second case, we change the line spacing to 180 m, so there are total five lines and 95 survey locations only. One per cent Gaussian noise is added to both the synthetic data. We compare inversion results of db2, db6 and db10 Daubechies wavelets with inversion results from traditional space-domain methods for both the data set. Same selection strategies for λ and ε is used as done for the synthetic model 1. The inversion results of two different synthetic data sets are shown in Figs 4 and 6, respectively, for layers 1−5. Fig. 4 reveals that when there is good data coverage, then all the five inversion methods recover the true model well and all of them converge to the target data misfit of rms 1 (Fig. 5a). L1-norm inversion could not recover the second body well in layer 5 and recovers lower resistivity of the bodies in all the layers. L2-norm and wavelet-based inversions recovers both bodies well in xy-plane for all the layers except db6 for layer 5. Wavelet-based inversions recovers a closer resistivity to the true resistivity for layers 2–4, while the L2- and L1-norm inversions recover much lower resistivity for the layer 3 than those for layers 2 and 4. The smoothness constraint in space domain always forces the model to change smoothly. But wavelet-domain inversions have constraints in multiscale on the model complexity to let it adapt better. If we have a good data density, then we can recover much sharper boundary in a smoother background through wavelet-domain inversions than the inversion methods based on smoothness constraint in space domain. Figure 5. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with good data coverage as shown in Fig. 3. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 5. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with good data coverage as shown in Fig. 3. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Fig. 6 shows the 3-D inversion result for larger flight line spacing data for same synthetic model 2. Since we use much less data than that in Fig. 4, the inversion is much more ill-posed. Figs 6(b) and (f) show although the resolution is not high in this situation, still these two inversion methods (L2-norm and db10) can recover the model better, which means the db10 wavelet-based inversion is as stable as the L2-norm inversion in the space domain. The other three inversion methods can also get a good data misfit (as shown in Fig. 7a), but the recovered models are different from the true model. From definition of norm, L1-norm prefers to generate a sparser model variation than the L2-norm, so L1-norm is not as stable as the L2-norm for seriously ill-posed inversion problems. When there is enough data as in inversion of synthetic model 1 (Fig. 2), the L1-norm-based inversion can give much sharper boundary than L2-norm-based inversion. Comparing Figs 6(d)–(f) reveals when the number N of dbN increases then the inversion becomes more stable and recovers more features. This is because when N is larger, the regularity of the wavelet is also larger and the basis function is smoother, which makes the recovered model smoother and the inversion more stable. We also observe from Fig. 6(c) that the inversion result of db2 wavelet is simpler with clear discontinuity. Theoretically, if the model is simple and has quite clear discontinuities, the db2 wavelet can give sharper boundaries because its basis function is not smooth. Figure 6. View largeDownload slide 3-D inversion result for synthetic model 2 with poor data coverage. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 6. View largeDownload slide 3-D inversion result for synthetic model 2 with poor data coverage. (a) True model and models constructed by (b) L2-norm inversion in the space domain, (c) L1-norm inversion in the space domain, (d) wavelet inversion using the db2 basis function, (e) wavelet inversion using db6 and (f) wavelet inversion using db10. The white dots are the locations of surveying points. Layers 1–5 correspond to the depth at 0, 5, 11, 18 and 27 m, respectively. Figure 7. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with poor data coverage as shown in Fig. 4. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 7. View largeDownload slide Plots of inversion parameters obtained from the inversion for the synthetic model 2 with poor data coverage as shown in Fig. 4. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. For better understanding of recovery of the model and sensitivity information in the wavelet-domain inversion, we show wavelet coefficients of the right-hand side of eq. (14) (the minus the gradient of the objective function) in Figs 8(b)–(d) for the first layer (at the surface) of the good data coverage case. Fig. 8(a) shows the minus gradient of objective function in space domain. The distribution in Fig. 8(a) resembles the shape of two anomalous bodies that means if we select a suitable step length then we can use this minus gradient to update the model directly. The wavelet representation of the minus gradient of the objective function transformed by db2, db6 and db10 wavelets are shown in Figs 8(b), (c) and (d), respectively. Since we use three times 1-D wavelet to transform it subsequently in three directions (x, y and z), the wavelet coefficients are mainly focused in one corner of the model. They mainly locate at the bottom left corner. Most of the coefficients are close to zero in Figs 8(b) and (c) for db2 and db6, respectively, except in Fig. 8(d) for db10 coefficients which has a more flat distribution than db2 and db6 because it has a larger vanishing moment. Figure 8. View largeDownload slide Minus gradient of objective function from the first layer of synthetic model 2 for good data coverage case. (a) In space domain, (b) in wavelet domain transformed by db2 wavelet, (c) in wavelet domain transformed by db6 wavelet and (d) in wavelet domain transformed by db10 wavelet. Figure 8. View largeDownload slide Minus gradient of objective function from the first layer of synthetic model 2 for good data coverage case. (a) In space domain, (b) in wavelet domain transformed by db2 wavelet, (c) in wavelet domain transformed by db6 wavelet and (d) in wavelet domain transformed by db10 wavelet. 4 FIELD DATA INVERSION In order to further test the resolution and stability of wavelet-based inversion technique in real world, we inverted a field FDAEM data set. This HEM data set was acquired by NGU in Byneset, Norway. Byneset is located in Trondheim municipality in Mid Norway. Fig. 9 shows the Quaternary geology map of Byneset. Most of the area is covered by marine sediments deposited during the deglaciation around 10 000 yr ago. Due to glacio-isostatic rebound, the sediments are now exposed on surface with a ground surface elevation around 160 m a.s.l. (Reite et al.1999). Central parts of Byneset are covered by ravines and about 100 landslide scars (Solberg et al.2016). A recent landslide in 2012 happened due to presence of the quick clay in the area and it is marked by a green curve in the middle of Fig. 9. The bedrock in the area is dominated by chlorite slate in west and phyllites in the east. The phyllite may contain graphite which is an electronically conductive mineral (Solberg et al.2016). Figure 9. View largeDownload slide Quaternary geological map of Byneset with location of HEM flight lines. Figure 9. View largeDownload slide Quaternary geological map of Byneset with location of HEM flight lines. This data set was collected mostly at 100 m line interval along 60 flight lines shown by black lines in Fig. 9. In-phase and quadrature components of induced secondary magnetic field in part per million (ppm) of transmitted primary field were recorded at 0.1 s interval resulting in 3 m data interval. The data were collected at five different frequencies, three of them in coplanar (34 kHz, 6.6 kHz and 880 Hz) and two of them in co-axial (7 kHz and 980 Hz) setting of transmitter and receiver coils. Details of the processing and spatially constrained inversion of the data are discussed by Baranwal et al. (2015, 2017). For inverting this field data set, we build a model with 231 × 261 × 21 cells. The cell sizes in the x- and y-directions are 30 m, while in the z-direction the thickness of the top layer is 5 m and the thicknesses of the underlying layers increases with a ratio of 1.2. The initial resistivity of background half-space is set to be 100 Ω⋅m. Since the 3-D inversion forward modeling speed is not very fast, we extract total 4148 data from the original data set with a surveying point space of about 50 m on the flight lines. The noise floor is set to be 5 ppm plus 1 per cent of the amplitude. Flight line spacing ranges from about 60 to 300 m (Fig. 9). In order to get a stable inversion, we only test db10 wavelet-based inversion and compare its result with the result of L2-norm-based inversion in the space domain. For both of these two inversion methods, we set the initial value of the trade-off parameter λ to be 1000. Data misfit reduces very slowly so we set λ to be 10 per cent of its previous value. After it reduces to 1, it is kept constant in the subsequent iterations. The inversion result of L2-norm and db10 wavelet methods are illustrated in Fig. 10 for top eight layers (at 0, 5, 11, 18, 27, 37, 50 and 65 m depths). We can see that the two inversion methods give overall similar results which has good agreement with the geology as shown in Fig. 9. The whole surveying area is surrounded by high resistivity outcropping bedrock (observed in the field). The central region of the area is conductive due to thick layers of marine deposits. It comprises a thin and moderately resistive (10–100 Ω⋅m) layer at the top (in layers 1 and 2 of Fig. 10) interpreted as leached clay or silty sediments with a thick unleached marine clay (< 10 Ω⋅m) underneath. Unleached marine sediment is interpreted thickest in the central region of Byneset from 3-D inversion of HEM data. Exposed to deep-seated high resistive bedrock is interpreted at different locations in the area. Upper part of the horizontal depth slices in Fig. 10 ( y > 7 027 500 m) recovered by L2-norm and db10 wavelet look very alike, while lower part ( y < 7 027 500 m) recovered by these two inversions are quite different. There is a good data density to ensure a stable inversion in the upper part where the flight line spacing is smaller. While in the lower part, the line spacing is quite large, therefore it is difficult for smoothing constraint in space domain to bridge the structures from two flight lines because of smaller footprint (about 70 m) of the AEM system in comparison to the line spacing (150–300 m). Figs 10(a)–(h) reveal that the recovered model from L2-norm in the lower part is not continuous. The db10 wavelet-based inversion (Figs 10i–p) gives a smoother result in the lower part because inherent constraint of wavelet inversion can recover the main structure correctly. Inversion result around y = 7 028 000 m is not reliable from both the inversions because the data have large noise. Figure 10. View largeDownload slide 3-D inversion result of field data acquired in Byneset, Norway. (a)–(h) are the horizontal depth slices of top eight layers of L2-norm-based 3-D inversion. (i) – (p) are the horizontal depth slices of top eight layers of db10 wavelet-based 3-D inversion. Layers–8 correspond to the depth at 0, 5, 11, 18, 27, 37, 50 and 65 m, respectively. Figure 10. View largeDownload slide 3-D inversion result of field data acquired in Byneset, Norway. (a)–(h) are the horizontal depth slices of top eight layers of L2-norm-based 3-D inversion. (i) – (p) are the horizontal depth slices of top eight layers of db10 wavelet-based 3-D inversion. Layers–8 correspond to the depth at 0, 5, 11, 18, 27, 37, 50 and 65 m, respectively. Plots of the inversion parameters during iterations are shown in Fig. 11. L2-norm-based inversion converges a little faster than the db10 wavelet-based inversion. Because we do not provide an accurate noise floor to the data, therefore the data misfit (rms) of these two methods reduces from initial 9.52 to about 2 only. However, we observe a good fit between field data and predicted data visually from both the inversions. As illustrated in Fig. 11, we have reset the trade-off parameter λ to be 100 after four iterations for db10 inversion. This is done because trade-off parameter λ is allowed to be smaller to get a suitable model update in the beginning of the inversion and later it is reset to a large λ again to ensure a stable inversion. Figure 11. View largeDownload slide Plots of inversion parameters obtained from the inversion for the field data from Byneset as shown in Fig. 9. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. Figure 11. View largeDownload slide Plots of inversion parameters obtained from the inversion for the field data from Byneset as shown in Fig. 9. (a) Data misfit rms, (b) objective function Φ, (c) model complexity ϕm and (d) trade-off parameter λ. To show the vertical resolution of the inverted models obtained from these two inversions of the field data, we illustrate vertical sections at y = 7 030 191 m from both the inversions in Fig. 12. This vertical section passes through 2012 landslide location. Two vertical slices obtained from L2-norm and db10 inversions have similar layered structure. The 2012 landslide took place at around x = 557 000 m, where we see a clear discontinuity of the conductive layers (unleached clay) in both the images. Presence of the quick clay at this location is confirmed by soil samples analysis (Solberg et al.2016). The result shown in Fig. 12 has a good agreement with HEM and ERT inversion results shown in Solberg et al. (2016). The inversion result of db10 wavelet (Fig. 12b) provides a clearer image of unleached clay layer than the L2-norm inversion (Fig. 12a). Due to the smoothing constraint in the space-domain method, the result of L2-norm inversion connects the surface and middle thick conductive layer of unleached clay at several locations. The field data and predicted data illustrated in Figs 12(c) and (d) show that the two inversion methods fit the data to a similar level. Figure 12. View largeDownload slide Vertical sections of inversion result at y = 7 030 191 m from (a) L2-norm-based method and (b) db10 wavelet-based method; (c) and (d) are data fitting for in-phase and quadrature components. Figure 12. View largeDownload slide Vertical sections of inversion result at y = 7 030 191 m from (a) L2-norm-based method and (b) db10 wavelet-based method; (c) and (d) are data fitting for in-phase and quadrature components. In Fig. 13, we illustrate the first (at 0 m depth) and fourth layer (at 18 m depth) inversion results at different iterations to show how the two inversion methods work. The two methods have different convergence process, so it is not suitable to compare the result from the same iteration, therefore we compare the results from the iterations with closer rms values. Since the AEM data are dense enough and one site relates to several nearby sites. The smoothing constraint in L2-norm-based inversion recovers the local structures quickly through the inversion (as shown in Figs 13a–d for the first layer and in Figs 13e–h for the fourth layer). In the wavelet-based inversion, we use a large lambda in the beginning, which enforces a very strong sparsity constraint and only allows few large-scale coefficients, so the inversion recovers the main structure of this model first. Figs 13(i)–(l) (for first layer) and Figs 13(m)–(p) (for fourth layer) show this property of the wavelet inversion by recovering the main structure of the entire model without any local detail in the first several iterations (Figs 13j and k and n and o). More local details appear in later iterations of the inversion (Figs 13l and p). Therefore, the wavelet-based inversion can recover structures of the model at different scales. Figure 13. View largeDownload slide Inversion results of first and fourth layers from two inversion approaches in different iterations. (a)–(d) depict first layer at the surface obtained from L2-norm inversion at rms values 8.51, 5.38, 4.06 and 1.90. (e)–(h) depict fourth layer at 18 m depth obtained from L2-norm-based inversion at rms values 8.51, 5.38, 4.06 and 1.90. (i)–(l) depict first layer at the surface obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. (m)–(p) depict fourth layer at 18 m depth obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. Figure 13. View largeDownload slide Inversion results of first and fourth layers from two inversion approaches in different iterations. (a)–(d) depict first layer at the surface obtained from L2-norm inversion at rms values 8.51, 5.38, 4.06 and 1.90. (e)–(h) depict fourth layer at 18 m depth obtained from L2-norm-based inversion at rms values 8.51, 5.38, 4.06 and 1.90. (i)–(l) depict first layer at the surface obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. (m)–(p) depict fourth layer at 18 m depth obtained from db10 wavelet-based inversion at rms values 8.94, 4.92, 3.92 and 2.30. 5 CONCLUSIONS We propose a new stable and well-resolved 3-D inversion method based on the wavelet transform and sparsity constraint to invert FDAEM data. This new method takes advantages of multiscale property of the wavelet representation which can be used to recover different scale structure of the model in different iterations. The wavelet-based inversion gives constraints in multiscale on the model complexity, while the smooth constraint in space domain is affected more by nearby cells. Synthetic and field data inversion tests show that convergence is stable and the inversion can produce as good resolution as L1-norm approaches in the space domain for sufficiently dense data. Both the wavelet-based inversion and L1-norm inversion give more focusing result than L2-norm inversion. The sparsity constraint on model complexity in wavelet domain can overcome the problem of sparse data distribution which always results in discontinuous structure, and generates a data-based resolution that adaptive to the data distribution rather than the model-based resolution that depends on the model discretization. The method has the potential to improve inversion resolution for many other EM methods without compromising the stability of the inversion. Acknowledgements We want to thank the reviewers and editors for their suggestions that helped to clarify this paper. This research is funded by the Key National Research Project of China (2016YFC0303100, 2017YFC0601900). This research is also funded by Key Program of National Natural Science Foundation of China (41530320) and Natural Science Foundation for Young Scientist (41404093). NGU is thanked for making the Byneset HEM data available. REFERENCES Abers G.A., Roecker S.W., 1991. Deep structure of an arc-continent collision: earthquake relocation and inversion for upper mantle P and S wave velocities beneath Papua New Guinea. J. Geophys. Res. , 96( B4), 6379– 6401. https://doi.org/10.1029/91JB00145 Aster R.C., Borchers B., Thurber C.H., 2013. Parameter Estimation and Inverse Problems , Academic Press. Baranwal V.C., Dalsegg E., Tønnesen J.F., Rønning J.S., Solberg I.L., Rodionov A., Dretvik H., 2015. Mapping of marine clay layers using airborne EM and ground geophysical surveys at Byneset, Trondheim municipality. Geological Survey of Norway (NGU) Report 2015.006 . Available at: https://www.ngu.no/en/publikasjon/mapping-marine-clay-layers-using-airborne-em-and-ground-geophysical-methods-byneset. Baranwal V.C., Rønning J.S., Solberg I.-L., Dalsegg E., Tønnesen J.F., Long M., 2017, Investigation of a sesnsitive clay landslide area using frequency-domain helicopter-borne EM and ground geophysical methods, in Landslides in Sensitive Clays, Series in Advances in Natural and Technological Hazards Research , eds Thakur V., L’Heureux J.-S., Locat A., p. 46. Springer International Publishing. Beamish D., 2003. Airborne EM footprints, Geophys. Prospect. , 51( 1), 49– 60. https://doi.org/10.1046/j.1365-2478.2003.00353.x Chiao L.Y., Kuo B.Y., 2001. Multiscale seismic tomography, Geophys. J. Int. , 145( 2), 517– 527. https://doi.org/10.1046/j.0956-540x.2001.01403.x Chiao L.Y., Liang W.T., 2003. Multiresolution parameterization for geophysical inverse problems, Geophysics , 68( 1), 199– 209. https://doi.org/10.1190/1.1543207 Cox L.H., Wilson G.A., Zhdanov M.S., 2010. 3D inversion of airborne electromagnetic data using a moving footprint, Explor. Geophys. , 41( 4), 250– 259. https://doi.org/10.1071/EG10003 Daubechies I., 1992. Ten Lectures on Wavelets , Society for Industrial and Applied Mathematics. Daubechies I., Defrise M., De Mol C., 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Comm. Pure appl. Math. , 57( 11), 1413– 1457. https://doi.org/10.1002/cpa.20042 Davis K., Li Y., 2011. Fast solution of geophysical inversion using adaptive mesh, space-filling curves and wavelet compression, Geophys. J. Int. , 185( 1), 157– 166. https://doi.org/10.1111/j.1365-246X.2011.04929.x Delost M., Virieux J., Operto S., 2008. First-arrival traveltime tomography using second generation wavelets, Geophys Prospect , 56( 4), 505– 526. https://doi.org/10.1111/j.1365-2478.2008.00710.x Donoho D.L., 2006. For most large underdetermined systems of linear equations the minimal?1-norm solution is also the sparsest solution, Comm. Pure appl. Math. , 59( 6), 797– 829. https://doi.org/10.1002/cpa.20132 Ekblom H., 1987, The L1-estimate as limiting case of an Lp-or Huber-estimate, in Statistical Data Analysis based on the L1-Norm and Related Methods , ed. Doge Y., pp. 109– 116. Elsevier Science Publ. Co., Inc. Fang H., Zhang H., 2014. Wavelet-based double-difference seismic tomography with sparsity regularization, Geophys. J. Int. , 199( 2), 944– 955. https://doi.org/10.1093/gji/ggu305 Farquharson C.G., 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions, Geophysics , 73( 1), K1– K9. https://doi.org/10.1190/1.2816650 Haber E., 2014, Computational Methods in Geophysical Electromagnetics , SIAM. Hung S.H., Chen W.P., Chiao L.Y., 2011. A data‐adaptive, multiscale approach of finite‐frequency, traveltime tomography with special reference to P and S wave data from central Tibet, J. geophys. Res. , 116, B06307, doi:10.1029/2010JB008190. Hung S.H., Chen W.P., Chiao L.Y., Tseng T.L., 2010. First multi‐scale, finite‐frequency tomography illuminates 3‐D anatomy of the Tibetan plateau, Geophys. Res. Lett. , 41, 7025– 7034, doi:10.1002/2014GL061644. Li Y., Oldenburg D.W., 1999. Rapid construction of equivalent sources using wavelets, in 69th Annual International Meeting, Expanded Abstracts , Vol. 18, pp. 374– 378. SEG, Houston TX. Li Y., Oldenburg D.W., 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method, Geophys. J. Int. , 152( 2), 251– 265. https://doi.org/10.1046/j.1365-246X.2003.01766.x Li Y., Oldenburg D.W., 2010. Rapid construction of equivalent sources using wavelets, Geophysics , 75( 3), L51– L59. https://doi.org/10.1190/1.3378764 Liu Y.H., Yin C.C., 2013. 3D inversion for frequency-domain HEM data, Chin. J. Geophys. , 56, 4278– 4287. Liu Y.H., Yin C.C., 2016. 3D inversion for multipulse airborne transient electromagnetic data, Geophysics , 81( 6), E401– E408. https://doi.org/10.1190/geo2015-0481.1 Loris I., Nolet G., Daubechies I., Dahlen F., 2007. Tomographic inversion using? 1 -norm regularization of wavelet coefficients, Geophys. J. Int. , 170( 1), 359– 370. https://doi.org/10.1111/j.1365-246X.2007.03409.x Nittinger C.G., Becken M., 2016. Inversion of magnetotelluric data in a sparse model domain, Geophys. J. Int. , 206( 2), 1398– 1409. https://doi.org/10.1093/gji/ggw222 Oldenburg D., McGillivray P., Ellis R., 1993. Generalized subspace methods for large-scale inverse problems, Geophys. J. Int. , 114( 1), 12– 20. https://doi.org/10.1111/j.1365-246X.1993.tb01462.x Portniaguine O., Zhdanov M.S., 1999. Focusing geophysical inversion images, Geophysics , 64( 3), 874– 887. https://doi.org/10.1190/1.1444596 Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1996. Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing , Cambridge University Press. Reite A.J., Sveian H., Erichsen E., 1999. Trondheim fra istid til nåtid – landskapshistorie og løsmasser , Geological Survey of Norway (NGU), Gråsteinen 5 (in Norwegian). Available at: https://www.ngu.no/publikasjon/gr%C3%A5steinen-5-trondheim-fra-istid-til-n%C3%A5tid. Simons F.J.et al.  , 2011. Solving or resolving global tomographic models with spherical wavelets, and the scale and sparsity of seismic heterogeneity, Geophys. J. Int. , 187( 2), 969– 988. https://doi.org/10.1111/j.1365-246X.2011.05190.x Siripunvaraporn W., Egbert G., Lenbury Y., Uyeshima M., 2005. Three-dimensional magnetotelluric inversion: data-space method, Phys. Earth planet. Inter. , 150( 1–3), 3– 14. https://doi.org/10.1016/j.pepi.2004.08.023 Solberg I.L., Long M., Baranwal V.C., Gylland A.S., Rønning J.S., 2016. Geophysical and geotechnical studies of geology and sediment properties at a quick-clay landslide site at Esp, Trondheim, Norway, Eng. Geol. , 208, 214– 230. https://doi.org/10.1016/j.enggeo.2016.04.031 Spakman W., Bijwaard H., 2001. Optimization of cell parameterizations for tomographic inverse problems, Pure appl. geophys. , 158( 8), 1401– 1423. https://doi.org/10.1007/PL00001227 Zhang H., Thurber C., 2005. Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: application to Parkfield, California, J. geophys. Res. , 110( B4), B04303, doi:10.1029/2004JB003186. https://doi.org/10.1029/2004JB003186 © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations