1 

Geophysical modelling of 3D electromagnetic diffusion with multigrid
The performance of a multigrid solver for timeharmonic electromagnetic problems in geophysical settings was investigated. With the low frequencies used in geophysical surveys for deeper targets, the lightspeed waves in the earth can be neglected. Diffusion of induced currents is the dominant physical effect. The governing equations were discretised by the FiniteIntegration Technique. The resulting set of discrete equation was solved by a multigrid method. The multigrid method provided excellent convergence with constant grid spacings, but not on stretched grids. The slower convergence rate of the multigrid method could be compensated by using bicgstab2, in which case multigrid acted as a preconditioner. Still, the overall performance was less than satisfactory with substantial grid stretching.

[PDF]
[Abstract]

2 

Subsurface offset behaviour in velocity analysis with extended reflectivity images
Migration velocity analysis with the wave equation can be accomplished by focusing of extended migration images, obtained by introducing a subsurface offset or shift. A reflector in the wrong velocity model will show up as a curve in the extended image. In the correct model, it should collapse to a point. The usual approach to obtain a focused image involves a cost functional that penalizes energy in the extended image at nonzero shift. Its minimization by a gradientbased method should then produce the correct velocity model. Here, asymptotic analysis and numerical examples show that this method may be too sensitive to amplitude peaks at large shifts at the wrong depth and to artifacts. A more robust alternative is proposed that can be interpreted as a generalization of stack power and maximizes the energy at zero subsurface shift. A realdata example is included.

[PDF]
[Abstract]

3 

Onder de grond kijken

[PDF]

4 

New triangular masslumped finite elements of degree six for wave propagation
Masslumped continuous finite elements allow for explicit time stepping with the secondorder wave equation if the resulting integration weights are positive and provide sufficient accuracy. To meet these requirements on triangular and tetrahedral meshes, the construction of continuous finite elements for a given polynomial degree on the edges involves polynomials of higher degree in the interior. The parameters describing the supporting nodes of the Lagrange interpolating polynomials and the integration weights are
the unknowns of a polynomial system of equations, which is linear in the integration weights. To find candidate sets for the nodes, it is usually required that the number of equations equals the number of
unknowns, although this may be neither necessary nor sufficient. Here, this condition is relaxed by requiring that the number of equations does not exceed the number of unknowns. This resulted in two new types elements of degree 6 for symmetrically placed nodes. Unfortunately, the first type is not unisolvent. There are many different elements of the second type with a large range in their associated timestepping stability limit. To assess the efficiency of the elements of various degrees, numerical tests on a simple problem with an exact solution were performed. Efficiency was measured by the computational time
required to obtain a solution at a given accuracy. For the chosen example, elements of degree 4 with fourthorder time stepping appear to be the most efficient.

[PDF]
[Abstract]

5 

Inversion of seismic reflection data through focusing
Since most of the easy hydrocarbon reservoirs have been found, accurate subsurface imaging for oil and gas exploration and production is crucial. Seismic data can provide a bandlimited reconstruction of impedance contrasts between different rock formations as well as a subsurface velocity model. Current compute power allows for the use of the full acoustic wave equation for seismic imaging and processing. Leastsquares data fitting is an obvious approach but may provide the wrong answer because of local minima in the cost function caused by the absence of low frequencies in the data. An alternative formulation is based on extended images that invokes action at distance to make up for errors in the estimated subsurface velocity model. The action at distance is accomplished by a subsurface shift and the correct model should focus amplitudes at zero shift. The related cost function penalizes image amplitudes at nonzero shift. It has a large basin of attraction but loses its sharpness closer to the minimum. A new formulation is proposed – motivated by highfrequency asymptotic analysis – that provides better results on synthetic and real seismic data and is less sensitive to imaging artifacts.

[PDF]
[Abstract]

6 

A realdata example of automatic migration velocity analysis with extended images in the presence of multiples
Migration velocity analysis for the twoway wave equation based on focusing extended migration images at zero subsurface offset are sensitive to the presence of multiples. Even after thorough multiple suppression, remnant multiple energy can lead to conflicting events that will focus at different velocities, one for the primary and another for the multiple. An automatic algorithm will produce a velocity in between these two. If the multiples are mainly caused by the presence of a water layer, a bias of the algorithm towards higher velocities may help. Here, an application of this approach to a marine data is presented and the result is compared to a well log.

[PDF]
[Abstract]

7 

A Comparison of Triangular Masslumped Finite Elements for 2D Wave Propagation
Masslumped continuous finite elements provide accurate solutions of the secondorder acoustic or elastic wave equation in complex geological settings, in particular in the presence of topography and large impedance contrasts. Elements of higher polynomial degree have better accuracy than those of lower degree but are also more costly. Here, the performance of elements of degree 1 to 6 is compared for a simple acoustic test problem. The element of degree 6 is new. The numerical test confirm the increase of accuracy with the polynomial degree of the basis functions. In terms of computational cost, the element of degree 4 performs best in the specific example. For higher degrees, the cost of having extra nodes and a more restrictive stability constraint on the time step can no longer be compensated by having a smaller number of larger triangular elements. Only at very high accuracy, the new element of degree 6 wins.

[PDF]
[Abstract]

8 

Subsurface offset behaviour in velocity analysis with extended reflectivity images
Migration velocity analysis with the constantdensity acoustic wave equation can be accomplished by the focusing of extended migration images, obtained by introducing a subsurface shift in the imaging condition. A reflector in a wrong velocity model will show up as a curve in the extended image. In the correct model, it should collapse to a point. The usual approach to obtain a focused image involves a cost functional that penalizes energy in the extended image at nonzero shift. Its minimization by a gradientbased method should then produce the correct velocity model. Here, asymptotic analysis and numerical examples show that this method may be too sensitive to amplitude peaks at large shifts at the wrong depth and to artefacts. A more robust alternative is proposed that can be interpreted as a generalization of stack power and maximizes the energy at zerosubsurface shift. A realdata example is included.

[PDF]
[Abstract]

9 

An ambiguity in attenuation scattering imaging

[PDF]

10 

Seismic attenuation imaging with causality
Seismic data enable imaging of the Earth, not only of velocity and density but also of attenuation contrasts. Unfortunately, the Born approximation of the constantdensity viscoacoustic wave equation, which can serve as a forward modelling operator related to seismic migration, exhibits an ambiguity when attenuation is included. Different scattering models involving velocity and attenuation perturbations may provide nearly identical data. This result was obtained earlier for scatterers that did not contain a correction term for causality. Such a term leads to dispersion when considering a range of frequencies. We demonstrate that with this term, linearized inversion or iterative migration will almost, but not fully, remove the ambiguity. We also investigate if attenuation imaging suffers from the same ambiguity when using nonlinear or full waveform inversion. A numerical experiment shows that nonlinear inversion with causality convergences to the true model, whereas without causality, a substantial difference with the true model remains even after a very large number of iterations. For both linearized and nonlinear inversion, the initial update in a gradientbased optimization scheme that minimizes the difference between modelled and observed data is still affected by the ambiguity and does not provide a good result. This first update corresponds to a classic migration operation. In our numerical experiments, the reconstructed model started to approximate the true model only after a large number of iterations.

[PDF]
[Abstract]

11 

A multigrid based iterative solver for the frequency domain elastic wave equation

[PDF]

12 

A multigridbased iterative solver for the frequencydomain elastic wave equation
Efficient numerical wave modelling is essential for imaging methods such as reversetime migration (RTM) and full waveform inversion (FWI). In 2D, frequencydomain modelling with LU factorization as a direct solver can outperform timedomain methods by one order. For 3D problems, the computational complexity of the LU factorization as well as its memory requirements are a disadvantage and the time domain becomes more attractive. Recently, it has been shown that with compute cores in abundance, a parallel frequencydomain iterative solver can be a competitive alternative to the timedomain approach for 3D acoustic RTM. The solver relies on a preconditioned Krylov subspace method, where the preconditioner involves the multigrid solution of a heavily damped wave equation. Here, we generalize this idea to the isotropic elastic wave equation in the frequency domain.

[PDF]
[Abstract]

13 

Resistivity imaging with controlledsource electromagnetic data: depth and data weighting
We discuss some computational aspects of resistivity imaging by inversion of offshore controlledsource electromagnetic data. We adopt the classic approach to imaging by formulating it as an inverse problem. A weighted leastsquares functional measures the misfit between synthetic and observed data. Its minimization by a quasiNewton algorithm requires the gradient of the functional with respect to the model parameters. We compute the gradient with the adjointstate technique. Preconditioners can improve the convergence of the inversion. Diagonal preconditioner based on a Born approximation are commonly used. In the context of CSEM inversion, the Born approximation is not really accurate, this limits the possibility of estimating a correct approximation of the Hessian in a smooth medium or, in fact, in any reference background that does not roughly account for the resistors. We hence rely on the limited memory BFGS approximation of the inverse of the Hessian and we improve the inversion convergence with the help of a heuristic data and depth weighting. Based on a numerical example, we show that a simple exponential depth weighting combined with an offset or frequency data weighting significantly improves the convergence rate of a deepwater controlledsource electromagnetic data inversion.

[PDF]
[Abstract]

14 

A comparison of seismic velocity inversion methods for layered acoustics
In seismic imaging, one tries to infer the medium properties of the subsurface from seismic reflection data. These data are the result of an active source experiment, where an explosive source and an array of receivers are placed at the surface. Due to the absence of low frequencies in the data, the corresponding inverse problem is strongly nonlinear in the slowly varying component of the velocity. The leastsquares misfit functional typically exhibits local minima and has a small basin of attraction. The usual approach of fitting the data in a leastsquares sense by employing a gradientbased optimisation method will therefore most likely result in a wrong velocity model. In the geophysical community, this problem has long been recognised and alternative formulations of the inverse problem have been developed. We review several of these formulations and analyse the sensitivity to the error in the smooth velocity component. This analysis is carried out for laterally homogeneous velocities using an asymptotic solution of the wave equation. The analysis suggests that formulations which are geared towards fitting the phases of the data, rather than the amplitudes, have smooth corresponding misfit functionals with a large basin of attraction.

[PDF]
[Abstract]

15 

The diagonalator: An alternative cost functional for waveequation inversion
The classic leastsquares cost functional for full waveform inversion suffers from local minima due to loop skipping in the absence of low frequencies in the seismic data. Velocity model building based on subsurface spatial or temporal shifts may break down in the presence of multiples in the data. Cost functionals that translate this idea to the data domain, with offset or timeshifts, can handle multiples.
An earlier datadomain formulation suffered from crosstalk between events. Here, we present a multishot extension that should be less sensitive to crosstalk. It has the property of an annihilator, similar to the functional used for velocity analysis with extended images based on subsurface shifs. However, since it operates in the data domain, it should be able to handle multiples.
For 2D models with line acquistion, the proposed functional applies a singularvalue decomposition on the observed data and uses the eigenvectors to build data panel that should be diagonal in the correct velocity model, but has significant offdiagonal entries in the wrong model. By minimizing these offdiagonal entries or maximizing the main diagonal, the correct model should be found. We therefore named it the diagonalator.
We present initial tests on a simple, horizontally layered velocity model.

[PDF]
[Abstract]

16 

A correlationbased misfit criterion for waveequation traveltime tomography
Waveequation traveltime tomography tries to obtain a subsurface velocity model from seismic data, either passive or active, that explains their traveltimes. A key step is the extraction of traveltime differences, or relative phase shifts, between observed and modelled finitefrequency waveforms. A standard approach involves a correlation of the observed and measured waveforms. When the amplitude spectra of the waveforms are identical, the maximum of the correlation is indicative of the relative phase shift. When the amplitude spectra are not identical, however, this argument is no longer valid. We propose an alternative criterion to measure the relative phase shift. This misfit criterion is a weighted norm of the correlation and is less sensitive to differences in the amplitude spectra. For practical application it is important to use a sensitivity kernel that is consistent with the way the misfit is measured. We derive this sensitivity kernel and show how it differs from the standard banana–doughnut sensitivity kernel. We illustrate the approach on a crosswell data set.

[PDF]
[Abstract]

17 

Decoupling of elastic parameters with iterative linearized inversion
Three model parameters as a function of position describe wave propagation in an isotropic elastic medium. Ideally, imaging of data for a point scatterer that consists of a perturbation in one of the elastic parameters should only provide a reconstruction of that perturbation, without crosstalk into the other parameters. This is not the case for seismic migration, where a perturbation of one elastic parameter contributes to the images of all three model parameters. For a reliable reconstruction of the true elastic reflectivity, one can apply iterative migration or linearized inversion, where the misfit cost function is minimized by the conjugategradient method. We investigated the decoupling of the three isotropic elastic medium parameters with the iterative linearized approach. Instead of iterating, the final result can be obtained directly by means of Newton's method, using the pseudoinverse of the Hessian matrix. Although the calculation of the Hessian for a realistic model is an extremely resourceintensive problem, it is feasible for the simple case of a point scatterer in a homogeneous medium, for which we present numerical results. We consider the iterative approach with the conjugategradient method and Newton's method with the complete Hessian. Experiments show that in both cases the elastic parameters are decoupled much better when compared to migration. The iterative approach achieves acceptable inversion results but requires a large number of iterations. For faster convergence, preconditioning is required. An optimal preconditioner, if found, can be used in other iterative methods including LBFGS.We considered two well known types of preconditioners, based on diagonal and on blockdiagonal Hessian approximations. Somewhat to our surprise, both preconditioners fail to improve the convergence rate. Hence, a more sophisticated preconditioning is required.

[PDF]
[Abstract]

18 

On the timestepping stability of continuous masslumped and discontinuous Galerkin finite elements for the 3D acoustic wave equation
We solve the threedimensional acoustic wave equation, discretized on tetrahedral meshes. Two methods are considered: masslumped continuous finite elements and the symmetric interiorpenalty discontinuous Galerkin method (SIPDG). Combining the spatial discretization with the leapfrog timestepping scheme, which is secondorder accurate and conditionally stable, leads to a fully explicit scheme. We provide estimates of its stability limit for simple cases, namely, the reference element with Neumann boundary conditions, its distorted version of arbitrary shape, the unit cube that can be partitioned into 6 tetrahedra with periodic boundary conditions, and its distortions. The CFL stability limit contains a length scale for which we considered different options. The one based on the sum of the eigenvalues of the spatial operator for the first degree masslumped element gives the best results. It resembles the diameter of the inscribed sphere but is slightly easier to compute. The stability estimates show that masslumped continuous and SIPDG finite elements have comparable stability conditions, with the exception of the elements of the first degree. The stability limit for the masslumped elements is less restrictive and allows for larger time steps.

[PDF]
[Abstract]

19 

Timestepping stability of continuous and discontinuous finiteelement methods for 3D wave propagation
We analyse the timestepping stability for the 3D acoustic wave equation, discretized on tetrahedral meshes. Two types of methods are considered: masslumped continuous finite elements and the symmetric interiorpenalty discontinuous Galerkin method. Combining the spatial discretization with the leapfrog timestepping scheme, which is secondorder accurate and conditionally stable, leads to a fully explicit scheme. We provide estimates of its stability limit for simple cases, namely, the reference element with Neumann boundary conditions, its distorted version of arbitrary shape, the unit cube that can be partitioned into six tetrahedra with periodic boundary conditions and its distortions. The Courant–Friedrichs–Lewy stability limit contains an element diameter for which we considered different options. The one based on the sum of the eigenvalues of the spatial operator for the firstdegree masslumped element gives the best results. It resembles the diameter of the inscribed sphere but is slightly easier to compute. The stability estimates show that the masslumped continuous and the discontinuous Galerkin finite elements of degree 2 have comparable stability conditions, whereas the masslumped elements of degree one and three allow for larger time steps.

[PDF]
[Abstract]

20 

The impact of surfacewave separation on seismic survey design
3D seismic survey design provides an acquisition geometry for obtaining seismic data that enable imaging and amplitudeversusoffset applications of target reflectors with sufficient quality under given economical and operational constraints. However, in land or shallow water environments, surface waves are often dominant and will lower the quality of the final output. The necessity to remove them from the seismic data imposes additional constraints on the acquisition parameters.
Here, we try to understand how the application of surfacewave separation affects the choice of survey parameters and the resulting data quality. Quality is quantified by the signaltonoise ratio of the seismic data after surfacewave separation or removal. In a case study, we applied surfacewave separation and signaltonoise ratio estimation to several data sets with different survey parameters. We found that the spatial sampling intervals of the basic subset are the most important ones among the various types of survey parameters. The resulting data quality as a function of the spatial sampling intervals follows a trend curve. Finer spatial sampling intervals lead to better data quality up to a point where the curve flattens and a plateau is reached. The actual shape of the trend curve depends on the method of surfacewave separation.

[PDF]
[Abstract]
