FWI model domain angle stacks with amplitude preservation

Abstract

A method, including: obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges; performing, with a computer, an acoustic full wavefield inversion process on each of the subsets, respectively, to invert for density and generate respective density models; generating acoustic impedances for each of the subsets, as a function of reflection angle, using the respective density models; and transforming, using a computer, the acoustic impedances for each of the subsets into reflectivity sections, wherein the transforming includes normalizing the reflectivity sections by their respective bandwidth.

Description

FIELD OF THE INVENTION

Exemplary embodiments described herein pertain to the field of geophysical prospecting, and more particularly to geophysical data processing. Specifically, embodiments described herein relate to a method for more efficiently generating FWI model domain angle stacks.

BACKGROUND

This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present invention. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present invention. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.

An important goal of seismic prospecting is to accurately image subsurface structures commonly referred to as reflectors. Seismic prospecting is facilitated by obtaining raw seismic data during performance of a seismic survey. During a seismic survey, seismic energy is generated at ground level by, for example, a controlled explosion, and delivered to the earth. Seismic waves are reflected from underground structures and are received by a number of sensors referred to as geophones. The seismic data received by the geophones is processed in an effort to create an accurate mapping of the underground environment. The processed data is then examined with a goal of identifying geological formations that may contain hydrocarbons.

Full Wavefield Inversion (FWI) is a geophysical method which is used to estimate subsurface properties (such as velocity or density). It is known to be advanced for the higher resolution and more accurate physics compared to conventional methods. The fundamental components of an FWI algorithm can be described as follows: using a starting subsurface physical properties model, synthetic seismic data are generated by solving a wave equation using a numerical scheme (e.g., finite-difference, finite-element etc.). The synthetic seismic data are compared with the field seismic data and using the difference between the two, the value of an objective function is calculated. To minimize the objective function, a modified subsurface model is generated which is used to simulate a new set of synthetic seismic data. This new set of synthetic seismic data is compared with the field data to recalculate the value of the objective function. The objective function optimization procedure is iterated by using the new updated model as the starting model for finding another search direction, which will then be used to perturb the model in order to better explain the observed data. The process continues until an updated model is found that satisfactorily explains the observed data. A global or local optimization method can be used to minimize the objective function and to update the subsurface model. Commonly used local objective function optimization methods include, but are not limited to, gradient search, conjugate gradients, quasi-Newton, Gauss-Newton and Newton's method. Commonly used global methods included, but are not limited to, Monte Carlo or grid search.

Although FWI is expected to provide the subsurface properties, it is difficult to extract the correct viscoelastic properties from the seismic data directly with FWI. As FWI estimates the properties by fitting the data with synthetic waveforms, it relies on how accurate the wave equation can explain the actual physics, and how well the optimization method can separate the effects from different properties. When an acoustic wave equation is used, FWI can generate P-wave velocity models based on the travel time information in the datasets. However, the amplitude information is not fully utilized because the real earth is visco-elastic, and an acoustic model cannot explain all the amplitudes in the acquired data. If FWI is expected to provide interpretable products like elastic impedances, elastic simulation is often needed but very expensive; in general it is 6 to 10 times the computation of acoustic FWI. In addition, the initial model for shear wave velocity is difficult to obtain due to the limited shear wave kinematic information and often poor signal to noise ratio in the acquisitions.

An alternative way of using the elastic amplitude information is to form angle stacks. Amplitude versus angle (AVA) analysis [5] can be performed on the angle stacks to extract the elastic properties. Traditional AVA stacks generated with Kirchhoff migration need geometric spreading corrections to account for the amplitude loss during propagation. However, it is not guaranteed that the amplitude after correction would reflect the true amplitude of the data. In addition, Kirchhoff migration is based on ray-tracing which favors smooth velocity models and would likely fail in high contrast medium. Angle calculations are under a 1-D assumption that is not accurate enough when subsurface structures are complex. Reverse time migration (RTM) based angle stacks [1, 2] are more advanced for making use of the high-resolution velocity models. Nonetheless, amplitude preservation is still difficult. Yu Zhang et al (2014) [3] reported a least-squares RTM to balance the image amplitudes; however, it has not been proved to be able to generate angle stacks.

SUMMARY

A method, including: obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges; performing, with a computer, an acoustic full wavefield inversion process on each of the subsets, respectively, to invert for density and generate respective density models; generating acoustic impedances for each of the subsets, as a function of reflection angle, using the respective density models; and transforming, using a computer, the acoustic impedances for each of the subsets into reflectivity sections, wherein the transforming includes normalizing the reflectivity sections by their respective bandwidth.

In the method, each of the full wavefield inversion processes start from a same velocity model.

In the method, each of full wavefield inversion processes are independently applied to the subsets.

In the method, the obtaining includes dividing a shot gather into the subsets by using a data mask that includes information of reflector dipping angles and P-wave velocity.

The method can further include using, for each of the reflectivity sections, a Fourier transform, discrete Fourier transform, or a fast Fourier transform to calculate an average spectrum within at least one local window that is applied at a same location to all of the reflectivity sections, and determining a bandwidth for each average spectrum.

In the method, the determining the bandwidth is based on a distance between 10-dB points.

In the method, the determining the bandwidth is based on a distance between points with steepest slope.

In the method, the average spectrum is calculated within a plurality of local windows, and is averaged.

The method can further include determining reflectivity values at a plurality of angles and constructing an angle-vs-amplitude curve by interpolation.

The method can further include managing hydrocarbon production using the reflectivity sections.

In the method, the managing hydrocarbon production includes drilling a well at a location determined at least in part by the reflectivity sections.

A non-transitory computer readable storage medium encoded with instructions, which when executed by a computer cause the computer to implement a method including: obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges; performing, with a computer, an acoustic full wavefield inversion process on each of the subsets, respectively, to invert for density and generate respective density models; generating acoustic impedances for each of the subsets, as a function of reflection angle, using the respective density models; and transforming, using a computer, the acoustic impedances for each of the subsets into reflectivity sections, wherein the transforming includes normalizing the reflectivity sections by their respective bandwidth.

BRIEF DESCRIPTION OF THE DRAWINGS

While the present disclosure is susceptible to various modifications and alternative forms, specific example embodiments thereof have been shown in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific example embodiments is not intended to limit the disclosure to the particular forms disclosed herein, but on the contrary, this disclosure is to cover all modifications and equivalents as defined by the appended claims. It should also be understood that the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating principles of exemplary embodiments of the present invention. Moreover, certain dimensions may be exaggerated to help visually convey such principles.

FIG. 1 illustrates an exemplary method for generating FWI AVA stacks.

FIG. 2 illustrates the spectrum spreading effect due to scattering angles.

FIG. 3A illustrates a single shot gather.

FIG. 3B illustrates the single shot gather that is muted into five different angle ranges.

FIG. 4A illustrates angle stacks with acoustic FWI.

FIG. 4B illustrates angle stacks with convolutions.

FIG. 5 illustrates the vertical lines from each angle stack in FIGS. 4A and 4B, overlaid to show consistency.

DETAILED DESCRIPTION

Exemplary embodiments are described herein. However, to the extent that the following description is specific to a particular embodiment, this is intended to be for exemplary purposes only and simply provides a description of the exemplary embodiments. Accordingly, the invention is not limited to the specific embodiments described below, but rather, it includes all alternatives, modifications, and equivalents falling within the true spirit and scope of the appended claims.

Exemplary embodiments described herein provide a method than can: 1) be robust with complex geology; 2) preserve the amplitude versus angle information; and 3) be less expensive than elastic FWI. The proposed FWI model domain angle stacks can be generated by inverting the datasets of different angle ranges for different acoustic models. Amplitude preservation can be achieved through the data fitting process, and the angle calculation is more accurate using Poynting vectors. The Poynting vector describes energy flow for body waves, interface waves, guided waves and inhomogeneous waves in isotropic and anisotropic media. Poynting vectors naturally take the advantage of the high-resolution FWI velocity models. When implemented as an integrated part of a FWI workflow, the present technological advantage can utilize the FWI products angle stacks without changing platforms. More importantly, exemplary method are not limited to the incomplete physics in the modeling engine.

Exemplary embodiments of the present technological advancement generates model domain amplitude preserved angle stacks using FWI. Advantageously, the present technological advancement can use only acoustic simulations, but can be applied to the full offsets of the acquired seismic data. It is impossible to use one acoustic model to fit all the data that contains all kinds of physics. However, if the datasets are separated by the reflection angles, for each angle, there is an acoustic model that can explain the data. With all the models combined, an impedance model is formed as a function of reflection angle: I(θ). Acoustic impedance is a measure of the ease with which seismic energy travels through a particular portion of the subsurface environment. Those of ordinary skill in the art will appreciate that acoustic impedance may be defined as a product of density and seismic velocity. From the impedance, the reflectivity can be derived as a function of angle: R(θ), which is exactly the definition of AVA.

In practice, it is not necessary to find a continuous form of R(θ). Instead, R(θ) can be determined at several angles, and the AVA curve can be reconstructed by interpolation.

FIG. 1 illustrates an exemplary method for generating FWI AVA stacks. In step 101, processed data is obtained, which can be a single shot gather generated from the collected seismic data. Such processed data is seismic data conditioned according to conventional techniques known to those of ordinary skill in the art. In step 102, the shot gather is divided into several subsets. Each subset is within a relatively small range of reflection angles. This can be achieved by using carefully designed data masks which include the information of reflector dipping angles and P-wave velocity. With the dipping angles and P-wave velocity, a ray-tracing method can be used to find the time and offset of the reflected wave in the data from each of the subsurface points at each reflection angle. Data masking is a process of hiding a portion of the original data so that only the desired portion (i.e., desired angle range) is available for subsequent processing. The angle ranges shown in the accompanying figures are only exemplary, and other angle ranges can be used. Moreover, the number of angle ranges is also exemplary as the shot gather could be divided in to more or less sections.

In step 103, on each of the data subsets generated in step 102, acoustic FWI is applied to obtain acoustic impedances independently. Those of ordinary skill in the art are familiar with acoustic FWI and further details of this process are omitted. All inversions can start from the same velocity model, and the kinematics are not updated in this process assuming that velocity model building is already finished and accurate enough. The model updates are meant to explain the data amplitude only. After the inversion, the synthetic waveforms simulated with these impedances fit the real data well so that the amplitude information is preserved in the model domain. Each impedance model can only explain the data of a certain range of reflection angles constrained by the data masks. Within one angle range, the mid angle can be chosen to be the nominal angle of the reflectivity. This is additionally guaranteed by using Poynting vectors [2] when forming the gradient in the inversion, and so the gradient is most sensitive to the reflections at the nominal angle. Poynting vectors are used to separate the wave propagation directions during the finite difference simulation and gradient calculation. The data separation can be conducted based on ray-theory. Since the Poynting vector is based on wave-theory, it can be a helpful check on the accuracy of the data separation.

The term velocity model, density model, or physical property model as used herein refers to an array of numbers, typically a 3-D array, where each number, which may be called a model parameter, is a value of velocity, density or another physical property in a cell, where a subsurface region has been conceptually divided into discrete cells for computational purposes.

After the inversion, both density and velocity are known and step 104 can determine impedance from the results of the acoustic FWI as impedance is a function of density and velocity.

In acoustic FWI, two parameters can be inverted for to fit the data amplitudes: P-wave velocity and density. P-wave velocity is often chosen to fit the amplitude and travel time at the same time when an L-2 norm type of objective function is used. However, density may be a better parameter for reflectivity inversion. Density has a much simpler AVA response than P-wave velocity. As described in the Aki-Richards equation:

$\begin{array}{cc}R\left(\theta \right)=\frac{1}{2}\frac{\Delta \rho }{\rho }+\frac{1}{2{\cos }^2\theta }\frac{\Delta \alpha }{\alpha }& \left[1\right]\end{array}$ where Δp is density perturbation and Act is the P-wave velocity perturbation, density has a constant AVA response which does not vary with angle [6]. This indicates that when a density perturbation is inverted for to fit the data amplitude at a certain angle θ, the value of the perturbation directly represents the reflectivity at that angle θ regardless of the actual value of θ. On the contrary, if a P-wave velocity perturbation is used, in order to obtain R(θ), we need to apply a correction of

$\frac{1}{2{\cos }^2\theta }.$ Moreover, Equation [1] is only valid when the perturbation is weak, and when the perturbation is strong, the correction does not have an explicit form. For density, the constant AVA response is valid for all cases.

After the acoustic impedances are obtained, they can be shaped or converted into reflectivity sections (P-P reflectivity) in step 105. The reflectivity sections can be approximately determined from the derivative of the acoustic impedance with respect to space (or more generally the vertical derivative, which in some cases could be time). However, there is one more step to balance the reflectivity spectrum across different angles, i.e., “stretch”. While FIG. 1 shows “shape” and “stretch” in the same step, these are not necessarily performed simultaneously.

Similar to the wavelet stretching effect in migration, the reflectivity's obtained from data of different reflection angles are of different resolutions. Because there is only bandlimited data (ω₀˜ω_f), where ω₀ and ω_f are the minimum and maximum frequencies present in the data, at each angle θ the reflectivity spectrum is only sampled from

$\frac{{\omega }_0}{\alpha }\cos \theta \mathrm{to}\frac{{\omega }_f}{\alpha }\cos \theta$ in the wavenumber domain as shown in FIG. 2. Different bandwidth in the wavenumber domain leads to different amplitude in the space domain. Assuming the true reflectivity values R(θ₁) and R(θ₂) are the same, the relation between the inverted reflectivities would be

$\frac{R\left({\theta }_1\right)}{R\left({\theta }_2\right)}=\frac{\cos {\theta }_1}{\cos {\theta }_2}.$ Therefore, to preserve the data AVO in the model domain, there is a need to compensate the spectrum stretch by dividing by a factor of cos θ. In practice, it is difficult to use data of a single reflection angle. Therefore, an alternative way of obtaining the compensation factor is to measure the bandwidth of the reflectivity. For all the reflectivity sections, a Fourier transform can be used to calculate the averaged spectrum within a local window that is applied at the same location to all sections. This can be performed at multiple locations and averaged, as long as the locations are the same in all sections. The Fourier transform may be preferred, but other transforms could be used, such as FFT or DFT. The bandwidth can be defined, for example, as the distance between the 10-dB points. However, other measures of bandwidth can be used (i.e., 3-dB points, full-width half maximum, points of steepest slope, etc.), but the distance between the 10-dB points may be preferred. Then, to complete step 105 and compensate for the “stretch,” each reflectivity section is normalized by its own bandwidth so that the spectrum stretching effect on the reflectivity amplitude is corrected for.

The final output of the method in FIG. 1 (step 106) are the reflectivity stacks for different angles.

The present technological advancement was applied on a synthetic dataset generated with a 2-D slice extracted from the SEG SEAM Phase I model. Pressure data are simulated with streamers of 4 km maximum offset. An absorbing boundary condition is used on the water surface. Therefore, no free surface related multiples are present in the data. As shown in FIG. 3A, a single shot gather is divided into five sections by the dashed lines 301, 302, 303, 304, and 305, which corresponds to the following angle ranges: 0-10; 10-20; 20-30; 30-40; and 40-50 degrees. The angle ranges do not necessarily need to overlap. Five data masks are designed based on these lines 301-305, each covering 10 degrees of reflection angle. Therefore, five data subsets are generated (see step 102) as shown in FIG. 3B. Starting from the same velocity model, acoustic FWI was performed (see step 103) using each of the subsets respectively. The resulting acoustic impedances (see step 104) after the inversion are “shaped” and “stretched” (see step 105) into reflectivity sections (see step 106), and shown in the panels in FIG. 4A. The gap area 401 in the last panel is because the large angle reflection from the deep part is outside of the acquisition offset. To verify the quality of the results, true reflectivity sections are generated by the convolution between the seismic wavelet and the reflectivity sequences, and are shown in FIG. 4B. The reflectivity sequences are calculated using the Zoeppritz equation [4] and the true synthetic model. The AVA in FIGS. 4A and 4B appears to be consistent. For a closer scrutiny, five vertical lines 402a, 403a, 404a, 405a, and 406a from the model domain stacks, and five from the convolution sections 402b, 403b, 404b, 405b, and 406b are overlaid in FIG. 5 to show the consistency. It is clear that at all angles and all depths the FWI model domain stacks have very similar amplitudes compared to the convolution stacks. It demonstrates that the workflow is reliable and accurate with a complicated geology.

The final reflectivity's are an example of a subsurface image that can be used for interpretation of the subsurface and/or management of hydrocarbon exploration. As used herein, hydrocarbon management includes hydrocarbon extraction, hydrocarbon production, hydrocarbon exploration, identifying potential hydrocarbon resources, identifying well locations, determining well injection and/or extraction rates, identifying reservoir connectivity, acquiring, disposing of and/or abandoning hydrocarbon resources, reviewing prior hydrocarbon management decisions, and any other hydrocarbon-related acts or activities.

In all practical applications, the present technological advancement must be used in conjunction with a computer, programmed in accordance with the disclosures herein. Preferably, in order to efficiently perform FWI, the computer is a high performance computer (HPC), known as to those skilled in the art, Such high performance computers typically involve clusters of nodes, each node having multiple CPU's and computer memory that allow parallel computation. The models may be visualized and edited using any interactive visualization programs and associated hardware, such as monitors and projectors. The architecture of system may vary and may be composed of any number of suitable hardware structures capable of executing logical operations and displaying the output according to the present technological advancement. Those of ordinary skill in the art are aware of suitable supercomputers available from Cray or IBM.

REFERENCES

The following references are hereby incorporated by reference in their entirety:

[1] Xu, S., Y. Zhang and B. Tang, 2011, 3D angle gathers from reverse time migration: Geophysics, 76:2, S77-S92. doi:10.1190/1.3536527;

[2] Thomas A. Dickens and Graham A. Winbow (2011) RTM angle gathers using Poynting vectors. SEG Technical Program Expanded Abstracts 2011: pp. 3109-3113;

[3] Yu Zhang, Lian Duan, and Yi Xie (2013) A stable and practical implementation of least-squares reverse time migration. SEG Technical Program Expanded Abstracts 2013: pp. 3716-3720;

[4] Encyclopedic Dictionary of Applied Geophysics, R. E. Sheriff, 4^th edition., SEG, 2002, p. 400;

[5] Encyclopedic Dictionary of Applied Geophysics, R. E. Sheriff, 4th edition., SEG, 2002, p. 12; and

[6] Aki, K and Richards, P (1980) Quantitative seismology, 2nd edition, University Science Books, p. 133-155.

Claims

1. A method, comprising: obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges;performing, with a computer, an acoustic full wavefield inversion process on each of the subsets, respectively, to invert for density and generate respective density models;generating acoustic impedances for each of the subsets, as a function of reflection angle, using the respective density models;transforming, using a computer, the acoustic impedances for each of the subsets into reflectivity sections, wherein the transforming includes normalizing the reflectivity sections by their respective bandwidth; andusing, for each of the reflectivity sections, a Fourier transform, discrete Fourier transform, or a fast Fourier transform to calculate an average spectrum within at least one local window that is applied at a same location to all of the reflectivity sections, and determining a bandwidth for each average spectrum.

2. The method of claim 1, wherein each of the full wavefield inversion processes start from a same velocity model.

3. The method of claim 1, wherein each of full wavefield inversion processes are independently applied to the subsets.

4. The method of claim 1, wherein the obtaining includes dividing a shot gather into the subsets by using a data mask that includes information of reflector dipping angles and P-wave velocity.

5. The method of claim 1, wherein the determining the bandwidth is based on a distance between 10-dB points.

6. The method of claim 1, wherein the determining the bandwidth is based on a distance between points with steepest slope.

7. The method of claim 1, wherein the average spectrum is calculated within a plurality of local windows, and is averaged.

8. The method of claim 1, further comprising determining reflectivity values at a plurality of angles and constructing an angle-vs-amplitude curve by interpolation.

9. The method of claim 1, further comprising managing hydrocarbon production using the reflectivity sections.

10. The method of claim 1, wherein the managing hydrocarbon production includes drilling a well at a location determined at least in part by the reflectivity sections.

11. A non-transitory computer readable storage medium encoded with instructions, which when executed by a computer cause the computer to implement a method comprising: obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges;performing, with a computer, an acoustic full wavefield inversion process on each of the subsets, respectively, to invert for density and generate respective density models;generating acoustic impedances for each of the subsets, as a function of reflection angle, using the respective density models;transforming, using a computer, the acoustic impedances for each of the subsets into reflectivity sections, wherein the transforming includes normalizing the reflectivity sections by their respective bandwidth; andusing, for each of the reflectivity sections, a Fourier transform, discrete Fourier transform, or a fast Fourier transform to calculate an average spectrum within at least one local window that is applied at a same location to all of the reflectivity sections, and determining a bandwidth for each average spectrum.