# An efficient sequential strategy for realizing cross-gradient joint inversion: method and its application to 2-D cross borehole seismic traveltime and DC resistivity tomography

An efficient sequential strategy for realizing cross-gradient joint inversion: method and its... Summary Cross-gradient joint inversion that enforces structural similarity between different models has been widely utilized in jointly inverting different geophysical data types. However, it is a challenge to combine different geophysical inversion systems with the cross-gradient structural constraint into one joint inversion system because they may differ greatly in the model representation, forward modelling and inversion algorithm. Here we propose a new joint inversion strategy that can avoid this issue. Different models are separately inverted using the existing inversion packages and model structure similarity is only enforced through cross-gradient minimization between two models after each iteration. Although the data fitting and structural similarity enforcing processes are decoupled, our proposed strategy is still able to choose appropriate models to balance the trade-off between geophysical data fitting and structural similarity. This is realized by using model perturbations from separate data inversions to constrain the cross-gradient minimization process. We have tested this new strategy on 2-D cross borehole synthetic seismic traveltime and DC resistivity data sets. Compared to separate geophysical inversions, our proposed joint inversion strategy fits the separate data sets at comparable levels while at the same time resulting in a higher structural similarity between the velocity and resistivity models. Electrical resistivity tomography (ERT), Joint inversion, Seismic tomography 1 INTRODUCTION Geophysical imaging using different types of data such as seismic, gravity, magnetic, electrical and electromagnetic data has been widely used in characterizing subsurface structures at various scales. However, due to inherent data uncertainties, limited data coverage and limited data sensitivities, geophysical imaging using an individual geophysical data set produces a subsurface image with inherent ambiguity and uncertainty. Generally different geophysical methods with disparate data sets have complementary strengths in resolving the Earth structures. Therefore, to improve the model resolution and reduce its non-uniqueness, it is desirable to exploit a variety of data types in order to better image the subsurface structures (Gallardo & Meju 2011; Moorkamp et al. 2011; Bortolozo et al. 2014; Bennington et al. 2015). Joint geophysical inversion generally falls into three categories. The first category involves using multiple geophysical data sets to solve for common geophysical attributes. For example, both seismic shear wave arrival times and surface wave dispersion data are sensitive to shear wave velocity (Vs), thus, it is straightforward to combine the two into a common inversion for Vs (Zhang et al. 2014; Fang et al. 2016). Other examples include joint inversion of direct current (DC) and transient electromagnetic (TEM) data (Albouy et al. 2001), and joint inversion of magnetotelluric (MT) and controlled source electromagnetic data (Abubakar et al. 2009), which are all related to the resistivity structure of the earth. The second category refers to using different geophysical data to solve for different geophysical models linked through an empirical relationship. In this way, different geophysical data can be either directly or indirectly sensitive to common geophysical model parameters (Jupp & Vozoff 1975; Lines et al. 1988; Bosch 1999; Maceira & Ammon 2009). For example, by using the empirical relationship between seismic velocity and density, Maceira & Ammon (2009) jointly inverted gravity anomaly and surface wave data to determine the 3-D Vs, and implicitly density models of the Tarim basin in China. However, it is not always the case that an empirical relationship between different geophysical attributes can be found. For example, it is a challenge to find an empirical relationship between seismic velocity and electrical resistivity. The third category of joint inversion uses multiple geophysical data sets to invert for multiple geophysical models that cannot be related through an empirical relationship. In this case, the joint inversion is realized via structure similarity among different geophysical models (e.g. Haber & Oldenburg 1997; Gallardo & Meju 2003; Günther & Rücker 2006, Günther et al.2010; Bennington et al. 2015; Ronczka et al. 2017). For this category of joint inversion, the inverted models are required to satisfy their respective data while at the same time achieving structural similarity. To characterize the structural similarity, Gallardo & Meju (2003) proposed the cross-gradient constraint. The model structures are characterized by its gradients and thus the structure similarity can be determined by checking whether the respective model gradients are consistent. If cross-gradients of two models (defined as cross product of two gradients) are close to zero in a zone, it indicates they are either similar in structures or one model does not have detectable structures. Otherwise, the two model structures are different when the cross gradients are large. The cross-gradient joint inversion is designed to simultaneously minimize data residuals and model cross gradients (Gallardo & Meju 2003, 2004). By enforcing the model cross-gradients close to zero, the two models are thus similar in structure. Instead of using the cross-gradient constraint, Günther & Rücker (2006) used one model to modulate the roughness matrix via correlation, which is applied to regularize the inversion system for another model. After the first demonstration of cross-gradient joint inversion scheme by Gallardo & Meju (2003, 2004), the technique has been applied to 2-D or 3-D problems with various geophysical data types over a range of spatial scales. Fregoso & Gallardo (2009) extended the cross-gradient joint inversion to 3-D gravity and magnetics data. Moorkamp et al. (2011) presented a joint inversion framework for multiple data sets containing 3-D MT, gravity and seismic refraction data. Peng et al. (2012) jointly inverted receiver functions and MT data in order to determine the crustal and mantle structure beneath central Namche Barwa, eastern Himalayan syntaxis. Ogunbo & Zhang (2014) used seismic traveltime and TEM data to jointly invert for the near surface velocity and resistivity structure. Bennington et al. (2015) developed the joint inversion of 3-D earthquake arrival times and 2-D MT data based on cross-gradient structural constraint. Scholl et al. (2016) applied the cross-gradient constraint to the joint inversion of 3-D airborne EM and gravity data. Shi et al. (2017) developed cross-gradient joint inversion of 3-D seismic refraction and DC resistivity data for archaeological investigation. However, inversion systems related to different geophysical data sets may differ greatly in the model representation, forward modelling, and inversion algorithm. In addition, the disparate geophysical data sets may have significantly different sensitivities to various physical properties. As a result, it is a challenge to combine these systems of equation with the cross-gradient system in a common joint inversion. Um et al. (2014) used a sequential strategy where the separate geophysical data sets were inverted followed by cross-gradient minimization (hereafter referred as strategy I). The above strategy to some extent follows the joint inversion scheme of Günther & Rücker (2006), which iteratively and alternately inverts different geophysical data where the structures from one model are projected onto the roughness matrix used to regularize the inversion for the other model. The strategy of Um et al. (2014) has the advantage of easily realizing cross-gradient minimization while still fitting the separate data sets since the two minimization processes are totally decoupled. Geophysical models have good structural similarity after cross-gradient minimization, however, the fit to the geophysical data may be worse. In the following iteration, the structural similarity may worsen due to separate data inversions. In this study, we propose a new cross-gradient joint inversion strategy that keeps the advantage of the strategy used by Um et al. (2014) while better balancing the trade-off between data fitting and structural similarity (hereafter referred as strategy II). The key point of our proposed strategy is to use the model perturbations from separate data inversions to better constrain the cross-gradient minimization. We test our joint inversion strategy using synthetic 2-D cross borehole DC resistivity and seismic traveltime data. 2 METHODOLOGY In this section, we describe the new joint inversion strategy based on the cross-gradient constraint. For two 3-D models such as velocity (ms) and resistivity (mr), the cross-gradient t is given by   $${{\bf t}}(x,y,z) = \nabla {{{\bf m}}_{\rm{s}}}(x,y,z) \times \nabla {{{\bf m}}_{\rm{r}}}(x,y,z),$$ (1)with$$\,{{{\bf t}}_x} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial y}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial z}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial z}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial y}}$$, $${{{\bf t}}_y} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial z}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial x}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial x}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial z}}$$ and $${{{\bf t}}_z} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial x}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial y}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial y}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial x}}$$, respectively. For the 2-D models, the cross-gradient vector is represented as   $${{\bf t}}(x,y,z) = \nabla {{{\bf m}}_{\rm{s}}}(x,z) \times \nabla {{{\bf m}}_{\rm{r}}}(x,z).$$ (2) If the cross-gradient value is very small, this indicates that the two models are structurally similar. The goal of joint inversion under the cross-gradient constraint is to obtain two models that fit their respective data sets while achieving structural similarity (Gallardo & Meju 2003). An overall objective function for the cross-gradient joint inversion can be represented as follows:   $$\phi ({{{\bf m}}_{\rm{s}}}{\rm{,}}{{{\bf m}}_{\rm{r}}}) = {\phi _1}({{{\bf m}}_{\rm{s}}}) + {\phi _2}({{{\bf m}}_{\rm{r}}}) + \left\| {t({{{\bf m}}_{\rm{s}}}{\rm{,}}{{{\bf m}}_{\rm{r}}})} \right\|.$$ (3) For the case of the joint inversion of seismic traveltime and DC resistivity data, the objective function related to minimizing seismic traveltime residuals is represented as   $$\phi ({{{\bf m}}_{\rm{s}}}) = \left\| {{{{\bf W}}_\sigma }(f({{{\bf m}}_{\rm{s}}}) - {{{\bf d}}_{\rm{s}}})} \right\|_2^2 + {\alpha _{\rm{s}}}^2\left\| {{{{\bf L}}_{\rm{s}}}{{{\bf m}}_{\rm{s}}}} \right\|_2^2 + {\beta _{\rm{s}}}^2\left\| {{{\bf I}}{{{\bf m}}_{\rm{s}}}} \right\|_2^2,$$ (4) where f(ms) is the theoretical traveltime calculated using the current velocity model ms, ds is the observed traveltime, Wσ is a diagonal matrix with the diagonal value as the inverse of data uncertainty, Ls is the smoothing operator, αs is the smoothing weight, I is the identity matrix, and βs is the damping weight, respectively. For calculating synthetic seismic traveltimes, first we discretize the velocity model into rectangular grid nodes and then calculate the 2-D traveltime field by using the finite difference method to solve the eikonal equation (Sambridge 1990; Rawlinson & Sambridge 2005; Hassouna & Farag 2007). For seismic traveltime tomography, we solve the objective function shown in eq. (4) by the least squares method. The smoothing and damping parameters are carefully selected through the trade-off analysis between model smoothness and data fitting (Aster et al. 2013). For DC resistivity tomography, the objective function is   $$\phi ({{{\bf m}}_{\rm{r}}}) = \frac{1}{2}\left\| {{{{\bf W}}_{\rm{d}}}\left[ {f({{{\bf m}}_{\rm{r}}}) - {{{\bf d}}_{\rm{r}}}} \right]} \right\|_2^2 + \frac{{{\beta _{\rm{r}}}}}{2}\left\| {{{{\bf W}}_{\rm{m}}}({{{\bf m}}_{\rm{r}}} - {{{\bf m}}_{\rm{r}}}{{_,}_{{\rm{ref}}}})} \right\|_2^2 ,\!\!\!\!$$ (5) where φ(mr) is the theoretically calculated potential according to the current resistivity model mr, dr is the observed potential, mref is the reference resistivity model, Wd is the weight for the data misfit, Wm is the weight for fitting the reference model, βr is the weight parameter to balance the trade-off between the data misfit $$\| {{{{\bf W}}_{\rm{d}}}[ {f( {{{{\bf m}}_{\rm{r}}}} ) - {{{\bf d}}_{\rm{r}}}} ]} \|_2^2$$ and the regularization term of $$\| {{{{\bf W}}_{\rm{m}}}( {{{{\bf m}}_{\rm{r}}} - {{{\bf m}}_{{\rm{r,ref}}}}} )} \|_2^2$$ (Tikhonov & Arsenin 1977). For calculating synthetic electrical potentials, we use the finite difference method to solve the electrostatic field equation to calculate the electric potential on the rectangular grid (Pidlisecky et al. 2007). For DC resistivity tomography, we applied the Gauss–Newton method to minimize the objective function shown in eq. (5). The Jacobian matrix and Hessian matrix are calculated based on the adjoint method (Mcgillivray & Oldenburg 1990; Pidlisecky et al.2007; Wu & Xu 2000). For the cross-gradient structural constraint, the related inversion system for minimizing the cross-gradients between two models (i.e. the third term in the objective function in eq. (3)) is as follows:   $$\left[ {\begin{array}{c@{\quad}c} {\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {\frac{{\partial {{{\bf t}}_y}}}{{\partial {{\bf m}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {\Delta {{{\bf m}}_{\rm{s}}}}\\ {\Delta {{{\bf m}}_{\rm{r}}}} \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{@{}*{1}{c}@{}} { - {{{\bf t}}_x}}\\ { - {{{\bf t}}_y}}\\ { - {{{\bf t}}_z}} \end{array}} \right],$$ (6)where $$\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}},\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{s}}}}},\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}$$ are the partial derivatives of cross-gradient t with respect to ms in the X, Y, and Z directions, $$\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}},\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}},\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}$$ are the partial derivatives of the cross-gradient t with respect to mr in X, Y, Z directions, tx, ty, tz are X-, Y- and Z-components of the cross-gradient vector between velocity and resistivity models from separate inversions. To minimize the objective function for cross-gradient joint inversion of seismic traveltime and DC resistivity data represented in eq. (3), it is extremely challenging to construct a single inversion system combining seismic traveltime tomography, DC resistivity tomography and cross-gradient structural inversion. This is because seismic traveltime tomography and DC resistivity tomography use very different model representations and inversion strategies. From the point of computer coding, it is difficult to merge the existing individual software packages because they may use different coding languages. To overcome this, Um et al. (2014) used an iterative inversion strategy to avoid combining the different geophysical inversions and the cross-gradient structural constraint. Essentially, the inversion consists of two loops. Initially, the two geophysical models are driven to have structural similarity via the cross-gradient objective function based on eq. (6), without regard for misfit to the geophysical data sets. Following this step, inversion of the separate geophysical data sets is carried out in order to obtain acceptable fits to the respective data measurements, independent of cross-gradient minimization. These iterative steps are repeated until the two geophysical models converge. In this way, the geophysical data constraints are completely decoupled and independent of the cross-gradient structural constraint. The advantage of this joint inversion strategy is that the separate inversion package can still be used. However, the disadvantage of this inversion strategy is that because geophysical data inversion is decoupled from the cross-gradient structure constraint, it is difficult to balance the trade-off between geophysical data fit and structural similarity. Here we propose a new joint inversion strategy that overcomes the disadvantages of Um et al. (2014). From separate seismic traveltime and DC resistivity inversions, we can determine the velocity model perturbation $$\Delta {{\bf m}}_{\rm{s}}^0$$ and resistivity model perturbation $$\Delta {{\bf m}}_{\rm{r}}^0$$. Thus, we can subsequently get the updated ms and mr. Instead of applying the cross-gradient constraint to the updated ms and mr based on eq. (6), to better balance the data misfit and model structural similarity, we modify the eq. (6) by incorporating model perturbations from the separate data inversions into the system of equations representing the cross-gradient constraint, as follows:   $$\left[ {\begin{array}{c@{\quad}c} {{w_{\rm{s}}}{{\bf I}}}&{{\bf 0}}\\ {{\bf 0}}&{{w_{\rm{r}}}{{\bf I}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}} \end{array}} \right]\left[ \begin{array}{@{}l@{}} \Delta {{{\bf m}}_s}\\ \Delta {{{\bf m}}_r} \end{array} \right] = \left[ {\begin{array}{@{}*{1}{c}@{}} {{w_{\rm{s}}}\Delta {{\bf m}}_{\rm{s}}^0}\\ {{w_{\rm{r}}}\Delta {{\bf m}}_{\rm{r}}^0}\\ { - {w_{\rm{t}}}{{{\bf t}}_x}}\\ { - {w_{\rm{t}}}{{{\bf t}}_y}}\\ { - {w_{\rm{t}}}{{{\bf t}}_z}} \end{array}} \right],$$ (7) where ws and wr are weights for velocity and resistivity model updates from separate inversions, wt is the weight parameter for the cross-gradient structural constraint. Δms is the new velocity model update and Δmr is the new resistivity model update. In the case of two dimensions, the eq. (7) can be simplified as   $$\left[ {\begin{array}{c@{\quad}c} {{w_{\rm{s}}}{{\bf I}}}&{{\bf 0}}\\ {{\bf 0}}&{{w_{\rm{r}}}{{\bf I}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{\bf m}_s}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_r}}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {\Delta {{{\bf m}}_{\rm{s}}}}\\ {\Delta {{{\bf m}}_{\rm{r}}}} \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{@{}*{1}{c}@{}} {{w_{\rm{s}}}\Delta {{\bf m}}_{\rm{s}}^0}\\ {{w_{\rm{r}}}\Delta {{\bf m}}_{\rm{r}}^0}\\ { - {w_{\rm{t}}}{{{\bf t}}_y}} \end{array}} \right].$$ (8) Fig. 1 shows the flow chart of the new joint inversion strategy, summarized as follows: Input initial velocity and resistivity models and calculate their cross-gradients. Conduct separate data inversions to minimize objective functions (4) and (5) corresponding to seismic traveltime data and DC resistivity data, respectively. Model pertuabations $$\Delta {{\bf m}}_{\rm{r}}^0$$ and $$\Delta {{\bf m}}_{\rm{s}}^0$$ are obtained. Set up eq. (7) or (8) using cross-gradients calculated in step (1) and model updates from step (2). By carrying out the inversion of eq. (7) or (8), we can solve for Δms and Δmr. Update the velocity and resistivity models using Δms and Δmr from step (3). Using the updated geophysical models, calculate the cross-gradients and the seismic traveltime and DC resistivity data residuals. Judge whether the data fitting and cross-gradient values meet the iteration termination criteria. If so, output the results; otherwise go back to step (2) to start another iteration. Figure 1. View largeDownload slide Proposed joint inversion flow for cross-gradient joint inversion. Figure 1. View largeDownload slide Proposed joint inversion flow for cross-gradient joint inversion. Compared to the fully coupled joint inversion system used by Gallardo & Meju (2003, 2004) as well as by Bennington et al. (2015), our proposed joint inversion scheme has the data inversion and structural constraint decoupled, which makes it very easy to combine any type of geophysical data into a cross-gradient joint inversion system. In eq. (7), it is very important to choose appropriate balance weight factors ws and wr because velocity and resistivity models have very different magnitudes. For seismic traveltime tomography, the slowness is typically inverted for instead of velocity. The magnitude of slowness is ∼10–3 s m-1 and the associated slowness update has the magnitude of ∼10–4 s m-1. For the DC resistivity tomography, the resistivity model is generally represented in a natural logarithm scale. Its magnitude is ∼101 and the associated model update magnitude is ∼100. Therefore, to ensure the two model updates are equivalent, weight factor ws should be a factor of 4 greater than weight factor wr. For the parameter wt, we use the L-curve criterion to select the appropriate weight to balance the model updates from the separate inversions (indirectly data fitting) and the application of the cross-gradient constraint (Aster et al. 2013). 3 SYNTHETIC DATA TEST To validate the effectiveness of our proposed joint inversion strategy, we setup a 2-D cross borehole observation system for seismic traveltime and DC resistivity tomography. The two geophysical methods are widely used in geological engineering exploration to understand the near surface, and each method has certain limitations (Rao et al. 2016; Lu et al. 2017). For example, seismic traveltime tomography has difficulty accurately resolving segregated low velocity anomalies. Similarly, DC resistivity tomography has difficulty resolving high resistivity anomalies. In this study, synthetic velocity and resistivity models have two resistivity and velocity anomalies having the same size and location (Figs 2a and c). With synthetic traveltime and electrical potential data, twenty iterations of the inversion are conducted to achieve model convergence. Fig. 2(b) shows the separately inverted velocity model using synthetic traveltimes. It can be seen that both high and low velocity anomalies are well resolved, but there are also smaller artefacts surrounding the two anomalies (Fig. 2b). Compared to the recovered velocity model, the separately inverted resistivity model resolves low and high resistivity anomalies with fewer artefacts but with poorly imaged boundaries (Fig. 2d). From the above synthetic tests, it can be clearly seen that seismic traveltime tomography and DC resistivity tomography have different strengths in recovering model anomalies. Figure 2. View largeDownload slide Synthetic and inverted velocity and resistivity models for cross borehole seismic traveltime and DC resistivity tomography. (a) Synthetic velocity model with traveltime observation system. The velocity model is constructed by adding two anomalies of 2200 and 1800 m s−1 to the background velocity of 2000 m s−1. There are two 64 m deep boreholes with the separation of 44 m between them. 60 receivers with an interval of 1 m were put into each borehole and 40 receivers were put onto the surface. Each borehole has 30 shots with an interval of 2 m. Each shot is recorded by all receivers in the borehole and on the surface. (b) Inverted velocity model from seismic traveltime tomography. (c) Synthetic resistivity model with DC resistivity observation system. The background resistivity is 100 Ω m. The two resistivity anomalies (1000 and 10 Ω m) have the same size as the two velocity anomalies. 30 electrodes with an interval of 2 m were put into each borehole and 21 electrodes were put onto the surface. The power supply mode is that for any two electrodes with power supply, all the rest of the electrodes collect the electrical potential. (d) Inverted resistivity model from DC resistivity tomography. Figure 2. View largeDownload slide Synthetic and inverted velocity and resistivity models for cross borehole seismic traveltime and DC resistivity tomography. (a) Synthetic velocity model with traveltime observation system. The velocity model is constructed by adding two anomalies of 2200 and 1800 m s−1 to the background velocity of 2000 m s−1. There are two 64 m deep boreholes with the separation of 44 m between them. 60 receivers with an interval of 1 m were put into each borehole and 40 receivers were put onto the surface. Each borehole has 30 shots with an interval of 2 m. Each shot is recorded by all receivers in the borehole and on the surface. (b) Inverted velocity model from seismic traveltime tomography. (c) Synthetic resistivity model with DC resistivity observation system. The background resistivity is 100 Ω m. The two resistivity anomalies (1000 and 10 Ω m) have the same size as the two velocity anomalies. 30 electrodes with an interval of 2 m were put into each borehole and 21 electrodes were put onto the surface. The power supply mode is that for any two electrodes with power supply, all the rest of the electrodes collect the electrical potential. (d) Inverted resistivity model from DC resistivity tomography. For comparison, we conducted cross-gradient joint inversions using the strategy of Um et al. (2014) and our proposed strategy in this study. Same as separate geophysical data inversions, twenty iterations are also performed to achieve the model convergence. For our joint inversion strategy, it is important to choose appropriate weighting parameter wt in eq. (7) or (8) in order to balance geophysical data fitting and the structural similarity constraint. Here we used the L-curve analysis to find the appropriate wt by varying the parameter from 10–10 to 1010 with an interval of 101 (Calvetti et al. 2000; Fig. 3). From the L-curve analysis, it can be seen that the optimal wt value is around 105 for the joint inversion. However, for the first 5 iterations, we choose a relatively smaller wt = 102 to weight more on the geophysical data fitting to update the models. In this way, the inverted models would not likely be trapped in structures of the initial models. After two models gain necessary structures from the geophysical data fitting, for iterations 6–20 we set wt = 105 to further enforce structural similarity. Figure 3. View largeDownload slide L-curve analysis for choosing appropriate weighting parameter for balancing data fitting and cross-gradient structure constraint. Figure 3. View largeDownload slide L-curve analysis for choosing appropriate weighting parameter for balancing data fitting and cross-gradient structure constraint. Fig. 4 shows the comparison of velocity and resistivity models after the first iteration of separate and joint inversions. Overall, for the first iteration, joint inversion with our proposed strategy behaves almost the same as separate data inversions, with very similar images of models (Fig. 4) and data residual as well as cross-gradient RMS values (Table 1). This is because initial iterations have the relative weighting factor wt set to be relatively low, thus, the cross-gradient constraint has very little effect on model updates. For joint inversion with the Um et al. (2014) strategy, due to the cross-gradient structural constraint, the model cross-gradients become smaller compared to separate inversion; the seismic traveltime fitting becomes worse; and the electrical potential fitting becomes better. This is expected because the inside loop of the Um et al. (2014) strategy only considers structure consistency. While the structural similarity between the two models is enforced, it is not guaranteed the two models can also fit the data well. As shown in Fig. 4 and Table 1, one model fits data better and the other model fits data worse after the inside loop of structural constraint. It may also occur that both models fit the data worse. In comparison, the inverted velocity and resistivity models using the joint inversion strategy proposed in this study can better balance data fitting and structure consistency (Table 1). Figure 4. View largeDownload slide Comparison of velocity (first row) and resistivity (second row) models after the first iteration of separate inversions (a and d), joint inversion using the strategy of Um et al. (2014) (b and e), and joint inversion using the strategy proposed in this study (c and f). Figure 4. View largeDownload slide Comparison of velocity (first row) and resistivity (second row) models after the first iteration of separate inversions (a and d), joint inversion using the strategy of Um et al. (2014) (b and e), and joint inversion using the strategy proposed in this study (c and f). Table 1. Data residual and cross-gradient RMS values for different inversion strategies after the first iteration.   Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  1.046  1.069  1.046  Electrical potential (mv)  21.098  20.918  21.098  Cross-gradient (10–6)  1.049  0.759  1.049    Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  1.046  1.069  1.046  Electrical potential (mv)  21.098  20.918  21.098  Cross-gradient (10–6)  1.049  0.759  1.049  View Large Fig. 5 and Table 2 show inverted velocity and resistivity models and the associated statistics after 20 iterations of the inversion, respectively. Overall, the performances of the joint and separate inversions after 20 iterations are similar to those after the first iteration (Table 2). Compared with the separately inverted velocity model (Fig. 2b), both jointly inverted velocity models significantly eliminate artefacts originally appearing in the separately inverted velocity model (Fig. 5). This is because the resistivity model is better resolved with fewer artefacts compared to separately inverted velocity model (Fig. 2), which can help to constrain the velocity model through the cross-gradient constraint. As a result, the two anomalies are well resolved without spurious artefacts surrounding the velocity anomalies. The model difference RMS values between inverted and true velocity models are 0.4910 and 0.5120 m s−1 for our proposed strategy and the Um et al. (2014) strategy, respectively, while it is 0.5134 m s−1 for the separate seismic traveltime inversion, indicating our proposed strategy performs the best. It can also be clearly seen that our proposed strategy better recovers the two velocity anomalies that closely resemble the true anomalies (Fig. 5). For the resistivity model, both high and low resistivity anomalies are much better resolved relative to the separate inversions. The resistivity model difference RMS values are 1.6142 and 1.6234 Ω m for our proposed strategy and the Um et al. (2014) strategy, respectively, while the RMS value is 1.6275 Ω m for the separate DC resistivity inversion. In comparison, our proposed joint inversion strategy better resolves the shapes of high and low resistivity anomalies, which are very close to true anomalies. This is because structures of two corresponding velocity anomalies help to resolve resistivity anomalies via the cross-gradient constraint (Figs 2 and 5). Figure 5. View largeDownload slide Comparison of velocity and resistivity models from separate inversions as well as joint inversions with different strategies after 20th iteration. (a) and (b) show velocity and resistivity models from separate inversions, (c) and (d) show velocity and resistivity models from joint inversion with the strategy of Um et al. (2014), and (e) and (f) show velocity and resistivity models from joint inversion using the proposed strategy, respectively. Figure 5. View largeDownload slide Comparison of velocity and resistivity models from separate inversions as well as joint inversions with different strategies after 20th iteration. (a) and (b) show velocity and resistivity models from separate inversions, (c) and (d) show velocity and resistivity models from joint inversion with the strategy of Um et al. (2014), and (e) and (f) show velocity and resistivity models from joint inversion using the proposed strategy, respectively. Table 2. Data residual and cross-gradient RMS values for different inversion strategies after 20 iterations.   Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  0.799  0.853  0.791  Electrical potential (mv)  10.454  10.906  11.431  Cross-gradient (10–6)  3.572  1.127  1.386    Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  0.799  0.853  0.791  Electrical potential (mv)  10.454  10.906  11.431  Cross-gradient (10–6)  3.572  1.127  1.386  View Large In terms of statistical analysis, velocity and resistivity models from separate inversions have the worse structural consistency, as indicated by the larger cross-gradient RMS value of 3.572. For cross-gradient joint inversions with either strategy, the two models are more structurally similar (RMS values of 1.127 for the Um et al. (2014) strategy and 1.386 for our proposed strategy; Table 2). For the geophysical data fitting, the joint inversion using the Um et al. (2014) strategy results in larger misfits for both seismic traveltime and electrical potential data relative to separate inversions (Table 2). In comparison, joint inversion with our proposed strategy results in similar seismic traveltime data fitting and slightly worse electrical potential data fitting than separate inversions (Table 2). Finally, the two joint inversion strategies produce velocity and resistivity models fitting the separate data sets at similar levels. Thus it is difficult to judge whether it is advantageous to use our proposed strategy over the Um et al. (2014) strategy solely from the point of data fitting. However, when we compare final velocity and resistivity models from the joint inversion with different strategies, it is clear that they better approximate the true models with our proposed strategy (Fig. 5). Compared to the joint inversion strategy of Um et al. (2014), our joint inversion strategy better resolves the anomalies’ shapes and amplitudes for both velocity and resistivity models. This can be seen from the comparison of the distribution of cross gradients for the final velocity and resistivity models from the three inversion schemes examined in this study (Fig. 6). The model cross gradients are close to zero around anomaly zones for the models obtained via our joint inversion strategy. Figure 6. View largeDownload slide Distribution of cross gradients for different inversions. (a) Separate inversions; (b) joint inversion using the strategy of Um et al. (2014); (c) Joint inversion using the proposed strategy. Figure 6. View largeDownload slide Distribution of cross gradients for different inversions. (a) Separate inversions; (b) joint inversion using the strategy of Um et al. (2014); (c) Joint inversion using the proposed strategy. Fig. 7 shows the data RMS residuals with progressive iterations for seismic traveltimes and electrical potentials as well as model cross-gradient RMS values for separate and joint inversions. From these curves, we can better see the differences between joint inversions with different strategies and the separate inversions. Similar to the separate inversions, the data RMS residual curves from the joint inversion with our proposed strategy decrease smoothly and monotonically over progressive iterations. In comparison, joint inversion with the Um et al. (2014) strategy may have the unstable issue during iterations. For example, as indicated in Fig. 7, at 13th iteration the cross-gradient minimization produces a resistivity model that bears similar structures to the velocity model (Fig. 8), but it poorly fits the data. For the following iteration of the separate DC resistivity tomography, the resistivity model has to be modified greatly to fit the data, thus leading to poorer structural similarity to the velocity model (Fig. 8). For this synthetic test, it actually takes several iterations for the inverted resistivity model to fit the data at comparable levels before the data fitting gets considerably worse (Fig. 7). This may not always happen but could be a potential issue for the Um et al. (2014) strategy. In comparison, our proposed strategy has a much smoother decrease in data residual because the model updates are partially constrained by geophysical data. Figure 7. View largeDownload slide RMS values of (a) cross gradients between velocity and resistivity models, (b) seismic traveltime residuals and (c) electrical potential residuals for separate and joint inversions. Figure 7. View largeDownload slide RMS values of (a) cross gradients between velocity and resistivity models, (b) seismic traveltime residuals and (c) electrical potential residuals for separate and joint inversions. Figure 8. View largeDownload slide Comparison of cross-gradient and resistivity models after 13th and 14th iterations of joint inversion using the strategy of Um et al. (2014) strategy. (a) and (b) show resistivity models after 13th and 14th iterations, respectively. (c and (d) show cross-gradients after 13th and 14th iterations, respectively. Figure 8. View largeDownload slide Comparison of cross-gradient and resistivity models after 13th and 14th iterations of joint inversion using the strategy of Um et al. (2014) strategy. (a) and (b) show resistivity models after 13th and 14th iterations, respectively. (c and (d) show cross-gradients after 13th and 14th iterations, respectively. In the above synthetic test, velocity and resistivity models have anomalies at the same position and of the same size (Fig. 2). To further test the ability of joint inversion strategy proposed in this study, we design another synthetic test with velocity and resistivity models having similar structures only in some regions of the model (Fig. 9). Anomaly V1 in the velocity model and anomaly R1 in the resistivity model have exactly the same structures. In comparison, anomaly V2 in the velocity model does not have the counterpart in the resistivity model. Similarly, anomaly R2 in the resistivity model does not have the counterpart in the velocity model either. For anomaly V3 in the velocity model and R3 in the resistivity model, their structures are only partially similar. Using the joint inversion strategy proposed in this study, the anomaly R1 is better resolved through the common structural constraint by the anomaly V1 than the separate inversion (Fig. 10). In comparison, for anomalies V2 and R2 that only exist in one model, they are also well recovered by cross-gradient joint inversion, which does not bring artefacts to the corresponding region in the other model (Fig. 10). Especially for anomaly V2, the related smearing effect is reduced in joint inversion. For anomalies V3 and R3 that only partially overlap, each anomaly is still well recovered by the joint inversion and the resistivity anomaly R3 is better resolved than the separate DC resistivity inversion (Fig. 10). Overall, as seen in the previous synthetic test, the artefacts around velocity anomalies are much smaller in the velocity model from joint inversion (Figs 9 and 10). Compared to the separate inversions, resistivity and velocity models from joint inversion are more structurally similar with a decrease of 56 per cent in the cross-gradient RMS value. For the geophysical data fitting, seismic traveltime RMS residual only slightly increases 3 per cent (from 0.832 to 0855 ms) while electrical potential RMS value decreases 26 per cent (from 16.69 to 12.39 ms). This test shows that the cross-gradient joint inversion strategy proposed in this study also works well even when two models have different structures in some regions. This is because when one model has anomalies and the other has no anomalies in some regions, the cross-gradients are still zero. As a result, cross-gradient structural constraint may not force structures in one model into another if data do not require so. Figure 9. View largeDownload slide Comparison of synthetic and inverted velocity and resistivity models from the separate and joint inversions using the strategy proposed in this study. The observation systems are the same as those in Fig. 2. (a) Synthetic velocity model with three anomalies: V1 (2200 m s−1), V2 (2200 m s−1) and V3 (1800 m s−1). The background velocity is 2000 m s−1. (b) Synthetic resistivity model with three anomalies: R1 (1000 Ω m), R2 (10 Ω m) and R3 (10 Ω m). The background resistivity is 100 Ω m. (c) Separately inverted velocity model. (d) Separately inverted resistivity model. (e) Jointly inverted velocity model. (f) Jointly inverted resistivity model. Figure 9. View largeDownload slide Comparison of synthetic and inverted velocity and resistivity models from the separate and joint inversions using the strategy proposed in this study. The observation systems are the same as those in Fig. 2. (a) Synthetic velocity model with three anomalies: V1 (2200 m s−1), V2 (2200 m s−1) and V3 (1800 m s−1). The background velocity is 2000 m s−1. (b) Synthetic resistivity model with three anomalies: R1 (1000 Ω m), R2 (10 Ω m) and R3 (10 Ω m). The background resistivity is 100 Ω m. (c) Separately inverted velocity model. (d) Separately inverted resistivity model. (e) Jointly inverted velocity model. (f) Jointly inverted resistivity model. Figure 10. View largeDownload slide Comparison of velocity and resistivity differences between inverted and true models for separate and joint inversions with our proposed strategy. (a) Velocity differences for separate inversion. (b) Velocity differences for joint inversion. (c) Resistivity differences for separate inversion. (d) Resistivity differences for joint inversion. Figure 10. View largeDownload slide Comparison of velocity and resistivity differences between inverted and true models for separate and joint inversions with our proposed strategy. (a) Velocity differences for separate inversion. (b) Velocity differences for joint inversion. (c) Resistivity differences for separate inversion. (d) Resistivity differences for joint inversion. 4 CONCLUSIONS Based on the joint inversion strategy of Um et al. (2014), we have proposed a new strategy to efficiently realize the cross-gradient joint inversion. In this strategy, different data sets are inverted separately and thus the existing inversion packages can be used in the joint inversion with trivial modifications. The model updates from separate inversions are then combined with the cross-gradient minimization system to make the two models similar in structure. Compared to the strategy of Um et al. (2014), our proposed strategy can better balance structural similarity and data fitting (via model updates from the separate inversions). Compared to the complete joint inversion system used by Bennington et al. (2015), the new joint inversion strategy avoids coupling between different systems and, can therefore, better account for varying data sensitivity amplitudes and different data coverage for different inversion systems. At the same time, it can produce models that are constrained by both geophysical data and cross-gradient minimization, efficiently realizing the cross-gradient structure constraint. We tested our proposed strategy on 2-D cross borehole seismic traveltime and DC resistivity tomography systems. Compared to the separate inversions, our strategy can better recover the resistivity anomalies while suppressing artefacts surrounding velocity anomalies. The synthetic test indicates that our proposed joint inversion strategy can be easily used to realize the cross-gradient joint inversion constraint while at the same time providing a better balance between data fitting and structural similarity than the strategy of Um et al. (2014). Acknowledgements We are grateful for the availability of the MSFM package (Hassouna & Farag 2007) for seismic traveltime calculation and the RESINVM3D package (Pidlisecky et al.2007) for DC resistivity tomography. We acknowledge the constructive comments from two anonymous reviewers that are very helpful for improving the manuscript. This work is supported by the National Natural Science Foundation of China (41704039) and by the Open Projects of State Key Laboratory of Coal Resources and Safe Mining, CUMT (SKLCRSM16KF09). REFERENCES Abubakar A., Li M., Liu J., Habashy T.M., 2009. Simultaneous joint inversion of MT and CSEM data using a multiplicative cost function, in SEG Technical Program Expanded Abstracts 2009, pp. 719– 723, Society of Exploration Geophysicists. Albouy Y., Andrieux P., Rakotondrasoa G., Ritz M., Descloitres M., Join J.L., Rasolomanana E., 2001. Mapping coastal aquifers by joint inversion of DC and TEM soundings—three case histories, Ground Water , 39( 1), 87– 97. https://doi.org/10.1111/j.1745-6584.2001.tb00354.x Google Scholar CrossRef Search ADS   Aster R.C., Borchers B., Thurber C.H. 2013. Parameter Estimation and Inverse Problems , 2nd edn, Elsevier Inc. Bennington N.L., Zhang H., Thurber C.H., Bedrosian P.A., 2015. Joint inversion of seismic and magnetotelluric data in the Parkfield Region of California using the normalized cross-gradient constraint, Pure appl. Geophys. , 172( 5), 1033– 1052. https://doi.org/10.1007/s00024-014-1002-9 Google Scholar CrossRef Search ADS   Bortolozo C.A., Couto M.A., Porsani J.L., Almeida E.R., Santos F.A.M., 2014. Geoelectrical characterization using joint inversion of VES/TEM data: a case study in Paraná Sedimentary Basin, São Paulo State, Brazil, J. Appl. Geophys. , 111: 33– 46. https://doi.org/10.1016/j.jappgeo.2014.09.009 Google Scholar CrossRef Search ADS   Bosch M., 1999. Lithologic tomography: from plural geophysical data to lithology estimation, J. geophys. Res. , 104( B1), 749– 766. https://doi.org/10.1029/1998JB900014 Google Scholar CrossRef Search ADS   Calvetti D., Morigi S., Reichel L., Sgallari F., 2000. Tikhonov regularization and the L-curve for large discrete ill-posed problems, J. Comput. Appl. Math. , 123( 1–2), 423– 446. https://doi.org/10.1016/S0377-0427(00)00414-3 Google Scholar CrossRef Search ADS   Fang H., Zhang H., Yao H., Allam A., Zigone D., Ben-Zion Y., Thurber C., van der Hilst R.D., 2016. A new algorithm for three-dimensional joint inversion of body wave and surface wave data and its application to the Southern California plate boundary region, J. geophys. Res. , 121( 5), 3557– 3569. Google Scholar CrossRef Search ADS   Fregoso E., Gallardo L.A., 2009. Cross-gradients joint 3D inversion with applications to gravity and magnetic data, Geophysics , 74( 4), L31– L42. https://doi.org/10.1190/1.3119263 Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data, Geophys. Res. Lett. , 30( 13), doi:10.1029/2003GL017370. Gallardo L.A., Meju M.A., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints, J. geophys. Res. , 109( B3), doi:10.1029/2003JB002716. Gallardo L.A., Meju M.A., 2011. Structure-coupled multiphysics imaging in geophysical sciences, Rev. Geophys. , 49( 1), doi: 10.1029/2010RG000330 Günther T., Rücker C., 2006. A new joint inversion approach applied to the combined tomography of DC resistivity and seismic refraction data, in Symposium on the Application of Geophysics to Engineering and Environmental Problems 2006, pp. 1196– 1202, doi.org/10.4133/1.2923578 Günther T., Dlugosch R., Holland R., Yaramanci U., 2010, Aquifer characterization using coupled inversion of DC/IP and MRS data on a hydrogeophysical test-site, in SAGEEP, Symposium on the Application of Geophysics to Engineering and Environmental Problems 2010 , pp. 302– 307. Haber E., Oldenburg D., 1997. Joint inversion: a structural approach, Inverse Probl. , 13( 1), 63– 77. https://doi.org/10.1088/0266-5611/13/1/006 Google Scholar CrossRef Search ADS   Hassouna M.S., Farag A.A., 2007. Multi-stencils fast marching methods: a highly accurate solution to the eikonal equation on cartesian domains, in IIEEE Trans. Pattern Anal. Mach. Intell. , 29( 9), 1563– 1574. Google Scholar CrossRef Search ADS   Jupp D.L.B., Vozoff K., 1975. Stable iterative methods for the inversion of geophysical data, Geophys. J. Int. , 42( 3), 957– 976. https://doi.org/10.1111/j.1365-246X.1975.tb06461.x Google Scholar CrossRef Search ADS   Lines L.R., Schultz A.K., Treitel S., 1988. Cooperative inversion of geophysical data, Geophysics , 53( 1), 8– 20. https://doi.org/10.1190/1.1442403 Google Scholar CrossRef Search ADS   Lu D.B., Wang F., Chen X.D., Qu J., Wang H., 2017. An improved ERT approach for the investigation of subsurface structures, Pure appl. Geophys. , 174( 1), 375– 386. https://doi.org/10.1007/s00024-016-1386-9 Google Scholar CrossRef Search ADS   Maceira M., Ammon C.J., 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure, J. geophys. Res. , 114( B2), doi:10.1029/2007JB005157. https://doi.org/10.1029/2007JB005157 McGillivray P.R., Oldenburg D.W., 1990. Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study, Geophys. Prospect. , 38( 5), 499– 524. https://doi.org/10.1111/j.1365-2478.1990.tb01859.x Google Scholar CrossRef Search ADS   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, Geophys. J. Int. , 184( 1), 477– 493. https://doi.org/10.1111/j.1365-246X.2010.04856.x Google Scholar CrossRef Search ADS   Ogunbo J.N., Zhang J., 2014. Joint seismic traveltime and TEM inversion for near surface imaging, in SEG Technical Program Expanded Abstracts 2014, pp. 2104– 2108, Society of Exploration Geophysicists. Peng M., Tan H.D., Jiang M., Wang W., Li Q.-Q., Zhang L.-S., 2012. Joint inversion of receiver functions and magnetotelluric data: application to crustal and mantle structure beneath central Namche Barwa, eastern Himalayan syntaxis, Chinese J. Geophys (in Chinese) , 55( 7), 2281– 2291. Pidlisecky A., Haber E., Knight R., 2007. RESINVM3D: a 3D resistivity inversion package, Geophysics , 72( 2), H1– H10. https://doi.org/10.1190/1.2402499 Google Scholar CrossRef Search ADS   Rao Y., Wang Y., Chen S., Wang J., 2016. Crosshole seismic tomography with cross-firing geometry, Geophysics , 81( 4), R139– R146. https://doi.org/10.1190/geo2015-0677.1 Google Scholar CrossRef Search ADS   Rawlinson N., Sambridge M., 2005. The fast marching method: an effective tool for tomographic imaging and tracking multiple phases in complex layered media, Explor. Geophys. , 36( 4), 341– 350. https://doi.org/10.1071/EG05341 Google Scholar CrossRef Search ADS   Ronczka M., Hellman K., Günther T., Wisén R., Dahlin T. ( 2017), Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory, Solid Earth , 8( 3), 671– 682. https://doi.org/10.5194/se-8-671-2017 Google Scholar CrossRef Search ADS   Sambridge M.S., 1990. Non-linear arrival time inversion: constraining velocity anomalies by seeking smooth models in 3-D, Geophys. J. Int. , 102( 3), 653– 677. https://doi.org/10.1111/j.1365-246X.1990.tb04588.x Google Scholar CrossRef Search ADS   Scholl C., Neumann J., Watts M.D., Hallinan S., Mule S., 2016. Geologically constrained 2D and 3D airborne EM inversion through cross-gradient regularization and multi-grid efficiency, in ASEG Extended Abstracts 2016 , pp. 54– 67. https://doi.org/10.1071/ASEG2016ab229 Shi Z., Hobbs R.W., Moorkamp M., Tian G., Jiang L., 2017. 3-D cross-gradient joint inversion of seismic refraction and DC resistivity data, J. Appl. Geophys. , doi:10.1016/j.jappgeo.2017.04.008. https://doi.org/10.1016/j.jappgeo.2017.04.008 Tikhonov A.N., Arsenin V.I.A.. 1977. Solutions of Ill-posed Problems , Vh Winston. Um E.S., Commer M., Newman G.A. 2014. A strategy for coupled 3D imaging of large-scale seismic and electromagnetic data sets: application to subsalt imaging, Geophysics , 79( 3), ID1– ID13. https://doi.org/10.1190/geo2013-0053.1 Google Scholar CrossRef Search ADS   Wu X.-P., Xu G.-M., 2000. Study on 3D resistivity inversion using conjugate gradient method, China J. Geophys . ( in Chinese), 43( 3), 450– 458. Google Scholar CrossRef Search ADS   Zhang H., Maceira M., Roux P., Thurber C., 2014. Joint inversion of body-wave arrival times and surface-wave dispersion for three-dimensional seismic structure around SAFOD, Pure appl. Geophys. , 171( 11), 3013– 3022. https://doi.org/10.1007/s00024-014-0806-y Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# An efficient sequential strategy for realizing cross-gradient joint inversion: method and its application to 2-D cross borehole seismic traveltime and DC resistivity tomography

, Volume 213 (2) – May 1, 2018
12 pages

Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy026
Publisher site
See Article on Publisher Site

### Abstract

Summary Cross-gradient joint inversion that enforces structural similarity between different models has been widely utilized in jointly inverting different geophysical data types. However, it is a challenge to combine different geophysical inversion systems with the cross-gradient structural constraint into one joint inversion system because they may differ greatly in the model representation, forward modelling and inversion algorithm. Here we propose a new joint inversion strategy that can avoid this issue. Different models are separately inverted using the existing inversion packages and model structure similarity is only enforced through cross-gradient minimization between two models after each iteration. Although the data fitting and structural similarity enforcing processes are decoupled, our proposed strategy is still able to choose appropriate models to balance the trade-off between geophysical data fitting and structural similarity. This is realized by using model perturbations from separate data inversions to constrain the cross-gradient minimization process. We have tested this new strategy on 2-D cross borehole synthetic seismic traveltime and DC resistivity data sets. Compared to separate geophysical inversions, our proposed joint inversion strategy fits the separate data sets at comparable levels while at the same time resulting in a higher structural similarity between the velocity and resistivity models. Electrical resistivity tomography (ERT), Joint inversion, Seismic tomography 1 INTRODUCTION Geophysical imaging using different types of data such as seismic, gravity, magnetic, electrical and electromagnetic data has been widely used in characterizing subsurface structures at various scales. However, due to inherent data uncertainties, limited data coverage and limited data sensitivities, geophysical imaging using an individual geophysical data set produces a subsurface image with inherent ambiguity and uncertainty. Generally different geophysical methods with disparate data sets have complementary strengths in resolving the Earth structures. Therefore, to improve the model resolution and reduce its non-uniqueness, it is desirable to exploit a variety of data types in order to better image the subsurface structures (Gallardo & Meju 2011; Moorkamp et al. 2011; Bortolozo et al. 2014; Bennington et al. 2015). Joint geophysical inversion generally falls into three categories. The first category involves using multiple geophysical data sets to solve for common geophysical attributes. For example, both seismic shear wave arrival times and surface wave dispersion data are sensitive to shear wave velocity (Vs), thus, it is straightforward to combine the two into a common inversion for Vs (Zhang et al. 2014; Fang et al. 2016). Other examples include joint inversion of direct current (DC) and transient electromagnetic (TEM) data (Albouy et al. 2001), and joint inversion of magnetotelluric (MT) and controlled source electromagnetic data (Abubakar et al. 2009), which are all related to the resistivity structure of the earth. The second category refers to using different geophysical data to solve for different geophysical models linked through an empirical relationship. In this way, different geophysical data can be either directly or indirectly sensitive to common geophysical model parameters (Jupp & Vozoff 1975; Lines et al. 1988; Bosch 1999; Maceira & Ammon 2009). For example, by using the empirical relationship between seismic velocity and density, Maceira & Ammon (2009) jointly inverted gravity anomaly and surface wave data to determine the 3-D Vs, and implicitly density models of the Tarim basin in China. However, it is not always the case that an empirical relationship between different geophysical attributes can be found. For example, it is a challenge to find an empirical relationship between seismic velocity and electrical resistivity. The third category of joint inversion uses multiple geophysical data sets to invert for multiple geophysical models that cannot be related through an empirical relationship. In this case, the joint inversion is realized via structure similarity among different geophysical models (e.g. Haber & Oldenburg 1997; Gallardo & Meju 2003; Günther & Rücker 2006, Günther et al.2010; Bennington et al. 2015; Ronczka et al. 2017). For this category of joint inversion, the inverted models are required to satisfy their respective data while at the same time achieving structural similarity. To characterize the structural similarity, Gallardo & Meju (2003) proposed the cross-gradient constraint. The model structures are characterized by its gradients and thus the structure similarity can be determined by checking whether the respective model gradients are consistent. If cross-gradients of two models (defined as cross product of two gradients) are close to zero in a zone, it indicates they are either similar in structures or one model does not have detectable structures. Otherwise, the two model structures are different when the cross gradients are large. The cross-gradient joint inversion is designed to simultaneously minimize data residuals and model cross gradients (Gallardo & Meju 2003, 2004). By enforcing the model cross-gradients close to zero, the two models are thus similar in structure. Instead of using the cross-gradient constraint, Günther & Rücker (2006) used one model to modulate the roughness matrix via correlation, which is applied to regularize the inversion system for another model. After the first demonstration of cross-gradient joint inversion scheme by Gallardo & Meju (2003, 2004), the technique has been applied to 2-D or 3-D problems with various geophysical data types over a range of spatial scales. Fregoso & Gallardo (2009) extended the cross-gradient joint inversion to 3-D gravity and magnetics data. Moorkamp et al. (2011) presented a joint inversion framework for multiple data sets containing 3-D MT, gravity and seismic refraction data. Peng et al. (2012) jointly inverted receiver functions and MT data in order to determine the crustal and mantle structure beneath central Namche Barwa, eastern Himalayan syntaxis. Ogunbo & Zhang (2014) used seismic traveltime and TEM data to jointly invert for the near surface velocity and resistivity structure. Bennington et al. (2015) developed the joint inversion of 3-D earthquake arrival times and 2-D MT data based on cross-gradient structural constraint. Scholl et al. (2016) applied the cross-gradient constraint to the joint inversion of 3-D airborne EM and gravity data. Shi et al. (2017) developed cross-gradient joint inversion of 3-D seismic refraction and DC resistivity data for archaeological investigation. However, inversion systems related to different geophysical data sets may differ greatly in the model representation, forward modelling, and inversion algorithm. In addition, the disparate geophysical data sets may have significantly different sensitivities to various physical properties. As a result, it is a challenge to combine these systems of equation with the cross-gradient system in a common joint inversion. Um et al. (2014) used a sequential strategy where the separate geophysical data sets were inverted followed by cross-gradient minimization (hereafter referred as strategy I). The above strategy to some extent follows the joint inversion scheme of Günther & Rücker (2006), which iteratively and alternately inverts different geophysical data where the structures from one model are projected onto the roughness matrix used to regularize the inversion for the other model. The strategy of Um et al. (2014) has the advantage of easily realizing cross-gradient minimization while still fitting the separate data sets since the two minimization processes are totally decoupled. Geophysical models have good structural similarity after cross-gradient minimization, however, the fit to the geophysical data may be worse. In the following iteration, the structural similarity may worsen due to separate data inversions. In this study, we propose a new cross-gradient joint inversion strategy that keeps the advantage of the strategy used by Um et al. (2014) while better balancing the trade-off between data fitting and structural similarity (hereafter referred as strategy II). The key point of our proposed strategy is to use the model perturbations from separate data inversions to better constrain the cross-gradient minimization. We test our joint inversion strategy using synthetic 2-D cross borehole DC resistivity and seismic traveltime data. 2 METHODOLOGY In this section, we describe the new joint inversion strategy based on the cross-gradient constraint. For two 3-D models such as velocity (ms) and resistivity (mr), the cross-gradient t is given by   $${{\bf t}}(x,y,z) = \nabla {{{\bf m}}_{\rm{s}}}(x,y,z) \times \nabla {{{\bf m}}_{\rm{r}}}(x,y,z),$$ (1)with$$\,{{{\bf t}}_x} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial y}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial z}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial z}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial y}}$$, $${{{\bf t}}_y} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial z}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial x}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial x}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial z}}$$ and $${{{\bf t}}_z} = \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial x}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial y}} - \frac{{\partial {{{\bf m}}_{\rm{r}}}}}{{\partial y}} \cdot \frac{{\partial {{{\bf m}}_{\rm{s}}}}}{{\partial x}}$$, respectively. For the 2-D models, the cross-gradient vector is represented as   $${{\bf t}}(x,y,z) = \nabla {{{\bf m}}_{\rm{s}}}(x,z) \times \nabla {{{\bf m}}_{\rm{r}}}(x,z).$$ (2) If the cross-gradient value is very small, this indicates that the two models are structurally similar. The goal of joint inversion under the cross-gradient constraint is to obtain two models that fit their respective data sets while achieving structural similarity (Gallardo & Meju 2003). An overall objective function for the cross-gradient joint inversion can be represented as follows:   $$\phi ({{{\bf m}}_{\rm{s}}}{\rm{,}}{{{\bf m}}_{\rm{r}}}) = {\phi _1}({{{\bf m}}_{\rm{s}}}) + {\phi _2}({{{\bf m}}_{\rm{r}}}) + \left\| {t({{{\bf m}}_{\rm{s}}}{\rm{,}}{{{\bf m}}_{\rm{r}}})} \right\|.$$ (3) For the case of the joint inversion of seismic traveltime and DC resistivity data, the objective function related to minimizing seismic traveltime residuals is represented as   $$\phi ({{{\bf m}}_{\rm{s}}}) = \left\| {{{{\bf W}}_\sigma }(f({{{\bf m}}_{\rm{s}}}) - {{{\bf d}}_{\rm{s}}})} \right\|_2^2 + {\alpha _{\rm{s}}}^2\left\| {{{{\bf L}}_{\rm{s}}}{{{\bf m}}_{\rm{s}}}} \right\|_2^2 + {\beta _{\rm{s}}}^2\left\| {{{\bf I}}{{{\bf m}}_{\rm{s}}}} \right\|_2^2,$$ (4) where f(ms) is the theoretical traveltime calculated using the current velocity model ms, ds is the observed traveltime, Wσ is a diagonal matrix with the diagonal value as the inverse of data uncertainty, Ls is the smoothing operator, αs is the smoothing weight, I is the identity matrix, and βs is the damping weight, respectively. For calculating synthetic seismic traveltimes, first we discretize the velocity model into rectangular grid nodes and then calculate the 2-D traveltime field by using the finite difference method to solve the eikonal equation (Sambridge 1990; Rawlinson & Sambridge 2005; Hassouna & Farag 2007). For seismic traveltime tomography, we solve the objective function shown in eq. (4) by the least squares method. The smoothing and damping parameters are carefully selected through the trade-off analysis between model smoothness and data fitting (Aster et al. 2013). For DC resistivity tomography, the objective function is   $$\phi ({{{\bf m}}_{\rm{r}}}) = \frac{1}{2}\left\| {{{{\bf W}}_{\rm{d}}}\left[ {f({{{\bf m}}_{\rm{r}}}) - {{{\bf d}}_{\rm{r}}}} \right]} \right\|_2^2 + \frac{{{\beta _{\rm{r}}}}}{2}\left\| {{{{\bf W}}_{\rm{m}}}({{{\bf m}}_{\rm{r}}} - {{{\bf m}}_{\rm{r}}}{{_,}_{{\rm{ref}}}})} \right\|_2^2 ,\!\!\!\!$$ (5) where φ(mr) is the theoretically calculated potential according to the current resistivity model mr, dr is the observed potential, mref is the reference resistivity model, Wd is the weight for the data misfit, Wm is the weight for fitting the reference model, βr is the weight parameter to balance the trade-off between the data misfit $$\| {{{{\bf W}}_{\rm{d}}}[ {f( {{{{\bf m}}_{\rm{r}}}} ) - {{{\bf d}}_{\rm{r}}}} ]} \|_2^2$$ and the regularization term of $$\| {{{{\bf W}}_{\rm{m}}}( {{{{\bf m}}_{\rm{r}}} - {{{\bf m}}_{{\rm{r,ref}}}}} )} \|_2^2$$ (Tikhonov & Arsenin 1977). For calculating synthetic electrical potentials, we use the finite difference method to solve the electrostatic field equation to calculate the electric potential on the rectangular grid (Pidlisecky et al. 2007). For DC resistivity tomography, we applied the Gauss–Newton method to minimize the objective function shown in eq. (5). The Jacobian matrix and Hessian matrix are calculated based on the adjoint method (Mcgillivray & Oldenburg 1990; Pidlisecky et al.2007; Wu & Xu 2000). For the cross-gradient structural constraint, the related inversion system for minimizing the cross-gradients between two models (i.e. the third term in the objective function in eq. (3)) is as follows:   $$\left[ {\begin{array}{c@{\quad}c} {\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {\frac{{\partial {{{\bf t}}_y}}}{{\partial {{\bf m}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {\Delta {{{\bf m}}_{\rm{s}}}}\\ {\Delta {{{\bf m}}_{\rm{r}}}} \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{@{}*{1}{c}@{}} { - {{{\bf t}}_x}}\\ { - {{{\bf t}}_y}}\\ { - {{{\bf t}}_z}} \end{array}} \right],$$ (6)where $$\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}},\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{s}}}}},\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}$$ are the partial derivatives of cross-gradient t with respect to ms in the X, Y, and Z directions, $$\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}},\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}},\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}$$ are the partial derivatives of the cross-gradient t with respect to mr in X, Y, Z directions, tx, ty, tz are X-, Y- and Z-components of the cross-gradient vector between velocity and resistivity models from separate inversions. To minimize the objective function for cross-gradient joint inversion of seismic traveltime and DC resistivity data represented in eq. (3), it is extremely challenging to construct a single inversion system combining seismic traveltime tomography, DC resistivity tomography and cross-gradient structural inversion. This is because seismic traveltime tomography and DC resistivity tomography use very different model representations and inversion strategies. From the point of computer coding, it is difficult to merge the existing individual software packages because they may use different coding languages. To overcome this, Um et al. (2014) used an iterative inversion strategy to avoid combining the different geophysical inversions and the cross-gradient structural constraint. Essentially, the inversion consists of two loops. Initially, the two geophysical models are driven to have structural similarity via the cross-gradient objective function based on eq. (6), without regard for misfit to the geophysical data sets. Following this step, inversion of the separate geophysical data sets is carried out in order to obtain acceptable fits to the respective data measurements, independent of cross-gradient minimization. These iterative steps are repeated until the two geophysical models converge. In this way, the geophysical data constraints are completely decoupled and independent of the cross-gradient structural constraint. The advantage of this joint inversion strategy is that the separate inversion package can still be used. However, the disadvantage of this inversion strategy is that because geophysical data inversion is decoupled from the cross-gradient structure constraint, it is difficult to balance the trade-off between geophysical data fit and structural similarity. Here we propose a new joint inversion strategy that overcomes the disadvantages of Um et al. (2014). From separate seismic traveltime and DC resistivity inversions, we can determine the velocity model perturbation $$\Delta {{\bf m}}_{\rm{s}}^0$$ and resistivity model perturbation $$\Delta {{\bf m}}_{\rm{r}}^0$$. Thus, we can subsequently get the updated ms and mr. Instead of applying the cross-gradient constraint to the updated ms and mr based on eq. (6), to better balance the data misfit and model structural similarity, we modify the eq. (6) by incorporating model perturbations from the separate data inversions into the system of equations representing the cross-gradient constraint, as follows:   $$\left[ {\begin{array}{c@{\quad}c} {{w_{\rm{s}}}{{\bf I}}}&{{\bf 0}}\\ {{\bf 0}}&{{w_{\rm{r}}}{{\bf I}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_x}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_{\rm{r}}}}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{s}}}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_z}}}{{\partial {{{\bf m}}_{\rm{r}}}}}} \end{array}} \right]\left[ \begin{array}{@{}l@{}} \Delta {{{\bf m}}_s}\\ \Delta {{{\bf m}}_r} \end{array} \right] = \left[ {\begin{array}{@{}*{1}{c}@{}} {{w_{\rm{s}}}\Delta {{\bf m}}_{\rm{s}}^0}\\ {{w_{\rm{r}}}\Delta {{\bf m}}_{\rm{r}}^0}\\ { - {w_{\rm{t}}}{{{\bf t}}_x}}\\ { - {w_{\rm{t}}}{{{\bf t}}_y}}\\ { - {w_{\rm{t}}}{{{\bf t}}_z}} \end{array}} \right],$$ (7) where ws and wr are weights for velocity and resistivity model updates from separate inversions, wt is the weight parameter for the cross-gradient structural constraint. Δms is the new velocity model update and Δmr is the new resistivity model update. In the case of two dimensions, the eq. (7) can be simplified as   $$\left[ {\begin{array}{c@{\quad}c} {{w_{\rm{s}}}{{\bf I}}}&{{\bf 0}}\\ {{\bf 0}}&{{w_{\rm{r}}}{{\bf I}}}\\ {{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{\bf m}_s}}}}&{{w_{\rm{t}}}\frac{{\partial {{{\bf t}}_y}}}{{\partial {{{\bf m}}_r}}}} \end{array}} \right]\left[ {\begin{array}{@{}*{1}{c}@{}} {\Delta {{{\bf m}}_{\rm{s}}}}\\ {\Delta {{{\bf m}}_{\rm{r}}}} \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{@{}*{1}{c}@{}} {{w_{\rm{s}}}\Delta {{\bf m}}_{\rm{s}}^0}\\ {{w_{\rm{r}}}\Delta {{\bf m}}_{\rm{r}}^0}\\ { - {w_{\rm{t}}}{{{\bf t}}_y}} \end{array}} \right].$$ (8) Fig. 1 shows the flow chart of the new joint inversion strategy, summarized as follows: Input initial velocity and resistivity models and calculate their cross-gradients. Conduct separate data inversions to minimize objective functions (4) and (5) corresponding to seismic traveltime data and DC resistivity data, respectively. Model pertuabations $$\Delta {{\bf m}}_{\rm{r}}^0$$ and $$\Delta {{\bf m}}_{\rm{s}}^0$$ are obtained. Set up eq. (7) or (8) using cross-gradients calculated in step (1) and model updates from step (2). By carrying out the inversion of eq. (7) or (8), we can solve for Δms and Δmr. Update the velocity and resistivity models using Δms and Δmr from step (3). Using the updated geophysical models, calculate the cross-gradients and the seismic traveltime and DC resistivity data residuals. Judge whether the data fitting and cross-gradient values meet the iteration termination criteria. If so, output the results; otherwise go back to step (2) to start another iteration. Figure 1. View largeDownload slide Proposed joint inversion flow for cross-gradient joint inversion. Figure 1. View largeDownload slide Proposed joint inversion flow for cross-gradient joint inversion. Compared to the fully coupled joint inversion system used by Gallardo & Meju (2003, 2004) as well as by Bennington et al. (2015), our proposed joint inversion scheme has the data inversion and structural constraint decoupled, which makes it very easy to combine any type of geophysical data into a cross-gradient joint inversion system. In eq. (7), it is very important to choose appropriate balance weight factors ws and wr because velocity and resistivity models have very different magnitudes. For seismic traveltime tomography, the slowness is typically inverted for instead of velocity. The magnitude of slowness is ∼10–3 s m-1 and the associated slowness update has the magnitude of ∼10–4 s m-1. For the DC resistivity tomography, the resistivity model is generally represented in a natural logarithm scale. Its magnitude is ∼101 and the associated model update magnitude is ∼100. Therefore, to ensure the two model updates are equivalent, weight factor ws should be a factor of 4 greater than weight factor wr. For the parameter wt, we use the L-curve criterion to select the appropriate weight to balance the model updates from the separate inversions (indirectly data fitting) and the application of the cross-gradient constraint (Aster et al. 2013). 3 SYNTHETIC DATA TEST To validate the effectiveness of our proposed joint inversion strategy, we setup a 2-D cross borehole observation system for seismic traveltime and DC resistivity tomography. The two geophysical methods are widely used in geological engineering exploration to understand the near surface, and each method has certain limitations (Rao et al. 2016; Lu et al. 2017). For example, seismic traveltime tomography has difficulty accurately resolving segregated low velocity anomalies. Similarly, DC resistivity tomography has difficulty resolving high resistivity anomalies. In this study, synthetic velocity and resistivity models have two resistivity and velocity anomalies having the same size and location (Figs 2a and c). With synthetic traveltime and electrical potential data, twenty iterations of the inversion are conducted to achieve model convergence. Fig. 2(b) shows the separately inverted velocity model using synthetic traveltimes. It can be seen that both high and low velocity anomalies are well resolved, but there are also smaller artefacts surrounding the two anomalies (Fig. 2b). Compared to the recovered velocity model, the separately inverted resistivity model resolves low and high resistivity anomalies with fewer artefacts but with poorly imaged boundaries (Fig. 2d). From the above synthetic tests, it can be clearly seen that seismic traveltime tomography and DC resistivity tomography have different strengths in recovering model anomalies. Figure 2. View largeDownload slide Synthetic and inverted velocity and resistivity models for cross borehole seismic traveltime and DC resistivity tomography. (a) Synthetic velocity model with traveltime observation system. The velocity model is constructed by adding two anomalies of 2200 and 1800 m s−1 to the background velocity of 2000 m s−1. There are two 64 m deep boreholes with the separation of 44 m between them. 60 receivers with an interval of 1 m were put into each borehole and 40 receivers were put onto the surface. Each borehole has 30 shots with an interval of 2 m. Each shot is recorded by all receivers in the borehole and on the surface. (b) Inverted velocity model from seismic traveltime tomography. (c) Synthetic resistivity model with DC resistivity observation system. The background resistivity is 100 Ω m. The two resistivity anomalies (1000 and 10 Ω m) have the same size as the two velocity anomalies. 30 electrodes with an interval of 2 m were put into each borehole and 21 electrodes were put onto the surface. The power supply mode is that for any two electrodes with power supply, all the rest of the electrodes collect the electrical potential. (d) Inverted resistivity model from DC resistivity tomography. Figure 2. View largeDownload slide Synthetic and inverted velocity and resistivity models for cross borehole seismic traveltime and DC resistivity tomography. (a) Synthetic velocity model with traveltime observation system. The velocity model is constructed by adding two anomalies of 2200 and 1800 m s−1 to the background velocity of 2000 m s−1. There are two 64 m deep boreholes with the separation of 44 m between them. 60 receivers with an interval of 1 m were put into each borehole and 40 receivers were put onto the surface. Each borehole has 30 shots with an interval of 2 m. Each shot is recorded by all receivers in the borehole and on the surface. (b) Inverted velocity model from seismic traveltime tomography. (c) Synthetic resistivity model with DC resistivity observation system. The background resistivity is 100 Ω m. The two resistivity anomalies (1000 and 10 Ω m) have the same size as the two velocity anomalies. 30 electrodes with an interval of 2 m were put into each borehole and 21 electrodes were put onto the surface. The power supply mode is that for any two electrodes with power supply, all the rest of the electrodes collect the electrical potential. (d) Inverted resistivity model from DC resistivity tomography. For comparison, we conducted cross-gradient joint inversions using the strategy of Um et al. (2014) and our proposed strategy in this study. Same as separate geophysical data inversions, twenty iterations are also performed to achieve the model convergence. For our joint inversion strategy, it is important to choose appropriate weighting parameter wt in eq. (7) or (8) in order to balance geophysical data fitting and the structural similarity constraint. Here we used the L-curve analysis to find the appropriate wt by varying the parameter from 10–10 to 1010 with an interval of 101 (Calvetti et al. 2000; Fig. 3). From the L-curve analysis, it can be seen that the optimal wt value is around 105 for the joint inversion. However, for the first 5 iterations, we choose a relatively smaller wt = 102 to weight more on the geophysical data fitting to update the models. In this way, the inverted models would not likely be trapped in structures of the initial models. After two models gain necessary structures from the geophysical data fitting, for iterations 6–20 we set wt = 105 to further enforce structural similarity. Figure 3. View largeDownload slide L-curve analysis for choosing appropriate weighting parameter for balancing data fitting and cross-gradient structure constraint. Figure 3. View largeDownload slide L-curve analysis for choosing appropriate weighting parameter for balancing data fitting and cross-gradient structure constraint. Fig. 4 shows the comparison of velocity and resistivity models after the first iteration of separate and joint inversions. Overall, for the first iteration, joint inversion with our proposed strategy behaves almost the same as separate data inversions, with very similar images of models (Fig. 4) and data residual as well as cross-gradient RMS values (Table 1). This is because initial iterations have the relative weighting factor wt set to be relatively low, thus, the cross-gradient constraint has very little effect on model updates. For joint inversion with the Um et al. (2014) strategy, due to the cross-gradient structural constraint, the model cross-gradients become smaller compared to separate inversion; the seismic traveltime fitting becomes worse; and the electrical potential fitting becomes better. This is expected because the inside loop of the Um et al. (2014) strategy only considers structure consistency. While the structural similarity between the two models is enforced, it is not guaranteed the two models can also fit the data well. As shown in Fig. 4 and Table 1, one model fits data better and the other model fits data worse after the inside loop of structural constraint. It may also occur that both models fit the data worse. In comparison, the inverted velocity and resistivity models using the joint inversion strategy proposed in this study can better balance data fitting and structure consistency (Table 1). Figure 4. View largeDownload slide Comparison of velocity (first row) and resistivity (second row) models after the first iteration of separate inversions (a and d), joint inversion using the strategy of Um et al. (2014) (b and e), and joint inversion using the strategy proposed in this study (c and f). Figure 4. View largeDownload slide Comparison of velocity (first row) and resistivity (second row) models after the first iteration of separate inversions (a and d), joint inversion using the strategy of Um et al. (2014) (b and e), and joint inversion using the strategy proposed in this study (c and f). Table 1. Data residual and cross-gradient RMS values for different inversion strategies after the first iteration.   Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  1.046  1.069  1.046  Electrical potential (mv)  21.098  20.918  21.098  Cross-gradient (10–6)  1.049  0.759  1.049    Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  1.046  1.069  1.046  Electrical potential (mv)  21.098  20.918  21.098  Cross-gradient (10–6)  1.049  0.759  1.049  View Large Fig. 5 and Table 2 show inverted velocity and resistivity models and the associated statistics after 20 iterations of the inversion, respectively. Overall, the performances of the joint and separate inversions after 20 iterations are similar to those after the first iteration (Table 2). Compared with the separately inverted velocity model (Fig. 2b), both jointly inverted velocity models significantly eliminate artefacts originally appearing in the separately inverted velocity model (Fig. 5). This is because the resistivity model is better resolved with fewer artefacts compared to separately inverted velocity model (Fig. 2), which can help to constrain the velocity model through the cross-gradient constraint. As a result, the two anomalies are well resolved without spurious artefacts surrounding the velocity anomalies. The model difference RMS values between inverted and true velocity models are 0.4910 and 0.5120 m s−1 for our proposed strategy and the Um et al. (2014) strategy, respectively, while it is 0.5134 m s−1 for the separate seismic traveltime inversion, indicating our proposed strategy performs the best. It can also be clearly seen that our proposed strategy better recovers the two velocity anomalies that closely resemble the true anomalies (Fig. 5). For the resistivity model, both high and low resistivity anomalies are much better resolved relative to the separate inversions. The resistivity model difference RMS values are 1.6142 and 1.6234 Ω m for our proposed strategy and the Um et al. (2014) strategy, respectively, while the RMS value is 1.6275 Ω m for the separate DC resistivity inversion. In comparison, our proposed joint inversion strategy better resolves the shapes of high and low resistivity anomalies, which are very close to true anomalies. This is because structures of two corresponding velocity anomalies help to resolve resistivity anomalies via the cross-gradient constraint (Figs 2 and 5). Figure 5. View largeDownload slide Comparison of velocity and resistivity models from separate inversions as well as joint inversions with different strategies after 20th iteration. (a) and (b) show velocity and resistivity models from separate inversions, (c) and (d) show velocity and resistivity models from joint inversion with the strategy of Um et al. (2014), and (e) and (f) show velocity and resistivity models from joint inversion using the proposed strategy, respectively. Figure 5. View largeDownload slide Comparison of velocity and resistivity models from separate inversions as well as joint inversions with different strategies after 20th iteration. (a) and (b) show velocity and resistivity models from separate inversions, (c) and (d) show velocity and resistivity models from joint inversion with the strategy of Um et al. (2014), and (e) and (f) show velocity and resistivity models from joint inversion using the proposed strategy, respectively. Table 2. Data residual and cross-gradient RMS values for different inversion strategies after 20 iterations.   Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  0.799  0.853  0.791  Electrical potential (mv)  10.454  10.906  11.431  Cross-gradient (10–6)  3.572  1.127  1.386    Separate inversion  Joint inversion strategy I  Joint inversion strategy II  Seismic traveltime (ms)  0.799  0.853  0.791  Electrical potential (mv)  10.454  10.906  11.431  Cross-gradient (10–6)  3.572  1.127  1.386  View Large In terms of statistical analysis, velocity and resistivity models from separate inversions have the worse structural consistency, as indicated by the larger cross-gradient RMS value of 3.572. For cross-gradient joint inversions with either strategy, the two models are more structurally similar (RMS values of 1.127 for the Um et al. (2014) strategy and 1.386 for our proposed strategy; Table 2). For the geophysical data fitting, the joint inversion using the Um et al. (2014) strategy results in larger misfits for both seismic traveltime and electrical potential data relative to separate inversions (Table 2). In comparison, joint inversion with our proposed strategy results in similar seismic traveltime data fitting and slightly worse electrical potential data fitting than separate inversions (Table 2). Finally, the two joint inversion strategies produce velocity and resistivity models fitting the separate data sets at similar levels. Thus it is difficult to judge whether it is advantageous to use our proposed strategy over the Um et al. (2014) strategy solely from the point of data fitting. However, when we compare final velocity and resistivity models from the joint inversion with different strategies, it is clear that they better approximate the true models with our proposed strategy (Fig. 5). Compared to the joint inversion strategy of Um et al. (2014), our joint inversion strategy better resolves the anomalies’ shapes and amplitudes for both velocity and resistivity models. This can be seen from the comparison of the distribution of cross gradients for the final velocity and resistivity models from the three inversion schemes examined in this study (Fig. 6). The model cross gradients are close to zero around anomaly zones for the models obtained via our joint inversion strategy. Figure 6. View largeDownload slide Distribution of cross gradients for different inversions. (a) Separate inversions; (b) joint inversion using the strategy of Um et al. (2014); (c) Joint inversion using the proposed strategy. Figure 6. View largeDownload slide Distribution of cross gradients for different inversions. (a) Separate inversions; (b) joint inversion using the strategy of Um et al. (2014); (c) Joint inversion using the proposed strategy. Fig. 7 shows the data RMS residuals with progressive iterations for seismic traveltimes and electrical potentials as well as model cross-gradient RMS values for separate and joint inversions. From these curves, we can better see the differences between joint inversions with different strategies and the separate inversions. Similar to the separate inversions, the data RMS residual curves from the joint inversion with our proposed strategy decrease smoothly and monotonically over progressive iterations. In comparison, joint inversion with the Um et al. (2014) strategy may have the unstable issue during iterations. For example, as indicated in Fig. 7, at 13th iteration the cross-gradient minimization produces a resistivity model that bears similar structures to the velocity model (Fig. 8), but it poorly fits the data. For the following iteration of the separate DC resistivity tomography, the resistivity model has to be modified greatly to fit the data, thus leading to poorer structural similarity to the velocity model (Fig. 8). For this synthetic test, it actually takes several iterations for the inverted resistivity model to fit the data at comparable levels before the data fitting gets considerably worse (Fig. 7). This may not always happen but could be a potential issue for the Um et al. (2014) strategy. In comparison, our proposed strategy has a much smoother decrease in data residual because the model updates are partially constrained by geophysical data. Figure 7. View largeDownload slide RMS values of (a) cross gradients between velocity and resistivity models, (b) seismic traveltime residuals and (c) electrical potential residuals for separate and joint inversions. Figure 7. View largeDownload slide RMS values of (a) cross gradients between velocity and resistivity models, (b) seismic traveltime residuals and (c) electrical potential residuals for separate and joint inversions. Figure 8. View largeDownload slide Comparison of cross-gradient and resistivity models after 13th and 14th iterations of joint inversion using the strategy of Um et al. (2014) strategy. (a) and (b) show resistivity models after 13th and 14th iterations, respectively. (c and (d) show cross-gradients after 13th and 14th iterations, respectively. Figure 8. View largeDownload slide Comparison of cross-gradient and resistivity models after 13th and 14th iterations of joint inversion using the strategy of Um et al. (2014) strategy. (a) and (b) show resistivity models after 13th and 14th iterations, respectively. (c and (d) show cross-gradients after 13th and 14th iterations, respectively. In the above synthetic test, velocity and resistivity models have anomalies at the same position and of the same size (Fig. 2). To further test the ability of joint inversion strategy proposed in this study, we design another synthetic test with velocity and resistivity models having similar structures only in some regions of the model (Fig. 9). Anomaly V1 in the velocity model and anomaly R1 in the resistivity model have exactly the same structures. In comparison, anomaly V2 in the velocity model does not have the counterpart in the resistivity model. Similarly, anomaly R2 in the resistivity model does not have the counterpart in the velocity model either. For anomaly V3 in the velocity model and R3 in the resistivity model, their structures are only partially similar. Using the joint inversion strategy proposed in this study, the anomaly R1 is better resolved through the common structural constraint by the anomaly V1 than the separate inversion (Fig. 10). In comparison, for anomalies V2 and R2 that only exist in one model, they are also well recovered by cross-gradient joint inversion, which does not bring artefacts to the corresponding region in the other model (Fig. 10). Especially for anomaly V2, the related smearing effect is reduced in joint inversion. For anomalies V3 and R3 that only partially overlap, each anomaly is still well recovered by the joint inversion and the resistivity anomaly R3 is better resolved than the separate DC resistivity inversion (Fig. 10). Overall, as seen in the previous synthetic test, the artefacts around velocity anomalies are much smaller in the velocity model from joint inversion (Figs 9 and 10). Compared to the separate inversions, resistivity and velocity models from joint inversion are more structurally similar with a decrease of 56 per cent in the cross-gradient RMS value. For the geophysical data fitting, seismic traveltime RMS residual only slightly increases 3 per cent (from 0.832 to 0855 ms) while electrical potential RMS value decreases 26 per cent (from 16.69 to 12.39 ms). This test shows that the cross-gradient joint inversion strategy proposed in this study also works well even when two models have different structures in some regions. This is because when one model has anomalies and the other has no anomalies in some regions, the cross-gradients are still zero. As a result, cross-gradient structural constraint may not force structures in one model into another if data do not require so. Figure 9. View largeDownload slide Comparison of synthetic and inverted velocity and resistivity models from the separate and joint inversions using the strategy proposed in this study. The observation systems are the same as those in Fig. 2. (a) Synthetic velocity model with three anomalies: V1 (2200 m s−1), V2 (2200 m s−1) and V3 (1800 m s−1). The background velocity is 2000 m s−1. (b) Synthetic resistivity model with three anomalies: R1 (1000 Ω m), R2 (10 Ω m) and R3 (10 Ω m). The background resistivity is 100 Ω m. (c) Separately inverted velocity model. (d) Separately inverted resistivity model. (e) Jointly inverted velocity model. (f) Jointly inverted resistivity model. Figure 9. View largeDownload slide Comparison of synthetic and inverted velocity and resistivity models from the separate and joint inversions using the strategy proposed in this study. The observation systems are the same as those in Fig. 2. (a) Synthetic velocity model with three anomalies: V1 (2200 m s−1), V2 (2200 m s−1) and V3 (1800 m s−1). The background velocity is 2000 m s−1. (b) Synthetic resistivity model with three anomalies: R1 (1000 Ω m), R2 (10 Ω m) and R3 (10 Ω m). The background resistivity is 100 Ω m. (c) Separately inverted velocity model. (d) Separately inverted resistivity model. (e) Jointly inverted velocity model. (f) Jointly inverted resistivity model. Figure 10. View largeDownload slide Comparison of velocity and resistivity differences between inverted and true models for separate and joint inversions with our proposed strategy. (a) Velocity differences for separate inversion. (b) Velocity differences for joint inversion. (c) Resistivity differences for separate inversion. (d) Resistivity differences for joint inversion. Figure 10. View largeDownload slide Comparison of velocity and resistivity differences between inverted and true models for separate and joint inversions with our proposed strategy. (a) Velocity differences for separate inversion. (b) Velocity differences for joint inversion. (c) Resistivity differences for separate inversion. (d) Resistivity differences for joint inversion. 4 CONCLUSIONS Based on the joint inversion strategy of Um et al. (2014), we have proposed a new strategy to efficiently realize the cross-gradient joint inversion. In this strategy, different data sets are inverted separately and thus the existing inversion packages can be used in the joint inversion with trivial modifications. The model updates from separate inversions are then combined with the cross-gradient minimization system to make the two models similar in structure. Compared to the strategy of Um et al. (2014), our proposed strategy can better balance structural similarity and data fitting (via model updates from the separate inversions). Compared to the complete joint inversion system used by Bennington et al. (2015), the new joint inversion strategy avoids coupling between different systems and, can therefore, better account for varying data sensitivity amplitudes and different data coverage for different inversion systems. At the same time, it can produce models that are constrained by both geophysical data and cross-gradient minimization, efficiently realizing the cross-gradient structure constraint. We tested our proposed strategy on 2-D cross borehole seismic traveltime and DC resistivity tomography systems. Compared to the separate inversions, our strategy can better recover the resistivity anomalies while suppressing artefacts surrounding velocity anomalies. The synthetic test indicates that our proposed joint inversion strategy can be easily used to realize the cross-gradient joint inversion constraint while at the same time providing a better balance between data fitting and structural similarity than the strategy of Um et al. (2014). Acknowledgements We are grateful for the availability of the MSFM package (Hassouna & Farag 2007) for seismic traveltime calculation and the RESINVM3D package (Pidlisecky et al.2007) for DC resistivity tomography. We acknowledge the constructive comments from two anonymous reviewers that are very helpful for improving the manuscript. This work is supported by the National Natural Science Foundation of China (41704039) and by the Open Projects of State Key Laboratory of Coal Resources and Safe Mining, CUMT (SKLCRSM16KF09). REFERENCES Abubakar A., Li M., Liu J., Habashy T.M., 2009. Simultaneous joint inversion of MT and CSEM data using a multiplicative cost function, in SEG Technical Program Expanded Abstracts 2009, pp. 719– 723, Society of Exploration Geophysicists. Albouy Y., Andrieux P., Rakotondrasoa G., Ritz M., Descloitres M., Join J.L., Rasolomanana E., 2001. Mapping coastal aquifers by joint inversion of DC and TEM soundings—three case histories, Ground Water , 39( 1), 87– 97. https://doi.org/10.1111/j.1745-6584.2001.tb00354.x Google Scholar CrossRef Search ADS   Aster R.C., Borchers B., Thurber C.H. 2013. Parameter Estimation and Inverse Problems , 2nd edn, Elsevier Inc. Bennington N.L., Zhang H., Thurber C.H., Bedrosian P.A., 2015. Joint inversion of seismic and magnetotelluric data in the Parkfield Region of California using the normalized cross-gradient constraint, Pure appl. Geophys. , 172( 5), 1033– 1052. https://doi.org/10.1007/s00024-014-1002-9 Google Scholar CrossRef Search ADS   Bortolozo C.A., Couto M.A., Porsani J.L., Almeida E.R., Santos F.A.M., 2014. Geoelectrical characterization using joint inversion of VES/TEM data: a case study in Paraná Sedimentary Basin, São Paulo State, Brazil, J. Appl. Geophys. , 111: 33– 46. https://doi.org/10.1016/j.jappgeo.2014.09.009 Google Scholar CrossRef Search ADS   Bosch M., 1999. Lithologic tomography: from plural geophysical data to lithology estimation, J. geophys. Res. , 104( B1), 749– 766. https://doi.org/10.1029/1998JB900014 Google Scholar CrossRef Search ADS   Calvetti D., Morigi S., Reichel L., Sgallari F., 2000. Tikhonov regularization and the L-curve for large discrete ill-posed problems, J. Comput. Appl. Math. , 123( 1–2), 423– 446. https://doi.org/10.1016/S0377-0427(00)00414-3 Google Scholar CrossRef Search ADS   Fang H., Zhang H., Yao H., Allam A., Zigone D., Ben-Zion Y., Thurber C., van der Hilst R.D., 2016. A new algorithm for three-dimensional joint inversion of body wave and surface wave data and its application to the Southern California plate boundary region, J. geophys. Res. , 121( 5), 3557– 3569. Google Scholar CrossRef Search ADS   Fregoso E., Gallardo L.A., 2009. Cross-gradients joint 3D inversion with applications to gravity and magnetic data, Geophysics , 74( 4), L31– L42. https://doi.org/10.1190/1.3119263 Google Scholar CrossRef Search ADS   Gallardo L.A., Meju M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data, Geophys. Res. Lett. , 30( 13), doi:10.1029/2003GL017370. Gallardo L.A., Meju M.A., 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints, J. geophys. Res. , 109( B3), doi:10.1029/2003JB002716. Gallardo L.A., Meju M.A., 2011. Structure-coupled multiphysics imaging in geophysical sciences, Rev. Geophys. , 49( 1), doi: 10.1029/2010RG000330 Günther T., Rücker C., 2006. A new joint inversion approach applied to the combined tomography of DC resistivity and seismic refraction data, in Symposium on the Application of Geophysics to Engineering and Environmental Problems 2006, pp. 1196– 1202, doi.org/10.4133/1.2923578 Günther T., Dlugosch R., Holland R., Yaramanci U., 2010, Aquifer characterization using coupled inversion of DC/IP and MRS data on a hydrogeophysical test-site, in SAGEEP, Symposium on the Application of Geophysics to Engineering and Environmental Problems 2010 , pp. 302– 307. Haber E., Oldenburg D., 1997. Joint inversion: a structural approach, Inverse Probl. , 13( 1), 63– 77. https://doi.org/10.1088/0266-5611/13/1/006 Google Scholar CrossRef Search ADS   Hassouna M.S., Farag A.A., 2007. Multi-stencils fast marching methods: a highly accurate solution to the eikonal equation on cartesian domains, in IIEEE Trans. Pattern Anal. Mach. Intell. , 29( 9), 1563– 1574. Google Scholar CrossRef Search ADS   Jupp D.L.B., Vozoff K., 1975. Stable iterative methods for the inversion of geophysical data, Geophys. J. Int. , 42( 3), 957– 976. https://doi.org/10.1111/j.1365-246X.1975.tb06461.x Google Scholar CrossRef Search ADS   Lines L.R., Schultz A.K., Treitel S., 1988. Cooperative inversion of geophysical data, Geophysics , 53( 1), 8– 20. https://doi.org/10.1190/1.1442403 Google Scholar CrossRef Search ADS   Lu D.B., Wang F., Chen X.D., Qu J., Wang H., 2017. An improved ERT approach for the investigation of subsurface structures, Pure appl. Geophys. , 174( 1), 375– 386. https://doi.org/10.1007/s00024-016-1386-9 Google Scholar CrossRef Search ADS   Maceira M., Ammon C.J., 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure, J. geophys. Res. , 114( B2), doi:10.1029/2007JB005157. https://doi.org/10.1029/2007JB005157 McGillivray P.R., Oldenburg D.W., 1990. Methods for calculating Frechet derivatives and sensitivities for the non-linear inverse problem: a comparative study, Geophys. Prospect. , 38( 5), 499– 524. https://doi.org/10.1111/j.1365-2478.1990.tb01859.x Google Scholar CrossRef Search ADS   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, Geophys. J. Int. , 184( 1), 477– 493. https://doi.org/10.1111/j.1365-246X.2010.04856.x Google Scholar CrossRef Search ADS   Ogunbo J.N., Zhang J., 2014. Joint seismic traveltime and TEM inversion for near surface imaging, in SEG Technical Program Expanded Abstracts 2014, pp. 2104– 2108, Society of Exploration Geophysicists. Peng M., Tan H.D., Jiang M., Wang W., Li Q.-Q., Zhang L.-S., 2012. Joint inversion of receiver functions and magnetotelluric data: application to crustal and mantle structure beneath central Namche Barwa, eastern Himalayan syntaxis, Chinese J. Geophys (in Chinese) , 55( 7), 2281– 2291. Pidlisecky A., Haber E., Knight R., 2007. RESINVM3D: a 3D resistivity inversion package, Geophysics , 72( 2), H1– H10. https://doi.org/10.1190/1.2402499 Google Scholar CrossRef Search ADS   Rao Y., Wang Y., Chen S., Wang J., 2016. Crosshole seismic tomography with cross-firing geometry, Geophysics , 81( 4), R139– R146. https://doi.org/10.1190/geo2015-0677.1 Google Scholar CrossRef Search ADS   Rawlinson N., Sambridge M., 2005. The fast marching method: an effective tool for tomographic imaging and tracking multiple phases in complex layered media, Explor. Geophys. , 36( 4), 341– 350. https://doi.org/10.1071/EG05341 Google Scholar CrossRef Search ADS   Ronczka M., Hellman K., Günther T., Wisén R., Dahlin T. ( 2017), Electric resistivity and seismic refraction tomography: a challenging joint underwater survey at Äspö Hard Rock Laboratory, Solid Earth , 8( 3), 671– 682. https://doi.org/10.5194/se-8-671-2017 Google Scholar CrossRef Search ADS   Sambridge M.S., 1990. Non-linear arrival time inversion: constraining velocity anomalies by seeking smooth models in 3-D, Geophys. J. Int. , 102( 3), 653– 677. https://doi.org/10.1111/j.1365-246X.1990.tb04588.x Google Scholar CrossRef Search ADS   Scholl C., Neumann J., Watts M.D., Hallinan S., Mule S., 2016. Geologically constrained 2D and 3D airborne EM inversion through cross-gradient regularization and multi-grid efficiency, in ASEG Extended Abstracts 2016 , pp. 54– 67. https://doi.org/10.1071/ASEG2016ab229 Shi Z., Hobbs R.W., Moorkamp M., Tian G., Jiang L., 2017. 3-D cross-gradient joint inversion of seismic refraction and DC resistivity data, J. Appl. Geophys. , doi:10.1016/j.jappgeo.2017.04.008. https://doi.org/10.1016/j.jappgeo.2017.04.008 Tikhonov A.N., Arsenin V.I.A.. 1977. Solutions of Ill-posed Problems , Vh Winston. Um E.S., Commer M., Newman G.A. 2014. A strategy for coupled 3D imaging of large-scale seismic and electromagnetic data sets: application to subsalt imaging, Geophysics , 79( 3), ID1– ID13. https://doi.org/10.1190/geo2013-0053.1 Google Scholar CrossRef Search ADS   Wu X.-P., Xu G.-M., 2000. Study on 3D resistivity inversion using conjugate gradient method, China J. Geophys . ( in Chinese), 43( 3), 450– 458. Google Scholar CrossRef Search ADS   Zhang H., Maceira M., Roux P., Thurber C., 2014. Joint inversion of body-wave arrival times and surface-wave dispersion for three-dimensional seismic structure around SAFOD, Pure appl. Geophys. , 171( 11), 3013– 3022. https://doi.org/10.1007/s00024-014-0806-y Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

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

### DeepDyve is your personal research library

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

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### DeepDyve Freelancer ### DeepDyve Pro Price FREE$49/month

\$360/year
Save searches from