Wavefront Construction Research Articles

Analytical and numerical construction of the optimal outcome function for a class of time-optimal problems

Analytical and numerical algorithms are proposed for constructing the optimal outcome function and its Lebesgue set for the time-optimal control problem with a circular velocity indicatrix. Our approach to the solution of the time-optimal control problem essentially utilizes the specific dynamic properties of the controlled system. The circular vectogram of possible velocities enables us to interpret the cross sections of the controllability set as wavefronts whose source is uniformly distributed on the boundary of the target set. Procedures have been developed for analytical and numerical construction of the evolution of wavefronts based on prior (given the geometry of the target set boundary) identification of their nonsmoothness sets. An essential feature of the construction is the point-to-set distance function. We investigate the differential properties of this function and identify the manifolds on which it loses its smoothness. The proposed wavefront construction algorithms are of independent interest in so far as they enable us to investigate the geometry of the sets and compute their nonconvexity measure. The results are useful not only when studying the evolution of reachability sets of controlled systems, but also for computing the eikonal in geometrical optics and investigating the solutions of the wave equation.

A practical implementation of wave front construction for 3-D isotropic media

SUMMARY Wave front construction (WFC) methods are a useful tool for tracking wave fronts and are a natural extension to standard ray shooting methods. Here we describe and implement a simple WFC method that is used to interpolate wavefield properties throughout a 3-D heterogeneous medium. Our approach differs from previous 3-D WFC procedures primarily in the use of a ray interpolation scheme, based on approximating the wave front as a ‘locally spherical’ surface and a ‘first arrival mode’, which reduces computation times, where only first arrivals are required. Both of these features have previously been included in 2-D WFC algorithms; however, until now they have not been extended to 3-D systems. The wave front interpolation scheme allows for rays to be traced from a nearly arbitrary distribution of take-off angles, and the calculation of derivatives with respect to take-off angles is not required for wave front interpolation. However, in regions of steep velocity gradient, the locally spherical approximation is not valid, and it is necessary to backpropagate rays to a sufficiently homogenous region before interpolation of the new ray. Our WFC technique is illustrated using a realistic velocity model, based on a North Sea oil reservoir. We examine wavefield quantities such as traveltimes, ray angles, source take-off angles and geometrical spreading factors, all of which are interpolated on to a regular grid. We compare geometrical spreading factors calculated using two methods: using the ray Jacobian and by taking the ratio of a triangular area of wave front to the corresponding solid angle at the source. The results show that care must be taken when using ray Jacobians to calculate geometrical spreading factors, as the poles of the source coordinate system produce unreliable values, which can be spread over a large area, as only a few initial rays are traced in WFC. We also show that the use of the first arrival mode can reduce computation time by ∼65 per cent, with the accuracy of the interpolated traveltimes, ray angles and source take-off angles largely unchanged. However, the first arrival mode does lead to inaccuracies in interpolated angles near caustic surfaces, as well as small variations in geometrical spreading factors for ray tubes that have passed through caustic surfaces.

Efficient level set methods for constructing wavefronts in three spatial dimensions

Wavefront construction in geometrical optics has long faced the twin difficulties of dealing with multi-valued forms and resolution of wavefront surfaces. A recent change in viewpoint, however, has demonstrated that working in phase space on bicharacteristic strips using eulerian methods can bypass both difficulties. The level set method for interface dynamics makes a suitable choice for the eulerian method. Unfortunately, in three-dimensional space, the setting of interest for most practical applications, the advantages of this method are largely offset by a new problem: the high dimension of phase space. In this work, we present new types of level set algorithms that remove this obstacle and demonstrate their abilities to accurately construct wavefronts under high resolution. These results propel the level set method forward significantly as a competitive approach in geometrical optics under realistic conditions.

An improved mesh generation scheme for the wavefront construction method

Wavefront construction is an effective tool for the rapid calculation of ray fields in anisotropic media. The method explicitly tracks the propagation of a wavefront through a model, mapping it to a computational mesh that is interpolated when accuracy criteria based on paraxial ray methods are violated. Takeoff angles are used often to define the initial ray directions, but uniform sampling in the two angles leads to oversampling of the ray field in the direction of the axis. Such sampling can lead also to numerical instability associated with vanishing derivatives with respect to the azimuthal angle. We suggest a new wavefront mesh definition using the cubed-sphere mesh, which is a coordinate system used to solve partial differential equations in spherical geometries. When using this mesh, ray directions are assigned by mapping points on a regular discretization of the faces of a cube surrounding the source to corresponding rays. This scheme produces a nearly uniform distribution of rays with minimal effort and using the cubed-sphere coordinates as ray parameters to calculate partial derivatives completely eliminates the singularities that arise when takeoff angles are used as ray parameters. Numerical results for quantities related to seismic amplitudes confirm that this new mesh does provide more stable and reliable results.

Amplitude, Fresnel zone, and NMO velocity for PP and SS normal-incidence reflections

Inspired by recent ray-theoretical developments, the theory of normal-incidence rays is generalized to accommodate P- and S-waves in layered isotropic and anisotropic media. The calculation of the three main factors contributing to the two-way amplitude — i.e., geometric spreading, phase shift from caustics, and accumulated reflection/transmission coefficients — is formulated as a recursive process in the upward direction of the normal-incidence rays. This step-by-step approach makes it possible to implement zero-offset amplitude modeling as an efficient one-way wavefront construction process. For the purpose of upward dynamic ray tracing, the one-way eigensolution matrix is introduced, having as minors the paraxial ray-tracing matrices for the wavefronts of two hypothetical waves, referred to by Hubral as the normal-incidence point (NIP) wave and the normal wave. Dynamic ray tracing expressed in terms of the one-way eigensolution matrix has two advantages: The formulas for geometric spreading, phase shift from caustics, and Fresnel zone matrix become particularly simple, and the amplitude and Fresnel zone matrix can be calculated without explicit knowledge of the interface curvatures at the point of normal-incidence reflection.

Modeling and velocity analysis with a wavefront-construction algorithm for anisotropic media

Wavefront construction methods allow for efficient simulation of seismic wave propagation in anisotropic media using ray theory. This approach models an individual wavefront as a mesh that is adaptively refined as the wave travels through an earth model. The algorithm, implemented here for quasi-compressional waves, provides an effective tool for modeling and analyzing a large vertical seismic-profile data set from the Gulf of Mexico. Seismic velocities in this region generally increase linearly with depth, but simulations show that an isotropic vertical-gradient model cannot accurately reproduce traveltimes of first arrivals. However, a simple velocity-analysis inversion based on the wavefront-construction method allows a straightforward estimation of improved isotropic and transversely isotropic models that significantly reduce the misfit between computed traveltimes and the observations. The goal of this analysis is to produce a simple model suitable for use as a starting point for migrations and for assessing the importance of anisotropy in the field site. In fact, the final anisotropic model is much more consistent with background geologic information than the best isotropic model, which also has a larger rms traveltime error. The results strongly support the presence of measurable anisotropy.

Traveltime computation by wavefront‐orientated ray tracing

ABSTRACTFor multivalued traveltime computation on dense grids, we propose a wavefront‐orientated ray‐tracing (WRT) technique. At the source, we start with a few rays which are propagated stepwise through a smooth two‐dimensional (2D) velocity model. The ray field is examined at wavefronts and a new ray might be inserted between two adjacent rays if one of the following criteria is satisfied: (1) the distance between the two rays is larger than a predefined threshold; (2) the difference in wavefront curvature between the rays is larger than a predefined threshold; (3) the adjacent rays intersect. The last two criteria may lead to oversampling by rays in caustic regions. To avoid this oversampling, we do not insert a ray if the distance between adjacent rays is smaller than a predefined threshold. We insert the new ray by tracing it from the source. This approach leads to an improved accuracy compared with the insertion of a new ray by interpolation, which is the method usually applied in wavefront construction. The traveltimes computed along the rays are used for the estimation of traveltimes on a rectangular grid. This estimation is carried out within a region bounded by adjacent wavefronts and rays. As for the insertion criterion, we consider the wavefront curvature and extrapolate the traveltimes, up to the second order, from the intersection points between rays and wavefronts to a gridpoint. The extrapolated values are weighted with respect to the distances to wavefronts and rays. Because dynamic ray tracing is not applied, we approximate the wavefront curvature at a given point using the slowness vector at this point and an adjacent point on the same wavefront. The efficiency of the WRT technique is strongly dependent on the input parameters which control the wavefront and ray densities. On the basis of traveltimes computed in a smoothed Marmousi model, we analyse these dependences and suggest some rules for a correct choice of input parameters. With suitable input parameters, the WRT technique allows an accurate traveltime computation using a small number of rays and wavefronts.

Fast migration/inversion with multivalued ray fields: Part 1—Method, validation test, and application in 2D to Marmousi

We address the problem of building an efficient algorithm for ray‐based migration/inversion in complex media characterized by multivalued traveltime. Since the multivalued ray maps can be computed by wavefront construction, the key point lies in the interpolation of these multivalued coarse maps. For this goal we propose two new strategies. The first one is based on an interpolation of the ray field after identifying the ray branches on the elementary cells of the coarse grid; the second is based on an extrapolation of the arrivals from the coarse grid points. Both strategies are evaluated in two dimensions with the Marmousi data set, demonstrating the CPU efficiency and accuracy of both approaches. Since they both consider multiarrivals locally, the associated overhead cost with respect to a single‐arrival implementation is quite low (twice the computing time of the Marmousi test case). The application also reveals that the new interpolation strategy is more efficient in terms of CPU time compared to the extrapolation strategy but at the expense of a more complex numerical implementation and less accurate results.

Fast migration/inversion with multivalued ray fields: Part 2—Applications to the 3D SEG/EAGE salt model

Three‐dimensional prestack depth migration is the convenient approach for seismic imaging in the case of strong lateral variations of the velocity. Because of computing limitations, it has been limited to single‐arrival kinematic Kirchhoff migration until recently. This approach fails in the case of complex media characterized by multiarrival traveltimes. We present numerical strategies for extending in three dimensions first‐arrival kinematic Kirchhoff migration to multiarrival quantitative ray‐based migration (preserved amplitude migration). We rely on wavefront construction in a smooth velocity model to compute the multivalued traveltime and amplitude maps, and the CPU efficiency of migration itself is ensured by efficient and robust interpolation or extrapolation strategies. We present an application to the synthetic 3D SEG/EAGE salt model. Taking into account multiarrivals clearly improves subsalt imaging at the price of quite limited computing costs (a 20% increase in our case, with respect to a preserved‐amplitude single‐arrival migration).

Efficient and accurate computation of seismic traveltimes and amplitudes

ABSTRACTWe 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.

Automatic velocity analysis by differential semblance optimization

We present a method for automatic velocity analysis of seismic data based on differential semblance optimization (DSO). The data are mapped for each offset from the time domain to the depth domain by a Born migration scheme using ray tracing with the efficient wavefront construction method. The DSO cost functional is evaluated by taking differences of the migration images for neighboring offsets. The gradient of this functional with respect to the underlying velocity model is obtained by a first‐order approximation of the adjoint‐state method, leading to an optimal complexity: the cost of evaluating the gradient is about the same as that of evaluating the functional. The method has been applied to a marine line. Multiples turned out to be a problem, but were handled effectively by incorporating a multiple filter inside the DSO cost functional.

Nonlinear Hyperbolic Waves in Multi-Dimensions. Monographs and Surveys in Pure and Applied Mathematics, Vol. 121.

7R10. Nonlinear Hyperbolic Waves in Multi-Dimensions.Monographs and Surveys in Pure and Applied Mathematics, Vol 121.- P Prasad (Dept of Math, Indian Inst of Sci, Bangalore, India). Chapman and Hall/CRC, Boca Raton FL. 2001. 338 pp. ISBN 1-58488-072-4. $89.95.Reviewed by Shi Tsan Wu (Dept of Mech and Aerospace Eng, Univ of Alabama, Sparkman Dr, Huntsville AL 35899).In this book, the authors have discussed the propagation of a curved nonlinear wavefront for the ordinary fluid dynamics (ie without considering magneto-fluids) according to two fundamental physical processes: i) at different points of the wavefront, it travels with different speeds depending on the local amplitude leading to a longitudinal stretching of rays; and ii) a lateral deviation of rays is produced due to non-uniform distribution of the amplitude of the wave front. This book consists of ten chapters. The first three chapters introduce all necessary mathematical concepts and some basic results. For example, the fundamentals of nonlinear hyperbolic waves are given in Chapter 1. In Chapter 2, the basic results of hyperbolic systems are presented. These results include i) discussion of hyperbolic system of first order equations in two independent variables with a canonical form of a system of linear and semi-linear equations, and ii) the hyperbolic system in more than two variables together with propagation of discontinuities of first order derivatives along rays is given. In Chapter 3, the results for simple wave high frequency approximation and ray theory are given. Specifically, the Huygen’s method of wavefront construction, Huygen’s method and ray theory, and Fermat’s principle are discussed in detail. The next seven chapters deal with specifics of the wave propagation of nonlinear hyperbolic systems; the derivation of weakly nonlinear ray theory (WNLRT) is given in Chapter 4. The stability of solutions for one-dimensional weakly nonlinear wave propagation as well as waves in a multi-dimensional steady transonic flow near a sonic singularity are discussed in Chapter 5. In Chapter 6, a special case of WNLRT for polytropic gas was presented in detail for two spatial dimensions. In Chapter 7, the authors presented a single conservation law approach to derive the compatibility conditions on a shock which is claimed as a New Theory of Shock Dynamics (NTSC). To establish this NTSC, the authors present a proof of the existence and uniqueness of the solution and made comparison of the numerical results deduced from NTSD with the exact solution which shows good agreement. Then, this NTSD is applied to a one-dimensional piston problem to illustrate its capabilities in Chapter 8. In Chapter 9, the authors presented a derivation of the compatibility conditions on a shock in multi-dimensions using shock ray theory. The final Chapter 10 includes the application of NTSD to the propagation of a curved weak shock. Overall, the book has three parts: i) The first three chapters have given all the necessary fundamental mathematical concepts, models and methods to calculate complete evolutionary history of a curved nonlinear wave front consisting of shock front. ii) The detailed discussions of WNLRT and NTSD are presented in Chapters 4-8; and iii) the specific applications of shock ray theory and NTSD are given in Chapters 9 and 10, respectively. In summary, Nonlinear Hyperbolic Waves in Multi-Dimensions includes numerical examples to illustrate theory, the quality of the graphic is good and an adequate subject index is provided. The author has achieved his goal as a textbook for advanced undergraduate applied mathematics. However, it could be used as a textbook for first year graduate students in engineering and physics if the author had presented some problems at the end of each chapter. This should be able to be remedied by the instructor. This reviewer strongly recommends that scientists and engineers who are interested in these topics read this book; it would be a useful reference book.

Geometric Optics in a Phase-Space-Based Level Set and Eulerian Framework

We introduce a level set approach for ray tracing and the construction of wavefronts in geometric optics. This is important in a wide variety of applications in wave propagation. Our approach automatically handles the multivalued solutions that appear and automatically resolves the wavefronts. This is achieved through solving for the bicharacteristic strips, whose projection to spatial space gives the wavefronts, in a reduced phase space under a Eulerian and partial differential-equation-based framework. The bicharacteristic strips are represented using a level set approach for handling higher codimensional objects and the partial differential equations responsible for the evolution are reduced forms of the Liouville equations. Results for the two-dimensional case for constant and variable indices of refraction are shown and compared to those of other current methods in the field. Results are also introduced to show the ability to handle reflection and to extend the method to the three-dimensional case.

Review of Ray Theory Applications in Modelling and Imaging of Seismic Data

Throughout the last twenty years, 3D seismic ray modelling has developed from a research tool to a more operational tool that has gained growing interest in the petroleum industry. Various areas of application have been established and new ones are under development. Many of these applications require a modelling system with flexible, robust and efficient modelling algorithms in the core. The present paper reviews the basic elements of such a system, based on the ‘open model’ concept and the ‘wavefront construction’ technique. In the latter, Cervený's dynamic ray tracing is an intrinsic part. The modelling system can be used for generating ray attributes and synthetic seismograms for realistic 3D surveys with tens of thousands of shots and receivers. Moreover, some other types of application areas are illustrated: Production of Green's functions for prestack depth migration and hybrid modelling (combined ray and finite-difference modelling), attribute mapping and illumination analysis, both for survey planning and interpretation. Finally, the concepts of ‘isochron rays’ and ‘velocity rays’ related to seismic isochrons have been introduced recently, with very interesting future applications.

Can we quantitatively image complex structures with rays?

Ray‐based prestack depth migration fails to image quantitatively complex structures when a single arrival—for example, the first or the strongest one—is taken into account. In this paper, we investigate whether accounting for multiple arrivals in ray‐based preserved amplitude prestack depth migration allows one to improve quantitative imaging of complex media. The asymptotic ray‐Born migration/inversion, originally designed to process one single arrival, is extended to the case of multiple arrivals by accounting for the cross‐contributions of all the source and receiver raypaths. Multiple arrivals in the folded ray fields are computed by a dynamic ray tracing based on a wavefront construction technique. With an application to the complex Marmousi model, we demonstrate that ray‐Born inversion can provide a reliable quantitative estimation of the relative impedance perturbation even in the complex deep part of the model, for which the amplitudes were underestimated drastically when only a single arrival was used, either the first or the strongest one. This 2-D case study shows that complex structures can be imaged quantitatively with rays. Future studies will require the optimizing of the implementation of ray‐Born migration/inversion with multiple arrivals before considering 3-D applications.

3‐D Multivalued Seismic wavefield Reconstruction and Computation for Green's Function

AbstractBased on wavefront construction, an efficient algorithm for multivalued traveltime wavefield is presented. This method overcomes the difficulties with ray series method in conventional configuration space due to complex velocity structure by implementing the calculation in phase space instead. It solves out the required multi‐valued traveltime and amplitude at the uniform dispersed grids efficiently when applied to 3‐D prestack depth migration. Two criterions for controlling and balancing the global computational accuracy and efficiency are also discussed in this paper. The numerical examples show that our method is applicable to complex geological models.

3-D ray modeling by wavefront construction in open models

A synthesis of two newly developed concepts in 3-D modeling is developed in this paper: (1) The open, (noncomplete) seismic model, and (2) the ray tracing based wavefront (WF) construction method. The open model may contain interfaces with holes and other missing parts, which simplifies model building considerably because the input horizon data from standard interpretation and processing packages are often incomplete. A set of volumes is defined in the model. A volume is a logical unit that points to a set of property functions, e.g., P-velocity, S-velocity, and density. The properties are represented either as constants or as B-spline functions of the spatial coordinates (x, y) or (x, y, z). The volumes of the model are assigned to opposing sides of each interface and not to specific spatial areas of the model, which is the case in most (blocky) model representations. The interfaces are given explicitly by triangular grids where the sizes of the triangles are determined locally by the curvature of the interface. We show how modeling by WF construction is both possible and computationally efficient in open models, but only after some modifications to deal with the ambiguity of the model representation. It is not possible to find a unique volume for the spatial positions in an open model. Instead, the volumes (with associated velocities, etc.) are determined from the last interface encountered by each ray in the WF. To find an arrival in a receiver, the volume associated to the receiver has to match the volume of the WF hitting the receiver.

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.

Comparison of six different methods for calculating traveltimes

We study six different methods for the calculation of seismic traveltimes. All methods yield traveltimes at all points of a regular grid.The methods examined comprise three different variants of finite‐difference (FD) eikonal solvers, the graph method, wavefront construction and a combined FD and Runge–Kutta method.The main points of investigation are computational time, accuracy and memory requirements. We took measures to obtain a high level of both general validity and clear understanding of the results. We used a profiling program to be able to measure the time that the actual core algorithm needs, thus avoiding any overhead of highly system‐dependent in‐/output operations.The comparison shows that no single method is the most appropriate but that the choice depends on the task to be fulfilled. The FD eikonal solver that uses expanding squares proves to be best suited for models which are not too complicated because it offers the best compromise between speed and accuracy, whereas wavefront construction should be applied to complex media because of its superior reliability which then justifies the much higher computational times.

3-D preserved amplitude prestack depth migration on a workstation

We present an algorithm based on the ray+Born approximation for 3-D preserved amplitude prestack depth migration (PAPsDM) of seismic reflection data. This ray+Born inversion scheme allows the quantitative recovery of model perturbations. The Green’s functions are estimated by dynamic ray tracing in 3-D heterogeneous smooth velocity fields with a wavefront construction (WFC) method. The PAPsDM algorithm was implemented on a single‐processor Sun Sparc 20 workstation. Special attention was paid to CPU efficiency and memory requirements. We present an application on a 3-D real marine data set (13 Gbytes). About one week of CPU time is needed to obtain a migrated image of 7 × 1 × 1 km.