Title: Method of seismic inversion using artificial kinematic constraints

Technical Field

[0001] The present invention relates to geophysics, and more particularly seismic prospecting or detecting and processing of seismic data for analysis and interpretation.

Background Art

[0002] In geophysical exploration, seismic data is used to image the subsurface where geological and civil structures might be identified, including minerals, hydrocarbons, and other fluids. In seismic exploration, one or more sources emit elastic waves in the form of modulation of pressure or ground motion. They propagate from specific locations (wavefield) that can be on or below the surface or sea surface, or in a borehole. The wavefield propagates away from the source(s) through the subsurface. Properties of the rocks can be estimated from the wavefield response of the media contrast (such as acoustic impedance). The response can be detected/recorded at the surface, borehole or other receiver locations through the reflected wavefield as pressure, partial motion or some derivative quantities.

[0003] There is no known solution for rock properties estimation from seismic trace at normal incident wavefield response - stacked seismic data (post-stack) or seismic trace that is generated or obtained by having source and receiver at the same position – that would allow to estimate independently velocity and density properties of rocks. Known solutions are the seismic inversion algorithms or methods for normal incident wavefield response. These methods and algorithms allow inverting or estimating velocity and density functions as a product of the two functions and is normally called acoustic impedance or p-impedance.

[0004] Other existing inversion solutions allow to get velocity and density performing an inversion of only seismic offset traces or angle traces also called as seismic gather (pre-stack data) resolving for angle dependent reflectivity or velocity and density of linearized approximations. There are two main types of pre-stack inversion, deterministic and stochastic. It is usually set as least squares problem that can be described as model parameters = generalized inverse x observations. Deterministic inversion tries to achieve production of a single “best” solution. The stochastic inversion produces many possible good solutions that are average to the deterministic solution. Density is resolved poorly because it is linked to impedance and solutions assumed that the seismic offset is long enough to reach 45 degrees angle between sours and receiver from a reflected point - often is not available due to seismic data acquisition or seismic data processing issues.

[0005] For example, the latest, most accurate seismic inversion algorithms that are used in the industry worldwide, is a Bayesian-based inversion. (Reference: Buland A. and H. Omre Bayesian linearized AVO inversion. Geophysics, 2003.63(1)). However, it is stated in the article, that density is not reliable: «With realistic noise levels, acoustic impedance was the best determined parameter, while the inversion provided practically no information about the density.” The major seismic inversion software providers supply deterministic, stochastic approaches, or approaches with Bayesian-based extension where density is linked to impedance. It can only be resolved when inverting pre-stack seismic data.

[0006] US patent nr. US10067264B2 by Blanchard et al. disclose a method a seismic inversion for petrophysical properties of a subsurface volume comprising the steps of: obtaining petrophysical data relating to valid geological and/or dynamical scenarios, converting this data into valid combinations of elastic parameters; projecting the valid combinations of elastic parameters onto a spherical plot; and determining a penalty term from the distances between each cell of the spherical plot and the nearest valid combination of elastic parameters within the subsurface volume. Valid geological and/or dynamical scenarios comprise those which are petrophysically possible. The penalty term is then used to constrain an inversion minimizing a cost function associated with seismic mismatch between two or more seismic surveys. This method is limited to prestack data.

[0007] Further, US Application published nr. US20170108602A1 by Di Yang discloses a Full Wavefield Inversion (FWI) method to estimate subsurface properties. The Full Wavefield Inversion (FWI) requires seismic data (traces) recorded at the different angles and uses a mathematical approach (a Fourier transform, discrete Fourier transform, or a fast Fourier transform) through misfit function. According to this method only P-wave velocity (Vp) can be estimated. Further this method of relies on obtaining a seismic dataset that is separated into subsets according to predetermined subsurface reflection angle ranges and inverting for density models. Where each of the data subsets generated, acoustic FWI is applied to obtain acoustic impedances. All inversions can start from the same velocity model, and the kinematics are not updated in this process.

[0008] It is vital to evaluate rock properties such as velocity and density to explore for mineral resources. Rock properties are estimated, for example, through acquiring seismic data, with further processing and inversion. Accordingly, there is need for a method that can estimate reliably velocity (Vp) and (Dn) separately from post-stack data as well as to being able to reliably compute non-coupled density for pre-stack data.

Brief summary of the disclosure

[0009] A computer implemented method for extracting or estimating rock properties from seismic traces is disclosed. This method assists in understanding the interpretation of post-stack and pre-stack seismic data by predicting density and velocity from subsurface seismic data.

[0010] The computer implemented method comprising the follow steps:

(1) preparing an initial model with initial functions of velocity and density and synthetic trace(s),

(2) generating synthetic trace(s) from the velocity and density function(s) and with added artificial kinematic constraints and randomly updating the initial functions of velocity and density in random start and length time or depth window(s);

(3) creating updated synthetic traces using randomly updated velocity and density functions; for each iteration, artificial wave(s) traveling from a source point to a reflection point and back to any receiver are simulated as a constraint; (4) performing a search of the misfit object function of any norm between original/or real trace or traces and the trace(s) generated in the following iteration, and

(5) using probabilistic techniques for approximating the global optimum and minimizing the cost function associated with seismic trace(s) mismatch.

[0011] Furthermore, a computer program is disclosed comprising instructions which, when the program is executed by a computer, cause the computer to carry out the steps of the computer implemented method.

[0012] The invention is defined by the claims.

Description of the drawings

[0013] Fig. 1 Outline of existing seismic post-stack and pre-stack inversion principal.

[0014] Fig. 2 Outline of computer implemented method.

Detailed description

[0015] The method and its different embodiments will be described in detail referring to the enclosed figures.

[0016] In the case of post-stack seismic data or seismic trace(s) at zero-incident (real trace): simulating offset seismic traces iteratively using artificial kinematic constraints and random updates of the velocity and density functions, and searching for an optimal misfit function (when some optimal minimum is reached) between simulated traces (gather) and the real trace, accepting the simulated parameters of density and velocity. In the case of pre-stack seismic data or gathers (real traces): perform the same steps as for post-stack, with the exception that the optimal misfit function is searched between artificial offset trace(s) and the real offset/or angle trace(s). The new method allows inverting seismic stacked data or seismic trace of a zero-incident (post-stack) for velocity, density and impedances. Before this method, inversion on stacked data was able to resolve impedance only, where impedance is velocity multiplied by density. This method is also applicable for seismic pre-stack data. Seismic pre-stack data is a gather that collects several traces and shares some common geometric attribute: common source or receiver, midpoint, image point, offset, azimuth, angle, dipping reflector gathers. Applying algorithm to the seismic pre-stack data, allows increasing quality on resulting velocity, density and impedance.

[0017] There are no known solutions for estimating Vp and Dn separately from poststack data, in other words, with all known physics it is not possible to solve this problem. The solution achieved here does not use the “classic” physics as we get used to understanding for a plane wave propagation and seismic inversion problem. It uses algorithms from artificial intelligence (AI) families and constructed assumption, that let us solve this problem.

[0018] The new approach does not use any of known methods: deterministic and stochastic, or geostatistical. It does not use any linearized form of any amplitude versus offset or amplitude versus angle equations. The basis of the existing algorithm are Zoepprits, Aki-Richards (1980s) and Fatti et al. (1994) equations and theirs improved approximations or modifications. (Aki, K., and Richards.1980. “Quantitative seismology: Theory and method.” (W. H. Freeman and Co.); Aki, K., B.Chouet.1975. “Origin of coda waves: source, attenuation, and scattering effects.” (J. Geophys. Res.) 80: 3322-3342; Zoeppritz, Karl.1991. “VIIb. Über Reflexion und Durchgang seismischer Wellen durch Unstetigkeitsflächen.” (Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-physikalische Klasse) 66-84; Fatti Jan L., George C. Smith, Peter J. Vail, Peter J. Strauss, and Philip R. Levitt Detection of gas in sandstone reservoirs using AVO analysis. GEOPHYSICS.1994.9(59).362-1376).

[0019] Extension of the proposed method to the non-normal incident environment allows estimate parameters: at least compressional wave velocity (Vp), density (Dn), and shear wave velocity (Vs) and their derivatives, with the difference that resulted or estimated density is estimated independently of other parameters. A coupling between velocity and density parameters can be only in some cases of initial models. Thus, the proposed method allows producing a result of rock properties estimation that is more accurate, computationally efficient, and does not require specific extended educational skills for an operator, which executes the algorithm on the suitable computer apparatus as it would be needed for any other seismic inversion process.

[0020] Seismic amplitudes contain information about the geological and physical properties of the rocks 100. A seismic is performed to extract potential properties. The mathematical methods allow inverting the seismic amplitudes to physical properties of the rocks, Fig 1. Thus, post-stack inversion methods (model-based, recursive, sparse spike, colored) allows inverting seismic data to rock properties of only p-impedance 105; pre-stack inversion methods with the same mathematical approaches allow extending inversion of seismic data to rock properties such as elastic impedance, Lame parameters, compressional-velocity (Vp), shear-velocity (Vs) and density (ρ) 110.

[0021] Pre-stack seismic inversion is an ill-posed problem that gives equivalent solutions for a large ensemble of rock properties. Other integrated and extended seismic inversion approaches of processes seismic data using multivariate or Bayesian statistic allow at some level to reduce the uncertainties. The parameter that is least accurately predicted is density, due to linearized dependency on impedance and sensitivity to the seismic noise.

[0022] To illustrate the basic definition of the problem formulation, the amplitudes of the seismic gather can be presented as a linearization of amplitude versus offset equation. In addition, pre-stack seismic data can be written as function of the reflectivity of acoustic impedance (RP_imp), shear impedance reflectivity (RS_imp), and density reflectivity (RD):

where

[0023] The solution for the parameters can be found by inverting separately each of the extracted reflectivity components at time t and trace N 110, using least-squares inversion and post-stack type scheme 100; or simultaneously, through derivative operator 115, for example, through a least-square problem: modelled parameters = generalised inverse x observation. The solution also can be presented deterministically – single best solutions, or stochastically – many possible good solutions that are averaged to the deterministic one. Bayesian stochastic inversion approach produces reasonable results however has difficulty in converging to an answer. A benefit of stochastic inversion based on multivariate Bayesian statistics is the possibility of getting the conditional probability for the answer choice; however, even a small level of noise density solution is not reliable. Thus, existing seismic inversion approaches resolve density only as a coupled parameter and from the inverting only a pre-stack data.

[0024] Fig. 2 shows a flowchart of the preferred embodiment of the method for estimating of Vp and ρ independently (the coupling can happen only beyond the frequency on the initial model), including the approach of these parameters estimation from the post-stack seismic data.

[0025] In the first aspect there is provided a method for compressional wave and density estimation at normal incident of wave field response or stacked seismic data (post-stack), or any given trace or seismic trace when seismic wave source and receiver at the same position. Here, the seismic source and receiver can be considered at nearby positions that the error when creating post-stack seismic trace from summing several traces or projecting trace(s) to the zero-incident location is acceptable.

[0026] Thus, it is necessary to generate seismic trace(s) 205 by simulating or modelling any artificial wave(s) travelling from any source point to any reflection point and back to any receiver point. It can be considered as simulating a seismic gather with any offset or angle. The travelling to destinations can be defined differently; for example, a receiver can be located at any angle from the travelling wave's reflection point; or the receiver and source can be propagated to the same position from any other offset or angle.

[0027] To create the artificial waves and further trace(s), it is necessary to define initial functions of seismic velocity and density for each location of the real seismic trace – the trace or traces participating in the inversion process. The function(s) can be constructed, taken as constants or measured 200.

[0028] The generated synthetic trace(s) are used to define the penalty term 225. The synthetic trace(s) is “compared” to the real or seismic trace through performing a search of the misfit object function, using probabilistic techniques for approximating the global optimum and minimizing the cost function - trace(s) mismatch. In optimization definition, a cost function is a function that maps an event or the values of one or more variables to a real number that intuitively represents some "cost" associated with the event. The optimization problem is aimed at minimizing the cost function. This can be done on the point-to-point-base, window-based of the whole trace-based. Or rather than computing minimum distance, an average distance to all valid points or to an average distance to the some or a group of nearest points can be used. When the penalty term is in compliance or the misfit function criteria is reached, the velocity and density function used in the synthetic trace(s) generating process output as found in the inversion process 230.

[0029] The generating synthetic trace(s) 215 is done iteratively 205, the first one or group of traces can be constructed from the initial functions 200 of at least velocity and density, the next synthetic trace(s) generated from the at least velocity and at least density function(s) randomly updated in random start and length time or depth window(s) and added artificial and/or kinematic constraints 220. Each generated synthetic trace or group of traces go through the penalty term 225. If penalty criteria 225 is not reached, the iteration is repeated. This is used as constraint in the inversion process for each iteration. Another option is that mismatch estimation can be associated with to the alignment offset or angle seismic traces at each or several time event (reflecting point(s)) of a seismic wave.

[0030] The update of initial functions of at least velocity and density can be random with a random start, and random length of time or depth window. The updated synthetic trace(s), created using randomly updated at least velocity and density functions, is performed for each iteration.

[0031] In another embodiment the method is extended to non-normal incident wave field response - the offset or angle trace(s) of pre-stack seismic data (often called seismic gather). In this case, all processes showed on Fig.2 per the first aspect is repeated. In this case, the penalty term is defined by “comparison” generated offset or angle trace(s) with the real offset or angle trace(s). The generated traces as well as per the first aspect can be the offset traces or in any way summed up group(s) of offset or angle traces. Shear wave velocity is included in the initial functions. It can be acquired or computed or constructed shear wave velocity.

[0032] In another example, the method is applied to seismic pre-stack data. Seismic pre-stack data is a gather that collects several traces and share some common geometric attribute: common source or receiver, midpoint, image point, offset, azimuth, angle, dipping reflector gathers. Applying the algorithm to seismic pre-stack data, allows increasing quality on resulting velocity, density and impedance.

[0033] One or more steps of the methods and the concepts described herein may be embodied in the form of computer readable instructions for running on suitable computer apparatus, or in the form of a computer system comprising at least a storage means for storing program instructions embodying the concepts described herein and a processing unit for performing the instructions. Such a computer system may also comprise a display unit and one or more input/output devices.