Decoupled joint inversion with variable splitting: example scheme for magnetotelluric, seismic and gravity data

SUMMARY

We present a general framework for multiphysics joint inversion of any number of geophysical data sets. Its main feature is the use of the variable splitting approach: an auxiliary multiparameter model space is introduced in which minimization of the coupling and stabilizing functionals is carried out. The use of rediscretization and interpolation to map between this auxiliary space and the model spaces allows the coupled models to have completely different parametrizations. Joint inversion is decoupled into the individual inversion and the coupling-regularization subproblems, each of which can be solved by a different optimization algorithm. For each subproblem, the linking term controlling the distance between the model and the corresponding auxiliary variable takes the form of a quadratic regularization with a reference model. As a result, any existing inversion code supporting such regularization can be integrated without modifications into the developed framework. As a concrete example scheme, we consider an application of the framework to 3-D joint inversion of magnetotelluric, seismic refraction and gravity data. We discuss different coupling functionals, mainly those corresponding to the more universal structural constraints: joint total variation, joint minimum support, cross-gradient, one-way cross-gradient and their combinations for a general multimodel case. The use of coupling based on explicit ‘petrophysical’ relationship between the properties is also considered. Performance of the developed framework is studied on three synthetic cases: a time-lapse joint inversion of full-tensor gravity gradiometry and seismic data, a joint inversion of magnetotelluric, seismic and gravity data and a joint inversion for electrical resistivity tomography and audio-magnetotellurics.

1 INTRODUCTION

Joint inversion is an actively developing approach to the integration of multiple geophysical methods that aims at reducing the uncertainties inherent to the geophysical inversions and subsequent geological interpretation. Related approach, when a model obtained with one geophysical method is incorporated in the inversion of another method, is termed ‘model fusion’ (Haber & Holtzman Gazit 2013) or ‘cooperative inversion’ (Moorkamp et al. 2016). Both approaches can be applied to methods sensitive to the same physical property, for example magnetotellurics and electrical resistivity tomography, which are both sensitive to the resistivity (see Vozoff & Jupp 1975), or to methods sensitive to different properties. In the latter case, some coupling between the models has to be assumed (the single-property joint/cooperative inversion can be viewed as a special case with trivial coupling).

Depending on the nature of the available prior information, coupling between the geophysical models can be based on structural similarity (see Crestel et al. 2018, for an overview), explicitly defined common structural parameters (e.g. Juhojuntti & Kamm 2015; Zheglova et al. 2018) or relationships between the physical properties, defined either analytically (e.g. Heincke et al. 2017; Colombo & Rovetta 2018) or using statistical data analysis techniques such as fuzzy clustering (e.g. Lelièvre et al. 2012; Sun & Li 2016) and Gaussian mixture models (Astic et al. 2021). A promising approach to model coupling is based on maximization of various measures of statistical dependence: Haber & Holtzman Gazit (2013) and Mandolesi & Jones (2014) experimented with the mutual information; as an alternative, the use of the variation of information (Moorkamp 2021) and the joint entropy (Zhdanov et al. 2022) have been suggested. Gramian constraints (Zhdanov et al. 2012; Malovichko et al. 2020) allow for efficient maximization of the correlation coefficients between the models or their transforms. Alternatively to the coupling of geophysical models, a joint inverse problem can be solved directly for lithologies (Bosch 1999) or for thermochemical (Khan et al. 2007; Afonso et al. 2013) or petrophysical (Chen et al. 2007; Gao et al. 2012) parameters.

Deterministic inversion methods based on local optimization have well-known limitations: the uncertainty estimation, when available, is restricted to linearized estimates, solutions depend on the choice of the regularization parameters (as a scalarization of a multiobjective optimization problem) and, for nonlinear problems, on the initial model. However, at the present time, for most 3-D inverse problems with a large number of unknowns, stochastic, global optimization methods are not computationally feasible. This situation may change with the advances in techniques like the reduced order modelling (Manassero et al. 2021) and deep learning (Colombo et al. 2021).

In this paper we consider a deterministic 3-D joint inversion framework, with the primary aim of applying it to large magnetotelluric (MT), gravity and seismic first-arrival traveltime (FATT) data sets, on both regional and local scales. During the preparation of this publication, magnetic data inversion was added to the framework. In the future, we plan to add inversion of seismic surface wave dispersion data. Thus, we need a framework to which geophysical methods can be easily added, preferably as existing inversion codes. Model discretization required by a particular method is dictated by the physics (resolution, boundary conditions, etc.) and by the numerical implementation of the forward problem. Although it is possible to use a common grid for the joint inverse problem and individual, problem-specific grids for the forward problems (e.g. see Moorkamp et al. 2011, Heincke et al. 2017), our choice here is to use an individual grid for each inverse problem—this makes it easier to adapt inversion codes and ensures that a model parametrization is capable of entirely explaining each data set. The ‘philosophy’ of our approach is that the priority is given to the individual inversions (i.e. the model parametrization, the forward solver and the optimization method are optimal for each data set in the sense that the same choices would have been made in a ‘separate’ inversion case), while the ‘joint’ part, that is the coupling, is added as a logically independent and easily modifiable block. We show how this can be achieved with a simple technique, often referred to as variable splitting (e.g. Wang et al. 2008): terms of an objective function are decoupled by introducing an auxiliary variable and penalizing the distance between this variable and the model vector. By combining this technique with interpolation between the grids, the joint inverse problem is decoupled into the separate inversion (with standard quadratic regularization) and the coupling (and regularization) subproblems that are solved in an alternating manner.

The paper is organized as follows. First, we derive a general split-variable formulation of joint inversion and discuss its practical implementation as an alternating minimization algorithm. Next, we describe several structure-coupling functionals, which are implemented in our joint inversion framework: the joint total variation (JTV), the joint minimum support (JMS), the cross-gradient (XG) and the one-way cross-gradient (OWXG). We then discuss the main aspects of the application of the proposed scheme to MT, FATT and gravity problems. Finally, performance of the framework is demonstrated on synthetic data: in the first example we consider a time-lapse joint inversion of FATT and gravity gradiometry data using JMS constraints, in the second example we compare performance of several combinations of JTV, XG, OWXG and ‘petrophysical’ constraints in a joint inversion of MT, FATT and gravity data, and in the third example we consider joint inversion of audio-magnetotelluric (AMT) and electrical resistivity tomography (ERT) data combining 2-D finite-element and 3-D finite-difference grids.

2 JOINT INVERSION BY VARIABLE SPLITTING AND ALTERNATING MINIMIZATION

Consider a multiphysics joint inverse problem—the same region of the Earth’s interior is imaged with N geophysical methods, each sensitive to a different physical property of the medium and modelled by a forward operator ${\mathbf{f}}_i$⁠; the goal is, given N observed data vectors ${\mathbf{d}}_i$⁠, to simultaneously recover the physical properties described by the model vectors ${\mathbf{m}}_i$⁠:

For each method, assuming the noise ${\mathit{ϵ}}_i$ is Gaussian with zero mean, the maximum likelihood solution can be found by minimizing the data misfit functional

where ${\mathbf{C}}_{d_i}$ is the covariance matrix of the observed data. Because this solution is generally non-unique due to ill-posedness of the inverse problem (1) and unattainable for local optimization due to nonlinearity of ${\mathbf{f}}_i({\mathbf{m}}_i)$⁠, prior information about the model must be introduced—in the form of regularization and, in the case of the multiphysics joint inversion, also by postulating some coupling between the individual models. This coupling, as well as the regularization, can generally be imposed as a penalty functional. Denoting by $\mathcal{R}$ the combination of all the ‘prior’ functionals, we can set the joint inversion as an unconstrained optimization problem:

Coefficients ${\lambda }_i>0$ weight the individual data misfit functionals relative to each other and to the ‘prior’ terms. Practical implementation of eq. (3) can be difficult due to the complexity and heterogeneity of the objective function—the individual inverse problems, as well as the ‘prior’ terms, often have different requirements to model parametrization and optimization techniques (and regardless of this, maximum reuse of the existing codes is always desirable). A straightforward workaround is based on the block coordinate descent (BCD) approach (Hu et al. 2009; Haber & Holtzman Gazit 2013) or alternating minimization: the idea is to cycle through the models ${\mathbf{m}}_i$ (i.e. blocks of $[{\mathbf{m}}_0^{\mathrm{T}},{\mathbf{m}}_1^{\mathrm{T}},\dots ,{\mathbf{m}}_{N-1}^{\mathrm{T}}{]}^{\mathrm{T}}$⁠) each time minimizing eq. (3) with respect to one model while keeping the other models fixed; this decouples the ${\mathrm{\Phi }}_{d_i}$ terms, significantly simplifying practical implementation of the joint inversion (at the cost of theoretically slower convergence). We develop a different alternating minimization approach, whose advantages over BCD will be demonstrated below.

Let us introduce for each ${\mathbf{m}}_i\in {\mathbb{R}}^{M_i}$ an auxiliary variable ${\mathbf{u}}_i\in {\mathbb{R}}^{M_i}$ and, instead of eq. (3), write:

With all ${\alpha }_i\to \mathrm{\infty }$⁠, this problem is equivalent to eq. (3). This technique, known as variable splitting, following the works of Wang et al. (2008) and Goldstein & Osher (2009) is widely used with ${\ell }_1$-regularized problems as a means of separating the non-differentiable ${\ell }_1$ and the ${\ell }_2$ terms of the objective function. Here it separates the ${\mathrm{\Phi }}_{d_i}$ terms from the regularization and coupling term $\mathcal{R}$⁠.

Assume that each vector ${\mathbf{m}}_i=[m_{i,1},m_{i,2},\dots ,m_{i,M_i}{]}^{\mathrm{T}}$ parametrizes a grid-based discretization of a physical property $m_i(\mathbf{x})$ of the medium in the following form:

where ${\phi }_{i,j}(\mathbf{x})$ are some basis functions defining a 3-D grid. Thus, we assume that individual models are discretized independently and the model regions ${\mathrm{\Omega }}_i$ are, in general, different (with a necessary limitation that the intersection of all ${\mathrm{\Omega }}_i$ is not empty). In contrast to ${\mathbf{m}}_i$⁠, vectors ${\mathbf{u}}_i$ in eq. (4) are coupled by the term $\mathcal{R}$ which implies that they must, at some point, be mapped to the same basis. Let us discretize all ${\mathbf{u}}_i$ by the same grid, that is analogously to eq. (5), $u_i(\mathbf{x})=\sum _{j=1}^M{u}_{i,j}{\phi }_j(\mathbf{x})$⁠, with $\mathbf{x}\in \mathrm{\Omega }:\bigcap _{i=0}^{N-1}{\mathrm{\Omega }}_i\subset \mathrm{\Omega }$⁠, and denote by matrix ${\mathbf{P}}_i$ interpolation onto this grid from the grid of ${\mathbf{m}}_i$⁠. Let us also define the multiparameter model vector $\mathbf{u}=[{\mathbf{u}}_0^{\mathrm{T}},{\mathbf{u}}_1^{\mathrm{T}},\dots ,{\mathbf{u}}_{N-1}^{\mathrm{T}}{]}^{\mathrm{T}}$⁠. Using this notation, we can rewrite problem (4) as follows:

The problem is thus decomposed into a system of N data-fitting (inversion) subproblems, each regularized by a simple quadratic functional, and the coupling-regularization subproblem. Derivatives of the linking terms with respect to ${\mathbf{m}}_i$ required by the inversion subproblems include transpose ${\mathbf{P}}_i^{\mathrm{T}}$ of the interpolation operators, which is often difficult to implement; an iterative solver such as the conjugate gradient method would require computationally expensive (for large 3-D grids) application of ${\mathbf{P}}_i$ and ${\mathbf{P}}_i^{\mathrm{T}}$ at each iteration. However, given that both ${\mathbf{m}}_i$ and ${\mathbf{u}}_i$ adequately represent the models (e.g. satisfy the Nyquist criterion), distance between the models can as well be calculated in the ${\mathbf{m}}_i$ space. Another difficulty with the inversion subproblems, especially when the gradient of ${\mathrm{\Phi }}_{d_i}({\mathbf{m}}_i)$ is ‘rough’ (as, e.g. in ray-based tomography), is that, even if the smoothness of ${\mathbf{u}}_i$ is controlled by $\mathcal{R}(\mathbf{u})$⁠, lack of smoothness in the updates to ${\mathbf{m}}_i$ at early stages, when ${\mathbf{m}}_i$ and ${\mathbf{u}}_i$ are far apart, may lead to instability and breakdown of the inversion. The required additional smoothing can be introduced by replacing the Euclidean norm in the linking terms by a Sobolev norm. Applying these changes to system (6), we get:

where ${\hat{\alpha }}_i={\alpha }_i/{\lambda }_i$⁠, ${\mathbf{Q}}_i$ is interpolation from the grid of $\mathbf{u}$ onto the grid of ${\mathbf{m}}_i$ and $W_2^1$ denotes (weighted) discretized Sobolev norm of order one:

where $\mathbf{L}$ is a negative discrete Laplacian and w is a weighting coefficient. In the coupling-regularization subproblem, we omit additional smoothing since in the cases considered here it is either redundant or conflicting with the regularization present in $\mathcal{R}$⁠. Note that if all the models and $\mathbf{u}$ use the same discretization and ${\alpha }_i\to \mathrm{\infty }$⁠, problem (7) is equivalent to eq. (3).

Discretizations of ${\mathbf{m}}_i$ are dictated by the inverse problems. Since the regularization and coupling should not create excessive details, discretizing $\mathbf{u}$ by a grid that is not coarser than the finest ${\mathbf{m}}_i$ grid is always safe. However, this is not always necessary and computationally feasible. For example, inversion of MT data often requires a very fine grid in the near surface, while its joint inversion with tomography and/or gravity aims at resolving the deeper (and larger scale) structures. In this case, discretizing $\mathbf{u}$ by a regular grid with the smallest step from the resistivity grid would be prohibitively expensive, but since $\mathbf{u}$ enters MT inversion as a reference model only, using a coarser grid will still allow for recovery of a fine resistivity structure in the near surface (see Section 5.2).

System (7) can be solved by alternating minimization: several (or one, for Newton-type methods) iterations of the inversion subproblems are followed by several iterations of the coupling-regularization subproblem, after which this cycle, which we will refer to as the outer iteration, is repeated. We use a continuation (see, e.g. Nocedal & Wright 1999) in ${\alpha }_i$ and ${\hat{\alpha }}_i$⁠: following Wang et al. (2008), ${\alpha }_i$ are gradually increased with the outer iterations as a geometric series ${\alpha }_{i,n}=q^n{\alpha }_{i,0}$⁠, where $q>1$⁠. On the other hand, ${\hat{\alpha }}_i={\alpha }_i/{\lambda }_i$⁠. Recalling eq. (3), we recognize ${\lambda }_i$ as the Lagrange multipliers of the constraints $2{\mathrm{\Phi }}_{d_i}({\mathbf{m}}_i)=N_{d_i}$⁠, where $N_{d_i}$ is the dimension of ${\mathbf{d}}_i$⁠; a common practical approach to determining these Lagrange multipliers, known as cooling or adaptive regularization, is to use continuation ${\lambda }_{i,n}=q_{\lambda }^n{\lambda }_{i,0}$ with $q_{\lambda }>1$⁠, starting from a “small” value ${\lambda }_{i,0}$⁠, until the constraint is satisfied (see, e.g. Haber et al. 2000; Zhdanov 2002). Assuming that $q_{\lambda }\approx q$⁠, we can keep ${\hat{\alpha }}_i$ constant during the inversion–this simplification was used in all our numerical examples and proved to be working satisfactorily. Overall, with ${\hat{\alpha }}_i$ fixed and ${\alpha }_i$ increasing, the balance is gradually shifted from the prior $\mathcal{R}$ towards the data terms, similar to the (separate inversion) cooling. A finer tuning of ${\alpha }_i$ and ${\hat{\alpha }}_i$ is usually needed at the final iterations in order to prevent overfitting and to assure simultaneous convergence of all the methods (see lines 11–15 of Algorithm 1 below). Convergence is controlled through the root-mean-square (RMS) values of the covariance-weighted data misfits:

and the relative mean misfits between the models and the auxiliary variables:

Parameters ${\alpha }_{i,0}$ and ${\hat{\alpha }}_i$ are problem-dependent and are selected in such a way that the effects of coupling and regularization are maximized while acceptable values of ${\text{RMS}}_i$ and $r_i$ are reached. Moreover, ${\alpha }_{i,0}$ control relative importance of the individual methods in the joint inversion. Practical aspects of this are discussed in Section 5.1.

The resulting inversion scheme is summarized below as Algorithm 1.

Decoupling minimization problems for the different data misfit functionals simplifies programming and allows each inverse subproblem to be solved using both numerical grid and optimization method tailored to the specific physics. Compared to the BCD, Algorithm 1 does not require implementation of $\mathcal{R}$ on the side of each inversion subproblem. Another advantage over the BCD is that the steps in different ${\mathbf{m}}_i$ can be taken simultaneously, that is the inverse problems can be run in parallel. The linking terms act as the standard quadratic ‘damping + smoothing’ stabilizers and the information exchange between the subproblems is performed by means of the reference models, thus most existing (separate) inversion codes can be used without any modification. With $N=1$⁠, Algorithm 1 effectively changes regularization in an existing inversion code. Similarly, isolation of the coupling-regularization subproblem simplifies its implementation and modification and allows the use of tailored optimization methods. Strategies of Um et al. (2014) and Gao & Zhang (2018) demonstrate a similar level of decoupling, but instead of variable splitting they essentially rely on sequential minimization of components of the joint objective function: decoupled inversions are restarted from the models that are ‘preconditioned’ by minimization of cross-gradient with damping. Another similar strategy, although restricted to the case of explicit petrophysical relationships between the models, is that of Heincke et al. (2017) and Kamm et al. (2015).

3 REGULARIZATION AND COUPLING

Here we give several illustrative examples of $\mathcal{R}(\mathbf{u})$⁠, focusing on the more universal structural priors.

For discretization of the regularization-coupling subproblem, we use a ‘voxel’ regular grid with step h and $M=M_x\times M_y\times M_z$ cells. With (⁠${\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_M$⁠) being lexicographically ordered coordinates ${\mathbf{x}}_{i,j,k}=(x_i,y_j,z_k)$ of the cell centres, the multiparameter model vector is defined as follows:

Partial derivatives are approximated by forward differences and zero Neumann boundary conditions are applied:

which we write in matrix form as ${\partial }_x{u}_i\approx {\mathbf{D}}_x{\mathbf{u}}_i$⁠, and analogously for y and z directions.

To update $\mathbf{u}$ at each inner iteration of Algorithm 1, the following system of linear equations is solved:

where ${\mathbf{H}}_{\mathcal{R}}$ is some approximation of the Hessian of $\mathcal{R}$⁠, $\mathbf{A}=\mathrm{diag}[{\alpha }_i]\otimes \mathbf{I}$⁠, $\mathbf{I}$ is $M\times M$ identity matrix and $\otimes$ denotes Kronecker product. Then, $\mathbf{u}\leftarrow \mathbf{u}+\mu \delta \mathbf{u}$ with the step length $\mu \le 1$ determined by backtracking line search. Here we focus on Gauss–Newton-like forms of eq. (13), as from our experience they combine ease of implementation with reasonable efficiency; the provided expressions for the gradient can be used to implement gradient-based and quasi-Newton schemes.

3.1 Joint total variation

Total variation (TV) regularization is efficient in many inverse problems because of its edge-preserving properties (Rudin et al. 1992). The TV is defined as ${\ell }_1$ norm of the magnitude of the model’s gradient; in order to make the absolute value differentiable at zero it is commonly ‘smoothed’ using a small positive constant $\beta$⁠:

TV of the multiparameter model $\mathbf{u}$ can be defined simply as a sum of TV functionals of the individual models: ${\mathcal{R}}_{\mathcal{\text{TV}}}(\mathbf{u})=\sum _i{\mathcal{R}}_{\mathcal{\text{TV}}}({\mathbf{u}}_i)$⁠.

Another way to generalize TV to the multiparameter model case is by using the ${\ell }_{1,2}$ mixed norm: taking ${\ell }_2$ norm of the vector $[|\mathrm{\nabla }{u}_0|,|\mathrm{\nabla }{u}_1|,\dots ,|\mathrm{\nabla }{u}_{N-1}|{]}^{\mathrm{T}}$ at each spatial position and then ${\ell }_1$ norm of the result gives functional known as colour, vectorial or joint total variation (JTV, Blomgren & Chan 1998; Haber & Holtzman Gazit 2013; Crestel et al. 2018). We use a slightly generalized, weighted version of JTV:

where ${\mathbf{W}}_i=\mathrm{diag}[w_{i,x},w_{i,y},w_{i,z}]$⁠, $w_{i,\xi }>0$⁠. This functional is convex (see the references above) and thus a stabilizing functional according to Tikhonov & Arsenin (1977), at the same time it introduces coupling between the models by aligning their edges, that is local maxima of the magnitude of the gradient. The coupling property of JTV can be illustrated by the following example: consider $\mathbf{u}$ consisting of two piecewise constant models with the jumps of unit amplitude and with unit values of the total variation; it is clear that if the jumps in the two models are disjoint, then ${\mathcal{R}}_{\mathcal{\text{JTV}}}(\mathbf{u})=2$⁠, but if instead they are fully coincident, ${\mathcal{R}}_{\mathcal{\text{JTV}}}(\mathbf{u})=\sqrt{2}$⁠.

A simple and relatively efficient way to minimize TV and JTV functionals is by using the iteratively reweighted least squares (IRLS) method (Vogel & Oman 1996; Molodtsov & Troyan 2017). It leads to an approximation of the Hessian of JTV that is block-diagonal, that is system (13) is decoupled into the following N systems:

where

and $g_{w;i,j}=\sqrt{w_{i,x}[{\mathbf{D}}_x{\mathbf{u}}_i{]}_j^2+w_{i,y}[{\mathbf{D}}_y{\mathbf{u}}_i{]}_j^2+w_{i,z}[{\mathbf{D}}_z{\mathbf{u}}_i{]}_j^2}$⁠.

Matrix ${\mathbf{L}}_{w;i}$ can be viewed as a finite-difference diffusion operator with anisotropic (due to the weighting ${\mathbf{W}}_i$⁠) diffusivity that depends on the $\mathbf{u}$ computed at the previous IRLS iteration. The edges (large $g_{w;i,j}$⁠) locally reduce all the diffusivities and in this way persist between the iterations, resulting in edge-preserving smoothing (similar to TV), and between the models, resulting in coupling. It follows that coincident edges are more resistant to smoothing, but new edges can only be introduced by ${\mathbf{u}}_{\text{ref}}$⁠, that is by the data.

Because systems (16) are symmetric positive definite, they can be efficiently solved by the preconditioned conjugate gradient (PCG) method. We use Jacobi preconditioner and, as usual in Gauss–Newton-like schemes (see, e.g. Nocedal & Wright 1999), terminate PCG early by setting relative tolerance to 0.001–0.01. The condition number of the system depends on $\beta$⁠; if $\beta$ is too small (and when ${\alpha }_i$ are small) PCG may stagnate. We find that, given the model gradients are reasonably scaled (see below), $\beta ={10}^{-7}$ provides a good balance between PCG convergence rate and the accuracy of ${\ell }_1$ norm approximation.

Scaling of the models, which represent different physical quantities, is important for efficient coupling; apart from the bounding transforms that are discussed below, it is achieved by defining ${\mathbf{W}}_i=s_i^{-1}\mathrm{diag}[w_x,w_y,w_z]$⁠, where $s_i$ is the RMS magnitude of the gradient of ${\mathbf{u}}_{\text{ref},i}$⁠. Similarly to the classical ${\ell }_2$-norm regularization case (e.g. Li & Oldenburg 1998), constants $w_x=w_y\approx 10$ and $w_z=1$ usually work well for ‘subhorizontal’ structures, while some scenarios with steeply dipping geology may benefit from $w_z>w_x=w_y$⁠.

3.2 Joint minimum support

For linear inverse problems, regularization based on the ${\ell }_1$ norm leads to sparse solutions (Loris et al. 2007). Similarly, the ${\ell }_{1,2}$ mixed norm promotes joint sparsity—same sparsity pattern between the individual model vectors (e.g. Gramfort et al. 2012). These approaches are particularly efficient because the corresponding stabilizers are convex. Non-convex stabilizers are more problematic from the optimization standpoint, but have stronger sparsity-promoting properties; a well-known example is the minimum support (MS) functional introduced by Last & Kubik (1983) and generalized by Portniaguine & Zhdanov (1999):

With $\beta \to 0$⁠, it approximates the volume of the regions of $\mathrm{\Omega }$ where $u_i$ deviates from a (known) background model $u_{b;i}$⁠, that is the volume of the support of $u_i(\mathbf{x})-u_{b;i}(\mathbf{x})$⁠. In the multiparameter case, by analogy with mixed norms, we arrive at the joint minimum support (JMS) functional (Molodtsov & Troyan 2017):

Similarly to the ${\ell }_{1,2}$ norm, JMS is a stabilizing functional and introduces coupling through joint sparsity. The coupling property follows from the fact that the integrand of eq. (19), the function $t\mapsto t/(t+\beta )$ with $t\ge 0$⁠, is subadditive. JMS regularization tends to ‘focus’ the models $u_i-u_{b;i}$ such that they have the same support, which can be especially useful in time-lapse problems. It was applied to magnetic, gravity and marine controlled-source electromagnetic data (Jorgensen & Zhdanov 2021; Tu & Zhdanov 2022).

To illustrate the behaviour of JMS, let us consider an ‘inverse problem’ with one data point $d=f({\mathbf{m}}_0)$ and two model parameters ${\mathbf{m}}_0=(m_{0,1},m_{0,2})$ representing two voxels and let ${\mathbf{m}}_1=(m_{1,1},m_{1,2})$ be a ‘fixed’ reference model. JMS of ${\mathbf{m}}_0$ and ${\mathbf{m}}_1$ (with zero background models and $\beta =0.001$⁠) is plotted in Fig. 1: at ${\mathbf{m}}_1=\mathbf{0}$⁠, JMS regularization selects two sparse solutions to $d=f({\mathbf{m}}_0)$⁠, but as we move away from zero along $m_{1,1}$⁠, only the solution having the same sparsity pattern as ${\mathbf{m}}_1$ remains.

Minimization of JMS can be carried out using IRLS (by analogy with MS, following Last & Kubik (1983) and Portniaguine & Zhdanov (1999)). In the vicinity of a model ${\mathbf{u}}^{(k)}$ found at IRLS iteration k, (eq. 19) is approximated by a quadratic functional:

with

Substituting it into the eq. (13) and taking unit step length, we get an explicit formula for the updated model vector components:

for $i=0,\dots ,N-1$⁠, $j=1,\dots ,M$⁠. JMS is straightforward to modify to include scaling of the models similarly to JTV (but with the coefficients $s_i$ equal to RMS values of ${\mathbf{u}}_{\text{ref},i}-{\mathbf{u}}_{b;i}$⁠).

3.3 Modifications of cross-gradient

The cross-gradient constraint of Gallardo & Meju (2003) can be imposed, for a pair of models, by minimizing ${\ell }_2$ norm of the cross product of the model gradients:

The cross product is zero if the gradients (and, consequently, the level sets of the models) are parallel or if a gradient vanishes. The latter means that the constraint is rather inefficient for piecewise-constant models and also demonstrates that (eq. 23) is not a stabilizing functional. JTV, on the other hand, does not efficiently couple smooth models since diffusivity in eq. (17) depends only on the magnitude of the gradients. Naturally occurring models are usually combinations of the smooth and the piecewise constant cases, so it seems reasonable to use JTV as a stabilizing functional for the cross-gradient. Considering all pairs of the models and introducing for each pair a positive weighting coefficient ${\gamma }_{ij}$⁠, we have

The gradient of eq. (24) consists of the blocks

where ${\mathbf{r}}_{ij,\xi }$ is the discretized $\xi$ component of $\mathrm{\nabla }{u}_j\times (\mathrm{\nabla }{u}_i\times \mathrm{\nabla }{u}_j)$⁠. Using Gauss–Newton approximation for the cross-gradient terms, the Hessian of eq. (24) is approximated by

where

and

The resulting linear system (13) is solved by the PCG analogously to the JTV case, the main difference is that the search directions $\delta {\mathbf{u}}_i$ for the individual models are now coupled.

Denoting by $\theta$ the angle between the two gradients, we can factorize the cross-gradient as follows:

The first two factors zero the cross-gradient when the corresponding gradient vanishes, the third factor—when the gradients are antiparallel and the fourth—parallel. In many cases, the sign ${\sigma }_{ij}$ of correlation between the two reconstructed properties is known a priori and this information should be incorporated into the coupling constraints. This can be done (following Molodtsov et al. 2011) by selecting either the ‘parallel’ or the ‘antiparallel’ factor of eq. (29) and taking ${\ell }_2$ norm of the result, which leads to the one-way cross-gradient functional:

where ${\sigma }_{ij}=\pm 1$⁠. Minimization of this functional with ${\sigma }_{ij}=1$ forces model gradients to be parallel, with ${\sigma }_{ij}=-1$ – antiparallel and in both cases the gradients are free to vanish.

The gradient and Gauss–Newton approximation to the Hessian of eq. (30) are given by the following expressions:

where ${\mathbf{r}}_{ij}$ is the discretized $|\mathrm{\nabla }{u}_i||\mathrm{\nabla }{u}_j|-{\sigma }_{ij}\mathrm{\nabla }{u}_i\cdot \mathrm{\nabla }{u}_j$⁠,

and $g_{i,k}=\sqrt{\sum _{\xi }[{\mathbf{D}}_{\xi }{\mathbf{u}}_i{]}_k^2+\beta }$⁠. These expressions are substituted into eqs (25) and (26) when some of the cross-gradient terms in the functional (24) are replaced by eq. (30).

Regularization and coupling can be easily generalized by appending fixed models to the vector $\mathbf{u}$⁠, for example in order to use geological or legacy geophysical models as a structural reference. Extra regularization (such as ‘damping’ towards a reference model) and coupling (such as ad hoc petrophysical relations, as demonstrated in the numerical example in Section 5.2) terms are often required when dealing with real data and are straightforward to implement.

4 APPLICATION TO MT, FATT AND GRAVITY

Here we outline the implementation of the inversion subproblems used in the numerical examples.

3-D MT inversion is performed using the nonlinear conjugate gradient (NLCG) method implemented in the ModEM software package (Kelbert et al. 2014). Resistivity is discretized on a non-uniform rectangular grid, with the cells increasing in size with the decay of the MT resolution with depth and laterally outside of the data-covered area; the model vector $\mathbf{m}$ is formed from the natural logarithm of the cells’ resistivities. Regularization in ModEM has the form $(\mathbf{m}-{\mathbf{m}}_{\text{ref}}{)}^{\mathrm{T}}{\mathbf{C}}_m^{-1}(\mathbf{m}-{\mathbf{m}}_{\text{ref}})$⁠, with the model covariance ${\mathbf{C}}_m$ implemented as a first-order autoregressive smoothing operator; it can be shown (Yanovskaya & Ditmar 1990; Xu 2005) that it is equivalent to a squared $W_2^1$ Sobolev norm. NLCG iterations are performed in transformed model space $\tilde{\mathbf{m}}={\mathbf{C}}_m^{-1/2}(\mathbf{m}-{\mathbf{m}}_{\text{ref}})$⁠. In $\mathbf{u}$⁠, a logarithmic resistivity vector is used, so $\mathbf{m}$ is transformed to $\tilde{\mathbf{m}}$ and back at each outer iteration of Algorithm 1. To make use of conjugacy, several NLCG iterations are needed at each outer iteration; at the same time, for the alternating minimization with continuation, this number should not be too large. In our experiments, five iterations were sufficient for a satisfactory convergence rate.

Seismic first-arrival traveltime tomography (FATT) is implemented using the Gauss–Newton method, a single iteration of which is performed at each outer iteration of Algorithm 1. Velocity $v_p$ is defined over a voxel grid. The forward problem is solved by the 3-D finite-difference method of Podvin & Lecomte (1991) and its Jacobian is computed using our implementation of the ‘backtracking’ ray tracing, principle of which is described in the same paper; normal equations are solved by PCG; computations are parallelized over the sources. Inversion (as well as regularization and coupling) is performed in the transformed model domain, with the transform

which guarantees that the slowness $s=1/v_p$ lies within the predefined interval $(s_{min},s_{max})$⁠. Details on the bounding transform eq. (34) can be found in Commer & Newman (2008).

Our implementation of 3-D gravity inversion is based on the works of Li & Oldenburg (1998, 2003). A single Gauss–Newton iteration with PCG as a linear solver is performed at each outer iteration of Algorithm 1. The model is discretized with a non-uniform rectangular grid and analytical expressions for the gravity field of a constant-density rectangular prism are used to precompute the sensitivity matrix. Each row of the sensitivity matrix is compressed using a 3-D discrete wavelet transform (DWT) with the D4 Daubechies wavelet. Trade-off between the compression ratio and the forward problem accuracy is controlled by setting the relative error of the sensitivity reconstruction. The depth (or distance) weighting (Li & Oldenburg 1998) is applied prior to the compression; in addition to improving the recovery of deep structures, it also reduces the dynamic range of the sensitivity thus reducing the compression artefacts. Inversion is performed in the depth/distance-weighted model space, but the vector ${\mathbf{u}}_i$ is unweighted, that is the weighting is only applied to the linking terms in the inversion subproblem. To constrain the density range, instead of a logarithmic barrier method used by Li & Oldenburg (2003), we use the bounding transform eq. (34).

Trilinear interpolation is used for ${\mathbf{P}}_i$⁠, ${\mathbf{Q}}_i$⁠, which is quite efficient when all the grids are rectangular. In case of a mismatch between the regions $\mathrm{\Omega }$ and ${\mathrm{\Omega }}_i$⁠, nearest-neighbour extrapolation is invoked and the interpolation-extrapolation result m is smoothly merged with the prior (or initial) model $m_0$ outside the interpolation region via a sigmoid function s centred on the regions boundary as $m\leftarrow sm+(1-s)m_0$⁠.

5 NUMERICAL EXAMPLES

5.1 Example 1: time-lapse seismic–gravity joint inversion using JMS regularization

Joint inversion attempting to reconstruct compact, highly (co)localized features can be efficiently performed with the help of the joint minimum support regularization. To demonstrate this, we consider a simple synthetic example combining time-lapse FATT and full-tensor gravity gradiometry (FTG). Seismic tomography is widely used for time-lapse applications at various scales (e.g. Patanè et al. 2006; Ajo-Franklin et al. 2007; Hilbich 2010). The feasibility of time-lapse applications of FTG was studied, for example by Droujinine et al. (2007) and by Reitz et al. (2015); in this example, we use six components of the gravity gradient tensor (⁠$U_{xx}$⁠, $U_{yy}$⁠, $U_{zz}$⁠, $U_{xy}$⁠, $U_{xz}$⁠, $U_{yz}$⁠) and the vertical gravity component $U_z$⁠.

Using notation introduced in Section 2, let us consider two experiments–FATT (⁠$i=0$⁠) and FTG (⁠$i=1$⁠), both performed twice—at time $t_0$ and $t_1$⁠:

In the context of time-lapse inversion, the states at $t_0$ and $t_1$ are usually referred to as baseline and monitor, respectively. If the data acquisition geometry stays constant, it is advantageous to work with the difference of monitor and baseline data, $\mathrm{\Delta }{\mathbf{d}}_i={\mathbf{d}}_i(t_1)-{\mathbf{d}}_i(t_0)$⁠. A common strategy is: (1) find ${\mathbf{m}}_i(t_0)$ by inverting baseline data and (2) find the model temporal perturbation $\mathrm{\Delta }{\mathbf{m}}_i$ by minimizing the double-difference data misfit functional (Waldhauser & Ellsworth 2000):

Given smallness of $\mathrm{\Delta }{\mathbf{m}}_i$⁠, the second step is relatively fast and both observational and modelling systematic errors are cancelled in eq. (36), increasing sensitivity to the temporal changes, compared with a straightforward sequential inversion of baseline and monitor data. The gravity forward problem is linear with respect to the density, so for the FTG inversion, eq. (36) simplifies to

that is there is no need for the baseline density model if we are only interested in the temporal changes of density. Thus, after the reconstruction of the baseline velocity, which we perform by stand-alone FATT, time-lapse joint inversion reduces to the problem (7), in which the unknown models ${\mathbf{m}}_i$ are replaced with $\mathrm{\Delta }{\mathbf{m}}_i$⁠. Nevertheless, baseline gravity data can still be used in a scenario where baseline velocity model can be improved by joint inversion with gravity.

The SEG/EAGE Overthrust Model (Aminzadeh et al. 1997) is used as the baseline velocity model. The time-lapse velocity perturbation is a 1 km s⁻¹ cube with a 600 m edge, positioned horizontally at the centre of the model, at 1825 m depth (Fig. 2a). The true models are discretized by a regular cubic grid with 50 m spacing. 50 sources and 900 receivers are placed over a regular grid at 25 m depth (Fig. 2a); every source is registered by all the receivers. In the FTG problem, we directly model the time-lapse data difference. The density perturbation is 0.1 g cm⁻³ cube coincident with the velocity perturbation. Because there is no baseline density and discretizations of FTG and FATT are independent, a smaller grid can be used for density as shown in Fig. 2(b). The grid consists of a uniform core with 100 m horizontal and 90 m vertical steps padded at the lateral edges by three padding layers expanding by a factor of two, a total of ${64}^3$ cells. Vertical gravity and the gravity gradient components are calculated at ${13}^2$ points at 10 m height (Fig. 2b). Synthetic data are contaminated with Gaussian noise with the standard deviation of 1 ms for the traveltimes, 1 μGal for the gravity and 1 mE for the gravity gradient.

Gravity inversion is performed on the same grid that is used for generating the synthetic data. However, in the inversion, unlike in the synthetic data calculation, the sensitivity is compressed using DWT with the relative reconstruction error of 1 per cent, resulting in the compression ratio of 58; the depth weighting is set proportional to $1/z$⁠. For the tomography, the model region that was used for the forward modelling is rediscretized with a 100 m step. Bounding transform eq. (34) is applied to the models such that $v_p\in (0.2,6.5)$ km s⁻¹ and the density $\sigma \in (-0.2,0.2)$ g cm⁻³. Baseline traveltimes are inverted using Algorithm 1 with the TV regularization; initial model is 1-D, velocity linearly increases with depth from 2 to 6 km s⁻¹; reconstructed baseline model is shown in Fig. 2(c). The subsequent time-lapse inversion is essentially a simple ‘spike test’ (in which FTG problem is linear, but FATT is, of course, nonlinear and depends on the baseline).

Regularization for the time-lapse inversion is chosen based on the assumption that the volume of the support of the temporal changes is much smaller than the whole model volume—this corresponds to the MS functional (eq. 18) for separate inversions and the JMS functional (eq. 19) for joint inversion. The regularization-coupling grid is identical to the FATT inversion grid (⁠$200\times 200\times 46$ cells, 100 m spacing).

To assess the coupling effect of JMS regularization separately from its regularizing effect and to quantify the degree of alignment of the reconstructed anomalies, we introduce the following quantity:

Given that the anomalous part of $\mathbf{u}$ does not vanish everywhere, it follows from the definitions of MS and JMS functionals that $1/N\le \gamma (\mathbf{u})\le 1$ (recall that ${\mathcal{R}}_{\text{MS}}$ approximates the volume of the support of an individual model and ${\mathcal{R}}_{\text{JMS}}$—that of the multiparameter model, i.e. volume of the union of the individual supports). The upper bound, $\gamma (\mathbf{u})=1$⁠, corresponds to fully disjoint anomalies, whose supports do not intersect, and the lower bound, $1/N$⁠,—to full alignment, when the supports of all N anomalies are identical.

To determine the regularization parameters ${\alpha }_i$ (to simplify the notation, instead of writing ${\alpha }_{i,0}$ we will understand ${\alpha }_i$ as the initial values) and $\widehat{{\alpha }_i}$⁠, we use the following strategy:

1. Determine $\widehat{{\alpha }_0}$⁠, $\widehat{{\alpha }_1}$⁠. For each method, run a series of separate inversions with fixed ${\mathbf{m}}_{\text{ref},i}$ (i.e. Algorithm 1 without the coupling-regularization step) and varying $\widehat{{\alpha }_i}$ and select its optimal value according to the Occam/discrepancy principle—maximum $\widehat{{\alpha }_i}$ for which ${\text{RMS}}_i\approx 1$⁠.

2. Determine ${\alpha }_0$⁠, ${\alpha }_1$⁠. A reasonable starting point is ${\alpha }_i=\widehat{{\alpha }_i}$⁠. Run a series of joint inversions with varying ${\alpha }_i$ and select those that return minimum values of $\mathcal{R}$ (of course, while the inversion converges in ${\text{RMS}}_i$ and $r_i$⁠). Such combinations of ${\alpha }_i$ are, in general, non-unique. At the same time, we usually want to equalize contributions of the methods/data sets to the joint inversion as much as possible. Since the inverse problems are decoupled, relative weights of the data sets are primarily controlled by ${\alpha }_i$ and not $\widehat{{\alpha }_i}$⁠: increasing ${\alpha }_i$ increases the influence of ith data set on the joint inversion. It follows that good balance between the data sets can be achieved through similar final values of $r_i$⁠.

The values of $\widehat{{\alpha }_i}$ selected with separate inversions ensure convergence of the joint inversion in the data RMS; however, since in joint inversion, ${\mathbf{m}}_{\text{ref},i}$ are no longer fixed but move towards ${\mathbf{m}}_i$⁠, it is usually possible to increase $\widehat{{\alpha }_i}$ slightly (usually less than by an order of magnitude) to enhance the regularization and coupling. In the current example we simply use the values found at the step 1: $\widehat{{\alpha }_0}=100$ and $\widehat{{\alpha }_1}={10}^8$⁠.

With the $\widehat{{\alpha }_i}$ fixed, we perform grid search for ${\alpha }_0$ and ${\alpha }_1$⁠. For the purpose of finding the optimal parameters, full grid search is redundant, here it serves the illustrative purpose. By changing ${\alpha }_i$ by an order of magnitude, we build a 15 $\times$ 15 logarithmic grid, centred roughly on (⁠$\widehat{{\alpha }_0}$⁠, $\widehat{{\alpha }_1}$⁠). At each grid point $({\alpha }_0,{\alpha }_1)$⁠, we run the joint inversion with the maximum number of outer iterations set to 30, $q=2$ and the number of the coupling-regularization iterations $n_{\text{itn}}=100$⁠. To slightly optimize the search, if joint inversion does not reach the target value of $r_0$⁠, then the current ${\alpha }_0$ is deemed too low and exploration of the current column of the grid is terminated, and analogously for $r_1$ and the respective rows. Fig. 3 shows the grid with the final values of JMS, $\gamma (\mathbf{u})$⁠, ${\text{RMS}}_0$⁠, ${\text{RMS}}_1$⁠, $r_0-r_1$ and the iteration count. For each i, as expected, as ${\alpha }_i$ increases, $r_i$ and ${\text{RMS}}_i$ decrease and JMS increases. We are interested in the minimum of JMS and, to balance FTG and FATT, $r_0\approx r_1$⁠. The region satisfying these criteria, outlined by the rectangle in Fig. 3, is explored with a finer grid, with a factor of two step (Fig. 4). This time we use $q=1.5$⁠, which increases the iteration count and, since the maximum number of iterations is unchanged, shifts the region of convergence towards larger ${\alpha }_i$⁠. At the expense of more iterations, reducing q usually improves the inversion results. As we see from the fine grid, the quantities vary smoothly enough that determining ${\alpha }_i$ to an order of magnitude should give satisfactory results.

The minimum of ${\mathcal{R}}_{\mathcal{\text{JMS}}}$ is associated with FATT dominating FTG (Figs 4a and d), which is consistent with the overall higher resolution of FATT. For the point $({10}^5,{10}^8)$ in this region, marked with ‘$\mathrm{\Delta }$’ in Fig. 4, the reconstructed velocity and density models are shown in Figs 5 and 6(e) and (f). Best balance between FTG and FATT, in terms of the smallest $|r_0-r_1|$⁠, along the ${\mathcal{R}}_{\mathcal{\text{JMS}}}$ valley bottom is reached at the point $(6\times {10}^4,{10}^{10})$ marked with ‘+’ in Fig. 4; the corresponding models are shown in Figs 6(c) and (d).

For separate inversion with MS regularization, grid search results in ${\alpha }_0={10}^6$ and ${\alpha }_1={10}^{10}$⁠, and the corresponding reconstructed models are shown in Figs 6(a) and (b). Not surprisingly, FTG inversion correctly positions the density anomaly in the horizontal plane but not vertically and its upward displacement is compensated by its greatly reduced amplitude. Joint inversion with ${\alpha }_0=6\times {10}^4$ and ${\alpha }_1={10}^{10}$ (the best FTG–FATT balance) moves the density anomaly downwards (Figs 6c and d), but both velocity and density anomalies are poorly focused, as a result of the trade-off in JMS. Increasing the relative weight of FATT by moving down the valley of ${\mathcal{R}}_{\mathcal{\text{JMS}}}$ to ${\alpha }_0={10}^5$⁠, ${\alpha }_1={10}^8$ results in both anomalies satisfactorily aligned and positioned (Figs 6e and f); note the stronger focusing of the velocity anomaly compared to the separate FATT inversion.

Fig. 7 shows the progress of the inversions corresponding to the models in Fig. 6, in terms of ${\text{RMS}}_i$⁠, $r_i$⁠, $\gamma (\mathbf{u})$ and ${\mathcal{R}}_{\mathcal{\text{JMS}}}$⁠. The target RMS values are reached after the first few iterations, after which the models and the auxiliary variables continue to be gradually brought together by increasing ${\alpha }_i$ until their relative difference $r_i$ reaches the target value of 0.1. Note that FATT is ahead of FTG in terms of $r_i$⁠, except in Fig. 7(b) where they converge simultaneously. While the stabilizer, ${\mathcal{R}}_{\mathcal{\text{JMS}}}$⁠, is being minimized, the structure of $\mathbf{u}$ develops (by the transfer from ${\mathbf{m}}_i$⁠), resulting in the overall growth of ${\mathcal{R}}_{\mathcal{\text{JMS}}}$⁠. At the same time, in the joint inversions, $\gamma (\mathbf{u})$ decreases as the support of the anomalies becomes more similar. In the MS-regularized separate inversion (Fig. 7a), the values of ${\mathcal{R}}_{\mathcal{\text{JMS}}}$ are similar, but $\gamma (\mathbf{u})$ stays close to 1 as the anomalies are almost disjoint (Figs 6a and b). Note that the points of the $\gamma (\mathbf{u})$ and ${\mathcal{R}}_{\mathcal{\text{JMS}}}$ curves correspond to the inner IRLS iterations of the coupling-regularization subproblem. The sawtooth-like appearance of these curves is a result of the continuation—at each following outer iteration, IRLS starts with new ${\mathbf{u}}_{\text{ref}}$ and higher values of ${\alpha }_i$⁠, the increased relative weight of the linking terms causes a jump in ${\mathcal{R}}_{\mathcal{\text{JMS}}}$⁠. With $n_{\text{itn}}=100$⁠, we can clearly see the convergence of the IRLS. Terminating the IRLS very early with $n_{\text{itn}}=5$ (Fig. 7d) results in higher final value of ${\mathcal{R}}_{\mathcal{\text{JMS}}}$⁠, which is manifested as somewhat smoother models in Figs 6(a) and (b), but the coupling remains satisfactory, both visually and as indicated by $\gamma (\mathbf{u})$⁠.

5.2 Example 2: MT–seismic–gravity joint inversion using different regularization and coupling

This example is based on a modification of MT checkerboard-type benchmark model presented in Egbert & Kelbert (2012). Electrical resistivity, P-wave velocity and density models are displayed in Fig. 8; the models have the same pattern except for the top checkerboard layer (with the anomalous blocks being coincident, up to the grid discretization error). The ‘checkers’ are bounded by the planes $x,y=\{-27.6,-11.7,-8.7,7.2,10.2,26.1\}$ km and $z=\{0,1,1.8,5,5.8,10.6\}$ km.

The resistivity values are 10, 100 and 1000 $\mathrm{\Omega }\cdot$m. The resistivity model grid contains $66\times 66\times 59$ cells. It consists of the core area with the horizontal cell size of 1 km $\times$ 1 km and four lateral padding layers. In the vertical direction, the cell size increases with depth from 50 m to 16 km. Impedance tensor and vertical magnetic transfer functions data are computed over 15 $\times$ 15 regular grid of stations at the surface (Figs 8a and 9) for 12 logarithmically equispaced periods in 0.1–100 s range. Gaussian noise is added to the data, its standard deviation is 3 per cent of $\sqrt{|Z_{xy}||Z_{yx}|}$ for all components of the impedance tensor and 0.03 for the tipper.

The background velocity is linearly increasing from 4 km s⁻¹ at the surface to 6.8 km s⁻¹ at 20 km depth, the checkerboard pattern imposes $\pm 0.4$ km s⁻¹ anomalies with respect to this trend; the model is discretized on a regular cubic grid with 500 m step, resulting in $140\times 140\times 40$ cells. 64 sources and 225 receivers are placed at the surface over the regular grid with 8 and 4 km spacing, respectively (Figs 8b and 9); for each source, first arrival traveltimes are computed at all receivers within 60 km radius and contaminated with Gaussian noise with the standard deviation of 12.5 ms.

The density model contains two layers with $\pm 0.1$ g cm⁻³ anomalies and is discretized into ${64}^3$ cells; the cells are approximately 1 km $\times$ 1 km $\times$ 0.5 km in the core area, which is surrounded by two padding layers expanding with a factor of two (Figs 8c and 9). Vertical component of the gravity field is computed at 10 m above the surface at 225 points horizontally coinciding with the MT locations and Gaussian noise with a standard deviation of 0.1 mGal is added.

The synthetic data are inverted using Algorithm 1 with different stabilizing and coupling functionals: (a) TV, corresponding to separate inversions, (b) JTV, (c) JTV combined with the cross-gradient (XG), (d) TV combined with XG and the one-way cross-gradient (OWXG), (e) JTV combined with XG and OWXG and (f) JTV combined with XG, OWXG and a ‘petrophysics’ coupling; the corresponding reconstructed models are shown in Fig. 11. Apart from the choice of the stabilizing and coupling functionals, all the inversions are run with the same settings, which are summarized below.

The MT inversion is started from the homogeneous 100 $\mathrm{\Omega }\cdot$m model. The number of NLCG iterations is set to five; the initial line-search step length is set to 20. The model covariance smoothing is applied twice with the factor 0.3 along each coordinate axis. FATT is started from the 1-D model corresponding to the true background velocity trend, with the model transform bounds 0.96 and 7.2 km s⁻¹; the relative tolerance of PCG is set to 0.01. Gravity inversion is started from a zero density model, with the model transform bounds −2.5 and 2.5 g cm⁻³. Depth weighting $\sim 1/z$ is used and the DWT reconstruction error is set to 5 per cent, resulting in the compression ratio of 77; PCG tolerance is set to ${10}^{-3}$⁠. All the inversion grids are the same as those of the true models used for computing the synthetic data.

In order to determine the regularization parameters ${\alpha }_i$ and $\widehat{{\alpha }_i}$⁠, we follow the same strategy as in Example 1: first find $\widehat{{\alpha }_i}$ using separate inversions, then find ${\alpha }_i$ that minimize JTV in the JTV-regularized joint inversion; however, because of the higher dimensionality and numerical complexity of the problem, we determine ${\alpha }_i$ by trial and error (using order-of-magnitude steps) without exhaustive grid search, understanding that the results may be suboptimal. The same values of ${\alpha }_i$ and $\widehat{{\alpha }_i}$ are then used in the rest of the inversions (with XG, OWXG and ‘petrophysics’ coupling), as we attempt to compare the effects of the coupling constraints. For all the cases, the results appear to be satisfactory, however, they probably can be improved by changing ${\alpha }_i$⁠. To determine ${\gamma }_{ij}$⁠, with ${\alpha }_i$ and $\widehat{{\alpha }_i}$ fixed, we choose the maximum values for which the convergence is maintained and the models are sufficiently ‘smooth’ (in terms of (J)TV).

Using the above procedure, the parameters are set as follows (with index 0 denoting FATT, 1–gravity, 2–MT): ${\alpha }_0=200$⁠, ${\hat{\alpha }}_0=20$⁠, ${\alpha }_1={10}^5$⁠, ${\hat{\alpha }}_1={10}^6$⁠, ${\alpha }_2=1$⁠, ${\hat{\alpha }}_2=10$ and $q=1.3$⁠. In the inversions with XG, ${\gamma }_{01}={\gamma }_{02}=200$⁠, ${\gamma }_{12}=20$ and when OWXG is applied to the FATT–gravity pair, ${\gamma }_{01}=2000$ and ${\sigma }_{01}=-1$⁠.

For FATT and gravity, the Sobolev norm weight $w=1000$⁠. For TV and JTV, $w_x=w_y=10$ and $w_z=1$⁠. The number of Gauss–Newton iterations $n_{\text{itn}}=4$ and PCG tolerance is set to 0.01. The regularization-coupling grid is the same as the FATT grid.

To measure the alignment of the models’ edges, similarly to the JMS case considered in Example 1, we normalize JTV and define

which is bounded by $1/\sqrt{N}$ and 1.

The change of RMS data misfits, relative distances $r_i$ between ${\mathbf{m}}_i$ and ${\mathbf{u}}_i$⁠, XG and OWXG functionals and ${\gamma }_{\text{TV}}(\mathbf{u})$ with the outer iterations is plotted in Fig. 10. Note that, except for the Fig. 10(c), the same functionals are plotted irrespective of whether or not they are part of the corresponding objective function. Table 1 summarizes the relative performance of the inversions in terms of the relative mean error with respect to the true models: $\text{mRMS}=\Vert {\mathbf{m}}_i-{\mathbf{m}}_{i,\text{true}}\Vert /\Vert {\mathbf{m}}_{i,\text{true}}\Vert \cdot 100\mathrm{\%}$⁠, as well as the RMS covariance-weighted data misfits and the number of iterations.

Fig. 11 presents a comparison of the models obtained with different combinations of the structural constraints. Using JTV instead of TV mainly increases the contrast of the anomalies, while the improvement in their structural similarity is minor, most visible in vertical sections of the density model. Adding XG further increases the structural similarity, as evidenced, for example by the shape of the low-velocity blocks (see 3 km depth slices) and the reduced ‘bleeding’ in the resistivity models. Postulating positivity of the velocity–density correlation via the use of OWXG appears to be critical for the reconstruction of density (cf. columns 3 and 4–6 of Fig. 11)—essentially, it brings depth resolution into the gravity inversion. This effect is similar to the effect of the positivity constraint on the susceptibility in magnetic data inversion (Li & Oldenburg 1998). However, also applying OWXG to the resistivity only marginally improves the models and does not result in reconstruction of the bottom-layer resistive blocks since their absence satisfies the OWXG constraint (see the irregular-geometry test below). Note that the pattern of the upper checkerboard layer is different for velocity, resistivity and density, that is structure present in one of the properties is absent in the others, and all the tested structural constraints work fine in this situation, producing no artefacts.

In the next test, we examine the effects of direct ‘petrophysics’ coupling, assuming such information is available for a certain subset of the model parameters. With the exception of the top checkerboard layer, true velocity $v_p$ and logarithm of resistivity $\rho$ are linearly related as follows:

where $u_0$ is the slowness model transformed by eq. (34), $u_2=\ln \rho$ and $v_{p,0}$ is the background velocity. This constraint is enforced for the bottom checkerboard layer only by adding to $\mathcal{R}(\mathbf{u})$ the term

where $\mathbf{W}$ is diagonal matrix whose entries are 1000 for the cells below $z=5$ km and zero otherwise. Adding it to the combination of JTV, XG and OWXG from the previous test somewhat accelerates convergence of the inversion and dramatically improves models of all the properties (see Table 1), in particular, resistive blocks of the bottom checkerboard layer are now recovered (last column of Fig. 11).

To simulate an irregular acquisition geometry often encountered in the real data, particularly in regional studies, we randomly draw from the regular-grid data sets 25 per cent of the MT and gravity locations and 50 per cent of seismic sources and receivers (Fig. 12). These subsamples of the data are inverted separately using TV (Fig. 12a) and jointly using JTV and OWXG (Fig. 12b). OWXG is applied to all pairs of the models with ${\gamma }_{01}={\gamma }_{02}=2000$⁠, ${\gamma }_{12}=200$ and ${\sigma }_{01}={\sigma }_{02}=-1$⁠, ${\sigma }_{12}=1$⁠. Data and model RMS errors and the number of iterations for separate and joint inversions are summarized in Table 1. As we can see from Fig. 12 (especially the resistivity model), joint inversion, by making use of the spatially complementary sensitivities, improves recovery of the features that are poorly resolved by the individual inversions due to the missing stations.

5.3 Example 3: joint inversion of MT and electrical resistivity tomography data

In this example we demonstrate how the proposed approach can be applied to the situation where the coupled inverse problems use radically different grids. It is inspired by field data from Furnas Volcano, Azores (Hogg et al. 2018), where a 1.5 km electrical resistivity tomography (ERT) profile was acquired across a fault-zone associated fumarole area together with audio-magnetotelluric (AMT) data. For ERT we use 2-D finite-element modelling and inversion implemented in pyGIMLI (Rücker et al. 2017). For AMT, the 3-D version of ModEM (Kelbert et al. 2014) is used. The true model is 2-D, consisting of a 1 $\mathrm{\Omega }\cdot$m ‘fault’, 10 $\mathrm{\Omega }\cdot$m ‘clay cap’ and a 1000 $\mathrm{\Omega }\cdot$m resistor embedded in 100 $\mathrm{\Omega }\cdot$m background with topography (Fig. 13). The ERT data are generated using 64 electrodes with 24 m spacing in the Wenner configuration and contaminated by Gaussian noise with the standard deviation of 2 per cent + 1 μV. The AMT data are calculated for four stations with 380 m spacing at nine frequencies in 1–100 Hz range and contaminated by Gaussian noise with 3 per cent standard deviation; only the off-diagonal components of the impedance tensor are inverted.

The AMT data are inverted using a $46\times 16\times 40$ model grid consisting of a core area with the horizontal cell size of 50 m $\times$ 50 m which is expanded by 5 padding cells in each horizontal direction, with the cell size increasing by a factor of 2.5. In the vertical direction, the cell size increases with depth from 5 m to 1.4 km. The model covariance smoothing is applied twice with a factor of 0.3 along the x- and z-directions and 0.9 along the y-direction. The ERT data are inverted using 2-D Delaunay triangulation grid with the cells expanding with depth starting from two edges between neighbouring electrodes at the surface. Both ERT and AMT inversions start with 100 $\mathrm{\Omega }\cdot$m homogeneous initial models.

For the joint inversion, using eq. (7) with $\mathcal{R}(\mathbf{u})=\Vert {\mathbf{u}}_0-{\mathbf{u}}_1{\Vert }^2=0$⁠, we get

so the reference models for the inverse subproblems are essentially the weighted average of the ERT and AMT models from the previous outer iteration. As both model vectors are logarithms of resistivity, this corresponds to the weighted geometric mean of the resistivities. We define $\mathbf{u}$ on a $340\times 130$ regular grid with a 5 m step and use bilinear interpolation between the cell centres for ${\mathbf{P}}_i$ and ${\mathbf{Q}}_i$⁠. In the AMT case, when mapping from the 3-D model to the 2-D $\mathbf{u}$⁠, we average along the y-direction over the core part of the grid, to minimize the effect of the model ‘fading’ away from the profile, and when mapping back, ‘broadcast’ along the y-direction over the whole grid.

At each outer iteration of the joint inversion, the ERT inversion is terminated after three Gauss–Newton iterations and the AMT inversion—after 10 NLCG iterations. The regularization parameters are set as follows for all the inversions: ${\hat{\alpha }}_0=40$ (ERT) and ${\hat{\alpha }}_1=1$ (AMT) and are kept constant. To analyse the effect of ERT/AMT weighting we test three values of ${\alpha }_0/{\alpha }_1$⁠: 0.1, 1 and 10.

The models resulting from the joint and separate inversions are shown in Fig. 14. As we can see from the separate inversion results, the near-surface part of the model is well resolved by ERT, but not by AMT, due to large spacing between the stations and the lack of high frequencies, while the deeper structure is resolved by AMT. Joint inversions achieve better reconstruction of the whole model, but at the expense of slightly higher RMS values in AMT: all joint inversions were terminated after 30 outer iterations without the RMS falling below 1.1. Fig. 15 shows the progress of the joint inversions in terms of RMS and $r_i$⁠. The high RMS values in AMT are largely due to the resistive artefact to the left of the ‘fault’, which is produced by the ERT inversion on the edge of its high-sensitivity region. Using, instead of scalar ${\alpha }_i$⁠, weighting matrices proportional to the integrated sensitivities should suppress such artefacts and further improve the joint inversion results. We see that the mapping between different grids based on a simple linear interpolation works reasonably well at least when the methods have comparable resolution lengths and the resistivity models are smooth. However, in more extreme cases accurate mapping of resistivities may require material averaging (see, e.g. Commer & Newman 2008) which will replace ${\mathbf{P}}_i$ and ${\mathbf{Q}}_i$ by nonlinear operators.

6 CONCLUSIONS

We have developed a flexible general framework for multiphysics joint inversion of any number of geophysical data sets. It is based on the variable splitting technique that introduces an auxiliary multiparameter model space in which minimization of the coupling and stabilizing functionals is carried out. By defining a mapping between this space and the model spaces via interpolation, completely different parametrizations can be applied to different models. Each of the inversion supbroblems and the coupling-regularization subproblem can be solved using a different optimization algorithm. For each subproblem, the linking term controlling the distance between the model and the corresponding auxiliary variable takes the form of a quadratic regularization with a reference model. Thus, any existing inversion code supporting such regularization can be integrated, without modifications, into this joint inversion framework. Using the developed framework, we have implemented 3-D joint inversion of magnetotelluric, seismic refraction and gravity data. We have implemented and tested on synthetic data different combinations of coupling constraints based on joint total variation, joint minimum support, cross-gradient and one-way cross-gradient, as well as an explicit relationship between the properties. The joint minimum support is a promising coupling regularization for multiphysics time-lapse inversions. The joint total variation and the cross-gradient are complementary, with JTV providing regularization and coupling edge-like features of the models by aligning their positions and XG coupling smoother features of the models by aligning gradient directions. The one-way cross-gradient can dramatically improve the effects of structural coupling in certain problems: in particular, for the gravity inversion it eliminates the non-uniqueness related to the lack of depth resolution. We expect that the developed approach will be especially useful for joint inversion with more complicated priors.

ACKNOWLEDGEMENTS

We would like to thank the editors, Colin Farquharson and two anonymous reviewers for their constructive comments that helped to improve the manuscript. This work has emanated from research supported in part by: a research grant to iCRAG from Science Foundation Ireland (SFI) under Grant Number 13/RC/2092 and co-funded under the European Regional Development Fund and by iCRAG industry partners; and a research grant to iCRAG2 from SFI under Grant Number 13/RC/2092_P2. Finally, we would like to thank Gary Egbert, Anna Kelbert and Naser Meqbel for making their ModEM code available to us.

DATA AVAILABILITY

The data underlying this paper will be shared on request to the corresponding author.

REFERENCES