A method for reconstructing seismic reflection signals

Abstract

The application discloses a seismic reflection signal reconstruction method and relates to the technical field of oil exploration. The method comprises the following steps: inputting a missing seismic signal; constructing a target function; solving the target function by using an adaptive accelerating sparse optimization algorithm; transforming a three-dimensional seismic data body containing 50% missing primary reflection waves and multiple waves to a hyperbolic sparse domain to obtain a three-dimensional sparse constraint hyperbolic Radon forward transformation model; converting the three-dimensional sparse constraint hyperbolic Radon forward transformation model to a seismic data domain to obtain a seismic reflection signal reconstruction result; and reconstructing all data bodies. The fast three-dimensional seismic data reconstruction method adopts a three-dimensional hyperbolic Radon transformation method based on an adaptive accelerating sparse optimization algorithm to reconstruct effective signals, better reconstructs missing data in the three-dimensional seismic data body, and significantly reduces the calculation cost.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil exploration technology, and particularly relates to a seismic reflection signal reconstruction method. BACKGROUND

[0002] The existing seismic acquisition is limited by high-density cost, and the actual data trace interval is much larger than the high-density standard. Moreover, the city buildings and obstacles often cause missing traces, resulting in insufficient imaging resolution. The conventional reflection signal reconstruction is based on sparse representation, low-rank completion or five-dimensional interpolation, which can partially repair the missing traces, but needs to iteratively solve large matrices, and the calculation amount increases in cubic with the data amount, which is difficult to meet the on-site rapid processing demand.

[0003] Therefore, a new method considering accuracy and efficiency is needed, which reduces the calculation complexity to linear or quasi-linear while maintaining the amplitude fidelity, so as to adapt to the low-cost acquisition and real-time decision-making scene. SUMMARY

[0004] To solve the above technical problems, the present application discloses a seismic reflection signal reconstruction method. Compared with the conventional seismic reflection signal reconstruction method, the present application introduces a dynamic weight adjustment mechanism to enhance the signal sparsity constraint, which significantly improves the inversion resolution. In combination with the spectral space momentum transfer technology, the algorithm convergence rate is improved from the sub-linear O(1 / k) of the weighted conjugate gradient method to O(1 / k²), and the Radon forward and inverse transformation calculation amount of a single iteration is reduced to 2 times, which improves the calculation efficiency by more than 3 times.

[0005] To achieve the above purpose, the present application adopts the following technical solutions:

[0006] A seismic reflection signal reconstruction method comprises the following steps:

[0007] Step 1. Input the missing seismic signal;

[0008] Step 2. Construct the objective function, and solve the objective function by using the adaptive accelerated sparse optimization algorithm. The three-dimensional seismic data body containing 50% missing primary reflection wave and multiple wave is transformed to the hyperbolic sparse domain to obtain a three-dimensional sparse constraint hyperbolic Radon forward transformation model;

[0009] Step 3. Convert the three-dimensional sparse constraint hyperbolic Radon forward transformation model obtained in step 2 into the seismic data domain to obtain the seismic reflection signal reconstruction result;

[0010] Step 4. Determine whether all data bodies have been reconstructed, if not, return to step 2; if yes, output the reconstructed three-dimensional seismic data body, and the data reconstruction is completed.

[0011] Optionally, in step 1, for the three-dimensional CMP record obtained by sequential excitation, the sparse constraint random subsampling process is represented as follows:

[0012] (1);

[0013] (2);

[0014] (3);

[0015] where x, y are the longitudinal and lateral coordinates in three-dimensional space, t is time, is the time in the new coordinate system, is the positive transformation, which belongs to the forward order, that is, another order is equal to , , is the positive transformation of the regular and complete high-density three-dimensional CMP record, is the regular and complete high-density three-dimensional CMP record after the positive transformation, is the incomplete three-dimensional CMP record obtained by random subsampling after the positive transformation, is the incomplete three-dimensional CMP record obtained by random subsampling after the positive transformation, N x and N y are the number of lateral and longitudinal lines, respectively, is the regular and complete high-density three-dimensional CMP record after the positive transformation, is the regular and complete high-density three-dimensional CMP record after the positive transformation, is a diagonal matrix for random subsampling. Optionally, step 2 specifically comprises:

[0016] Step 21. The time-domain continuous form of the three-dimensional hyperbolic Lagrangian transformation is expressed as:

[0017] Step 21. The time-domain continuous form of the three-dimensional hyperbolic Lagrangian transformation is expressed as:

[0018] (4);

[0019] where D h (t, x, y) is the value of the input seismic signal at time t coordinate (x, y), wherein x and y are along the longitudinal and lateral line directions, respectively; is the HRT model, represents the Dirac function, is the basis function of HRT; variables , and represent the intercept, slowness along the longitudinal and lateral lines, respectively;

[0020] Step 22. Perform three-dimensional hyperbolic Lagrangian transformation in the frequency domain, and transform formula (4) by using formula (5) to obtain formula (6):

[0021] ​(5);

[0022] (6);

[0023] where, is the intercept in the new coordinate system;

[0024] From equation (4) to equation (6) The forward transformation process is symbolically abbreviated as:

[0025] (7);

[0026] According to equation (5) definition The inverse transformation process is as follows:

[0027] (8);

[0028] Therefore, equation (6) is converted to equation (4) using equation (7) The inverse transformation process is symbolically abbreviated as:

[0029] (9);

[0030] where, is The inverse transformation belongs to the inverse order, that is, another order t is equal to ;

[0031] Step 23. Through the coordinate transformation of the time axis, the three-dimensional hyperbolic Lagrangian transformation is converted into a three-dimensional parabolic Lagrangian transformation in the new coordinate system, and the Fourier transform is used on equation (6) to obtain:

[0032] (10);

[0033] where, i is the imaginary unit, is the frequency domain seismic signal, is the frequency, is the frequency domain HRT model;

[0034] Equation (10) is discretized into a discrete form shown in equation (11):

[0035] (11);

[0036] where, and are the number of seismic traces in the inline and crossline directions, respectively; and are the number of slowness sampling points in the inline and crossline directions, respectively;

[0037] Equation (11) is in the frequency domain For a complex linear system of variables, it can be written in matrix form as follows:

[0038] (12);

[0039] wherein, and are vectors and matrices with frequency as variable; for each frequency, and are the HRT operators in the frequency domain along the NMO and the crossline directions respectively, is the HRT model; is the 3D hyperbolic Radon transform operator in the frequency domain; the matrix and have the forms of equation (13) and equation (14) respectively:

[0040] (13);

[0041] (14);

[0042] Step 24. Obtain the 3D mixed-domain hyperbolic Radon transform according to equation (12) for the purpose of balancing the computational efficiency in the frequency domain and the sparsity in the time domain, as follows:

[0043] (15);

[0044] wherein, and are the forward and inverse Fourier transforms respectively, is the 3D CMP record obtained using the 3D hyperbolic Radon transform;

[0045] Let the regular and complete high-density 3D CMP record after the forward transform can be sparsely expressed in the 3D hyperbolic sparse domain, and combining equation (2) and equation (15) gives:

[0046] (16);

[0047] Step 25. Substitute equation (16) into equation (2) and construct the objective function based on the sparsity constraint theory as follows:

[0048] (17);

[0049] wherein, is the objective function, is a trade-off parameter for adjusting the proportion between the L2 norm and the L1 norm;

[0050] Step 26. Solve equation (17) by using the adaptive accelerating sparse optimization algorithm to obtain the sparse constraint 3D hyperbolic Radon forward transform model where the algorithm iteration formula is as follows:

[0051] (18);

[0052] (19);

[0053] (20);

[0054] (21);

[0055] wherein r represents a residual, k represents an iteration number, is a step length, c is an iteration step length, represents a current solution of the algorithm at the kth iteration, is an auxiliary variable of the algorithm at the kth iteration, used to accelerate the convergence of the algorithm, and is obtained by updating a momentum term; represents a soft threshold function, and the expression is as follows:

[0056] (22);

[0057] wherein is a sign function, is a threshold parameter, and the empirical expression is as follows:

[0058] (23).

[0059] Optionally, in step 3, a three-dimensional hyperbolic Radon inverse transform operator is used with inverse transform to restore complete high-density three-dimensional CMP records The processes are respectively expressed as follows:

[0060] (24);

[0061] (25).

[0062] The present application has the beneficial effects that the sparse constraint seismic reflection signal reconstruction method based on the fast three-dimensional hyperbolic Radon transform is used to solve the sparse expression base in the seismic reflection signal reconstruction method by using an adaptive accelerated sparse optimization algorithm; compared with the conventional reweighted conjugate gradient method, the dynamic weight adjustment mechanism is introduced in the present application to enhance the signal sparsity constraint, and the inversion resolution is significantly improved; in combination with the spectral space momentum transfer technology, the algorithm convergence rate is improved from the sub-linear O(1 / k) of the reweighted conjugate gradient method to O(1 / k²), and the Radon forward and inverse transform calculation amount of a single iteration is reduced to 2 times, the calculation efficiency is improved by more than 3 times, and thus the seismic reflection signal in the three-dimensional CMP record is more accurately reconstructed. Attached Figure Description

[0063] Figure 1 This is a schematic flowchart of a seismic reflection signal reconstruction method according to the present invention;

[0064] Figure 2 This is a 3D seismic data volume containing 50% missing primary and secondary reflections;

[0065] Figure 3 To utilize the three-dimensional hyperbolic sparse domain zero-intercept time slice obtained by this invention;

[0066] Figure 4 To obtain a zero-intercept time slice for a three-dimensional hyperbolic sparse domain using the traditional reweighted conjugate gradient method;

[0067] Figure 5 The reconstructed three-dimensional data volume using this invention;

[0068] Figure 6 This is a 3D data volume reconstructed using the traditional reweighted conjugate gradient method. Detailed Implementation

[0069] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0070] A method for reconstructing seismic reflection signals, such as Figure 1 As shown, it includes the following steps:

[0071] Step 1. Input the missing seismic signal; for the 3D CMP record obtained by sequential excitation, the sparse-constrained random subsampling process is represented as follows:

[0072] (1);

[0073] (2);

[0074] (3);

[0075] where x, y are the longitudinal and lateral coordinates in three-dimensional space, and t is time, is the time in the new coordinate system, is forward transform, is the regular and complete high-density three-dimensional CMP record after forward transform, is the regular and complete high-density three-dimensional CMP record after forward transform, is is the regular and complete high-density three-dimensional CMP record after forward transform, , N x and N y are the number of lateral and longitudinal traces, respectively, is is the incomplete three-dimensional CMP record obtained by random subsampling after forward transform, is is the regular and complete high-density three-dimensional CMP record after forward transform, is a diagonal matrix used for random subsampling.

[0076] Step 2. Construct the objective function, and use the adaptive accelerated sparse optimization algorithm to solve the objective function to transform the three-dimensional seismic data volume containing 50% missing primary and multiple waves to the hyperbolic sparse domain to obtain a three-dimensional sparse constraint hyperbolic Radon forward transform model; specifically including:

[0077] Step 21. The continuous form of the three-dimensional hyperbolic Radon transform in the time domain is expressed as:

[0078] (4);

[0079] where D h (t, x, y) is the value of the input seismic signal at time t coordinate (x, y), wherein x and y are along the longitudinal and lateral directions, respectively; is the HRT model, represents Delta (Delta function), is the basis function of HRT; variables , and represent the intercept, slowness along the longitudinal and lateral directions, respectively;

[0080] Step 22. The three-dimensional hyperbolic Radon transform in the frequency domain has the advantage of high computational efficiency. In order to be able to perform three-dimensional hyperbolic Radon transform in the frequency domain, the radical in equation (4) needs to be eliminated, and equation (4) is transformed using equation (5) to obtain equation (6):

[0081] (5);

[0082] (6);

[0083] In the formula, The intercept in the new coordinate system;

[0084] From equation (4) to equation (6) The forward transformation process is abbreviated as:

[0085] (7);

[0086] Accordingly, according to the definition in equation (5) The inverse transformation process is as follows:

[0087] (8);

[0088] Therefore, equation (6) can be transformed into equation (4) using equation (7). The inverse transformation process is abbreviated as:

[0089] (9);

[0090] In the formula, yes Inverse transform;

[0091] Step 23. By transforming the coordinates of the time axis, the three-dimensional hyperbolic Radon transform is transformed into a three-dimensional parabolic Radon transform in the new coordinate system. Applying a Fourier transform to equation (6) yields:

[0092] (10);

[0093] In the formula, i is the imaginary unit. It is a frequency domain seismic signal. For frequency, It is a frequency domain HRT model;

[0094] Equation (10) is discretized into the discrete form shown in Equation (11):

[0095] (11);

[0096] In the formula, and These are the number of seismic traces along the longitudinal and transverse directions, respectively. and These represent the number of slow sampling points in the longitudinal and transverse survey lines, respectively.

[0097] Equation (11) is based on frequency For a complex linear system with variables, the matrix form can be used as follows:

[0098] (12);

[0099] wherein, and are vectors and matrices with frequency as variable; for each frequency, and are HRT operators in frequency domain along the NMO and crossline directions respectively, is the HRT model; is the 3D hyperbolic Radon transform operator; the matrix and have the forms of equation (13) and equation (14) respectively:

[0100] (13);

[0101] (14);

[0102] Step 24. Obtain the 3D mixed domain hyperbolic Radon transform for the purpose of balancing the computational efficiency in frequency domain and sparsity in time domain according to equation (12) as follows:

[0103] (15);

[0104] wherein, and are forward and inverse Fourier transforms respectively, is the 3D CMP record obtained using the 3D hyperbolic Radon transform;

[0105] Let the regular and complete high-density 3D CMP record after forward transform can be sparsely expressed in 3D hyperbolic sparse domain, and equation (2) and equation (15) are combined to obtain:

[0106] (16);

[0107] Step 25. Substitute equation (16) into equation (2) and construct the objective function based on the sparsity constraint theory as follows:

[0108] (17);

[0109] wherein, is the objective function, is a trade-off parameter for adjusting the proportion between L2 norm and L1 norm, and in the embodiment = 0.5;

[0110] Step 26. Solve equation (17) by using the adaptive accelerating sparse optimization algorithm to obtain the sparse constraint 3D hyperbolic Radon forward transform model , wherein the algorithm iteration formula is as follows:

[0111] (18);

[0112] (19);

[0113] (20);

[0114] (21);

[0115] where r represents residual, k represents iteration number, is a step size, c is an iteration step size, represents the current solution of the algorithm at the kth iteration, is an auxiliary variable of the algorithm at the kth iteration, which is used to accelerate the convergence of the algorithm, and is obtained by updating the momentum term; represents a soft threshold function, and its expression is:

[0116] (22);

[0117] where, is a sign function, is a threshold parameter, and its empirical expression is:

[0118] (23).

[0119] Step 3. Convert the three-dimensional sparse constraint hyperbolic Radon forward transform model obtained in step 2 into a seismic data domain to obtain a seismic reflection signal reconstruction result; use a three-dimensional hyperbolic Radon inverse transform operator and inverse transform to restore complete high-density three-dimensional CMP records , and the process is expressed as follows, respectively:

[0120] (24);

[0121] (25).

[0122] Step 4. Determine whether all data bodies have been reconstructed, if not, return to step 2; if yes, output the reconstructed three-dimensional seismic data body, and the data reconstruction is completed.

[0123] The fast three-dimensional seismic data reconstruction method of the application is applied to a three-dimensional seismic data body, and ideal calculation effect is obtained.

[0124] Figure 2 is a three-dimensional seismic data body with 50% missing of primary reflection wave and multiple wave, and it can be seen that the effective signal is seriously missing. Figure 3 is a zero-intercept time slice obtained by using the three-dimensional seismic data reconstruction optimization algorithm of the application; in order to be used as a comparison, Figure 4The zero-intercept time slice obtained by the traditional reweighted conjugate gradient method is shown, and it can be seen by comparison that the method has higher resolution than the traditional method. Figure 5 The reconstructed three-dimensional data volume obtained by the three-dimensional seismic data reconstruction optimization algorithm of the application is used as a comparison, Figure 6 The reconstructed three-dimensional data volume obtained by the traditional reweighted conjugate gradient method is shown, and it can be seen by comparison that the method has higher resolution than the traditional method.

[0125] Of course, the above description is not a limitation of the application, and the application is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the scope of the application should also be within the scope of the application.

Claims

1. A method of seismic reflection signal reconstruction, characterized by, The method comprises the following steps: Step 1. inputting a missing seismic signal; Step 2. constructing an objective function, solving the objective function by using an adaptive accelerating sparse optimization algorithm, transforming a three-dimensional seismic data body containing 50% missing primary reflection wave and multiple wave to a hyperbolic sparse domain to obtain a three-dimensional sparse constraint hyperbolic Radon forward transformation model; Step 3. converting the three-dimensional sparse constraint hyperbolic Radon forward transformation model obtained in step 2 into a seismic data domain to obtain a seismic reflection signal reconstruction result; Step 4. judging whether all data bodies have been reconstructed, if not, returning to step 2; if yes, outputting the reconstructed three-dimensional seismic data body, and ending the data reconstruction; In step 1, for a three-dimensional CMP record obtained by sequential excitation, the sparse constraint random subsampling process is as follows: （1）； （2）； （3）； where x, y are the longitudinal and lateral coordinates in three-dimensional space, and t is time, is the time in the new coordinate system, is the forward transform, is the regular and complete high-density 3D CMP record after forward transform, is the regular and complete high-density 3D CMP record after forward transform, , N x and N y are the number of seismic traces in the longitudinal and lateral directions, respectively, is the incomplete 3D CMP record after forward transform and random sub-sampling, is the regular and complete high-density 3D CMP record after forward transform, is a diagonal matrix for random sub-sampling; Step 2 specifically comprises: Step 21. a continuous form of three-dimensional hyperbolic Radon transformation in a time domain is expressed as: （4）； where D h (t, x, y) is the value of the input seismic signal at time t and coordinates (x, y), where x and y are along the in-line and cross-line directions, respectively; is the HRT model, denotes the Dirac function, is the basis function of HRT; the variable , and denote the intercept, the slowness along the in-line and cross-line directions, respectively. Step 22. performing three-dimensional hyperbolic Radon transformation in a frequency domain, transforming formula (4) into formula (6) by using formula (5): （5）； （6）； In the formula, is the intercept in the new coordinate system; From equation (4) to equation (6) The forward transform process is symbolically abbreviated as: （7）； According to the definition of formula (5) The inverse transform process is as follows: （8）； Thus, equation (6) is converted to equation (4) using equation (7) The inverse transform process is symbolically abbreviated as: （9）； In the formula, is inverse transform; Step 23. by coordinate transformation on a time axis, the three-dimensional hyperbolic Radon transformation is converted into three-dimensional parabolic Radon transformation in a new coordinate system, and formula (6) is subjected to Fourier transformation to obtain formula (7): （10）； where i is the imaginary unit, is a frequency domain seismic signal, is a frequency, is a frequency domain HRT model; Formula (10) is discretized into a discrete form shown in formula (11): （11）； In the formula, N x N y respectively are the number of seismic traces in the in-line and cross-line directions, respectively are the number of seismic traces in the in-line and cross-line directions, respectively are the number of seismic traces in the in-line and cross-line directions, Equation (11) is a complex linear system in the frequency variable, written in matrix form as follows: （12）； where are vectors and matrices in frequency; for each frequency, are vectors and matrices in frequency; for each frequency, are the HRT operators along the N and T directions in frequency domain, are the HRT operators along the N and T directions in frequency domain, is the HRT model; is the 3D hyperbolic Radon transform operator; the matrix are in the form of equation (13) and equation (14), respectively: are in the form of equation (13) and equation (14), respectively: （13）； （14）； Step 24. according to formula (12), a three-dimensional mixed domain hyperbolic Radon transformation is obtained, which is used for giving consideration to the calculation efficiency in the frequency domain and the sparsity in the time domain, as follows: （15）； wherein and are forward and inverse Fourier transforms, respectively, for obtaining a 3D CMP record using a 3D hyperbolic Radon transform; Let A regular and complete high-density 3D CMP record after positive transformation can be sparsely represented in a 3D hyperbolic sparse domain, combining equation (2) and equation (15) gives: （16）； Step 25. substituting formula (16) into formula (2) and constructing an objective function based on the sparse constraint theory as follows: （17）； In the formula, is a target function, is a trade-off parameter for adjusting the proportion between the L2 norm and the L1 norm; Step 26. Solve equation (17) by the adaptive accelerating sparse optimization algorithm to obtain the sparse constrained three-dimensional hyperbolic Rayleigh normal transform model where the iterative formula of the algorithm is as follows: （18）； （19）； （20）； （21）； where r denotes the residual, k denotes the iteration number, is the step size, c is the iteration step size, denotes the current solution of the algorithm at the kth iteration, is an auxiliary variable of the algorithm at the kth iteration, denotes the soft threshold function, which is expressed as: （22）； wherein is a symbol function, is a threshold parameter, whose empirical expression is: （23）。 2. A method of seismic reflection signal reconstruction as claimed in claim 1, characterized in that, In step 3, the three-dimensional inverse Radon transform operator is used with inverse transform to recover the full high-density 3D CMP record The process is expressed as follows, respectively: （24）； （25）。

Images