Seismic data high-resolution processing method

Abstract

The invention provides a seismic data high-resolution processing method for solving the problem that existing seismic data is insufficient in resolution for underground structures in complex regions. According to the method, a sparse seismic reflection coefficient is obtained by solving an optimization problem of sparse constraint to carry out frequency expansion on seismic data, and the method comprises the following steps: obtaining a low-frequency wavelet function according to an average spectrum of the seismic data, and obtaining an initial seismic reflection coefficient through channel-by-channel deconvolution seismic data and the low-frequency wavelet function; and iteratively extracting numerical reconstruction sparse seismic reflection coefficients with dominant values in the initial reflection coefficients by using the sparsity of the reflection coefficients, and carrying out convolution on the sparse seismic reflection coefficients and the high-frequency seismic wavelets to obtain seismic data after frequency expansion. According to the method, a seismic section with higher resolution can be provided through signal processing, and description and interpretation of a thin-layer oil and gas reservoir structure are facilitated.

Description

Technical field:

[0001] The present invention relates to the technical field of seismic surveying, and in particular to a method for high-resolution processing of seismic data. Background technology:

[0002] Seismic data is an important surveying technology that uses seismic waves reflected underground to image underground structures. It has application significance in the fields of petroleum geophysical prospecting and engineering geological survey. With the deepening of underground exploration, engineering practice has put forward higher requirements for fine exploration of thin-layer structures. Affected by the source excitation technology and stratum propagation attenuation, the high-frequency components in seismic data are often relatively limited, resulting in the corresponding seismic data resolution being insufficient to meet the needs.

[0003] In order to fully explore the ability of seismic data to resolve underground structures, high-resolution processing technology has been well developed, including but not limited to deconvolution, inverse Q filtering, Hilbert transform, wavelet analysis, spectral inversion and other frequency extension techniques. However, these methods usually do not fully utilize the sparsity of seismic reflection coefficients and are still limited in amplitude-preserving and high-resolution processing of complex structures.

[0004] For fine exploration, using high-resolution processing technology that preserves amplitude and fidelity to improve the resolution of seismic data for underground structures in complex areas is a key issue that needs to be solved urgently. To this end, the present invention proposes a high-resolution seismic processing method that uses sparse constrained inversion to accurately extract seismic reflection coefficients and broadens the frequency band of seismic data by preserving amplitude and fidelity through the full-band nature of the reflection coefficients. Summary of the invention:

[0005] The present invention aims to provide a high-resolution seismic data processing method to solve the problem that the existing seismic data in the background technology have insufficient resolution for the underground structure in complex areas. The high-resolution seismic data processing method can robustly obtain accurate sparse seismic reflection coefficients, and based on this, combined with continuous wavelet transform high-frequency wavelets, obtain seismic profile data with higher resolution while preserving amplitude and fidelity.

[0006] The present invention can solve the problem through the following technical solution: the high-resolution processing method of seismic data comprises the following steps:

[0007] S1. Determine the average frequency spectrum range of seismic data;

[0008] S2, designing a low-frequency wavelet w for extracting the reflection coefficient based on the determined average spectrum range;

[0009] S3, based on the selected low-frequency wavelet, process the seismic data one by one; by minimizing the objective function, using the sparse constrained pulse deconvolution inversion method, obtain the optimal sparse seismic reflection coefficient corresponding to each seismic data;

[0010] S4, obtaining a high-frequency wavelet h according to the obtained optimal sparse seismic reflection coefficient r;

[0011] S5. Based on the obtained optimal sparse seismic reflection coefficient r and high-frequency wavelet h, high-resolution seismic data after frequency extension is obtained.

[0012] Preferably, the method for determining the average frequency spectrum range of the seismic data in step S1 is: determining the average frequency spectrum range of the seismic data by performing a one-dimensional fast Fourier transform on the seismic data along the time / depth direction.

[0013] Preferably, in step S2, based on the determined average spectrum range, the method for designing a low-frequency sub-wave w for extracting the reflection coefficient is:

[0014] Based on the determined average spectrum range of seismic data, the seismic reflection coefficient r obtained by joint logging is 测井 , through the seismic data corresponding to the well 井震 , using the seismic data convolution model, the low-frequency wavelet of the reflection coefficient is extracted.

[0015] Preferably, the convolution model is:

[0016] d 井震 =R 测井 w;

[0017] The expression of the low-frequency wavelet for extracting the reflection coefficient is:

[0018] w=(R 测井 T R 测井 +∈I) -1 R 测井 T d 井震 ;

[0019] Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, R 测井 w means r 测井 The convolution result with w, ∈ is the regularization factor, I is the convolution result with R T R is the identity matrix of the same size, and T is the matrix rank conversion operator.

[0020] Preferably, the specific method of obtaining the optimal sparse seismic reflection coefficient by minimizing the objective function in S3 includes:

[0021] a. Define the target function J(r) about the sparse seismic reflection coefficient r; define the continuous reflection coefficient cache variable r tmp = r, initialize the optimal sparse reflection coefficient r to zero; define the data cache variable d tmp=d; where d is single-channel seismic data;

[0022] The expression of the objective function J(r) is:

[0023] J(r)=||Wr-d|| 2 +λ||r|| 0 ;

[0024] Where: W is the Toeplitz matrix corresponding to the low-frequency wavelet w, and Wr represents the convolution result of w and r;

[0025] b. Obtain the initial continuous reflection coefficient through deconvolution calculation

[0026] r init =(W T W) -1 W T d tmp ;

[0027] Where: T is the matrix rank operator; calculate the initial objective function value J init =||Wr init -d|| 2 +λ||r init || 0 ;

[0028] c. Make the sparse reflection coefficient r tmp The value of is zero, and then r tmp (maxloc(r init ))=maxval(r init );Get the updated sparse reflection coefficient r tmp ;

[0029] Among them: maxloc is the position of the maximum value of the sequence, maxval is the maximum value of the sequence;

[0030] d. Through convolution Wr tmp Get synthetic data f and calculate the normalization coefficient

[0031] e. Update the optimal sparse reflection coefficient r = r + μr tmp ;

[0032] f. Determine the current objective function J value based on the current sparse reflection coefficient r, and compare it with the initial objective function J init The values ​​are compared to determine whether the iteration has converged; if it has converged, it is terminated to obtain the optimal sparse seismic reflection coefficient r; if it has not converged, d is updated tmp , so that d tmp =d tmp -μWr tmp , repeat bf until convergence.

[0033] Preferably, the method for judging iterative convergence in step f is: the current objective function value J is less than the initial objective function value J init When a certain proportion of is reached, the iteration is judged to have converged.

[0034] Preferably, the ratio is a specified value in the range of 0.1%-10%.

[0035] Preferably, the method of obtaining the high-frequency wavelet h according to the obtained optimal sparse seismic reflection coefficient r in step S4 comprises the following steps:

[0036] 1.1. Obtain seismic profile observation data d;

[0037] 1.2. Based on the optimal sparse seismic reflection coefficient r and observation data d, based on the well logging data constraints, according to the convolution model d = Rw, the low-frequency sub-wave w = (R T R+∈I) -1 R T d;

[0038] Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, and Rw represents the convolution result of r and w;

[0039] 1.3. Perform Morlet wavelet-based transformation on the extracted low-frequency wavelet w Continuous wavelet transform Get its time-frequency spectrum response Xw(τ,a);

[0040] Where: * is the conjugate operator, t is time, τ is the time axis of the time-frequency spectrum, a is the frequency axis of the time-frequency spectrum, ω ψ Association, relationship is

[0041] 1.4. In the time-frequency spectrum, the wavelet coefficients are supplemented according to the required high-frequency wavelet main frequency and frequency band to obtain Using inverse continuous wavelet transform Get the high-frequency sub-wave h;

[0042] Where: ψ(ω) is the FFT transform of ψ(t).

[0043] Preferably, the method for obtaining high-resolution seismic data after spectrum extension based on the obtained optimal sparse seismic reflection coefficient r and high-frequency wavelet h is:

[0044] The obtained optimal sparse seismic reflection coefficient r is convolved with the high-frequency wavelet h to obtain high-resolution seismic data after frequency extension;

[0045] The expression of the high-resolution seismic data after frequency extension is:

[0046] d hf =Hr,

[0047] Where H is the Toeplitz matrix corresponding to the high-frequency wavelet h, and Hr represents the convolution result of the high-frequency wavelets h and r.

[0048] Preferably, the method for verifying the high-resolution seismic data after frequency extension is:

[0049] The underground structure information of a large area obtained through ground observations in the seismic profile after frequency extension needs to be compared with the well logging data in a small area before frequency extension after being improved by signal processing methods, in order to analyze the improvement of the local thin-layer structure characterization and conduct quality control on the well-to-seismic ratio through well logging information to confirm whether the high-resolution processing is reasonable. If so, the high-resolution processing of the seismic data is confirmed to be effective frequency extension.

[0050] Compared with the above background technology, the present invention has the following beneficial effects:

[0051] When analyzing underground structures with seismic data, insufficient resolution can lead to incorrect interpretation of thin target layers, which can cause difficulties in subsequent development. However, conventional data acquisition modes are unlikely to provide high-frequency seismic signals with sufficient resolution to identify thin layers, and conventional high-resolution seismic processing technology is unlikely to meet the needs of actual development and production.

[0052] Based on sparse constrained deconvolution, this paper proposes a seismic high-resolution frequency extension processing technology, which can robustly obtain accurate sparse seismic reflection coefficients, and based on this combined with continuous wavelet transform high-frequency wavelets, it can effectively improve the main frequency and bandwidth of seismic data, and obtain seismic profile data with higher resolution while maintaining amplitude fidelity. This method can effectively improve the accuracy of using seismic profiles to finely depict underground target structures, enhance the ability to interpret and identify underground thin layers, and thus assist the process of oil and gas resource reservoir exploration and prediction. Description of the drawings:

[0053] Attached Figure 1 is a flow chart of the present invention;

[0054] Attached Figure 2 is a flow chart of optimal sparse reflection coefficient extraction according to an embodiment of the present invention;

[0055] Attached Figure 3 is a sparse reflection coefficient map extracted by an embodiment of the present invention;

[0056] Attached Figure 4 is a seismic data profile after high-frequency wavelet replacement in an embodiment of the present invention;

[0057] Attached Figure 5It is a comparison diagram of seismic profile details before and after high-resolution processing according to an embodiment of the present invention;

[0058] Attached Figure 6 It is a comparison diagram of well-seismic matching before and after high-resolution processing in an embodiment of the present invention. Specific implementation method:

[0059] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0060] like Figure 1 As shown, a high-resolution processing method for seismic data comprises the following steps: S1, determining an average frequency spectrum range of seismic data;

[0061] The average spectrum range of the seismic data is determined by performing a one-dimensional fast Fourier transform on the seismic data along the time / depth direction.

[0062] S2. Based on the determined average spectrum range, a low-frequency sub-wave w is designed for extracting the reflection coefficient; the specific method is:

[0063] Based on the determined average spectrum range of seismic data, the seismic reflection coefficient r obtained by joint logging is 测井 , through the seismic data corresponding to the well 井震 , using the seismic data convolution model, the low-frequency wavelet of the reflection coefficient is extracted.

[0064] The convolution model is:

[0065] d 井震 =R 测井 w;

[0066] The expression of the low-frequency wavelet for extracting the reflection coefficient is:

[0067] w=(R 测井 T R 测井 +∈I) -1 R 测井 T d 井震 ;

[0068] Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, R 测井 w means r 测井 The convolution result with w, ∈ is the regularization factor, I is the convolution result with R T R is the identity matrix of the same size, and T is the matrix rank conversion operator.

[0069] S3. Based on the selected low-frequency wavelet, the seismic data are processed one by one; by minimizing the objective function, the sparse constrained pulse deconvolution inversion method is used to obtain the optimal sparse seismic reflection coefficient corresponding to each seismic data; the specific method includes:

[0070] a. Define the target function J(r) about the sparse seismic reflection coefficient r; define the continuous reflection coefficient cache variable r tmp = r, initialize the optimal sparse reflection coefficient r to zero; define the data cache variable d tmp =d; where d is single-channel seismic data;

[0071] The expression of the objective function J(r) is:

[0072] J(r)=||Wr-d|| 2 +λ||r|| 0 ;

[0073] Where: W is the Toeplitz matrix corresponding to the low-frequency wavelet w, and Wr represents the convolution result of w and r;

[0074] b. Through deconvolution calculation, the initial continuous reflection coefficient r is obtained init =(W T W) -1 W T d tmp ;

[0075] Where: T is the matrix rank operator; calculate the initial objective function value J init =||Wr init -d|| 2 +λ||r init || 0 ;

[0076] c. Make the sparse reflection coefficient r tmp The value of is zero, and then r tmp (maxloc(r init ))=maxval(r init );Get the updated sparse reflection coefficient r tmp ;

[0077] Among them: maxloc is the position of the maximum value of the sequence, maxval is the maximum value of the sequence;

[0078] d. Through convolution Wr tmp Get synthetic data f and calculate the normalization coefficient

[0079] e. Update the optimal sparse reflection coefficient r = r + μr tmp ;

[0080] f. Determine the current objective function J value based on the current sparse reflection coefficient r, and compare it with the initial objective function J init The values ​​are compared to determine whether the iteration has converged; if it has converged, it is terminated to obtain the optimal sparse seismic reflection coefficient r; if it has not converged, d is updated tmp , so that d tmp =d tmp -μWr tmp , repeat bf until convergence.

[0081] The method for judging iterative convergence is: the current objective function J value is less than the initial objective function value J init When a certain proportion of , the iteration is judged to have converged. The proportion is a specified value between 0.1% and 10%.

[0082] S4. Obtaining a high-frequency wavelet according to the obtained optimal sparse seismic reflection coefficient; the specific method includes:

[0083] 1.1. Obtain seismic profile observation data d;

[0084] 1.2. Based on the optimal sparse seismic reflection coefficient r and observation data d, based on the well logging data constraints, according to the convolution model d = Rw, the low-frequency sub-wave w = (R T R+∈I) -1 R T d;

[0085] Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, and Rw represents the convolution result of r and w;

[0086] 1.3. Perform Morlet wavelet-based transformation on the extracted low-frequency wavelet w Continuous wavelet transform Get its time-frequency spectrum response Xw(τ,a);

[0087] Where: * is the conjugate operator, t is time, τ is the time axis of the time-frequency spectrum, a is the frequency axis of the time-frequency spectrum, ω ψ Association, relationship is

[0088] 1.4. In the time-frequency spectrum, the wavelet coefficients are supplemented according to the required high-frequency wavelet main frequency and frequency band to obtain Using inverse continuous wavelet transform Get the high-frequency sub-wave h;

[0089] Where: ψ(ω) is the FFT transform of ψ(t).

[0090] S5. Based on the obtained optimal sparse seismic reflection coefficient and high-frequency wavelet, high-resolution seismic data after spectrum extension is obtained. The specific method is as follows.

[0091] The obtained optimal sparse seismic reflection coefficient r is convolved with the high-frequency wavelet h to obtain high-resolution seismic data after frequency extension;

[0092] The expression of the high-resolution seismic data after frequency extension is:

[0093] d hf =Hr;

[0094] Where: H is the Toeplitz matrix corresponding to the high-frequency wavelet h, and Hr represents the convolution result of the high-frequency wavelets h and r.

[0095] The high-resolution seismic data after frequency extension needs to be verified, and the specific verification method is as follows:

[0096] The underground structure information of a large area obtained through ground observations in the seismic profile after frequency extension needs to be compared with the well logging data in a small area before frequency extension after being improved by signal processing methods, in order to analyze the improvement of the local thin-layer structure characterization and conduct quality control on the well-to-seismic ratio through well logging information to confirm whether the high-resolution processing is reasonable. If so, the high-resolution processing of the seismic data is confirmed to be effective frequency extension.

[0097] Example 1

[0098] The thin dolomite reservoirs in the Lower Permian in Sichuan and Chongqing are thin, the seismic resolution is low, and the response characteristics of the thin reservoirs are unclear. This invention proposes a specific implementation example to improve the resolution of thin reservoirs in central Sichuan. Note that this example is only used to illustrate the application method of the present invention and is not used to limit specific application scenarios. The specific steps of this case are as follows:

[0099] 1. Perform one-dimensional fast Fourier transform on the seismic data along the time / depth direction to determine the average spectrum range of the seismic data;

[0100] 2. Based on the average spectrum range of seismic data, the seismic reflection coefficient r obtained by joint logging 测井 , through the seismic data corresponding to the well 井震 , using the convolution model d 井震 =R 测井 w, extract the low-frequency wavelet of reflection coefficient w = (R 测井 T R 测井 +∈I) - 1 R 测井 T d 井震; Where R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, R 测井 w means r 测井 The convolution result with w, ∈ is the regularization factor, I is the convolution result with R T R is the same size as the identity matrix, and T is the matrix rank conversion operator;

[0101] 3. Based on the obtained low-frequency wavelet, process the seismic data one by one; by minimizing the objective function and using the sparse constrained pulse deconvolution inversion method, obtain the sparse seismic reflection coefficient corresponding to each seismic data:

[0102] Assume that the single-channel seismic data is d, the corresponding reflection coefficient is r, and the low-frequency wavelet sequence is w. Expressed in matrix form, the size of these three variables is N t ×1, where N t is the number of sampling points of seismic single-channel data. It can be constructed that the size of w is N t ×N t The Toeplitz matrix W of , expresses the convolution of w and r through Wr.

[0103] The expression of the subsequent objective function for the sparse seismic reflection coefficient r by minimization is:

[0104] r = argmin r J(r)=||Wr-d|| 2 +λ||r|| 0

[0105] Among them, argmin is the subsequent objective function J to be minimized, ||·|| p is the LP norm, and λ is the weight of the hyperparameter balanced sparsification regularization term in the optimization. If the L1 norm approximation is used, the gradient of the objective function with respect to the reflection coefficient can be obtained as:

[0106]

[0107] Where sgn represents the symbolic operation. Based on the objective function J(r) and the gradient The optimal sparse seismic reflection coefficient r can be solved iteratively. To maximize the sparsity of the reflection coefficient, Figure 2 As shown, the specific implementation process of the present invention is as follows:

[0108] aDefine the objective function J(r) = ||Wr-d|| about the sparse seismic reflection coefficient r 2 +λ||r|| 0 , where W is the Toeplitz matrix corresponding to the low-frequency wavelet w, and Wr represents the convolution result of w and r; define the continuous reflection coefficient cache variable r tmp= r, initialize the optimal sparse reflection coefficient r to zero; define the data cache variable d tmp =d;

[0109] b. Through deconvolution calculation, the initial continuous reflection coefficient r is obtained init =(W T W) -1 W T d tmp , where T is the matrix rank operator; calculate the initial objective function value J init =||Wr init -d|| 2 +λ||r init || 0 ;

[0110] c. Define the reflection coefficient r with an initial value of zero tmp , then make r tmp (maxloc(r init ))=maxval(r init ), where maxloc is the position of the maximum value of the sequence, and maxval is the maximum value of the sequence; the initial sparse reflection coefficient r is obtained tmp ;

[0111] d. Through convolution Wr tmp Get the synthetic data f and calculate the normalization coefficient

[0112] e. Update the optimal sparse reflection coefficient r = r + μr tmp ;

[0113] f. Determine whether the iteration converges: Determine the current objective function J value based on the current sparse reflection coefficient r and compare it with the initial objective function J init If the current objective function J value is less than a certain proportion (for example, 1%) of the initial objective function J value, the iteration is judged to be converged. If it converges, it is terminated to obtain the optimal sparse seismic reflection coefficient r; if it does not converge, update d tmp , so that d tmp =d tmp -μWr tmp , repeat bf until convergence. The extracted sparse reflection coefficients are shown in Figure 3 .

[0114] 4. According to the obtained optimal sparse seismic reflection coefficient r and observation data d, based on the well logging data constraints, the shaping regularization method is used to extract the initial wavelet, and then the frequency band of the initial wavelet is expanded by the continuous wavelet transform method. Here, Morlet wavelet replacement is used to obtain the high-frequency wavelet h and its corresponding Toeplitz matrix H; the specific method is as follows:

[0115] Obtain seismic profile observation data d;

[0116] According to the optimal sparse seismic reflection coefficient r and observation data d, based on the logging data constraints, according to the convolution model d = Rw, the low-frequency sub-wave w = (R T R+∈I) -1 R T d; where R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, and Rw represents the convolution result of r and w;

[0117] The extracted low-frequency wavelet w is subjected to Morlet wavelet-based Continuous wavelet transform Get its time-frequency spectrum response Xw(τ,a); where * is the conjugate operator, τ is the time axis of the time-frequency spectrum, and a is the frequency axis ω of the time-frequency spectrum. ψ Association, relationship is In the time-frequency spectrum, the wavelet coefficients are supplemented according to the required high-frequency wavelet main frequency and frequency band to obtain Using inverse continuous wavelet transform Get the high-frequency sub-wave h.

[0118] 5. Finally, the optimal sparse seismic reflection coefficient r is convolved with the high-frequency wavelet h to obtain the high-resolution seismic data d after frequency extension. hf =hr, such as Figure 4 shown.

[0119] 6. Comparison of seismic data before and after the frequency extension of the present invention to improve resolution shows that the seismic frequency band is increased from 10 Hz-45 Hz before the frequency extension to 5 Hz-75 Hz, and the main frequency is increased from 30 Hz to 45 Hz; the seismic resolution is significantly improved, and the thin layer characteristics inside the reservoir are clear, such as Figure 5 On the other hand, from the comparison of synthetic record calibration, the correlation coefficient before frequency extension is 0.65 and after frequency extension is 0.74. The thin layer of dolomite at 1800ms is obvious and the well-seismic matching relationship is good, confirming that the high-resolution processing of seismic data is effective frequency extension, as shown in Figure 2. Figure 6 shown.

Claims

1. A method for high-resolution processing of seismic data, characterized in that: The following steps are involved: S1. Determine the average frequency spectrum range of seismic data; S2, designing a low-frequency wavelet w for extracting the reflection coefficient based on the determined average spectrum range; S3, based on the selected low-frequency wavelet, process the seismic data one by one; by minimizing the objective function, using the sparse constrained pulse deconvolution inversion method, obtain the optimal sparse seismic reflection coefficient corresponding to each seismic data; S4, obtaining a high-frequency wavelet according to the obtained optimal sparse seismic reflection coefficient; S5. Based on the obtained optimal sparse seismic reflection coefficient and high-frequency wavelet, high-resolution seismic data after spectrum extension is obtained.

2. The method for high-resolution processing of seismic data according to claim 1, characterized in that: The method for determining the average frequency spectrum range of the seismic data in step S1 is: determining the average frequency spectrum range of the seismic data by performing a one-dimensional fast Fourier transform on the seismic data along the time / depth direction.

3. The method for high-resolution processing of seismic data according to claim 1, characterized in that: The method of step S2, designing a low-frequency sub-wave w for extracting the reflection coefficient based on the determined average spectrum range, is: Based on the determined average spectrum range of seismic data, the seismic reflection coefficient r obtained by joint logging is 测井 , through the seismic data corresponding to the well 井震 , using the seismic data convolution model, the low-frequency wavelet of the reflection coefficient is extracted.

4. The method for high-resolution processing of seismic data according to claim 3, characterized in that: The convolution model is: d 井震 =R 测井 w; The expression of the low-frequency wavelet for extracting the reflection coefficient is: w=(R 测井 T R 测井 +∈I) -1 R 测井 T d 井震 ; Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, R 测井 w means r 测井 The convolution result with w, ∈ is the regularization factor, I is the convolution result with R T R is the identity matrix of the same size, and T is the matrix rank conversion operator.

5. The method for high-resolution processing of seismic data according to claim 4, characterized in that: The specific method of obtaining the optimal sparse seismic reflection coefficient by minimizing the objective function in S3 includes: a. Define the target function J(r) about the sparse seismic reflection coefficient r; define the continuous reflection coefficient cache variable r tmp = r, initialize the optimal sparse reflection coefficient r to zero; define the data cache variable d tmp =d; where d is single-channel seismic data; The expression of the objective function J(r) is: J(r)=||Wr-d|| 2 +λ||r|| 0 ; Where: W is the Toeplitz matrix corresponding to the low-frequency wavelet w, and Wr represents the convolution result of w and r; b. Through deconvolution calculation, the initial continuous reflection coefficient r is obtained init =(W T W) -1 W T d tmp ; Where: T is the matrix rank operator; calculate the initial objective function value J init =||Wr init -d|| 2 +λ||r init || 0 ; c. Make the sparse reflection coefficient r tmp The value of is zero, and then r tmp (maxloc(r init ))=maxval(r init );Get the updated sparse reflection coefficient r tmp ; Among them: maxloc is the position of the maximum value of the sequence, maxval is the maximum value of the sequence; d. Through convolution Wr tmp Get synthetic data f and calculate the normalization coefficient e. Update the optimal sparse reflection coefficient r = r + μr tmp ; f. Determine the current objective function J value based on the current sparse reflection coefficient r, and compare it with the initial objective function J init The value is compared to determine whether the iteration converges; if it converges, it is terminated to obtain the optimal sparse seismic reflection coefficient r; if it does not converge, d is updated tmp , so that d tmp =d tmp -μWr tmp , repeat bf until convergence.

6. The method for high-resolution processing of seismic data according to claim 5, characterized in that: The method for judging iterative convergence in step f is: the current objective function value J is less than the initial objective function value J init When a certain proportion of is reached, the iteration is judged to have converged.

7. The method for high-resolution processing of seismic data according to claim 6, characterized in that: The current objective function value J is less than the initial objective function value J init The ratio is 0.1%-10% of the specified value.

8. The method for high-resolution processing of seismic data according to claim 6, characterized in that: The method of obtaining the high-frequency wavelet h according to the obtained optimal sparse seismic reflection coefficient r in step S4 comprises the following steps: 1.

1. Obtain seismic profile observation data d; 1.

2. Based on the optimal sparse seismic reflection coefficient r and observation data d, based on the well logging data constraints, according to the convolution model d = Rw, the low-frequency sub-wave w = (R T R+∈I) -1 R T d; Where: R is the Toeplitz matrix corresponding to the sparse seismic reflection coefficient r, and Rw represents the convolution result of r and w; 1.

3. Perform Morlet wavelet-based transformation on the extracted low-frequency wavelet w Continuous wavelet transform Get its time-frequency spectrum response Xw(τ,a); Where: * is the conjugate operator, t is time, τ is the time axis of the time-frequency spectrum, a is the frequency axis of the time-frequency spectrum, ω ψ Association, relationship is 1.

4. In the time-frequency spectrum, the wavelet coefficients are supplemented according to the required high-frequency wavelet main frequency and frequency band to obtain Using inverse continuous wavelet transform Get the high-frequency sub-wave h; Where: ψ(ω) is the FFT transform of ψ(t).

9. The method for high-resolution processing of seismic data according to claim 8, characterized in that: The method for obtaining high-resolution seismic data after spectrum extension based on the obtained optimal sparse seismic reflection coefficient r and high-frequency wavelet h is: The obtained optimal sparse seismic reflection coefficient r is convolved with the high-frequency wavelet h to obtain high-resolution seismic data after frequency extension; The expression of the high-resolution seismic data after frequency extension is: the hf = Mr.; Where: H is the Toeplitz matrix corresponding to the high-frequency wavelet h, and Hr represents the convolution result of the high-frequency wavelets h and r.

10. The method for high-resolution processing of seismic data according to claim 1, characterized in that: The method for verifying the high-resolution seismic data after frequency extension is as follows: The underground structure information of a large area obtained through ground observations in the seismic profile after frequency extension needs to be compared with the well logging data in a small area before frequency extension after being improved by signal processing methods, in order to analyze the improvement of the local thin-layer structure characterization and conduct quality control on the well-to-seismic ratio through well logging information to confirm whether the high-resolution processing is reasonable. If so, the high-resolution processing of the seismic data is confirmed to be effective frequency extension.

Images