Showing papers in "Geophysical Prospecting in 2004"

A numerical comparison of 2D resistivity imaging with 10 electrode arrays

Abstract: Numerical simulations are used to compare the resolution and efficiency of 2D resistivity imaging surveys for 10 electrode arrays. The arrays analysed include polepole (PP), pole-dipole (PD), half-Wenner (HW), Wenner-α (WN), Schlumberger (SC), dipole-dipole (DD), Wenner-β (WB), γ -array (GM), multiple or moving gradient array (GD) and midpoint-potential-referred measurement (MPR) arrays. Five synthetic geological models, simulating a buried channel, a narrow conductive dike, a narrow resistive dike, dipping blocks and covered waste ponds, were used to examine the surveying efficiency (anomaly effects, signal-to-noise ratios) and the imaging capabilities of these arrays. The responses to variations in the data density and noise sensitivities of these electrode configurations were also investigated using robust (L1-norm) inversion and smoothness-constrained least-squares (L2-norm) inversion for the five synthetic models. The results show the following. (i) GM and WN are less contaminated by noise than the other electrode arrays. (ii) The relative anomaly effects for the different arrays vary with the geological models. However, the relatively high anomaly effects of PP, GM and WB surveys do not always give a high-resolution image. PD, DD and GD can yield better resolution images than GM, PP, WN and WB, although they are more susceptible to noise contamination. SC is also a strong candidate but is expected to give more edge effects. (iii) The imaging quality of these arrays is relatively robust with respect to reductions in the data density of a multi-electrode layout within the tested ranges. (iv) The robust inversion generally gives better imaging results than the L2-norm inversion, especially with noisy data, except for the dipping block structure presented here. (v) GD and MPR are well suited to multichannel surveying and GD may produce images that are comparable to those obtained with DD and PD. Accordingly, the GD, PD, DD and SC arrays are strongly recommended for 2D resistivity imaging, where the final choice will be determined by the expected geology, the purpose of the survey and logistical considerations.

Wave-equation migration velocity analysis. I. Theory

Abstract: We present a migration velocity analysis (MVA) method based on wavefield extrapolation. Similarly to conventional MVA, our method aims at iteratively improving the quality of the migrated image, as measured by the flatness of angle-domain common-image gathers (ADCIGs) over the aperture-angle axis. However, instead of inverting the depth errors measured in ADCIGs using ray-based tomography, we invert ‘image perturbations’ using a linearized wave-equation operator. This operator relates perturbations of the migrated image to perturbations of the migration velocity. We use prestack Stolt residual migration to define the image perturbations that maximize the focusing and flatness of ADCIGs. Our linearized operator relates slowness perturbations to image perturbations, based on a truncation of the Born scattering series to the first-order term. To avoid divergence of the inversion procedure when the velocity perturbations are too large for Born linearization of the wave equation, we do not invert directly the image perturbations obtained by residual migration, but a linearized version of the image perturbations. The linearized image perturbations are computed by a linearized prestack residual migration operator applied to the background image. We use numerical examples to illustrate how the backprojection of the linearized image perturbations, i.e. the gradient of our objective function, is well behaved, even in cases when backprojection of the original image perturbations would mislead the inversion and take it in the wrong direction. We demonstrate with simple synthetic examples that our method converges even when the initial velocity model is far from correct. In a companion paper, we illustrate the full potential of our method for estimating velocity anomalies under complex salt bodies.

Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study

Abstract: An integrated multiscale seismic imaging flow is applied to dense onshore wideaperture seismic data recorded in a complex geological setting (thrust belt). An initial P-wave velocity macromodel is first developed by first-arrival traveltime tomography. This model is used as an initial guess for subsequent full-waveform tomography, which leads to greatly improved spatial resolution of the P-wave velocity model. However, the application of full-waveform tomography to the high-frequency part of the source bandwidth is difficult, due to the non-linearity of this kind of method. Moreover, it is computationally expensive at high frequencies since a finitedifference method is used to model the wave propagation. Hence, full-waveform tomography was complemented by asymptotic prestack depth migration to process the full-source bandwidth and develop a sharp image of the short wavelengths. The final traveltime tomography model and two smoothed versions of the final full-waveform tomography model were used as a macromodel for the prestack depth migration. In this study, wide-aperture multifold seismic data are used. After specific preprocessing of the data, 16 frequency components ranging from 5.4 Hz to 20 Hz were inverted in cascade by the full-waveform tomography algorithm. The full-waveform tomography successfully imaged SW-dipping structures previously identified as highresistivity bodies. The relevance of the full-waveform tomography models is demonstrated locally by comparison with a coincident vertical seismic profiling (VSP) log available on the profile. The prestack depth-migrated images, inferred from the traveltime, and the smoothed full-waveform tomography macromodels are shown to be, on the whole, consistent with the final full-waveform tomography model. A more detailed analysis, based on common-image gather computations, and local comparison with the VSP log revealed that the most accurate migrated sections are those obtained from the full-waveform tomography macromodels. A resolution analysis suggests that the asymptotic prestack depth migration successfully migrated the wideaperture components of the data, allowing medium wavelengths in addition to the short wavelengths of the structure to be imaged.

Adaptive subtraction of multiples using the l1-norm

Abstract: A strategy for multiple removal consists of estimating a model of the multiples and then adaptively subtracting this model from the data by estimating shaping filters. A possible and efficient way of computing these filters is by minimizing the difference or misfit between the input data and the filtered multiples in a least-squares sense. Therefore, the signal is assumed to have minimum energy and to be orthogonal to the noise. Some problems arise when these conditions are not met. For instance, for strong primaries with weak multiples, we might fit the multiple model to the signal (primaries) and not to the noise (multiples). Consequently, when the signal does not exhibit minimum energy, we propose using the L 1 -norm, as opposed to the L 2 -norm, for the filter estimation step. This choice comes from the well-known fact that the L 1 -norm is robust to ‘large’ amplitude differences when measuring data misfit. The L 1 -norm is approximated by a hybrid L 1 /L 2 -norm minimized with an iteratively reweighted least-squares (IRLS) method. The hybrid norm is obtained by applying a simple weight to the data residual. This technique is an excellent approximation to the L 1 -norm. We illustrate our method with synthetic and field data where internal multiples are attenuated. We show that the L 1 -norm leads to much improved attenuation of the multiples when the minimum energy assumption is violated. In particular, the multiple model is fitted to the multiples in the data only, while preserving the primaries.

On anelliptic approximations for qP velocities in VTI media

Abstract: A unified approach to approximating phase and group velocities of qP seismic waves in a transversely isotropic medium with vertical axis of symmetry (VTI) is developed. While the exact phase‐velocity expressions involve four independent parameters to characterize the elastic medium, the proposed approximate expressions use only three parameters. This makes them more convenient for use in surface seismic experiments, where the estimation of all four parameters is problematic. The three‐parameter phase‐velocity approximation coincides with the previously published ‘acoustic’ approximation of Alkhalifah. The group‐velocity approximation is new and noticeably more accurate than some of the previously published approximations. An application of the group‐velocity approximation for finite‐difference computation of traveltimes is shown.

Euler deconvolution of the analytic signal and its application to magnetic interpretation

Abstract: Euler deconvolution and the analytic signal are both used for semi-automatic interpretation of magnetic data. They are used mostly to delineate contacts and obtain rapid source depth estimates. For Euler deconvolution, the quality of the depth estimation depends mainly on the choice of the proper structural index, which is a function of the geometry of the causative bodies. Euler deconvolution applies only to functions that are homogeneous. This is the case for the magnetic field due to contacts, thin dikes and poles. Fortunately, many complex geological structures can be approximated by these simple geometries. In practice, the Euler equation is also solved for a background regional field. For the analytic signal, the model used is generally a contact, although other models, such as a thin dike, can be considered. It can be shown that if a function is homogeneous, its analytic signal is also homogeneous. Deconvolution of the analytic signal is then equivalent to Euler deconvolution of the magnetic field with a background field. However, computation of the analytic signal effectively removes the background field from the data. Consequently, it is possible to solve for both the source location and structural index. Once these parameters are determined, the local dip and the susceptibility contrast can be determined from relationships between the analytic signal and the orthogonal gradients of the magnetic field. The major advantage of this technique is that it allows the automatic identification of the type of source. Implementation of this approach is demonstrated for recent high-resolution survey data from an Archean granite-greenstone terrane in northern Ontario, Canada.

Effective elastic properties of randomly fractured soils: 3D numerical experiments

Abstract: This paper is concerned with numerical tests of several rock physical relationships. The focus is on effective velocities and scattering attenuation in 3D fractured media. We apply the so-called rotated staggered finite-difference grid (RSG) technique for numerical experiments. Using this modified grid, it is possible to simulate the propagation of elastic waves in a 3D medium containing cracks, pores or free surfaces without applying explicit boundary conditions and without averaging the elastic moduli. We simulate the propagation of plane waves through a set of randomly cracked 3D media. In these numerical experiments we vary the number and the distribution of cracks. The synthetic results are compared with several (most popular) theories predicting the effective elastic properties of fractured materials. We find that, for randomly distributed and randomly orientated non-intersecting thin penny-shaped dry cracks, the numerical simulations of P- and S-wave velocities are in good agreement with the predictions of the self-consistent approximation. We observe similar results for fluid-filled cracks. The standard Gassmann equation cannot be applied to our 3D fractured media, although we have very low porosity in our models. This is explained by the absence of a connected porosity. There is only a slight difference in effective velocities between the cases of intersecting and non-intersecting cracks. This can be clearly demonstrated up to a crack density that is close to the connectivity percolation threshold. For crack densities beyond this threshold, we observe that the differential effective-medium (DEM) theory gives the best fit with numerical results for intersecting cracks. Additionally, it is shown that the scattering attenuation coefficient (of the mean field) predicted by the classical Hudson approach is in excellent agreement with our numerical results.

Wave‐equation migration velocity analysis. II. Subsalt imaging examples

Abstract: Subsalt imaging is strongly dependent on the quality of the velocity model. However, rugose salt bodies complicate wavefield propagation and lead to subsalt multipathing, illumination gaps and shadow zones, which cannot be handled correctly by conventional traveltime-based migration velocity analysis (MVA). We overcome these limitations by the wave-equation MVA technique, introduced in a companion paper, and demonstrate the methodology on a realistic synthetic data set simulating a salt-dome environment and a Gulf of Mexico data set. We model subsalt propagation using wave paths created by one-way wavefield extrapolation. Those wave paths are much more accurate and robust than broadband rays, since they inherit the frequency dependence and multipathing of the underlying wavefield. We formulate an objective function for optimization in the image space by relating an image perturbation to a perturbation of the velocity model. The image perturbations are defined using linearized prestack residual migration, thus ensuring stability, relative to the first-order Born approximation assumptions. Synthetic and real data examples demonstrate that wave-equation MVA is an effective tool for subsalt velocity analysis, even when shadows and illumination gaps are present.

Effect of grain scale alignment on seismic anisotropy and reflectivity of shales

Abstract: The elastic properties and anisotropy of shales are strongly influenced by the degree of alignment of the grain scale texture. In general, an orientation distribution function (ODF) can be used to describe this alignment, which, in practice, can be characterized by two Legendre coefficients. We discuss various statistical ODFs that define the alignment by spreading from a mean value; in particular, the Gaussian, Fisher and Bingham distributions. We compare the statistical models with an ODF resulting from pure vertical compaction (no shear strain) of a sediment. The compaction ODF may be used to estimate how the elastic properties and anisotropy evolve due to burial of clayey sediments. Our study shows that the three statistical ODFs produce almost identical correspondence between the two Legendre coefficients as a function of the spreading parameter, so that the spreading parameter of one ODF can be converted to the spreading parameter of another ODF. In most cases it is then sufficient to apply the spreading parameter for the ODF instead of the two Legendre coefficients. The effect of compaction on the ODF gives a slightly different correspondence between the two Legendre coefficients from that for the other models. In principle, this opens up the possibility of distinguishing anisotropy effects due to compaction from those due to other processes. We also study reflection amplitudes versus angle of incidence (AVA) for all wave modes, where shales having various ODFs overlie an isotropic medium. The AVA responses are modelled using both exact and approximation formulae, and their intercepts and gradients are compared. The modelling shows that the S-wave velocity is sensitive to any perturbation in the spreading parameter, while the P-wave velocity becomes increasingly sensitive to a perturbation of a less ordered system. Similar observations are found for the AVA of the P-P and P-SV waves. Modelling indicates that a combined use of the amplitude versus offset of P-P and P-SV reflected waves may reveal certain grain scale alignment properties of shale-like rocks.

Modelling of GPR waves for lossy media obeying a complex power law of frequency for dielectric permittivity

Abstract: The attenuation of ground-penetrating radar (GPR) energy in the subsurface decreases and shifts the amplitude spectrum of the radar pulse to lower frequencies (absorption) with increasing traveltime and causes also a distortion of wavelet phase (dispersion). The attenuation is often expressed by the quality factor Q . For GPR studies, Q can be estimated from the ratio of the real part to the imaginary part of the dielectric permittivity. We consider a complex power function of frequency for the dielectric permittivity, and show that this dielectric response corresponds to a frequency-independent- Q or simply a constant- Q model. The phase velocity (dispersion relationship) and the absorption coefficient of electromagnetic waves also obey a frequency power law. This approach is easy to use in the frequency domain and the wave propagation can be described by two parameters only, for example Q and the phase velocity at an arbitrary reference frequency. This simplicity makes it practical for any inversion technique. Furthermore, by using the Hilbert transform relating the velocity and the absorption coefficient (which obeys a frequency power law), we find the same dispersion relationship for the phase velocity. Both approaches are valid for a constant value of Q over a restricted frequency-bandwidth, and are applicable in a material that is assumed to have no instantaneous dielectric response. Many GPR profiles acquired in a dry aeolian environment have shown a strong reflectivity inside dunes. Changes in water content are believed to be the origin of this reflectivity. We model the radar reflections from the bottom of a dry aeolian dune using the 1D wavelet modelling method. We discuss the choice of the reference wavelet in this modelling approach. A trial-and-error match of modelled and observed data was performed to estimate the optimum set of parameters characterizing the materials composing the site. Additionally, by combining the complex refractive index method (CRIM) and/or Topp equations for the bulk permittivity (dielectric constant) of moist sandy soils with a frequency power law for the dielectric response, we introduce them into the expression for the reflection coefficient. Using this method, we can estimate the water content and explain its effect on the reflection coefficient and on wavelet modelling.

3D angle‐domain common‐image gathers for migration velocity analysis

Abstract: Angle-domain common-image gathers (ADCIGs) are an essential tool for migration velocity analysis (MVA). We present a method for computing ADCIGs in 3D from the results of wavefield-continuation migration. The proposed methodology can be applied before or after the imaging step in a migration procedure. When computed before imaging, 3D ADCIGs are functions of the offset ray parameters ; we derive the geometric relationship that links the offset ray parameters to the aperture angle γ and the reflection azimuth φ. When computed after imaging, 3D ADCIGs are directly produced as functions of γ and φ. The mapping of the offset ray parameters into the angles (γ, φ) depends on both the local dips and the local interval velocity; therefore, the transformation of ADCIGs computed before imaging into ADCIGs that are functions of the actual angles is difficult in complex structure. By contrast, the computation of ADCIGs after imaging is efficient and accurate even in the presence of complex structure and a heterogeneous velocity function. On the other hand, the estimation of the offset ray parameters is less sensitive to velocity errors than the estimation of the angles (γ, φ). When ADCIGs that are functions of the offset ray parameters are adequate for the application of interest (e.g. ray-based tomography), the computation of ADCIGs before imaging might be preferable. Errors in the migration velocity cause the image point in the angle domain to shift along the normal to the apparent geological dip. By assuming stationary rays (i.e. small velocity errors), we derive a quantitative relationship between this normal shift and the traveltime perturbation caused by velocity errors. This relationship can be directly used in an MVA procedure to invert depth errors measured from ADCIGs into migration velocity updates. In this paper, we use it to derive an approximate 3D residual moveout (RMO) function for measuring inconsistencies between the migrated images at different γ and φ. We tested the accuracy of our kinematic analysis on a 3D synthetic data set with steeply dipping reflectors and a vertically varying propagation velocity. The tests confirm the accuracy of our analysis and illustrate the limitations of the straight-ray approximation underlying our derivation of the 3D RMO function.

Stereotomography: a semi-automatic approach for velocity macromodel estimation

Abstract: Most methods for velocity macromodel estimation require considerable operator input, mainly concerning the regularization and the picking of events in the data set or in the migrated images. For both these aspects, slope tomography methods offer interesting solutions. They consider locally coherent events characterized by their slopes in the data cube. Picking is then much easier and consequently denser than in standard traveltime tomography. Stereotomography is the latest slope tomography method. In recent years it has been improved significantly, both from an algorithmic point of view and in terms of practical use. Robust and fast procedures are now available for 2D stereotomographic picking and optimization. Concerning the picking, we propose simple criteria for the selection of relevant data among the automatically picked events. This enables an accurate smooth velocity macromodel to be estimated quite rapidly and with very limited operator intervention. We demonstrate the method using a 2D line extracted from the Oseberg NH8906 data set.

3D tomographic velocity model estimation with kinematic wavefield attributes

Abstract: A tomographic inversion method is presented that uses kinematic information in the form of zero-offset traveltimes and kinematic wavefield attributes (first and second spatial traveltime derivatives) to determine smooth, laterally inhomogeneous 3D subsurface velocity models for depth imaging. The kinematic wavefield attributes can be extracted from the seismic prestack data by means of the common reflection surface (CRS) stack. The input for the tomography is then taken from the resulting attribute volumes at a number of pick locations in the CRS stacked zero-offset volume. As a smooth model description based on B-splines is used and reflection points are treated independently of each other, only locally coherent events in the stacked volume are required and very few picks are needed. Thus, picking is considerably simplified. During the iterative inversion process, the required forward-modelled quantities are obtained by dynamic ray tracing along normal rays pertaining to the input data points. Frechet derivatives for the tomographic matrix are calculated with ray perturbation theory. The inversion algorithm is demonstrated on a 3D synthetic data example, where the kinematic wavefield attributes have directly been obtained by forward modelling.

Comparison of methods for modelling the behaviour of bubbles produced by marine seismic airguns

Abstract: Three models for the dynamics of seismic airgun-generated bubbles and their associated far-field signals are developed and compared with geophysical data. The first model of an airgun-generated bubble uses a spherical approximation, the second is an approximate Lagrangian model which allows for small deformations from a spherical shape, whilst the final model is an axisymmetric boundary-integral method which permits the bubble to evolve into highly non-spherical geometries. The boundary-integral method also allows both geometric interference and strong dynamic interactions in multi-bubble studies. When comparing the spherical model to experimental data there are three apparent, significant differences: the magnitude of the primary pressure peak, which is greater in the model; the subsequent decay of the pressure peaks and motion – the experimental data demonstrating greater decay and a slower rise rate; and the frequency of oscillation, which is slower in the experimental data. It is believed that the first discrepancy is due to the initial stages of expansion where the compressed air is forced to sparge through the airgun ports. The other differences indicate that there is some other energy-loss mechanism which is not accounted for in the spherical bubble model. Non-spherical bubble behaviour is investigated through the use of two different deformable many-bubble codes and their predictions are compared with the spherical model and experimental data. The Lagrangian model predicts the formation of a buoyancy-driven liquid jet on the first collapse of a typical airgun bubble; however, the model breaks down when the bubble becomes significantly deformed, due to a low-order spherical-harmonic approximation for the potential. The axisymmetric boundary-integral code models the jet shape accurately and it is found that these bubbles evolve to toroidal geometries when the jet impacts on the opposite surface of the bubble. This highly non-spherical behaviour is readily observed on high-speed films of airgun bubbles, and is one key source of energy loss; it damps the pulsations of the bubble and slows its rise speed. Inter-bubble interactions are investigated using the two deformable bubble models, and the predictions are compared to field data. It was found that as the bubbles approach each other, their periods of oscillation increase in accordance with observations, and jets are formed in the direction of motion upon collapse.

Dispersion and radial depth of investigation of borehole modes

Abstract: Sonic techniques in geophysical prospecting involve elastic wave velocity measurements that are performed by placing acoustic transmitters and receivers in a fluid-filled borehole. The signals recorded at the receivers are processed to obtain compressional- and shear-wave velocities in the surrounding formation. These velocities are generally used in seismic surveys for the time-to-depth conversion and other formation parameters, such as porosity and lithology. Depending upon the type of transmitter used (e.g. monopole or dipole) and as a result of eccentering, it is possible to excite axisymmetric (n = 0), flexural (n = 1) and quadrupole (n = 2) families of modes propagating along the borehole. We present a study of various propagating and leaky modes that includes their dispersion and attenuation characteristics caused by radiation into the surrounding formation. A knowledge of propagation characteristics of borehole modes helps in a proper selection of transmitter bandwidth for suppressing unwanted modes that create problems in the inversion for the compressional- and shear-wave velocities from the dispersive arrivals. It also helps in the design of a transmitter for a preferential excitation of a given mode in order to reduce interference with drill-collar or drilling noise for sonic measurements-while-drilling. Computational results for the axisymmetric family of modes in a fast formation with a shear-wave velocity of 2032 m/s show the existence of Stoneley, pseudo-Rayleigh and anharmonic cut-off modes. In a slow formation with a shear-wave velocity of 508 m/s, we find the existence of the Stoneley mode and the first leaky compressional mode which cuts in at approximately the same normalized frequency ωa/V S = 2.5 (a is the borehole radius) as that of the fast formation. The corresponding modes among the flexural family include the lowest-order flexural and anharmonic cut-off modes. For both the fast and slow formations, the first anharmonic mode cuts in at a normalized frequency ωa/V S = 1.5 approximately. Cut-off frequencies of anharmonic modes are inversely proportional to the borehole radius in the absence of any tool. The borehole quadrupole mode can also be used for estimating formation shear slownesses. The radial depth of investigation with a quadrupole mode is marginally less than that of a flexural mode because of its higher frequency of excitation.

Airborne EM skin depths

Abstract: Skin depth is an electromagnetic (EM) scale length that provides a measure of the degree of attenuation experienced by a particular frequency of an EM system. As has been discussed in the literature, skin depth is not a complete measure of the depth of investigation, but the two may be related. Frequency-domain airborne EM systems employ pairs of transmitter and receiver coils that use a frequency range from several hundred hertz to over 100 kHz. For elevated dipoles, both geometrical and frequency-dependent attenuation of the induced fields must be considered. For airborne EM systems it is possible to define a skin depth based only on the electric field induced by the transmitter. A vertical decay scale length, here defined from the at-surface position of maximum electric field, enables the same skin-depth estimate to be obtained for both cases of vertical and horizontal dipolar excitation. Such dipolar skin depths associated with towed-bird and fixed-wing airborne systems are studied in relation to frequency, conductivity and sensor elevation. Dipolar skin depths are found to be much smaller than their plane-wave counterparts except at high frequency (>50 kHz) and in combination with high conductivity. For the majority of airborne systems the influence of altitude on skin depth is highly significant. Dipolar skin depths increase with increasing sensor elevation. Low frequencies display the greatest sensitivity. At low elevation (<40 m), geometrical attenuation dominates the behaviour of the skin depth. The study indicates that typical low-altitude airborne surveys provide vertically compact assessments of subsurface conductivity, well suited to near-surface, environmental applications.

Characterization of dipping fractures in a transverselyisotropic background

Abstract: Although it is believed that natural fracture sets predominantly have near-vertical orientation, oblique stresses and some other mechanisms may tilt fractures away from the vertical Here, we examine an effective medium produced by a single system of obliquely dipping rotationally invariant fractures embedded in a transversely isotropic with a vertical symmetry axis (VTI) background rock This model is monoclinic with a vertical symmetry plane that coincides with the dip plane of the fractures Multicomponent seismic data acquired over such a medium possess several distinct features that make it possible to estimate the fracture orientation For example, the vertically propagating fast shear wave (and the fast converted PS-wave) is typically polarized in the direction of the fracture strike The normal-moveout (NMO) ellipses of horizontal reflection events are co-orientated with the dip and strike directions of the fractures, which provides an independent estimate of the fracture azimuth However, the polarization vector of the slow shear wave at vertical incidence does not lie in the horizontal plane – an unusual phenomenon that can be used to evaluate fracture dip Also, for oblique fractures the shear-wave splitting coefficient at vertical incidence becomes dependent on fracture infill (saturation) A complete medium-characterization procedure includes estimating the fracture compliances and orientation (dip and azimuth), as well as the Thomsen parameters of the VTI background We demonstrate that both the fracture and background parameters can be obtained from multicomponent wide-azimuth data using the vertical velocities and NMO ellipses of PP-waves and two split SS-waves (or the traveltimes of PS-waves) reflected from horizontal interfaces Numerical tests corroborate the accuracy and stability of the inversion algorithm based on the exact expressions for the vertical and NMO velocities

Adjustment of regularization in ill-posed linear inverse problems by the empirical Bayes approach

Abstract: Regularization is the most popular technique to overcome the null space of model parameters in geophysical inverse problems, and is implemented by including a constraint term as well as the data‐misfit term in the objective function being minimized. The weighting of the constraint term relative to the data‐fitting term is controlled by a regularization parameter, and its adjustment to obtain the best model has received much attention. The empirical Bayes approach discussed in this paper determines the optimum value of the regularization parameter from a given data set. The regularization term can be regarded as representing a priori information about the model parameters. The empirical Bayes approach and its more practical variant, Akaike's Bayesian Information Criterion, adjust the regularization parameter automatically in response to the level of data noise and to the suitability of the assumed a priori model information for the given data. When the noise level is high, the regularization parameter is made large, which means that the a priori information is emphasized. If the assumed a priori information is not suitable for the given data, the regularization parameter is made small. Both these behaviours are desirable characteristics for the regularized solutions of practical inverse problems. Four simple examples are presented to illustrate these characteristics for an underdetermined problem, a problem adopting an improper prior constraint and a problem having an unknown data variance, all frequently encountered geophysical inverse problems. Numerical experiments using Akaike's Bayesian Information Criterion for synthetic data provide results consistent with these characteristics. In addition, concerning the selection of an appropriate type of a priori model information, a comparison between four types of difference‐operator model – the zeroth‐, first‐, second‐ and third‐order difference‐operator models – suggests that the automatic determination of the optimum regularization parameter becomes more difficult with increasing order of the difference operators. Accordingly, taking the effect of data noise into account, it is better to employ the lower‐order difference‐operator models for inversions of noisy data.

A comparison of cross‐hole electrical and seismic data in fractured rock

Abstract: Cross-hole anisotropic electrical and seismic tomograms of fractured metamorphic rock have been obtained at a test site where extensive hydrological data were available. A strong correlation between electrical resistivity anisotropy and seismic compressional-wave velocity anisotropy has been observed. Analysis of core samples from the site reveal that the shale-rich rocks have fabric-related average velocity anisotropy of between 10% and 30%. The cross-hole seismic data are consistent with these values, indicating that observed anisotropy might be principally due to the inherent rock fabric rather than to the aligned sets of open fractures. One region with velocity anisotropy greater than 30% has been modelled as aligned open fractures within an anisotropic rock matrix and this model is consistent with available fracture density and hydraulic transmissivity data from the boreholes and the cross-hole resistivity tomography data. However, in general the study highlights the uncertainties that can arise, due to the relative influence of rock fabric and fluid-filled fractures, when using geophysical techniques for hydrological investigations.

Determination of the regularization level of truncated singular-value decomposition inversion: The case of 1D inversion of MT data

Abstract: The total mean-square error (MSE) of the estimated model, defined as the sum of the standard model variance and the bias variance, is used to define the truncation level of the singular-value decomposition to give a reasonable balance between model resolution and model variance. This balance is determined largely by the data and no further assumptions are necessary except that the bias terms are estimated sufficiently well. This principle has been tested on the 1D magnetotelluric inverse problem with special emphasis on high-frequency radio magnetotelluric (RMT) data. Simulations clearly demonstrate that the method provides a good balance between resolution and variance. Starting from a homogeneous half-space, the best solution is sought for a fixed set of singular values. The model variance is estimated from the sum of the inverse eigenvalues squared, up to a certain threshold, and the bias variance is estimated from the model projections on the remaining eigenvectors. By varying the threshold, the minimum of the MSE is found for an increasing number of fixed singular values until the number of active singular values becomes greater than or equal to the estimated number. As a side-effect, the depth of penetration of a given set of measurements can be estimated very efficiently by simply noting at which depth the final model deviates little from the starting homogeneous half-space model. A suite of synthetic data is inverted and an example of inversion of one site is shown to illustrate how the truncation is carried out as the non-linear inversion process proceeds. A field example with a profile across a plume of contaminated groundwater in the Netherlands shows good agreement with the electrical resistivity obtained in a nearby borehole.

The effect of the electrical anisotropy on the response of helicopter‐borne frequency‐domain electromagnetic systems

Abstract: Helicopter electromagnetic (HEM) systems are commonly used for conductivity mapping and the data are often interpreted using an isotropic horizontally layered earth model. However, in regions with distinct dipping stratification, it is useful to extend the model to a layered earth with general anisotropy by assigning each layer a symmetrical 3 × 3 resistivity tensor. The electromagnetic (EM) field is represented by two scalar potentials, which describe the poloidal and toroidal parts of the magnetic field. Via a 2D Fourier transform, we obtain two coupled ordinary differential equations in the vertical coordinate. To stabilize the numerical calculation, the wavenumber domain is divided into two parts associated with small and large wavenumbers. The EM field for small wavenumbers is continued from layer to layer with the continuity conditions. For large wavenumbers, the EM field behaves like a DC field and therefore cannot be sensed by airborne EM systems. Thus, the contribution from the large wavenumbers is simply ignored. The magnetic fields are calculated for the vertical coaxial (VCX), horizontal coplanar (HCP) and vertical coplanar (VCP) coil configurations for a helicopter EM system. The apparent resistivities defined from the VCX, VCP and HCP coil responses, when plotted in polar coordinates, clearly identify the principal anisotropic axes of an anisotropic earth. The field example from the Edwards Aquifer recharge area in Texas confirms that the polar plots of the apparent resistivities identify the principal anisotropic axes that coincide well with the direction of the underground structures.

Efficient and accurate computation of seismic traveltimes and amplitudes

Abstract: We describe two practicable approaches for an efficient computation of seismic traveltimes and amplitudes The first approach is based on a combined finite-difference solution of the eikonal equation and the transport equation (the ‘FD approach’) These equations are formulated as hyperbolic conservation laws; the eikonal equation is solved numerically by a third-order ENO–Godunov scheme for the traveltimes whereas the transport equation is solved by a first-order upwind scheme for the amplitudes The schemes are implemented in 2D using polar coordinates The results are first-arrival traveltimes and the corresponding amplitudes The second approach uses ray tracing (the ‘ray approach’) and employs a wavefront construction (WFC) method to calculate the traveltimes Geometrical spreading factors are then computed from these traveltimes via the ray propagator without the need for dynamic ray tracing or numerical differentiation With this procedure it is also possible to obtain multivalued traveltimes and the corresponding geometrical spreading factors Both methods are compared using the Marmousi model The results show that the FD eikonal traveltimes are highly accurate and perfectly match the WFC traveltimes The resulting FD amplitudes are smooth and consistent with the geometrical spreading factors obtained from the ray approach Hence, both approaches can be used for fast and reliable computation of seismic first-arrival traveltimes and amplitudes in complex models In addition, the capabilities of the ray approach for computing traveltimes and spreading factors of later arrivals are demonstrated with the help of the Shell benchmark model

Unsupervised separation of seismic waves using the watershed algorithm on time-scale images

Abstract: This paper illustrates the use of image processing techniques for separating seismic waves. Because of the non-stationarity of seismic signals, the continuous wavelet transform is more suitable than the conventional Fourier transforms for the representation, and thus the analysis, of seismic processes. It provides a 2D representation, called a scalogram, of a 1D signal where the seismic events are well localized and isolated. Supervised methods based on this time-scale representation have already been used to separate seismic events, but they require strong interactions with the geophysicist. This paper focuses on the use of the watershed algorithm to segment time-scale representations of seismic signals, which leads to an automatic estimation of the wavelet representation of each wave separately. The computation of the inverse wavelet transform then leads to the reconstruction of the different waves. This segmentation, tracked over the different traces of the seismic profile, enables an accurate separation of the different wavefields. This method has been successfully validated on several real data sets.

3D imaging of a reservoir analogue in point bar deposits in the Ferron Sandstone, Utah, using ground-penetrating radar

Abstract: Most existing reservoir models are based on 2D outcrop studies; 3D aspects are inferred from correlation between wells, and so are inadequately constrained for reservoir simulations. To overcome these deficiencies, we have initiated a multidimensional characterization of reservoir analogues in the Cretaceous Ferron Sandstone in Utah. Detailed sedimentary facies maps of cliff faces define the geometry and distribution of reservoir flow units, barriers and baffles at the outcrop. High-resolution 2D and 3D ground-penetrating radar (GPR) images extend these reservoir characteristics into 3D to allow the development of realistic 3D reservoir models. Models use geometric information from mapping and the GPR data, combined with petrophysical data from surface and cliff-face outcrops, and laboratory analyses of outcrop and core samples. The site of the field work is Corbula Gulch, on the western flank of the San Rafael Swell, in east-central Utah. The outcrop consists of an 8–17 m thick sandstone body which contains various sedimentary structures, such as cross-bedding, inclined stratification and erosional surfaces, which range in scale from less than a metre to hundreds of metres. 3D depth migration of the common-offset GPR data produces data volumes within which the inclined surfaces and erosional surfaces are visible. Correlation between fluid permeability, clay content, instantaneous frequency and instantaneous amplitude of the GPR data provides estimates of the 3D distribution of fluid permeability and clay content.

The influence of induced polarization on long‐offset transient electromagnetic data

Abstract: The diffusion of electromagnetic fields is dependent not only on conductivity, but also on magnetic permeability, dielectric permittivity and polarizability, i.e. dispersive conductivity. The long‐offset transient electromagnetic (LOTEM) method is mainly used to determine the spatial distribution of conductivity in the subsurface. However, earlier work on loop‐loop TEM suggests that transient EM methods can also be affected by induced polarization (IP). Numerous 1D forward calculations were carried out to study the IP effect on LOTEM data, using the Cole‐Cole relaxation model to simulate the polarizability of the ground. Besides the polarizability of each layer, the IP effect depends on the LOTEM field set‐up and the spatial distribution of conductivity in the ground. In particular, near‐surface layers with high chargeabilities can significantly distort the late time transients of the electric field components in the vicinity of the transmitter. The influence of polarizable layers on the magnetic field components can be neglected under normal circumstances. In 1997 and 1999, LOTEM measurements were carried out at Mt. Vesuvius in Italy to explore the geological structure of the volcano. Sensitivity studies on the effect of polarizable layers suggest that high chargeabilities in connection with conductive layers at greater depths would result in a detectable distortion of the electric field transients. Although the simultaneous IP measurements revealed high chargeabilities in a near‐surface layer, no evidence of IP effects could be found in the measured LOTEM data. We conclude that the observed chargeabilities are local and that 3D effects are probably present in the data. Another aspect is the measurement of the system response, which is usually measured by placing a receiver very close to the transmitter. Therefore, large distortions can be expected if near‐surface polarizable layers exist. This was verified in practice by field measurements in an area with high chargeabilities in Longerich, Cologne.

Migration velocity analysis by depth image‐wave remigration: first results

Abstract: The image-wave equation for depth remigration is a partial differential equation that is similar to the acoustic wave equation. In this work, we study its finite-difference solution and possible applications. The conditions for stability, dispersion and dissipation exhibit a strong wavenumber dependence. Where higher horizontal than vertical wavenumbers are present in the data to be remigrated, stability may be difficult to achieve. Grid dispersion and dissipation can only be reduced to acceptable levels by the choice of very small grid intervals. Numerical tests demonstrate that, upon reaching the true medium velocity, remigrated images of curved reflectors propagate to the correct depth and those of diffractions collapse to single points. The latter property points towards the method's potential for use as a tool for migration velocity analysis. A first application to inhomogeneous media shows that in a horizontally layered medium, the reflector images reach their true depth when the remigration velocity equals the inverse of the mean medium slowness.

Velocity modelling and prestack depth imaging below complex salt structures: a case history from on-shore Germany

Abstract: In recent years, the advances in velocity model building and depth imaging have provided a better understanding of complex subsalt plays. The tomographic approach to subsurface velocity modelling, using interpretive processes, has led to significant progress in solving subsalt imaging problems, which were once considered to be impenetrable barriers. We show how gravity data, as an alternative data source, can be integrated into iterative velocity–depth model building to constrain the overburden velocity model and delineate the shape of the salt body above the target reflector. In this way, a structurally accurate image of subsalt reflectors is achieved.

Optimal use of PP and PS time-lapse stacks for fluid–pressure discrimination

Abstract: The combined use of time-lapse PP and PS seismic data is analysed for optimal discrimination between pressure and saturation changes. The theory is based on a combination of the well-known Gassmann model and the geomechanical grain model derived by Hertz and Mindlin. A key parameter in the discrimination process is the opening angle between curves representing constant changes in PP and PS reflectivity plotted against pressure and saturation changes. The optimal discrimination angle in the pressure-saturation space is 90° and this is used to determine optimal offset ranges for both PP and PS data. For typical production scenarios, we find an optimal offset range corresponding to an angle of incidence of 25-30°, for both PP and PS data. For gas we find slightly different results. This means that conventional survey parameters used in marine multicomponent acquisition should be sufficient for the purpose of estimating pressure and fluid saturation changes during production.

Joint cross-well and single-well seismic studies of CO2 injection in an oil reservoir

Abstract: A series of time-lapse seismic cross-well and single-well experiments were conducted in a diatomite reservoir to monitor the injection of CO2 into a hydrofracture zone, based on P- and S-wave data. A high-frequency piezo-electric P-wave source and an orbital-vibrator S-wave source were used to generate waves that were recorded by hydrophones as well as 3-component geophones. During the first phase the set of seismic experiments was conducted after the injection of water into the hydrofractured zone. The set of seismic experiments was repeated after a time period of seven months during which CO2 was injected into the hydrofractured zone. The questions to be answered ranged from the detectability of the geological structure in the diatomic reservoir to the detectability of CO2 within the hydrofracture. Furthermore, it was intended to determine which experiment (cross-well or single-well) is best suited to resolve these features. During the pre-injection experiment, the P-wave velocities exhibited relatively low values between 1700 and 1900 m/s, which decreased to 1600‐1800 m/s during the post-injection phase (−5%). The analysis of the pre-injection S-wave data revealed slow S-wave velocities between 600 and 800 m/s, while the post-injection data revealed velocities between 500 and 700 m/s (−6%). These velocity estimates produced high Poisson’s ratios between 0.36 and 0.46 for this highly porous (∼50%) material. Differencing post- and pre-injection data revealed an increase in Poisson’s ratio of up to 5%. Both velocity and Poisson’s ratio estimates indicate the dissolution of CO2 in the liquid phase of the reservoir accompanied by an increase in pore pressure. The single-well data supported the findings of the cross-well experiments. P- and S-wave velocities as well as Poisson’s ratios were comparable to the estimates of the cross-well data. The cross-well experiment did not detect the presence of the hydrofracture but appeared to be sensitive to overall changes in the reservoir and possibly the presence of a fault. In contrast, the single-well reflection data revealed an arrival that could indicate the presence of the hydrofracture between the source and receiver wells, while it did not detect the presence of the fault, possibly due to out-of-plane reflections.

3D sparse‐data Kirchhoff redatuming

Abstract: Wave-equation redatuming can be a very efficient method of overcoming the overburden imprint on the target area. Owing to the growing amount of 3D data, it is increasingly important to develop a feasible method for the redatuming of 3D prestack data. Common 3D acquisition designs produce relatively sparse data sets, which cannot be redatumed successfully by applying conventional wave-equation redatuming. We propose a redatuming approach that can be used to perform wave-equation redatuming of sparse 3D data. In this new approach, additional information about the medium velocity below the new datum is included, i.e. redatumed root-mean-square (RMS) velocities, which can be extracted from the input data set by conventional velocity analysis, are used. Inclusion of this additional information has the following implications: (i) it becomes possible to simplify the 4D redatuming integral into a 2D integral such that the number of traces needed to calculate one output time sample and the computational effort are both reduced; (ii) the information about the subsurface enables an infill of traces which are needed for the integral calculation but which are missing in the sparse input data set. Two tests applying this new approach to fully sampled 2D data show satisfactory results, implying that this method can certainly be used for the redatuming of sparse 3D data sets.