A viscoelastic medium compensation migration method, system and electronic device for multi-component seismic records

Abstract

The present invention discloses a viscoelastic medium compensation migration method, system and electronic device for multi-component seismic records. The method includes: obtaining the phase-delay amplitude-decayed seismic source wavefield corresponding to the shot point by using the corresponding shot point position information and known conditions; performing reverse wavefield continuation on the preprocessed multi-component reflected wave seismic records according to the phase dispersion amplitude compensation wave equation to obtain the phase-delay amplitude-decayed geophone wavefield and the phase-delay amplitude-constant geophone wavefield, and calculating the compensation coefficient and the geophone compensation wavefield respectively; calculating the P-wave and S-wave separated seismic source wavefield and the geophone compensation wavefield; calculating the PP component, PS component, SS component and SP component of the single-shot migration result respectively; superimposing all the single-shot compensation migration results to obtain the underground compensation profile, calculating the multi-shot migration superimposition result, and obtaining the multi-component compensation migration record. Compared with other viscoelastic medium migration techniques, the calculation process of the present invention is more stable.

Description

Technical Field

[0001] The present invention relates to the technical field of reflection wave seismic exploration, and in particular to a viscoelastic medium compensation migration method, system and electronic device for multi-component seismic records. Background Art

[0002] The actual underground medium is closer to non-perfect elasticity. The multi-component seismic records collected in seismic exploration are affected by absorption attenuation, resulting in amplitude attenuation and phase distortion. If only the elastic wave migration method is used, there will be problems of amplitude and phase distortion in the imaging results. Migration considering the viscoelasticity of the medium can compensate the amplitude of seismic waves, correct the phase, and obtain higher-resolution results.

[0003] In the prior art, Zhang et al. (2010) added a pseudo-difference operator to the visco-acoustic wave equation of the constant Q model to perform visco-acoustic reverse time migration imaging to eliminate amplitude attenuation and velocity dispersion. Since the seismic energy increases with the increase of frequency during the wave field reverse continuation process, a regularization step was introduced to avoid numerical instability, thus realizing the reverse propagation of the viscous wave field. Suh et al. (2012) developed this operator to the VTI medium. Bai et al. (2013) adopted a similar method for attenuation compensation during the reverse time migration imaging process, using a visco-acoustic wave equation without memory variables. Such algorithms can perform a certain degree of compensation in terms of energy amplitude, but further research is needed for the problems of multi-component and phase correction. Zhu et al. (2014) proposed a new decoupled fractional Laplacian constant Q visco-acoustic wave equation, which divides the attenuation term into two terms: amplitude attenuation and velocity dispersion, and applied it to the reverse time migration imaging of visco-acoustic media with Q compensation, and measurements were carried out on theoretical data and actual data. For the stability of the wave equation continuation in the Q compensation migration method, Fletcher et al. (2012) proposed to design filtering operators independent of both amplitude and phase, and through two acoustic wave field continuations, before Q compensation imaging, the filtering operators are used for the source and receiver wave fields. The filtering method is based on the travel time attenuation of the propagation path, so it may not be applicable to the cases of multiple waves and drastic changes in the Q model. Zhao et al. (2018) used the excitation amplitude imaging condition for viscoelastic medium migration imaging, and some scholars realized the migration imaging problem of viscous media by modifying the wave equation (Yang et al., 2018). Sun et al. (2018) compensated for attenuation in the least squares reverse time migration. The least squares reverse time migration can obtain high-precision imaging results and will not have stability problems (Guo et al., 2018), but it inevitably brings a greater computational burden. At present, the migration for viscous media mainly focuses on single-component data, only considering the propagation process of longitudinal waves, or the viscoelastic model used cannot handle complex absorption attenuation regions. There is a need to improve the high-precision viscoelastic medium compensation migration method for multi-component data.

[0004] In the prior art, for multi-component seismic records, only elastic media are mainly considered, but it is impossible to cope with areas with strong absorption attenuation. The existing viscous medium compensation imaging technology mainly compensates and migrates based on viscoacoustic media, unable to make full use of the properties of elastic waves. Only the longitudinal wave data component is used for imaging, and multi-component data cannot be fully utilized. Although some public materials consider the viscoelasticity of the medium and construct wave field imaging based on the cross-correlation imaging condition, there are unstable phenomena in the application process, and it is impossible to cope with areas with relatively complex lateral variations. Summary of the Invention

[0005] To solve the above technical problems, the present invention proposes a viscoelastic medium compensation migration method for multi-component seismic records. Based on a stable imaging condition, stable forward and reverse wave fields are calculated and constructed respectively to realize multi-component data compensation viscoelastic migration imaging, and four component migration results can be obtained simultaneously, which can improve the imaging resolution.

[0006] On the one hand, to achieve the above object, the present invention provides a viscoelastic medium compensation migration method for multi-component seismic records, including:

[0007] S1. Obtain multi-component reflection wave seismic records of multiple shots. For single-shot data, use the corresponding shot point position information and known conditions to obtain the phase-delay amplitude attenuation source wave field of the corresponding shot point, where the known conditions include the longitudinal wave velocity model, transverse wave velocity model, longitudinal wave quality factor, and transverse wave quality factor obtained after velocity modeling is completed in the target area;

[0008] S2. Use the preprocessed multi-component reflection wave seismic records to perform reverse wave field continuation according to the phase dispersion amplitude compensation wave equation to obtain the phase-delay amplitude attenuation geophone wave field and the phase-delay amplitude invariant geophone wave field respectively;

[0009] S3. Calculate the compensation coefficient and the geophone compensation wave field respectively through the phase-delay amplitude attenuation geophone wave field and the phase-delay amplitude invariant geophone wave field;

[0010] S4. Calculate the source wave field with longitudinal and transverse waves separated and the geophone compensation wave field with longitudinal and transverse waves separated based on the Helmholtz decomposition method to obtain the longitudinal wave component of the source wave field, the transverse wave component of the source wave field, the longitudinal wave component of the geophone compensation wave field, and the transverse wave component of the geophone compensation wave field;

[0011] S5. Under the stable imaging condition, use the longitudinal wave component of the source wave field, the transverse wave component of the source wave field, the longitudinal wave component of the geophone compensation wave field, and the transverse wave component of the geophone compensation wave field at N moments to calculate the single-shot migration results of the PP component, PS component, SS component, and SP component respectively;

[0012] S6. Repeat S2 - S5, stack all the single-shot compensation migration results to obtain the subsurface compensation profile, calculate the multi-shot migration stack result, and obtain the multi-component compensation migration record.

[0013] Preferably, obtaining the phase-delay amplitude attenuated source wavefield corresponding to the shot point includes:

[0014] Using the shot point position information and the known conditions, calculate the source wavefield, and substitute the source wavelet and the single-shot acquisition geometry into the phase-delay amplitude compensation wave equation for wavefield extrapolation to obtain the phase-delay amplitude attenuated source wavefield corresponding to the shot point; where the method for performing the wavefield extrapolation is:

[0015]

[0016] In equations (1) and (2), A xx , A zz , A xz are the stress components of the wavefield respectively, B x , B z are the horizontal velocity component and the vertical velocity component of the wavefield respectively, S(x, z, t) is the loaded source function, ρ is the known density model, V p , V s , Q p , Q s correspond to the known P-wave velocity model, S-wave velocity model, P-wave quality factor model, and S-wave quality factor model respectively, t is time, x is distance, z is depth, is the Laplace operator, e p , f p , d p , d s , e s , f s are intermediate variables respectively, and ω0 is the dominant frequency of the wavelet.

[0017] Preferably, in the wavefield extrapolation formula, the integer-order time partial derivative is solved using second-order finite differences, and the fractional-order time partial derivative and the spatial partial derivative terms are solved using Fourier transforms. The specific discrete iteration formula is:

[0018]

[0019] In equation (3), Δt is the calculation time interval, A xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1They are the wave field values of each component at the time of (n + 1)Δt; A xx n and A zz n and A xz n and B x n and B z n They are the wave field values of each component at the time of nΔt; B x n-1 and B z n-1 They are the wave field values of each component at the time of (n - 1)Δt; e p and f p and d p and d s and e s and f s are intermediate variables, F is the forward Fourier transform, and F -1 is the inverse Fourier transform; k x and k z are the wave numbers in each direction, is the average of f p over the entire model, i is the imaginary unit, and S n is the source amplitude value at the time of nΔt; Using the wave field values of each component at the times of T = nΔt and T = (n - 1)Δt (A xx n and A zz n and A xz n and B x n and B z n and B x n-1 and B z n-1 ), the wave field values at the time of T = (n + 1)Δt (A xx n+1 and A zz n+1 and A xz n+1 and B x n+1 and B z n+1 ) are iteratively calculated;

[0020] Substitute the source function S(x, z, t), calculate a total of N times, start loading the source function from T = 0, and finally obtain the source wave field values at N times, including B x and B z and A xx and Azz , component A xz forms the phase-delay amplitude-decay seismic source wavefield W corresponding to the shot point shot (x, z).

[0021] Preferably, the method for obtaining the phase-delay amplitude-decay geophone wavefield is as follows:

[0022] Based on the vertical velocity component and the horizontal velocity component of the seismic record, inverse-time iterative calculation is performed, specifically:

[0023]

[0024] In equation (4), R x (x, z, t) and R z (x, z, t) are the x and z components of the seismic record respectively;

[0025] The specific calculation iteration method is as follows:

[0026]

[0027] In equation (5), R x n (x, z) and R z n (x, z) are the x component and the z component of the seismic record at nΔt, Δt is the calculation time interval. Using the wavefield values of each component at T = nΔt and T = (n + 1)Δt (A xx n , A zz n , A xz n , B x n , B z n , B x n+1 , B z n +1 ), and the seismic record values R x n (x, z) and R z n (x, z) at T = nΔt, the wavefield values at T = (n - 1)Δt are obtained by iterative calculation (A xx n-1 , A zz n-1 , A xz n-1 , B x n-1 , B z n-1 );

[0028] Calculate the wave field at the geophone points at N moments, and substitute it into the vertical and horizontal components R of the seismic record x (x, z, t), R z (x, z, t), start loading the seismic record from T = N for iteration, and finally obtain the wave field values at the geophone points at N moments, including B x 、B z 、A xx 、A zz 、A xz components, that is, form the phase delay amplitude attenuation wave field W of the corresponding shot point rec1 (x, z).

[0029] Preferably, the method for obtaining the phase delay amplitude invariant wave field W rec2 (x, z) is as follows:

[0030]

[0031] The specific calculation iteration method is as follows:

[0032]

[0033] In equations (6) and (7), use the wave field values of each component at T = nΔt and T = (n + 1)Δt (A xx n 、A zz n 、A xz n 、B x n 、B z n 、B x n+1 、B z n+1 ), the seismic record value R at T = nΔt x n (x, z)、R z n (x, z), and iteratively calculate the wave field values at T = (n - 1)Δt (A xx n-1 、A zz n-1 、A xz n-1 、B x n-1 、B z n-1 );

[0034] Calculate the wave field at the geophone points at N moments, and substitute it into the vertical and horizontal components R of the seismic record x (x, z, t), Rz (x, z, t), starting from T = N, load the seismic records for iteration, and finally obtain the wavefield values of the geophone points at N moments, including B x 、B z 、A xx 、A zz 、A xz The components, namely, constitute the phase-delay amplitude-invariant geophone point wavefield W rec2 (x, z).

[0035] Preferably, the method for calculating the compensation coefficient is as follows:

[0036]

[0037] In the formula, ε is white noise, is the compensation coefficient, W rec1 (x, z) is the phase-delay amplitude-decaying geophone point wavefield, and W rec2 (x, z) is the phase-delay amplitude-invariant geophone point wavefield;

[0038] The method for calculating the compensated wavefield of the geophone point is as follows:

[0039]

[0040] In the formula, is the compensated wavefield of the geophone point.

[0041] Preferably, the method for calculating the PP component, PS component, SS component, and SP component of the single-shot migration result is as follows:

[0042]

[0043] In the formula, ε is white noise, is the compensated wavefield of the geophone point, W shot (x, z) is the phase-delay amplitude-decaying source wavefield, I(x, z) is the component of the single-shot migration result, and N is the number of time intervals.

[0044] Preferably, the method for calculating the multi-shot migration stacking result is as follows:

[0045]

[0046] In the formula, I stack (x, z) is the stacking result, and M is the total number of calculated shots.

[0047] On the other hand, to achieve the above object, the present invention also provides a viscoelastic medium compensation migration system for multi-component seismic records, which is applied to the viscoelastic medium compensation migration method for multi-component seismic records, including:

[0048] Phase delay amplitude attenuation seismic source wavefield acquisition module: For single-shot data, using the corresponding shot point position information and known conditions, obtain the phase delay amplitude attenuation seismic source wavefield of the corresponding shot point, where the known conditions include multi-component reflection wave seismic records of multiple shots, the P-wave velocity model, S-wave velocity model, P-wave quality factor, and S-wave quality factor obtained after velocity modeling of the target area;

[0049] Reverse wavefield continuation processing module: Use the preprocessed multi-component reflection wave seismic records to perform reverse wavefield continuation according to the phase dispersion amplitude compensation wave equation, and obtain the phase delay amplitude attenuation geophone wavefield and the phase delay amplitude invariant geophone wavefield respectively;

[0050] Compensation coefficient and geophone compensation wavefield calculation module: Calculate the compensation coefficient and the geophone compensation wavefield respectively through the phase delay amplitude attenuation geophone wavefield and the phase delay amplitude invariant geophone wavefield;

[0051] P-S wave separation module: Calculate the seismic source wavefield of P-S wave separation and the geophone compensation wavefield of P-S wave separation based on the Helmholtz decomposition method, and obtain the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield;

[0052] Single-shot migration result component acquisition module: Under stable imaging conditions, use the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield at N time moments to calculate the PP component, PS component, SS component, and SP component of the single-shot migration result respectively;

[0053] Multi-component compensation migration output module: Superimpose all single-shot compensation migration results to obtain an underground compensation profile, calculate the multi-shot migration stack result, and obtain a multi-component compensation migration record.

[0054] This embodiment also provides an electronic device, including a processor, a memory, and a program or instruction stored on the memory and executable on the processor. When the program or instruction is executed by the processor, it implements the steps of the viscoelastic medium compensation migration method for multi-component seismic records.

[0055] Compared with the prior art, the present invention has the following advantages and technical effects:

[0056] Compared with the existing multi-component migration imaging technology, the present invention simultaneously considers the medium absorption attenuation characteristics and elastic characteristics. After compensating the energy, the imaging position under the gas cloud area is accurate and the resolution is higher; compared with other viscoelastic medium migration technologies, the calculation process of the present invention is more stable. Description of the Drawings

[0057] The accompanying drawings, which form a part of this application, are used to provide a further understanding of this application. The schematic embodiments of this application and their descriptions are used to explain this application and do not constitute an improper limitation to this application. In the accompanying drawings:

[0058] Figure 1 It is a flowchart of a viscoelastic medium compensation migration method for multi-component seismic records according to an embodiment of the present invention;

[0059] Figure 2 It is a schematic diagram of the longitudinal and transverse wave velocities and Q-value parameters of the model according to an embodiment of the present invention. Among them, (a) is the longitudinal wave velocity model, (b) is the transverse wave velocity model, (c) is the longitudinal wave quality factor model, and (d) is the transverse wave quality factor model;

[0060] Figure 3 It is a schematic diagram of a typical single-shot seismic record according to an embodiment of the present invention. Among them, (a) is the horizontal velocity component, and (b) is the vertical velocity component;

[0061] Figure 4 It is a schematic diagram of the source wave field at the 500 ms moment of a single shot according to an embodiment of the present invention. Among them, (a) is the B x component, and (b) is the B z component;

[0062] Figure 5 It is a schematic diagram of the wave field of the phase delay amplitude attenuation geophone at the 500 ms moment of a single shot according to an embodiment of the present invention. Among them, (a) is the B x component, and (b) is the B z component;

[0063] Figure 6 It is a schematic diagram of the wave field of the phase delay amplitude invariant geophone at the 500 ms moment of a single shot according to an embodiment of the present invention. Among them, (a) is the B x component, and (b) is the B z component;

[0064] Figure 7 It is a schematic diagram of the longitudinal and transverse wave components of the source wave field and the longitudinal and transverse wave components of the compensated geophone wave field at the 500 ms moment of a single shot according to an embodiment of the present invention. Among them, (a) is the W shot (x,z) P component, and (b) is the W shot (x,z) S component, (c) is the component, and (d) is the component;

[0065] Figure 8 It is a schematic diagram of the single-shot migration result according to an embodiment of the present invention. Among them, (a) is the PP component, (b) is the PS component, (c) is the SP component, and (d) is the SS component;

[0066] Figure 9 Schematic diagram of multi-component stacking profile of an embodiment of the present invention, where (a) is the PP component, (b) is the PS component, (c) is the SP component, and (d) is the SS component. Detailed implementation manners

[0067] It should be noted that, without conflict, the embodiments in this application and the features in the embodiments can be combined with each other. The following will detail this application with reference to the accompanying drawings and in combination with the embodiments.

[0068] It should be noted that the steps shown in the flowchart of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although the logical order is shown in the flowchart, in some cases, the steps shown or described can be executed in a different order than here.

[0069] The present invention proposes a viscoelastic medium compensation migration method for multi-component seismic records, as Figure 1 , including:

[0070] S1. For single-shot data, using the corresponding shot point position information and known conditions, obtain the source wave field with phase delay and amplitude attenuation of the corresponding shot point, where the known conditions include multi-component reflection wave seismic records of multiple shots, the longitudinal wave velocity model, the transverse wave velocity model, the longitudinal wave quality factor, and the transverse wave quality factor obtained after velocity modeling is completed in the target area;

[0071] S2. Using the preprocessed multi-component reflection wave seismic records, perform reverse wave field continuation according to the phase dispersion amplitude compensation wave equation to obtain the phase delay amplitude attenuation geophone wave field and the phase delay amplitude invariant geophone wave field respectively;

[0072] S3. Through the phase delay amplitude attenuation geophone wave field and the phase delay amplitude invariant geophone wave field, calculate the compensation coefficient and the geophone compensation wave field respectively;

[0073] S4. Based on the Helmholtz decomposition method, calculate the source wave field with longitudinal and transverse waves separated and the geophone compensation wave field with longitudinal and transverse waves separated to obtain the longitudinal wave component of the source wave field, the transverse wave component of the source wave field, the longitudinal wave component of the geophone compensation wave field, and the transverse wave component of the geophone compensation wave field;

[0074] S5. Under stable imaging conditions, use the longitudinal wave component of the source wave field, the transverse wave component of the source wave field, the longitudinal wave component of the geophone compensation wave field, and the transverse wave component of the geophone compensation wave field at N moments to calculate the single-shot migration results of the PP component, the PS component, the SS component, and the SP component respectively;

[0075] S6. Repeat S2 - S5, stack all the single - shot compensation migration results to obtain the subsurface compensation profile, calculate the multi - shot migration stack result, and obtain the multi - component compensation migration record.

[0076] Compared with the existing multi - component migration imaging technology, this embodiment simultaneously considers the characteristics of medium absorption attenuation and elastic characteristics. After compensating the energy, the imaging position under the gas cloud area is accurate and the resolution is higher. Compared with other viscous medium migration technologies, the calculation process of this embodiment is more stable.

[0077] Further, obtaining the phase - delay amplitude - attenuated seismic source wavefield corresponding to the shot point includes:

[0078] Using the shot - point position information and the previously known conditions, calculate the seismic source wavefield, and substitute the seismic source wavelet and the single - shot observation system into the phase - delay amplitude compensation wave equation for wavefield extrapolation to obtain the phase - delay amplitude - attenuated seismic source wavefield corresponding to the shot point.

[0079] Specifically, in this embodiment, the known data includes the multi - component reflection wave seismic records (vertical component and horizontal component) of multiple shots, the P - wave velocity model, S - wave velocity model, P - wave quality factor, and S - wave quality factor obtained after velocity modeling of the target area.

[0080] For the calculated single - shot data, use the shot - point position information and the previously known conditions to calculate the seismic source wavefield. Substitute the seismic source wavelet and the single - shot observation system into the phase - delay amplitude compensation wave equation for wavefield extrapolation to obtain the phase - delay amplitude - attenuated seismic source wavefield corresponding to the shot point. The equations for wavefield extrapolation are shown in Formulas (1) - (2).

[0081]

[0082]

[0083] In the formula, A xx 、A zz 、A xz are the stress components of the wavefield respectively, B x 、B z are the horizontal velocity component and vertical velocity component of the wavefield respectively (3 stress components and 2 velocity components constitute the phase - delay amplitude - attenuated seismic source wavefield W shot (x, z)), S(x, z, t) is the loaded seismic source function, ρ is the known density model, V p 、V s 、Q p 、Q s correspond to the known P - wave velocity model, S - wave velocity model, P - wave quality factor model, and S - wave quality factor model respectively, t is time, x is distance, z is depth, is the Laplace operator, ep , f p , d p , d s , e s , f s are intermediate variables, and ω0 is the dominant frequency of the wavelet.

[0084] Furthermore, in the wavefield extrapolation formula, the integer-order time partial derivative is solved using second-order finite differences, and the fractional-order time partial derivative and the spatial partial derivative terms are solved using Fourier transforms. The specific discrete iteration formula is:

[0085]

[0086] In Equation (3), Δt is the calculation time interval, and A xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1 are the wavefield values of each component at time (n + 1)Δt, and A xx n , A zz n , A xz n , B x n , B z n are the wavefield values of each component at time nΔt; B x n-1 , B z n-1 are the wavefield values of each component at time (n - 1)Δt; e p , f p , d p , d s , e s , f s are intermediate variables, F is the forward Fourier transform, and F -1 is the inverse Fourier transform; k x , k z are the wave numbers in each direction, is the average of f p for the entire model, i is the imaginary unit, S n is the source amplitude value at time nΔt. Using the wavefield values of each component at times T = nΔt and T = (n - 1)Δt (A xx n , A zz n , A xzn , B x n , B z n , B x n-1 , B z n-1 ) The wave field value at time T = (n + 1)Δt is obtained by iterative calculation (A xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1 );

[0087] Substitute the source function S(x, z, t). A total of N time instants are calculated. Starting from T = 0, the source function is loaded, and finally the source wave field values at N time instants are obtained, including B x , B z , A xx , A zz , A xz The components, namely, constitute the phase delay amplitude attenuation source wave field W shot (x, z) of the corresponding shot point.

[0088] Further, the method for obtaining the phase delay amplitude attenuation receiver wave field is as follows:

[0089] Use the multi-component seismic record after conventional noise suppression, and obtain the phase delay amplitude attenuation receiver wave field W rec1 (x, z) by backpropagating the wave field according to the phase dispersion amplitude compensation wave equation. The calculation is shown in formula (4). Substitute the vertical velocity component and horizontal velocity component of the seismic record into the equation and perform inverse time iterative calculation:

[0090]

[0091] The specific iterative calculation method is shown in formula (5):

[0092]

[0093] In formula (4), R x (x, z, t), R z (x, z, t) are the x and z components of the seismic record respectively; the meanings of other components are the same as those in (1) and (2); in formula (5), R x n (x, z), R z n(x, z) are the x - component and z - component of the seismic record at nΔt time, and the meanings of other components are the same as those in Equation (3); using the wave - field values of each component at T = nΔt and T=(n + 1)Δt (A xx n 、A zz n 、A xz n 、B x n 、B z n 、B x n +1 、B z n+1 ), the seismic record values R x n (x, z)、R z n (x, z) at T = nΔt, and iteratively calculate the wave - field values at T=(n - 1)Δt (A xx n-1 、A zz n-1 、A xz n-1 、B x n-1 、B z n-1 );

[0094] In this embodiment, it is assumed that the wave - field of the geophone points at N times is calculated, and the vertical and horizontal components R x (x, z, t)、R z (x, z, t) of the seismic record are substituted. Starting from T = N, the seismic record is loaded for iteration, and finally the wave - field values of the geophone points at N times are obtained, including B x 、B z 、A xx 、A zz 、A xz components, that is, the phase - delay amplitude - attenuation geophone - point wave - field W rec1 (x, z) corresponding to the shot point is formed.

[0095] Furthermore, the method for obtaining the phase - delay amplitude - invariant geophone - point wave - field W rec2 (x, z) is as follows:

[0096] Using the pre - processed multi - component seismic record, the phase - delay amplitude - invariant geophone - point wave - field is obtained by reverse wave - field continuation according to the phase - dispersion amplitude - compensation wave equation:

[0097]

[0098] The specific iterative calculation method is shown in Equation (7):

[0099]

[0100] The components in equations (6) and (7) have the same meanings as those in equations (1), (2), and (3). Using the wave field values of each component at times T = nΔt and T = (n + 1)Δt (A xx n 、A zz n 、A xz n 、B x n 、B z n 、B x n+1 、B z n+1 ), and the seismic record values R x n (x,z)、R z n (x,z) at time T = nΔt, the wave field values at time T = (n - 1)Δt (A xx n-1 、A zz n-1 、A xz n-1 、B x n-1 、B z n-1 ) are obtained through iterative calculation;

[0101] Calculate the wave fields at the geophone points for N time instances, substitute them into the vertical and horizontal components R x (x,z,t)、R z (x,z,t) of the seismic record. Starting from T = N, load the seismic record for iteration, and finally obtain the wave field values at the geophone points for N time instances, including B x 、B z 、A xx 、A zz 、A xz components, which constitute the phase-delay amplitude-invariant geophone point wave field W rec2 (x,z) corresponding to the shot point.

[0102] Furthermore, the method for calculating the compensation coefficient is as follows:

[0103]

[0104] where ε is white noise, is the compensation coefficient, W rec1 (x,z) is the phase-delay amplitude-decaying geophone point wave field, W rec2(x, z) is the wave field of the phase-delay amplitude-invariant detection point;

[0105] The method for calculating the compensated wave field of the detection point is as follows:

[0106]

[0107] In the formula, is the compensated wave field of the detection point.

[0108] Furthermore, the method for calculating the source wave field separated into longitudinal and transverse waves and the compensated wave field of the detection point separated into longitudinal and transverse waves based on the Helmholtz decomposition method is as follows:

[0109]

[0110] Obtain the longitudinal wave component and transverse wave component (W shot (x, z) P 、W shot (x, z) S ) of the source wave field, and compensate the longitudinal wave component and transverse wave component of the detection point wave field

[0111] Furthermore, calculate the single-shot imaging result based on the stable imaging condition, specifically:

[0112] The imaging condition is shown in formula (11). Substitute the longitudinal and transverse wave components of the source wave field and the compensated detection point wave field at N moments into the formula respectively to obtain the single-shot migration results of the PP component, PS component, SS component, and SP component (I(x, z) PP 、I(x, z) PS 、I(x, z) SP 、I(x, z) SS ):

[0113]

[0114] In the formula, ε is white noise, is the compensated wave field of the detection point, W shot (x, z) is the source wave field with phase-delay amplitude attenuation, I(x, z) is the component of the single-shot migration result, and N is the total number of time intervals.

[0115] Furthermore, stack all the single-shot compensated migration results to obtain the underground compensated profile:

[0116] Repeat S2 - S5 for the seismic records of multiple shots, and then stack the migration results of the entire observed multiple shots (formula 12) to obtain the compensated multi-component stacked profile I stack (x, z) PP 、I stack (x, z) PS 、I stack(x,z) SP 、I stack (x,z) SS

[0117]

[0118] where I stack (x,z) is a multi-component superposition profile, and M is the total number of shots.

[0119] This embodiment also provides a viscoelastic medium compensation migration system for multi-component seismic records, which is applied to the viscoelastic medium compensation migration method for multi-component seismic records, and includes:

[0120] Phase-delay amplitude-decay seismic source wavefield acquisition module: For single-shot data, using the corresponding shot-point position information and known conditions, obtain the phase-delay amplitude-decay seismic source wavefield of the corresponding shot point, where the known conditions include multi-component reflection wave seismic records of multiple shots, the P-wave velocity model, S-wave velocity model, P-wave quality factor, and S-wave quality factor obtained after velocity modeling in the target area;

[0121] Reverse wavefield continuation processing module: Use the preprocessed multi-component reflection wave seismic records to perform reverse wavefield continuation according to the phase dispersion amplitude compensation wave equation, and respectively obtain the phase-delay amplitude-decay geophone wavefield and the phase-delay amplitude-invariant geophone wavefield;

[0122] Compensation coefficient and geophone compensation wavefield calculation module: Calculate the compensation coefficient and the geophone compensation wavefield respectively through the phase-delay amplitude-decay geophone wavefield and the phase-delay amplitude-invariant geophone wavefield;

[0123] P-S wave separation module: Based on the Helmholtz decomposition method, calculate the seismic source wavefield of P-S wave separation and the geophone compensation wavefield of P-S wave separation to obtain the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield;

[0124] Single-shot migration result component acquisition module: Under stable imaging conditions, use the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield at N moments to calculate the PP component, PS component, SS component, and SP component of the single-shot migration result respectively;

[0125] Multi-component compensation migration output module: Superimpose all single-shot compensation migration results to obtain an underground compensation profile, calculate the multi-shot migration stacking result, and obtain a multi-component compensation migration record.

[0126] This embodiment also provides an electronic device, including a processor, a memory, and a program or instruction stored on the memory and executable on the processor. When the program or instruction is executed by the processor, the steps of the viscoelastic medium compensation migration method for multi-component seismic records are implemented.

[0127] To more clearly express the technical solution of the present invention, specific embodiments are provided below for scheme introduction:

[0128] Select the multi-component seismic record of a classic gas cloud area model for migration imaging.

[0129] S2.1. Basic parameters of the model:

[0130] It is known that the basic parameters of this model are shown in Figure 2 (a)- Figure 2 (d). The longitudinal wave velocity, transverse wave velocity, longitudinal wave quality factor, and transverse wave quality factor. The model length is 5000m, the depth is 2000m, and the model grid size is 12.5m in both length and depth.

[0131] S2.2. Parameters of the multi-component seismic record:

[0132] It is known that the multi-component seismic record based on this model has a total of 40 shots. The shot point depth is 0, the starting position of the shot point is 0, and the shot interval is 125m; each shot is received by 398 geophones. The receiving point depth is 0, the trace interval is 12.5m, the geophone coverage ranges from 0 to 4962.5m, the time sampling interval is 1ms, and the record length is 3.5s. The source is a 25Hz Ricker wavelet. A typical single-shot multi-component seismic record (the 21st shot) is shown in Figure 3 (a)- Figure 3 (b).

[0133] S2.3. Calculate the phase-delay amplitude-decay source wavefield:

[0134] Using the known model-related parameters in S2.1, the Ricker wavelet source, and the recurrence formula (3), the source wavefield at N moments is calculated. Among them, the B x 、B z components at 500ms are shown in Figure 4 (a)- Figure 4 (b).

[0135] S2.4. Calculate the phase-delay amplitude-decay geophone wavefield:

[0136] Using the single-shot seismic record in S2.2, the known model-related parameters in S2.1, and the recurrence formula (5), the phase-delay amplitude-decay geophone wavefield at N moments is calculated. Among them, the B x 、B z components at 500ms are shown inFigure 5 (a)- Figure 5 (as shown in (b)).

[0137] S2.5. Calculate the phase-delay amplitude-invariant geophone-point wavefield:

[0138] Using the single-shot seismic record in S2.2, the known model-related parameters in S2.1, and the recurrence formula (7), the phase-delay amplitude-invariant geophone-point wavefield at N time instances is calculated. Among them, the B x and B z components are as shown in Figure 6 (a)- Figure 6 (b).

[0139] S2.6. Calculate the geophone-point wavefield amplitude compensation coefficient and compensate the geophone-point wavefield:

[0140] (1) Using the two geophone-point wavefields W rec1 (x,z) and W rec2 (x,z) obtained in S2.4 and S2.5, the compensation coefficients are calculated through formula (8).

[0141] (2) Using the phase-delay amplitude-invariant geophone-point wavefield W rec2 (x,z) and the compensation coefficients, substitute them into formula (9) to calculate the compensated geophone-point wavefield at N time instances.

[0142] S2.7. Calculate the compensated wavefields of the longitudinal and transverse wave separations of the source and geophones based on Helmholtz decomposition:

[0143] Based on formula (10), for each component of the source wavefield and the compensated geophone-point wavefield at N time instances, the longitudinal wave component, transverse wave component of the source wavefield, and the longitudinal wave component and transverse wave component of the compensated geophone-point wavefield are calculated Among them, the components at 500 ms are as shown in Figure 7 (a)- Figure 7 (d).

[0144] S2.8. Calculate the single-shot imaging result based on the stable imaging condition:

[0145] Substitute the longitudinal and transverse wave components of the source wavefield and the compensated geophone-point wavefield at N time instances into formula (11) respectively to obtain the single-shot migration results of the PP, PS, SS, and SP components, as shown in Figure 8 (a)- Figure 8 (d).

[0146] S2.9. Calculate the multi-shot migration stacking result to obtain the multi-component compensated migration record:

[0147] This model has 40-shot seismic records. The calculation steps from S2.2 to S2.8 are repeated for each shot to obtain the migration imaging results for the 40 shots. Finally, the multi-component results of the entire stacked section are obtained by performing stacking processing using Equation (12), as shown in Figure 9 (a)- Figure 9 (d).

[0148] The above is only a preferred specific embodiment of the present application. However, the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present application should be covered by the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A viscoelastic medium compensation migration method for multi-component seismic records, characterized in that, Including: S1. Obtain multi-component reflection wave seismic records of multiple shots. For single-shot data, use the corresponding shot point position information and known conditions to obtain the phase-delay amplitude-decayed seismic source wavefield of the corresponding shot point, where the known conditions include the P-wave velocity model, S-wave velocity model, P-wave quality factor, and S-wave quality factor obtained after velocity modeling of the target area; S2. Use the preprocessed multi-component reflection wave seismic records to perform reverse wavefield continuation according to the phase dispersion amplitude compensation wave equation to obtain the phase-delay amplitude-decayed geophone wavefield and the phase-delay amplitude-constant geophone wavefield respectively; S3. Calculate the compensation coefficient and the geophone compensation wavefield respectively through the phase-delay amplitude-decayed geophone wavefield and the phase-delay amplitude-constant geophone wavefield; S4. Calculate the seismic source wavefield separated into P-wave and S-wave and the geophone compensation wavefield separated into P-wave and S-wave based on the Helmholtz decomposition method to obtain the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield; S5. Under stable imaging conditions, use the P-wave component of the seismic source wavefield, the S-wave component of the seismic source wavefield, the P-wave component of the geophone compensation wavefield, and the S-wave component of the geophone compensation wavefield at N moments to calculate the single-shot migration results of PP component, PS component, SS component, and SP component respectively; S6. Repeat S2 - S5, stack all single-shot compensation migration results to obtain the underground compensation profile, calculate the multi-shot migration stack result, and obtain the multi-component compensation migration record.

2. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, wherein Obtaining the phase-delay amplitude-decayed seismic source wavefield of the corresponding shot point includes: Using the corresponding shot point position information and the known conditions, calculate the seismic source wavefield, and bring the seismic source wavelet and the single-shot observation system into the phase-delay amplitude compensation wave equation for wavefield continuation to obtain the phase-delay amplitude-decayed seismic source wavefield of the corresponding shot point; where the method for performing the wavefield continuation is: In equations (1) and (2), A xx , A zz , A xz are the stress components of the wave field respectively, B x , B z are the horizontal and vertical velocity components of the wave field respectively, S(x, z, t) is the loaded source function, ρ is the known density model, V p , V s , Q p , Q s correspond to the known P-wave velocity model, S-wave velocity model, P-wave quality factor model, and S-wave quality factor model respectively, t is time, x is distance, z is depth, is the Laplace operator, e p , f p , d p , d s , e s , f s are intermediate variables respectively, and ω0 is the dominant frequency of the wavelet.

3. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 2, characterized in that In the wavefield continuation formula, the integer-order partial derivative with respect to time is solved using second-order finite difference, and the fractional-order partial derivative with respect to time and the spatial partial derivative term are solved using Fourier transform. The specific discrete iteration formula is: In Equation (3), Δt is the calculation time interval, A xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1 are the wave field values of each component at the time of (n + 1)Δt; A xx n , A zz n , A xz n , B x n , B z n are the wave field values of each component at the time of nΔt; B x n-1 , B z n-1 are the wave field values of each component at the time of (n - 1)Δt; e p , f p , d p , d s , e s , f s are intermediate variables, F is the forward Fourier transform, F -1 is the inverse Fourier transform; k x , k z are the wave numbers in each direction, is the average of f p for the entire model, i is the imaginary unit, S n is the source amplitude value at the time of nΔt; Using the wave field values of each component at the times of T = nΔt and T = (n - 1)Δt (A xx n , A zz n , A xz n , B x n , B z n , B x n-1 , B z n-1 ) to iteratively calculate the wave field values at the time of T = (n + 1)Δt (A xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1 )； Substitute the source function \(S(x,z,t)\), calculate a total of \(N\) time instants, start loading the source function from \(T = 0\), and finally obtain the source wavefield values at \(N\) time instants, including \(B\) x and \(B\) z and \(A\) xx and \(A\) zz and \(A\) xz components, which together form the phase-delay amplitude-decay source wavefield \(W\) shot (x,z).

4. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 3, characterized in that The method for obtaining the phase-delay amplitude-decayed geophone wavefield is: Based on the vertical velocity component and horizontal velocity component of the seismic record, perform reverse-time iterative calculation, specifically: In Equation (4), R x (x, z, t) and R z (x, z, t) are the x and z components of the seismic record, respectively; The specific calculation iteration method is: In Equation (5), R x n (x, z) and R z n (x, z) are the x-component and z-component of the seismic record at nΔt time respectively, Δt is the calculation time interval. Using the wave field values of each component at T = nΔt and T = (n + 1)Δt (A xx n , A zz n , A xz n , B x n , B z n , B x n+1 , B z n+1 ), and the seismic record values R x n (x, z) and R z n (x, z) at T = nΔt, the wave field values at T = (n - 1)Δt (A xx n-1 , A zz n -1 , A xz n-1 , B x n-1 , B z n-1 ) are obtained by iterative calculation; Calculate the wave field at the geophone points at N moments and substitute it into the vertical and horizontal components R of the seismic record x (x,z,t) and R z (x,z,t). Start loading the seismic record from T = N for iteration, and finally obtain the wave field values at the geophone points at N moments, including B x and B z and A xx and A zz and A xz components, which constitute the wave field W of the phase delay amplitude attenuation at the geophone points corresponding to the shot point rec1 (x,z).

5. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 4, characterized in that, Method for obtaining phase-delay amplitude-invariant detected wavefield W rec2 (x, z) is as follows: The specific calculation iteration method is: In equations (6) and (7), using the component wave field values (A xx n 、A zz n 、A xz n 、B x n 、B z n 、B x n+1 、B z n+1 ) at times T = nΔt and T = (n + 1)Δt, and the seismic record values R x n (x,z)、R z n (x,z) at time T = nΔt, the wave field values (A xx n-1 、A zz n-1 、A xz n-1 、B x n-1 、B z n-1 ) at time T = (n - 1)Δt are obtained through iterative calculation; Calculate the wave field at the geophone point at N moments and substitute it into the vertical and horizontal components R x (x,z,t) and R z (x,z,t). Start loading the seismic record from T = N for iteration, and finally obtain the wave field values at the geophone point at N moments, including B x and B z and A xx and A zz and A xz components, which together form the phase-delay amplitude-invariant geophone wave field W rec2 (x,z).

6. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, wherein The method for calculating the compensation coefficient is: where ε is white noise, is the compensation coefficient, and W rec1 (x, z) is the wave field of the phase-delay amplitude attenuation geophone point, and W rec2 (x, z) is the wave field of the phase-delay amplitude invariant geophone point; The method for calculating the geophone compensation wavefield is: In the formula, is the compensated wave field at the detection point.

7. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, wherein The method for calculating the single-shot migration results of PP component, PS component, SS component, and SP component is: where ε is white noise, is the compensated wavefield at the detection point, and W shot (x, z) is the source wavefield with phase delay and amplitude attenuation, I(x, z) is the single-shot migration result component, and N is the number of time intervals.

8. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 7, characterized in that The method for calculating the multi-shot migration stack result is: Where, I stack (x, z) is the stacking result, and M is the total number of calculated shots.

9. A viscoelastic medium compensation migration system for multi-component seismic records, which is applied to the viscoelastic medium compensation migration method for multi-component seismic records according to any one of claims 1-8, and is characterized in that, Including: Phase-delay amplitude-decayed seismic source wavefield acquisition module: used for single-shot data, using the corresponding shot point position information and known conditions to obtain the phase-delay amplitude-decayed seismic source wavefield of the corresponding shot point, where the known conditions include multi-component reflection wave seismic records of multiple shots, the P-wave velocity model, S-wave velocity model, P-wave quality factor, and S-wave quality factor obtained after velocity modeling of the target area; Reverse wavefield continuation processing module: It is used to perform reverse wavefield continuation according to the phase dispersion amplitude compensation wave equation by using the preprocessed multi-component reflected wave seismic records, and respectively obtain the phase-delay amplitude-attenuated geophone wavefield and the phase-delay amplitude-constant geophone wavefield; Compensation coefficient and geophone compensated wavefield calculation module: It is used to calculate the compensation coefficient and the geophone compensated wavefield respectively through the phase-delay amplitude-attenuated geophone wavefield and the phase-delay amplitude-constant geophone wavefield; P-S wave separation module: It is used to calculate the source wavefield of P-S wave separation and the geophone compensated wavefield of P-S wave separation based on the Helmholtz decomposition method, and obtain the longitudinal wave component of the source wavefield, the transverse wave component of the source wavefield, the longitudinal wave component of the geophone compensated wavefield, and the transverse wave component of the geophone compensated wavefield; Single-shot migration result component acquisition module: It is used to calculate the single-shot migration result PP component, PS component, SS component and SP component respectively by using the longitudinal wave component of the source wavefield, the transverse wave component of the source wavefield, the longitudinal wave component of the geophone compensated wavefield, and the transverse wave component of the geophone compensated wavefield at N moments under the stable imaging condition; Multi-component compensated migration output module: It is used to superimpose all single-shot compensated migration results to obtain the underground compensated profile, calculate the multi-shot migration stack result, and obtain the multi-component compensated migration record.

10. An electronic device, characterized in that, It includes a processor, a memory, and a program or instruction stored on the memory and executable on the processor. When the program or instruction is executed by the processor, it implements the steps of the viscoelastic medium compensated migration method for multi-component seismic records according to any one of claims 1-8.

Images