This invention generally relates to the field of seismic prospecting and reservoir delineation and, more particularly, to subsurface rock property prediction in a fractured media using geophysical seismic data. More specifically, the invention provides a system and method to remove the effects of fracture-induced anisotropy on the seismic amplitude of converted-wave data seismic data thereby allowing the data to be inverted for at least one subsurface rock property using traditional seismic inversion schemes.
It is well know that seismic reflection amplitudes vary with the incidence angle or the acquisition offset. This variation in the seismic amplitude is governed by subsurface rock properties, such as P-wave velocity or impedance, S-wave velocity or impedance, and density. It is a fairly common practice in the geophysics community to predict or invert for these subsurface properties by exploiting this property of the seismic reflection amplitudes. Such inversion methodologies are conventionally referred to as Amplitude-vs-Angle (AVA) or Amplitude-vs-Offset (AVO) inversion. The most common types of reflection modes used in the geophysics community are PP and PS modes. A reflected PP mode is excited when both the incident and the reflected waves are P-wave, whereas a PS mode is defined by an incident P-wave but a reflected S-wave.
Reflectivity at the interface between two elastic media can be calculated by Zoeppritz's equations. However, Zoeppritz's equations are complex and, in their original form, do not provide any insight to the physics of the wave propagation. Due to their nonlinear form, they are difficult to use in linear inverse problems and may cause nonuniqueness in the inversion solutions. Hence, a number of its linearized forms have been proposed in the last thirty years for both isotropic and anisotropic media. Aki and Richard [1] derived linearized forms of reflectivity between two isotropic media for all possible modes of wave propagation. Their equations have been further simplified by several authors and have served as a basis for most of the AVA-based inversion.
As known by those of ordinary skill in the art, if the medium is isotropic or has polar anisotropy, the seismic reflection amplitudes vary only with the incidence angle. However, if the medium is azimuthally anisotropic, the seismic amplitudes vary both with the incidence angle and the azimuth of the incidence plane. It is well known that at seismic wavelengths, fractured reservoirs exhibit azimuthal anisotropy. More specifically, one set of parallel vertical fractures in isotropic rocks causes HTI anisotropy. Ruger and Tsvankin [2] proposed a seismic inversion method for HTI media. Their method, however, is limited to PP data only. They derived new linearized equations for the reflection coefficients for HTI media and they showed that the reflectivity between two HTI media is governed by the additional anisotropic or Thomsen's parameter [3] along with P-wave velocity and density.
S-waves while travelling through fractured rocks split into the fast (S1) and slow (S2) S-waves. This phenomenon, along with the anisotropic nature of reflection coefficients, compounds the problem of seismic inversion of PS modes. Due to azimuthally varying elastic properties in HTI media, seismic reflectivity also varies with azimuth of wave propagation. The extra anisotropic constants in HTI media and complex wave propagation render the equations for reflection coefficients very complex and cause the AVA-based inverse problem to be intractable. Due to these problems, inversion for PS data in azimuthally anisotropic media is not very popular.
Jilek [4] proposed a method to invert for isotropic and anisotropic subsurface parameters in azimuthally anisotropic media. His method simultaneously inverts for both isotropic and anisotropic parameters, which render his approach unstable and unsuitable for application on field data which usually have low signal-to-noise ratio.
Thus, there is a need for improvement in this field.
The present invention provides a system and method for removal of fracture-induced anisotropy from converted wave seismic amplitudes.
In one embodiment, the invention is a method for inferring one or more physical property parameters of a subsurface media by inverting converted wave data acquired during a seismic survey, the method comprising: generating a composite seismic trace at a plurality of survey azimuths, said composite traces being composed such that their amplitudes are free of effects of subsurface anisotropy; and inverting at least one composite seismic trace by isotropic inversion to determine a property parameter of a subsurface media.
The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings.
The invention will be described in connection with example embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims. Persons skilled in the technical field will readily recognize that in practical applications of the present inventive method, it must be performed on a computer, typically a suitably programmed digital computer.
As noted above, Jilek [4] discloses a variant of the linearized form of reflectivity in anisotropic media for converted waves (PS mode). However, his equations are very complex and it is difficult for them to be used to set up a working inverse problem. The present disclosure provides an alternate form of the PS reflection coefficients between two HTI media and demonstrates how the reflection coefficients may be utilized.
The linearized form of reflection coefficients between two HTI media for PSV and PSH modes, which are the vertical and horizontal components of the reflected PS mode, can be written as follows:
R
PSV
HTI (θ, φ)=(DPSVISO,S+DPSVANISO,S)sin θ+(EPSVISO,S+EPSVANISO,S)sin3 θ (1)
and
R
PSH
HTI (θ, φ)=DPSHANISO sin θ+EPSHANISO sin3 θ, (2)
where θ is the incidence angle of the slowness vector. It is important to note that RPSVHTI is controlled by both isotropic and anisotropic terms, whereas RPSHHTI is controlled by anisotropic terms only. In deriving the above equations, it was assumed that the contrasts between the two surfaces are small. It was also assumed that the azimuth of the fracture normal is same in both media. Barring anisotropic terms, equation 1 is very similar to the PS reflection coefficient given by Khare et al. [5].
As will be demonstrated below, the anisotropic terms are dependent on anisotropic constants and the azimuth φ of the incidence plane with respect to the fracture normal. In the given form, it is assumed that the fracture normal is pointing in the x-direction. By using a perturbation approach, the following expressions can be derived:
In DPSVISO,S and EPSVISO,S, S signifies slow S-wave based attributes, βS, IS, and ρ are averages of the slow S-wave velocities, slow S-wave impedances and densities of the surfaces. ΔβS, ΔIsS and ρ are the differences of slow S-wave velocities, slow S-wave impedances and densities between the surfaces. g is the ratio of the average P- and S-wave velocities. The anisotropic components of RPSVHTI can be written as follows:
where δ(v) and ε(v) are Tsvankin's Thomsen-style anisotropic parameters for HTI media [6]. γ0 is the generic Thomsen's parameter [3]. Similarly, we can write the following expressions for PSH mode:
Unlike isotropic media, where an incident P-wave on a flat reflector can only generate a converted PSV mode, both PSV and PSH modes can be generated in HTI media. If the incidence medium is isotropic and the reflecting medium is HTI, no S-wave splitting occurs and both PSV and PSH modes are recorded concurrently at the receiver. If the incident medium is HTI, S-wave splitting occurs and the reflected S-wave travels as fast (PS1) and slow (PS2) modes and a time-lag is introduced between both modes.
Bansal et al. [7], which is incorporated by reference in its entirety, presented a method to generate PS1 and PS2 modes from radial and transverse PS data. This process can be referred to as 2-component rotation, or 2C rotation. Bansal et al. [7] also taught how to remove the time-lag in the PS2 mode in order to align it to the PS1 mode. The time-aligned PS2 mode may also be referred to as the time-compensated PS2 mode. In other embodiments, the PS1 mode may be shifted to align with the PS2 mode. In either case, the resulting PS1 and PS2 mode data has been time-adjusted.
It can be shown that PSV and PSH modes can be derived from PS1 and time-compensated PS2 data as follows:
where φ is the azimuth of the survey with respect to the fractures normal (symmetry axis). Equation 9 is strictly valid only for vertical propagation. However, in weakly anisotropic medium, where Vp/Vs is reasonably large, it remains valid for most practical purposes.
As shown in
But as shown in
but also by anisotropic parameters (Δδ(v), Δγ, Δε(v)). If the seismic amplitudes of PSV modes from orthogonal azimuths are added together, the following is obtained:
In DPSVISO,F and EPSVISO,F, F signifies fast S-wave based attributes (see equations 3 and 4). The fast S-wave based attributes were utilized in order to allow some of the anisotripic terms to cancel out. Numerical modeling was performed and it was found that the fast S-wave based attributes introduce smaller error into the approximation than slow S-wave based attributes.
In reviewing equation 10, it is first noted that the multiplier of sinθ is independent of azimuth angle φ. Moreover, in fractured rocks, δ(v) is either negative (in dry fractures) or close to zero (in fully-saturated fractures) and γ is positive. Hence in equation 10,
is very small for most practical purposes. Similarly, anisotropic terms in the multiplier of sin3 θ are also usually small. As a result, the sum of the reflection coefficients of orthogonal PSV modes can be expected to be independent of azimuth φ.
If the anisotropic coefficients are ignored, equation 10 takes the following form:
In equations 11, 12 and 13, βF, IsF and ρ are averages of the fast S-wave velocities, fast S-wave impedances and densities of the media. For example, βF=(β1F+β2F)/2 where β1F and β2F are the fast S-wave velocities of the incident and reflecting media, respectively. ΔβF , ΔIsF and Δρ are the differences of fast S-wave velocities, fast S-wave impedances and densities between the two media. It is important to note that the quantities in equations 11, 12 and 13 are free of anisotropic parameters.
When reviewing these figures, the difference between the two approximations initially appears to be large. However, after considering the error introduced during data acquisition, processing, and generating angleazimuth stacks, the difference in the AVA pattern is negligible. Moreover, the maximum difference between the two curves occurs when the fractures are fully-saturated with water, which is rarely the case. As a result, the difference between the anisotropic and isotropic approximations is very small and can be ignored for most practical purposes.
In the above discussion, the PSV and PSH modes were determined from PS1 and PS2 mode data. However, fracture parameters may also be estimated by a data rotation scheme first described by Alford [8] for vertical seismic profile (VSP) surveys and later modified by Gaiser [9, 10] for surface seismic survey. This rotations scheme is referred to herein as 4-component (4C) rotation.
The method of 4C rotation performs tensor rotation on the radial and transverse components from orthogonal azimuths. The following equation describes the tensor rotation process:
In the above equation, φ is the azimuth of the survey. PSRφ and PSTφ are radial and transverse components, respectively, of the PS mode at survey azimuth φ. PSRφ+90 and PSTφ+90 are radial and transverse components from the orthogonal azimuth φ+90. This tensor rotation moves the seismic energy from the off-diagonal components (PSTφ and PSTφ+90) to diagonal components. The new diagonal components, PS11 and PS22, carry most of the seismic energy and off-diagonal components, PS12 and PS21, modes have typically very little energy left. PS11 is represents the pseudo-fast mode, while PS22 represents the pseudo-slow mode.
After performing the 4C rotation, PS22 is aligned with PS11 mode by applying a positive time-shift on PS22 mode. The new time-aligned PS22 mode may be referred to as the time-compensated PS22 mode. Alternatively, PS11 is aligned with PS22 mode by applying a negative time-shift on PS11 mode. In either case, the resulting PS11 and PS22 mode data has been time-adjusted.
After matrix multiplication of equation 14, it can be shown that PS11 and time-compensated PS22 modes and PSV are related as follow:
PS
11
φ
+PS
22
φ
=PSV
φ
+PSV
φ+90 (15)
Similar to the PS1 and PS2 modes, the PSV and PSH modes can be thought of as re-computed radial and transverse components if no time-lag existed between PS11 and PS22 modes. Since the time-lag between PS11 and PS22 has been removed, PSV and PSH are same as PSR and PST, respectively. As a result, equation 15 can be derived from equation 14.
The flowchart of
In this type of acquisition, any type of P-wave source is used, which may be, but is not limited to, dynamite, a vertical vibrator, or air gun. The vertical component of the data mostly contains P-wave energy and the two horizontal components carry converted-wave PS energy. The PS energy is defined as the P-wave energy reflected back from a reflector as S-wave energy; i.e., P-wave goes down and some of that energy is reflect back up in an S-wave mode.
The process continues as the acquired PS energy is rotated or resolved into radial (direction along the line connecting source to receiver) and transverse (direction perpendicular to the line connecting source to receiver) components (603). Free-surface related seismic noise, such as surface-waves and free-surface multiples, may then be removed from the radial and transverse components (605). After noise correction, normal moveout (NMO) correction is applied (607) on the data to flatten the reflections. In other embodiments, the reflections may be flattened by pre-stack time migration. As appreciated by those skilled in the relevant art, steps 601-607 are standard processing steps and are routinely applied in seismic data processing.
The process continues by generating composite seismic traces at all azimuths whose amplitudes are free of the effects of subsurface anisotropy (step 609). These composite traces may then be further used to estimate subsurface rock properties. The next steps of the illustrated embodiment depend on how the PS energy is processed. In one instance, PS11 and time-compensated PS22 gathers are generated (step 611). In one embodiment, the PS11 and PS22 modes are generated by the 4C rotation process described herein. The PS11 and time-compensated PS22 gathers at common survey aziumuths are then summed to create composite traces (step 613).
In another instance, PS1 and time-compensated PS2 gathers are generated (step 615). In one embodiment, the PS1 and PS2 are generated by the 2C rotation process described herein. Next, PS1 and time-compensated PS2 gathers are used to generate PSV traces at all survey azimuths using equation 9 (step 617). PSV traces from orthogonal azimuths are then summed to create composite traces (step 619).
Regardless of how the composite traces are generated, process 600 continues by stacking the composite traces from all azimuths to generate a full-azimuth composite trace (step 621). In other embodiments, the composite traces from some, but not all, azimuths are stacked to generate a multi-azimuth composite trace. Due to stacking of the traces from the different azimuths, the signal-to-noise ratio of the full-azimuth composite trace improves. Because the seismic amplitudes of the composite traces are substantially free of subsurface anisotropy, isotropic inversion is performed on the full-azimuth composite trace (step 623). In one embodiment, subsurface rock properties are inverted for by using the equation 11.
Another embodiment of the present disclosure includes only steps 609 and 623. In additional embodiments, other steps may also be included to meet system and design objectives.
Also depicted are the seismic amplitudes, in polar plots, of the PS11 and PS22 traces (
However, as seen in
Equation 16 allows subsurface rock properties to be calculated based upon known isotropic inversion techniques.
Also depicted are the seismic amplitudes, in polar plots, of the PS1 and PS2 traces (
Equation 17 allows subsurface rock properties to be calculated based upon isotropic inversion techniques.
Subsurface rock properties can be solved in a variety of ways. In one embodiment, forwarding modeling is performed using equation 16 or 17 and using a set of rock properties
The modeled synthetic traces are compared with the field data, i.e., the generated composite traces. If there is a good match between the two, it is assumed that the correct rock properties were used to do the modeling. If there is not a good match, then modeling is performed again using a new set of rock properties. This process is iterated until there is a good match between the synthetics and the composite traces. The final rock properties used to do the modeling is assumed to be the correct subsurface properties.
Another way to estimate the rock properties are by directly inverting for them from the seismic data. In such approach, a least-squares inverse problem is set up using equation 16 or 17. The seismic data (or, in the case of the present invention, composite traces) at multiple incidence angles is used as the input to the inverse problem. The least-squares inverse problem is then solved to estimate at least one rock property at each time sample.
While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiments have been shown and described and that all changes and modifications that come within the wording of the claims are desired to be protected. It is also contemplated that structures and features embodied in the present examples can be altered, rearranged, substituted, deleted, duplicated, combined, or added to each other. The articles “the”, “a” and “an” are not necessarily limited to mean only one, but rather are inclusive and open ended so as to include, optionally, multiple such elements.
This application claims the benefit of U.S. Provisional Patent Application 61/567,516, filed Dec. 6, 2011, entitled REMOVAL OF FRACTURE-INDUCED ANISOTROPY FROM CONVERTED-WAVE SEISMIC AMPLITUDES, the entirety of which is incorporated by reference herein.
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/US2012/059376 | 10/9/2012 | WO | 00 | 4/23/2014 |
Number | Date | Country | |
---|---|---|---|
61567516 | Dec 2011 | US |