METHOD OF OPERATING A DATA-PROCESSING SYSTEM FOR THE SIMULATION OF THE ACOUSTIC WAVE PROPAGATION IN THE TRANSVERSELY ISOTROPIC MEDIA COMPRISING AN HYDROCARBON RESERVOIR

Abstract

A computer implemented method including a numerical model of a region of the earth modeling the acoustic behavior of that region. The method implements an acoustic wave propagator allowing the simulation of the propagation of pure P-waves in transversal isotropic media. The propagator can be applied to applications such as seismic forward modeling, reverse time migration and other two-way wave-equation based applications.

Claims

1. A method of operating a data-processing system for the simulation of an acoustic wave propagation in a transversely isotropic media comprising an hydrocarbon reservoir, the method comprising the steps of: a) determining a domain (Ω) comprising the reservoir, wherein said domain (Ω) comprises subsurface points and surface points (δ.sub.1Ω); b) populating the domain (Ω) with properties of an earth model, wherein said earth model at least comprises velocity field V.sub.z along a direction of a symmetry axis of the transversely isotropic media; c) determining boundary conditions for boundaries (δ.sub.2Ω) corresponding to the subsurface points and boundary conditions (δ.sub.1Ω) comprising one or more predetermined acoustic sources emitting at certain frequency located at the surface (δΩ) of the domain (Ω), and initial conditions; d) automatically solving, via the data-processing system, forward in time the earth model in the domain by using the wavefield propagator defined by: $\frac{1}{V_z^2}.Math.\frac{d^2.Math.u}{{\mathrm{dt}}^2}=\left(1+2.Math.\mathrm{.Math.}\right).Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}+\frac{{\partial }^2.Math.u}{\partial y^2}\right)+\frac{{\partial }^2.Math.u}{\partial z^2}+\left(\mathrm{.Math.}-\delta \right).Math.\frac{\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right).Math.\frac{{\partial }^2.Math.u}{\partial z^2}+{\left(\frac{\partial u}{\partial z}\right)}^2.Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}+\frac{{\partial }^2.Math.u}{\partial y^2}\right)}{\left(1+2.Math.\mathrm{.Math.}\right).Math.\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right)+{\left(\frac{\partial u}{\partial z}\right)}^2}$ wherein the u is the wavefield, ε and δ are the anisotropic parameters, (x, y) are the spatial coordinates at the surface (δ.sub.1Ω) of the domain (Ω) and z de spatial coordinate in the depth direction and t the time coordinate; and, e) making available the wavefield determined on the domain (Ω).

2. The method according to claim 1, wherein a plurality of receivers are located in the domain (Ω) and, after step d) the values calculated at the receiver locations are provided.

3. The method according to claim 1, wherein the boundary conditions for the boundaries corresponding to the subsurface points are absorbing free conditions.

4. The method according to claim 1, wherein: a) the initial conditions comprises a plurality of receivers located at predetermined locations of the surface of the domain being simulated; b) the wavefield generated by the acoustic sources are being interpolated and gathered in each receiver; and, c) the output of the receivers are made available.

5. A method of operating a data-processing system for a reverse time migration simulation, the method comprising the steps of: i. implementing the method according to claim 1 wherein: step b) is carried out by providing a recorded seismic field data and with a given estimation of the earth model in the domain (Ω); step c) is carried out by automatically determining, via the data-processing system, the acoustic source from the recorded seismic field data; step e) provides a time migrated seismic source wavefield; ii. implementing the method according to claim 1 wherein: step b) is determined by using the same estimated earth model as that used in step i); step c) is carried out by automatically providing, via the data-processing system, the time reversed seismic data as boundary conditions or initial conditions; step e) provides a reconstructed receiver wavefield, iii. automatically generating, via the data-processing system, a cross-correlation image between the source wavefield of step ii) and the receiver wavefield of step iii) and provide an estimated structure of the earth model in the domain (Ω).

6. The method according to claim 5 wherein the cross-correlation image is used as the gradient of the velocity field being used to update the estimated earth model; and steps i)-iii) are iteratively applied until convergence of the new estimated earth structure in the domain (Ω).

7. The method according to claim 1, wherein the integration of the partial differential equations is by using a space discretization based method.

8. The method according to claim 1, wherein the reservoir is tilted and the equation of step d) is being rotated; being equation of step d) being formulated as: $.Math.\frac{1}{V^2}.Math.\frac{d^2.Math.u}{{\mathrm{dt}}^2}=\left[{\mathrm{pa}}_{\mathrm{xy}}.Math..Math.p_{\mathrm{az}}\right]\left[\begin{array}{cc}{\mathrm{xy}}_{11}& z_{13}\\ {\mathrm{xy}}_{21}& z_{23}\\ {\mathrm{xy}}_{31}& z_{33}\\ {\mathrm{xy}}_{41}& z_{43}\\ {\mathrm{xy}}_{51}& z_{53}\\ {\mathrm{xy}}_{61}& z_{63}\end{array}\right]\left[\begin{array}{c}\frac{{\partial }^2.Math.u}{\partial x^2}\\ \frac{{\partial }^2.Math.u}{\partial y^2}\\ \frac{{\partial }^2.Math.u}{\partial z^2}\\ \frac{{\partial }^2.Math.u}{\partial x.Math.\partial y}\\ \frac{{\partial }^2.Math.u}{\partial y.Math.\partial z}\\ \frac{{\partial }^2.Math.u}{\partial z.Math.\partial x}\end{array}\right].Math..Math.\mathrm{where}$ $.Math.{\mathrm{pa}}_{\mathrm{xy}}=\left(1+2.Math.\mathrm{.Math.}\right).Math.V^2+\left(\mathrm{.Math.}-\delta \right).Math.V^2.Math.\frac{{\left(\frac{\partial u}{\partial z}\right)}^2}{\left(1+2.Math.\mathrm{.Math.}\right).Math.\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right)+{\left(\frac{\partial u}{\partial z}\right)}^2};$ ${\mathrm{pa}}_{\mathrm{az}}=\left(1+2.Math.\mathrm{.Math.}\right).Math.V^2+\left(\mathrm{.Math.}-\delta \right).Math.V^2.Math.\frac{{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2}{\left(1+2.Math.\mathrm{.Math.}\right).Math.\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right)+{\left(\frac{\partial u}{\partial z}\right)}^2};$ $.Math.z_{13}={\cos }^2\left(\phi \right).Math.\left(1-{\cos }^2\left(\theta \right)\right);$ $.Math.{\mathrm{xy}}_{21}={\cos }^2\left(\theta \right)*\left(1-{\cos }^2\left(\phi \right)\right)+{\cos }^2\left(\phi \right);$ $.Math.z_{33}={\cos }^2\left(\theta \right);$ $.Math.z_{43}=2.Math.\left(1-{\cos }^2\left(\theta \right)\right).Math.\sin \left(\phi \right).Math.\cos \left(\phi \right);$ $.Math.{\mathrm{xy}}_{51}=2.Math..Math.\sin \left(\theta \right).Math.\cos \left(\theta \right).Math.\sin \left(\phi \right);$ $.Math.{\mathrm{xy}}_{61}=2.Math..Math.\sin \left(\theta \right).Math.\cos \left(\theta \right).Math.\cos \left(\phi \right);$ $.Math.{\mathrm{xy}}_{11}=1-z_{13};$ $.Math.z_{23}=1-{\mathrm{xy}}_{21};$ $.Math.{\mathrm{xy}}_{31}=1-z_{33};$ $.Math.{\mathrm{xy}}_{41}=1-z_{43};$ $.Math.z_{53}=1-{\mathrm{xy}}_{51};$ $.Math.z_{63}=1-{\mathrm{xy}}_{61}.$

9. The method according to claim 1, wherein the wavefield propagator of step d) is based on the elastic stiffness tensor element further comprising the orthorhombic parameters, wherein said propagator is defined by: $.Math.\frac{d^2.Math.u}{{\mathrm{dt}}^2}=N.Math.\frac{{\partial }^2.Math.u}{\partial x^2}+P.Math.\frac{{\partial }^2.Math.u}{\partial y^2}+Q.Math.\frac{{\partial }^2.Math.u}{\partial z^2}.Math..Math.\mathrm{where}$ $N=V_1^2\left(1+2.Math.{\eta }_1\right)+\frac{\frac{1}{2}.Math.\left(V_1^4.Math.{\gamma }^2.Math.{\xi }_1^2-V_1^2.Math.V_2^2.Math.{\xi }_1.Math.{\xi }_2\right).Math.{\left(\frac{\partial u}{\partial y}\right)}^2+V_1^2.Math.V_z^2.Math.{{\eta }_1\left(\frac{\partial u}{\partial z}\right)}^2}{\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2}+\frac{\frac{1}{3}.Math.{\mathrm{CC}\left(\frac{\partial u}{\partial y}\right)}^2.Math.{\left(\frac{\partial u}{\partial z}\right)}^2}{{\left(\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2\right)}^2}$ $P=V_2^2\left(1+2.Math.{\eta }_2\right)+\frac{\frac{1}{2}.Math.\left(V_1^4.Math.{\gamma }^2.Math.{\xi }_1^2-V_1^2.Math.V_2^2.Math.{\xi }_1.Math.{\xi }_2\right).Math.{\left(\frac{\partial u}{\partial x}\right)}^2+V_2^2.Math.V_z^2.Math.{{\eta }_2\left(\frac{\partial u}{\partial z}\right)}^2}{\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2}+\frac{\frac{1}{3}.Math.{\mathrm{CC}\left(\frac{\partial u}{\partial x}\right)}^2.Math.{\left(\frac{\partial u}{\partial z}\right)}^2}{{\left(\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2\right)}^2}$ $Q=V_z^2+\frac{\frac{1}{2}.Math.V_z^2\left(V_1^2.Math.{{\eta }_1\left(\frac{\partial u}{\partial x}\right)}^2+V_2^2.Math.{{\eta }_2\left(\frac{\partial u}{\partial z}\right)}^2\right)}{\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2}+\frac{\frac{1}{3}.Math.{\mathrm{CC}\left(\frac{\partial u}{\partial x}\right)}^2.Math.{\left(\frac{\partial u}{\partial y}\right)}^2}{{\left(\left(1+2.Math.{\eta }_1\right).Math.{V_1^2\left(\frac{\partial u}{\partial x}\right)}^2+\left(1+2.Math.{\eta }_2\right).Math.{V_2^2\left(\frac{\partial u}{\partial y}\right)}^2+{V_z^2\left(\frac{\partial u}{\partial z}\right)}^2\right)}^2}$ wherein, the “1” index represents the “x” coordinate and “2” the “y” coordinate; being the horizontal velocity in the “x” direction and the horizontal velocity in the “y” direction:
V.sub.1=v.sub.1√{square root over (1+2η.sub.1)} and
V.sub.2=v.sub.2√{square root over (1+2η.sub.2)} respectively, being ${\eta }_1=\frac{c_{11}\left(c_{33}-c_{55}\right)}{2.Math.c_{13}\left(c_{13}+2.Math.c_{55}\right)+2.Math.c_{33}.Math.c_{55}}-\frac{1}{2}$ ${\eta }_2=\frac{c_{22}\left(c_{33}-c_{44}\right)}{2.Math.c_{23}\left(c_{23}+2.Math.c_{44}\right)+2.Math.c_{33}.Math.c_{44}}-\frac{1}{2}$ and, the quantity
y=√{square root over (1+2δ)} being $\delta =\frac{{\left(c_{12}+c_{66}\right)}^2-{\left(c_{11}-c_{66}\right)}^2}{2.Math.c_{11}\left(c_{11}-c_{66}\right)}$ ${\xi }_1=1+2.Math.{\eta }_1$ ${\xi }_2=1+2.Math.{\eta }_2$ $\mathrm{and}$ $\mathrm{CC}=2.Math.V_1^3.Math.V_2.Math.V_z^2.Math.{\mathrm{\gamma \xi }}_1-V_1^4.Math.V_z^2.Math.{\gamma }^2.Math.{\xi }_1^2-V_1^2.Math.V_2^2.Math.V_z^2\left(1-4.Math.{\eta }_1.Math.{\eta }_2\right)$

10. A computer program product adapted to carry out a method according to claim 1.

11. An electronic device for data analysis, the device comprising: a processor; and a computer readable medium storing computer-executable instructions which, when executed by the processor, result in the method according to claim 1.

Description

DESCRIPTION OF THE DRAWINGS

[0049] These and other features and advantages of the invention will be seen more clearly from the following detailed description of a preferred embodiment provided only by way of illustrative and non-limiting example in reference to the attached drawings.

[0050] FIG. 1 This figure shows a data processing system for carrying out a method according to the invention.

[0051] FIG. 2 This figure shows the domain of a numerical model being simulated by an example of method according to the invention.

[0052] FIG. 3a This figure shows a wavefield being simulated according to a two-equations based propagator in an example of the invention showing the s-wave noise and instabilities.

[0053] FIG. 3b This figure shows the same example being simulated with a propagator according to step d) wherein neither S-wave noise nor instabilities are being generated.

[0054] FIG. 4a This figure shows a wavefield being simulated according to a two-equations based propagator in another example when a single shot generates a wavefield. In this figure the scale of the image is not clearly shown as it is not relevant for the purpose of the teachings disclosed.

[0055] FIG. 4b This figure shows the same example being simulated with a propagator according to step d) wherein neither S-wave noise nor instabilities are being generated.

[0056] FIG. 5 This figure is an scheme showing that the reverse time migration RTM and forward modeling share the same wavefield propagator.

DETAILED DESCRIPTION OF THE INVENTION

[0057] As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon.

[0058] Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

[0059] 30

[0060] A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device.

[0061] Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.

[0062] Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The program code may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider).

[0063] Aspects of the present invention are described below with reference to illustrations and/or diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each illustration can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

[0064] These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks.

[0065] The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

[0066] Turning now to the drawings and more particularly, FIG. 1 shows an example of a data processing system 100 for the simulation of the acoustic wave propagation in the transversely isotropic media, according to a preferred embodiment of the present invention.

[0067] The preferred data processing system 100 determines a wavefield in an efficient manner free of noise, in particular without shear wave content.

[0068] A preferred data processing system 100 includes one or more computers 102, 104, 106 (3 in this example), coupled together, e.g., wired or wirelessly over a network 108. The network 108 may be, for example, a local area network (LAN), the Internet, an intranet or a combination thereof. Typically, the computers 102, 104, 106 include one or more processors, e.g., central processing unit (CPU) 110, memory 112, local storage 114 and some form of input/output device 116 providing a user interface. The local storage 114 may store and/or include earth models, in particular velocity fields and acoustic impedance data, computed wavefields; and raw data obtained from experiments carried out in specifics fields; being accessible by the plurality of computers 102, 104, 106, processing data by using the claimed propagator.

[0069] FIG. 2 shows an embodiment of the invention wherein a region of the earth comprising an oil reservoir is being limited by a domain (Ω). The domain (Ω) is also the numerical domain where the equations are solved. Given an earth model, the domain (Ω) is populated with at least the vertical velocity field V.sub.z. Additional rock properties or velocity field values may be populated.

[0070] The domain (Ω) is limited by a first boundary (δ.sub.1Ω) comprising the points of the surface of the earth region, and a second boundary (δ.sub.2Ω) corresponding to the subsurface points of the boundary of the domain (Ω).

[0071] Once the boundary conditions and the initial condition are set, the wave propagation is determined by solving forward in time the earth model in the domain by using the wavefield propagator, that is,

[00008]$\frac{1}{V_z^2}.Math.\frac{d^2.Math.u}{{\mathrm{dt}}^2}=\left(1+2.Math.\mathrm{.Math.}\right).Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}+\frac{{\partial }^2.Math.u}{\partial y^2}\right)+\frac{{\partial }^2.Math.u}{\partial z^2}+\left(\mathrm{.Math.}-\delta \right).Math.\frac{\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right).Math.\frac{{\partial }^2.Math.u}{\partial z^2}+{\left(\frac{\partial u}{\partial z}\right)}^2.Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}+\frac{{\partial }^2.Math.u}{\partial y^2}\right)}{\left(1+2.Math.\mathrm{.Math.}\right).Math.\left({\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial y}\right)}^2\right)+{\left(\frac{\partial u}{\partial z}\right)}^2},$

[0072] This equation is solved only for one wavefield component variable u reducing the computational cost of solving two coupled equations.

[0073] FIG. 2 shows a 2D domain (Ω) wherein only the x horizontal coordinate and the depth z coordinate are relevant.

[0074] Al terms comprising the y variable are neglected, being the equation simplified as:

[00009]$\frac{1}{V_z^2}.Math.\frac{d^2.Math.u}{{\mathrm{dt}}^2}=\left(1+2.Math.\mathrm{.Math.}\right).Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}\right)+\frac{{\partial }^2.Math.u}{\partial z^2}+\left(\mathrm{.Math.}-\delta \right).Math.\frac{{\left(\frac{\partial u}{\partial x}\right)}^2.Math.\frac{{\partial }^2.Math.u}{\partial z^2}+{\left(\frac{\partial u}{\partial z}\right)}^2.Math.\left(\frac{{\partial }^2.Math.u}{\partial x^2}\right)}{\left(1+2.Math.\mathrm{.Math.}\right).Math.{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial z}\right)}^2}.$

[0075] For a subsurface point I(x, y) the two way propagation between a determined acoustic source (S.sub.4) is depicted by a dashed line.

[0076] The examples shown in pictures 3a, 3b, 4a and 4b has been integrated by using a finite-differences method by imposing a 10 hz of dominant frequency Ricker wavelet as source and wherein the size of the spatial grid is dx=dy=dz=25 m.

[0077] FIG. 3a shows a wavefield being simulated according to two-equations based propagator showing the s-wave noise and instabilities. The upper part of the wavefield has been magnified allowing to show the crimped region due to the noise generation. The arrows points to the locations where the shear-waves noise exists.

[0078] The lower part also shows a distorted area within the wave fronts that have been magnified in the lower part of FIG. 3a.

[0079] The largest amplitudes do not satisfy geometry spreading and transmission/reflection amplitude rules. This issue is due to the existence of much lower shear-waves leading to a wavefield which is not stable when using the same time step used to integrate numerically higher acoustic p-waves.

[0080] The same example has been simulated by using the equation according to step d) according to the invention. The same parts of the new wavefield are magnified on FIG. 3b showing that the same regions are noise free and without instabilities.

[0081] A new embodiment is being simulated by using a single shot in the middle of the surface. FIG. 4a shows a wavefield being simulated according to a two-equation based propagator showing a region identified with an oval where the S-wave noise generates a distorted image. The arrows points to the locations where the shear-waves noise exists.

[0082] FIG. 4b shows an RTM result for an impulse response with marine data. It shows the same example being simulated with a propagator according to step d) of the invention wherein neither S-wave noise nor instabilities are being generated. The plot does not have the shear-waves noise and appears very clean.

[0083] The method according to the invention may also be used for other applications such as reversed time migration of images. In this particular example the propagator used for the source wavelet of the forward modeling and the propagator used for the time reversed seismic common shot gather are the same, that identified in step d) of the present invention.

[0084] FIG. 5 shows the revers time migration workflow. The left part of the figure identifies the forward modeling wherein the source wavelet is determined in the domain (Ω), and the earth model provides the properties being populated in said domain (Ω).

[0085] The boundary and initial conditions allows to integrate the initial value problem stated by equation of step d), then obtaining as a result the source wavefield.

[0086] This process is identified as a particular embodiment of the method of the invention deeming that (step i): [0087] step b) is carried out by providing a recorded seismic field data and with a given estimation of the earth model in the domain (Ω); [0088] step c) is carried out by determining the acoustic source from the recorded seismic field data; [0089] step e) provides a time migrated seismic source wavefield.

[0090] The right part of figure identifies the time reversed seismic common shot gather which provides the boundary and initial conditions to be used when solving the wavefield propagator.

[0091] The result is a receiver wavefield as it is shown at the end of the right branch of the scheme.

[0092] This process is also identified as a particular embodiment of the method of the invention deeming that: [0093] step b) is determined by using the same estimated earth model as that used in step i); [0094] step c) is carried out by providing the time reversed seismic data as boundary conditions or initial conditions; [0095] step e) provides a reconstructed receiver wavefield.

[0096] A cross-correlation image between the source wavefield and the receiver wavefield provides an estimated earth structure in the domain (Ω), the RTM Image output.