DEVICE AND METHOD FOR ESTIMATING PRE-STACK WAVELET MODEL FROM SEISMIC GATHERS

Abstract

Computing device, computer instructions and method for calculating a characteristic of the earth based on recorded seismic data related to a subsurface. The method includes receiving well log data; calculating a statistical model based on the well log data; calculating a reflectivity R of the subsurface based on an intercept A and a gradient B of the statistical model; calculating a reconstructed gather D based on the reflectivity R and a wavelet model W; calculating parameters w.sub.n of the wavelet model W with a first inversion function C that depends on (i) the reconstructed gather D, (ii) a recorded gather D.sub.0, and (iii) a probability density function (pdf) associated with the statistical model; and calculating the characteristic of the earth for the subsurface, based on the parameters w.sub.n of the wavelet model W and the recorded gather D.sub.0.

Claims

1. A method for calculating a characteristic of the earth based on recorded seismic data related to a subsurface, the method comprising: receiving well log data; calculating a statistical model based on the well log data; calculating a reflectivity R of the subsurface based on an intercept A and a gradient B of the statistical model; calculating a reconstructed gather D based on the reflectivity R and a wavelet model W; calculating parameters w.sub.n of the wavelet model W with a first inversion function C that depends on (i) the reconstructed gather D, (ii) a recorded gather D.sub.0, and (iii) a probability density function (pdf) associated with the statistical model; and calculating the characteristic of the earth for the subsurface, based on the parameters w.sub.n of the wavelet model W and the recorded gather D.sub.0.

2. The method of claim 1, wherein the characteristic of the earth is one of a density, P-impedance or S-impedance, the recorded gather D.sub.0 is obtained from the recorded seismic data and the recorded seismic data is acquired with seismic sensors.

3. The method of claim 1, wherein the step of calculating a reflectivity R comprises: calculating plural values of the intercept A and gradient B based on the well log data; and calculating an anticorrelation, a blue spectrum and sparseness of the plural values to obtain the statistical model.

4. The method of claim 1, wherein the step of calculating a reconstructed gather D comprises: selecting a structure of the wavelet model that has a scale term a, a phase-only term W.sub.P(z) and a zero-phase term W.sub.A(z).

5. The method of claim 4, wherein the phase-only term W.sub.P(z) and the zero-phase term W.sub.A(z) are each an all-pass autoregressive moving average structure.

6. The method of claim 1, wherein the first inversion function is a Bayesian function.

7. The method of claim 1, wherein the first inversion function has (i) a first term C.sub.P that depends with a probability density function (pdf) of the statistical model and (ii) a second term C.sub.D that depends on a difference between the reconstructed gather D and the recorded gather D.sub.0.

8. The method of claim 7, wherein the first term C.sub.P is corrected with a Jacobian term C.sub.J to define a scaling of the first inversion function.

9. The method of claim 7, wherein the second term C.sub.D is modified with a noise-whitening operator for a noise that is non-white.

10. The method of claim 1, wherein the wavelet model W is continuous in a time and an angle that characterize the recorded gather.

11. The method of claim 1, wherein the recorded gather is a common image gather.

12. The method of claim 1, further comprising: using a second inversion function and the wavelet model W to process the seismic data to generate an image of the subsurface.

13. A computing device for calculating a characteristic of the earth based on recorded seismic data related to a subsurface, the computing device comprising: an input/output interface for receiving well log data; and a processor connected to the input/output interface and configured to, calculate a statistical model based on the well log data; calculate a reflectivity R of the subsurface based on an intercept A and a gradient B of the statistical model; calculate a reconstructed gather D based on the reflectivity R and a wavelet model W; calculate parameters w.sub.n of the wavelet model W with a first inversion function C that depends on (i) the reconstructed gather D, (ii) a recorded gather D.sub.0, and (iii) a probability density function (pdf) associated with the statistical model; and calculate the characteristic of the earth for the subsurface, based on the parameters w.sub.n of the wavelet model W and the recorded gather D.sub.0.

14. The computing device of claim 13, wherein the characteristic of the earth is one of a density, P-impedance or S-impedance, the recorded gather D.sub.0 is obtained from the recorded seismic data and the recorded seismic data is acquired with seismic sensors.

15. The computing device of claim 13, wherein the processor is further configured to, calculate plural values of the intercept A and gradient B based on the well log data; and calculate an anticorrelation, a blue spectrum and sparseness of the plural values to obtain the statistical model.

16. The computing device of claim 13, wherein the processor is further configured to, select a structure of the wavelet model that has a scale term , a phase-only term W.sub.P(z) and a zero-phase term W.sub.A(z).

17. The computing device of claim 16, wherein the phase-only term W.sub.P(z) and the zero-phase term W.sub.A(z) are each an all-pass autoregressive moving average structure.

18. The computing device of claim 13, wherein the first inversion function is a Bayesian function, and the first inversion function has (i) a first term C.sub.P that depends with a probability density function (pdf) of the statistical model and (ii) a second term C.sub.D that depends on a difference between the reconstructed gather D and the recorded gather D.sub.0.

19. The computing device of claim 13, wherein the processor is further configured to use a second inversion function and the wavelet model W to process the seismic data to generate an image of the subsurface.

20. A non-transitory computer readable medium, including computer executable instructions, wherein the instructions, when executed by a processor, implement a method for calculating a characteristic of the earth based on recorded seismic data related to a subsurface, the method comprising: receiving well log data; calculating a statistical model based on the well log data; calculating a reflectivity R of the subsurface based on an intercept A and a gradient B of the statistical model; calculating a reconstructed gather D based on the reflectivity R and a wavelet model W; calculating parameters w.sub.n of the wavelet model W with a first inversion function C that depends on (i) the reconstructed gather D, (ii) a recorded gather D.sub.0, and (iii) a probability density function (pdf) associated with the statistical model; and calculating the characteristic of the earth for the subsurface, based on the parameters w.sub.n of the wavelet model W and the recorded gather D.sub.0.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0018] The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings:

[0019] FIG. 1 is a schematic diagram of a marine seismic data acquisition system;

[0020] FIG. 2 is a schematic diagram of a land seismic data acquisition system;

[0021] FIG. 3 illustrates a seismic wave's path from a source to a receiver underground;

[0022] FIG. 4 illustrates a migrated angle gather;

[0023] FIG. 5 illustrates a well reflectivity;

[0024] FIG. 6 is a flowchart of a method for calculating a characteristic of the earth using a wavelet model;

[0025] FIG. 7 illustrates the intercept and gradient values for a given well;

[0026] FIG. 8 illustrates a coloring wavelet;

[0027] FIG. 9 illustrates the coloring wavelet's spectrum;

[0028] FIG. 10 illustrates a probability density function for a statistical model of the well;

[0029] FIG. 11 is a flowchart of a method for calculating a wavelet model;

[0030] FIG. 12 illustrates a reconstruction error of an image gather calculated based on the wavelet model;

[0031] FIG. 13 illustrates an estimated reflectivity calculated based on the reconstructed image gather;

[0032] FIG. 14 illustrates wavelets calculated using a blind separation method;

[0033] FIG. 15 illustrates wavelets calculated using a non-blind separation method;

[0034] FIG. 16 is a flowchart of a method for calculating a characteristic of the earth based on the wavelet model; and

[0035] FIG. 17 is a schematic diagram of a computing device configured to implement the methods discussed above.

DETAILED DESCRIPTION

[0036] The following description of the embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are discussed, for simplicity, with regard to common image gathers (CIG). However, the embodiments to be discussed next are not limited to CIG, but may be also applied to other type of gathers.

[0037] Reference throughout the specification to one embodiment or an embodiment means that a particular feature, structure, or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases in one embodiment or in an embodiment in various places throughout the specification is not necessarily referring to the same embodiment. Further, the particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0038] As previously discussed, the seismic elastic inversion is a process that takes as input a pre-stack processed seismic dataset. A pre-stack seismic dataset is a set for which the traces were not previously stacked. This dataset may take the form of CIG gathers computed on a dense set of surface locations (x,y) by a process called migration. Although in the following embodiments only CIG gathers are used to exemplify the novel method, other type of gathers may be used, as the common midpoint gathers, or the common shot gathers, or the common receiver gathers, etc. The CIG gather is a subset of the entire data seismic set collected during the acquisition phase. The traces present in the CIG gather are selected to take into account the dipping reflector geometry of the subsurface.

[0039] For a given surface location (x,y), as illustrated in FIG. 3, the measured CIG is denoted D.sub.0(t,), where t is the migrated time and the local reflection angle. FIG. 3 shows a wavenumber vector k.sub.s for the wavefields from the source S, a wavenumber vector k.sub.r for the wavefields from the receiver R, and a migration vector k. Gather D.sub.0 is shown in FIG. 4. The migrated time t can be replaced by the migrated depth z and/or the angle may be replaced by the offset h (i.e., a distance between the source S and the receiver R in the plane that includes the source S and receiver R) without changing the nature of the described method.

[0040] During the migration process, a velocity model is used. The velocity model is correct if the gathers exhibit only horizontal events. When this is not the case, the gathers are not totally horizontal, as shown in FIG. 4 (note the up curving of the black and white lines at large amplitudes).

[0041] The seismic elastic inversion process takes CIGs as input and outputs, for every surface location (x,y), quantities related to the local elastic properties of the earth, as for example, the P-impedance I.sub.P(t), the S-impedance I.sub.S(t) and, optionally, the density (t).

[0042] However, in order to do so, the seismic elastic inversion also needs a pre-stack wavelet model W(t,). More precisely, the seismic elastic inversion considers that each gather at a given surface location can be modeled as described by the following equation:

D.sub.0(t,)=W(t,)*R(t,)+N(t,)(1)

where D.sub.0(t,) is the input pre-stack gather, W(t,) is the wavelet model, R(t,) is the reflectivity, and N(t,) is the noise. The symbol * denotes convolution in time. The inputs to the seismic elastic inversion are the gather and the wavelet model. The output of the elastic inversion are the elastic parameters such as P-impedance I.sub.P(t), the S-impedance I.sub.S(t) and the density (t).

[0043] The seismic elastic inversion process estimates/calculates, according to this embodiment, from the elastic parameters (t), V.sub.P(t), and V.sub.S(t) (or alternatively, uses the Zoeppritz equation), an intercept A(t), a gradient B(t) and optionally a curvature C(t) based on equation (2) below, and then the process estimates a reflectivity R(t,) by using the two-term or three term Aki-Richards equation (3) (or other equivalent methods) as follows:

[00001]$\begin{array}{cc}A=\frac{1}{2}.Math.\left(\frac{.Math..Math.V_P}{V_P}+\frac{}{}\right),B=\frac{.Math..Math.V_P}{2.Math..Math.V_P}-4.Math.{\left(\frac{V_S}{V_P}\right)}^2.Math.\frac{.Math..Math.V_S}{V_S}-2.Math.{\left(\frac{V_S}{V_P}\right)}^2.Math.\frac{}{},.Math.C=\frac{.Math..Math.V_P}{2.Math..Math.V_P},& \left(2\right)\\ R\left(t,\right)=A\left(t\right)+B\left(t\right).Math.{\sin }^2.Math.+C\left(t\right).Math.{\sin }^2.Math..Math..Math.{\tan }^2.Math..& \left(3\right)\end{array}$

[0044] A well reflectivity calculated based on equations (2) and (3) is shown in FIG. 5. Then, the modeled gather D(t,)=W(t,)*R(t,) is computed by using the input wavelet model W(t,) and the seismic inversion process minimizes the misfit between the modeled gather D(t,) and the measured gather D.sub.0(t,). However, the results of this process can be improved as discussed below.

[0045] According to an embodiment, the pre-stack wavelet model W(t,) is estimated with a continuous variation in time and angle, with an amplitude and a phase term that uses the available wells information only in a statistical sense. In other words, the novel pre-stack wavelet model is based on (1) estimating a statistical model for the intercept A and gradient B (and optionally the curvature C discussed with regard to equation (2)) from the well logs, and (2) using this statistical model in an inversion (e.g., Bayesian inversion) in which the cost function has (i) a first term representing the conformity of the estimated intercept and gradient (and optionally the curvature) with the statistical model and (ii) a second term representing the conformity of the gather computed from the estimated wavelet model, with the recorded pre-stack seismic gathers.

[0046] A statistical model for intercept A and gradient B has three statistical features. First, there is an anti-correlation between A and B which is equivalent to stating that usually B and A are of opposite sign, so that the amplitude decreases with angle, while keeping the possibility of an unusual case where A and B have the same sign and the amplitude increases versus angle. Second, the reflectivities do not have a white spectrum, but rather a blue spectrum. Third, the reflectivities do not have a Gaussian distribution, but rather a long-tail distribution. Well data, when available, can be used to extract these statistical properties in a quantitative way.

[0047] In one application, the wavelet model uses a structure which is autoregressive moving-average (ARMA) parametrized in lattice parameters.

[0048] In the following, the pre-stack wavelet model is discussed based only on the intercept A term and the gradient B term of the Aki-Richards equation (i.e., the two terms Shuey equation). Those skilled in the art would understand that the same concepts apply to the full three-term equation, by incorporating the curvature C term.

[0049] The method for constructing the pre-stack wavelet model is now discussed with regard to FIG. 6. In step 600, data from one or more well logs is received. This data may include, but is not limited to, density (t), P-wave velocity V.sub.P(t), and S-wave velocity V.sub.S(t). This data is received either directly from real measurements made at one or more wells or from processed well logs. In case that no well is available in the area of interest, a generic statistical model may be used to simulate these measurements. For example, the generic statistical model may use a zero-phase wavelet that can be statistically estimated for a given time window and a given partial angle stack (stacking of traces based on a common angle) by assuming a white reflectivity. Thus, in the following, well log recorded data is understood to mean either physically recorded data or a generic statistical model.

[0050] In step 602, a statistical model for the intercept A(t) and gradient B(t) is calculated. The statistical model is based on equation (2), where the third term is ignored, i.e.,

[00002]$\begin{array}{cc}A=\frac{1}{2}.Math.\left(\frac{.Math..Math.V_P}{V_P}+\frac{}{}\right),B=\frac{.Math..Math.V_P}{2.Math..Math.V_P}-4.Math.{\left(\frac{V_S}{V_P}\right)}^2.Math.\frac{.Math..Math.V_S}{V_S}-2.Math.{\left(\frac{V_S}{V_P}\right)}^2.Math.\frac{}{}.& \left(4\right)\end{array}$

[0051] It is noted that each quantity in equation (4) depends on time, but for simplicity, this has been omitted. This means that as the log well data is measured over time, the density and various velocities change in time. These changes in time are reflected in the intercept A and gradient B values. When the values of the intercept A and gradient B are calculated based on equation (4) and then plotted on a graph, as in FIG. 7 (the intercept is plotted on the X axis and the gradient is plotted on the Y axis), a cross-plot of A(t) and B(t) is obtained for a given well. FIG. 7 shows that the dominant diagonal corresponds to the anticorrelation of these two quantities. The information shown in FIG. 7 illustrates the statistical model.

[0052] Next, in step 604, a blind source separation of the intercept A(t) and gradient B(t) of the statistical model is calculated. Prior to actually applying the blind source separation to the intercept and gradient of the statistical model, a few general considerations about source separation are now presented. Note that source separation problems in digital signal processing are those in which several signals have been mixed together into a combined signal and the objective is to recover the original component signals from the combined signal. Source separation is a process that estimates N signals s.sub.i(t) from N input independent components e.sub.j(t) and a correlation matrix c.sub.ij, such that these quantities respect the following equation:

[00003]$\begin{array}{cc}{s}_i\left(t\right)=\underset{j=1}{\overset{N}{.Math.}}.Math..Math.c_{\mathrm{ij}}.Math.e_j\left(t\right),\mathrm{with}.Math..Math.i=1,\mathrm{.Math.}.Math.,N.& \left(5\right)\end{array}$

[0053] In blind source separation, the input independent components e.sub.j(t) are unknown, but they are supposed to be independent with a known probability density function (pdf) p(x). If the pdf is Gaussian, then only a matrix V=C.sup.TC can be recovered, where matrix C has components c.sub.ij. However, if the pdf is non-Gaussian, then the matrix C can be estimated, and the independent components e.sub.j(t) may be computed by applying the inverse matrix C.sup.1 to the signals s.sub.i(t). Matrix C is selected such that after applying the matrix C.sup.1 to the signals s.sub.i(t), the components e.sub.j(t) respect the following relation:

[00004]$\begin{array}{cc}\frac{1}{N_t}.Math.\underset{t}{.Math.}.Math..Math.e_i\left(t\right).Math.f^{}\left[e_j\left(t\right)\right]={}_{\mathrm{ij}},& \left(6\right)\end{array}$

where N.sub.t is the number of t values, function f is defined as f(x)=log p(x), the prime sign next to function f indicates a space derivative, and .sub.ij=1 if i=j, and 0 otherwise.

[0054] Once the independent components e.sub.j(t) are calculated, they can be whitened by computing coloring wavelets. FIG. 8 shows a coloring wavelet 800 obtained by using a non-gaussian blind deconvolution algorithm, and FIG. 9 shows the amplitude spectrum 900 of the coloring wavelet. It can be seen that the amplitude spectrum of the wavelet increases with frequency (blue spectrum). Combining these two steps results in a convolutional source separation as expressed by the following equation:

[00005]$\begin{array}{cc}{s}_i\left(t\right)=\underset{j=1}{\overset{N}{.Math.}}.Math..Math.c_{\mathrm{ij}}*e_j\left(t\right),\mathrm{with}.Math..Math.i=1,\mathrm{.Math.}.Math.,N,& \left(7\right)\end{array}$

where c.sub.ij are the wavelets which are convolved with the independent whitened components e.sub.j(t).

[0055] Returning to step 604, the blind source separation process discussed above is applied to the statistical model that includes the intercept A and the gradient B. In this case, N=2, s.sub.1(t)=A(t), s.sub.2(t)=B(t) and the statistical model includes the wavelets c.sub.ij(t) obtained by applying equation (7), i.e., the convolutional source separation. The statistical model also includes the probability density function p(x) of the independent whitened components e.sub.j(t). Different pdfs can be used and the function p(x) can be adapted to the statistics of each well, once the independent whitened components e.sub.j(t) have been computed. A pdf corresponding to a logistic distribution is given by:

[00006]$\begin{array}{cc}p\left(x\right)=\frac{1}{{\left(2.Math..Math.\cosh .Math.\frac{x}{2}\right)}^2}.& \left(8\right)\end{array}$

The function of equation (8) is a good candidate for the pdf when the number of samples of the well is not sufficient to estimate precisely the exact pdf of the components. This pdf is illustrated in FIG. 10 by curve 1000. Curve 1002 is a Gaussian distribution. It can be seen that the logistic distribution 1000 has a longer tail than the Gaussian distribution 1002, which corresponds to the observation that reflectivities are sparser than Gaussian distributions (they have more events that are much larger than the background).

[0056] Extracting the statistical model from the well logs has allowed to precisely quantify three known features of the intercept A(t) and gradient B(t) components: their anticorrelation, their blue spectrum and their sparseness.

[0057] Returning to the method of FIG. 6, after performing the separation of the statistical model in step 604, the method advances to step 606 for calculating the wavelet model W(t,). The structure of the wavelet model is best explained in the z-transform notation. The wavelet W(f) in the spectral domain is transformed to the z-transform W(z) by using z=exp(2jft), where t is the time sampling interval, j is the complex number given by square root of minus one, and f is the frequency. The model used for the wavelet model W(z) is a product of (1) a scale term , (2) a normalized phase-only term W.sub.P(z) and (3) a normalized zero-phase term W.sub.A(z). Thus, the wavelet model can be written as:

W(z)=W.sub.P(z)W.sub.A(z).(9)

[0058] To construct the phase-only term W.sub.P(z), it is possible to use an all-pass autoregressive moving average (ARMA) structure as described by the following equations:

[00007]$\begin{array}{cc}S\left(z\right)=\frac{z^P.Math.G\left(z^{-1}\right)}{G\left(z\right)},\mathrm{with}.Math..Math.G\left(z\right)=\underset{n=0}{\overset{P}{.Math.}}.Math..Math.g_n.Math.z^n.Math..Math.\mathrm{and}.Math..Math.g_0=1,& \left(10\right)\end{array}$

where G(z) is a normalized minimum phase filter of length P.

[0059] The normalized phase-only term W.sub.P(z) is then defined as the product of a casual all-pass wavelet and an anticausal all-pass wavelet as described by equation:

[00008]$\begin{array}{cc}{W}_P\left(z\right)=\frac{z^P.Math.G\left(z^{-1}\right)}{G\left(z\right)}.Math.\frac{z^{-Q}.Math.H\left(z\right)}{H\left(z^{-1}\right)},\mathrm{with}.Math.\text{.Math.}.Math.G\left(z\right)=\underset{n=0}{\overset{P}{.Math.}}.Math..Math.g_n.Math.z^n,H\left(z\right)=\underset{n=0}{\overset{Q}{.Math.}}.Math..Math.h_n.Math.z^n,g_0=1.Math..Math.\mathrm{and}.Math..Math.h_0=1,& \left(11\right)\end{array}$

where H(z) is a normalized minimum phase filters of length Q.

[0060] The zero-phase term W.sub.A(z) is an ARMA structure with a zero-phase numerator and a zero-phase denominator as described by the following equation:

[00009]$\begin{array}{cc}{W}_A\left(z\right)=\frac{U\left(z\right).Math.U\left(z^{-1}\right)}{V\left(z\right).Math.V\left(z^{-1}\right)},\mathrm{with}.Math..Math.U\left(z\right)=\underset{n=0}{\overset{R}{.Math.}}.Math..Math.u_n.Math.z^n,V\left(z\right)=\underset{n=0}{\overset{s}{.Math.}}.Math..Math.v_n.Math.z^n,u_0=1.Math..Math.\mathrm{and}.Math..Math.v_0=1,& \left(12\right)\end{array}$

where U(z) and V(z) are normalized minimum phase filters of lengths R and S, respectively.

[0061] The four normalized minimum phase filters G(z), H(z), U(z) and V(z) can be parametrized by their lattice parameters k.sub.n. With this lattice structure, the filters are minimum phase if all their lattice parameters are between 1 and 1. In order to make the wavelet model smooth in terms of time and angle variations while having a value between 1 and 1, the lattice parameters k.sub.n are selected as follows:

[00010]$\begin{array}{cc}{k}_n\left(t,\right)=\sin \left[\underset{\mathrm{pq}}{\overset{}{.Math.}}.Math.w_{\mathrm{npq}}.Math.T_p\left(t\right).Math.T_q\left(\right)\right],& \left(13\right)\end{array}$

where T.sub.p(t) and T.sub.q() can be taken as orthogonal polynomials.

[0062] The scale term of the wavelet W(z) can be taken as follows:

[00011]$\begin{array}{cc}\left(t,\right)=\exp \left[-\underset{\mathrm{pq}}{\overset{}{.Math.}}.Math.w_{\mathrm{pq}}.Math.T_p\left(t\right).Math.T_q\left(\right)\right].& \left(14\right)\end{array}$

[0063] The coefficients W.sub.npq for the filters G, H, U, and V and the coefficients w.sub.pg for the scaling are denoted w.sub.n and they describe the wavelet model W(t,).

[0064] The method now advances to step 608 for applying an inversion to the wavelet model W(t,). This inversion is different than the seismic elastic inversion previously discussed. The inversion may include an a priori term that describes the conformity of the estimated intercept and gradient with the statistical model and a posteriori term that describes the conformity of the gather computed from the estimated wavelet model with the measured pre-stack seismic data. Thus, the a priori term describes parameters of the earth acquired independent of the seismic data from which the gathers are determined, and the a posteriori term describes the acquired seismic data. In other words, the a priori term relates to events before the seismic data acquisition takes place and the a posteriori term relates to events of the seismic survey.

[0065] In the embodiment discussed now, a Bayesian inversion having the a priori term and the a posteriori term is applied to the wavelet model. Again, this inversion is not to be mixed up with an inversion that is applied later to the CIG gathers for determining various characteristics of the earth or an image of the surveyed subsurface.

[0066] The unknowns for the Bayesian inversion are the independent whitened components e.sub.1(t) and e.sub.2(t) as well as the wavelet parameters w.sub.n. The independent whitened components cannot be computed from the well logs because it is assumed that the desired calculations are performed at a location (x,y) where no well is available. However, it is still possible to consider that the statistical model computed from a nearby well is still valid.

[0067] The Bayesian inversion process, which is illustrated in FIG. 11, selects in step 1100 a cost function C(e.sub.1,e.sub.2,w.sub.n) (Bayesian cost function) that needs to be minimized or maximized for determining the unknown quantities e.sub.1, e.sub.2, and w.sub.n. Before doing this, the method computes, in step 1102, based on the independent components e.sub.1 and e.sub.2, the intercept and gradient of the statistical model of the well, which are given by the following equations:

A(t)=c.sub.11(t)*e.sub.1(t)+c.sub.12(t)*e.sub.2(t)

B(t)=c.sub.21(t)*e.sub.1(t)+c.sub.22(t)*e.sub.2(t).(15)

[0068] Then, in step 1104, the method calculates the angle dependent reflectivity using Shuey equation:

R(t,)=A(t)+B(t) sin.sup.2.(16)

[0069] In step 1106, a reconstructed gather is calculated by applying the time and angle variant lattice filter W to the reflectivity R as illustrated in the following equation:

D(t,)=W(t,)*R(t,).(17)

[0070] The Bayesian cost function C selected in step 1100 has one a priori term C.sub.P, which is the negative logarithm of the a priori probability density function p(x) of the independent components e.sub.1 and e.sub.2. This term can be written, after normalization, as follows:

[00012]$\begin{array}{cc}{C}_P=-\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[\log .Math..Math.p\left[e_1\left(t\right)\right]+\log .Math..Math.p\left[e_2\left(t\right)\right]\right],& \left(18\right)\end{array}$

where Nt is the number of t values.

[0071] The Bayesian cost function C also has an a posteriori term C.sub.D that depends on the measured seismic data D.sub.0(t,), which is the negative logarithm of the probability density function of the noise N(t,). The a posteriori term can be written, if the noise is white and after normalization, as follows:

[00013]$\begin{array}{cc}{C}_D=\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.{\left[D\left(t,\right)-D_0\left(t,\right)\right]}^2,& \left(19\right)\end{array}$

where D.sub.0(t,) is the angle gather obtained from processing the recorded seismic data and D(t,) is the gather reconstructed in step 1106 from (1) components e.sub.1 and e.sub.2 and (2) the wavelet model W(t,) reconstructed from the wavelet parameters w.sub.n based on equations (11) to (14).

[0072] Thus, the Bayesian cost function C can be written as:

[00014]$\begin{array}{cc}C\left[e_1\left(t\right),e_2\left(t\right),w_n\right]=.Math..Math.C_P+C_D=.Math.\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[f\left[e_1\left(t\right)\right]+f\left[e_2\left(t\right)\right]\right]+\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.{\left[D\left(t,\right)-D_0\left(t,\right)\right]}^2,& \left(20\right)\end{array}$

where is the hyperparameter governing the compromise between the two terms of the cost function and may account for the minus sign in equation (18).

[0073] Having the cost function, the method minimizes or maximizes it in step 1108, by using known processes, for example, the conjugate gradient method. However, there is a scale ambiguity associated with this cost function. If all scale terms (t,) are divided by a factor .sub.0 and the independent components e.sub.1(t) and e.sub.2(t) are multiplied by .sub.0, then the angle gather D(t,) does not change, which means that the a posteriori term C.sub.D of the cost function also does not change. Contrary to this, the a priori term C.sub.P decreases. Thus, the cost function needs to be modified to prevent the scale term (t,) from going to infinity.

[0074] One approach, which is applied in step 1110, is to enforce a normalization of the scale term (t,) during the minimization of the cost function C. For example, the geometrical mean of the scale term (t,) can be forced to unity (or a constant value) as follows:

[00015]$\begin{array}{cc}\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.\log .Math..Math.{}^N\left(t,\right)=0& \left(21\right)\end{array}$

and then compute, after the minimization, a scaling quantity .sub.0 of the normalized values .sup.N(t,), e.sub.1.sup.N(t) and e.sub.2.sup.N(t) to obtain the final values (t,), e.sub.1(t) and e.sub.2(t) by using the following equation:

[00016]$\begin{array}{cc}\left(t,\right)=\frac{{}^N\left(t,\right)}{{}_0},e_1\left(t\right)={}_0.Math.e_1^N\left(t\right),\mathrm{and}.Math..Math.e_2\left(t\right)={}_0.Math.e_2^N\left(t\right).& \left(22\right)\end{array}$

[0075] The scaling quantity Go can be computed by using equation (6) for i=j:

[00017]$\begin{array}{cc}\frac{1}{N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[{}_0.Math.e_1^N\left(t\right).Math.f^{}\left[{}_0.Math.e_1^N\left(t\right)\right]+{}_0.Math.e_2^N\left(t\right).Math.f^{}\left[{}_0.Math.e_2^N\left(t\right)\right]\right]=1,& \left(23\right)\end{array}$

which is equivalent to

[00018]$\begin{array}{cc}\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[e_1^N\left(t\right).Math.f^{}\left[{}_0.Math.e_1^N\left(t\right)\right]+e_2^N\left(t\right).Math.f^{}\left[{}_0.Math.e_2^N\left(t\right)\right]\right]-\frac{1}{{}_0}=0,& \left(24\right)\end{array}$

which is equivalent to minimizing in .sub.0 the function:

[00019]$\begin{array}{cc}c\left({}_0\right)=\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[f\left[{}_0.Math.e_1^N\left(t\right)\right]+f\left[{}_0.Math.e_2^N\left(t\right)\right]\right]-\log \left({}_0\right)=0.& \left(25\right)\end{array}$

[0076] By replacing the independent components with their normalized values, i.e.,

[00020]$\begin{array}{cc}.Math.e_1\left(t\right)={}_0.Math.e_1^N\left(t\right).Math..Math.\mathrm{and}.Math..Math.e_2\left(t\right)={}_0.Math.e_2^N\left(t\right)& \left(26\right)\\ .Math.\mathrm{and}& \\ \frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.\log .Math..Math.\left(t,\right)=\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.\log .Math..Math.{}^N.Math.\left(t,\right)-\log .Math..Math.{}_0=-\log .Math..Math.{}_0,& \left(27\right)\end{array}$

it is possible to solve the scaling ambiguity of the cost function during the minimization by adding a Jacobian term C.sub.J to the a priori term C.sub.P as follows:

[00021]$\begin{array}{cc}{C}_P+C_J=\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[f\left[e_1\left(t\right)\right]+f\left[e_2\left(t\right)\right]\right]+\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.\log .Math..Math.\left(t,\right).& \left(28\right)\end{array}$

[0077] Thus, as a result of minimizing the cost function in step 1108 with the normalization noted in equation (28), which is applied in step 1110, the method generates in step 1112 the wavelet model, which is continuous in time and angle. FIG. 12 shows the reconstruction error D(t,)-D.sub.0(t,) for the CIG D.sub.0(t,) shown in FIG. 4 and FIG. 13 shows the estimated well reflectivity, which can be compared to the actual well reflectivity shown in FIG. 5. FIG. 14 shows the estimated blind wavelet D(t,) (for the middle of the time window shown in the figure) and FIG. 15 shows the deterministic wavelet D(t,) computed in non-blind mode (using the actual reflectivity of FIG. 5). One skilled in the art would note that the blind wavelets are very close to the non-blind wavelets.

[0078] If the noise N (t,) of the angle gathers is non-white, it is possible to use a noise-whitening operator A.sub.N(t) for modifying the data dependent (a posteriori) term C.sub.D, as follows:

[00022]$\begin{array}{cc}{C}_D=\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.{\left[A_N\left(t\right)*\left(D\left(t,\right)-D_0\left(t,\right)\right)\right]}^2,& \left(29\right)\end{array}$

[0079] The noise-whitening operator A.sub.N(t) has a power spectral density of the noise NN(f) given by:

[00023]$\begin{array}{cc}\mathrm{NN}\left(f\right)=\frac{{}_N^2}{{.Math.A_n\left(f\right).Math.}^2}.& \left(30\right)\end{array}$

[0080] Such a whitening operator can be computed during the iterative minimization step 1108 illustrated in FIG. 11, using, for example, a Levinson algorithm applied to the reconstruction error D(t,)-D.sub.0(t,) illustrated in FIG. 12. If this is the case, then the Bayesian cost function C, after adding the Jacobian term discussed with regard to equation (28) and the noise-whitening operator discussed with regard to equation (29) is given by:

[00024]$\begin{array}{cc}C\left[e_1\left(t\right),e_2\left(t\right),w_n\right]=.Math.\frac{1}{2.Math.N_t}.Math.\underset{t}{\overset{}{.Math.}}.Math.\left[f\left[e_1\left(t\right)\right]+f\left[e_2\left(t\right)\right]\right]+.Math.\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.\log .Math..Math.\left(t,\right)+\frac{1}{N_t.Math.N_{}}.Math.\underset{t,}{\overset{}{.Math.}}.Math.{\left[A_N\left(t\right)*\left(D\left(t,\right)-D_0\left(t,\right)\right)\right]}^2.& \left(31\right)\end{array}$

[0081] The whitening operator A.sub.N(t) have been shown in equation (29) to depend only on time t, however, it is possible to also depend on angle .

[0082] Returning to the method discussed with regard to FIG. 6, after the wavelet model has been calculated in step 608, as discussed with regard to the method of FIG. 11, the process advances to step 610 for processing the recorded data, for example, by applying the wavelet model to the CIG gather for obtaining processed seismic data. The processed seismic data may then be used in step 612 for calculating one or more characteristics of the earth (e.g., density, speed, etc.) and/or an image of the subsurface that was seismically surveyed.

[0083] The methods discussed above have used only the intercept and gradient (two term Shuey equation) for the statistical model. However, the methods may be extended to the three-term Aki-Richards equation (see equation 2), thus using the intercept A(t), the gradient B(t), and the curvature C(t). This means that the source separation is performed on three components, resulting in a statistical model with a three-by-three matrix. Then, the Bayesian inversion discussed in FIG. 11 would use a statistical model with 3 independent whitened components e.sub.1(t), e.sub.2(t) and e.sub.3(t).

[0084] A method for calculating a characteristic of the earth with an improved wavelet model is now discussed with regard to FIG. 16. It is noted that this method not only improves the technological field of imaging the subsurface of the earth for better exploring the oil and gas reserves, but also improves the functionality of a computer that runs such methods by using a wavelet model that is continuous in time and angle. Because the previous methods use wavelet models that are not continuous in time and angle, the mathematical difficulties to be handled by such models require human intervention due to the model discontinuities.

[0085] Further, for conventional narrowband acquired seismic data, the wavelet is weakly varying, as a Ricker-type wavelet with a dominant frequency decreasing with time and angle. For broadband data, which may be used in the methods discussed above, the wavelet becomes strongly dependent on time and angle. This is because the method attempts to push the maximum frequency f.sub.max(t,) of the wavelet as far as possible: this means f.sub.max decrease quite fast with time and angle, due to both the stretch effect and the absorption. Note that for broadband data, the wavelets exhibit more variations, due to the fact that the maximum frequency is very high for small time and angle values, but much smaller for large time and angle values.

[0086] FIG. 16 shows a method 1600 that receives well log data in step 1602. The method calculates in step 1604 the intercept and gradient based on equation (4). Note that the intercept and gradient depend on time and for this reason, these two parameters take many values during a given time interval, as illustrated in FIG. 7. Based on these multiple values, in step 1606, the method computes, for example, by convolutional source separation as described by equation (7), the statistical model of the intercept and gradient in terms of sparseness, coloring wavelet (see FIGS. 8 and 9), and correlation matrix c.sub.ij. The statistical model includes the correlation matrix and a probability density function p(x). In one application, the method used the pdf introduced in equation (8), e.g., the logistic distribution (shown in FIG. 10), which is not a Gaussian distribution. The statistical model includes the anticorrelation of the intercept and gradient components, their blue spectrum (illustrated in FIG. 9) and their sparseness.

[0087] In step 1608, the intercept and gradient are calculated from the statistical model. Note that by calculating the intercept and gradient from the statistical model, which is different from step 1604 in which the intercept and gradient were calculated from an equation, the method makes use of the well in a statistical sense, which means that the well does not have to be at the exact location of the seismic gather. The statistics of a well can be considered to be valid regionally, and not only locally. If no well is available, a generic statistical model can be used and steps 1602 and 1604 may be omitted.

[0088] The method then advances to step 1610 where a wavelet model structure is selected, for example, a parametrically time and angle ARMA model, and the parameters of the ARMA model are the lattice parameters (see, e.g., equation (13)). In step 1612, the time and angle dependent reflectivity is calculated based on equation (16) or (3), using the intercept and gradient calculated from the statistical model in step 1608. Then, in step 1614, the reconstructed gather D(t,) is calculated based on the wavelet model from step 1610 and the reflectivity from step 1612. In step 1616, seismic data acquired during a seismic survey is received and processed to obtain recorded gather D.sub.0(t,). As previously discussed, the gathers noted in this method are common image gathers. However, other type of gathers can be used.

[0089] The method advances to step 1618 in which a first cost function is selected for calculating the parameters w.sub.n of the wavelet function. In one example, the first cost function is a Bayesian cost function, that has (i) a first term (a priori) that is related to the probability density function of the statistical model (see equation (18) and (ii) a second term (a posteriori) that is related to a difference between the recorded gather D.sub.0(t,) and the reconstructed gather D(t,). In step 1620, the parameters w.sub.n of the wavelet model are calculated by Bayesian inversion. The Bayesian inversion may be solved with a conjugate gradient method. Optionally, in step 1622, a noise-whitening operator is calculated during the iterative process of Bayesian inversion. The noise-whitening operator may be calculated with a Levinson algorithm. In step 1624, the scale term of the cost function is normalized, as discussed above with regard to FIG. 11, i.e., a Jacobian term is added to the cost function, for preventing a scale ambiguity (see equation (28)).

[0090] As a result of step 1620, the parameters w.sub.n of the wavelet model W(t,) are calculated and thus, the reconstructed gather D(t,) is known. A second cost function (associated with an elastic inversion of the acquired seismic data) may then be applied in step 1626 to the reconstructed gather for generating the processed data. This processed data may then be used to generate an image of the subsurface. Alternatively, the wavelet model may be used in step 1626 for processing the acquired seismic data for generating one or more characteristics of the earth. For example, the phase only part of the wavelet W may be applied to the acquired seismic data to obtain the processed seismic data and then one or more characteristics of the earth are calculated based on the processed seismic data (using, for example, an inversion process).

[0091] The method discussed above can be used in a variety of other ways. For example, in one application, the independent whitened components e.sub.1(t) and e.sub.2(t) of the statistical model can be used to compute the intercept and gradient A(t) and B(t), which can be used to invert for the elastic parameters. Alternatively, the wavelet model W(t,) can be forwarded to any elastic inversion together with the gathers D.sub.0(t,).

[0092] In addition, the wavelet model can be forwarded together with the modeled gathers D(t,) for generating the image of the subsurface, which is equivalent to performing a noise attenuation or data conditioning. Another alternative is to apply the phase part of the wavelet model to either D.sub.0(t,) or D(t,), which is equivalent to performing a residual moveout (that is a flattening of the gathers), and provide these processed gathers to the seismic inversion with the zero-phase part of the wavelet model. The wavelet model can be estimated in a multichannel way, where multiples CIGs corresponding to surface positions (x,y) covering a given surface contribute to the estimation of a single common wavelet model.

[0093] The embodiments discussed above make uses of the well in a statistical sense, which means the well does not have to be at the exact location of the seismic gather. This means that the statistics of a well can be considered to be valid regionally, and not only locally. If no well is available, a generic statistical model can be used.

[0094] An advantage of one or more of the embodiments discussed above is the continuous variation in time and angle of the gather, which allows to have a more precise wavelet model. The continuous variation in time allows the use of large windows, and therefore a better estimation of the low frequencies, while taking into account the time variability of the wavelet. As the method separates the signal from the noise, no partial angle stacking is necessary and the angle dependency can also be precisely modeled.

[0095] In one embodiment, the wavelet model can estimate the scaling, the normalized amplitude spectrum and the normalized phase spectrum, while existing statistical methods estimate only a normalized amplitude spectrum.

[0096] The seismic data that may be processed with one or more of the above embodiments may include both hydrophone (i.e., pressure) data (H) and particle motion data (V). The particle motion data (V) may be used to directly or indirectly determine a vertical velocity (Vz) for particles proximate to the detectors that provided the particle motion data. The additional particle velocity data may be kinematically the same as hydrophone data (i.e., arrivals at the same times for the same traces), but with a polarity reversal on the migration-ghost event. It should be noted that some sensors, such as particle motion sensors, may be directional and contain only a component of the pressure recordings and such sensors may only be sensitive to acoustic energy travelling along the orientation of the sensor.

[0097] The above-discussed procedures and methods provide improved results over existing techniques and may be implemented in a computing device illustrated in FIG. 17. Hardware, firmware, software, or a combination thereof may be used to perform the various steps and operations described herein. The computing device 1700 of FIG. 17 is an exemplary computing structure that may be used in connection with such a system.

[0098] The exemplary computing device 1700 suitable for performing the activities described in the exemplary embodiments may include a server 1701. Such a server 1701 may include a central processor (CPU) 1702 coupled to a random access memory (RAM) 1704 and to a read only memory (ROM) 1706. The ROM 1706 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. The processor 1702 may communicate with other internal and external components through input/output (I/O) circuitry 1708 and bussing 1710, to provide control signals and the like. The processor 1702 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

[0099] The server 1701 may also include one or more data storage devices, including hard drives 1712, CDROM drives 1714, and other hardware capable of reading and/or storing information such as DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CDROM or DVD 1716, a USB storage device 1718 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as the CDDROM drive 1714, the disk drive 1712, etc. The server 1701 may be coupled to a display 1720, which may be any type of known display or presentation screen, such as LCD displays, plasma display, cathode ray tubes (CRT), etc. A user input interface 1722 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

[0100] The server 1701 may be coupled to other devices, such as sources, detectors, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 1728, which allows ultimate connection to the various landline and/or mobile computing devices.

[0101] The disclosed exemplary embodiments provide a computing device, a method for seismic data processing, and systems corresponding thereto. For example, the disclosed computing device and method could be integrated into a variety of seismic survey and processing systems including land, marine, ocean bottom, and transitional zone systems with either cabled or autonomous data acquisition nodes. It should be understood that this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications, and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

[0102] Although the features and elements of the present exemplary embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

[0103] This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.