SEISMIC WAVEFIELD MODELING HONORING AVO/AVA WITH APPLICATIONS TO FULL WAVEFORM INVERSION AND LEAST-SQUARES IMAGING

Abstract

A method for modelling and migrating seismic data, that includes using an acoustic wave equation and a spatial distribution of one or more earth-model parameters. The acoustic wave equation is modified by including at least one secondary source term, and based on a seismic acquisition configuration, either calculating the seismic signals that would be detected from the modelled wavefield or migrating observed seismic signals or migrating residual signals as part of an inversion.

Claims

1. A method for deducing subsurface physical property parameters for seismic exploration, comprising, in a computer: calculating seismic signals that would be observed based on an arrangement of at least one seismic source and at least one seismic receiver and an initial model of Earth's subsurface, the initial model comprising spatial distribution of at least one parameter related to propagation of elastic waves, the calculating comprising entering the initial model into an acoustic wave equation modified by at least one secondary source term, the at least one secondary source term enabling estimation of amplitude versus offset and/or angle (AVO/A) phenomena in the calculated seismic signals, the acoustic wave equation capable of the estimation only when modified by the at least one secondary source term; comparing the calculated seismic signals with measured seismic signals obtained using the arrangement; determining differences between the calculated seismic signals and the measured seismic signals; and using an inversion method to adjust the initial model and repeating the calculating, comparing and determining to reduce the determined differences.

2. The method of claim 1 wherein the adjusting comprises changing at least one of a value of the at least one parameter and the spatial distribution.

3. The method of claim 2 wherein the at least one secondary source term represents reflectivity.

4. The method of claim 3 wherein the reflectivity results from changes in at the least one parameter.

5. The method of claim 4 wherein the at least one parameter comprises compressional (P) wave velocity.

6. The method of claim 4 wherein the at least one parameter comprises shear (S) wave velocity.

7. The method of claim 4 wherein the at least one parameter comprises density.

8. The method of claim 4 wherein the at least one parameter comprises compressional (P) wave impedance.

9. The method of claim 4 wherein the at least one parameter comprises shear (S) wave impedance.

10. The method of claim 4 wherein the at least one parameter comprises elastic impedance.

11. The method of claim 4 wherein the at least one parameter comprises at least one component of an elastic tensor.

12. The method of claim 4 wherein the at least one parameter comprises attenuation.

13. The method of claim 4 further comprising producing a reflectivity image from the adjusted initial earth model.

14. The method of claim 4 wherein the reflectivity is expressed as the grad of a spatially-variable scalar parameter field.

15. The method of claim 13 further comprising selectively modifying the at least one parameter and/or a strength of the at least one secondary source term in order to control production of diving wave-paths (bananas) and reflection wave-paths (rabbit ears) as substantially separate terms in a gradient of an objective function of a full waveform inversion.

16. The method of claim 13 further comprising selectively modifying the at least parameter and/or a strength of the at least one secondary source term in order to substantially produce only diving wave-paths (bananas) in a gradient of an objective function of a full waveform inversion.

17. The method of claim 13 further comprising selectively modifying the at least one parameter and/or a strength of the at least one secondary source term in order to substantially produce only reflection wave-paths (rabbit ears) in a gradient of an objective function of a full waveform inversion.

18. The method of claim 1 wherein the at least one secondary source term effectively reflects or scatters an incident wavefield with directional variations.

19. The method of claim 1 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A intercept attribute.

20. The method of claim 1 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A gradient attribute.

21. The method of claim 1 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A curvature attribute.

22. The method of claim 1 further comprising injecting resulting secondary source terms into a different wavefield from that which their values depend on.

23. The method of claim 1 further comprising Born modelling.

24. The method of claim 1 wherein a termination criterion of the inversion method is given by a predetermined number of repetitions.

25. The method of claim 1 wherein a termination criterion of the inversion is given by a predetermined similarity between the calculated and measured signals.

26. A computer program stored in a non-transitory computer readable medium, the program having logic operable to cause a programmable computer to perform a process for deducing subsurface physical property parameters for seismic exploration, the process comprising: calculating seismic signals that would be observed based on an arrangement of at least one seismic source and at least one seismic receiver and an initial model of Earth's subsurface, the initial model comprising spatial distribution of at least one parameter related to propagation of elastic waves, the calculating comprising entering the initial model into an acoustic wave equation modified by at least one secondary source term, the at least one secondary source term enabling estimation of amplitude versus offset and/or angle (AVO/A) phenomena in the calculated seismic signals, the acoustic wave equation capable of the estimation only when modified by the at least one secondary source; comparing the calculated seismic signals with measured seismic signals obtained using the arrangement; determining differences between the calculated seismic signals and the measured seismic signals; and using an inversion method to adjust the initial model and repeating the calculating, comparing and determining to reduce the determined differences.

27. The computer program of claim 26 wherein the adjusting comprises changing at least one of a value of the at least one parameter and the spatial distribution.

28. The computer program of claim 27 wherein the at least one secondary source term represents reflectivity.

29. The computer program of claim 28 wherein the reflectivity results from changes in the at least one parameter.

30. The computer program of claim 29 wherein the at least one parameter comprises compressional (P) wave velocity.

31. The computer program of claim 29 wherein the at least one parameter comprises shear (S) wave velocity.

32. The computer program of claim 29 wherein the at least one parameter comprises density.

33. The computer program of claim 29 wherein the at least one parameter comprises compressional (P) wave impedance.

34. The computer program of claim 29 wherein the at least one parameter comprises shear (S) wave impedance.

35. The computer program of claim 29 wherein the at least one parameter comprises elastic impedance.

36. The computer program of claim 29 wherein the at least one parameter comprises at least one component of an elastic tensor.

37. The computer program of claim 29 wherein the at least one parameter comprises attenuation.

38. The computer program of claim 29 wherein the logic is further operable to cause the computer to perform producing a reflectivity image from the adjusted initial earth model.

39. The method of claim 38 wherein the reflectivity is expressed as the grad of a spatially-variable scalar parameter field.

40. The computer program of claim 39 wherein the logic is further operable to cause the computer to perform selectively modifying the at least one parameter and/or a strength of the at least one secondary source term in order to control production of diving wave-paths (bananas) and reflection wave-paths (rabbit ears) as substantially separate terms in a gradient of an objective function of a full waveform inversion.

41. The computer program of claim 39 wherein the logic is further operable to cause the computer to perform selectively modifying the at least parameter and/or a strength of the at least one secondary source term in order to substantially produce only diving wave-paths (bananas) in a gradient of an objective function of a full waveform inversion.

42. The computer program of claim 39 wherein the logic is further operable to cause the computer to perform selectively modifying the at least one parameter and/or a strength of the at least one secondary source term in order to substantially produce only reflection wave-paths (rabbit ears) in a gradient of an objective function of a full waveform inversion.

43. The computer program of claim 26 wherein the at least one secondary source term effectively reflects or scatters an incident wavefield with directional variations.

44. The computer program of claim 26 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A intercept attribute.

45. The computer program of claim 26 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A gradient attribute.

46. The computer program of claim 26 wherein the at least one secondary source term represents reflectivity corresponding to an AVO/A curvature attribute.

47. The computer program of claim 26 further comprising injecting resulting secondary source terms into a different wavefield from that which their values depend on.

48. The computer program of claim 26 wherein the logic is further operable to cause the computer to perform Born modelling.

49. The computer program of claim 26 wherein a termination criterion of the inversion is given by a predetermined number of repetitions.

50. The computer program of claim 26 wherein a termination criterion of the inversion is given by a predetermined similarity between the calculated and measured signals.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0076] FIG. 1 shows parts of a full waveform inversion (FWI) gradient: The diving wave-path (banana), the reflection wave-path (rabbit ears) and the reflection.

[0077] FIG. 2 shows that at a boundary between two media, an incident P-wave will be partially reflected and transmitted and will also undergo partial mode conversion into S-waves.

[0078] FIG. 3 shows an example embodiment of acquiring seismic signals that may be used in accordance with the present disclosure.

[0079] FIG. 4 shows the virtual secondary source wavefield reflecting off a boundary. The orientation of the incident wavefront normal, and the normal to the boundary, may be used to determine the sine of the angle of incidence using the cross product.

[0080] FIG. 5 shows an example embodiment of a computer system that may be used in accordance with the present disclosure.

DETAILED DESCRIPTION

[0081] Methods according to the present disclosure may be used in connection with seismic data acquired using techniques known in the art. Such techniques comprise deploying one or more seismic energy sources, e.g., air guns or arrays of such guns, vibrators or arrays of vibrators or explosives, at selected locations, or such techniques may use naturally occurring sources. Arrays of seismic sensors or receivers are deployed in selected arrangements and detect seismic energy that ultimately originates from the source(s). An example embodiment of acquiring seismic signals (data) is shown in FIG. 3.

[0082] FIG. 3 shows a vertical section view of an ocean bottom cable (OBC) seismic survey being conducted using, for example, two different “source” vessels for towing seismic energy sources. In some embodiments, only one source vessel may be used. According to the present disclosure, any number of source vessels may be used and the following description is not intended to limit the scope of the present disclosure. The source vessels move along the surface 16A of a body of water 16 such as a lake or the ocean. In the present example, a vessel referred to as a “primary source vessel” 10 may include equipment, shown generally at 14, that comprises components or subsystems (none of which is shown separately) for navigation of the primary source vessel 10, for actuation of seismic energy sources and for retrieving and processing seismic signal recordings. The primary source vessel 10 is shown towing two, spaced apart seismic energy sources 18, 18A.

[0083] The equipment 14 on the primary source vessel 10 may be in signal communication with corresponding equipment 13 (including similar components to the equipment on the primary source vessel 10) disposed on a vessel referred to as a “secondary source vessel” 12. The secondary source vessel 12 in the present example also tows spaced apart seismic energy sources 20, 20A near the water surface 16A. In the present example, the equipment 14 on the primary source vessel 10 may, for example, send a control signal to the corresponding equipment 13 on the secondary source vessel 12, such as by radio telemetry, to indicate the time of actuation (firing) of each of the sources 18, 18A towed by the primary source vessel 10. The corresponding equipment 13 may, in response to such signal, actuate the seismic energy sources 20, 20A towed by the secondary source vessel 12.

[0084] The seismic energy sources 18, 18A, 20, 20A may be air guns, water guns, marine vibrators, or arrays of such devices. The seismic energy sources are shown as discrete devices in FIG. 3 to illustrate the general principle of seismic signal acquisition. The type of and number of seismic energy sources that can be used in any embodiment is not intended to limit the scope of the disclosure.

[0085] In FIG. 3, an OBC 22 is deployed on the bottom 16B of the water 16 such that spaced apart seismic receiver modules 24 are disposed on the water bottom 16B in a preselected pattern. The receiver modules 24 may include a pressure or pressure time gradient responsive seismic sensor, and one or more seismic particle motion sensors, for example, one-component or three-component geophones, or one- or three-component accelerometers (none of the sensors are shown separately). The type of and the number of seismic sensors in each module 24 is not intended to limit the scope of the disclosure. The seismic sensors in each module 24 generate electrical and/or optical signals (depending on the sensor type) in response to, in particular, detected seismic energy resulting from actuations of the seismic energy sources 18, 18A, 20, 20A. In some embodiments, the sensors may comprise pressure or pressure time derivative sensors such as hydrophones, and one or more particle motion responsive sensors such as geophones or accelerometers. In some embodiments, such particle motion responsive sensors may be oriented so as to be sensitive principally (ignoring any cross-component coupling for purposes of explanation) to vertical particle motion. The signals generated by the various sensors may be conducted to a device near the water surface 16A such as a recording buoy 23, which may include a data recorder (not shown separately) for storing the signals for later retrieval and processing by the equipment 14 on the primary source vessel 10, or other processing equipment to be described further below. The data storage functions performed by the recording buoy 23 may be performed by different types of equipment, such as a data storage unit on a recording vessel (not shown) or a recording module (not shown) deployed on the water bottom 16B, e.g., proximate each sensor module 24, or even on the primary or secondary source vessels. Accordingly, the disclosure is not limited in scope to use with a recording buoy or any other specific recording device(s).

[0086] Although the description of acquiring signals explained with reference to FIG. 3 is for sensors deployed on the water bottom, it will be appreciated that it is possible to obtain corresponding measurements at any selected depth in the water, using, for example, seismic sensors disposed in a towed cable such as described in U.S. Pat. No. 7,239,577 issued to Tenghamn, et al. Further, the present disclosure is not limited to use with seismic signals acquired in or on the surface of a body of water; methods described herein may also be used with seismic signals acquired on land.

[0087] Seismic energy detected by the sensors or receivers generates signals that are recorded for later processing, including by methods according to this disclosure. Such methods may comprise generating a model of materials (e.g., formations) in the Earth, more specifically, the values of one or more parameters related to the physical properties of the materials and their spatial distribution in the Earth. The model may be used to calculate modeled or synthetic seismic signals, that is, signals that would be detected by each of the sensors or receivers such as in FIG. 3, in response to the energy imparted by the sources if the Earth had properties as defined in the model. Differences between the calculated, modeled or synthetic seismic signals, and actual measured (observed) seismic signals may be used to adjust the model and subsequently to repeat performing calculation of the modeled seismic signals. Adjusting the model may comprise changing a value of one or more parameters and/or the spatial distribution of one or more parameters. In this sense, “parameters” may mean any physical property of the formations that will affect propagation of the seismic wavefield, including without limitation, compressional wave (P) velocity, shear wave (S) velocity, density, and fluid content in porous media. The forgoing actions may be repeated until differences between the modeled seismic signals and the observed (measured) seismic signals are reduced to below a predetermined difference limit, or minimized, as explained in the Background section herein. In some embodiments, the foregoing process, which may be referred to as an “inversion process”, may be repeated for a predetermined number of cycles each comprising model adjustment, recalculation of the modeled seismic signals and comparison to the observed seismic signals.

[0088] Reflections of the seismic wavefield from impedance boundaries in real (measured) seismic data exhibit amplitude vs. offset and angle (AVO/A) phenomena, i.e., the amplitudes R of the reflection events in the seismic data vary with the reflection angle (θ). Amplitude variation for angles to within approximately 10° of the critical angle as explained, for example, in Shuey (1985) and may be represented by the expression:

R(θ)=I+G sin.sup.2 θ+C(tan.sup.2 θ−sin.sup.2 θ)  (3)

[0089] In Eq. (3), I, G and C are the AVO/A intercept, gradient and curvature terms, respectively. These terms are functions of the physical properties of the media (i.e., the subsurface formations). For example, the intercept term may relate to certain physical properties of formations by the expression:

[00003]$\begin{array}{cc}I\approx \frac{1}{2}\left(\frac{\Delta c}{c}+\frac{\mathrm{\Delta \rho }}{\rho }\right)& \left(3a\right)\end{array}$

wherein c represents the P-wave (compressional wave) velocity and ρ represents bulk density. Seismic reflections modeled with an acoustic wave equation, as explained in the Background section herein, do not correctly exhibit these AVO/A phenomena. For more modest reflection angles, (e.g., 0°<θ<30°), which apply in many practical cases, the curvature term may be neglected. However, in other cases, accurate modelling of seismic data would be enhanced by suitably representing AVO/A phenomena.

[0090] An explanation of a method according to the present disclosure may begin with the variable density acoustic wave equation and may then show how dynamic behavior that is not captured by the P-wave velocity field alone can be modeled by the introduction of additional terms in the variable density acoustic wave equation that act as secondary sources. One of these terms will produce AVO/A intercept behavior that is not captured by the P-wave velocity field alone and a second term will capture the AVO/A gradient behavior. It follows that a third term to capture AVO/A curvature behavior may be readily included, although the present disclosure may omit certain details in the interest of brevity of the explanation. Once it is established how to formulate the foregoing reflectivity terms in a modified isotropic variable density acoustic wave equation (as explained in the Background section herein, an isotropic wave equation does not use shear velocity as a parameter), it then becomes a relatively simple matter to selectively supplement the background velocity field in full waveform inversion (FWI) to generate bananas and/or rabbit ears, e.g., as explained with reference to FIG. 1. For illustrative purposes, the approach may be described using an isotropic wave equation, however, this method readily extends to the anisotropic case. Note that the anisotropic acoustic case will naturally require use the shear velocity as a parameter but will still inadequately describe AVO/AVA without the method herein.

[0091] 1. The Intercept Term

[0092] Beginning with the isotropic variable-density, acoustic, inhomogeneous wave-equation:

[00004]$\begin{array}{cc}\frac{1}{\rho c^2}\frac{{\partial }^{2}u}{\partial t^2}=\nabla .Math.\left(\frac{1}{\rho }\nabla u\right)+s& \left(4\right)\end{array}$

in which, the scalar wavefield, the source term, the bulk density, the P-wave velocity and time are denoted by u, s, ρ, c and t respectively, the above wave equation may be rearranged, and making use of the product rule it may be observed that:

[00005]$\begin{array}{cc}\begin{array}{c}\frac{1}{c^2}\frac{{\partial }_u^2}{\partial t^2}=\rho \nabla .Math.\left(\frac{1}{\rho }\nabla u\right)+\rho s\\ ={\nabla }^{2}u+\rho \nabla \frac{1}{\rho }.Math.\nabla u+\rho s\\ ={\nabla }^{2}u-\nabla \left(\log \rho \right).Math.\nabla u+\rho s\\ ={\nabla }^{2}u-\frac{1}{\rho }\nabla \rho .Math.\nabla u+\rho s\end{array}& \left(5\right)\end{array}$

[0093] Bulk density contrasts can be modeled using an additional secondary source term, −∇(log ρ).Math.∇u, (or alternatively,

[00006]$-\frac{1}{\rho }\nabla \rho .Math.\nabla u)$

in a constant density wave equation. The second interpretation shows a direct correspondence with the density contribution to the intercept term in the AVO/A equation, as shown by the expression:

[00007]$\begin{array}{cc}R\left(\theta =0\right)=\frac{1}{2}\left(\frac{\Delta c}{c}+\frac{\mathrm{\Delta \rho }}{\rho }\right)& \left(6\right)\end{array}$

[0094] As a result it may be observed that this modified acoustic wave equation may be written in terms of a dot product between a spatially-varying vector parameter field r.sub.1 and with the gradient of the scalar wavefield, ∇u, to produce intercept-like scattering as in the expression:

[00008]$\begin{array}{cc}\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-r_I.Math.\nabla u+\rho s& \left(7\right)\end{array}$

in which the new term, r.sub.1.Math.∇u has been emphasized. It is important to note that although the above Eq. (7) was developed to model scattering due to density variations, any other suitable parameter field may be chosen to produce scattering behavior in modeled seismic data. For example, in the case that the background velocity field, c, is smooth, the additional term in Eq. (7) could be used to generate intercept-like scattering. In such case, the intercept parameter field r.sub.1 can take on a role similar to ∇(log ρ+log c)=∇ log ρc=∇ log Z. Note here that Z is the acoustic impedance. Similar reasoning may be applied to any earth formation parameters that generate scattering due to abrupt changes in value (e.g., change with respect to spatial position), for example P-wave velocity, shear-wave velocity, bulk density, P-wave impedance and shear-wave impedance.

[0095] It follows that the introduction of other secondary source terms could serve other purposes, for example, to introduce AVO/A gradient-like and/or AVO/A curvature-like scattering behavior. However, for these a directional aspect is required in the secondary source term.

[0096] 2. The Gradient Term

[0097] For the AVO/A gradient term, one may introduce an additional secondary source term to the modified acoustic wave equation, e.g., Eq. (7), that has sin.sup.2 θ directivity. This additional secondary source term depends upon a gradient-related vector parameter field r.sub.G. The modified wave equation, Eq. (7) then becomes:

[00009]$\begin{array}{cc}\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-{\stackrel{.fwdarw.}{r}}_I.Math.\nabla u-{\stackrel{.fwdarw.}{r}}_G.Math.{\sin }^2\theta \nabla u+\rho s& \left(8\right)\end{array}$

in which the new gradient term has been emphasized. Unless there is any ambiguity, in what follows one may omit the subscript G. A useful mathematical approach may be to produce the directivity behavior using a finite difference stencil.

[0098] The sine of the angle between the incident wavefront normal and the normal to a parameter boundary, e.g., as shown in FIG. 2, can be determined from the vector cross-product, ∇r×∇u, since |∇r×∇u|=|∇r||∇u| sin θ.

[0099] FIG. 4 shows a virtual secondary source wavefield reflecting off a boundary. The orientation of the incident wavefront normal, and the normal to the boundary, ∇r may be used to determine the sine of the angle of incidence using the cross product, |∇r×∇u|=|∇r||∇u| sin θ.

[0100] Consider the gradient of a scalar wavefield u,

[00010]$\begin{array}{cc}\nabla u=\left(\frac{\partial }{\partial x}\hat{x}+\frac{\partial }{\partial y}\hat{y}+\frac{\partial }{\partial z}\hat{z}\right)u& \left(9\right)\end{array}$

and assume that u is the plane monochromatic wave represented by the expression,

u=e.sup.−i(ωt-k.sup.x.sup.x-k.sup.y.sup.y-k.sup.z.sup.z)  (10)

then the gradient of the scalar wavefield is given by the expression:

∇u=i(k.sub.x{circumflex over (x)}+k.sub.yŷ+k.sub.z{circumflex over (z)})u  (11)

[0101] The normal to the reflecting boundary can be estimated by taking the gradient of a suitable physical field, r (e.g., density, S-wave velocity or P-wave velocity):

[00011]$\begin{array}{cc}\stackrel{.fwdarw.}{r}=\nabla r=\left(\frac{\partial }{\partial x}\hat{x}+\frac{\partial }{\partial y}\hat{y}+\frac{\partial }{\partial y}\hat{z}\right)r& \left(12\right)\end{array}$

whose values are different in different physical media. This is typically the field that will be used to induce AVO/A effects into modeled reflections. For convenience, the notation, r.sub.x=∂r/∂x may be adopted, so that:

{right arrow over (r)}=∇r=r.sub.x{circumflex over (x)}+r.sub.yŷ+r.sub.z{circumflex over (z)}  (13)

[0102] The cross product of the wavefront normal and the normal to the reflecting boundary may be given by the expression:

|∇u×{right arrow over (r)}|=|∇u||{right arrow over (r)}| sin θ  (14)

and the directivity required may be given by the expression:

[00012]$\begin{array}{cc}\begin{array}{c}{\sin }^2\theta =\frac{{.Math.\nabla u\times \stackrel{.fwdarw.}{r}.Math.}^2}{{.Math.\nabla u.Math.}^2{.Math.\stackrel{.fwdarw.}{r}.Math.}^2}\\ =\frac{\begin{array}{c}{k}_x^2\left(r_y^2+r_z^2\right)+k_y^2\left(r_x^2+r_z^2\right)+k_z^2\left(r_x^2+r_y^2\right)-\\ 2k_x{k}_y{r}_x{r}_y-2k_x{k}_z{r}_x{r}_z-2k_y{k}_z{r}_y{r}_z\end{array}}{\left(k_x^2+k_y^2+k_z^2\right)\left(r_x^2+r_y^2+r_z^2\right)}\end{array}& \left(15\right)\end{array}$

[0103] This means that the directivity of the secondary source term with sin.sup.2 θ behavior can be constructed by applying an operator that is a combination of the components of the parameter field to the gradient of the wavefield, ∇u:

[00013]$\begin{array}{cc}\begin{array}{c}{\sin }^2\theta =\frac{{.Math.\nabla u\times \stackrel{.fwdarw.}{r}.Math.}^2}{{.Math.\nabla u.Math.}^2{.Math.\stackrel{.fwdarw.}{r}.Math.}^2}\\ =\frac{\begin{array}{c}{k}_x^2\left(r_y^2+r_z^2\right)+k_y^2\left(r_x^2+r_z^2\right)+k_z^2\left(r_x^2+r_y^2\right)-\\ 2k_x{k}_y{r}_x{r}_y-2k_x{k}_z{r}_x{r}_z-2k_y{k}_z{r}_y{r}_z\end{array}}{\left(k_x^2+k_y^2+k_z^2\right)\left(r_x^2+r_y^2+r_z^2\right)}\end{array}& \left(15\right)\end{array}$

[0104] Define:

[00014]$\begin{array}{cc}v=\frac{u}{k_x^2+k_y^2+k_z^2}=\frac{u}{k^2},& \left(17\right)\end{array}$

so that ∇v=∇u/k.sup.2 and then:

[00015]$\begin{array}{cc}{\sin }^2\theta \nabla u={\sin }^2\theta k^2\nabla v=\frac{\begin{array}{c}{k}_x^2\left(r_y^2+r_z^2\right)+k_y^2\left(r_x^2+r_z^2\right)+k_z^2\left(r_x^2+r_y^2\right)-\\ 2k_x{k}_y{r}_x{r}_y-2k_x{k}_z{r}_x{r}_z-2k_y{k}_z{r}_y{r}_z\end{array}}{r_x^2+r_y^2+r_z^2}\nabla v& \left(18\right)\end{array}$

[0105] It remains to complete the new AVO/A term by applying the negative gradient of the parameter field:

[00016]$\begin{array}{cc}-{\stackrel{.fwdarw.}{r}}_G.Math.{\sin }^2\theta k^2-\nabla v=-{\stackrel{.fwdarw.}{r}}_G.Math.\frac{\begin{array}{c}{k}_x^2\left(r_y^2+r_z^2\right)+k_y^2\left(r_x^2+r_z^2\right)+k_z^2\left(r_x^2+r_y^2\right)-\\ 2k_x{k}_y{r}_x{r}_y-2k_x{k}_z{r}_x{r}_z-2k_y{k}_z{r}_y{r}_z\end{array}}{r_x^2+r_y^2+r_z^2}\nabla v& \left(19\right)\end{array}$

[0106] Noting the Fourier duals k.sub.x.sup.2=−∂.sup.2/∂x.sup.2 and k.sub.xk.sub.y=−∂.sup.2/∂x∂y, re-introducing the partial differential operators in the numerator provides the expression:

[00017]$\begin{array}{cc}-{\stackrel{.fwdarw.}{r}}_G.Math.{\sin }^2\theta k^2\nabla v=.Math.{\stackrel{.fwdarw.}{r}}_G.Math.\frac{\begin{array}{c}\left(r_y^2+r_z^2\right)\frac{{\partial }^2}{\partial x^2}+\left(r_x^2+r_z^2\right)\frac{{\partial }^2}{\partial y^2}+\left(r_x^2+r_y^2\right)\frac{{\partial }^2}{\partial z^2}-\\ 2r_x{r}_y\frac{{\partial }^2}{\partial x\partial y}-2r_x{r}_z\frac{{\partial }^2}{\partial x\partial z}-2r_y{r}_z\frac{{\partial }^2}{\partial y\partial z}\end{array}}{r_x^2+r_y^2+r_z^2}\nabla v& \left(20\right)\end{array}$

[0107] It should be noted that the product with ∇v generates triple partial derivatives of the wavefield.

[0108] There remains the calculation of the auxiliary field, v. One approach is to assume that the isotropic dispersion relation, ω.sup.2/c.sup.2=k.sub.x.sup.2+k.sub.y.sup.2+k.sub.z.sup.2, holds, in which case it is possible to re-arrange Eq. (8) so that ω.sup.2v=c.sup.2u. Recognizing the Fourier transform pair, ∂.sup.2/∂t.sup.2⇔−ω.sup.2, one may write that:

[00018]$\begin{array}{cc}\frac{{\partial }^2}{\partial t^2}v=-c^{2}u& \left(21\right)\end{array}$

which admits a solution for v by double time integration of −c.sup.2u.

[0109] Having shown herein how the intercept and gradient terms can be introduced as secondary sources if it is denoted that a scatterer generating parameter field, r.sub.j that has a directivity, described by an operator, D.sub.j, in principle it is possible to generalize the modified wave-equation for n secondary source terms as follows:

[00019]$\begin{array}{cc}\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-\underset{j=1}{\overset{n}{.Math.}}{\stackrel{.fwdarw.}{r}}_j.Math.D_j\nabla v_j+\rho s,& \left(22\right)\end{array}$

in which it is assumed that in the most general case n auxiliary fields are required. A summary of the different forms of the wave equation discussed can be seen in Table 2 below

TABLE-US-00002 TABLE 2 Forms of the Acoustic Wave Equation for Modelling AVO/A Phenomena Secondary source Form of wave physical phenomena equation Equation represented Isotropic acoustic wave equation with constant density [00020]$\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u+\rho s$ N/A Isotropic acoustic wave equation with variable density [00021]$\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-\nabla \left(\log \rho \right).Math.\nabla u+\rho s$ Bulk density contrast Isotropic acoustic augmented wave equation 1 [00022]$\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-r_I.Math.\nabla u+\rho s$ Change in earth formation parameter (e.g. AVO/A intercept) Isotropic acoustic augmented wave equation 2 [00023]$\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-r_I.Math.\nabla u-r_G.Math.{\sin }^2\theta \nabla u+\rho s$ Change in earth formation parameter (e.g. AVO/A intercept) and AVO/A gradient Isotropic acoustic augmented wave equation 3 [00024]$\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-\underset{j=1}{\overset{n}{.Math.}}{r}_j.Math.D_j\nabla v_j+\rho s$ This is the sum of many secondary sources representing, for example, intercept, gradient and curvature

[0110] Summary of the Method

[0111] Herein is described a method in which the acoustic wave equation (see Eq. 5 for the isotropic case), that is, a wave equation that uses compressional velocity as a velocity and bulk density as parameters, may be augmented with secondary source terms to represent scattering of various types, in particular so that AVO/A behavior may be represented in synthetic seismic signals or data. Because the secondary source terms can be designed to have appropriate directivity, they may be used to model aspects of wavefield behavior normally only possible with an elastic wave equation. As an example, it has been shown how the augmented acoustic wave equation given by Eq. (8), along with an auxiliary equation, e.g., Eq. (22):

[00025]$$\begin{gathered} \begin{array}{cc}\frac{1}{c^2}\frac{{\partial }^{2}u}{\partial t^2}={\nabla }^{2}u-{\stackrel{.fwdarw.}{r}}_I.Math.\nabla u-{\stackrel{.fwdarw.}{r}}_G.Math.{\sin }^2\theta k^2\nabla u+\rho s \\ \frac{{\partial }^{2}v}{\partial t^2}=-c^{2}u& \left(23\right)\end{array} \end{gathered}$$

may be used to model the AVO/A behavior suggested by the 2-term approximation described in Shuey (1985). The intercept term may be given by −{right arrow over (r)}.sub.1.Math.∇u and the AVO/A gradient term may be given by −{right arrow over (r)}.sub.G.Math.k.sup.2 sin.sup.2 θ∇ν. This augmented acoustic wave equation requires an auxiliary parameter field, v, which is the solution of the auxiliary equation, Eq. (23). Using the augmented acoustic wave equation to synthesize seismic data or its adjoint to backproject seismic data into the image space can yield results more similar to those that may be obtained using the full elastic wave equation but at much reduced computational cost.

[0112] Since it has been shown how to determine the sine of the incident angle using the vector cross-product, it will be appreciated that it is possible to derive the cosine using the vector dot-product. Consequently, it follows that it is possible to include a secondary source term representing AVO/A curvature. The directivity of these secondary source terms is not limited to those described. Indeed, a wide range of other forms could also be considered for a variety of purposes.

[0113] This type of augmented wave equation may be used to model seismic wavefields, produce seismic images and perform full waveform inversion (FWI) much more efficiently than using an elastic wave equation. Using this approach it may be possible to produce estimates of a range of parameters such as intercept, gradient, curvature, density, attenuation and velocity-related fields.

[0114] By using such reflectivity terms in a FWI setting, the reflectivity terms may be used selectively to supplement the velocity field in order to control the generation of reflections and reflection wave-paths in the FWI gradient. Since this technique permits a separation of the propagation and scattering terms, not only does it naturally lend itself to Born modelling, but it has also been reported to make, “the inversion problem significantly more linear” (Verschuur et al., 2016).

[0115] The foregoing process may be implemented on any form of computer or computer system programmed to perform the actions described. FIG. 5 shows an example computing system 100 in accordance with some embodiments. The computing system 100 may be an individual computer system 101A or an arrangement of distributed computer systems. The individual computer system 101A may include one or more analysis modules 102 that may be configured to perform various tasks according to some embodiments, such as the tasks explained with reference to FIG. 5. To perform these various tasks, the analysis module 102 may operate independently or in coordination with one or more processors 104, which may be connected to one or more storage media 106. A display device 105 such as a graphic user interface of any known type may be in signal communication with the processor 104 to enable user entry of commands and/or data and to display results of execution of a set of instructions according to the present disclosure.

[0116] The processor(s) 104 may also be connected to a network interface 108 to allow the individual computer system 101A to communicate over a data network 110 with one or more additional individual computer systems and/or computing systems, such as 101B, 101C, and/or 101D (note that computer systems 101B, 101C and/or 101D may or may not share the same architecture as computer system 101A, and may be located in different physical locations, for example, computer systems 101A and 101B may be at a well drilling location, while in communication with one or more computer systems such as 101C and/or 101D that may be located in one or more data centers on shore, aboard ships, and/or located in varying countries on different continents).

[0117] A processor may include, without limitation, a microprocessor, microcontroller, processor module or subsystem, programmable integrated circuit, programmable gate array, or another control or computing device.

[0118] The storage media 106 may be implemented as one or more computer-readable or machine-readable storage media. Note that while in the example embodiment of FIG. 5 the storage media 106 are shown as being disposed within the individual computer system 101A, in some embodiments, the storage media 106 may be distributed within and/or across multiple internal and/or external enclosures of the individual computing system 101A and/or additional computing systems, e.g., 101B, 101C, 101D. Storage media 106 may include, without limitation, one or more different forms of memory including semiconductor memory devices such as dynamic or static random access memories (DRAMs or SRAMs), erasable and programmable read-only memories (EPROMs), electrically erasable and programmable read-only memories (EEPROMs) and flash memories; magnetic disks such as fixed, floppy and removable disks; other magnetic media including tape; optical media such as compact disks (CDs) or digital video disks (DVDs); or other types of storage devices. Note that computer instructions to cause any individual computer system or a computing system to perform the tasks described above may be provided on one computer-readable or machine-readable storage medium, or may be provided on multiple computer-readable or machine-readable storage media distributed in a multiple component computing system having one or more nodes. Such computer-readable or machine-readable storage medium or media may be considered to be part of an article (or article of manufacture). An article or article of manufacture can refer to any manufactured single component or multiple components. The storage medium or media can be located either in the machine running the machine-readable instructions, or located at a remote site from which machine-readable instructions can be downloaded over a network for execution.

[0119] It should be appreciated that computing system 100 is only one example of a computing system, and that any other embodiment of a computing system may have more or fewer components than shown, may combine additional components not shown in the example embodiment of FIG. 5, and/or the computing system 100 may have a different configuration or arrangement of the components shown in FIG. 5. The various components shown in FIG. 5 may be implemented in hardware, software, or a combination of both hardware and software, including one or more signal processing and/or application specific integrated circuits.

[0120] Further, the acts of the processing methods described above may be implemented by running one or more functional modules in information processing apparatus such as general purpose processors or application specific chips, such as ASICs, FPGAs, PLDs, GPUs, coprocessors or other appropriate devices. These modules, combinations of these modules, and/or their combination with general hardware are all included within the scope of the present disclosure.

REFERENCES CITED IN THE DISCLOSURE

[0121] Aki, K., and Richards, P. G., 1980, Quantitative seismology theory and methods, Volume 1, W.H. Freeman and Company, New York. [0122] Plessix, R.-E., 2006, A review of the adjoint-state method for computing the gradient of a functional with geophysical application, Geophysical Journal International, 167, 2, 495-503 [0123] Shuey, R. T., 1985, A simplification of the Zoeppritz equations, Geophysics, 50(4), 609-614. [0124] Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation, Geophysics, 49, 8, 1259-1266 [0125] Verschuur, D. J., Staal, X. R. and A. J. Berkhout, 2016, Joint migration inversion: Simultaneous determination of velocity fields and depth images using all orders of scattering, The Leading Edge 35: 1037-1046. [0126] Zoeppritz, Karl, 1919, VIIb. Über Reflexion and Durchgang seismischer Wellen durch Unstetigkeitsflächen, [VIIb. On reflection and transmission of seismic waves by surfaces of discontinuity], Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-physikalische Klasse, 66-84.

[0127] In light of the principles and example embodiments described and illustrated herein, it will be recognized that the example embodiments can be modified in arrangement and detail without departing from such principles. The foregoing discussion has focused on specific embodiments, but other configurations are also contemplated. In particular, even though expressions such as in “an embodiment,” or the like are used herein, these phrases are meant to generally reference embodiment possibilities, and are not intended to limit the disclosure to particular embodiment configurations. As used herein, these terms may reference the same or different embodiments that are combinable into other embodiments. As a rule, any embodiment referenced herein is freely combinable with any one or more of the other embodiments referenced herein, and any number of features of different embodiments are combinable with one another, unless indicated otherwise. Although only a few examples have been described in detail above, those skilled in the art will readily appreciate that many modifications are possible within the scope of the described examples. Accordingly, all such modifications are intended to be included within the scope of this disclosure as defined in the following claims.