WO2014117284A2

WO2014117284A2 - Wave propagation and imaging method

Info

Publication number: WO2014117284A2
Application number: PCT/CH2014/000012
Authority: WO
Inventors: Shravan Hanasoge; Tarje Nissen-Meyer; Johan Robertsson
Original assignee: Shravan Hanasoge; Tarje Nissen-Meyer; Johan Robertsson
Priority date: 2013-01-30
Filing date: 2014-01-30
Publication date: 2014-08-07
Also published as: GB201301620D0; GB2510338A; WO2014117284A3

Abstract

According to a first preferred embodiment of the present invention, there is provided an efficient method of computing synthetic wavefields. In the preferred arrangement, wavefields are computed by means of the Lattice Boltzmann Method (LBM) by increasing the number of degrees of freedom to be greater than the order of the lattice. This new degree of freedom is referred to as colour. Thus, particles in the present invention are fully defined by position, speed, direction and colour. In a preferred embodiment, the present invention is used to image of seismic or electromagnetic wavefield data through forward propagating and back propagating source and receiver wavefields followed by applying an imaging condition. In another preferred embodiment, the present invention is used to invert seismic or electromagnetic wavefield data through an iterative scheme where the difference between modelled and recorded data is minimized by updating an Earth model until the misfit is sufficiently small.

Description

Wave propagation and imaging method

Field of the Invention

The present invention relates to a system and method of modeling acoustic, elastic and electromagnetic waves as well as imaging and inversion, particularly for Earth properties.

Background

The invention relates to modeling and imaging of geophysical data.

Some exploration involves surveying subterranean geological formations for hydrocarbon deposits. A survey can involve seismic and electromagnetic methods. In the following, the invention is described in the context of a seismic survey. However, the invention applies equally well to

electromagnetic surveys (e.g., using so-called controlled source

electromagnetic survey methods). A survey typically involves deploying seismic source(s) and seismic sensors at predetermined locations. The sources generate seismic waves, which propagate into the geological formations by means of pressure changes and vibrations. Variations in the elastic properties of the geological formation scatter the seismic waves, changing their direction of propagation and other properties. Part of the energy emitted by the sources reaches the seismic sensors. Some seismic sensors are sensitive to pressure changes (hydrophones) whereas others are sensitive to either vertical (e.g., geophones) or three-component particle motion (e.g., accelerometers), and industrial surveys may deploy only one type of sensor or both. In response to the detected seismic events, the sensors generate electrical signals to produce time series of seismic data. Analysis of the seismic data can indicate the presence or absence of potential locations of hydrocarbon deposits.

So-called forward modeling is part of many data processing workflows of seismic data. The forward modeling step includes the ability to generate synthetic seismic data at predetermined locations in a synthetic Earth model after exciting one or several sources at predetermined locations in the Earth model. Many forward modeling methods for wave simulations such as seismic, acoustic or electromagnetic data are for instance based on taking discrete differences to approximate the partial differential equations of motion via finite temporal and spatial derivatives such as finite differences, finite elements, or pseudospectral methods (Robertsson, J. O. A., J. O. Blanch, K. Nihei, and J. Tromp (Editors), 2012, Numerical modeling of seismic wave propagation: Reprint volume published by the Society of Exploration Geophysicists ).

These systems describe the macroscopic behavior of waves interacting with a discretized representation of medium parameters such as seismic

wavespeeds and density. Alternatively, discrete particle methods work on a microscopic representation of the physical system, in which the continuum behavior is achieved by interacting individual particles according to specific rules. Many such examples exist based on different types of microscopic approaches such as the Hamiltonian particle method (Takekawa, J.,

Madariaga, R., Mikada, H., Tada-nori, G. 2012, Numerical simulation of seismic wave propagation produced by earthquake by using a particle method, Geophysical Journal International, doi: 10.1111/j.1365- 246X.2012.05676.x), and the discrete particle method (O'Brien, G. 2008, Discrete visco-elastic lattice methods for seismic wave propagation,

Geophysical Research Letters, 35). Further seismic imaging methods are described in the international application WO 2012/128873 and by Walsh et al, Accelerating geoscience and engineering system simulations on graphic hardware, Computers and Geosciences, Vol. 35, 2009, 2353-2364.

The Lattice Boltzmann (LB) Method (LBM) has considerable precedence in the realm of fluid mechanics and multi-phase flow, where it is used to simulate non-linear compressible flow. Recent years have seen the extension of LB methods to systems such as linear wave propagation. The extension of LBM to solve the linear electro-magnetic wave equation can be traced back to Dellar (Dellar P. 2010, European Physics Letters, 90, 50002), although it did not address the general problem, whereas Hanasoge, Succi & Orszag (Hanasoge, S. M., Succi, S., and Orszag, S. A. 2011 , European Physics Letters, 96, 14002) described a general means of solving the Maxwell equations for electromagnetic wave propagation using LB techniques.

Hanasoge et al. numerically solve the electromagnetic wave equation applying the Lattice Boltzmann method, specifically studying the propagation of light through complex media. They used a limited definition of the

distribution function, considering a specific pseudo-vector distribution function, and applied it to Maxwell's equations of electromagnetism.

Earlier attempts related to the Boltzmann equation for acoustic and elastic waves include, among others, Mora (Mora P. 1992, The Lattice Boltzmann Photonic Lattice Solid, Journal of Statistical Physics, 68, 3) and Guangwu (Guangwu, Y. 2000, A Lattice Boltzmann Equation for Waves, Journal of Computational Physics, 161 , 61). O'Brien et al. (O'Brien, G., Nissen-Meyer, T. and Bean, C. 201 1 , A Lattice Boltzmann Method for Elastic Wave

Propagation in a Poisson Solid, Bulletin of the Seismologica! Survey of America, 102, 3) presented a first application of the LBM to elastodynamics, however only for isotropic Poisson media (i.e., elastic media where the two Lame constants are equal), thus failing to accommodate the full range of elastic parameters relevant to heterogeneous media.

Wave propagation can be computed using 2D, 2.5D or 3D assumptions. For 2D wave propagation computations, line sources and line receivers are assumed and a model symmetry around the plane of computation. In 2.5D wave propagation simulations, point sources and point receivers can be assumed but the model is symmetric around the plane of interest. In 3D wave propagation both models and sources/receivers are general and best approximate real experimental data. The term 4D refers to repeat

experiments when a 3D seismic survey is carried out at two different stages to monitor changes in the sub-surface for instance caused by hydrocarbon production.

Finally, waves can also propagate in media described by acoustic,

viscoacoustic, viscoelastic or poroelastic wave propagation (Moczo, P., J. O. A. Robertsson and L. Eisner, 2007, The Finite-Difference Time-Domain Method for Modeling of Seismic Wave Propagation, pp. 421-516, in Advances in wave propagation in heterogeneous Earth (eds. R.S. Wu and V. Maupin) Vol 48, Advances in Geophysics (ed. R. Dmowska), Elsevier-Pergamon, Oxford; Carcione, J. M., and Quiroga-Goode, G., 1995, Some aspects of the physics and numerical modeling of Biot compressional waves: J. Comp.

Acoust., 3, 261 ). In addition, a particularly useful method for simulating wave propagation in anisotropic media is based on the pseudo-acoustic

approximation for tilted transversely anisotropic media introduced by

Alkhalifah (Alkhalifah, T., 2000, An acoustic wave equation for anisotropic media: Geophysics, 65, 1239-1250). Based on Alkhalifah's approximation, different space/time-domain wave equations have been proposed (e.g., Zhou H., Zhang, G., and Bloor, R., 2006, An anisotropic acoustic wave equation for VTI media: 68th EAGE conference and exhibition, extended abstracts, H033) that for instance can be used to image seismic data in anisotropic media without using a full anisotropic elastic wave equation.

Summary of the Invention

According to a first preferred aspect of the present invention, there is provided an efficient method of computing synthetic wavefields. In the preferred arrangement, wavefields are computed by means of the LBM by increasing the number of degrees of freedom to be greater than the order of the lattice. This new class of degrees of freedom is referred to as "colour". Thus mesoscopic particles in the present invention are fully defined by position, speed, direction and colour.

The colour parameter can be, for example, created through reverse engineering from the PDE of interest such that when the LB equation is averaged over the lattice, it produces the desired continuum equations.

Preferably, a general tensorial approach for distribution functions is applied, enabling the use of LBM to solve a range of problems in physics and engineering that has not been possible in the prior art. Specifically, the application of this method to the problem of wave propagation in elastic media is unique and its success in solving a broad range of wave equations relevant to exploration seismology is seen as surprising since prior art in this area claims Lattice Boltzmann fails at the general terrestrial seismology problem (O'Brien, Nissen-Meyer and Bean 201 1 , Bulletin of the Seismological Survey of America).

The method is best used to propagate waves through a representation or model of the physical body through which the waves travel. The waves are injected into such a physical body using controlled wave sources such as vibratory or explosive sources, sonic or ultrasonic transducers and the like.

By way of summary, in a preferred embodiment, the present invention is used to image sub-surface properties using seismic or electromagnetic data by means of wavefield computations. In one particular preferred embodiment, an image of the sub-surface is produced through a sequence of steps where:

1. a background Earth model is defined,

2. source wavefields are forward propagated through an Earth model,

3. recorded wavefields are back propagated through the same Earth model,

4. an image of the sub-surface is produced through a correlation between back propagated recorded wavefield and forward propagated source wavefield,

wherein the LBM is used as the wavefield propagator in at least one of the steps 2 and/or 3.

By way of summary, in another preferred embodiment, the present invention is used to invert seismic or electromagnetic wavefield data for sub-surface properties. In one particular preferred embodiment, a sub-surface property volume is produced through a sequence of steps where:

1. a starting (i.e., initial background) Earth model is defined,

2. synthetic geophysical wave information is computed using the LBM and the Earth model, 3. the difference between recorded and synthetic geophysical wave information is computed using an appropriate measure,

4. the Earth model is updated based on the difference between recorded and synthetic data,

5. the process is repeated from steps (2)-(4) until the difference between synthetic and recorded data is sufficiently small and a new Earth model found.

In the present invention, the use of a so-called Lattice-Boltzmann numerical method is described to simulate linear (and non-linear) elastic, viscoelastic, poroelastic, acoustic, viscoacoustic and electromagnetic wave propagation in heterogeneous, isotropic and anisotropic media with or without external or internal topography. Applications relate to heterogeneous domains at the scale of the Earth's crust or smaller, relevant for hydrocarbon and mineral exploration. Specifically, the method is intended to replace the forward modeling step in, for instance, acoustic and elastic or electromagnetic imaging and inversion methods for the Earth's subsurface. To obtain seismic images of the subsurface, one needs to solve the wave equation (synthetics) and combine this solution with actual observations. In imaging, one can, for instance, propagate the wavefield due to an assumed source, and back- propagate recorded data to satisfy an imaging condition wherever those two wavefields interact constructively. Full-waveform inversion (FWI) generalizes the back-propagating source as a misfit between data and synthetics, and allows for iterative improvements of the sub-surface image by gradient techniques. In the general framework, this includes modeling the full seismic wave equation. In practice, this modeling step needs to be undertaken 10M- 10^Λ7 times (number of sources times iteration count) in modern

imaging/inversion techniques before a reliable image is obtained. Thus, it is of utmost importance to tune and optimize these forward modeling approaches towards the intended application in terms of the necessary and desired complexity in structure and modeling, including uncertainties. The method described (the Lattice-Boltzmann method for wave propagation) is based on the motion of mesoscopic particles along a lattice, whose continuum properties correspond accurately to the desired wave equation. It runs efficiently and highly scalable on parallel computers, and suits GPU architectures particularly well due to the nature of its bulk operations akin to multi-threading. The Lattice-Boltzmann (LB) Method (LBM) is especially useful for coupling domains of different physical systems (e.g. solid-fluid domains, nonlinear and linear wave propagation, elastostatics, brittle failure, interaction between waves and porous fluid flow).

The present invention introduces the use of special vector and tensorial distribution functions in conjunction with methods of LB. It does not suffer from limitations in the prior art such as only being applicable to limited classes of sub-surface materials (e.g., Poisson solids), offering a complete

computational engine for the hydrocarbon exploration applications mentioned above.

The foregoing has outlined in broad terms features of preferred embodiments of the invention disclosed herein so that the detailed description that follows may be more clearly understood, and so that the contribution of the present invention to the art may be better appreciated. The present invention is not to be limited in its applications to the details of the construction and to the arrangements of the components set forth in the following description or illustrated in the drawings. Rather, the invention is capable of other

embodiments and of being practiced and carried out in various other ways not specifically enumerated herein. Finally, it should be understood that the phraseology and terminology employed herein are for the purpose of description and should not be regarded as limiting, unless the specification specifically so limits the invention. For example, the term "backward

propagation" can include adjoint propagation or to compare measured and propagated waves can include comparing observed ground or body motion and the propagated/modelled wave. Further the invention can have applications to fields other than seismic imaging, for example in the examination of structural damages or other internal properties of industrially manufactured bodies, imaging of the interior of the human body, using sound, ultrasound, or other types of waves.

Brief Description of the Drawings

Exemplary embodiments of the invention will now be described, with reference to the accompanying drawing, in which:

FIG.1 shows a flowchart illustrating the role of a modeling method to solve geophysical modeling (forward problem) and imaging/inversion problems (inverse problem);

FIG. 2 is a flow chart with steps in accordance with an example of the

invention, and

FIG. 3 is a comparison between waves propagated in accordance with an example of the invention and in accordance with a known method.

Detailed Description of the Invention

The present invention is related to modeling, imaging and inversion of seismic data. With reference to Figure 1 , the relation between these different applications is illustrated. "Earth model space" in step 10 refers to a computer representation of the Earth's subsurface structure and/or properties. "Seismic data space" in step 11 refers to seismic data that would have been acquired or has been acquired over said corresponding Earth model in step 10. The process of going from the Earth model space in step 10 to the seismic data space in step 11 is referred to as the forward problem or forward modeling in step 12. The process of going from the seismic data space in step 11 to the Earth model space in step 10 is referred to as the inverse problem or imaging problem in step 13. At the heart of both the forward problem in step 12 and the inverse problem in step 13 is a computational engine involving simulating wave propagation problems in step 4, enabled by the present invention. It is therefore clear that the present invention has applications both in forward modeling of seismic data and imaging and inversion of seismic data.

The task of predicting signals observed at various recording stations due to specific externally applied stimuli (explosions) is computationally very expensive, especially under demands of high accuracy and fidelity in modeling a complex 3-D heterogeneous Earth. Currently, the only viably fast methods of solution, finite differences, are subject sometimes to significant errors and the cost of computation therefore becomes very high to mitigate for these errors. Further, complex topography and varying mesh spacing is hard to incorporate in these methods.

In the present invention these issues are addressed by using the LBM with innovations described here. Classical LB methods have been successfully designed to solve the Navier-Stokes equations (i.e., fluid mechanics). The Navier-Stokes equations may be viewed as a continuum representation of a kinetic theory of microscopic particle motion. The classical LBM solves the kinetic equation albeit with mesoscopic particle units, moving on a user- specified lattice. Lattices consist of nodes and paths radiating from these nodes along which particles are constrained to move and are denoted by terms such as "D3Q27", "D2Q4" etc. The lattice D3Q27 is 3 dimensional and has 27 paths connecting each node to a neighbouring set of other 27 nodes. Similarly, D2Q4 is a 2 dimensional lattice with 4 paths radiating out from each node. The "order" of a lattice increases with the number of paths; high-order lattices contain more paths than low-order ones.

The emergence of the continuum equations is assured given the right choice for the mesoscopic canonical ensemble. The classical LBM is second-order accurate, maintained everywhere including at the free-surface boundary. It can account for variable mesh spacing and external as well as internal topography; mesh generation time is trivial. Recent developments have

towards it. Correspondingly, the velocities that the particles can assume are denoted by c_r

A departure from previous approaches involves broadening the class of distribution functions g. to a more general g_ia, where a is a vector index, to allow for the simulation of the full anisotropic seismic wave equation. The new class of distribution functions is tensorial, vectorial or may equivalently be termed as containing multiple components, more than allowed by the classical choice of scalar functions for each lattice connector. Thus there are more degrees of freedom than the number of lattice connectors radiating from each node. These additional degrees of freedom is referred to as colour, an additional descriptor of the particles. Thus, populations of particles are characterized by four properties: position, speed, direction, and colour.

Generally, one or more colour(s) can be derived by

1. Posing the relevant wave equation as a transport equation (see for example eq. (1) above and (2) below)

2. Construct a tensor equilibrium distribution function (see for example, equation (3) below)

3. Use the standard Lattice Boltzmann methodology (see for example the sample code below).

As an example of the above applied to the field of seismic imaging, the general seismic wave equation can be written as

(2) where Q_a contains various wavefield properties, Τ_αβγ (χ) is a user-specified wave-propagation- physics encapsulating tensor and f_a x, t) is a forcing function, specified also externally (see for instance Dumbser, M., and Kaeser, M., 2006, An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes - II. The three-dimensional isotropic case, Geophysical Journal International, 167, 319-336).

Then the distribution function is postulated to obey the conservation equation (1 ). Furthermore, lattices are chosen so that they obey the following properties:

If one also assumes that the distribution function were to satisfy

(8)D one obtains

(9)D

This produces the desired continuum equation: d,Q_a + djr_apyQy) = fa a >

(10)

A quantity of central importance in the LBM is the equilibrium distribution function. It is a descriptor of the interaction of the mesoscopic particle units and relates continuum properties to the local density of particles. In prior art, although not always, it has been derived based on a kinetic theory, limiting the applicability of LBM. In the present invention, these constraints are eliminated. The present invention is not restricted to scalar functions; instead, user- prescribed vector and tensorial definitions of the equilibrium distribution functions are employed. When used in conjunction with classical LB

methodology, these vectorial/tensorial equilibrium distribution functions, result in the solution of various wave equations, of direct utility and applicability to imaging methodologies such as reverse time migration (RTM; see Etgen, J., Gray, S. H., and Zhang, Y., 2009, An overview of depth imaging in exploration geophysics: Geophysics, 74, WCA5-WCA17, and references therein), oneway wave equation migration (OWEM; Etgen et a!., 2009, and references therein), and full-waveform inversion.

Typically in the prior art, one scalar distribution function is defined for each connector of the lattice. Thus, for instance, the D2Q4 lattice (simplest lattice in 2D), that contains one connector respectively from North, East, West and South (i.e., directions along the Cartesian axes), would have a scalar distribution function with four components for each of the 4 connectors from each node. O'Brien et al. found that this form of distribution is able only to reproduce seismic wave propagation in Poisson solids. In the present invention, the number of degrees of freedom is increased and a multi- component scalar distribution function or a vectorial or tensorial distribution function is used, defined for each lattice connector. In other words, more distribution functions are defined than the order of the lattice i.e., more distribution functions than the number of connections from each node on the lattice. This new degree of freedom is referred to as "colour". Thus,

mesoscopic particles in the present invention are fully defined by position, speed, direction and colour. In contrast, in prior art for the seismic wave equation (O'Brien et al. 2012), particles have had only three descriptors, position, speed and direction.

The classic LB algorithm is employed to compute the evolution of the populations of the coloured mesoscopic particles, specifically the steps of streaming and collision. These algorithmic steps are well known to

practitioners of LB methods (e.g., Guangwu, Y. 2000, A Lattice Boltzmann Equation for Waves, Journal of Computational Physics, 161 , 61). Particles are streamed (or moved) based on their current direction and speed. Only particles of the same colour, direction and position are allowed to collide. Alternately, every component of the special tensorial distribution function is subjected to the steps of streaming and collision in the same identical manner, as would be done using a scalar function (e.g., O'Brien et al., 2012). The continuum properties may be easily recovered from the distribution function at any point in space and time.

An example of a realization of the LB method for a special distribution function involves the steps shown in Fig 2. In step 21 , particles are streamed, i.e., moved from one lattice node to another. This is done for each component of the tensorial distribution function. In step 22, the continuum variables are recovered using the definition of the distribution function. The equilibrium distribution function is reconstructed from these continuum variables in step 23 and the particles collided at each lattice node in step 24. Finally, a special set of vectorial source terms is used to force the equations in step 25. These source terms may comprise of not just wave sources but also terms whose function is to ensure that the resultant continuum equation is the desired one. Following steps 21-25 results in evolving the system forward by one time step. Thus to evolve it temporally forward by a desired amount, the process must be iterated, as is indicated in step 26.

Two example MATLAB codes with documentation follow below. The first code shows the classical Lattice Boltzmann approach, where there are as many distribution functions as lattice connectors from each node. The second code shows a realization of a code with the innovation described here.

% Iterate over required number of time steps

for it = l:nsteps,

% FOR A SIMPLE 4-direction LATTICE (particles are

% constrained to move left, right, up and down, i.e.,

% the Cartesian directions) .

% g a three dimensional array is the distribution of

% particles, two spatial coordinates and the lattice

% direction (1, 2, 3, 4) representing (East, North, West,

% South) respectively, nx and ny grid points are used in

% the x and y directions respectively.

% STREAMING - PARTICLES MOVE ALONG THE LATTICE - STEP 21 % (PART OF CLASSICAL LATTICE-BOLTZMANN METHODOLOGY)

% TRANSFERRING PARTICLES FROM ONE LATTICE POINT TO

% ANOTHER BASED ON THEIR DIRECTION OF MOVEMENT.

% MOVING PARTICLES NORTHWARD

for j=ny:-l:l,

g ( : , j , 2) = g ( : , j-1, 2) ;

end;

% MOVING PARTICLES SOUTHWARD

for j=l:ny,

g(:,j, 4 ) = g(:, j+1, 4 ) ;

end;

% MOVING PARTICLES WESTWARD

g(l:nx,:,3) = g (2 : nx+1, : , 3) ;

% MOVING PARTICLES EASTWARD

g(nx:-l:l, : ,1) = g ( (nx-1) : -1 : 0, : , 1) ;

REPRESENTING VARIABLES OF THE PHYSICAL SYSTEM

ux ( : , : ) , uy ( : , : ) , sigma ( : , : ) > AS PARTICLES

; By construction, sigma ( : , : ) = sum_i ( g(:,:,i) ) ; By construction, ux ( : , : ) = sum_i ( cx_i * g(:,:,i) ) ; By construction, uy ( : , : ) = sum_i ( cy_i * g(:,:,i) ) ; i denotes the directions that particles can move sigma (:,:) = g(:,:,l) + g(:,:,2) + ...

g ( : , : , 3 ) + g ( : , : , 4 ) ;

,

% DIRECTION

% Set omega = 2.0 for no dissipation, g = omega* geq + (1. - omega) *g; % ADD SOURCE TERMS

% TO MANIPULATE THE CONTINUUM AVERAGE temp = force. ^source (it) . *w; % ADDING SEISMIC SOURCE

% IN THIS CASE the same seismic source is added to all ribution function i) + temp;

This second code shows a realization of a code with the innovation described here, i.e., where vectorial (or tensorial) distribution functions are used, with more components than lattice connectors radiating from each node.

We consider simulating waves through a vertically transversely isotropic medium (VTI) using the equation described by Du et al. (Du, Fletcher, and Fowler 2008, A New Pseudo-Acoustic Equation for VTI Media, EAGE

Conference & Exhibition, Rome, article H033) :

d q 7 d p d²q ,_{Λ Λ \}

For the pseudo-acoustic VTI equation described above (eq. 1 1), there exists a transformation mapping this system of equations to the tensor form in equation (2). In 2-D, the tensor has 5 X 2 X 5 components. The first and third dimension have a size of 5 because there are two components of

displacement and three independent components of stress in 2-D. And since the problem is in 2-D, the second dimension of matrix T has a size of 2. In 3- D, the size of matrix will be 9 X 3 X 9, since there are three components of displacement and 6 independent components of stress (see Dumbser and Kaeser, 2006, Geophysical Journal International). In this 2-D problem, 5 colours are therefore introduced, i.e., 5 types of particles, each with a different colour, representing two displacements and three stresses. Defining ncol = 5, a sample code can be shown that can be used to simulate wave propagation in this pseudo-acoustic VTI medium.

% IMBUING PARTICLES WITH COLOUR (INNOVATIVE)

% ALTERNATELY: USING TENSORIAL/VECTORIAL

% DISTRIBUTION FUNCTION

% ncol >=2 is number of colours

% THE DISTRIBUTION FUNCTION HAS FOUR INDICES,

% 2 spatial, one for direction and one for colour

% STREAMING

% ALL THE DIFFERENT SETS OF PARTICLES ARE STREAMED

% IN THE SAME MANNER

% MOVING PARTICLES NORTHWARD

for i=l:4,

geq( : , : = (cx(i) * ux ( : , : , j ) + ...

cy(i) * uy(:,:,j) + sigma ( : , : , j ) ) . *w; end;

end;

% INTERACTING PARTICLES WHICH POSSESS THE SAME MOVEMENT

DIRECTION AND SAME COLOUR

% Set omega = 2.0 for no dissipation. g = omega* geq + (1. - omega) *g;

% ADD SOURCE TERMS

% TO MANIPULATE THE CONTINUUM AVERAGE temp = force. ^source (it) . *w; % ADDING SEISMIC SOURCE temp = temp + ux * dx_rho0 + uy * dy_rho0;

% ADDING TERMS TO THE SOURCE SO AS TO MODIFY THE

% DIFFERENTIAL EQUATION (INNOVATIVE) for j=l:ncol,

for i=l:4,

g(:,:,i,j) = g(:,:,i,j) + temp;

end;

A modified version of the above code is presented below to emphasize the introduction of the colour parameter: IMBUING PARTICLES WITH COLOUR (INNOVATIVE) ALTERNATELY: USING TENSORIAL/VECTORIAL

DISTRIBUTION FUNCTION ncol >-2 is number of colours

THE DISTRIBUTION FUNCTION HAS FOUR INDICES,

2 spatial, one for direction and one for colour

STREAMING

ALL THE DIFFERENT SETS OF PARTICLES ARE STREAMED IN THE SAME MANNER

MOVING PARTICLES NORTHWARD

or col=l:ncol,

for j =ny : -1 : 1 ,

g( : ,j,2,col) = g ( : , j -1, 2, col) ;

end; end;

% MOVING PARTICLES SOUTHWARD

for col=l:ncol,

for j=l:ny,

g(:,j,4,col) = g ( : , j +1 , 4 , col) ;

end;

% MOVING PARTICLES WESTWARD

g ( 1 : nx, : , 3 , 1 : ncol) = g (2 : nx+1, : , 3, 1 : ncol) ;

% MOVING PARTICLES EASTWARD

g (nx: -1 : 1, : , 1 : ncol) = g ( (nx-1) : -1 : 0, : , 1, 1 : ncol) ;

% REPRESENTING VARIABLES OF THE PHYSICAL SYSTEM

% Q_l, Q_2, .... ,Q_ncol;

% AS PARTICLES

% By construction, Q_j = sum_i ( cx_i * g(:,:,i,j) ) % i denotes the directions that particles can move

% j denotes the colour of the particles for j=l:ncol,

Q(:, :,j) = g(:, :,l,j) + g(:, :,2, j) + ...

g ( : , : , 3 , j ) + g ( : , : , 4 , j ) ; end;

% COMPUTING PARTICLE DISTRIBUTION GIVEN VARIABLES OF THE % PHYSICAL SYSTEM (Q_j).

% THE MATRIX T_a/fyIS REWRITTEN AS TWO MATRICES Tx_{jk}, Ty_{jk}, where the indices j,k run between 1 through 5. a,

% By construction, geq__{ij} = [sum_k (cx__i * Tx__{jk} + cy_i * Ty__{jk})* Q__k) + Q_j ] * w_i,

% where i denotes the directions that particles can move

% and j denotes the colour of the particles

% w = weights .

% See equation (3) .

% The sum over k in the definition of geq_{ij} mixes the % different coloured particles geq (:,:,:,:) = 0.0;

for j=l:ncol,

for i=l:4,

for k=l:ncol,

geq(:, :,i, j) = (cx(i) * Tx ( : , : , j , k) * Q ( : , : , k) + ... cy(i) * Ty(:,:,j,k) * Q(:,:,k) + Q ( : , : , j ) ) . *w + ... geq( : , : , i, j ) end;

end; end;

% INTERACTING PARTICLES WHICH POSSESS THE SAME MOV:

DIRECTION AND SAME COLOUR

% Set omega = 2.0 for no dissipation. g = omega* geq + (1. - omega) *g;

% ADD SOURCE TERMS

% TO MANIPULATE THE CONTINUUM AVERAGE temp = force. *source (it) . *w; % ADDING SEISMIC SOURCE for j=l:ncol,

for i=l: 4 ,

g ( : , : , i , ) = g ( : , : , i , j ) + temp ;

end;

In the comparative plots of FIG. 3, which compare wave propagation modelled using the above method (LB) and a finite-difference method, vx refers to p and vy refers to q as defined in equation (11). The agreement is excellent, demonstrating the ability of LBM to model a broad range of wave equations. FD refers to finite differences and LB to Lattice Boltzmann. Using the method, TTI media or fully anisotropic elastic wave equations can be modelled and it is possible to jointly treat such different regions containing respective medium properties due to the generality of the method which simply requires different distribution functions for each of the domains

Depth imaging of seismic data as illustrated in Fig. 1 , step 13, is often thought of as a three-step process (Etgen et al., 2009). First, a velocity model of the subsurface is estimated by some means (e.g., tomography or full waveform inversion). Secondly, the seismic data are propagated into the subsurface in this velocity model using some kind of numerical method for wave

propagation as in Fig. 1 , step 14. This is the heart of the imaging algorithm and one important application of the present invention. Thirdly, an imaging condition is carried out to form the image. To image seismic data, a source wavefield is combined with a receiver wavefield. The source wavefield is sometimes computed on a velocity model using a synthetic or estimated source wavelet. Other times, the source wavefield is estimated by decomposing the recorded data into an up-going and a down-going wavefield at the recording locations (e.g., Muijs, R., J. O. A. Robertsson and K. Holliger, 2007, Prestack depth migration of primary and surface-related multiple reflections, Part !: imaging: Geophysics, 72, S59- S69). The down-going wavefield would then correspond to the source wavefields. The receiver wavefield contains as a "source" most commonly the data recorded at the receivers. However, in other cases it may constitute the up-going wavefield at the receiver locations (Muijs et al., 2007).

Source and receiver wavefields need to be downward continued into the subsurface. This is sometimes done using the one-way wave equation (OWEM; e.g., see Etgen et al., 2009, and references therein), in such a case, the wavefield is propagated to successively deeper depth levels. Other times, the wavefield is propagated into the subsurface using the two-way wave equation sometime referred to as Reverse Time Migration (RTM; e.g., see Etgen et al., 2009, and references therein), !n such a case, wavefields propagate in all spatial directions at each time step. Forward propagation of the source wavefield is done in the forward time direction. For back propagation of the receiver wavefield, the wavefield is propagated backwards (i.e., opposite to the actual propagation direction) in negative time. For a two- way wave equation it would then ideally focus back on the source from which it originated.

The forward- and back-propagated wavefields are combined to produce a seismic image using various imaging conditions. Sometimes, the imaging condition is as simple as the center-lag of the cross-correlation of the two wavefields. Alternatively, the imaging condition can be more sophisticated, involving temporal and/or spatial deconvolutions of one wavefield from the other. In the present invention, the term correlation is used to refer to all types of imaging conditions. Full waveform inversion (FWI) can be seen as a generalization of the imaging procedure described above. The key concept is the definition of a misfit functional between recorded and simulated waveforms, to be iteratively minimized by updating Earth properties (e.g. velocity model). The misfit functional is usually given by a well-chosen time window of a L2 norm between the (possibly frequency-filtered) waveforms, a phase or traveltime representation. As the iterations progress, the definition of the misfit may change, e.g. altering the time window, or increasing the frequency content. The model update requires knowledge of the gradient of the misfit functional with respect to parameter variations (such as different types of wave speed, impedance). This gradient is then used for iterative gradient schemes such as steepest descent, conjugate gradient, or Gauss-Newton methods to minimize the misfit functional. Obtaining the gradient is the most computationally demanding aspect of FWI: At each iteration, one needs to simulate the entire dataset for each shot to all geophones and backwards, similar to reverse time migration. Thus, an FWI requires at least 2^*Nit^*Nshots simulations. A typical dataset may consist of tens or hundreds of thousands of shots (Nshots), and achieving an acceptable convergence level may take up to hundreds of iterations Nit. For each shot, forward and backward wavefields are convolved with the corresponding backward wavefields, and summed over all shots to gradually assemble the gradient. Choices for the misfit functional, filtering, and model parameterization are subjective and may vary not only from one application to another, but within one FWI framework. FWI is a powerful technique, which attempts to exploit the entire information in recorded data to derive sharper images of the subsurface than, for example, by conventional methods such as ray tomography.

It is important to note that any inversion procedure such as FWI is non- unique, and convergence to a local minimum poses a danger. It is therefore paramount to have as accurate a starting model as possible before conducting FWI. The starting model may for instance be a state-of-the-art tomography model of the region of interest, and FWI is mainly used to sharpen the image.

Since the advent of FWI, most applications have relied on frequency-domain solvers to the 2D acoustic wave equation. This is a feasible approach especially if one iterates towards higher and higher frequencies as the inversion procedure refines its imaging capabilities with each step towards lower misfits. However, frequency-domain solvers suffer from the caveat of requiring global communication when inverting a matrix system. For elastic applications, and any applications in 3D, frequency domain solvers are incapable of resolving the frequencies of interest even on modern

supercomputers. Thus, most 3D applications have turned towards time- domain forward solvers to compute the gradient of the misfit functional, which can be conveniently accomplished by the adjoint method.

Full waveforms encapsulate effects of wave propagation such as diffraction, caustics, 3D elastic focusing or evanescent waves which may have been ignored by approximate techniques (e.g. ray theory, acoustic media). They contain interface waves sensitive to solid-fluid boundaries or salt flanks, and wave effects due to undulating surface or bathymetric topography. To fully replicate high-frequency waveforms via modeling and project their nature back to their physical origin within the subsurface, it is crucial to have a modeling engine at hand which can handle such complexities, while at the same time runs sufficiently fast to tackle the computational burden outlined above.

In summary, 3D FWI benefits from a fast forward solver for complex media including sharp interfaces, solid-fluid boundaries, diverse rheologies (e.g. poroelastic media), surface topography and bathymetry. Each iteration of the FWI requires a complete re-definition of model parameters in the forward solver. Thus, an efficient method to re-mesh and re-populate 3D wavespeeds is also of central concern for feasible FWI applications. The present invention satisfies all of these constraints, and is quite uniquely positioned as such: Re-meshing of 3D models is a critical issue for any sophisticated numerical technique such as finite or spectral elements, whereas complexities such as multiple rheologies or surface topography pose a formidable challenge to conventional methods such as finite differences.

As the present invention has been described above purely by way of example, the above modifications or others may be made within the scope of the invention. The invention may also comprise any individual features described or implicit herein or shown or implicit in the drawings, or any combination of any such features or any generalisation of any such features or combination, which extends to equivalents thereof. Thus, the breadth and scope of the present invention should not be limited by any of the above-described exemplary embodiments. Alternative features serving the same, equivalent or similar purposes may replace each feature disclosed in the specification, including the drawings, unless expressly stated otherwise, for example using the principles as described above to acoustic waves propagating in an acoustic medium or electromagnetic waves.

Unless explicitly stated herein, any discussion of the prior art throughout the specification is not an admission that such prior art is widely known or forms part of the common general knowledge in the field.

Claims

1. A method for representing the propagation of waves injected by

controlled wave-generating sources into a physical body in a medium being a representation of the body through which the injected waves travel , comprising the steps of:

a. representing the wave in the medium as a discrete collection of particles at locations defining a lattice,

b. defining paths interconnecting nodes, and specifying possible directions and speeds between nodes, in the lattice,

c. classifying the collection of particles according to location, speed, direction and at least one additional parameter representing an additional degree of freedom for the particles,

d. moving particles according to their instantaneous direction of movement to neighboring lattice locations,

e. identifying at least one set of particles that share the same location, speed, direction and at least one additional parameter representing an additional degree of freedom for the particles,

f. interacting this set of particles through a specified interaction rule.

2. The method of claim 1 comprising the further step of locating

inhomogeneities within the physical body by propagating the waves forward or backward to a location where the waves are accessible for measurement by sensors and comparing measured and propagated waves.

3. The method of claim 1 or 2 where the medium is selected from a group consisting of an elastic medium, a non-Poisson medium, an anisotropic elastic medium, an acoustic medium, a poroelastic medium, a

viscoacoustic medium, a viscoelastic medium, and a pseudo anisotropic acoustic medium, whereby any medium can be of varying density and/or wave propagation speed.

4. The method any of the preceding claims wherein the propagation is determined by a hyperbolic systems of wave equations and wherein the body is a subsection of the earth.

5. The method of any of the preceding claims for propagating waves in two dimensions (2D) and/or in three dimensions (3D)

6. The method of any of the preceding claims further comprising the step of determining the difference in wave propagation at essentially the same locations caused by changes in sub-surface properties (4D).

7. The method of any of the preceding claims, using lattice orders of

accuracy greater than 2.

8. The method of any of the preceding claims, using lattice-decomposed domains with different mesh spacings.

9. The method of any of the preceding claims, using lattices that are

designed to possess isotropic errors or on the basis of entropic arguments

10. A method for depth-migrating common source or common receiver

gathers, or both, of geophysical wave information using an appropriate Earth model, comprising:

(a) forward-propagate a source wavefield through an Earth model,

(b) back-propagate the recorded wavefield through an Earth model,

(c) produce an image of the sub-surface through a correlation between the back propagated recorded wavefield and the forward propagated source wavefield,

wherein the Lattice-Boltzmann method is used in at least one of steps (a) and (b).

1 1. The method of claim 10 using a Lattice-Boltzmann method in

accordance with any of claims 1 - 9.

12. The method of claim 10 or 11 being applied to one or more of reverse- time migration imaging method and a one-way wave equation migration.

13. The method of any of claims 10 to 12 where geophysical waves

are seismic waves or electromagnetic waves.

14. The method of any of claims 10 to 13 for propagating waves in two

dimensions (2D) and/or in three dimensions (3D) .

15. The method of any of claims 10 to 14 further comprising the step of

determining the difference in wave propagation at essentially the same locations caused by changes in sub-surface properties (4D).

16. The method of claim 10 wherein the Earth model comprises at least one of an elastic medium, a non-Poisson medium, an anisotropic elastic medium, an acoustic medium, an acoustic medium with variable or constant density and variable sound-speed, a pseudo anisotropic acoustic medium, a poroelastic medium, a viscoacoustic medium, a viscoelastic medium, whereby any medium can be of varying density and/or speed of sound.

17. The method of any of claims 10 to 16 wherein the waves are forward or backward propagated based on a hyperbolic systems of wave equations and used for hydrocarbon exploration.

18. The method of any of claims 10 to 17, using lattice orders of accuracy greater than 2.

19. The method of any of claims 10 to 18, using lattice decomposed

domains with different mesh spacings.

20. The method of any of claims 10 to 19 using lattices that are designed to possess isotropic errors or on the basis of entropic arguments.

21 . A process of determining an Earth property model by means of full

waveform inversion applied to geophysical wave information,

comprising:

(a) defining an Earth model,

(b) computing synthetic geophysical wave information using the Lattice- Boltzmann Method and the Earth model,

(c) defining a difference between recorded and synthetic geophysical wave information using an appropriate misfit measure,

(d) updating the Earth model based on the misfit between recorded and synthetic data,

(e) repeating the process from (b)-(d) until the difference between synthetic and recorded data is sufficiently small.

22. The method of claim 21 using a Lattice-Boltzmann method in

accordance with any of claims 1 -9.

23. The method of claim 21 or 22 being applied to one or more of reverse- time migration imaging method and a one-way wave equation migration.

24. The method of any of claims 21 to 23 where geophysical waves

are seismic waves or electromagnetic waves.

25. The method of any of claims 21 to 24 for propagating waves in two

dimensions (2D) and/or in three dimensions (3D) .

26. The method of any of claims 21 to 25 further comprising the step of

27. The method of claim 21 wherein the Earth model comprises at least one of an elastic medium, a non-Poisson medium, an anisotropic elastic medium, an acoustic medium, an acoustic medium with variable or constant density and variable sound-speed, a pseudo anisotropic acoustic medium, a poroelastic medium, a viscoacoustic medium, a viscoelastic medium, whereby any medium can be of varying density and/or speed of sound.

28. The method of any of claims 21 to 27 wherein the waves are forward or backward propagated based on a hyperbolic systems of wave equations and used for hydrocarbon exploration.

29. The method of any of claims 21 to 28, using lattice orders of accuracy greater than 2.

30. The method of any of claims 21 to 29, using lattice decomposed

domains with different mesh spacings.

31 . The method of any of claims 21 to 30, using lattices that are designed to possess isotropic errors or on the basis of entropic arguments.