Generalized L p -norm joint inversion of gravity and magnetic data using cross-gradient constraint
Saeed Vatankhah, Shuang Liu, Rosemary A. Renaut, Xiangyun Hu, Mostafa Gharloghi
GGeneralized L p -norm joint inversion of gravity andmagnetic data using cross-gradient constraint Saeed Vatankhah , , Shuang Liu , Rosemary A. Renaut , Xiangyun Hu ,Mostafa Gharloghi Hubei Subsurface Multi-scale Imaging Key Laboratory, Institute of Geophysics and Geomatics,China University of Geosciences, Wuhan, China Institute of Geophysics, University of Tehran, Tehran, Iran School of Mathematical and Statistical Sciences, Arizona State University, Tempe, AZ, USA
SUMMARY
A generalized unifying approach for L p -norm joint inversion of gravity and mag-netic data using the cross-gradient constraint is presented. The presented frameworkincorporates stabilizers that use L , L , and L -norms of the model parameters,and/or the gradient of the model parameters. Furthermore, the formulation is de-veloped from standard approaches for independent inversion of single data sets, and,thus, also facilitates the inclusion of necessary model and data weighting matricesthat provide, for example, depth weighting and imposition of hard constraint data.The developed efficient algorithm can, therefore, be employed to provide physically-relevant smooth, sparse, or blocky target(s) which are relevant to the geophysicalcommunity. Here, the nonlinear objective function, that describes the inclusion ofall stabilizing terms and the fit to data measurements, is minimized iteratively byimposing stationarity on the linear equation that results from applying linearization a r X i v : . [ phy s i c s . g e o - ph ] J a n S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi of the objective function about a starting model. To numerically solve the result-ing linear system, at each iteration, the conjugate gradient (CG) algorithm is used.The general framework is then validated for three-dimensional synthetic models forboth sparse and smooth reconstructions, and the results are compared with thoseof individual gravity and magnetic inversions. It is demonstrated that the presentedjoint inversion algorithm is practical and significantly improves reconstructed mod-els obtained by independent inversion.
Key words:
Potential field surveys, gravity and magnetic, have been reported as effective strategies for de-lineating subsurface geological targets. They are applied in a wide range of studies including,for example, oil and gas exploration, mining applications, and mapping the basement topogra-phy (Nabighian et al. 2005; Blakely 1996). These surveys are relatively cheap, non-destructivepassive remote sensing methods, and can provide valuable information of the subsurface tar-gets. Yet, they only require the measurement of variations in the Earth’s natural fields thatare caused by changes in the physical properties of the subsurface rocks. In the interpreta-tion process, the acquired survey data can be used in an automatic inversion algorithm forestimation of specific parameters of a subsurface target, for example its geometry or physicalproperties. It is well known, however, that the potential field inversion problem is ill posed.Identifiability of stable and physically-relevant solutions is then obtained by application ofsuitable regularization strategies. An independent solution of the inverse problem for eithergravity, or magnetic, data for the survey area will only provide information about the densityor susceptibility, respectively, of the subsurface. On the other hand, complementary solutionof the inverse problem for both data sets can be used to reveal both density and magneti-zation variations present in a subsurface target. Thus, it is more appropriate to perform asimultaneous joint inversion that uses both data sets. Combined with regularization, this isan effective strategy for yielding a reliable subsurface geological model, that simplifies theinterpretation of the subsurface target(s). Thus, the development of efficient and stable jointinversion algorithms has received increased attention in the geophysical community.Many different techniques have been developed for the simultaneous joint inversion ofgeophysical data sets. Generally, these techniques can be categorized into two main groups: eneralized L p -norm joint inversion p formulation for efficient simultaneous joint inversion of gravityand magnetic data sets. S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Different types of stabilizers have been adopted for the inversion of potential field data,dependent on the desired smoothness, or otherwise, of model features that are to be recovered.For example, it may be appropriate to reconstruct a model which only represents the large-scale features of the subsurface under the survey area without any arbitrary discontinuities.This is achieved with the maximum smoothness stabilizer which uses a L -norm I of the gra-dient of the model parameters, (Constable et al. 1987; Li & Oldenburg 1996; Pilkington 1997;Li & Oldenburg 1998). On the contrary, when it is anticipated that the subsurface structureexhibits discontinuities, stabilization can be achieved by imposing the L -norm, or L -norm,on the gradient of the model parameters (Farquharson & Oldenburg 1998; Portniaguine &Zhdanov 1999; Bertete-Aguirre et al. 2002; Vatankhah et al. 2018a; Fournier & Oldenburg2019). Alternatively, when it can be assumed that the subsurface targets are localized andcompact, it is more appropriate to apply the L or L -norms directly on the model param-eters (Last & Kubik 1983; Portniaguine & Zhdanov 1999; Barbosa & Silva 1994; Zhdanov& Tolstaya 2004; Ajo-Franklin et al. 2007; Vatankhah et al. 2014a; Vatankhah et al. 2014b;Sun & Li 2014; Zhdanov 2015; Vatankhah et al. 2017). In the potential field literature, sta-bilization by application of the L -norm on the model parameters is usually referred to asthe compactness constraint, whereas application of the L -norm, or L -norm, on the gradientis referred to as total variation (TV) and minimum gradient support (MGS), stabilization,respectively. A unifying approach for application of these constraints for single potential fieldinversion is presented in Vatankhah et al. (2019a). This approach also includes the modifica-tion of the stabilizers to account for additional model and data weighting matrices, such asrequired for imposition of depth weighting and hard constraint conditions, (Li & Oldenburg1996; Boulanger & Chouteau 2001). Note that the use of depth weighting counteracts thenatural rapid decay of the kernels, dependent on the specific kernel, whether magnetic orgravity. This then facilitates the contribution of all prisms at depth to the surface measure-ments with an approximately equal probability through the inversion algorithm. The hardconstraint weighting is used to impose available geological information on the reconstructedmodel. Consequently, any practical algorithm for the joint inversion of gravity and magneticdata should also incorporate such weighting schemes. Here, we extend this unifying frameworkfor simultaneous joint inversion of gravity and magnetic data sets in conjunction with the useof the cross-gradient constraint.The well-known and widely-used formula of Fregoso & Gallardo (2009) for the joint in- I The L p -norm of a vector x ∈ R n is defined as (cid:107) x (cid:107) p = ( (cid:80) ni | x i | p ) /p , p ≥
1, and (cid:107) x (cid:107) counts the number of nonzeroentries in x . eneralized L p -norm joint inversion To formulate the problem, we use the well-known strategy for linear inversion of potentialfield data in which the subsurface is divided into a set of rectangular prisms with fixed sizebut unknown physical properties (Li & Oldenburg 1996; Boulanger & Chouteau 2001). Here,it is assumed that there is no remanent magnetization, and that self-demagnetization effectsare also negligible. For ease of exposition, we first introduce some basic notation for stackingof vectors (matrices) and generation of block diagonal matrices. We use block stack ( · , · ) toindicate the stacking of vectors (or matrices) with the same number of columns in one vector(or matrix). Further, block diag ( A, B ) indicates a block diagonal matrix of size ( m A + m B ) × ( n A + n B ) when A and B are of sizes ( m A × n A ) and ( m B × n B ), respectively. Both definitionsextend immediately for more than two entities.We suppose that m measurements are taken for two sets of potential field data II . Theseare the vertical components of the gravity and total magnetic fields, and they are stacked invectors d obs1 , and d obs2 , each of length m , respectively. The unknown physical parameters, the II We could assume different numbers of measurements for each field, m and m but for simplicity of the discussionwe immediately assume m = m = m . S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi density and the susceptibility, of n prisms are also stacked in vectors m , and m , respec-tively. The data vectors and model parameters are then stacked consistently in vectors d obs =block stack (cid:0) d obs1 , d obs2 (cid:1) ∈ R m , and m = block stack ( m , m ) ∈ R n . The measurements areconnected to the model parameters via G m = d obs where G = block diag ( G , G ) ∈ R m × n ,and G and G are the linear forward modeling operators for gravity and magnetic kernels,respectively. There are different alternative formulas which can be used to compute the entriesof matrices G and G . Here, we use the formulas developed by Ha´az (1953), for computingthe vertical gravitational component, and Rao & Babu (1991), for the total magnetic fieldanomaly, of a right rectangular prism, respectively.The goal of the inversion is to find geologically plausible models m and m that predict d obs1 and d obs2 , respectively, via a simultaneous joint algorithm that also facilitates the incor-poration of relevant weighting matrices in the algorithm. We formulate the joint inversionfor the determination of the model parameters m and m as the minimization of the globalobjective function, in which parameters α and λ are relative weighting parameters for therespective terms, P ( α,λ ) ( m ) = (cid:107) W d ( d obs − G m ) (cid:107) + α (cid:107) W D ( m − m apr ) (cid:107) + λ (cid:107) t ( m ) (cid:107) . (1)The data misfit term, (cid:107) W d ( d obs − G m ) (cid:107) , measures how well the calculated data reproducethe observed data. Diagonal matrix W d = block diag ( W d , W d ) ∈ R m × m , where W d and W d are diagonal weighting matrices for the gravity and magnetic data, respectively. Herewe suppose that these diagonal elements are the inverses of the standard deviations of theindependent, but potentially colored, noise in the data. Zhang & Wang (2019) consideredan alternative weighting based on the individual row norms which could also be used here.The stabilizer, (cid:107) W D ( m − m apr ) (cid:107) , controls the growth of the solution with respect to theweighted norm, and is especially significant as it determines the structural qualities of thedesired solution. Here, this stabilizer is presented through a general L -norm formulation, butwe will discuss how different choices of W and D lead to different L p -norm stabilizations.In (1) the vector m apr = block stack ( m apr1 , m apr2 ) ∈ R n is an initial starting model thatmaybe known from previous investigations. It is also possible to set m apr = . The linkbetween the gravity and magnetic models in this inversion algorithm is the cross-gradientfunction t ( m ) ∈ R n . For this study, we assume that the model structures for m and m areapproximately the same, and thus that it is important to measure the structural similarityusing the cross-gradient constraint which will be approximately zero for models with similar eneralized L p -norm joint inversion t ( m ) = ∇ m ( x, y, z ) × ∇ m ( x, y, z ) , (2)where ∇ indicates the gradient operator (Gallardo & Meju 2003) and structural similarityis achieved when t ( m ) = , see Appendix A for details. As noted already in Section 1, thiscorresponds to the case in which the gradient vectors are in the same or opposite direction, or,alternatively, one of them is zero (Gallardo & Meju 2003; Gallardo & Meju 2004; Tryggvason& Linde 2006; Gallardo 2007; Fregoso & Gallardo 2009; Haber & Gazit 2013; Fregoso et al.2015). From a geological viewpoint this means that if a boundary exists then it must be sensedby both methods in a common orientation regardless of how the amplitude of the physicalproperty changes (Gallardo & Meju 2003). This means that information that is contained inone model is relevant to the other model and vice versa. Therefore, structures determinedby one model can assist with the identification of structures in the other model, and, as aconsequence, the structures of the two models can correct each other throughout the jointinversion process (Gallardo & Meju 2003; Haber & Gazit 2013). On the other hand, while itis assumed that both models have similar structures at similar locations, it is also possiblefor one model to have a structure in a location where the other model has none (Haber &Gazit 2013). Further details about characteristic properties of the cross-gradient constraintare provided in Gallardo & Meju (2003), Tryggvason & Linde (2006), and Fregoso & Gallardo(2009).The stabilizing term, (cid:107) W D ( m − m apr ) (cid:107) , in (1) has a very significant impact on thesolution that is obtained by minimizing (1). Depending on the type of desired model to berecovered, there are many choices that can be considered for this stabilization, and that havebeen extensively adopted by the geophysical community. Here, we show how it is possibleto use the given weighted L -norm regularizer to approximate different L p -norm stabilizers,0 ≤ p ≤
2, (Vatankhah et al. 2019a). Suppose D is the identity matrix, D = I n , and W isselected as W = block diag ( W , W ) ∈ R n × n , in which, W i = ( W depth ) i ( W h ) i ( W L p ) i ∈ R n × n , i = 1 , . (3)Here, the diagonal weighting matrix ( W L p ) i ∈ R n × n is defined, assuming entries are calculatedelementwise, by ( W L p ) i = diag (cid:16) / (( m i − m apri ) + (cid:15) ) (2 − p ) / (cid:17) , i = 1 , . (4)When p = 0 and p = 1, compact L -norm and L -norm solutions are obtained, respectively.The choice p = 2 provides the L -norm solution of the model parameters. The parameter (cid:15) is S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi a small positive number, 0 < (cid:15) (cid:28)
1, which is added to avoid the possibility of division by zero,and has an important effect on the solution. When (cid:15) is very small the solutions are sparse,while for large values the solutions are smooth. Further discussion on the impact of (cid:15) is given,for example, in Last & Kubik (1983), Farquharson & Oldenburg (1998), Zhdanov & Tolstaya(2004), Ajo-Franklin et al. ( 2007), Fiandaca et al. (2015), and Fournier & Oldenburg (2019).On the other hand, it is also possible to chose D to provide an approximation to the gradientof the model parameters. Suppose, for example, that D = block stack ( D x , D y , D z ) ∈ R n × n ,where D x , D y , and D z are square and provide discrete approximations for derivative operatorsin x , y , and z -directions. Then, defining n × n to be the zero matrix of size 3 n × n , and setting D = D , we can use the matrix D = D n × n n × n D ∈ R n × n , (5)so that D m yields the approximate gradient of m . More details about the structure of thematrices D x , D y , and D z can be found in Li & Oldenburg (2000), Leli´evre & Oldenburg(2009), and Leli´evre & Oldenburg (2013). Then, with this definition for D , and again usingelement-wise calculations, ( W L p ) i in (4) is replaced by( W L p ) i = diag (cid:16) / (cid:0) ( D x ( m − m apri )) + ( D y ( m − m apri )) + ( D z ( m − m apri )) + (cid:15) (cid:1) (2 − p ) / (cid:17) . (6)But now, for the multiplication in equation (1) to be dimensionally consistent, W is replacedby W = block diag ( W , W , W , W , W , W ) ∈ R n × n . (7)In the definition of ( W L p ) i , given in (6), picking p = 2 produces a solution with minimumstructure yielding a smooth model without arbitrary discontinuities (Constable et al. 1987; Li& Oldenburg 1996; Pilkington 1997). If we anticipate that there are true discontinuous jumps,then, it is possible to take p = 1 or p = 0 for a total variation (TV) or minimum gradientsupport (MGS) stabilization, respectively. In summary, all aforementioned definitions indicatehow, dependent on the choices of p and D in W , it is possible to use the objective function(1) with a desired stabilizer. Well-known stabilizers, including TV, MGS, minimum structure,compactness, and L -norm can all be incorporated in a joint inversion methodology. Moreover,this unifying framework allows the use of additional stabilizers, which are not common inpotential field inversion, simply by changing the choice of p .In (3), the diagonal depth weighting matrix, ( W depth ) i = diag(1 / ( z j + z ) β i ) is used to eneralized L p -norm joint inversion z j is the mean depth of prism j , z depends both upon prism size and the height of theobserved data. With application of appropriate depth weighting, as determined by parameter β i , all prisms participate with an approximately equal probability in the inversion process.The diagonal hard constraint matrix ( W h ) i , is generally an identity matrix. If geological in-formation, or prior investigations in the survey area, can provide the values of the modelparameter for some prisms, then, the information is included in m apri , and the correspondingdiagonal entries of ( W h ) i are set to a large value (Boulanger & Chouteau 2001; Vatankhahet al. 2014a; Vatankhah et al. 2018b). These known parameters are kept fixed during theiterative minimization of (1). Equivalently, the inversion algorithm searches only for unknownmodel parameters. As an important aside, note that all matrices D x , D y , D z , W depth , W h ,and W L p , are sparse and can therefore be saved using a MATLAB sparse format, with verylimited demand on the memory.Now, in (1) both W L p and t ( m ) depend on the model parameters m . Thus, the objectivefunction (1) is nonlinear with respect to m . We use a simple iterative strategy to convert P ( α,λ ) ( m ) to a linear form, in which to linearize the cross-gradient constraint, the first orderTaylor expansion is used (Gallardo & Meju 2003; Gallardo & Meju 2004; Tryggvason & Linde2006; Gallardo 2007; Fregoso & Gallardo 2009). First, suppose that the superscript (cid:96) appliedto any variable indicates the value of that variable at iteration (cid:96) , so that m ( (cid:96) ) is the estimateof the model parameters at iteration (cid:96) . Then, we suppose that m (1) = m apr and rewrite theobjective function (1) as P ( α,λ ) ( m ) = (cid:107) W d ( d obs − G m ) (cid:107) + α (cid:107) W ( (cid:96) − D ( m − m ( (cid:96) − ) (cid:107) (8)+ λ (cid:107) t ( (cid:96) − + B ( (cid:96) − ( m − m ( (cid:96) − ) (cid:107) , (cid:96) = 2 , , ... Note that W ( (cid:96) − indicates W estimated at iteration (cid:96) − W L p as given in (4) or (6). Here, t ( (cid:96) − and B ( (cid:96) − = ( ∇ m t ( (cid:96) − ) are the cross-gradient and theJacobian matrix of the discrete approximation for the cross-gradient function, respectively,evaluated at m ( (cid:96) − , consistent with the linear Taylor expansion for t around m ( (cid:96) − . Theformulae used are given in Appendix A. Taking ∇ m P ( α,λ ) ( m ) = defines the update m ( (cid:96) ) asthe solution of − G T W T d W d ( d obs − G m ) + α D T ( W ( (cid:96) − ) T W ( (cid:96) − D ( m − m ( (cid:96) − )+ (9) λ (cid:16) ( B ( (cid:96) − ) T (cid:110) t ( (cid:96) − + B ( (cid:96) − ( m − m ( (cid:96) − ) (cid:111)(cid:17) = 0 . S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Then, after some algebraic manipulation, the desired update for m is given by m ( (cid:96) ) = m ( (cid:96) − + ∆ m ( (cid:96) − , (10)where ∆ m ( (cid:96) − is the solution of the equation (cid:16) G T W T d W d G + α D T ( W ( (cid:96) − ) T W ( (cid:96) − D + λ ( B ( (cid:96) − ) T B ( (cid:96) − (cid:17) ∆ m ( (cid:96) − = (cid:16) G T W T d W d ( d obs − G m ( (cid:96) − ) − λ ( B ( (cid:96) − ) T t ( (cid:96) − (cid:17) . (11)Equivalently, defining E ( (cid:96) ) = (cid:16) G T W T d W d G + α D T ( W ( (cid:96) − ) T W ( (cid:96) − D + λ ( B ( (cid:96) − ) T B ( (cid:96) − (cid:17) , (12)and f ( (cid:96) ) = (cid:16) G T W T d W d ( d obs − G m ( (cid:96) − ) − λ ( B ( (cid:96) − ) T t ( (cid:96) − (cid:17) , (13)then ∆ m ( (cid:96) − given by (11) solves the linear system E ( (cid:96) ) ∆ m ( (cid:96) − = f ( (cid:96) ) . Numerically theCG algorithm can be used to solve (11) and a line search can also be included in the update(10), replacing ∆ m ( (cid:96) − by γ ( (cid:96) − ∆ m ( (cid:96) − where 0 < γ ( (cid:96) − is chosen to speed convergence,see for example Fournier & Oldenburg (2019). We should note that at each iteration of thealgorithm lower and upper bounds on density and susceptibility are imposed. During theinversion process if an estimated physical property falls outside the specified bounds, it willbe returned back to the nearest bound. Further, to test the convergence of the solution ateach iteration (cid:96) , we calculate a χ measure for the respective data fit term at each iteration,( χ ) ( (cid:96) ) = (cid:107) W d i ( d obsi − G i m ( (cid:96) )i ) (cid:107) , i = 1 , . (14)The iteration will terminate at convergence only when ( χ ) ( (cid:96) ) ≤ m + √ m , for both i = 1,and i = 2. Otherwise, the iteration is allowed to proceed to a maximum number of iterationsMAXIT. The steps of the joint inversion algorithm are summarized in Algorithm 1.In (1) the important parameters α and λ are the regularization parameters which giverelative weights to the stabilizer and the cross-gradient term, respectively. To be more precise,we define α as block diag ( α I n , α I n ) ∈ R n × n , where α and α are the relative weightsfor the gravity and magnetic terms, respectively. We should note that, although α is a diag-onal matrix and can be used inside the stabilizer, we prefer to put α outside in order for theformulation to be consistent with the conventional Tikhonov objective functional. Weightingparameters α , α , and λ have an important effect on the estimated solution. Thus, they needto be determined carefully. But application of an automatic parameter-choice method for de-termining α , α , and λ is difficult, or potentially impossible, and is outside the scope of this eneralized L p -norm joint inversion α and α . In subsequent iterations the parameters are reduced slowly dependent on parameters q and q , respectively, using α ( (cid:96) )1 = α ( (cid:96) − q and α ( (cid:96) )2 = α ( (cid:96) − q , where q and q are smallnumbers, 0 (cid:28) q , q <
1. The process continues until the predicted data of one of the recon-structed models satisfies the observed data at the noise level. For that data set, the relevantparameter is then kept fixed during the following iterations. The parameter λ is held fixedin the implementation, although it is quite feasible that it is also iteration dependent. Theamount of structural similarity obtained through the joint inversion algorithm can be adjustedusing different choices of λ . The impact of selecting the parameters α , α , (equivalently q and q ), and λ is demonstrated in Section 3. The validity of the presented joint inversion algorithm is evaluated for synthetic examplesand the results are compared with those of the separate gravity and magnetic inversions.We select a model that consists of two dikes, one is a small vertical dike and the other isa larger dipping dike. In the first study, we suppose that the targets have both a densitycontrast and a susceptibility distribution. The density contrast and the susceptibility of thetargets are selected as ρ = 0 . − and κ = 0 .
06 (SI unit), respectively, embeddedin a homogeneous non-susceptible background. For the second study, we assume that thedipping dike has a density and susceptibility distribution, but that the vertical dike is notmagnetic. For all inversions, the bound constraints 0 = ρ min ≤ m ≤ ρ max = 0 .
6, gr cm − , and0 = κ min ≤ m ≤ κ max = 0 .
06, SI unit, are imposed. To generate the total field anomaly, theintensity of the geomagnetic field, the inclination, and the declination are selected as 50000 nT,45 ◦ and 45 ◦ , respectively. To simulate realistic noise-contaminated observed data, we supposethat Gaussian noise with zero mean and standard deviations ( τ | ( d exacti ) j | + τ max( | d exacti | )), j = 1 ...n , are added to each exact true measurement j , ( d exacti ) j . The parameters pairs ( τ , τ )are selected as (0 . , . . , .
01) for gravity and magnetic data, respectively. Thesestandard deviations are selected such that the signal to noise ratios (SNR), as given bySNR i = 20 log ( (cid:107) d exacti (cid:107) (cid:107) d obsi − d exacti (cid:107) ) , i = 1 , , (15)2 S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Algorithm 1
Generalized L p -norm joint inversion of gravity and magnetic data. Input: d obs , m apr , G , W d , ( W h ) , ( W h ) , D x , D y , D z , (cid:15) , ρ min , ρ max , κ min , κ max , MAXIT, β , β , α (1)1 , α (1)2 , q , q and λ . Calculate ( W depth ) and ( W depth ) as determined by β and β , respectively. Set W = ( W depth ) ( W h ) and W = ( W depth ) ( W h ) . If D is the identity form W = block diag ( W , W ). Otherwise set W =block diag ( W , W , W , W , W , W ). Initialize m (1) = m apr , ( W L p ) (1)1 = I , ( W L p ) (1)2 = I , and (cid:96) = 1. while Not converged, noise level not satisfied, and (cid:96) <
MAXIT do (cid:96) = (cid:96) + 1. Compute t ( (cid:96) − and B ( (cid:96) − , as given in Appendix A. Use CG to solve E ( (cid:96) ) ∆ m ( (cid:96) − = f ( (cid:96) ) for m ( (cid:96) − , defined by (12) and (13). Update m ( (cid:96) ) = m ( (cid:96) − + ∆ m ( (cid:96) − . Impose constraints on m ( (cid:96) ) to force ρ min ≤ m ( (cid:96) )1 ≤ ρ max and κ min ≤ m ( (cid:96) )2 ≤ κ max . Test convergence criteria, (14), for χ and χ . Exit loop if both satisfied. Set α ( (cid:96) )1 = α ( (cid:96) − q and α ( (cid:96) )2 = α ( (cid:96) − q . Update α ( (cid:96) ) . Determine ( W L p ) ( (cid:96) )1 and ( W L p ) ( (cid:96) )2 , dependent on D and p , equations (4) or (6). Set W ( (cid:96) ) = block diag (cid:16) W ( (cid:96) )1 , W ( (cid:96) )2 (cid:17) when D is the identity. Otherwise set W ( (cid:96) ) =block diag (cid:16) W ( (cid:96) )1 , W ( (cid:96) )1 , W ( (cid:96) )1 , W ( (cid:96) )2 , W ( (cid:96) )2 , W ( (cid:96) )2 (cid:17) end whileOutput: Solution ρ = m ( (cid:96) )1 , κ = m ( (cid:96) )2 , IT = (cid:96) .become nearly the same. The resulting SNRs are indicated in the captions of the respectivefigures associated with the results for the given data set. In all simulations we use β = 0 . β = 1 . W depth ) and ( W depth ) , respectively. The maximum number of iterationsof the algorithm is selected as MAXIT = 100 for all inversions. Further, in the following, wepresent the results for two different stabilizers; the L -norm of the model parameters imposedusing D = I and p = 1 in (4), and the minimum structure stabilizer, imposed using p = 2and D (cid:54) = I in (6). This then demonstrates the validity of the algorithm for both sparse andsmooth reconstructions of the subsurface targets. Focusing parameter (cid:15) is held fixed in allinversions, (cid:15) = 1 e − . A summary of the parameters chosen for the simulations discussed inSections 3.1-3.2 are provided in Table 1. eneralized L p -norm joint inversion Model Case Figure p D α (1)1 q α (1)2 q λ MAXIT W h Results1 1 1 1 I ,
000 0 . ,
000 0 .
95 0 100 - 31 2 1 1 I ,
000 0 . ,
000 0 .
95 10
100 - 61 2 1 1 I ,
000 0 . ,
000 0 .
95 10
100 - 91 3 1 1 I ,
000 0 . ,
000 0 .
95 10
100 Yes 121 4 1 2 D ,
000 0 . ,
000 0 .
95 10
100 - 132 1 16 1 I ,
000 0 . ,
000 0 .
95 10
100 - 182 2 16 2 D ,
000 0 . ,
000 0 .
95 10
100 - 19
Table 1.
Summary of the parameters chosen for each of the synthetic simulations that are discussedin Sections 3.1 and 3.2. The true cross section of the models are provided in the figure indicated by thenumber in the column with heading Figure. The obtained cross sections by the inversion are illustratedin the figure with number in the column with heading Results.
The first model is illustrated in Fig. 1. The top depths of the bodies are located at 50 m. Thedipping dike is extended to 300 m, but the maximum depth of the small vertical dike is 200 m.The data on the surface are generated on a grid with 25 ×
20 = 500 points and grid spacing50 m. The noise-contaminated gravity and magnetic data are illustrated in Figs. 2(a) and2(b), respectively. To perform the inversion, the subsurface volume is discretized into 4000prisms of sizes 50 m in each dimension. The initial regularization parameters are selected as α (1)1 = 20 ,
000 and α (1)2 = 50 , G and G , have different spectral properties and,therefore, the regularization parameter should be much larger for the inversion of magneticdata as compared to that used for the inversion of gravity data. Thus, it is appropriate tocontrol the speed of convergence for each model with parameters q = 0 . q = 0 .
95, thatare different. -norm inversion We first implement the algorithm using the L -norm stabilizer but without using the cross-gradient constraint. This is easy to do by selecting λ = 0. Hence the algorithm proceeds exactlyas in Algorithm 1, with the same termination criteria, but without the cross-gradient term.The inversion is initiated using m apr = . After IT = 61 iterations the convergence criteria, χ and χ , are satisfied and the inversion terminates. In this simulation the χ termination is4 S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 1.
Cross-section of the synthetic model that consists of a vertical and a dipping dike. (a)Density distribution; (b) Susceptibility distribution. reached at iteration 43. The susceptibility model is recovered more quickly than the densitymodel, which requires additional iterations until both termination criteria are satisfied.The reconstructed density and susceptibility models are illustrated in Figs. 3(a) and 3(b),respectively. They are in good agreement with the original models, and are acceptable. Theresults are very close to those that can be obtained using another single L -norm algorithm, seeVatankhah et al. (2019b), in which the unbiased predictive risk estimator is used as parameter-choice method and the singular value decomposition (SVD) is used computationally. It is clearthat sharp and focused images of the subsurface are obtained, but the recovered susceptibilitymodel shows more extension in depth for both targets. Further, in the density model, theextension of the vertical dike is underestimated. For comparison, we present the data producedby the models in Figs. 4(a) and 4(b), respectively. In addition, the the progression of the datamisfit, the regularization term, and the regularization parameter, with iteration (cid:96) , for bothgravity and magnetic problems, are illustrated in Fig. 5. (a) Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -2000200400600 nT Figure 2.
Noise contaminated anomaly produced by the model shown in Fig. 1. (a) Vertical componentof gravity; (b) Total magnetic field. The SNR for gravity and magnetic data, respectively, are 25 . . eneralized L p -norm joint inversion Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 3.
Cross-section of the reconstructed model using individual inversions, Case 1. (a) Densitydistribution; (b) Susceptibility distribution. (a)
Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -2000200400600 nT Figure 4.
The data produced by the models shown in Fig. 3, Case 1. (a) Vertical component of gravity;(b) Total magnetic field.
10 20 30 40 50 60
Iteration number -10 -5 (a) Data misfitStabilizer
10 20 30 40 50 60
Iteration number -10 -5 (b) Data misfitStabilizer Figure 5.
The progression of the data misfit, the regularization term, and the regularization parameter,with iteration (cid:96) , for the models shown in Fig. 3, Model 1 and Case 1. (a) Gravity inversion; (b) Magneticinversion. S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 6.
Reconstructed models using joint inversion with λ = 10 and with the L -norm of themodel parameters as the stabilizer, Model 1 and Case 2. (a) Density distribution; (b) Susceptibilitydistribution. -norm inversion with cross-gradient constraint We now implement the inversion algorithm using the cross-gradient constraint with λ =10 . This selection is based on an analysis of entries of matrix B , which are very small.If we want to give enough weight to the cross-gradient term, it is necessary to use a largevalue for λ . We will also show, on the other hand, that if λ is too large the results canbe unsatisfactory. All other parameters are selected as for the simulation for Model 1 andCase 1, in Section 3.1.1. The inversion terminates at iteration IT = 59, two iterations lessthan required for the presented individual inversions. The results are presented in Figs. 6, 7,and 8, respectively. Now, the reconstructed density and susceptibility models are quite similar,and are close to the original models. Indeed, the results indicate that the application of thecross-gradient function significantly improves the quality of the solutions. Both targets are inbetter agreement with the original models.We next perform the same simulation but select the weighting on the cross-gradient to be (a) Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -2000200400600 nT Figure 7.
The data produced by the models shown in Fig. 6, Model 1 and Case 2. (a) Verticalcomponent of gravity; (b) Total magnetic field. eneralized L p -norm joint inversion
10 20 30 40 50
Iteration number -10 -5 (a) Data misfitStabilizer
10 20 30 40 50
Iteration number -10 -5 (b) Data misfitStabilizer Figure 8.
The progression of the data misfit, the regularization term, and the regularization parameter,with iteration (cid:96) , for the models shown in Fig. 6, Model 1 and Case 2. (a) Gravity; and (b) Magnetic. very large, with λ = 10 . In this case, the inversion terminates at MAXIT = 100, withoutsatisfying the noise levels. That means that it is not possible to satisfy the data misfit criteriawith a large λ . The results are presented in Figs. 9, 10, and 11. The reconstructed modelsare not at all consistent with the original models, either with respect to the shape or tothe maximum values of the physical properties, especially for reconstruction of the density.Although, the main reconstructed bodies are quite similar to each other, as expected due tothe strong requirement imposed by the use of the cross-gradient constraint, some additionalunrealistic structures appear in the density model. Clearly the selection of λ is very important.But, this is not a difficult task. It is sufficient to consider the entries of B , or to run thealgorithm once, to determine a suitable λ . Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 9.
Reconstructed models using joint inversion with λ = 10 and with the L -norm of themodel parameters as the stabilizer, Model 1 and Case 2. (a) Density distribution; (b) Susceptibilitydistribution. S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi (a)
Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -200-1000100200300400500 nT Figure 10.
The data produced by the models shown in Fig. 9, Model 1 and Case 2. (a) Verticalcomponent of gravity; (b) Total magnetic field. -norm inversion with cross-gradient constraint and hard constraintmatrix The algorithm is implemented again with λ = 10 , but now, based on available information,we suppose that the physical properties for one prism in the vertical dike, located in theupper left side, are known. These known values are communicated to the algorithm via m apr and by setting the corresponding entries of W h to 100. The inversion terminates at iterationIT = 59 < MAXIT. The reconstructed models are presented in Figs. 12(a) and 12(b). Theresults are even better than those obtained by the previous reconstruction, as demonstratedby the model illustrated in Fig. 6. Here, the known model parameters are kept fixed during
20 40 60 80 100
Iteration number -15 -10 -5 (a) Data misfitStabilizer
20 40 60 80 100
Iteration number -15 -10 -5 (b) Data misfitStabilizer Figure 11.
The progression of the data misfit, the regularization term, and the regularization pa-rameter, with iteration (cid:96) , for the models shown in Fig. 9, Model 1 and Case 2. (a) Gravity; and (b)Magnetic. eneralized L p -norm joint inversion Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 12.
Reconstructed models using joint inversion with λ = 10 , with the L -norm of the modelparameters as the stabilizer, and with the application of the hard constraint matrix, Model 1 andCase 3. (a) Density distribution; (b) Susceptibility distribution. the iterations. The results demonstrate how the incorporation of available information forsome model parameters can increase the reliability of the model obtained using the inversealgorithm. For this example, we implement the joint inversion algorithm using the minimum structurestabilizer and initial regularization parameters, α (1)1 = 200 ,
000 and α (1)2 = 500 , -norm implementations discussed for Model 1 withCases 1 to 3, in Sections 3.1.1-3.1.3. This increase in the regularization parameters is dueto the significant change in the minimization function when changing the stabilizer from theL -norm to the L -norm. Moreover, we also use λ = 10 in order to increase the weightingon the cross-gradient constraint, consistently with increasing the weight on the stabilizer. Allother parameters remain the same as selected for Cases 1 to 3. It is also important to notethat for this simulation, even though α and α decrease as determined by q and q , the datamisfit does not necessarily decrease monotonically with the iterations. Because this increasein the data misfit can occur, the algorithm is modified to hold α or α fixed for any iterationin which an increase in the calculated data misfit occurs. Moreover, the related model updatefor that iteration is rejected. Specifically, this means that the step is repeated without thedecrease in the regularization parameter. This experience demonstrates that it is importantto determine a reliable automatic parameter-choice strategy during the inversion, which is acomplicated problem when three parameters are involved. The approach adopted here, whichis quite simple but probably not optimal, leads to acceptable solutions, as compared withthose presented by Fregoso & Gallardo (2009). The inversion terminates at MAXIT = 100.Equivalently, this means that the χ tests on the noise level are not satisfied for both datasets. The results are presented in Figs. 13, 14, and 15. As expected from the use of the0 S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 13.
Reconstructed models using joint inversion with λ = 10 and with the L -norm of thegradient of the model parameters as the stabilizer, Model 1 and Case 4. (a) Density distribution; (b)Susceptibility distribution. minimum structure constraint, the reconstructed models are smooth. Moreover, as imposedby the use of the cross-gradient constraint, the models have a similar structure. They indicatea large dipping dike along with a small vertical structure in the subsurface. We note that theminimum structure inversion has its own advantages, which should not be disregarded, andthat make the algorithm a safe strategy for the reconstruction of low-frequency subsurfacestructures. The progression of the data misfit, the regularization term, and the regularizationparameters, with iteration (cid:96) , are presented in Fig. 15. These plots demonstrate that it isonly at the initial iterations where there are significant changes in the parameters betweensuccessive iterative steps. Once the parameters are effectively fixed, the changes are very small,which suggests that the algorithm can be safely terminated for a smaller value of MAXIT.For direct comparison with the presented results given for Model 1 with Cases 1 to 3 we keepMAXIT = 100. (a) Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -2000200400600 nT Figure 14.
The data produced by the models shown in Fig. 13, Model 1 and Case 4. (a) Verticalcomponent of gravity; (b) Total magnetic field. eneralized L p -norm joint inversion
20 40 60 80 100
Iteration number -5 (a) Data misfitStabilizer
20 40 60 80 100
Iteration number -5 (b) Data misfitStabilizer Figure 15.
The progression of the data misfit, the regularization term, and the regularization pa-rameter, with iteration (cid:96) , for the models shown in Fig. 13, Model 1 and Case 4. (a) Gravity; and (b)Magnetic.
For the second model we assume that the dipping dike has both a density and a susceptibilitydistribution, but that the vertical dike is not magnetic and has, therefore, only a densitydistribution. Figs. 16(a) and 16(b) illustrate the models. Here, the aim is to test the jointinversion algorithm, when one model has a structure in a location where the other model doesnot. The data produced by the models are presented in Figs. 17(a) and 17(b), respectively. -norm inversion with cross-gradient constraint The joint inversion algorithm is implemented using the L -norm stabilizer and with the pa-rameters the same as were selected for the Case 2 simulation of Model 1 in Section 3.1.2. Theinversion terminates at IT = 58, one iteration less than its counterpart in Section 3.1.2. Thereconstructed models are shown in Figs. 18(a) and 18(b). The dipping dikes, in both densityand susceptibility distributions, are reconstructed very well, and they are completely similar. Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 16.
Cross-section of the synthetic model 2. (a) Density distribution; (b) Susceptibility distri-bution. Here, it is assumed that the vertical dike is not magnetic. S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi (a)
Easting (m) N o r t h i ng ( m ) mGal (b) Easting (m) N o r t h i ng ( m ) -2000200400600 nT Figure 17.
Noise contaminated anomaly produced by the model shown in Fig. 16. (a) Vertical com-ponent of gravity; (b) Total magnetic field. The SNR for gravity and magnetic data, respectively, are25 . . The density distribution of the vertical dike is almost reconstructed and no susceptibility dis-tribution is obtained. This confirms, as noted by the theory of the cross-gradient constraint,that it is possible to reconstruct models which do not share all structures. Only the structureswhich are supported by the data will be similar.
For the final validation of the algorithm, we implemented the joint inversion algorithm usingthe minimum structure stabilizer for the data of Fig 17. All the parameters are selected asfor the simulation for Model 1 with Case 4 in Section 3.1.4. The algorithm terminated atthe maximum iteration, IT = 100. The reconstructed models are presented in Figs. 19(a)and 19(b). As for the simulation with the L -norm of the model parameters in Case 1, thesusceptibility model exhibits none of the structure of the vertical dike. Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 18.
Reconstructed models using joint inversion with λ = 10 and with the L -norm of themodel parameters as the stabilizer, Model 2 and Case 1. (a) Density distribution; (b) Susceptibilitydistribution. eneralized L p -norm joint inversion Easting (m) D e p t h ( m ) (a) g/cm Easting (m) D e p t h ( m ) (b) SI Figure 19.
Reconstructed models using joint inversion with λ = 10 and with the L -norm of thegradient of the model parameters as the stabilizer, Model 2 and Case 2. (a) Density distribution; (b)Susceptibility distribution. A unified framework for the incorporation of the L p -norm constraint in an algorithm forjoint inversion of gravity and magnetic data, in which the cross-gradient constraint providesthe link between the two models, has been developed. This unifying framework shows howit is possible to incorporate all well-known and widely used stabilizers, that are used forpotential field inversion, within a joint inversion algorithm with the cross-gradient constraint.By suitable choices of the parameter p and the weighting matrix, that define the L p -normconstraint, it is possible to reconstruct a subsurface target exhibiting smooth, sparse, or blockycharacteristics. The global objective function for the joint inversion consists of a data misfitterm, a general form for the stabilizer, and the cross-gradient constraint. Their contributionsto the global objective function are obtained using three different regularization parameters. Asimple iterative strategy is used to convert the global non-linear objective function to a linearform at each iteration, and the regularization weights can be adjusted at each iteration. Depthweighting and hard constraint matrices are also used in the presented inversion algorithm.These make it possible to weight prisms at depth, and to include the known values of someprisms in the reconstructed model. Bound constraints on the model parameters may also beimposed at each iteration.Results presented for two synthetic three-dimensional models illustrate the performance ofthe developed algorithm. These results indicate that, when suitable regularization parameterscan be estimated, the joint inversion algorithm yields suitable reconstructions of the subsurfacestructures. These reconstructions are improved in comparison with reconstructions obtainedusing independent gravity and magnetic inversions. The structures of the subsurface targets,for both density and susceptibility distributions, are quite similar and are close to the originalmodels. A simple but practical strategy for the estimation of the regularization parametersis provided, by which large values are used at the initial step of the iteration, with a gradual4 S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi decrease in subsequent iterations, dependent on selected scaling parameters for each of theimposed gravity and magnetic constraint terms. The weight on the cross-gradient linkageconstraint is chosen to balance to the three regularization terms. The results shows that thisstrategy is effective, particularly given the lack of any known robust methods for automaticallyestimating these parameters. The latter is a topic for our future study, as is the developmentof an improved implementation for the practical and efficient solution of three-dimensionallarge-scale problems and its application for real data.
ACKNOWLEDGMENTS
Rosemary Renaut acknowledges the support of NSF grant DMS 1913136: “Approximate Sin-gular Value Expansions and Solutions of Ill-Posed Problems”.
REFERENCES
Ajo-Franklin, J. B., Minsley, B. J. & Daley, T. M., 2007. Applying compactness constraints to differ-ential traveltime tomography,
Geophysics ,
72 (4) , R67-R75.Barbosa, V. C. F. & Silva, J. B. C., 1994. Generalized compact gravity inversion,
Geophysics ,
59 (1) ,57-68, doi: 10.1190/1.1443534 .Bertete-Aguirre, H., Cherkaev, E. & Oristaglio, M., 2002. Non-smooth gravity problem with to-tal variation penalization functional,
Geophysical Journal International ,
149 (2) , 499-507, doi:10.1046/j.1365-246X.2002.01664.x.Blakely, R. J., 1996.
Potential Theory in Gravity & Magnetic Applications , Cambridge University Press,Cambridge, United Kingdom.Boulanger, O. & Chouteau, M., 2001. Constraint in 3D gravity inversion,
Geophysical Prospecting ,
49 (2) , 265-280, doi:10.1046/j.1365-2478.2001.00254.x.Constable, S. C., Parker, R. L. & Constable, C. G., 1987. Occam’s inversion: A practical algorithm forgenerating smooth models from electromagnetic sounding data,
Geophysics ,
52 (3) , 289-300.Fiandaca, G., Doetsch, J., Vignoli, G. & Auken, E, 2015. Generalized focusing of time-lapse changeswith applications to direct current and time-domain induced polarization inversions,
GeophysicalJournal International ,
203 (2) , 1101-1112, doi: 10.1093/gji/ggv350.Fournier, D. & Oldenburg, D. W., 2019. Inversion using spatially variable mixed (cid:96) p -norms, GeophysicalJournal International ,
218 (1) , 268-282, doi: 10.1093/gji/ggz156.Farquharson, C. G. & Oldenburg, D. W., 1998. Nonlinear inversion using general measure of data misfitand model structure,
Geophysical Journal International ,
134 (1) , 213-227, doi: 10.1046/j.1365-246x.1998.00555.x.Farquharson, C. G. & Oldenburg, D. W., 2004. A comparison of automatic techniques for estimating eneralized L p -norm joint inversion the regularization parameter in non-linear inverse problems, Geophysical Journal International ,
156 (3) , 411-425, doi: 10.1111/j.1365-246X.2004.02190.x.Fregoso, E. & Gallardo, L. A., 2009. Cross-gradients joint 3D inversion with applications to gravityand magnetic data,
Geophysics ,
74 (4) , L31-L42.Fregoso, E., Gallardo, L. A. & Garcia-Abdeslem, J., 2015. Structural joint inversion coupled with Eulerdeconvolution of isolated gravity and magnetic anomalies,
Geophysics ,
80 (2) , G67-G79.Gallardo, L. A., 2007. Multiple cross-gradient joint inversion for geospectral imaging,
GeophysicalResearch Letters , , L19301, doi: 10.1029/2007GL030409.Gallardo, L. A. & Meju, M. A., 2003. Characterization of heterogeneous near-surface materials byjoint 2D inversion of DC resistivity and seismic data, Geophys. Res. Lett. ,
30 (13) , 1658,doi:10.1029/2003GL017370.Gallardo, L. A. & Meju, M. A., 2004. Joint two-dimensional DC resistivity and seismic travel timeinversion with cross-gradients constraints,
Journal of Geophysical Research , , B03311.Gross, L., 2019. Weighted cross-gradient function for joint inversion with the application to regional3-D gravity and magnetic anomalies, Geophysical Journal International ,
217 (3) , 2035-2046.doi:10.1093/gji/ggz134.Ha´az, I. B., 1953. Relations between the potential of the attraction of the mass contained in a fi-nite rectangular prism and its first and second derivatives,
Geofizikai Kozlemenyek , , (inHungarian).Haber, E. & Oldenburg, D. W., 1997. Joint inversion: a structural approach, Inverse Problems , ,6377.Haber, E. & Gazit, M. H., 2013. Model Fusion and Joint Inversion, Surv. Geophys. , , 675-695.Jorgensen, M. & Zhdanov, M. S., 2019. Imaging Yellowstone magmatic system by the joint Gramianinversion of gravity and magnetotelluric data, Physics of the Earth and Planetary Interiors , ,12-20. doi: 10.1016/j.pepi.2019.05.003.Last, B. J. & Kubik, K., 1983. Compact gravity inversion, Geophysics ,
48 (6) , 713-721, doi:10.1190/1.1441501.Leli´evre, P. G. & Oldenburg, D. W., 2009. A comprehensive study of including structural orientationinformation in geophysical inversions,
Geophysical Journal International ,
178 (2) , 623-637, doi:10.1111/j.1365-246X.2009.04188.x.Leli´evre, P. G. & Farquharson, C. G., 2013. Gradient and smoothness regularization operators forgeophysical inversion on unstructured meshes,
Geophysical Journal International ,
195 (1) , 330-341, doi: 10.1093/gji/ggt255.Li, Y. & Oldenburg, D. W., 1996. 3-D inversion of magnetic data,
Geophysics ,
61 (2) , 394-408.Li, Y. & Oldenburg, D. W., 1998. 3-D inversion of gravity data,
Geophysics ,
63 (1) , 109-119.Li, Y. & Oldenburg, D. W., 2000. Incorporating geologic dip information into geophysical inversions,
Geophysics ,
65 (1) , 148-157, doi: 10.1190/1.1444705. S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi
Lin, W. & Zhdanov, M. S., 2018. Joint multinary inversion of gravity and magnetic data using Gramianconstraints,
Geophysical Journal International ,
215 (3) ,1540-1577. doi: 10.1093/gji/ggy351.Moorkamp, M., Heincke, B., Jegen, M., Roberts A. W. & Hobbs, R. W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data,
Geophysical Journal International ,
184 (1) , 477-493.Moorkamp, M., Roberts, A. W., Jegen, M., Heincke, B., & Hobbs, R. W., 2013. Verification of velocity-resistivity relationships derived from structural joint inversion with borehole data,
GeophysicalResearch Letters , , 3596-3601.Nabighian, M. N., Ander, M. E., Grauch, V. J. S., Hansen, R. O., Lafehr, T. R., Li, Y., Pearson, W. C.,Peirce, J. W., Philips, J. D. & Ruder, M. E., 2005. Historical development of the gravity methodin exploration, Geophysics , , 63-89.Nielsen, L. & Jacobsen, B. H, 2000. Integrated gravity and wide-angle seismic inversion for two-dimensional crustal modelling, Geophysical Journal International ,
140 (1) , 222-232.Pilkington, M., 1997. 3-D magnetic imaging using conjugate gradients,
Geophysics ,
62 (4) , 1132-1142.Portniaguine, O. & Zhdanov, M. S., 1999. Focusing geophysical inversion images,
Geophysics ,
64 (3) ,874-887, doi:10.1190/1.1444596.Rao, D. B. & Babu, N. R., 1991. A rapid method for three-dimensional modeling of magnetic anomalies,
Geophysics ,
56 (11) , 1729-1737. doi:10.1190/1.1442985.Sun, J. & Li, Y., 2014. Adaptive L p inversion for simultaneous recovery of both blocky andsmooth features in a geophysical model, Geophysical Journal International ,
197 (2) , 882-899,doi: 10.1093/gji/ggu067.Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least squarescriterion,
Reviews of Geophysics and Space Physics ,
20 (2) , 219-232.Tryggvason, A & Linde, N., 2006. Local earthquake (LE) tomography with joint inversion for P-and S-wave velocities using structural constraints,
Geophysical Research Letters , , L07303,doi:10.1029/2005GL025485.Vatankhah, S., Ardestani, V. E. & Renaut, R. A., 2014a. Automatic estimation of the regularizationparameter in 2D focusing gravity inversion: application of the method to the Safo manganese minein northwest of Iran, Journal Of Geophysics and Engineering ,
11 (4) , 045001, doi:10.1088/1742-2132/11/4/045001.Vatankhah, S., Renaut, R. A. & Ardestani, V. E., 2014b. Regularization parameter estimation forunderdetermined problems by the χ principle with application to 2D focusing gravity inversion, Inverse Problems , , 085002, doi: 10.1088/0266-5611/30/8/085002.Vatankhah, S., Renaut, R. A. & Ardestani, V. E., 2017. 3-D Projected L inversion of gravity data usingtruncated unbiased predictive risk estimator for regularization parameter estimation, GeophysicalJournal International ,
210 (3) , 1872-1887, doi:10.1093/gji/ggx274.Vatankhah, S., Renaut, R. A. & Ardestani, V. E., 2018a. Total variation regularization of the 3-D eneralized L p -norm joint inversion gravity inverse problem using a randomized generalized singular value decomposition, GeophysicalJournal International ,
213 (1) , 695-705, doi: 10.1093/gji/ggy014.Vatankhah, S., Renaut, R. A. & Ardestani, V. E., 2018b. A fast algorithm for regularized focused3D inversion of gravity data using randomized singular-value decomposition,
Geophysics ,
83 (4) ,G25-G34, doi: 10.1190/geo2017-0386.1.Vatankhah, S., Liu, S., Renaut, R. A., Hu, X. & Baniamerian J., 2019b. Improving the use of therandomized singular value decomposition for the inversion of gravity and magnetic data,
Submitted ,, -, arXiv:1906.11221.Vatankhah, S., Renaut, R. A. & Liu, S., 2019a. A unifying framework for widely-used stabiliza-tion of potential field inverse problems,
Geophysical Prospecting , Under preparation , 0-0, doi:10.1111/1365-2478.12926.Wohlberg, B. & Rodr´ıguez, P. 2007. An iteratively reweighted norm algorithm for minimization oftotal variation functionals,
IEEE Signal Processing Letters ,
14 (12) , 948–951.Zhang, Y. & Wang, Y., 2019.
Three-Dimensional Gravity-Magnetic Cross-Gradient Joint InversionBased on Structural Coupling and a Fast Gradient Method , Journal of Computational Mathematics ,
37 (6) , 758-777. doi:10.4208/jcm.1905-m2018-0240.Zhdanov, M. S., 2015.
Inverse Theory and Applications in Geophysics , ElsevierZhdanov, M. S. & Tolstaya, E., 2004. Minimum support nonlinear parameterization in the solu-tion of 3D magnetotelluric inverse problem,
Inverse Problems , , 937-952, doi: 10.1088/0266-5611/20/3/017.Zhdanov, M. S., Gribenko, A. & Wilson, G., 2012. Generalized joint inversion of multi-modal geophysical data using Gramian constraints, Geophysical Research Letters , , L09301,doi:10.1029/2012GL051233. APPENDIX A: CROSS-GRADIENT FORMULATION
The components of the cross-gradient function (2) are given by t x ( m ( x, y, z ) , m ( x, y, z )) = ∂ m ∂y ∂ m ∂z − ∂ m ∂z ∂ m ∂y ∈ R n (A.1) t y ( m ( x, y, z ) , m ( x, y, z )) = ∂ m ∂z ∂ m ∂x − ∂ m ∂x ∂ m ∂z ∈ R n (A.2) t z ( m ( x, y, z ) , m ( x, y, z )) = ∂ m ∂x ∂ m ∂y − ∂ m ∂y ∂ m ∂x ∈ R n , (A.3)yielding t ( m ( x, y, z )) = block stack ( t x , t y , t z ) ∈ R n . As illustrated in Fig. A1, the subsurfaceis commonly divided into right rectangular prisms. Here, we suppose all prisms have samedimensions and that m ( i, j, k ) represents the value of the current estimate for m at ( x, y, z ) =( i ∆ x, j ∆ y, k ∆ z ) where i ≥ j ≥ k ≥
0, with the origin, m (0 , , S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi at the given grid point of the volume. We use forward difference operators for each of thederivatives in (A.1) to give . t x ( i, j, k ) = (cid:18) m ( i, j + 1 , k ) − m ( i, j, k )∆ y (cid:19) (cid:18) m ( i, j, k + 1) − m ( i, j, k )∆ z (cid:19) − (cid:18) m ( i, j, k + 1) − m ( i, j, k )∆ z (cid:19) (cid:18) m ( i, j + 1 , k ) − m ( i, j, k )∆ y (cid:19) t y ( i, j, k ) = (cid:18) m ( i, j, k + 1) − m ( i, j, k )∆ z (cid:19) (cid:18) m ( i + 1 , j, k ) − m ( i, j, k )∆ x (cid:19) − (cid:18) m ( i + 1 , j, k ) − m ( i, j, k )∆ x (cid:19) (cid:18) m ( i, j, k + 1) − m ( i, j, k )∆ z (cid:19) t z ( i, j, k ) = (cid:18) m ( i + 1 , j, k ) − m ( i, j, k )∆ x (cid:19) (cid:18) m ( i, j + 1 , k ) − m ( i, j, k )∆ y (cid:19) − (cid:18) m ( i, j + 1 , k ) − m ( i, j, k )∆ y (cid:19) (cid:18) m ( i + 1 , j, k ) − m ( i, j, k )∆ x (cid:19) . These simplify as t x ( i, j, k ) = 1∆ y ∆ z (cid:18) m ( i, j, k ) (cid:16) m ( i, j + 1 , k ) − m ( i, j, k + 1) (cid:17) + m ( i, j + 1 , k ) (cid:16) m ( i, j, k + 1) − m ( i, j, k ) (cid:17) + m ( i, j, k + 1) (cid:16) m ( i, j, k ) − m ( i, j + 1 , k ) (cid:17)(cid:19) (A.4) t y ( i, j, k ) = 1∆ z ∆ x (cid:18) m ( i, j, k ) (cid:16) m ( i, j, k + 1) − m ( i + 1 , j, k ) (cid:17) + m ( i, j, k + 1) (cid:16) m ( i + 1 , j, k ) − m ( i, j, k ) (cid:17) + m ( i + 1 , j, k ) (cid:16) m ( i, j, k ) − m ( i, j, k + 1) (cid:17)(cid:19) (A.5) t z ( i, j, k ) = 1∆ x ∆ y (cid:18) m ( i, j, k ) (cid:16) m ( i + 1 , j, k ) − m ( i, j + 1 , k ) (cid:17) + m ( i + 1 , j, k ) (cid:16) m ( i, j + 1 , k ) − m ( i, j, k ) (cid:17) + m ( i, j + 1 , k ) (cid:16) m ( i, j, k ) − m ( i + 1 , j, k ) (cid:17)(cid:19) . (A.6)The Jacobian matrix for the cross-gradient function is given by B = B x B x B y B y B z B z = ( B , B ) ∈ R n × n (A.7)where B i x = ∂ t x ∂ m i , B i y = ∂ t y ∂ m i , B i z = ∂ t z ∂ m i , i = 1 , . (A.8)We illustrate the discrete derivative for t x with respect to m and m , first noting from (A.4)that t x is a nonlinear combination of m ( i, j, k ) , m ( i, j + 1 , k ) , m ( i, j, k + 1) and m ( i, j, k ) , m ( i, j + 1 , k ) , m ( i, j, k + 1) . eneralized L p -norm joint inversion m and m that are nonzero; in total there are only six non zero column entries for each row of the firstrow block matrix ( B x , B x ). Specifically, we only have for row ijk the column entries pqr asgiven by( B x ) ijk,pqr = 1∆ y ∆ z m ( i, j + 1 , k ) − m ( i, j, k + 1) p = i q = j r = k m ( i, j, k + 1) − m ( i, j, k ) p = i q = j + 1 r = k m ( i, j, k ) − m ( i, j + 1 , k ) p = i q = j r = k + 1(A.9)( B x ) ijk,pqr = 1∆ y ∆ z m ( i, j, k + 1) − m ( i, j + 1 , k ) p = i q = j r = k m ( i, j, k ) − m ( i, j, k + 1) p = i q = j + 1 r = k m ( i, j + 1 , k ) − m ( i, j, k ) p = i q = j r = k + 1(A.10)Here ( B ) ijk,pqr indicates the row associated with grid point ( x i , y y , z k ), and the nonzeroentries are in the relevant columns indexed by pqr . Therefore, the nonzero elements on eachrow are given by( B ) ijk, ··· = 1∆ y ∆ z (cid:18) · · · (cid:16) m ( i, j + 1 , k ) − m ( i, j, k + 1) (cid:17) · · · (cid:16) m ( i, j, k + 1) − m ( i, j, k ) (cid:17) · · · (cid:16) m ( i, j, k ) − m ( i, j + 1 , k ) (cid:17)(cid:19) ( B ) ijk, ··· = 1∆ y ∆ z (cid:18) · · · (cid:16) m ( i, j, k + 1) − m ( i, j + 1 , k ) (cid:17) · · · (cid:16) m ( i, j, k ) − m ( i, j, k + 1) (cid:17) · · · (cid:16) m ( i, j + 1 , k ) − m ( i, j, k ) (cid:17)(cid:19) . These equations are consistent with (A.9) and (A.10). Furthermore, the non zero entries fortwo row block matrices associated with derivatives of t y and t z are obtained similarly from(A.5) and (A.6).0 S. Vatankhah, S. Liu, R. A. Renaut, X. Hu, M. Gharloghi