Abstracts

3D raytracing through homogeneous anisotropic media with smooth interfaces

Jessé Costa^I; Michael Schoenberg; Jaime Urban^II

^IDepartamento de Geofísica. Universidade Federal do Pará (UFPA), Rua Augusto Correa, nº 1, Bairro do Guamá, CEP: 66070-110, Belém-PA. Fax: 91-3184-1693, Tel: 91-3184-1692. E-mail: jesse@ufpa.br

^IIDepartamento de Física. Universidade Federal do Pará (UFPA), Rua Augusto Correa, nº 1, Bairro do Guamá, CEP: 66070-110, Belém-PA. Fax: 91-3184-1693, Tel: 91-3183-1403. E-mail: urban@ufpa.br

ABSTRACT

Two-point raytracing problem is solved for events in a piecewise homogeneous and laterally varying 3D anisotropic media by continuation techniques. In conjunction with the shooting method the algorithm can be used for computation of qP, qS1, and qS2 events. The algorithm has the same performance and robustness as previous implementations of the continuation method for tracing rays in isotropic models. Routines based on our algorithm have several useful applications. First, an efficient forward problem solver for traveltime inversion of elastic parameters in the presence of anisotropy. Second, Newton-Raphson iterations during two-point raytracing produce wavefront attributes, slowness and wavefront curvature. These attributes allows the computation of geometrical spreading and second order approximations for traveltimes. Therefore it can be used to investigate the effects of anisotropy on CRS, in simple velocity models.

Keywords: Ray theory, continuation method, anisotropy

RESUMO

O raio conectando dois pontos em um meio anisotrópico, homogêneo por partes e com variação lateral, é calculado utilizando-se técnicas de continuação em 3D. Se combinado com algoritmos para solução do problema de valor inicial, o método pode ser estendido para o cálculo de eventos qS₁ e qS₂. O algoritmo apresenta a mesma eficiência e robustez que implementações de técnicas de continuação em meios isotrópicos. Rotinas baseadas neste algoritmo têm várias aplicações de interesse. Primeiramente, na modelagem e inversão de parâmetros elásticos na presença de anisotropia. Em segundo lugar, as iterações de Newton-Raphson produzem atributos da frente de onda como vetor vagarosidade e a matrix hessiana do tempo de trânsito, quantidades que permitem determinar o espalhamento geométrico e aproximações de segunda ordem para o tempo de trânsito. Estes atributos permitem calcular as amplitudes ao longo do raio e investigar os efeitos da anisotropia no empilhamento CRS em modelos de velocidade simples.

Palavras-chave: Teoria do Raio, Método de Continuação, Anisotropia

INTRODUCTION

A computational scheme is proposed to solve two-point raytracing, in a piecewise homogeneous and laterally varying 3D anisotropic layered media, based on the continuation method . The algorithm is an extension to the anisotropic case of previous works for isotropic media (KELLER; PEROZZI, 1983; DOCHERTY; BLEISTEIN, 1984).

The continuation method permits to trace rays through a complex model by continuous transformations from simpler intermediary configurations. This procedure warrant robustness to the iterative techniques used to solve the nonlinear system of equations that arises from Fermat's principle. A sequence of continuation steps is performed. Starting with a vertical ray in an isotropic layered medium, the isotropic slowness surfaces are deformed to the intended anisotropic slowness surfaces at each layer. The next step is to transform the flat interfaces to the desired curved ones. Finally we move source and receiver positions to the configuration we want to achieve.

The algorithm has the same perfomance described by Docherty and Bleistein (1984) for isotropic models when computing qP raypaths. The computation of qS₁ and qS₂ raypaths presents difficulties at the singular regions where the slowness surfaces for these waves intersects one another and regions where they are not convex. A shooting approach can remedy most of these problems, but the performance degrades if many restarts are required.

This procedure is well adapted to compute raypath and traveltimes for surface seismics and VSP experiments in anisotropic models. The algorithm can be a forward problem routine for a traveltime inversion scheme which produces estimatives of the elastic constants structure for subsurface, most like the works of Whitmore and Lines (1986), Chiu and Stewart (1987) and Guiziou and Haas (1988), without the limitations of the elliptical anisotropy assumption. Other possible applications include VSP-CDP mapping and the computation of wavefront attributes,slowness vector and wavefront curvature, along the raypath.

FERMAT'S PRINCIPLE AND RAYPATH COMPUTATION

We assume a piecewise homogeneous, anisotropic, layered medium where the interfaces are described by functions

f(x) = x₃ – z(x₁, x₃) = 0

where x₁ represents Cartesian coordinates. The two-point raytracing problem in layered anisotropic medium through N interfaces can be stated, by Fermat's principle, as the following variational problem in phase space,

subjected to,

where, t is the traveltime along the ray, , represents the ray intersection with the a-th interface, is the source position, is the receiver position, is the slowness vector at the a-th layer, is the dispersion relation at the a-th layer. Summation convention is assumed on subscript index j which identifies the cartesian coordinates directions.

Equation (1) arises from the relations (MUSGRAVE, 1970),

Here s is the slowness , v is the group velocity and t traveltime along the ray. The last identity results from the polar reciprocal relationship between slowness and group velocity s · v º 1. For a piecewise homogeneous layered media equation (3) reduces to (1).

The solution of (1) subject to (2) is straightforward using Lagrange multipliers,

following standard procedures we obtain the nonlinear system of 5N + 3 equations on 5N + 3 unknowns,namely, :

These five equations express the basic geometrical relations for rays and slownesses through a homogeneous layered medium. Equations for and can be written as

where n^a is the unit vector normal to the a-th interface.This equation is a statement of the Snell's law. It requires that the variation of the slowness vector across an interface be parallel to the interface normal, i.e., tangential slowness component is continuous across the interface. Equations for and , express the relation between ray direction and slowness surface, and can be written as

this equation states that ray direction should be normal to the slowness surface. The last equation for F₅, defines the slowness surfaces at each layer. The system of equations (5) in vector notation is,

where,

The system (8) can be solved iteratively using Newton-Raphson method. At the K-th iteration this corresponds to:

1. solve

Du^(K) · Ñ_uF = F(u^(K)) for Du ,

2. update

CONTINUATION METHODS

Continuation techniques are used to solve the system (8) in order to add robustness to the Newton-Raphson iterations (9). The continuation method consists in solving the following nonlinear equation

instead of (8).

The parameter g is intended to control the smooth transformation from problem F(u,0) = 0, whose solution can be readily found, to F(u,1) = F(u), which is the problem to be solved.

The solution of the two-point raytracing in 3D layered anisotropic medium can be constructed imbedding a sequential application of the continuation technique for the solution of (8).

Continuation of slowness surfaces For each layer, the dispersion relation may be written in the form

where and are the dispersion relations for the anisotropic layers and for initial isotropic layers, respectively. The starting layered isotropic model (g = 0), has compressional velocities = max(a₁₁, a₂₂, a₃₃), at each layer, were are the density normalized elastic moduli in condensed notation (MUSGRAVE, 1970).

Begining with a vertical ray in the initial flat layered isotropic model, so u(g = 0) is trivially computed, the algorithm proceeds as follows. The initial guess for Newton iterations (9) for a new g = g + D g is computed (DOCHERTY; BLEISTEIN, 1984):

1. solve

2. start Newton iterations from

Then iterate using Newton-Raphson equations (9). If the Newton method does not converge in 5 iterations, halve Dg and start again. After convergence double the value Dg and continue in the same way until g = 1. Usually, we can go from g = 0 to g = 1 in one step Dg = 1. When g = 1, the model is the desired set of anisotropic layers with flat interfaces.

Continuation of interfaces The next step is to deform the interfaces from flat layers to the desired ones keeping source and receiver vertically aligned. This is done using,

so, when g = 0, the interfaces are flat, and, when g = 1, the desired interfaces are achieved. The interfaces are defined by specifying a set of points z(x₁, x₂) over a grid for each interface. B-spline interpolation are used to define the interfaces and evaluation of the required derivatives.

The procedure follows as before. The starting values for Newton iterations is computeted similarly from

and

where u(g = 0) is the solution of the slowness continuation procedure.

Receiver and Source continuation Now that the intended media structure was achieved, the next step is to move the receiver to the desired configuration applying the continuation procedure to the relation,

where x_rec = (x_rec, y_rec, z_rec). Here g = 0 means a previous computed receiver position and g = 1 the new receiver position for raypath computation.

The starting guess for Newton iterations comes from the solution of

with

here u(g = 0) is the solution for the previous receiver position. The same procedure is applied for source continuation. The corresponding expressions can be obtained just changing the labels of receiver position for the source position at the above equations.

SHOOTING STILL REQUIRED

Two-point raytracing has its drawbacks, though. Firstly, Newton-Raphson iterations fail when matrix Ñ_uF is singular. This occurs in several instances during raytracing: over an interface saddle point, if ray is tangent to an interface, or when receivers are at shadow zones. Secondly, when more than one ray connects source and receiver, two-point raytracing converges to a single trajectory depending on initial guess, in other words, only a single branch of traveltime surface is computed and the method fails at the borders of these branches. To overcome these problems the algorithm has to be combined with a shooting strategy (PRESS et al., 1989). Whenever Newton iterations fail a ray is computed using the shooting algorithm. Since raypath is a straight line in each layer, raytracing is reduced to compute the intersection of the ray with the interfaces and apply Snell's law to proceed across an interface until we reach the layer containing the receivers. Once a new ray is successfully computed in the neighborhood of a singular point, two-point raytracing can be tried again to compute rays to nearby receivers. If the ray identifies a shadow zone the procedure is interrupted.

NUMERICAL RESULTS

The algorithm was applied to compute raypaths for a model consisting of five layers, having strongly anisotropic elastic properties, separated by smooth interfaces. The density normalized elastic parameters for each layer are presented in ^{appendix A} APPENDIX A . The layers are indexed from 1 to 5 from top to bottom. The qP group velocity surface for each medium are shown in Figure 1, Figure 2 and Figure 3.

We only present computation of rays associated with qP events, although qS events can be computed if an initial ray is determined using the shooting method. Figures 4(a) and 4(b) show the result of raypath computation for receivers along a line in the surface. The raypaths include a multiple in the third layer.

The second example, Figure 5, includes an unusual configuration of source and receivers at the surface just to point out the possibilities of this implementation. Again a multiple was computed in the third layer. The third example, Figure 6, has the same distribution of receivers at the surface but now the source is located in the third layer as in a reverse 3-D VSP experiment.

One can always find complex models, which present shadow zones or rough interfaces, where Newton iterations fail and many restarts using the shooting method are required. This algorithm is not well suited for these models. For applications which do not require such complex models and where qP waves are the main concern the algorithm is quite efficient.

CONCLUSION

An extension of the continuation method proposed by Keller & Perozzi (1983) to anisotropic models in 3-D was developed. The algorithm is efficient when computing qP events. The computation of the qS₁ and qS₂ trajectories requires the use of shooting method to find, besides the starting ray, new initial approximations whenever two-point raytracing fails due to singular regions in the associated slowness surface or other possible causes. Although the algorithm do not handle general inhomogeneous anisotropic models, it allows a simple model specification for interfaces and arbitrary anisotropy. The performance of the algorithm permits its application to procedures requiring the computation of a large number of raypaths, as on inversion algorithms and VSP-CDP mapping in 3-D. It can also be helpful, when wavefront attributes as slowness vector and wavefront curvature are required as occurs in geometrical spreading computations or in CRS studies.

Acknowledgements

This work was kindly supported by the sponsors of the Wave Inversion Technology (WIT) Consortium, Karlsruhe, Germany.

Recebido em 16 set., 2002 / Aceito em 11 dez., 2003

Received Set. 16, 2002 / Accepted Dec. 11, 2003

NOTES ABOUT THE AUTHORS

Jessé Carvalho Costa é formado em Física (UFPA/1983), Mestre e Doutor em Geofísica (UFPA/1987, 1993 respectivamente). Summer Student na Schlumberger Cambridge Research em 1991 e 1992. Estágio de pós-doutoramento no Departamento de Geofísica na Universidade de Stanford (1994-1996) e Visiting Assistent Professor no departamento de Geofísica da Universidade de Stanford (1995). É Professor da Universidade Federal do Pará desde 1989, no Departamento de Física de 1989/2003, atualmente é Professor do Departamento de Geofísica desta mesma Universidade. Áreas de interesse: anisotropia, modelagem sísmica e tomografia.

Michael Schoenberg se aposentou em 1999, após 21 anos de carreira como pesquisador para a Schlumberger em Ridgefield, CT, Cambridge, UK, e Tókio. Desde então, tem ministrado cursos, coordenado seminários e prestado consultoria em várias instituições pelo mundo do Lawrence Berkeley Lab, na Califórniam ao CSIRO em Perth, com estadas breves em Deft, na Holanda e Edimburgo, na Escócia. Seus interesses se concentram em propagação de ondas elásticas e processamento de dados sísmicos em meios anisotrópicos, teoria de meio efetivo for maios estratificados e fraturados, e física de rochas. As aplicações principais de sua pesquisa estão na interpretação de dados de monitoramento sísmicos com o tempo de reservatórios para produção ou injeção de fluidos.

Jaime A. Urban é bacharel em Física (1995) e mestre em Geofísica (1999) pela UFPA. Em 1997 Jaime foi professor substituto do departamento de Matemática da UFPA e, desde 1998, é professor Assistente de Física da UFPA. Atualmente Jaime é aluno do programa de doutorado em Geofísica da Stanford University (USA). Suas principais áreas de pesquisa são imageamento e monitoramento sísmico de reservatórios de petroleo e monitoramento da produção de gás em reservatórios de carvão mineral (coalbed methane). Jaime é membro da SBGf e da SEG.

The density normalized elastic tensor for the three different materials forming the layers in the model used for numerical test. The units are in (km/s)².

APPENDIX A 

Publication Dates

History