Estimating multiple subsurface parameters by cascaded inversion of wavefield components

Abstract

A method, including: obtaining initial estimates of a plurality of subsurface parameters; obtaining a recorded wavefield decomposed into a plurality of discrete components; performing, with a computer, a cascaded inversion where the initial estimates of the subsurface parameters are individually updated, wherein each of the subsurface parameters are updated using a different discrete component of the recorded wavefield of the plurality of discrete components; and generating, with the computer, updated subsurface models from the cascaded inversion for each of the subsurface parameters.

Claims

1. A method, comprising: obtaining initial estimates of a plurality of subsurface parameters, wherein the subsurface parameters comprise both anisotropic parameters and physical parameters; obtaining a dataset representative of a wavefield from a seismic survey decomposed into a plurality of discrete offset ranges; determining for each subsurface parameter a discrete offset range to which each subsurface parameter is most sensitive; iteratively performing, with a computer, a cascaded inversion, wherein the iteratively performing the cascaded inversion includes: iteratively updating, for each iteration of the cascaded inversion, each subsurface parameter; expanding for each iteration of the cascaded inversion the plurality of discrete offset ranges to provide a plurality of expanded discrete offset ranges, wherein during each iteration of the cascaded inversion, an update to each iteratively updated subsurface parameter is constructed with the expanded discrete offset range that each iteratively updated subsurface parameter is sensitive to; and performing each iteration of the cascaded inversion until the plurality of expanded discrete offset ranges are each expanded to include all data within the dataset; and generating and displaying, with the computer, one or more subsurface models from the cascaded inversion for each of the subsurface parameters.

2. The method of claim 1, wherein each subsurface parameter is updated by performing a sub-iterative inversion process that comprises individually performing a predetermined number of mono-parameter full wavefield inversion updates for each subsurface parameter.

3. The method of claim 2, wherein the mono-parameter full wavefield inversion updates for each subsurface parameter are performed in series.

4. The method of claim 1, further comprising: iteratively updating each subsurface parameter until a predetermined level of convergence between the dataset and the cascaded inversion for each subsurface parameter is achieved.

5. The method of claim 1, wherein the cascaded inversion is iteratively performed until a global convergence is obtained.

6. The method of claim 1, further comprising: managing hydrocarbons based on the one or more subsurface models.

7. The method of claim 6, wherein the managing includes causing a well to be drilled at a location determined from the one or more subsurface models.

8. The method of claim 7, wherein the managing includes producing hydrocarbons from the well.

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 is a non-limiting example of a method for multi-parameter estimation via cascaded inversion of wavefield components.

(3) FIG. 2 is a non-limiting example of an application of the present technological advancement to anisotropic acoustic full wavefield inversion.

(4) FIG. 3 is a non-limiting example of the inversion of seismic data.

(5) FIGS. 4A and 4B are non-limiting examples of models useable with the present technological advancement.

(6) FIGS. 5A and 5B are non-limiting examples of inverted models derived from the present technological advancement.

(7) FIGS. 6A and 6B are non-limiting examples of pre-stack depth migration image gathers derived from using the models of FIGS. 3 and 4.

(8) FIG. 7 is a non-limiting example of a computer that can execute methods that embody the present technological advancement.

DETAILED DESCRIPTION

(9) Exemplary embodiments are described herein. However, to the extent that the following description is specific to a particular, 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) The present technological advancement described herein relates to a practical solution to multi-parameter full wavefield inversion (FWI). For example, the method depicted in FIG. 1 reformats multi-parameter FWI into a cascaded inversion of systematically selected wavefield components. As explained in more detail below, the inverted parameters can be readily constrained to parts of the recorded wavefield to which they are most sensitive while ensuring consistent data-fit of the full recorded wavefield.

(11) FIG. 1 is a non-limiting example of a method for multi-parameter estimation via cascaded inversion of wavefield components. The non-limiting examples provided herein pertain to anisotropic inversion. However, the present technological advancement can be applied to estimate any number of subsurface parameters. As used herein, subsurface parameter includes geophysical properties such as p-wave velocity, shear wave velocity, vertical velocity, normal moveout velocity, anisotropic parameters, attenuation parameters, porosity, permeability, and resistivity. This is not an exhaustive list of parameters useable with the present technological advancement.

(12) As described above, in production-scale applications, one way to account for anisotropy in FWI is by utilizing anisotropic models derived from other methods (e.g. travel-time tomography) for FWI forward simulations while keeping such models fixed throughout the FWI workflow. Here, distinctions are made between physical subsurface parameters (e.g., velocity or impedance) and anisotropic parameters (e.g., Thomsen parameters). If the anisotropic parameters derived outside the FWI engine are sufficiently accurate, it is possible to obtain physical subsurface models which are within acceptable levels of accuracy (e.g., Warner et al., 2013). However, where such anisotropic parameters are inaccurate or anisotropy is neglected, the estimated physical parameters are strongly degraded (e.g. Prieux et al., 2011). To overcome this, some authors have proposed ad-hoc methods that utilize a different engine (such as ray-based tomography) to update anisotropic parameters using FWI physical models and then repeating FWI for re-estimation of the physical parameters while keeping the updated anisotropic parameters fixed. Such alternating approach (for physical parameters through FWI and for anisotropic parameters through a different engine) may provide improvements over the conventional fixed-anisotropy approach. However, because of the disjoint between the estimation engines, it is difficult to simultaneously retrieve consistent anisotropic and physical parameters that simultaneously predict the full seismic wavefield.

(13) By posing anisotropy estimation as an integral part of the FWI formulation in the present technological advancement, reliable and consistent estimates can be obtained of both anisotropic and physical models simultaneously. Utilizing known characteristics of the seismic wavefield, along with sensitivity to changes in various physical parameters, allows multi-parameter anisotropic FWI to be reformatted into coupled cascades of mono-parameter FWI updates. By doing this, the high cost associated with full multi-parameter inversion, which requires reliable second-order (Hessian) information, can be avoided without sacrificing consistency or robustness between the estimated subsurface parameters.

(14) In step 101, a full recorded data set is obtained. The data set can be obtained by either the generation of the data using a source and receivers as is well known in the art, or by operations by a computer to select, access, download, or receive the electronic data set.

(15) In step 103, the parameters of interest are obtained. This can include the determination of discriminating properties of parameters of interest and the components of the full recorded data to which such parameters are most sensitive.

(16) In this anisotropic inversion example, the parameters could be any two or more of vertical velocity V.sub.z, normal moveout velocity V.sub.NMO, the Thomsen parameters and/or , , and horizontal velocity V.sub.H. For the present example, the parameters under consideration are V.sub.NMO and .

(17) Normal moveout or NMO is the effect of the separation between receiver and source on the arrival time of a reflection that does not dip, abbreviated NMO. A reflection typically arrives first at the receiver nearest the source. The offset between the source and other receivers induces a delay in the arrival time of a reflection from a horizontal surface at depth. A plot of arrival times versus offset has a hyperbolic shape. V.sub.NMO can be determined from the following equation:

(18) $\begin{array}{cc}t=\sqrt{t_o^2+\frac{x^2}{V_{\mathrm{NMO}}^2}}& \left(1\right)\end{array}$
wherein t is the time that the receiver at offset x receives a reflected signal, and t.sub.o is the time for a receiver at zero offset to receive a reflected signal.

(19) is an anellipicity coefficient that is defined as follows:

(20) $\begin{array}{cc}=\frac{1}{2}\left[{\left(\frac{V_H}{V_{\mathrm{NMO}}}\right)}^{-1}-1\right]& \left(2\right)\end{array}$

(21) In step 105, the full recorded data set is decomposed into predefined components. FIG. 2 provides a non-limiting example of the data set decomposition. In FIG. 2, the recorded wavefield 201 is decomposed into discrete offset ranges 203 and 205 labeled near offset and far offset, respectively, and are defined by the vertical dashed lines. Offset is the horizontal distance from the source-position (or from the center of a source group) to a seismic detector, or more commonly to the center of a detector group measured in the horizontal (x-y) plane.

(22) The near offset region 203 and far offset region 205 intentionally do not encompass the full recorded wavefield 201. As explained below, the near offset region 203 and the far offset region 205 will expand for subsequent cascades, and can eventually encompass the entire full recorded wavefield 201. Those of ordinary skill in the art will be able to determine the initial sizes of the near offset region 203 and far offset region 205 based on the particular parameterization being employed, the size of the data set, and the depth of penetration.

(23) Far offset and near offset are examples that are useful in anisotropic inversion. The full recorded data set could be decomposed differently, depending upon the particular application.

(24) In step 107a, component(s) of the recorded wavefield 201 are selected for parameter 1. This can include defining selection criteria for components of the recorded wavefield 201 to utilize for each parameter of interest.

(25) In the present example, parameter 1 is V.sub.NMO and the corresponding component of the recorded wavefield is the near offset region 203. The present technological advancement can use contributions from the corresponding components of the recorded wavefield to which a parameter is most sensitive (i.e., V.sub.NMO is more sensitive to the near offset than the far offset, or small- to mid-reflection angles). By definition, V.sub.NMO provides information about the normal-moveout characteristics of the recorded seismic wavefield 201 at short- to mid-offset ranges.

(26) In step 107b, corresponding component(s) of the recorded wavefield 201 are selected for parameter 2. In the present example, parameter 2 is and the corresponding component of the recorded wavefield is the far offset region 205. is more sensitive to the far offset (or wide reflection angles) than the near offset. By definition, provides information about the non-hyperbolic characteristics of the recorded seismic wavefield 201.

(27) In step 107n, component(s) could be selected for a nth parameter in a manner similar to what was described in steps 107a and 107b.

(28) In step 109a, iterative mono-parameter FWI updates are performed for parameter 1 (V.sub.NMO). The FWI updates can be accomplished through an iterative cost function optimization (the iterative nature being indicating the curved arrow in the box labeled 109a). Cost function optimization involves iterative minimization of the value, with respect to the model, of a cost function which is a measure of the misfit between the calculated and observed data (this is also sometimes referred to as the objective function), where the calculated data is simulated with a computer using the current geophysical properties model and the physics governing propagation of the source signal in a medium represented by a given geophysical properties model. The simulation computations may be done by any of several numerical methods including but not limited to finite difference, finite element or ray tracing. A geophysical properties model gives one or more subsurface properties as a function of location in a region. Cost functions are well known to those of ordinary skill in the art, and an example of a suitable cost function is discussed in U.S. Pat. No. 8,121,823, the entire content of which is hereby incorporated by reference.

(29) Optionally, within each cascade, the inversions for each parameter can use a frequency continuation technique (e.g., Bunks et. al, 1995; Sirgue and Pratt, 2004; Prieux et. al, 2009). The inversion for each parameter can use different frequency pass bands. For the particular example described here, because typically only the long wavelength background anisotropic model is required, one can utilize lower frequency pass bands for updates and higher frequency pass bands for V.sub.NMO updates.

(30) FIG. 3 provides a non-limiting example of iterative cost function optimization. In step 301, an initial geophysical property model for the particular parameter of interest is obtained. Examples of such initial geophysical property models are provide in FIGS. 4A and 4B, and are discussed further below. In step 303, the initial geophysical property model is used to compute synthetic data. Generating and/or obtaining synthetic data based on an initial geophysical property model is well known to those of ordinary skill in the art. In step 305, an objective function is computed using observed geophysical data and corresponding synthetic data obtained in step 303.

(31) As part of step 303, a gradient of the cost function can be taken with respect to only parameter 1 (V.sub.NMO); other parameters and properties are held constant. Where m is the parameter, f is the cost function, and is the gradient, this calculation can be expressed as follows:

(32) $\begin{array}{cc}{m}_i=\frac{f}{m_i}& \left(3\right)\end{array}$

(33) In step 307, the gradient (which provides the rate of the change of the cost function in a given direction) is then used to update parameter 1 (V.sub.NMO) in order to minimize the cost function. Step 307 can include searching for an updated geophysical property model that is a perturbation of the initial geophysical property model in the gradient direction that better explains the observed data.

(34) In step 309, the preceding steps are repeated by using the new updated model as the starting model for another gradient search. The process can continue until an updated model is found which satisfactorily explains the observed data or until a predetermined number of iterations have been completed. Commonly used cost function inversion methods include gradient search, conjugate gradients and Newton's method. However, this is not necessarily an exhaustive list, and the present technological advancement can be used with any cost function inversion method.

(35) Steps 109b to 109n are the same as step 109a, except they are performed for their respective parameters and use the model obtained from the preceding step that was updated for parameter n1. For example, the model obtained from step 109a, having been updated for only parameter 1, serves as an input for step 109b. At the end of step 109n, a set of updated models, one for each parameter, is generated.

(36) By decomposing the recorded wavefield as discussed above, it can be ensured that initial updates to each inversion parameter are constructed only using parts of the wavefield to which such parameter is most sensitive. The present technological advancement has an advantage that fewer crosstalk artifacts are present in the gradient search direction for the each parameter. This enables the avoidance of possible local minima at early stages of FWI iterations.

(37) Once the FWI mono-parameter updates are completed for at least two parameters, the process can proceed to step 111. In step 111, it is determined whether the desired convergence is obtained. Step 111 can include defining convergence criteria for each parameter within each cascade.

(38) In step 111, the desired convergence can be judged with a cost function based on the actual data used from the recorded data set 201 and the observed data. Cost function optimization methods involve computing the cost function for a population of models {M.sub.1, M.sub.2 . . . M.sub.n} and selecting a set of one or more models from the population that approximately minimizes the cost function. For example, synthetic data generated from the geophysical property models resulting from steps 109a to 109n can be compared to the subset of observed data used in steps 109a to 109n according to a similar process described in FIG. 3, except that the gradient is now taken for all parameters at the same time. This cost function optimization technique is well known to those of ordinary skill in the art. If the result of step 111 is that the desired level of convergence is not obtained, another iteration of steps 109a to 109n are repeated. Steps 109a to 109n are repeated until the desired convergence is obtained.

(39) Convergence refers to a condition that occurs during an iterative data modeling procedure when predictive output data, or a function measuring error (e.g., cost function), remains substantially the same between iterations. Convergence may be used to determine an end point for the iterative process by indicating an acceptable level of correspondence between predictive data with known data at a given point in space or time.

(40) Once the desired level of convergence is obtained, the process proceeds to step 113. Step 113 can include defining global convergence criteria for the full data set. A global cost function optimization method, referenced above, can be used here as well.

(41) The model analyzed in step 111 with respect to actual data used from the data set 200 can now be analyzed with respect to the entire data set 200. For example, synthetic data is generated from the models and is compared to observed data of the entire data set 200 according to a process that is similar to what is described in FIG. 3, except that the gradient is taken for all parameters at the same time. If a cost function based on this synthetic data and observed data has not converged, then the answer to step 113 is no and the process proceeds to step 115 in order for another cascade to be performed.

(42) In step 115, the near offset region 203 and the far offset region 205 are expanded (i.e., the dashed vertical lines in FIG. 2 are systematically moved to II in the direction of the respective arrows) in order for each region to include more of the recorded data 201. Then steps 107 to 111 are repeated for this second cascade, which will use the expanded regions 203/205. Once convergence is achieved, the process will proceed again to step 113 in order to test for global convergence.

(43) Cascades can be iteratively performed until the offset regions 203/205 are expanded from those used in a preceding cascade until each region 203/205 includes the full recorded wavefield data set 201 (i.e., until the offset regions 203/204 are systematically expanded to regions III and IV in respective iterative cascades). As shown in FIG. 2, at the time the answer to step 113 is yes, the full range of offsets can be utilized to update all inversion parameters. However, it is possible that convergence is deemed to have been satisfied before the offset regions 203/205 encompass the entire data set 203/205.

(44) FIG. 2 divides the data set 201 in a manner that corresponds to four cascades before the offset regions 203/205 each encompass the entire data set 201. Four cascades is merely an example, and the data set 201 can be divided to correspond to any number of cascades.

(45) Once global convergence is achieved, the answer to step 113 is yes and the process can proceed to step 117. The updated subsurface parameters that are obtained after step 113 are robust and consistent models that ensure a good fit to the recorded data at all offsets, and generates flattened image gathers (discussed below).

(46) In step 117, hydrocarbons are managed according to the updated subsurface parameters. 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.

(47) The following is an exemplary application of the present technological advancement to anisotropic acoustic inversion of a field 3D Narrow azimuth streamer data set. The maximum offset is approximately 5 km and the maximum (cut-off) frequency for the inversion is 15 Hz. The target area is approximately 4 km12 km. A total of 20 FWI iterations are performed comprising four cascades (See FIG. 1). Each cascade is made up of 3 iterations for V.sub.NMO and 2 iterations for . The initial models for V.sub.NMO and are shown in FIGS. 4A and 4B, respectively. In FIGS. 4A and 4B, and in similar figure, the inline (left) and cross-line (right) panels are from positions indicated by the white broken lines on the opposite opposite panels. Note that whereas the initial V.sub.NMO model was derived from conventional ray-ray-base tomography, the starting model was set to zero (consequently, FIG. 4B is entirely black).

(48) The inverted V.sub.NMO and models are shown in FIGS. 5A and 5B, respectively. The estimated V.sub.NMO model (FIG. 5A) is of much higher resolution than the starting model (FIG. 4A) and the model (FIG. 5B) is geologic and of within known realistic values.

(49) One measure of the accuracy of the estimated models is the flatness of the migrated image gathers. FIGS. 6A and 6B show the migrated gathers within the target area derived using the initial models (FIGS. 4A and 4B) and the inverted models derived using the present technological method (FIGS. 5A and 5B). Note that the gathers derived using the present technological advancement produce image gathers that are flatter at all offsets (FIG. 6B), indicating that these models are of superior quality to the starting models (FIG. 4A) derived from ray-based tomography.

(50) While the above-embodiments described geophysical inversion based on cost function optimization, the present technological advancement can be used with other techniques, such as iterative series inversion. Iterative series inversion, based on iteration of the Lippmann-Schwinger equation is discussed in U.S. Pat. No. 8,121,823.

(51) FIG. 7 is a block diagram of a computer system 2400 that can be used to execute the present geophysical inversion techniques. A central processing unit (CPU) 2402 is coupled to system bus 2404. The CPU 2402 may be any general-purpose CPU, although other types of architectures of CPU 2402 (or other components of exemplary system 2400) may be used as long as CPU 2402 (and other components of system 2400) supports the operations as described herein. Those of ordinary skill in the art will appreciate that, while only a single CPU 2402 is shown in 7, additional CPUs may be present. Moreover, the computer system 2400 may comprise a networked, multi-processor computer system that may include a hybrid parallel CPU/CPU system. The CPU 402 may execute the various logical instructions according to various teachings disclosed herein. For example, the CPU 2402 may execute machine-level instructions for performing processing according to the operational flow described.

(52) The computer system 2400 may also include computer components s as non-transitory, computer-readable media. Examples of computer-readable media include a random access memory (RAM) 2406, which may be SRAM, DRAM, SDRAM, or the like. The computer system 2400 may also include additional non-transitory, computer-readable media such as a read-only memory (ROM) 2408, which may be PROM, EPROM, EEPROM, or the like. RAM 2406 and ROM 2408 hold user and system data and programs, as is known in the art. The computer system 2400 may also include an input/output (I/O) adapter 2410, a communications adapter 2422, a user interface adapter 2424, and a display adapter 2418.

(53) The I/O adapter 2410 may connect additional non-transitory, computer-readable media such as a storage device(s) 2412, including, for example, a hard drive, a compact disc (CD) drive, a floppy disk drive, a tape drive, and the like to computer system 2400. The storage device(s) may be used when RAM 2406 is insufficient for the memory requirements associated with storing data for operations of the present techniques. The data storage of the computer system 2400 may be used for storing information and/or other data used or generated as disclosed herein. For example, storage device(s) 2412 may be used to store configuration information or additional plug-ins in accordance with an the present techniques. Further, user interface adapter 2424 couples user input devices, such as a keyboard 2428, a pointing device 2426 and/or output devices to the computer system 400. The display adapter 2418 is driven by the CPU 2402 to control the display on a display device 2420 to, for example, present information to the user regarding available plug-ins.

(54) The architecture of system 2400 may be varied as desired. For example, any suitable processor-based device may be used, including without limitation personal computers, laptop computers, computer workstations, and multi-processor servers. Moreover, the present technological advancement may be implemented on application specific integrated circuits (ASICs) or very large scale integrated (VLSI) circuits. In fact, persons of ordinary skill in the art may use any number of suitable hardware structures capable of executing logical operations according to the present technological advancement. The term processing circuit encompasses a hardware processor (such as those found in the hardware devices noted above), ASICs, and VLSI circuits, Input data to the computer system 2400 may include various plug-ins and library files. Input data may additionally include configuration information.

(55) The present techniques may be susceptible to various modifications and alternative forms, and the examples discussed above have been shown only by way of example. However, the present techniques are not intended to be limited to the particular examples disclosed herein. Indeed, the present techniques include all alternatives, modifications, and equivalents falling within the spirit and scope of the appended claims.

REFERENCES

(56) Each of the following references are hereby incorporated by reference in their entirety: 1. Bunks, C., F. M. Saleck, and S. Zaleski and G. Chavent (1995). Multiscale seismic waveform inversion: Geophysics, 60, 1457-1473. 2. Prieux, V., Operto, S., Brossier, R., and Virieux, J. (2009). Application of acoustic Full waveform inversion to the synthetic Valhall velocity model, SEG Expanded abstracts, 2268-2272. 3. Prieux, V., R. Brossier, Y. Gholami, S. Operto, J. Virieux, O. I. Barkved, and J. H. Kommedal (2011). On the footprint of anisotropy on isotropic full waveform inversion: The Valhall case study: Geophysical Journal International, 187, 1495-1515. 4. Sirgue, L., and Pratt, R. G. (2004). Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies: Geophysics, 69, 231-248. 5. Warner, M., Ratcliffe, A., Nangoo, T., Morgan, J., Umpleby, A., Shah, N., Vinje, V., tekl, I., Guasch, L., Win, C., Conroy, G., and Bertrand, A. (2013). Anisotropic 3D full-waveform inversion. GEOPHYSICS, 78(2), R59-R80.