Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain

Abstract

A method for estimating the depth-domain seismic wavelets from depth-domain seismic data and synthesizing depth-domain seismic records. The method includes: obtaining depth coordinates and P-wave velocity v and density from well log, calculating a corresponding reflectivity series r; performing constant-velocity depth conversion for a seismic trace S and a reflectivity series r by using a velocity v.sub.c as a reference velocity to obtain the converted seismic trace S.sub.1 and the converted reflectivity series r.sub.1; and estimating a depth-domain seismic wavelet w based on the Gibbs sampling method; synthesizing depth-domain seismic record by using the P-wave v, the reflectivity series r and the estimated depth-domain seismic wavelet w.

Claims

1. An integrated method for estimation of depth-domain seismic wavelets and synthesis of seismic records in depth domain, the method comprising: a) drilling a borehole to perform well logging such that a well log is obtained; b) obtaining depth coordinates and P-wave velocity v and density from the well log, and calculating a corresponding reflectivity series r; c) performing constant-velocity depth conversion for a seismic trace S and a reflectivity series r by using a velocity v.sub.c as a reference velocity to obtain the converted seismic trace S.sub.1 and the converted reflectivity series r.sub.1; and estimating the depth domain seismic wavelet w based on the Gibbs sampling method; and d) synthesizing depth-domain seismic records by using the P-wave velocity v and the reflectivity series r obtained in step 1) and the estimated depth-domain seismic wavelet w in step 2).

2. The integrated method of claim 1, wherein estimating the depth-domain seismic wavelet w based on the Gibbs sampling method in the step c) further comprises: performing the constant-velocity depth conversion according to following formula: $\Delta d_c=\Delta d\times \frac{v_c}{v_{\max }};$ where, Δd is a depth sampling interval of well log, v.sub.c is a user-defined constant velocity, v.sub.max, is a maximum velocity recorded in well log, and Δd.sub.c is a depth sampling interval after the constant-velocity depth conversion; resampling the seismic trace S and the reflectivity series r at the converted sampling interval Δd.sub.c to obtain the converted seismic trace S.sub.1 and the converted reflectivity series r.sub.1; wherein a convolution model for estimating the seismic wavelet is as follows:
S.sub.1=r.sub.1*w+e; where, S.sub.1 and r.sub.1 are the seismic trace and the reflectivity series in the constant-velocity depth domain, respectively, “*” represents the convolution operation; w is the seismic wavelet in the constant-velocity depth domain; e is the error term; where, the error term e meets the multivariate normal distribution with a mean vector of 0 and a covariance matrix of Σ.sub.s; the seismic trace S.sub.1 in the constant-velocity depth domain satisfies a following normal distribution:
S.sub.1|μ.sub.s,Σ.sub.s˜N.sub.m(μ.sub.s,Σ.sub.s); where, N represents a normal distribution, m is the number of sampling points of S.sub.1, μ.sub.s=r.sub.1*w, Σ.sub.s=σ.sub.s.sup.2M.sub.s; where M.sub.s is a correlation coefficient matrix constructed by S.sub.1, σ.sub.s.sup.2 satisfies an inverse-gamma distribution:
σ.sub.s.sup.2˜IGa(α.sub.s,λ.sub.s); where, α.sub.s is a shape parameter, λ.sub.s is a scale parameter; wherein the seismic wavelet w estimated in the constant-velocity depth domain satisfies following normal distribution:
w|μ.sub.w,Σ.sub.w˜N.sub.o(μ.sub.w,Σ.sub.w); wherein, N represents the normal distribution and o is the number of sampling points of w, μ.sub.w meets a multivariate normal distribution with a mean vector of 0 and a covariance matrix of Σ.sub.wa=10000M.sub.w, Σ.sub.w=σ.sub.w.sup.2M.sub.w; where M.sub.w is a correlation coefficient matrix constructed from a known wavelet; σ.sub.w.sup.2 also satisfies an inverse-gamma distribution:
σ.sub.w.sup.2˜IGa(α.sub.w,λ.sub.w); where, α.sub.w is a shape parameter, λ.sub.w is a scale parameter; unknown parameters involved in the depth-domain seismic wavelet estimation model are: w, μ.sub.w, σ.sub.w.sup.2 and σ.sub.s.sup.2; according to the Bayesian theorem, a posterior distribution of the unknown parameters is expressed as:
p(w,μ.sub.w,σ.sub.w.sup.2,σ.sub.s.sup.2|s,r)∝p(s|r,w,σ.sub.s.sup.2)p(w|μ.sub.w,σ.sub.w.sup.2)p(μ.sub.w)p(σ.sub.w.sup.2)p(σ.sub.s.sup.2); estimating the depth-domain seismic wavelet w based on the Gibbs sampling method further comprises: 1) giving the initial values of w.sup.(0), μ.sub.w.sup.(0), σ.sub.w.sup.2(0), σ.sub.s.sup.2(0) to w, μ.sub.w, σ.sub.w.sup.2, σ.sub.s.sup.2 respectively; 2) repeating the following steps at a given number of cycles t=1, 2, . . . : i) generating a candidate value w.sup.(t) from following full conditional distribution:
f(w|μ.sub.w.sup.(t-1),σ.sub.w.sup.2(t-1),σ.sub.s.sup.2(t-1),s,r):N.sub.o(μ.sub.wf,Σ.sub.wf); ii) generating a candidate value μ.sub.w.sup.(t) from following full conditional distribution:
f(μ.sub.w|w.sup.(t),σ.sub.w.sup.2(t-1),σ.sub.s.sup.2(t-1),s,r):N.sub.o(μ.sub.μ.sub.w.sub.f,Σ.sub.μ.sub.w.sub.f); iii) generating a candidate value σ.sub.w.sup.2(t) from following full conditional distribution:
f(σ.sub.w.sup.2|w.sup.(t),μ.sub.w.sup.(t),σ.sub.s.sup.2(t-1),s,r):IGa(α.sub.wf,λ.sub.wf); and iv) generating a candidate value σ.sub.s.sup.2(t) from following full conditional distribution:
f(σ.sub.s.sup.2|w.sup.(t),μ.sub.w.sup.(t),σ.sub.w.sup.2(t),s,r):IGa(α.sub.sf,λ.sub.sf); and obtaining the estimated depth-domain seismic wavelet if the candidate value w.sup.(t) meets the convergence condition.

3. The integrated method of claim 1, wherein the step d) further comprises: 1) using the calculated reflectivity series r as the initial input data p(z;0)=[r(0), . . . , r(n), . . . , r(D)].sup.T; performing the Fourier transform to the depth-domain wavelet w to obtain the wave number domain wavelet W; 2) calculating e.sup.−ik.sup.z.sup.Δd; where, $k_z=\frac{4\pi \times \left(j-1\right)}{0.002\times v\left(n\right)\times t},$ j is the j.sub.th frequency sampling point, v(n) is a velocity at the n.sub.th depth sampling point, t is the number of frequency sampling points, Δd is a depth sampling interval; multiplying e.sup.−ik.sup.z.sup.Δd with p(n−1;j) to obtain p(n;j):
p(n;j)=e.sup.−ik.sup.z.sup.Δd×p((n−1);j); 3) returning to step 2) and continuing the depth cycle until a maximum depth is reached, so that the temporary variable p(j) at current frequency j is obtained as follows:
p(j)=p(D;j)×W(j); where D is the number of the depth sampling points; 4) taking another frequency j+1, and repeating steps 2) and 3) until all frequencies have been repeated, so that following temporary array is obtained:
p′=[p(1), . . . ,p(j), . . . ,p(t)]; 5) initializing the values of frequency and depth; taking a frequency of j, when the depth is n, calculating e.sup.ik.sup.z.sup.Δd and multiplying e.sup.ik.sup.z.sup.Δd with p(j), so that a temporary variable P(n;j) is obtained:
P(n;j)=e.sup.ik.sup.z.sup.Δd×p(j); 6) taking another frequency j+1, returning to step 5) and continuing the frequency cycle until the maximum frequency is reached; so that the value of synthetic seismogram at depth n is obtained as follows: $P\left(n\right)=\underset{j=1}{\overset{j=t}{.Math.}}P\left(n;j\right);$ and 7) taking another depth n+1, repeating steps 5) and 6) until a maximum depth is reached; a final synthetic seismogram P is obtained:
P=[P(1), . . . ,P(n), . . . ,P(D)].

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a diagram of one-dimensional multi-layer impedance model; wherein the abscissa represents the impedance in units of m/s×g/cc, and the ordinate represents depth in units of meters; the depth sampling interval is 1 meter.

(2) FIG. 2 is a diagram of depth-domain seismic records synthesized by using the conventional approach; wherein the ordinate represents depth in units of meters; and the abscissa represents the impedance in units of m/s×g/cc.

(3) FIG. 3 is a diagram of depth-domain synthetic seismic records synthesized by the method of the present disclosure; wherein the ordinate represents depth in units of meters; and the abscissa represents the impedance in units of m/s×g/cc.

(4) FIG. 4 is a diagram of Ricker wavelet in the constant-velocity depth domain, which is used for making depth-domain synthetic seismograms in FIG. 1, with the velocity of 3000 m/s and a dominant frequency of 50 Hz under theoretical conditions.

(5) FIG. 5 is a diagram of depth-domain seismic wavelet estimated from a borehole-side seismic trace; wherein the abscissa represents the length in the unit of meters, and the ordinate represents the amplitude.

(6) FIG. 6 is a diagram of comparison of the depth-domain synthetic seismic records with the borehole-side seismic traces; wherein the ordinate represents the depth in the unit of meters, and the range of the depth is from 3550 meters to 3975 meters.

DETAILED DESCRIPTION OF EMBODIMENTS

(7) (2) Calculating a reflectivity series form the well log as follow:

(8) $r_i=\frac{{\rho }_{i+1}{v}_{i+1}-{\rho }_i{v}_i}{{\rho }_{i+1}{v}_{i+1}+{\rho }_i{v}_i};$

(9) where, r.sub.1 is a depth-domain reflection coefficient between the i.sub.th sampling point and the (i+1).sub.th sampling point, ρ.sub.i is a density of the i.sub.th sampling point and v.sub.i is a velocity of the i.sub.th sampling point.

(10) (2) Selecting an appropriate borehole-side seismic trace from prestack depth migration seismic data and defining the velocity v.sub.c as a standard velocity, performing constant-velocity depth conversion according to the following formula to satisfy the linear time-invariant condition:

(11) $\Delta d_c=\Delta d\times \frac{v_c}{v_{\max }};$

(12) where Δd is a depth sampling interval of well log, v.sub.max is the maximum velocity recorded in well log, and Δd.sub.c is a depth sampling interval after the constant-velocity depth conversion. The borehole-side seismic trace S and reflectivity series r are resampled at the converted sampling interval Δd.sub.c to obtain the converted seismic trace S.sub.1 and reflectivity series r.sub.1.

(13) Before performing constant-velocity depth conversion, the true depth h is:

(14) $h=\underset{i=1}{\overset{i=l}{.Math.}}\Delta d;$

(15) where, l is the number of depth sampling points of well log.

(16) The depth h after constant-velocity depth conversion is h.sub.c:

(17) $h_c=\underset{i=1}{\overset{i=l}{.Math.}}\Delta d_c;$

(18) (3) A convolution model in constant-velocity depth domain is as follows:
S.sub.1=r.sub.1*w+e;

(19) where, S.sub.1 and r.sub.1 are the seismic trace and reflectivity series in the constant-velocity depth domain, respectively; “*” represents the convolution operation; w is the seismic wavelet to be estimated in constant-velocity depth domain; e is an error term.

(20) For the error term e in the above equation, we assume that it meets a multivariate normal distribution with a mean vector of 0 and a covariance matrix of Σ.sub.s. For seismic trace S.sub.1 in the constant-velocity depth domain, we assume that it satisfies following normal distribution:
S.sub.1|μ.sub.s,Σ.sub.s˜N.sub.m(μ.sub.s,Σ.sub.s);

(21) where, N represents the normal distribution, m is the number of sampling points of S.sub.1, μ.sub.s=r.sub.1*w, Σ.sub.s=σ.sub.s.sup.2M.sub.s, where M.sub.s is a correlation coefficient matrix constructed by S.sub.1 and is a positive definite matrix; σ.sub.s.sup.2 is a variance factor which satisfies an inverse-gamma distribution:
σ.sub.s.sup.2˜IGa(α.sub.s,λ.sub.s)

(22) in the formula, α.sub.s is the shape parameter, λ.sub.s is a scale parameter.

(23) For the seismic wavelet w estimated in the constant-velocity depth domain, we assume that it satisfies following normal distribution:
w|μ.sub.w,Σ.sub.w˜N.sub.o(μ.sub.w,Σ.sub.w)

(24) where, N represents the normal distribution and o is the number of sampling points of w, μ.sub.w meets a multivariate normal distribution with a mean vector of 0 and a covariance matrix of Σ.sub.wa=10000M.sub.w. Σ.sub.w=σ.sub.w.sup.2M.sub.w, where M.sub.w is a correlation coefficient matrix constructed from a known wavelet (e.g. the Ricker wavelet), σ.sub.w.sup.2 also satisfies an inverse gamma distribution:
σ.sub.w.sup.2˜IGa(α.sub.w,λ.sub.w);

(25) where, α.sub.w is a shape parameter, λ.sub.w is a scale parameter.

(26) Unknown parameters involved in the depth-domain seismic wavelet estimation model of the method of the present disclosure are: w, μ.sub.w, σ.sub.w.sup.2 and σ.sub.s.sup.2. According to the Bayesian theorem, a posterior distribution of the unknown parameters is expressed as:
p(w,μ.sub.w,σ.sub.w.sup.2,σ.sub.s.sup.2|s,r)∝p(s|r,w,σ.sub.s.sup.2)p(w|μ.sub.w,σ.sub.w.sup.2)p(μ.sub.w)p(σ.sub.w.sup.2)p(σ.sub.s.sup.2);

(27) The steps to estimate the depth-domain seismic wavelet based on the Gibbs sampling method are as follows:

(28) 1) giving the initial values of w.sup.(0), μ.sub.w.sup.(0), σ.sub.w.sup.2(0), σ.sub.s.sup.2(0) to w, μ.sub.w, σ.sub.w.sup.2, σ.sub.s.sup.2 respectively;

(29) 2) repeating the following steps at a given number of cycles t=1, 2, . . . :

(30) a) generating a candidate value w.sup.(t) from following full conditional distribution:
f(w|μ.sub.w.sup.(t-1),σ.sub.w.sup.2(t-1),σ.sub.s.sup.2(t-1),s,r):N.sub.o(μ.sub.wf,Σ.sub.wf);

(31) where, μ.sub.wf=μ.sub.w+Σ.sub.wR.sup.T(RΣ.sub.wR.sup.T+Σ.sub.s).sup.−1(S.sub.1−Rμ.sub.w), T stands for the transpose, R is a Toeplitz matrix constructed from the reflectivity series r.sub.1, Σ.sub.wf=Σ.sub.w−Σ.sub.wR.sup.T(RΣ.sub.wR.sup.T+Σ.sub.s).sup.−1RΣ.sub.s;

(32) b) generating a candidate value μ.sub.w.sup.(t) from following full conditional distribution:
f(μ.sub.w|w.sup.(t),σ.sub.w.sup.2(t-1),σ.sub.s.sup.2(t-1),s,r):N.sub.o(μ.sub.μ.sub.w.sub.f,Σ.sub.μ.sub.w.sub.f);

(33) where, μ.sub.μ.sub.w.sub.f=Σ.sub.wa(Σ.sub.wa+Σ.sub.w).sup.−1w, Σ.sub.μ.sub.w.sub.f=Σ.sub.wa−Σ.sub.wa(Σ.sub.wa+Σ.sub.w).sup.−1Σ.sub.wa;

(34) c) generating candidate value σ.sub.w.sup.2(t) from following full conditional distribution:

(35) $f\left({\sigma }_w^2|w^{\left(t\right)},{\mu }_w^{\left(t\right)},{\sigma }_s^{2\left(t-1\right)},s,r\right)\text{:}\mathrm{IGa}\left({\alpha }_{\mathrm{wf}},{\lambda }_{\mathrm{wf}}\right);$$\mathrm{where},{\alpha }_{\mathrm{wf}}={\alpha }_w+\frac{o}{2},{\lambda }_{\mathrm{wf}}={\lambda }_w+\frac{{\left(w-{\mu }_w\right)}^T{M}_w^{-1}\left(w-{\mu }_w\right)}{2};$

(36) d) generating a candidate value σ.sub.s.sup.2(t) from following full conditional distribution:

(37) $f\left({\sigma }_s^2|w^{\left(t\right)},{\mu }_w^{\left(t\right)},{\sigma }_w^{2\left(t\right)},s,r\right)\text{:}\mathrm{IGa}\left({\alpha }_{\mathrm{sf}},{\lambda }_{\mathrm{sf}}\right);$$\mathrm{where},{\alpha }_{\mathrm{sf}}={\alpha }_s+\frac{m}{2},{\lambda }_{\mathrm{sf}}={\lambda }_s+\frac{{\left(S_1-r_1*w\right)}^T{M}_s^{-1}\left(S_1-r_1*w\right)}{2};$

(38) Obtaining the estimated depth-domain seismic wavelet if the candidate value w meets the convergence condition.

(39) (4) using the calculated reflectivity series r as the initial input data p(z;0)=[r(0),L,r(n),L,r(D)].sup.T;

(40) performing a Fourier transform to the depth-domain wavelet w to obtain the wave number domain wavelet W;

(41) (5) calculating e.sup.−ik.sup.z.sup.Δd;

(42) where,

(43) 0$k_z=\frac{4\pi \times \left(j-1\right)}{0.002\times v\left(n\right)\times t},$
j is the j.sub.th frequency sampling point, v(n) is a velocity at n.sub.th depth sampling point, t is the number of frequency sampling points, Δd is a depth sampling interval;

(44) multiplying e.sup.−ik.sup.z.sup.Δd with p(n−1;j) to obtain p(n;j); where
p(n;j)=e.sup.−ik.sup.z.sup.Δd×p((n−1);j);

(45) (6) Returning to step (5) and continuing the depth cycle until a maximum depth is reached; so that the temporary variable p(j) at current frequency j is obtained as follow:
p(j)=p(D;j)×W(j); where D is the number of the depth sampling points;

(46) (7) taking another frequency of j+1, repeating steps (5) and (6) until all frequencies have been repeated, so that the following temporary array is obtained:
p′=[p(1); . . . ;p(j); . . . ;p(t)];

(47) (8) initializing the values of frequency and depth; taking a frequency of j, when the depth is n, calculating e.sup.ik.sup.z.sup.Δd and multiplying e.sup.ik.sup.z.sup.Δd with p(j), so that a temporary variable P(n; j) is obtained:
P(n;j)=e.sup.ik.sup.z.sup.Δd×p(j);

(48) (9) taking another frequency j+1, returning to step (8) and continuing the frequency cycle until the maximum frequency is reached; so that the value of synthetic seismogram at depth n is obtained as follow:

(49) $P\left(n\right)=\underset{j=1}{\overset{j=t}{.Math.}}P\left(n;j\right);$

(50) (10) taking another depth n+1, repeating steps (8) and (9) until the maximum depth is reached; a final synthetic seismic record P is obtained:
P=[P(1), . . . ,P(n), . . . ,P(D)].

(51) It should be understood that for those of ordinary skills in the art, improvements or variations can be made based on the above descriptions, and such improvements and variations fall within the scope of the appended claims.

(52) The embodiments are only illustrative of the present disclosure, and apparently the implementations are not limited by the above modes. The embodiments described herein and various modifications based on the ideas and technical solutions of the present disclosure fall within the scope of the present application.