A three-dimensional reflection coefficient inversion method with jitter artifact suppression function

Abstract

The application provides a three-dimensional reflection coefficient inversion method with a jitter artifact suppression function, which first obtains a seismic wavelet by using a damped least square inversion method in combination with logging data and a seismic trace beside a well; then constructs a constraint term by using a horizontal second-order difference operator to add a three-dimensional reflection coefficient inversion objective function; and finally splits the inversion objective function into several sub-problems in a frequency domain according to a split Bregman optimization framework, and solves the sub-problems alternately to obtain a three-dimensional reflection coefficient inversion result. The method adds the correlation between adjacent seismic traces in a three-dimensional space to the inversion objective function through a second-order difference, can suppress jitter artifacts in a main survey line direction and a contact survey line direction at the same time, can significantly enhance the traceability of top and bottom interfaces of thin strata and small structure features in a profile and a plane in strong noise-containing seismic data, and further improves the precision of geological interpretation by using three-dimensional reflection coefficients.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas seismic exploration data processing, and particularly relates to a three-dimensional reflection coefficient inversion method with all-round jitter artifact suppression function. Background Technology

[0002] Reflectance coefficient data volumes have wide applications in interpreting tectonic boundaries and identifying the top and bottom interfaces of thin strata; therefore, three-dimensional reflection coefficients are considered a bridge between seismic data and subsurface geological structures. Currently, there are two methods for obtaining three-dimensional reflection coefficient data volumes: trace-by-trace one-dimensional inversion and line-by-line two-dimensional inversion. The trace-by-trace one-dimensional inversion method only uses sparse constraints in the longitudinal time axis direction to establish the inversion objective function. Because it does not utilize the high similarity between adjacent seismic traces in the main seismic line and connecting seismic lines, when the signal-to-noise ratio is low, the same interface exhibits slight time shifts in adjacent seismic traces, causing obvious fluctuations and instability in the overall shape of the reflection coefficients. In reality, subsurface interfaces generally exhibit a smooth and undulating state. These fluctuation artifacts can make the boundaries of subsurface structures and the top and bottom interfaces of thin strata very blurry, causing significant difficulties in geological interpretation.

[0003] To suppress jitter artifacts in 3D reflection coefficient data volumes, various line-by-line inversion methods have been proposed. These methods generally incorporate the correlation information between adjacent seismic traces along the main seismic line into the inversion objective function to achieve jitter suppression. Then, the final 3D reflection coefficient volume is obtained by performing 2D reflection coefficient inversion on each seismic line. Among these line-by-line inversion methods, some utilize Markov random fields to simulate the correlation between adjacent seismic traces, while others utilize the local linearity of phase axes in the time-space domain. f - k The sparsity in the domain reflects this correlation. Some methods use the covariance matrix to describe the lateral correlation between seismic traces, while others use the dip information of the phase axis extracted from the seismic data to establish the correlation between adjacent traces. However, these methods for describing the correlation between seismic traces only involve the main survey line direction. The correlation between seismic traces in the connecting survey line direction is not included in the inversion objective function. Moreover, these methods are difficult to directly extend to three-dimensional space. Therefore, when the signal-to-noise ratio is low, the jitter problem of the reflection coefficient is still very prominent, reducing the accuracy of geological interpretation. Summary of the Invention

[0004] To address the aforementioned problems in the prior art, this invention provides a three-dimensional reflection coefficient inversion method with jitter artifact suppression functionality. This method utilizes a lateral second-order difference operator to incorporate the correlation between adjacent seismic traces in three-dimensional space as a lateral constraint into the inversion objective function, thereby simultaneously achieving jitter suppression functionality for both the survey line direction and the connecting survey line direction, thus solving the problems existing in the prior art. Specifically, it includes:

[0005] Acquire seismic data volumes and P-wave velocity and density logging curves;

[0006] The well logging reflection coefficient is calculated using the P-wave velocity and density logging curves. The wavelet inversion objective function is constructed using the wellside traces and well logging reflection coefficients in the seismic data volume. The seismic wavelet sequence is calculated using the damped least squares method.

[0007] Based on the dimensions of the earthquake data volume and the calculated earthquake wavelet sequence, construct a wavelet three-dimensional convolution operator and a transverse second-order difference three-dimensional convolution operator, and calculate the three-dimensional Fourier transform of the above three-dimensional convolution operator;

[0008] The three-dimensional reflection coefficient inversion objective function with lateral constraint terms is constructed using the aforementioned lateral second-order difference three-dimensional convolution operator;

[0009] The objective function for the inversion of the three-dimensional reflection coefficients is split in the frequency domain using the split Bregman algorithm framework, and the sub-problems are alternately optimized to obtain the inversion results of the three-dimensional reflection coefficients. Attached Figure Description

[0010] The accompanying drawings, which form part of the specification, provide an embodiment of the method described in this patent. They are intended to illustrate the principle and effects of the method more intuitively. Obviously, those skilled in the art can obtain other drawings based on these drawings without any inventive effort. In the drawings:

[0011] Figure 1 A flowchart illustrating a three-dimensional reflection coefficient inversion method with jitter artifact suppression function provided by the present invention;

[0012] Figure 2 This is a schematic diagram illustrating the construction of a three-dimensional wavelet convolution operator using seismic wavelet sequences, provided in an embodiment of the present invention.

[0013] Figure 3 This is a schematic diagram illustrating the construction of a three-dimensional transverse second-order difference convolution operator provided in an embodiment of the present invention;

[0014] Figure 4 A schematic diagram of a three-dimensional seismic data volume provided in an embodiment of the present invention;

[0015] Figure 5 This is a schematic diagram of the three-dimensional reflection coefficient inversion results provided in an embodiment of the present invention;

[0016] Figure 6 This is a cross-sectional view of the reflection coefficient inversion results provided in an embodiment of the present invention;

[0017] Figure 7 This is a time slice diagram of the reflection coefficient inversion results provided in an embodiment of the present invention. Detailed Implementation

[0018] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments selected herein are only used to explain the present invention and are not intended to limit the present invention; that is, the provided embodiments are merely one embodiment of the present invention, and not all embodiments.

[0019] The following is combined Figure 1-7 The present invention will be described in detail below.

[0020] Figure 1 This is a flowchart of a three-dimensional reflection coefficient inversion method with jitter artifact suppression function provided by an embodiment of the present invention, as shown below. Figure 1 As shown, it includes:

[0021] Step S101: Acquire seismic data volume and P-wave velocity and density logging curves;

[0022] Step S102: Calculate the well logging reflection coefficient using the well logging curve, construct the wavelet inversion objective function using the well bypass channel and well logging reflection coefficient in the seismic data volume, and calculate the seismic wavelet sequence using the damped least squares method;

[0023] Step S103: Construct the frequency domain wavelet three-dimensional convolution operator and the transverse second-order difference three-dimensional convolution operator based on the dimension of the seismic data volume and the calculated seismic wavelet;

[0024] Step S104: Construct a three-dimensional reflection coefficient inversion objective function with lateral constraint terms using the aforementioned lateral second-order difference three-dimensional convolution operator;

[0025] Step S105: Combine the split Bregman algorithm to split the objective function of the three-dimensional reflection coefficient inversion in the frequency domain, and alternately optimize each sub-problem to obtain the three-dimensional reflection coefficient inversion result.

[0026] Specifically, the acquisition of seismic data volumes and P-wave velocity and density logging curves includes:

[0027] Step 1.1: Convert the seismic data volume Corresponding data file and P-wave velocity ,density The data files corresponding to the two logging curves were read into the computer memory respectively.

[0028] Furthermore, the steps of calculating the well logging reflection coefficient using the well logging curve, constructing the wavelet inversion objective function using the wellside traces and well logging reflection coefficients in the seismic data volume, and calculating the seismic wavelet sequence using the damped least squares method include:

[0029] Step 2.1: Calculate the well logging reflection coefficient sequence using well logging impedance. The formula is as follows:

[0030]

[0031] in, , indicating the first i Well logging impedance values ​​at each time point N This represents the number of time sampling points in the well logging sequence.

[0032] Step 2.2: Utilize the aforementioned well logging reflection coefficient sequence Construct Topritz matrices for well-side seismic traces. and column vectors :

[0033] ,

[0034] in, M This represents the number of time sampling points for the seismic trace. The first earthquake tunnel near the well i The amplitude value of each sampling point.

[0035] Step 2.3: Utilize the aforementioned Topritz matrix The vector formed by the well-side seismic trace Constructing sub-wave vectors Inversion objective function:

[0036]

[0037] in, Describing the L2 norm of a vector, This means finding the expression that minimizes the value of the formula within the curly braces. , This is the regularization parameter.

[0038] Step 2.4: The wavelet sequence is calculated using the damped least squares method. The formula is:

[0039]

[0040] Among them, the superscript " T " indicates matrix transpose, "-1" indicates matrix inversion, and I is the identity matrix.

[0041] Furthermore, the three-dimensional convolution operator for the frequency domain wavelet and the transverse second-order difference three-dimensional convolution operator calculated based on the dimension of the seismic data volume include:

[0042] Step 3.1: Constructing the original seismic data volume Four zero-value data volumes with the same scale: , , and The first dimension X represents the direction of the survey line, and the sampling point number is 1~ G The second dimension Y represents the direction of the connecting survey line, and the sampling point numbers are 1~ F The third dimension Z represents the time direction, and the sampling point numbers are 1~ M ;

[0043] Step 3.2: Change The value of an element at a specific position in the middle: , ,in p This refers to the sampling point number corresponding to time zero of the wavelet; change , and The value of an element at a specific position in the middle: , , , , ; , , , ; , , , .

[0044] When the scale of the 3D seismic data volume G = F = M When =10, Figure 2 The sampling point number at time zero of the wavelet is given. p An example of constructing a wavelet three-dimensional convolution operator W from the wavelet vector w when =5; Figure 3 The construction of a transverse second-order difference three-dimensional convolution operator is given. , and Examples of other scales , , and This can be derived from this embodiment without any additional creative effort.

[0045] Step 3.3: Calculation , , and 3D Fourier Transform , , and and 3D seismic data volume 3D Fourier Transform .

[0046] Furthermore, constructing a three-dimensional reflection coefficient inversion objective function with lateral constraints using the aforementioned second-order difference three-dimensional convolution operator includes:

[0047] Step 4.1:

[0048] in, For three-dimensional reflectance, Denotes the Frobenius norm. Describes the norm 1. Represents three-dimensional temporal convolution. This means finding the expression that minimizes the value of the formula within the curly braces. , μ For time-direction sparsity regularization parameters. λ These are the parameters for lateral constraint regularization.

[0049] Furthermore, the step of combining the split Bregman algorithm to split the objective function for the three-dimensional reflection coefficient inversion in the frequency domain, and alternately optimizing each sub-problem to obtain the three-dimensional reflection coefficient inversion result includes:

[0050] Step 5.1: Given the maximum number of iterations K Given an index number k =0, given eight three-dimensional zero-value data volumes with the same scale as the original seismic data: Given a three-dimensional fast Fourier transform operator and .

[0051] Step 5.2: Split the objective function for retrieving the three-dimensional reflection coefficients into the following frequency domain optimization sub-problem:

[0052] ①

[0053] ②

[0054] ③

[0055] ④

[0056] ⑤

[0057] Among them, the symbol " "" indicates the dot product operation between two three-dimensional data volumes.

[0058] Step 5.3: According to the following formula and description , , , and Reflection coefficient in the three-dimensional Fourier domain and each intermediate variable , , and Alternate updates:

[0059] ①

[0060] in, and They are and Conjugate;

[0061] ②

[0062] in, , It is a symbolic function;

[0063] ③

[0064] in,

[0065] ④

[0066] ⑤ .

[0067] Step 5.4: Let ,like Repeat step 5.3 for the alternating update; otherwise, end the update process and calculate the three-dimensional reflection coefficient. .

[0068] Figure 4 This is a schematic diagram of a three-dimensional seismic data volume provided in an embodiment of the present invention. The data contains significant random noise. Figure 5 The three-dimensional reflection coefficient inversion results obtained using the method of this invention are smooth and natural at the top and bottom interfaces of each stratum, without any jitter artifacts. Figure 6 for Figure 5 The profile corresponding to line 100 in the inversion results shown by the arrow indicates that the extension of the thin layer's top and bottom interfaces conforms to geological laws and does not contain any shaking artifacts. Figure 7 for Figure 5 The time slice at 530 milliseconds in the inversion results shown can be clearly traced and interpreted on the plane, as indicated by the arrow. This further illustrates that the method of the present invention can effectively suppress the artifact of reflection coefficient jitter caused by noise.

Claims

1. A three-dimensional reflection coefficient inversion method with a function of suppressing a jitter artifact, characterized by, The method comprises the following steps: obtaining a seismic data volume S and P-wave velocity alpha and density rho logging curves; Using well logs reflection coefficients r i and seismic traces a i Constructing a wavelet inversion objective function and calculating seismic wavelet vector w by damped least squares method; Constructing a frequency domain wavelet 3-D convolution operator of the same dimension as seismic data volume and a lateral second difference 3-D convolution operator and constructing a three-dimensional reflectivity inversion objective function with a lateral constraint term by using the lateral second-order difference three-dimensional convolution operator; combining a split Bregman algorithm to perform frequency domain splitting on the objective function, and alternately optimizing each split sub-problem to obtain a three-dimensional reflectivity inversion result R.

2. The three-dimensional reflectivity inversion method with a function of suppressing a jitter artifact according to claim 1, characterized by, The method for constructing a wavelet inversion objective function by using logging reflectivity and constructing a wavelet from a seismic trace beside a well, and calculating a seismic wavelet sequence by using a damped least square method comprises: Using the well-logging reflection coefficients r i and the seismic trace beside the well to construct the topatriz matrix R J and the column vector S J : wherein M is the number of time sampling points of the seismic trace, a i (i = 1, 2,..., M) is the amplitude value of the i-th sampling point in the seismic trace near the well. The matrix R J and the column vector S J Construct the inversion objective function of the wavelet vector w: where || | | |2denotes the vector two-norm, denotes the w that minimizes the expression in the braces, and η is a regularization parameter. The method for solving formula (2) by using a damped least square method is: Wherein, the superscript "T" represents matrix transposition, "-1" represents matrix inversion, and I is a unit matrix.

3. The three-dimensional reflectivity inversion method with a function of suppressing a jitter artifact according to claim 1, characterized by, The method for constructing a frequency domain wavelet three-dimensional convolution operator and a lateral second-order difference three-dimensional convolution operator with the same dimension as a seismic data volume comprises: four zero-value data volumes W, Ψ1, Ψ2 and Ψ3 with the same scale as the original seismic data volume S are constructed, the first dimension X of these data volumes is a survey line direction, the sampling point number is 1, 2,..., G, the second dimension Y is a tie line direction, the sampling point number is 1, 2,..., F, and the third dimension Z is a time direction, the sampling point number is 1, 2,..., M; the element values are changed as follows: W(1, 1, 1:M-p+1) = W(p:M), W(1, 1, M-p+2:M) = W(1:p-1), wherein p is a sampling point number corresponding to a zero time of a wavelet; the element values are changed as follows: Ψ1(1, 1, 1) =-4, Ψ1(2, 1, 1) = 1, Ψ1(G, 1, 1) = 1, Ψ1(1, 2, 1) = 1, Ψ1(1, F, 1) = 1, Ψ2(2, 1, 1) = 1, Ψ2(G, 1, 1) = 1, Ψ2(1, 2, 1) =-1, Ψ2(1, F, 1) =-1, Ψ3(1, 1, 1) = 2, Ψ3(2, 1, 1) =-2, Ψ3(1, 2, 1) =-2, Ψ3(2, 2, 1) = 2; A three-dimensional Fourier transform is performed on S, W, Ψ1, Ψ2, and Ψ3 to obtain corresponding frequency domain data and 4. The three-dimensional reflectivity inversion method with a function of suppressing a jitter artifact according to claim 1, characterized by, The method for constructing a three-dimensional reflectivity inversion objective function with a lateral constraint term by using the lateral second-order difference three-dimensional convolution operator comprises: an objective function is established: where R is the three-dimensional reflection coefficient, || ||1is the one-norm, F denotes the Frobenius norm, || ||1denotes the one-norm, denotes the three-dimensional time-domain convolution, denotes the R that minimizes the expression in the braces, μ is a time- direction sparsity regularization parameter, and λ is a lateral constraint regularization parameter.

5. The three-dimensional reflectivity inversion method with a function of suppressing a jitter artifact according to claim 1, characterized by, The method for combining a split Bregman algorithm to perform frequency domain splitting on the objective function, and alternately optimizing each split sub-problem to obtain a three-dimensional reflectivity inversion result R comprises: Given a maximum number of iterations K, given an index number k = 0, given eight three-dimensional zero-value data volumes of the same dimensions as the original seismic data: Given a three-dimensional fast Fourier transform operator FFT 3D and the three-dimensional reflectivity inversion objective function is split into the following frequency domain optimization sub-problems: ① ② ③ ④ ⑤ wherein the symbol represents a point multiplication operation of two three-dimensional data volumes; The reflectivity in the three-dimensional Fourier domain is updated alternately according to the following equations and the intermediate variables and are updated alternately. ① wherein and are respectively and conjugates of ② wherein sign is the sign function; ③ wherein ④ ⑤ The iterations are stopped when the maximum number of iterations K is reached, and the final 3D reflectivity inversion results are obtained using the formula The final 3D reflectivity inversion results are obtained.

Images