Applications of Directional Wavefield Decomposition, Phase Space, and Path Integral Methods to Seismic Wave Propagation and Inversion

— Recently, de Hoop and coworkers developed an asymptotic, seismic inversion formula for application in complex environments supporting multi-pathed and multi-mode wave propagation (de Hoop et al., 1999; de Hoop and Brandsberg-Dahl, 2000; Stolk and de Hoop, 2000). This inversion is based on the Born/Kirchhoff approximation, and employs the global, uniform asymptotic extension of the geometrical method of “tracing rays” to account for caustic phenomena. While this approach has successfully inverted the multicomponent, ocean-bottom data from the Valhall field in Norway, accounting for severe focusing effects (de Hoop and Brandsberg-Dahl, 2000), it is not able to account properly for wave phenomena neglected in the “high-frequency” limit (i.e., diffraction effects) and strong scattering effects. To proceed further and incorporate wave effects in a nonlinear inversion scheme, the theory of directional wavefield decomposition and the construction of the generalized Bremmer coupling series are combined with the application of modern phase space and path (functional) integral methods to, ultimately, suggest an inversion algorithm which can be interpreted as a method of “tracing waves.” This paper is intended to provide the seismic community with an introduction to these approaches to direct and inverse wave propagation and scattering, intertwining some of the most recent new results with the basic outline of the theory, and culminating in an outline of the extended, asymptotic, seismic inversion algorithm. Modeling at the level of the fixed-frequency (elliptic), scalar Helmholtz equation, exact and uniform asymptotic constructions of the well-known, and fundamentally important, square-root Helmholtz operator (symbol) provide the most important results.     

T-matrix approach to seismic forward modelling in the acoustic approximation

Forward seismic modelling in the acoustic approximation, for variable velocity but constant density, is dealt with. The wave equation and the boundary conditions are represented by a volume integral equation of the Lippmann-Schwinger (LS) or Fredholm type. A T-matrix (or transition operator) approach from quantum mechanical potential scattering theory is used to derive a family of linear and nonlinear approximations (cluster expansions), as well as an exact numerical solution of the LS equation. For models of 4D anomalies involving small or moderate contrasts, the Born approximation gives identical numerical results as the first-order t-matrix approximation, but the predictions of an exact T-matrix solution can be quite different (depending on spatial extention of the perturbations). For models of fluid-saturated cavities involving large or huge contrasts, the first-order t-matrix approximation is much more accurate than the Born approximation, although it does not lead to significantly more time-consuming computations. If the spatial extention of the perturbations is not too large, it is practical to use the exact T-matrix solution which allows for arbitrary contrasts and includes all the effects of multiple scattering.     

2.5维地震波场褶积微分算子法数值模拟

早期的褶积微分算子都是基于正反傅立叶变换而实现的,其精度比四阶有限差分的精度稍高,本文将计算数学中的Forsyte广义正交多项式微分算子与褶积算子相结合,构建了一个新的快速、高精度褶积微分算子,其计算结果非常接近实验函数微分的精确值,精度与16阶有限差分的精度相当,远优于错格伪谱法的精确度.另外,2.5维数值模拟比二维模拟可以更真实地模拟三维介质的臬个剖面的波场,并且2.5维地震波模拟的计算量比三维模拟的计算量及计算耗时要大大减少.本文利用基于Forsyte广义正交多项式褶积微分算子法计算2.5维非均匀介质地震波场,模拟结果表明,该算法的计算速度快,计算精度高,能够直观、高效地反映复杂介质中波场的传播规律,并且2.5维波场数值模拟具有更高的计算效率,是一种非常值得深入研究并广泛应用的方法.     

3-D acoustic modelling of edge diffractions — revisited

An accurate, fast, and simple algorithm for 3-D acoustic modelling of seismic edge diffractions, originally developed in the 1980s, is revisited in this paper. The main objective is to reintroduce this simple approach to edge-diffraction modelling and for the first time give the details of the theory in the open literature. The method is based on a combination of Kirchhoff theory and uniform asymptotic techniques developed within a high-frequency assumption. The diffraction contributions are then computed at stationary edge points only, by analogy with the geometrical ray contributions associated with internal stationary points or specular points. To be able to handle sampling inaccuracies of the critical edge points, a modified algorithm is proposed. Its robustness is verified in case of scattering from a circular edge. Also the extension from rigid or free boundary conditions to the case of edges defined by two penetrable surfaces is discussed in this paper. Both experimental and synthetic 3-D data are presented to demonstrate the potential of this edge-diffraction modelling technique. Since all parameters needed in the computations are obtained from dynamic ray tracing, the algorithm can readily be incorporated in existing software packages for 3-D seismic ray modelling.     

Time evolution of the first‐order linear acoustic/elastic wave equation using Lie product formula and Taylor expansion

We propose a new numerical solution to the first‐order linear acoustic/elastic wave equation. This numerical solution is based on the analytic solution of the linear acoustic/elastic wave equation and uses the Lie product formula, where the time evolution operator of the analytic solution is written as a product of exponential matrices where each exponential matrix term is then approximated by Taylor series expansion. Initially, we check the proposed approach numerically and then demonstrate that it is more accurate to apply a Taylor expansion for the exponential function identity rather than the exponential function itself. The numerical solution formulated employs a recursive procedure and also incorporates the split perfectly matched layer boundary condition. Thus, our scheme can be used to extrapolate wavefields in a stable manner with even larger time‐steps than traditional finite‐difference schemes. This new numerical solution is examined through the comparison of the solution of full acoustic wave equation using the Chebyshev expansion approach for the matrix exponential term. Moreover, to demonstrate the efficiency and applicability of our proposed solution, seismic modelling results of three geological models are presented and the processing time for each model is compared with the computing time taking by the Chebyshev expansion method. We also present the result of seismic modelling using the scheme based in Lie product formula and Taylor series expansion for the first‐order linear elastic wave equation in vertical transversely isotropic and tilted transversely isotropic media as well. Finally, a post‐stack migration results are also shown using the proposed method.     

Diffraction imaging by uniform asymptotic theory and double exponential fitting

Seismic diffracted waves carry valuable information for identifying geological discontinuities. Unfortunately, the diffraction energy is generally too weak, and standard seismic processing is biased to imaging reflection. In this paper, we present a dynamic diffraction imaging method with the aim of enhancing diffraction and increasing the signal‐to‐noise ratio. The correlation between diffraction amplitudes and their traveltimes generally exists in two forms, with one form based on the Kirchhoff integral formulation, and the other on the uniform asymptotic theory. However, the former will encounter singularities at geometrical shadow boundaries, and the latter requires the computation of a Fresnel integral. Therefore, neither of these methods is appropriate for practical applications. Noting the special form of the Fresnel integral, we propose a least‐squares fitting method based on double exponential functions to study the amplitude function of diffracted waves. The simple form of the fitting function has no singularities and can accelerate the calculation of diffraction amplitude weakening coefficients. By considering both the fitting weakening function and the polarity reversal property of the diffracted waves, we modify the conventional Kirchhoff imaging conditions and formulate a diffraction imaging formula. The mechanism of the proposed diffraction imaging procedure is based on the edge diffractor, instead of the idealized point diffractor. The polarity reversal property can eliminate the background of strong reflection and enhance the diffraction by same‐phase summation. Moreover,the fitting weakening function of diffraction amplitudes behaves like an inherent window to optimize the diffraction imaging aperture by its decaying trend. Synthetic and field data examples reveal that the proposed diffraction imaging method can meet the requirement of high‐resolution imaging, with the edge diffraction fully reinforced and the strong reflection mostly eliminated.     

Hydrodynamic effects in a reservoir with semicircular cross-section and absorptive bottom

A three-dimensional dam-reservoir system under seismic load is analysed. The dam is assumed to be rigid. The reservoir is an infinite channel with semi-circular cross-section. The exact analytical solution, based on the assumption of potential fluid motion is presented, as well as numerical results for selected parameters.The most significant parameters are: the direction and frequency content of the seismic input; the radiation damping at the reservoir bottom; and the compressibility of the fluid. The response of the system depends strongly on the direction of the input ground motion. This is shown by the transfer functions as well as by the pressure time histories due to two earthquakes with different frequency content. The energy absorption at the reservoir bottom is important. A simple plane-wave model shows, that even for a rock foundation, the amount of transmitted energy can reach up to 80%. For comparison the case without bottom absorption is also shown. Compressbility has to be included to capture the resonance effects. The exact analytical solution is also used to verify numerical results obtained by a new method that combines a finite element model with a rigorous radiation boundary for the infinite channel in the time domain.     

How rough sea affects marine seismic data and deghosting procedures

Most seismic processing algorithms generally consider the sea surface as a flat reflector. However, acquisition of marine seismic data often takes place in weather conditions where this approximation is inaccurate. The distortion in the seismic wavelet introduced by the rough sea may influence (for example) deghosting results, as deghosting operators are typically recursive and sensitive to the changes in the seismic signal. In this paper, we study the effect of sea surface roughness on conventional (5–160 Hz) and ultra‐high‐resolution (200–3500 Hz) single‐component towed‐streamer data. To this end, we numerically simulate reflections from a rough sea surface using the Kirchhoff approximation. Our modelling demonstrates that for conventional seismic frequency band sea roughness can distort results of standard one‐dimensional and two‐dimensional deterministic deghosting. To mitigate this effect, we introduce regularisation and optimisation based on the minimum‐energy criterion and show that this improves the processing output significantly. Analysis of ultra‐high‐resolution field data in conjunction with modelling shows that even relatively calm sea state (i.e., 15 cm wave height) introduces significant changes in the seismic signal for ultra‐high‐frequency band. These changes in amplitude and arrival time may degrade the results of deghosting. Using the field dataset, we show how the minimum‐energy optimisation of deghosting parameters improves the processing result.     

Application of the nearly perfectly matched layer for seismic wave propagation in 2D homogeneous isotropic media

Numerical modelling plays an important role in helping us understand the characteristics of seismic wave propagation. The presence of spurious reflections from the boundaries of the truncated computational domain is a prominent problem in finite difference computations. The nearly perfectly matched layer has been proven to be a very effective boundary condition to absorb outgoing waves in both electromagnetic and acoustic media. In this paper, the nearly perfectly matched layer technique is applied to elastic isotropic media to further test the method's absorbing ability. The staggered‐grid finite‐difference method (fourth‐order accuracy in space and second‐order accuracy in time) is used in the numerical simulation of seismic wave propagation in 2D Cartesian coordinates. In the numerical tests, numerical comparisons between the nearly perfectly matched layer and the convolutional perfectly matched layer, which is considered the best absorbing layer boundary condition, is also provided. Three numerical experiments demonstrate that the nearly perfectly matched layer has a similar performance to the convolutional perfectly matched layer and can be a valuable alternative to other absorbing layer boundary conditions.     

Improved seismic hazard model with application to probabilistic seismic demand analysis

An improved seismic hazard model for use in performance‐based earthquake engineering is presented. The model is an improved approximation from the so‐called ‘power law’ model, which is linear in log–log space. The mathematics of the model and uncertainty incorporation is briefly discussed. Various means of fitting the approximation to hazard data derived from probabilistic seismic hazard analysis are discussed, including the limitations of the model. Based on these ‘exact’ hazard data for major centres in New Zealand, the parameters for the proposed model are calibrated. To illustrate the significance of the proposed model, a performance‐based assessment is conducted on a typical bridge, via probabilistic seismic demand analysis. The new hazard model is compared to the current power law relationship to illustrate its effects on the risk assessment. The propagation of epistemic uncertainty in the seismic hazard is also considered. To allow further use of the model in conceptual calculations, a semi‐analytical method is proposed to calculate the demand hazard in closed form. For the case study shown, the resulting semi‐analytical closed form solution is shown to be significantly more accurate than the analytical closed‐form solution using the power law hazard model, capturing the ‘exact’ numerical integration solution to within 7% accuracy over the entire range of exceedance rate. Copyright © 2007 John Wiley & Sons, Ltd.     

Complex frequency-shifted multi-axial perfectly matched layer for frequency-domain seismic wavefield simulation in anisotropic media

Seismic anisotropy has an important influence on seismic data processing and interpretation. Although the frequency-domain seismic wavefield simulation has a problem of solving the large scale linear sparse matrix due to the computational limitations, it has some advantages over the time-domain seismic wavefield simulation including efficient inversion using only a limited number of frequency components and easy implementation of multiple sources. To accurately simulate seismic wave propagation in the frequency domain, we also need to choose the absorbing boundary conditions to absorb artificial reflections from edges of the model as we do in the time domain. Compared with the classical boundary conditions including the perfectly matched layer and complex frequency-shifted perfectly matched layer, the complex frequency-shifted multi-axial perfectly matched layer has been proven to effectively suppress the unwanted reflections at grazing incidence and solve the instability problem in the time-domain seismic numerical modelling in anisotropic elastic media. In this paper, we propose to extend the complex frequency-shifted multi-axial perfectly matched layer absorbing boundary condition to the frequency-domain seismic wavefield simulation in anisotropic elastic media. To test the validity of our proposed algorithm, we compare the results (snapshots and seismograms) of the frequency-domain seismic wavefield simulation with those of the time-domain modelling. The model studies indicate that the complex frequency-shifted multi-axial perfectly matched layer absorbing boundary condition is stable in the frequency-domain seismic wavefield simulation in anisotropic media, and provides better absorbing performance than the complex frequency-shifted perfectly matched layer boundary condition.     

Perfectly matched layer boundary conditions for the second-order acoustic wave equation solved by the rapid expansion method

We derive a governing second-order acoustic wave equation in the time domain with a perfectly matched layer absorbing boundary condition for general inhomogeneous media. Besides, a new scheme to solve the perfectly matched layer equation for absorbing reflections from the model boundaries based on the rapid expansion method is proposed. The suggested scheme can be easily applied to a wide class of wave equations and numerical methods for seismic modelling. The absorbing boundary condition method is formulated based on the split perfectly matched layer method and we employ the rapid expansion method to solve the derived new perfectly matched layer equation. The use of the rapid expansion method allows us to extrapolate wavefields with a time step larger than the ones commonly used by traditional finite-difference schemes in a stable way and free of dispersion noise. Furthermore, in order to demonstrate the efficiency and applicability of the proposed perfectly matched layer scheme, numerical modelling examples are also presented. The numerical results obtained with the put forward perfectly matched layer scheme are compared with results from traditional attenuation absorbing boundary conditions and enlarged models as well. The analysis of the numerical results indicates that the proposed perfectly matched layer scheme is significantly effective and more efficient in absorbing spurious reflections from the model boundaries.     

Exact time-domain solutions for acoustic diffraction by a half plane

We derive exact time-domain solutions for scattering of acoustic waves by a half plane by inverse Fourier transforming the frequency-domain integral solutions. The solutions consist of a direct term, a reflected term and two diffraction terms. The diffracting edge induces step function discontinuities in the direct and reflected, terms at two shadow boundries. At each boundary, the associated diffraction term reaches a maximum amplitude of half the geometrical optics term and has a signum function discontinuity so that the total field remains continuous. We evaluate solutions for practical point source configurations by numerically convolving the impulse diffraction responses with a wavelet. We solve the associated problems of convolution with a singular, truncated diffraction operator by analytically derived correction techniques. We produce a zero offset section and compare it to a Kirchhoff integral solution. Our exact diffraction hyperbola exhibits noticeable asymmetry, with higher amplitudes on the reflector side of the edge. Near the apex of the hyperbola the Kirchhoff solution approximates the exact diffraction term symmetric in amplitude about the reflection shadow boundary, but omits the other low amplitude term necessary to ensure continuity at the direct shadow boundary.     

Using the pseudospectral technique on curved grids for 2D acoustic forward modelling1

The Fourier pseudospectral method has been widely accepted for seismic forward modelling because of its high accuracy compared to other numerical techniques. Conventionally, the modelling is performed on Cartesian grids. This means that curved interfaces are represented in a ‘staircase fashion‘causing spurious diffractions. It is the aim of this work to eliminate these non-physical diffractions by using curved grids that generally follow the interfaces. A further advantage of using curved grids is that the local grid density can be adjusted according to the velocity of the individual layers, i.e. the overall grid density is not restricted by the lowest velocity in the subsurface. This means that considerable savings in computer storage can be obtained and thus larger computational models can be handled. One of the major problems in using the curved grid approach has been the generation of a suitable grid that fits all the interfaces. However, as a new approach, we adopt techniques originally developed for computational fluid dynamics (CFD) applications. This allows us to put the curved grid technique into a general framework, enabling the grid to follow all interfaces. In principle, a separate grid is generated for each geological layer, patching the grid lines across the interfaces to obtain a globally continuous grid (the so-called multiblock strategy). The curved grid is taken to constitute a generalised curvilinear coordinate system, where each grid line corresponds to a constant value of one of the curvilinear coordinates. That means that the forward modelling equations have to be written in curvilinear coordinates, resulting in additional terms in the equations. However, the subsurface geometry is much simpler in the curvilinear space. The advantages of the curved grid technique are demonstrated for the 2D acoustic wave equation. This includes a verification of the method against an analytic reference solution for wedge diffraction and a comparison with the pseudospectral method on Cartesian grids. The results demonstrate that high accuracies are obtained with few grid points and without extra computational costs as compared with Cartesian methods.     

Numerical-analytical interfacing in two dimensions with applications to modeling NTS seismograms

A new method for interfacing numerical and integral techniques allows greater flexibility in seismic modeling. Specifically, numerical calculations in laterally varying structure are interfaced with analytic methods that enable propagation to great distances. Such modeling is important for studying situations containing localized complex regions not easily handled by analytic means. The calculations involved are entirely two-dimensional, but the use of an appropriate source in combination with a filter applied to the resulting seismograms produces synthetic seismograms which are point-source responses in three dimensions. The integral technique is called two-dimensional Kirchhoff because its form is similar to the classical three-dimensional Kirchhoff. Data from Yucca Flat at the Nevada Test Site are modeled as a demonstration of the usefulness of the new method. In this application, both local and teleseismic records are modeled simultaneously from the same model with the same finite-difference run. This application indicates the importance of locally scattered Rayleigh waves in the production of teleseismic body-wave complexity and coda.     

Local and controlled prestack depth migration in complex areas

Prestack depth migrations based on wavefield extrapolation may be computationally expensive, especially in 3D. They are also very dependent on the acquisition geometry and are not flexible regarding the geometry of the imaging zone. Moreover, they do not deal with all types of wave, considering only primary reflection events through the model. Integral approaches using precalculated Green's functions, such as Kirchhoff migration and Born-based imaging, may overcome these problems. In the present paper, both finite-difference traveltimes and wavefront construction are used to obtain asymptotic Green's functions, and a generalized diffraction tomography is applied as an example of Born-based acoustic imaging. Target-orientated imaging is easy to perform, from any type of survey and subselection of shot/receiver pairs. Multifield imaging is possible using Green's functions that take into account, for instance, reflections at model boundaries. This may help to recover parts of complex structures which would be missing using a paraxial wave equation approach. Finally, a numerical evaluation of the resolution, or point-spread, function at any point of the depth-migrated section provides valuable information, either at the survey planning stage or for the interpretation.     

From the Hagedoorn imaging technique to Kirchhoff migration and inversion

The seminal 1954 paper by J.G. Hagedoorn introduced a heuristic for seismic reflector imaging. That heuristic was a construction technique – a 'string construction' or 'ruler and compass' method – for finding reflectors as an envelope of equal traveltime curves defined by events on a seismic trace. Later, Kirchhoff migration was developed. This method is based on an integral representation of the solution of the wave equation. For decades Kirchhoff migration has been one of the most popular methods for imaging seismic data. Parallel with the development of Kirchhoff wave-equation migration has been that of Kirchhoff inversion, which has as its objectives both structural imaging and the recovery of angle-dependent reflection coefficients. The relationship between Kirchhoff migration/inversion and Hagedoorn's constructive technique has only recently been explored. This paper addresses this relationship, presenting the mathematical structure that the Kirchhoff approach adds to Hagedoorn's constructive method and showing the relationship between the two.     

TIME MIGRATION—SOME RAY THEORETICAL ASPECTS*

Using an elementary theory of migration one can consider a reflecting horizon as a continuum of scattering centres for seismic waves. Reflections arising at interfaces can thus be looked upon as the sum of energy scattered by interface points. The energy from one point is distributed among signals upon its reflection time surface. This surface is usually well approximated by a hyperboloid in the vicinity of its apex. Migration aims at focusing the scattered energy of each depth point into an image point upon the reflection time surface. To ensure a complete migration the image must be vertical above the depth point. This is difficult to achieve for subsurface interfaces which fall below laterally in-homogeneous velocity media. Migration is hence frequently performed for these interfaces as well by the Kirchhoff summation method which systematically sums signals into the apex of the approximation hyperboloid even though the Kirchhoff integral is in this case not strictly valid. For a multilayered subsurface isovelocity layer model with interfaces of a generally curved nature this can only provide a complete migration for the uppermost interface. Still there are various advantages gained by having a process which sums signals consistently into the minimum of the reflection time surface. The position of the time surface minimum is the place where a ray from the depth point emerges vertically to the surface. The Kirchhoff migration, if applied to media with laterally inhomogeneous velocity, must necessarily be followed by a further time-to-depth migration if the true depth structure is to be recovered. Primary normal reflections and their respective migrated reflections have a complementary relationship to each other. Normal reflections relate to rays normal to the reflector and migrated reflections relate to rays normal to the free surface. Ray modeling is performed to indicate a new approach for simulating seismic reflections. Commonly occuring situations are investigated from which lessons can be learned which are of immediate value for those concerned with interpreting time migrated reflections. The concept of the ‘image ray’ is introduced.     

论绕射波机制

本文讨论一项用以研究折射绕射波的三度空间超声波模拟实验的结果。实验说明,在三度空间中,在某一点观测到的绕射波包含着来自无穷个绕射源的波涟,这些绕射源沿着绕射楔或断层而分布。这点在观测点位于垂直绕射楔或断层的方向上也是同样真实。 结果也指出绕射波的动力学性质随着绕射楔或断层的几何形状而变化,因此根据这种结果可以设想,若详细研究绕射波的动力学特性,可以对绕射楔或断层的几何形状得出有用的证据。 此外,根据模拟实验的结果,可用“最小走程原理”来解释波涟到达观测点的行程和求得绕射点的位置。