Multi-parameter inversion through offset dependent elastic FWI

Abstract

Method for multi-parameter inversion using elastic inversion. This method decomposes data into offset/angle groups and performs inversion on them in sequential order. This method can significantly speed up convergence of the iterative inversion process, and is therefore most advantageous when used for full waveform inversion (FWI). The present inventive approach draws upon relationships between reflection energy and reflection angle, or equivalently, offset dependence in elastic FWI. The invention uses recognition that the amplitudes of small angle (near offset) reflections are largely determined by acoustic impedance alone (1), independent for the most part of Vp/Vs. Large angle (middle and far offset) reflections are affected by Ip, Vp/Vs (2) and other earth parameters such as density (3) and anisotropy. Therefore, the present inventive method decomposes data into angle or offset groups in performing multi-parameter FWI to reduce crosstalk between the different model parameters being determined in the inversion.

Claims

1. A computer-implemented method for full wavefield inversion of seismic data to infer subsurface physical property parameters including P-wave (pressure wave) velocity, S-wave (shear wave) velocity, and density, comprising: extracting only PP mode (P-wave down/P-wave up) from the seismic data, and inverting, with a full wavefield inversion algorithm, the PP mode data sequentially in two or more different offset ranges, each offset range full wavefield inversion determining at least one physical property parameter, wherein in a second and subsequent full wavefield inversions, parameters determined in a previous inversion are held fixed, and further wherein all full wavefield inversions are performed using a computer; and generating, with a computer, a subsurface physical property model that includes the P-wave velocity, S-wave velocity, and the density from the full wavefield inversion algorithm, which transforms the seismic data into the subsurface physical property model, wherein the subsurface physical property model is a quantitative rock-property description of a hydrocarbon reservoir, and using the subsurface physical property model for geophysical prospecting.

2. The method of claim 1, wherein a near offset range is sequentially first to be inverted, and said first full wavefield inversion infers P-wave acoustic impedance Ip, using a computer programmed with an acoustic full wavefield inversion algorithm.

3. The method of claim 2, wherein a mid-offset range is sequentially second to be inverted, and said second full wavefield inversion infers S-wave acoustic impedance I.sub.S, or P-wave velocity V.sub.P divided by S-wave velocity V.sub.S, with I.sub.P fixed at its value from the first full wavefield inversion, said second full wavefield inversion using an elastic inversion algorithm.

4. The method of claim 3, wherein a far-offset range is sequentially third to be inverted and said third full wavefield inversion infers density or V.sub.P, using an elastic full wavefield inversion algorithm, with I.sub.P fixed at its value from the first full wavefield inversion and V.sub.P/V.sub.S fixed at a value determined from the second full wavefield inversion.

5. The method of claim 4, wherein V.sub.P and Vs are computed from I.sub.P and I.sub.S using definition of acoustic impedance, and using density as inferred in the third full wavefield inversion.

6. The method of claim 4, wherein V.sub.P is inferred in the third full wavefield inversion, and density is computed from the relationship I.sub.P=V.sub.P and I.sub.P is as determined in the first full wavefield inversion.

7. The method of claim 4, wherein one or both of the relationships I.sub.P=V.sub.P and I.sub.P=V.sub.s are used in performing the method.

8. The method of claim 4, further comprising repeating the sequential full wavefield inversions at least one time to update the inferred physical property parameters.

9. A computer-implemented method for inversion of seismic data to infer at least P-wave (pressure wave) velocity, S-wave (shear wave) velocity, and density, comprising: (a) taking only PP-mode (P-wave down/P-wave up) data from the seismic data, and dividing the seismic data into a near-offset range, a mid-offset range, and a far offset range, which ranges may or may not overlap; (b) inverting the near offset range for P-wave acoustic impedance IP, using a computer programmed with an acoustic full wavefield inversion algorithm; (c) inverting the mid-offset range for S-wave acoustic impedance Is, or for P-wave velocity V.sub.P divided by S-wave velocity V.sub.S, with I.sub.P fixed at its value from (b), using an elastic full wavefield inversion algorithm; (d) inverting the far-offset range for density, using an elastic full wavefield inversion algorithm, with I.sub.P fixed at its value from (b) and V.sub.P/V.sub.S fixed at a value determined from the value of I.sub.S from (c); (e) computing V.sub.P and V.sub.S from I.sub.P and I.sub.S using definition of acoustic impedance and density as determined in (d); and (f) generating, with a computer, a subsurface physical property model that includes the P-wave velocity, S-wave velocity, and the density, wherein the steps (b), (c), and (d) transform the seismic data into the subsurface physical property model, wherein the subsurface physical property model is a quantitative rock-property description of a hydrocarbon reservoir, and using the subsurface physical property model for geophysical prospecting.

10. The method of claim 1, wherein at least some of the two or more different offset ranges overlap.

11. The method of claim 1, wherein at least some of the two or more different offset ranges do not overlap.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The advantages of the present invention are better understood by referring to the following detailed description and the attached drawings, in which:

(2) FIG. 1 is a flowchart showing basic steps in one embodiment of the seismic processing method of the present invention;

(3) FIG. 2 shows the true Vp, Vs and density profiles used to generate a synthetic gather, and one of the shot gathers;

(4) FIG. 3 shows inversion of I.sub.P using near offset data and data misfit, compared with true I.sub.P and synthetic data;

(5) FIG. 4 shows I.sub.p alone without knowledge of V.sub.p/V.sub.s is not able to explain middle offset data;

(6) FIG. 5 shows inversion of V.sub.P/V.sub.S with I.sub.p fixed from FIG. 2 explains seismic data up to the middle offsets; and

(7) FIG. 6 shows the results of inversion of density from far offset data, with I.sub.P and V.sub.P/V.sub.S fixed from FIG. 2 and FIG. 4.

(8) Many of the drawings are color originals converted to gray scale because of patent law restrictions on the use of color.

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

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

(10) In the elastic FWI method presented by (SSB for short) Sears, Singh and Barton (2008), a three-stage workflow was proposed to estimate Vp and Vs from P-wave and S-wave seismic data: stage one, inversion for short and intermediate scale Vp using normal-incidence and wide-angle P-wave data; stage two, inversion for intermediate Vs using wide-angle P-wave data; and stage three, inversion for short-scale Vs using PS-wave data. Short and intermediate scale are terms used in the SSB paper. General speaking, short-scale refers to spatial scales that can be inferred directly from high frequency reflection energy in seismic data, and large-scale refers to spatial scales whose reflected frequencies are below typical seismic sources (e.g., 4-6 Hz in marine acquisition). Therefore, the large-scale is typically inferred from migration velocity analysis. The gap between large-scale and short-scale is usually called intermediate-scale.

(11) While the SSB method may at first appear similar to the 3-step inventive method that is disclosed herein, there are important features that distinguish them. First, the SSB method uses different wave modes through the 3 stages. The present inventive method uses the same wave mode (PP-wave) but different reflection angle/offset through the 3 steps. It is well known that PP-wave data represent most of the recorded energy in a typical seismic survey, and therefore most of the value in marine streamer acquisition. Second, the SSB method does not separate normal-incidence and wide-angle P-wave data in stage 1 and uses them simultaneously. The present inventive method uses only small angle reflection data in step 1, which is the critical step of speeding up convergence.

(12) A synthetic example is used to demonstrate that this method is very robust and effective in retrieving Ip and Vp/Vs. The total number of iterations needed to get Ip and Vp/Vs is 10. Retrieving density information in step 3 (see the FIG. 1 flow chart) may require an additional 10-15 iterations in the synthetic example. Tests on field data show that accurate and robust estimate of I.sub.p and V.sub.p/V.sub.s can be obtained within 10 iterations as well. However, in the field data case, the reliability of the density inversion is strongly subject to the accuracy of the velocity model, including anisotropy, and data quality at far-offsets.

(13) The synthetic example follows the embodiment of the present inventive method illustrated in the flow chart of FIG. 1. Synthetic (computer simulated) data are used in this test example to demonstrate the invention. The data set is generated by isotropic elastic finite difference modeling on a layered (1D) earth model shown in FIG. 2, where V.sub.P, V.sub.S and density are plotted vs. depth in the subsurface. The units for velocity and density are m/s and kg/m.sup.3. A common-shot gather of the synthetic measured data is also shown at 8 in FIG. 2. Time in seconds is plotted on the vertical axis, and offset in meters is plotted on the horizontal axis. The maximum depth of the earth model is 2.3 km and the maximum offset available is 5 km. Due to patent law restrictions on the use of color, the depicted shot gather 8 is a gray scale conversion of a colored data display, where color is used to represent the magnitude of seismic amplitudes. The same is true of the comparisons of simulated to measured data, and the misfits, shown in FIGS. 3-6.

(14) Step 1: Inversion of Ip from Near Offset Data.

(15) First, acoustic FWI is performed using near offset PP data (offset <500 m) to get an estimate of I.sub.P, which is plotted in FIG. 3. As explained above, PP-wave data at small reflection angles (equivalently, small offsets in this example) are determined by acoustic impedance I.sub.p. Elastic parameters have very little effect on small angle PP reflection data. Initial V.sub.P and density models are needed to perform acoustic FWI. The initial V.sub.P model can be derived from traditional migration velocity analysis, and for this synthetic test, a smoothed version of the true V.sub.P profile (used to forward model the synthetic data) in FIG. 2 was used. The initial density model can be derived from an empirical relationship between density and V.sub.P. For simplicity, a constant density (1,000 kg/m.sup.3) model was used to start with. From the mathematical definition
I.sub.p=V.sub.p,(1)
it is clear that inverted I.sub.P with known density can be directly translated to V.sub.p after dividing I.sub.P by density . The results at iteration 5 of I.sub.P and V.sub.P are shown in both time and depth domain in FIG. 3, where the dark lines are the inverted model and the lighter shaded lines are the synthetic model. The inverted unknown is I.sub.p in this case. An estimate of V.sub.p may then be obtained by dividing inverted I.sub.p by according to equation (1). In FIG. 3, the inverted models are overlaid with the true synthetic models for comparison. All inversions are performed in depth domain (meters); the results are shown at 11 and 12. For comparison in certain frequency range, inversion results are converted to time (seconds) by depth-to-time conversion using the smoothed version of true V.sub.p in FIG. 2. The comparisons in time domain (9 and 10) are limited within 5-40 Hz after applying band-pass filter. From 9 and 11, it can be seen that the inverted Ip matches synthetic model very well. Since V.sub.p was derived from the inverted I.sub.p based on an assumed constant according to Eqn. (1), a good match between derived V.sub.p and the true V.sub.p is not expected (no updated estimation of has been performed yet). Thus, the initial density model (constant) is very different from the synthetic density model (7 in FIG. 2), and this difference is reflected in V.sub.P due to equation (1). This is particularly indicated in 10 by the mismatch in time domain at about 1.75 s and a similar mismatch in depth domain (12) at about 1800 m. It can be seen in 9 and 11 that the mismatch for I.sub.P is much less at that particular time and depth.

(16) Data misfit 15, i.e. the difference between measured data 13 (from synthetic models) and simulated data 14 (from inverted I.sub.p, constant density and derived V.sub.p according to (1)), is shown in FIG. 3. The difference is actually negligible. Data misfit is a very important criterion for convergence check during inversion of field (actual) data because in a field data application, a true model is seldom known. Generally speaking, when other conditions are similar, better data misfit usually, but not always, indicates higher confidence in the inversion product. The negligible amount of misfit indicates that near offset data can be well explained by Ip alone.

(17) Step 2: Inversion of I.sub.S or V.sub.P/V.sub.S from Middle Offset (<2 km) Data with I.sub.P Fixed from the Previous Step.

(18) The following are known, simple relationships:

(19) $\begin{array}{cc}{I}_s=V_s& \left(2\right)\\ I_s=\frac{V_s}{V_p}{I}_p,& \left(3\right)\end{array}$
where Eqn. (3) results directly from Eqs. (1) and (2). In this step 2, the inversion needs to be elastic and the inversion unknown was V.sub.p/V.sub.s. Since I.sub.p is fixed from the previous step, inverting for V.sub.p/V.sub.s is equivalent to inverting for I.sub.s in this step according to (3). Alternatively, the inversion unknown could be I.sub.S. FIG. 4 shows difference between initial Vs model (dark line, constant) and synthetic model (lighter shaded line) in 18, and the ratio Vp/Vs is shown in 19. With this initial Vs model, and Vp (shown in 17) and density (constant) from step 1, a large data misfit may be observed in panel 22 when extending the offset to 2 km, as indicated in FIG. 4. This is because I.sub.p alone is not adequate to explain middle reflection angle (offset) data. A good estimate for the second parameter, which is V.sub.p/V.sub.s is needed to explain middle offset data. However, the data misfit at near offset is still as small as in FIG. 3 (15) because I.sub.P is fixed (16, 9) from step 1.

(20) Following the same layout as FIG. 3 used in displaying step 1 inversion results, FIG. 5 shows the inverted V.sub.P/V.sub.S (dark line, 26) after 5 iterations, overlaid with the synthetic model (lighter shaded line, 26). The inverted model matches the synthetic model very well. As indicated at panel 29, the data misfit at the middle offset range (500 m to 2 km, scale not shown in the drawing) is greatly reduced by having the benefit of the inverted V.sub.P/V.sub.S model. In the step 2 inversion, I.sub.P (23) and V.sub.P (24) are fixed from step 1. From equation (3), an accurate I.sub.s can be derived from accurate inversion results of I.sub.p and V.sub.p/V.sub.s. But V.sub.s from Eqn. (2) will not be as accurate if information on density is missing or inaccurate. This is indicted at time1.75 s in panel 25 in FIG. 5, where it can be seen that Vs derived from V.sub.p/Vs does not match synthetic model to the same degree as V.sub.p/V.sub.s.

(21) Step 3: Inversion of Density from Far Offset (Up to 5 km) Data with I.sub.P and V.sub.P/V.sub.S Fixed from the Previous Two Steps.

(22) The mathematical relations (1)-(3) indicate that any update of density with I.sub.p and V.sub.P/V.sub.S being fixed results in an update to V.sub.P and V.sub.S. Therefore, inversion of density with I.sub.P and V.sub.P/V.sub.S fixed is equivalent to inversion of V.sub.p. In step 3, all available offsets up to 5 km (in this example) are used to perform an elastic inversion for density, with I.sub.P and V.sub.P/V.sub.S fixed from steps 1 and 2. FIG. 6 shows inverted density (dark line, 33) after 10 iterations, overlaid with the synthetic model (lighter shaded line, 33), where the synthetic model is 7 in FIG. 2, converted to time domain. At the same time, step 3 results in an improved prediction of V.sub.P (31, dark line) compared with that of FIG. 3 (10, dark line) due to the updated density profile 33. Data misfit is mostly at far offsets (2 km to 5 km) as is shown at 36 in FIG. 3.

(23) The foregoing description is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. All such modifications and variations are intended to be within the scope of the present invention, as defined by the appended claims.

REFERENCES

(24) Aki and Richards, Quantitative Seismology, Theory and Methods, chapter 5.20, W. H. Freeman & Co. (1980). Lazaratos, S., Chikichev, I. and Wang, K., 2011, Improving convergence rate of Full Wavefield Inversion (FWI) using spectral shaping, PCT patent application publication WO2012/134621. Hampson, Russell, and Bankhead, Simultaneous inversion of pre-stack seismic data, 75.sup.th Annual International Meeting, SEG, Expanded Abstracts, 1633-1637 (2005). Sears, Singh and Barton, Elastic full waveform inversion of multi-component OBC seismic data, Geophysical Prospecting 56, 843-862 (2008)