An Elastic Reverse Time Migration Imaging Method Based on Transformation Sensing Gated

Abstract

This invention discloses an elastic reverse time migration imaging method based on conversion-sensor gating, relating to the field of seismic exploration technology. First, the elastic wavefield is separated into P-waves and S-waves to obtain wavefields of different modes. Then, based on the spatial gradient of the medium's Lamé parameters and density, an acceleration term driven solely by medium inhomogeneity is constructed, and curl-driven and divergence-driven coupled sources are obtained through Poisson projection. Combining the instantaneous power density, spatial smoothing results, and elastic parameter gradient information corresponding to the coupled sources, a material gating function and conversion intensity ratio are constructed, forming a conversion-sensor gating weight that varies with space and time. This weight is then introduced into the elastic reverse time migration imaging conditions to weight and control the imaging contribution. This method can improve the continuity, amplitude stability, and physical rationality of multi-mode imaging results under complex medium conditions, and has good practical application value.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration technology, and in particular to an elastic reverse time migration imaging method based on conversion sensing gating. Background Technology

[0002] With the development of multi-component seismic exploration and multi-wave imaging technologies, elastic reverse-time migration imaging has become an important means of high-precision imaging in complex structural areas. Multi-component seismic data simultaneously contains both P-wave and S-wave information, providing richer wavefield information for fault identification, reservoir prediction, and detailed interpretation of complex geological bodies. Therefore, it has significant application value in oil and gas exploration and geophysical imaging. Existing technologies have developed various P-wave and S-wave separation and imaging methods for multi-component seismic data processing and elastic reverse-time migration imaging, including Helmholtz decomposition, polarization direction projection, pseudo-differential operator separation, and vector wavefield decomposition, which have improved the quality of multi-wave imaging to a certain extent.

[0003] However, in practical applications, traditional elastic reverse-time migration imaging still faces significant crosstalk problems between P-waves and S-waves. This problem stems from two main causes: firstly, incomplete decoupling of P-waves and S-waves during imaging conditions easily leads to non-physical energy leakage and spurious reflections at steeply dipping interfaces, highly inhomogeneous media, and regions with abrupt velocity changes; secondly, the injection of seismic record boundaries and the reverse-time migration process introduce additional PS and SP conversion components, increasing the complexity of wavefield reconstruction. Furthermore, actual acquired data often includes prior P-wave and S-wave coupling due to factors such as imperfect source coupling and near-surface inhomogeneity, making it difficult to fundamentally eliminate crosstalk using conventional P-wave and S-wave separation methods alone. These problems further result in defects in the imaging results, such as artifact interference, poor phase axis continuity, amplitude distortion, and unclear converted wave responses, thus affecting the imaging accuracy and the reliability of geological interpretation under complex media conditions.

[0004] Therefore, there is an urgent need for an elastic reverse time migration imaging method applicable to complex media conditions, which can effectively constrain the P-wave and S-wave conversion relationship during wave field propagation, reverse time propagation and imaging, reduce non-physical crosstalk and artifact interference, preserve the true converted wave energy, and thus improve the continuity, stability and accuracy of multi-mode imaging results. Summary of the Invention

[0005] To address the aforementioned problems, this invention discloses an elastic reverse-time migration imaging method based on conversion-sensing gating, which effectively suppresses P-wave and S-wave conversion crosstalk and non-physical artifacts, improving the continuity, stability, and accuracy of multi-mode imaging results under complex media conditions. This method first separates the elastic wavefield into P-waves and S-waves, obtaining the wavefields of different wave modes after separation. Then, based on the spatial gradient of the medium's Lamé parameters and density, it constructs an acceleration term driven solely by medium inhomogeneity, and obtains the curl-driven acceleration term through Poisson projection. Coupled source and divergence driven The coupling source is used to characterize the location and type of wave mode conversion. Furthermore, by combining the instantaneous power density, spatial smoothing results, and elastic parameter gradient information corresponding to the coupling source, a material gating function and conversion intensity ratio are constructed to form a conversion-aware gating weight that varies spatially and temporally. This weight is then introduced into the elastic reverse-time migration imaging condition to weight and control the imaging contribution. This method can retain effective imaging energy only at physically reasonable conversion locations and times, automatically suppressing non-physical conversion responses in homogeneous medium regions. This effectively reduces PS crosstalk, co-mode reflection interference, multiple wave effects, and artifacts caused by existing wave mode coupling, improving the continuity, amplitude stability, and physical rationality of multi-mode imaging results under complex medium conditions, and has good practical application value.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] An elastic reverse time migration imaging method based on conversion sensing gating includes the following steps:

[0008] (1) Obtain multi-component earthquake records, source information, and background elastic parameter model, wherein the background elastic parameter model includes at least Lamé parameters. , and density ;

[0009] (2) Based on the elastic wave equation, the source wave field is propagated in the forward direction and the received wave field is propagated in the reverse time direction. The source wave field and the received wave field are separated into P-wave and S-wave respectively to obtain the separated source wave field and the received wave field.

[0010] (3) Based on the spatial gradient of the separated wave field and the background elastic parameter model, construct a non-uniform driving acceleration that only contains the gradient term of the medium parameter;

[0011] (4) Perform a Poisson projection on the non-uniform driving acceleration to obtain Coupled source and Coupled source;

[0012] (5) Construct the converted power density and the same-mode power density based on the coupled source and the separated wave field, and perform spatial smoothing on the power density;

[0013] (6) Construct a material gating function based on the parameter gradient of the background elastic parameter model, and construct the conversion intensity ratio by combining the smoothed conversion power density and the same-mode power density;

[0014] (7) Construct spatiotemporal transformation sensing gating weights based on the material gating function and the transformation intensity ratio;

[0015] (8) Introduce the spatiotemporal conversion sensing gating weight into the elastic reverse time-shift imaging condition, weight the imaging contribution and accumulate it over time to obtain one or more imaging results among PP, PS, SP and SS.

[0016] Optionally, the elastic wave equation used in step (2) satisfies:

[0017] (1);

[0018] in, Indicates spatial location Density at that location Represents a spatial position vector. Indicates time, Represents the particle vibration velocity vector wave field. Represents the stress tensor. Indicates the focal term. This represents the first derivative of the particle's vibrational velocity with respect to time. It represents the divergence of the stress tensor, i.e., the body force term generated by stress changes;

[0019] The corresponding wavefield separation relationship is:

[0020] (2);

[0021] in, This represents the P-wave particle vibration velocity vector wave field. Let S represent the particle vibration velocity vector wave field, with superscripts P and S indicating P-wave mode and S-wave mode respectively. The equation indicates that the particle vibration velocity vector wave field is decomposed into the sum of P-wave component and S-wave component.

[0022] Optionally, in step (3), the non-uniform driving acceleration retains only the term caused by spatial variations in the medium parameters, used to characterize the mode transition drive caused by medium non-uniformity. Non-uniform driving acceleration for:

[0023] (3);

[0024] Among them, wave mode index Take the P wave or S wave. ; For spatial components or spatial direction indices. ; , , Indicates by The first acceleration response generated by the velocity vector wave field of separated particles One portion, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first Partial derivatives in spatial directions, express In the wave field of the vibration velocity vector of separated particles, the first The component is related to the first... Partial derivatives in each spatial direction, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first The partial derivatives in the spatial direction; the Lamé parameter and the spatial gradient of density are used to characterize the driving effect of medium inhomogeneity on wave mode conversion.

[0025] Optionally, the Poisson projection in step (4) includes the following process:

[0026] right The coupling source first constructs the scalar curl of the P-wave non-uniform driving acceleration. Solve the stream function using the scalar curl as the source term. Used to recover the projected acceleration in the S-wave subspace:

[0027] (4);

[0028] in, The scalar curl representing the non-uniform driving acceleration of the P-wave is used to extract shear-type cross-mode information generated by the non-uniform driving acceleration of the incident P-wave and serves as the right-hand term of the subsequent Poisson projection to construct the coupling source of the P-wave to S-wave conversion. and These represent the non-uniform driving accelerations generated by the P-wave separated wavefield, respectively. direction and Components in direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the stream function obtained from the right-hand side of the P-wave curl, used for subsequent recovery of the S-wave subspace. Coupled source;

[0029] The stream function obtained from the Poisson equation is combined with its spatial derivatives to reconstruct the cross-mode projection components, thus obtaining... Coupled source :

[0030] (5);

[0031] in, This represents the S-wave subspace coupling source obtained by Poisson projection of P-wave non-uniform driving acceleration. The subscript S indicates projection onto the S-wave subspace, and the superscript P indicates that the coupling source is constructed by P-wave non-uniform driving acceleration. and Representing the stream function right direction and First-order partial derivative in the direction;

[0032] right Coupled source, first construct the divergence of S-wave non-uniform driving acceleration. Solve the scalar potential function using this divergence as the source term. Used to recover the projected acceleration in the P-wave subspace:

[0033] (6);

[0034] in, This represents the divergence of the non-uniform S-wave driven acceleration, which is used to extract the information of the compressed cross-wave mode generated by the incident non-uniform S-wave driven acceleration and serves as the right-hand term of the subsequent Poisson projection to construct the coupling source of the S-wave to P-wave conversion. and These represent the non-uniform driving accelerations generated by the S-wave separated wavefield in... direction and Components in direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the scalar potential function obtained from the right-hand side of the S-wave divergence, used for subsequent recovery of the P-wave subspace. Coupled source;

[0035] Based on the potential function obtained from the Poisson equation, its spatial derivatives are combined to reconstruct the cross-mode projection components, thus obtaining... Coupled source :

[0036] (7);

[0037] in, This represents a P-wave subspace coupled source obtained by Poisson projection of S-wave non-uniform driven acceleration. The subscript P indicates projection onto the P-wave subspace, and the superscript S indicates that the coupled source is constructed by S-wave non-uniform driven acceleration. and They represent scalar potential functions respectively. right direction and The first-order partial derivative in the direction.

[0038] Optionally, the conversion power density and the same-mode power density in step (5) are defined as follows:

[0039] (8a);

[0040] (8b);

[0041] (8c);

[0042] (8d);

[0043] in, This is used to extract local positive power and avoid mutual cancellation of oscillation phases; in equation (8a), express Instantaneous conversion power density of the conversion channel The S-wave velocity vector wave field is represented by the symbol ·, which indicates the vector dot product; in equation (8b), Represents the S-wave co-mode reference power density. This represents the same-mode projection acceleration obtained by projecting the S-wave non-uniform driving acceleration onto the S-wave subspace, used as... The S-wave co-mode reference term when the conversion intensity is normalized; in equation (8c), express Instantaneous conversion power density of the conversion channel This expression represents the P-wave velocity vector wave field, obtained by relating the target P-wave velocity field to... The dot product of coupled sources measures locality. Transformation response intensity; in equation (8d), This represents the P-wave co-mode reference power density. This represents the same-mode projected acceleration obtained by projecting the P-wave non-uniform driving acceleration onto the P-wave subspace; this quantity is used as... P-wave co-mode reference term when the conversion intensity is normalized.

[0044] Optionally, the spatial smoothing operator in step (5) is:

[0045] (9);

[0046] in, Represents grid points Power density after smoothing Represented by grid points Unsmoothed power density in the neighborhood centered on the target; For grid indexing; To smooth the window half-width, the grid spacing can be selected based on the dominant wavelength, grid spacing, and noise level. Generally, 1 to 5 grid points are used, with 2 or 3 grid points being preferred. The window size is... ; and To smooth the inner edge of the window direction and Integer offset of direction, all values ​​are in the range of - to Both are merely summation indices and do not represent P / S wave modes; and They are respectively direction and The unit grid vector in the direction;

[0047] After smoothing, the result is , , and .

[0048] Optionally, in step (6), the material gating function approaches zero in locally uniform regions and increases in regions with abrupt parameter changes, thus defining physically reasonable transition locations. The material gating function is constructed from the gradient of elastic parameters. First, the spatial gradient magnitudes of the Lamé parameters and density are calculated, then normalized using the corresponding parameter magnitudes, and the weighted summation of the normalized gradient contributions is truncated to the range of 0 to 1, satisfying:

[0049] (10);

[0050] (11);

[0051] in, , and They represent , and Spatial gradient magnitudes of the three elastic parameters; , They represent right direction and Partial derivatives in direction, , They represent right direction and Partial derivatives in direction, , They represent right direction and The larger the partial derivative of the direction and the larger the gradient magnitude, the stronger the change in medium parameters at that location, and the more likely it is to correspond to the physical wave mode conversion location; Represents the material gating function; , and They represent , and The weighting coefficients of the gradient term are used to adjust the relative contributions of the Lamé parameter and density gradient to the material gating function, and can be taken as values ​​between 0 and 1; in the absence of special priors, they can be taken as equal weights, for example, all three can be taken as 0.33; , and This indicates the magnitude or absolute value of the corresponding elastic parameter; As a stability constant, it is used to avoid the normalization denominator from approaching zero. It is a positive number that is much smaller than the representative value of the corresponding elastic parameter, for example, one-thousandth of the maximum value of the corresponding parameter model, so as to avoid the denominator from approaching zero without changing the gating response at the real interface. and Used to restrict the gating function to the range of 0 to 1.

[0052] Equation (10) is used to calculate the spatial gradient magnitudes of the Lamé parameters and density, respectively. The larger the gradient magnitude, the more likely the location is to be an interface or a strongly non-uniform region, and the more likely it is to support physically reasonable wave mode conversion. Equation (11) is used to weightedly synthesize the normalized gradient contributions of each elastic parameter into a material gating function. The material gating function approaches 0 in locally uniform regions, increases in parameter abrupt changes or strongly non-uniform regions, and is restricted to between 0 and 1.

[0053] Optionally, the conversion intensity ratio in step (6) satisfies:

[0054] (12a);

[0055] (12b);

[0056] (12c);

[0057] (12d);

[0058] In equation (12a), This represents the normalized intensity ratio of the P-wave co-mode channel, which measures the relative proportion of the P-wave co-mode response in the P-wave correlated response. The smoothed P-wave co-mode reference power density, For smoothed Conversion power density, To prevent a positive stability constant with a denominator of zero, its value should be smaller than the characteristic order of the denominator; it can be taken as 10 times the maximum value of the corresponding denominator. -6 Up to 10 -3 10 are preferred -6 To ensure numerical stability in the low-energy region and without affecting the intensity ratio in the effective conversion region; in equation (12b), express The normalized conversion intensity ratio of the conversion channel; in equation (12c), This represents the normalized intensity ratio of the S-wave co-mode channel, which measures the relative proportion of the S-wave co-mode response in the S-wave correlated response. For S-wave co-mode reference power density, for Conversion power density; in equation (12d), express The normalized conversion intensity ratio of the conversion channel; the larger this ratio, the more supportive the local wavefield. Transform it into something like a contribution.

[0059] The conversion intensity ratio in equations (12a) to (12d) is a normalized conversion index used to characterize the ratio of the cross-mode coupling intensity to the same-mode reference intensity; it does not represent a strict physical probability, but is used to measure the degree of support of the local wave field for the corresponding conversion channel in imaging conditions.

[0060] Optionally, the spatiotemporal transformation sensing gating weights in step (7) satisfy:

[0061] (13);

[0062] in, This indicates that the source-detector mode combination is... Spatiotemporal transformation sensing gating weights, , Both are P or S; and These represent the corresponding conversion intensity ratios calculated from the source-end wavefield and the receiver-end wavefield, respectively. Indicates the source end, The terminator represents the detector end; this formula requires that the gating weights be increased only when the model parameter gradient, the source end wavefield conversion index, and the detector end wavefield conversion index are all supported simultaneously.

[0063] Equation (13) multiplies the material gating function with the conversion intensity ratio of the same imaging channel at the source and detector ends to form a conversion-aware gating weight that changes with both space and time; the imaging contribution is preserved or enhanced only when the model parameter gradient and the wave field conversion index at both ends support the channel.

[0064] Optionally, step (8) includes elastic reverse time migration imaging conditions and imaging results. satisfy:

[0065] (14);

[0066] (15);

[0067] in, express The passage in spatial location and time Cross-correlation imaging terms at the location; Indicates the detector end Separate vector wave field Indicates the source end Separate vector wave field; express The final imaging result of the channel, , , and These represent the imaging results for the PP, PS, SP, and SS channels, respectively.

[0068] Equation (14) represents the introduction of conversion sensing gating weights into the cross-correlation imaging conditions of the vector wave fields at the source and detector ends, weighting the candidate imaging contributions at each time step. The imaging results of Equation (15) are uniformly denoted as... This means that the cross-correlation imaging contributions at each time step are summed after gating and weighting to obtain the final imaging result.

[0069] When the background medium is locally homogeneous, the gradient of the medium parameters satisfies and Therefore, Corresponding Coupled source, The coupling source, conversion power density, and conversion intensity ratio are all zero to suppress non-physical conversion contributions in local uniform regions.

[0070] The beneficial effects of the present invention are: 1) The present invention proposes an elastic reverse time migration imaging method based on conversion sensing gating. By performing reverse time back propagation on the P-wave and S-wave received wave fields, only physically reasonable converted wave contributions are retained during the imaging process. This can effectively suppress non-physical PS crosstalk, boundary injection artifacts and mode leakage, and improve the accuracy of multi-mode imaging results.

[0071] 2) This invention utilizes the gradient of medium parameters to construct a non-uniform acceleration term, and combines Poisson projection and spatiotemporal gating to constrain the conversion position. This can enhance the real converted wave response under complex medium conditions, reduce the influence of same-mode reflection, multiple waves and pre-stored crosstalk, and make the imaging results have better continuity, stability and amplitude preservation capabilities.

[0072] 3) The present invention has been verified by model experiments and actual multi-component seismic data, and has good applicability and noise resistance. It can obtain high-quality imaging results under different imaging source types and complex structural conditions, and has good engineering application value. Attached Figure Description

[0073] Figure 1 This is a schematic diagram of the elastic reverse time migration imaging method based on conversion sensing gating of the present invention;

[0074] Figure 2 This is a simple horizontal two-layer model shown in one embodiment of the present invention;

[0075] Figure 3 It uses the finite difference method to... Figure 2 The mixed-source seismic record obtained from the simple horizontal two-layer model shown;

[0076] Figure 4 It is by utilizing the present invention from Figure 2 The intermediate quantities of the incident P-wave at the source end obtained in the simple horizontal two-layer model shown are: where, Figure 4 (a) in the image is a snapshot of the wavefield of the separated P-wave components obtained using the finite difference method. Figure 4 (b) in the figure is a material gating method obtained using the present invention. Figure 4 (c)-(f) in the figure represent the Poisson projection components obtained using this invention. , , , , Figure 4 In this paper, (g)-(j) represent the instantaneous power density after spatial smoothing obtained using this invention. , , , ;

[0077] Figure 5 It is by utilizing the present invention from Figure 2 The longitudinal wave source migration imaging results obtained in the simple horizontal two-layer model shown;

[0078] Figure 6 This is an embodiment of the Overthrust model shown in this invention;

[0079] Figure 7 It is by utilizing the present invention from Figure 6 The relevant intermediate quantities obtained in the Overthrust model shown are: where, Figure 7 (a) A snapshot of the wavefield of the separated P-wave components obtained using the finite difference method. Figure 7 (b) A snapshot of the wave field of the separated S-wave components obtained using the finite difference method. Figure 7 (c)-(f) in the figure represent the instantaneous power density after spatial smoothing obtained using the present invention. , , , ;

[0080] Figure 8 It is by utilizing the present invention from Figure 6 The longitudinal wave source migration imaging results obtained from the Overthrust model are shown below;

[0081] Figure 9 This is a migration velocity model of actual multi-component seismic data shown in an embodiment of the present invention;

[0082] Figure 10 It is by utilizing the present invention from Figure 9 The results of longitudinal wave source migration imaging obtained from the actual data velocity model are shown. Detailed Implementation

[0083] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0084] An elastic reverse time migration imaging method based on conversion sensing gating, such as Figure 1 As shown, it includes the following steps:

[0085] (1) Obtain multi-component earthquake records, source information, and background elastic parameter model, wherein the background elastic parameter model includes at least Lamé parameters. , and density .

[0086] (2) Based on the elastic wave equation, the source wave field is propagated forward, and the received wave field is propagated backward. P-wave and S-wave separation is performed on the source and received wave fields respectively to obtain the separated source and received wave fields. The P-wave and S-wave separation is achieved using a separation method based on Helmholtz-type vector wave field decomposition. In a two-dimensional isotropic medium, the P-wave component is extracted using divergence correlation components, and the S-wave component is extracted using curl correlation components. The elastic wave equation satisfies:

[0087] (1);

[0088] in, Indicates spatial location Density at that location Represents a spatial position vector. Indicates time, Represents the particle vibration velocity vector wave field. Represents the stress tensor. Indicates the focal term. This represents the first derivative of the particle's vibrational velocity with respect to time. The divergence of the stress tensor, i.e., the body force term generated by stress changes, is represented by the wave field separation relation:

[0089] (2);

[0090] in, This represents the P-wave particle vibration velocity vector wave field. Let S represent the particle vibration velocity vector wave field, with superscripts P and S indicating P-wave mode and S-wave mode respectively. The equation indicates that the particle vibration velocity vector wave field is decomposed into the sum of P-wave component and S-wave component.

[0091] (3) Based on the spatial gradient of the separated wavefield and the background elastic parameter model, a non-uniform driving acceleration containing only the medium parameter gradient term is constructed to characterize the mode conversion drive caused by medium non-uniformity. Non-uniform driving acceleration for:

[0092] (3);

[0093] Among them, wave mode index Take the P wave or S wave. ; For spatial components or spatial direction indices. ; , , Indicates by The first acceleration response generated by the velocity vector wave field of separated particles One portion, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first Partial derivatives in spatial directions, express In the wave field of the vibration velocity vector of separated particles, the first The component is related to the first... Partial derivatives in each spatial direction, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first The partial derivatives in the spatial direction; the Lamé parameter and the spatial gradient of density are used to characterize the driving effect of medium inhomogeneity on wave mode conversion.

[0094] (4) Perform a Poisson projection on the non-uniform driving acceleration to obtain Coupled source and Coupled source; Poisson projection includes the following process:

[0095] right The coupling source first constructs the scalar curl of the P-wave non-uniform driving acceleration. Solve the stream function using the scalar curl as the source term. Used to recover the projected acceleration in the S-wave subspace:

[0096] (4);

[0097] in, The scalar curl representing the non-uniform driving acceleration of the P-wave is used to extract shear-type cross-mode information generated by the non-uniform driving acceleration of the incident P-wave and serves as the right-hand term of the subsequent Poisson projection to construct the coupling source of the P-wave to S-wave conversion. and These represent the non-uniform driving accelerations generated by the P-wave separated wavefield, respectively. direction and Components in direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the stream function obtained from the right-hand side of the P-wave curl, used for subsequent recovery of the S-wave subspace. Coupled source;

[0098] The stream function obtained from the Poisson equation is combined with its spatial derivatives to reconstruct the cross-mode projection components, thus obtaining... Coupled source :

[0099] (5);

[0100] in, This represents the S-wave subspace coupling source obtained by Poisson projection of P-wave non-uniform driving acceleration. The subscript S indicates projection onto the S-wave subspace, and the superscript P indicates that the coupling source is constructed by P-wave non-uniform driving acceleration. and Representing the stream function right direction and First-order partial derivative in the direction;

[0101] right Coupled source, first construct the divergence of S-wave non-uniform driving acceleration. Solve the scalar potential function using this divergence as the source term. Used to recover the projected acceleration in the P-wave subspace:

[0102] (6);

[0103] in, This represents the divergence of the non-uniform S-wave driven acceleration, which is used to extract the information of the compressed cross-wave mode generated by the incident non-uniform S-wave driven acceleration and serves as the right-hand term of the subsequent Poisson projection to construct the coupling source of the S-wave to P-wave conversion. and These represent the non-uniform driving accelerations generated by the S-wave separated wavefield in... Components in the z-direction and z-direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the scalar potential function obtained from the right-hand side of the S-wave divergence, used for subsequent recovery of the P-wave subspace. Coupled source;

[0104] Based on the potential function obtained from the Poisson equation, its spatial derivatives are combined to reconstruct the cross-mode projection components, thus obtaining... Coupled source :

[0105] (7);

[0106] in, This represents a P-wave subspace coupled source obtained by Poisson projection of S-wave non-uniform driven acceleration. The subscript P indicates projection onto the P-wave subspace, and the superscript S indicates that the coupled source is constructed by S-wave non-uniform driven acceleration. and They represent scalar potential functions respectively. right direction and The first-order partial derivative in the direction.

[0107] (5) Construct the converted power density and the same-mode power density based on the coupled source and the separated wave field, and perform spatial smoothing on the power density; the converted power density and the same-mode power density are defined as follows:

[0108] (8a);

[0109] (8b);

[0110] (8c);

[0111] (8d);

[0112] in, This is used to extract local positive power and avoid mutual cancellation of oscillation phases; in equation (8a), express Instantaneous conversion power density of the conversion channel The S-wave velocity vector wave field is represented by the symbol ·, which indicates the vector dot product; in equation (8b), Represents the S-wave co-mode reference power density. This represents the same-mode projection acceleration obtained by projecting the S-wave non-uniform driving acceleration onto the S-wave subspace, used as... The S-wave co-mode reference term when the conversion intensity is normalized; in equation (8c), express Instantaneous conversion power density of the conversion channel This expression represents the P-wave velocity vector wave field, obtained by relating the target P-wave velocity field to... The dot product of coupled sources measures locality. Transformation response intensity; in equation (8d), This represents the P-wave co-mode reference power density. This represents the same-mode projected acceleration obtained by projecting the P-wave non-uniform driving acceleration onto the P-wave subspace; this quantity is used as... P-wave co-mode reference term when the conversion intensity is normalized.

[0113] The spatial smoothing operator is:

[0114] (9);

[0115] in, Represents grid points Power density after smoothing Represented by grid points Unsmoothed power density in the neighborhood centered on the target; For grid indexing; To smooth the window half-width, the grid spacing can be selected based on the dominant wavelength, grid spacing, and noise level. Generally, 1 to 5 grid points are used, with 2 or 3 grid points being preferred. The window size is... ; and To smooth the inner edge of the window direction and Integer offset of direction, all values ​​are in the range of - to Both are merely summation indices and do not represent P / S wave modes; and They are respectively direction and The unit grid vector in the direction;

[0116] After smoothing, the result is , , and .

[0117] (6) Construct a material gating function based on the parameter gradient of the background elastic parameter model, and construct the conversion intensity ratio by combining the smoothed conversion power density and the same-mode power density; the material gating function approaches zero in the local uniform region and increases in the parameter abrupt region, which is used to limit the physically reasonable conversion position. The material gating function is constructed from the elastic parameter gradient. First, calculate the spatial gradient magnitude of the Lamé parameter and density, then normalize them with the corresponding parameter magnitudes, and then weight and superimpose the normalized gradient contributions and cut them to the range of 0 to 1, satisfying:

[0118] (10);

[0119] (11);

[0120] in, , and They represent , and Spatial gradient magnitudes of the three elastic parameters; , They represent right direction and Partial derivatives in direction, , They represent right direction and Partial derivatives in direction, , They represent right direction and The larger the partial derivative of the direction and the larger the gradient magnitude, the stronger the change in medium parameters at that location, and the more likely it is to correspond to the physical wave mode conversion location; Represents the material gating function; , and They represent , and The weighting coefficients of the gradient term are used to adjust the relative contributions of the Lamé parameter and density gradient to the material gating function, and can be taken as values ​​between 0 and 1; in the absence of special priors, they can be taken as equal weights, for example, all three can be taken as 0.33; , and This indicates the magnitude or absolute value of the corresponding elastic parameter; As a stability constant, it is used to avoid the normalization denominator from approaching zero. It is a positive number that is much smaller than the representative value of the corresponding elastic parameter, for example, one-thousandth of the maximum value of the corresponding parameter model, so as to avoid the denominator from approaching zero without changing the gating response at the real interface. and Used to restrict the gating function to the range of 0 to 1.

[0121] Equation (10) is used to calculate the spatial gradient magnitudes of the Lamé parameters and density, respectively. The larger the gradient magnitude, the more likely the location is to be an interface or a strongly non-uniform region, and the more likely it is to support physically reasonable wave mode conversion. Equation (11) is used to weightedly synthesize the normalized gradient contributions of each elastic parameter into a material gating function. The material gating function approaches 0 in locally uniform regions, increases in parameter abrupt changes or strongly non-uniform regions, and is restricted to between 0 and 1.

[0122] The conversion intensity ratio satisfies:

[0123] (12a);

[0124] (12b);

[0125] (12c);

[0126] (12d);

[0127] In equation (12a), This represents the normalized intensity ratio of the P-wave co-mode channel, which measures the relative proportion of the P-wave co-mode response in the P-wave correlated response. The smoothed P-wave co-mode reference power density, For smoothed Conversion power density, To prevent a positive stability constant with a denominator of zero, its value should be smaller than the characteristic order of the denominator; it can be taken as 10 times the maximum value of the corresponding denominator. -6 Up to 10 -3 10 are preferred -6 To ensure numerical stability in the low-energy region and without affecting the intensity ratio in the effective conversion region; in equation (12b), express The normalized conversion intensity ratio of the conversion channel; in equation (12c), This represents the normalized intensity ratio of the S-wave co-mode channel, which measures the relative proportion of the S-wave co-mode response in the S-wave correlated response. For S-wave co-mode reference power density, for Conversion power density; in equation (12d), express The normalized conversion intensity ratio of the conversion channel; the larger this ratio, the more supportive the local wavefield. Transform it into something like a contribution.

[0128] The conversion intensity ratio in equations (12a) to (12d) is a normalized conversion index used to characterize the ratio of the cross-mode coupling intensity to the same-mode reference intensity; it does not represent a strict physical probability, but is used to measure the degree of support of the local wave field for the corresponding conversion channel in imaging conditions.

[0129] (7) Construct spatiotemporal conversion sensing gating weights based on the material gating function and the conversion intensity ratio; the spatiotemporal conversion sensing gating weights satisfy:

[0130] (13);

[0131] in, This indicates that the source-detector mode combination is... Spatiotemporal transformation sensing gating weights, , Both are P or S; and These represent the corresponding conversion intensity ratios calculated from the source-end wavefield and the receiver-end wavefield, respectively. Indicates the source end, The terminator represents the detector end; this formula requires that the gating weights be increased only when the model parameter gradient, the source end wavefield conversion index, and the detector end wavefield conversion index are all supported simultaneously.

[0132] Equation (13) multiplies the material gating function with the conversion intensity ratio of the same imaging channel at the source and detector ends to form a conversion-aware gating weight that changes with both space and time; the imaging contribution is preserved or enhanced only when the model parameter gradient and the wave field conversion index at both ends support the channel.

[0133] (8) The spatiotemporal transformation sensing gating weights are introduced into the elastic reverse time migration imaging condition to weight the imaging contribution and accumulate it over time, obtaining one or more imaging results among PP, PS, SP, and SS. Elastic reverse time migration imaging condition and imaging results. satisfy:

[0134] (14);

[0135] (15);

[0136] in, express The passage in spatial location and time Cross-correlation imaging terms at the location; Indicates the detector end Separate vector wave field Indicates the source end Separate vector wave field; express The final imaging result of the channel, , , and These represent the imaging results for the PP, PS, SP, and SS channels, respectively.

[0137] Equation (14) represents the introduction of conversion sensing gating weights into the cross-correlation imaging conditions of the vector wave fields at the source and detector ends, weighting the candidate imaging contributions at each time step. The imaging results of Equation (15) are uniformly denoted as... This means that the cross-correlation imaging contributions at each time step are summed after gating and weighting to obtain the final imaging result.

[0138] When the background medium is locally homogeneous, the gradient of the medium parameters satisfies and Therefore, Corresponding Coupled source, The coupling source, conversion power density, and conversion intensity ratio are all zero to suppress non-physical conversion contributions in local uniform regions.

[0139] The feasibility and effectiveness of the present invention are further illustrated through three embodiments.

[0140] Example 1

[0141] Figure 2 It is a simple horizontal two-layer model, namely the P-wave velocity model and the S-wave velocity model. The model contains 401×401 grids, with a spatial interval of 10m in each direction and a time sampling interval of 1ms. An explosive source is used to generate seismic waves, and the wavelet used is the Rick wavelet with a dominant frequency of 25Hz. Figure 3 It uses the finite difference method to... Figure 2 The mixed-source seismic record obtained in the simple horizontal two-layer model shown is the X component and Z component of the seismic record obtained using the finite difference method; this seismic record will be used for subsequent migration imaging. Figure 4 (a)-(j) in the present invention are derived from Figure 2 The intermediate quantities of the incident P-wave at the source end obtained in the simple horizontal two-layer model shown can reflect the location and intensity of the wave transition at that moment, thereby providing a gating function. Figure 5 It is by utilizing the present invention from Figure 2 The longitudinal wave source migration imaging results obtained in the simple horizontal two-layer model shown are the PP imaging results and the PS imaging results, respectively. It can be seen that the imaging results obtained by the method of the present invention significantly suppress crosstalk artifacts and enhance the continuity of the in-phase axis, proving the effectiveness of the present invention.

[0142] Example 2

[0143] Figure 6 The Overthrust model provided by this invention includes the P-wave velocity model. Shear wave velocity model The model contains 1201×413 grids with a spatial spacing of 10m in each direction. It uses an explosion source to excite seismic waves, and the wavelet used is the Rick wavelet with a dominant frequency of 25Hz. The time sampling interval is 1ms. Figure 7 (a)-(f) in the present invention are derived from Figure 6 The relevant intermediate quantities obtained from the Overthrust model shown accurately reflect the location and intensity of the wave transition at that moment. Figure 8 It is by utilizing the present invention from Figure 6 The P-wave source migration imaging results obtained in the Overthrust model shown are the PP imaging result and the PS imaging result, respectively. It can be seen that for complex models, the elastic reverse time migration imaging method proposed in this invention can still obtain high-precision imaging results. The above results verify the effectiveness of the method in complex models.

[0144] Example 3

[0145] Figure 9 The present invention provides actual data migration velocity models, namely, P-wave velocity models. Shear wave velocity model The spatial interval in each direction is 10m. Seismic waves are excited by an explosive source. The wavelet used is the Rick wavelet with a main frequency of 35Hz. The time sampling interval is 1ms. Figure 10 It is by utilizing the present invention from Figure 9 The P-wave source migration imaging results obtained from the actual data velocity model are shown below, representing PP and PS imaging results, respectively. It can be seen that the elastic reverse-time migration imaging method proposed in this invention can still obtain high-precision imaging results for actual data. These results verify the effectiveness of the method in actual data.

[0146] In summary, the method of the present invention has good feasibility and practicality in complex geological and geophysical models.

[0147] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for elastic reverse-time migration imaging based on conversion sensing gating, characterized in that, Includes the following steps: (1) Obtain multi-component earthquake records, source information, and background elastic parameter model, wherein the background elastic parameter model includes at least Lamé parameters. , and density ; (2) Based on the elastic wave equation, the source wave field is propagated in the forward direction and the received wave field is propagated in the reverse time direction. The source wave field and the received wave field are separated into P-wave and S-wave respectively to obtain the separated source wave field and the received wave field. (3) Based on the spatial gradient of the separated wave field and the background elastic parameter model, construct a non-uniform driving acceleration that only contains the gradient term of the medium parameter; (4) Perform a Poisson projection on the non-uniform driving acceleration to obtain Coupled source and Coupled source; (5) Construct the converted power density and the same-mode power density based on the coupled source and the separated wave field, and perform spatial smoothing on the power density; (6) Construct a material gating function based on the parameter gradient of the background elastic parameter model, and construct the conversion intensity ratio by combining the smoothed conversion power density and the same-mode power density; (7) Construct spatiotemporal transformation sensing gating weights based on the material gating function and the transformation intensity ratio; (8) Introduce the spatiotemporal conversion sensing gating weight into the elastic reverse time-shift imaging condition, weight the imaging contribution and accumulate it over time to obtain one or more imaging results among PP, PS, SP and SS.

2. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, The elastic wave equation used in step (2) satisfies: (1)； in, Indicates spatial location Density at that location Represents a spatial position vector. Indicates time, Represents the particle vibration velocity vector wave field. Represents the stress tensor. Indicates the focal term. This represents the first derivative of the particle's vibrational velocity with respect to time. It represents the divergence of the stress tensor, i.e., the body force term generated by stress changes; The corresponding wavefield separation relationship is: (2)； in, This represents the P-wave particle vibration velocity vector wave field. Let S represent the particle vibration velocity vector wave field, with superscripts P and S indicating P-wave mode and S-wave mode respectively. The equation indicates that the particle vibration velocity vector wave field is decomposed into the sum of P-wave component and S-wave component.

3. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, Non-uniform driving acceleration in step (3) for: (3)； Among them, wave mode index Take the P wave or S wave. ; For spatial components or spatial direction indices. ; , , Indicates by The first acceleration response generated by the velocity vector wave field of separated particles One portion, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first Partial derivatives in spatial directions, express In the wave field of the vibration velocity vector of separated particles, the first The component is related to the first... Partial derivatives in each spatial direction, express In the wave field of the vibration velocity vector of separated particles, the first Components for the first The partial derivatives in the spatial direction; the Lamé parameter and the spatial gradient of density are used to characterize the driving effect of medium inhomogeneity on wave mode conversion.

4. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, Step (4) of the Poisson projection includes the following process: right The coupling source first constructs the scalar curl of the P-wave non-uniform driving acceleration. Solve the stream function using the scalar curl as the source term. : (4)； in, Scalar curl representing the non-uniform driving acceleration of the P-wave; and These represent the non-uniform driving accelerations generated by the P-wave separated wavefield, respectively. direction and Components in direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the stream function obtained from the right-hand side of the P-wave curl, used for subsequent recovery of the S-wave subspace. Coupled source; The stream function obtained from the Poisson equation is combined with its spatial derivatives to reconstruct the cross-mode projection components, thus obtaining... Coupled source : (5)； in, This represents the S-wave subspace coupling source obtained by Poisson projection of P-wave non-uniform driving acceleration. The subscript S indicates projection onto the S-wave subspace, and the superscript P indicates that the coupling source is constructed by P-wave non-uniform driving acceleration. and Representing the stream function right direction and First-order partial derivative in the direction; right Coupled source, first construct the divergence of S-wave non-uniform driving acceleration. Solve the scalar potential function using this divergence as the source term. : (6)； in, This represents the divergence of the S-wave non-uniform driving acceleration. and These represent the non-uniform driving accelerations generated by the S-wave separated wavefield in... direction and Components in direction; and These represent the spatial partial derivatives of the corresponding acceleration components; Represents the two-dimensional Laplace operator; This represents the scalar potential function obtained from the right-hand side of the S-wave divergence, used for subsequent recovery of the P-wave subspace. Coupled source; Based on the potential function obtained from the Poisson equation, its spatial derivatives are combined to reconstruct the cross-mode projection components, thus obtaining... Coupled source : (7)； in, This represents a P-wave subspace coupled source obtained by Poisson projection of S-wave non-uniform driven acceleration. The subscript P indicates projection onto the P-wave subspace, and the superscript S indicates that the coupled source is constructed by S-wave non-uniform driven acceleration. and They represent scalar potential functions respectively. right The first-order partial derivatives in the z-direction and z-direction.

5. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, In step (5), the conversion power density and the same-mode power density are defined as follows: (8a)； (8b)； (8c)； (8d)； in, This is used to extract local positive power and avoid mutual cancellation of oscillation phases; express Instantaneous conversion power density of the conversion channel; Represents the S-wave velocity vector wave field; This represents the S-wave co-mode reference power density; This represents the same-mode projection acceleration obtained by projecting the S-wave non-uniform driving acceleration onto the S-wave subspace, used as... S-wave co-mode reference term when the conversion intensity is normalized; express Instantaneous conversion power density of the conversion channel; Represents the P-wave velocity vector wave field; This represents the P-wave co-mode reference power density. This represents the same-mode projected acceleration obtained by projecting the P-wave non-uniform driving acceleration onto the P-wave subspace; this quantity is used as... P-wave co-mode reference term when the conversion intensity is normalized.

6. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1 or 5, characterized in that, The spatial smoothing operator in step (5) is: (9)； in, Represents grid points Power density after smoothing Represented by grid points Unsmoothed power density in the neighborhood centered on the target; For grid indexing; To smooth the window's half-width, the window size is... ; and To smooth the inner edge of the window direction and Integer offset of direction, all values ​​are in the range of - to Both are merely summation indices and do not represent P / S wave modes; and They are respectively direction and The unit grid vector in the direction; After smoothing, the result is , , and .

7. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, In step (6), the material gating function is constructed from the elastic parameter gradient. First, the spatial gradient magnitudes of the Lamé parameter and density are calculated, and then normalized using the corresponding parameter magnitudes. The weighted summation of the normalized gradient contributions is then truncated to the range of 0 to 1, satisfying the following: (10)； (11)； in, , and They represent , and Spatial gradient magnitudes of the three elastic parameters; , They represent right direction and Partial derivatives in direction, , They represent right direction and Partial derivatives in direction, , They represent right direction and Partial derivatives in direction; Represents the material gating function; , and They represent , and The weighting coefficients of the gradient term; , and This indicates the magnitude or absolute value of the corresponding elastic parameter; It is a stability constant used to prevent the normalized denominator from approaching zero; and Used to restrict the gating function to the range of 0 to 1.

8. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, The conversion intensity ratio in step (6) satisfies: (12a)； (12b)； (12c)； (12d)； in, This represents the normalized intensity ratio of the P-wave co-mode channel; The smoothed P-wave co-mode reference power density, For smoothed Conversion power density, To prevent positive stability constants with a denominator of zero; express Normalized conversion intensity ratio of the conversion channel; This represents the normalized intensity ratio of the S-wave co-mode channel; For S-wave co-mode reference power density, for Conversion power density; express Normalized conversion intensity ratio of the conversion channel.

9. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, The spatiotemporal transformation sensing gating weights in step (7) satisfy: (13)； in, This indicates that the source-detector mode combination is... Spatiotemporal transformation sensing gating weights, , Both are P or S; and These represent the corresponding conversion intensity ratios calculated from the source-end wavefield and the receiver-end wavefield, respectively. Indicates the source end, Indicates the detector end.

10. The elastic reverse-time migration imaging method based on conversion sensing gating as described in claim 1, characterized in that, Step (8) Elastic reverse time migration imaging conditions and imaging results satisfy: (14)； (15)； in, express The passage in spatial location and time Cross-correlation imaging terms at the location; Indicates the detector end Separate vector wave field Indicates the source end Separate vector wave field; express The final imaging result of the channel, , , and These represent the imaging results for the PP, PS, SP, and SS channels, respectively.

Images