Viscoelastic medium offset compensation method and system for multi-component seismic record and electronic equipment

Abstract

The invention discloses a viscoelastic medium offset compensation method and system for a multi-component seismic record and electronic equipment, and the method comprises the steps: obtaining a phase delay amplitude attenuation seismic source wave field of a corresponding shot point through the position information of the corresponding shot point and a known condition; performing reverse wave field continuation by using the preprocessed multi-component reflected wave seismic record according to a phase frequency dispersion amplitude compensation wave equation to obtain a phase delay amplitude attenuation detection point wave field and a phase delay amplitude invariant detection point wave field, and respectively calculating a compensation coefficient and a detection point compensation wave field; calculating a seismic source wave field and a detection point compensation wave field of longitudinal and transverse wave separation; respectively calculating a PP component, a PS component, an SS component and an SP component of a single-shot offset result; and superposing all single-shot offset compensation results to obtain an underground compensation profile, and calculating a multi-shot offset superposition result to obtain a multi-component offset compensation record. Compared with other viscous medium offset technologies, the method has the advantage that the calculation process is more stable.

Description

Technical Field

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

[0002] The actual underground medium is closer to non-perfect elasticity. The multi-component seismic records collected by seismic exploration are subject to absorption and attenuation, and the amplitude is attenuated and the phase is distorted. If only the elastic wave migration method is used, the imaging results will have the problem of amplitude and phase distortion. Considering the offset of the medium viscoelasticity can compensate for the amplitude of the seismic wave, correct the phase, and obtain a higher resolution result.

[0003] In the prior art, Zhang et al. (2010) added a pseudo-difference operator to the constant-Q model viscoacoustic wave equation to perform viscoacoustic reverse time migration imaging to eliminate amplitude attenuation and velocity dispersion. Since the seismic energy increases with the frequency during the reverse extension of the wave field, a regularization step was introduced to avoid numerical instability, thereby realizing the reverse propagation of the viscous wave field. Suh et al. (2012) developed this operator into VTI media. Bai et al. (2013) used a similar method to compensate for attenuation during reverse time migration imaging, using a viscoacoustic wave equation that does not require memory variables. This type of algorithm can compensate for energy amplitude to a certain extent, but further research is needed on the issues of multi-component and phase correction. Zhu et al. (2014) proposed a new decoupled fractional-order Laplace constant-Q viscoacoustic wave equation, which divides the attenuation term into two terms: amplitude attenuation and velocity dispersion. It was applied to Q-compensated viscoacoustic medium reverse time migration imaging and was measured on theoretical and actual data. In order to ensure the stability of the wave equation extension in the Q compensation migration method, Fletcher et al. (2012) proposed to design independent filter operators for both amplitude and phase. By performing two acoustic wave field extensions, the filter operators were applied to the source and receiver wave fields before Q compensation imaging. The filtering method is based on the travel time attenuation of the propagation path, so it may not be applicable to the situation where multiple waves and Q models change drastically. Zhao et al. (2018) used the excitation amplitude imaging condition to perform viscoelastic medium migration imaging. Some scholars have realized the migration imaging problem of viscous media by transforming 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 without stability problems (Guo et al., 2018), but it inevitably brings a larger amount of calculation. At present, the migration for viscous media is mainly for single-component data, only considering the propagation process of longitudinal waves, or the viscoelastic model used cannot cope with complex absorption and attenuation areas. The high-precision viscoelastic medium compensation migration method for multi-component data needs to be improved.

[0004] In the existing technology, only elastic media are considered for multi-component seismic recordings, but they cannot cope with areas with strong absorption attenuation. Existing viscous medium compensation imaging technology is mainly based on viscoacoustic media for compensation offset, which cannot fully utilize the properties of elastic waves. It only uses longitudinal wave data components for imaging and cannot fully utilize multi-component data. Although some public materials consider the viscoelasticity of the medium and construct wave field imaging based on cross-correlation imaging conditions, there are instabilities in the application process and it cannot cope with areas with more complex lateral changes. Summary of the invention

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

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

[0007] S1. Acquire multi-component reflection wave seismic records of multiple shots, and for single shot data, obtain the phase-delayed amplitude-attenuated source wave field of the corresponding shot point using the corresponding shot point position information and known conditions, wherein 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 is completed;

[0008] S2, using the pre-processed multi-component reflection wave seismic record to perform reverse wave field extension according to the phase dispersion amplitude compensation wave equation, to obtain a phase-delayed amplitude attenuated detection point wave field and a phase-delayed amplitude unchanged detection point wave field;

[0009] S3, respectively calculating a compensation coefficient and a detection point compensation wavefield through the phase-delayed amplitude-attenuated detection point wavefield and the phase-delayed amplitude-unchanged detection point wavefield;

[0010] S4. Calculate the source wave field separated by P- and S-waves and the compensation wave field of the detection point separated by P- and S-waves based on the Helmholtz decomposition method to obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point;

[0011] S5, under stable imaging conditions, using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component, detection point compensation wave field shear wave component at N moments, respectively calculating the single shot migration result PP component, PS component, SS component and SP component;

[0012] S6. Repeat S2-S5, superimpose all single-shot compensation offset results to obtain an underground compensation profile, calculate multi-shot offset superposition results, and obtain a multi-component compensation offset record.

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

[0014] The source wave field is calculated by using the corresponding shot point position information and the known conditions, and the source wavelet and the single shot observation system are introduced into the phase-delayed amplitude compensation wave equation to perform wave field extension, so as to obtain the phase-delayed amplitude attenuated source wave field of the corresponding shot point; wherein the method for performing the wave field extension is:

[0015]

[0016] In formula (1) and formula (2), A xx , A zz , A xz are the stress components of the wave field, B x , B z are the horizontal velocity component and the vertical velocity component of the wave field, S(x,z,t) is the loaded source function, ρ is the known density model, V p 、V s , Q p , Q s They correspond to the known P-wave velocity model, S-wave velocity model, P-wave quality factor model, S-wave quality factor model, 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, ω 0 is the main frequency of the wavelet.

[0017] Preferably, in the wave field extension formula, the time integer-order partial derivative is solved by second-order finite difference, and the time fractional-order partial derivative and the space partial derivative are solved by Fourier transform. The specific discrete iterative formula is:

[0018]

[0019] In formula (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+1are the wave field values ​​of each component at (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 time nΔt; B x n-1 , B z n-1 are the wave field values ​​of each component at (n-1)Δt; e p 、f p d p d s 、e s 、f s is the intermediate variable, F is the Fourier transform, F -1 is the inverse Fourier transform; k x , k z are the wave numbers in each direction, f p The average of the entire model, i is the imaginary basic unit, S n is the source amplitude value at time nΔt; using the wave field values ​​of each component at time T = nΔt, 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 ) Iterative calculation to obtain the wave field value (A) at time T = (n + 1) Δt xx n+1 , A zz n+1 , A xz n+1 , B x n+1 , B z n+1 );

[0020] Substitute the source function S(x,z,t) and calculate N time points in total. Load the source function starting from T=0 and finally obtain the source wave field values ​​at N time points, including B x , B z , A xx , Azz , A xz The components constitute the phase-delayed amplitude-attenuated source wave field W corresponding to the shot point. shot (x,z).

[0021] Preferably, the method for obtaining the phase-delayed amplitude-attenuated detection point wave field is:

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

[0023]

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

[0025] The specific calculation iteration method is:

[0026]

[0027] In formula (5), R x n (x,z),R z n (x,z) are the x component and z component of the earthquake record at time nΔt, Δt is the calculation time interval, and the wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt xx n-1 , A zz n-1 , A xz n-1 , B x n-1 , B z n-1 );

[0028] Calculate the wave field of the N detection points at each moment and substitute the vertical and horizontal components R of the seismic record into x (x,z,t),R z (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz Component, that is, the phase delay amplitude attenuation detection point wave field W of the corresponding shot point rec1 (x,z).

[0029] Preferably, the phase-delayed amplitude-invariant detection point wave field W is obtained rec2 The method for (x,z) is:

[0030]

[0031] The specific calculation iteration method is:

[0032]

[0033] In formula (6) and formula (7), the wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt xx n-1 , A zz n-1 , A xz n-1 , B x n-1 , B z n-1 );

[0034] Calculate the wave field of the N detection points at each moment and substitute the vertical and horizontal components R of the seismic record into x (x,z,t),Rz (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz The components constitute the phase-delayed amplitude-invariant detection point wave field W corresponding to the shot point rec2 (x,z).

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

[0036]

[0037] In the formula, ε is white noise, is the compensation coefficient, W rec1 (x,z) is the phase-delayed amplitude-attenuated detection point wave field, W rec2 (x,z) is the phase-delayed amplitude-invariant detection point wave field;

[0038] The method for calculating the compensation wave field of the detection point is:

[0039]

[0040] In the formula, Compensate the wave field for the detection point.

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

[0042]

[0043] In the formula, ε is white noise, is the compensation wave field of the detection point, W shot (x,z) is the phase-delayed amplitude-attenuated 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 offset stacking result is:

[0045]

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

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

[0048] Phase-delayed amplitude-attenuated source wavefield acquisition module: used for obtaining the phase-delayed amplitude-attenuated source wavefield of the corresponding shot point for single-shot data using the corresponding shot point position information and known conditions, wherein 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 extension processing module: used to use the pre-processed multi-component reflection wave seismic records to perform reverse wavefield extension according to the phase dispersion amplitude compensation wave equation, and obtain the phase-delayed amplitude attenuation detection point wavefield and the phase-delayed amplitude unchanged detection point wavefield respectively;

[0050] A compensation coefficient and detection point compensation wave field calculation module is used to calculate the compensation coefficient and the detection point compensation wave field respectively through the phase delay amplitude attenuation detection point wave field and the phase delay amplitude unchanged detection point wave field;

[0051] P-wave and S-wave separation module: used to calculate the source wave field of P-wave and S-wave separation and the compensation wave field of the detection point of P-wave and S-wave separation based on the Helmholtz decomposition method, and obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point;

[0052] Single shot migration result component acquisition module: used to calculate the PP component, PS component, SS component and SP component of the single shot migration result respectively by using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component and detection point compensation wave field shear wave component at N moments under stable imaging conditions;

[0053] Multi-component compensation offset output module: used to superimpose all single-shot compensation offset results to obtain underground compensation profiles, calculate multi-shot offset superposition results, and obtain multi-component compensation offset records.

[0054] This embodiment also provides an electronic device, including a processor, a memory, and a program or instruction stored in the memory and executable on the processor, wherein the program or instruction, when executed by the processor, implements the steps of a viscoelastic medium compensation offset 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 takes into account the medium absorption attenuation characteristics and elastic characteristics at the same time. After compensating the energy, the imaging position in the gas cloud area is accurate and the resolution is higher. Compared with other viscous medium migration technologies, the calculation process of the present invention is more stable. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

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

[0059] Figure 2 Schematic diagram of longitudinal and transverse wave velocities and Q value parameters of the model of the embodiment of the present invention, wherein (a) is a longitudinal wave velocity model, (b) is a transverse wave velocity model, (c) is a longitudinal wave quality factor model, and (d) is a transverse wave quality factor model;

[0060] Figure 3 Schematic diagram of a typical single-shot seismic record of an embodiment of the present invention, wherein (a) is the horizontal velocity component, and (b) is the vertical velocity component;

[0061] Figure 4 Schematic diagram of the source wave field of a single shot at 500ms in an embodiment of the present invention, where (a) is B x Component, (b) is B z Quantity;

[0062] Figure 5 Schematic diagram of the wave field of the phase-delayed amplitude-attenuated detection point of a single shot at 500ms in an embodiment of the present invention, where (a) is B x Component, (b) is B z Quantity;

[0063] Figure 6 Schematic diagram of the wave field of the phase-delayed amplitude-invariant detection point of a single shot at 500ms in an embodiment of the present invention, where (a) is B x Component, (b) is B z Quantity;

[0064] Figure 7 Schematic diagram of the longitudinal and transverse wave components of the source wave field and the longitudinal and transverse wave components of the compensation detection point wave field at 500ms for a single shot according to an embodiment of the present invention, where (a) is W shot (x,z) P component, (b) is W shot (x,z) S Component, (c) is Component, (d) is Quantity;

[0065] Figure 8 Schematic diagram of single shot migration results of an embodiment of the present invention, wherein (a) is the PP component, (b) is the PS component, (c) is the SP component, and (d) is the SS component;

[0066] Fig. 9 It is a schematic diagram of a multi-component superposition cross-section of an embodiment of the present invention, wherein (a) is the PP component, (b) is the PS component, (c) is the SP component, and (d) is the SS component. DETAILED DESCRIPTION

[0067] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

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

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

[0070] S1. For single shot data, using the corresponding shot point location information and known conditions, obtain the phase-delayed amplitude-attenuated source wave field of the corresponding shot point, wherein 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;

[0071] S2, using the pre-processed multi-component reflection wave seismic record to perform reverse wave field extension according to the phase dispersion amplitude compensation wave equation, to obtain a phase-delayed amplitude attenuated detection point wave field and a phase-delayed amplitude unchanged detection point wave field;

[0072] S3, respectively calculating a compensation coefficient and a detection point compensation wavefield through the phase-delayed amplitude-attenuated detection point wavefield and the phase-delayed amplitude-unchanged detection point wavefield;

[0073] S4. Calculate the source wave field separated by P- and S-waves and the compensation wave field of the detection point separated by P- and S-waves based on the Helmholtz decomposition method to obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point;

[0074] S5, under stable imaging conditions, using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component, detection point compensation wave field shear wave component at N moments, respectively calculating the single shot migration result PP component, PS component, SS component and SP component;

[0075] S6. Repeat S2-S5, superimpose all single-shot compensation offset results to obtain an underground compensation profile, calculate multi-shot offset superposition results, and obtain a multi-component compensation offset record.

[0076] Compared with the existing multi-component migration imaging technology, this embodiment takes into account both the absorption and attenuation characteristics of the medium and the elastic characteristics. After compensating the energy, the imaging position in 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] Furthermore, the phase-delayed amplitude-attenuated source wave field corresponding to the shot point is obtained, including:

[0078] The source wavefield is calculated using the corresponding shot point location information and the previously known conditions, and the source wavelet and single-shot observation system are introduced into the phase-delayed amplitude-compensated wave equation for wavefield extension to obtain the phase-delayed amplitude-attenuated source wavefield of the corresponding shot point.

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

[0080] For the calculated single shot data, the source wave field is calculated using the corresponding shot point location information and the previously known conditions. The source wavelet and the single shot observation system are introduced into the phase-delayed amplitude compensation wave equation for wave field extension to obtain the phase-delayed amplitude attenuated source wave field corresponding to the shot point. The equations for wave field extension are shown in formulas (1)-(2).

[0081]

[0082]

[0083] In the formula, A xx , A zz , A xz are the stress components of the wave field, B x , B z They are the horizontal velocity component and the vertical velocity component of the wave field (the three stress components and the two velocity components constitute the phase-delayed amplitude-attenuated source wave field W shot (x,z)), S(x,z,t) is the source function of the load, ρ is the known density model, V p 、V s , Q p , Q s They correspond to the known P-wave velocity model, S-wave velocity model, P-wave quality factor model, S-wave quality factor model, 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, ω 0 is the main frequency of the wavelet.

[0084] Furthermore, in the wave field extension formula, the time integer-order partial derivatives are solved using the second-order finite difference method, and the time fractional-order partial derivatives and spatial partial derivatives are solved using Fourier transform. The specific discrete iterative formula is:

[0085]

[0086] In formula (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 (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 time nΔt; B x n-1 , B z n-1 are the wave field values ​​of each component at (n-1)Δt; e p 、f p d p d s 、e s 、f s is the intermediate variable, F is the Fourier transform, F -1 is the inverse Fourier transform; k x , k z are the wave numbers in each direction, f p The average of the entire model, i is the imaginary basic unit, S n is the source amplitude value at time nΔt, and the wave field values ​​of each component at time 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 ) Iterative calculation to obtain the wave field value (A) at time T = (n + 1) Δt 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) and calculate N time points in total. Load the source function starting from T=0 and finally obtain the source wave field values ​​at N time points, including B x , B z , A xx , A zz , A xz The components constitute the phase-delayed amplitude-attenuated source wave field W corresponding to the shot point. shot (x,z).

[0088] Furthermore, the method for obtaining the phase-delayed amplitude-attenuated detection point wave field is:

[0089] The phase-delayed amplitude-attenuated detection point wave field W is obtained by using the multi-component seismic records after conventional noise suppression and the reverse wave field continuation according to the phase-dispersion amplitude compensation wave equation. rec1 (x,z), the calculation is shown in formula (4), the vertical velocity component and horizontal velocity component of the seismic record are substituted into the equation and reverse time iteration is performed:

[0090]

[0091] The specific calculation iterative 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 earthquake record respectively; the other components have the same meanings as (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 earthquake record at time nΔt, respectively. The other components have the same meaning as in formula (3). Using the wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt 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 fields of the detection points at N moments are calculated and the vertical and horizontal components R of the seismic records are substituted into x (x,z,t),R z (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz Component, that is, the phase delay amplitude attenuation detection point wave field W of the corresponding shot point rec1 (x,z).

[0095] Furthermore, the phase-delayed amplitude-invariant detection point wave field W is obtained. rec2 The method for (x,z) is:

[0096] The phase-delayed amplitude-invariant detection point wave field is obtained by reverse wavefield extension based on the phase-dispersion amplitude compensation wave equation using the preprocessed multi-component seismic records:

[0097]

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

[0099]

[0100] The components in equations (6) and (7) have the same meaning as those in equations (1), (2), and (3). The wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt xx n-1 , A zz n-1 , A xz n-1 , B x n-1 , B z n-1 );

[0101] Calculate the wave field of the N detection points at each moment and substitute the vertical and horizontal components R of the seismic record into x (x,z,t),R z (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz The components constitute the phase-delayed amplitude-invariant detection point wave field W corresponding to the shot point rec2 (x,z).

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

[0103]

[0104] In the formula, ε is white noise, is the compensation coefficient, W rec1 (x,z) is the phase-delayed amplitude-attenuated detection point wave field, W rec2(x,z) is the phase-delayed amplitude-invariant detection point wave field;

[0105] The method for calculating the compensation wave field of the detection point is:

[0106]

[0107] In the formula, Compensate the wave field for the detection point.

[0108] Furthermore, the method for calculating the source wavefield of P-wave and S-wave separation and the compensation wavefield of the detection point of P-wave and S-wave separation based on the Helmholtz decomposition method is as follows:

[0109]

[0110] The longitudinal wave component and the transverse wave component (W shot (x,z) P , W shot (x,z) S ), compensate the longitudinal and transverse components of the wave field at the detection point

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

[0112] The imaging conditions are shown in formula (11). Substituting the longitudinal and transverse wave components of the source wave field and the compensation detection point wave field at N moments into the formula, the PP component, PS component, SS component and SP component of the single shot migration result (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 compensation wave field of the detection point, W shot (x,z) is the phase-delayed amplitude-attenuated source wavefield, I(x,z) is the component of the single-shot migration result, and N is the total number of time intervals.

[0115] Furthermore, all the single-shot compensation migration results are superimposed to obtain the underground compensation profile:

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

[0117]

[0118] In the formula, I stack (x, z) is the multi-component stacked section, and M is the total number of shots.

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

[0120] Phase-delayed amplitude-attenuated source wavefield acquisition module: used for obtaining the phase-delayed amplitude-attenuated source wavefield of the corresponding shot point for single-shot data using the corresponding shot point position information and known conditions, wherein 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;

[0121] Reverse wavefield extension processing module: used to use the pre-processed multi-component reflection wave seismic records to perform reverse wavefield extension according to the phase dispersion amplitude compensation wave equation, and obtain the phase-delayed amplitude attenuation detection point wavefield and the phase-delayed amplitude unchanged detection point wavefield respectively;

[0122] A compensation coefficient and detection point compensation wave field calculation module is used to calculate the compensation coefficient and the detection point compensation wave field respectively through the phase delay amplitude attenuation detection point wave field and the phase delay amplitude unchanged detection point wave field;

[0123] P-wave and S-wave separation module: used to calculate the source wave field of P-wave and S-wave separation and the compensation wave field of the detection point of P-wave and S-wave separation based on the Helmholtz decomposition method, and obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point;

[0124] Single shot migration result component acquisition module: used to calculate the PP component, PS component, SS component and SP component of the single shot migration result respectively by using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component and detection point compensation wave field shear wave component at N moments under stable imaging conditions;

[0125] Multi-component compensation offset output module: used to superimpose all single-shot compensation offset results to obtain underground compensation profiles, calculate multi-shot offset superposition results, and obtain multi-component compensation offset records.

[0126] This embodiment also provides an electronic device, including a processor, a memory, and a program or instruction stored in the memory and executable on the processor, wherein the program or instruction, when executed by the processor, implements the steps of a viscoelastic medium compensation offset method for multi-component seismic records.

[0127] In order to more clearly express the technical solution of the present invention, the following specific embodiments are provided to introduce the solution:

[0128] Multi-component seismic records of a classic gas cloud region model are selected for migration imaging.

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

[0130] The basic parameters of the model are known. Figure 2 (a)- Figure 2 (d) shows the P-wave velocity, S-wave velocity, P-wave quality factor, and S-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. Multi-component seismic recording parameters:

[0132] The multi-component seismic record based on this model is known to have a total of 40 shots, with a shot depth of 0, a shot start position of 0, and a shot spacing of 125m; each shot has 398 receivers, with a receiving point depth of 0, a track spacing of 12.5m, a detector coverage range from 0 to 4962.5m, a time sampling interval of 1ms, and a record length of 3.5s. The source is a 25Hz Ricker wavelet. A typical single-shot multi-component seismic record (shot 21) is shown below: Figure 3 (a)- Figure 3 (b) as shown.

[0133] S2.3. Calculate the phase-delayed amplitude-attenuated source wave field:

[0134] Using the known model-related parameters in S2.1, the Ricker wavelet source and the recursive formula (3), the source wavefield at N times is calculated, among which the 500ms time B x , B z Quantity Figure 4 (a)- Figure 4 (b) as shown.

[0135] S2.4. Calculate the phase delay amplitude attenuation detection point wave field:

[0136] Using the single-shot seismic record in S2.2, the known model-related parameters in S2.1 and the recursive formula (5), the phase-delayed amplitude-attenuated detection point wave fields at N times are calculated, among which the 500ms time B x , B z Quantity Figure 5 (a)- Figure 5 (b) as shown.

[0137] S2.5. Calculate the phase-delayed amplitude-invariant detection point wave field:

[0138] Using the single-shot seismic record in S2.2, the known model-related parameters in S2.1 and the recursive formula (7), the phase-delayed amplitude-invariant detection point wave fields at N times are calculated, among which the 500ms time B x , B z Quantity Figure 6 (a)- Figure 6 (b) as shown.

[0139] S2.6. Calculate the detection point wave field amplitude compensation coefficient and compensate the detection point wave field:

[0140] (1) Using S2.4 and S2.5, we can obtain the wave fields W of two detection points at N moments. rec1 (x,z),W rec2 The compensation coefficients of the components (x, z) are calculated by formula (8).

[0141] (2) Using the phase-delayed amplitude-invariant detection point wave field W rec2 Substitute (x, z) and the compensation coefficient into formula (9) to obtain the wave fields of the compensated detection points at N moments.

[0142] S2.7. Compensation wave fields of seismic sources and receivers based on Helmholtz decomposition for separation of P- and S-waves:

[0143] Based on formula (10), the components of the source wave field and the compensation detection point wave field at N moments are calculated to obtain the longitudinal wave component and the transverse wave component of the source wave field, and the longitudinal wave component and the transverse wave component of the compensation detection point wave field (W shot (x,z) P , W shot (x,z) S , Among them, the components at 500ms are as follows Figure 7 (a)- Figure 7 (d) as shown.

[0144] S2.8. Calculate single shot imaging results based on stable imaging conditions:

[0145] Substitute the longitudinal and transverse wave components of the source wave field and the compensation detection point wave field at N moments into formula (11) to obtain the PP, PS, SS, and SP components of the single-shot migration results, as follows: Figure 8 (a)- Figure 8 (d) as shown.

[0146] S2.9. Calculate the multi-shot offset superposition results and obtain the multi-component compensation offset record:

[0147] This model has a total of 40 shot seismic records. The calculation steps from S2.2 to S2.8 are repeated to obtain the migration imaging results of 40 shots. Finally, the stacking process is performed using formula (12) to obtain the multi-component results of the entire stacked section, as shown in Fig. 9 (a)- Fig. 9 (d) as shown.

[0148] The above are only preferred specific implementations of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A viscoelastic medium compensation migration method for multi-component seismic records, characterized in that: include: S1. Acquire multi-component reflection wave seismic records of multiple shots, and for single shot data, obtain the phase-delayed amplitude-attenuated source wave field of the corresponding shot point using the corresponding shot point position information and known conditions, wherein 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 is completed; S2, using the pre-processed multi-component reflection wave seismic record to perform reverse wave field extension according to the phase dispersion amplitude compensation wave equation, to obtain a phase-delayed amplitude attenuated detection point wave field and a phase-delayed amplitude unchanged detection point wave field; S3, respectively calculating a compensation coefficient and a detection point compensation wavefield through the phase-delayed amplitude-attenuated detection point wavefield and the phase-delayed amplitude-unchanged detection point wavefield; S4. Calculate the source wave field separated by P- and S-waves and the compensation wave field of the detection point separated by P- and S-waves based on the Helmholtz decomposition method to obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point; S5, under stable imaging conditions, using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component, detection point compensation wave field shear wave component at N moments, respectively calculating the single shot migration result PP component, PS component, SS component and SP component; S6. Repeat S2-S5, superimpose all single-shot compensation offset results to obtain an underground compensation profile, calculate multi-shot offset superposition results, and obtain a multi-component compensation offset record.

2. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, characterized in that: Obtain the phase-delayed amplitude-attenuated source wave field corresponding to the shot point, including: The source wave field is calculated by using the corresponding shot point position information and the known conditions, and the source wavelet and the single shot observation system are introduced into the phase-delayed amplitude compensation wave equation to perform wave field extension, so as to obtain the phase-delayed amplitude attenuated source wave field of the corresponding shot point; wherein the method for performing the wave field extension is: In formula (1) and formula (2), A xx , A zz , A xz are the stress components of the wave field, B x , B z are the horizontal velocity component and the vertical velocity component of the wave field, S(x,z,t) is the loaded source function, ρ is the known density model, V p 、V s , Q p , Q s They correspond to the known P-wave velocity model, S-wave velocity model, P-wave quality factor model, S-wave quality factor model, 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 main 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 wave field extension formula, the time integer-order partial derivative is solved by second-order finite difference, and the time fractional-order partial derivative and the spatial partial derivative are solved by Fourier transform. The specific discrete iterative formula is: In formula (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 (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 time nΔt; B x n-1 , B z n-1 are the wave field values ​​of each component at (n-1)Δt; e p 、f p ,d p ,d s 、e s 、f s is the intermediate variable, F is the Fourier transform, F -1 is the inverse Fourier transform; k x , k z are the wave numbers in each direction, f p The average of the entire model, i is the imaginary basic unit, S n is the source amplitude value at time nΔt; using the wave field values ​​of each component at time T = nΔt, 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 ) Iterative calculation to obtain the wave field value (A) at time T = (n + 1) Δt 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) and calculate N time points in total. Load the source function starting from T=0 and finally obtain the source wave field values ​​at N time points, including B x , B z , A xx , A zz , A xz The components constitute the phase-delayed amplitude-attenuated source wave field W corresponding to the shot point. 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-delayed amplitude-attenuated detection point wave field is: Based on the vertical velocity component and horizontal velocity component of the seismic record, reverse time iteration calculation is performed, specifically: In formula (4), R x (x,z,t),R z (x,z,t) are the x and z components of the earthquake record, respectively; The specific calculation iteration method is: In formula (5), R x n (x,z),R z n (x,z) are the x component and z component of the earthquake record at time nΔt, Δt is the calculation time interval, and the wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt xx n-1 , A zz n -1 , A xz n-1 , B x n-1 , B z n-1 ); Calculate the wave field of the N detection points at each moment and substitute the vertical and horizontal components R of the seismic record into x (x,z,t),R z (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz Component, that is, the phase delay amplitude attenuation detection point wave field W of the corresponding shot point rec1 (x,z).

5. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 4, characterized in that: Get the phase-delayed amplitude-invariant detection point wave field W rec2 The method for (x,z) is: The specific calculation iteration method is: In formula (6) and formula (7), the wave field values ​​of each component at time 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 ), T = earthquake record value R at time nΔt x n (x,z),R z n (x, z), iterative calculation to obtain the wave field value (A) at time T = (n-1) Δt xx n-1 , A zz n-1 , A xz n-1 , B x n-1 , B z n-1 ); Calculate the wave field of the N detection points at each moment and substitute the vertical and horizontal components R of the seismic record into x (x,z,t),R z (x,z,t), load the seismic records from T = N and iterate, and finally obtain the wave field values ​​of the detection points at N moments, including B x , B z , A xx , A zz , A xz The components constitute the phase-delayed amplitude-invariant detection point wave field W corresponding to the shot point rec2 (x,z).

6. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, characterized in that: The method for calculating the compensation coefficient is: In the formula, ε is white noise, is the compensation coefficient, W rec1 (x,z) is the phase-delayed amplitude-attenuated detection point wave field, W rec2 (x,z) is the phase-delayed amplitude-invariant detection point wave field; The method for calculating the compensation wave field of the detection point is: In the formula, Compensate the wave field for the detection point.

7. The viscoelastic medium compensation migration method for multi-component seismic records according to claim 1, characterized in that: The method for calculating the PP component, PS component, SS component and SP component of the single shot offset result is: In the formula, ε is white noise, is the compensation wave field of the detection point, W shot (x,z) is the phase-delayed amplitude-attenuated source wavefield, I(x,z) is the component of the single-shot migration result, 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 offset stacking results is: In the formula, I stack (x, z) is the superposition result, and M is the total number of calculated shots.

9. A viscoelastic medium compensation migration system for multi-component seismic records, applied to the viscoelastic medium compensation migration method for multi-component seismic records as claimed in any one of claims 1 to 8, characterized in that: include: Phase-delayed amplitude-attenuated source wavefield acquisition module: used for obtaining the phase-delayed amplitude-attenuated source wavefield of the corresponding shot point for single-shot data using the corresponding shot point position information and known conditions, wherein 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 extension processing module: used to use the pre-processed multi-component reflection wave seismic records to perform reverse wavefield extension according to the phase dispersion amplitude compensation wave equation, and obtain the phase-delayed amplitude attenuation detection point wavefield and the phase-delayed amplitude unchanged detection point wavefield respectively; A compensation coefficient and detection point compensation wave field calculation module is used to calculate the compensation coefficient and the detection point compensation wave field respectively through the phase delay amplitude attenuation detection point wave field and the phase delay amplitude unchanged detection point wave field; P-wave and S-wave separation module: used to calculate the source wave field of P-wave and S-wave separation and the compensation wave field of the detection point of P-wave and S-wave separation based on the Helmholtz decomposition method, and obtain the P-wave component of the source wave field, the S-wave component of the source wave field, the P-wave component of the compensation wave field of the detection point, and the S-wave component of the compensation wave field of the detection point; Single shot migration result component acquisition module: used to calculate the PP component, PS component, SS component and SP component of the single shot migration result respectively by using the source wave field longitudinal wave component, source wave field shear wave component, detection point compensation wave field longitudinal wave component and detection point compensation wave field shear wave component at N moments under stable imaging conditions; Multi-component compensation offset output module: used to superimpose all single-shot compensation offset results to obtain underground compensation profiles, calculate multi-shot offset superposition results, and obtain multi-component compensation offset records.

10. An electronic device, characterized in that: The method comprises a processor, a memory, and a program or instruction stored in the memory and executable on the processor, wherein the program or instruction, when executed by the processor, implements the steps of the viscoelastic medium compensation migration method for multi-component seismic records as described in any one of claims 1 to 8.

Images