Reciprocity theorems for one-way wavefields

Acoustic reciprocity theorems have proved their usefulness in the study of forward and inverse scattering problems. The reciprocity theorems in the literature apply to the two-way (i.e. total) wavefield, and are thus not compatible with one-way wave theory, which is often applied in seismic exploration. By transforming the two-way wave equation into a coupled system of one-way wave equations for downgoing and upgoing waves it appears to be possible to derive 'one-way reciprocity theorems" along the same lines as the usual derivation of the 'two-way reciprocity theorems'. However, for the one-way reciprocity theorems it is not directly obvious that the 'contrast term' vanishes when the medium parameters in the two different states are identical. By introducing a modal expansion of the Helraholtz operator, its square root can be derived, which appears to have a symmetric kernel. This symmetry property appears to be sufficient to let the contrast term vanish in the above-mentioned situation.
The one-way reciprocity theorem of the convolution type is exact, whereas the one-way reciprocity theorem of the correlation type ignores evanescent wave modes. The extension to the elastodynamic situation is not trivial, but it can be shown relatively easily that similar reciprocity theorems apply if the (non-unique) decomposition of the elastodynamic two-way operator is done in such a way that the elastodynamic one-way operators satisfy similar symmetry properties to the acoustic one-way operators.     

Teleseismic waveform modelling with a one-way wave equation

It is now widely accepted that elastic properties of the continental lithosphere and the underlying sublithospheric mantle are both anisotropic and laterally heterogeneous at a range of scales. To fully exploit modern three-component broad-band array data sets requires the use of comprehensive modelling tools. In this work, we investigate the use of a wide-angle, one-way wave equation to model variations in teleseismic 3-D waveforms due to 2-D elastic heterogeneity and anisotropy. The one-way operators are derived based on a high-frequency approximation of the square-root operator and include the effects of wave propagation as well as multiple scattering. Computational cost is reduced through a number of physically motivated approximations. We present synthetic results from simple 1-D (layer over a half-space) and 2-D (subduction zone) models that are compared with reference solutions. The algorithm is then used to model data from an array of broad-band seismograph stations deployed in northwestern Canada as part of the IRIS-PASSCAL/LITHOPROBE CANOE experiment. In this region radial-component receiver functions show a clear continental Moho and the presence of crustal material dipping into the mantle at the suture of two Palaeo-Proterozoic terranes. The geometry of the suture is better defined on the transverse component where subduction is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. The modelling reproduces the main features of the receiver functions, including the effects of anisotropy, heterogeneity and finite-frequency scattering.     

Quadrangle-grid velocity-stress finite-difference method for elastic-wave-propagation simulation

I present a 2-D numerical-modelling algorithm based on a first-order velocity-stress hyperbolic system and a non-rectangular-grid finite-difference operator. In this method the velocity and stress are defined at different nodes for a staggered grid. The scheme uses non-orthogonal grids, thereby surface topography and curved interfaces can be easily modelled in the seismic-wave-propagation stimulation. The free-surface conditions of complex geometry are achieved by using integral equilibrium equations on the surface, and the stability of the free-surface conditions is improved by introducing local filter modification. The method incorporates desirable qualities of the finite-element method and the staggered-grid finite-difference scheme, which is of high accuracy and low computational cost.     

Three-dimensional massively parallel electromagnetic inversion—I. Theory

An iterative solution to the non-linear 3-D electromagnetic inverse problem is obtained by successive linearized model updates using the method of conjugate gradients. Full wave equation modelling for controlled sources is employed to compute model sensitivities and predicted data in the frequency domain with an efficient 3-D finite-difference algorithm. Necessity dictates that the inverse be underdetermined, since realistic reconstructions require the solution for tens of thousands of parameters. In addition, large-scale 3-D forward modelling is required and this can easily involve the solution of over several million electric field unknowns per solve. A massively parallel computing platform has therefore been utilized to obtain reasonable execution times, and results are given for the 1840-node Intel Paragon. The solution is demonstrated with a synthetic example with added Gaussian noise, where the data were produced from an integral equation forward-modelling code, and is different from the finite difference code embedded in the inversion algorithm     

One-way representations of seismic data

A general one-way representation of seismic data can be obtained by substituting a Green's one-way wavefield matrix into a reciprocity theorem of the convolution type for one-way wavefields. From this general one-way representation, several special cases can be derived.
By introducing a Green's one-way wavefield matrix for primaries , a generalized Bremmer series representation is obtained. Terminating this series after the first-order term yields a primary representation of seismic reflection data. According to this representation, primary seismic reflection data are proportional to a reflection operator, 'modified' by primary propagators for downgoing and upgoing waves. For seismic imaging, these propagators need to be inverted. Stable inverse primary propagators can easily be obtained from a one-way reciprocity theorem of the correlation type.
By introducing a Green's one-way wavefield matrix for generalized primaries , an alternative representation is obtained in which multiple scattering is organized quite differently (in comparison with the generalized Bremmer series representation). According to the generalized primary representation, full seismic reflection data are proportional to a reflection operator, 'modified' by generalized primary propagators for downgoing and upgoing waves. Internal multiple scattering is fully included in the generalized primary propagators {either via a series expansion or in a parametrized way). Stable inverse generalized primary propagators can be obtained from the one-way reciprocity theorem of the correlation type. These inverse propagators are the nucleus for seismic imaging techniques that take the angle-dependent dispersion effects due to fine-layering into account.     

Surface wave tomography: global membrane waves and adjoint methods

We implement the wave equation on a spherical membrane, with a finite-difference algorithm that accounts for finite-frequency effects in the smooth-Earth approximation, and use the resulting 'membrane waves' as an analogue for surface wave propagation in the Earth. In this formulation, we derive fully numerical 2-D sensitivity kernels for phase anomaly measurements, and employ them in a preliminary tomographic application. To speed up the computation of kernels, so that it is practical to formulate the inverse problem also with respect to a laterally heterogeneous starting model, we calculate them via the adjoint method, based on backpropagation, and parallelize our software on a Linux cluster. Our method is a step forward from ray theory, as it surpasses the inherent infinite-frequency approximation. It differs from analytical Born theory in that it does not involve a far-field approximation, and accounts, in principle, for non-linear effects like multiple scattering and wave front healing. It is much cheaper than the more accurate, fully 3-D numerical solution of the Earth's equations of motion, which has not yet been applied to large-scale tomography. Our tomographic results and trade-off analysis are compatible with those found in the ray- and analytical-Born-theory approaches.     

Modelling of scattering of seismic waves from a corrugated rough sea surface: a comparison of three methods

We compare three numerical methods to model the sea surface interaction in a marine seismic reflection experiment (the frequencies considered are in the band 10–100 Hz): the finite-difference method (FDM), the spectral element method (SEM) and the Kirchhoff method (KM). A plane wave is incident at angles of 0° and 30° with respect to the vertical on a rough Pierson–Moskowitz surface with 2 m significant wave height and the response is synthesized at 6, 10 and 50 m below the average height of the sea surface. All three methods display an excellent agreement for the main reflected arrival. The FDM and SEM also agree very well all through the scattered coda. The KM shows some discrepancies, particularly in terms of amplitudes.     

Comparison of global synthetic seismograms calculated using the spherical 2.5-D finite-difference method with observed long-period waveforms including data from the intra-Antarctic region

We have been developing an accurate and efficient numerical scheme, which uses the finite-difference method (FDM) in spherical coordinates, for the computation of global seismic wave propagation through laterally heterogeneous realistic Earth models. In the field of global seismology, traditional axisymmetric modeling has been used widely as an efficient approach since it can solve the 3-D elastodynamic equation in spherical coordinates on a 2-D cross-section of the Earth, assuming structures to be invariant with respect to the axis through the seismic source. However, it has the severe disadvantages that asymmetric structures about the axis cannot be incorporated and the source mechanisms with arbitrary shear dislocation have not been attempted for a long time. Our scheme is based on the framework of axisymmetric modeling but has been extended to treat asymmetric structures, arbitrary moment-tensor point sources, anelastic attenuation, and the Earth center which is a singularity of wave equations in spherical coordinates. All these types of schemes which solve 3-D wavefields on a 2-D model cross-section are classified as 2.5-D modeling, so we have named our scheme the spherical 2.5-D FDM. In this study, we compare synthetic seismograms calculated using our FDM scheme with three-component observed long-period seismograms including data from stations newly installed in Antarctica in conjunction with the International Polar Year (IPY) 2007–2008. Seismic data from inland Antarctica are expected to reveal images of the Earth's deep interior with enhanced resolution because of the high signal-to-noise ratio and wide extent of this region, in addition to the rarity of sampling paths along the rotation axis of the Earth. We calculate synthetic seismograms through the preliminary reference earth model (PREM) including attenuation using a moment-tensor point source for the November 9, 2009 Fiji earthquake. Our results show quite good agreement between synthetic and observed seismograms, which indicates the accuracy of observations in the Antarctica, as well as the feasibility of the spherical 2.5-D modeling scheme.     

A data-mapping approach to 3-D wavefield redatuming

In case of a complex overburden, the seismic data can be greatly improved by applying a full wavefield redatuming procedure. In practice, the application of the redatuming process to 3-D data acquired by conventional acquisition designs is non-trivial. Because of the large amount of data involved in the 3-D redatuming process and because of the sparseness of these data, it is impossible to apply conventional wave equation datuming directly.
We present a data mapping approach to redatuming (DMR), which follows the concept of Kirchhoff data mapping. A simplified background medium where no ray bending occurs is assumed for the medium below the datum in order to map an input data set referenced to the acquisition surface to an output data set referenced to the new datum level. The DMR method can be interpreted as a simplified version of the Kirchhoff summation redatuming (KSR) method, where one of the 2-D integrals over the acquisition coordinates can be solved analytically. Consequently, in this approach fewer traces are involved in the computation of one time sample (a 2-D integral is computed instead of a 4-D integral), which makes it particularly attractive for the application to 3-D data sets.
In this paper the theory underlying data mapping redatuming is discussed and the proposed approach is tested on fully sampled 2-D and 3-D synthetic data from models with both simple and complex velocity distributions in the subsurface.
The tests clearly show that the objective of producing results that are comparable to the conventional KSR has been achieved. The redatumed traces are dynamically and kinematically correct. Furthermore, these results confirm that the dependency of the new approach on the assumed medium below the datum level is, indeed, weak because the assumption of a velocity medium where no ray bending occurs is already sufficient to produce correct results.     

Inversion of a normally incident reflection seismogram by the Gopinath–Sondhi integral equation

Summary. The one-dimensional acoustic wave equation has been transformed to two coupled first-order equations whose inverse solution is obtained through application of the Gopinath and Sondhi integral equation. A scattering solution of the Schrödinger wave equation for an explosive source leads us to express the kernel of the Gopinath–Sondhi integral equation in terms of a seismic reflection response. A numerical solution of the integral equation obtained by a trapezoidal rule yields a continuous impedance profile whose derivative has step-like discontinuities. The method is illustrated with computer model studies.     

Finite difference seismograms for laterally varying marine models

Summary. The method of finite differences is applied to the elastic wave equation to generate synthetic seismograms for laterally varying seafloor structures. The results are compared with borehole seismic data from the Gulf of California (Deep Sea Drilling Project Site 485) in which lines were shot over flat and rough topography. The significant new phenomenon observed in both the synthetic seismograms and the field data is the generation of a 'double head wave' due to the interaction of the incident wavefront with the side of a hill and the flat seafoor adjacent to the hill.
In these models the hills are on the order of a seismic wavelength in height and steep velocity gradients occur over distances comparable to wavelengths. Ray theoretical methods would not be suitable for studying such structures. True amplitude record sections are obtained by the finite difference method, which show for these models that the head wave generated at the flat seafloor adjacent to the hill is lower in amplitude than if the hill were not present and is lower in amplitude than the head wave generated at the hill.
A second feature which is important for borehole receivers is the existence of the 'direct wave root' in the upper basement. This energy occurs below the sharp interface when the direct wave impinges on the interface from above. There is no corresponding Snell's law ray path for this energy and the energy is evanescent with depth in the lower medium.
The properties of both the double head wave and the direct wave root are clearly demonstrated in the finite difference 'snapshot' displays.     

Modelling of ground motion in the Higashinada (Kobe) area for an aftershock of the 1995 January 17 Hyogo-ken Nanbu, Japan, earthquake

We model the ground motion from an aftershock of the 1995 January 17 Hyogo-ken Nanbu (Kobe) earthquake to investigate basin edge effects on wave propagation in Higashinada ward, downtown Kobe. Point-source finite-difference seismograms calculated using a double-couple solution and a 2-D basin structure are compared with the ground-motion velocity seismograms recorded in a small station array deployed at sites within and outside the heavily damage zone in Higashinada ward. The comparison suggests that in the frequency range 0.1-2 Hz that was analysed, the observed spatial amplitude variation of the aftershock ground motion is attributable mainly to the basin edge effects. We found that the basin edge effect, caused by the superposition of the direct S wave and the basin-edge-diffracted waves, amplified the ground motion in a narrow zone that is offset by about 0.7 km from the basin edge.     

The construction of effective methods for electromagnetic modelling

Summary. This paper deals with the further development of finite-difference methods for electromagnetic field modelling in two-and three-dimensional cases. The main feature of the approach suggested here is the application of generalized asymptotic boundary conditions valid with the accuracy (1/ρ^N), where ρ is the distance from the heterogeneities. The finite-difference approximation of problems under solution is made using the balance method, which results in 5-point difference schemes in the 2-D case and 7-point difference schemes in the 3-D case. To solve the linear system of difference equations the successive over-relaxation (SOR) method is used, the relaxation factor being chosen during the iteration procedure. In view of the vectorial character of the problem for the 3-D case, a successive blocked over-relaxation method (SBOR) is applied.
The model's validity is based on the comparison of the fields accounted at the ground surface with those computed by the integral transformation of excessive currents, determined in the heterogeneity region using the finite-difference scheme.     

Inverse kinematical problems of reflection seismology—II. Problems of downward field continuation

Summary. This paper discusses the inverse (downward) continuation of a wavefield into a medium, as a means of estimating seismic velocities and the imaging of the geometry of reflecting objects. The notion of a kinematically equivalent (K-equivalent) operator of the field continuation is introduced, which allows both the possibility of replacing the initial Lamé equation by a simpler one (for example, a wave equation) and also the existence of an infinite set of continuation operators which permit us to construct a wavefield with proper characteristics. Any of the K-equivalent operators can be used in the task of imaging of reflectors given the known velocities. If velocities are unknown and reflectors are irregular, it is possible to improve the regularity of reflections by continuing the field downward to a fictitious surface in the media with a velocity V_o, that differs from the true velocity V . The conditions needed to untie loops in the true patterns in the continued field are also given.
Another approach is connected with calculations of dynamic sections using different values of velocity (say, an extremal velocity) under which the amplitudes of reflectors on sections are brought to a maximum due to the phenomena of focusing. Extremal velocities can be also used for solving inverse kinematical problems.     

Numerical modelling of propagation of elastic waves in anisotropic inhomogeneous media for the half-space and the sphere

Summary. A method based on a combination of partial separation of variables and finite-difference method is used for the calculation of complete theoretical seismograms for inhomogeneous anisotropic media. Examples of theoretical seismograms for several anisotropic models are presented.     

Numerical scheme for the inversion of acoustical impedance profile based on the Gelfand-Levitan method

Summary. An exact method for the solution of the inverse problem in plane wave propagation modelled after the Gelfand-Levitan technique is reviewed and refined. A numerical scheme for the solution of the integral equation that arises in the method is proposed. A discussion on the stability and an error analysis of the numerical approximation are presented. The applicability of the inversion algorithm is demonstrated in a numerical experiment.     

Fast transform eigenfunction methods for boundary value problems

An extension of the classic eigenfunction expansion method is presented for solving boundary value problems on boundaries of essentially arbitrary shape in R ³. the algorithm yields approximate solutions in terms of simple eigenfunctions of the underlying differential operator. the coefficients of these eigenfunctions are calculated with two-dimensional FFTs. the algorithm has many applications in geophysics, such as potential theory, scattering theory, earthquake source dynamics, and wave propagation, where non-separable boundaries are common.     

Regional tomographic inversion of the amplitude and phase of Rayleigh waves with 2-D sensitivity kernels

In this study, we test the adequacy of 2-D sensitivity kernels for fundamental-mode Rayleigh waves based on the single-scattering (Born) approximation to account for the effects of heterogeneous structure on the wavefield in a regional surface wave study. The calculated phase and amplitude data using the 2-D sensitivity kernels are compared to phase and amplitude data obtained from seismic waveforms synthesized by the pseudo-spectral method for plane Rayleigh waves propagating through heterogeneous structure. We find that the kernels can accurately predict the perturbation of the wavefield even when the size of anomaly is larger than one wavelength. The only exception is a systematic bias in the amplitude within the anomaly itself due to a site response.
An inversion method of surface wave tomography based on the sensitivity kernels is developed and applied to synthesized data obtained from a numerical simulation modelling Rayleigh wave propagation over checkerboard structure. By comparing recovered images to input structure, we illustrate that the method can almost completely recover anomalies within an array of stations when the size of the anomalies is larger than or close to one wavelength of the surface waves. Surface wave amplitude contains important information about Earth structure and should be inverted together with phase data in surface wave tomography.     

A bottom‐up approach to discover transition rules of cellular automata using ant intelligence

This paper presents a new method to discover transition rules of geographical cellular automata (CA) based on a bottom‐up approach, ant colony optimization (ACO). CA are capable of simulating the evolution of complex geographical phenomena. The core of a CA model is how to define transition rules so that realistic patterns can be simulated using empirical data. Transition rules are often defined by using mathematical equations, which do not provide easily understandable explicit forms. Furthermore, it is very difficult, if not impossible, to specify equation‐based transition rules for reflecting complex geographical processes. This paper presents a method of using ant intelligence to discover explicit transition rules of urban CA to overcome these limitations. This ‘bottom‐up’ ACO approach for achieving complex task through cooperation and interaction of ants is effective for capturing complex relationships between spatial variables and urban dynamics. A discretization technique is proposed to deal with continuous spatial variables for discovering transition rules hidden in large datasets. The ACO–CA model has been used to simulate rural–urban land conversions in Guangzhou, Guangdong, China. Preliminary results suggest that this ACO–CA method can have a better performance than the decision‐tree CA method.