Characterisation of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data

Introduction

• The study focuses on establishing the relationship between electrical resistivity and seismic compressional (P) wave velocity in heterogeneous near-surface materials.
• This relationship is crucial for improved petrophysical characterization, which is necessary for understanding flow and transport processes in shallow-depth environments.
• Correlations between resistivity and velocity are observed at different spatial scales in collocated experiments.
• Conventional separate inversions of DC resistivity and seismic refraction data can lead to conflicting models.
• Joint inversion is a better approach, allowing objective testing of resistivity-velocity interrelationships derived from separate interpretations.
• Joint multidimensional resistivity-velocity inversion is challenging due to the lack of an established analytical relationship between resistivity and velocity.

Two Philosophical Approaches to Joint 2D Inversion

• Petrophysical Approach:
  • Assumes resistivity and velocity are functions of porosity and water saturation.
  • The validity of existing correlation models in general geological media is debated.
  • There is no single relationship for accurately predicting variations in clay content, shape, or interconnections of pores on geophysical measurements in field situations.
• Structural Approach:
  • Assumes both methods sense the same underlying geology, which structurally controls petrophysical property distribution.
  • This study favors the structural approach, deriving petrophysical information from resultant models.
• Zhang and Dale Morgan's study:
  • The only reported account of two-dimensional joint structural inversion of non-invasive dc resistivity and seismic sounding data in the literature.
  • Assume that the Laplacian of a DC resistivity image should resemble that of the seismic image.
  • Required a scaling factor to weight the relevance of each part (i.e., DC or seismic Laplacians).
  • The reported synthetic example showed a direct correlation between the resistivities and seismic velocities, which may not always be the case in field situations.

Questions Addressed

• How to recover a reliable joint inversion model over porous and non-porous heterogeneous natural materials where predictive models may not always hold.
• Is it possible to develop a generalized quantitative criterion for evaluating resistivity-velocity images for congruency or similarity in complex environments with heterogeneous materials?
• Can joint 2D inversion offer a means for improved characterization of heterogeneous materials?

Methodology for Joint 2D Resistivity-Seismic Inversion

Problem Formulation: Cross-Gradients Function

• A joint problem formulation is adopted with a novel resistivity-velocity cross-gradients function to link the DC resistivity model and the seismic velocity model.
• The cross-product of the gradients is defined as: $!t(x, y, z) = \nabla m<em>r(x, y, z) \times \nabla m</em>s(x, y, z)$
  • Where $m<em>r$ refers to the logarithm of resistivity, and $m</em>s$ refers to the P-wave slowness.
• In the two-dimensional case, $!t(x, y, z)$ is treated as a scalar $t$.
• Objective function to minimize: $\underset{m<em>r, m</em>s}{Minimize} \space ||d - f(m)||<em>{C</em>{dd}^{-1}}^2 + a<em>r ||D m</em>r||^2 + a<em>s ||D m</em>s||^2 + ||m<em>r - m</em>{Rr}||<em>{C</em>{RR}^{-1}}^2 + ||m<em>s - m</em>{Rs}||<em>{C</em>{RR}^{-1}}^2$
  • Subject to $t(m) = 0$
• Where:
  • $d$ is the vector of observed data (logarithm of apparent resistivity, $d<em>r$ and seismic traveltimes, $d</em>s$).
  • $m = [m<em>r: m</em>s]^T$ is the vector of the model parameters.
  • $f$ is the theoretical model response.
  • $D$ is the discrete version of a smoothing operator acting on $m$.
  • $m<em>R = [m</em>{Rr}: m_{Rs}]^T$ is an a priori model.
  • $C_{dd}$ is the covariance of the field data (assumed diagonal, i.e., fully uncorrelated data).
  • $C_{RR}$ is the covariance of the a priori model (assumed diagonal).
  • $a<em>r$ and $a</em>s$ are weighting factors to control the level of smoothing of the resistivity and seismic models.
• The cross-gradients criterion requires that spatial changes in both resistivity and velocity point in the same or opposite direction, regardless of amplitude.
• If a boundary exists, both methods must sense it in a common orientation, regardless of the physical property changes' amplitude.
• The technique allows for geological boundaries with significant changes in only electrical resistivity or seismic velocity.

Description of the Inversion Scheme

• The subsurface model is discretized into rectangular cells of variable sizes, optimized according to the sensitivity of resistivity and seismic measurements.
• The DC forward problem is solved using the fast approximate scheme of Pe´rez-Flores et al. [2001].
• The seismic forward problem uses the technique of progressive finite differences of Vidale [1990], incorporating features from Zelt and Barton [1998] to compute accurate travel times efficiently.
• Discrete version of equation (1) simplification:
  $t \approx \frac{1}{4 \Delta x \Delta z} ((m<em>{rc} + m</em>{rb} - m<em>{rr} - m</em>{rc})(m<em>{sc}- m</em>{sb}) - (m<em>{rb} - m</em>{rr})(m<em>{sr} - m</em>{sc}))$
• Using a first-order Taylor series expansion, equation (2) is equivalent to:
  $min \frac{1}{2} m^T N<em>1 m + n</em>2^T m$
  $subject \space to \space t(m<em>0) + B(m - m</em>0) = 0$
• Where:
  $N<em>1 = \begin{bmatrix} A</em>r^T C<em>{dd}^{-1} A</em>r + a<em>r^2 D^T D + C</em>{RRr}^{-1} &amp; 0 \ 0 &amp; A<em>s^T C</em>{dd}^{-1} A<em>s + a</em>s^2 D^T D + C<em>{RRs}^{-1} \end{bmatrix}$$n</em>2 = \begin{bmatrix} A<em>r^T C</em>{dd}^{-1} (d<em>r - f</em>r(m<em>{0r})) + C</em>{RRr}^{-1} (m<em>{0r} - m</em>{Rr}) \ A<em>s^T C</em>{dd}^{-1} (d<em>s - f</em>s(m<em>{0s})) + C</em>{RRs}^{-1} (m<em>{0s} - m</em>{Rs}) \end{bmatrix}$
• $A<em>r$, $A</em>s$, and $B$ are the respective partial derivatives of $f<em>r$, $f</em>s$, and $t$ evaluated at the initial model, $m<em>0 = [m</em>{0r}: m_{0s}]^T$.
• The Jacobian matrix for seismics, $A_s$, is computed using ray tracing as suggested by Vidale [1990] and Zelt and Barton [1998].
• The Jacobian matrix for DC resistivity, $A_r$, is computed using the code of Pe´rez-Flores et al. [2001].
• Solution to (4) used in the iterative scheme:
  $m = N<em>1^{-1} n</em>2 - N<em>1^{-1} B^T (B N</em>1^{-1} B^T)^{-1} (B N<em>1^{-1} n</em>2 - B m<em>0 + t(m</em>0) )$
• Weighting factors are initially assigned large values and gradually reduced in subsequent iterations until the data are fitted to the required level.
• The joint inversion is initiated using a half-space model in the absence of reliable a priori information.
• The cross-gradients criterion and regularization measures ensure that the resolution characteristics of the individual data sets are fully exploited in the search for structurally linked models.
• With the cross-gradients criterion, there is no need to define or assume any interdependence of resistivity and seismic velocity ab initio, avoiding bias in the inverse solution.

Near-Surface Characterization by Joint Data Inversion

• Field data sets from collocated DC resistivity and seismic refraction experiments were used and have been separately inverted, and the resulting models correlated to define a resistivity-velocity relation [see Meju et al., 2003]. The same dataset is used for simultaneous inversion using the presented algorithm.
• The data consists of first arrivals from 4 shot points recorded at 2m intervals over a 166m spread.
• The DC resistivity data consists of six in-line Schlumberger soundings with half-current electrode spacings (AB/2) ranging from 1.5 to 90 m.
• The test area is a buried hillside composed of fractured granodiorite, overlain by a heterogeneous mudstone sequence and glacial drift deposits.
• The optimal seismic velocity and resistivity images are derived by joint inversion using initial half-space models.
• Structural similarities and good spatial correlation of velocity and resistivity are observed in the reconstructed distributions of model parameters.
• The adopted cross-gradients criterion serves for geologic structural control but does not force the two models into conformity if not justified by the field data.
• Observations:
  • A zone of coincident low velocities and high resistivities in the top 3 m where there are glacial drift deposits.
  • An intermediate zone (3 – 20 m depth) showing a heterogeneous character and undulatory base, possibly consisting of subcropping basement blocks and minor troughs filled by mudstone or clayey materials.
  • A basal section of high velocity and high resistivity, possibly corresponding to hard granodiorite bedrock.
  • The distribution of possible rock-contacts is concordant in both models as a consequence of the cross-gradients structural constraint.
• The resistivity-velocity relationship for the reconstructed models suggests distinct sub-groupings within these near-surface materials.
• The trend reversal yielding the v-shaped curves may be a consequence of water saturation (water table) or a natural divide between unconsolidated and consolidated materials [cf. Meju et al., 2003].
• Seismic rays are critically refracted at the top of the structural units mapped as zones 3 and 4, defining the possible boundary between consolidated and unconsolidated materials.

Conclusion

• The incorporation of the cross-gradients criterion in 2D inversion improves the structural conformity between velocity and resistivity images without forcing or assuming the form of the relationship between them, leading to geologically meaningful solutions.
• The cross-gradients criterion also allows the delineation of those subsurface features for which only one of the geophysical techniques is sensitive, leading to a better structural characterization.
• The application of joint 2D inversion with cross-gradients to field data from collocated seismic and DC resistivity experiments has led to an improved characterization of near-surface heterogeneous materials.
• Unconsolidated and consolidated (or possibly unsaturated and saturated) materials may be sub-classified on the basis of their resistivity-velocity relationship evinced from joint inversion, a feature that was not observed in a previous study using separate inversion models.
• The cross-gradients approach adopted here can also be used for 3D problems and for any combination of independent geophysical methods.