WO2012087176A2

WO2012087176A2 - Method for computing wavefields

Info

Publication number: WO2012087176A2
Application number: PCT/RU2010/000770
Authority: WO
Inventors: Vladimir Albertovich Cheverda; Galina Vitalyevna RESHETOVA; Mikhail Andreevich BELONOSOV; Vadim Viktorovich LISITSA; Clement Kostov
Original assignee: Schlumberger Canada Ltd; Services Petroliers Schlumberger SA; Schlumberger Technology BV; Schlumberger Holdings Ltd; Prad Research and Development Ltd
Current assignee: Schlumberger Canada Ltd; Services Petroliers Schlumberger SA; Schlumberger Technology BV; Schlumberger Holdings Ltd; Prad Research and Development Ltd
Priority date: 2010-12-21
Filing date: 2010-12-21
Publication date: 2012-06-28
Anticipated expiration: 2013-06-21

Description

METHOD FOR COMPUTING WAVEFIELDS

Field of the invention

This invention relates to methods for numerical simulations of different kinds of wavefields namely, acoustic, elastic, visco-elastic, poro-elastic, or electro-magnetic wavefields. More particularly, the invention relates to methods for solving the acoustic or elastic wave equations for applications of survey design, imaging, interpretation and wavefield inversion of seismic and sonic data.

Background art

Currently, the state-of-the-art in industrial simulations of seismic data is represented by computations of 3D acoustic wavefields (at relatively low frequencies) for applications of survey design and seismic data interpretation, including imaging and wavefield inversion of seismic data. Increasing the bandwidth of the data and solving elastic wave equations, as required for more realistic simulations, increases significantly the computational requirements. The development of algorithms and computer implementations for solving efficiently the elastic wave equations is a very active field of research.

The conventional approaches to solving the acoustic or elastic wave equations are finite-difference methods applied in the time-space domain. Typically, the problem is too large to fit in a single memory space, and therefore various strategies need to be implemented to deal with memory limitations, while solving the wave equation. The two most common approaches are domain decomposition and out-of-core computations (Etgen, J., O'Brien, M., 2007, Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial, Geophysics, vol 72, pp SM223- SM230.). Such approaches complicate and slow down the algorithms (dealing with memory limitations) and require either intensive communication between processors, or large amounts of disk access during the computations. In the other major approach, leading to solution of Helmholtz equations after Fourier transform, parallelization is also difficult to achieve (Operto, S., Virieux, J., Amestoy, P., L'Excellent, J-Y., Giraud, L., Hafedh Ben Hadj, A., 2007, 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study; Geophysics, 72.).

It is an object of the invention to provide a method for improving the efficiency of the calculations and the accuracy of the computed wavefields.

The object of the invention is achieved by a method as set forth in the appended claims.

Summary of the invention

The present invention provides a method for computing wave fields based on the Laguerre transform applied to the time variable in the elastic (or acoustic, or visco-elastic) wave equations which leads to Helmholtz equations. After the Laguerre transform, the associated Helmholtz operators are negative definite, which is a favorable situation for further parallelization by domain decomposition. Here however, the domain decomposition occurs before the solution of the Helmholtz wave equations, and the solver of the Helmholtz equation doesn't need to deal with the domain decomposition. Data exchanges occur between domains only after relatively long independent computations, and therefore, they can be done efficiently via simple disk access, decoupled from the computations. The combination of the Laguerre transform and an effective domain decomposition method such as the Schwartz iterative method is applied to the numerical simulations of sonic and seismic data surveys, for applications of optimizing the survey design and data interpretation workflows.

So the invention comprises the steps of specifying the governing wave equations, setting up the corresponding material parameters, the dimensions of a computational domain, the conditions at the boundaries and the initial conditions. Then the parameters of the Laguerre transform are selected by examining the transform for the initial conditions. Furthermore the invention comprises the steps of decomposing the computational domain into at least two sub-domains, applying Laguerre transform along the time axis for each sub-domain and obtaining Helmholtz equations, solving Helmholtz equations for each sub-domain, exchanging data between sub-domains as per Schwartz iteration method, iterating previous two steps until sufficient level of accuracy is reached, applying inverse Laguerre transform and obtaining the solution of wave equations.

The wave equations can be for acoustic, elastic, visco-elastic, poro- elastic, or electro-magnetic wavefields

Helmholtz equations can be solved using standard finite-difference or finite-element methods.

Helmholtz equations for different sub-domains can be solved by different numerical methods.

Those and other features of the invention will be understood by those skilled in the art from the following detailed description and drawings.

Brief description of the drawings

Fig. 1 shows computational domain for a simple layered medium.

Fig. 2 shows computational domain for more complicated model. Fig. 3 shows the wave field computed for a simple layered medium shown on Fig.1.

Fig. 4 shows the wave field computed for more complicated model shown on Fig. 2

Detailed description of the invention

The detailed steps for applying the proposed method are as follows. At first the task of computing seismic wave fields is defined - e.g. simulation of data, or wavefield propagation as needed for imaging and inversion. This step is common to all computational approaches, and includes specifying the governing equations (e.g. acoustic wave equation, or elastic, or visco-elastic, etc wave equations), the corresponding material parameters (compressional and shear velocities, density, etc), the dimensions of a computational domain and the conditions at the boundaries, the initial conditions (e.g. source functions), what results from computations are saved.

Then the parameters of the Laguerre transform are selected by examining the transform for the initial conditions (e.g. source time function). The computational domain is decomposed into sub-domains, taking into account variations in the complexity of the material parameters, target sizes for the computational tasks in each sub-domain, as well characteristics of the available computing resources. The Laguerre transform is applied along the time axis for each sub-domain and Helmholtz equations obtained after Laguerre transform are solved for each sub-domain, using standard finite- difference or finite-element methods. Data exchanges is made between sub- domains as per Schwartz iteration method. Last two steps (solving Helmholtz equations and data exchange) are iterated until sufficient level of accuracy is reached for the solution. Then the inverse Laguerre transform is applied and the solution to the defined task is obtained. Further, it is proposed to apply different numerical schemes within each of the sub-domains (e.g. finite-difference schemes of different orders in space, or finite-element schemes), thus adapting the computational cost to the accuracy requirements and the earth model properties.

The integral Laguerre transforms (forward and backward) are provided by the formulae:

where ⁼ _«n ^'—(ht)² e ² L"(ht) , andL^a _n (ht) are the classical Laguerre polynomials of degree n, is an auxiliary parameter. The most important variable appearing in the formulae is h which is real and positive. As it follows from the representation of the Laguerre transform, h is a scaling parameter, the higher h is the narrower the Laguerre polynomials are if considered as functions of time. This means that the higher h is the less polynomials are needed to approximate a source wavelet.

Having applied Laguerre transform to scalar wave equation one achieves so-called negative defined Helmholtz equation, i.e. all the eigenvalues of the operator are strictly negative:

K =∑ .(*) (AO

n=0 The main advantages of the obtained equation are: the spectra is well separated from zero, hence - good conditionality; the operator does not depend on degree of Laguerre polynomial, hence can be inverted only ones and used for different right-hand sides; right-hand side depends on the linear combination of the previously computed coefficients, but fi_k do not depend on current number of polynomial (n), hence no need to store all the previous coefficients but only the linear combination.

As it was mentioned above negative-definite Helmholtz equation, obtained by means of Laguerre transform, is elliptic one, hence the well- developed theory of Schwartz domain decomposition method can be applied in straightforward way. Taking into account a number of different ways designed to improve convergence of Schwartz method for elliptical operators one may expect high convergence rate for negative-definite Helmholtz equation.

To illustrate the efficiency of the proposed method a set of numerical experiments was done. For the very beginning ID problem was considered and the comparison of Schwartz iterations was done for both classical and negative-definite Helmholtz equation. It was observed that for variable velocity Schwartz method converged within 45 iterations for each frequency if applied for classical Helmholtz and within only 5 iterations if done for negative-definite one.

2D simulations were done with the help of Laguerre transform only. Wave fields were computed inside the domain of 1 over 1 kilometer. Two models were considered. The first one was rather simple, i.e. horizontally layered - it is shown on Fig. l . The other one was more complicate and designed to observe refracted waves (see Fig. 2). In order to apply Schwartz method the domain was divided into four subdomains - one division in each spatial direction. In Figure 3 one can see the wave field computed for the simple layered medium of Fig. l, and the synthetic wave filed calculated for the second model (Fig.2) is presented in Figure 3. Reflected waves as well as refracted ones are clearly seen.

It is worth mentioning that for both experiments it took Schwartz algorithm only 7 iterations to converges, while the simples conjugation conditions for were used. Applying more specific conditions one may expect significant increase of the convergence rate (number of iterations close to number of divisions along a single spatial direction).

Claims

1. A method for computing wave fields comprising the steps of:

- specifying the governing wave equations,

- setting up the corresponding material parameters, the dimensions of a computational domain, the conditions at the boundaries and the initial conditions,

- selecting the parameters of the Laguerre transform by examining the transform for the initial conditions,

- decomposing the computational domain into at least two sub-domains,

- applying the Laguerre transform along the time axis for each sub-domain and obtaining Helmholtz equations,

- solving Helmholtz equations for each sub-domain,

- exchanging data between sub-domains as per Schwartz iteration method,

- iterating previous two steps until sufficient level of accuracy is reached,

- applying inverse Laguerre transform and obtaining the solution of wave equations.

2. A method of claim 1 wherein wave equations are acoustic or elastic wave equations.

3. A method of claim 1 wherein wave equations are visco-elastic or poro- elastic wave equations.

4. A method of claim 1 wherein wave equations are electromagnetic wave equations.

5. A method of claim 1 wherein Helmholtz equations are solved using standard finite-difference method.

6. A method of claim 1 wherein Helmholtz equations are solved using standard finite-element method.

7. A method of claim 1 wherein Helmholtz equations for different sub- domains are solved by different numerical methods.