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.

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; and 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.

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; and 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.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) 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.

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

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

(4) FIG. 3A illustrates a single shot gather.

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

(6) FIG. 4A illustrates angle stacks with acoustic FWI.

(7) FIG. 4B illustrates angle stacks with convolutions.

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

DETAILED DESCRIPTION

(9) 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.

(10) 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.

(11) 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.

(12) 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.

(13) 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.

(14) 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.

(15) 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.

(16) 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.

(17) 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:

(18) $\begin{array}{cc}R\left(\right)=\frac{1}{2}\frac{}{}+\frac{1}{2{\cos }^2}\frac{}{}& \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

(19) $\frac{1}{2{\cos }^2}.$
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.

(20) 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.

(21) 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 (.sub.0.sub.f), where .sub.0 and .sub.f are the minimum and maximum frequencies present in the data, at each angle the reflectivity spectrum is only sampled from

(22) $\frac{{}_0}{}\cos \mathrm{to}\frac{{}_f}{}\cos$
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(.sub.1) and R(.sub.2) are the same, the relation between the inverted reflectivities would be

(23) $\frac{R\left({}_1\right)}{R\left({}_2\right)}=\frac{\cos {}_1}{\cos {}_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.

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

(25) 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.

(26) 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.

(27) 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

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

(29) [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;

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

(31) [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;

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

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

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