Computer-implemented method and system employing nonlinear direct prestack seismic inversion for poisson impedance

Abstract

A computer-implemented method, and system implementing the method, are disclosed for computing a final model of elastic properties, using nonlinear direct prestack seismic inversion for Poisson impedance. User inputs and earth-model data is obtained over points of incidence of a survey region, at various angles of incidence. Various models are then computed that serve for lithology identification and fluid discrimination and take part in preliminary seismic exploration and reservoir characterization. Therefore, further refinement of these models is required due to changes in burial depths, compaction and overburden pressure, as they provide limitations for reservoirs on porous media. The further refinement using nonlinear direct prestack seismic model is performed on a system computer, which produces a final model of elastic properties. This model can then be applied for lithology prediction and fluid detection to identify potential targets of oil and gas exploration and estimating spots in unconventional shale gas applications.

Claims

1. A computer-implemented method for estimating Poisson impedance in a survey region having at least one well location and for generating a final model of elastic properties of the survey region, the method comprising: positioning seismic data recording sensors in the survey region at different locations including a wellbore; blasting at points of incidence in the survey region; retrieving, using the seismic data recording sensors, seismic data and transmitting the seismic data to a database; retrieving, by a computer, an upscaled well log data represented in time domain, including compressional velocity, shear velocity, and density data over various points of incidence in the survey region from the database; retrieving, by the computer, a set of angle image gathers, preconditioned to preserve a signal amplitude information of various angles of incidence for the survey region from the database; retrieving, by the computer, an earth model data that includes horizons information data for the survey region from the database; retrieving, by the computer, seismic velocity data over time for the survey region from the database; generating, by the computer, a set of low frequency models through interpolation based on the earth model data, and the upscaled well log data over various points of incidence for the survey region which are retrieved from the database; generating, by the computer, a set of angle dependent wavelets from a plurality of angles of incidence for the survey region which is retrieved from the database, by applying a well-to-seismic-tie calculation to the set of angle image gathers, and the upscaled well log data which is retrieved from the database; combining, by the computer, a first angle image gather from the set of angle image gathers, with the corresponding first low frequency model from the set of low frequency models, and to the set of angle wavelets; executing, by the computer, a nonlinear direct prestack seismic inversion model by employing the combined first angle image gather, with the first low frequency model, and the set of angle wavelets; repeating the combining using another angle image gather from the set of angle image gathers with the corresponding first low frequency model and the executing a nonlinear direct prestack seismic inversion model operations for the survey region until all gathers from the set of angle image gathers, and all corresponding low frequency models from the set of low frequency models have been combined to the set of angle wavelets; verifying that the repeating operation executed a final nonlinear direct prestack seismic inversion model by combining a final angle image gather from the set of angle image gathers, with a corresponding final low frequency model from the set of low frequency models to the set of angle wavelets; and generating a final model of elastic properties from the results of said nonlinear direct prestack seismic inversion model in the survey region by calculating the exponential integration of the elastic reflectiveness to characterize the survey region.

2. The computer-implemented method of claim 1, wherein the retrieving, by the computer, the set of angle image gathers includes: retrieving seismic survey data for the survey region; retrieving seismic velocity data for the survey region; executing a prestack time migration calculation using the retrieved seismic survey data, and the retrieved seismic velocity data for the survey region; generating offset image gathers; stacking the retrieved offset image gathers over time at various offsets of incidence; converting the offset image gathers to angle image gathers; and calculating the set of angle image gathers.

3. The computer-implemented method of claim 1, wherein the generating, by the computer, the set of low frequency models further includes: retrieving seismic velocity data for the survey region; retrieving interpreted horizons information from an earth-model data for the survey region; retrieving the upscaled well log data represented in time domain including compressional velocity, shear velocity, and density data over various points of incidence, for the survey region; smoothing the retrieved seismic velocity data, and the retrieved upscaled well log data; converting from the smoothed seismic velocity data, using the retrieved interpreted horizons information to an interval velocity data; interpolating the converted interval velocity data, with the interpreted horizons information, and the smoothed upscaled well log data; and generating the set of low frequency models.

4. The computer-implemented method of claim 1, wherein the generating, by the computer, the set of angle dependent wavelets from the plurality of angles of incidence over the survey region includes: retrieving angle image gathers for the survey region; retrieving upscaled well log data for the survey region; extracting an initial set of angle wavelets from the retrieved angle image gathers; computing reflectivities represented in time domain, from the retrieved upscaled well log data; executing a well-to-seismic calculation to the initial set of angle wavelets for the survey region; and generating a set of angle-dependent wavelets.

5. The computer-implemented method of claim 1, wherein the generating, by the computer, the final model of elastic properties includes: inverting the elastic reflectiveness of all nonlinear direct prestack seismic inversion models comprising the expression: d=Wf(m.sub.0)+e; generating volume models for compression velocity, Poisson impedance, and density; exporting the volume models to the database; and generating the final model of elastic properties.

6. The computer-implemented method of claim 1, wherein the set of low frequency models equal, in quantity, to the same amount of angle image gathers from the set of angle image gathers.

7. A digital computing system, comprising: a computer system output device; a system computer, coupled to the memory resource and to the computer system output device, comprising at least one memory resource comprising instructions and at least one hardware processor to execute the instructions within the at least one memory resource to implement: receiving seismic survey data from a plurality of seismic data recorders located over a plurality of offset distances from the plurality of points of incidence; executing the nonlinear direct prestack seismic inversion model of a seismic survey region associated with the plurality of points of incidence with at least one well location; inputting a user defined tolerance value; inputting an initial estimate model equal to a user input value, or zero; computing a set of initial data residuals using the initial estimate model; setting up the initial estimate model equal to an input estimate model; computing a set of model gradients; computing a set of data residual changes; updating the input estimate model using the set of model gradients, with the set of data residual changes; computing a set of data residuals using the updated input estimate model; generating a matrix norm from the computed set of data residuals; repeating the operations of computing a set of model gradients, computing a set of data residual changes, updating the input estimate model, and computing a set of data residuals until the generated norm equals the user defined tolerance value; and generating a final set of data residuals.

8. The digital computing system of claim 7, wherein the inputting the initial estimate model further comprises the expression m.sub.0.

9. The digital computing system of claim 7, wherein the computing flail the set of initial data residuals further comprises the expression r=d−Wf(m.sub.0).

10. The digital computing system of claim 7, wherein the setting the initial estimate model equal to the input estimate model further comprises the expression m.sub.0=m.sub.i.

11. The digital computing system of claim 7, wherein the set of model gradients further comprises the expression: $g\left(m_i\right)=\frac{\partial }{\partial m^T}J\left(m_i\right)=F^T{W}^T{C}_e^{-1}r+C_m^{-1}{m}_i+C_p^{-1}\left(m_i-\hat{m}\right)+C_{\ell }^{-1}ℒ\left(m_i\right).$

12. The digital computing system of claim 7, wherein the computing the set of residual changes further comprises the expression Δr=WFg(m.sub.i).

13. The digital computing system of claim 12, wherein the term F further comprises the expression: $F=\frac{\partial d}{\partial m}=\left[\begin{array}{ccc}{A}_1+2{\Psi }^T{D}_1{\Psi }_{\alpha }& B_1+2{\Psi }^T{D}_1{\Psi }_{\mathrm{PI}}& C_1+2{\Psi }^T{D}_1{\Psi }_{\rho }\\ A_2+2{\Psi }^T{D}_2{\Psi }_{\alpha }& B_2+2{\Psi }^T{D}_2{\Psi }_{\mathrm{PI}}& C_2+2{\Psi }^T{D}_2{\Psi }_{\rho }\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}\\ A_M+2{\Psi }^T{D}_M{\Psi }_{\alpha }& B_M+2{\Psi }^T{D}_M{\Psi }_{\mathrm{PI}}& C_M+2{\Psi }^T{D}_M{\Psi }_{\rho }\end{array}\right].$

14. The digital computing system of claim 7, wherein the updating of the input estimate model further comprises the expression m.sub.i+1=m.sub.i+Δm.

15. The digital computing system of claim 7, wherein the computing the set of data residuals further comprises the expression r=d−Wf(m.sub.i+1).

16. The digital computing system of claim 7 wherein the repeating, further increases the iterator i by one every time said operation is executed, until the generated norm is less than or equal to the user defined tolerance values.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The teachings of the present invention can be readily understood by considering the following detailed description in conjunction with the accompanying drawings.

(2) FIG. 1. Is a schematic diagram showing a top view of a survey region with the various points of incidence of seismic sources according to an embodiment of the present disclosure;

(3) FIG. 2. Is a schematic diagram showing a cross-sectional view of an illustrative environment with points of incidence of seismic sources, seismic receivers, a well location, a wellbore, the various transmission rays, and the various angles of incidence, according to certain embodiments of the present disclosure;

(4) FIG. 3. Is a flow chart showing the computer-implemented method employing nonlinear direct prestack seismic inversion for Poisson impedance, according to an embodiment of the present disclosure;

(5) FIG. 4. Is a flow diagram illustrating a method for generating a set of angle image gathers, according to an embodiment of the present disclosure;

(6) FIG. 5. Is a flow chart showing the operation of generating a set of low frequency models, according to an embodiment of the present disclosure;

(7) FIG. 6. Is a flow diagram illustrating an operation of generating a set of angle wavelets, according to an embodiment of the present disclosure;

(8) FIG. 7. Illustrates in graphic representation, the outputs generated by the steps of generating a set of angle image gathers (FIG. 7A), the steps of generating a set of low frequency models (FIG. 7B), the steps of generating a set of angle dependent wavelets (FIG. 7C), according to an embodiment of the present disclosure;

(9) FIG. 8. Illustrates in graphic representation, the test results of a well-to-seismic approach for comparing the synthetic to measured seismic data at well location, according to an embodiment of the present disclosure;

(10) FIG. 9. Illustrates in graphic representation, a comparison between the inverted Poisson impedance (PI) using the present computer-implemented method and system illustrated in FIG. 9A, versus the indirect Aki-Richards approximation, illustrated in FIG. 9B, according to an embodiment of the present disclosure;

(11) FIG. 10. Illustrates in graphic representation, a final model of elastic properties produced by the present computer-implemented method and system. FIG. 10A shows the P-wave velocity, FIG. 10B shows the Poisson impedance, and FIG. 10C shows the density, according to an embodiment of the present disclosure;

(12) FIG. 11. Is an electric diagram, in block form of the digital high performance computing system performing the computer-implemented method, according to an embodiment of the present disclosure; and

(13) FIG. 12. Is a flow chart showing the sub-routine of generating a nonlinear direct prestack seismic inversion model as performed by the system computer of the high performance digital computing system, according to an embodiment of the present disclosure.

DETAILED DESCRIPTION OF THE INVENTION

(14) Reference will now be made in detail, to several embodiments of the present disclosures, examples of which, are illustrated in the accompanying figures. It is noted that wherever practicable similar or like reference numbers may be used in the figures and may indicate similar or like functionality. The figures depict embodiments of the present disclosure for purposes of illustration only. One skilled in the art will readily recognize from the following description that alternative embodiments of the structures, systems, and methods illustrated herein may be employed without departing from the principles of the disclosure described herein.

(15) FIG. 1 illustrates a seismic survey region, 101, in which the preferred embodiment of the present invention is useful. It is important to note, that the survey region of FIG. 1 is a land-based region represented as 102. Persons of ordinary skill in the art, will recognize that seismic survey regions produce detailed images of local geology in order to determine the location and size of possible hydrocarbon (oil and gas) reservoirs, and therefore a well location 103. In these survey regions, sound waves bounce off underground rock formations during blasts at various points of incidence 104, and the waves that reflect back to the surface are captured by seismic data recording sensors, 105, transmitted by data transmission systems, 1101, wirelessly, 1102, from said sensors, 105, then stored for later processing, and analysis by the digital high performance computing system of FIG. 11.

(16) In FIG. 2 a cross-sectional view of a portion of the earth over the survey region is illustrated as 201, showing different types of earth formation, 102, 203, 204, which will comprise the Oseismic survey data in the present invention. In particular, persons having ordinary skill in the art will soon realize that the present example shows a common midpoint-style gather, where seismic data traces are sorted by surface geometry to approximate a single reflection point in the earth. In this example, data from several shots and receivers may be combined into a single image gather, or used individually depending upon the type of analysis to be performed. Although the present example may illustrate a flat reflector and a respective image gather class, other types or classes of image gathers known in the art maybe used, and its selection may depend upon the presence of various earth conditions or events.

(17) As shown on FIG. 2, seismic energy from multiple points of incidence 104, will be reflected from the interface between the different earth formations. These reflections will be captured by multiple seismic data recording sensors 105, each of which will be placed at different location offsets 210 from each other, and the well 103. Because all points of incidences 104, and all seismic data recording sensors 105 are placed at different offsets 210, the survey seismic data or traces, also known in the art as gathers, will be recorded at various angles of incidence represented by 208. The points of incidence 104 generate downward transmission rays 205, in the earth that are captured by their upward transmission reflection through the recording sensors 105. Well location 103, in this example, is illustrated with an existing drilled well attached to a wellbore, 209, along which multiple measurements are obtained using techniques known in the art. This wellbore 209, is used to obtain well log data, that includes P-wave velocity, S-wave velocity, Density, among others. Other sensors, not depicted in FIG. 2, are placed within the survey region to also capture horizons data information required for interpreters and persons of ordinary skilled in the art to perform various geophysical analysis. In the present example, the gathers will be sorted from field records in order to examine the dependence of amplitude, signal-to-noise, move-out, frequency content, phase, and other seismic attributes, on incidence angles 208, offset measurements 210, azimuth, and other geometric attributes that are important for data processing and imaging and known by persons having ordinary skills in the art.

(18) However, to perform elastic reflectivity analysis it is often useful in the art, to exploit the redundancy of seismic data and produce images with more dimensions than the two coordinates of the physical space. In the present invention, the final model of elastic properties 314 results in a three-dimensional cube made by computing a wide array of data, mainly represented in time domain. When the final model of elastic properties 314 is created, no additional dimensions get included which in other situations, could have either be the absolute offset and azimuth. In the preferred embodiment of this invention, because the final model of elastic properties 314 is created according to a three-dimensional coordinate, the cube will end up comprising the in-line, cross-line, and time input variables.

(19) With regards to FIG. 3, it illustrates a flow chart with an overview of the preferred embodiment of the invention. The retrieving phase 301 initiates the process, and comprises obtaining and preparing data and information from the survey region. In particular, four different types of inputs are retrieved from seismic data over a survey region: upscaled well log data represented in time domain 302, a set of angle image gathers 310, horizons information 304, and seismic velocity data 305. The well log data may be obtained directly from the well location and remotely transferred onto a database, 1103, for further compilation, such as that contemplated on FIG. 11. It may also be retrieved from database 1103, by system computer 1105. For the present invention, the well log data comprises the P-wave velocity, S-wave velocity, and density over the survey region. Further, the set of angle image gathers 310 (further exemplified by FIG. 7A) comprise a plurality of image gathers, constant with the amount of points of incidence 104 within a survey region 101 represented in the time domain at the various angles of incidence 208. Nonetheless, since this seismic survey data 402 data is too raw, noisy, or from various points of incidence 104, it needs to be further processed. This further refinement of the set of angle image gathers is sub-routine procedure 303, within the present computer-implement method is illustrated in FIG. 4. Procedure 303 retrieves 401 seismic survey data 402, and seismic velocity data 305 for survey region 101, all of which are to be stored in a database, such as storage device 1104, and retrieved therefrom by system computer 1105 for further processing. The sub-routine, using high performance computing of FIG. 11, performs a prestack time migration approach 403 which solves the problem of conflicting dips with different stacking velocities, and data. The offset image gathers are then generated 404, stacked 405, and converted back to the time domain 406 with a radial-trace mapping originally introduced for shot-geophone migration. The final set of angle image gathers are generated 310 for the various points of incidence 104 over the survey region 101. FIG. 7A illustrates an example of what comprises a set of angle image gathers in the present invention. Further, the high performance computing system of FIG. 11, contemplate these set of angle image gathers data 310 to be stored in a database, such as storage device 1104, and retrieved therefrom by system computer 1105.

(20) Continuing with the computer-implemented method of FIG. 3, in particular the retrieved horizons information data 304, may be remotely transferred onto a database, 1103, for further compilation, such as that contemplated on FIG. 11. Further, data 304 may also be retrieved from database 1103, by system computer 1105. The retrieved horizons information data 304, are typically represented for purposes of the present invention as a layer, a surface, an interface, in the earth, and its manifestation in the seismic data. In its purest form, 304 is unattainable because a person having skills in the art, can rarely know exactly where the horizon is in the subsurface. Yet, the horizons information data 304 comprises a matrix of discrete samples in a three dimensional plan (x,y,z) that can be stored in an ASCII file. As such, a horizon is something that can be unambiguously drawn on a map, and treated like a raster image. Some software tools even call attribute maps horizons, blurring the definition further; and therefore creating it as an artifact of software.

(21) Another set of data gathered by the present invention is seismic velocity, 305. Seismic velocity (both compressional and shear) 305 is a fundamental property of materials, a property that varies with changes in conditions both external (stress, temperature) and internal (fluid saturation, crack density). Monitoring of changes in these external or internal conditions is a goal of geophysical investigations such as earthquake prediction (via stress change monitoring) and reservoir exploitation (via fluid saturation monitoring).

(22) Once all data variables have been retrieved (302, 310, 304, and 305), the computer-implemented method proceeds to execute two sub-routine, using system computer 1105.

(23) The first sub-routine, using system computer 1105, involves generating a set of low frequency models 307 through interpolation based on the earth model data, and the upscaled well log data over various points of incidence 104 for the survey region 101. Sub-routine 307 then outputs, and stores on device 1104, a set of low frequency models 309, better exemplified by FIG. 7B. In tandem, sub-routine 306 is performed, which involves generating a set of angle dependent wavelets 308 (exemplified by FIG. 7C) from a plurality of angles of incidence 208 for the survey region 101, by computing a well-to-seismic-tie approach 605 to the set of angle image gathers 312, and the upscaled well log data 302.

(24) In connection with the first-subroutine of generating a set of low frequency models 307, FIG. 5 describes it further. Performing another retrieving operation, 500, of the upscaled well log data which is represented in time domain 302, in combination with the interpreted horizons information data 304, and the seismic velocity data 305, the present invention generates a set of low frequency models (further illustrated by FIG. 7B), using interpolation process, 506, based also on 302 and 305, over various points of incidence 204 for the survey region 101. However, due to the fact that the retrieved seismic data from geo-sensors 206, like those depicted in the present invention are bandlimited, the retrieved information lacks the low frequencies that are essential for a relative-to-absolute property conversion. These unfiltered or noisy set of low frequency models can be regarded as the mean trend of the data and the higher frequencies are the deviation from that mean trend. As such, the present invention filters or smooths-out 302 and 305, to generate a set low frequency models FIG. 7B that are a critical step for the present invention as they are used to account for missing information below the seismic data. Further, the present invention solves the problem of how to go about interpolating the logs between the different points of incidence 205.

(25) Because interpreters are usually uncertain about how to obtain a set of low frequency models like FIG. 7B, and often uncomfortable with the results, the present invention discloses a method that is compatible with the preferred embodiment of the present invention and produces better results. To achieve said better results, the present invention discloses the steps of generating a set of low frequency models that further comprise digitally smoothing 501 the retrieved seismic velocity data 305, and digitally smoothing 502 the upscaled well log data 302, by running the retrieved data over a set of low pass filters. This filtering computes an amplitude spectrum to meet specifications in the frequency domain, thereby cutting out high frequencies and significantly reducing the noise from retrieved upscaled log well data 302, as well as that of the retrieved seismic velocity 305. Once the seismic velocity 305 has been smoothed to 503, it is then converted to internal velocity data 504, using the retrieved interpreted horizons information 304 and following the Dix transform inversion method (Durbaum, 1954; and Dix, 1955). Said method relies primarily in computing interval velocities from picked root mean squares (rms), or stacking, velocities, and corresponding travel times by the well-known formula of:

(26) $\begin{array}{cc}=\sqrt{\frac{V_{\mathrm{rm}s,i}^2{t}_i-V_{\mathrm{rm}s,i-1}^2{t}_{i-1}}{\Delta t_i}}& \left(36\right)\end{array}$
where is the local internal interval velocity over the time interval Δt.sub.i=t.sub.i−t.sub.i−1, where i is an iterator starting from 1 to nth-value. Parameters V.sub.rms,i.sup.2 and V.sub.rms,i−1.sup.2 are rms velocities at the top and bottom interfaces of the interval. Instantaneous velocity V.sub.int.sup.i is assumed as a piecewise constant, V.sub.int.sup.i=, with discontinuities at the interfaces.

(27) The converted internal velocity 504, is then combined with the interpreted horizons information 304, and the smoothed upscaled well log data 505 to perform the step of interpolation 506. Step 506 is performed due to the fact that data, even after going through a smoothening process, the input data is still noisy and requires further refinement. As such, the present computer-implemented method uses an interpolation method known in the art as cokriging. The special characteristic of this method, and why it is of relevant use in the present invention, is because it can utilize data of different nature to model, optimize and interpolate particular attributes. In general terms, the computer-implemented method of the present invention applies an algorithm to an interpolated program coded to accept s set of variables (as oppose to just 2 sets found on typical cokriging algorithms) producing a linear system following the equation of:

(28) $\begin{array}{cc}\left(\begin{array}{cc}C& E\\ E^T& 0\end{array}\right)\left(\begin{array}{c}T\\ \mu \end{array}\right)=\left(\begin{array}{c}{C}_0\\ I_0\end{array}\right)& \left(37\right)\end{array}$
where C is the covariance (or its estimate) matrix of all know variables' pair, and C.sub.0 is the vector of pairwise covariances between the unknown variable and all other variables.
where, μ is the vector of all LaGrange multiples μ.sub.1 . . . μ.sub.s. E is a vector of matrices of I.sub.1 . . . I.sub.s. Each matrix I.sub.i, i∈{1 . . . s} is of size (number of points in i.sup.th variable set)×s. All elements in the i.sup.th column of I.sub.i, are one and all other entries are zero.
where, T is the vector of all coefficients, and I.sub.0 is a column vector of size s×1 of all elements under C.sub.0 on the right-hand side of the equation.

(29) Once the computer-implemented method has interpolated all variables for 304, 504, and 505 following the above algorithm, it generates a set of low frequency models 309 comprising of various individual low frequency models, of the same amount as there were angle image gathers.

(30) Similarly, sub-routine 306 involves the steps of generating a set of angle wavelets (illustrated on FIG. 7C), using system computer 1105, from a plurality of angles of incidence 208 over the survey region 101, and is represented by FIG. 6. The steps required to fully execute this subroutine further involve retrieving 601 the set of angle image gathers 310, and the upscaled well log data 302. An initial set of angle wavelets 602 is then extracted from the retrieved angle image gathers 312, which acts as an initial estimate of what a final angle dependent wavelet should look like. The extraction is computed within the present invention, using statistical wavelet extraction procedures which use the seismic traces alone to extract the angle dependent wavelets 604. However, in order to extract each angle wavelet 604, the corresponding trace in the set of angle image gathers 310 is used. Then an analysis window is extracted, and the start and end of the extracted analysis window are tapered in lengths equal to the lesser of 10 samples, or ¼ of the analysis window. An autocorrelation process is computed in a data window, wherein the length of the autocorrelation is equal to ½ of the desire wavelength. The amplitude spectrum of the autocorrelation is computed, and the square root of the autocorrelation is calculated. This operation is performed to approximate the amplitude spectrum of the wavelet. Thereafter, the desired phase (e.g. 0°, 5°, 10°, 15°, etc.) is added, and the inverse Fast Fourier Transform (FFT) is computed at 602 by 1105 to produce or extract the initial angle wavelet 604. The extraction procedure 602 verifies if there were other wavelets produced before, and generates the initial wavelet. If other wavelets were created, the new wavelet is generated next to the previously generated wavelet from other traces of the same angle image gathers 312. The extraction procedure continues until all angles from the angle image gathers 312 have been analyzed. Once all angle wavelets have been statistically extracted, an extracted set of initial angle wavelets is generated 604.

(31) According to the preferred embodiment of the invention, the next step in sub-routine 306 is computing angle dependent reflectivity 603 using the quadratic approximation of equation (38) that combines the angles of incidence 208 represented by θ.sub.i, with respect to the: (a) upscaled well log data 302 (b) P-wave velocity represented as α, (c) S-wave velocity represented as β (d) the Poisson impedance represented as PI, and (e) the density represented as ρ, over the survey region 101.

(32) $\begin{array}{cc}{R}_{\mathrm{PP}}\left({\theta }_i\right)=\left[\frac{1}{2}{\sec }^2{\theta }_i-\frac{4}{c}\left(\frac{\beta }{\alpha }\right){\sin }^2{\theta }_i\right]\frac{\mathrm{\Delta \alpha }}{\alpha }+\left[\frac{4}{c}\left(\frac{\beta }{\alpha }\right)-4{\left(\frac{\beta }{\alpha }\right)}^2\right]{\sin }^2{\theta }_i\frac{\Delta \mathrm{PI}}{\mathrm{PI}}+\left[\frac{1}{2}-\frac{4}{c}\left(\frac{\beta }{\alpha }\right){\sin }^2{\theta }_i+2{\left(\frac{\beta }{\alpha }\right)}^2{\sin }^2{\theta }_i\right]\frac{\mathrm{\Delta \rho }}{\rho }+{\left[\frac{2}{c}\frac{\Delta \alpha }{\alpha }-\left(\frac{2}{c}-2\frac{\beta }{\alpha }\right)\frac{\Delta \mathrm{PI}}{\mathrm{PI}}+\left(\frac{2}{c}-\frac{\beta }{\alpha }\right)\frac{\mathrm{\Delta \rho }}{\rho }\right]}^2\frac{\beta }{\alpha }\cos {\mathrm{\theta sin}}^2{\theta }_i& \left(38\right)\end{array}$

(33) Next, by performing a well-to-seismic approach 605 (graphically illustrated in FIG. 8) to the extracted initial set of angle wavelets 604 for the survey region, a set of angle-dependent wavelets 308 is generated, and stored on device 1104. Said set of angle-dependent wavelets is further exemplify in FIG. 7A.

(34) According to the preferred embodiment of the invention, once the set of angle image gathers 310, the set of angle wavelets 308, and the set of low frequency models 309 are generated, process 311 is next performed. This process combines the three different sets of data as a precursor to sub-routine 312, using system computer 1105 and storing them in device 1104.

(35) Sub-routine 312 (illustrated in FIG. 12) it is then initiated by a user inputting 1201 through a combination of user input devices like 1108 (keyboard), and 1109 (mouse) a user defined tolerance value 1202; that a person having ordinary skills in the art would understand it to be less than 1 (i.e. 10.sup.−6). The user also inputs an initial estimate model 1203 comprising of variable m.sub.0, which persons of ordinary skill in the art also regard it as an initial model prediction. Both inputs are then sent from a user computer system to be stored in database device 1104 for processing by system computer 1105 following step 1204. Process 1204 computes a set of initial data residuals using the initial estimate model 1203, and the user defined tolerance value 1202. In geophysics term, a data residual is made by the difference between a measured data and predicted data. During this process, the computed initial data residuals are calculated using the following expression, which contemplates that those skilled in the art will readily recognize it as a term derived from other known data residual expression for computing the data residuals between observed data and predicted data:
r=d−Wf(m.sub.0)  (39)

(36) Thereafter, the computer-implemented method of the present invention begins preparing the system computer 1105 to run a system of loops, in order to obtain the final non-linear direct prestack seismic inversion 313. However, a user must first set up the initial estimate model 1205 in the personal computer from value m.sub.0 to m.sub.i. This is interim computing step, messages the system computer 1105 that model m.sub.i, as defined in the system computer 1154, will get reiterated until a certain requirement has been met. According to the preferred embodiment of the present invention, i in expression m.sub.i is an iterator stored in device 1104, used as a pointer by system computer 1105, for abstracting addresses of elements of different structures, in other words multiple linear objects. As a result, i allows data structures to interface with algorithms of the present invention, that aren't aware of the structure they are operation on. Further, said iterator is used by system computer 1105 to generate a new iterator, based on the iterator that it is currently available, and upon passing it through the various required process (see 1206, 1207, 1208, 1209 and 1210) it returns the next iterator in the sequence in for the model variable m.

(37) As all the values have inputted into system computer 1105, as well as the initial data residuals 1204 already computed; a set of model gradients is then computed, 1206, using expression:

(38) $\begin{array}{cc}g\left(m_i\right)=\frac{\partial }{\partial m^T}J\left(m_i\right)=F^T{W}^T{C}_e^{-1}r+C_m^{-1}{m}_i+C_p^{-1}\left(m_i-\hat{m}\right)+C_{\ell }^{-1}ℒ\left(m_i\right)& \left(40\right)\end{array}$

(39) The digital computing system of FIG. 11, using system computer 1105 then computes a set of data residuals changes 1207 for the linear and non-linear reflectivities of the approximation to the Zoeppritz equation which further comprises the expression:
Δr=WFg(m.sub.i)  (41)

(40) According to the preferred embodiment of the invention, the process of updating the input estimate model 1208 is performed by system computer 1105 which calculates the changes from any previous model (i.e. m.sub.0) to the current model (m.sub.i), and adds it to the current model. An updated input model m.sub.i+1 is then generated and the iterator gets a value of one added.

(41) The step of computing the data residuals 1209 using expression r=d−Wf(m.sub.i+1) is performed next. With step 1209 completed, system computer 1105 generates a matrix norm 1210 for the computed data residuals 1209. System computer 1105 then performs a series of iterations 1211, using the norm computed data residuals 1209 with the inputted user defined tolerance value 1202, until it confirms that the computed data residual 1209 is less than the tolerance value 1202. Once system computer 1209 confirms such inequality, the iteration process 1211 is concluded and the inverted results, in the form of a final set of updated models is generated, 1212, and outputted onto storage device 1106, alongside the entire set of data generated during sub-routine 312. A user may then retrieve the data using a personal computer 1107 by selecting an export function.

(42) With sub-routine 312 concluded, the computer-implemented method of the present invention illustrated in FIG. 3 performs iteration 313, wherein it repeats the steps of combining the set of low frequency models 309 with the set of angle dependent wavelets 308, until all gathers from the set of angle image gathers 310, and all corresponding low frequency models from the set of low frequency models 309 have been combined to the set of angle wavelets 308. System computer 1105 then verifies that the repeating step performed a final nonlinear direct prestack seismic inversion model, in other words, there are no more final angle image gathers (from the set of angle image gathers), to combine with corresponding final low frequency models (from the set of low frequency models) to the set of angle wavelets. If the system computer 1105 verifies 313 that the entirety set of models have been combined, it produces a final model of elastic properties, 314, as shown on FIG. 10. If the system computer 1105 is not able to verify 313 that it is ready to produce a final model of elastic properties 314, it will initiate again the step of combining all set of models and data, as well as sub-routine 312, until it is able to produce the final model of elastic properties 313, which will then get outputted by system computer 1105 onto storage device 1106. The user of personal computer 1107 may then retrieved and display the final model of elastic properties 314 for further analysis or recommendations.

(43) For any of the foregoing embodiments, the computer system 1107 may include any of the following elements, alone or in combination with each other. The functions performed by the system computer 1105 include functions to input, retrieve, convert, combine, generate, produce, calibrate, perform, extract, filter, and output; an array of different values, including well log data, horizons information, seismic velocity, P-wave velocity, Poisson Impedance, density, reflectivity, angle image gathers, low frequency models, angle dependent wavelets, as well as more complex functions like performing the nonlinear direct prestack seismic inversion model of the present invention based upon information related, and obtained from, survey region containing at least one well location and wellbore penetrating the subterranean earth formation.

(44) Turning over to FIG. 7, it illustrates in graphical form, the different sets of data (310, 308, and 309) generated as an output of the aforementioned sub-routines (303, 306, and 307). Similarly, FIG. 8 illustrates in graphical form the process of performing a well-to-tie approach by comparing synthetic seismogram to survey seismic data using angle-dependent wavelets, and a computed angle-dependent reflectivity.

(45) In FIG. 9, the graphical representation shows a comparison between the results using the computer-implemented method of the present invention, employing nonlinear direct prestack seismic inversion for Poisson (PI) impedance, illustrated by FIG. 9A, against another computer-implemented method employing the indirect Aki-Richards approximation, illustrated by FIG. 9B. On both graphic representations, line 901 represents the true model, and line 902 represents the results obtained from each computer-implemented method. Clearly, graphical representation of FIG. 9A shows more accurate results than FIG. 9B as represented by gap created by the differences between line 901 and line 902.

(46) In FIG. 10 the graphical representation of the computer-implemented final model of elastic properties 315 is shown. When dividing up the parameters according to their individual representations, FIG. 10A shows the P-wave velocity, FIG. 10B shows the Poisson impedance, and FIG. 10C shows the density. The 1001 lines represent the true parameters, the 1002 lines represent the initial parameters, and the 1003 lines represent the parameters as produced in the final model of elastic properties 315 by the present invention. As expected, the parameters produced by the present invention, match well with the true ones in velocity, Poisson impedance, and density.

(47) As used herein, the term “computing” encompasses a wide variety of actions, including calculating, determining, processing, deriving, investigation, look ups (e.g. looking up in a table, a database or another data structure), ascertaining and the like. It may also include receiving (e.g. receiving information), accessing (e.g. accessing data in a memory) and the like. Also, “computing” may include resolving, selecting, choosing, establishing, and the like.

(48) According the preferred embodiment of the present invention, certain hardware, and software descriptions were detailed, merely as example embodiments and are not to limit the structure of implementation of the disclosed embodiments. For example, although many internal, and external components of the computer system of FIG. 12 have been described, those with ordinary skills in the art will appreciate that such components and their interconnection are well known. Additionally, certain aspects of the disclosed invention may be embodied in software that is executed using one or more processing units as 1105. Program aspects of the technology may be thought of as “products” or “articles of manufacture” typically in the form of executable code and/or associated data that is carried on, or embodied in, a type of machine readable medium. Tangible non-transitory “storage” type media and devices include any or all memory or other storage for the computers, process or the like, or associated modules thereof such as various semiconductor memories, tape drives, disk drives, optical or magnetic disks, and the like which may provide storage at any time for the software programming.

(49) Additionally, the flowcharts and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present disclosure. It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the Figures. For examples, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowcharts illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified hardware functions or acts, or combinations of special purpose hardware and computer instructions.

(50) While in the foregoing specification this disclosure has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. On the contrary, a wide variety of modifications and alternative embodiments will be apparent to a person skilled in the art, without departing from the true scope of the invention, as defined in the claims set forth below. Additionally, it should be appreciated that structural features or method steps shown or described in any one embodiment herein can be used in other embodiments as well.

(51) TABLE-US-00001 Symbols Table Symbol Brief Definition R.sub.PP P-wave reflection coefficient on a horizontally elastic interface R.sub.PP (p) P-wave reflection coefficient as a function of ray parameter R.sub.PP (θ) P-wave reflection coefficient as a function of angle of incidence R.sub.f Fluid-to-fluid reflection coefficient on a horizontally acoustic interface p Ray parameter q Vertical slowness as a function of velocity and ray parameter α P-wave velocity β S-wave velocity ρ Density μ Shear moduli PI Poisson impedance Δ α Contrast in P-wave velocity across a horizontally elastic interface Δ β Contrast in S-wave velocity across a horizontally elastic interface Δ ρ Contrast in density across a horizontally elastic interface Δμ Contrast in shear moduli across a horizontally elastic interface ΔPI Contrast in Poisson impedance across a horizontally elastic interface θ Angle, Angle of incidence between the reflecting ray and the vertical c Rotation optimization factor of Poisson impedance R.sub.PP Linear reflectivity of the approximation to the Zoeppritz equation {tilde over (R)}.sub.PP Nonlinear reflectivity of the approximation to the Zoeppritz equation A, B, C, D, Other variables used for formula derivation E, F, G and Ψ m.sub.α Model vector of elastic parameter reflectivity for P-wave velocity m.sub.PI Model vector of elastic parameter reflectivity for S-wave velocity m.sub.ρ Model vector of elastic parameter reflectivity for density m Model vector of three elastic parameter reflectivities m.sub.0 Initial model vector for three elastic parameter reflectivities {circumflex over (m)} Low frequency trend model vector from well-logging or stacking velocity Δm Model perturbation vector between actual model and initial model r Linear reflectivity vector for R.sub.PP {tilde over (r)} Nonlinear reflectivity vector for {tilde over (R)}.sub.PP d Observed seismic data vector e Vector containing noise W Given wavelet as a function of angle of incidence f Nonlinear forward modeling operator that predicts data from model F Fréchet derivatives indicating the sensitivity of data to model perturbation C Covariance matrix ∂d/∂m Matrix of partial derivatives of data vector with respect to model vector J (m) Objective function as a function of model vector g (m) Model gradient as a function of model vector r Data residual between observed data and predicted data A, B, C, D Other matrices used for formula derivation and Ψ