Pre-stack Q value attenuation compensation method and device based on L2 norm regularization constraint

Abstract

The invention relates to the technical field of oil-gas exploration and pre-stack seismic inverse Q compensation, discloses a pre-stack Q value attenuation compensation method and device based on L2 norm regularization constraint, and aims to solve the problem of signal attenuation caused by Q value attenuation in pre-stack seismic data. According to the main scheme, the method comprises the steps of preprocessing pre-stack seismic data and obtaining a pre-stack Q value and incident angle information; according to a frequency domain wave field continuation formula of a plane wave in a viscoelastic medium, constructing a pre-stack Q value filtering method related to the offset by using the obtained pre-stack Q value and incident angle information, and establishing a pre-stack attenuation model related to the offset; introducing regularization to construct an inversion objective function, iteratively updating the seismic record of the initial frequency domain for the objective function by using a conjugate gradient algorithm, judging whether the seismic record of the updated frequency domain meets an iteration termination condition or not, if so, outputting the seismic record, and if not, outputting the seismic record; and iteratively updating the seismic record of the updated frequency domain again until an iteration termination condition is met.

Description

Technical Field

[0001] The 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 underground propagation of seismic waves, due to the viscoelastic properties of the underground medium, their amplitude and phase will be affected by attenuation, thereby reducing the resolution of seismic records. As an effective compensation method, inverse Q filtering plays a vital role in improving the resolution of seismic records. At present, the inverse Q filtering method is mainly used in post-stack seismic records. The inverse Q filtering processing of pre-stack seismic records is conducive to improving the resolution and fidelity of seismic records to meet the requirements of today's seismic exploration for high resolution and high fidelity of seismic records.

[0003] The inverse Q filtering method needs to consider the influence of offset distance when compensating prestack seismic records. The inverse Q filtering method of amplitude gain limitation and stability factor proposed by Wang can stably compensate prestack seismic records to a certain extent, but it is insufficient for deep layers and low Q values. In summary, the existing prestack Q value attenuation compensation method technology has large errors in compensating seismic records in deep layers and low Q values. Summary of the invention

[0004] The purpose of the present invention is to solve the signal loss problem caused by Q value attenuation in pre-stack seismic records, and to achieve high-precision and stable compensation for seismic records by constructing a Q value filtering formula related to offset distance and a regularized inversion method.

[0005] In order 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, comprising the following steps:

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

[0008] Step 2: Based on the frequency domain wave field extension formula of plane waves in viscoelastic media, the prestack Q value and incident angle information are obtained to construct a prestack Q value filtering formula related to the offset, and a prestack attenuation model related to the offset is established;

[0009] Step 3: Based on the pre-stack Q-value filter constructed in step 2, regularization is introduced to construct the inversion objective function, and the conjugate gradient algorithm is used to iteratively update the initial frequency domain seismic records.

[0010] Step 4: Determine whether the updated seismic record in the frequency domain meets 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 and return to step 3 to iterate again until the iteration termination condition is met.

[0011] In the above scheme, step 1 comprises the following steps:

[0012] Step 1.1: Using the relationship between the time, offset and Q of the prestack seismic record, the Q value of the prestack seismic record is extracted 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 ray tracing technology i .

[0014] In the above scheme, step 2 comprises the following steps:

[0015] Step 2.1: According to the frequency domain wave field extension formula of plane waves in viscoelastic media, a propagation model of seismic waves in the medium is established. The frequency domain wave field extension formula of plane waves is expressed as:

[0016]

[0017] Among them, B(T,ω) represents the wave field with propagation time T and angular frequency ω, ΔT i is the propagation time increment of the seismic wave in the i-th layer, k(ω) is the frequency-dependent wave number, v r is the speed of 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 acquired pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of the seismic wave in the i-th layer according to Snell's law i , and ΔT i Substituting into the frequency domain wave field extension formula, we get the wave field extension:

[0019]

[0020] ΔT i Expressed as:

[0021]

[0022] The ΔT i,0 is the two-way propagation time increment of the seismic wave vertically propagating 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;

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

[0024]

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

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

[0027] The time-shift properties of Fourier transform are:

[0028]

[0029] Substituting formula (5) into formula (4) yields the following wave field extension formula:

[0030]

[0031] When the seismic wave propagates through n layers of media, the wave field extension formula becomes:

[0032]

[0033] Use the sampling interval dt instead of the propagation time △T i，0 , and discretize the frequency ω into ω m , formula (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 the seismic wave passes through the i-th layer of the medium, the equivalent quality factor Q is used from the top formation to the i-th layer. e，i The formula is:

[0036]

[0037] Among them, T i represents the accumulated travel time of layer i. When i=n, ​​formula (8) can be expressed as follows by replacing the inter-layer Q value with the equivalent Q value:

[0038]

[0039] Among them, Qe .n represents the equivalent quality factor of the nth layer. At this time, the prestack positive Q filter formula related to the offset distance is written as:

[0040]

[0041] In the above scheme, step 3 comprises the following steps:

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

[0043] The positive Q filtering process is represented by a matrix:

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

[0045] Among them, u, L and U are the spectra of the effective frequencies of the attenuated seismic record, the positive Q filter operator and the unattenuated seismic record, respectively. The stabilized target functional is constructed using Tikhonov regularization:

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

[0047] Among them, λ is a regularization factor, W is a unit operator, a first-order or second-order differential operator, and the analytical solution of the target functional G is expressed as:

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

[0049] In step 3.2, the inversion objective function G constructed in step 3.1 is iteratively updated using the conjugate gradient algorithm to solve the spectrum U of the effective frequency of the unattenuated seismic record, which specifically includes the following steps:

[0050] The conjugate gradient iteration method is used to solve equation (14), all U values ​​are assigned to 0 as the initial term of the first iteration, and the data fitting error r is calculated. n , and use the data to fit the difference r n Find the direction of fastest ascent and conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows:

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

[0052]

[0053] Among them, F T is the transpose of the Frechet differential operator of L, To determine the coefficient of the conjugate direction, the value in the first iteration is 0, so as to iteratively update the compensated data U:

[0054]

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

[0056] Determine whether the iteratively updated data U meets the condition || r n || 2 ≤Tol, where Tol is the given error term value. If yes, the last iterative inversion data U is output as the final prestack compensation data. Otherwise, the last iterative inversion data U is used 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, comprising:

[0058] A preprocessing module is used to preprocess the prestack seismic records and obtain the prestack Q value and incident angle information;

[0059] The pre-stack Q-value filter construction module constructs a pre-stack Q-value filter formula related to the offset distance based on the frequency-domain wave field extension formula of plane waves in viscoelastic media and uses the pre-stack Q-value and incident angle information obtained by the pre-processing module to build a pre-stack Q-value filter formula related to the offset distance, and establishes a pre-stack forward model related to the offset distance;

[0060] The inversion objective function construction module introduces regularization to construct the inversion objective function based on the prestack Q-value filter formula constructed by the prestack Q-value filter formula construction module, and uses the conjugate gradient algorithm to iteratively update the seismic records in the initial frequency domain;

[0061] The iteration termination condition judgment module is used to judge whether the frequency domain seismic record updated by the inversion objective function construction module meets the iteration termination condition. If so, the updated frequency domain seismic record is output as the final inversion result. Otherwise, the updated frequency domain seismic record is used as the new initial model and returned to the inversion objective function construction module for another iterative update until the iteration termination condition is met.

[0062] In the above device, the preprocessing module is implemented including the following steps:

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

[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 pre-stack Q value filter formula construction module is implemented by the following steps:

[0066] Step 2.1: According to the frequency domain wave field extension formula of plane waves in viscoelastic media, a propagation model of seismic waves in the medium is established. The frequency domain wave field extension formula of plane waves is expressed as:

[0067]

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

[0069] Step 2.2: Using the acquired pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of the seismic wave in the i-th layer according to Snell's law i , and ΔT i Substituting into the frequency domain wave field extension formula, we get the wave field extension:

[0070]

[0071] ΔT i Expressed as:

[0072]

[0073] Among them, △T i，0 is the two-way propagation time increment of the seismic wave vertically propagating 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;

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

[0075]

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

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

[0078] The time-shift properties of Fourier transform are:

[0079]

[0080] Substituting formula (5) into formula (4) yields the following wave field extension formula:

[0081]

[0082] When the seismic wave propagates through n layers of media, the wave field extension formula becomes:

[0083]

[0084] Use the sampling interval dt instead of the propagation time △T i，0 , and discretize the frequency ω into ω m , formula (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 the seismic wave passes through the i-th layer of the medium, the equivalent quality factor Q is used from the top formation to the i-th layer. e，i The formula is:

[0087]

[0088] Among them, T i represents the accumulated travel time of layer i. When i=n, ​​formula (8) can be expressed as follows by replacing the inter-layer Q value with the equivalent Q value:

[0089]

[0090] Among them, Q e，n represents the equivalent quality factor of the nth layer. At this time, the prestack positive Q filter formula related to the offset distance can be written as:

[0091]

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

[0093] In step 3.1, the present invention further introduces an L2 norm regularization constraint based on the prestack Q value filter constructed in the prestack Q value filter formula construction module to construct an inversion objective function. The specific operations are as follows:

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

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

[0096] Among them, u, L and U are the spectra of the effective frequencies of the attenuated seismic record, the positive Q filter operator and the unattenuated seismic record, respectively. The stabilized target functional is constructed using Tikhonov regularization:

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

[0098] Among them, λ is a regularization factor, W is a unit operator, a first-order or second-order differential operator, and the analytical solution of the target functional G is expressed as:

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

[0100] In step 3.2, the inversion objective function G constructed in step 3.1 is iteratively updated using the conjugate gradient algorithm to solve the spectrum U of the effective frequency of the unattenuated seismic record, which specifically includes the following steps:

[0101] The conjugate gradient iteration method is used to solve equation (14), all U values ​​are assigned to 0 as the initial term of the first iteration, and the data fitting error r is calculated. n , and use the data to fit the difference r n Find the direction of fastest ascent and conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows:

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

[0103]

[0104]

[0105] Among them, F T is the transpose of the Frechet differential operator of L, To determine the coefficient of the conjugate direction, the value in the first iteration is 0, so as to iteratively update 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 iteratively updated data U meets the condition || r n || 2 ≤Tol, where Tol is the given error term value. If yes, the last iterative inversion data U is output as the final prestack compensation data. Otherwise, the last iterative inversion data U is used as the initial value and goes to the inversion objective function construction module.

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

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

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

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

[0113] 4. By judging whether the seismic records in the frequency domain after iterative updating meet the iteration termination conditions and performing cyclic iterative updating, the problem of insufficient convergence that may occur in the compensation process is solved, thereby achieving the effect of ensuring the accuracy of the final inversion results.

[0114] 5. The Q value of the pre-stack seismic record is extracted by using the relationship between the pre-stack seismic record time, offset distance and Q, and the incident angle of each layer is obtained through ray tracing technology, which solves the problem of inaccurate extraction of Q value and incident angle information and achieves the effect of providing reliable parameters for constructing the pre-stack Q value filtering formula.

[0115] 6. An attenuation compensation model based on the modified Kolsky model is constructed, and the attenuation and velocity dispersion characteristics of seismic waves are taken into consideration. This solves the problem that the traditional attenuation model cannot accurately simulate the seismic wave propagation process, and achieves the effect of more accurate compensation for the attenuation of seismic records.

[0116] 7. The positive Q filtering process is represented by a matrix, and Tikhonov regularization is used to construct a stabilized target functional, which solves the problem of complex target function construction in the inversion process and achieves the effect of simplifying the inversion calculation process.

[0117] 8. The conjugate gradient algorithm is used to iteratively update the inversion objective function and solve the spectrum of the effective frequency of the unattenuated seismic record, which solves the problem of low computational efficiency in the inversion process and improves the inversion speed and effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] Figure 1 Flow chart of the implementation plan

[0119] Figure 2 Prestack seismic record test without noise influence, wherein (a) is an unattenuated prestack synthetic seismic record, wherein (b) is an attenuated prestack synthetic seismic record, wherein (c) is a prestack seismic record compensated by the present method;

[0120] Figure 3 Test of noisy prestack seismic records, including (a) noisy prestack synthetic seismic records, (b) seismic records compensated by stability factor, (c) prestack seismic records compensated by amplitude gain limitation method, and (d) prestack seismic records compensated by this method. DETAILED DESCRIPTION

[0121] The following is a detailed description of the embodiments of the present invention. Although the present invention will be described and illustrated 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 substitutions made to the present invention should all be included in the scope of the claims of the present invention.

[0122] In addition, in order to better illustrate the present invention, numerous specific details are given in the following specific embodiments. It will be understood by those skilled in the art that the present invention can also be implemented without these specific details.

[0123] In view of the above research problems: the present invention provides a prestack Q-value attenuation compensation method based on L2 norm regularization constraint. The method is based on the prestack Q-value filtering method of forward wave field propagation. The most important step is to build a prestack Q-value model according to the wave equation propagation mechanism, use the filtering model to treat compensation as an inversion problem, and use the regularization constraint strategy to invert to ensure the numerical stability of amplitude compensation. This method can effectively improve the problem of insufficient compensation of seismic records in deep layers and with low Q values ​​by the inverse Q filtering method.

[0124] In order to achieve the above object, the present invention adopts the following technical solution:

[0125] A prestack Q-value attenuation compensation method based on L2 norm regularization constraint, the technical scheme of which is as follows:

[0126] Step 1: Preprocess the prestack seismic records and obtain the prestack Q value and incident angle information.

[0127] Step 2: Based on the frequency domain wave field extension formula of plane waves in viscoelastic media, the prestack Q value and incident angle information are obtained to construct a prestack Q value filtering formula related to the offset distance, and a prestack forward model related to the offset distance is established.

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

[0129] Step 4: Determine whether the updated seismic record in the frequency domain meets 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 and return to step 3 to iterate again until the iteration termination condition is met.

[0130] Further, the step 1 comprises the following steps:

[0131] Step 1.1: Using the relationship between the time, offset and Q of prestack seismic records, the Q value of prestack seismic records is extracted using the prestack Q value extraction method.

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

[0133] Further, the step 2 comprises the following steps:

[0134] Step 2.1: Construct a decay compensation decay model based on the modified Kolsky model:

[0135] The backward wave field extension of a plane wave in the frequency domain can be expressed as:

[0136]

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

[0138]

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

[0140] Substituting equation (2) into equation (1), the wave field extension is:

[0141]

[0142] The improved Kolsky model is used to describe velocity dispersion and seismic attenuation. At this time, the wave field extension can be expressed as:

[0143]

[0144] Where ω0 is the tuning frequency, Compared with the traditional post-stack inverse Q decay model, both exponential operators include the incident angle θ i The information of the wave field can be used to correlate the wave field extension with the offset distance, which can then be used for pre-stack record compensation.

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

[0146]

[0147] Substituting formula (5) into formula (4) yields the following wave field extension formula:

[0148]

[0149] When the seismic wave propagates through n layers of media, the wave field extension formula becomes:

[0150]

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

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

[0153]

[0154] At this time, the wave field extension formula contains an additional cosθ i , making recursive calculation difficult. At this time, it is considered that the medium in which the seismic wave propagates is approximately uniform and equivalent. In this approximate case, the interlayer quality factor Q n Effective quality factor Q e,n Expressed as:

[0155]

[0156] Where Q n represents the interval Q, T of the nth layer n represents the accumulated travel time of n layers. Therefore, formula (8) can be approximated as:

[0157]

[0158] At this time, the prestack positive Q formula related to the offset distance can be written as:

[0159]

[0160]

[0161] Further, the step 3 comprises the following steps:

[0162] Step 3.1: Use the matrix to represent the positive Q filtering process:

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

[0164] Among them, u, L and U are the spectra of the effective frequencies of the attenuated seismic records, the positive Q filter operator and the unattenuated seismic records, respectively. The stabilized target functional is constructed using Tikhonov regularization:

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

[0166] Among them, λ is a regularization factor, W is a unit operator, a first-order or second-order differential operator, and 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. All U values ​​are assigned to 0 as the initial term of the first iteration, and the data fitting error r is calculated. n , and use the data to fit the difference r n Find the direction of fastest ascent and conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows:

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

[0170]

[0171] Among them, F T is the transpose of the Frechet differential operator of L, To determine the coefficient of the conjugate direction, the value in the first iteration is 0, so as to iteratively update the compensated data U:

[0172]

[0173] Further, the step 4 comprises 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 the given error term value. If yes, the last iterative inversion data U is output as the final prestack compensation data. Otherwise, the last iterative inversion data U is used as the initial value and go to step 3.

[0175] Example 1

[0176] In order to verify the feasibility of the method in this embodiment, the present invention is described using synthetic CMP seismic records as a column, such as Figure 2 (a) is an unattenuated prestack synthetic seismic record, which is obtained by convolution of 50Hz Ricker wavelet. The model Q value parameters are 40, 60, and 100, the detector spacing is 100m, 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) Attenuated pre-stack CMP seismic record. The pre-stack seismic record inverted and compensated using this method without noise influence is shown in Figure 2. Figure 2 As shown in (c), it can be seen that in the absence of noise, this method has a good compensation for seismic attenuation and is basically consistent with the original seismic record. Furthermore, in order to verify the compensation effect of this method under the influence of noise, this method is compared with the conventional stability factor method and amplitude gain limitation method. Figure 3 (a) is a pre-stack synthetic seismic record with Gaussian white noise of 30% added to it. Figure 3 (b) Stability factor compensated earthquake records, Figure 3 (c) Prestack seismic records compensated by amplitude gain limitation method, Figure 3 (d) Prestack earthquakes compensated by this method. By comparison, it can be seen that the stability factor method can stably compensate prestack seismic records, but the amplitude compensation of prestack seismic records is insufficient. The amplitude gain limitation method can effectively avoid distortion during deep amplitude compensation by limiting the amplitude compensation of high-frequency components, but the amplitude compensation of seismic records is also limited. At the same time, improper selection of high-frequency components can easily amplify noise. Compared with the above two methods, this method has a better improvement in amplitude and noise suppression, which verifies the noise resistance of the proposed method and its compensation effect on deep seismic records.

[0177] In summary, the present invention proposes a pre-stack Q-value attenuation compensation method based on L2 norm regularization constraint, applies the inverse Q filtering method to pre-stack seismic records by deriving the wave field extension formula, and uses L2 norm regularization for inversion compensation, thereby improving the stability and high precision of pre-stack seismic record compensation.

[0178] The present invention applies the inverse Q filtering method to the prestack by constructing a prestack Q value attenuation model, and avoids the instability of the inverse Q filtering method by using the L2 norm regularization constrained inversion method, thus overcoming the influence of insufficient compensation of the inverse Q filtering method. In the case of deep layers and low Q values, the compensation of seismic records has good stability and high precision.

Claims

1. A prestack Q-value attenuation compensation method based on L2 norm regularization constraint, characterized in that: The following steps are involved: Step 1: Preprocess the prestack seismic records and obtain the prestack Q value and incident angle information; Step 2: Based on the frequency domain wave field extension formula of plane waves in viscoelastic media, the prestack Q value formula related to the offset is constructed using the obtained prestack Q value and incident angle information, and the prestack attenuation model related to the offset is established; Step 3: Based on the pre-stack attenuation model constructed in step 2, regularization is introduced to construct the inversion objective function, and the conjugate gradient algorithm is used to iteratively update the seismic records in the initial frequency domain. Step 4: Determine whether the updated seismic record in the frequency domain meets 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 and return to step 3 to iterate again until the iteration termination condition is met.

2. The prestack Q-value attenuation compensation method based on L2 norm regularization constraint according to claim 1 is characterized in that: The step 1 comprises the following steps: Step 1.1: Using the relationship between the time, offset and Q of the prestack seismic record, the Q value of the prestack seismic record is extracted using the prestack Q value extraction method; Step 1.2: Obtain the incident angle θ of each layer in the prestack seismic record by using ray tracing technology i .

3. The prestack Q-value attenuation compensation method based on L2 norm regularization constraint according to claim 2 is characterized in that: The step 2 comprises the following steps: Step 2.1: According to the frequency domain wave field extension formula of plane waves in viscoelastic media, a propagation model of seismic waves in the medium is established. The frequency domain wave field extension formula of plane waves is expressed as: Among them, B(T, ω) represents the wave field with propagation time T and angular frequency ω, ΔT i is the propagation time increment of the seismic wave in the i-th layer, k(ω) is the frequency-dependent wave number, v r is the speed of the reference frequency, the variable j represents the imaginary unit, and e is the base of the natural logarithm; Step 2.2: Using the acquired pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of the seismic wave in the i-th layer according to Snell's law i , and ΔT i Substituting into the frequency domain wave field extension formula, we get the wave field extension: ΔT i It is expressed as: The ΔT i，0 is the two-way propagation time increment of the seismic wave vertically propagating 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, the attenuation and velocity dispersion characteristics of seismic waves are introduced into the wave field extension formula. The modified Kolsky model describes the viscoelastic properties of the medium by introducing the quality factor Q and the tuning frequency ω0, thereby more accurately simulating the propagation process of seismic waves. At this time, the wave field extension is expressed as: Where ω0 is the tuning frequency, Q i represents the quality factor of the i-th layer of medium, i = 1, 2, 3, ... n; Step 2.4: Using the time-shift property of Fourier transform, the wave field continuation formula is converted into a pre-stack attenuation model related to the offset. Specifically: The time-shift properties of Fourier transform are: Substituting formula (5) into formula (4) yields the following wave field extension formula: When the seismic wave propagates through n layers of media, the wave field extension formula becomes: Replace the propagation time ΔT with the sampling interval dt i，0 , and discretize the frequency ω into ω m , formula (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 the medium, the equivalent quality factor Q is used from the top formation to the i-th layer. e，i The formula is: Among them, T i represents the accumulated travel time of layer i. When i=n, ​​formula (8) can be expressed as follows by replacing the inter-layer Q value with the equivalent Q value: Among them, Q e，n represents the equivalent quality factor of the nth layer. At this time, the prestack positive Q attenuation formula related to the offset distance is written as:

4. The prestack Q-value attenuation compensation method based on L2 norm regularization constraint according to claim 1, characterized in that: The step 3 comprises the following steps: In step 3.1, the present invention further introduces an L2 norm regularization constraint based on the prestack Q value attenuation formula constructed in step 2 to construct an inversion objective function. The specific operations are as follows: The positive Q filtering process is represented by 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 filter operator and the unattenuated seismic record, respectively. The stabilized target functional is constructed using Tikhonov regularization: G=||LU-u|| 2 +λ||WU|| 2 (13) Among them, λ is a regularization factor, W is a unit operator, a first-order or second-order differential operator, and the analytical solution of the target functional G is expressed as: U=(L T L+λW T W) -1 L T u (14); In step 3.2, the inversion objective function G constructed in step 3.1 is iteratively updated using the conjugate gradient algorithm to solve the spectrum U of the effective frequency of the unattenuated seismic record, which specifically includes the following steps: The conjugate gradient iteration method is used to solve equation (14), all U values ​​are assigned to 0 as the initial term of the first iteration, and the data fitting error r is calculated. n , and use the data to fit the difference r n Find the direction of fastest ascent and conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows: r n =LU-d (15) Among them, F T is the transpose of the Frechet differential operator of L, To determine the coefficient of the conjugate direction, the value in the first iteration is 0, so as to iteratively update the compensated data U:

5. The prestack Q-value attenuation compensation method based on L2 norm regularization constraint according to claim 1, characterized in that: The step 4 comprises the following steps: Determine whether the iteratively updated data U meets the condition || r n || 2 ≤Tol, where Tol is the given error term value. If yes, the last iterative inversion data U is output as the final prestack compensation data. Otherwise, the last iterative inversion data U is used 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: include: A preprocessing module is used to preprocess the prestack seismic records and obtain the prestack Q value and incident angle information; The pre-stack Q value filter construction module, based on the frequency domain wave field extension formula of plane waves in viscoelastic media, uses the pre-stack Q value and incident angle information obtained by the pre-processing module to construct a pre-stack Q value filter related to the offset, and establishes a pre-stack attenuation model related to the offset; The inversion objective function building module introduces regularization to build the inversion objective function based on the prestack Q-value filter built by the prestack Q-value filter building module, and uses the conjugate gradient algorithm to iteratively update the seismic records in the initial frequency domain; The iteration termination condition judgment module is used to judge whether the frequency domain seismic record updated by the inversion objective function construction module meets the iteration termination condition. If so, the updated frequency domain seismic record is output as the final inversion result. Otherwise, the updated frequency domain seismic record is used as the new initial model and returned to the inversion objective function construction module for another iterative update until the iteration termination condition is met.

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

8. The prestack Q-value attenuation compensation device based on L2 norm regularization constraint according to claim 7, characterized in that: The pre-stack Q-value filter building block implementation includes the following steps: Step 2.1: According to the frequency domain wave field extension formula of plane waves in viscoelastic media, a propagation model of seismic waves in the medium is established. The frequency domain wave field extension formula of plane waves is expressed as: Among them, B(T, ω) represents the wave field with propagation time T and angular frequency ω, ΔT i is the propagation time increment of the seismic wave in the i-th layer, k(ω) is the frequency-dependent wave number, v r is the speed of the reference frequency, the variable j represents the imaginary unit, and e is the base of the natural logarithm; Step 2.2: Using the acquired pre-stack Q value and incident angle information, calculate the propagation time increment ΔT of the seismic wave in the i-th layer according to Snell's law i , and ΔT i Substituting into the frequency domain wave field extension formula, we get the wave field extension: ΔT i It is expressed as: The ΔT i，0 is the two-way propagation time increment of the seismic wave vertically propagating 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, the attenuation and velocity dispersion characteristics of seismic waves are introduced into the wave field extension formula. The modified Kolsky model describes the viscoelastic properties of the medium by introducing the quality factor Q and the tuning frequency ω0, thereby more accurately simulating the propagation process of seismic waves. At this time, the wave field extension is expressed as: Where ω0 is the tuning frequency, Q i represents the quality factor of the i-th layer of medium, i = 1, 2, 3, ... n; Step 2.4: Using the time-shift property of Fourier transform, the wave field extension formula is converted into a pre-stack Q value filtering formula related to the offset. Specifically: The time-shift properties of Fourier transform are: Substituting formula (5) into formula (4) yields the following wave field extension formula: When the seismic wave propagates through n layers of media, the wave field extension formula becomes: Replace the propagation time ΔT with the sampling interval dt i，0 , and discretize the frequency ω into ω m , formula (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 the medium, the equivalent quality factor Q is used from the top formation to the i-th layer. e，i The formula is: Among them, T i represents the accumulated travel time of layer i. When i=n, ​​formula (8) can be expressed as follows by replacing the inter-layer Q value with the equivalent Q value: Among them, Q e，n represents the equivalent quality factor of the nth layer. At this time, the prestack positive Q filter formula related to the offset distance is written as:

9. The prestack Q-value attenuation compensation device based on L2 norm regularization constraint according to claim 1, characterized in that: The inversion objective function building module implementation includes the following steps: In step 3.1, the present invention further introduces an L2 norm regularization constraint based on the prestack Q value filter constructed in the prestack Q value filter construction module to construct an inversion objective function. The specific operations are as follows: The positive Q filtering process is represented by 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 filter operator and the unattenuated seismic record, respectively. The stabilized target functional is constructed using Tikhonov regularization: G=||LU-u|| 2 +λ||WU|| 2 (13) Among them, λ is a regularization factor, W is a unit operator, a first-order or second-order differential operator, and the analytical solution of the target functional G is expressed as: U=(L T L+λW T W) -1 L T u (14); In step 3.2, the inversion objective function G constructed in step 3.1 is iteratively updated using the conjugate gradient algorithm to solve the spectrum U of the effective frequency of the unattenuated seismic record, which specifically includes the following steps: The conjugate gradient iteration method is used to solve equation (14), all U values ​​are assigned to 0 as the initial term of the first iteration, and the data fitting error r is calculated. n , and use the data to fit the difference r n Find the direction of fastest ascent and conjugate direction Then calculate the conjugate gradient iteration step size The formula is as follows: r n =LU-d (15) Among them, F T is the transpose of the Frechet differential operator of L, To determine the coefficient of the conjugate direction, the value in the first iteration is 0, so as to iteratively update the compensated data U:

10. The prestack Q-value attenuation compensation method based on L2 norm regularization constraint according to claim 7, characterized in that: The implementation of the iteration termination condition judgment module includes the following steps: Determine whether the iteratively updated data U meets the condition || r n || 2 ≤Tol, where Tol is the given error term value. If yes, the last iterative inversion data U is output as the final prestack compensation data. Otherwise, the last iterative inversion data U is used as the initial value and the module is transferred to the inversion objective function construction module.

Images