CN110146923B - An efficient and high-precision depth-domain seismic wavelet extraction method - Google Patents

Abstract

The invention provides an efficient and high-precision depth domain seismic wavelet extraction method. The method can accurately extract the depth domain seismic wavelets in a short depth window in the constant-speed depth domain, and accords with the stable hypothesis condition of the convolution theory. According to the method, the short depth window and the short limiting matrix are adopted, so that the influence of high-frequency noise in the logging information can be reduced, the data volume participating in operation can be reduced, the requirement on the memory of a computer is reduced, the calculation time is shortened, and the calculation efficiency is improved.

Description

An efficient and high-precision depth-domain seismic wavelet extraction method

technical field

The invention belongs to the field of petroleum seismic exploration, and relates to a method for rapidly and accurately extracting depth-domain seismic wavelets from depth-domain seismic data.

Background technique

In recent years, thanks to the development of computer technology, the depth migration imaging technology of seismic data has been widely used in the field of petroleum seismic exploration. Depth migration imaging is a structural inversion method. Compared with time-migration imaging results, depth-migration imaging results are more accurate in determining subsurface complex structures. Therefore, there is an increasing demand for other depth-domain inversion methods for reservoir description, such as impedance inversion, elastic parameter inversion, etc. Depth-domain seismic wavelets largely control the results of depth-domain inversion. Therefore, the accurate extraction of seismic wavelets in the depth domain is the key. However, the waveform of the seismic wavelet in the depth domain changes with the change of the medium velocity when it propagates underground, and has strong non-stationarity. This makes the extraction of seismic wavelets in the depth domain not as easy as in the time domain. Some domestic scholars have done a lot of exploration in this area. In a limited space, the medium with different layer velocities in the depth domain is converted into a medium with the same layer velocity through the velocity replacement method. This depth domain space after velocity replacement is called constant velocity depth domain space. In the constant velocity depth domain space, the waveform of the seismic wavelet in the depth domain can remain stable, so it can satisfy the "linear time-invariant" assumption of the convolution model and eliminate the influence of non-stationarity. On this basis, we can extract seismic wavelets in constant velocity depth domain based on the convolution model.

However, the sampling interval of the depth domain seismic data is much smaller than that of the time domain. For example, in the logging process, the depth sampling interval will usually be 0.1 m or 0.15 m. If the propagation speed of the seismic wave is 3000 m/s, then its corresponding time sampling interval is about 0.033 msec or 0.05 msec, while the normal high resolution The rate sampling rate is only 1 ms. Therefore, a large amount of data needs to be processed in the process of seismic wavelet extraction in the depth domain, which occupies a large amount of computer memory resources. Although today's computer memory can meet these needs, mathematical operations on large matrices are time-consuming. For example, the time required for inversion of a matrix is not a simple linear relationship but an exponential relationship with the size of the matrix.

SUMMARY OF THE INVENTION

The object of the present invention is to avoid the deficiencies of the prior art, and provides a new method for rapidly and accurately extracting depth-domain seismic wavelets from depth-domain seismic data, the method comprising the following main steps:

(1) Given a velocity value, and transform the velocity according to the velocity, respectively map the logging data and the side-hole seismic records from the depth domain to the constant velocity depth domain;

(2) Using the density and P-wave velocity information in the logging data in the constant velocity depth domain, the reflection coefficient r in the constant velocity depth domain is calculated;

(3) According to the number of sampling points M of seismic wavelets, an integer value a is given, and a depth window for seismic wavelet extraction is selected, and the window length is N=aM;

(4) Construct a reflection coefficient matrix R with a Toeplitz structure using the constant velocity depth domain reflection coefficient r _w in the depth window, and the size of the matrix R is N×N;

⑸ Construct the restriction matrix P according to the following formula:

In the formula, the size of the matrix P is N×M;

(6) Given the initial seismic wavelet w, given the initial coefficient matrix Φ, and given the threshold ε, the coefficient matrix Φ has the following structure:

Φ=diag[φ ₁ φ ₂ … φ _i … φ _N ],

Among them, diag[ ] represents the diagonal matrix symbol, and φ _i is the i-th diagonal element in the coefficient matrix Φ;

⑺According to the following formula, the synthetic seismic record of constant velocity depth domain is obtained

Wherein, R _p =RP;

与常速度深度域井旁地震记录s的残差绝对值η：⑻ Computational synthetic seismic records

The absolute value η of the residual error with the seismic record s beside the wellbore in the constant velocity depth domain:

If the L ₂ norm of η is less than or equal to the termination condition, the final seismic wavelet is output to terminate the cycle; if the L ₂ norm of η is greater than the termination condition, the coefficient matrix needs to be updated according to the following formula:

where η _i is the ith element in the residual vector η;

At the same time, the seismic wavelet needs to be updated according to the following formula:

Among them, T represents the transpose operation of the matrix;

⑼ If the set number of cycles is reached, terminate the cycle, otherwise return to step ⑺.

Description of drawings

FIG. 1 is a one-dimensional forward modeling model according to an embodiment of the present invention. Among them, Figure 1(a) is the density logging curve in the actual logging data in a certain area, the abscissa is the density, the unit is g/cm3, the ordinate is the depth, the unit is meters, and the depth is from 3550 meters to 3975 meters. . Figure 1(b) is the P-wave velocity logging curve, the abscissa is the velocity, the unit is km/s, and the ordinate is the depth, the unit is m. Fig. 1(c) is the reflection coefficient calculated according to the logging density and the logging longitudinal wave velocity. The ordinate is the depth, and the unit is meters. Figure 1(d) is the reflection coefficient transformed to the constant velocity depth domain, the velocity value used for constant velocity transformation is 3000 m/s, and its ordinate is the constant velocity depth, in meters. Fig. 1(e) is the constant-velocity depth-domain seismic record synthesized by the reflection coefficient of Fig. 1(d) and the constant-velocity depth-domain seismic wavelet shown in Fig. 1(f), and the ordinate is the constant-velocity depth, The unit is meters. The data in the black dotted rectangles in Fig. 1(d) and Fig. 1(e) are used for seismic wavelet extraction. Figure 1(f) shows the constant velocity depth domain Ricker wavelet when the constant velocity is 3000 m/s, the corresponding time domain Ricker wavelet has a dominant frequency of 43 Hz, and the ordinate is the length in meters.

FIG. 2 is a result of extracting seismic wavelets in a constant velocity depth domain on the forward modeling model shown in FIG. 1 according to an embodiment of the present invention. Among them, the black line is the known Ricker wavelet, the gray line is the extracted seismic wavelet, the abscissa is the length, the unit is meters, and the ordinate is the amplitude.

FIG. 3 is a result of extracting seismic wavelets in a constant velocity depth domain by using the actual logging data and seismic data of Well A in a certain work area according to an embodiment of the present invention. Among them, Fig. 3(a) is the reflection coefficient calculated according to the logging data of Well A, the ordinate is the depth, the unit is meters, and the depth is from 3435 meters to 3845 meters. Fig. 3(b) is the comparison between the seismic record gathers (black) corresponding to the well A and the synthetic seismic records in the depth domain (gray, indicated by arrows) made by using the extracted seismic wavelets. The ordinate is the depth, and the unit is meters. . Figures 3(c) and 3(d) are the reflection coefficients transformed to the constant velocity depth domain and the wellbore seismic record gathers (among them, the wellbore seismic records used for seismic wavelet extraction are gray, indicated by arrows), using The velocity value transformed at constant velocity is 3000 m/s, and the ordinate is the constant velocity depth, in meters. Seismic wavelets are extracted using the data in the black dotted rectangles in Fig. 3(c) and Fig. 3(d). Figure 3(e) is the extracted seismic wavelet in the constant velocity depth domain, the ordinate is the length, and the unit is meters.

Detailed ways

(1) Given a velocity value, and transform the velocity according to the velocity, respectively map the logging data and the side-hole seismic records from the depth domain to the constant velocity depth domain;

(2) Using the density ρ and P-wave velocity v information in the logging data in the constant velocity depth domain, the reflection coefficient r in the constant velocity depth domain is calculated according to the following formula:

(3) According to the wavenumber spectrum of the seismic records in the constant velocity depth domain or by the trial and error method, the sampling points M of the seismic wavelets in the constant velocity depth domain are estimated; the sampling points of the seismic wavelets are multiplied by the depth sampling interval in the constant velocity depth domain to calculate the length of the seismic wavelet;

(4) According to the number of sampling points M of seismic wavelets, an integer value a is given, and a depth window for seismic wavelet extraction is selected, and the window length is N=aM; the given integer value a needs to satisfy: a≥3 and aM≤L, where L is the total number of sampling points of the reflection coefficient r in the constant velocity depth domain;

(5) Construct a reflection coefficient matrix R with a Toeplitz structure using the constant velocity depth domain reflection coefficient r _w in the depth window, and the size of the matrix R is N×N;

(6) Construct the restriction matrix P according to the following formula:

In the formula, the size of the matrix P is N×M;

⑺ Use the least squares method to determine the initial seismic wavelet w:

表示矩阵的广义逆运算；In the formula, R _p = RP,

Represents the generalized inverse operation of a matrix;

Given a threshold ε and an initial coefficient matrix Φ:

⑻ According to the following formula, the synthetic seismic records of constant velocity depth domain are obtained

与常速度深度域井旁地震记录s的残差绝对值η：⑼ Computational synthetic seismic records

The absolute value η of the residual error with the seismic record s beside the wellbore in the constant velocity depth domain:

If the L ₂ norm of η is less than or equal to the termination condition, the final seismic wavelet is output to terminate the cycle; if the L ₂ norm of η is greater than the termination condition, the coefficient matrix needs to be updated according to the following formula:

Among them, φ _i is the ith diagonal element in the coefficient matrix Φ; η _i is the ith element in the residual vector η;

At the same time, the seismic wavelet needs to be updated according to the following formula:

Among them, T represents the transpose operation of the matrix;

⑽ If the set number of cycles is reached, terminate the cycle, otherwise return to step ⑻.

In the one-dimensional forward model shown in Figure 1, the sampling points of the reflection coefficient r in the constant velocity depth domain are L=3176, the sampling points of the Ricker wavelet in the constant velocity depth domain are M=681, and the given integer value a=3 , the length of the depth window shown by the black dotted rectangle is aM=N=2043;

Figure 2 shows the result of extracting seismic wavelets in constant velocity depth domain for the forward modeling model shown in The seismic wavelet extracted by the embodiment of the present invention is consistent with the known seismic wavelet.

In the actual data shown in Figure 3, the number of sampling points of the reflection coefficient r in the constant velocity depth domain is L=2394, the number of sampling points of the Ricker wavelet in the constant velocity depth domain is M=437, the given integer value a=4, the black dotted line The length of the depth window shown in the rectangular box is aM=N=1748; given the threshold ε=0.00001, the termination condition is 0.0001, and the set number of cycles is 100, the results of seismic wavelet extraction in the constant velocity depth domain are shown in Figure 3 ( e) shown. The synthetic seismic record in the depth domain obtained by the inverse transformation of the synthetic seismic record made with the constant velocity depth domain seismic wavelet is shown in Fig. 3(b). The depth-domain synthetic seismic records produced by the seismic wavelet of , are in good agreement with the wellbore depth-domain seismic records, and the correlation coefficient between the two is 0.85.

The advantages of the present invention are: (1) in the constant velocity depth domain, the accurate extraction of seismic wavelets in the depth domain can be realized with the data in the shorter depth window, which conforms to the stationary assumption of the convolution theory; (2) by using The shorter depth window and restriction matrix P can not only reduce the influence of high-frequency noise in the logging information, but also reduce the amount of data involved in the operation (for example, reduce the size of the reflection coefficient matrix R), thereby reducing the demand for computer memory, Shorten computing time and improve computing efficiency.

The above embodiments are only used to illustrate the present invention, and each implementation step of the method can be changed to some extent, and all equivalent transformations and improvements carried out on the basis of the technical solutions of the present invention should not be excluded from the present invention. outside the scope of protection.

Claims (1)

1. An efficient high-precision depth domain seismic wavelet extraction method mainly comprises the following steps:

the method includes the steps of giving a speed value, carrying out speed transformation according to the speed, and mapping logging data and well side seismic records from a depth domain to a constant speed depth domain;

secondly, calculating a constant velocity depth domain reflection coefficient r by utilizing density and longitudinal wave velocity information in the logging data of the constant velocity depth domain;

thirdly, extracting seismic wavelets in the constant velocity depth domain by using the reflection coefficient r in the constant velocity depth domain and the well side seismic record s; in the step three, extracting the seismic wavelet in the constant velocity depth domain is carried out according to the following steps:

(a) according to the number M of sampling points of the seismic wavelet, giving an integer value a, and selecting a depth window for extracting the seismic wavelet, wherein the window length is N-aM; a given integer value a needs to satisfy: a is more than or equal to 3, aM is less than or equal to L, wherein L is the total sampling point number of the reflection coefficient r of the constant-speed depth domain;

(b) by constant velocity depth-domain reflection coefficient r within a depth window_wConstructing a reflection coefficient matrix R with a Toeplitz structure, wherein the size of the matrix R is N×N；

(c) The restriction matrix P is constructed as follows:

in the formula, the size of the matrix P is NxM;

(d) giving an initial seismic wavelet w, a threshold value and an initial coefficient matrix phi, wherein the coefficient matrix phi has the following structure:

Φ＝diag[φ₁φ₂… φ_i… φ_N]，

wherein, diag [. cndot]Representing diagonal matrix symbols, phi_iIs the ith diagonal element in the coefficient matrix Φ;

(e) obtaining a constant velocity depth domain synthetic seismic record according to the formula

In the formula, R_p＝RP；

(f) Computing synthetic seismic records

Residual absolute value η from the constant velocity depth domain well side seismic record s:

if L of η₂Outputting the final seismic wavelet and stopping the cycle if the norm is less than or equal to the stop condition, and if the L of η is less than or equal to the stop condition₂If the norm is greater than the termination condition, the coefficient matrix needs to be updated according to the following formula:

wherein, η_iIs the ith element in residual vector η;

meanwhile, the seismic wavelets need to be updated according to the following formula:

wherein T represents a transpose operation of a matrix;

(g) if the set circulation times are reached, the circulation is terminated, otherwise, the step (e) is returned to.

Images