Joint full wavefield inversion of P-wave velocity and attenuation using an efficient first order optimization

Abstract

A method for iteratively inverting seismic data to jointly infer a model for at least P-wave velocity and attenuation parameters of the subsurface, the method including: jointly inverting the P-wave velocity and attenuation parameters with an iterative visco-acoustic full wavefield inversion process, wherein the iterative visco-acoustic full wavefield inversion process includes computing a gradient of an objective function, the objective function measuring a misfit between all or part of the seismic data and corresponding model-simulated seismic data; for each of the P-wave velocity and attenuation parameters, computing a search direction in model space from the gradient; determining line search step sizes and for the search directions for the P-wave velocity and attenuation parameters, respectively, wherein a ratio of the step sizes is a function of the P-wave velocity parameter; and using the step sizes and and the search directions for each of the P-wave velocity and attenuation parameters, computing a new search direction in model space, then performing a line search along the new search direction to arrive at a new step size, and using the new step size and the new search direction to generate an updated model for a current iteration of the iterative visco-acoustic full wavefield inversion process.

Claims

1. A method for iteratively inverting seismic data to jointly infer a model for at least P-wave velocity and attenuation parameters of the subsurface, said method comprising: jointly inverting the P-wave velocity and attenuation parameters with an iterative visco-acoustic full wavefield inversion process, wherein the iterative visco-acoustic full wavefield inversion process includes computing a gradient of an objective function, the objective function measuring a misfit between all or part of the seismic data and corresponding model-simulated seismic data; for each of the P-wave velocity and attenuation parameters, computing a search direction in model space from the gradient; determining line search step sizes and for the search directions for the P-wave velocity and attenuation parameters, respectively, wherein a ratio of the step sizes is a function of the P-wave velocity parameter; and using the step sizes and and the search directions for each of the P-wave velocity and attenuation parameters, computing a new search direction in model space, then performing a line search along the new search direction to arrive at a new step size, and using the new step size and the new search direction to generate an updated model for a current iteration of the iterative visco-acoustic full wavefield inversion process; wherein the step sizes satisfy $\frac{}{} - \frac{}{V_{p} \log (\frac{}{_{0}})},$ V.sub.p is the P-wave velocity parameter, is frequency and .sub.o is a reference frequency.

2. The method of claim 1, wherein the obtaining the starting P-wave velocity model includes determining, with a computer, the starting P-wave velocity model from an iterative acoustic full wavefield inversion process, wherein the iterative acoustic full wavefield inversion process includes using a cross-correlation objective function to minimize a misfit between all or part of the seismic data and corresponding model-simulated seismic data.

3. The method of claim 1, wherein the objective function is an L2 objective function.

4. The method of claim 1, wherein the objective function is a cross-correlation objective function.

5. The method of claim 1, further comprising: generating a subsurface image from the updated model; and using the subsurface image to prospect for hydrocarbons.

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) The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

(3) FIG. 1 illustrates an exemplary method for jointly inverting P-wave velocity and attenuation;

(4) FIG. 2A illustrates an example of an inverted P-wave velocity model obtained from an acoustic FWI process;

(5) FIG. 2B illustrates the true P-wave velocity model for FIG. 2A;

(6) FIG. 2C illustrates the difference between FIGS. 2A and 2B;

(7) FIG. 3A illustrates P-wave velocity model obtained from a viscoacoustic FWI process;

(8) FIG. 3B illustrates the true P-wave velocity model for FIG. 3A;

(9) FIG. 3C illustrates the difference between FIGS. 3A and 3B;

(10) FIG. 4A illustrates an example of an inverted 1/Qp model obtained from a viscoacoustic FWI process;

(11) FIG. 4B illustrates the true 1/Qp model for FIG. 4A; and

(12) FIG. 5 illustrates an example of the determination of a new search direction in model space.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

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

(14) The present technological advancement can invert for both Vp and Qp with mitigated crosstalk between the two parameters. This can be accomplished with an efficient first-order optimization method for joint FWI of Vp and Qp. However, the present technological advancement is not limited to first-order optimization methods.

(15) FIG. 1 illustrates an exemplary method utilizing the present technological advancement. Step 101 includes inverting for Vp and/or anisotropic parameters (e.g., Thomsen parameters , , and ) using acoustic FWI with a cross-correlation objective function. Step 101 does not have to take attenuation into account.

(16) The delta parameter, the short offset effect, captures the relationship between the velocity required to flatten gathers (the normal moveout velocity) and the zero-offset average velocity as recorded by checkshots. The epsilon parameter, the long offset effect, is the fractional difference between vertical and horizontal P velocities; i.e., it is the parameter usually referred to as the anisotropy of a rock. The gamma parameter, the shear wave effect, relates a horizontal shear wave with horizontal polarization to a vertical shear wave.

(17) FWI is a computer-implemented geophysical method that is used to invert for subsurface properties such as velocity or acoustic impedance. The crux of any FWI algorithm can be described as follows: using a starting subsurface physical property model, synthetic seismic data are generated, i.e. modeled or simulated, by solving the wave equation using a numerical scheme (e.g., finite-difference, finite-element etc.). The term velocity 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 or another physical property in a cell, where a subsurface region has been conceptually divided into discrete cells for computational purposes. The synthetic seismic data are compared with the field seismic data and using the difference between the two, an error or objective function is calculated. Using 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 generate a new objective function. This process is repeated until the objective function is satisfactorily minimized and the final subsurface model is generated. A global or local optimization method is used to minimize the objective function and to update the subsurface model.

(18) The cross-correlation objective function is commonly used in seismic inversion to match the phase of the data, and it is often considered to be robust when precise amplitudes cannot fit the simulation physics. A typical normalized cross-correlation objective function is given by:

(19) $= \frac{_{obs} .Math._{calc}}{.Math._{obs} .Math. .Math._{calc} .Math.}$
where .sub.obs is the measured data, .sub.calc is the calculated data, and custom character is the non-zero lag cross-correlation operator.

(20) The cross-correlation objective function is used in step 101 with the aim of fitting the travel time of the data. With step 101, the process is able to fit the travel time of the data to a large extent, but with an incorrect velocity model. FIG. 2A illustrates an example of the inverted Vp model after step 101, together with the true Vp model (FIG. 2B) and the difference between the left and middle panels (FIG. 2C). Compared with the true Vp model, the inverted velocity is lower in the strong attenuation area. This is because FWI puts all attenuation efforts into the P-wave velocity. In FIG. 2A, the inverted Vp has the correct geometry, but the values are lower than the true velocity model.

(21) In step 102, starting from the velocity model obtained in step 101 or another starting velocity model, joint inversion for Vp and Qp is carried out using a viscoacoustic FWI process with or without updates to or inverting for the anisotropic parameters. In step 102, an L2 objective function can be used in order to fit the amplitude of the data. The L2 objective function may be preferred when the initial Qp model is merely background (i.e., the initial Qp model does not include an indication of an anomaly). If the initial Qp model does include an indication of an anomaly, then step 101 can be omitted and step 102 can be implemented with the cross-correlation objective function.

(22) FIG. 3A shows the inverted Vp model after step 102, the true Vp model (FIG. 3B) and the difference between the left and middle panels (FIG. 3C). It is observed from FIG. 3A that the inverted Vp from step 102 is close to the true VP after the joint inversion of Vp and Qp. The inverted Qp from step 102 is also close to the true Qp as is shown in FIGS. 4A and 4B, wherein FIG. 4A shows the inverted 1/Qp model after step 102 and FIG. 4B shows the true 1/Qp.

(23) Local cost function optimization of seismic data in the acoustic approximation is a common geophysical inversion task, and is generally illustrative of other types of geophysical inversion. When inverting seismic data in the acoustic approximation, the cost function can be written as:
custom character ()=.sub.g=1.sup.N.sup.g.sub.r=1.sup.N.sup.r.sub.t=1.sup.N.sup.tW(.sub.calc(,r,t,w.sub.g).sub.obs(r,t,w.sub.g)),
where
()=cost function,
=vector of N parameters, (.sub.1, .sub.2, . . . .sub.N) describing the subsurface model,
g=gather index,
w.sub.g=source function for gather g which is a function of spatial coordinates and time, for a point source this is a delta function of the spatial coordinates,
N.sub.g=number of gathers,
r=receiver index within gather,
N.sub.r=number of receivers in a gather,
t=time sample index within a trace,
N.sub.t=number of time samples,
W=norm function (minimization function, e.g. for least squares function (x)=x.sup.2),
.sub.calc=calculated seismic data from the model ,
.sub.obs=measured seismic data (pressure, stress, velocities and/or acceleration).

(24) The gathers, data from a number of sensors that share a common geometry, can be any type of gather (common midpoint, common source, common offset, common receiver, etc.) that can be simulated in one run of a seismic forward modeling program. Usually the gathers correspond to a seismic shot, although the shots can be more general than point sources. For point sources, the gather index g corresponds to the location of individual point sources. This generalized source data, .sub.obs, can either be acquired in the field or can be synthesized from data acquired using point sources. The calculated data .sub.calc on the other hand can usually be computed directly by using a generalized source function when forward modeling.

(25) FWI attempts to update the discretized model such that custom character () is a minimum. This can be accomplished by local cost function optimization which updates the given model .sup.(k) as follows:
.sup.(i+1)=.sup.(i)+.sup.(i)S(.sup.(i)),
where i is the iteration number, is the scalar step size of the model update, and S() is the search direction. For steepest descent, S()=.sub. custom character (), which is the negative of the gradient of the misfit function taken with respect to the model parameters. In this case, the model perturbations, or the values by which the model is updated, are calculated by multiplication of the gradient of the objective function with a step length , which must be repeatedly calculated.

(26) For second-order optimization techniques, the gradient is scaled by the Hessian (second-order derivatives of objective function with respect to the model parameters). The computation of .sub. custom character () requires computation of the derivative of () with respect to each of the N model parameters. N is usually very large in geophysical problems (more than one million), and this computation can be extremely time consuming if it has to be performed for each individual model parameter. Fortunately, the adjoint method can be used to efficiently perform this computation for all model parameters at once (Tarantola, 1984).

(27) FWI generates high-resolution property models for prestack depth migration and geological interpretation through iterative inversion of seismic data (Tarantola, 1984; Pratt et al., 1998). With increasing computer resources and recent technical advances, FWI is capable of handling much larger data sets and has gradually become affordable in 3D real data applications. However, in conventional FWI, the data being inverted are often treated as they were collected in an acoustic subsurface medium, which is inconsistent with the fact that the earth is always attenuating. When gas clouds exist in the medium, the quality factor (Q) which controls the attenuation effect plays an important role in seismic wave propagation, leading to distorted phase, dim amplitude and lower frequency. Therefore, conventional acoustic FWI does not compensate the Q effect and cannot recover amplitude and bandwidth loss beneath gas anomalies.

(28) Viscoacoustic FWI, on the other hand, uses both the medium velocity and the Q values in wave field propagation. Thus, the Q-effect is naturally compensated while wave-front proceeds. In some cases, where shallow gas anomalies overlay the reservoir, severe Q-effect screens off signals and causes cycle skipping issue in acoustic FWI implementation.

(29) Step 102 can include sub-steps 102.1, 102.2, 102.3, and 102.4. Sub-step 102.1 includes computing the gradient of the objective function, which is discussed above. Sub-step 102.2 includes determining the search direction, which is discussed above.

(30) Sub-step 102.3 includes appropriately scaling the search directions for Vp and Qp using the present technological advancement. The present technological advancement can parameterize the model by =(V.sub.p, 1/Q.sub.p). Anisotropy parameters have been omitted for ease of notation. In gradient based optimization methods, the models are updated along the normalized search direction S.sub.V.sub.p and S.sub.1/Q.sub.p
V.sub.p=V.sub.p+S.sub.V.sub.p,
1/Q.sub.p=1/Q.sub.p+S.sub.1/Q.sub.p,
where and are step sizes of the line search.

(31) The observed data and the current subsurface model are used in computing the gradient of the objective function for each parameter class, resulting in a search direction for each parameter class. The relative scale / decides the relative parameter updates between Vp and 1/Qp. In the conventional first-order optimization method, /=1.

(32) The line search directions and the objective function contours are illustrated in FIG. 5, wherein S.sub.new=S.sub.V.sub.p+S.sub.1/Q.sub.p. It has been observed that the combined search direction S, obtained with conventional technology, does not point to the location of the minimal objective function value which is indicated by a star in the figure. This implies that the conventional first-order optimization methods converge slowly. After applying an appropriate value of /, the new combined search direction S.sub.new points to the location of the minimal objective function value, leading to a much faster convergence than the conventional method.

(33) The appropriate relative scale of /, can be determined as follows. According to constant Q model, the dispersion relation is:

(34) $V_{p}^{att} () = V_{p} (_{0}) [1 + \frac{1}{Q_{p}} \log (\frac{}{_{0}})],$
where V.sub.p.sup.att() is the attenuated phase velocity at frequency , V.sub.p(.sub.0) is the reference P-wave velocity at a reference frequency .sub.0. We would like to maintain V.sub.p.sup.att() while updating V.sub.p and 1/Q.sub.p, because after the first step of FWI, we have already fitted the travel time of the data to a large extent. The updates of V.sub.p and 1/Q.sub.p should thus follow the relation

(35) $(\frac{1}{Q_{p}}) = - \frac{V_{p}}{V_{p} A} (1 + \frac{A}{Q_{p}}),$
where

(36) $A = \log (\frac{}{_{0}}) .$
Therefore, the relative scales between the line search step sizes are determined by

(37) $\frac{}{} = - \frac{1}{V_{p} A} (1 + \frac{A}{Q_{p}}) - \frac{}{V_{p} \log (\frac{}{_{0}})} .$

(38) Advantageously, the ratio / reduces down to a representation based only on Vp, and not Qp.

(39) A preferred approach to determine the direction of S.sub.new may be to assume that either or has a predetermined value (e.g., equal to 1) and then using the above-noted relationship to determine the other.

(40) Sub-step 102.3 could alternatively include the efficient line search methods for general multi-parameter FWI (Tang and Ayeni, 2015), in which the relative scales between the line search step sizes / are computed using trial line searches. Although the results are similar, for each iteration of FWI, the computational cost of the present technological advancement is at least 30% lower than that of the method of Tang and Ayeni (2015).

(41) Sub step 102.4 includes performing a line search and updating the models. As step 102.3 determined the direction of S.sub.new, the size of the step along that direction can be determined using a line search method. A line search involves model-simulating seismic data for various different step sizes along the search direction (e.g., S.sub.new) in the highly multi-dimensional model parameter space. The step size that minimizes the misfit between the model-simulated data and the measured data is selected.

(42) At step 103, the process is repeated for another iteration unless updated models have converged according to some predetermined convergence criterion or has reached another stopping point.

(43) At step 104, the final inverted models are generated.

(44) At step 105, the final models can be used to generate images of the subsurface that are useable to prospect for and/or manage hydrocarbons. 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.

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

(46) The foregoing description is directed to particular embodiments of the present technological advancement 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.

(47) The following references are incorporated by reference in their entirety: Bai, J., Yingst, D., Bloor, R. and Leveille, J., 2014. Viscoacoustic waveform inversion of velocity structures in the time domain. Geophysics, 79(3), pp. R103-R119. Cheng, X., Jiao; Cheng, X., Jiao, K., Sun, D. and Vigh, D., 2015, December. A New Approach of Visco-Acoustic Waveform Inversion in the Time Domain. In 2015 SEG Annual Meeting. Society of Exploration Geophysicists; Denli, H., Akcelik, V., Kanevsky, A., Trenev, D., White, L. and Lacasse, M., 2013. Full-wavefield inversion for acoustic wave velocity and attenuation. In 2013 SEG Annual Meeting. Society of Exploration Geophysicists; Malinowski, M., Operto, S. and Ribodetti, A., 2011. High-resolution seismic attenuation imaging from wide-aperture onshore data by visco-acoustic frequency-domain full-waveform inversion. Geophysical Journal International, 186(3), pp. 1179-1204; Plessix, R. E., Stopin, A., Kuehl, H., Goh, V. and Overgaag, K., 2016. Visco-acoustic full waveform inversion. In 78th EAGE Conference and Exhibition 2016. European Association of Geoscientists and Engineers; Prieux, V., Operto, S., Virieux, J. and Brossier, R., 2012. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from Valhall. Part 1: imaging compressional wavespeed, density and attenuation. Geophysical Journal International, 194(3), pp. 1640-1664; Tang, Y. and Ayeni, G., 2015. Efficient line search methods for multi-parameter full wavefield inversion. US patent application, US2015/0323689; Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 1259-1266; US Patent Publication 2015/0362622; and U.S. Provisional Patent Application 62/236,190.

Joint full wavefield inversion of P-wave velocity and attenuation using an efficient first order optimization

Assignee

Inventors

Cpc classification

Classification Explorer

G01V1/368

PHYSICS

Classification Explorer

G01V1/303

PHYSICS

Classification Explorer

G01V1/282

PHYSICS

Classification Explorer

G06F17/11

PHYSICS

International classification

Classification Explorer

G01V1/28

PHYSICS

Classification Explorer

G01V1/36

PHYSICS

Classification Explorer

G01V1/30

PHYSICS

Classification Explorer

G06F17/11

PHYSICS

Abstract

Claims

Description