A prestack Q-value attenuation compensation method and device based on L2-norm regularization constraint

Abstract

The present invention relates to the technical field of oil and gas exploration and prestack seismic anti-Q compensation, and discloses a prestack Q-value attenuation compensation method and device based on L2 norm regularization constraint, aiming to solve the problem of signal attenuation caused by Q-value attenuation in prestack seismic data. The main solutions include preprocessing prestack seismic data and obtaining prestack Q-value and incident angle information; according to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, using the obtained prestack Q-value and incident angle information to construct a prestack Q-value filtering method related to offset, and establishing a prestack attenuation model related to offset; introducing regularization to construct an inversion objective function, using the conjugate gradient algorithm to iteratively update the seismic record in the initial frequency domain for the objective function, and determining whether the updated seismic record in the frequency domain meets the iteration termination condition. If so, output the seismic record. Otherwise, iteratively update the updated seismic record in the frequency domain again until the iteration termination condition is met.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas exploration and prestack seismic inversion, and discloses a prestack Q - value attenuation compensation method and device based on L2 - norm regularization constraint. Background Art

[0002] During the propagation of seismic waves in the subsurface, due to the visco - elastic properties of the subsurface medium, both its amplitude and phase will be affected by attenuation, thereby reducing the resolution of seismic records. As an effective compensation means, inverse Q - filtering plays a crucial role in improving the resolution of seismic records. Currently, inverse Q - filtering methods are mainly applied to post - stack seismic records. The inverse Q - filtering processing of prestack seismic records is beneficial to improving the resolution and fidelity of seismic records to meet the requirements of high - resolution and high - fidelity seismic records in current seismic exploration.

[0003] When compensating prestack seismic records using inverse Q - filtering methods, the influence of offset needs to be considered. The inverse Q - filtering method of amplitude gain limitation and stability factor proposed by Wang can, to a certain extent, stably compensate prestack seismic records, but the compensation is insufficient in the case of deep layers and low Q - values. To sum up, the existing prestack Q - value attenuation compensation method has a large error in compensating seismic records in the case of deep layers and low Q - values. Summary of the Invention

[0004] The purpose of the present invention is to solve the problem of signal loss caused by Q - value attenuation in prestack seismic records. By constructing a Q - value filtering formula related to offset and a regularization inversion method, high - precision and stable compensation of seismic records is achieved.

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

[0006] The present invention provides a prestack Q - value attenuation compensation method based on L2 - norm regularization constraint, including the following steps:

[0007] Step 1: Pre - process the prestack seismic records and obtain prestack Q - value and incident angle information;

[0008] Step 2: According to the frequency - domain wave - field continuation formula of plane waves in visco - elastic media, use the obtained prestack Q - value and incident angle information to construct a prestack Q - value filtering formula related to offset, and establish a prestack attenuation model related to offset;

[0009] Step 3: On the basis of the prestack Q - value filtering constructed in Step 2, introduce regularization to construct an inversion objective function, and use the conjugate gradient algorithm to iteratively update the seismic records in the initial frequency domain;

[0010] Step 4: Determine whether the updated seismic record in the frequency domain meets the iterative termination condition. If so, output the updated seismic record in the frequency domain as the final inversion result. Otherwise, use the updated seismic record in the frequency domain as the new initial model, return to Step 3, and iterate and update again until the iterative termination condition is met.

[0011] In the above solution, Step 1 includes the following steps:

[0012] Step 1.1: Using the relationship between the prestack seismic record time, offset, and Q, extract the Q value of the prestack seismic record using the prestack Q value extraction method;

[0013] Step 1.2: Obtain the incident angle θ of each layer in the prestack seismic record by using the ray tracing technique i 。

[0014] In the above solution, Step 2 includes the following steps:

[0015] Step 2.1: Establish a propagation model of seismic waves in the medium according to the frequency-domain wavefield continuation formula of plane waves in a viscoelastic medium. The frequency-domain wavefield continuation formula of plane waves is expressed as:

[0016]

[0017] where B(T,ω) represents the wavefield with a propagation time of T and an angular frequency of ω, ΔT i is the propagation time increment of seismic waves in the i-th layer, k(ω) is the frequency-dependent wave number, v r is the velocity at the reference frequency, the variable j represents the imaginary unit, and e is the base of the natural logarithm;

[0018] Step 2.2: Using the obtained prestack Q value and incident angle information, calculate the propagation time increment ΔT of seismic waves in the i-th layer according to Snell's law i , and substitute ΔT i into the frequency-domain wavefield continuation formula to obtain wavefield continuation:

[0019]

[0020] ΔT i is expressed as:

[0021]

[0022] where the ΔT i,0 is the two-way propagation time increment of seismic waves vertically propagating in the i-th layer, θ i is the incident angle of each layer of the medium, and h i represents the thickness of the i-th path or layer;

[0023] Step 2.3: Then, according to the modified Kolsky model, the attenuation and velocity dispersion characteristics of seismic waves are introduced into the wavefield extrapolation formula. The modified Kolsky model describes the viscoelastic properties of the medium by introducing the quality factor Q and the tuning frequency ω0, so as to more accurately simulate the propagation process of seismic waves. At this time, the wavefield extrapolation is expressed as:

[0024]

[0025] where ω0 is the tuning frequency, and Q i represents the quality factor of the i-th layer of the medium, i = 1, 2, 3…n;

[0026] Step 2.4: Using the time-shift property of the Fourier transform, the wavefield extrapolation formula is converted into pre-stack Q-value filtering related to the offset. Specifically:

[0027] Using the time-shift property of the Fourier transform, we have:

[0028]

[0029] Substituting Equation (5) into Equation (4), we obtain the following wavefield extrapolation formula:

[0030]

[0031] When seismic waves propagate through n layers of the medium, the wavefield extrapolation formula at this time becomes:

[0032]

[0033] Replace the propagation time △T with the sampling interval dt i，0 , and discretize the frequency ω into ω m , Equation (7) is rewritten as:

[0034]

[0035] The equivalent quality factor represents the accumulation of the Q-filtering effect of the formation, that is, from the top formation to the current depth formation. When seismic waves pass through the i-th layer of the medium, from the top formation to the i-th layer, it is represented by the equivalent quality factor Q e，i , and the formula is expressed as:

[0036]

[0037] where T i represents the cumulative travel time of the i-th layer. When i = n, Equation (8) can be expressed by replacing the interlayer Q-value with the equivalent Q-value as:

[0038]

[0039] where Qe .n represents the equivalent quality factor of the nth layer. At this time, the pre-stack positive Q filtering formula related to the offset is written as:

[0040]

[0041] In the above solution, step 3 includes the following steps:

[0042] In step 3.1, based on the pre-stack Q value filtering formula constructed in step 2, the present invention further introduces the L2 norm regularization constraint to construct the inversion objective function. The specific operation is as follows:

[0043] Represent the positive Q filtering process using a matrix:

[0044] u = real(LU) (12)

[0045] Where u, L, and U are the spectra of the attenuated seismic record, the positive Q filtering operator, and the effective frequency of the non-attenuated seismic record respectively. Use Tikhonov regularization to construct the stabilized objective functional:

[0046] G = ||LU - u|| 2 + λ||WU|| 2 (13)

[0047] Where λ is a regularization factor and W is an identity operator, a first-order or second-order differential operator. The analytical solution of the objective functional G is expressed as:

[0048] U = (L T L + λW T W) -1 L T u (14);

[0049] In step 3.2, use the conjugate gradient algorithm to iteratively update the inversion objective function G constructed in step 3.1 to solve the spectrum U of the effective frequency of the non-attenuated seismic record. Specifically, it includes the following steps:

[0050] Solve equation (14) by the conjugate gradient iteration method. Assign all U values to 0 as the initial term of the first iteration, and calculate the data fitting difference r n , and use the data fitting difference r n To obtain the steepest ascent direction And the conjugate direction Furthermore, calculate the conjugate gradient iteration step size Its formula is as follows:

[0051] r n = LU - d (15)

[0052]

[0053] Among them, F T is the transpose of the Fréchet differential operator of L, is the coefficient for determining the conjugate direction, and its value in the first iteration is 0, so as to iteratively update the compensated data U:

[0054]

[0055] In the above solution, the said step 4 includes the following steps:

[0056] Judge whether the iteratively updated data U satisfies the condition ||r n || 2 ≤ Tol, where Tol is the given error term value. If so, output the final iteratively inverted data U as the final pre-stack compensation data. Otherwise, use the final iteratively inverted data U as the initial value and go to step 3.

[0057] The present invention also provides a pre-stack Q-value attenuation compensation device based on L2-norm regularization constraint, including:

[0058] A preprocessing module for preprocessing the pre-stack seismic record and obtaining the pre-stack Q-value and incident angle information;

[0059] A pre-stack Q-value filtering construction module, according to the frequency-domain wavefield continuation formula of plane waves in a viscoelastic medium, constructs a pre-stack Q-value filtering formula related to the offset using the pre-stack Q-value and incident angle information obtained by the preprocessing module, and establishes a pre-stack forward model related to the offset;

[0060] An inversion objective function construction module, based on the pre-stack Q-value filtering formula constructed by the pre-stack Q-value filtering formula construction module, introduces regularization to construct an inversion objective function, and uses the conjugate gradient algorithm to iteratively update the seismic record in the initial frequency domain;

[0061] An iteration termination condition judgment module for judging whether the frequency-domain seismic record updated by the inversion objective function construction module satisfies the iteration termination condition. If so, output the updated frequency-domain seismic record as the final inversion result. Otherwise, use the updated frequency-domain seismic record as a new initial model and return to the inversion objective function construction module for re-iterative update until the iteration termination condition is satisfied.

[0062] In the above device, the implementation of the preprocessing module includes the following steps:

[0063] Step 1.1: Use the relationship between the time, offset and Q of the pre-stack seismic record, and use the pre-stack Q-value extraction method to extract the Q-value of the pre-stack seismic record;

[0064] Step 1.2: Obtain the incident angle θ of each layer in the prestack seismic record by using ray tracing technology i .

[0065] In the above device, the prestack Q-value filtering formula construction module realizes the following steps:

[0066] Step 2.1: Establish a propagation model of seismic waves in the medium according to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media. The frequency-domain wavefield continuation formula of plane waves is expressed as:

[0067]

[0068] where b(T,ω) represents the wavefield with propagation time T and angular frequency ω, and ΔT i is the propagation time increment of seismic waves in the i-th layer, k(ω) is the frequency-dependent wave number, v r is the velocity at the reference frequency, the variable j represents the imaginary unit, and e is the base of the natural logarithm;

[0069] Step 2.2: Use the obtained prestack Q-value and incident angle information to calculate the propagation time increment ΔT of seismic waves in the i-th layer according to Snell's law i , and substitute ΔT i into the frequency-domain wavefield continuation formula to obtain wavefield continuation:

[0070]

[0071] ΔT i is expressed as:

[0072]

[0073] where △T i，0 is the two-way propagation time increment of seismic waves vertically propagating in the i-th layer, θ i is the incident angle of each layer of the medium, and h i represents the thickness of the i-th path or layer;

[0074] Step 2.3: Then, according to the modified Kolsky model, introduce the attenuation and velocity dispersion characteristics of seismic waves into the wavefield continuation formula. The modified Kolsky model describes the viscoelastic characteristics of the medium by introducing the quality factor Q and the tuning frequency ω0, so as to more accurately simulate the propagation process of seismic waves. At this time, the wavefield continuation is expressed as:

[0075]

[0076] where ω0 is the tuning frequency, and Q i represents the quality factor of the i-th layer of the medium, i = 1, 2, 3...n;

[0077] Step 2.4: Using the time-shift property of Fourier transform, convert the wavefield continuation formula into pre-stack Q-value filtering related to the offset. Specifically:

[0078] Using the time-shift property of Fourier transform, we have:

[0079]

[0080] Substitute Equation (5) into Equation (4) to obtain the following wavefield continuation formula:

[0081]

[0082] When seismic waves propagate through n layers of media, the wavefield continuation formula at this time becomes:

[0083]

[0084] Replace the propagation time △T with the sampling interval dt i，0 and discretize the frequency ω into ω m , Equation (7) is rewritten as:

[0085]

[0086] The equivalent quality factor represents the accumulation of the Q-filtering effect of the formation, that is, from the top formation to the current depth formation. When seismic waves pass through the i-th layer of media, from the top formation to the i-th layer, the equivalent quality factor Q e，i is used to represent, and the formula is expressed as:

[0087]

[0088] where, T i represents the cumulative travel time of the i-th layer. When i = n, Equation (8) can be expressed by replacing the interlayer Q-value with the equivalent Q-value as:

[0089]

[0090] where, Q e，n represents the equivalent quality factor of the n-th layer. At this time, the pre-stack positive Q-filtering formula related to the offset can be written as:

[0091]

[0092] In the above device, the inversion objective function construction module realizes the following steps:

[0093] In step 3.1, the present invention further introduces the L2-norm regularization constraint based on the pre-stack Q-value filter constructed in the pre-stack Q-value filtering formula construction module to construct the inversion objective function. The specific operation is as follows:

[0094] The positive Q filtering process is represented using a matrix:

[0095] u = real(LU) (12)

[0096] where u, L, and U are the spectra of the attenuated seismic record, the positive Q filtering operator, and the effective frequency of the non-attenuated seismic record, respectively. A stabilized objective functional is constructed using Tikhonov regularization:

[0097] G = ||LU - u|| 2 + λ||WU|| 2 (13)

[0098] where λ is a regularization factor and W is an identity operator, a first-order or second-order differential operator. The analytical solution of the objective functional G is expressed as:

[0099] U = (L T L + λW T W) -1 L T u (14);

[0100] In step 3.2, the conjugate gradient algorithm is used to iteratively update the inversion objective function G constructed in step 3.1 to solve for the spectrum U of the effective frequency of the non-attenuated seismic record. This specifically includes the following steps:

[0101] Solve equation (14) using the conjugate gradient iteration method. Set the U value to 0 as the initial term for the first iteration and calculate the data fitting error r n , and use the data fitting error r n to calculate the steepest ascent direction and the conjugate direction and then calculate the conjugate gradient iteration step size The formula is as follows:

[0102] r n = LU - d (15)

[0103]

[0104]

[0105] where F T is the transpose of the Fréchet differential operator of L, is the coefficient for determining the conjugate direction, with a value of 0 in the first iteration, thereby iteratively updating the compensated data U:

[0106]

[0107] In the above device, the iteration termination condition judgment module is implemented including the following steps:

[0108] Determine whether the data U after iterative update satisfies the condition ||r n || 2 ≤ Tol, where Tol is the given error term value. If so, output the final iterative inversion data U as the final pre-stack compensation data. Otherwise, use the final iterative inversion data U as the initial value and go to the inversion objective function construction module.

[0109] The present invention adopts the above technical means, and thus has the following beneficial effects:

[0110] 1. By means of pre-stack seismic record preprocessing and obtaining Q-value and incident angle information, the problem of signal loss caused by Q-value attenuation in seismic records is solved, and the effect of providing an accurate data basis for subsequent attenuation compensation is achieved.

[0111] 2. By constructing a pre-stack Q-value filtering formula related to the offset and establishing a pre-stack forward modeling, the problem that the traditional inverse Q-filtering method cannot be effectively applied to pre-stack seismic records is solved, and the effect of improving the compensation accuracy of pre-stack seismic records is achieved.

[0112] 3. Based on the pre-stack Q-value filtering formula, introducing L2-norm regularization constraint to construct an inversion objective function and using the conjugate gradient algorithm to iteratively update the seismic record in the initial frequency domain, the problem of numerical instability in the inverse Q-filtering process is solved, and the effect of ensuring the stability of amplitude compensation is achieved.

[0113] 4. By judging whether the seismic record in the frequency domain after iterative update satisfies the iterative termination condition and performing cyclic iterative update, the problem of insufficient convergence that may occur in the compensation process is solved, and the effect of ensuring the accuracy of the final inversion result is achieved.

[0114] 5. Using the relationship between the time, offset and Q of the pre-stack seismic record to extract the Q-value of the pre-stack seismic record and obtaining the incident angle of each layer through ray tracing technology, the problem of inaccurate extraction of Q-value and incident angle information is solved, and the effect of providing reliable parameters for constructing the pre-stack Q-value filtering formula is achieved.

[0115] 6. Constructing an attenuation compensation model based on the modified Kolsky model and considering the attenuation and velocity dispersion characteristics of seismic waves, the problem that the traditional attenuation model cannot accurately simulate the propagation process of seismic waves is solved, and the effect of more accurately compensating the attenuation of seismic records is achieved.

[0116] 7. Representing the forward Q-filtering process using a matrix and constructing a stabilized objective functional using Tikhonov regularization, the problem of complex construction of the objective function in the inversion process is solved, and the effect of simplifying the inversion calculation process is achieved.

[0117] 8. The conjugate gradient algorithm is used to iteratively update the inversion objective function, which is a means to solve the spectrum of the effective frequency of the non-attenuated seismic record, solving the problem of low computational efficiency in the inversion process and achieving improved inversion speed and effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] Figure 1 Flow chart of the implementation scheme

[0119] Figure 2 Pre-stack seismic record test without noise influence, where (a) is the non-attenuated pre-stack synthetic seismic record, (b) is the attenuated pre-stack synthetic seismic record, and (c) is the pre-stack seismic record compensated by this method;

[0120] Figure 3 Pre-stack seismic record test with added noise, where (a) is the pre-stack synthetic seismic record with added noise, (b) is the seismic record compensated by the stability factor, (c) is the pre-stack seismic record compensated by the amplitude gain limit method, and (d) is the pre-stack seismic record compensated by this method. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0121] The following will give a detailed description of the embodiments of the present invention. Although the present invention will be described and explained in conjunction with some specific embodiments, it should be noted that the present invention is not limited to these embodiments. On the contrary, modifications or equivalent replacements made to the present invention should be covered within the scope of the claims of the present invention.

[0122] In addition, for a better illustration of the present invention, numerous specific details are given in the following detailed description. Those skilled in the art will understand that the present invention can be implemented without these specific details.

[0123] For the problems studied above: The present invention provides a pre-stack Q-value attenuation compensation method based on L2-norm regularization constraint. This method is based on the pre-stack Q-value filtering method of forward wavefield propagation. The most important step is to construct a pre-stack Q-value model according to the wave equation propagation mechanism, use this filtering model to regard compensation as an inversion problem, and use the regularization constraint strategy for inversion to ensure the numerical stability of amplitude compensation. This method can effectively improve the problem of insufficient compensation of seismic records in the case of deep layers and low Q-values by the inverse Q-filtering method.

[0124] To achieve the above object, the present invention adopts the following technical solutions:

[0125] A pre-stack Q-value attenuation compensation method based on L2-norm regularization constraint, the technical solution is as follows:

[0126] Step 1: Preprocess the pre-stack seismic record and obtain the pre-stack Q-value and incident angle information.

[0127] Step 2: According to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, use the obtained pre-stack Q value and incident angle information to construct a pre-stack Q value filtering formula related to offset, and establish a pre-stack forward modeling model related to offset.

[0128] Step 3: Introduce regularization to construct an inversion objective function, and use the conjugate gradient algorithm to iteratively update the seismic record in the initial frequency domain for the objective function.

[0129] Step 4: Determine whether the updated seismic record in the frequency domain satisfies the iteration termination condition. If so, output the updated seismic record in the frequency domain as the final inversion result. Otherwise, use the updated seismic record in the frequency domain as the new initial model, return to Step 3, and iterate and update again until the iteration termination condition is satisfied.

[0130] Furthermore, Step 1 includes the following steps:

[0131] Step 1.1: Use the relationship between the time, offset, and Q of the pre-stack seismic record, and use the pre-stack Q value extraction method to extract the Q value of the pre-stack seismic record.

[0132] Step 1.2: Obtain the incident angle θ of each layer in the pre-stack seismic record by using ray tracing technology i .

[0133] Furthermore, Step 2 includes the following steps:

[0134] Step 2.1: Construct an attenuation compensation attenuation model based on the modified Kolsky model:

[0135] The backward wavefield continuation of plane waves in the frequency domain can be expressed as:

[0136]

[0137] where b(T, ω) represents the wavefield with propagation time T and angular frequency ω, ΔT i is the propagation time increment of seismic waves in the i-th layer, k(ω) is the frequency-dependent wave number, and v r is the velocity at the reference frequency. The variable j represents the imaginary unit. According to Snell's law, ΔT i can be expressed as:

[0138]

[0139] where the ΔT i,0 is the two-way propagation time increment of seismic waves vertically propagating in the i-th layer, and θ i is the incident angle of each layer of medium.

[0140] The wavefield continuation obtained by substituting Equation (2) into Equation (1) is:

[0141]

[0142] Among them, the improved Kolsky model is used to describe the velocity dispersion and seismic attenuation. At this time, the wavefield continuation can be expressed as:

[0143]

[0144] where ω0 is the tuning frequency, Compared with the traditional post-stack inverse Q attenuation model, the incident angle θ is included in both exponential operators, i so that the wavefield continuation is related to the offset, and then it is applied to pre-stack record compensation.

[0145] Using the time-shift property of the Fourier transform, we have:

[0146]

[0147] Substituting Equation (5) into Equation (4) gives the following wavefield continuation formula:

[0148]

[0149] When the seismic wave propagates through n layers of media, the wavefield continuation formula at this time becomes:

[0150]

[0151] In actual calculations, the sampling interval dt is used to replace the propagation time ΔT i,0 , and the frequency ω is discretized into ω m .

[0152] At this time, Equation (7) is rewritten as:

[0153]

[0154] At this time, the wavefield continuation formula contains an additional cosθ i , making the recursive calculation difficult. At this time, it is considered that the medium through which the seismic wave propagates is approximately homogeneous and equivalent. In this approximation, the interlayer quality factor Q n is represented by the effective quality factor Q e,n as:

[0155]

[0156] where Q n represents the interval Q of the nth layer, and T n represents the cumulative travel time of the nth layer. Therefore, Equation (8) can be approximated as:

[0157]

[0158] At this time, the pre-stack positive Q formula related to the offset can be written as:

[0159]

[0160]

[0161] Furthermore, step 3 includes the following steps:

[0162] Step 3.1: Represent the positive Q filtering process using a matrix:

[0163] u = real(LU) (12)

[0164] where u, L, and U are the spectra of the attenuated seismic record, the positive Q filtering operator, and the effective frequency of the non-attenuated seismic record respectively. A stabilized objective functional is constructed using Tikhonov regularization:

[0165] G = ||LU - u|| 2 + λ||WU|| 2 (13)

[0166] where λ is a regularization factor and W is an identity operator, a first-order or second-order differential operator. Its analytical solution can be expressed as:

[0167] U = (L T L + λW T W) -1 L T u (14)

[0168] Step 3.2: Equation (14) is solved by the conjugate gradient iteration method. The U value is fully assigned 0 as the initial term for the first iteration, and the data fitting difference r n , and the data fitting difference r n is used to obtain the steepest ascent direction and the conjugate direction Then, the conjugate gradient iteration step size is calculated. Its formula is as follows:

[0169] r n = LU - d (15)

[0170]

[0171] where F T is the transpose of the Fréchet differential operator of L, is the coefficient for determining the conjugate direction, and its value in the first iteration is 0. Thus, the compensated data U is iteratively updated:

[0172]

[0173] Further, step 4 includes the following steps:

[0174] Step 4.1: Further, step 4 is to determine whether the iteratively updated data U satisfies the condition ||r n || 2 ≤ Tol, where Tol is a given error term value. If so, the final iteratively inverted data U is output as the final prestack compensation data. Otherwise, the final iteratively inverted data U is used as the initial value and the process goes to step 3.

[0175] Embodiment 1

[0176] To verify the feasibility of the method in this embodiment, the present invention is described by taking the synthetic CMP seismic record as an example, as Figure 2 (a) is the unattenuated prestack synthetic seismic record, which is obtained by convolving with a 50 Hz Ricker wavelet. The model Q - value parameters are 40, 60, and 100 respectively, the geophone spacing is 100 m, and there are 12 single - channel data in total. The model parameters are introduced into the prestack Q - value attenuation to obtain Figure 2 (b) the attenuated prestack CMP seismic record. The prestack seismic record inversely compensated by using this method under the condition of no noise influence is as shown in Figure 2 (c). It can be seen that under the condition of no noise, this method compensates the seismic attenuation well and basically coincides with the original seismic record. Further, to verify the compensation benefit of this method under the influence of noise, this method is compared with the conventional stable factor method and amplitude gain limiting method. Figure 3 (a) is the prestack synthetic seismic record with added Gaussian white noise with a signal - to - noise ratio of 30%, Figure 3 (b) is the seismic record compensated by the stable factor, Figure 3 (c) is the prestack seismic record compensated by the amplitude gain limiting method, Figure 3 (d) is the prestack seismic compensated by this method. By comparison, it can be seen that the stable factor method can stably compensate the prestack seismic record, but the amplitude compensation of the prestack seismic record is insufficient. The amplitude gain limiting method can effectively avoid the distortion during the amplitude compensation of the deep layer by restricting the amplitude compensation of the high - frequency components, but the amplitude compensation of the seismic record is also limited. At the same time, if the selection of high - frequency components is inappropriate, it is easy to amplify the noise. This method has a better improvement in terms of amplitude and noise suppression compared with the above two methods, verifying the noise resistance and compensation benefit of the proposed method for deep - layer seismic records.

[0177] In summary, the present invention proposes a prestack Q - value attenuation compensation method based on L2 - norm regularization constraint. By deriving the wave - field continuation formula, the inverse Q - filtering method is applied to the prestack seismic record, and L2 - norm regularization is used for inverse compensation, improving the stability and high precision of the compensation for the prestack seismic record.

[0178] The present invention applies the inverse Q filtering method to the prestack by constructing a prestack Q-value attenuation model. Meanwhile, the instability of the inverse Q filtering method is avoided through the method of L2-norm regularization constraint inversion, overcoming the influence of insufficient compensation of the inverse Q filtering method. In the case of deep layers and low Q values, it has good stability and high precision in compensating seismic records.

Claims

1. A pre-stack Q-value attenuation compensation method based on L2-norm regularization constraint, characterized in that, Including the following steps: Step 1: Preprocess the prestack seismic records and obtain prestack Q value and incident angle information; Step 2: According to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, use the obtained prestack Q value and incident angle information to construct a prestack Q value formula related to offset, and establish a prestack attenuation model related to offset; Step 3: Based on the prestack attenuation model constructed in Step 2, introduce regularization to construct an inversion objective function, and use the conjugate gradient algorithm to iteratively update the seismic records in the initial frequency domain for the objective function; Step 4: Determine whether the updated seismic records in the frequency domain satisfy the iteration termination condition. If so, output the updated seismic records in the frequency domain as the final inversion result. Otherwise, use the updated seismic records in the frequency domain as a new initial model, return to Step 3, and iterate and update again until the iteration termination condition is satisfied.

2. The prestack Q - value attenuation compensation method based on L2 - norm regularization constraint according to claim 1, wherein, The said Step 1 includes the following steps: Step 1.1: Utilize the relationship among the time, offset, and Q of the prestack seismic records, and use the prestack Q value extraction method to extract the Q value of the prestack seismic records; Step 1.2: Obtain the incident angle θ of each layer in the pre-stack seismic record by using ray tracing technology i .

3. A prestack Q-value attenuation compensation method based on L2-norm regularization constraint according to claim 2, characterized in that, The said Step 2 includes the following steps: Step 2.1: According to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, establish a propagation model of seismic waves in the medium. The frequency-domain wavefield continuation formula of plane waves is expressed as: where \(B(T,\omega)\) represents the wave field with propagation time \(T\) and angular frequency \(\omega\), \(\Delta T\) i is the propagation time increment of seismic waves in the \(i\)-th layer, \(k(\omega)\) is the frequency-dependent wave number, \(v\) r is the velocity at the reference frequency, the variable \(j\) represents the imaginary unit, and \(e\) is the base of the natural logarithm; Step 2.2: Using the obtained pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of seismic waves in the i-th layer according to Snell's law i , and substitute ΔT i into the frequency-domain wavefield continuation formula to obtain wavefield continuation: ΔT i Expressed as: where ΔT i,0 is the two-way propagation time increment of the seismic wave propagating vertically in the i-th layer, θ i is the incident angle of each layer of medium, h i represents the thickness of the i-th path or layer; Step 2.3: Then, according to the modified Kolsky model, introduce the attenuation and velocity dispersion characteristics of seismic waves into the wavefield continuation formula. The modified Kolsky model describes the viscoelastic characteristics of the medium by introducing the quality factor Q and the tuning frequency ω0, so as to more accurately simulate the propagation process of seismic waves. At this time, the wavefield continuation is expressed as: where ω0 is the tuning frequency, Q i represents the quality factor of the i-th layer of dielectric, Step 2.4: Utilize the time-shift property of the Fourier transform to convert the wavefield continuation formula into a prestack attenuation model related to offset. Specifically: Using the time-shift property of the Fourier transform, we have: Substitute Equation (5) into Equation (4) to obtain the following wavefield continuation formula: When the seismic wave propagates through n layers of media, the wavefield continuation formula at this time becomes: Replace the propagation time ΔT with the sampling interval dt i,0 , and discretize the frequency ω into ω m , Equation (7) is rewritten as: The equivalent quality factor represents the accumulation of the Q-filtering effect of the formation, that is, from the top formation to the current depth formation. When the seismic wave passes through the i-th layer of medium, the equivalent quality factor Q is used from the top formation to the i-th layer e,i is expressed, and the formula is expressed as: Among them, T i represents the cumulative travel time of the i-th layer. When i = n, formula (8) can be expressed by replacing the interlayer Q value with the equivalent Q value as follows: Among them, Q e,n represents the equivalent quality factor of the nth layer. At this time, the prestack forward Q attenuation formula related to the offset is written as:

4. A pre-stack Q-value attenuation compensation method based on L2-norm regularization constraint according to claim 1, characterized in that The said Step 3 includes the following steps: In Step 3.1, based on the prestack Q value attenuation formula constructed in Step 2, further introduce the L2 norm regularization constraint to construct an inversion objective function. The specific operation is as follows: Represent the positive Q filtering process using a matrix: u = real(LU) (12) Where u, L, and U are the spectra of the effective frequencies of the attenuated seismic records, the positive Q filtering operator, and the non-attenuated seismic records respectively. Use Tikhonov regularization to construct a stabilized objective functional: G = ||LU - u|| 2 + λ||WU|| 2 (13) Where λ is a regularization factor, W is an identity operator, a first-order or second-order differential operator, and the analytical solution of the objective functional G is expressed as: U = (L T L + λW T W) -1 L T u (14); In Step 3.2, adopt the conjugate gradient algorithm to iteratively update the inversion objective function G constructed in Step 3.1 to solve the spectrum U of the effective frequencies of the non-attenuated seismic records. Specifically, it includes the following steps: Solve Equation (14) by the conjugate gradient iteration method, assign all U values to 0 as the initial term for the first iteration, and calculate the data fitting error r n , and use the data fitting error r n to calculate the steepest ascent direction and the conjugate direction Furthermore, calculate the conjugate gradient iteration step size The formula is as follows: r n = LU-d (15) where, F T is the transpose of the Fréchet differential operator of L, is the coefficient for determining the conjugate direction, with a value of 0 in the first iteration, thereby iteratively updating the compensated data U:

5. A prestack Q-value attenuation compensation method based on L2-norm regularization constraint according to claim 1, characterized in that The said Step 4 includes the following steps: Determine whether the data U after iterative update satisfies the condition ||r n || 2 ≤ Tol, where Tol is the given error term value. If so, output the final iterative inversion data U as the final prestack compensation data. Otherwise, use the final iterative inversion data U as the initial value and go to step 3.

6. A pre-stack Q-value attenuation compensation device based on L2-norm regularization constraint, characterized in that, Including: A preprocessing module for preprocessing the prestack seismic records and obtaining prestack Q value and incident angle information; Pre-stack Q-value filtering construction module, according to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, uses the pre-stack Q-value and incident angle information obtained by the preprocessing module to construct pre-stack Q-value filtering related to offset, and establishes a pre-stack attenuation model related to offset; Inversion objective function construction module, based on the pre-stack Q-value filter constructed by the pre-stack Q-value filtering construction module, introduces regularization to construct an inversion objective function, and uses the conjugate gradient algorithm to iteratively update the seismic record in the initial frequency domain; Iterative termination condition judgment module, used to judge whether the frequency-domain seismic record updated by the inversion objective function construction module meets the iterative termination condition. If it meets, the updated frequency-domain seismic record is output as the final inversion result. Otherwise, the updated frequency-domain seismic record is used as a new initial model and returned to the inversion objective function construction module for iterative update again until the iterative termination condition is met.

7. An apparatus for prestack Q-value attenuation compensation based on L2-norm regularization constraint according to claim 6, wherein The preprocessing module implementation includes the following steps: Step 1.1: Using the relationship between the time, offset, and Q of the pre-stack seismic record, the pre-stack Q-value extraction method is used to extract the Q-value of the pre-stack seismic record; Step 1.2: Obtain the incident angle θ of each layer in the pre-stack seismic record by using ray tracing technology i .

8. The pre-stack Q-value attenuation compensation device based on L2-norm regularization constraint according to claim 7, characterized in that, The pre-stack Q-value filtering construction module implementation includes the following steps: Step 2.1: According to the frequency-domain wavefield continuation formula of plane waves in viscoelastic media, establish a propagation model of seismic waves in the medium. The frequency-domain wavefield continuation formula of plane waves is expressed as: where \(B(T,\omega)\) represents the wave field with propagation time \(T\) and angular frequency \(\omega\), \(\Delta T\) i is the propagation time increment of the seismic wave in the \(i\)-th layer, \(k(\omega)\) is the frequency-dependent wavenumber, \(v\) r is the velocity at the reference frequency, the variable \(j\) represents the imaginary unit, and \(e\) is the base of the natural logarithm; Step 2.2: Using the obtained pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of seismic waves in the i-th layer according to Snell's law i , and substitute ΔT i into the frequency-domain wavefield continuation formula to obtain wavefield continuation: ΔT i Expressed as: where ΔT i,0 is the two-way propagation time increment of the seismic wave vertically propagating through the i-th layer, and θ i is the incident angle of each layer of medium, and h i represents the thickness of the i-th path or layer; Step 2.3: Then, according to the modified Kolsky model, introduce the attenuation and velocity dispersion characteristics of seismic waves into the wavefield continuation formula. The modified Kolsky model describes the viscoelastic characteristics of the medium by introducing the quality factor Q and the tuning frequency ω0, so as to more accurately simulate the propagation process of seismic waves. At this time, the wavefield continuation is expressed as: where ω0 is the tuning frequency, and Q i represents the quality factor of the i-th layer of the medium, Step 2.4: Using the time-shift property of the Fourier transform, convert the wavefield continuation formula into a pre-stack Q-value filtering formula related to offset. Specifically: Using the time-shift property of the Fourier transform, we have: Substitute Equation (5) into Equation (4) to obtain the following wavefield continuation formula: When seismic waves propagate through n layers of media, the wavefield continuation formula at this time becomes: Replace the propagation time ΔT with the sampling interval dt i,0 , and discretize the frequency ω into ω m , Equation (7) is rewritten as: The equivalent quality factor represents the accumulation of the Q-filtering effect of the formation, that is, from the top formation to the current depth formation. When the seismic wave passes through the i-th layer of medium, the equivalent quality factor Q is used from the top formation to the i-th layer. e,i It is expressed by the formula as follows: Among them, T i represents the cumulative travel time of the i-th layer. When i = n, formula (8) can be expressed by replacing the interlayer Q value with the equivalent Q value as follows: where Q e,n represents the equivalent quality factor of the nth layer. At this time, the prestack forward Q-filtering formula related to the offset is written as:

9. The pre-stack Q-value attenuation compensation device based on L2-norm regularization constraint according to claim 6, characterized in that, The inversion objective function construction module implementation includes the following steps: In Step 3.1, further introduce the L2 norm regularization constraint to the pre-stack Q-value filter constructed in the pre-stack Q-value filtering construction module to construct an inversion objective function. The specific operation is as follows: Represent the positive Q filtering process using a matrix: u = real(LU) (12) Among them, u, L, and U are the spectra of the effective frequencies of the attenuated seismic record, the positive Q filtering operator, and the non-attenuated seismic record respectively. Use Tikhonov regularization to construct a stabilized objective functional: G = ||LU - u|| 2 + λ||WU|| 2 (13) Among them, λ is a regularization factor, W is an identity operator, a first-order or second-order differential operator, and the analytical solution of the objective functional G is expressed as: U = (L T L + λW T W) -1 L T u (14); In Step 3.2, use the conjugate gradient algorithm to iteratively update the inversion objective function G constructed in Step 3.1 to solve the spectrum U of the effective frequencies of the non-attenuated seismic record. Specifically, it includes the following steps: Solve equation (14) by the conjugate gradient iteration method. Assign all U values to 0 as the initial term for the first iteration, and calculate the data fitting error r n , and use the data fitting error r n to calculate the steepest ascent direction and the conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows: r n = LU - d (15) where F T is the transpose of the Fréchet differential operator of L, is the coefficient for determining the conjugate direction, with a value of 0 in the first iteration, and thus iteratively updates the compensated data U:

10. A prestack Q - value attenuation compensation device based on L2 - norm regularization constraint according to claim 7, characterized in that, The iterative termination condition judgment module implementation includes the following steps: Determine whether the data U after iterative update satisfies the condition ||r n || 2 ≤ Tol, where Tol is the given error term value. If so, output the final iterative inversion data U as the final prestack compensation data; otherwise, use the final iterative inversion data U as the initial value and go to the inversion objective function construction module.

Images